Support:Documents:Examples:Estimate Input Parameters Simultaneously with Flow

Estimating Input Function Parameters Simultaneously with Flow

This example demonstrates an approach to simultaneously estimating input functions parameters along with blood flow in the context of the PET blood flow model.

Tissue Uptake

The flow model has one tissue compartment. Radioactive water (or another freely diffusible molecule) in the blood is assumed to rapidly exchange with the water in the tissue. The tissue model has two rate constants: K₁ and k₂ and can be depicted in a diagram:

It can also be described by the differential equation:

dC_T/dt = K₁ C_p - k₂ C_T

where C_T is the tissue concentration and C_p is the plasma concentration of radioactive water. C_T and C_p are interpreted as either molar concentrations or decay-corrected activity concentration (Salinas 2007).

Instead of estimating K₁ and k₂, it sometimes is advantageous to estimate a different set of parameters defined as

lambda = K₁/k₂ and F = K₁.

F is blood flow (perfusion) per unit volume of tissue (for simplicity we're assuming 100% extraction) and lambda is the ratio of tissue to blood concentrations at equilibrium.

To use this parameterization in COMKAT, one needs to express the rate constants in terms of lambda and F: K₁ = F, k₂ = K₁/lambda

This example demonstrates how to estimate F (simultaneously with input function parameters). The example could be easily modified to estimate lambda as well.

Input Function

This example also demonstrates how to use a parameterized input function. Specifically, the plasma concentration vs. time t is modeled as

C_p = (p₁(t-d) - p₂ - p₃) exp(p₄ (t-d)) + p₂ exp(p₅ (t-d)) + p₃ exp(p₆(t-d)) [aka Feng input]

d is a time-delay and C_p is defined to be zero for t < d.

In this example, the values of the parameters p₁, p₂, ... p₆ will be estimated.

Note: The Feng input function has been implemented in fin.m in the COMKAT examples folder.

Pixel Values

In a PET image, pixel values represent the radioactivity concentration. There are contributions from radioactivity in blood and tissue. The arterial blood concentration includes radioactive water in plasma and the blood cells and is designated C_a. Here we assume C_a can be described using the Feng input equation but with different values for p₁, p₂, ... p₆ than were used for C_p. The tissue activity concentration is C_T x S where S is the exponentially-decaying specific activity and C_T is given above. The pixel activity concentration is calculated as the activity concentration in the tissue plus the pixel's fractional blood space F_v times the arterial blood activity concentration: PET = C_TS + F_vC_a.

COMKAT Implementation

Implement the Compartment Model for Calculating Output Concentrations

First create an new model blood flow compartment model (bfcm) and add compartments. The J compartment is not shown in the figure above but is needed to accept the flux from the k₂ arrow.

% create new, empty model
bfcm = compartmentModel; 
 
% add compartments
bfcm = addCompartment(bfcm, 'CT');
bfcm = addCompartment(bfcm, 'J');

The specific activity (ratio of activity to molar concentration) is an exponentially decaying function S(t) = S₀ exp(-0.341 t) where t is expressed in units of minutes and 0.341 = ln(2)/half-life of ¹⁵O and S₀ is the initial specific activity:

% define the initial specific activity and decay constant
sa0 = 1; % specific activity at t=0
dk = 0.341; % O-15 decay constant

Next specify the input functions C_p and C_a for the compartment model. fin is the name of the function that implements Feng input. pCp and pCa are parameter vectors that hold the values of p₁, p₂, ... p₆ for C_p and C_a.

bfcm = addInput(bfcm, 'Cp', sa0, dk, 'fin', 'pCp');
bfcm = addParameter(bfcm, 'pCp', 2*[851.1 21.88 20.81 -4.134 -0.1191 -0.01043 .5]);

bfcm = addInput(bfcm, 'Ca', sa0, dk, 'fin', 'pCa');
bfcm = addParameter(bfcm, 'pCa', [851.1 21.88 20.81 -4.134 -0.1191 -0.01043 .5]);

Define the connections (arrows) between the compartment and input. The link type 'L' indicates linear kinetics from an input to a compartment and type 'K' indicates linear kinetics from one compartment to another.

% connect compartments and inputs
bfcm = addLink(bfcm, 'L', 'Cp', 'CT', 'K1');
bfcm = addLink(bfcm, 'K', 'CT', 'J', 'k2');
 
% define values for K1, k2 ...
bfcm = addParameter(bfcm, 'K1', 'F');
bfcm = addParameter(bfcm, 'k2', 'K1/lambda');
bfcm = addParameter(bfcm, 'F', 0.5);
bfcm = addParameter(bfcm, 'lambda', 1.0);
bfcm = addParameter(bfcm, 'Fv', 0.05);

Define two model outputs. The output called TissuePixel predicts the activity in a tissue pixel and the output called ArterialPixel predicts the activity that would be observed in a pixel in the aorta. TissuePixel corresponds to C_TS + F_vC_a whereas AterialPixel is 100% C_a.

% define the activity concentration in tissue pixel
bfcm = addOutput(bfcm, 'TissuePixel', {'CT', '1'}, {'Ca', 'Fv'});
 
% define the activity concentration in arterial blood pixel
bfcm = addOutput(bfcm, 'ArterialPixel', {'CT', '0'}, {'Ca', '1'});

Note that COMKAT calculates the values of the model outputs as time-averaged concentrations over the scan frames. These are defined as sequence of contiguous 10-second frames (frame start times 0, 10, 20, ... 590 sec and frame end times 10, 20, 30, ... 600).

% specify scan frame times
bfcm = set(bfcm, 'ScanTime', [[0:10:590]' [10:10:600]']/60);

Finally, tell COMKAT to solve the model to obtain simulation data that will be used in fitting.

[PET, PETIndex, Output, OutputIndex] = solve(bfcm);

% plot activity in tissue PET(:,3) and blood PET(:,4)
midTime = (PET(:,1) + PET(:,2)) / 2;
plot(midTime, PET(:,3), 'b-', midTime, PET(:,4), 'r.')
legend(PETIndex{3}, PETIndex{4})

% make the figure for this wiki
set(1,'PaperPosition',[0.25 2.5 3 2])
print -dpng c:\temp\ExampleFigureBloodFlowInputModelDataToFit.png

Note, different rows of PET correspond to different frames. The columns of PET hold various quantities as indicated in PETIndex. That is, the contents of column i are described in PETIndex{i}.

Fit Data

Modify the COMKAT model a bit so that it can be used to fit data created above.

Specify the Data to Fit

Just use the result from above

data = PET(:, [3 4]);
bfcm = set(bfcm, 'ExperimentalData', data);

Specify parameters to estimate along with their initial guesses and range of feasible values.

 
% specifyparameters to estimate
bfcm = addSensitivity(bfcm, 'pCp', 'F');

% define initial guess and lower and upper bounds
pCpInit = 0.9*2*[851.1 21.88 20.81 -4.134 -0.1191 -0.01043 .5];
FInit = 0.3;
pinit = [pCpInit FInit];
lb = [2*[85 2 2 -8 -3 -0.1 .05] 0.1];
ub = [2*[2000 100 200 -1 -0.01 -0.001 2] 1.];

Note the values of the initial guesses for the input function parameters are 0.9 times the true values. (We're are not cheating by starting with the correct values.)

Fit the data to estimate values of the parameters

[pEst, modelFit] = fit(bfcm, pInit, lb,ub);

COMKAT makes fitting easy!

Plot the results and make the figure below for this wiki

plot(midTime, modelFit(:,1), 'b', midTime, PET(:,3), 'b.', midTime, modelFit(:,2), 'r', midTime, PET(:,4), 'r.')
legend('Tissue Fit', PETIndex{3}, 'Arterial Fit', PETIndex{4})
xlabel('Time')
ylabel('Activity')
print -dpng c:\temp\ExampleFigureBloodFlowInputModelFit.png

Ta Da

It looks like a pretty good fit.

Text

Write some introduction to your example.

Code

You can put code in a different format using <pre> tags like:

For text that you want to appear as part of a MATLAB function or program, encase it inside a "pre" block.
Line 2 of pre block
    Line 3 (indented) of pre block.

Edit the current page to see how this is done.

Figures

Remember to include figures and plots in your example to make it look pretty and convey information. To insert a figure, create a figure file in png (preferred) or other graphic format.

figure
plot(1:10)
title('Example Figure')
ylabel('y-axis label')
xlabel('x-axis label')
set(1,'PaperPosition',[0.25 2.5 3 2])
print -dpng ExampleFigure.png

Next, upload the image file ("upload file" at the bottom of any wiki page) giving it a unique name.

Use double-square brackets to insert the figure into the page

Edit the example index page to include your example in the list

Add your new example to the list of examples on the example index page by adding a line like this to the index page:

[http://comkat.case.edu/comkat/comkat_wiki/index.php?title=Support:Documents:Examples:Instructions_for_Adding_Your_Own_Example The title of you example]

Within the brackets, the first part is the URL address to your example page. Then you type your example title following that address. On the index page it will look like:

The title of you example