BrainStates

BrainStates - Determine brain state using LFP, EMG and movement.

  USAGE

    [exploration,sws,rem] = BrainStates(s,t,f,q,emg,<options>)

    s              spectrogram
    t              time bins for spectrogram
    f              frequency bins for spectrogram
    q              'quiescence' vector, obtained using <a href="matlab:help QuietPeriods">QuietPeriods</a> (see NOTE)
    emg            electromyogram (pass [] if missing)
    <options>      optional list of property-value pairs (see table below)

    =========================================================================
     Properties    Values
    -------------------------------------------------------------------------
     'nClusters'   number of clusters for K-means clustering (default = 2)
     'method'      see below (default = 'hippocampus')
     'nComponents' number of principal components (for 'pca' method only)
                   (default = automatically computed to account for 85% of
                   the variance)
     'show'        plot K-means clustering in feature space ('kmeans'),
                   clusters on spectrogram ('clusters') or both ('all')
                   (default = 'none')
    =========================================================================

  METHODS

    Option 'method' can take several possible values:

     'pca'                     PCA on the spectrogram
     'hippocampus', 'direct'   theta / delta ratio
     'cortex', 'ratios'        ratios described in Gervasoni et al. (2004)
     'amygdala'                gamma / low ratio

  OUTPUT

    exploration    at each time t, whether the rat was exploring
    sws            at each time t, whether the rat was in slow wave sleep
    rem            at each time t, whether the rat was in REM sleep

  NOTE

    Quiescence q is a time series of boolean values specifying for each
    timestamp whether the animal was quiet (as determined by <a href="matlab:help QuietPeriods">QuietPeriods</a>).

  SEE

    See also QuietPeriods, SpectrogramBands.

0001 function [exploration,sws,rem] = BrainStates(s,t,f,q,emg,varargin)
0002 
0003 %BrainStates - Determine brain state using LFP, EMG and movement.
0004 %
0005 %  USAGE
0006 %
0007 %    [exploration,sws,rem] = BrainStates(s,t,f,q,emg,<options>)
0008 %
0009 %    s              spectrogram
0010 %    t              time bins for spectrogram
0011 %    f              frequency bins for spectrogram
0012 %    q              'quiescence' vector, obtained using <a href="matlab:help QuietPeriods">QuietPeriods</a> (see NOTE)
0013 %    emg            electromyogram (pass [] if missing)
0014 %    <options>      optional list of property-value pairs (see table below)
0015 %
0016 %    =========================================================================
0017 %     Properties    Values
0018 %    -------------------------------------------------------------------------
0019 %     'nClusters'   number of clusters for K-means clustering (default = 2)
0020 %     'method'      see below (default = 'hippocampus')
0021 %     'nComponents' number of principal components (for 'pca' method only)
0022 %                   (default = automatically computed to account for 85% of
0023 %                   the variance)
0024 %     'show'        plot K-means clustering in feature space ('kmeans'),
0025 %                   clusters on spectrogram ('clusters') or both ('all')
0026 %                   (default = 'none')
0027 %    =========================================================================
0028 %
0029 %  METHODS
0030 %
0031 %    Option 'method' can take several possible values:
0032 %
0033 %     'pca'                     PCA on the spectrogram
0034 %     'hippocampus', 'direct'   theta / delta ratio
0035 %     'cortex', 'ratios'        ratios described in Gervasoni et al. (2004)
0036 %     'amygdala'                gamma / low ratio
0037 %
0038 %  OUTPUT
0039 %
0040 %    exploration    at each time t, whether the rat was exploring
0041 %    sws            at each time t, whether the rat was in slow wave sleep
0042 %    rem            at each time t, whether the rat was in REM sleep
0043 %
0044 %  NOTE
0045 %
0046 %    Quiescence q is a time series of boolean values specifying for each
0047 %    timestamp whether the animal was quiet (as determined by <a href="matlab:help QuietPeriods">QuietPeriods</a>).
0048 %
0049 %  SEE
0050 %
0051 %    See also QuietPeriods, SpectrogramBands.
0052 
0053 % Copyright (C) 2008-2016 by Michaël Zugaro, Gabrielle Girardeau
0054 %
0055 % This program is free software; you can redistribute it and/or modify
0056 % it under the terms of the GNU General Public License as published by
0057 % the Free Software Foundation; either version 3 of the License, or
0058 % (at your option) any later version.
0059 
0060 % Defaults
0061 nClusters = 2;
0062 window = 5*1250;
0063 show = 'none';
0064 nComponents = 0;
0065 method = 'hippocampus';
0066 
0067 % Check number of parameters
0068 if nargin < 5,
0069     error('Incorrect number of parameters (type ''help <a href="matlab:help BrainStates">BrainStates</a>'' for details).');
0070 end
0071 
0072 if mod(length(varargin),2) ~= 0,
0073   error('Incorrect number of parameters (type ''help <a href="matlab:help BrainStates">BrainStates</a>'' for details).');
0074 end
0075 
0076 % Parse options
0077 for i = 1:2:length(varargin),
0078     if ~ischar(varargin{i}),
0079         error(['Parameter ' num2str(i) ' is not a property (type ''help <a href="matlab:help BrainStates">BrainStates</a>'' for details).']);
0080     end
0081     switch(lower(varargin{i})),
0082         case 'nclusters',
0083             nClusters = varargin{i+1};
0084             if ~isiscalar(nClusters,'>0'),
0085                 error('Incorrect value for property ''nClusters'' (type ''help <a href="matlab:help BrainStates">BrainStates</a>'' for details).');
0086             end
0087         case 'method',
0088             method = lower(varargin{i+1});
0089             % 'direct' and 'ratios' are for backward compatibility (deprecated)
0090             if ~isastring(method,'pca','hippocampus','cortex','amygdala','direct','ratios'),
0091                 error('Incorrect value for property ''method'' (type ''help <a href="matlab:help BrainStates">BrainStates</a>'' for details).');
0092             end
0093         case 'ncomponents',
0094             nComponents = varargin{i+1};
0095             if ~isiscalar(nComponents,'>0'),
0096                 error('Incorrect value for property ''nComponents'' (type ''help <a href="matlab:help BrainStates">BrainStates</a>'' for details).');
0097             end
0098         case 'show',
0099             show = lower(varargin{i+1});
0100             if ~isastring(show,'kmeans','clusters','all','none'),
0101                 error('Incorrect value for property ''show'' (type ''help <a href="matlab:help BrainStates">BrainStates</a>'' for details).');
0102             end
0103         otherwise,
0104             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help BrainStates">BrainStates</a>'' for details).']);
0105     end
0106 end
0107 
0108 % Compute and interpolate EMG and quiescence at spectrogram timestamps
0109 emg0 = [];
0110 if ~isempty(emg),
0111     emg0 = Interpolate(emg,t,'trim','off');
0112     emg0 = Smooth(abs(emg0(:,2)),5);
0113 end
0114 q0 = [];
0115 if all(q(:,2)==0|q(:,2)==1),
0116     q0 = Interpolate(double(q),t,'trim','off');
0117     q0 = round(q0(:,2));
0118 else
0119     q0 = Interpolate(q,t,'trim','off');
0120     q0 = q0(:,2);
0121 end
0122 
0123 % No movement info => undetermined
0124 usable = ~isnan(q0);
0125 
0126 % Exploration corresponds to movement periods
0127 exploration = zeros(size(t));
0128 exploration(usable) = ~q0(usable);
0129 
0130 % Sleep/rest corresponds to quiescence
0131 sleep = zeros(size(t));
0132 sleep(usable) = q0(usable);
0133 sleep = logical(sleep);
0134 
0135 % Get theta/delta ratio
0136 bands = SpectrogramBands(s,f);
0137 
0138 % Determine features for automatic data clustering
0139 switch(method),
0140     case 'pca',
0141         % Compute PCA on spectrogram and reduce dimensionality
0142         S = s(:,sleep)';
0143         S(:,f>30) = 0;
0144         [eigenvectors,projected,lambda] = princomp(S,'econ');
0145         if nComponents == 0,
0146             nComponents = find(cumsum(lambda)/sum(lambda)>0.85);
0147             nComponents = nComponents(1);
0148         end
0149         eigenvectors = eigenvectors(:,1:nComponents);
0150         lambda = lambda(1:nComponents);
0151         projected = projected(:,1:nComponents);
0152         features = projected;
0153     case {'hippocampus','direct'},
0154         % Use theta/delta ratio
0155         features = bands.ratios.hippocampus(sleep);
0156     case {'cortex','ratios'},
0157         % Use heuristic ratios of Gervasoni et al. (2004)
0158         features = bands.ratios.cortex(sleep,:);
0159     case 'amygdala',
0160         % Use gamma / low ratio
0161         features = bands.ratios.amygdala(sleep);
0162 end
0163 
0164 % Add EMG to feature list
0165 if ~isempty(emg0),
0166     features = [features emg0(sleep)];
0167 end
0168 
0169 % Normalize columns (we need this because K-means clustering uses Euclidean
0170 % distances)
0171 features = zscore(features);
0172 
0173 % Cluster (K-means)
0174 cluster = kmeans(features,nClusters);
0175 
0176 % Identify clusters
0177 sws = logical(zeros(size(t)));
0178 rem = logical(zeros(size(t)));
0179 % 1) Measure mean ratio for each state
0180 switch(method),
0181     case 'amygdala',
0182         ratio = bands.ratios.amygdala(sleep);
0183     otherwise,
0184         ratio = bands.ratios.hippocampus(sleep);
0185 end
0186 for i = 1:nClusters,
0187     meanRatio(i) = mean(ratio(cluster==i));
0188 end
0189 % 2) REM is the quiet state with highest mean ratio
0190 [~,i] = sort(meanRatio(:),1,'descend');
0191 h = i(1);
0192 rem(sleep) = cluster==h;
0193 highest = meanRatio(h);
0194 % 3) SWS is the quiet state with lowest mean ratio
0195 [~,i] = sort(meanRatio(:),1,'ascend');
0196 l = i(1);
0197 sws(sleep) = cluster==l;
0198 lowest = meanRatio(l);
0199 % 4) REM must have a ratio at least twice as high as SWS
0200 if highest < 2*lowest,
0201     sws(sleep) = 1;
0202     rem(sleep) = 0;
0203 end
0204 
0205 % Show K-means clustering in feature space
0206 if strcmp(show,'kmeans') | strcmp(show,'all'),
0207     figure;
0208     % Determine number N of subplots (depends on number of principal components kept; clip at 15)
0209     n = size(features,2);
0210     N = min([n*(n-1)/2 21]);
0211     if N == 0, N = 1; end
0212     m0 = round(sqrt(N));
0213     n0 = ceil(N/m0);
0214     colors = jet(nComponents);
0215     if strcmp(method,'pca'),
0216         % Plot principal components
0217         subplot(m0+1,1,1);hold on;
0218         l = 'h = legend(';
0219         for i = 1:nComponents,
0220             plot(f,eigenvectors(:,i),'color',colors(i,:));
0221             l = [l '''' int2str(i) ''','];
0222         end
0223         l = [l(1:end-1) ',''location'',''northoutside'',''orientation'',''horizontal'');'];
0224         eval(l);
0225         set(h,'color',[1 1 1]*.7);
0226         shift = 1;
0227     else
0228         shift = 0;
0229     end
0230     % List pairs of principal component IDs to display in successive 2D projections
0231     colors = jet(nClusters);
0232     if N > 1,
0233         feature1 = [];
0234         feature2 = [];
0235         i = 1;
0236         while length(feature1) < N,
0237             feature1 = [feature1 repmat(i,1,n-i)];
0238             feature2 = [feature2 i+1:n];
0239             i = i+1;
0240         end
0241         % Plot 2D projections of data in PC space, colored by cluster ID
0242         for i = 1:N,
0243             subplot(m0+shift,n0,i+shift*n0);hold on;
0244             for j = 1:size(features,1),
0245                 plot(features(j,feature1(i)),features(j,feature2(i)),'.','color',colors(cluster(j),:));
0246             end
0247             xlabel([int2str(feature1(i)) ' vs ' int2str(feature2(i))]);
0248         end
0249     else
0250         subplot(shift+1,1,shift+1);hold on;
0251         r = rand(size(features));
0252         for j = 1:size(features,1),
0253             plot(r(j),features(j),'.','color',colors(cluster(j),:));
0254         end
0255     end
0256 end
0257 
0258 % Show clusters on spectrogram
0259 if strcmp(show,'clusters') | strcmp(show,'all'),
0260     figure;
0261     % 1) Plot spectrogram and ratio
0262     subplot(2,1,1);hold on;
0263     PlotColorMap(log(s),1,'cutoffs',[0 12],'x',t,'y',f);
0264     plot(t(sleep),ratio,'k');
0265     l = ['h = legend(''' method ' ratio'',''q'''];
0266     ylim([0 30]);
0267     yLim = ylim;
0268     dy = yLim(2) - yLim(1);
0269     y0 = yLim(1);
0270     % ... quiescence...
0271     plot(t,5*ZeroToOne(q0)+y0+0.5*dy,'b');
0272     % ... and EMG
0273     if ~isempty(emg0),
0274         plot(t,ZeroToOne(emg0)*5+y0+0.65*dy,'r');
0275         l = [l ',''emg'''];
0276     end
0277     l = [l ');'];
0278     eval(l);
0279     set(h,'color',[1 1 1]);
0280     xlim([t(1) t(end)]);
0281     % 2) Plot spectrogram and clusters
0282     subplot(2,1,2);hold on;
0283     PlotColorMap(log(s),1,'cutoffs',[0 12],'x',t,'y',f);
0284     ylim([0 30]);
0285     l = 'h = legend(';
0286     yLim = ylim;
0287     dy = yLim(2) - yLim(1);
0288     y0 = yLim(1);
0289     % Plot clustering output
0290     c = nan(size(t));
0291     c(sleep) = cluster-1;
0292     plot(t,5*ZeroToOne(c)+y0+0.35*dy,'k');
0293     l = [l '''cluster'',''location'',''north'',''orientation'',''horizontal'');'];
0294     eval(l);
0295     % Show SWS/REM
0296     if ~isempty(sws),
0297         z = nan(size(t));
0298         z(sws) = t(sws);
0299         plot(z,15*~isnan(z),'b','linewidth',5);
0300     end
0301     if ~isempty(rem),
0302         z = nan(size(t));
0303         z(rem) = t(rem);
0304         plot(z,15*~isnan(z),'r','linewidth',5);
0305     end
0306     xlim([t(1) t(end)]);
0307     set(h,'color',[1 1 1]);
0308 end