diff --git a/Cortical_Column.h b/Cortical_Column.h index 4948214..b196f62 100644 --- a/Cortical_Column.h +++ b/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; diff --git a/Plots.m b/Plots.m index 65be3a8..4cb45c5 100644 --- a/Plots.m +++ b/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' ); % diff --git a/TC.cpp b/TC.cpp index bab0f8f..04edf75 100644 --- a/TC.cpp +++ b/TC.cpp @@ -78,10 +78,12 @@ 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; @@ -89,7 +91,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 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; } } @@ -97,7 +99,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { /* Output of the simulation */ plhs[0] = Ve; plhs[1] = Vt; - plhs[2] = Stimulation.get_marker(); + plhs[2] = ah; + plhs[3] = Stimulation.get_marker(); return; } /****************************************************************************************************/ diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp index 15154a0..7524d81 100644 --- a/Thalamic_Column.cpp +++ b/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 MTRands; diff --git a/saves.h b/saves.h index 3f95a57..bcf9292 100644 --- a/saves.h +++ b/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 */