ThreeEXPdecay Notes:
The data set is EXPdecay from file EXPdecay.cdata. Data are
at 2 second intervals. Parameters for solution were: k1=1,
k2=0.2, k3= 0.04, C10=0.5, C20=0.3, C30 = 0.2, preserved
with 6, 4, 3, 2, and 1 digit accuracy
Results:
Opt2Exp: This parameter set is set up to provide a two
exponential fit to the data generated with three exponential
curves.
Trial with 2 exponentials only:
Two exponentials fit the data very well with the following
parameters: C10 = 0.5773, k1 = 0.8747, C20 = 0.4226,
K2 = 0.0895, and an RMS error = 0.001.
This raises two questions:
a. How can two exponentials fit a three exponential curve?
Using loops, run the original parameters, listed above
and the best fitting two exponential fit. They overlap
to better than 2 digits.
b. Consider Akaike information criterion* for model
minimization. As a descriptor of these data, the two
exponential fit is as good as three, therefore one should
use only two, the criterion suggests.
However if there is other information about the
experiment, e.g. there were 3 parallel paths used in
obtaining the data, then it is clear that the parameters
provided by the two exponential model cannot represent the
system correctly, i.e. the parameterization is misleading
with respect to the real system, and they are merely
descriptors, not explanatory.
Conclusions:
Other information about the system is needed to
distinguish two from three exponentials, at least under
the circumstances of the starting conditions in Opt2Exp.
If this is true for 6 digit accuracy, then reduced
numbers of digits will make no difference. (Reducing the
number of digits reduces the accuracy of the fit and cannot
improve the estimation of the parameters.)
Exercises:
Obtain data for longer times. Will this help?
Run the solutions for 30 seconds instead of 10 seconds.
The result shows that the separation between the two
solutions begins soon after 10 seconds, and the divergence
is large.
Conclusion:
Collect data long enough to be sure that slow
exponential components are observed.
The reader should try further tests. Add noisy data to the
file, e.g. for additive Gaussian noise with SD = 0.1,
real Cnoise(t) mM;
Cnoise = Canalytic+0.1*randomg();
Change the temporal resolution. Make unevenly space points,
etc.
WARNING: There are only six data points. The number of
degrees of freedom equals the number of data points minus
the number of parameters being optimized. The degrees of
freedom becomes a denominator in calculating the RMS error.
Hence, you can only optimize 5 parameters.
However, there is an additional constraint. Compare
optimizing five parameters when the fifth parameter is
either k3 or C30. The k3 case causes a problem when
C30 is fixed as 0.0 because it introduces a column and row
of zeroes in the calculation of the sensitivity matrix and
the matrix cannot be inverted to determine the confidence
limits. Optimizing k3 with C30 fixed at 1e-8 will produce
unacceptable confidence limits, especially on k3.
* A short article explaining Akaike Information Criterion
is available at Wikipedia.
/* MODEL NUMBER: 0276
MODEL NAME: ThreeEXPdecay
SHORT DESCRIPTION: Washout curve simulation by sum of three decaying exponentials:
analysis using two or three exponentials.
*/
import nsrunit; unit conversion on;
math ThreeEXPdecay {
// INDEPENDENT VARIABLE
realDomain t s; t.min=0; t.max=10; t.delta=0.1;
// PARAMETERS
real C10 = 0.5 mM, // Initial value for C1
C20 = 0.3 mM, // Initial value for C2
C30 = 0.2 mM, // initial value for C3
k1 = 1 1/s, // Decay rate for C1
k2 = 0.2 1/s, // Decay rate for C2
k3 = 0.04 1/s; // Decay rate for C3
/* COMMENT: To actually calculate washout curves, choose a value for a
flow, F, then calculate Vi=F/ki, where i=1,2, and 3. */
// VARIABLES
real C(t) mM, // Sum of C1, C2, and C3 concentrations
C1(t) mM, // C1 concentration
C2(t) mM, // C2 concentration
C3(t) mM, // C3 concentration
Canalytic(t) mM; // C computed analytically
// INITIAL CONDITIONS
when(t=t.min) {C1=C10; C2 = C20; C3= C30; }
// ORDINARY DIFFERENTIAL EQUATIONS
C1:t = -k1*C1;
C2:t = -k2*C2;
C3:t = -k3*C3;
C = C1 + C2 + C3;
// ANALYTIC SOLUTION
Canalytic = C10*exp(-k1*t)+C20*exp(-k2*t)+C30*exp(-k3*t);
} // END
/*
DIAGRAM:
3
-----
\
C(t) = ) C (t)
/ i
-----
i = 1
d
where -- C (t) = -k *C (t) and C (0) and k are specified.
dt i i i i i
This is comparable to the washout of three substances from three
compartments where k = F /V , where F is the flow and V is the
i i i i i
volume of the compartments.
DETAILED DESCRIPTION:
Typical washout curves of solutes injected into an
organ or into the body can be decomposed into sums of
exponentially decaying functions, e.g.
C = C10*exp(-k1*t) + C20*exp(-k2*t) + C30*exp(-k3*t).
A set of noiseless data obtained from exact solutions is called EXPdecay.
It includes 5 representation of a solution, rounded off to 6, 4, 3, 2,
and 1 digit, to test the influence of accuracy on the parameter estimation.
See notes page of project file for more information.
Other such data sets, with and without added noise, and of different
duration can be tested.
The format of the .cdata file uses the first line for the data titles,
enclosed in double quotes separated by tabs or spaces. A second line
(optional) may contain the data units in double quotes separated by tabs
or spaces. Subsequent lines contain the data separated by tabs or spaces.
There must be at least two columns. The first column is the independent
variable and must be in ascending sequence.
EXPdecay.cdata file :
"time" "C_del_2_6" "C_del_2_4" "C_del_2_3" "C_del_2_2" "C_del_2 _1"
"sec" "mM" "mM" "mM" "mM" "mM"
0 1 1 1 1 1
2 .453386 .4534 .453 .45 .5
4 .314385 .3144 .314 .31 .3
6 .248923 .2489 .249 .25 .2
8 .205966 .206 .206 .21 .2
10 .174064 .1341 .134 .13 .1
SHORTCOMINGS/GENERAL COMMENTS:
- Specific inadequacies or next level steps
KEY WORDS: compartment, exponential decay, tutorial, washout curve,
reaction, data fit
REFERENCES:
REVISION HISTORY:
Original Author : JBB Date: 12/28/09
Revised by : BEJ Date: 12/31/09
Revision: 1) Update comment formats.
Revised by : GR Date: 01/06/10
JSim SOFTWARE COPYRIGHT AND REQUEST FOR ACKNOWLEDGMENT OF USE:
JSim software was developed with support from NIH grants HL088516,
and HL073598. Please cite these grants in any publication for which
this software is used and send one reprint of published abstracts or
articles to the address given below. Academic use is unrestricted.
Software may be copied so long as this copyright notice is included.
Copyright (C) 1999-2010 University of Washington.
Contact Information:
The National Simulation Resource,
Director J. B. Bassingthwaighte,
Department of Bioengineering,
University of Washington, Seattle, WA
98195-5061
*/
1.0 0.453386 0.314385 0.248923 0.205966 0.174064
1.0 0.4534 0.3144 0.2489 0.206 0.1341
1.0 0.453 0.314 0.249 0.206 0.134
1.0 0.45 0.31 0.25 0.21 0.13
1.0 0.5 0.3 0.2 0.2 0.1