%% Lesson 13: Event-Related fMRI and 'Deconvolution'
% 
% Since the expected fMRI response is the convolution of the impulse
% response with the neuronal response, it's possible work backwards and
% estimate the hemodynamic response that best predicts a measured fMRI
% signal.  This is essentially undoing the convolution process, or
% 'deconvolution'.
%
% Let 'h' be be the vector representing the hemodynamic response function,
% and 's' be the vector of zeros and ones containing the time-points when
% 'events' occured.  A linear system predicts that the fMRI response, y, is
% the convolution of the h and s.  Here's an example:

dt = 1;  %step size (seconds)
maxt = 15;  %ending time (seconds)
th = 0:dt:(maxt-dt);
n = length(th);  %n is the length of the hdr.

k = 1; %seconds
nCascades = 3;
delay = 2; %seconds
h = Gamma(nCascades,k,th-delay)';  %make it a column vector

figure(1)
clf
plot(th,h);
xlabel('Time (s)');

% Here's the event sequence:
maxt = 255;  %ending time (seconds)
t = 0:dt:(maxt-dt);
m = length(t);  %m is the length of the event sequence

prob = .5;  %probability of an event for each time point.
s = floor(rand(1,m)+prob)';
y = conv(s,h)*dt;
y = y(1:m);

%% Deconvolution
%
% Suppose, now, that we only know the response, y, and our event sequence,
% s.  How can we reconstruct the hemodynamic response?  Specifically, if
% the hdr has, say, 15 time points, then we need to find the best 15
% numbers that will predict our fMRI response. A brute-force method would
% be to use a search algorithm to find the best 15 numbers that minimizes
% the error between the predicted and the measured fMRI response.
% Fortunately, there's an easier and more efficient way using linear
% algebra.
%
% The trick is to think of this is a giant linear algebra problem with 15
% unknowns and n knowns, where n is the length of the fMRI response.  We
% can see that this is a linear algebra problem by noting that convolution
% can be conducted through matrix multiplication (indeed, this is how
% matlab does it).  This is done by transforming our stimulus sequence, s,
% into a 'design matrix' X that, when multipied by the hdr, performs the
% convolution.
%
% Remember that convolution is the process of multiplying shifted copies of
% the hdr with the input and adding.  Equivalently (since convolution is
% commutative) this can be done with shifted copies of the input instead.
%
% The columns of the design matrix, X, are generated by successive shifts of
% the stimulus, like this:

X = zeros(m,n);
temp = s;
for i=1:n
    X(:,i) = temp;
    temp = [0;temp(1:end-1)];
end

figure(1)
clf
imagesc(X)
axis equal
axis off
colormap(gray)

%%
% convolution of s and h is done by multiplying the design matrix X and h
% (as a column vector):

r = X*h;

%%
% we can compare this to matlab's convolution by taking the norm of the
% difference between the two vectors.  A small number means they're the
% same (the norm of a vector is the square root of the sum of squares).

rconv = conv(s,h); rconv = rconv(1:m);
norm(r-rconv)

%%
% Here's where the linear algebra comes in.  if r = X*h, then if this were
% normal algebra we could solve for h by dividing r by X (h = r/X).  When r
% and h are vectors and X is a matrix, this 'division' is done by finding
% the least-squared solution with the 'psueudo-inverse'.  How this works is
% beyond our scope (but elegant). Details can be found in most engineering
% and linear algebra textbooks.

PX = inv(X'*X)*X';

%%
% Matlab has its own function that does this called 'pinv'

PX = pinv(X);

hest = PX*r;

%%
% We'll use 'norm' again to compare hest to h:

norm(hest-h)

%%
% What did we just do?  To recap, we made up a hemodynamic response
% function, h, convolved it with a random stimulus sequence, s, to get an
% fMRI response, r.  We then did the whole thing in reverse and
% reconstructed h from r and s by deconvolving via the pseudo-inverse.
%
% Real fMRI data is noisy.  Let's add independently and identically
% distributed ('iid') noise to the response and see how it affects our
% ability to reconstruct the hemodynamic response

noiseSD = .2;
rnoise = r+noiseSD*randn(m,1);

figure(1)
clf

plot(t,r,'b-');
hold on
plot(t,rnoise,'bo');
legend({'actual response','response+noise'});

hest = PX*rnoise;

figure(2)
clf
plot(th,h,'b-');
hold on
plot(th,hest,'bo','MarkerFaceColor','k')
xlabel('Time');
legend({'actual hdr','reconstructed hdr'});

norm(hest-h)

%%
% You might be surprised how robust our estimate of the hdr is considering
% how much noise we added to the response.  This is because we have so many
% knowns to help constrain the number of unknowns.  As you can guess, the
% longer the event sequence, the better we will be at estimating the true
% hdr.
%
% Go ahead and play with the magnitude of the noise (noiseSD) and see how
% it affects the estimate of the hdr.
%
% It turns out that the average value of norm(h-hest) grows in proportion
% to the standard deviation of the added noise (or equivalently the sums of
% squared error grows in proportion to the variance of the noise)  We can
% demonstrate this with a simulation:

noiseVar = linspace(0,1,21);  %list of noise variances
nReps = 5000;  %number of repetitions per noise variance level

err=  zeros(nReps,length(noiseVar));
for i=1:length(noiseVar)
    for j=1:nReps
        rnoise = r+sqrt(noiseVar(i))*randn(m,1);
        hest = PX*rnoise;
        err(j,i) = norm(hest-h)^2;
    end
end

meanErr=  mean(err);

figure(1)
clf
plot(noiseVar,meanErr,'ro','MarkerFaceColor','r');
xlabel('Noise Variance');
ylabel('SSE');

%% Efficiency of a stimulus sequence
%
% The slope of this line tells us something about the way the ability to
% reconstruct the hdr falls apart with increasing noise.  A lower slope is
% better.  Different event sequences will have different slopes. Try it out
% yourself by running the program with a different kind of sequence.
%
% Since low slopes are good, the inverse of the slope has a special
% meaning, called the 'efficiency' of the event sequence.  Let's calculate
% this number from our simulated data:
%
%Fit a line to the data
p = polyfit(noiseVar,meanErr,1);
bestLine = polyval(p,noiseVar);
hold on
plot(noiseVar,bestLine,'k-');
legend({'Average of norm(h-hest)','Best-fitting line'});


%%
% It also turns out that there's a simple way of calculating the efficiency
% without a simulation.  Here's how:

E = 1/trace(inv(X'*X));

%%
% The trace of a matrix is the sum of its diagonals.  The derivation of
% this equation is, once again, beyond the scope of this lesson (and beyond
% the scope of my full understanding).  But let's just assume it's true.
% This number should be close to the inverse of the slope from the
% simulation:

title(sprintf('1/Slope = %5.3g, Efficiency = %5.3g',1/p(1),E));

%%
% So we now have an easy way to assess how well a given event sequence will
% be able to reconstruct the hemodynamic response (assuming that the
% 'hemodynamic coupling' process is linear, and that noise is additive,
% which can't be completely right).
%
% Here's something interesting: the calculation of efficiency doesn't
% depend on the shape of the hdr.  It's kind of an amazing fact if you
% think about it.  We could have chosen any crazy shape for the hdr and the
% ability to reconstruct it from noise data will be exactly the same.  This
% is a very convenient fact, too, since we don't know the actual shape of
% the hdr (it's what we're trying to find out).
%
%% m-sequences
%
% There is one kind of sequence that has the greatest efficiency - the
% m-sequence.  M-sequences (or maximum-length shift register sequences)
% were originally developed about 50 years ago and have been used in
% encryption, error-correcting codes - any time you want to generate a
% signal that is minimally corruptible by noise.  Generation of m-sequences
% is not intuitive - it involves prime numbers and modulo arithmetic.  But
% the result is a sequence that has a zero autocorrelation function.  I've
% provided a function 'mseq' that generates these sequences for you. 'mseq'
% takes in a minimum of two parameters.  The first is the number of event
% types (here we've only talked about zeros and ones, so the number is 2,
% since blank trials count as an event type).  The second is a 'power
% value' that determines the length of the sequence, which will be
% nTypes^(powerVal-1).
%
% A more thorough discussion of m-sequences applied to event-related fMRI
% can be found here:
%
% G.T. Buracas and G.,M. Boynton (2002), "Efficient Design of Event-Related
% fMRI Experiments Using M-Sequences" , NeuroImage 16, 801-813

powerVal = 8;
ms = mseq(2,powerVal);

X = zeros(m,n);
temp = ms;

%Calculate the efficiency:
for i=1:n
    X(:,i) = temp;
    temp = [0;temp(1:end-1)];
end
Emseq = 1/trace(inv(X'*X));

disp(sprintf('Efficiency of m-sequence: %5.3g',Emseq));

%%
% This number should be greater than our example random sequence (it had
% better be anyway).

%% Multiple event types
%
% Event-related fMRI often involves more than one type of stimulus in a
% given scan.  Stimuli may include, for example, a parametric range of
% stimulus values, like contrast.  Or they may be different cognitive or
% attentive tasks.  There's an easy way to incorporate multiple event types
% into your design and analysis.  The trick is to simply concatenate the
% separate design matrices horizontally (with zeros and ones in each
% matrix).  The resulting estimate of the hemodynamic response will be a
% single vector with the estimate of the response to each condition in
% succession.  We'll demonstrate this by creating four 'fake' hdrs,
% generating a response with convolution, adding noise, and reconstructing
% the hdrs.

amps = [1,2,3,4];  %amplitudes of the hdr for the four event types (e.g. for increasing stimulus contrast).

h = zeros(n,4);
for i=1:4
    delay = 2; %seconds
    h(:,i) = amps(i)*Gamma(nCascades,k,th-delay)';  %make it a column vector
end

figure(1)
clf
plot(th,h);
xlabel('Time (s)');

%%
% get an m-sequence with four event types - plus null events.  This vector
% will have values from 0 to 4 corresponding to what happens on each time
% step.

s = mseq(5,3);
m = length(s);
t = 0:dt:(m-1);
%%
% Generate a concatenated design matrix
X = [];
for j=1:4
    Xj = zeros(m,n);
    temp = s==j;

    for i=1:n
        Xj(:,i) = temp;
        temp = [0;temp(1:end-1)];
    end
    X = [X,Xj];
end

figure(1)
clf
imagesc(X)
axis equal
axis off
colormap(gray)

r = X*h(:);

%add noise
noiseSD = .25;
rNoise = r+noiseSD*randn(size(r));

hest = pinv(X)*rNoise;

%%
% The resulting hdr estimates are all in one long vector.  To show them in
% the form of the original hdr we can reshape this vector back into an nx4
% matrix:

hest = reshape(hest,n,4);

figure(1)
clf
plot(th,h);
hold on
plot(th,hest,'o');
xlabel('Time (s)');

%% Other regressors:
%
% In addition to the additive-type noise, fMRI signals are contaminated
% with low-frequency noise, often characterized as linear trend either up
% or down, along with a large y-intercept.  The y-intercept is just the
% mean of the image which doesn't carry any functional information (the
% mean is a T2-weighted image that contains structural information,
% however).  The linear trend is more troublesome - it can dominate the
% time-course of a voxel when viewed by eye, yet its causes aren't well
% understood.  So we just subtract it out.  Scary?  Yes.
%
% It's easy to remove these trends in the same step as the deconvolution
% process.  This is because each column in the design matrix 'X' represents
% the time-course to be explained by each free parameter.  The first column
% is for the first time-point in the hdr, for example.  To regress out a
% DC, then, we simply add a column of 1's to the end of the design matrix.
% The resulting last parameter will be the best-fitting parameter for the
% DC (which will be the mean). Similarly, if we add a ramp to X, the
% regressor will be the best fitting scale factor that makes that ramp fit
% the data.  We don't use these parameters, instead by adding them to the
% matrix their influences get factored, or regressed out.  
%
% First, let's see how a DC and trend screws up our estimation of the hdr
% before regressing them out:

rtrend = r+100 + linspace(0,10,m)';

figure(1)
clf
plot(t,rtrend);
xlabel('Time (s)');

%%
% Now THAT looks more like fMRI data.  Here's the estimated hdr:

hest = pinv(X)*rtrend;
hest = reshape(hest,n,4);

figure(1)
clf
plot(th,h);
hold on
plot(th,hest,'o');
xlabel('Time (s)');

%%
% What a mess.  Let's add a column of ones to X to deal with the DC:

X(:,n*4+1) = ones(m,1);
hest = pinv(X)*rtrend;

%%
% The DC regressor is the last element in hest.  It should be about 100
dc = hest(end);

hest = reshape(hest(1:n*4),n,4);

figure(1)
clf
plot(th,h);
hold on
plot(th,hest,'o');
xlabel('Time (s)');

%%
% Better, but we have to get rid of the linear trend, too

X(:,n*4+2) = linspace(0,1,m)';

hest = pinv(X)*rtrend;
lowFac = hest(end-1:end);
hest = reshape(hest(1:end-2),n,4);

figure(1)
clf
h1=  plot(th,h);
hold on
h2= plot(th,hest,'o');
xlabel('Time (s)');

%% 
% It's back.  You can see how important it is to get rid of these factors -
% the DC is essential because the deconvolution process only works if you
% either remove the DC via regression, or if you remove it ahead of time by
% subtracting out the mean.  
%

%% Estimating a voxel's sensitivity to a condition type
%
% The deconvolution method described above provides an estimate of the
% shape and amplitude of the hemodyanmic response function that generated a
% voxel's response.  How can we translate this into a single number that
% represents how sensitive a voxel is to a stimulus?  In other words, how
% do we collapse the 15 or so numbers that represent the estimated HDR into
% a single value representing amplitude?
%
% This is the same thing as finding the 'beta weight' in many fMRI data
% analysis software packages.
%
% Almost any solution will require some assumptions about the shape of the
% underlying hemodynamic response function.  For example, we could take the
% average of the hdr estimate between two time points, say between 4 and 8
% seconds.  This assumes that most of the information about the hemodyamic
% response lies within this time period.  
%
% Or, we could go further and assume a specific hemodynamic response
% function.  This could either be estimated empirically in a separate
% experiment, or more commonly, be chosen based on other work (like the
% Boynton 1996 paper).  
%
% It works like this.  Let's generate a fake, noisy fMRI data set using a single
% event type and an m-sequence stimulus:

% Stimulus

s = mseq(2,7);
m = length(s);  %2^7-1
t = 0:dt:dt*(m-1);

k = 1; %seconds
nCascades = 3;
delay = 2; %seconds
amp = 2;  %amplitude of hdr
h = amp*Gamma(nCascades,k,th-delay)';  %make it a column vector 

% event matrix:

X = zeros(m,n);
temp = s;
for i=1:n
    X(:,i) = temp;
    temp = [0;temp(1:end-1)];
end

y = X*h;
y = y+noiseSD*rand(size(y));

%%
% Originally we estimated the hdr using the pseudoinverse, which found the
% best scale parameters for each time point to minimize the sums of squared
% error between the predicted hdr and the actual time course.

hest = pinv(X)*y;


%% 
% Instead, if we multiply our event matrix by a guess at the hdr and use
% that instead of the event matrix: (here we'll use the parameters from the
% actual hdr that generated the fake data)

hGuess = Gamma(nCascades,k,th)';  %note, no scaling by 'amp'

XX = X*hGuess;
ampEst = pinv(XX)*y;

disp(sprintf('estimated amp: %5.2f, actual amp: %5.2f',ampEst,amp));

%% 
% The result is a single number that is our best estimate of the amplitude
% parameter.  Why?  It's because the matrix XX is just the convolution of
% the stimulus with our best (unscaled) estimate of the hdr, so it's our
% best (unscaled) guess at the fMRI time-course.  The 'pinv' process then
% simply finds the best scale factor for this time course to best match our
% data in a least squares sense.  


%% Beta weights for multiple event types
%
% The process of estimating 'beta weights' extend easily to multiple event
% types.  Let's go back to our earlier example of four event types and four
% scaled hdrs:

amps = [1,2,3,4];  %amplitudes of the hdr for the four event types (e.g. for increasing stimulus contrast).

h = zeros(n,4);
for i=1:4
    delay = 2; %seconds
    h(:,i) = amps(i)*Gamma(nCascades,k,th-delay)';  %make it a column vector
end

%%
% Generate a concatenated design matrix using an m-sequence

s = mseq(5,3);
m = length(s);
t = 0:dt:(m-1);

X = [];
for j=1:4
    Xj = zeros(m,n);
    temp = s==j;

    for i=1:n
        Xj(:,i) = temp;
        temp = [0;temp(1:end-1)];
    end
    X = [X,Xj];
end

%%
% Use the design matrix and our hdr's to generate a fake, noisy time-course
y = X*h(:);

%%
%add noise, dc and linear ramp
y = y+noiseSD*randn(size(y))+ rand(1) +randn(1)*linspace(0,1,m)';

%% 
% We can recover estimates of our four amplitudes by multiplying our design
% matrix a matrix with four columns, each containing a shifted, unscaled
% estimate of the hdr:

hMat = zeros(n*4,4);
id = 1:n;
for i=1:4
    hMat(id,i) = Gamma(nCascades,k,th-delay)';
    id = id+n;
end

XX = X*hMat;
XX = [XX,ones(m,1),linspace(0,1,m)'];

ampEst = pinv(XX)*y;

%% 
% The first four parameters of ampEst are our best estimates of the scale
% factors needed to scale each hdr so that when convolved with the event
% matrix gives us our best estimate of the fMRI time course.  The last two
% parameters are the dc and linear ramp factors that we can throw away.

figure(3)
clf
bar([ampEst(1:4),amps']);
legend({'Predicted','actual'},'Location','NorthWest')
xlabel('Event type');
ylabel('Amplitude (beta-weight)');

%%
% This works because multiplying the design matrix X by the four-column
% vector hMat is doing the same thing as convolving the stimulus sequence
% for each event type with the estimated hdr.  Each column of XX is
% therefore the estimated fMRI response for each isolated event type. Under
% linearity, the resulting fMRI time-course should be a linear combination
% of these four time-courses.  Deconvolving using pinv provides the four
% weights (or beta weights) that provide the best fit for the linear
% combination.  
%
% Often when estimating beta weights, something closer to a blocked-design
% is used instead of m-sequences.  For example, if there are four event
% types, we might present each stimulus for 20 seconds at a time in a
% shuffled order with a 10 second gap in between.  Each cycle through the
% four event types will take two minutes, and we'll repeat the whole thing
% twice for a four-minute scan.

dt =1;
maxt = 240;  %six minutes
t = (0:dt:(maxt-dt))';
m = length(t);

seq = shuffle(1:4);

s = zeros(size(t));
eventDur = 20; %seconds
eventGap = 10;
t0 = 0;
for i=1:2 % 4 cycles
    for j=1:4 % 4 event types
        s(t>=t0 & t<t0+eventDur) = seq(j);
        t0 = t0 + eventDur+eventGap;
    end
end

%%
% We'll display this design like some software packages do:

figure(1)
clf
subplot(2,1,1);
hold on
colorList = {[.7,.7,1],[.7,1,.7],[1,.7,.7],[.7,1,1]};
for i=1:(m-1)
    if s(i)
        patch([t(i),t(i+1),t(i+1),t(i)],[0,0,1,1],colorList{s(i)},'EdgeColor','none');
    end
end

set(gca,'XLim',[0,maxt]);
xlabel('Time (sec)');
ylabel('Event type');


%%
% We'll generate our event matrix like before:
X = []; 
for j=1:4
    Xj = zeros(m,n);
    temp = s==j;

    for i=1:n
        Xj(:,i) = temp;
        temp = [0;temp(1:end-1)];
    end
    X = [X,Xj];
end

figure(2)
clf
imagesc(X)
axis equal
colormap(gray)
axis off
    
%%
% Again, we'll generate a fake voxel response to this by multiplying our
% event matrix by h:

y = X*h(:);

% And add noise, dc and ramp:

y = y+.5*randn(size(y))+ 100 +10*linspace(0,1,m)';

figure(1)
subplot(2,1,2)
plot(t,y,'b-');
set(gca,'XLim',[0,maxt]);
xlabel('Time (sec)');
ylabel('Response');

%% 
% And just like before, we can estimate the four amplitude parameters
% through deconvolution with our event matrix mulitplied by our special
% four-column matrix containing our template hdr:

XX = X*hMat;

%%
% This matrix XX contains four columns, each containing the unscaled
% predicted response to each isolated stimulus condition.  We can draw
% these time-courses on top of our experimental design like the software
% packages do:

figure(1)
subplot(2,1,1)
hold on
plot(t,XX,'LineWidth',2);

%%
% Now we can add dc and linear trend columns to XX so that we have size
% 'regressor'.

XX = [XX,ones(m,1),linspace(0,1,m)'];

%%
% pinv(XX)*y will give us the six scale values for the best linear combination of
% the columns of XX to fit the data, y.

ampEst = pinv(XX)*y;

figure(3)
clf
bar([ampEst(1:4),amps']);
legend({'Predicted','actual'},'Location','NorthWest')
xlabel('Event type');
ylabel('Amplitude (beta-weight)');

%%
% Finally we can get our best estimate of the fMRI response by multiplying
% our design matrix by our estimated amplitudes (and dc and ramp values)

predY = XX*ampEst;

%%
% and overlay it on top of the original time-course
figure(1)
subplot(2,1,2)
hold on
plot(t,predY,'r-');

%%
% In this discussion of deconvolution, we assumed that the fMRI response is
% a linear transformation of the neuronal response with iid noise.  It's
% known that the noise in the signal has temporal correlations in it.  For
% a discussion on how to deal with that via deconvolution, see an excellent
% paper by Anders Dale:
%
% Dale, A. (1999) "Optimal Experimental Design for Event-Related fMRI",
% Human Brain Mapping 8:109–114

%% Exercises
%
% 1) The efficiency of a random sequence depends on the probability of an
% event occurring on a given trial (the variable 'prob' above).  Plot a
% graph of efficiency as a function of prob (between 0 and 1) to see what
% probability yields the most efficient sequence.
%
% 2) Generate a bunch of random 50% probability event sequences to see if
% you get anything close to the efficiency of an m-sequence.  Do you think
% this random procedure is worthwhile?  What are its advantages?
%
% 3) Removing the DC and linear trend is one way of removing
% 'low-frequency' factors in the fMRI signal.  Using what we learned from
% Lesson 12, how would you remove low-frequency stuff using the frequency
% domain?  Many fMRI software packages do this sort of thing.  
%
% 4) The optimal design matrix for reconstructing the hdr may be different
% from the design for estimating the amplitudes (or beta-weights) given a
% template HDR.  For example, an on-off blocked design is terrible for
% estimating the shape of the hdr, but optimal for estimating the
% amplitude.  Look at the ability to reconstruct hdrs and amplitudes by
% generating noisy data and re-estimating your parameters to see how
% efficient different designs (e.g. m-sequences vs. blocked designs) are
% for estimating hdrs vs. amplitudes.