From 2a0c42c0878d21dcb644836252bae3cd751d42c1 Mon Sep 17 00:00:00 2001 From: Schellenberger Date: Fri, 7 Aug 2015 11:39:46 +0200 Subject: [PATCH] =?UTF-8?q?Update=20of=20RK=20iteration=20Updated=20to=20n?= =?UTF-8?q?ew=20scheme=20after=20R=C3=B6=C3=9Fler2010=20No=20more=20macro?= =?UTF-8?q?=20shenanigans=20Replaced=20pow(*,=20N)=20function=20with=20N?= =?UTF-8?q?=20times=20product?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Cortical_Column.cpp | 117 ++++++++++++++++------------------ Cortical_Column.h | 111 ++++++++++++++++++--------------- ODE.h | 18 +----- Plots.m | 34 ++++------ TC.cpp | 33 +++++----- TC_model.pro | 2 +- Test_Stimulation.m | 176 +++++++++++++++++++++++++--------------------------- Thalamic_Column.cpp | 176 ++++++++++++++++++++++------------------------------ Thalamic_Column.h | 164 +++++++++++++++++++++++++----------------------- macros.h | 83 ------------------------- saves.h | 24 +++---- 11 files changed, 405 insertions(+), 533 deletions(-) delete mode 100644 macros.h diff --git a/Cortical_Column.cpp b/Cortical_Column.cpp index 947eb6d..5443cff 100644 --- a/Cortical_Column.cpp +++ b/Cortical_Column.cpp @@ -29,16 +29,21 @@ /* Initialization of RNG */ /****************************************************************************************************/ void Cortical_Column::set_RNG(void) { - /* Number of independent streams */ - int N = 4; + extern const double dt; + /* Number of independent random variables */ + int N = 2; /* Create RNG for each stream */ for (int i=0; i B1, B2; - double n = 1 / h * (B1[N-1] * Rand_vars[2*M] + B2[N-1] * Rand_vars[2*M+1]); - return n; + return gamma_e * gamma_e * (Rand_vars[2*M] + Rand_vars[2*M+1]/std::sqrt(3))*B[N]; +} + +double Cortical_Column::noise_aRK(int M) const{ + return gamma_e * gamma_e * (Rand_vars[2*M] - Rand_vars[2*M+1]*std::sqrt(3))/4; } /****************************************************************************************************/ /* end */ @@ -162,21 +158,19 @@ double Cortical_Column::noise_xRK(int N, int M) const{ /****************************************************************************************************/ void Cortical_Column::set_RK (int N) { extern const double dt; - _SWITCH((Phi_ee)(Phi_ei)(Phi_ie)(Phi_ii)(phi) - (x_ee) (x_ei) (x_ie) (x_ii) (y)) - Ve [N] = dt*(-(I_L_e(N) + I_ee(N) + I_ie(N))/tau_e - I_KNa(N)); - Vi [N] = dt*(-(I_L_i(N) + I_ei(N) + I_ii(N))/tau_i); - Na [N] = dt*(alpha_Na * get_Qe(N) - Na_pump(N))/tau_Na; - Phi_ee [N] = dt*(var_x_ee); - Phi_ei [N] = dt*(var_x_ei); - Phi_ie [N] = dt*(var_x_ie); - Phi_ii [N] = dt*(var_x_ii); - phi [N] = dt*(var_y); - x_ee [N] = dt*(pow(gamma_e, 2) * (N_ee * get_Qe(N) + noise_xRK(N, 0) + N_te * Thalamus->get_phi(N) - var_Phi_ee) - 2 * gamma_e * var_x_ee); - x_ei [N] = dt*(pow(gamma_e, 2) * (N_ei * get_Qe(N) + noise_xRK(N, 1) + N_ti * Thalamus->get_phi(N) - var_Phi_ei) - 2 * gamma_e * var_x_ei); - x_ie [N] = dt*(pow(gamma_i, 2) * (N_ie * get_Qi(N) - var_Phi_ie) - 2 * gamma_i * var_x_ie); - x_ii [N] = dt*(pow(gamma_i, 2) * (N_ii * get_Qi(N) - var_Phi_ii) - 2 * gamma_i * var_x_ii); - y [N] = dt*(pow(nu, 2) * ( get_Qe(N) - var_phi) - 2 * nu * var_y); + Ve [N+1] = Ve [0] + A[N] * dt*(-(I_L_e(N) + I_ee(N) + I_ie(N))/tau_e - I_KNa(N)); + Vi [N+1] = Vi [0] + A[N] * dt*(-(I_L_i(N) + I_ei(N) + I_ii(N))/tau_i); + Na [N+1] = Na [0] + A[N] * dt*(alpha_Na * get_Qe(N) - Na_pump(N))/tau_Na; + y_ee[N+1] = y_ee[0] + A[N] * dt*(x_ee[N]); + y_ei[N+1] = y_ei[0] + A[N] * dt*(x_ei[N]); + y_ie[N+1] = y_ie[0] + A[N] * dt*(x_ie[N]); + y_ii[N+1] = y_ii[0] + A[N] * dt*(x_ii[N]); + y [N+1] = y [0] + A[N] * dt*(x [N]); + x_ee[N+1] = x_ee[0] + A[N] * dt*(pow(gamma_e, 2) * (N_ee * get_Qe(N) + N_te * Thalamus->y[N] - y_ee[N]) - 2 * gamma_e * x_ee[N]) + noise_xRK(N, 0); + x_ei[N+1] = x_ei[0] + A[N] * dt*(pow(gamma_e, 2) * (N_ei * get_Qe(N) + N_ti * Thalamus->y[N] - y_ei[N]) - 2 * gamma_e * x_ei[N]) + noise_xRK(N, 1) ; + x_ie[N+1] = x_ie[0] + A[N] * dt*(pow(gamma_i, 2) * (N_ie * get_Qi(N) - y_ie[N]) - 2 * gamma_i * x_ie[N]); + x_ii[N+1] = x_ii[0] + A[N] * dt*(pow(gamma_i, 2) * (N_ii * get_Qi(N) - y_ii[N]) - 2 * gamma_i * x_ii[N]); + x [N+1] = x [0] + A[N] * dt*(pow(nu, 2) * ( get_Qe(N) - y [N]) - 2 * nu * x [N]); } /****************************************************************************************************/ /* end */ @@ -187,24 +181,23 @@ void Cortical_Column::set_RK (int N) { /* Function that adds all SRK terms */ /****************************************************************************************************/ void Cortical_Column::add_RK(void) { - extern const double h; - Ve [0] += (Ve [1] + Ve [2] * 2 + Ve [3] * 2 + Ve [4])/6; - Vi [0] += (Vi [1] + Vi [2] * 2 + Vi [3] * 2 + Vi [4])/6; - Na [0] += (Na [1] + Na [2] * 2 + Na [3] * 2 + Na [4])/6; - Phi_ee [0] += (Phi_ee [1] + Phi_ee[2] * 2 + Phi_ee[3] * 2 + Phi_ee[4])/6; - Phi_ei [0] += (Phi_ei [1] + Phi_ei[2] * 2 + Phi_ei[3] * 2 + Phi_ei[4])/6; - Phi_ie [0] += (Phi_ie [1] + Phi_ie[2] * 2 + Phi_ie[3] * 2 + Phi_ie[4])/6; - Phi_ii [0] += (Phi_ii [1] + Phi_ii[2] * 2 + Phi_ii[3] * 2 + Phi_ii[4])/6; - phi [0] += (phi [1] + phi [2] * 2 + phi [3] * 2 + phi [4])/6; - x_ee [0] += (x_ee [1] + x_ee [2] * 2 + x_ee [3] * 2 + x_ee [4])/6 + pow(gamma_e, 2) * h * Rand_vars[0]; - x_ei [0] += (x_ei [1] + x_ei [2] * 2 + x_ei [3] * 2 + x_ei [4])/6 + pow(gamma_e, 2) * h * Rand_vars[2]; - x_ie [0] += (x_ie [1] + x_ie [2] * 2 + x_ie [3] * 2 + x_ie [4])/6; - x_ii [0] += (x_ii [1] + x_ii [2] * 2 + x_ii [3] * 2 + x_ii [4])/6; - y [0] += (y [1] + y [2] * 2 + y [3] * 2 + y [4])/6; - - /* Generat noise for the next iteration */ + Ve [0] = (-3*Ve [0] + 2*Ve [1] + 4*Ve [2] + 2*Ve [3] + Ve [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; + y_ee[0] = (-3*y_ee[0] + 2*y_ee[1] + 4*y_ee[2] + 2*y_ee[3] + y_ee[4])/6; + y_ei[0] = (-3*y_ei[0] + 2*y_ei[1] + 4*y_ei[2] + 2*y_ei[3] + y_ei[4])/6; + y_ie[0] = (-3*y_ie[0] + 2*y_ie[1] + 4*y_ie[2] + 2*y_ie[3] + y_ie[4])/6; + y_ii[0] = (-3*y_ii[0] + 2*y_ii[1] + 4*y_ii[2] + 2*y_ii[3] + y_ii[4])/6; + y [0] = (-3*y [0] + 2*y [1] + 4*y [2] + 2*y [3] + y [4])/6; + x_ee[0] = (-3*x_ee[0] + 2*x_ee[1] + 4*x_ee[2] + 2*x_ee[3] + x_ee[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_ie[0] = (-3*x_ie[0] + 2*x_ie[1] + 4*x_ie[2] + 2*x_ie[3] + x_ie[4])/6; + x_ii[0] = (-3*x_ii[0] + 2*x_ii[1] + 4*x_ii[2] + 2*x_ii[3] + x_ii[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 "Thalamic_Column.h" -#include "macros.h" using std::vector; class Thalamic_Column; @@ -46,6 +45,17 @@ typedef boost::variate_generator GEN; /* Variate generator */ /****************************************************************************************************/ +/* 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 { @@ -62,56 +72,55 @@ class Cortical_Column { /* Connect to the thalamic module */ void get_Thalamus(Thalamic_Column& T) {Thalamus = &T;} - /* Return axonal flux */ - double get_phi (int N) const {_SWITCH((phi)); return var_phi;} - /* Data storage access */ - friend void get_data (int, Cortical_Column&, Thalamic_Column&, _REPEAT(double*, 4)); + friend void get_data (int, Cortical_Column&, Thalamic_Column&, vector); - /* ODE functions */ - void set_RK (int); - void add_RK (void); + /* ODE functions */ + void set_RK (int); + void add_RK (void); /* Stimulation protocol access */ - friend class Stim; + friend class Stim; + friend class Thalamic_Column; private: - /* Initialize the RNGs */ - void set_RNG (void); - - /* Firing rates */ - double get_Qe (int) const; - double get_Qi (int) const; - - /* Currents */ - double I_ee (int) const; - double I_ei (int) const; - double I_ie (int) const; - double I_ii (int) const; - double I_L_e (int) const; - double I_L_i (int) const; - double I_KNa (int) const; - - /* Potassium pump */ - double Na_pump (int) const; - - /* Noise function */ - double noise_xRK (int, int) const; - - /* Population variables */ + /* Initialize the RNGs */ + void set_RNG (void); + + /* Firing rates */ + double get_Qe (int) const; + double get_Qi (int) const; + + /* Currents */ + double I_ee (int) const; + double I_ei (int) const; + double I_ie (int) const; + double I_ii (int) const; + double I_L_e (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; + + /* Population variables */ vector Ve = _INIT(E_L_e), /* excitatory membrane voltage */ Vi = _INIT(E_L_i), /* inhibitory membrane voltage */ Na = _INIT(Na_eq), /* Na concentration */ - Phi_ee = _INIT(0.0), /* PostSP from excitatory to excitatory population */ - Phi_ei = _INIT(0.0), /* PostSP from excitatory to inhibitory population */ - Phi_ie = _INIT(0.0), /* PostSP from inhibitory to excitatory population */ - Phi_ii = _INIT(0.0), /* PostSP from inhibitory to inhibitory population */ - phi = _INIT(0.0), /* axonal flux */ - x_ee = _INIT(0.0), /* derivative of Phi_ee */ - x_ei = _INIT(0.0), /* derivative of Phi_ei */ - x_ie = _INIT(0.0), /* derivative of Phi_ie */ - x_ii = _INIT(0.0), /* derivative of Phi_ii */ - y = _INIT(0.0); /* derivative of phi */ + y_ee = _INIT(0.0), /* PostSP from excitatory to excitatory population */ + y_ei = _INIT(0.0), /* PostSP from excitatory to inhibitory population */ + y_ie = _INIT(0.0), /* PostSP from inhibitory to excitatory population */ + y_ii = _INIT(0.0), /* PostSP from inhibitory to inhibitory population */ + y = _INIT(0.0), /* axonal flux */ + x_ee = _INIT(0.0), /* derivative of y_ee */ + x_ei = _INIT(0.0), /* derivative of y_ei */ + x_ie = _INIT(0.0), /* derivative of y_ie */ + x_ii = _INIT(0.0), /* derivative of y_ii */ + x = _INIT(0.0); /* derivative of y */ /* Random number generators */ vector MTRands; @@ -141,9 +150,9 @@ class Cortical_Column { /* 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 tau_Na = 1.7; /* Sodium time constant in ms */ - const double R_pump = 0.09; /* Na-K pump constant in mM/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 */ @@ -166,7 +175,7 @@ class Cortical_Column { const double E_GABA = -70; /* Leak */ - const double E_L_e = -64; + const double E_L_e = -64; const double E_L_i = -64; /* Potassium */ @@ -174,19 +183,23 @@ class Cortical_Column { /* Noise parameters in ms^-1 */ const double mphi = 0E-3; - const double dphi = 60E-3; + const double dphi = 20E-1; double input = 0.0; /* Connectivities (dimensionless) */ - const double N_ee = 115; + const double N_ee = 115; const double N_ei = 72; - const double N_ie = 90; - const double N_ii = 90; + const double N_ie = 90; + const double N_ii = 90; const double N_te = 2.5; const double N_ti = 2.5; /* Pointer to thalamic column */ Thalamic_Column* Thalamus; + + /* Interation matrix for SRK4 */ + const vector A = {0.5, 0.5, 1.0, 1.0}; + const vector B = {0.75, 0.75, 0.0, 0.0}; }; /****************************************************************************************************/ /* end */ diff --git a/ODE.h b/ODE.h index 44f647c..a659131 100644 --- a/ODE.h +++ b/ODE.h @@ -32,7 +32,7 @@ /****************************************************************************************************/ inline 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=1; i<=4; ++i) { + for (int i=0; i<4; ++i) { Cortex.set_RK(i); Thalamus.set_RK(i); } @@ -44,19 +44,3 @@ inline void ODE(Cortical_Column& Cortex, Thalamic_Column& Thalamus) { /****************************************************************************************************/ /* end */ /****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Parameters for SRK4 integration */ -/****************************************************************************************************/ -extern const vector B1 = {0, - 0.626708569400000081728308032325, - 1.7296310295000001389098542858846, - 1.2703689705000000831347506391467}; -extern const vector B2 = {0, - 0.78000033203198970710445792065002, - 1.28727807507536762265942797967, - 0.44477273249350995909523476257164}; -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ diff --git a/Plots.m b/Plots.m index 0fedd7e..c43fd57 100644 --- a/Plots.m +++ b/Plots.m @@ -6,20 +6,22 @@ function Plots(type) type = 2; end +%mex CXXFLAGS="\$CXXFLAGS -std=c++11 -O3 -lgopm" TC.cpp Cortical_Column.cpp Thalamic_Column.cpp + if type == 1 Param_Cortex = [4.7; % sigma_e 1.33; % g_KNa - 120E-3]; % dphi + 20E-1]; % dphi - Param_Thalamus = [0.051; % g_h + Param_Thalamus = [0.048; % g_h 0.024]; % g_LK else - Param_Cortex = [6; % sigma_e - 2.05; % g_KNa - 120E-3]; % dphi + Param_Cortex = [6; % sigma_e + 2.05; % g_KNa + 20E-1]; % dphi - Param_Thalamus = [0.051; % g_h - 0.020]; % g_LK + Param_Thalamus = [0.048; % g_h + 0.02]; % g_LK end Connectivity = [ 2.6; % N_et @@ -49,7 +51,7 @@ function Plots(type) L = length(Vt); timeaxis = linspace(0,T,L); -figure(1) +figure(2) subplot(411), plot(timeaxis,Ve) title('Pyramidal membrane voltage'), xlabel('Time in s'), @@ -67,12 +69,9 @@ function Plots(type) xlabel('Time in s'), ylabel('V_{t} in mV') xlim([0,30]); -ylim([-70,-35]); +ylim([-80,-35]); % vertical line for markers -for i=1:var_stim(6) - hx = graph2d.constantline(Marker_Stim/1E2+(i-1)*var_stim(7)/1E3,'ydata', get(gca,'ylim'),'xdata', get(gca,'xlim'), 'color', 'black', 'LineStyle', ':'); - changedependvar(hx,'x'); -end + subplot(413), plot(timeaxis,ah) title('Thalamic relay I_h activation'), @@ -81,10 +80,7 @@ function Plots(type) xlim([0,30]); ylim(get(gca,'ylim')); % vertical line for markers -for i=1:var_stim(6) - hx = graph2d.constantline(Marker_Stim/1E2+(i-1)*var_stim(7)/1E3,'ydata', get(gca,'ylim'),'xdata', get(gca,'xlim'), 'color', 'black', 'LineStyle', ':'); - changedependvar(hx,'x'); -end + subplot(414), plot(timeaxis,Ca) title('Thalamic relay ca concentration'), @@ -93,10 +89,6 @@ function Plots(type) xlim([0,30]); ylim(get(gca,'ylim')); % vertical line for markers -for i=1:var_stim(6) - hx = graph2d.constantline(Marker_Stim/1E2+(i-1)*var_stim(7)/1E3,'ydata', get(gca,'ylim'),'xdata', get(gca,'xlim'), 'color', 'black', 'LineStyle', ':'); - changedependvar(hx,'x'); -end % [Pxx,f] = pwelch(Ve-mean(Ve),hamming(10*L/T), 2*L/T, [], L/T); % n = find(f<=30, 1, 'last' ); diff --git a/TC.cpp b/TC.cpp index 99cbda4..3c0c825 100644 --- a/TC.cpp +++ b/TC.cpp @@ -65,8 +65,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double* var_stim = mxGetPr (prhs[4]); /* Parameters of stimulation protocol */ /* Initialize the populations */ - Cortical_Column Cortex (Param_Cortex, Connections); - Thalamic_Column Thalamus (Param_Thalamus, Connections); + Cortical_Column Cortex = Cortical_Column(Param_Cortex, Connections); + Thalamic_Column Thalamus = Thalamic_Column(Param_Thalamus, Connections); /* Link both modules */ Cortex.get_Thalamus(Thalamus); @@ -76,16 +76,16 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { Stim Stimulation(Cortex, Thalamus, var_stim); /* Create data containers */ - mxArray* Ve = SetMexArray(1, T*res/red); - mxArray* Vt = SetMexArray(1, T*res/red); - mxArray* Ca = SetMexArray(1, T*res/red); - mxArray* ah = SetMexArray(1, T*res/red); + 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 */ - double* Pr_Ve = mxGetPr(Ve); - double* Pr_Vt = mxGetPr(Vt); - double* Pr_Ca = mxGetPr(Ca); - double* Pr_ah = mxGetPr(ah); + vector pData(Data.size(), NULL); + for(unsigned i=0; i=onset*res && t%red==0){ - get_data(count, Cortex, Thalamus, Pr_Ve, Pr_Vt, Pr_Ca, Pr_ah); + get_data(count, Cortex, Thalamus, pData); ++count; } } /* Output of the simulation */ - plhs[0] = Ve; - plhs[1] = Vt; - plhs[2] = Ca; - plhs[3] = ah; - plhs[4] = Stimulation.get_marker(); + /* Return the data containers */ + for(unsigned i=0; i + +Option_Data_ERP = {Data_Range_ERP; + -80:40:40; + -80:40:40; + 'black'; + xRange; + [xRange(1),xRange(end)]}'; + +Option_Model_FSP = {Model_Range_FSP; + linspace(Model_Range_FSP(1), Model_Range_FSP(2), 4); + linspace(Model_Range_FSP(1), Model_Range_FSP(2), 4); + 'black'; + xRange; + [xRange(1),xRange(end)]}'; + +Option_Data_FSP = {Data_Range_FSP; + linspace(Data_Range_FSP(1), Data_Range_FSP(2), 4); + linspace(Data_Range_FSP(1), Data_Range_FSP(2), 4); + 'black'; + xRange; + [xRange(1),xRange(end)]}'; + +[Ve, Vt, Ca, ah, Marker_Stim] = TC(T, Param_Cortex, Param_Thalamus, Connectivity, var_stim); Fs = length(Ve)/T; Ve_FSP = ft_preproc_hilbert(ft_preproc_bandpassfilter(Ve, Fs, [12, 15], 513, 'fir'), 'abs').^2; +xRange = [-1, 3]; % Search for peaks x_SO = Marker_Stim; @@ -93,14 +105,14 @@ function Test_Stimulation(type) % Set the variables N_Stim = length(x_SO); -time_event = linspace(xRange(1), xRange(end), (xRange(end)-xRange(1))*Fs); +time_event = linspace(xRange(1), xRange(end), (xRange(end)-xRange(1))*Fs+1); Events = zeros(length(time_event), N_Stim); Events_FSP = zeros(length(time_event), N_Stim); % Segmentation for i=1:N_Stim - Events(:,i) = Ve ((x_SO(i)+xRange(1)*Fs)+1:(x_SO(i)+xRange(end)*Fs)); - Events_FSP(:,i) = Ve_FSP((x_SO(i)+xRange(1)*Fs)+1:(x_SO(i)+xRange(end)*Fs)); + Events(:,i) = Ve ((x_SO(i)+xRange(1)*Fs):(x_SO(i)+xRange(end)*Fs)); + Events_FSP(:,i) = Ve_FSP((x_SO(i)+xRange(1)*Fs):(x_SO(i)+xRange(end)*Fs)); end mean_ERP_model= mean(Events, 2); %#ok<*NASGU> @@ -111,67 +123,45 @@ function Test_Stimulation(type) % Define handle for plotting BL_model =@(y,x) boundedline(y,x(:,1), x(:,2), 'alpha', 'transparency', 0.1, 'r'); BL_data =@(y,x) boundedline(y,x(:,1), x(:,2), 'alpha', 'transparency', 0.1, 'black'); - -% Option array for set -Option_Name = { 'ylim'; - 'ytick'; - 'yticklabel'; - 'ycolor'; - 'xtick'}'; - -Option_Model_ERP = {Model_xRange; - linspace(Model_xRange(1), Model_xRange(2), 5); - linspace(Model_xRange(1), Model_xRange(2), 5); - 'black'; - xRange}'; %#ok<*NBRAK> - -Option_Data_ERP = {Data_xRange; - linspace(Data_xRange(1), Data_xRange(2), 5); - linspace(Data_xRange(1), Data_xRange(2), 5); - 'black'; - xRange}'; - -Option_Model_FSP = {Model_Range_FSP; - linspace(Model_Range_FSP(1), Model_Range_FSP(2), 5); - linspace(Model_Range_FSP(1), Model_Range_FSP(2), 5); - 'black'; - xRange}'; - -Option_Data_FSP = {Data_Range_FSP; - linspace(Data_Range_FSP(1), Data_Range_FSP(2), 5); - linspace(Data_Range_FSP(1), Data_Range_FSP(2), 5); - 'black'; - xRange}'; - + figure(1) subplot(411) plot(linspace(0,30,3000),Ve(101:3100)); title(['Ve with a mean of :',num2str(mean(Ve))]); subplot(412) -plot(linspace(0,30,3000),Vi(101:3100)); -title(['Vi with a mean of :',num2str(mean(Vi))]); -subplot(413) plot(linspace(0,30,3000),Vt(101:3100)); title(['Vt with a mean of :',num2str(mean(Vt))]); +subplot(413) +plot(linspace(0,30,3000),Ca(101:3100)); +title(['Ca with a mean of :',num2str(mean(Ca))]); subplot(414) -plot(linspace(0,30,3000),Vr(101:3100)); -title(['Vr with a mean of :',num2str(mean(Vr))]); +plot(linspace(0,30,3000),ah(101:3100)); +title(['ah with a mean of :',num2str(mean(ah))]); % Create figure figure(2) clf subplot(211) -[AX1, ~, ~] = plotyy(time_events,[mean_ERP_data, sd_ERP_data],time_events,[mean_ERP_model, sd_ERP_model], BL_data, BL_model); -set(AX1(1), Option_Name, Option_Data_ERP); -set(AX1(2), Option_Name, Option_Model_ERP); +[AX1, ~, ~] = plotyy(time_events,[mean_ERP, sem_FSP],time_events,[mean_ERP_model, sd_ERP_model], BL_data, BL_model); +set(AX1(1),Option_Name, Option_Data_ERP, 'box', 'off'); +set(AX1(2),Option_Name, Option_Model_ERP); ylabel(AX1(1),'EEG [$\mu$V]'); -ylabel(AX1(2),'$V_{e}$ [mV]'); -title([num2str(N_Stim), ' Events']) +ylabel(AX1(2),'V$_{p}$ [mV]'); + subplot(212) -[AX2, ~, ~] = plotyy(time_events,[mean_FSP_data, sd_FSP_data],time_events,[mean_FSP_model, sd_FSP_model], BL_data, BL_model); -set(AX2(1), Option_Name, Option_Data_FSP); -set(AX2(2), Option_Name, Option_Model_FSP); -ylabel(AX2(1),'FSP data [a.u.]'); -ylabel(AX2(2),'$FSP model$ [a.u.]'); +[AX2, ~, ~] = plotyy(time_events,[mean_FSP, sem_FSP],time_events,[mean_FSP_model, sd_FSP_model], BL_data, BL_model); +set(AX2(1),Option_Name, Option_Data_FSP); +set(AX2(2),Option_Name, Option_Model_FSP); +ylabel(AX2(1),'Spindle Power [$\mu$V$^{2}$]'); +ylabel(AX2(2),'Spindle Power [mV$^{2}$]'); title([num2str(N_Stim), ' Events']) + +% Marker for stimulation +for i=1:2 + hx1 = graph2d.constantline((i-1)*1.05+0.125*(i-1)*(i-2)/2,'ydata', get(AX1(1),'ylim'), 'parent', AX1(1), 'color', 'black', 'LineStyle', ':'); + hx2 = graph2d.constantline((i-1)*1.05+0.125*(i-1)*(i-2)/2,'ydata', get(AX2(1),'ylim'), 'parent', AX2(1), 'color', 'black', 'LineStyle', ':'); + changedependvar(hx1,'x'); + changedependvar(hx2,'x'); +end + end \ No newline at end of file diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp index b730519..f56b06e 100644 --- a/Thalamic_Column.cpp +++ b/Thalamic_Column.cpp @@ -29,16 +29,19 @@ /* Initialization of RNG */ /****************************************************************************************************/ void Thalamic_Column::set_RNG(void) { - /* Number of independent streams */ - int N = 2; + extern const double dt; + /* Number of independent random variables */ + int N = 1; /* Create RNG for each stream */ 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; +double Thalamic_Column::noise_xRK(int N, int M) const{ + return gamma_e * gamma_e * (Rand_vars[2*M] + Rand_vars[2*M+1]/std::sqrt(3))*B[N]; +} + +double Thalamic_Column::noise_aRK(int M) const{ + return gamma_e * gamma_e * (Rand_vars[2*M] - Rand_vars[2*M+1]*std::sqrt(3))/4; } /****************************************************************************************************/ /* end */ @@ -262,28 +242,23 @@ double Thalamic_Column::noise_xRK(int N) const{ /****************************************************************************************************/ void Thalamic_Column::set_RK (int N) { extern const double dt; - _SWITCH((Ca) - (Phi_tt)(Phi_tr)(Phi_rt)(Phi_rr)(phi) - (x_tt) (x_tr) (x_rt) (x_rr) (y) - (h_T_t) (h_T_r) (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))); - Ca [N] = dt*(alpha_Ca * I_T_t(N) - (var_Ca - Ca_0)/tau_Ca); - 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_h [N] = dt*((m_inf_h(N) * (1 - var_m_h2) - var_m_h)/tau_m_h(N) - k3 * P_H(N) * var_m_h + k4 * var_m_h2); - m_h2 [N] = dt*(k3 * P_H(N) * var_m_h - k4 * var_m_h2); - P_h [N] = dt*(k1 * pow(var_Ca, n_P) - 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); - Phi_rr [N] = dt*(var_x_rr); - phi [N] = dt*(var_y); - x_tt [N] = dt*(pow(gamma_e, 2) * (noise_xRK(N) + N_et * Cortex->get_phi(N) - var_Phi_tt) - 2 * gamma_e * var_x_tt); - x_tr [N] = dt*(pow(gamma_e, 2) * (N_tr * get_Qt(N) + N_er * Cortex->get_phi(N) - var_Phi_tr) - 2 * gamma_e * var_x_tr); - x_rt [N] = dt*(pow(gamma_i, 2) * (N_rt * get_Qr(N) - var_Phi_rt) - 2 * gamma_i * var_x_rt); - x_rr [N] = dt*(pow(gamma_i, 2) * (N_rr * get_Qr(N) - var_Phi_rr) - 2 * gamma_i * var_x_rr); - y [N] = dt*(pow(nu, 2) * ( get_Qt(N) - var_phi) - 2 * nu * var_y); + Vt [N+1] = Vt [0] + A[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+1] = Vr [0] + A[N]*dt*(-(I_L_r(N) + I_er(N) + I_ir(N))/tau_r - (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); + y_tt [N+1] = y_tt [0] + A[N]*dt*(x_tt[N]); + y_tr [N+1] = y_tr [0] + A[N]*dt*(x_tr[N]); + y_rt [N+1] = y_rt [0] + A[N]*dt*(x_rt[N]); + y_rr [N+1] = y_rr [0] + A[N]*dt*(x_rr[N]); + y [N+1] = y [0] + A[N]*dt*(x [N]); + x_tt [N+1] = x_tt [0] + A[N]*dt*(pow(gamma_e, 2) * ( + N_et * Cortex->y[N] - y_tt[N]) - 2 * gamma_e * x_tt[N]) + noise_xRK(N,0); + x_tr [N+1] = x_tr [0] + A[N]*dt*(pow(gamma_e, 2) * (N_tr * get_Qt(N) + N_er * Cortex->y[N] - y_tr[N]) - 2 * gamma_e * x_tr[N]); + x_rt [N+1] = x_rt [0] + A[N]*dt*(pow(gamma_i, 2) * (N_rt * get_Qr(N) - y_rt[N]) - 2 * gamma_i * x_rt[N]); + x_rr [N+1] = x_rr [0] + A[N]*dt*(pow(gamma_i, 2) * (N_rr * get_Qr(N) - y_rr[N]) - 2 * gamma_i * x_rr[N]); + x [N+1] = x [0] + A[N]*dt*(pow(nu, 2) * ( get_Qt(N) - y [N]) - 2 * nu * x [N]); + 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]); } /****************************************************************************************************/ /* end */ @@ -294,26 +269,23 @@ void Thalamic_Column::set_RK (int N) { /* 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; - 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; - Phi_rr [0] += (Phi_rr [1] + Phi_rr [2] * 2 + Phi_rr [3] * 2 + Phi_rr [4])/6; - phi [0] += (phi [1] + phi [2] * 2 + phi [3] * 2 + phi [4])/6; - x_tt [0] += (x_tt [1] + x_tt [2] * 2 + x_tt [3] * 2 + x_tt [4])/6 + pow(gamma_e, 2) * h * Rand_vars[0]; - x_tr [0] += (x_tr [1] + x_tr [2] * 2 + x_tr [3] * 2 + x_tr [4])/6; - x_rt [0] += (x_rt [1] + x_rt [2] * 2 + x_rt [3] * 2 + x_rt [4])/6; - x_rr [0] += (x_rr [1] + x_rr [2] * 2 + x_rr [3] * 2 + x_rr [4])/6; - y [0] += (y [1] + y [2] * 2 + y [3] * 2 + y [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_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; - + 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; + y_tt [0] =(-3*y_tt [0] + 2*y_tt [1] + 4*y_tt [2] + 2* y_tt [3] + y_tt [4])/6; + y_tr [0] =(-3*y_tr [0] + 2*y_tr [1] + 4*y_tr [2] + 2* y_tr [3] + y_tr [4])/6; + y_rt [0] =(-3*y_rt [0] + 2*y_rt [1] + 4*y_rt [2] + 2* y_rt [3] + y_rt [4])/6; + y_rr [0] =(-3*y_rr [0] + 2*y_rr [1] + 4*y_rr [2] + 2* y_rr [3] + y_rr [4])/6; + y [0] =(-3*y [0] + 2*y [1] + 4*y [2] + 2* y [3] + y [4])/6; + x_tt [0] =(-3*x_tt [0] + 2*x_tt [1] + 4*x_tt [2] + 2* x_tt [3] + x_tt [4])/6 + noise_aRK(0); + x_tr [0] =(-3*x_tr [0] + 2*x_tr [1] + 4*x_tr [2] + 2* x_tr [3] + x_tr [4])/6; + x_rt [0] =(-3*x_rt [0] + 2*x_rt [1] + 4*x_rt [2] + 2* x_rt [3] + x_rt [4])/6; + x_rr [0] =(-3*x_rr [0] + 2*x_rr [1] + 4*x_rr [2] + 2* x_rr [3] + x_rr [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 -#include "macros.h" #include "Cortical_Column.h" using std::vector; @@ -47,6 +46,17 @@ typedef boost::variate_generator GEN; /* Variate generator */ /****************************************************************************************************/ +/* 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 { @@ -57,86 +67,87 @@ class Thalamic_Column { /* Constructor for simulation */ Thalamic_Column(double* Param, double* Con) - : g_LK_t (Param[1]), g_LK_r (Param[1]), g_h (Param[0]), + : g_LK (Param[1]), g_h (Param[0]), N_et (Con[0]), N_er (Con[1]) {set_RNG();} - /* Get the pointer to the cortical module */ - void get_Cortex (Cortical_Column& C) {Cortex = &C;} - - /* Get axonal flux */ - double get_phi (int N) const {_SWITCH((phi)); return var_phi;} + /* Get the pointer to the cortical module */ + void get_Cortex (Cortical_Column& C) {Cortex = &C;} - /* ODE functions */ - void set_RK (int); - void add_RK (void); + /* ODE functions */ + void set_RK (int); + void add_RK (void); /* Data storage access */ - friend void get_data (int, Cortical_Column&, Thalamic_Column&, _REPEAT(double*, 4)); + friend void get_data (int, Cortical_Column&, Thalamic_Column&, vector); + friend class Cortical_Column; - /* Set strength of input */ - void set_input (double I) {input = I;} + /* Set strength of input */ + void set_input (double I) {input = I;} private: - - /* 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_it (int) const; - double I_er (int) const; - double I_ir (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; - - /* 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) const; - - /* Population variables */ + /* 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_it (int) const; + double I_er (int) const; + double I_ir (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; + + /* 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 */ - 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 */ - phi = _INIT(0.0), /* axonal flux */ - 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 */ - y = _INIT(0.0), /* derivative of phi */ + y_tt = _INIT(0.0), /* PostSP from TC population to TC population */ + y_tr = _INIT(0.0), /* PostSP from TC population to RE population */ + y_rt = _INIT(0.0), /* PostSP from RE population to TC population */ + y_rr = _INIT(0.0), /* PostSP from RE population to RE population */ + y = _INIT(0.0), /* axonal flux */ + x_tt = _INIT(0.0), /* derivative of y_tt */ + x_tr = _INIT(0.0), /* derivative of y_tr */ + x_rt = _INIT(0.0), /* derivative of y_rt */ + x_rr = _INIT(0.0), /* derivative of y_rr */ + 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 */ - P_h = _INIT(0.0); /* messenger protein bound with calcium */ + m_h2 = _INIT(0.0); /* activation of h channel bound with protein */ /* Random number generators */ vector MTRands; @@ -173,12 +184,10 @@ class Thalamic_Column { /* Conductivities in mS/cm^-2 */ /* Leak current */ - const double g_L_t = 1; - const double g_L_r = 1; + const double g_L = 1; /* Potassium leak current */ - const double g_LK_t = 0.02; - const double g_LK_r = 0.02; + const double g_LK = 0.02; /* T current */ const double g_T_t = 3; @@ -220,19 +229,22 @@ class Thalamic_Column { /* Noise parameters in ms^-1 */ const double mphi = 0E-3; - const double dphi = 10E-3; + const double dphi = 20E-1; double input = 0.0; - /* Connectivities (dimensionless) */ const double N_tr = 3; - const double N_rt = 5; - const double N_rr = 19; - const double N_et = 2.5; - const double N_er = 2.5; + const double N_rt = 5; + const double N_rr = 19; + const double N_et = 2.6; + const double N_er = 2.6; /* Pointer to cortical column */ Cortical_Column* Cortex; + + /* SRK integration parameters */ + const vector A = {0.5, 0.5, 1.0, 1.0}; + const vector B = {0.75, 0.75, 0.0, 0.0}; }; /****************************************************************************************************/ /* end */ diff --git a/macros.h b/macros.h deleted file mode 100644 index 072533d..0000000 --- a/macros.h +++ /dev/null @@ -1,83 +0,0 @@ -/* -* Copyright (c) 2014 Michael Schellenberger Costa -* -* Permission is hereby granted, free of charge, to any person obtaining a copy -* of this software and associated documentation files (the "Software"), to deal -* in the Software without restriction, including without limitation the rights -* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -* copies of the Software, and to permit persons to whom the Software is -* furnished to do so, subject to the following conditions: -* -* The above copyright notice and this permission notice shall be included in -* all copies or substantial portions of the Software. -* -* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN -* THE SOFTWARE. -*/ - -/****************************************************************************************************/ -/* Definition of all macros used */ -/****************************************************************************************************/ -#ifndef MACROS_H -#define MACROS_H -#include -#include -#include -#include -#include -#include - - -/****************************************************************************************************/ -/* Macro for vector initialization */ -/****************************************************************************************************/ -#define _INIT(x) {x, 0.0, 0.0, 0.0, 0.0} -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* 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); -#define _SET_RK3(r, data, elem) BOOST_PP_CAT(data, elem) = (elem[0] + elem[2] * 0.5); -#define _SET_RK4(r, data, elem) BOOST_PP_CAT(data, elem) = (elem[0] + elem[3]); - -#define _SWITCH(...) \ - double BOOST_PP_SEQ_FOR_EACH_I(_GET_VARNAMES, var_, __VA_ARGS__); \ - switch(N) { \ - default: \ - BOOST_PP_SEQ_FOR_EACH(_SET_RK1, var_, __VA_ARGS__) \ - break; \ - case 2: \ - BOOST_PP_SEQ_FOR_EACH(_SET_RK2, var_, __VA_ARGS__) \ - break; \ - case 3: \ - BOOST_PP_SEQ_FOR_EACH(_SET_RK3, var_, __VA_ARGS__) \ - break; \ - case 4: \ - BOOST_PP_SEQ_FOR_EACH(_SET_RK4, var_, __VA_ARGS__) \ - break; \ - } -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* 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) -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ -#endif // MACROS_H diff --git a/saves.h b/saves.h index 60b2494..e9f0f62 100644 --- a/saves.h +++ b/saves.h @@ -30,11 +30,11 @@ /****************************************************************************************************/ /* Save data */ /****************************************************************************************************/ -inline void get_data(int counter, Cortical_Column& Cortex, Thalamic_Column& Thalamus, double* Ve, double* Vt, double* Ca, double* ah) { - Ve [counter] = Cortex.Ve [0]; - Vt [counter] = Thalamus.Vt [0]; - Ca [counter] = Cortex.Phi_ee [0]; - ah [counter] = Thalamus.act_h (); +inline void get_data(int counter, Cortical_Column& Cortex, Thalamic_Column& Thalamus, vector Data) { + Data[0][counter] = Cortex.Ve [0]; + Data[1][counter] = Thalamus.Vt [0]; + Data[2][counter] = Thalamus.Ca [0]; + Data[3][counter] = Thalamus.act_h (); } /****************************************************************************************************/ /* end */ @@ -44,13 +44,13 @@ inline void get_data(int counter, Cortical_Column& Cortex, Thalamic_Column& Thal /****************************************************************************************************/ /* Create MATLAB data container */ /****************************************************************************************************/ -mxArray* SetMexArray(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* 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; } /****************************************************************************************************/ /* end */