diff --git a/Data_Storage.h b/Data_Storage.h
index 9e5d34b..25b14be 100644
--- a/Data_Storage.h
+++ b/Data_Storage.h
@@ -20,28 +20,17 @@
  *	THE SOFTWARE.
  *
  *	AUTHORS:	Michael Schellenberger Costa: mschellenbergercosta@gmail.com
- *
- *	Based on:	A thalamocortical neural mass model of the EEG during NREM sleep and its response
- *				to auditory stimulation.
- *				M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz,
- *				JC Claussen.
- *				PLoS Computational Biology (in review).
  */
 
-/****************************************************************************************************/
-/*									Functions for data storage										*/
-/****************************************************************************************************/
+/******************************************************************************/
+/*                        Functions for data storage                          */
+/******************************************************************************/
 #pragma once
+#include <vector>
 #include "Thalamic_Column.h"
 
-/****************************************************************************************************/
-/*											Save data												*/
-/****************************************************************************************************/
-void get_data(int counter, Thalamic_Column& Col,	double* Vt,	double* Vr,	double* ah) {
-	Vt 		[counter] = Col.Vt 		[0];
-	Vr 		[counter] = Col.Vr 		[0];
-	ah 		[counter] = Col.act_h 	();
+void get_data(unsigned counter, Thalamic_Column& Thalamus, std::vector<double*> pData) {
+    pData[0][counter] = Thalamus.Vt 		[0];
+    pData[0][counter] = Thalamus.Vr 		[0];
+    pData[0][counter] = Thalamus.act_h 	();
 }
-/****************************************************************************************************/
-/*										 		end													*/
-/****************************************************************************************************/
diff --git a/NM_Thalamus.pro b/NM_Thalamus.pro
index f6c53f8..0e7cabf 100644
--- a/NM_Thalamus.pro
+++ b/NM_Thalamus.pro
@@ -6,14 +6,17 @@ CONFIG -= qt
 TARGET = Thalamus.cpp
 
 SOURCES +=  Thalamic_Column.cpp \
-	    Thalamus.cpp	\
-	    Thalamus_mex.cpp
+			Thalamus.cpp	\
+			Thalamus_mex.cpp
 
 HEADERS +=  Data_Storage.h	\
-	    ODE.h		\
-	    Random_Stream.h	\
-	    Thalamic_Column.h
-
-QMAKE_CXXFLAGS += -std=c++11 -O3
+			ODE.h		\
+			Random_Stream.h	\
+			Thalamic_Column.h
 
 SOURCES -= Thalamus_mex.cpp
+
+QMAKE_CXXFLAGS += -std=c++11
+QMAKE_CXXFLAGS_RELEASE -= -O1
+QMAKE_CXXFLAGS_RELEASE -= -O2
+QMAKE_CXXFLAGS_RELEASE *= -O3
diff --git a/Random_Stream.h b/Random_Stream.h
index f5fb19a..29d013e 100644
--- a/Random_Stream.h
+++ b/Random_Stream.h
@@ -23,58 +23,35 @@
  *              Stefanie Gareis: gareis@inb.uni-luebeck.de
  */
 
-/****************************************************************************************************/
-/*                                       Random number streams                                      */
-/****************************************************************************************************/
+/******************************************************************************/
+/*                          Random number streams                             */
+/******************************************************************************/
 #pragma once
 #include <random>
 
-/****************************************************************************************************/
-/*									Struct for normal distribution                                  */
-/****************************************************************************************************/
-struct random_stream_normal
-{
-    /* Random number engine: Mersenne-Twister */
+class random_stream_normal {
+public:
+    explicit random_stream_normal(double mean, double stddev)
+    : mt(rand()), norm_dist(mean, stddev) {}
+    explicit random_stream_normal(double mean, double stddev, double seed)
+    : mt(seed), norm_dist(mean, stddev) {}
+
+    double operator ()(void) { return norm_dist(mt); }
+private:
     std::mt19937_64                     mt;
-    /* Random number generator: Normal-distribution */
     std::normal_distribution<double>    norm_dist;
-
-    /* Constructors */
-    random_stream_normal(){}
-    random_stream_normal(double mean, double stddev)
-    : mt(rand()) , norm_dist(mean, stddev)
-    {}
-
-    /* Overwrites the function-call operator "( )" */
-    double operator( )(void) {
-        return norm_dist(mt);
-    }
 };
-/****************************************************************************************************/
-/*										 		end													*/
-/****************************************************************************************************/
 
-/****************************************************************************************************/
-/*									Struct for uniform int distribution                             */
-/****************************************************************************************************/
-struct random_stream_uniform_int
-{
-    /* Random number engine: Mersenne-Twister */
+class random_stream_uniform_int {
+public:
+    explicit random_stream_uniform_int(int lower_bound, int upper_bound)
+    : mt(rand()), uniform_dist(lower_bound, upper_bound) {}
+    explicit random_stream_uniform_int(int lower_bound, int upper_bound, double seed)
+    : mt(seed), uniform_dist(lower_bound, upper_bound) {}
+
+    double operator ()(void) { return uniform_dist(mt); }
+private:
     std::mt19937_64                     mt;
-    /* Random number generator: Uniform integer-distribution */
     std::uniform_int_distribution<>     uniform_dist;
-
-    /* Constructors */
-    random_stream_uniform_int(){}
-    random_stream_uniform_int(double lower_bound, double upper_bound)
-    : mt(rand()) , uniform_dist(lower_bound, upper_bound)
-    {}
-
-    /* Overwrites the function-call operator "( )" */
-    double operator( )(void) {
-        return uniform_dist(mt);
-    }
 };
-/****************************************************************************************************/
-/*										 		end													*/
-/****************************************************************************************************/
+
diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp
index 29de11b..f9dd008 100644
--- a/Thalamic_Column.cpp
+++ b/Thalamic_Column.cpp
@@ -25,299 +25,228 @@
  *				to auditory stimulation.
  *				M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz,
  *				JC Claussen.
- *				PLoS Computational Biology (in review).
+ *				PLoS Computational Biology http://dx.doi.org/10.1371/journal.pcbi.1005022
  */
 
-/****************************************************************************************************/
-/*									Functions of the thalamic module								*/
-/****************************************************************************************************/
+/******************************************************************************/
+/*							Functions of the thalamic module				  */
+/******************************************************************************/
 #include "Thalamic_Column.h"
-/****************************************************************************************************/
-/*										 Initialization of RNG 										*/
-/****************************************************************************************************/
+
+/******************************************************************************/
+/*							Initialization of RNG 							  */
+/******************************************************************************/
 void Thalamic_Column::set_RNG(void) {
-	extern const double dt;
-	/* Number of independent random variables */
-	int N = 1;
-
-	/* Create RNG for each stream */
-	for (int i=0; i<N; ++i){
-		/* Add the RNG for I_{l}*/
-		MTRands.push_back(random_stream_normal(0.0, dphi*dt));
-
-		/* Add the RNG for I_{l,0}
-		 * In the case of a noise free system (dphi=0) set I_{l,0}
-		 * to zero.
-		 */
-		if(dphi==0)
-			MTRands.push_back(random_stream_normal(0.0, 0));
-		else
-			MTRands.push_back(random_stream_normal(0.0, dt));
-
-		/* Get the random number for the first iteration */
-		Rand_vars.push_back(MTRands[2*i]());
-		Rand_vars.push_back(MTRands[2*i+1]());
-	}
+    extern const double dt;
+    unsigned numRandomVariables = 1;
+
+    MTRands.reserve(2*numRandomVariables);
+    Rand_vars.reserve(2*numRandomVariables);
+    for (unsigned i=0; i < numRandomVariables; ++i){
+        /* Add the RNG for I_{l}*/
+        MTRands.push_back(random_stream_normal(0.0, dphi*dt));
+
+        /* Add the RNG for I_{l,0} */
+        MTRands.push_back(random_stream_normal(0.0, dt));
+
+        /* Get the random number for the first iteration */
+        Rand_vars.push_back(MTRands[2*i]());
+        Rand_vars.push_back(MTRands[2*i+1]());
+    }
+}
+/******************************************************************************/
+/*                          RK noise scaling 								  */
+/******************************************************************************/
+double Thalamic_Column::noise_xRK(int N, int M) const{
+    return gamma_e * gamma_e * (Rand_vars[2*M] + Rand_vars[2*M+1]/std::sqrt(3))*B[N];
 }
-/****************************************************************************************************/
-/*										 		end			 										*/
-/****************************************************************************************************/
 
+double Thalamic_Column::noise_aRK(int M) const{
+    return gamma_e * gamma_e * (Rand_vars[2*M] - Rand_vars[2*M+1]*std::sqrt(3))/4;
+}
 
-/****************************************************************************************************/
-/*										 Firing Rate functions 										*/
-/****************************************************************************************************/
-/* Thalamic relay (TC) firing rate */
+/******************************************************************************/
+/*                          Firing Rate functions 							  */
+/******************************************************************************/
 double Thalamic_Column::get_Qt	(int N) const{
-	double q = Qt_max/ (1 + exp(-C1 * (Vt[N] - theta_t) / sigma_t));
-	return q;
+    return Qt_max/ (1 + exp(-C1 * (Vt[N] - theta_t) / sigma_t));
 }
 
-/* Thalamic reticular (RE) firing rate */
 double Thalamic_Column::get_Qr	(int N) const{
-	double q = Qr_max / (1 + exp(-C1 * (Vr[N]- theta_r) / sigma_r));
-	return q;
+    return Qr_max / (1 + exp(-C1 * (Vr[N]- theta_r) / sigma_r));
 }
-/****************************************************************************************************/
-/*										 		end			 										*/
-/****************************************************************************************************/
-
 
-/****************************************************************************************************/
-/*										Synaptic currents											*/
-/****************************************************************************************************/
+/******************************************************************************/
+/*							Synaptic currents								  */
+/******************************************************************************/
 /* Excitatory input to TC population */
 double Thalamic_Column::I_et	(int N) const{
-	double I = g_AMPA * s_et[N]* (Vt[N]- E_AMPA);
-	return I;
+    return g_AMPA * s_et[N]* (Vt[N]- E_AMPA);
 }
 
 /* Inhibitory input to TC population */
 double Thalamic_Column::I_gt	(int N) const{
-	double I = g_GABA * s_gt[N]* (Vt[N]- E_GABA);
-	return I;
+    return g_GABA * s_gt[N]* (Vt[N]- E_GABA);
 }
 /* Excitatory input to RE population */
 double Thalamic_Column::I_er	(int N) const{
-	double I = g_AMPA * s_er[N]* (Vr[N]- E_AMPA);
-	return I;
+    return g_AMPA * s_er[N]* (Vr[N]- E_AMPA);
 }
 
 /* Inhibitory input to RE population */
 double Thalamic_Column::I_gr	(int N) const{
-	double I = g_GABA * s_gr[N]* (Vr[N]- E_GABA);
-	return I;
+    return g_GABA * s_gr[N]* (Vr[N]- E_GABA);
 }
-/****************************************************************************************************/
-/*										 		end			 										*/
-/****************************************************************************************************/
-
 
-/****************************************************************************************************/
-/*										 	I_T gating	 											*/
-/****************************************************************************************************/
+/******************************************************************************/
+/*                          I_T gating functions 							  */
+/******************************************************************************/
 /* Activation in TC population after Destexhe 1996 */
 double Thalamic_Column::m_inf_T_t	(int N) const{
-	double m = 1/(1+exp(-(Vt[N]+59)/6.2));
-	return m;
+    return 1/(1+exp(-(Vt[N]+59)/6.2));
 }
 
 /* Activation in RE population after Destexhe 1996 */
 double Thalamic_Column::m_inf_T_r	(int N) const{
-	double m = 1/(1+exp(-(Vr[N]+52)/7.4));
-	return m;
+    return 1/(1+exp(-(Vr[N]+52)/7.4));
 }
 
 /* Deactivation in TC population after Destexhe 1996 */
 double Thalamic_Column::h_inf_T_t	(int N) const{
-	double h = 1/(1+exp( (Vt[N]+81)/4));
-	return h;
+    return 1/(1+exp( (Vt[N]+81)/4));
 }
 
 /* Deactivation in RE population after Destexhe 1996 */
 double Thalamic_Column::h_inf_T_r	(int N) const{
-	double h = 1/(1+exp( (Vr[N]+80)/5));
-	return h;
+    return 1/(1+exp( (Vr[N]+80)/5));
 }
 
 /* Deactivation time in RE population after Destexhe 1996 */
 double Thalamic_Column::tau_h_T_t	(int N) const{
-	double tau =  (30.8 + (211.4 + exp((Vt[N]+115.2)/5))/(1 + exp((Vt[N]+86)/3.2)))/3.7371928;
-	return tau;
+    return (30.8 + (211.4 + exp((Vt[N]+115.2)/5))/(1 + exp((Vt[N]+86)/3.2)))/3.7371928;
 }
 
 /* Deactivation time in RE population after Destexhe 1996 */
 double Thalamic_Column::tau_h_T_r	(int N) const{
-	double tau =  (85 + 1/(exp((Vr[N]+48)/4) + exp(-(Vr[N]+407)/50)))/3.7371928;
-	return tau;
+    return (85 + 1/(exp((Vr[N]+48)/4) + exp(-(Vr[N]+407)/50)))/3.7371928;
 }
-/****************************************************************************************************/
-/*										 		end			 										*/
-/****************************************************************************************************/
-
 
-/****************************************************************************************************/
-/*										 	I_h gating 												*/
-/****************************************************************************************************/
+/******************************************************************************/
+/*                          I_h gating functions 							  */
+/******************************************************************************/
 /* Activation in TC population after Destexhe 1993 */
 double Thalamic_Column::m_inf_h	(int N) const{
-	double h = 1/(1+exp( (Vt[N]+75)/5.5));
-	return h;
+    return 1/(1+exp( (Vt[N]+75)/5.5));
 }
 
 /* Activation time for slow components in TC population after Chen2012 */
 double Thalamic_Column::tau_m_h	(int N) const{
-	double tau = (20 + 1000/(exp((Vt[N]+ 71.5)/14.2) + exp(-(Vt[N]+ 89)/11.6)));
-	return tau;
+    return (20 + 1000/(exp((Vt[N]+ 71.5)/14.2) + exp(-(Vt[N]+ 89)/11.6)));
 }
 
 /* Instantaneous calcium binding onto messenger protein after Chen2012 */
 double Thalamic_Column::P_h	(int N) const{
-	//double P_h = k1 * pow(Ca[N], n_P)/(k1*pow(Ca[N], n_P)+k2);
-	double P_h = k1 * Ca[N] * Ca[N] * Ca[N] * Ca[N]/(k1 * Ca[N] * Ca[N] * Ca[N] * Ca[N]+k2);
-	return P_h;
+    //return k1 * pow(Ca[N], n_P)/(k1*pow(Ca[N], n_P)+k2);
+    return k1 * Ca[N] * Ca[N] * Ca[N] * Ca[N]/(k1 * Ca[N] * Ca[N] * Ca[N] * Ca[N]+k2);
 }
 
 /* Return I_h activation */
 double Thalamic_Column::act_h	(void) const{
-	double activation = m_h[0] + g_inc * m_h2[0];
-	return activation;
+    return m_h[0] + g_inc * m_h2[0];
 }
-/****************************************************************************************************/
-/*										 		end			 										*/
-/****************************************************************************************************/
-
 
-/****************************************************************************************************/
-/*										 Current functions 											*/
-/****************************************************************************************************/
+/******************************************************************************/
+/*							Intrinsic currents                                */
+/******************************************************************************/
 /* Leak current of TC population */
 double Thalamic_Column::I_L_t	(int N) const{
-	double I = g_L * (Vt[N]- E_L_t);
-	return I;
+    return g_L * (Vt[N]- E_L_t);
+
 }
 
 /* Potassium leak current of TC population */
 double Thalamic_Column::I_LK_t	(int N) const{
-	double I = g_LK * (Vt[N]- E_K);
-	return I;
+    return g_LK * (Vt[N]- E_K);
 }
 
 /* Leak current of RE population */
 double Thalamic_Column::I_L_r	(int N) const{
-	double I = g_L * (Vr[N]- E_L_r);
-	return I;
+    return g_L * (Vr[N]- E_L_r);
 }
 
 /* Potassium leak current of RE population */
 double Thalamic_Column::I_LK_r	(int N) const{
-	double I = g_LK * (Vr[N]- E_K);
-	return I;
+    return g_LK	* (Vr[N]- E_K);
 }
 
 /* T-type current of TC population */
 double Thalamic_Column::I_T_t	(int N) const{
-	double I = g_T_t * m_inf_T_t(N) * m_inf_T_t(N) * h_T_t[N] * (Vt[N]- E_Ca);
-	return I;
+    return g_T_t * m_inf_T_t(N) * m_inf_T_t(N) * h_T_t[N] * (Vt[N]- E_Ca);
 }
 
 /* T-type current of RE population */
 double Thalamic_Column::I_T_r	(int N) const{
-	double I = g_T_r * m_inf_T_r(N) * m_inf_T_r(N) * h_T_r[N] * (Vr[N]- E_Ca);
-	return I;
+    return g_T_r * m_inf_T_r(N) * m_inf_T_r(N) * h_T_r[N] * (Vr[N]- E_Ca);
 }
 
 /* h-type current of TC population */
 double Thalamic_Column::I_h		(int N) const{
-	double I = g_h * (m_h[N] + g_inc * m_h2[N]) * (Vt[N]- E_h);
-	return I;
+    return g_h * (m_h[N] + g_inc * m_h2[N]) * (Vt[N]- E_h);
 }
-/****************************************************************************************************/
-/*										 		end			 										*/
-/****************************************************************************************************/
-
 
-/****************************************************************************************************/
-/*										 RK noise scaling 											*/
-/****************************************************************************************************/
-double Thalamic_Column::noise_xRK(int N, int M) const{
-	return gamma_e * gamma_e * (Rand_vars[2*M] + Rand_vars[2*M+1]/std::sqrt(3))*B[N];
-}
-
-double Thalamic_Column::noise_aRK(int M) const{
-	return gamma_e * gamma_e * (Rand_vars[2*M] - Rand_vars[2*M+1]*std::sqrt(3))/4;
-}
-/****************************************************************************************************/
-/*										 		end			 										*/
-/****************************************************************************************************/
-
-
-/****************************************************************************************************/
-/*										Calculate the Nth SRK term									*/
-/****************************************************************************************************/
+/******************************************************************************/
+/*                              SRK iteration                                 */
+/******************************************************************************/
 void Thalamic_Column::set_RK (int N) {
-	extern const double dt;
-	Vt	  	[N+1] = Vt   [0] + A[N]*dt*(-(I_L_t(N) + I_et(N) + I_gt(N))/tau_t - C_m * (I_LK_t(N) + I_T_t(N) + I_h(N)));
-	Vr	  	[N+1] = Vr   [0] + A[N]*dt*(-(I_L_r(N) + I_er(N) + I_gr(N))/tau_r - C_m * (I_LK_r(N) + I_T_r(N)));
-	Ca      [N+1] = Ca   [0] + A[N]*dt*(alpha_Ca * I_T_t(N) - (Ca[N] - Ca_0)/tau_Ca);
-	h_T_t   [N+1] = h_T_t[0] + A[N]*dt*(h_inf_T_t(N) - h_T_t[N])/tau_h_T_t(N);
-	h_T_r 	[N+1] = h_T_r[0] + A[N]*dt*(h_inf_T_r(N) - h_T_r[N])/tau_h_T_r(N);
-	m_h 	[N+1] = m_h  [0] + A[N]*dt*((m_inf_h(N) * (1 - m_h2[N]) - m_h[N])/tau_m_h(N) - k3 * P_h(N) * m_h[N] + k4 * m_h2[N]);
-	m_h2 	[N+1] = m_h2 [0] + A[N]*dt*(k3 * P_h(N) * m_h[N] - k4 * m_h2[N]);
-	s_et	[N+1] = s_et [0] + A[N]*dt*(x_et[N]);
-	s_er	[N+1] = s_er [0] + A[N]*dt*(x_er[N]);
-	s_gt	[N+1] = s_gt [0] + A[N]*dt*(x_gt[N]);
-	s_gr	[N+1] = s_gr [0] + A[N]*dt*(x_gr[N]);
-	x_et  	[N+1] = x_et [0] + A[N]*dt*(pow(gamma_e, 2) * (                 - s_et[N]) - 2 * gamma_e * x_et[N]) + noise_xRK(N,0);
-	x_er  	[N+1] = x_er [0] + A[N]*dt*(pow(gamma_e, 2) * (N_tr * get_Qt(N)	- s_er[N]) - 2 * gamma_e * x_er[N]);
-	x_gt  	[N+1] = x_gt [0] + A[N]*dt*(pow(gamma_g, 2) * (N_rt * get_Qr(N) - s_gt[N]) - 2 * gamma_g * x_gt[N]);
-	x_gr  	[N+1] = x_gr [0] + A[N]*dt*(pow(gamma_g, 2) * (N_rr * get_Qr(N)	- s_gr[N]) - 2 * gamma_g * x_gr[N]);
+    extern const double dt;
+    Vt	  	[N+1] = Vt   [0] + A[N]*dt*(-(I_L_t(N) + I_et(N) + I_gt(N))/tau_t - C_m * (I_LK_t(N) + I_T_t(N) + I_h(N)));
+    Vr	  	[N+1] = Vr   [0] + A[N]*dt*(-(I_L_r(N) + I_er(N) + I_gr(N))/tau_r - C_m * (I_LK_r(N) + I_T_r(N)));
+    Ca      [N+1] = Ca   [0] + A[N]*dt*(alpha_Ca * I_T_t(N) - (Ca[N] - Ca_0)/tau_Ca);
+    h_T_t   [N+1] = h_T_t[0] + A[N]*dt*(h_inf_T_t(N) - h_T_t[N])/tau_h_T_t(N);
+    h_T_r 	[N+1] = h_T_r[0] + A[N]*dt*(h_inf_T_r(N) - h_T_r[N])/tau_h_T_r(N);
+    m_h 	[N+1] = m_h  [0] + A[N]*dt*((m_inf_h(N) * (1 - m_h2[N]) - m_h[N])/tau_m_h(N) - k3 * P_h(N) * m_h[N] + k4 * m_h2[N]);
+    m_h2 	[N+1] = m_h2 [0] + A[N]*dt*(k3 * P_h(N) * m_h[N] - k4 * m_h2[N]);
+    s_et	[N+1] = s_et [0] + A[N]*dt*(x_et[N]);
+    s_er	[N+1] = s_er [0] + A[N]*dt*(x_er[N]);
+    s_gt	[N+1] = s_gt [0] + A[N]*dt*(x_gt[N]);
+    s_gr	[N+1] = s_gr [0] + A[N]*dt*(x_gr[N]);
+    x_et  	[N+1] = x_et [0] + A[N]*dt*(gamma_e*gamma_e * (                 - s_et[N]) - 2 * gamma_e * x_et[N]) + noise_xRK(N,0);
+    x_er  	[N+1] = x_er [0] + A[N]*dt*(gamma_e*gamma_e * (N_rt * get_Qt(N)	- s_er[N]) - 2 * gamma_e * x_er[N]);
+    x_gt  	[N+1] = x_gt [0] + A[N]*dt*(gamma_g*gamma_g * (N_tr * get_Qr(N)	- s_gt[N]) - 2 * gamma_g * x_gt[N]);
+    x_gr  	[N+1] = x_gr [0] + A[N]*dt*(gamma_g*gamma_g * (N_rr * get_Qr(N)	- s_gr[N]) - 2 * gamma_g * x_gr[N]);
 }
 /****************************************************************************************************/
 /*										 		end			 										*/
 /****************************************************************************************************/
 
 
-/****************************************************************************************************/
-/*									Function that adds all SRK terms								*/
-/****************************************************************************************************/
 void Thalamic_Column::add_RK(void) {
-	Vt	  	[0] =(-3*Vt   [0] + 2*Vt   [1] + 4*Vt	[2] + 2* Vt     [3] + Vt	[4])/6;
-	Vr	  	[0] =(-3*Vr   [0] + 2*Vr   [1] + 4*Vr	[2] + 2* Vr     [3] + Vr	[4])/6;
-	Ca	  	[0] =(-3*Ca   [0] + 2*Ca   [1] + 4*Ca	[2] + 2* Ca     [3] + Ca	[4])/6;
-	s_et	[0] =(-3*s_et [0] + 2*s_et [1] + 4*s_et	[2] + 2* s_et	[3] + s_et	[4])/6;
-	s_er	[0] =(-3*s_er [0] + 2*s_er [1] + 4*s_er	[2] + 2* s_er	[3] + s_er	[4])/6;
-	s_gt	[0] =(-3*s_gt [0] + 2*s_gt [1] + 4*s_gt	[2] + 2* s_gt	[3] + s_gt	[4])/6;
-	s_gr	[0] =(-3*s_gr [0] + 2*s_gr [1] + 4*s_gr	[2] + 2* s_gr	[3] + s_gr	[4])/6;
-	x_et  	[0] =(-3*x_et [0] + 2*x_et [1] + 4*x_et	[2] + 2* x_et	[3] + x_et	[4])/6 + noise_aRK(0);
-	x_er  	[0] =(-3*x_er [0] + 2*x_er [1] + 4*x_er	[2] + 2* x_er	[3] + x_er	[4])/6;
-	x_gt  	[0] =(-3*x_gt [0] + 2*x_gt [1] + 4*x_gt	[2] + 2* x_gt	[3] + x_gt	[4])/6;
-	x_gr  	[0] =(-3*x_gr [0] + 2*x_gr [1] + 4*x_gr	[2] + 2* x_gr	[3] + x_gr	[4])/6;
-	h_T_t 	[0] =(-3*h_T_t[0] + 2*h_T_t[1] + 4*h_T_t[2] + 2* h_T_t  [3] + h_T_t	[4])/6;
-	h_T_r 	[0] =(-3*h_T_r[0] + 2*h_T_r[1] + 4*h_T_r[2] + 2* h_T_r  [3] + h_T_r	[4])/6;
-	m_h		[0] =(-3*m_h  [0] + 2*m_h  [1] + 4*m_h	[2] + 2* m_h    [3] + m_h	[4])/6;
-	m_h2	[0] =(-3*m_h2 [0] + 2*m_h2 [1] + 4*m_h2	[2] + 2* m_h2	[3] + m_h2	[4])/6;
-
-	/* Generate noise for the next iteration */
-	for (unsigned i=0; i<Rand_vars.size(); ++i) {
-		Rand_vars[i] = MTRands[i]() + input;
-	}
+    add_RK(Vt);
+    add_RK(Vr);
+    add_RK(Ca);
+    add_RK(s_et);
+    add_RK(s_er);
+    add_RK(s_gt);
+    add_RK(s_gr);
+    add_RK_noise(x_et, 0);
+    add_RK(x_er);
+    add_RK(x_gt);
+    add_RK(x_gr);
+    add_RK(h_T_t);
+    add_RK(h_T_r);
+    add_RK(m_h);
+    add_RK(m_h2);
+
+    /* Generate noise for the next iteration */
+    for (unsigned i=0; i<Rand_vars.size(); ++i) {
+        Rand_vars[i] = MTRands[i]() + input;
+    }
 }
-/****************************************************************************************************/
-/*										 		end			 										*/
-/****************************************************************************************************/
-
 
-/****************************************************************************************************/
-/*										Evaluation of SRK4											*/
-/****************************************************************************************************/
 void Thalamic_Column::iterate_ODE(void) {
-	/* First calculating every ith RK moment. Has to be in order, 1th moment first */
-	for (int i=0; i<4; ++i) {
-		set_RK(i);
-	}
-	add_RK();
+    /* First calculating every ith RK moment. Has to be in order, 1th moment first */
+    for (unsigned i=0; i<4; ++i) {
+        set_RK(i);
+    }
+    add_RK();
 }
-/****************************************************************************************************/
-/*										 		end													*/
-/****************************************************************************************************/
diff --git a/Thalamic_Column.h b/Thalamic_Column.h
index 498b2d9..364711c 100644
--- a/Thalamic_Column.h
+++ b/Thalamic_Column.h
@@ -25,207 +25,196 @@
  *				to auditory stimulation.
  *				M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz,
  *				JC Claussen.
- *				PLoS Computational Biology (in review).
+ *				PLoS Computational Biology http://dx.doi.org/10.1371/journal.pcbi.1005022
  */
 
-/************************************************************************************************/
-/*								Header file of a thalamic module								*/
-/************************************************************************************************/
+/******************************************************************************/
+/*						Implementation of a thalamic module					  */
+/******************************************************************************/
 #pragma once
 #include <cmath>
 #include <vector>
-#include "Random_Stream.h"
-using std::vector;
-
-/****************************************************************************************************/
-/*									Macro for vector initialization									*/
-/****************************************************************************************************/
-#ifndef _INIT
-#define _INIT(x)	{x, 0.0, 0.0, 0.0, 0.0}
-#endif
-/****************************************************************************************************/
-/*										 		end			 										*/
-/****************************************************************************************************/
 
+#include "Random_Stream.h"
 
-/****************************************************************************************************/
-/*								Implementation of the thalamic module 								*/
-/****************************************************************************************************/
 class Thalamic_Column {
 public:
-	/* Constructors for compile check */
-	Thalamic_Column(void)
-	{set_RNG();}
-
-	/* Constructor for simulation */
-	Thalamic_Column(double* Par)
-	: 	g_LK		(Par[0]),	g_h 		(Par[1])
-	{set_RNG();}
-
-	/* Iterate one time step through SRK4 */
-	void	iterate_ODE (void);
-
-	/* Data storage  access */
-	friend void get_data (int, Thalamic_Column&, double*, double*, double*);
+    Thalamic_Column(double* Par)
+    : g_LK (Par[0]), g_h (Par[1]) {set_RNG();}
 
+    /* Iterate one time step through SRK4 */
+    void	iterate_ODE (void);
 private:
-	/* Declaration of private functions */
-	/* Initialize the RNGs */
-	void 	set_RNG		(void);
-
-	/* Firing rates */
-	double 	get_Qt		(int) const;
-	double 	get_Qr		(int) const;
-
-	/* Synaptic currents */
-	double 	I_et		(int) const;
-	double 	I_gt		(int) const;
-	double 	I_er		(int) const;
-	double 	I_gr		(int) const;
-
-	/* Activation functions */
-	double  m_inf_T_t	(int) const;
-	double  m_inf_T_r	(int) const;
-	double  m_inf_h		(int) const;
-	double  tau_m_h		(int) const;
-	double  P_h			(int) const;
-	double  act_h		(void)const;
-
-	/* 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 */
-	double 	I_L_t		(int) const;
-	double 	I_L_r		(int) const;
-	double 	I_LK_t		(int) const;
-	double 	I_LK_r		(int) const;
-	double 	I_T_t		(int) const;
-	double 	I_T_r		(int) const;
-	double 	I_h			(int) const;
-
-	/* Noise functions */
-	double 	noise_xRK 	(int,int) const;
-	double 	noise_aRK 	(int) const;
-
-	/* ODE functions */
-	void 	set_RK		(int);
-	void 	add_RK	 	(void);
-
-	/* Declaration and Initialization of parameters */
-	/* Membrane time in ms */
-	const double 	tau_t 		= 20.;
-	const double 	tau_r 		= 20.;
-
-	/* Maximum firing rate in ms^-1 */
-	const double 	Qt_max		= 400.E-3;
-	const double 	Qr_max		= 400.E-3;
-
-	/* Sigmoid threshold in mV */
-	const double 	theta_t		= -58.5;
-	const double 	theta_r		= -58.5;
-
-	/* Sigmoid gain in mV */
-	const double 	sigma_t		= 6.;
-	const double 	sigma_r		= 6.;
-
-	/* Scaling parameter for sigmoidal mapping (dimensionless) */
-	const double 	C1          = (M_PI/sqrt(3));
-
-	/* PSP rise time in ms^-1 */
-	const double 	gamma_e		= 70E-3;
-	const double 	gamma_g		= 100E-3;
-
-	/* Membrane capacitance in muF/cm^2 */
-	const double	C_m			= 1.;
-
-	/* Weights/ conductivities */
-	/* Leak  in aU */
-	const double 	g_L    		= 1.;
-
-	/* Synaptic conductivity in ms */
-	const double 	g_AMPA 		= 1.;
-	const double 	g_GABA 		= 1.;
-
-	/* Potassium leak current in mS/m^2 */
-	double			g_LK 		= 0.033;
-
-	/* T current in mS/m^2 */
-	const double	g_T_t		= 3;
-	const double	g_T_r		= 2.3;
-
-	/* h current in mS/m^2 */
-	double			g_h			= 0.02;
-
-	/* 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 		= -70.;
-
-	/* Potassium */
-	const double 	E_K    		= -100.;
-
-	/* I_T current */
-	const double 	E_Ca    	= 120.;
-
-	/* I_h current */
-	const double 	E_h    		= -40.;
-
-	/* Calcium parameters */
-	const double	alpha_Ca	= -51.8E-6;			/* influx per spike in nmol		*/
-	const double	tau_Ca		= 10.;				/* calcium time constant in ms	*/
-	const double	Ca_0		= 2.4E-4;			/* resting concentration 		*/
-
-	/* I_h activation parameters */
-	const double 	k1			= 2.5E7;
-	const double 	k2			= 4E-4;
-	const double 	k3			= 1E-1;
-	const double 	k4			= 1E-3;
-	const double 	n_P			= 4.;
-	const double 	g_inc		= 2.;
-
-	/* Noise parameters in ms^-1 */
-	const double 	mphi		= 0E-3;
-	const double	dphi		= 0E-3;
-	double			input		= 0.0;
-
-	/* Connectivities (dimensionless) */
-	const double 	N_tr		= 3.;
-	const double 	N_rt		= 5.;
-	const double 	N_rr		= 25.;
-
-	/* Parameters for SRK4 iteration */
-	const vector<double> A = {0.5,  0.5,  1.0, 1.0};
-	const vector<double> B = {0.75, 0.75, 0.0, 0.0};
-
-	/* Random number generators */
-	vector<random_stream_normal> MTRands;
-
-	/* Container for noise */
-	vector<double>	Rand_vars;
-
-	/* 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			*/
-					s_et	= _INIT(0.0),		/* PostSP from TC population to TC population		*/
-					s_er	= _INIT(0.0),		/* PostSP from TC population to RE population		*/
-					s_gt	= _INIT(0.0),		/* PostSP from RE population to TC population		*/
-					s_gr	= _INIT(0.0),		/* PostSP from RE population to RE population		*/
-					x_et	= _INIT(0.0),		/* derivative of s_et								*/
-					x_er	= _INIT(0.0),		/* derivative of s_er								*/
-					x_gt	= _INIT(0.0),		/* derivative of s_gt								*/
-					x_gr	= _INIT(0.0),		/* derivative of s_gr								*/
-					h_T_t	= _INIT(0.0),		/* inactivation of T channel						*/
-					h_T_r	= _INIT(0.0),		/* inactivation of T channel						*/
-					m_h		= _INIT(0.0),		/* activation 	of h   channel						*/
-					m_h2	= _INIT(0.0);		/* activation 	of h   channel bound with protein 	*/
+    /* Declaration of private functions */
+    /* Initialize the RNGs */
+    void 	set_RNG		(void);
+
+    /* Firing rates */
+    double 	get_Qt		(int) const;
+    double 	get_Qr		(int) const;
+
+    /* Synaptic currents */
+    double 	I_et		(int) const;
+    double 	I_gt		(int) const;
+    double 	I_er		(int) const;
+    double 	I_gr		(int) const;
+
+    /* Activation functions */
+    double  m_inf_T_t	(int) const;
+    double  m_inf_T_r	(int) const;
+    double  m_inf_h		(int) const;
+    double  tau_m_h		(int) const;
+    double  P_h			(int) const;
+    double  act_h		(void)const;
+
+    /* 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 */
+    double 	I_L_t		(int) const;
+    double 	I_L_r		(int) const;
+    double 	I_LK_t		(int) const;
+    double 	I_LK_r		(int) const;
+    double 	I_T_t		(int) const;
+    double 	I_T_r		(int) const;
+    double 	I_h			(int) const;
+
+    /* Noise functions */
+    double 	noise_xRK 	(int,int) const;
+    double 	noise_aRK 	(int) const;
+
+    /* ODE functions */
+    void 	set_RK		(int);
+    void 	add_RK	 	(void);
+
+    /* Helper functions */
+    inline std::vector<double> init (double value)
+    {return {value, 0.0, 0.0, 0.0, 0.0};}
+
+    inline void add_RK (std::vector<double>& var)
+    {var[0] = (-3*var[0] + 2*var[1] + 4*var[2] + 2*var[3] + var[4])/6;}
+
+    inline void add_RK_noise (std::vector<double>& var, unsigned noise)
+    {var[0] = (-3*var[0] + 2*var[1] + 4*var[2] + 2*var[3] + var[4])/6 + noise_aRK(noise);}
+
+    /* Declaration and Initialization of parameters */
+    /* Membrane time in ms */
+    const double 	tau_t 		= 20.;
+    const double 	tau_r 		= 20.;
+
+    /* Maximum firing rate in ms^-1 */
+    const double 	Qt_max		= 400.E-3;
+    const double 	Qr_max		= 400.E-3;
+
+    /* Sigmoid threshold in mV */
+    const double 	theta_t		= -58.5;
+    const double 	theta_r		= -58.5;
+
+    /* Sigmoid gain in mV */
+    const double 	sigma_t		= 6.;
+    const double 	sigma_r		= 6.;
+
+    /* Scaling parameter for sigmoidal mapping (dimensionless) */
+    const double 	C1          = (M_PI/sqrt(3));
+
+    /* PSP rise time in ms^-1 */
+    const double 	gamma_e		= 70E-3;
+    const double 	gamma_g		= 100E-3;
+
+    /* Membrane capacitance in muF/cm^2 */
+    const double	C_m			= 1.;
+
+    /* Weights/ conductivities */
+    /* Leak  in aU */
+    const double 	g_L    		= 1.;
+
+    /* Synaptic conductivity in ms */
+    const double 	g_AMPA 		= 1.;
+    const double 	g_GABA 		= 1.;
+
+    /* Potassium leak current in mS/m^2 */
+    double			g_LK 		= 0.033;
+
+    /* T current in mS/m^2 */
+    const double	g_T_t		= 3;
+    const double	g_T_r		= 2.3;
+
+    /* h current in mS/m^2 */
+    double			g_h			= 0.02;
+
+    /* 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 		= -70.;
+
+    /* Potassium */
+    const double 	E_K    		= -100.;
+
+    /* I_T current */
+    const double 	E_Ca    	= 120.;
+
+    /* I_h current */
+    const double 	E_h    		= -40.;
+
+    /* Calcium parameters */
+    const double	alpha_Ca	= -51.8E-6;			/* influx per spike in nmol		*/
+    const double	tau_Ca		= 10.;				/* calcium time constant in ms	*/
+    const double	Ca_0		= 2.4E-4;			/* resting concentration 		*/
+
+    /* I_h activation parameters */
+    const double 	k1			= 2.5E7;
+    const double 	k2			= 4E-4;
+    const double 	k3			= 1E-1;
+    const double 	k4			= 1E-3;
+    const double 	n_P			= 4.;
+    const double 	g_inc		= 2.;
+
+    /* Noise parameters in ms^-1 */
+    const double 	mphi		= 0E-3;
+    const double	dphi		= 0E-3;
+    double			input		= 0.0;
+
+    /* Connectivities (dimensionless) */
+    const double 	N_tr		= 3.;
+    const double 	N_rt		= 5.;
+    const double 	N_rr		= 25.;
+
+    /* Parameters for SRK4 iteration */
+    const std::vector<double> A = {0.5,  0.5,  1.0, 1.0};
+    const std::vector<double> B = {0.75, 0.75, 0.0, 0.0};
+
+    /* Random number generators */
+    std::vector<random_stream_normal> MTRands;
+
+    /* Container for noise */
+    std::vector<double>	Rand_vars;
+
+    /* Population variables																			*/
+    std::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			*/
+                        s_et	= init(0.0),	/* PostSP from TC population to TC population		*/
+                        s_er	= init(0.0),	/* PostSP from TC population to RE population		*/
+                        s_gt	= init(0.0),	/* PostSP from RE population to TC population		*/
+                        s_gr	= init(0.0),	/* PostSP from RE population to RE population		*/
+                        x_et	= init(0.0),	/* derivative of s_et								*/
+                        x_er	= init(0.0),	/* derivative of s_er								*/
+                        x_gt	= init(0.0),	/* derivative of s_gt								*/
+                        x_gr	= init(0.0),	/* derivative of s_gr								*/
+                        h_T_t	= init(0.0),	/* inactivation of T channel						*/
+                        h_T_r	= init(0.0),	/* inactivation of T channel						*/
+                        m_h		= init(0.0),	/* activation 	of h   channel						*/
+                        m_h2	= init(0.0);	/* activation 	of h   channel bound with protein 	*/
+
+    /* Data storage  access */
+    friend void get_data (unsigned, Thalamic_Column&, std::vector<double*>);
 };
 /****************************************************************************************************/
 /*										 		end			 										*/
diff --git a/Thalamus.cpp b/Thalamus.cpp
index edec717..2ee9e3a 100644
--- a/Thalamus.cpp
+++ b/Thalamus.cpp
@@ -25,52 +25,46 @@
  *				to auditory stimulation.
  *				M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz,
  *				JC Claussen.
- *				PLoS Computational Biology (in review).
+ *				PLoS Computational Biology http://dx.doi.org/10.1371/journal.pcbi.1005022
  */
 
-/****************************************************************************************************/
-/*		Main file for compilation tests																*/
-/****************************************************************************************************/
+/******************************************************************************/
+/*                  Main file for compilation and runtime tests				  */
+/******************************************************************************/
 #include <iostream>
 #include <chrono>
+
 #include "Thalamic_Column.h"
 
-/****************************************************************************************************/
-/*										Fixed simulation settings									*/
-/****************************************************************************************************/
+/******************************************************************************/
+/*                          Fixed simulation settings						  */
+/******************************************************************************/
 typedef std::chrono::high_resolution_clock::time_point timer;
-extern const int T		= 30;								/* Simulation length s					*/
-extern const int res 	= 1E4;								/* number of iteration steps per s		*/
-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			 										*/
-/****************************************************************************************************/
-
+extern const int T      = 30;		/* Time until data is stored in  s		  */
+extern const int res 	= 1E4;		/* Number of iteration steps per s		  */
+extern const double dt 	= 1E3/res;	/* Duration of a time step in ms		  */
+extern const double h	= sqrt(dt); /* Square root of dt for SRK iteration	  */
 
-/****************************************************************************************************/
-/*										Main simulation routine										*/
-/****************************************************************************************************/
+/******************************************************************************/
+/*                              Main simulation routine						  */
+/******************************************************************************/
 int main(void) {
-	/* Initialize the populations */
-	Thalamic_Column Thalamus = Thalamic_Column();
+    std::vector<double> param = {0.2, 0.06};
+    Thalamic_Column Thalamus (param.data());
 
-	/* Take the time of the simulation */
-	timer start,end;
+    /* Take the time of the simulation */
+    timer start,end;
 
-	/* Simulation */
-	start = std::chrono::high_resolution_clock::now();
-	for (int t=0; t< T*res; ++t) {
-		Thalamus.iterate_ODE();
-	}
-	end = std::chrono::high_resolution_clock::now();
+    /* Simulation */
+    start = std::chrono::high_resolution_clock::now();
+    for (unsigned t=0; t< T*res; ++t) {
+        Thalamus.iterate_ODE();
+    }
+    end = std::chrono::high_resolution_clock::now();
 
-	/* Time consumed by the simulation */
-	double dif = 1E-3*std::chrono::duration_cast<std::chrono::milliseconds>( end - start ).count();
-	std::cout << "simulation done!\n";
-	std::cout << "took " << dif 	<< " seconds" << "\n";
-	std::cout << "end\n";
+    /* Time consumed by the simulation */
+    double dif = 1E-3*std::chrono::duration_cast<std::chrono::milliseconds>( end - start ).count();
+    std::cout << "simulation done!\n";
+    std::cout << "took " << dif 	<< " seconds" << "\n";
+    std::cout << "end\n";
 }
-/****************************************************************************************************/
-/*										 		end			 										*/
-/****************************************************************************************************/
diff --git a/Thalamus_mex.cpp b/Thalamus_mex.cpp
index ad1ef58..218aa9c 100644
--- a/Thalamus_mex.cpp
+++ b/Thalamus_mex.cpp
@@ -28,86 +28,87 @@
  *				PLoS Computational Biology (in review).
  */
 
-/****************************************************************************************************/
-/* 		Implementation of the simulation as MATLAB routine (mex compiler)							*/
-/* 		mex command is given by:																	*/
-/* 		mex CXXFLAGS="\$CXXFLAGS -std=c++11 -O3" Thalamus_mex.cpp Thalamic_Column.cpp				*/
-/****************************************************************************************************/
+/******************************************************************************/
+/* Implementation of the simulation as MATLAB routine (mex compiler)		  */
+/* mex command is given by:													  */
+/* mex CXXFLAGS="\$CXXFLAGS -std=c++11 -O3" Thalamus_mex.cpp                  */
+/*                                          Thalamic_Column.cpp               */
+/******************************************************************************/
 #include "mex.h"
 #include "matrix.h"
+
+#include <iterator>
+#include <vector>
+
 #include "Data_Storage.h"
-mxArray* SetMexArray(int N, int M);
+#include "Thalamic_Column.h"
+mxArray* GetMexArray(int N, int M);
 
-/****************************************************************************************************/
-/*										Fixed simulation settings									*/
-/****************************************************************************************************/
-extern const int onset	= 15;								/* time until data is stored in  s		*/
-extern const int res 	= 1E3;								/* 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			 										*/
-/****************************************************************************************************/
+/******************************************************************************/
+/*                          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 	= 1E2;		/* Number of iterations steps not saved	  */
+extern const double dt 	= 1E3/res;	/* Duration of a time step in ms		  */
+extern const double h	= sqrt(dt); /* Square root of dt for SRK iteration	  */
 
-/****************************************************************************************************/
-/*										Simulation routine	 										*/
-/*										lhs defines outputs											*/
-/*										rhs defines inputs											*/
-/****************************************************************************************************/
+/******************************************************************************/
+/*                              Simulation routine	 						  */
+/*								lhs defines outputs							  */
+/*								rhs defines inputs							  */
+/******************************************************************************/
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
-	/* Initialize the seeder */
-	srand(time(NULL));
+    /* 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 cortical module		*/
+    /* 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		*/
 
-	/* Initialize the populations */
-	Thalamic_Column Thalamus(Param_Thalamus);
+    /* Initialize the populations */
+    Thalamic_Column Thalamus(Param_Thalamus);
 
-	/* Create data containers */
-	mxArray* Vt		= SetMexArray(1, T*res/red);
-	mxArray* Vr		= SetMexArray(1, T*res/red);
-	mxArray* ah		= SetMexArray(1, T*res/red);
+    /* Create data containers */
+    std::vector<mxArray*> dataArray;
+    dataArray.reserve(3);
+    dataArray.push_back(GetMexArray(1, T*res/red));	// Vt
+    dataArray.push_back(GetMexArray(1, T*res/red));	// Vr
+    dataArray.push_back(GetMexArray(1, T*res/red));	// act_h
 
-	/* Pointer to the actual data block */
-	double* Pr_Vt	= mxGetPr(Vt);
-	double* Pr_Vr	= mxGetPr(Vr);
-	double* Pr_ah	= mxGetPr(ah);
+    /* Pointer to the data blocks */
+    std::vector<double*> dataPointer;
+    dataPointer.reserve(dataArray.size());
+    for (auto &dataptr : dataArray) {
+        dataPointer.push_back(mxGetPr(dataptr));
+    }
 
-	/* Simulation */
-	int count = 0;
-	for (int t=0; t<Time; ++t) {
-		Thalamus.iterate_ODE();
-		if(t>=onset*res && t%red==0){
-			get_data(count, Thalamus, Pr_Vt, Pr_Vr, Pr_ah);
-			++count;
-		}
-	}
+    /* Simulation */
+    int count = 0;
+    for (unsigned t=0; t < Time; ++t) {
+        Thalamus.iterate_ODE();
+        if(t >= onset*res && t%red == 0){
+            get_data(count, Thalamus, dataPointer);
+            ++count;
+        }
+    }
 
-	/* Output of the simulation */
-	plhs[0] = Vt;
-	plhs[1] = Vr;
-	plhs[2] = ah;
-return;
+    /* Return the data containers */
+    nlhs = dataArray.size()+1;
+    for (auto &dataptr : dataArray) {
+        plhs[std::distance(&dataptr, dataArray.data())] = dataptr;
+    }
+    return;
 }
-/****************************************************************************************************/
-/*												end													*/
-/****************************************************************************************************/
-
 
-/****************************************************************************************************/
-/*									Create MATLAB data container									*/
-/****************************************************************************************************/
-mxArray* SetMexArray(int N, int M) {
-	mxArray* Array	= mxCreateDoubleMatrix(0, 0, mxREAL);
-	mxSetM(Array, N);
-	mxSetN(Array, M);
-	mxSetData(Array, mxMalloc(sizeof(double)*M*N));
-	return Array;
+/******************************************************************************/
+/*                          Create MATLAB data containers					  */
+/******************************************************************************/
+mxArray* GetMexArray(int N, int M) {
+    mxArray* Array	= mxCreateDoubleMatrix(0, 0, mxREAL);
+    mxSetM(Array, N);
+    mxSetN(Array, M);
+    mxSetData(Array, mxMalloc(sizeof(double)*M*N));
+    return Array;
 }
-/****************************************************************************************************/
-/*										 		end													*/
-/****************************************************************************************************/