Support:Documents:Manual:COMKAT Command Line Interface

COMKAT Overview

Summary

COMKAT allows a user to construct a compartment model of his/her own design without resorting to writing equations. Models may then be used for a simulation and fitting data. COMKAT uses an object-oriented representation of models which is illustrated in this article.

Background

COMKAT enables a user to specify a model configuration without having to resort to writing and programming differential equations. This is all handled behind the scenes using state-of-the art differential equations solvers. On several platforms a solver written in C is used and transparently accessed from MATLAB via a mexfile interface created by the COMKAT authors. This means that solving the equations will be efficient and reliable. On platforms for which the compiled solver is not available, COMKAT will automatically use a differential equation solver (e.g. ode15s) from the MATLAB ODE Suite. This too will be reliable but may take more computer time.

Interfaces

There are two interfaces to COMKAT: a graphical user interface (GUI) and the command line interface (CLI). While the GUI has snazzy appeal and easy of use, the CLI is of particular utility to people doing modeling research. This is because the CLI enables one to set up simulations that loop and create thousands of simulated data sets. Because this is all done within the MATLAB environment, results of such simulations may be analyzed, graphed, and otherwise be prepared for publication. Alternatively, results could be saved to EXCEL spreadsheet or pasted into documents.

Principles

Modeling generally entails analyzing experimental data using mathematical rules that predict relationships between input functions, tissue responses and physiologic parameters. Modeling is used in two forms: 1 simulation or prediction and 2 fitting or characterization of function. In practice, to do fitting entails repeating simulations and adjusting values of parameters until the model-predicted output matches measured data. Consequently, 2 can be thought of as a superset of 1.

CompartmentModel Object

COMKAT uses an object-oriented approach to represent models. A compartmentModel object can be thought of as the type for a variable that holds all the information needed to describe the kinetics, to setup and solve the equations, and even to fit data. Special functions called methods are used to configure the model, perform the solving and do the fitting. What makes something a method and not just a function is that a method acts on information stored in an object. While this might sound abstract and obtuse, hopefully the following example will convince you that this is a convenient, intuitive approach.

compartmentModel Object

Using the MATLAB command line interface (CLI), type this to create a new model object:

cm = compartmentModel;

Next use the whos command to see what was created:

>> whos cm
Name Size Bytes Class
cm     1x1 18830 compartmentModel object
Grand total is 428 elements using 18830 bytes

Notice that cm is not a double array, not an integer array and not even a structure. cm is a special data type created using the special function compartmentModel. This function is called a constructor because it creates new instance of an object. Actually, the function is a method. It is the constructor method for a compartmentModel object. (MATLAB "knows" that is a constructor since the name of the method, compartmentModel is the same as the name of the object type.)

At this point cm is a model with no compartments, no links connecting compartments, no data... as is revealed when examining its contents

>> cm
untitled1
Inputs: Name Spec. Act. Decay Const.
Compartments: Name (Saturation)
Outputs 
Sensitivities
Links
Parameters: Name=Value

Note that nothing is special about the variable name cm. One could have just as easily created a compartmentModel object and called it Fred (Fred = compartmentModel;).

Compartments

Not much can be demonstrated with an empty model so we might as well create something useful. We start by adding a couple of compartments. Lets create a model that describes a hypothetical model of blood flow to the brain. The model will have two compartments one representing the brain and one representing the rest of the body.

>> brainModel = compartmentModel; % create a new, empty model
>> brainModel = addCompartment(brainModel, 'Brain');
>> brainModel = addCompartment(brainModel, 'Body');

[For those not familiar with MATLAB, % create a new, empty model is a comment. It is something that MATLAB ignores and is only for the convenience of the person reading the commands.]

Inputs

Lets assume that we are doing an imaging experiment and that something like radioactive water or gadolinium is administered intravenously to the subject. Blood in the carotid arteries will carry the material to the brain. Hence we consider the carotids as input function. That is, the time-course of the concentration of such a material in the blood would be experimentally measured. For the sake of simplicity lets assume that a MATLAB function called fengInput is created which can provide the blood concentration at time t. [In a more experimentally realistic case the concentration of the input function at time t may be determined by interpolating (e.g. linear interpolation) measured data. This can be done too and is, for the sake of keeping the present example simple, shown elsewhere.]

fengInput.m calculates the blood concentration at time t as a sum of gamma variate and exponential terms. The amplitudes and the exponential constants are all specified in input p to fengInput. For the sake of example, let p be defined as

>> p=[2 0 851.1 21.88 20.81 -4.134 -0.1191 -0.01043];

[The details of what each element of p means is not important to the current discussion.]

The following commands plot the input function concentration vs time [0., 0.1; 0.2, .... 100.0] minutes

>> p=[2 .5 16000 2188 2081 -7.5 -10 -0.01043];
>> t = [0:0.1:100];
>> plot(t,fengInput(p,t))

To include this input function in the model, use this command

>> brainModel = addInput(brainModel, 'Carotid', 1, 0, 'fengInput', p);

[For now don't worry about the 1 and 0 in the command. These are the specific activity (e.g. Bq/pmol) and tracer halflife (e.g. minutes). Because of the values used in this example the model can be interpreted as using a non-radioactive tracer or a tracer with infinite half life.]

Links

At this point the model has two compartments and one input but there is no specification as how tracer moves between the input and compartments. Herein we assume compartments are connected as depicted in the figure

To create the k1 link from the corotid input as the source to the brain as the destination, use this command

>> brainModel = addLink(brainModel, 'L', 'Carotid', 'Brain', 'k1');

[COMKAT is case-sensitive. Since the input was called Carotid in addInput(), it must be called Carotid in addLink() and not carotid and not CaRoTiD... You may use whatever mixture of upper and lower case that you prefer so long as you are consistent!]

The parameter 'L' in the statement identifies the link type. 'L' means the link source is an input and the kinetics are linear. To include the k2 link in the model the same type of command is used except that the linktype is 'K' because the link source is a compartment and not an input. [A list of the possible link types is given in the help for addLink.]

>> brainModel = addLink(brainModel, 'K', 'Brain', 'Body', 'k2');

For the mathematically inclined readers, we have just set up a model described by these equations:

 dBrain/dt = k1 • Corotid – k2 • Brain
 dBody/dt = k2 • Brain

Parameters

The next step is to assign values for the rate constants k1 and k2. These can be adjusted later so the values used here are not so critical. For now just assign k1 = 0.8, k2 = 0.5.

>> brainModel=addParameter(brainModel, 'k1', 0.8);
>> brainModel=addParameter(brainModel, 'k2', 0.5);

Model Output

At this point we are almost at the point where the model equations can be solved. What is missing is a specification of what is the model output. The output is not restricted to being the concentration in a single compartment but rather can be a weighted sum of the concentration in more than one compartment. Also there is the issue of defining the times at which solutions to the equations are desired. Model outputs are often defined as corresponding to the measured data. For example, pixel values may correspond to brain concentration averaged over the duration of the image acquisition.

Ignoring this time-averaging for a moment, an output might be thought of as one times the brain compartment concentration plus zero times the body compartment concentration. In COMKAT only the non-zero contributions of a compartment to outputs need be specified. Here is an example of how it might be done.

>> wlist = {'Brain', 1};
>> xlist = {};
>> brainModel=addOutput(brainModel, 'pixels', wlist, xlist);

The wlist defines the compartments and weights with which they contribut to outputs. (E.g. to define an output as 0.5*Brain + 0.4*Body, wlist = {'Brain', 0.5; 'Body', 0.4 }; xlist is similar to wlist except that xlist defines how input functions contribute directly to outputs. If for example, one would like the model the pixel as having 4% contribution from vascular space, do this xlist = {'Carotid', 0.04 };.)

To account for time-averaging inherent in data collection processes, we need only specify the the start and stop time of each image in the dynamic sequence. Lets assume the are 14 images and the start time (column 1) and end times (column 2 ) of each image is given in a MATLAB variable

>> st = [...
0    0.5
0.5  1
1    2
2    3
3    4
4    5
5   10
10  15
15  20
20  25
25  30
30  40
40  50
50  60
];

For example, image 3 (row 3 in st) begins at 1 minutes and ends at 2 minutes. To tell specify this in the model description, use the following command

>> brainModel = set(brainModel, 'ScanTime', st);

Solving the Model

At this point the model is completely defined and can be solved to obtain model-predicted concentration curves using the solve method:

>> [result, resultIndex] = solve(brainModel);

That is all there is to it. To see what we got, use the whos command

>> whos result
Name Size  Bytes Class
result 14x3 336    double array

The array result has as many rows as there images in the dyanmic image sequence. In this simple case the first column holds the scan start times, the second holds the scan end times and the third holds the output which was called 'pixels'. The number of columns and their contents depends on the model specification and all of this is described resultIndex

>> resultIndex
resultIndex =
'Scan Start'
'Scan End'
'pixels'

In general resultIndex{i} tells what is in column i of result.

Standard MATLAB commands may be used to plot the pixel concentration versus the mid-time of each scan.

>> t=0.5*(result(:,1)+result(:,2)); % t is the mid-scan time
>> plot(t,result(:,3));
>> xlabel('MInutes')
>> ylabel('Concentration')

Closing Remarks

COMKAT allows a great deal of flexibility in defining, solving models and fitting models to data. This article just scratches the surface of the possiblities. This article focused on the command line interface (CLI) of COMKAT. Another article provides an introduction to the graphical user interface (GUI). Importantly, the GUI is just a "front end" that – behind the scenes – uses CLI.

COMKAT Command Line Function Reference

Command	Syntax	Comment
compartmentModel	cm = compartmentModel;	Creates a new compartment model object
addCompartment	cm2 = addCompartment(cm1, compartmentName[ ,saturationParameter])	Adds a compartment to a model
addLink	cm2 = addLink(cm1, linkType, linkSource, linkDestination, linkName)	Adds a connection between two compartments or an input and a compartment
addInput	cm2 = addInput(cm1, inputName, specificActivity, decayConstant, functionName, inputParameter [, extraParameter])	Adds an input to a model
addOutput	cm2 = addOutput(cm1, outputName, wSpec, xSpec)	Adds an output to a model (multiple outputs are supported)
addParameter	cm2 = addParameter(cm1, parameterName, parameterValue[, parameterMin, parameterMax])	Adds a parameter to a model and specified initial value and optional limits on range
addSensitivity	cm2 = addSensitivity(cm1, parameterName)	Tells COMKAT that sensitivty functions are to be solved for the specified parameter(s). These parameters will be adjusted when fitting.
analyzeFit	[standardDeviation, covariance, correlation] = analyzeFit(cm1, jacobian, parameterFitValues, residuals)	Estimates precision of parameter estimates
cIdx	index = cIdx(cm1, compartmentName)	Returns compartment number associated with named compartment
dependsOn	yesNo = dependsOn(cm1, parmeter1, parameter2)	Returns 1 of parameter1 depends on parameter 2, 0 otherwise. parameter1 may be cell string array
exportAsScript	status = exportAsScript(cm1, fileName, variableName)	Creates a new file that contains COMKAT commands for creating the model in cm1. In the commands the model is put into the specified variable name.
fit	[newParmVal,modelfit,resnorm,residual,exitFlag,output,lambda,jacobian] = fit(cm, initParmVal, lb, ub)	Fit model to data returning parameter estimates and other optional information
fitGen	[pfit, qfit, modelfit, exitflag, output, lambda, grad, hessian] = fitGen(cm, pinit, plb, pub)	Generalized fitting. Supports a number of different objective functions.
get	value = (cm1, property)	Retrieve information from compartment model object
getInputType	inputType = getInputType(cm, inputName)	Returns information about whether named input is used as a plasma concentration or as whole blood activity input
iIdx	index = iIdx(cm, inputName)	Retuns index number of named input.
load	cm2 = load(cm1, fileName)	Loads into memory a compartment model object that is stored in the specified mat file
pDeriv	dp1dp2 =pDeriv(cm, p1, p2)	Evaluates the partial derivative of parameter p1 with respect to parameter p2
removeCompartment	cm2 = removeCompartment(cm1, idx)	Makes a copy of a model with compartment (specified by index) deleted
removeInput	cm2 = removeInput(cm1, input)	Returns a copy of the model with input deleted. Input may be specified by name or index
removeLink	cm2 = removeLink(cm1, idx)	Returns a copy of the model with link deleted. Link may be specified by name or index
save	cm2 = save(cm1, fileName)	Save specified compartment model object to the named mat file
set	cm2 = set(cm1, property, value[, one or more optional inputs])	Set value of specified model property
solve	[PET, PETindex, component, componentIndex] = solve(cm)	Solve model equations

Model construction basics

Summary

COMKAT allows a user to construct a model of his/her own design without resorting to writing equations.

Background

COMKAT enables a user to specify a model configuration without having to resort to writing and solving differential equations. This is all handled behind the scenes using state-of-the art differential equations solvers. On several platforms, the solver is written in C and accessed from MATLAB via a mexfile interface. This means that solving the equations will be efficient and reliable.

There are two interfaces to COMKAT: a graphical user interface (GUI) and the command line interface (CLI). While the GUI has snazzy appeal and easy of use, the CLI is the workhorse for people doing modeling research. This is because the CLI enables one to setup simulations that loop and create thousands of simulated data sets all within the MATLAB programming environment.

Principles

Modeling generally entails analyzing experimental data using mathematical rules that predict relationships between input functions, tissue responses and physiologic parameters. Values of physiologic parameters, such as perfusion, are estimated by adjusting their values until the model-predicted output agrees with the experimental data. Input function fitting is no different. One proposes a mathematical form as describing the input and then adjusts values of its parameters to achieve the best match between model-predicted and experimental data.

COMKAT Implementation

In this section an example is presented of using COMKAT to estimate parameters of the FDG input function. The first form of the example is intentionally kept simple to best illustrate concepts. The second example is more complex and perhaps a good starting point for you to research strategies for simultaneously estimating input functions and tissue kinetic parameters.

Example 1: Estimate Input Parameters, Rate Constants Known

Here we assume that the input function has the form

Cp = (p1t - p2 - p3) exp(-p4t) + p2 exp(-p5t) + p3 exp(-p6t)

where Cp is the plasma concentration of 18FDG, pi are the input function parameters to be estimated, and t is time. (This form comes from Feng, Huang, Wang. "Models for Computer simulation Studies of Input Functions for Tracer Kinetic Modeling with Positron Emission Tomography", International Journal of Biomedical Computing, 32(2):95-110, March 1993.)

Step 1 To set this up in MATLAB for simulation, create the following function and put it in a file called refCp.m:

function [fval, jac] = refCp(p, t)
 if (nargout > 0),
       t = t(:);  % make it a column vector
       fval = (p(1)*t - p(2) - p(3)) .* exp(-p(4)*t) + p(2) * exp(-p(5)*t) + p(3) * exp(-p(6)*t);
       if (nargout > 1),
           jac = [t .* exp(-p(4)*t) ...
                  -exp(-p(4)*t) + exp(-p(5)*t) ...
                  -exp(-p(4)*t) + exp(-p(6)*t) ...
                  -t .* (p(1)*t - p(2) - p(3)) .* exp(-p(4)*t) ...
                  -t * p(2) .* exp(-p(5)*t) ...
                  -t * p(3) .* exp(-p(6)*t)];
       end
   end
return

This function can be used to plot an example input curve using these commands

   t=0:.1:60;
   pin = [28 0.75 0.70 4.134 0.1191 0.01043];
   plot(t,refCp(pin,t))
   set(gca,'FontSize',6);
   axis([0 60 0 5]);
   xlabel('Minutes');
   ylabel('pmol/ml');

Note that he above function, refCp.m not only calculates values Cp, it also optionally calculates values for the derivative of Cp with respect to parameter vector p. These will be needed later since to fit the input function, it is important to quantify how changes in values of p will effects Cp values which, in turn, effect model-predicted tissue concentration reflected in voxels or regions of interest.

Step 2 Create an FDG model

   cm = compartmentModel;

   %        k1     k2      k3     k4
   ktrue=[0.1020 0.1300 0.0620 0.0068];

   % define the parameters
   cm = addParameter(cm, 'k1',    ktrue(1));     % 1/min
   cm = addParameter(cm, 'k2',    ktrue(2));     % 1/min
   cm = addParameter(cm, 'k3',    ktrue(3));     % ml/(pmol*min)
   cm = addParameter(cm, 'k4',    ktrue(4));     % 1/min
   cm = addParameter(cm, 'sa',    75);           % specific activity of injection, kBq/pmol
   cm = addParameter(cm, 'dk',    log(2)/109.8); % radioactive decay
   cm = addParameter(cm, 'PV',    1);            % (none)

   % define input function parameter vector
   cm = addParameter(cm, 'pin', [28; 0.75; 0.70; 4.134; 0.1191; 0.01043]);

   % define compartments
   cm = addCompartment(cm, 'Ce');
   cm = addCompartment(cm, 'Cm');
   cm = addCompartment(cm, 'Junk');

   % define plasma input function
   % specifying function as refCp with parameters pin
   cm = addInput(cm, 'Cp', 'sa', 'dk', 'refCp', 'pin'); % plamsa pmol/ml

   % connect inputs and compartments
   cm = addLink(cm, 'L', 'Cp', 'Ce', 'k1');
   cm = addLink(cm, 'K', 'Ce', 'Junk', 'k2');
   cm = addLink(cm, 'K', 'Ce', 'Cm', 'k3');
   cm = addLink(cm, 'K', 'Cm', 'Ce', 'k4');

   % specify scan begin and end times
   ttt=[ ones(6,1)*5/60; ...    %  6 frames x  5   sec
         ones(2,1)*15/60; ...   %  2 frames x 15   sec
         ones(6,1)*0.5;...      %  6 frames x  0.5 min
         ones(3,1)*2;...        %  3 frames x  2   min
         ones(2,1)*5;...        %  2 frames x  5   min
         ones(10,1)*10];        % 10 frames x 10   min
   scant = [[0;cumsum(ttt(1:(length(ttt)-1)))] cumsum(ttt)];
   cm = set(cm, 'ScanTime', scant);
 
   % define an output which is sum of Ce and Cm adjusted to account 
   % for partial volume (PV).  Ignore vascular activity.
   Wlist = {...
       'Ce', 'PV';
       'Cm', 'PV'};
   Xlist = {};
   cm = addOutput(cm, 'PET', Wlist, Xlist);

   % solve model and generate example output
   [PET, PETindex]=solve(cm);
   plot(0.5*(PET(:,1)+PET(:,2)),PET(:,3), '.');
   set(gca,'FontSize',6);
   axis([0 60 0 100]);
   xlabel('Minutes');
   ylabel('KBq/ml');

Step 3 Use model output as perfect "data"

   data = PET(:,3);

Step 4 Fit the "data"

   % tell model abbout data to be fit
   cm = set(cm, 'ExperimentalData', data);

   % specify parameters to be adjusted in fitting
   cm = addSensitivity(cm, 'pin');

   % set parameter values initial guess, lower and upper bounds
   pinit = [ 10; 0.4;  0.4;  3;  0.05; 0.01];
   plb =   [ 10; 0.1;  0.1;  1;  0.05; 0.001];
   pub =   [100; 2. ;  2. ; 10;  1.  ; 0.05];

   % actually do the fitting
   pfit = fit(cm, pinit, plb, pub);

Values of parameter pfit estimates are [28.0167; 0.7652; 0.7059; 4.1586; 0.1257; 0.0105] and are reasonably close to the true values (pxEval(cm, 'pin')) [28; 0.75; 0.70; 4.134; 0.1191; 0.01043]

Step 5 Examine model outputs for "data", initial guess, and fit

   PETfit = solve(set(cm, 'ParameterValue', 'pin', pfit));
   PETinit = solve(set(cm, 'ParameterValue', 'pin', pinit));
   t = 0.5*(PET(:,1)+PET(:,2));
   plot(t,data,'.', t, PETinit(:,3), t, PETfit(:,3));
   set(gca,'FontSize',8);
   axis([0 60 0 75]);
   xlabel('Minutes');
   ylabel('KBq/ml');
   ah=legend('Data', 'Initial Guess', 'Fit','Location','SouthEast');
   set(ah,'FontSize',6,'Box','off');

Example 2: Simultaneously Estimate Input Parameters and Rate Constants

Aside from the obvious difference of estimatting values for rate constants as well as input parameters, this example simultaneoulsy fits three tissue regions. To do this, three different outputs (regA, regB, regC) are used and these are based on three different sets of compartments. We assume a single, common input.

Another important detail is that because of the way the K1 parameters appear in the differential equations - always multiplied by Cp, there is insufficient information to estimate both the K1 values and the scale of the input function. Think about it. It is what is called an identifiability problem. If you double K1 and halve Cp, their product - and hence model output - is unchange. Consequently, here we take the approach of fixing K1 values at their true value.

Step 1 Do nothing, just use the input from step 1 above.

Step 2 Create an FDG model for three regions all sharing the same input. (This example is a lot of typing or just one giant Copy and Paste. Nevertheless, it should be quite straightforward if you understood Example 1.

   cm = compartmentModel;  % start with a new, empty model

   %        k1     k2      k3     k4
   ktrueA=[0.1 ;  0.13 ; 0.06 ; 0.0068];
   ktrueB=[0.2 ;  0.25 ; 0.04 ; 0.0068];
   ktrueC=[0.05;  0.13 ; 0.02 ; 0.0068];

   % define the parameters
   cm = addParameter(cm, 'sa',    75);           % specific activity of injection, kBq/pmol
   cm = addParameter(cm, 'dk',    log(2)/109.8); % radioactive decay
   cm = addParameter(cm, 'PV',    1);            % (none)

   % region A
   cm = addParameter(cm, 'k1A',    ktrueA(1));     % 1/min
   cm = addParameter(cm, 'k2A',    ktrueA(2));     % 1/min
   cm = addParameter(cm, 'k3A',    ktrueA(3));     % ml/(pmol*min)
   cm = addParameter(cm, 'k4A',    ktrueA(4));     % 1/min

   % region B
   cm = addParameter(cm, 'k1B',    ktrueB(1));     % 1/min
   cm = addParameter(cm, 'k2B',    ktrueB(2));     % 1/min
   cm = addParameter(cm, 'k3B',    ktrueB(3));     % ml/(pmol*min)
   cm = addParameter(cm, 'k4B',    ktrueB(4));     % 1/min

   % region C
   cm = addParameter(cm, 'k1C',    ktrueC(1));     % 1/min
   cm = addParameter(cm, 'k2C',    ktrueC(2));     % 1/min
   cm = addParameter(cm, 'k3C',    ktrueC(3));     % ml/(pmol*min)
   cm = addParameter(cm, 'k4C',    ktrueC(4));     % 1/min

   % define input function parameter vector
   cm = addParameter(cm, 'pin', [28; 0.75; 0.70; 4.134; 0.1191; 0.01043]);

   % define compartments for regions A, B, C (they can share the same junk or sink compartment)
   cm = addCompartment(cm, 'Junk');
   cm = addCompartment(cm, 'CeA');
   cm = addCompartment(cm, 'CmA');
   cm = addCompartment(cm, 'CeB');
   cm = addCompartment(cm, 'CmB');
   cm = addCompartment(cm, 'CeC');
   cm = addCompartment(cm, 'CmC');

   % define plasma input function
   % specifying function as refCp with parameters pin
   cm = addInput(cm, 'Cp', 'sa', 'dk', 'refCp', 'pin'); % plamsa pmol/ml

   % connect inputs and compartments
   % region A
   cm = addLink(cm, 'L', 'Cp',  'CeA', 'k1A');
   cm = addLink(cm, 'K', 'CeA', 'Junk','k2A');
   cm = addLink(cm, 'K', 'CeA', 'CmA', 'k3A');
   cm = addLink(cm, 'K', 'CmA', 'CeA', 'k4A');

   % region B
   cm = addLink(cm, 'L', 'Cp',  'CeB', 'k1B');
   cm = addLink(cm, 'K', 'CeB', 'Junk','k2B');
   cm = addLink(cm, 'K', 'CeB', 'CmB', 'k3B');
   cm = addLink(cm, 'K', 'CmB', 'CeB', 'k4B');

   % region C
   cm = addLink(cm, 'L', 'Cp',  'CeC', 'k1C');
   cm = addLink(cm, 'K', 'CeC', 'Junk','k2C');
   cm = addLink(cm, 'K', 'CeC', 'CmC', 'k3C');
   cm = addLink(cm, 'K', 'CmC', 'CeC', 'k4C');

   % specify scan begin and end times
   ttt=[ ones(6,1)*5/60; ...    %  6 frames x  5   sec
         ones(2,1)*15/60; ...   %  2 frames x 15   sec
         ones(6,1)*0.5;...      %  6 frames x  0.5 min
         ones(3,1)*2;...        %  3 frames x  2   min
         ones(2,1)*5;...        %  2 frames x  5   min
         ones(10,1)*10];        % 10 frames x 10   min
   scant = [[0;cumsum(ttt(1:(length(ttt)-1)))] cumsum(ttt)];
   cm = set(cm, 'ScanTime', scant);

   % define an outputs, one for each region
   cm = addOutput(cm, 'RegA', {'CeA', 'PV'; 'CmA', 'PV'}, {});
   cm = addOutput(cm, 'RegB', {'CeB', 'PV'; 'CmB', 'PV'}, {});
   cm = addOutput(cm, 'RegC', {'CeC', 'PV'; 'CmC', 'PV'}, {});

   % solve model and generate example output
   [PET, PETindex]=solve(cm);
   idx = [3 4 5];
   plot(0.5*(PET(:,1)+PET(:,2)), PET(:,idx), '.');  
   set(gca,'FontSize',6);
   xlabel('Minutes');
   ylabel('KBq/ml');
   axis([0 60 0 120]);
   ah=legend(PETindex(idx));
   set(ah,'FontSize',5,'Box','off');

Step 3 Use model output as perfect "data"

   data = PET(:,idx);  % data will have 3 columns, one for each region

Step 4 Fit the "data"

   % tell model abbout data to be fit
   cm = set(cm, 'ExperimentalData', data);

   % specify parameters to be adjusted in fitting
   cm = addSensitivity(cm, 'pin', 'k2A', 'k3A', 'k4A', 'k2B', 'k3B', 'k4B',  'k2C', 'k3C', 'k4C');

   % set parameter values initial guess, lower and upper bounds.  values are in same order as ensitivities
   %        _____________pin_________________  ______RegA______    ______RegB______  ______RegC______
   pinit = [ 10; 0.4;  0.4;  3;  0.05; 0.01;   0.1;  0.03; 0.001; 0.1;  0.03; 0.001; 0.1;  0.03; 0.001];
   plb =   [ 10; 0.1;  0.1;  1;  0.05; 0.001;  1e-3; 1e-3; 1e-3 ; 1e-3; 1e-3; 1e-3;  1e-3; 1e-3; 1e-3];
   pub =   [100; 2. ;  2. ; 10;  1.  ; 0.05;   1.;   1.;   1.;    1.;   1.;   1.;    1.;   1.;   1.];

   % actually do the fitting
   pfit = fit(cm, pinit, plb, pub);

Here are the results

         True       Fit
p1      28.0000  28.1555
p2      0.7500    0.9002
p3      0.7000    0.5201
p4      4.1340    4.1324
p5      0.1191    0.1240
p6      0.0104    0.0072
k2A    0.1300    0.1141
k3A    0.0600    0.0739
k4A    0.0068    0.0114
k2B    0.2500    0.2316
K3B    0.0400    0.0502
k4B    0.0068    0.0109
k2C    0.1300    0.1122
k3C    0.0200    0.0236
k4C    0.0068    0.0125

and, yes that is right, we just estimated values of 15 parameters. Wasn't that easy. Now lets see how well we did. Some or the parameters are estimated well; others are not.

Lets see how well well the curves fit the "data". For convenience, make a copy of the model setting the parameter to estimated values.

   cmfit = cm; % make a copy
   cmfit = set(cmfit, 'ParameterValue', 'pin', pfit(1:6));
   cmfit = set(cmfit, 'ParameterValue', 'k2A', pfit(7));
   cmfit = set(cmfit, 'ParameterValue', 'k3A', pfit(8));
   cmfit = set(cmfit, 'ParameterValue', 'k4A', pfit(9));
   cmfit = set(cmfit, 'ParameterValue', 'k2B', pfit(10));
   cmfit = set(cmfit, 'ParameterValue', 'k3B', pfit(11));
   cmfit = set(cmfit, 'ParameterValue', 'k4B', pfit(12));
   cmfit = set(cmfit, 'ParameterValue', 'k2C', pfit(13));
   cmfit = set(cmfit, 'ParameterValue', 'k3C', pfit(14));
   cmfit = set(cmfit, 'ParameterValue', 'k4C', pfit(15));
   set(gca,'FontSize', 6);
   PETfit = solve(cmfit);
   plot(0.5*(PET(:,1)+PET(:,2)),PET(:,idx),'.',0.5*(PET(:,1)+PET(:,2)),PETfit(:,idx));
   axis([0 60 0 120]);
   ylabel('KBq/ml');
   xlabel('Minutes');
   ah=legend('Data A', 'Data B', 'Data C', 'Fit A', 'FIt B', 'Fit C','Location','EastOutside');
   set(ah,'FontSize',5,'Box','off');

The fits look pretty good so it appears the imperfect parameter estimates are a consequence of identifiability problems. This is described in another lesson.