function [BOLD,q,dCBV] = caoBalloonModel(dCMRO2,dt)

% Usage:   [BOLD,q,dCBV] = caoBalloonModel(dCMRO2,dt)
%
% Based on Buxton's (famous) balloon model, which predicts local
% BOLD signal strength as a function of CBF, CBV, and CMRO2.
% Equations taken from Obata et al. (2004) NeuroImage 21:144.
%   and Mildner et al. (2001), MRM46:819.
% Inputs:
% dCMRO2 is a 1 x n (row) vector showing the value of CMRO2 at 
% each time point; dt indicates temporal resolution of model.
%   e.g. dCMRO2 = [1 1.05 1.05 1.05 1.05 1 1 1 1 1 1 1 1 1 1 1 1 1]; 
%        dt = 2.0 for 8s stim with 
%
% Calculations:
% dCBF is calculated directly from dCMRO2 (model could be much
%   more flexible if CMRO2 were varied independently), multiplying
%   input normalized dCMRO2 by 3 to get dCBF, e.g. 5% increase in
%   CMRO2 at timepoint t:  dCMRO2(t) = 1.05 will correspond to
%   dCBF(t) = 1.30;
% dCBV, normalized change in CBV, is calculated from inflow/outflow
%   relationships reported in the literature:
%   dv/dt = 1/tau * [flow_in(t) - flow_out(v,t)] 
%   flow_out(v,t) = dCBV^(1/alpha) + tau_v*dv/dt 
% finally, q(t), the total deoxyhemoglobin concentration, is 
%   calculated from dCMRO2, dCBF, and flow_out
%
% Output:
% BOLD = V0*(A1*(1-q)-A2*(1-dCBV));
%   where V0 is resting state venous blood volume (2.0%)
%           A1 is a field-dependant constant (3.4 at 1.5T)
%           A2 is a field-dependant constant (1.0 at 1.5T)
%
% a hack in progress ... cao 1/24/05

% First interpolate input data to get temporal resolution of 0.01s.
timebase = dt*(0:(length(dCMRO2)-1));  % input time
dt = 0.1; % new temporal resolution
dCMRO2 = interp1(timebase,dCMRO2,0:dt:max(timebase)); % new CMRO2
timebase = dt*(0:(length(dCMRO2)-1)); % new time, for plotting

% Calculate CBF from CMRO2.  It will be seen later that the model
% allows room to vary oxygen extraction rate separately from CBF
% and CMRO2, so this is an oversimplification, but for now ...
dCBF = (dCMRO2-1)*3 + 1;  % just define CBF from CMRO2
                          % taking a "standard" assumption of
                          % CBF = 3*CMRO2.
                          
% Model assumes that blood flow is the driving force, and 
% determines volume.  Basic equation is 
%   dv/dt = 1/tau * [f_in(t) - f_out(v,t)]
% f_in(t) is the specified dCBF function
% One solution that Buxton has worked out is:
%   f_out(v,t) = dCBV^(1/alpha) + tau_v*dv/dt
%   (Equation from Mildner et al. (2001), MRM46:819.)
%   where tau_0 is the transit time through the balloon
%           under resting conditions
%         tau_v describes the resistance of the balloon 
%           to changes (stiffness of balloon?).  

% Substituted back into the original equation, this gives
%   dv/dt = ______1______{f_in(t) - v(t)^1/alpha}
%           tau_0 + tau_v
% which describes how the volume increases as flow increases, but
% also how the balloon tends to shrink when volume is large
%
% Parameters for modeling blood volume and flow changes
alpha = .4; % volume = flow^alpha.
tau_v = 10; % activated state transit time???
tau_0 = 1.5; % resting state transit time through balloon

% Numerically calculate dv/dt from the given dCBF, assuming
% that v(0) = 1.
dCBV(1) = 1; dvdt(1) = 0;
for it = 2:length(dCBF)
    dvdt(it) = dt/(tau_0 + tau_v)* ...
        (dCBF(it-1) - dCBV(it-1)^(1/alpha));
    dCBV(it) = dCBV(it-1) + dvdt(it);
end
%
%just as a reality check to see that my flow-out vs. volume curve
% looks like the hysteresis plots in the book ...

%flow_out = dCBV.^(1/alpha) + tau_v*dvdt;
flow_out = dCBF - (tau_0/tau_v)*(dCBF-dCBV.^(1/alpha));
subplot(2,2,1); plot(dCBV,flow_out); xlabel('volume'); ylabel('flow_o_u_t')
  title(['tau_v = ' num2str(tau_v,2)]); axis tight
subplot(2,2,2); plot(timebase,[dCBF;dCBV;flow_out])
  xlabel('time'); ylabel('normalized parameter')
  legend('CBF','CBV','flow_o_u_t')
  %% something's obviously wrong with my equations if I end up
  % with an initial negative flow_out!!!
  
% constants to be used in converting dHb and volume to BOLD
V0 = 0.02; % fractional volume of blood in voxel at rest
A1 = 3.4; % a field-dependent constant including effects of 
  % oxygen extraction (resting), frequency shift of fully
  % deoxygenated blood, echo time of image acquisition, ratio
  % of intra- to extravascular signal, and dependence of 
  % relaxation rates on oxygenation of blood. Unitless.
  % Given value is appropriate for 1.5T.
A2 = 1.0; % constant similar to A1, including all of the above
  % except frequency shift of fully deoxygenated blood. Also 
  % field dependent, and current value is set for 1.5T.

% To calculate local deoxyhemoglobin concentration, include
% oxygen extraction fraction (what percent of available
% oxygen is taken out of blood by tissue), rate at which blood
% enters voxel (flow_in), rate at which blood leaves voxel
% (flow_out), and volume of blood in voxel, dCBV.
% Equation is:
% 
% dq/dt = (1/tau) * [flow_in(t)*E(t)/E0 - q(t)/v(t)*flow_out(v,t)]
%     where tau is the mean transit time through balloon at rest
%           E(t) is time-dependent oxygen extraction fraction
%           E0 is resting OEF
%           q(t) is normalized deoxyhemoglobin concentration
%           v(t) is dCBV
%
% It is worth noting that flow_in(t)*E(t)/E0 is CMRO2, normalized
% to resting values.  This model is therefore built to specify
% CBF, and model CMRO2 through a variable E(t).  I don't like that
% but will go with it for now.

subplot(2,2,3)
plot(timebase,[dCMRO2;dCBF;dCBV]); legend('CMRO_2','CBF','CBV')
    
% Then:
% dq/dt = (1/tau) * [dCMRO2 - q(t)/dCBV*flow_out)]
q(1) = 1; 
for it = 2:length(dCBF)
    q(it) = q(it-1) + (1/tau_0)*(dCMRO2(it-1) - ...
        q(it-1)/dCBV(it-1)*flow_out(it-1));
end

BOLD = V0*(A1*(1-q)-A2*(1-dCBV));
subplot(2,2,4); plot(timebase,[100*BOLD],'k'); axis tight
 xlabel('time'); ylabel('percent change, BOLD')
 set(gcf,'Color',[1 1 1])