%%% This code produces coherent hemodynamics spectroscopy spectra in the frequency domain.
%%% The outputs, which are stored in variable F_cmro2, are frequency dependent phase differences Arg(D)-Arg(O) and Arg(O)-Arg(T) as well as amplitude ratios of D/O and O/T
%%% Specifically, F_cmro2 have the structure of [f phi A phi2 A2];
%%% f ... frequency in Hz
%%% phi ... Arg(D)-Arg(O) in radians
%%% A ... |D|/|O|
%%% phi2 ... Arg(O)-Arg(T) in radians
%%% A2 ... |O|/|T|
%%% F_cmro2 is based on assuming that the cerebral metabolic rate of oxygen stays constant(cmro2=0) - the equations are based on the Physiological Measurements article
%%% For reference, equations, and citations please use:
%%% Fantini S. Dynamic model for the tissue concentration and oxygen saturation of hemoglobin in relation to blood volume,
%%% flow velocity, and oxygen consumption: Implications for functional neuroimaging and coherent hemodynamics spectroscopy (CHS).
%%% Neuroimage 2014; 85: 202-221.
%%% Fantini S. A new hemodynamic model shows that temporal perturbations of cerebral blood flow and metabolic rate of oxygen
%%% cannot be measured individually using functional near-infrared spectroscopy.
%%% Physiol Meas 2014; 35(1): N1-N9.
%%% Kainerstorfer JM. et al. Practical steps for applying a new dynamic model to near-infrared spectroscopy measurements of hemodynamic oscillations and transient changes: implications for cerebrovascular and functional brain studies.
%%% Acad. Radiol. 21(2), 185-196 (2014).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Parameters which can be changed:
%%% Parameters of the model (frequency independent):
Tc = 1.1; %%% Capillary transit time (s) - range: 0 - 2s
Tv = 5.6; %%% Venous transit time (s) - range: 0 - 10s
fc = 0.03;%0.03; %%% Autoregulation cut-off frequency (Hz) - range: 0 - 0.15Hz
k = 3; %%% Inverse of the modified Grubbs exponent - typical range: 2 - 5
CBV_a = 0.025*0.2; %0.005;%0.001; %%% Ratio of arterial blood volume to tissue volume - CBV0_a
CBV_c = 0.025*0.6; %0.015;%0.011; %%% Ratio of capillary blood volume to tissue volume - CBV0_c
CBV_v = 0.025*0.2; %0.005;%0.011; %%% Ratio of venous blood volume to tissue volume - CBV0_v
Far = 0.8; %%% Fahraeus factor
CBV=(CBV_a+CBV_c+CBV_v); %%% Total blood volume
ctHb = 2300; %%% Hemoglobin concentration in blood, expressed in units of micromoles per unit blood volume
%%%%%% (CBV_a+CBV_c+CBV_v)*ctHb equals the total hemoglobin content in tissue, expressed in uM
va = 0.02*exp(1i*0); %%% Arterial volume phasor ... 0.02 means 2% change in volume oscillations
vc = 0; %%% Capillary volume phasor. No capillary recruitment, hence no volume change in the capillaries.
vv = 0.02*exp(1i*0); %%% Venous volume phasor - typically the same as arterial volume change
vt = va*(CBV_a/CBV)+vv*(CBV_v/CBV); %%% cbv(w) with weighted volume contributions from the other compartments
Sa = 0.98; %%% Arterial saturation
alpha = 0.8; %%% Oxygen diffusion rate (1/sec)
fmax = 0.2; %%% Maximum frequency to be evaluated here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Sc = Sa*(1/Tc/alpha)*(1-exp(-alpha*Tc)); %%% Spatial average capillary saturation
Sv = Sa*exp(-alpha*Tc); %%% Venous saturation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Frequency dependent parameters:
idx=1;
for f=0.001:0.001:fmax; %%% Frequency vector for plotting the CHS spectra
H_RC_HP = f./(f-1i*fc); %%% Frequency response of RC high pass filter
cbf = k*H_RC_HP*(vt); %%% Blood flow phasor
flow(idx)=cbf;
H_RC_LP = 1./sqrt(1+(2*pi*f*Tc/exp(1)).^2).*exp(-1i*atan(2*pi*f*Tc/exp(1))); %%% frequency response of RC Low pass filter (capillary)
H_G_LP = exp(-0.5*log(2)*(1.765*(Tc+Tv)*f).^2).*exp(-1i*pi*(Tc+Tv)*f); %%% frequency response of Gaussian Low pass filter (venous)
%%% Absolute hemoglobin concentration oscillation:
%%% Physiological Measurements equation - what is used here is the relationship of oxy and deoxy with respect to the cerebral metabolic rate of oxygen (cmro2)
%%% The assumption is that the cerebral metabolic rate of oxygen stays constant (cmro2=0)
oxy = ctHb*((CBV_a*Sa*va+CBV_c*Far*Sc*vc+CBV_v*Sv*vv)*ones(length(f),1)+((Sc-Sv)*(Sc/Sv)*Far*CBV_c*H_RC_LP+(Sa-Sv)*CBV_v*H_G_LP).*cbf);
deoxy = ctHb*((CBV_a*(1-Sa)*va+CBV_c*Far*(1-Sc)*vc+CBV_v*(1-Sv)*vv)*ones(length(f),1)-((Sc-Sv)*(Sc/Sv)*Far*CBV_c*H_RC_LP+(Sa-Sv)*CBV_v*H_G_LP).*cbf);
total = oxy+deoxy;
A = abs(deoxy./oxy); %%% Ratio of Oxy and Deoxy hemoglobin
phi=unwrap(angle(deoxy./oxy)); %%% Phase difference between Oxy and Deoxy
phi(find(phi>0)) = phi(find(phi>0))-2*pi; %%% by definition the phase difference between Oxy and Deoxy is negative for us
A2 = abs(oxy./total); %%% Ratio of Oxy and Total hemoglobin
phi2=unwrap(angle(oxy./total)); %%% Phase difference between Oxy and Total hemoglobin
F_cmro2(idx,:) = [f phi A phi2 A2];
idx=idx+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2);
subplot(2,2,1); hold all; plot(F_cmro2(:,1),F_cmro2(:,4)*180./pi,'LineWidth',2); xlim([0 fmax]); ylabel('Arg(O)-Arg(T) (degrees)'); xlabel('Frequency (Hz)');
subplot(2,2,2); hold all; plot(F_cmro2(:,1),F_cmro2(:,5),'LineWidth',2); xlim([0 fmax]); ylabel('|O|/|T|'); xlabel('Frequency (Hz)');
subplot(2,2,3); hold all; plot(F_cmro2(:,1),F_cmro2(:,2)*180./pi,'LineWidth',2); xlim([0 fmax]); ylabel('Arg(D)-Arg(O) (degrees)'); xlabel('Frequency (Hz)');
subplot(2,2,4); hold all; plot(F_cmro2(:,1),F_cmro2(:,3),'LineWidth',2); xlim([0 fmax]); ylabel('|D|/|O|'); xlabel('Frequency (Hz)');