diff --git a/Cortical_Column.h b/Cortical_Column.h index 8d4972d..562240c 100644 --- a/Cortical_Column.h +++ b/Cortical_Column.h @@ -65,40 +65,40 @@ class Cortical_Column { /* Return axonal flux */ double get_phi (int N) const {_SWITCH((phi)); return var_phi;} - /* Initialize the RNGs */ - void set_RNG (void); + /* Data storage access */ + friend void get_data (int, Cortical_Column&, Thalamic_Column&, _REPEAT(double*, 4)); - /* Firing rates */ - double get_Qe (int) const; - double get_Qi (int) const; + /* ODE functions */ + void set_RK (int); + void add_RK (void); - /* 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; + /* Stimulation protocol access */ + friend class Stim; - /* Potassium pump */ - double Na_pump (int) const; +private: + /* Initialize the RNGs */ + void set_RNG (void); - /* Noise function */ - double noise_xRK (int, int) const; + /* Firing rates */ + double get_Qe (int) const; + double get_Qi (int) const; - /* ODE functions */ - void set_RK (int); - void add_RK (void); + /* 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; - /* Data storage access */ - friend void get_data (int, Cortical_Column&, Thalamic_Column&, _REPEAT(double*, 4)); + /* Potassium pump */ + double Na_pump (int) const; - /* Stimulation protocol access */ - friend class Stim; + /* Noise function */ + double noise_xRK (int, int) const; -private: - /* Population variables */ + /* 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 */ @@ -141,9 +141,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 = 2; /* 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 */ @@ -174,14 +174,14 @@ class Cortical_Column { /* Noise parameters in ms^-1 */ const double mphi = 0E-3; - const double dphi = 60E-3;; + const double dphi = 60E-3; double input = 0.0; /* Connectivities (dimensionless) */ - const double N_ee = 120; + const double N_ee = 125; const double N_ei = 72; - const double N_ie = 90; - const double N_ii = 90; + const double N_ie = 100; + const double N_ii = 100; const double N_te = 0; const double N_ti = 0; diff --git a/Plots.m b/Plots.m index cb10ba0..24bf3f8 100644 --- a/Plots.m +++ b/Plots.m @@ -1,31 +1,31 @@ % mex command is given by: -% mex CXXFLAGS="\$CXXFLAGS -std=c++11 -O3" TC.cpp Cortical_Column.cpp Thalamic_Column.cpp +% mex CXXFLAGS="\$CXXFLAGS -std=c++11 -O3 -lgopm" TC.cpp Cortical_Column.cpp Thalamic_Column.cpp function Plots(type) if nargin == 0 - type = 2; + type = 1; end if type == 1 Param_Cortex = [4.7; % sigma_e - 1.33; % g_KNa + 1.43; % g_KNa 120E-3]; % dphi Param_Thalamus = [0.051; % g_h 0.024]; % g_LK else - Param_Cortex = [6.; % sigma_e - 2.1; % g_KNa + Param_Cortex = [6; % sigma_e + 2.15; % g_KNa 120E-3]; % dphi Param_Thalamus = [0.051; % g_h - 0.02]; % g_LK + 0.02]; % g_LK end -Connectivity = [ 2.5; % N_et - 2.5; % N_er +Connectivity = [ 2.5; % N_et + 2.5; % N_er 6; % N_te - 15]; % N_ti + 15]; % N_ti % stimulation parameters % first number is the mode of stimulation @@ -34,11 +34,11 @@ function Plots(type) % 2 == phase dependend var_stim = [ 0; % mode of stimulation - 300; % strength of the stimulus in Hz (spikes per second) + 60; % strength of the stimulus in Hz (spikes per second) 120; % duration of the stimulus in ms 5; % time between stimulation events in s (ISI) 0; % range of ISI in s [ISI-range,ISI+range] - 1; % Number of stimuli per event + 3; % Number of stimuli per event 1050; % time between stimuli within a event in ms 450]; % time until stimuli after minimum in ms @@ -49,7 +49,7 @@ function Plots(type) L = length(Vt); timeaxis = linspace(0,T,L); -figure(1) +%figure(1) subplot(411), plot(timeaxis,Ve) title('Pyramidal membrane voltage'), xlabel('Time in s'), @@ -58,8 +58,8 @@ function Plots(type) ylim([-80, -40]); % 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'); + %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(412), plot(timeaxis,Vt) diff --git a/Stimulation.h b/Stimulation.h index 9810470..e372121 100644 --- a/Stimulation.h +++ b/Stimulation.h @@ -105,8 +105,11 @@ class Stim { /* If a stimulation event has occurred and there is a minimal time (pause) until the next one */ bool stimulation_paused = false; + /* If burst mode is enabled */ + bool burst_enabled = false; + /* In case of bursted stimulation start the bursts */ - bool burst_started = false; + bool burst_started = true; /* Length of a burst stimulus */ int burst_length = 10; @@ -183,8 +186,8 @@ void Stim::setup (double* var_stim) { time_between_stimuli = (int) var_stim[6] * res / 1000; /* Scale the length of burst_length and burst_ISI from ms to dt*/ - burst_length = (int) 6 * res / 1000; - burst_ISI = (int) 44 * res / 1000; + burst_length = (int) 2 * res / 1000; + burst_ISI = (int) 28 * res / 1000; /* If ISI is fixed do not create RNG */ if (mode == 1) { @@ -310,19 +313,21 @@ void Stim::check_stim (int time) { count_bursts++; /* Switch stimulation on and off wrt burst parameters */ - if(burst_started) { - if(count_bursts%burst_length==0) { - count_bursts = 0; - burst_started = false; - Thalamus->set_input(0.0); - } - } else { - if(count_bursts%burst_ISI==0) { - count_bursts = 0; - burst_started = true; - Thalamus->set_input(strength); - } - } + if(burst_enabled) { + if(burst_started) { + if(count_bursts%burst_length==0) { + count_bursts = 0; + burst_started = false; + Thalamus->set_input(0.0); + } + } else { + if(count_bursts%burst_ISI==0) { + count_bursts = 0; + burst_started = true; + Thalamus->set_input(strength); + } + } + } } /* Wait if there is a pause between stimulation events */ diff --git a/TC.cpp b/TC.cpp index b2a4737..99cbda4 100644 --- a/TC.cpp +++ b/TC.cpp @@ -103,8 +103,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { plhs[1] = Vt; plhs[2] = Ca; plhs[3] = ah; - plhs[4] = Stimulation.get_marker(); -return; + plhs[4] = Stimulation.get_marker(); + + return; } /****************************************************************************************************/ /* end */ diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp index bbde587..b730519 100644 --- a/Thalamic_Column.cpp +++ b/Thalamic_Column.cpp @@ -51,7 +51,7 @@ void Thalamic_Column::set_RNG(void) { /****************************************************************************************************/ /* Thalamic relay (TC) firing rate */ double Thalamic_Column::get_Qt (int N) const{ - _SWITCH((Vt)) + _SWITCH((Vt)) double q = Qt_max/ (1 + exp(-C1 * (var_Vt - theta_t) / sigma_t)); return q; } diff --git a/Thalamic_Column.h b/Thalamic_Column.h index 53f8988..9706e77 100644 --- a/Thalamic_Column.h +++ b/Thalamic_Column.h @@ -61,63 +61,64 @@ class Thalamic_Column { 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;} - - /* Set strength of input */ - void set_input (double I) {input = I;} - - /* Get axonal flux */ - double get_phi (int N) const {_SWITCH((phi)); return var_phi;} - - /* 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; - - /* ODE functions */ - void set_RK (int); - void add_RK (void); + /* Get the pointer to the cortical module */ + void get_Cortex (Cortical_Column& C) {Cortex = &C;} + + /* Get axonal flux */ + double get_phi (int N) const {_SWITCH((phi)); return var_phi;} + + /* 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)); + /* Set strength of input */ + void set_input (double I) {input = I;} + private: - /* 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; + + /* 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 */ 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 */ diff --git a/macros.h b/macros.h index 1baa7a1..072533d 100644 --- a/macros.h +++ b/macros.h @@ -23,7 +23,8 @@ /****************************************************************************************************/ /* Definition of all macros used */ /****************************************************************************************************/ -#pragma once +#ifndef MACROS_H +#define MACROS_H #include #include #include @@ -79,3 +80,4 @@ /****************************************************************************************************/ /* end */ /****************************************************************************************************/ +#endif // MACROS_H diff --git a/saves.h b/saves.h index 3450a8d..12a9284 100644 --- a/saves.h +++ b/saves.h @@ -32,8 +32,8 @@ /****************************************************************************************************/ 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] = Thalamus.phi [0]; + Vt [counter] = Cortex.N_te*Thalamus.phi [0]; + Ca [counter] = Cortex.Phi_ee [0]; ah [counter] = Thalamus.act_h (); } /****************************************************************************************************/ @@ -48,7 +48,8 @@ mxArray* SetMexArray(int N, int M) { mxArray* Array = mxCreateDoubleMatrix(0, 0, mxREAL); mxSetM(Array, N); mxSetN(Array, M); - mxSetData(Array, mxMalloc(sizeof(double)*M*N)); + #pragma omp critical + {mxSetData(Array, mxMalloc(sizeof(double)*M*N));} return Array; } /****************************************************************************************************/