Skip to content

Commit

Permalink
make the tutorial more accessible
Browse files Browse the repository at this point in the history
  • Loading branch information
@raltodo
raltodo committed Nov 25, 2024
1 parent 9d22944 commit f2a851d
Showing 1 changed file with 48 additions and 22 deletions.
70 changes: 48 additions & 22 deletions tutorials/pipelineReactivation.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,14 +14,18 @@
% I recommend using theta cycles during running for (a) and activity bursts
% during sws for (b), but other solutions are possible.


% This script first loads all the data that the tutorial will require. If
% your data is saved differently, skip the next section, at the end of
% which there is a list of all the variables the tutorial will need.

%% Load session data:

basepath = 'Z:\Data\Can\OML22\day6'; % we will be using this session for the sake of this tutorial
basename = basenameFromBasepath(basepath);
spikeStructure = importSpikes('basepath',basepath,'CellType',"Pyramidal Cell",'brainRegion','CA1'); % load CA1 pyramidal cells

% We will need a spikes matrix, which all the functions making up this
% tutorial will expect:
% We will need a spikes matrix, which all the functions making up this tutorial will expect:
try
spikes = Group(spikeStructure.times{:},'sort',true);
catch
Expand All @@ -35,8 +39,6 @@
spikes = sortrows(spikes); % sort spikes accoring to their time
end

%% Detect key intervals

% As stated in the note above, this tutorial uses theta cycles during running
% as the task period defining the patterns. We load them from here:
behavior = getStruct(basepath,'animal.behavior'); % load behavior structure
Expand Down Expand Up @@ -66,6 +68,9 @@
% Therefore make sure to define pre-task sleep and post-task sleep in a way
% that makes sense for your project.

ripples = getStruct(basepath,'ripples'); rippleTimestamps = ripples.timestamps(:,1);
activePeriods = behavior.run;

% Truncate preSleep and postSleep to one hour:
preSleep = TruncateIntervals(preSleep,3600);
postSleep = TruncateIntervals(postSleep,3600);
Expand All @@ -76,13 +81,12 @@
% happen to record for longer would appear to have poorer reactivation).

simple = true;

if simple,
if simple
% The simplest possible way to define our intervals would be the following:
% This is ALTERNATIVE 1:
taskIntervals = SplitIntervals(behavior.run,'pieceSize',0.02); % 20 ms bins in behavior
preIntervals = SplitIntervals(preSleep,'pieceSize',0.02); % 20 ms bins in pre-task sleep
postIntervals = SplitIntervals(postSleep,'pieceSize',0.02); % 20 ms bins in post-task sleep
binsActive = SplitIntervals(activePeriods,'pieceSize',0.02); % 20 ms bins in behavior
binsPre = SplitIntervals(preSleep,'pieceSize',0.02); % 20 ms bins in pre-task sleep
binsPost = SplitIntervals(postSleep,'pieceSize',0.02); % 20 ms bins in post-task sleep
% In general, I do recommend restricting the task bins to running periods,
% because if you include immobility, that also will end up including whatever
% is being reactivated during awake ripples, which may include remote replay
Expand All @@ -102,7 +106,7 @@
catch % if they don't exist, compute them
error(['No "thetacycles" for .' basepath '. Please run auto_theta_cycles']);
end
taskIntervals = Restrict(thetacycles.timestamps,behavior.run);
binsActive = Restrict(thetacycles.timestamps,behavior.run);

% Detect bursts of activity in NREM pre-task and post-task
% If you have already performed replay analysis, these are already detected:
Expand All @@ -129,10 +133,34 @@
inRunning = IntervalsIntersect(bursts(:,[1 3]),behavior.run);
bursts(inRunning,:);
end
preIntervals = Restrict(bursts(:,[1 3]), preSleep);
postIntervals = Restrict(bursts(:,[1 3]), postSleep);
binsPre = Restrict(bursts(:,[1 3]), preSleep);
binsPost = Restrict(bursts(:,[1 3]), postSleep);
end


%% At this stage, you are expected to have the following variables in your workspace:
%
% spikes: a matrix of [timestamp id] pairs, one row for each spike
% binsPre: a matrix of [start stop] timestamps of bins within
% pre-task sleep epochs. These can be e.g. 20-ms bins
% within all slow-wave sleep periods, or specifically
% events of interests (e.g. bursts of activity,
% ripples, etc.).
% binsActive: a matrix of [start stop] timestamps of bins within
% of active behavior. These can be e.g. 20-ms bins
% within active behavior, or specifically events of
% interest (e.g. exploration of an object, theta
% cycles, etc.)
% binsPost: a matrix of [start stop] timestamps of bins within
% post-task sleep epochs. These can be e.g. 20-ms bins
% within all slow-wave sleep periods, or specifically
% events of interests (e.g. bursts of activity,
% ripples, etc.).
% rippleTimestamps: a vector of timestamps of detected ripples. Feel free
% to use bursts or any other events of interest
%
% Feel free to load them from your own data before executing the rest of the tutorial.

%% Explained Variance (EV)

% This is the measure introduced by Kudrimoti et al. (1999):
Expand All @@ -147,7 +175,7 @@
% are reversed, yielding the Reverse Explained Variance. If we have more
% task-related activity in post-task sleep than pre-task sleep, we'd expect
% EV>REV.
[EV, REV] = ExplainedVariance(spikes, preIntervals, taskIntervals, postIntervals);
[EV, REV] = ExplainedVariance(spikes, binsPre, binsActive, binsPost);

% Plot the result:
figure(1); clf;
Expand All @@ -167,14 +195,14 @@
%% Cell pair correlations

% First, we compute the correlation matrix in behavior:
rTask = CorrInIntervals(spikes, taskIntervals);
rTask = CorrInIntervals(spikes, binsActive);

% We implement this measure as described by Giri et al (2019):
% The correlation matrices in pre- and post-task sleep are very similar.
% We can extract the part of the correlation matrix in post-task sleep
% that cannot be explained by the correlation matrix in pre-task sleep
% (residual correlation):
residuals = CorrResiduals(spikes, preIntervals, postIntervals);
residuals = CorrResiduals(spikes, binsPre, binsPost);
% Note that if the pre- and post- correlations were linearly correlated
% with slope=1 and intercept=0, the the residuals would simply be
% corrPost - corrPre*1 +0. However, in reality, the best fit might not be
Expand Down Expand Up @@ -204,7 +232,6 @@
title(['r= ' num2str(r) ', p=' num2str(p)]);
set(gca,'box','off','TickDir','out','fontsize',12);


% Here, we have used firing correlations in the task for the x-axis.
% Alternative approaches include using overlap between place fields,
% correlations between firing maps, etc.
Expand All @@ -219,20 +246,18 @@
% previous run may now be component 6). To make sure we can reproduce our
% results, we set the random number generator with a value (here, 0), which
% will result in the same output every time we execute the following line:
templates = ActivityTemplates(spikes,'bins',taskIntervals,'mode','ica');
templates = ActivityTemplates(spikes,'bins',binsActive,'mode','ica');
% Then, detect the (re)activation of each template/component over the whole session
binSize = 0.02; step = 0.001;
bins = Bins(0,MergePoints.timestamps(end),binSize,step);
bins = Bins(0,postSleep(end),binSize,step);
% We don't necessarity need the bins outside of pre-task sleep and post-task sleep
bins = Restrict(bins,[preSleep; postSleep]); % this saves memory when computing the strength
strength = ReactivationStrength(spikes,templates,'bins',bins);
% The first column of "strength" is timestamps and the following columns are the activation strength of each component for that timestamp

% Compute the activation strength around pre-task and post-task ripples:
% If you have not yet detected bursts, choosing the "simple" option above, detect them now:
ripples = getStruct(basepath,'ripples');
preRipples = Restrict(ripples.timestamps(:,1),preSleep);
postRipples = Restrict(ripples.timestamps(:,1),postSleep);
preRipples = Restrict(rippleTimestamps,preSleep);
postRipples = Restrict(rippleTimestamps,postSleep);
nComponents = size(templates,3);
nBins = 101; durations = [-1 1]*0.5;
[mPre,mPost] = deal(nan(nComponents,nBins)); % pre-allocate matrices
Expand All @@ -244,6 +269,7 @@
x = linspace(durations(1),durations(2),nBins);

%% Plot the result:

figure(3); clf;
set(gcf,'position',[300 200 1500 700]);
matrices = {mPre,mPost,mPost-mPre};
Expand Down

0 comments on commit f2a851d

Please sign in to comment.