Added the matlab functions to create the plots/data as well as the xp…

…paut ode file

Figures/Data_Thalamus.m

@@ -0,0 +1,63 @@
+function Data_Thalamus()
+% This function creates the data depicted in Figure 3 of
+%
+% A thalamocortical neural mass model of the EEG during NREM sleep and its
+% response to auditory stimulation.
+% M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz,
+% JC Claussen.
+% PLoS Computational Biology (in review).
+
+% Move to source folder(assuming it contains the figures folder
+cd ..;
+
+% Check if the executable exists and compile if needed
+if(exist('Thalamus_mex.mesa64', 'file')==0)
+    mex CXXFLAGS="\$CXXFLAGS -std=c++11 -O3" Thalamus_mex.cpp Thalamic_Column.cpp;
+end
+
+% Add the path to the simulation routine
+addpath(pwd);
+
+% Go back into figures folder
+cd Figures;
+
+for i=1:6
+    switch (i)
+        case 1
+            Param_Thalamus = [0.018;      % g_LK
+                              0.062];     % g_h
+            Name           = 'I';
+            
+        case 2
+            Param_Thalamus = [0.032;      % g_LK
+                              0.062];     % g_h
+            Name           = 'II';
+            
+        case 3
+            Param_Thalamus = [0.025;      % g_LK
+                              0.025];     % g_h
+            Name           = 'III';
+            
+        case 4
+            Param_Thalamus = [0.04;       % g_LK
+                              0.066];     % g_h
+            Name           = 'IV';
+            
+        case 5
+            Param_Thalamus = [0.052;      % g_LK
+                              0.066];     % g_h
+            Name           = 'V';
+            
+        case 6
+            Param_Thalamus = [0.052;      % g_LK
+                              0.04];      % g_h
+            Name           = 'VI';
+    end
+    
+    T            = 15; % Duration of the simulation
+    [Vt, Vr, ah] = Thalamus_mex(T, Param_Thalamus); %#ok<ASGLU>
+    
+    % Save the data
+    save(['Data/Timeseries_Thalamus_',Name], 'Vt', 'Vr', 'ah', 'T');
+end
+end

Figures/Plot_Bifurcation_2D.m

@@ -0,0 +1,77 @@
+function Plot_Bifurcation_2D()
+% This function creates Figure 2 of
+%
+% A thalamocortical neural mass model of the EEG during NREM sleep and its
+% response to auditory stimulation.
+% M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz,
+% JC Claussen.
+% PLoS Computational Biology (in review).
+%
+% Please note that the coloring of the two spindle regions was done
+% manually
+
+% Load the data from xppaut
+fname_Hopf  = 'Data/Bifurcation_Thalamus_Hopf.dat';
+fname_Torus = 'Data/Bifurcation_Thalamus_Torus.dat';
+fname_LP    = 'Data/Bifurcation_Thalamus_LP.dat';
+fname_PD    = 'Data/Bifurcation_Thalamus_PD.dat';
+
+% Create the figure
+close all;
+figure(1)
+clf, shg
+
+% Plot the data
+plotxppaut(fname_PD,'Green','-');
+plotxppaut(fname_Torus,'Blue','-');
+plotxppaut(fname_Hopf,'Red','-');
+plotxppaut(fname_LP,'Black','-');
+xlabel('\bar{g}_{LK} [mS/cm^2]')
+ylabel('\bar{g}_{h} [mS/cm^2]')
+set(gca, 'XLim',[0, 0.08], ...
+         'XTick',0:0.02:0.08, ...
+         'YLim', [0, 0.09],...
+         'YTick',0:0.02:0.08);
+     
+% Clean up the plot a bit for the legend
+lines = get(gca, 'Children');
+id = 1:numel(lines);
+% Keep the real lines
+id([1 3 5 7]) = [];
+delete(lines(id))
+
+% Add legends
+legend('Period-doubling', 'Saddle-node', 'Torus', 'Hopf', 'location', 'northoutside', 'orientation','horizontal');
+
+% Set the line style
+lines = get(gca, 'Children');
+for i=1:numel(lines)
+    set(lines(i), 'Marker', 'none');
+end
+
+% Plot the lines indicating the spindle regime
+hold on
+plot([0, 0.02932],[0.06357,0.06357],'--','Color',[0.7,0.7,0.7],'LineWidth',0.5)
+plot([0.0367, 0.08],[0.06357,0.06357],'--','Color',[0.7,0.7,0.7],'LineWidth',0.5)
+plot([0.02932, 0.03613],[0.06357,0.06357],'-','Color',[0.7,0.7,0.7],'LineWidth',0.5)
+
+plot([0.03613,0.03613], [0, 0.0412],'--','Color',[0.7,0.7,0.7],'LineWidth',0.5)
+plot([0.03613,0.03613], [0.064, 0.1],'--','Color',[0.7,0.7,0.7],'LineWidth',0.5)
+plot([0.03613,0.03613], [0.04146, 0.06357],'-','Color',[0.7,0.7,0.7],'LineWidth',0.5)
+
+% Delta regime
+plot([0.044,0.044], ylim,'--','Color',[0.7,0.7,0.7],'LineWidth',0.5)
+
+% Indicate the parameter settings for the example time series
+g_LK = [0.018, 0.0316, 0.052,  0.052, 0.025, 0.04];
+g_h  = [0.062,  0.062,  0.066, 0.04, 0.025, 0.066];
+plot(g_LK, g_h,'*', 'Color', 'black');
+text(g_LK(1), g_h(1), 'S_{I}', 'horizontal','right', 'vertical','bottom')
+text(g_LK(2), g_h(2), 'S_{II}', 'horizontal','center', 'vertical','bottom')
+text(g_LK(3), g_h(3), 'D_{I}', 'horizontal','right', 'vertical','bottom')
+text(g_LK(4), g_h(4), 'D_{II}', 'horizontal','right', 'vertical','bottom')
+text(g_LK(5), g_h(5), 'C_{I}', 'horizontal','right', 'vertical','bottom')
+text(g_LK(6), g_h(6), 'C_{II}', 'horizontal','center', 'vertical','bottom')
+hold off
+
+end

Figures/Plot_Timeseries_Thalamus.m

@@ -0,0 +1,70 @@
+function Plot_Timeseries_Thalamus()
+% This function creates Figure 2 of
+%
+% A thalamocortical neural mass model of the EEG during NREM sleep and its
+% response to auditory stimulation.
+% M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz,
+% JC Claussen.
+% PLoS Computational Biology (in review).
+
+% Path to panel tool
+if(isempty(strfind(path, [pwd, '/Tools'])))
+    addpath([pwd, '/Tools']);
+end
+
+% Names of the individual files
+Names    = {'I', 'II', 'III', 'IV', 'V', 'VI'};
+Labels   = {'S_{I}', 'S_{II}', 'C_{I}', 'C_{II}', 'D_{I}', 'D_{II}'};
+
+T        = 15;
+timeaxis = linspace(0,T,T*100);
+
+% Create figure
+figure(1);
+clf
+
+% Create panel
+p = panel('no-manage-font');
+p.pack(6,2);
+
+% set margins
+p.de.margintop = 1;
+p.de.marginbottom = 6;
+p.de.marginleft= 5;
+p.de.marginright= 5;
+
+for i=1:numel(Names)
+    load(['Data/Timeseries_Thalamus_',Names{i}]);
+    % Thalamic relay membrane voltage
+    p(i,1).select();
+    plot(timeaxis, Vt, 'Color', 'black');
+    set(gca, 'XTickLabel', [], 'YTick', [-80, -50, -20], 'YLim', [-80, -20]);
+    if(i==1)
+        title('Thalamic relay population');
+    end
+    if(i==6)
+    xlabel('Time [s]'); 
+    set(gca, 'XTick', 0:5:15); 
+    set(gca, 'XTick', 0:5:15, 'XTickLabel', 0:5:15);
+    end
+    
+    % Reticular membrane voltage
+    ylabel('V_{t} [mV]');
+    p(i,2).select();
+    plot(timeaxis, Vr, 'Color', 'black');
+    set(gca, 'XTickLabel', [], 'YAxisLocation', 'right', 'YTick', [-80, -50, -20], 'YLim', [-80, -20]);
+    if(i==1)
+        title('Thalamic reticular population');
+    end
+    if(i>=5)
+    set(gca, 'YTick', [-100, 0, 100], 'YLim', [-100, 100]);
+    end        
+    if(i==6)
+    xlabel('Time [s]'); 
+    set(gca, 'XTick', 0:5:15, 'XTickLabel', 0:5:15);
+    end     
+    ylabel('V_{r} [mV]');   
+    % Add the title
+    p(i).title(Labels(i));
+end
+end

Figures/Thalamus.ode

@@ -0,0 +1,123 @@
+#Thalamic model
+
+param g_LK=0.0
+param g_h=0.025
+
+#Time constants for excitatory and inhibitory neurons in ms
+!tau=20
+
+#Conductances 
+!g_L=1
+
+!g_T_t=3
+!g_T_r=2.3
+
+!g_AMPA=1
+!g_GABA=1
+
+#Reversal potential of excitatory and inhibitory synapses in mV
+!E_AMPA=0
+!E_GABA=-70
+
+# Reversal potential of other currents
+!E_L=-70
+!E_K=-100 
+!E_Ca=120
+!E_h=-40
+
+#Maximum firing rate in ms^-1
+!Q_max=400E-3
+
+#Sigmoid threshold in mV
+!theta=-58.5
+
+#Slope for sigmoid in mV
+!sigma=6
+
+#Scaling parameter for sigmoidal mapping(dimensionless)
+!C=(pi/sqrt(3))
+
+#PSC rise time in ms^-1
+!gamma_t=70E-3
+!gamma_r=100E-3
+
+#Connectivities
+!N_tr=3
+!N_rt=5
+!N_rr=25
+
+#Calcium dynamic
+!alpha_Ca=-51.8E-6
+!tau_Ca=10
+!Ca_0=2.4E-4
+
+#Parameters of h-current
+!k1= 2.5E7
+!k2= 4E-4
+!k3= 1E-1
+!k4= 1E-3
+!n_P= 4
+!g_inc= 2
+
+#Functions
+#Firing rates
+Qt(Vt)=Q_max/(1+exp(-C*(Vt-theta)/sigma))
+Qr(Vr)=Q_max/(1+exp(-C*(Vr-theta)/sigma))
+
+#Activation functions
+m_inf_T_t(Vt)=1/(1+exp(-(Vt+59)/6.2))
+m_inf_T_r(Vr)=1/(1+exp(-(Vr+52)/7.4))
+h_inf_T_t(Vt)=1/(1+exp((Vt+81)/4))
+h_inf_T_r(Vr)=1/(1+exp((Vr+80)/5))
+
+tau_h_T_t(Vt)=(30.8+(211.4+exp((Vt+115.2)/5))/(1+exp((Vt+86)/3.2)))/3.7371928
+tau_h_T_r(Vr)=(85+1/(exp((Vr+48)/4)+exp(-(Vr+407)/50)))/3.7371928
+
+m_inf_h(Vt)=1/(1+exp((Vt+75)/5.5))
+tau_m_h(Vt)=(20+1000/(exp((Vt+71.5)/14.2)+exp(-(Vt+89)/11.6)))
+
+#Synaptic currents
+I_tr(Vr,s_er)=s_er*(Vr-E_AMPA)
+I_rt(Vt,s_rt)=s_rt*(Vt-E_GABA)
+I_rr(Vr,s_rr)=s_rr*(Vr-E_GABA)
+
+#Leak currents
+I_L_t(Vt)=g_L*(Vt-E_L)
+I_L_r(Vr)=g_L*(Vr-E_L)
+I_LK_t(Vt)=g_LK*(Vt-E_K)
+I_LK_r(Vr)=g_LK*(Vr-E_K)
+
+#Calicum current
+I_T_t(Vt,h_T_t)=g_T_t*m_inf_T_t(Vt)*m_inf_T_t(Vt)*h_T_t*(Vt-E_Ca)
+I_T_r(Vr,h_T_r)=g_T_r*m_inf_T_r(Vr)*m_inf_T_r(Vr)*h_T_r*(Vr-E_Ca)
+
+#h-current
+I_h(Vt,m_h,m_h2)=g_h*(m_h+g_inc*m_h2)*(Vt-E_h)
+
+# Protein binding
+P_h(Ca)=(k1*Ca*Ca*Ca*Ca/(k1*Ca*Ca*Ca*Ca+k2))
+
+#System equation
+Vt'=-(I_L_t(Vt)+I_rt(Vt,s_rt))/tau-(I_LK_t(Vt)+I_T_t(Vt,h_T_t)+I_h(Vt,m_h,m_h2))
+Vr'=-(I_L_r(Vr)+I_tr(Vr,s_er)+I_rr(Vr,s_rr))/tau-(I_LK_r(Vr)+I_T_r(Vr,h_T_r))
+Ca'=(alpha_Ca*I_T_t(Vt,h_T_t)-(Ca-Ca_0)/tau_Ca)
+h_T_t'=(h_inf_T_t(Vt)-h_T_t)/tau_h_T_t(Vt)
+h_T_r'=(h_inf_T_r(Vr)-h_T_r)/tau_h_T_r(Vr)
+m_h'=((m_inf_h(Vt)*(1-m_h2)-m_h)/tau_m_h(Vt)-k3*P_h(Ca)*m_h+k4*m_h2)
+m_h2'=(k3*P_h(Ca)*m_h-k4*m_h2)
+s_er'=x_er
+s_rt'=x_rt
+s_rr'=x_rr
+x_er'=gamma_t*gamma_t*(N_tr*Qt(Vt)-s_er)-2*gamma_t*x_er
+x_rt'=gamma_r*gamma_r*(N_rt*Qr(Vr)-s_rt)-2*gamma_r*x_rt
+x_rr'=gamma_r*gamma_r*(N_rr*Qr(Vr)-s_rr)-2*gamma_r*x_rr
+
+# define the initial states
+init Vt=-68,Vr=-68,Ca=2.4E-4
+
+# simulation settings
+@ maxstor=1000000,total=3000,dt=0.1,meth=rungekutta,njmp=10,bound=200
+@ YLO=-80,YHI=0,XLO=0,XHI=3000
+@ nmax=80000,npr=0,ntst=100,epsl=1e-9,epss=1e-7,epsu=1e-7,parmin=0,parmax=0.1,
+@ autoxmin=0,autoxmax=0.1,autoymin=-100,autoymax=0, ds=1e-3,dsmin=1e-5,dsmax=0.005
+done

Figures/Tools/license_panel.txt

@@ -0,0 +1,24 @@
+Copyright (c) 2015, Ben Mitch
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in
+      the documentation and/or other materials provided with the distribution
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.