Homec4science
Diffusion Slepians (master) Browse

Slepi/Code/Utilities/spm12/spm_hdm_priors.m
8a186cd82567master

Last Change
Raw File

spm_hdm_priors.m
View Options

function [pE,pC] = spm_hdm_priors(m,h)
% returns priors for a hemodynamic dynamic causal model
% FORMAT [pE,pC] = spm_hdm_priors(m,[h])
% m - number of inputs
% h - number of hemodynamic modes (default = 3)
%
% pE - prior expectations
% pC - prior covariances
%
% (5) biophysical parameters
% P(1) - signal decay d(ds/dt)/ds)
% P(2) - autoregulation d(ds/dt)/df)
% P(3) - transit time (t0)
% P(4) - exponent for Fout(v) (alpha)
% P(5) - resting oxygen extraction (E0)
% P(6) - ratio of intra- to extra- (epsilon)
% vascular components of the
% gradient echo signal
%
% plus (m) efficacy priors
% P(7) - ....
%
%___________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
% Karl Friston
% $Id: spm_hdm_priors.m 4579 2011-12-02 20:21:07Z karl $
% default: 3 hemodynamic [eigen]modes
%---------------------------------------------------------------------------
if nargin < 2
h = 3;
end
% biophysical parameters with prior expectation and
%---------------------------------------------------------------------------
pE = [ 0.65 0.41 0.98 0.32 0.34 ];
% covariance restricted to h modes (v) scaled by eigenvales (e) {see below)
%---------------------------------------------------------------------------
v = [
-0.0535 0.0095 -0.1117 -0.0040 -0.0026
-0.0604 -0.0319 0.0430 -0.0077 0.0026
0.1116 -0.0347 -0.2539 -0.0169 -0.0115
0.1985 0.1698 0.4984 -0.4493 0.4434
0.0029 0.2081 1.9582 -0.5209 -1.1634]';
e = [2.1225 1.2006 0.3519 0.0039 0.0012];
% set variance of minor modes to zero
%---------------------------------------------------------------------------
i = (h + 1):5;
e(i) = 0;
pC = v*diag(e)*v'/32;
% append scaling parameter for epsilon: prior variance = 1/32,
%---------------------------------------------------------------------------
pE = [pE -1];
pC = blkdiag(pC,1/128);
% append m efficacy priors
%---------------------------------------------------------------------------
pE = [pE(:); zeros(m,1)];
pC = blkdiag(pC,eye(m,m)*32);
return
% NOTES: sample covariances from Friston et al (2000)
%---------------------------------------------------------------------------
qC = [0.0150 0.0052 0.0283 0.0002 -0.0027
0.0052 0.0020 0.0104 0.0004 -0.0013
0.0283 0.0104 0.0568 0.0010 -0.0069
0.0002 0.0004 0.0010 0.0013 -0.0010
-0.0027 -0.0013 -0.0069 -0.0010 0.0024];
% NOTES: Reduce rank of prior covariances for computational expediancy
%---------------------------------------------------------------------------
% assume independent priors in parameter space
%---------------------------------------------------------------------------
qC = diag(diag(qC));
% model specification (single node DCM)
%---------------------------------------------------------------------------
M.f = 'spm_fx_hdm';
M.g = 'spm_gx_hdm';
M.x = [0 0 0 0]';
M.pE = [pE(1:6) 1];
M.m = 1;
M.n = 4;
M.l = 1;
M.N = 32;
M.dt = 1/2;
% compute partial derivatives w.r.t. hemodynamic parameters [J] dy(t)/dp
%---------------------------------------------------------------------------
P = M.pE;
p = length(P);
dp = 1e-6;
[k J] = spm_nlsi(M);
for i = 1:5
M.pE = P;
M.pE(i) = M.pE(i) + dp;
[k q] = spm_nlsi(M);
Jq(:,i) = (q - J)/dp;
end
% implied covariance of impulse response
%---------------------------------------------------------------------------
Cq = Jq*qC*Jq';
% reduce to h hemodynamic modes in measurement space
%---------------------------------------------------------------------------
[v e] = spm_svd(Cq);
e = diag(e);
v = pinv(Jq)*v;
qC = v*diag(e)*v';
% NOTES: graphics - eigenvalues of qC
%---------------------------------------------------------------------------
subplot(2,2,1)
bar(e)
xlabel('eigen mode')
title('eigenvalue')
set(gca,'XLim',[0 6])
axis square
grid on
% graphics - response differentials
%---------------------------------------------------------------------------
subplot(2,2,2)
plot([1:M.N]*M.dt,Jq*v(:,1),[1:M.N]*M.dt,Jq*v(:,2),'-.')
xlabel('PST {secs}')
title('hemodynamic modes')
axis square
grid on
c4science · Help