/NM_TC

Permalink
Browse files

Imported the project into QTcreator, 

More parameter fixes and first tries with bursting stimulation.
  • Loading branch information...
1 parent 8ae1ecc commit 8385cf5dfd392108fa420004885606bcec91004d @miscco committed Dec 4, 2014
Unified Split
Showing with 102 additions and 33 deletions.
  1. +3 −3 Cortical_Column.cpp
  2. +1 −1 Main.cpp
  3. +25 −14 Plots.m
  4. +44 −6 Stimulation.h
  5. +20 −0 TC_model.pro
  6. +5 −5 Thalamic_Column.h
  7. +4 −4 saves.h
View
6 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;
}
View
2 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 */
/****************************************************************************************************/
View
39 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' );
%
View
50 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 */
View
20 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
View
10 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;
View
8 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));

0 comments on commit 8385cf5

Please sign in to comment.