From 790a3f1e58701edf36443947129bb5e15d0e1d7d Mon Sep 17 00:00:00 2001
From: Schellenberger <schellenberger@inb.uni-luebeck.de>
Date: Wed, 1 Oct 2014 10:39:36 +0200
Subject: [PATCH] Updated implementation of the Thalamo-cortical model
 Stimulation protocol has been reworked Multiple code cleanups and spelling
 fixes

Signed-off-by: Schellenberger <schellenberger@inb.uni-luebeck.de>
---
 .cproject         |   6 +-
 Cortical_Column.h |  12 +-
 Main.cpp          |   8 +-
 Plots.m           |  63 +++++----
 Stimulation.h     | 413 +++++++++++++++++++++++++++++++-----------------------
 TC.cpp            |  21 +--
 Thalamic_Column.h |  11 +-
 saves.h           |   1 -
 8 files changed, 307 insertions(+), 228 deletions(-)

diff --git a/.cproject b/.cproject
index 4957c40..bcd7e28 100644
--- a/.cproject
+++ b/.cproject
@@ -23,7 +23,7 @@
 							<tool id="cdt.managedbuild.tool.gnu.cpp.compiler.exe.debug.1024163887" name="GCC C++ Compiler" superClass="cdt.managedbuild.tool.gnu.cpp.compiler.exe.debug">
 								<option id="gnu.cpp.compiler.exe.debug.option.optimization.level.787584377" name="Optimization Level" superClass="gnu.cpp.compiler.exe.debug.option.optimization.level" value="gnu.cpp.compiler.optimization.level.most" valueType="enumerated"/>
 								<option id="gnu.cpp.compiler.exe.debug.option.debugging.level.312080211" name="Debug Level" superClass="gnu.cpp.compiler.exe.debug.option.debugging.level" value="gnu.cpp.compiler.debugging.level.max" valueType="enumerated"/>
-								<option id="gnu.cpp.compiler.option.other.other.74128570" name="Other flags" superClass="gnu.cpp.compiler.option.other.other" value="-std=c++0x -c -fmessage-length=0" valueType="string"/>
+								<option id="gnu.cpp.compiler.option.other.other.74128570" name="Other flags" superClass="gnu.cpp.compiler.option.other.other" value="-std=c++11 -c -fmessage-length=0" valueType="string"/>
 								<inputType id="cdt.managedbuild.tool.gnu.cpp.compiler.input.1054982019" superClass="cdt.managedbuild.tool.gnu.cpp.compiler.input"/>
 							</tool>
 							<tool id="cdt.managedbuild.tool.gnu.c.compiler.exe.debug.900081259" name="GCC C Compiler" superClass="cdt.managedbuild.tool.gnu.c.compiler.exe.debug">
@@ -33,7 +33,7 @@
 							</tool>
 							<tool id="cdt.managedbuild.tool.gnu.c.linker.exe.debug.2088031771" name="GCC C Linker" superClass="cdt.managedbuild.tool.gnu.c.linker.exe.debug"/>
 							<tool id="cdt.managedbuild.tool.gnu.cpp.linker.exe.debug.590332858" name="GCC C++ Linker" superClass="cdt.managedbuild.tool.gnu.cpp.linker.exe.debug">
-								<option id="gnu.cpp.link.option.flags.218353989" name="Linker flags" superClass="gnu.cpp.link.option.flags" value="-std=c++0x" valueType="string"/>
+								<option id="gnu.cpp.link.option.flags.218353989" name="Linker flags" superClass="gnu.cpp.link.option.flags" value="-std=c++11" valueType="string"/>
 								<inputType id="cdt.managedbuild.tool.gnu.cpp.linker.input.1214797218" superClass="cdt.managedbuild.tool.gnu.cpp.linker.input">
 									<additionalInput kind="additionalinputdependency" paths="$(USER_OBJS)"/>
 									<additionalInput kind="additionalinput" paths="$(LIBS)"/>
@@ -101,7 +101,7 @@
 		</cconfiguration>
 	</storageModule>
 	<storageModule moduleId="cdtBuildSystem" version="4.0.0">
-		<project id="TC_v2.cdt.managedbuild.target.gnu.exe.517165918" name="Executable" projectType="cdt.managedbuild.target.gnu.exe"/>
+		<project id="TC_model.cdt.managedbuild.target.gnu.exe.1203717936" name="Executable" projectType="cdt.managedbuild.target.gnu.exe"/>
 	</storageModule>
 	<storageModule moduleId="scannerConfiguration">
 		<autodiscovery enabled="true" problemReportingEnabled="true" selectedProfileId=""/>
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 <boost/random/uniform_int_distribution.hpp>
 #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<int>		marker_minimum;
+	vector<int>		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; i<marker_stimulation.size(); ++i) {
-			/* Division by res transforms marker time from sample rate to sec */
-			Pr_Marker[0+i*3] = marker_threshold[i]/res;
-			Pr_Marker[1+i*3] = marker_minimum[i]/res;
-			Pr_Marker[2+i*3] = marker_stimulation[i]/res;
-		}
-		return Marker;
-	}
+	/* Random number generator in case of semi-periodic stimulation */
+	ENG				Generator;
+	DIST_int		Uniform_Distribution;
+};
+/****************************************************************************************************/
+/*										 		end													*/
+/****************************************************************************************************/
 
-private:
-	/* Stimulation strength */
-	double strength 	= 0.0;
 
-	/* Duration of the stimulation */
-	int duration 		= 500;
+/****************************************************************************************************/
+/*										 Function definitions										*/
+/****************************************************************************************************/
+void Stim::setup (double* var_stim) {
+	extern const int onset;
+	extern const int res;
 
-	/* Inter-stimulus-interval in case of periodic stimulation */
-	int ISI				= 5E4;
+	/* Set the onset onset_correction for the marker */
+	onset_correction 		= onset * res;
 
-	/* Threshold for phase dependent stimulation */
-	double threshold	= -70;
+	/* Mode of stimulation */
+	mode					= (int) var_stim[0];
 
-	/* Time until stimulus after minimum was found */
-	int	time_to_stim	= 900;
+	/* Scale the stimulation strength from s^-1 (Hz) to ms^-1 */
+	strength 				= 		var_stim[1] / 1000;
 
-	/* Mode of stimulation 				*/
-	/* 0 == none 						*/
-	/* 1 == periodic 					*/
-	/* 2 == phase dependent up state 	*/
-	/* 3 == phase dependent down state 	*/
-	int mode			= 0;
+	/* Scale duration from ms to dt */
+	duration 				= (int) var_stim[2] * res / 1000;
 
-	/* Internal variables */
-	/* Simulation on for TRUE and off for FALSE */
-	bool stimulation_started= false;
+	/* Scale the inter stimulus event interval from s to dt */
+	ISI 					= (int) var_stim[3] * res;
 
-	/* Threshold has been reached */
-	bool threshold_crossed	= false;
+	/* Scale inter stimulus event interval range from s to dt */
+	ISI_range 				= (int) var_stim[4] * res;
 
-	/* Minimum found */
-	bool minimum_found		= false;
+	/* Number of stimuli per Stimulus event */
+	number_of_stimuli		= (int) var_stim[5];
 
-	/* Onset in time steps to correct the given time of the markers */
-	int correction			= 10000;
+	/* Scale time_between_stimuli from ms to dt */
+	time_between_stimuli 	= (int) var_stim[6] * res / 1000;
 
-	/* Counter for stimulation duration */
-	int count_duration		= 0;
+	/* If ISI is fixed do not create RNG */
+	if (mode == 1) {
+		/* Set first time_to_stimuli to 1 sec after onset */
+		time_to_stimuli = (int) (onset+1) * res;
 
-	/* Counter after minimum */
-	int count_to_start 		= 0;
+		/* If ISI is random create RNG */
+		if (ISI_range != 0){
+			/* Create the generator */
+			Generator = ENG(rand());
 
-	/* Old voltage value for minimum detection */
-	double Ve_old			= 0;
+			/* Combine RNG with distribution */
+			Uniform_Distribution = DIST_int(ISI-ISI_range, ISI+ISI_range);
+		}
+	} else {
+		/* In case of phase dependent stimulation, time_to_stim is the time from minimum detection to start of stimulation */
+		/* Scale time_to_stimuli from ms to dt */
+		time_to_stimuli = (int) var_stim[7] * res / 1000;
+	}
+};
 
-	/* Pointer to cortical column */
-	Cortical_Column* Cortex;
+void Stim::check_stim	(int time) {
+	/* Check if stimulation should start */
+	switch (mode) {
+
+	/* No stimulation */
+	default:
+		break;
+
+	/* Semi-periodic stimulation */
+	case 1:
+		/* Check if stimulation time is reached */
+		if(time == time_to_stimuli) {
+			/* Switch stimulation on */
+			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++;
+			}
+			/* 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<double>	marker_threshold;
-	std::vector<double>	marker_minimum;
-	std::vector<double>	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; i<number_of_stimuli; ++i) {
+					marker_minimum.push_back(time - onset_correction);
+				}
+				Ve_old = 0;
+			} else {
+				Ve_old = Cortex->Ve[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<marker_stimulation.size(); ++i) {
+		/* Division by res transforms marker time from dt to sampling rate */
+		Pr_Marker[0+i*2] = marker_minimum	  [i]/red;
+		Pr_Marker[1+i*2] = marker_stimulation [i]/red;
+	}
+	return Marker;
 };
+
 /****************************************************************************************************/
 /*										 		end													*/
 /****************************************************************************************************/
diff --git a/TC.cpp b/TC.cpp
index 54ae42e..bab0f8f 100644
--- a/TC.cpp
+++ b/TC.cpp
@@ -23,7 +23,7 @@
 /****************************************************************************************************/
 /* 	Implementation of the simulation as MATLAB routine (mex compiler)								*/
 /* 	mex command is given by:																		*/
-/* 	mex CXXFLAGS="\$CXXFLAGS -std=c++0x" TC.cpp Cortical_Column.cpp Thalamic_Column.cpp				*/
+/* 	mex CXXFLAGS="\$CXXFLAGS -std=c++11" TC.cpp Cortical_Column.cpp Thalamic_Column.cpp				*/
 /*	The Simulation requires the following boost libraries:	Preprocessor							*/
 /*															Random									*/
 /****************************************************************************************************/
@@ -39,11 +39,11 @@
 /****************************************************************************************************/
 /*										Fixed simulation settings									*/
 /****************************************************************************************************/
-extern const int onset	= 20;								/* time until data is stored in  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 onset	= 20;								/* Time until data is stored in  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 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			 										*/
 /****************************************************************************************************/
@@ -60,12 +60,13 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 	const int T				= (int) (mxGetScalar(prhs[0]));	/* Duration of simulation in s			*/
 	const int Time 			= (T+onset)*res;				/* Total number of iteration steps		*/
 	double* Param_Cortex	= mxGetPr (prhs[1]);			/* Parameters of cortical module		*/
-	double* Connections		= mxGetPr (prhs[2]);			/* Parameters of thalamic module		*/
-	double* var_stim	 	= mxGetPr (prhs[3]);			/* Parameters of stimulation protocol	*/
+	double* Param_Thalamus	= mxGetPr (prhs[2]);			/* Parameters of thalamic module		*/
+	double* Connections		= mxGetPr (prhs[3]);			/* Connectivity values C <-> 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 <iostream>
 #include "Cortical_Column.h"
 #include "Thalamic_Column.h"