Difference between revisions of "Support:Documents:Examples:FDG with Time-varying Rate Constants"

From COMKAT wiki
Jump to navigation Jump to search
(New page: == FDG Model with Time-varying Rate Constants == This example demonstrates how to implement the two-tissue compartment model with extended such that K1 and k2 vary in time. This models th...)
 
 
(7 intermediate revisions by the same user not shown)
Line 1: Line 1:
 
== FDG Model with Time-varying Rate Constants ==
 
== FDG Model with Time-varying Rate Constants ==
−
This example demonstrates how to implement the two-tissue compartment model with extended such that K1 and k2 vary in time.
+
 
 +
==== Overview ====
 +
This example demonstrates how to implement the two-tissue compartment model with extended such that K<sub>1</sub> and k<sub>2</sub> 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.
 
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.
  
−
==== Tissue uptake ====
+
The standard FDG model was described by
−
This model was described by Huang and Phelps.
+
[http://www.ncbi.nlm.nih.gov/pubmed/117743?ordinalpos=2&itool=EntrezSystem2.PEntrez.Pubmed.Pubmed_ResultsPanel.Pubmed_RVDocSum Phelps]
 +
and
 +
[http://www.ncbi.nlm.nih.gov/pubmed/6965568?ordinalpos=4&itool=EntrezSystem2.PEntrez.Pubmed.Pubmed_ResultsPanel.Pubmed_RVDocSum Huang]
 +
.
  
−
MAKE FIGURE
 
 
[[Image:ExampleFigureFDGModel.png]]
 
[[Image:ExampleFigureFDGModel.png]]
 +
  
 
It can also be described by the differential equations:
 
It can also be described by the differential equations:
Line 25: Line 30:
  
  
−
Here we allow K<sub>1</sub> and k<sub>2</sub> to be a linear function of time:
+
Here we allow K<sub>1</sub> and k<sub>2</sub> to be linear functions of time:
  
 
K<sub>1</sub>(t) = K<sub>1</sub><sup>0</sup> + a t
 
K<sub>1</sub>(t) = K<sub>1</sub><sup>0</sup> + a t
  
−
k<sub>2</sub>(t) = K<sub>2</sub><sup>0</sup> + b t
+
k<sub>2</sub>(t) = k<sub>2</sub><sup>0</sup> + b t
 +
 
 +
K<sub>1</sub><sup>0</sup> and k<sub>2</sub><sup>0</sup> are the initial (0 superscript) values of the rate constants
 +
and a and b are the derivatives (slopes) of the rate constant with respect to time.
 +
 
 +
 
 +
To implement the time-varying rate constants, create a MATLAB function in a .m file.
 +
By making the function general, the same function can be used for K<sub>1</sub> and k<sub>2</sub>.
 +
 
 +
<pre>
 +
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
 +
</pre>
 +
 
 +
 
 +
This function has a similar pattern to any function one might write to implement customized kinetic rules.
 +
Note the function is called numerous times as the model equations are solved.  The value of t (time) takes on values from 0 to the end of the last frame and many values in betweeen - even values that do not correspond to the frame begin and end times.
 +
 
 +
The first part of the function retrieves values that are automatically sent by COMKAT when it calls this function.
 +
<pre>
 +
K0Name = pName{1};
 +
dkdtName = pName{2};
 +
K0Value = pValue{1};
 +
dkdtValue = pValue{2};
 +
K0SensIndex = pSensIdx{1};
 +
dkdtSensIndex = pSensIdx{2};
 +
</pre>
 +
 
 +
To get COMKAT to call this function, a link (connection between compartments) would be included in the model.  For example
 +
<pre>
 +
cm = addLink(cm, 'LSpecial', 'Cp', 'Ce',  'K1', 'kinLinearTime', {'K10', 'dK1dt'});
 +
</pre>
 +
 
  
 +
When COMKAT calls the custom kinetic function, the last argument of the addLink command ({'K10', 'dK1dt'}) will be passed in as the pName argument.
 +
pValue will be a cell array with the value of pValue{i} set to the value of the parameter named in pName{i}.  (Although you could use pxEval to obtain the value of the parameter from pName{i}, it is more efficient for COMKAT to do this internally (only once and cache the result).
 +
pSensIdx will be a cell array with the value of pSensIdx{i} being the position (index) of parameter named in pName{i} in the list of parameters to be estimated.
 +
If the parameter named in pName{i} was not included in an addSensitivity() call, then pSensIdx{i} will be an empty matrix.
  
−
To implement the time-varying rate constants, it is necessary to
+
The emphasis of the rest of the description of kinLinearTime.m is on the mathematical details.
 +
 
 +
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.
 +
<pre>
 +
v = K0Value + dkdtValue * t
 +
</pre>
  
 +
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.
  
−
== Implement the Compartment Model ==
+
The line
 +
<pre>
 +
varargout{2} = dkdc; 
 +
</pre>
 +
 
 +
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
 +
<pre>
 +
dkdp(K0SensIndex) = 1;
 +
dkdp(dkdtSensIndex) = t;
 +
</pre>
 +
 
 +
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
 +
<pre>
 +
varargout{4} = ddkdpdc; % ddkdpdc
 +
</pre>
 +
 
 +
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<sub>1</sub> and k<sub>2</sub>.
 +
 
 +
 
 +
== Implementing the Compartment Model ==
  
 
==== Create a new model ====
 
==== Create a new model ====
−
First create an new model stored in variable bfcm (blood flow compartment model) and add compartments. 
+
First define some MATLAB variables for clarity and to facilitate exploring what happens if values are changed
−
The J compartment is not shown in the figure above but is needed to accept the flux from the k<sub>2</sub> arrow.
+
 
 
<pre>
 
<pre>
−
% create new, empty model
+
K10 =  0.15;  % K1 = K10 + dK1dt * t
−
bfcm = compartmentModel;  
+
dK1dt = -K10 * 0.5 / 60; % 50% decrease per 60 minutes
−
   
+
k20 = 0.325; % k2 = k20 + dk2dt * t
−
% add compartments
+
dk2dt = -k20 * 0.5 / 60;  % 50% decrease per 60 minutes
−
bfcm = addCompartment(bfcm, 'CT');
+
k3 = 0.05;
−
bfcm = addCompartment(bfcm, 'J');
+
k4 = 0.007;
 
</pre>
 
</pre>
  
  
−
==== Specify the input ====
+
Next, create a compartment model object and define the parameters within the model object
−
The specific activity (ratio of activity to molar concentration) is an exponentially decaying function S(t) = S<sub>0</sub> exp(-0.341 t) where t is expressed in units of minutes and 0.341 = ln(2)/half-life of <sup>15</sup>O and S<sub>0</sub> is the specific activity at time t = 0.
+
 
 
<pre>
 
<pre>
−
% define the initial specific activity and decay constant
+
% create empty compartmentModel object
−
sa0 = 1; % specific activity at t=0
+
cm = compartmentModel;
−
dk = 0.341; % O-15 decay constant
+
 
 +
% 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)
 
</pre>
 
</pre>
  
  
−
Next specify the input functions C<sub>p</sub> and C<sub>a</sub> 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<sub>1</sub>, p<sub>2</sub>, ... p<sub>6</sub> for C<sub>p</sub> and C<sub>a</sub>.
+
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.
 +
 
 
<pre>
 
<pre>
−
bfcm = addInput(bfcm, 'Cp', sa0, dk, 'fin', 'pCp');
+
cm = addCompartment(cm, 'Ce');
−
bfcm = addParameter(bfcm, 'pCp', 2*[851.1 21.88 20.81 -4.134 -0.1191 -0.01043 .5]);
+
cm = addCompartment(cm, 'Cm');
 +
cm = addCompartment(cm, 'Junk');
 +
</pre>
  
−
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]);
+
For the plasma concentration time-course of FDG or input function, we'll use the
 +
[http://www.ncbi.nlm.nih.gov/pubmed/8449593?ordinalpos=3&itool=EntrezSystem2.PEntrez.Pubmed.Pubmed_ResultsPanel.Pubmed_RVDocSum 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.
 +
 
 +
<pre>
 +
x0 = [2 0 500 5 5 -7 -0.1 -0.015];
 +
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
 
</pre>
 
</pre>
  
  
−
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. 
+
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<sub>e</sub> and C<sub>m</sub> compartments
 +
 
 
<pre>
 
<pre>
−
% connect compartments and inputs
+
Wlist = {...
−
bfcm = addLink(bfcm, 'L', 'Cp', 'CT', 'K1');
+
    'Ce', 'PV';
−
bfcm = addLink(bfcm, 'K', 'CT', 'J', 'k2');
+
    'Cm', 'PV'};
−
+
Xlist = {};
−
% define values for K1, k2 ...
+
cm = addOutput(cm, 'PET', Wlist, Xlist);
−
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);
 
 
</pre>
 
</pre>
  
  
−
==== Define tissue and blood outputs ====
+
Next define the scan frames as the start and stop time of each image in the dynamic sequence.
−
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<sub>T</sub>S + F<sub>v</sub>C<sub>a</sub> whereas ArterialPixel is 100% C<sub>a</sub>.
+
 
 
<pre>
 
<pre>
−
% define the activity concentration in tissue pixel
+
ttt=[ ones(12,1)*10/60; ...  % 10 sec
−
bfcm = addOutput(bfcm, 'TissuePixel', {'CT', '1'}, {'Ca', 'Fv'});
+
      ones(10,1)*0.5; ...    %  0.5 min
−
   
+
      ones(10,1)*2;...      %  2 min
−
% define the activity concentration in arterial blood pixel
+
      ones(10,1)*5;...      % 5 min
−
bfcm = addOutput(bfcm, 'ArterialPixel', {'CT', '0'}, {'Ca', '1'});
+
      ones(4,1)*10];        % 10 min
 +
scant=[[0;cumsum(ttt(1:(length(ttt)-1)))] cumsum(ttt)];
 +
cm = set(cm, 'ScanTime', scant);
 
</pre>
 
</pre>
  
−
Note that COMKAT even accounts for the time-averaging of concentrations over the durations of the scan frames.
+
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:
−
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).
 
 
<pre>
 
<pre>
−
% specify scan frame times
+
scant =
−
bfcm = set(bfcm, 'ScanTime', [[0:10:590]' [10:10:600]']/60); % division by 60 converts sec to min
+
 
 +
        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
 +
...
 
</pre>
 
</pre>
  
  
−
==== Solve for the model output ====
+
Notice thus far the links (arrows) that connect the compartments have not been defined.
−
Finally, tell COMKAT to solve the model to obtain simulation data that will be used in fitting.
+
This has been delayed since we want to make two models: one will be the standard or reference model with time-invariant rate constants
 +
and one with the time-varying K<sub>1</sub> and k<sub>2</sub>. To make two models, we now make a copy of what we have so far
 
<pre>
 
<pre>
−
[PET, PETIndex, Output, OutputIndex] = solve(bfcm);
+
cmRef = cm;
 +
</pre>
  
−
% 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
+
For the reference model, define the links as time-invariant
−
set(1,'PaperPosition',[0.25 2.5 3 2])
+
<pre>
−
print -dpng c:\temp\ExampleFigureBloodFlowInputModelDataToFit.png
+
% 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');
 +
cmRef = addLink(cmRef, 'K',        'Cm', 'Ce',  'k4');
 
</pre>
 
</pre>
  
−
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}.
 
  
−
[[Image:ExampleFigureBloodFlowInputModelDataToFit.png]]
+
For the fancy model, define the links with the time-varying K<sub>1</sub> and k<sub>2</sub> rate constants
 +
<pre>
 +
cm = addLink(cm, 'LSpecial', 'Cp', 'Ce',  'K1', 'kinLinearTime', {'K10', 'dK1dt'});
 +
cm = addLink(cm, 'KSpecial', 'Ce', 'Junk', 'k2', 'kinLinearTime', {'k20', 'dk2dt'});
 +
</pre>
 +
 
 +
noting that these make use of the kinLinearTime.m function defined above.
  
 +
Also, define the time-invariant links
 +
<pre>
 +
cm = addLink(cm, 'K',        'Ce', 'Cm',  'k3');
 +
cm = addLink(cm, 'K',        'Cm', 'Ce',  'k4');
 +
</pre>
  
−
== Fit Data ==
 
−
Modify the COMKAT model a bit so that it can be used to fit data created above.
 
  
−
==== Specify the data to fit ====
+
OK, the models are ready for use.  Begin by solving the model to obtain the time-courses of FDG activity that these models predict would be observed in PET images.
−
Just use the result from above.
 
 
<pre>
 
<pre>
−
data = PET(:, [3 4]);
+
[PET, PETIndex] = solve(cm);
−
bfcm = set(bfcm, 'ExperimentalData', data);
+
[PETRef, PETRefIndex] = solve(cmRef);
 
</pre>
 
</pre>
  
−
Specify parameters to estimate along with their initial guesses and range of feasible values.
+
 
 +
Now plot the time-courses.  Column 1 of PET (and PETRef) holds the frame start and end times.  For convenience, compute the mid-frame time
 +
<pre>
 +
tmid = (PET(:,1) + PET(:,2)) / 2;
 +
</pre>
 +
 
 +
And finally do the plotting
 +
<pre>
 +
plot(tmid, PET(:,3), 'b', tmid, PETRef(:,3), 'b:');
 +
xlabel('Time (Minutes)');
 +
ylabel('Activity (uCi/ml)');
 +
legend('Time-varying', 'Reference', 'Location', 'SouthEast');
 +
</pre>
 +
 
 +
 
 +
Well, we can get a bit fancy and plot the time-courses of the K<sub>1</sub> and k<sub>2</sub> rate constants above the model outputs and add other annotations.
 +
<pre>
 +
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');
 +
</pre>
 +
 
 +
[[Image:ExampleFigureFDGTimeVaryingK1k2.png]]
 +
 
 +
 
 +
These results are a bit surprising.  '''Why does reducing the both K<sub>1</sub> and k<sub>2</sub> 50% per hour have such a small impact on the model output?'''
 +
I speculate that the majority of the sensitivity to these parameters occurs when the plasma concentration is high.  Also, in this example, the changes in K<sub>1</sub> and k<sub>2</sub> were proportional such that the ratio was not changing.  Why does this matter?  Consider the FDG model if there were no phosphorylation (e.g. k<sub>3</sub> were zero) so that there was only a single compartment.  The differential equation for the compartment would be
 +
 
 +
dC<sub>e</sub>/dt = K<sub>1</sub> C<sub>p</sub> - k<sub>2</sub> C<sub>e</sub>
 +
 
 +
 
 +
At equilibrium, dC<sub>e</sub>/dt = 0 which implies K<sub>1</sub> C<sub>p</sub> = k<sub>2</sub> C<sub>e</sub> or
 +
C<sub>e</sub>/C<sub>p</sub> = K<sub>1</sub> / k<sub>2</sub>
 +
 
 +
This means, that the example that the equilibrium concentration ratio of tissue to plasma (in absence of phosphorylation) depends only on the ratio K<sub>1</sub> / k<sub>2</sub>.
 +
Perhaps in the example once the time gets a bit after the injection and things settle down, the model output is not so sensitive to the individual values K<sub>1</sub>  and k<sub>2</sub> but instead more dependent on their ratio.  To test this hypothesis, I can revise the above simulation to consider the case of k<sub>2</sub> decreasing 10% per hour so that the ratio K<sub>1</sub> / k<sub>2</sub> is altered over time. To do this, simply change one line:
 +
 
 +
Original
 
<pre>
 
<pre>
−
% specify parameters to estimate
+
dk2dt = -k20 * 0.5 / 60;  % 50% decrease per 60 minutes
−
bfcm = addSensitivity(bfcm, 'pCp', 'F');
+
</pre>
  
−
% define initial guess and lower and upper bounds
+
Revised
−
pCpInit = 0.9*2*[851.1 21.88 20.81 -4.134 -0.1191 -0.01043 .5];
+
<pre>
−
FInit = 0.3;
+
dk2dt = -k20 * 0.1 / 60; % 10% decrease per 60 minutes
−
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.];
 
 
</pre>
 
</pre>
 +
 +
When the simulation is performed with the new value of k<sub>2</sub>, the effect of the time-varying rate constants is more evident.
 +
 +
[[Image:ExampleFigureFDGTimeVaryingK1k2Diffk2.png]]
  
  
−
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.)
+
On the other hand, if the rate constants are changing proportionally (as modeled initially), then neglecting the time-variation while fitting data might have a minimal impact.
 +
This is examined next.
  
  
−
==== Do the actual fitting ====
+
 
−
Fit the data to estimate values of the parameters
+
==== Fitting the standard (time-invariant rate constants) model to data with changing physiology ====
 +
Lets examine the consequences of fitting the standard (time-invariant rate constants) model to data wherein the
 +
the physiology (rate constants) are changing.  Use the fancy model to create "data" and fit the standard model to that "data".
 +
 
 
<pre>
 
<pre>
−
[pEst, modelFit] = fit(bfcm, pInit, lb,ub);
+
data = PET(:,3);
 +
cmRef = set(cmRef, 'ExperimentalData', data);
 
</pre>
 
</pre>
−
COMKAT makes fitting easy!
 
  
  
−
Plot the results and make the figure below for this wiki
+
Set the initial guess and lower and upper bounds on the rate constants
 
<pre>
 
<pre>
−
plot(midTime, modelFit(:,1), 'b', midTime, PET(:,3), 'b.', midTime, modelFit(:,2), 'r', midTime, PET(:,4), 'r.')
+
pinit = [0.07; 0.1;  0.04;  0.005]; % inital guess
−
legend('Tissue Fit', PETIndex{3}, 'Arterial Fit', PETIndex{4})
+
plb =  [0.01; 0.001; 0.001; 0.0001];  % lower bound
−
xlabel('Time')
+
pub =  [0.5;  0.5;  0.2;  0.01]; % upper bound
−
ylabel('Activity')
+
cmRef = addSensitivity(cmRef, 'K10', 'k20', 'k3', 'k4');
−
print -dpng c:\temp\ExampleFigureBloodFlowInputModelFit.png
 
 
</pre>
 
</pre>
  
−
[[Image:ExampleFigureBloodFlowInputModelFit.png]]
+
Now do the fitting
 +
<pre>
 +
[pfit, modelFit] = fit(cmRef, pinit, plb, pub);
 +
</pre>
 +
 
 +
 
 +
Plot the fit over the "data"
 +
<pre>
 +
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)');
 +
legend('Time-varying Data', 'Time-invariant Fit', 'Location', 'SouthEast');
 +
</pre>
 +
 
 +
 
 +
[[Image:ExampleFigureFDGTimeVaryingK1k2Fit.png]]
 +
 
 +
 
 +
Examine estimated values of the rate constants
 +
 
 +
  Parm.  True (Init->Final)  Est      Est  vs. Init    Units
 +
    K<sub>1</sub>      0.150 -> 0.041    0.149    -0.4%            1/min
 +
    k<sub>2</sub>      0.325 -> 0.089    0.324    -0.3%            1/min
 +
    k<sub>3</sub>      0.050                  0.050    -0.9%            1/min
 +
    k<sub>4</sub>      0.007                  0.007    1.1%            1/min
 +
    k<sub>i</sub>        0.020 -> 0.020    0.020    -0.9%            1/min
 +
 
 +
It turns out that the estimated values are in close agreement with the values of the parameters at time zero (about the time of injection).
 +
Thus, in the circumstances simulated here, the model output is pretty insensitive to the significant decreases in K<sub>1</sub> and k<sub>2</sub>.
  
−
====Ta Da====
+
==== Get the full example ====
−
It looks like a pretty good fit.
+
'''Note'''  This code for this example can be found in '''fdgTimeVarying.m''' in the COMKAT examples folder.  That version has a few extra bells and whistles that allow us, for example. to use the FDG model with or without k<sub>4</sub>.

Latest revision as of 15:39, 31 March 2008

FDG Model with Time-varying Rate Constants

Overview

This example demonstrates how to implement the two-tissue compartment model with extended such that K1 and k2 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 .

ExampleFigureFDGModel.png


It can also be described by the differential equations:

dCe/dt = K1 Cp - k2 Ce - k3 Ce + k4 Cm

dCm/dt = k3 Ce - k4 Cm


Ce is the tissue concentration of FDG and Cm as the tissue concentration of metabolized FDG (FDG-6-phosphate)

Ce and Cm are interpreted as molar concentrations (Salinas, Muzic and Saidel 2007).


In the standard model. the rate constants K1, k2, k3 and k4 are assumed to be independent of time.


Here we allow K1 and k2 to be linear functions of time:

K1(t) = K10 + a t

k2(t) = k20 + b t

K10 and k20 are the initial (0 superscript) values of the rate constants and a and b are the derivatives (slopes) of the rate constant with respect to time.


To implement the time-varying rate constants, create a MATLAB function in a .m file. By making the function general, the same function can be used for K1 and k2. 

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 has a similar pattern to any function one might write to implement customized kinetic rules. Note the function is called numerous times as the model equations are solved. The value of t (time) takes on values from 0 to the end of the last frame and many values in betweeen - even values that do not correspond to the frame begin and end times.

The first part of the function retrieves values that are automatically sent by COMKAT when it calls this function. 

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

To get COMKAT to call this function, a link (connection between compartments) would be included in the model. For example 

cm = addLink(cm, 'LSpecial', 'Cp', 'Ce',   'K1', 'kinLinearTime', {'K10', 'dK1dt'});


When COMKAT calls the custom kinetic function, the last argument of the addLink command ({'K10', 'dK1dt'}) will be passed in as the pName argument. pValue will be a cell array with the value of pValue{i} set to the value of the parameter named in pName{i}. (Although you could use pxEval to obtain the value of the parameter from pName{i}, it is more efficient for COMKAT to do this internally (only once and cache the result). pSensIdx will be a cell array with the value of pSensIdx{i} being the position (index) of parameter named in pName{i} in the list of parameters to be estimated. If the parameter named in pName{i} was not included in an addSensitivity() call, then pSensIdx{i} will be an empty matrix.

The emphasis of the rest of the description of kinLinearTime.m is on the mathematical details.

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 K1 and k2.


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 500 5 5 -7 -0.1 -0.015];
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 Ce and Cm 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
...


Notice thus far the links (arrows) that connect the compartments have not been defined. This has been delayed since we want to make two models: one will be the standard or reference model with time-invariant rate constants and one with the time-varying K1 and k2. To make two models, we now make a copy of what we have so far 

cmRef = cm;


For the reference model, define the links as time-invariant 

% 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');
cmRef = addLink(cmRef, 'K',        'Cm', 'Ce',   'k4');


For the fancy model, define the links with the time-varying K1 and k2 rate constants 

cm = addLink(cm, 'LSpecial', 'Cp', 'Ce',   'K1', 'kinLinearTime', {'K10', 'dK1dt'});
cm = addLink(cm, 'KSpecial', 'Ce', 'Junk', 'k2', 'kinLinearTime', {'k20', 'dk2dt'});

noting that these make use of the kinLinearTime.m function defined above.

Also, define the time-invariant links 

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


OK, the models are ready for use. Begin by solving the model to obtain the time-courses of FDG activity that these models predict would be observed in PET images. 

[PET, PETIndex] = solve(cm);
[PETRef, PETRefIndex] = solve(cmRef);


Now plot the time-courses. Column 1 of PET (and PETRef) holds the frame start and end times. For convenience, compute the mid-frame time 

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

And finally do the plotting 

plot(tmid, PET(:,3), 'b', tmid, PETRef(:,3), 'b:');
xlabel('Time (Minutes)');
ylabel('Activity (uCi/ml)');
legend('Time-varying', 'Reference', 'Location', 'SouthEast');


Well, we can get a bit fancy and plot the time-courses of the K1 and k2 rate constants above the model outputs and add other annotations. 

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');

ExampleFigureFDGTimeVaryingK1k2.png


These results are a bit surprising. Why does reducing the both K1 and k2 50% per hour have such a small impact on the model output? I speculate that the majority of the sensitivity to these parameters occurs when the plasma concentration is high. Also, in this example, the changes in K1 and k2 were proportional such that the ratio was not changing. Why does this matter? Consider the FDG model if there were no phosphorylation (e.g. k3 were zero) so that there was only a single compartment. The differential equation for the compartment would be

dCe/dt = K1 Cp - k2 Ce


At equilibrium, dCe/dt = 0 which implies K1 Cp = k2 Ce or Ce/Cp = K1 / k2

This means, that the example that the equilibrium concentration ratio of tissue to plasma (in absence of phosphorylation) depends only on the ratio K1 / k2. Perhaps in the example once the time gets a bit after the injection and things settle down, the model output is not so sensitive to the individual values K1 and k2 but instead more dependent on their ratio. To test this hypothesis, I can revise the above simulation to consider the case of k2 decreasing 10% per hour so that the ratio K1 / k2 is altered over time. To do this, simply change one line:

Original 

dk2dt = -k20 * 0.5 / 60;  % 50% decrease per 60 minutes

Revised 

dk2dt = -k20 * 0.1 / 60;  % 10% decrease per 60 minutes

When the simulation is performed with the new value of k2, the effect of the time-varying rate constants is more evident.

ExampleFigureFDGTimeVaryingK1k2Diffk2.png


On the other hand, if the rate constants are changing proportionally (as modeled initially), then neglecting the time-variation while fitting data might have a minimal impact. This is examined next.


Fitting the standard (time-invariant rate constants) model to data with changing physiology

Lets examine the consequences of fitting the standard (time-invariant rate constants) model to data wherein the the physiology (rate constants) are changing. Use the fancy model to create "data" and fit the standard model to that "data". 

data = PET(:,3);
cmRef = set(cmRef, 'ExperimentalData', data);


Set the initial guess and lower and upper bounds on the rate constants 

pinit = [0.07; 0.1;   0.04;  0.005]; % inital guess
plb =   [0.01; 0.001; 0.001; 0.0001];  % lower bound
pub =   [0.5;  0.5;   0.2;   0.01]; % upper bound
cmRef = addSensitivity(cmRef, 'K10', 'k20', 'k3', 'k4');

Now do the fitting 

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


Plot the fit over the "data" 

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)');
legend('Time-varying Data', 'Time-invariant Fit', 'Location', 'SouthEast');


ExampleFigureFDGTimeVaryingK1k2Fit.png


Examine estimated values of the rate constants 

  Parm.   True (Init->Final)  Est       Est  vs. Init     Units
   K1       0.150 -> 0.041     0.149    -0.4%            1/min
   k2       0.325 -> 0.089     0.324    -0.3%            1/min
   k3       0.050                   0.050    -0.9%            1/min
   k4       0.007                   0.007    1.1%             1/min
   ki        0.020 -> 0.020     0.020    -0.9%            1/min

It turns out that the estimated values are in close agreement with the values of the parameters at time zero (about the time of injection). Thus, in the circumstances simulated here, the model output is pretty insensitive to the significant decreases in K1 and k2.

Get the full example

Note This code for this example can be found in fdgTimeVarying.m in the COMKAT examples folder. That version has a few extra bells and whistles that allow us, for example. to use the FDG model with or without k4.

Retrieved from "https://comkat.cwru.edu/index.php?title=Support:Documents:Examples:FDG_with_Time-varying_Rate_Constants&oldid=583"