diff --git a/Cortical_Column.cpp b/Cortical_Column.cpp index 4fe4ef9..89467e1 100644 --- a/Cortical_Column.cpp +++ b/Cortical_Column.cpp @@ -178,9 +178,9 @@ void Cortical_Column::set_RK (int N) { s_gp[N+1] = s_gp[0] + A[N] * dt*(x_gp[N]); s_gi[N+1] = s_gi[0] + A[N] * dt*(x_gi[N]); y [N+1] = y [0] + A[N] * dt*(x [N]); - x_ep[N+1] = x_ep[0] + A[N] * dt*(pow(gamma_e, 2) * (N_ep * get_Qp(N) + N_tp * Thalamus->y[N] - s_ep[N]) - 2 * gamma_e * x_ep[N]) + noise_xRK(N, 0); - x_ei[N+1] = x_ei[0] + A[N] * dt*(pow(gamma_e, 2) * (N_ei * get_Qp(N) + N_ti * Thalamus->y[N] - s_ei[N]) - 2 * gamma_e * x_ei[N]) + noise_xRK(N, 1) ; - x_gp[N+1] = x_gp[0] + A[N] * dt*(pow(gamma_g, 2) * (N_ip * get_Qi(N) - s_gp[N]) - 2 * gamma_g * x_gp[N]); + x_ep[N+1] = x_ep[0] + A[N] * dt*(pow(gamma_e, 2) * (N_pp * get_Qp(N) + N_pt * Thalamus->y[N] - s_ep[N]) - 2 * gamma_e * x_ep[N]) + noise_xRK(N, 0); + x_ei[N+1] = x_ei[0] + A[N] * dt*(pow(gamma_e, 2) * (N_ip * get_Qp(N) + N_it * Thalamus->y[N] - s_ei[N]) - 2 * gamma_e * x_ei[N]) + noise_xRK(N, 1) ; + x_gp[N+1] = x_gp[0] + A[N] * dt*(pow(gamma_g, 2) * (N_pi * get_Qi(N) - s_gp[N]) - 2 * gamma_g * x_gp[N]); x_gi[N+1] = x_gi[0] + A[N] * dt*(pow(gamma_g, 2) * (N_ii * get_Qi(N) - s_gi[N]) - 2 * gamma_g * x_gi[N]); x [N+1] = x [0] + A[N] * dt*(pow(nu, 2) * ( get_Qp(N) - y [N]) - 2 * nu * x [N]); } diff --git a/Cortical_Column.h b/Cortical_Column.h index f2eea2c..90c7420 100644 --- a/Cortical_Column.h +++ b/Cortical_Column.h @@ -64,8 +64,8 @@ class Cortical_Column { {set_RNG();} Cortical_Column(double* Param, double* Con) - :sigma_p (Param[0]), g_KNa (Param[1]), dphi (Param[2]), - N_tp (Con[2]), N_ti (Con[3]) + :sigma_p (Param[0]), g_KNa (Param[1]), dphi (Param[2]), + N_pt (Con[2]), N_it (Con[3]) {set_RNG();} /* Connect to the thalamic module */ @@ -170,12 +170,12 @@ class Cortical_Column { double input = 0.0; /* Connectivities (dimensionless) */ - const double N_ep = 115; - const double N_ei = 72; - const double N_ip = 90; + const double N_pp = 115; + const double N_ip = 72; + const double N_pi = 90; const double N_ii = 90; - const double N_tp = 2.5; - const double N_ti = 2.5; + const double N_pt = 2.5; + const double N_it = 2.5; /* Pointer to thalamic column */ Thalamic_Column* Thalamus; @@ -192,18 +192,18 @@ class Cortical_Column { /* Population variables */ vector Vp = _INIT(E_L_p), /* excitatory membrane voltage */ - Vi = _INIT(E_L_g), /* inhibitory membrane voltage */ - Na = _INIT(Na_eq), /* Na concentration */ - s_ep = _INIT(0.0), /* PostSP from excitatory to excitatory population */ - s_ei = _INIT(0.0), /* PostSP from excitatory to inhibitory population */ - s_gp = _INIT(0.0), /* PostSP from inhibitory to excitatory population */ - s_gi = _INIT(0.0), /* PostSP from inhibitory to inhibitory population */ - y = _INIT(0.0), /* axonal flux */ - x_ep = _INIT(0.0), /* derivative of s_ep */ - x_ei = _INIT(0.0), /* derivative of s_ei */ - x_gp = _INIT(0.0), /* derivative of s_gp */ - x_gi = _INIT(0.0), /* derivative of s_gi */ - x = _INIT(0.0); /* derivative of y */ + Vi = _INIT(E_L_g), /* inhibitory membrane voltage */ + Na = _INIT(Na_eq), /* Na concentration */ + s_ep = _INIT(0.0), /* PostSP from excitatory to excitatory population */ + s_ei = _INIT(0.0), /* PostSP from excitatory to inhibitory population */ + s_gp = _INIT(0.0), /* PostSP from inhibitory to excitatory population */ + s_gi = _INIT(0.0), /* PostSP from inhibitory to inhibitory population */ + y = _INIT(0.0), /* axonal flux */ + x_ep = _INIT(0.0), /* derivative of s_ep */ + x_ei = _INIT(0.0), /* derivative of s_ei */ + x_gp = _INIT(0.0), /* derivative of s_gp */ + x_gi = _INIT(0.0), /* derivative of s_gi */ + x = _INIT(0.0); /* derivative of y */ }; /****************************************************************************************************/ /* end */ diff --git a/Data_Storage.h b/Data_Storage.h index 06e79c6..6d2caa3 100644 --- a/Data_Storage.h +++ b/Data_Storage.h @@ -39,7 +39,7 @@ /* Save data */ /****************************************************************************************************/ void get_data(int counter, Cortical_Column& Cortex, Thalamic_Column& Thalamus, vector Data) { - Data[0][counter] = Cortex.Ve [0]; + Data[0][counter] = Cortex.Vp [0]; Data[1][counter] = Thalamus.Vt [0]; Data[2][counter] = Thalamus.Ca [0]; Data[3][counter] = Thalamus.act_h (); diff --git a/Plots.m b/Plots.m index 6334f90..dfcf2da 100644 --- a/Plots.m +++ b/Plots.m @@ -9,25 +9,25 @@ function Plots(type) %mex CXXFLAGS="\$CXXFLAGS -std=c++11 -O3" TC_mex.cpp Cortical_Column.cpp Thalamic_Column.cpp if type == 1 - Param_Cortex = [4.7; % sigma_e + Param_Cortex = [4.7; % sigma_p 1.33; % g_KNa - 20E-1]; % dphi + 2.]; % dphi - Param_Thalamus = [0.048; % g_h - 0.024]; % g_LK + Param_Thalamus = [0.034; % g_LK + 0.052]; % g_h else - Param_Cortex = [6; % sigma_e - 2.05; % g_KNa - 20E-1]; % dphi + Param_Cortex = [5.8; % sigma_p + 1.8; % g_KNa + 2.]; % dphi - Param_Thalamus = [0.048; % g_h - 0.02]; % g_LK + Param_Thalamus = [0.038; % g_LK + 0.052]; % g_h end -Connectivity = [ 2.6; % N_et - 2.6; % N_er - 5; % N_te - 10]; % N_ti +Connectivity = [ 2.5; % N_tp + 2.5; % N_rp + 5; % N_pt + 10]; % N_it % stimulation parameters % first number is the mode of stimulation @@ -35,27 +35,27 @@ function Plots(type) % 1 == semi-periodic % 2 == phase dependend -var_stim = [ 0; % mode of stimulation - 60; % strength of the stimulus in Hz (spikes per second) - 120; % duration of the stimulus in ms +var_stim = [ 2; % mode of stimulation + 200; % strength of the stimulus in Hz (spikes per second) + 100; % duration of the stimulus in ms 5; % time between stimulation events in s (ISI) 0; % range of ISI in s [ISI-range,ISI+range] 3; % Number of stimuli per event 1050; % time between stimuli within a event in ms - 450]; % time until stimuli after minimum in ms + 400]; % time until stimuli after minimum in ms T = 30; % duration of the simulation -[Ve, Vt, Ca, ah, Marker_Stim] = TC_mex(T, Param_Cortex, Param_Thalamus, Connectivity, var_stim); +[Vp, Vt, Ca, ah, Marker_Stim] = TC_mex(T, Param_Cortex, Param_Thalamus, Connectivity, var_stim); L = length(Vt); timeaxis = linspace(0,T,L); figure(2) -subplot(411), plot(timeaxis,Ve) +subplot(411), plot(timeaxis,Vp) title('Pyramidal membrane voltage'), xlabel('Time in s'), -ylabel('V_{e} in mV') +ylabel('V_{p} in mV') xlim([0,30]); ylim([-80, -40]); % vertical line for markers @@ -69,7 +69,7 @@ function Plots(type) xlabel('Time in s'), ylabel('V_{t} in mV') xlim([0,30]); -ylim([-80,-35]); +%ylim([-80,-35]); % vertical line for markers diff --git a/Stimulation.h b/Stimulation.h index b6e74fe..aa1ef64 100644 --- a/Stimulation.h +++ b/Stimulation.h @@ -89,7 +89,7 @@ class Stim { int time_between_stimuli = 1050E1; /* Threshold for phase dependent stimulation */ - double threshold = -68.5; + double threshold = -68; /* Internal variables */ /* Simulation on for TRUE and off for FALSE */ @@ -135,7 +135,7 @@ class Stim { int count_pause = 0; /* Old voltage value for minimum detection */ - double Ve_old = 0.0; + double Vp_old = 0.0; /* Pointer to columns */ Cortical_Column* Cortex; @@ -194,11 +194,8 @@ void Stim::setup (double* var_stim) { /* If ISI is random create RNG */ if (ISI_range != 0){ - /* Create the generator */ - Generator = ENG(rand()); - - /* Combine RNG with distribution */ - Uniform_Distribution = DIST_int(ISI-ISI_range, ISI+ISI_range); + /* Generate uniform distribution */ + Uniform_Distribution = random_stream_uniform_int(ISI-ISI_range, ISI+ISI_range); } } else { /* In case of phase dependent stimulation, time_to_stim is the time from minimum detection to start of stimulation */ @@ -236,7 +233,7 @@ void Stim::check_stim (int time) { } /* After last stimulus in event update the timer with respect to (random) ISI*/ else { - time_to_stimuli += (ISI_range==0)? ISI : Uniform_Distribution(Generator); + time_to_stimuli += (ISI_range==0)? ISI : Uniform_Distribution(); /* Reset the stimulus counter for next stimulation event */ count_stimuli = 1; @@ -248,19 +245,19 @@ void Stim::check_stim (int time) { case 2: /* Search for threshold */ if(!stimulation_started && !minimum_found && !threshold_crossed && time>onset_correction && !stimulation_paused) { - if(Cortex->Ve[0]<=threshold) { + if(Cortex->Vp[0]<=threshold) { threshold_crossed = true; } } /* Search for minimum */ if(threshold_crossed) { - if(Cortex->Ve[0]>Ve_old) { + if(Cortex->Vp[0]>Vp_old) { threshold_crossed = false; minimum_found = true; - Ve_old = 0; + Vp_old = 0; } else { - Ve_old = Cortex->Ve[0]; + Vp_old = Cortex->Vp[0]; } } diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp index 62e7f29..d021df3 100644 --- a/Thalamic_Column.cpp +++ b/Thalamic_Column.cpp @@ -89,18 +89,18 @@ double Thalamic_Column::I_et (int N) const{ /* Inhibitory input to TC population */ double Thalamic_Column::I_gt (int N) const{ - double I = g_AMPA * s_rt[N]* (Vt[N]- E_GABA); + double I = g_GABA * s_gt[N]* (Vt[N]- E_GABA); return I; } /* Excitatory input to RE population */ double Thalamic_Column::I_er (int N) const{ - double I = g_GABA * s_er[N]* (Vr[N]- E_AMPA); + double I = g_AMPA * s_er[N]* (Vr[N]- E_AMPA); return I; } /* Inhibitory input to RE population */ double Thalamic_Column::I_gr (int N) const{ - double I = g_GABA * s_rr[N]* (Vr[N]- E_GABA); + double I = g_GABA * s_gr[N]* (Vr[N]- E_GABA); return I; } /****************************************************************************************************/ @@ -135,15 +135,15 @@ double Thalamic_Column::h_inf_T_r (int N) const{ return h; } -/* 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{ - double tau = (30.8 + (211.4 + exp((Vt[N]+115.2)/5))/(1 + exp((Vt[N]+86)/3.2)))/pow(3, 1.2); + double tau = (30.8 + (211.4 + exp((Vt[N]+115.2)/5))/(1 + exp((Vt[N]+86)/3.2)))/3.7371928; return tau; } /* 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)))/pow(3, 1.2); + double tau = (85 + 1/(exp((Vr[N]+48)/4) + exp(-(Vr[N]+407)/50)))/3.7371928; return tau; } /****************************************************************************************************/ @@ -169,7 +169,7 @@ 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 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]/(k1* Ca[N]* Ca[N]* Ca[N]+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; } @@ -252,23 +252,23 @@ double Thalamic_Column::noise_aRK(int M) const{ /****************************************************************************************************/ 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 - (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 - (I_LK_r(N) + I_T_r(N))); + 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); - 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_rt [N+1] = s_rt [0] + A[N]*dt*(x_rt[N]); - s_rr [N+1] = s_rr [0] + A[N]*dt*(x_rr[N]); - y [N+1] = y [0] + A[N]*dt*(x [N]); - x_et [N+1] = x_et [0] + A[N]*dt*(pow(gamma_e, 2) * ( + N_pt * Cortex->y[N] - 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) + N_pr * Cortex->y[N] - s_er[N]) - 2 * gamma_e * x_er[N]); - x_rt [N+1] = x_rt [0] + A[N]*dt*(pow(gamma_g, 2) * (N_rt * get_Qr(N) - s_rt[N]) - 2 * gamma_g * x_rt[N]); - x_rr [N+1] = x_rr [0] + A[N]*dt*(pow(gamma_g, 2) * (N_rr * get_Qr(N) - s_rr[N]) - 2 * gamma_g * x_rr[N]); - x [N+1] = x [0] + A[N]*dt*(pow(nu, 2) * ( get_Qt(N) - y [N]) - 2 * nu * x [N]); 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]); + y [N+1] = y [0] + A[N]*dt*(x [N]); + x_et [N+1] = x_et [0] + A[N]*dt*(pow(gamma_e, 2) * ( + N_tp * Cortex->y[N] - 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_rt * get_Qt(N) + N_rp * Cortex->y[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_tr * 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]); + x [N+1] = x [0] + A[N]*dt*(pow(nu, 2) * ( get_Qt(N) - y [N]) - 2 * nu * x [N]); } /****************************************************************************************************/ /* end */ @@ -284,13 +284,13 @@ void Thalamic_Column::add_RK(void) { 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_rt [0] =(-3*s_rt [0] + 2*s_rt [1] + 4*s_rt [2] + 2* s_rt [3] + s_rt [4])/6; - s_rr [0] =(-3*s_rr [0] + 2*s_rr [1] + 4*s_rr [2] + 2* s_rr [3] + s_rr [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; y [0] =(-3*y [0] + 2*y [1] + 4*y [2] + 2* y [3] + y [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_rt [0] =(-3*x_rt [0] + 2*x_rt [1] + 4*x_rt [2] + 2* x_rt [3] + x_rt [4])/6; - x_rr [0] =(-3*x_rr [0] + 2*x_rr [1] + 4*x_rr [2] + 2* x_rr [3] + x_rr [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; x [0] =(-3*x [0] + 2*x [1] + 4*x [2] + 2* x [3] + x [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; diff --git a/Thalamic_Column.h b/Thalamic_Column.h index cd4421c..d67529e 100644 --- a/Thalamic_Column.h +++ b/Thalamic_Column.h @@ -61,8 +61,8 @@ class Thalamic_Column { /* Constructor for simulation */ Thalamic_Column(double* Param, double* Con) - : g_LK (Param[1]), g_h (Param[0]), - N_pt (Con[0]), N_pr (Con[1]) + : g_LK (Param[0]), g_h (Param[1]), + N_tp (Con[0]), N_rp (Con[1]) {set_RNG();} /* Get the pointer to the cortical module */ @@ -139,8 +139,8 @@ class Thalamic_Column { const double theta_r = -58.5; /* Sigmoid gain in mV */ - const double sigma_t = 6; - const double sigma_r = 6; + const double sigma_t = 6.; + const double sigma_r = 6.; /* Scaling parameter for sigmoidal mapping (dimensionless) */ const double C1 = (M_PI/sqrt(3)); @@ -152,6 +152,9 @@ class Thalamic_Column { /* Axonal flux time constant in ms^-1*/ const double nu = 120E-3; + /* Membrane capacitance in muF/cm^2 */ + const double C_m = 1.; + /* Leak weight in aU */ const double g_L = 1.; @@ -203,17 +206,17 @@ class Thalamic_Column { /* Noise parameters in ms^-1 */ const double mphi = 0E-3; - const double dphi = 20E-1; + const double dphi = 20E-3; double input = 0.0; /* Connectivities (dimensionless) */ - const double N_tr = 3; - const double N_rt = 5; - const double N_rr = 19; + const double N_rt = 3.; + const double N_tr = 5.; + const double N_rr = 25.; /* Connectivities from cortex (dimensionless) */ - const double N_pt = 2.6; - const double N_pr = 2.6; + const double N_tp = 2.6; + const double N_rp = 2.6; /* Pointer to cortical column */ Cortical_Column* Cortex; @@ -230,22 +233,22 @@ class Thalamic_Column { /* Population variables */ vector 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_rt = _INIT(0.0), /* PostSP from RE population to TC population */ - s_rr = _INIT(0.0), /* PostSP from RE population to RE population */ - y = _INIT(0.0), /* axonal flux */ - x_et = _INIT(0.0), /* derivative of s_et */ - x_er = _INIT(0.0), /* derivative of s_er */ - x_rt = _INIT(0.0), /* derivative of s_rt */ - x_rr = _INIT(0.0), /* derivative of s_rr */ - x = _INIT(0.0), /* derivative of y */ - 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 */ + 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 */ + y = _INIT(0.0), /* axonal flux */ + 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 */ + x = _INIT(0.0), /* derivative of y */ + 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 */ }; /****************************************************************************************************/ /* end */