Added a tracking function for I_h activation.

Put P_h back for testing.

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

Cortical_Column.h

@@ -92,7 +92,7 @@ class Cortical_Column {
 	void 	add_RK	 	(void);
 
 	/* Data storage access */
-	friend void get_data (int, Cortical_Column&, Thalamic_Column&, _REPEAT(double*, 2));
+	friend void get_data (int, Cortical_Column&, Thalamic_Column&, _REPEAT(double*, 3));
 
 	/* Stimulation protocol access */
 	friend class 	Stim;

Plots.m

@@ -21,9 +21,13 @@ function Plots(type)
                            2.1;          % g_KNa
                            120E-3];     % dphi
                        
-    Param_Thalamus      = [0.055;       % g_h
-                           0.024;       % g_LK_t
-                           0.024];       % g_LK_r                    
+    Param_Thalamus      = [0.042;       % g_h
+                           0.022;       % g_LK_t
+                           0.022;       % g_LK_r  
+                           2.5;         % k1
+                           4;           % k2
+                           1;           % k3
+                           1];          % k4
 end
                         
 Connectivity            = [ 4;          % N_et
@@ -48,13 +52,13 @@ function Plots(type)
 
 T       	= 30;           % duration of the simulation
 
-[Ve, Vt, Marker_Stim] = TC(T, Param_Cortex, Param_Thalamus, Connectivity, var_stim);
+[Ve, Vt, ah, Marker_Stim] = TC(T, Param_Cortex, Param_Thalamus, Connectivity, var_stim);
 
 L        = length(Vt);
 timeaxis = linspace(0,T,L);
 
 figure(1)
-subplot(211), plot(timeaxis,Ve)
+subplot(311), plot(timeaxis,Ve)
 title('Pyramidal membrane voltage'), 
 xlabel('Time in s'), 
 ylabel('V_{e} in mV')
@@ -66,7 +70,7 @@ function Plots(type)
     changedependvar(hx,'x');
 end
 
-subplot(212), plot(timeaxis,Vt)
+subplot(312), plot(timeaxis,Vt)
 title('Thalamic relay membrane voltage'), 
 xlabel('Time in s'), 
 ylabel('V_{t} in mV')
@@ -77,6 +81,18 @@ function Plots(type)
     hx = graph2d.constantline(Marker_Stim/1E2+(i-1)*var_stim(7)/1E3,'ydata', get(gca,'ylim'),'xdata', get(gca,'xlim'), 'color', 'black', 'LineStyle', ':');
     changedependvar(hx,'x');
 end
+
+subplot(313), plot(timeaxis,ah)
+title('Thalamic relay I_h activation'), 
+xlabel('Time in s'), 
+ylabel('m_h in mV')
+xlim([0,30]);
+ylim(get(gca,'ylim'));
+% vertical line for markers
+for i=1:var_stim(6)
+    hx = graph2d.constantline(Marker_Stim/1E2+(i-1)*var_stim(7)/1E3,'ydata', get(gca,'ylim'),'xdata', get(gca,'xlim'), 'color', 'black', 'LineStyle', ':');
+    changedependvar(hx,'x');
+end
 % [Pxx,f]   = pwelch(Ve-mean(Ve),hamming(5*L/T), 2*L/T, [], L/T);
 % n         = find(f<=30, 1, 'last' );
 % 

TC.cpp

@@ -78,26 +78,29 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 	/* Create data containers */
 	mxArray* Ve		= SetMexArray(1, T*res/red);
 	mxArray* Vt		= SetMexArray(1, T*res/red);
+	mxArray* ah		= SetMexArray(1, T*res/red);
 
 	/* Pointer to the actual data block */
 	double* Pr_Ve	= mxGetPr(Ve);
 	double* Pr_Vt	= mxGetPr(Vt);
+	double* Pr_ah	= mxGetPr(ah);
 
 	/* Simulation */
 	int count = 0;
 	for (int t=0; t<Time; ++t) {
 		ODE (Cortex, Thalamus);
 		Stimulation.check_stim(t);
 		if(t>=onset*res && t%red==0){
-			get_data(count, Cortex, Thalamus, Pr_Ve, Pr_Vt);
+			get_data(count, Cortex, Thalamus, Pr_Ve, Pr_Vt, Pr_ah);
 			++count;
 		}
 	}
 
 	/* Output of the simulation */
 	plhs[0] = Ve;
 	plhs[1] = Vt;
-	plhs[2] = Stimulation.get_marker();
+	plhs[2] = ah;
+	plhs[3] = Stimulation.get_marker();
 return;
 }
 /****************************************************************************************************/

Thalamic_Column.cpp

@@ -171,11 +171,17 @@ double Thalamic_Column::tau_m_h	(int N) const{
 }
 
 /* Instantaneous calcium binding onto messenger protein after Chen2012 */
-double Thalamic_Column::P_h	(int N) const{
+double Thalamic_Column::P_H	(int N) const{
 	_SWITCH((Ca))
 	double P_h = k1 * pow(var_Ca, n_P)/(k1*pow(var_Ca, n_P)+k2);
 	return P_h;
 }
+
+/* I_h activation */
+double Thalamic_Column::act_h	(void) const{
+	double activation = m_h[0] + g_inc * m_h2[0];
+	return activation;
+}
 /****************************************************************************************************/
 /*										 		end			 										*/
 /****************************************************************************************************/
@@ -259,14 +265,15 @@ void Thalamic_Column::set_RK (int N) {
 	_SWITCH((Ca)
 			(Phi_tt)(Phi_tr)(Phi_rt)(Phi_rr)(phi)
 			(x_tt)	(x_tr)	(x_rt)	(x_rr)	(y)
-			(h_T_t)	(h_T_r)	(m_h)	(m_h2))
+			(h_T_t)	(h_T_r)	(m_h)	(m_h2) 	(P_h))
 	Vt	  	[N] = dt*(-(I_L_t(N) + I_et(N) + I_it(N))/tau_t - (I_LK_t(N) + I_T_t(N) + I_h(N)));
 	Vr	  	[N] = dt*(-(I_L_r(N) + I_er(N) + I_ir(N))/tau_r - (I_LK_r(N) + I_T_r(N)));
 	Ca		[N] = dt*(alpha_Ca * I_T_t(N) - (var_Ca - Ca_0)/tau_Ca);
 	h_T_t 	[N] = dt*(h_inf_T_t(N) - var_h_T_t)/tau_h_T_t(N);
 	h_T_r 	[N] = dt*(h_inf_T_r(N) - var_h_T_r)/tau_h_T_r(N);
-	m_h 	[N] = dt*((m_inf_h(N) * (1 - var_m_h2) - var_m_h)/tau_m_h(N) - k3 * P_h(N) * var_m_h + k4 * var_m_h2);
-	m_h2 	[N] = dt*(k3 * P_h(N) * var_m_h - k4 * var_m_h2);
+	m_h 	[N] = dt*((m_inf_h(N) * (1 - var_m_h2) - var_m_h)/tau_m_h(N) - k3 * P_H(N) * var_m_h + k4 * var_m_h2);
+	m_h2 	[N] = dt*(k3 * P_H(N) * var_m_h - k4 * var_m_h2);
+	P_h 	[N] = dt*(k1 * pow(var_Ca, n_P) - k2 * var_P_h);
 	Phi_tt	[N] = dt*(var_x_tt);
 	Phi_tr	[N] = dt*(var_x_tr);
 	Phi_rt	[N] = dt*(var_x_rt);
@@ -305,6 +312,7 @@ void Thalamic_Column::add_RK(void) {
 	h_T_r 	[0] += (h_T_r	[1] + h_T_r		[2] * 2 + h_T_r		[3] * 2 + h_T_r		[4])/6;
 	m_h		[0] += (m_h		[1] + m_h		[2] * 2 + m_h		[3] * 2 + m_h		[4])/6;
 	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;
 
 	/* Generate noise for the next iteration */
 	for (unsigned i=0; i<Rand_vars.size(); ++i) {

Thalamic_Column.h

@@ -58,6 +58,8 @@ class Thalamic_Column {
 	/* Constructor for simulation */
 	Thalamic_Column(double* Param, double* Con)
 	: g_LK_t	(Param[1]),	g_LK_r	(Param[2]),	g_h (Param[0]),
+	  k1		(Param[3]*1E7),		k2		(Param[4]*1E-4),
+	  k3		(Param[5]*1E-1),	k4		(Param[6]*1E-3),
 	  N_et 		(Con[0]),	N_er	(Con[1])
 	{set_RNG();}
 
@@ -88,7 +90,8 @@ class Thalamic_Column {
 	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  P_H			(int) const;
+	double  act_h		(void)const;
 
 	/* Deactivation functions */
 	double  h_inf_T_t	(int) const;
@@ -113,7 +116,7 @@ class Thalamic_Column {
 	void 	add_RK	 	(void);
 
 	/* Data storage  access */
-	friend void get_data (int, Cortical_Column&, Thalamic_Column&, _REPEAT(double*, 2));
+	friend void get_data (int, Cortical_Column&, Thalamic_Column&, _REPEAT(double*, 3));
 
 private:
 	/* Population variables */
@@ -133,7 +136,8 @@ class Thalamic_Column {
 					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 	*/
+					m_h2	= _INIT(0.0),		/* activation 	of h   channel bound with protein 	*/
+					P_h		= _INIT(0.0);		/* messenger protein bound with calcium			 	*/
 
 	/* Random number generators */
 	vector<GEN>		MTRands;

saves.h

@@ -30,9 +30,10 @@
 /****************************************************************************************************/
 /*											Save data												*/
 /****************************************************************************************************/
-inline void get_data(int counter, Cortical_Column& Cortex, Thalamic_Column& Thalamus, double* Ve, double* Vt) {
+inline void get_data(int counter, Cortical_Column& Cortex, Thalamic_Column& Thalamus, double* Ve, double* Vt, double* ah) {
 	Ve 	[counter] = Cortex.Ve		[0];
 	Vt 	[counter] = Thalamus.Vt		[0];
+	ah 	[counter] = Thalamus.act_h	();
 }
 /****************************************************************************************************/
 /*										 		end													*/