From b9603ae6fa3bf2ba2c923c8cb48bf0f4afb39ea0 Mon Sep 17 00:00:00 2001 From: Schellenberger Date: Fri, 7 Nov 2014 10:33:43 +0100 Subject: [PATCH] More tuning, changed Ca0 to 2.4. alpha_Ca to 51.8 Signed-off-by: Schellenberger --- Cortical_Column.h | 2 +- Plots.m | 58 +++++++++++++++++++++++++++-------------------------- Stimulation.h | 2 +- Thalamic_Column.cpp | 23 ++++++++++++--------- Thalamic_Column.h | 20 +++++++++--------- 5 files changed, 55 insertions(+), 50 deletions(-) diff --git a/Cortical_Column.h b/Cortical_Column.h index af7af58..4948214 100644 --- a/Cortical_Column.h +++ b/Cortical_Column.h @@ -141,7 +141,7 @@ class Cortical_Column { /* parameters of the firing adaption */ const double alpha_Na = 2; /* Sodium influx per spike in mM ms */ - const double tau_Na = 1.5; /* Sodium time constant in ms */ + const double tau_Na = 1.7; /* 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 */ diff --git a/Plots.m b/Plots.m index 658be16..65be3a8 100644 --- a/Plots.m +++ b/Plots.m @@ -1,35 +1,35 @@ % mex command is given by: -% mex CXXFLAGS="\$CXXFLAGS -std=c++11" TC.cpp Cortical_Column.cpp Thalamic_Column.cpp +% mex CXXFLAGS="\$CXXFLAGS -std=c++11 -O3" TC.cpp Cortical_Column.cpp Thalamic_Column.cpp function Plots(type) if nargin == 0 - type = 1; + type = 2; end if type == 1 Param_Cortex = [4.7; % sigma_e - 1.33; % g_KNa + 1.35; % g_KNa 120E-3]; % dphi - Param_Thalamus = [0.048; % g_h - 0.0245; % g_LK_t - 0.0245]; % g_LK_r + Param_Thalamus = [0.037; % g_h + 0.024; % g_LK_t + 0.024]; % g_LK_r else - Param_Cortex = [6.05; % sigma_e + Param_Cortex = [6.; % sigma_e 2.1; % g_KNa 120E-3]; % dphi - Param_Thalamus = [0.052; % g_h - 0.02; % g_LK_t - 0.02]; % g_LK_r + Param_Thalamus = [0.055; % g_h + 0.024; % g_LK_t + 0.024]; % g_LK_r end -Connectivity = [ 3; % N_et - 3; % N_er +Connectivity = [ 4; % N_et + 4; % N_er 5; % N_te - 10]; % N_ti + 8]; % N_ti % stimulation parameters % first number is the mode of stimulation @@ -37,14 +37,14 @@ function Plots(type) % 1 == semi-periodic % 2 == phase dependend -var_stim = [ 0; % mode of stimulation - 50; % strength of the stimulus in Hz (spikes per second) - 120; % duration of the stimulus in ms +var_stim = [ 2; % mode of stimulation + 80; % strength of the stimulus in Hz (spikes per second) + 150; % duration of the stimulus in ms 5; % time between stimulation events in s (ISI) 0; % range of ISI in s [ISI-range,ISI+range] - 1; % Number of stimuli per event - 950; % time between stimuli within a event in ms - 5]; % time until stimuli after minimum in ms + 3; % Number of stimuli per event + 1050; % time between stimuli within a event in ms + 450]; % time until stimuli after minimum in ms T = 30; % duration of the simulation @@ -58,10 +58,11 @@ function Plots(type) title('Pyramidal membrane voltage'), xlabel('Time in s'), ylabel('V_{e} in mV') -ylim([-80, -40]) +xlim([0,30]); +ylim([-80, -40]); % 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'), 'color', 'black', 'LineStyle', ':'); + 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 @@ -69,17 +70,18 @@ function Plots(type) title('Thalamic relay membrane voltage'), xlabel('Time in s'), ylabel('V_{t} in mV') +xlim([0,30]); ylim([-70,-35]); % 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'), 'color', 'black', 'LineStyle', ':'); + 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' ); - -figure(2) -plot(f(1:n),log(Pxx(1:n))) -title('Powerspectrum with pwelch'), xlabel('frequency in Hz'), ylabel('Power (log)') +% [Pxx,f] = pwelch(Ve-mean(Ve),hamming(5*L/T), 2*L/T, [], L/T); +% n = find(f<=30, 1, 'last' ); +% +% figure(2) +% plot(f(1:n),log(Pxx(1:n))) +% title('Powerspectrum with pwelch'), xlabel('frequency in Hz'), ylabel('Power (log)') %save('Timeseries', 'Ve', 'Vt'); end \ No newline at end of file diff --git a/Stimulation.h b/Stimulation.h index 349fe7b..5b23547 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 = -72; + double threshold = -68.5; /* Internal variables */ /* Simulation on for TRUE and off for FALSE */ diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp index 2b3d421..15154a0 100644 --- a/Thalamic_Column.cpp +++ b/Thalamic_Column.cpp @@ -160,17 +160,22 @@ double Thalamic_Column::m_inf_h (int N) const{ return h; } -/* Activation time for slow components in TC population*/ +/* Activation time for slow components in TC population after Chen2012 */ double Thalamic_Column::tau_m_h (int N) const{ _SWITCH((Vt)) - /* Destexhe 1993 */ - //double tau = 1. / (exp(-14.59 - 0.086 * var_Vt) + exp(-1.87 + 0.07 * var_Vt)); /* Bazhenov1998 */ - double tau = (5.3 + 267/(exp((var_Vt + 71.5)/14.2) + exp(-(var_Vt + 89)/11.6))); + //double tau = (5.3 + 267/(exp((var_Vt + 71.5)/14.2) + exp(-(var_Vt + 89)/11.6))); /* Chen2012 */ - //double tau = (20 +1000/(exp((var_Vt + 71.5)/14.2) + exp(-(var_Vt + 89)/11.6))); + double tau = (20 + 1000/(exp((var_Vt + 71.5)/14.2) + exp(-(var_Vt + 89)/11.6))); return tau; } + +/* Instantaneous calcium binding onto messenger protein after Chen2012 */ +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; +} /****************************************************************************************************/ /* end */ /****************************************************************************************************/ @@ -254,15 +259,14 @@ 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) (P_h)) + (h_T_t) (h_T_r) (m_h) (m_h2)) 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 * var_P_h * var_m_h + k4 * var_m_h2); - m_h2 [N] = dt*(k3 * k1 * pow(var_Ca, n_P)/(k1*pow(var_Ca, n_P)+k2) * var_m_h - k4 * var_m_h2); - P_h [N] = dt*(k1 * pow(var_Ca, n_P) * (1 - var_P_h) - k2 * var_P_h); + 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); Phi_tt [N] = dt*(var_x_tt); Phi_tr [N] = dt*(var_x_tr); Phi_rt [N] = dt*(var_x_rt); @@ -301,7 +305,6 @@ 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; @@ -174,15 +174,15 @@ class Thalamic_Column { const double g_L_r = 1; /* Potassium leak current */ - const double g_LK_t = 0.024; - const double g_LK_r = 0.024; + const double g_LK_t = 0.02; + const double g_LK_r = 0.02; /* T current */ const double g_T_t = 3; const double g_T_r = 2.3; /* h current */ - const double g_h = 0.048; + const double g_h = 0.07; /* Reversal potentials in mV */ /* Synaptic */ @@ -203,9 +203,9 @@ class Thalamic_Column { const double E_h = -40; /* Calcium parameters */ - const double alpha_Ca = -50E-6; /* influx per spike in nmol */ + 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 = 2E-4; /* resting concentration */ + const double Ca_0 = 2.4E-4; /* resting concentration */ /* I_h activation parameters */ const double k1 = 2.5E7; @@ -223,8 +223,8 @@ class Thalamic_Column { /* Connectivities (dimensionless) */ const double N_tr = 3; - const double N_rt = 4; - const double N_rr = 22; + const double N_rt = 5; + const double N_rr = 16; const double N_et = 5; const double N_er = 5;