Simulation and estimation - Biomedical Engineering Desk Reference - page 46

using the constants k 01 ¼ 0.0041, k 12 ¼ 0.0585, k 21 ¼
0.0498, and V 1 ¼ 7.3, c (0) ¼ p (0) ¼ 0.
This problem is designed to illustrate the use of global
variables in MATLAB. The mathematical problem is
formulated as a conservation of mass with the isotope
injection beginning at time zero; the length of the bolus
injection varies with trial.
The parameters of the model and the parameters of
the bolus injection are set as global variables, which can
then be used by the script that generates the graph and
the script that evaluates the system of differential
equations, called renal.m.
First, the system of equations (2.2.2) is coded as
a MATLAB function as per the MATLAB convention,
with arguments time t and the vector of function values,
y. All other model parameters are global variables.

Figure 2.2-16 A two-compartment model of renal clearance,
from Estelberger and Popper (2002) .

central volume and a less-perfused peripheral compart-
ment. The marker kinetics, measured indirectly from the
time course of radionuclide concentration in the two
compartments, is the result of the isotope injection, the
transport between the two compartments, and the renal
elimination process. Figure 2.2-16 is a diagram of the
Estelberger/Popper model.
The transport model can be formulated by a set of two
simultaneous differential equations describing the rates
of change of radionuclide concentration in the two
compartments:

function yprime ¼ renal(t,y)
% separate the components of the state:
% central compartment concentration and
% peripheral compartment concentration
c ¼ y(1);
p ¼ y(2);
%
global D tau;
global k01 k12 k21 V1;
% let time t be in minutes
% now compute f(t)
%
if t <¼ tau
f ¼ D/tau;
else f ¼ 0;
end
%
dcdt ¼ (f-(k01 þ k21)*c þ k12*p) /V1;
dpdt ¼ k21*c-k12*p;
yprime ¼ [dcdt, dpdt]';

dc
dt ¼ dðtÞðk 01 þ k 21 Þc þ k 12 p
dp
dt ¼ k 21 c k 12 p

(2.2.2)

where

dðtÞ¼ D

T ;

0 t < T (2.2.3)

There is a single script that solves the system of
equations three timesdonce for each of the different
injections of isotopes. A single figure with the three
concentrations in the central compartment is plotted
with different symbols. The succeeding plots are added
to the same figure by using the hold on command.

that is, the bolus isotope injection is delivered over T
seconds. For purposes of this exercise, assume that the
isotope is delivered uniformly during this period of time.
Assume also that the central volume V 1 is unknown.
The brief paper by Estelberger and Popper shows how
the group estimates the parameters of the model from
empirical data. The optimization procedure is beyond
the scope of this text and rather than detract from the
example, the reader is referred to Estelberger and Popper
(2002) for details.
In comparison to the simulation given in Section 2.2.3,
this set of differential equations will be solved using
MATLAB, rather than Simulink.

% Example 2.2.5 Renal Clearance
global D tau;
global k01 k12 k21 V1;
% parameters
k01 ¼ 0.0041;
k12 ¼ 0.0585;
k21 ¼ 0.0498;
V1 ¼ 7.3;
% initial value conditions
c0 ¼ 0;
p0 ¼ 0;

Example 2.2.5 Renal clearance.
Solve the system of two differential equations (2.2.2)
and the initial value condition (2.2.3) for renal clearance