Support:Documents:Examples:FDG with Time-varying Rate Constants

FDG Model with Time-varying Rate Constants

Overview

This example demonstrates how to implement the two-tissue compartment model with extended such that K₁ and k₂ vary in time.

This models the case where there are physiologic changes during the course of the PET scanning. This might happen if tumors are being treated during scanning.

The standard FDG model was described by Phelps and Huang .

It can also be described by the differential equations:

dC_e/dt = K₁ C_p - k₂ C_e - k₃ C_e + k₄ C_m

dC_m/dt = k₃ C_e - k₄ C_m

C_e is the tissue concentration of FDG and C_m as the tissue concentration of metabolized FDG (FDG-6-phosphate)

C_e and C_m are interpreted as molar concentrations (Salinas, Muzic and Saidel 2007).

In the standard model. the rate constants K₁, k₂, k₃ and k₄ are assumed to be independent of time.

Here we allow K₁ and k₂ to be linear functions of time:

K₁(t) = K₁⁰ + a t

k₂(t) = K₂⁰ + b t

a and b are the derivatives (slopes) of the rate constant with respect to time.

To implement the time-varying rate constants, we create a MATLAB function in a .m file. By making the function general, the same function can be used for K₁ and k₂.

function varargout = kinLinearTime(t, c, ncomp, nsens, pName, pValue, pSensIdx, cm, xtra)
% rate constant varies linearly with time: k = k0 + m*t  where m = dk/dt

K0Name = pName{1};
dkdtName = pName{2};
K0Value = pValue{1};
dkdtValue = pValue{2};
K0SensIndex = pSensIdx{1};
dkdtSensIndex = pSensIdx{2};

if (nargout > 0), % return effective rate constant
    v = K0Value + dkdtValue * t;
    if (v < 0), % don't allow negative values for rate constant
        varargout{1} = 0;
    else
        varargout{1} = v;
    end        
    if (nargout > 1),  % return dk/dc 
        dkdc = zeros([1 ncomp]); 
        varargout{2} = dkdc;  
        if (nargout > 2), % return dk/dp
            dkdp = zeros([1 nsens]);
            if (v >= 0),
                dkdp(K0SensIndex) = 1;
                dkdp(dkdtSensIndex) = t;
            end        
            varargout{3} = dkdp; % dkdp
            if (nargout > 3),
                ddkdpdc = zeros([ncomp nsens]);
                varargout{4} = ddkdpdc; % ddkdpdc
            end
        end
    end
end
return

This function is pretty straightforward. The emphasis of our description here is on the mathematical aspects. The first output argument varargout{1} holds the value of the rate constant specified at the current time t. This value is first calculated and stored in the variable v.

v = K0Value + dkdtValue * t

k0 is the rate constant at time zero (t=0) and dkdtValue is the increase in the rate constant per unit time (time-derivative or slope). k0Value and dkdtValue are the MATLAB variables that hold the values of these parameters. The value of v is checked to see if it is negative and, if it is, the value zero is used instead since negative values for rate constants are non-physiologic.

The line

varargout{2} = dkdc;  

returns the derivative of the rate constant with respect to all the compartment concentrations (a vector). In this case, the rate constant is not dependent on the concentrations so the derivatives are all zeros.

The lines

dkdp(K0SensIndex) = 1;
dkdp(dkdtSensIndex) = t;

calculate the derivatives of the rate constant with respect to the model parameters k0 (the initial value of the rate constant) and dkdt (the increase in k per unit time). Since a model can have other parameters besides these, the function stores these derivatives in the appropriate element of the derivative vector. Derivatives of this rate constant with respect to all other parameters are zeros.

The line

varargout{4} = ddkdpdc; % ddkdpdc

returns the values of the mixed derivatives (derivative of rate constant with respect to concentrations and with respect to time), which, in this case are all zeros.

The function is stored in kinLinearTime.m in the comkat examples folder.

This function is called to evaluate both K₁ and k₂.

Implementing the Compartment Model

Create a new model

First define some MATLAB variables for clarity and to facilitate exploring what happens if values are changed

K10 =  0.15;  % K1 = K10 + dK1dt * t
dK1dt = -K10 * 0.5 / 60;  % 50% decrease per 60 minutes
k20 =  0.325; % k2 = k20 + dk2dt * t
dk2dt = -k20 * 0.5 / 60;  % 50% decrease per 60 minutes
k3 = 0.05;
k4 = 0.007;

Next, create a compartment model object and define the parameters within the model object

% create empty compartmentModel object
cm = compartmentModel;
  
% define the parameters
cm = addParameter(cm, 'K10',   K10);    % 1/min
cm = addParameter(cm, 'k20',   k20);    % 1/min
cm = addParameter(cm, 'k3',    k3);     % 1/min
cm = addParameter(cm, 'k4',    k4);     % 1/min
cm = addParameter(cm, 'dK1dt', dK1dt);  % 1/min/min
cm = addParameter(cm, 'dk2dt', dk2dt);  % 1/min/min

cm = addParameter(cm, 'sa',    1);            % specific activity of injection
cm = addParameter(cm, 'dk',    log(2)/109.8); % F-18 radioactive decay
cm = addParameter(cm, 'PV',    1);            % (none)

Now, define the model compartments. The first two are Ce and Cm as described above and Junk is a "sink" that collects the FDG that is cleared by the venous circulation.

cm = addCompartment(cm, 'Ce');
cm = addCompartment(cm, 'Cm');
cm = addCompartment(cm, 'Junk');

For the plasma concentration time-course of FDG or input function, we'll use the Feng input which is an analytic expression. Since piecewise polynomials (linear interpolation, cubic splines, etc...) are widely used in COMKAT and their implementation has been optimized, here we construct a cubic spline representation of the input function by sampling the Feng input.

x0=[2 0.5 16.000 2.188 2.081 -7.5 -10 -0.01043]; % Feng parameters
t = [0:.1:3 3.5 4 4.5 5:1:10 12 15:5:(expDuration+5)];
Cp = fengInput(x0, t);
ppCp = spline(t, Cp);
cm = addInput(cm, 'Cp', 'sa', 'dk', 'ppval', ppCp);  % Cp has units of pmol/ml

The PET scanner measurement is assumed to represent the sum of FDG and FDG-6-Phosphate so the model output is calculated as the sum of the C_e and C_m compartments

Wlist = {...
    'Ce', 'PV';
    'Cm', 'PV'};
Xlist = {};
cm = addOutput(cm, 'PET', Wlist, Xlist);

Next define the scan frames as the start and stop time of each image in the dynamic sequence.

ttt=[ ones(12,1)*10/60; ...  % 10 sec
      ones(10,1)*0.5; ...    %  0.5 min
      ones(10,1)*2;...       %  2 min
      ones(10,1)*5;...       %  5 min
      ones(4,1)*10];         % 10 min
scant=[[0;cumsum(ttt(1:(length(ttt)-1)))] cumsum(ttt)];
cm = set(cm, 'ScanTime', scant);

Note the acquisition begins with 12 frames of 10 sec duration, then has 10 frames of 0.5 minute duration, etc... The cumsum() function calculates the cumulative summation. Thus, the first column of scant holds the scan start times and the second column holds the scan end times. These are expressed in minutes:

scant =

         0    0.1667
    0.1667    0.3333
    0.3333    0.5000
    0.5000    0.6667
    0.6667    0.8333
    0.8333    1.0000
    1.0000    1.1667
    1.1667    1.3333
    1.3333    1.5000
    1.5000    1.6667
    1.6667    1.8333
    1.8333    2.0000
...

% make copy of model thus far cmRef = cm;

% define rate constants cm = addLink(cm, 'LSpecial', 'Cp', 'Ce', 'K1', 'kinLinearTime', {'K10', 'dK1dt'}); cm = addLink(cm, 'KSpecial', 'Ce', 'Junk', 'k2', 'kinLinearTime', {'k20', 'dk2dt'}); cm = addLink(cm, 'K', 'Ce', 'Cm', 'k3'); if strcmp(modelConfig, '4k'),

   cm = addLink(cm, 'K',        'Cm', 'Ce',   'k4');

end

% define reference model with time-invariant rate constants cmRef = addLink(cmRef, 'L', 'Cp', 'Ce', 'K10'); cmRef = addLink(cmRef, 'K', 'Ce', 'Junk', 'k20'); cmRef = addLink(cmRef, 'K', 'Ce', 'Cm', 'k3'); if strcmp(modelConfig, '4k'),

   cmRef = addLink(cmRef, 'K',        'Cm', 'Ce',   'k4');

end

% solve both models and compare outputs [PET, PETIndex] = solve(cm); [PETRef, PETRefIndex] = solve(cmRef);

tmid = (PET(:,1) + PET(:,2)) / 2;

figure clf axK1k2 = axes('position', [0.2 0.75 0.7 0.2]); phK = plot(axK1k2, tmid, K10 + dK1dt * tmid, 'b', tmid, k20 + dk2dt * tmid, 'r', ...

            tmid, K10 + 0 * tmid, 'b:',                tmid, k20 + 0 * tmid, 'r:');

set(phK, 'LineWidth', 1.5); ax = axis; ax(3) = 0; ax(4) = 0.5; axis(ax); ylabel('Rate Const. (1/min)'); lh = legend('K1', 'k2', 'K1', 'k2', 'Location', 'SouthEast'); set(lh, 'fontsize', 8);

axPET = axes('position', [0.2 0.2 0.7 0.4]); phPET = plot(axPET, tmid, PET(:,3), 'b', tmid, PETRef(:,3), 'b:'); set(phPET, 'LineWidth', 1.5); ax = axis; ax(3) = 0; axis(ax); xlabel('Time (Minutes)'); ylabel('Activity (uCi/ml)'); legend('Time-varying', 'Reference', 'Location', 'SouthEast');

% What would be the value of rate constants estimated while neglecting the time-variation? % Fit the time-invariant model to the time-varying model data. data = PET(:,3); cmRef = set(cmRef, 'ExperimentalData', data); pinit = [0.07; 0.1; 0.04; 0.005]; plb = [0.01; 0.001; 0.001; 0.0001]; pub = [0.5; 0.5; 0.2; 0.01]; if strcmp(modelConfig, '4k'),

   cmRef = addSensitivity(cmRef, 'K10', 'k20', 'k3', 'k4');

else

   cmRef = addSensitivity(cmRef, 'K10', 'k20', 'k3');
   pinit(4) = [];
   plb(4) = [];
   pub(4) = [];

end

[pfit, modelFit] = fit(cmRef, pinit, plb, pub);

figure phFit = plot(tmid, data, 'bo', tmid, modelFit, 'b:'); set(phFit, 'LineWidth', 1.5); ax = axis; ax(3) = 0; axis(ax); xlabel('Time (Minutes)'); ylabel('Activity (uCi/ml)');

KiRange = [K10*k3/(k20 + k3); (K10 + dK1dt)*k3/(k20 + k3); K10*k3/((k20 + dk2dt) + k3); (K10 + dK1dt)*k3/((k20 + dk2dt) + k3)]; KiInit = K10*k3/(k20 + k3); if strcmp(modelConfig, '4k'),

   c = {'Parm.', 'True', 'Est', 'vs. Init', 'Units'; ...
       'K1', sprintf('%1.3f -> %1.3f', K10, K10 + dK1dt * PET(end,2)), sprintf('%1.3f', pfit(1)), sprintf('%1.1f%%', 100*(pfit(1)-K10)/K10), '1/min';
       'k2', sprintf('%1.3f -> %1.3f', k20, k20 + dk2dt * PET(end,2)), sprintf('%1.3f', pfit(2)), sprintf('%1.1f%%', 100*(pfit(2)-k20)/k20), '1/min';
       'k3', sprintf('%1.3f', k3),                                     sprintf('%1.3f', pfit(3)), sprintf('%1.1f%%', 100*(pfit(3)-k3)/k3),   '1/min';
       'k4', sprintf('%1.3f', k4),                                     sprintf('%1.3f', pfit(4)), sprintf('%1.1f%%', 100*(pfit(4)-k4)/k4),   '1/min';
       'Ki', sprintf('%1.3f - %1.3f', min(KiRange), max(KiRange)),     sprintf('%1.3f', pfit(1)*pfit(3)/(pfit(2)+pfit(3))), sprintf('%1.1f%%', 100*(pfit(1)*pfit(3)/(pfit(2)+pfit(3))-KiInit)/KiInit), '1/min';
      };

else

   c = {'Parm.', 'True', 'Est', 'vs. Init', 'Units'; ...
       'K1', sprintf('%1.3f -> %1.3f', K10, K10 + dK1dt * PET(end,2)), sprintf('%1.3f', pfit(1)), sprintf('%1.1f%%', 100*(pfit(1)-K10)/K10), '1/min';
       'k2', sprintf('%1.3f -> %1.3f', k20, k20 + dk2dt * PET(end,2)), sprintf('%1.3f', pfit(2)), sprintf('%1.1f%%', 100*(pfit(2)-k20)/k20), '1/min';
       'k3', sprintf('%1.3f', k3),                                     sprintf('%1.3f', pfit(3)), sprintf('%1.1f%%', 100*(pfit(3)-k3)/k3),   '1/min';
       'Ki', sprintf('%1.3f - %1.3f', min(KiRange), max(KiRange)),     sprintf('%1.3f', pfit(1)*pfit(3)/(pfit(2)+pfit(3))), sprintf('%1.1f%%', 100*(pfit(1)*pfit(3)/(pfit(2)+pfit(3))-KiInit)/KiInit), '1/min';
      };

end

table(c, [0.4, 0.3] ); legend('Time-varying Data', 'Time-invariant Fit', 'Location', 'SouthEast');

First create an new model stored in variable bfcm (blood flow compartment model) 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');

Specify the input

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 specific activity at time t = 0.

% 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 tissue and blood outputs

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 ArterialPixel 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 even accounts for the time-averaging of concentrations over the durations of the scan frames. For simplicity here a sequence of contiguous 10-second frames is used. (Frame start times 0, 10, 20, ... 590 sec and frame end times 10, 20, 30, ... 600 sec).

% specify scan frame times
bfcm = set(bfcm, 'ScanTime', [[0:10:590]' [10:10:600]']/60); % division by 60 converts sec to min

Solve for the model output

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.

% specify parameters 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.)

Do the actual fitting

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.