diff --git a/Main.cpp b/Main.cpp
index 1ae8668..9095b5e 100644
--- a/Main.cpp
+++ b/Main.cpp
@@ -1,3 +1,25 @@
+/*
+*	Copyright (c) 2014 Michael Schellenberger Costa
+*
+*	Permission is hereby granted, free of charge, to any person obtaining a copy
+*	of this software and associated documentation files (the "Software"), to deal
+*	in the Software without restriction, including without limitation the rights
+*	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+*	copies of the Software, and to permit persons to whom the Software is
+*	furnished to do so, subject to the following conditions:
+*
+*	The above copyright notice and this permission notice shall be included in
+*	all copies or substantial portions of the Software.
+*
+*	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+*	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+*	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+*	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+*	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+*	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+*	THE SOFTWARE.
+*/
+
 /****************************************************************************************************/
 /*		Main file for compilation tests																*/
 /*		The Simulation requires the following boost libraries:	Preprocessor						*/
@@ -11,10 +33,11 @@
 /****************************************************************************************************/
 /*										Fixed simulation settings									*/
 /****************************************************************************************************/
-extern const int T 		= 60;
-extern const int res 	= 1E4;
-extern const double dt 	= 1E3/res;
-extern const double h	= sqrt(dt);
+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	*/
 /****************************************************************************************************/
 /*										 		end			 										*/
 /****************************************************************************************************/
@@ -24,20 +47,20 @@ extern const double h	= sqrt(dt);
 /*										Main simulation routine										*/
 /****************************************************************************************************/
 int main(void) {
-	// Initializing the populations;
+	/* Initialize the populations */
 	Thalamic_Column Thalamus;
 
-	// takes the time of the simulation
+	/* Takes the time of the simulation */
 	time_t start,end;
 	time (&start);
 
-	// simulation
+	/* Simulation */
 	for (int t=0; t< T*res; ++t) {
 		ODE (Thalamus);
 	}
 	time (&end);
 
-	// time consumed by the simulation
+	/* Time consumed by the simulation */
 	double dif = difftime(end,start);
 	std::cout << "simulation done!\n";
 	std::cout << "took " << dif 	<< " seconds" << "\n";
diff --git a/ODE.h b/ODE.h
index 607eb48..cb2165f 100644
--- a/ODE.h
+++ b/ODE.h
@@ -1,3 +1,25 @@
+/*
+*	Copyright (c) 2014 Michael Schellenberger Costa
+*
+*	Permission is hereby granted, free of charge, to any person obtaining a copy
+*	of this software and associated documentation files (the "Software"), to deal
+*	in the Software without restriction, including without limitation the rights
+*	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+*	copies of the Software, and to permit persons to whom the Software is
+*	furnished to do so, subject to the following conditions:
+*
+*	The above copyright notice and this permission notice shall be included in
+*	all copies or substantial portions of the Software.
+*
+*	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+*	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+*	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+*	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+*	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+*	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+*	THE SOFTWARE.
+*/
+
 /****************************************************************************************************/
 /*									Implementation of the ODE solver								*/
 /****************************************************************************************************/
@@ -8,11 +30,12 @@
 /*										Evaluation of SRK4											*/
 /****************************************************************************************************/
 inline void ODE(Thalamic_Column& Thalamus) {
-	// first calculating every ith RK moment
-	// has to be in order, 1th moment first
+	/* First calculating every ith RK moment. Has to be in order, 1th moment first */
 	for (int i=1; i<=4; ++i) {
 		Thalamus.set_RK(i);
 	}
+
+	/* Add moments */
 	Thalamus.add_RK();
 }
 /****************************************************************************************************/
diff --git a/Stimulation.h b/Stimulation.h
index fccdf15..b00fa70 100644
--- a/Stimulation.h
+++ b/Stimulation.h
@@ -1,3 +1,25 @@
+/*
+*	Copyright (c) 2014 Michael Schellenberger Costa
+*
+*	Permission is hereby granted, free of charge, to any person obtaining a copy
+*	of this software and associated documentation files (the "Software"), to deal
+*	in the Software without restriction, including without limitation the rights
+*	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+*	copies of the Software, and to permit persons to whom the Software is
+*	furnished to do so, subject to the following conditions:
+*
+*	The above copyright notice and this permission notice shall be included in
+*	all copies or substantial portions of the Software.
+*
+*	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+*	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+*	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+*	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+*	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+*	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+*	THE SOFTWARE.
+*/
+
 /****************************************************************************************************/
 /*								Implementation of the stimulation protocol							*/
 /****************************************************************************************************/
@@ -9,32 +31,32 @@
 /****************************************************************************************************/
 class Stim {
 public:
-	/* empty constructor for compiling */
+	/* Empty constructor for compiling */
 	Stim(void);
 
 	Stim(Thalamic_Column& T, double* var)
 	{ Thalamus = &T;
 	  setup(var);}
 
-	/* setup with respect to stimulation mode */
+	/* Setup with respect to stimulation mode */
 	void setup		(double* var_stim) {
 		extern const int onset;
 		extern const int res;
 		extern const int dt;
 
-		/* mode of stimulation */
+		/* Mode of stimulation */
 		mode		= (int) var_stim[0];
 
-		/* scale the stimulation strength from s^-1 to ms^-1 */
+		/* Scale the stimulation strength from s^-1 to ms^-1 */
 		strength 	= 		var_stim[1] / 1000;
 
-		/* scale duration from ms to dt */
+		/* Scale duration from ms to dt */
 		duration 	= (int) var_stim[2] * res / 1000;
 
-		/* scale the ISI from s to ms^-1 */
+		/* Scale the ISI from s to ms^-1 */
 		ISI 		= (int) var_stim[3] * res;
 
-		/* scale time to stimulus from ms to dt */
+		/* Scale time to stimulus from ms to dt */
 		time_to_stim= (int) var_stim[4] * res / 1000;
 
 		if(mode==1) {
@@ -47,33 +69,35 @@ class Stim {
 
 	void check_stim	(int time) {
 
-		/* check if stimulation should start */
+		/* Check if stimulation should start */
 		switch (mode) {
 
-		/* no stimulation */
+		/* No stimulation */
 		default:
 			break;
 
-		/* periodic stimulation */
+		/* Periodic stimulation */
 		case 1:
-			/* check if time is reached */
+			/* Check if time is reached */
 			if(time == time_to_stim) {
-				/* switch stimulation on */
+				/* Switch stimulation on */
 				stimulation_started 	= true;
 				Thalamus->set_input(strength);
 
-				/* update the timer */
+				/* Update the timer */
 				time_to_stim += ISI;
+
+				/* Update the marker */
+				marker_stimulation.push_back(time - correction);
 			}
 			break;
-
 		}
 
-		/* wait to switch the stimulation off */
+		/* Wait to switch the stimulation off */
 		if(stimulation_started) {
 			count_duration++;
 
-			/* switch stimulation off */
+			/* Switch stimulation off */
 			if(count_duration==duration) {
 				stimulation_started 	= false;
 				count_duration			= 0;
@@ -82,22 +106,38 @@ class Stim {
 		}
 	}
 
+	/* Create MATLAB container for marker storage */
+	mxArray* get_marker(void) {
+		mxArray* Marker	= mxCreateDoubleMatrix(0, 0, mxREAL);
+	    mxSetM(Marker, 1);
+	    mxSetN(Marker, marker_stimulation.size());
+	    mxSetData(Marker, mxMalloc(sizeof(double)*marker_stimulation.size()));
+		double* Pr_Marker	= mxGetPr(Marker);
+		for(unsigned i=0; i<marker_stimulation.size(); ++i) {
+			Pr_Marker[i] = marker_stimulation[i];
+		}
+		return Marker;
+	}
+
 private:
 
 	/* Stimulation parameters */
-	/* stimulation strength */
+	/* Stimulation strength */
 	double strength 	= 0.0;
 
-	/* duration of the stimulation */
+	/* Duration of the stimulation */
 	int duration 		= 500;
 
-	/* inter stimulus intervall in case of periodic stimulation */
+	/* Inter stimulus intervall in case of periodic stimulation */
 	int ISI				= 5E4;
 
-	/* time until stimulus after minimum was found */
+	/* Threshold for phase dependent stimulation */
+	double threshold	= -80;
+
+	/* Time until stimulus after minimum was found */
 	int	time_to_stim	= 5500;
 
-	/* mode of stimulation 				*/
+	/* Mode of stimulation 				*/
 	/* 0 == none 						*/
 	/* 1 == periodic 					*/
 	int mode			= 0;
@@ -106,18 +146,23 @@ class Stim {
 	/* Simulation on for TRUE and off for FALSE */
 	bool stimulation_started= false;
 
-	/* onset in timesteps to correct the given time of the markers */
+	/* Threshold has been reached */
+	bool threshold_crossed	= false;
+
+	/* Minimum found */
+	bool minimum_found		= false;
+
+	/* Onset in timesteps to correct the given time of the markers */
 	int correction			= 10000;
 
-	/* counter for stimulation duration */
+	/* Counter for stimulation duration */
 	int count_duration		= 0;
 
-	/* counter after minimum */
-	int count_to_start 		= 0;
-
-	/* pointer to thalamic column */
+	/* Pointer to thalamic column */
 	Thalamic_Column* Thalamus;
 
+	/* Data containers */
+	std::vector<int>	marker_stimulation;
 };
 /****************************************************************************************************/
 /*										 		end													*/
diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp
index 2fc8c32..62c41b2 100644
--- a/Thalamic_Column.cpp
+++ b/Thalamic_Column.cpp
@@ -1,3 +1,25 @@
+/*
+*	Copyright (c) 2014 Michael Schellenberger Costa
+*
+*	Permission is hereby granted, free of charge, to any person obtaining a copy
+*	of this software and associated documentation files (the "Software"), to deal
+*	in the Software without restriction, including without limitation the rights
+*	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+*	copies of the Software, and to permit persons to whom the Software is
+*	furnished to do so, subject to the following conditions:
+*
+*	The above copyright notice and this permission notice shall be included in
+*	all copies or substantial portions of the Software.
+*
+*	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+*	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+*	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+*	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+*	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+*	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+*	THE SOFTWARE.
+*/
+
 /****************************************************************************************************/
 /*									Functions of the cortical module								*/
 /****************************************************************************************************/
@@ -7,15 +29,15 @@
 /*										 Initialization of RNG 										*/
 /****************************************************************************************************/
 void Thalamic_Column::set_RNG(void) {
-	// the number of independent streams
+	/* Number of independent streams */
 	int N = 2;
 
-	// get the RNG
+	/* Create RNG for each stream */
 	for (int i=0; i<N; ++i){
-		// add the RNG
+		/* Add the RNG */
 		MTRands.push_back({ENG(rand()), DIST (mphi, dphi)});
 
-		// get the random number for the first iteration
+		/* Get the random number for the first iteration */
 		Rand_vars.push_back(MTRands[i]());
 	}
 }
@@ -27,14 +49,14 @@ void Thalamic_Column::set_RNG(void) {
 /****************************************************************************************************/
 /*										 Firing Rate functions 										*/
 /****************************************************************************************************/
-// thalamic relay (TC) firing rate
+/* Thalamic relay (TC) firing rate */
 double Thalamic_Column::get_Qt	(int N) const{
 	_SWITCH((Vt))
 	double q = Qt_max/ (1 + exp(-C1 * (var_Vt - theta_t) / sigma_t));
 	return q;
 }
 
-// thalamic reticular (RE) firing rate
+/* Thalamic reticular (RE) firing rate */
 double Thalamic_Column::get_Qr	(int N) const{
 	_SWITCH((Vr))
 	double q = Qr_max / (1 + exp(-C1 * (var_Vr - theta_r) / sigma_r));
@@ -48,27 +70,27 @@ double Thalamic_Column::get_Qr	(int N) const{
 /****************************************************************************************************/
 /*										Synaptic currents											*/
 /****************************************************************************************************/
-// excitatory input to TC population
+/* Excitatory input to TC population */
 double Thalamic_Column::I_et	(int N) const{
 	_SWITCH((Vt)(Phi_tt))
 	double psi = var_Phi_tt * (var_Vt - E_AMPA);
 	return psi;
 }
 
-// inhibitory input to TC population
+/* Inhibitory input to TC population */
 double Thalamic_Column::I_it	(int N) const{
 	_SWITCH((Vt)(Phi_rt))
 	double psi = var_Phi_rt * (var_Vt - E_GABA);
 	return psi;
 }
-// excitatory input to RE population
+/* Excitatory input to RE population */
 double Thalamic_Column::I_er	(int N) const{
 	_SWITCH((Vr)(Phi_tr))
 	double psi = var_Phi_tr * (var_Vr - E_AMPA);
 	return psi;
 }
 
-// inhibitory input to RE population
+/* Inhibitory input to RE population */
 double Thalamic_Column::I_ir	(int N) const{
 	_SWITCH((Vr)(Phi_rr))
 	double psi = var_Phi_rr * (var_Vr - E_GABA);
@@ -82,64 +104,56 @@ double Thalamic_Column::I_ir	(int N) const{
 /****************************************************************************************************/
 /*										 	I_T gating	 											*/
 /****************************************************************************************************/
-// instant activation in TC population
-// after Destexhe 1996
+/* Activation in TC population after Destexhe 1996 */
 double Thalamic_Column::m_inf_T_t	(int N) const{
 	_SWITCH((Vt))
 	double m = 1/(1+exp(-(var_Vt+59)/6.2));
 	return m;
 }
 
-// instant activation in RE population
-// after Destexhe 1996
+/* Activation in RE population after Destexhe 1996 */
 double Thalamic_Column::m_inf_T_r	(int N) const{
 	_SWITCH((Vr))
 	double m = 1/(1+exp(-(var_Vr+52)/7.2));
 	return m;
 }
 
-// instant deactivation in TC population
-// after Destexhe 1996
+/* Deactivation in TC population after Destexhe 1996 */
 double Thalamic_Column::h_inf_T_t	(int N) const{
 	_SWITCH((Vt))
 	double h = 1/(1+exp( (var_Vt+81)/4));
 	return h;
 }
 
-// instant deactivation in RE population
-// after Destexhe 1996
+/* Deactivation in RE population after Destexhe 1996 */
 double Thalamic_Column::h_inf_T_r	(int N) const{
 	_SWITCH((Vr))
 	double h = 1/(1+exp( (var_Vr+80)/5));
 	return h;
 }
 
-// activation time in TC population
-// after Bazhenov 1998
+/* Activation time in TC population after Bazhenov 1998 */
 double Thalamic_Column::tau_m_T_t	(int N) const{
 	_SWITCH((Vt))
 	double tau = (0.612 + 1 /(exp(-(var_Vt + 132) / 16.7) + exp((var_Vt + 16.8)/18.2)))/pow(3,1.2);
 	return tau;
 }
 
-// activation time in RE population
-// after Destexhe 1996
+/* Activation time in RE population after Destexhe 1996 */
 double Thalamic_Column::tau_m_T_r	(int N) const{
 	_SWITCH((Vr))
 	double tau = (1 + 0.33/( exp((var_Vr+27)/10.0) + exp(-(var_Vr+102)/15.0)))/pow(2.5, 1.2);
 	return tau;
 }
 
-// deactivation time in RE population
-// after Destexhe 1996
+/* deactivation time in RE population after Destexhe 1996 */
 double Thalamic_Column::tau_h_T_t	(int N) const{
 	_SWITCH((Vt))
 	double tau =  (30.8 + (211.4 + exp((var_Vt+115.2)/5))/(1 + exp((var_Vt+86)/3.2)))/pow(3, 1.2);
 	return tau;
 }
 
-// deactivation time in RE population
-// after Destexhe 1996
+/* Deactivation time in RE population after Destexhe 1996 */
 double Thalamic_Column::tau_h_T_r	(int N) const{
 	_SWITCH((Vr))
 	double tau =  (85 + 1/(exp((var_Vr+48)/4) + exp(-(var_Vr+407)/50)))/pow(3, 1.2);
@@ -153,16 +167,14 @@ double Thalamic_Column::tau_h_T_r	(int N) const{
 /****************************************************************************************************/
 /*										 	I_h gating 												*/
 /****************************************************************************************************/
-// instant activation in TC population
-// after Destexhe 1993
+/* Activation in TC population after Destexhe 1993 */
 double Thalamic_Column::m_inf_h	(int N) const{
 	_SWITCH((Vt))
 	double h = 1/(1+exp( (var_Vt+75)/5.5));
 	return h;
 }
 
-// activation time for slow components in TC population
-// after Destexhe 1993
+/* Activation time for slow components in TC population after Destexhe 1993 */
 double Thalamic_Column::tau_m_h	(int N) const{
 	_SWITCH((Vt))
 	double tau = 1. / (exp(-14.59 - 0.086 * var_Vt) + exp(-1.87 + 0.07 * var_Vt));
@@ -176,51 +188,51 @@ double Thalamic_Column::tau_m_h	(int N) const{
 /****************************************************************************************************/
 /*										 Current functions 											*/
 /****************************************************************************************************/
-// Leak current of TC population
+/* Leak current of TC population */
 double Thalamic_Column::I_L_t	(int N) const{
 	_SWITCH((Vt))
 	double I = g_L_t * (var_Vt - E_L_t);
 	return I;
 }
 
-// Potassium leak current of TC population
+/* Potassium leak current of TC population */
 double Thalamic_Column::I_LK_t	(int N) const{
 	_SWITCH((Vt))
 	double I = g_LK_t * (var_Vt - E_K);
 	return I;
 }
 
-// Leak current of RE population
+/* Leak current of RE population */
 double Thalamic_Column::I_L_r	(int N) const{
 	_SWITCH((Vr))
 	double I = g_L_r * (var_Vr - E_L_r);
 	return I;
 }
 
-// Potassium leak current of RE population
+/* Potassium leak current of RE population */
 double Thalamic_Column::I_LK_r	(int N) const{
 	_SWITCH((Vr))
 	double I = g_LK_r * (var_Vr - E_K);
 	return I;
 }
 
-// T-type current of TC population
+/* T-type current of TC population */
 double Thalamic_Column::I_T_t	(int N) const{
-	_SWITCH((Vt)(h_T_t))
-	//double I = g_T_t * pow(var_m_T_t, 2) * var_h_T_t * (var_Vt - E_Ca);
-	double I = g_T_t * pow(m_inf_T_t(N), 2) * var_h_T_t * (var_Vt - E_Ca);
+	_SWITCH((Vt)(h_T_t)(m_T_t))
+	double I = g_T_t * pow(var_m_T_t, 2) * var_h_T_t * (var_Vt - E_Ca);
+	/* double I = g_T_t * pow(m_inf_T_t(N), 2) * var_h_T_t * (var_Vt - E_Ca); */
 	return I;
 }
 
-// T-type current of RE population
+/* T-type current of RE population */
 double Thalamic_Column::I_T_r	(int N) const{
 	_SWITCH((Vr)(h_T_r)(m_T_r))
 	double I = g_T_r * pow(var_m_T_r, 2) * var_h_T_r * (var_Vr - E_Ca);
-	//double I = g_T_r * pow(m_inf_T_r(N), 2) * var_h_T_r * (var_Vr - E_Ca);
+	/* double I = g_T_r * pow(m_inf_T_r(N), 2) * var_h_T_r * (var_Vr - E_Ca); */
 	return I;
 }
 
-// h-type current of TC population
+/* h-type current of TC population */
 double Thalamic_Column::I_h		(int N) const{
 	_SWITCH((Vt)(m_h)(m_h2))
 	double I = g_h * (var_m_h + g_inc * var_m_h2) * (var_Vt - E_h);
@@ -269,7 +281,7 @@ void Thalamic_Column::set_RK (int N) {
 	Phi_tr	[N] = dt*(var_x_tr);
 	Phi_rt	[N] = dt*(var_x_rt);
 	Phi_rr	[N] = dt*(var_x_rr);
-	x_tt  	[N] = dt*(pow(gamma_e, 2) * (noise_xRK(N) 		- var_Phi_tt) - 2 * gamma_e * var_x_tt);
+	x_tt  	[N] = dt*(pow(gamma_e, 2) * (noise_xRK(N) 	 	- var_Phi_tt) - 2 * gamma_e * var_x_tt);
 	x_tr  	[N] = dt*(pow(gamma_e, 2) * (N_tr * get_Qt(N)	- var_Phi_tr) - 2 * gamma_e * var_x_tr);
 	x_rt  	[N] = dt*(pow(gamma_i, 2) * (N_rt * get_Qr(N) 	- var_Phi_rt) - 2 * gamma_i * var_x_rt);
 	x_rr  	[N] = dt*(pow(gamma_i, 2) * (N_rr * get_Qr(N)	- var_Phi_rr) - 2 * gamma_i * var_x_rr);
@@ -303,7 +315,7 @@ void Thalamic_Column::add_RK(void) {
 	m_h2	[0] += (m_h2	[1] + m_h2		[2] * 2 + m_h2		[3] * 2 + m_h2		[4])/6;
 	P_h	 	[0] += (P_h		[1] + P_h		[2] * 2 + P_h		[3] * 2 + P_h		[4])/6;
 
-	// generating the noise for the next iteration
+	/* Generate noise for the next iteration */
 	for (unsigned i=0; i<Rand_vars.size(); ++i) {
 		Rand_vars[i] = MTRands[i]() + input;
 	}
diff --git a/Thalamic_Column.h b/Thalamic_Column.h
index 796c608..b2e3742 100644
--- a/Thalamic_Column.h
+++ b/Thalamic_Column.h
@@ -1,3 +1,25 @@
+/*
+*	Copyright (c) 2014 Michael Schellenberger Costa
+*
+*	Permission is hereby granted, free of charge, to any person obtaining a copy
+*	of this software and associated documentation files (the "Software"), to deal
+*	in the Software without restriction, including without limitation the rights
+*	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+*	copies of the Software, and to permit persons to whom the Software is
+*	furnished to do so, subject to the following conditions:
+*
+*	The above copyright notice and this permission notice shall be included in
+*	all copies or substantial portions of the Software.
+*
+*	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+*	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+*	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+*	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+*	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+*	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+*	THE SOFTWARE.
+*/
+
 /************************************************************************************************/
 /*								Header file of a thalamic module								*/
 /************************************************************************************************/
@@ -13,9 +35,9 @@ using std::vector;
 /****************************************************************************************************/
 /*										Typedefs for RNG											*/
 /****************************************************************************************************/
-typedef boost::mt11213b                    	ENG;    // Mersenne Twister
-typedef boost::normal_distribution<double>	DIST;   // Normal Distribution
-typedef boost::variate_generator<ENG,DIST> 	GEN;    // Variate generator
+typedef boost::mt11213b                    	ENG;    /* Mersenne Twister		*/
+typedef boost::normal_distribution<double>	DIST;   /* Normal Distribution	*/
+typedef boost::variate_generator<ENG,DIST> 	GEN;    /* Variate generator	*/
 /****************************************************************************************************/
 /*										 		end			 										*/
 /****************************************************************************************************/
@@ -26,32 +48,32 @@ typedef boost::variate_generator<ENG,DIST> 	GEN;    // Variate generator
 /****************************************************************************************************/
 class Thalamic_Column {
 public:
-	// Constructors
+	/* Constructors for compile check */
 	Thalamic_Column(void)
 	{set_RNG();}
 
-	// Constructors
-	Thalamic_Column(double* Par)
-	: N_tr		(Par[0]),	N_rt		(Par[1]),	N_rr		(Par[2])
+	/* Constructor for simulation */
+	Thalamic_Column(double* Con)
+	: N_tr		(Con[0]),	N_rt		(Con[1]),	N_rr		(Con[2])
 	{set_RNG();}
 
-	// change the strength of input
+	/* Set strength of input */
 	void	set_input	(double I) {input = I;}
 
-	// Initialize the RNGs
+	/* Initialize the RNGs */
 	void 	set_RNG		(void);
 
-	// Firing rates
+	/* Firing rates */
 	double 	get_Qt		(int) const;
 	double 	get_Qr		(int) const;
 
-	// Synaptic currents
+	/* Synaptic currents */
 	double 	I_et		(int) const;
 	double 	I_it		(int) const;
 	double 	I_er		(int) const;
 	double 	I_ir		(int) const;
 
-	// Activation functions
+	/* Activation functions */
 	double  m_inf_T_t	(int) const;
 	double  m_inf_T_r	(int) const;
 	double  tau_m_T_t	(int) const;
@@ -59,13 +81,13 @@ class Thalamic_Column {
 	double  m_inf_h		(int) const;
 	double  tau_m_h		(int) const;
 
-	// Deactivation functions
+	/* 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
+	/* Currents */
 	double 	I_L_t		(int) const;
 	double 	I_L_r		(int) const;
 	double 	I_LK_t		(int) const;
@@ -74,113 +96,108 @@ class Thalamic_Column {
 	double 	I_T_r		(int) const;
 	double 	I_h			(int) const;
 
-	// Noise functions
+	/* Noise functions */
 	double 	noise_xRK 	(int) const;
 
-	// ODE functions
+	/* ODE functions */
 	void 	set_RK		(int);
 	void 	add_RK	 	(void);
 
-	// Data storage
+	/* Data storage  access */
 	friend void get_data (int, Thalamic_Column&, _REPEAT(double*, 2));
 
 private:
-	// Population variables
-	vector<double> 	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
-					Phi_tt	= _INIT(0.0),		// PostSP from TC population to TC population
-					Phi_tr	= _INIT(0.0),		// PostSP from TC population to RE population
-					Phi_rt	= _INIT(0.0),		// PostSP from RE population to TC population
-					Phi_rr	= _INIT(0.0),		// PostSP from RE population to RE population
-					x_tt	= _INIT(0.0),		// derivative of Phi_tt
-					x_tr	= _INIT(0.0),		// derivative of Phi_tr
-					x_rt	= _INIT(0.0),		// derivative of Phi_rt
-					x_rr	= _INIT(0.0),		// derivative of Phi_rr
-					h_T_t	= _INIT(0.0),		// inactivation of T channel
-					h_T_r	= _INIT(0.0),		// inactivation of T channel
-					m_T_t	= _INIT(0.0),		// activation 	of T channel
-					m_T_r	= _INIT(0.0),		// activation 	of T channel
-					m_h		= _INIT(0.0),		// activation 	of h   channel
-					m_h2	= _INIT(0.0),		// activation 	of h   channel bound with protein
-					P_h		= _INIT(0.0);		// fraction of protein bound with calcium
-
-	// Random number generators
+	/* Population variables */
+	vector<double> 	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			*/
+					Phi_tt	= _INIT(0.0),		/* PostSP from TC population to TC population		*/
+					Phi_tr	= _INIT(0.0),		/* PostSP from TC population to RE population		*/
+					Phi_rt	= _INIT(0.0),		/* PostSP from RE population to TC population		*/
+					Phi_rr	= _INIT(0.0),		/* PostSP from RE population to RE population		*/
+					x_tt	= _INIT(0.0),		/* derivative of Phi_tt								*/
+					x_tr	= _INIT(0.0),		/* derivative of Phi_tr								*/
+					x_rt	= _INIT(0.0),		/* derivative of Phi_rt								*/
+					x_rr	= _INIT(0.0),		/* derivative of Phi_rr								*/
+					h_T_t	= _INIT(0.0),		/* inactivation of T channel						*/
+					h_T_r	= _INIT(0.0),		/* inactivation of T channel						*/
+					m_T_t	= _INIT(0.0),		/* activation 	of T channel						*/
+					m_T_r	= _INIT(0.0),		/* activation 	of T channel						*/
+					m_h		= _INIT(0.0),		/* activation 	of h   channel						*/
+					m_h2	= _INIT(0.0),		/* activation 	of h   channel bound with protein 	*/
+					P_h		= _INIT(0.0);		/* fraction of protein bound with calcium 			*/
+
+	/* Random number generators */
 	vector<GEN>		MTRands;
 
-	// Container for noise
+	/* Container for noise */
 	vector<double>	Rand_vars;
 
-	// Declaration and Initialization of parameters
-	// Membrane time in ms
-	const double 	tau_t 		= 1;
-	const double 	tau_r 		= 1;
+	/* Declaration and Initialization of parameters */
+	/* Membrane time in ms */
+	const double 	tau_t 		= 30;
+	const double 	tau_r 		= 30;
 
-	// Maximum firing rate in ms^-1
+	/* Maximum firing rate in ms^-1 */
 	const double 	Qt_max		= 400.E-3;
 	const double 	Qr_max		= 400.E-3;
 
-	// Sigmoid threshold in mV
+	/* Sigmoid threshold in mV */
 	const double 	theta_t		= -45;
 	const double 	theta_r		= -45;
 
-	// Sigmoid gain in mV
-	const double 	sigma_t		= 9;
-	const double 	sigma_r		= 9;
+	/* Sigmoid gain in mV */
+	const double 	sigma_t		= 3;
+	const double 	sigma_r		= 3;
 
-	// Scaling parameter for sigmoidal mapping (dimensionless)
+	/* Scaling parameter for sigmoidal mapping (dimensionless) */
 	const double 	C1          = (3.14159265/sqrt(3));
 
-	// PSP rise time in ms^-1
+	/* PSP rise time in ms^-1 */
 	const double 	gamma_e		= 70E-3;
 	const double 	gamma_i		= 58.6E-3;
 
 	/* axonal flux time constant */
-	const double 	nu			= 60E-3;
+	const double 	nu			= 120E-3;
 
-	// Conductivities in mS/cm^-2
-	// Leak current
-	const double 	g_L_t  		= 0.02;
-	const double 	g_L_r  		= 0.05;
+	/* Conductivities in mS/cm^-2 */
+	/* Leak current */
+	const double 	g_L_t  		= 0.6;
+	const double 	g_L_r  		= 1.5;
 
-	// Potassium leak current
+	/* Potassium leak current */
 	const double 	g_LK_t 		= 0.02;
 	const double 	g_LK_r 		= 0.02;
 
-	// T current
-	const double	g_T_t		= 3;
+	/* T current */
+	const double	g_T_t		= 2.2;
 	const double	g_T_r		= 2;
 
 	/* h current */
 	const double	g_h			= 0.07;
 
-	// Connectivities (dimensionless)
-	const double 	N_tr		= 2;
-	const double 	N_rt		= 5.5;
-	const double 	N_rr		= 5;
-
-	// Reversal potentials in mV
-	// synaptic
+	/* Reversal potentials in mV */
+	/* Synaptic */
 	const double 	E_AMPA  	= 0;
 	const double 	E_GABA  	= -70;
 
-	// Leak
-	const double 	E_L_t 		= -70;
-	const double 	E_L_r 		= -55;
+	/* Leak */
+	const double 	E_L_t 		= -64;
+	const double 	E_L_r 		= -64;
 
-	// Potassium
+	/* Potassium */
 	const double 	E_K    		= -100;
 
-	// I_T current
+	/* I_T current */
 	const double 	E_Ca    	= 120;
 
-	// I_h current
+	/* I_h current */
 	const double 	E_h    		= -40;
 
 	/* Calcium parameters */
-	const double	alpha_Ca	= -50E-6;			// influx per spike in nmol
-	const double	tau_Ca		= 5;				// calcium time constant in ms
-	const double	Ca_0		= 2E-4;				// resting concentration
+	const double	alpha_Ca	= -50E-6;			/* influx per spike in nmol		*/
+	const double	tau_Ca		= 5;				/* calcium time constant in ms	*/
+	const double	Ca_0		= 2E-4;				/* resting concentration 		*/
 
 	/* I_h activation parameters */
 	const double 	k1			= 2.5E7;
@@ -190,11 +207,16 @@ class Thalamic_Column {
 	const double 	n_P			= 4;
 	const double 	g_inc		= 2;
 
-	// Noise parameters in ms^-1
+	/* Noise parameters in ms^-1 */
 	const double 	mphi		= 0E-3;
 	const double	dphi		= 2E-3;;
 	double			input		= 0.0;
 
+	/* Connectivities (dimensionless) */
+	const double 	N_tr		= 2;
+	const double 	N_rt		= 5.5;
+	const double 	N_rr		= 5;
+
 	friend class 	Stim;
 };
 /****************************************************************************************************/
diff --git a/Thalamus.cpp b/Thalamus.cpp
index 795a86c..f1ed9d0 100644
--- a/Thalamus.cpp
+++ b/Thalamus.cpp
@@ -1,3 +1,25 @@
+/*
+*	Copyright (c) 2014 Michael Schellenberger Costa
+*
+*	Permission is hereby granted, free of charge, to any person obtaining a copy
+*	of this software and associated documentation files (the "Software"), to deal
+*	in the Software without restriction, including without limitation the rights
+*	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+*	copies of the Software, and to permit persons to whom the Software is
+*	furnished to do so, subject to the following conditions:
+*
+*	The above copyright notice and this permission notice shall be included in
+*	all copies or substantial portions of the Software.
+*
+*	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+*	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+*	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+*	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+*	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+*	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+*	THE SOFTWARE.
+*/
+
 /****************************************************************************************************/
 /* 		Implementation of the simulation as MATLAB routine (mex compiler)							*/
 /* 		mex command is given by:																	*/
@@ -16,11 +38,11 @@
 /****************************************************************************************************/
 /*										Fixed simulation settings									*/
 /****************************************************************************************************/
-extern const int onset	= 5;								// 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	= 5;								/* 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	*/
 /****************************************************************************************************/
 /*										 		end			 										*/
 /****************************************************************************************************/
@@ -32,30 +54,30 @@ extern const double h	= sqrt(dt);							// squareroot of dt for SRK iteration
 /*										rhs defines inputs											*/
 /****************************************************************************************************/
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
-	// Initializing the seeder.
+	/* Initialize the seeder */
 	srand(time(NULL));
 
-	// Fetch inputs
-	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_Thalamus	= mxGetPr (prhs[1]);			// parameters of thalamic module
-	double* var_stim	 	= mxGetPr (prhs[2]);			// parameters of stimulation protocol
+	/* Fetch inputs */
+	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_Thalamus	= mxGetPr (prhs[1]);			/* Parameters of cortical module		*/
+	double* var_stim	 	= mxGetPr (prhs[2]);			/* Parameters of stimulation protocol	*/
 
-	// Initializing the populations;
+	/* Initialize the populations */
 	Thalamic_Column Thalamus(Param_Thalamus);
 
-	// Initialize the stimulation protocol
-	Stim	Stimulation(Thalamus, var_stim);
+	/* Initialize the stimulation protocol */
+	Stim Stimulation(Thalamus, var_stim);
 
-	// setting up the data containers
+	/* Create data containers */
 	mxArray* Vt		= SetMexArray(1, T*res/red);
 	mxArray* Vr		= SetMexArray(1, T*res/red);
 
-	// Pointer to the actual data block
+	/* Pointer to the actual data block */
 	double* Pr_Vt	= mxGetPr(Vt);
 	double* Pr_Vr	= mxGetPr(Vr);
 
-	// simulation
+	/* Simulation */
 	int count = 0;
 	for (int t=0; t<Time; ++t) {
 		ODE (Thalamus);
@@ -66,7 +88,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 		}
 	}
 
-	// output of the simulation
+	/* Output of the simulation */
 	plhs[0] = Vt;
 	plhs[1] = Vr;
 return;
diff --git a/macros.h b/macros.h
index 8fcf5b4..1baa7a1 100644
--- a/macros.h
+++ b/macros.h
@@ -1,3 +1,25 @@
+/*
+*	Copyright (c) 2014 Michael Schellenberger Costa
+*
+*	Permission is hereby granted, free of charge, to any person obtaining a copy
+*	of this software and associated documentation files (the "Software"), to deal
+*	in the Software without restriction, including without limitation the rights
+*	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+*	copies of the Software, and to permit persons to whom the Software is
+*	furnished to do so, subject to the following conditions:
+*
+*	The above copyright notice and this permission notice shall be included in
+*	all copies or substantial portions of the Software.
+*
+*	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+*	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+*	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+*	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+*	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+*	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+*	THE SOFTWARE.
+*/
+
 /****************************************************************************************************/
 /*									Definition of all macros used									*/
 /****************************************************************************************************/
@@ -9,6 +31,7 @@
 #include <boost/preprocessor/seq/for_each.hpp>
 #include <boost/preprocessor/facilities/empty.hpp>
 
+
 /****************************************************************************************************/
 /*									Macro for vector initialization									*/
 /****************************************************************************************************/
diff --git a/saves.h b/saves.h
index aa3bf7c..5d846ff 100644
--- a/saves.h
+++ b/saves.h
@@ -1,3 +1,25 @@
+/*
+*	Copyright (c) 2014 Michael Schellenberger Costa
+*
+*	Permission is hereby granted, free of charge, to any person obtaining a copy
+*	of this software and associated documentation files (the "Software"), to deal
+*	in the Software without restriction, including without limitation the rights
+*	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+*	copies of the Software, and to permit persons to whom the Software is
+*	furnished to do so, subject to the following conditions:
+*
+*	The above copyright notice and this permission notice shall be included in
+*	all copies or substantial portions of the Software.
+*
+*	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+*	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+*	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+*	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+*	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+*	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+*	THE SOFTWARE.
+*/
+
 /****************************************************************************************************/
 /*									Functions for data storage										*/
 /****************************************************************************************************/