/NM_TC

Permalink
Browse files

More parameter tuning. 

Changed the conduction speed nu to 60E-3, so that one round is around 30ms
  • Loading branch information...
1 parent 1648694 commit b7992b7fb6289c0eba9c8504615e2186ba592287 @miscco committed Jan 20, 2015
Unified Split
Showing with 196 additions and 19 deletions.
  1. +6 −6 Cortical_Column.h
  2. +6 −6 Plots.m
  3. +6 −6 Test_Parameters.m
  4. +177 −0 Test_Stimulation.m
  5. +1 −1 saves.h
View
12 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 = 2; /* 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 */
@@ -151,7 +151,7 @@ class Cortical_Column {
const double gamma_i = 58.6E-3;
/* Axonal flux time constant */
- const double nu = 120E-3;
+ const double nu = 60E-3;
/* Conductivities in mS/cm^-2 */
/* Leak */
@@ -166,7 +166,7 @@ class Cortical_Column {
const double E_GABA = -70;
/* Leak */
- const double E_L_e = -66;
+ const double E_L_e = -64;
const double E_L_i = -64;
/* Potassium */
@@ -178,10 +178,10 @@ class Cortical_Column {
double input = 0.0;
/* Connectivities (dimensionless) */
- const double N_ee = 125;
+ const double N_ee = 115;
const double N_ei = 72;
- const double N_ie = 100;
- const double N_ii = 100;
+ const double N_ie = 90;
+ const double N_ii = 90;
const double N_te = 0;
const double N_ti = 0;
View
12 Plots.m
@@ -3,19 +3,19 @@
function Plots(type)
if nargin == 0
- type = 1;
+ type = 2;
end
if type == 1
Param_Cortex = [4.7; % sigma_e
- 1.43; % g_KNa
+ 1.33; % g_KNa
120E-3]; % dphi
Param_Thalamus = [0.051; % g_h
0.024]; % g_LK
else
Param_Cortex = [6; % sigma_e
- 2.15; % g_KNa
+ 2.05; % g_KNa
120E-3]; % dphi
Param_Thalamus = [0.051; % g_h
@@ -24,8 +24,8 @@ function Plots(type)
Connectivity = [ 2.5; % N_et
2.5; % N_er
- 6; % N_te
- 15]; % N_ti
+ 5; % N_te
+ 10]; % N_ti
% stimulation parameters
% first number is the mode of stimulation
@@ -49,7 +49,7 @@ function Plots(type)
L = length(Vt);
timeaxis = linspace(0,T,L);
-%figure(1)
+figure(1)
subplot(411), plot(timeaxis,Ve)
title('Pyramidal membrane voltage'),
xlabel('Time in s'),
View
12 Test_Parameters.m
@@ -2,7 +2,7 @@
function Test_Parameters(type)
if nargin == 0
- type = 1;
+ type = 2;
end
@@ -23,7 +23,7 @@ function Test_Parameters(type)
1.33; % g_KNa
120E-3]; % dphi
- Param_Thalamus = [0.052; % g_h
+ Param_Thalamus = [0.051; % g_h
0.024]; % g_LK
fn_data = '/nfshome/schellen/Documents/MATLAB/TC_model/Data/SO_Average_N3';
@@ -34,7 +34,7 @@ function Test_Parameters(type)
xRange = -0.5:0.25:1.5;
else
Param_Cortex = [6; % sigma_e
- 1.92; % g_KNa
+ 2.0; % g_KNa
120E-3]; % dphi
Param_Thalamus = [0.051; % g_h
@@ -51,7 +51,7 @@ function Test_Parameters(type)
Connectivity = [ 2.5; % N_et
2.5; % N_er
5; % N_te
- 5]; % N_ti
+ 15]; % N_ti
load(fn_data);
@@ -62,8 +62,8 @@ function Test_Parameters(type)
% 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
+ 40; % 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
View
177 Test_Stimulation.m
@@ -0,0 +1,177 @@
+% mex command is given by:
+
+function Test_Stimulation(type)
+if nargin == 0
+ type = 1;
+end
+
+
+mex CXXFLAGS="\$CXXFLAGS -std=c++11 -O3 -lgopm" TC.cpp Cortical_Column.cpp Thalamic_Column.cpp
+
+% Path to fieltrip preprocessing function
+if(isempty(strfind(path, '/nfshome/schellen/Documents/MATLAB/Tools/fieldtrip/preproc')))
+ addpath('~/Documents/MATLAB/Tools/fieldtrip/preproc');
+end
+
+% Path to helper function
+if(isempty(strfind(path, '/nfshome/schellen/Documents/MATLAB/Tools/boundedline')))
+ addpath('~/Documents/MATLAB/Tools/boundedline');
+end
+
+Param_Cortex = [6; % sigma_e
+ 2.05; % g_KNa
+ 120E-3]; % dphi
+
+Param_Thalamus = [0.052; % g_h
+ 0.02]; % g_LK
+
+Connectivity = [ 2.5; % N_et
+ 2.5; % N_er
+ 5; % N_te
+ 5]; % N_ti
+
+switch type
+ case 1;
+ fn_data = '/nfshome/schellen/Documents/MATLAB/TC_model/Data/Open_Loop_3Klick.mat';
+ Model_xRange = [-70, -45];
+ Model_Range_FSP = [-0.25, 1.25];
+ Data_xRange = [-80, 25];
+ Data_Range_FSP = [-4, 8];
+ xRange = -0.5:0.5:4;
+ case 2;
+ fn_data = '/nfshome/schellen/Documents/MATLAB/TC_model/Data/Closed_Loop_2Klick.mat';
+ Model_xRange = [-75, -45];
+ Model_Range_FSP = [-0.25, 1.25];
+ Data_xRange = [-130, 80];
+ Data_Range_FSP = [-2, 4];
+ xRange = -1:0.5:3.5;
+ case 3;
+ fn_data = '/nfshome/schellen/Documents/MATLAB/TC_model/Data/Closed_Loop_3Klick.mat';
+ Model_xRange = [-75, -45];
+ Model_Range_FSP = [-0.25, 1.25];
+ Data_xRange = [-130, 80];
+ Data_Range_FSP = [-2, 4];
+ xRange = -1:0.5:3.5;
+ case 4;
+ fn_data = 'Data/ClosedLoopStim.mat';
+ Model_xRange = [-75, -45];
+ Model_Range_FSP = [-0.25, 1.25];
+ Data_xRange = [-80, 50];
+ Data_Range_FSP = [2., 6];
+ xRange = -1:0.5:3;
+end
+
+load(fn_data);
+
+% stimulation parameters
+% first number is the mode of stimulation
+% 0 == none
+% 1 == semi-periodic
+% 2 == phase dependend
+
+var_stim = [ 1; % mode of stimulation
+ 40; % 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
+
+T = 300; % duration of the simulation
+
+[Ve, Vi, Vt, Vr, Marker_Stim] = TC(T, Param_Cortex, Param_Thalamus, Connectivity, var_stim);
+Fs = length(Ve)/T;
+Ve_FSP = ft_preproc_hilbert(ft_preproc_bandpassfilter(Ve, Fs, [12, 15], 513, 'fir'), 'abs').^2;
+
+% Search for peaks
+x_SO = Marker_Stim;
+
+% Remove those events, that are too close to begin/end
+x_SO = x_SO(x_SO<(T-2)*Fs);
+x_SO = x_SO(x_SO> 2*Fs);
+
+% Set the variables
+N_Stim = length(x_SO);
+time_event = linspace(xRange(1), xRange(end), (xRange(end)-xRange(1))*Fs);
+Events = zeros(length(time_event), N_Stim);
+Events_FSP = zeros(length(time_event), N_Stim);
+
+% Segmentation
+for i=1:N_Stim
+ Events(:,i) = Ve ((x_SO(i)+xRange(1)*Fs)+1:(x_SO(i)+xRange(end)*Fs));
+ Events_FSP(:,i) = Ve_FSP((x_SO(i)+xRange(1)*Fs)+1:(x_SO(i)+xRange(end)*Fs));
+end
+
+mean_ERP_model= mean(Events, 2); %#ok<*NASGU>
+mean_FSP_model= mean(Events_FSP,2);
+sd_ERP_model = std (Events, 0, 2);
+sd_FSP_model = std (Events_FSP,0, 2);
+
+% Define handle for plotting
+BL_model =@(y,x) boundedline(y,x(:,1), x(:,2), 'alpha', 'transparency', 0.1, 'r');
+BL_data =@(y,x) boundedline(y,x(:,1), x(:,2), 'alpha', 'transparency', 0.1, 'black');
+
+% Option array for set
+Option_Name = { 'ylim';
+ 'ytick';
+ 'yticklabel';
+ 'ycolor';
+ 'xtick'}';
+
+Option_Model_ERP = {Model_xRange;
+ linspace(Model_xRange(1), Model_xRange(2), 5);
+ linspace(Model_xRange(1), Model_xRange(2), 5);
+ 'black';
+ xRange}'; %#ok<*NBRAK>
+
+Option_Data_ERP = {Data_xRange;
+ linspace(Data_xRange(1), Data_xRange(2), 5);
+ linspace(Data_xRange(1), Data_xRange(2), 5);
+ 'black';
+ xRange}';
+
+Option_Model_FSP = {Model_Range_FSP;
+ linspace(Model_Range_FSP(1), Model_Range_FSP(2), 5);
+ linspace(Model_Range_FSP(1), Model_Range_FSP(2), 5);
+ 'black';
+ xRange}';
+
+Option_Data_FSP = {Data_Range_FSP;
+ linspace(Data_Range_FSP(1), Data_Range_FSP(2), 5);
+ linspace(Data_Range_FSP(1), Data_Range_FSP(2), 5);
+ 'black';
+ xRange}';
+
+figure(1)
+subplot(411)
+plot(linspace(0,30,3000),Ve(101:3100));
+title(['Ve with a mean of :',num2str(mean(Ve))]);
+subplot(412)
+plot(linspace(0,30,3000),Vi(101:3100));
+title(['Vi with a mean of :',num2str(mean(Vi))]);
+subplot(413)
+plot(linspace(0,30,3000),Vt(101:3100));
+title(['Vt with a mean of :',num2str(mean(Vt))]);
+subplot(414)
+plot(linspace(0,30,3000),Vr(101:3100));
+title(['Vr with a mean of :',num2str(mean(Vr))]);
+
+% Create figure
+figure(2)
+clf
+subplot(211)
+[AX1, ~, ~] = plotyy(time_events,[mean_ERP_data, sd_ERP_data],time_events,[mean_ERP_model, sd_ERP_model], BL_data, BL_model);
+set(AX1(1), Option_Name, Option_Data_ERP);
+set(AX1(2), Option_Name, Option_Model_ERP);
+ylabel(AX1(1),'EEG [$\mu$V]');
+ylabel(AX1(2),'$V_{e}$ [mV]');
+title([num2str(N_Stim), ' Events'])
+subplot(212)
+[AX2, ~, ~] = plotyy(time_events,[mean_FSP_data, sd_FSP_data],time_events,[mean_FSP_model, sd_FSP_model], BL_data, BL_model);
+set(AX2(1), Option_Name, Option_Data_FSP);
+set(AX2(2), Option_Name, Option_Model_FSP);
+ylabel(AX2(1),'FSP data [a.u.]');
+ylabel(AX2(2),'$FSP model$ [a.u.]');
+title([num2str(N_Stim), ' Events'])
+end
View
2 saves.h
@@ -32,7 +32,7 @@
/****************************************************************************************************/
inline void get_data(int counter, Cortical_Column& Cortex, Thalamic_Column& Thalamus, double* Ve, double* Vt, double* Ca, double* ah) {
Ve [counter] = Cortex.Ve [0];
- Vt [counter] = Cortex.N_te*Thalamus.phi [0];
+ Vt [counter] = Thalamus.Vt [0];
Ca [counter] = Cortex.Phi_ee [0];
ah [counter] = Thalamus.act_h ();
}

0 comments on commit b7992b7

Please sign in to comment.