diff --git a/Cortical_Column.cpp b/Cortical_Column.cpp index 630b214..947eb6d 100644 --- a/Cortical_Column.cpp +++ b/Cortical_Column.cpp @@ -52,14 +52,14 @@ void Cortical_Column::set_RNG(void) { /* Pyramidal firing rate */ double Cortical_Column::get_Qe (int N) const{ _SWITCH((Ve)) - double q = Qe_max / (1 + exp(-C1 * (var_Ve - theta_e) / sigma_e)); + double q = Qe_max / (1 + exp(-C1 * (var_Ve - theta_e) / sigma_e)); return q; } /* Inhibitory firing rate */ double Cortical_Column::get_Qi (int N) const{ _SWITCH((Vi)) - double q = Qi_max / (1 + exp(-C1 * (var_Vi - theta_i) / sigma_i)); + double q = Qi_max / (1 + exp(-C1 * (var_Vi - theta_i) / sigma_i)); return q; } /****************************************************************************************************/ @@ -121,7 +121,7 @@ double Cortical_Column::I_L_i (int N) const{ /* Sodium dependent potassium current */ double Cortical_Column::I_KNa (int N) const{ _SWITCH((Ve)(Na)) - double w_KNa = 0.37/(1+pow(38.7/var_Na, 3.5)); + double w_KNa = 0.37/(1+pow(38.7/var_Na, 3.5)); double I_KNa = g_KNa * w_KNa * (var_Ve - E_K); return I_KNa; } diff --git a/Main.cpp b/Main.cpp index 7d26e1a..c288b3e 100644 --- a/Main.cpp +++ b/Main.cpp @@ -37,7 +37,7 @@ /****************************************************************************************************/ extern const int T = 30; /* Simulation length s */ extern const int res = 1E4; /* Number of iteration steps per s */ -extern const int red = 1E2; /* Fraction of iterations that is saved */ +extern const int red = 1E2; /* Fraction of iterations that is saved */ extern const double dt = 1E3/res; /* Duration of a iteration step in ms */ extern const double h = sqrt(dt); /* Square root of dt for SRK iteration */ /****************************************************************************************************/ diff --git a/Plots.m b/Plots.m index f1983c1..cb10ba0 100644 --- a/Plots.m +++ b/Plots.m @@ -6,28 +6,26 @@ function Plots(type) type = 2; end -if type == 1 - +if type == 1 Param_Cortex = [4.7; % sigma_e - 1.35; % g_KNa + 1.33; % g_KNa 120E-3]; % dphi - Param_Thalamus = [0.053; % g_h - 0.024]; % g_LK - + Param_Thalamus = [0.051; % g_h + 0.024]; % g_LK else Param_Cortex = [6.; % sigma_e 2.1; % g_KNa 120E-3]; % dphi - Param_Thalamus = [0.055; % g_h + Param_Thalamus = [0.051; % g_h 0.02]; % g_LK end Connectivity = [ 2.5; % N_et 2.5; % N_er - 5; % N_te - 10]; % N_ti + 6; % N_te + 15]; % N_ti % stimulation parameters % first number is the mode of stimulation @@ -36,11 +34,11 @@ function Plots(type) % 2 == phase dependend var_stim = [ 0; % mode of stimulation - 80; % strength of the stimulus in Hz (spikes per second) + 300; % 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] - 3; % Number of stimuli per event + 1; % Number of stimuli per event 1050; % time between stimuli within a event in ms 450]; % time until stimuli after minimum in ms @@ -52,7 +50,7 @@ function Plots(type) timeaxis = linspace(0,T,L); figure(1) -subplot(311), plot(timeaxis,Ve) +subplot(411), plot(timeaxis,Ve) title('Pyramidal membrane voltage'), xlabel('Time in s'), ylabel('V_{e} in mV') @@ -64,7 +62,7 @@ function Plots(type) changedependvar(hx,'x'); end -subplot(312), plot(timeaxis,Vt) +subplot(412), plot(timeaxis,Vt) title('Thalamic relay membrane voltage'), xlabel('Time in s'), ylabel('V_{t} in mV') @@ -76,7 +74,7 @@ function Plots(type) changedependvar(hx,'x'); end -subplot(313), plot(timeaxis,ah) +subplot(413), plot(timeaxis,ah) title('Thalamic relay I_h activation'), xlabel('Time in s'), ylabel('m_h in mV') @@ -87,6 +85,19 @@ function Plots(type) 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'), +xlabel('Time in s'), +ylabel('[Ca] in \mu mol') +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/Stimulation.h b/Stimulation.h index 5b23547..9810470 100644 --- a/Stimulation.h +++ b/Stimulation.h @@ -44,6 +44,7 @@ class Stim { /* Empty constructor for compiling */ Stim(void); + /* Constructor with references and stimulation variables */ Stim(Cortical_Column& C, Thalamic_Column& T, double* var) { Cortex = &C; Thalamus = &T; @@ -63,7 +64,7 @@ class Stim { /* 0 == none */ /* 1 == semi-periodic */ /* 2 == phase dependent */ - int mode = 0; + int mode = 0; /* Default values already in dt: E1==ms, E4==s */ /* Stimulation strength */ @@ -104,12 +105,24 @@ class Stim { /* If a stimulation event has occurred and there is a minimal time (pause) until the next one */ bool stimulation_paused = false; + /* In case of bursted stimulation start the bursts */ + bool burst_started = false; + + /* Length of a burst stimulus */ + int burst_length = 10; + + /* Length of silence between burst stimuli */ + int burst_ISI = 10; + /* Onset in time steps where no data is recorded, so stimulation has to be delayed */ int onset_correction = 10E4; /* Counter for number of stimuli that occurred within a stimulation event */ int count_stimuli = 1; + /* Counter for number of bursts per stimuli */ + int count_bursts = 0; + /* Counter for stimulation duration since started*/ int count_duration = 0; @@ -120,7 +133,7 @@ class Stim { int count_pause = 0; /* Old voltage value for minimum detection */ - double Ve_old = 0; + double Ve_old = 0.0; /* Pointer to columns */ Cortical_Column* Cortex; @@ -169,6 +182,10 @@ void Stim::setup (double* var_stim) { /* Scale time_between_stimuli from ms to dt */ time_between_stimuli = (int) var_stim[6] * res / 1000; + /* Scale the length of burst_length and burst_ISI from ms to dt*/ + burst_length = (int) 6 * res / 1000; + burst_ISI = (int) 44 * res / 1000; + /* If ISI is fixed do not create RNG */ if (mode == 1) { /* Set first time_to_stimuli to 1 sec after onset */ @@ -187,7 +204,7 @@ void Stim::setup (double* var_stim) { /* Scale time_to_stimuli from ms to dt */ time_to_stimuli = (int) var_stim[7] * res / 1000; } -}; +} void Stim::check_stim (int time) { /* Check if stimulation should start */ @@ -272,19 +289,40 @@ void Stim::check_stim (int time) { count_stimuli = 1; } } + /* Update counter */ count_to_start++; } break; } - /* Wait to switch the stimulation off */ + /* Actual stimulation protocols */ if(stimulation_started) { + /* Wait to switch the stimulation off */ if(count_duration==duration) { stimulation_started = false; + burst_started = true; count_duration = 0; + count_bursts = 0; Thalamus->set_input(0.0); } + count_duration++; + count_bursts++; + + /* Switch stimulation on and off wrt burst parameters */ + if(burst_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 */ @@ -295,7 +333,7 @@ void Stim::check_stim (int time) { } count_pause++; } -}; +} mxArray* Stim::get_marker(void) { extern const int red; @@ -309,7 +347,7 @@ mxArray* Stim::get_marker(void) { Pr_Marker[i] = marker_stimulation [i]/red; } return Marker; -}; +} /****************************************************************************************************/ /* end */ diff --git a/TC_model.pro b/TC_model.pro new file mode 100644 index 0000000..6ed6fe0 --- /dev/null +++ b/TC_model.pro @@ -0,0 +1,20 @@ +TEMPLATE = app +CONFIG += console +CONFIG -= app_bundle +CONFIG -= qt + +SOURCES += Main.cpp \ + Cortical_Column.cpp \ + Thalamic_Column.cpp \ + TC.cpp + +HEADERS += macros.h \ + ODE.h \ + saves.h \ + Cortical_Column.h \ + Thalamic_Column.h \ + Stimulation.h + +QMAKE_CXXFLAGS += -std=c++11 -O3 + +SOURCES -= TC.cpp diff --git a/Thalamic_Column.h b/Thalamic_Column.h index 373cdf3..53f8988 100644 --- a/Thalamic_Column.h +++ b/Thalamic_Column.h @@ -153,8 +153,8 @@ class Thalamic_Column { const double Qr_max = 400.E-3; /* Sigmoid threshold in mV */ - const double theta_t = -58.6; - const double theta_r = -58.6; + const double theta_t = -58.5; + const double theta_r = -58.5; /* Sigmoid gain in mV */ const double sigma_t = 6; @@ -219,14 +219,14 @@ class Thalamic_Column { /* Noise parameters in ms^-1 */ const double mphi = 0E-3; - const double dphi = 10E-3;; + const double dphi = 10E-3; double input = 0.0; /* Connectivities (dimensionless) */ const double N_tr = 3; - const double N_rt = 6; - const double N_rr = 19; + const double N_rt = 5; + const double N_rr = 19; const double N_et = 5; const double N_er = 5; diff --git a/saves.h b/saves.h index 29f447e..3450a8d 100644 --- a/saves.h +++ b/saves.h @@ -32,9 +32,9 @@ /****************************************************************************************************/ 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.Ca [0]; - ah [counter] = Thalamus.act_h (); + Vt [counter] = Thalamus.Vt [0]; + Ca [counter] = Thalamus.phi [0]; + ah [counter] = Thalamus.act_h (); } /****************************************************************************************************/ /* end */ @@ -45,7 +45,7 @@ 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); + mxArray* Array = mxCreateDoubleMatrix(0, 0, mxREAL); mxSetM(Array, N); mxSetN(Array, M); mxSetData(Array, mxMalloc(sizeof(double)*M*N));