- Email: [email protected]
Compartmental models with uncertain flow rates
GARY W. HARRISON* Department of Mathematics, Uniwrsi@ of Georgia, Athens, Georgia 30402 Received 1I May 1978; revised 14 Ju& 1978
ABSTRACT It is often extremely difficult to obtain precise values for the flow rates in compartmental models. Hence, in this paper it is assumed that the rate coefficients he within known intervals, though they may fluctuate in time. A method is given to compute upper and lower bounds at any time t for the concentration xi(t) in each compartment of any system whose coefficients remain in the given intervals. The bounds computed are the best possible, since there is some such system with xi(t) equal to the bound.
1.
INTRODUCTION
Linear compartmental models are often used to model biological phenomena such as pharmacokinetics or nutrient cycling in ecosystems. These models can be expressed as a vector differential equation dx dt
=Ax+r,
where xi, for i = 1,. . . , n, is the concentration of material in compartment i, ri is the flow rate into compartment i from outside the system, and, for i#j, uiixj is the flow rate from compartmentj into compartment i. The term uiixi is the total flow rate out of compartment i; hence ui,=-aOi-
2
a/c,,
(2)
k#i
where uoixi is the flow rate from compartment i to outside the system. In applications of these models it is often extremely difficult to determine the rate coefficients uti and ri. A frequent source of this uncertainty is the fact that each compartment may be an aggregation of subcompartments [l, 21. Other sources of uncertainty are inability to measure the flow rates precisely and incomplete knowledge of the process controlling the flow *This work was partially supported by NSF. grant No. DEB-02942. MATHEMATICAL
BIOSCIENCES
BElsevier North Holland, Inc., 1979
43: 13I- 139(1979)
131 00255564/79/010131+9$02.25
132
GARY
W. HARRISON
rates. (The linear model is usually an approximation for a nonlinear system.) Considerable work has been done on compartmental models with stochastic rate coefficients [e.g. 3-71. This paper takes a different approach to dealing with uncertain flow rates, by assuming that a priori bounds are known for A and r of the form
and computing intervals [x,(t),&(t)] which contain the components x,(t) of the solution. It is not assumed that the au and r, are constant; they may vary over time as long as they remain in the intervals given in (3). In practice, these intervals would be determined from physical and biological considerations or else would be statistically determined confidence intervals. Traditional methods for differential inequalities [8, 91 can be used to obtain upper and lower bounds, x(t) and x(t), for the solution of (1) under the restrictions (3) [l, lo], but these bounds may be very poor either when the entries of the matrix are not independent [as is the case for compartmental models because of Eq. (2)] or when the system is not quasimonotone as defined in [9] and [ 1 I] (which could happen in compartmental models if there were “recipient controlled flows” leading to negative off diagonal entries in A [12]). In the latter case some improvement is obtained by using Moore’s [13] correction for the wrapping effect [lo], but there is still no guarantee that the bounds are sharp. This paper gives a method to obtain the sharpest bounds possible, since there are solutions which reach the bounding functions x(t) and x(t) at each time t. An example of each type of difficulty and the improvement obtained is given in the last section. 2.
THEORY The equation dx x
(4)
=A(u(t))x+r(t),
where the matrix A is a function
of the parameter u
vector u satisfying (5)
and r(t) is a vector of inputs satisfying r
(6)
describes a family of differential equations, one for each possible pair u(t) and r(t) satisfying (5) and (6). A solution of (4) subject to the restrictions (5)
133
UNCERTAIN FLOW RATES
and (6) will mean a solution of any differential equation in this family. (In compartmental models the entries of the parameter u are the rate coefficients uij for i+j themselves.) THEOREM Let p(t) be the solution of the adjoint equations
-dpdt=(t) = -p=A(u*(t)),
(7)
where u*(t) satisfies p’(t)[A(u)-A(u*(t))]x for all
x >
GO
(8)
0 and all u satisfying (5) and let r*(t) satisfy, p=(t)r*(t)>p=(t)r
(9)
for alI r satis&ing (6). Zf x*(t) is a solution of (4) with u(t)=u*(t) and r(t)= r*(t), and x(t) is any nonnegative solution of (4) subject to the restrictions (5) and (6) with p =(0)x(O) G P =(0)x*(O),
(10)
pT(t)x(t)
(11)
then for all t > 0
Proof.
Equations
QJT(t)x*(t).
(4) and (7) imply that
4p ‘xl =pT[-A(u*)+A(u)]x+pTr. dt
(‘2)
Letting u = u* and r = r* in (12), it follows that
4p =x*1=p =r*. dt
Since it is assumed that x(t) > 0, applying that
which along with the initial inequality
the inequalities
(10) implies (11).
(‘3) (8) and (9) shows
GARY W. HARRISON
134
Remarks. The use of the adjoint equations (7) and the maximum condition (8) is in the spirit of optimal control theory. Indeed, the problem of finding bounds on a component x,(t) of the vector x(t) could be viewed as the problem of finding controls u*(t) and r*(t) that maximize or minimize x,(t), and Pontrjagin’s maximum principle could be applied to find necessary conditions for u* and r *. The simple proof which is given shows that the conditions (7)-(9) are both necessary and sufficient to maximize p Tx, and unlike traditional linear optimal control problems allow controls which act through the matrix A as well as controls on the input r(t). The assumption that x(t)> 0 is met by any realistic compartmental model. COROLLARY
Let dx z =A(t)x+r(t),
(‘5)
where the diagonal entries aii(t) of A are given by (2), and where r,(t) and a,(t) for i#j, i=O ,..., n, j= l,..., n, satisfv the bounds given in (3). Let p(t) be the solution of
dPT
dt
-p TA*(t),
=
(‘6)
where A *(t) = [a;(t)] with a&(t)=
a;(t) = for i # 0, i zj,
if pi(‘) > 0,
F i a,,
if pi(t)
aU
if p,(t) --p,(t) > 0,
a,
if
1
Pi(t)-Pi(t)
(‘7)
(‘8)
and r:(t)=
r’ i ‘i
if if
p;(t) > 0, pi(t)
(‘9)
and ai: is computed from a,$ and the a$‘s by Eq. (2). If x*(t) is a solution of dx* -=A*(t)x+r*(t), dt then any nonnegative solution of (15) subject to the restrictions (2) and (3) with p T(0)x(O)
0.
135
UNCERTAIN FLOW RATES
Proof: Apply the theorem with the parameters u(t) being the rate coefficients au for i#j, i=O ,..., n, j= l,..., n. Equation (19) clearly implies Eq. (9). Since Pr(A--*)x=
2 -pi(a,i-a$)xi+ I
2 i
2 (pi-pj)(ai,-a;)xj,
(21)
jti
Eq. (17) and (18) describe the proper choice for ai; to make Eq. (8) hold for all x > 0. 3.
COMPUTATION
OF BOUNDS
FOR SOLUTIONS
Equations (7)-(9) [or Eqs. (16)-(19)] and Eq. (13) can be integrated to calculate p(t) and p r( t)x*(t), which according to (11) gives an upper bound for p =(2)x(t). [Equation (4) with u = u* and r = r* or Eq. (20) is only needed if one wants to compute the full vector x*(r), not just p rx*.] But the stated objective was to compute upper and lower bounds for the components of x(t). Since there are no initial conditions specified for p, one possible approach would be to compute several different solutions p(‘), each with a different initial condition. The resulting set of inequalities p(“r( t)x( t) < ,~(~)~(t)x(‘)*(t) defines a hypeivolume in n-space which contains x(t), and finding the maximum or minimum of any component x,(t) of the vector x(t) subject to these inequalities is a linear programming problem. The difficulty with this approach is that the vectors p w all eventually approach the left eigenvector associated with the largest eigenvalue of the matrix A = [& + au)], making the region described by the inequalities long and needle-shaped and only giving good bounds on x(t) in one direction. Stewart [ 141 discusses the same difficulty in a slightly different setting. A better approach to find bounds for xi(T) at a specified time T is to use the boundary condition p(T)= e,, where e, is the ith column of an n dimensional identity matrix. Equations (7)-(9) [or (16)-(19)] can be integrated backwards in time to obtain p(O) and p ‘(0)x*(O), and then (13) [or (20)] can be integrated forward to obtain p ‘( T)x*( T). But since p(T) = e,, the inequality (11) reduces to xi(T) < xF( T). The same procedure with boundary condition p(T) = - e, yields the inequality - xi(T) < - xF( T), which gives a lower bound for x(T). Since x*(t) is itself a solution of (4) subject to the restrictions (5) and (6), these are the best possible bounds for x,(T). The conditionp(T)= e, or -e, is similar to the transversality conditions of optimal control theory. 4.
EXAMPLES
Example I is a four compartment system modeling the flows of nitrogen in an oak-hickory forest, which was obtained in [l] by lumping together
136
GARY
W. HARRISON
compartments of a seven compartment model constructed by Waide and Swank [ 151. The system has the form given by Eqs. (1) and (2), but, even assuming that the seven compartment model was exact, the process of combining compartments introduced uncertainty in the rate coefficients; hence the coefficients, to the extent known, are as follows: az, = 0.534,
0.158
0.3 10 < ad3 d 0.340,
ua, = 0.350,
0.017
az4 = 0.025,
ad, = 0.001,
0.054 < u4* < 0.088,
ao4= 0.005,
r, = 1.200,
(2’)
r2= 13.100.
All other entries of A and r are zero except the diagonal entries of A, which are given by Eq. (2). Figure 1 shows bounds for x4(f), soil nitrogen. The outer bounds were computed in [ 1] using Muller’s theorem [8, p. 961 and ignoring the depen-
X4( t)
3200
10
TIME
FIG. 1. explanation.
Bounds
20
(years)
for x4(l), concentration of soil nitrogen,
in Example
1. See text for
UNCERTAIN
137
FLOW RATES
X,(t) 25 20 15
10
5 0 -5
FIG. 2.
Bounds for x,(t) in Example 2. See text for explanation.
dence of the diagonal entries of A on the off diagonal entries. The inner bounds are the exact bounds computed from Eqs. (16)-( 19) and (13) with boundary condition p(T) = e, as described in the previous section, for a sequence of values of T. Example 2 is the equation
[
-0.5 0.95
I’d?]
dx/dt
= A(t)x + r(t), where
[ 14.925, - 9.045] < r’(t)
Q [ -;i5
-;T;],
< [ 15.075, - 8.9551.
(22)
Here the coefficients do not satisfy any relationship such as Eq. (2); they may vary independently within the prescribed intervals. The difficulty in obtaining bounds for the solution is that the system is not quasimonotone, due to the off diagonal entries of A having opposite sign. Figure 2 shows bounds for x,(t). The outer bounds were computed using Muller’s theorem [8, p. 961; the middle bounds were computed in [lo] using a transformation suggested by Moore to reduce the “wrapping effect” [13],
138
GARY W. HARRISON
and the inner bounds (7)-( 11) with boundary 5.
are the exact bounds computed by applying Eqs. condition&T) = e, for a sequence of values of T.
CONCLUSION
It is seen that the exact bounds for solutions of Eq. (4) subject to restrictions of the form (5) and (6) can be computed by the application of the theorem of this paper, and are much sharper than bounds computed by previous methods. Although, for ease of notation, it was not indicated in Eq. (5) that the intervals [c, il] and [ci, FJ vary over time, the results still apply as long as the functions s(t), r(t), etc. are smooth enough for the appropriate differential equations to have unique solutions. It is hoped that these results will provide a useful tool, in addition to the traditional stochastic process methods, for dealing with uncertain flow rates in compartmental nodels. The most obvious connection between the two is through the use of confidence intervals. Suppose that the coefficients u,(r) and r,(t) are random functions with a joint probability distribution, and for each choice of a set of coefficient functions {u,(t),r,(t); i= 1,. . . ,n, j= 1, . . . ,n} there is a corresponding “sample trajectory” x(t) which is the solution of Eq. (1). If the intervals in (3) are, say, 95% confidence intervals in the sense that, with probability 95%, any set of randomly chosen coefficient functions satisfy the bounds in (3) for all O
lumping in nutrient cycling processes, in Enoironmentat Chemistry and Cycling Processes (D. C. Adrians and I. L. Brisbin, Eds.), DOE
Symposium Series CONF760429, 1978, pp. 121-137. 2 R. V. O’Neill and B. Rust, Aggregation error in ecological models, in press. 3 J. A. Jacquez, Compartmental Ana&sis in Biology and Medicine, Elsevier, New York, 1972. 4 J. P. Lehocsky and D. P. Gaver, A diffusion-approximation analysis of a general n-compartment system, Math. Biosci. 36: 127-148 (1977). 5 T. T. Soong, Pharmacokinetics with uncertainties in rate constants, Math. Biosci. 12: 235-243 6
(1971).
T. T. Soong, Random Differential Equations in Science and Engineering, Academic, New York, 1973.
UNCERTAIN
FLOW RATES
139
7 J. 0. Tsokos and C. P. Tsokos, Statistical modeling of pharmacokinetic systems, J. Dyn. Syst. Meas. Control 98: 37-43 (1976). 8 E. F. Beckenbach and R. Bellman, Inequalities, Springer, New York, 1965. 9 W. Walter, Differential and Integral Inequalities, Springer, New York, 1970. 10 G. W. Harrison, Dynamic models with uncertain parameters, in Proceedings of the First International Conference on Mathematical Modeling (X. J. R. Avula, Ed.), University of Missouri, Rolla, 1977, Vol. 1, pp. 295-304. 11 V. Lakshmikantham and S. Leela, Differential and Integral Inequalities, Vol. 1, Academic, New York, 1969. 12 J. B. Waide and J. R. Webster, Engineering systems analysis, applicability to ecosystems, in Systems Analysis and Simulation in Ecology, Vol. 4 (B. C. Patten, Ed.), Academic, New York, 1975, pp. 329-371. 13 R. E. Moore, Inter& Analysis, Prentice-Hall, Englewood Cliffs, N. J., 1966. 14 N. E. Stewart, A heuristic to reduce the wrapping effect, BIT 11: 328-337 (1971). 15 J. B. Waide and W. T. Swank, Nutrient recycling and the stability of ecosystems: implications for forest management in the southeastern U. S., Sot. Amer. For. Proc., 404-424 (1975).