From 3495bacb8fc43e93fa923974221ec9cca5570281 Mon Sep 17 00:00:00 2001 From: Schellenberger Date: Fri, 31 Oct 2014 11:37:41 +0100 Subject: [PATCH] More parameter tuning Signed-off-by: Schellenberger --- Plots.m | 105 +++++++++++++++++++++++++++------------------------- Thalamic_Column.cpp | 11 ++++-- Thalamic_Column.h | 26 ++++++------- 3 files changed, 75 insertions(+), 67 deletions(-) diff --git a/Plots.m b/Plots.m index 5b9902b..658be16 100644 --- a/Plots.m +++ b/Plots.m @@ -1,82 +1,85 @@ % mex command is given by: % mex CXXFLAGS="\$CXXFLAGS -std=c++11" TC.cpp Cortical_Column.cpp Thalamic_Column.cpp -function Plots(T) - +function Plots(type) if nargin == 0 + type = 1; +end + +if type == 1 - Param_Cortex_N2 = [4.7; % sigma_e - 1.43; % g_KNa + Param_Cortex = [4.7; % sigma_e + 1.33; % g_KNa 120E-3]; % dphi - - Param_Cortex_N3 = [6.3; % sigma_e - 2.; % g_KNa - 120E-3]; % dphi - - - Param_Thalamus_N2 = [0.025; % g_LK_t - 0.025; % g_LK_r - 0.08]; % g_h - + + Param_Thalamus = [0.048; % g_h + 0.0245; % g_LK_t + 0.0245]; % g_LK_r + +else + Param_Cortex = [6.05; % 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 +end - Param_Thalamus_N3 = [0.021; % g_LK_t - 0.021; % g_LK_r - 0.08]; % g_h - - Connectivity = [2.4; % N_et - 2.6; % N_er - 5; % N_te - 10]; % N_ti +Connectivity = [ 3; % N_et + 3; % N_er + 5; % N_te + 10]; % N_ti - % stimulation parameters - % first number is the mode of stimulation - % 0 == none - % 1 == semi-periodic - % 2 == phase dependend +% stimulation parameters +% first number is the mode of stimulation +% 0 == none +% 1 == semi-periodic +% 2 == phase dependend - var_stim = [ 2; % mode of stimulation - 25; % strength of the stimulus in Hz (spikes per second) - 120; % duration of the stimulus in ms - 5; % time between stimulation events in s (ISI) - 0; % range of ISI in s [ISI-range,ISI+range] - 2; % Number of stimuli per event - 950; % time between stimuli within a event in ms - 500]; % time until stimuli after minimum in ms +var_stim = [ 0; % mode of stimulation + 50; % strength of the stimulus in Hz (spikes per second) + 120; % 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 - T = 30; % duration of the simulation -end +T = 30; % duration of the simulation -[Ve, Vt, Marker_Stim] = TC(T, Param_Cortex_N2, Param_Thalamus_N2, Connectivity, var_stim); -%[Ve, Vt, Marker_Stim] = TC(T, Param_Cortex_N3, Param_Thalamus_N3, Connectivity, var_stim); +[Ve, Vt, Marker_Stim] = TC(T, Param_Cortex, Param_Thalamus, Connectivity, var_stim); -L = max(size(Vt)); +L = length(Vt); timeaxis = linspace(0,T,L); figure(1) subplot(211), plot(timeaxis,Ve) -title('Pyramidal membrane voltage'), xlabel('time in s'), ylabel('Ve in mV') +title('Pyramidal membrane voltage'), +xlabel('Time in s'), +ylabel('V_{e} in mV') 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', ':'); changedependvar(hx,'x'); end -hx2 = graph2d.constantline((Marker_Stim/1E2 -var_stim(8)/1E3), 'color', 'red'); -changedependvar(hx2,'x'); subplot(212), plot(timeaxis,Vt) -title('Thalamic relay membrane voltage'), xlabel('time in s'), ylabel('Vt in mV') -ylim(get(gca,'ylim')); +title('Thalamic relay membrane voltage'), +xlabel('Time in s'), +ylabel('V_{t} in mV') +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', ':'); changedependvar(hx,'x'); end -% [Pxx,f] = pwelch(Ve-mean(Ve),hamming(L/30), 4*L/T, 2048, 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/Thalamic_Column.cpp b/Thalamic_Column.cpp index 35acbf5..2b3d421 100644 --- a/Thalamic_Column.cpp +++ b/Thalamic_Column.cpp @@ -160,10 +160,15 @@ double Thalamic_Column::m_inf_h (int N) const{ return h; } -/* Activation time for slow components in TC population after Destexhe 1993 */ +/* Activation time for slow components in TC population*/ double Thalamic_Column::tau_m_h (int N) const{ _SWITCH((Vt)) - double tau = 1. / (exp(-14.59 - 0.086 * var_Vt) + exp(-1.87 + 0.07 * var_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))); + /* Chen2012 */ + //double tau = (20 +1000/(exp((var_Vt + 71.5)/14.2) + exp(-(var_Vt + 89)/11.6))); return tau; } /****************************************************************************************************/ @@ -256,7 +261,7 @@ void Thalamic_Column::set_RK (int N) { 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 * 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); Phi_tt [N] = dt*(var_x_tt); Phi_tr [N] = dt*(var_x_tr); diff --git a/Thalamic_Column.h b/Thalamic_Column.h index abd2f6a..632086f 100644 --- a/Thalamic_Column.h +++ b/Thalamic_Column.h @@ -57,8 +57,8 @@ class Thalamic_Column { /* Constructor for simulation */ Thalamic_Column(double* Param, double* Con) - : g_LK_t (Param[0]), g_LK_r(Param[1]), g_h(Param[2]), - N_et (Con[0]), N_er (Con[1]) + : g_h (Param[0]), g_LK_t (Param[1]), g_LK_r (Param[2]), + N_et (Con[0]), N_er (Con[1]) {set_RNG();} /* Get the pointer to the cortical module */ @@ -155,8 +155,8 @@ class Thalamic_Column { const double theta_r = -58.6; /* Sigmoid gain in mV */ - const double sigma_t = 4; - const double sigma_r = 4; + const double sigma_t = 6; + const double sigma_r = 6; /* Scaling parameter for sigmoidal mapping (dimensionless) */ const double C1 = (3.14159265/sqrt(3)); @@ -174,15 +174,15 @@ class Thalamic_Column { const double g_L_r = 1; /* Potassium leak current */ - const double g_LK_t = 0.02; - const double g_LK_r = 0.02; + const double g_LK_t = 0.024; + const double g_LK_r = 0.024; /* T current */ const double g_T_t = 3; const double g_T_r = 2.3; /* h current */ - const double g_h = 0.05; + const double g_h = 0.048; /* Reversal potentials in mV */ /* Synaptic */ @@ -203,9 +203,9 @@ class Thalamic_Column { const double E_h = -40; /* Calcium parameters */ - const double alpha_Ca = -52E-6; /* influx per spike in nmol */ + const double alpha_Ca = -50E-6; /* influx per spike in nmol */ const double tau_Ca = 10; /* calcium time constant in ms */ - const double Ca_0 = 2.4E-4; /* resting concentration */ + const double Ca_0 = 2E-4; /* resting concentration */ /* I_h activation parameters */ const double k1 = 2.5E7; @@ -222,11 +222,11 @@ class Thalamic_Column { /* Connectivities (dimensionless) */ - const double N_tr = 4; + const double N_tr = 3; const double N_rt = 4; - const double N_rr = 20; - const double N_et = 10; - const double N_er = 10; + const double N_rr = 22; + const double N_et = 5; + const double N_er = 5; /* Pointer to cortical column */ Cortical_Column* Cortex;