diff --git a/.cproject b/.cproject
index 4957c40..bcd7e28 100644
--- a/.cproject
+++ b/.cproject
@@ -23,7 +23,7 @@
-
+
@@ -33,7 +33,7 @@
-
+
@@ -101,7 +101,7 @@
-
+
diff --git a/Cortical_Column.h b/Cortical_Column.h
index 1196be9..af7af58 100644
--- a/Cortical_Column.h
+++ b/Cortical_Column.h
@@ -54,9 +54,9 @@ class Cortical_Column {
Cortical_Column(void)
{set_RNG();}
- Cortical_Column(double* Par, double* Con)
- :sigma_e (Par[0]), g_KNa (Par[1]), dphi (Par[2]),
- N_te (Con[2]), N_ti (Con[3])
+ Cortical_Column(double* Param, double* Con)
+ :sigma_e (Param[0]), g_KNa (Param[1]), dphi (Param[2]),
+ N_te (Con[2]), N_ti (Con[3])
{set_RNG();}
/* Connect to the thalamic module */
@@ -94,7 +94,7 @@ class Cortical_Column {
/* Data storage access */
friend void get_data (int, Cortical_Column&, Thalamic_Column&, _REPEAT(double*, 2));
- /* Stimulation protocoll access */
+ /* Stimulation protocol access */
friend class Stim;
private:
@@ -141,7 +141,7 @@ 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; /* Sodium time constant in ms */
+ const double tau_Na = 1.5; /* Sodium time constant in 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 */
@@ -161,7 +161,7 @@ class Cortical_Column {
const double g_KNa = 1.33;
/* Reversal potentials in mV */
- /* synaptic */
+ /* Synaptic */
const double E_AMPA = 0;
const double E_GABA = -70;
diff --git a/Main.cpp b/Main.cpp
index 5a72229..7d26e1a 100644
--- a/Main.cpp
+++ b/Main.cpp
@@ -36,10 +36,10 @@
/* Fixed simulation settings */
/****************************************************************************************************/
extern const int T = 30; /* Simulation length s */
-extern const int res = 1E4; /* number of iteration steps per s */
-extern const int red = res/100; /* number of iterations that is saved */
-extern const double dt = 1E3/res; /* duration of a timestep in ms */
-extern const double h = sqrt(dt); /* squareroot of dt for SRK iteration */
+extern const int res = 1E4; /* Number of iteration steps per s */
+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 */
/****************************************************************************************************/
/* end */
/****************************************************************************************************/
diff --git a/Plots.m b/Plots.m
index d1b52f1..7879ab8 100644
--- a/Plots.m
+++ b/Plots.m
@@ -4,19 +4,29 @@
function Plots(T)
if nargin == 0
- Input_N3 = [ 6.5; % sigma_e
- 2.1; % g_KNa
- 120E-3]; % dphi
+
+ Param_Cortex_N2 = [4.7; % sigma_e
+ 1.43; % g_KNa
+ 120E-3]; % dphi
+
+ Param_Cortex_N3 = [6.3; % sigma_e
+ 2.; % g_KNa
+ 120E-3]; % dphi
+
+
+ Param_Thalamus_N2 = [0.025; % g_LK_t
+ 0.025; % g_LK_r
+ 0.08]; % g_h
- Input_N2 = [ 4.7; % sigma_e
- 1.5; % g_KNa
- 120E-3]; % dphi
+ Param_Thalamus_N3 = [0.021; % g_LK_t
+ 0.021; % g_LK_r
+ 0.08]; % g_h
- Connectivity= [2.4; % N_et
- 2.5; % N_er
- 5; % N_te
- 5]; % N_ti
+ Connectivity = [2.4; % N_et
+ 2.6; % N_er
+ 5; % N_te
+ 10]; % N_ti
% stimulation parameters
% first number is the mode of stimulation
@@ -27,14 +37,15 @@ function Plots(T)
var_stim = [ 0; % mode of stimulation
60; % strength of the stimulus in Hz (spikes per second)
- 100; % duration of the stimulus in ms
+ 200; % duration of the stimulus in ms
5; % time between stimuli in s
- 650]; % time until stimuli after negativ peak in ms
+ 300]; % time until stimuli after negativ peak in ms
- T = 30; % duration of the simulation
+ T = 600; % duration of the simulation
end
-[Ve, Vt, Marker_Stim] = TC(T, Input_N2, Connectivity, var_stim);
+%[Ve, Vt, Marker_Stim] = TC(T, Param_Cortex_N2, Param_Thalamus_N2, Connectivity, var_stim);
+[Ve, Vt, Marker_Stim] = TC(T, Param_Cortex_N3, Param_Thalamus_N3, Connectivity, var_stim);
L = max(size(Vt));
timeaxis = linspace(0,T,L);
@@ -42,22 +53,24 @@ function Plots(T)
figure(1)
subplot(211), plot(timeaxis,Ve)
title('Pyramidal membrane voltage'), xlabel('time in s'), ylabel('Ve in mV')
+ylim([-80, -40])
% vertical line for markers
-%hx = graph2d.constantline(Marker_Stim(3,:));
-%changedependvar(hx,'x');
+hx1 = graph2d.constantline(Marker_Stim(2,:), 'color', 'black');
+hx2 = graph2d.constantline(Marker_Stim(1,:), 'color', 'red');
+changedependvar(hx1,'x');
+changedependvar(hx2,'x');
subplot(212), plot(timeaxis,Vt)
title('Thalamic relay membrane voltage'), xlabel('time in s'), ylabel('Vt in mV')
% vertical line for markers
-%hx = graph2d.constantline(Marker_Stim(3,:));
-%changedependvar(hx,'x');
+hx = graph2d.constantline(Marker_Stim(2,:));
+changedependvar(hx,'x');
+[Pxx,f] = pwelch(Ve-mean(Ve),hamming(L/30), 4*L/T, 2048, L/T);
+n = find(f<=30, 1, 'last' );
-% [Pxx,f] = pwelch(Ve-mean(Ve),[], [], [], L/T);
-% n = find(f<=30, 1, 'last' );
-%
-% figure(2)
-% plot(f(1:n),log(Pxx(1:n)))
-% title('Powerspectrum with pwelch'), xlabel('frequency in Hz'), ylabel('Power (log)')
-% save('Timeseries', 'Ve', 'Vt');
+figure(2)
+plot(f(1:n),log(Pxx(1:n)))
+title('Powerspectrum with pwelch'), xlabel('frequency in Hz'), ylabel('Power (log)')
+%save('Timeseries', 'Ve', 'Vt');
end
\ No newline at end of file
diff --git a/Stimulation.h b/Stimulation.h
index 6e973b2..662d75e 100644
--- a/Stimulation.h
+++ b/Stimulation.h
@@ -24,10 +24,19 @@
/* Implementation of the stimulation protocol */
/****************************************************************************************************/
#pragma once
+#include
#include "Cortical_Column.h"
#include "Thalamic_Column.h"
/****************************************************************************************************/
+/* Typedefs for RNG */
+/****************************************************************************************************/
+typedef boost::random::uniform_int_distribution<> DIST_int; /* Uniform integer distribution */
+/****************************************************************************************************/
+/* end */
+/****************************************************************************************************/
+
+/****************************************************************************************************/
/* Stimulation object */
/****************************************************************************************************/
class Stim {
@@ -40,221 +49,277 @@ class Stim {
Thalamus = &T;
setup(var);}
- /* Setup with respect to stimulation mode */
- void setup (double* var_stim) {
- extern const int onset;
- extern const int res;
- extern const int dt;
+ /* Initialize stimulation class with respect to stimulation mode */
+ void setup (double* var_stim);
- /* Mode of stimulation */
- mode = (int) var_stim[0];
+ /* Check whether stimulation should be started/stopped */
+ void check_stim (int time);
- /* Scale the stimulation strength from s^-1 (Hz) to ms^-1 */
- strength = var_stim[1] / 1000;
+ /* Create MATLAB container for marker storage */
+ mxArray* get_marker(void);
- /* Scale duration from ms to dt */
- duration = (int) var_stim[2] * res / 1000;
+private:
+ /* Mode of stimulation */
+ /* 0 == none */
+ /* 1 == semi-periodic */
+ /* 2 == phase dependent */
+ int mode = 0;
- /* Scale the ISI from s to ms */
- ISI = (int) var_stim[3] * res;
+ /* Default values already in dt: E1==ms, E4==s */
+ /* Stimulation strength */
+ double strength = 0.0;
- /* Scale time to stimulus from ms to dt */
- time_to_stim= (int) var_stim[4] * res / 1000;
+ /* Duration of the stimulation */
+ int duration = 120E1;
- /* Set the onset correction for the marker */
- correction = onset * res;
+ /* Interval between different Stimulus events */
+ int ISI = 5E4;
- if(mode==1) {
- time_to_stim = (int) (onset+1) * res;
- }
+ /* Range of Inter Stimulus Interval */
+ int ISI_range = 1E4;
- correction = onset * res;
+ /* Number of stimuli in case of multiple stimuli per event */
+ int number_of_stimuli = 1;
- }
+ /* Time until next stimulus */
+ /* Function varies between different stimulation modes */
+ int time_to_stimuli = 350E1;
- void check_stim (int time) {
+ /* Time between stimuli in case of multiple stimuli per event */
+ int time_between_stimuli = 1050E1;
- /* Check if stimulation should start */
- switch (mode) {
+ /* Threshold for phase dependent stimulation */
+ double threshold = -72;
- /* No stimulation */
- default:
- break;
+ /* Internal variables */
+ /* Simulation on for TRUE and off for FALSE */
+ bool stimulation_started = false;
- /* Periodic stimulation */
- case 1:
- /* Check if stimulation time is reached */
- if(time == time_to_stim) {
- /* Switch stimulation on */
- stimulation_started = true;
- Thalamus->set_input(strength);
+ /* If threshold has been reached */
+ bool threshold_crossed = false;
- /* Update the timer */
- time_to_stim += ISI;
+ /* If minimum was found */
+ bool minimum_found = false;
- /* Add marker */
- marker_threshold.push_back(0);
- marker_minimum.push_back(0);
- marker_stimulation.push_back(time - correction);
- }
- break;
-
- /* Phase dependent up state stimulation */
- case 2:
- /* Search for threshold */
- if(!stimulation_started && !minimum_found && !threshold_crossed && time>correction) {
- if(Cortex->Ve[0]<=threshold) {
- threshold_crossed = true;
- marker_threshold.push_back(time - correction);
- }
- }
+ /* If a stimulation event has occurred and there is a minimal time (pause) until the next one */
+ bool stimulation_paused = false;
- /* Search for minimum */
- if(threshold_crossed) {
- if(Cortex->Ve[0]>Ve_old) {
- threshold_crossed = false;
- minimum_found = true;
- marker_minimum.push_back(time - correction);
- Ve_old = 0;
- } else {
- Ve_old = Cortex->Ve[0];
- }
- }
+ /* Onset in time steps where no data is recorded, so stimulation has to be delayed */
+ int onset_correction = 10E4;
- /* Wait until the stimulation should start */
- if(minimum_found) {
- count_to_start++;
+ /* Counter for number of stimuli that occurred within a stimulation event */
+ int count_stimuli = 1;
+ /* Counter for stimulation duration since started*/
+ int count_duration = 0;
- /* Start stimulation after time_to_stim has passed */
- if(count_to_start==time_to_stim) {
- minimum_found = false;
- stimulation_started = true;
- count_to_start = 0;
- marker_stimulation.push_back(time - correction);
- Thalamus->set_input(strength);
- }
- }
- break;
-
- /* Phase dependent down state stimulation */
- case 3:
- /* Search for threshold */
- if(!stimulation_started && !minimum_found && !threshold_crossed && time>correction) {
- if(Cortex->Ve[0]<=threshold) {
- threshold_crossed = true;
- marker_threshold.push_back(time - correction);
- }
- }
+ /* Counter after minimum was found */
+ int count_to_start = 0;
- /* Search for minimum */
- if(threshold_crossed) {
- if(Cortex->Ve[0]>Ve_old) {
- threshold_crossed = false;
- minimum_found = true;
- marker_minimum.push_back(time - correction);
- Ve_old = 0;
- } else {
- Ve_old = Cortex->Ve[0];
- }
- }
+ /* Counter for time between two stimulation events (with multiple tones) */
+ int count_pause = 0;
- /* Start the stimulation */
- if(minimum_found) {
- minimum_found = false;
- stimulation_started = true;
- marker_stimulation.push_back(time - correction);
- Thalamus->set_input(strength);
- }
- break;
- }
+ /* Old voltage value for minimum detection */
+ double Ve_old = 0;
- /* Wait to switch the stimulation off */
- if(stimulation_started) {
- count_duration++;
+ /* Pointer to columns */
+ Cortical_Column* Cortex;
+ Thalamic_Column* Thalamus;
- /* Switch stimulation off */
- if(count_duration==duration) {
- stimulation_started = false;
- count_duration = 0;
- Thalamus->set_input(0.0);
- }
- }
- }
+ /* Data containers */
+ vector marker_minimum;
+ vector marker_stimulation;
- /* Create MATLAB container for marker storage */
- mxArray* get_marker(void) {
- extern const int res;
- mxArray* Marker = mxCreateDoubleMatrix(0, 0, mxREAL);
- mxSetM(Marker, 3);
- mxSetN(Marker, marker_stimulation.size());
- mxSetData(Marker, mxMalloc(sizeof(double)*3*marker_stimulation.size()));
- double* Pr_Marker = mxGetPr(Marker);
- for(unsigned i=0; iset_input(strength);
+
+ /* Check if multiple stimuli should be applied */
+ if (count_stimuli < number_of_stimuli) {
+ /* Update the timer with respect to time between stimuli */
+ time_to_stimuli += time_between_stimuli;
+ count_stimuli++;
+ }
+ /* After last stimulus in event update the timer with respect to (random) ISI*/
+ else {
+ time_to_stimuli += (ISI_range==0)? ISI : Uniform_Distribution(Generator);
- /* Pointer to thalamic column */
- Thalamic_Column* Thalamus;
+ /* Reset the stimulus counter for next stimulation event */
+ if (number_of_stimuli != 1) {
+ count_stimuli = 1;
+ }
+ }
- /* Data containers */
- std::vector marker_threshold;
- std::vector marker_minimum;
- std::vector marker_stimulation;
+ /* Add marker */
+ marker_minimum.push_back(0);
+ marker_stimulation.push_back(time - onset_correction);
+ }
+ break;
+
+ /* Phase dependent stimulation */
+ case 2:
+ /* Search for threshold */
+ if(!stimulation_started && !minimum_found && !threshold_crossed && time>onset_correction && !stimulation_paused) {
+ if(Cortex->Ve[0]<=threshold) {
+ threshold_crossed = true;
+ }
+ }
+
+ /* Search for minimum */
+ if(threshold_crossed) {
+ if(Cortex->Ve[0]>Ve_old) {
+ threshold_crossed = false;
+ minimum_found = true;
+ /* add one marker for every stimuli in the event */
+ for(int i=0; iVe[0];
+ }
+ }
+
+ /* Wait until the stimulation should start */
+ if(minimum_found) {
+ /* Start stimulation after time_to_stimuli has passed */
+ if(count_to_start==time_to_stimuli) {
+ stimulation_started = true;
+ Thalamus->set_input(strength);
+
+ /* Check if multiple stimuli should be applied */
+ if (count_stimuli < number_of_stimuli) {
+ /* Update the timer with respect to time between stimuli */
+ time_to_stimuli += time_between_stimuli;
+ count_stimuli++;
+ } else {
+ /* After last stimulus in event pause the stimulation */
+ minimum_found = false;
+ stimulation_paused = true;
+ count_to_start = 0;
+
+ /* Reset the stimulus counter for next stimulation event */
+ if (number_of_stimuli != 1) {
+ count_stimuli = 1;
+ }
+ }
+
+ /* Add marker */
+ marker_stimulation.push_back(time - onset_correction);
+ }
+
+ count_to_start++;
+ }
+ break;
+ }
+
+ /* Wait to switch the stimulation off */
+ if(stimulation_started) {
+ if(count_duration==duration) {
+ stimulation_started = false;
+ count_duration = 0;
+ Thalamus->set_input(0.0);
+ }
+ count_duration++;
+ }
+
+ /* Wait if there is a pause between stimulation events */
+ if(stimulation_paused) {
+ if(count_pause == ISI) {
+ stimulation_paused = false;
+ count_pause = 0;
+ }
+ count_pause++;
+ }
+};
+
+mxArray* Stim::get_marker(void) {
+ extern const int red;
+ mxArray* Marker = mxCreateDoubleMatrix(0, 0, mxREAL);
+ mxSetM(Marker, 2);
+ mxSetN(Marker, marker_stimulation.size());
+ mxSetData(Marker, mxMalloc(sizeof(double)*2*marker_stimulation.size()));
+ double* Pr_Marker = mxGetPr(Marker);
+ for(unsigned i=0; i T */
+ double* var_stim = mxGetPr (prhs[4]); /* Parameters of stimulation protocol */
/* Initialize the populations */
- Cortical_Column Cortex(Param_Cortex, Connections);
- Thalamic_Column Thalamus(Connections);
+ Cortical_Column Cortex (Param_Cortex, Connections);
+ Thalamic_Column Thalamus (Param_Thalamus, Connections);
/* Link both modules */
Cortex.get_Thalamus(Thalamus);
diff --git a/Thalamic_Column.h b/Thalamic_Column.h
index 446d2a4..2c03175 100644
--- a/Thalamic_Column.h
+++ b/Thalamic_Column.h
@@ -56,8 +56,9 @@ class Thalamic_Column {
{set_RNG();}
/* Constructor for simulation */
- Thalamic_Column(double* Con)
- : N_et (Con[0]), N_er (Con[1])
+ Thalamic_Column(double* Param, double* Con)
+ : g_LK_t (Param[0]), g_LK_r(Param[1]), g_h(Param[2]),
+ N_et (Con[0]), N_er (Con[1])
{set_RNG();}
/* Get the pointer to the cortical module */
@@ -173,15 +174,15 @@ class Thalamic_Column {
const double g_L_r = 1;
/* Potassium leak current */
- const double g_LK_t = 0.025;
- const double g_LK_r = 0.025;
+ const double g_LK_t = 0.02;
+ const double g_LK_r = 0.02;
/* T current */
const double g_T_t = 3;
const double g_T_r = 2;
/* h current */
- const double g_h = 0.08;
+ const double g_h = 0.09;
/* Reversal potentials in mV */
/* Synaptic */
diff --git a/saves.h b/saves.h
index 852374c..3f95a57 100644
--- a/saves.h
+++ b/saves.h
@@ -24,7 +24,6 @@
/* Functions for data storage */
/****************************************************************************************************/
#pragma once
-#include
#include "Cortical_Column.h"
#include "Thalamic_Column.h"