MATLAB Functions

B_CURVE a two-compartment rCBF model used for delay correction

      integral = b_curve (args, shifted_g_even, ts_even, A, ... 
                          fstart, flengths)

Used by blood delay correction, this function implements a two-compartment rCBF model used for fitting the blood curve data to the head curve data. CORRECTBLOOD perform delay and dispersion corrections on blood curve

   [new_ts_even, Ca_even, delta] = correctblood (A, FrameTimes, ...
                                   FrameLengths, g_even, ts_even, progress)

The required input parameters are:

A - brain activity, averaged over all gray matter in a slice. This should be in units of decay / (gram-tissue * sec), and should just be a vector - one value per frame.

FrameTimes - the start time of every frame, in seconds

FrameLengths - the length of every frame, in seconds

g_even - the (uncorrected) arterial input function, resampled at some *evenly spaced* time domain. Should be in units of decay / (mL-blood * sec)

ts_even - the time domain at which g_evenis resampled

The returned variables are:

new_ts_even - generally the same as the old time scale, with some points missing from the end.

Ca_even - g_evenwith dispersion and delay hopefully corrected, in units of decay / (mL-blood * sec).

delay - the delay time (ie. shift) in seconds

A, FrameTimes, and FrameLengths must all be vectors with the same number of elements (presumably the number of frames in the study). g_evenand ts_evenmust also be vectors with the same number of elements, but their size should be much larger, due to the resampling at half-second intervals performed by resampleblood.

correctblood corrects for dispersion in blood activity by calculating g(t) + tau * dg/dt, where tau (the dispersion time constant) is taken to be 4.0 seconds.

It then attempts to correct for delay by fitting a theoretical blood curve to the observed brain activity A(t). This curve depends on the parameters alpha, beta, gamma (these correspond to K1, k2, and V0, although for the entire slice rather than pixel-by-pixel) and delta (which is the delay time). correctblood steps through a series of delta values (currently -5 to +10 sec), and performs a three- parameter fit with respect to alpha, beta, and gamma; the value of delta that results in the best fit is chosen as the delay time.

options is an entirely optional vector meant for debugging purposes. If options(1) is non-zero, then correctblood will show its progress, by printing out the results of progressive delay-correction fits. If it is at least 2, then correctblood will also show progress graphically, by displaying a graph of A(t) and the fits corresponding to every value of delta tried. If options(2) is zero, then no delay correction will be performed; if options(3) is supplied, then it will be used as delta to do delay correction without the time-consuming fitting. FINDINTCONVO calculate tables of the integrated convolutions commonly used

    [int1,int2,int3] = findintconvo (Ca_even, ts_even, k2_lookup,...
                                     midftimes, flengths, w1[, w2[, w3]])

given a table of k2 values, generates tables of weighted integrals that commonly occur in RCBF analysis. Namely, int_convois a table of the same size as k2_lookupcontaining

int ( conv (Ca(t), exp(-k2*t)) * weight )

where the integration is carried out across frames. weight is one of w1, w2, or w3, each of which will generally be some simple function of midftimes. findintconvo will return int2 if and only if w2 is supplied, and int3 if and only if w3 is supplied. w1 is required, and int1 will always be returned. Normally, the weight functions should be vectors with the same number of elements as midftimes; however, if w1 is empty then the weighting function is taken to be unity.

Note that in order to correctly calculate the convolution, Ca(t) must be resampled at evenly spaced time intervals, and this resampled blood activity should be passed as Ca_even. The times at which it is sampled should be passed as ts_even. (These can be calculated by resampleblood before calling findconvints.)

Then, the convolution of Ca(t) and exp(-k2*t) is resampled at the mid-frame times (passed as midftimes) and integrated across frames using flengths as dt. RCBF1 a one-compartment (double-weighted integral) rCBF model.

         [K1,k2] = rcbf1 (filename, slice)

A one-compartment rCBF model (without V0 or blood delay and dispersion) implemented as a MATLAB function. The compartmental equation is solved by integrating it across the entire study, and then weighting this integral with two different weights. When these two integrals are divided by each other, K1 is eliminated, leaving only k2. A lookup table is calculated, relating values of k2 to values of the integral. From this, k2 and be calculated. From k2, K1 is easily found by substitution into the original compartmental equation. See the document "rCBF Analysis Using Matlab" for further details of both the compartmental equations themselves, and the method of solution. RCBF2 a two-compartment (triple-weighted integral) rCBF model.

        [K1,k2,V0,delta] = rcbf2 (filename, slice)

rcbf2 implements the three-weighted integral method of calculating k2, K1, and V0 (in that order) for a particular slice. This function also returns the delay value calculated for blood correction. It first reads in a great mess of data (viz., the brain activity for every frame of the slice, frame start times and lengths, plasma activity, and blood sample times). Then, a simple mask is created and used to filter out roughly all points outside the head.

The actual calculations follow the procedure outlined in the document "RCBF Analysis Using Matlab". Occasionally, comments in the source code or documentation for various functions involved in the analysis will refer to equations in this document. The most relevant functions in this respect are rcbf2 itself, correctblood and findintconvos.

The starting point of the three-weighted integration method is Eq. 10 of the RCBF document. The left hand side of this equation, rL, is calculated for every pixel. Then, a lookup table relating a series of k2 values to the rR (the right-hand side of Eq. 10) is calculated. This lookup table should have a few hundred elements, as calculating rR is considerably more expensive than calculating rL. Since rL and rR are equal, we use the pixel-wise values of rL to lookup k2 for every pixel.

Then, Eq. XXX is used to calculate K1. This requires calculating the moderately complicated int_activity(left hand side of Eq. XXX) and the extremely complicated k2_conv_ints(right hand side). However, the expression for k2_conv_intsappeared already in the numerator of rR, so we preserve that lookup table as conv_int1and use it to lookup k2_conv_ints. These two long vectors (with one number for every pixel) are then divided to get K1. Finally, V0 is calculated via Eq. YYY. RCBFDEMO Demonstrate the RCBF blood analysis package.

    rcbfdemo(slice_number, frame_number)

Hard-coded to use the yates_19445data file, but allows input of a slice and frame to display.