Updated the TC model to latest code changes.

Parameters are now locally defined as const double within class
definition.

Therefore no Parameters.h file anymore

Updated the Stimulation protocol to reflect periodic and phase-dependent
stimulation

Signed-off-by: Schellenberger <schellenberger@inb.uni-luebeck.de>

Cortical_Column.cpp

@@ -13,7 +13,7 @@ void Cortical_Column::set_RNG(void) {
 	// get the RNG
 	for (int i=0; i<N; ++i){
 		// add the RNG
-		MTRands.push_back({ENG(rand()), DIST (mphi_c, dphi_c)});
+		MTRands.push_back({ENG(rand()), DIST (mphi, dphi)});
 
 		// get the random number for the first iteration
 		Rand_vars.push_back(MTRands[i]());
@@ -85,14 +85,14 @@ double Cortical_Column::I_ii	(int N) const{
 // Leak current of pyramidal population
 double Cortical_Column::I_L_e	(int N) const{
 	_SWITCH((Ve))
-	double I = gL_e * (var_Ve - E_L_e);
+	double I = g_L * (var_Ve - E_L_e);
 	return I;
 }
 
 // Leak current of inhibitory population
 double Cortical_Column::I_L_i	(int N) const{
 	_SWITCH((Vi))
-	double I = gL_i * (var_Vi - E_L_i);
+	double I = g_L * (var_Vi - E_L_i);
 	return I;
 }
 

Cortical_Column.h

@@ -8,9 +8,9 @@
 #include <boost/random/normal_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
 #include "macros.h"
-#include "parameters.h"
 #include "Thalamic_Column.h"
 using std::vector;
+
 class Thalamic_Column;
 
 /****************************************************************************************************/
@@ -31,26 +31,19 @@ class Cortical_Column {
 public:
 	// Constructors
 	Cortical_Column(void)
-	: Ve		(_INIT(E_L_e)),	Vi 	   	(_INIT(E_L_i)),	Na	 	(_INIT(Na_eq)),
-	  Phi_ee	(_INIT(0.0)), 	Phi_ei 	(_INIT(0.0)), 	Phi_ie 	(_INIT(0.0)), 	Phi_ii	(_INIT(0.0)), 	phi_e	(_INIT(0.0)),
-	  x_ee 		(_INIT(0.0)), 	x_ei   	(_INIT(0.0)),	x_ie   	(_INIT(0.0)), 	x_ii	(_INIT(0.0)), 	y_e		(_INIT(0.0)),
-	  alpha_Na 	(0), 			tau_Na	(0),			g_KNa	(0),	  		theta_e	(0),
-	  sigma_e 	(0), 			dphi_c	(0)
 	{set_RNG();}
 
-	Cortical_Column(double* Par)
-	: Ve		(_INIT(E_L_e)),	Vi 	   	(_INIT(E_L_i)),	Na	 	(_INIT(Na_eq)),
-	  Phi_ee	(_INIT(0.0)), 	Phi_ei 	(_INIT(0.0)), 	Phi_ie 	(_INIT(0.0)), 	Phi_ii	(_INIT(0.0)), 	phi_e	(_INIT(0.0)),
-	  x_ee 		(_INIT(0.0)), 	x_ei   	(_INIT(0.0)),	x_ie   	(_INIT(0.0)), 	x_ii	(_INIT(0.0)), 	y_e		(_INIT(0.0)),
-	  alpha_Na 	(Par[0]), 		tau_Na	(Par[1]),		g_KNa	(Par[2]),	  	theta_e	(Par[3]),
-	  sigma_e 	(Par[4]), 		dphi_c	(Par[5])
+	Cortical_Column(double* Par, double* Con)
+	 :tau_e		(Par[0]), 		theta_e		(Par[1]),		sigma_e 	(Par[2]),		alpha_Na 	(Par[3]),
+	  tau_Na	(Par[4]),	  	g_KNa		(Par[5]),		N_te		(Con[2]), 		N_ti		(Con[3]),
+	  dphi		(Par[6])
 	{set_RNG();}
 
 	// get the pointer to the cortical module
 	void	get_Thalamus(Thalamic_Column& T) {Thalamus = &T;}
 
 	// get axonal flux
-	double get_phi		(int N) const {_SWITCH((phi_e)); return var_phi_e;}
+	double 	get_phi		(int N) const {_SWITCH((phi_e)); return var_phi_e;}
 
 	// Initialize the RNGs
 	void 	set_RNG		(void);
@@ -83,40 +76,97 @@ class Cortical_Column {
 
 private:
 	// Population variables
-	vector<double> 	Ve,			// excitatory membrane voltage
-					Vi,			// inhibitory membrane voltage
-					Na,			// Na concentration
-					Phi_ee,		// PostSP from excitatory to excitatory population
-					Phi_ei,		// PostSP from excitatory to inhibitory population
-					Phi_ie,		// PostSP from inhibitory to excitatory population
-					Phi_ii,		// PostSP from inhibitory to inhibitory population
-					phi_e,		// axonal flux
-					x_ee,		// derivative of Phi_ee
-					x_ei,		// derivative of Phi_ei
-					x_ie,		// derivative of Phi_ie
-					x_ii,		// derivative of Phi_ii
-					y_e;		// derivative of phi_e
-
-	// Adaption parameters
-	double			alpha_Na,	// Sodium influx per spike
-					tau_Na,		// Sodium time constant
-					g_KNa;		// KNa conductance
-
-	// Firing rate parameters
-	double			theta_e,	// pyramidal firing threshold
-					sigma_e;	// pyramidal gain
-
-	// Noise parameters
-	double 			dphi_c;
-
-	// pointer to the thalamic module
-	Thalamic_Column*	Thalamus;
+	vector<double> 	Ve		= _INIT(E_L_e),		// excitatory membrane voltage
+					Vi		= _INIT(E_L_i),		// inhibitory membrane voltage
+					Na		= _INIT(Na_eq),		// Na concentration
+					Phi_ee	= _INIT(0.0),		// PostSP from excitatory to excitatory population
+					Phi_ei	= _INIT(0.0),		// PostSP from excitatory to inhibitory population
+					Phi_ie	= _INIT(0.0),		// PostSP from inhibitory to excitatory population
+					Phi_ii	= _INIT(0.0),		// PostSP from inhibitory to inhibitory population
+					phi_e	= _INIT(0.0),		// axonal flux
+					x_ee	= _INIT(0.0),		// derivative of Phi_ee
+					x_ei	= _INIT(0.0),		// derivative of Phi_ei
+					x_ie	= _INIT(0.0),		// derivative of Phi_ie
+					x_ii	= _INIT(0.0),		// derivative of Phi_ii
+					y_e		= _INIT(0.0);		// derivative of phi_t
 
 	// Random number generators
 	vector<GEN>		MTRands;
 
 	// Container for noise
 	vector<double>	Rand_vars;
+
+	// Declaration and Initialization of parameters
+	// Membrane time in ms
+	const double 	tau_e 		= 30;
+	const double 	tau_i 		= 30;
+
+	// Maximum firing rate in ms^-1
+	const double 	Qe_max		= 30.E-3;
+	const double 	Qi_max		= 60.E-3;
+
+	// Sigmoid threshold in mV
+	const double 	theta_e		= -58.5;
+	const double 	theta_i		= -58.5;
+
+	// Sigmoid gain in mV
+	const double 	sigma_e		= 4;
+	const double 	sigma_i		= 6;
+
+	// Scaling parameter for sigmoidal mapping (dimensionless)
+	const double 	C1          = (3.14159265/sqrt(3));
+
+	// 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 	R_pump   	= 0.09;        	// Na-K pump  constant              in mM/ms
+	const double 	Na_eq    	= 9.5;         	// Na-eq concentration              in mM
+
+	// 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;
+
+
+	// Conductivities in mS/cm^-2
+	// Leak current
+	const double 	g_L    		= 1;
+
+	// KNa current
+	const double	g_KNa		= 1.33;
+
+	// Connectivities (dimensionless)
+	const double 	N_ee		= 120;
+	const double 	N_ei		= 72;
+	const double 	N_ie		= 90;
+	const double 	N_ii		= 90;
+	const double 	N_te		= 10;
+	const double 	N_ti		= 10;
+
+	// Reversal potentials in mV
+	// synaptic
+	const double 	E_AMPA  	= 0;
+	const double 	E_GABA  	= -70;
+
+	// Leak
+	const double 	E_L_e 		= -66;
+	const double 	E_L_i 		= -64;
+
+	// Potassium
+	const double 	E_K    		= -100;
+
+	// Noise parameters in ms^-1
+	const double 	mphi		= 0E-3;
+	const double	dphi		= 30E-3;;
+	double			input		= 0.0;
+
+	/* Pointer to thalamic column */
+	Thalamic_Column* Thalamus;
+
+	friend class 	Stim;
 };
 /****************************************************************************************************/
 /*										 		end			 										*/

Plots.m

@@ -1,37 +1,49 @@
 % mex command is given by: 
 % mex CXXFLAGS="\$CXXFLAGS -std=gnu++0x -fpermissive" TC.cpp Cortical_Column.cpp Thalamic_Column.cpp
 
-function Plots(T, onset)
+function Plots(T)
 
 if nargin == 0
-    Input_N3    = [ 2.6;        % alpha_Na
-                    3;          % tau_Na
-                    1.33	% g_KNa
-                    -63;        % theta_e
-                    8;          % sigma_e
-                    30E-3];     % dphi
+    % fittet input
+    Input_N3    = [ 24;         % tau_e
+                    -64;        % theta_e
+                    8.3;          % sigma_e
+                    3.6;        % alpha_Na
+                    2.7;          % tau_Na
+                    0.63;        % g_KNa
+                    50E-3];      % dphi
                         
                         
-    Input_N2    = [ 2;          % alpha_Na
-                    1;          % tau_Na
-                    1.33;	% g_KNa
+    Input_N2    = [ 30;         % tau_e
                     -58.5;      % theta_e
-                    4.25;       % sigma_e
+                    4.5;          % sigma_e
+                    2;          % alpha_Na
+                    1;          % tau_Na
+                    1.33;       % g_KNa
                     30E-3];     % dphi
 
-    Con     	= [10;		% N_tr
-               	   20;		% N_rt
-                   40];		% N_rr    
-
-    var_stim    = [ 0;          % strength of the stimulus 	in Hz (spikes per second)
-                    0;          % time between stimuli 		in s    
-                    0;          % time until first stimuli 	in s
-                    0];        	% duration of the stimulus 	in ms
-
-    T       	= 30;  		% duration of the simulation
+    Con     	= [0;		% N_et
+               	   0;		% N_er
+                   0;       % N_te
+                   0];		% N_ti    
+
+    % stimulation parameters
+    % first number is the mode of stimulation
+    % 0 == none
+    % 1 == periodic
+    % 2 == phase dependend up state
+    % 3 == phase dependend down state
+    
+    var_stim    = [ 0           % mode of stimulation
+                    100.0;      % strength of the stimulus      in Hz (spikes per second)
+                    100;       	% duration of the stimulus      in ms
+                    8;          % time between stimuli          in s    
+                    550];       % time until stimuli after min 	in ms
+
+    T       	= 60;  		% duration of the simulation
 end
 
-[Ve, Vt] = TC(T, Input_N2, Con, var_stim);
+[Ve, Vt] = TC(T, Input_N3, Con, var_stim);
 
 L        = max(size(Vt));
 timeaxis = linspace(0,T,L);