Equilibrium points for nonlinear compartmental models.

Equilibrium Points for Nonlinear Compartmental Models DAVID H . ANDERSON Department of Mathematics, Southern Methodist University, Dallas, Texas 75275 AND

TOWANNA ROLLER Department of Mathematics, Dallas Baptist University, Dallas, Texas 75211 Received 6 December 1989; revised 30 July 1990

ABSTRACT Equilibrium points for nonlinear autonomous compartmental models with constant input are discussed . Upper and lower bounds for the steady states are derived . Theorems guaranteeing existence and uniqueness of equilibrium points for a large collection of system are proved . New information relating to mean residence times is developedAsymptotic results and a section on stability are included . A recursive process is discussed that generates iterates that converge to steady states for certain types of models . An interesting range of models are included as examples . An attempt is made to provide general qualitative theory for such nonlinear compartmental systems.

1 . INTRODUCTION Compartmental analysis is a useful tool in the study of epidemiology, drug kinetics, lipoprotein kinetics, chemical reaction kinetics, metabolic systems, ecosystems, and other systems . In compartmental analysis the dynamic system is separated into a finite number of component parts, called compartments, and flow between the compartments is described by ordinary differential equations usually arising from mass-balance considerations [3, 9, 11, 18, 21] . In this paper we study equilibrium points and qualitative analysis for continuous dynamic nonlinear autonomous models x(t) = f(x) which can be factored into the form x=A(x)x+b,

(1 .1)

where x= x(t) > 0 and b ~> 0 are n X I nonnegative vectors and A(x)=[a,t(x)] is an n x n "compartmental matrix." (In this article, all inequalities, MATHEMATICAL BIOSCIENCES 103 :159-201 (1991)

©Elsevier Science Publishing Co ., Inc ., 1991 655 Avenue of the Americas, New York, NY 10010

159 0025-5564/91/$03 .50

160

DAVID H . ANDERSON AND TOWANNA ROLLER

including those between matrices or vectors, are taken componentwise .) An individual element a jx), i ~ j, is the fractional transfer coefficient for material from the jth compartment to the ith compartment, while a o;(x) 3 0 is the fractional excretion coefficient to the environment from compartment i . The diagonal element a„(x) is by definition the negative of the sum of all coefficients for material leaving compartment l, n a„(x)=-

a,,(x))

i=o i*i

Physical considerations involving mass balance require that for appropriate x =[x,] 3 0, A(x) has

(1) Nonnegative off-diagonal elements :

a, (x) 3

0;

i * j

(2) Nonpositive diagonal elements : a ;;(x)c0 ;

i=1,2, . . .,n

(1 .3)

(3) Nonpositive column sums : n

~aii(x)c0 ;

1=1,2, . ., n .

i-1

From (1 .2) and (1 .3) observe that the jth column sum of A(x) is - a o (x) and A(x) is a column diagonally dominant matrix . Any square matrix A(x) that has the properties of (1 .3) when evaluated at a given vector x 3 0 will be called a compartmental matrix . The components of A are assumed to be continuous functions of the state variable x and to satisfy a Lipschitz condition also . This is no real restriction, because usually the functions a,,(x) are either polynomials or rational functions such that the denominator never vanishes . The purpose, of course, is to establish existence and uniqueness of a solution to system (1 .1) for time t 3 0 [4, 151. Here the matrix A(x) is not dependent on t except by way of the state vector, and the input vector b---[b,] >- 0 is constant over time . As is always the case in compartmental analysis, x remains a nonnegative vector over time due to condition (1 .3) and x(0) 3 0 . For, choose i to be an arbitrary index in {1, 2, . . ., n] . In vector x, let x ; = 0 and all x i 3 0 for j * i . Then the

161

NONLINEAR COMPARTMENTAL MODELS

ith component equation of (1 .1) yields n

n

,i i = F, a il (x)x t +b,= L a i,(x)x l +b1 i=l i=] I I

0.

Hence x i = 0 is in general a repulsive hyperplane for any trajectory in the nonnegative orthant of R". Thus if x(t) is on the boundary of the nonnegative orthant, its direction of movement cannot take it outside that region . An equilibrium point or steady state of 0 .1) is a nonnegative vector x that is found by solving the simultaneous nonlinear algebraic equations 0= x= A(x)x+h .

(1 .4)

Some questions of interest include existence and uniqueness of x, bounds on x, any information on computing z, stability of x, convergence of trajectories to z, and how x depends on x(0), b. Let R', =(x E R" : x 30) . For any given model (1 .1), our work will be restricted to the subset of R'defined by Sl=-{xE R'+ : A(x) is a compartmental matrix} .

(1 .5)

We will often be working with the function g(x) - - A(x) -16

(1 .6)

and thus will need A(x) to be invertible . For the given model (1 .1), define fl-,

(x a fl : A(x) is an invertible matrix),

(1J)

which is the set of all x E fl such that the compartmental model contains no traps [16] . Since A(x) is compartmental, then - A(x) -1 contains nonnegative entries, that is, - A(x) - '3 0 for each x E fl_ 1 [2, 3, 20] . With b 0 we have g : f2_ 1 -j R" . Also any equilibrium point of compartmental system (1 .1) in ,fl_, is a fixed point of g : x=g(z) . Example 1.1 . The domain fl on which system (1 .1) is compartmental often depends directly on the parameters of the model. A good illustration of this dependence is the famous three-dimensional Lorenz system from meteorology [5, p . 217; 17 ; 19, p . 92 ; 29] . This model factors into form (1 .1)

162


with b = 0, -a A(x)= P - x3 0

a -1 x,

0 0 - v

and parameters p, a, v > 0 . This matrix is compartmental if 0 < x, < I - a (which requires that parameter a < 1), if x 2 3 0, if 0 5x 3 < p, and, to ensure a 01(x) = a - p + x 3 is nonnegative, if a 3 p . Note that det A(x) =(p-X3-1)av, which is zero if and only if x 3 =p-1 . We must have p < 1, because if p >_ 1, then o > p 31, a contradiction to a < 1 . Domain SI is defined by these conditions . The point z = 0 is in db_, =f1 and is the only equilibrium point in 12 . The system appears to be well behaved on f . However, for a much different selection of parameters (e .g ., Lorenz [29] fixed a=10, v=8/3, and p=28), the model is not compartmental and illustrates strange attractors and chaotic behavior 119, p . 95] . Models that factor into the form of Equation (1 .1), where A(x) is compartmental (and where b is constant or nonconstant) have proliferated in the literature in biology, ecology, physiology, pharmacology, chemistry, and engineering. Recent developments in compartmental analysis are found in [31, 191, 111], [181, and [211 . Articles devoted to nonlinear compartmental systems are [4], [15], [22], [25], [26], [311, [321, [36], [37], and [40] . Particular examples that have the nonlinear compartmental structure of (1 .1) are the SIRS epidemiology model [10, 14]; AIDS models [12, 23] ; chemical reaction models [3, 14, 21, 24, 30, 42] ; chemostat models 114, 38] ; air-flow and weather models [17, 291 ; guerilla/conventional forces combat models [5, 8] ; predator-prey models [8] ; helminth infections [33] ; enzymatic systems [11, 18, 28, 381 ; metabolic control systems [9, p . 75); and blood sugar/insulin models [30] . For a variety of examples, see [34] . This paper is organized as follows . Section 2 presents bounds and other basic information on any equilibrium point x and input vector b of (1 .1) . Also given are bounds on x and the function g(x) of (1 .6). Sections 3-5 are about existence and uniqueness of equilibrium points . Certain models of type (1 .1) have easily described hyperplanes in R'in which all equilibrium points lie . Section 6 discusses these hyperplanes and the convergence of the trajectories to equilibrium points on such hyperplanes . Model stability via Lyapunov functions and linearization is the subject of Section 7 . Section 8 gives new information on the relation between fractional transfer coeffi-

NONLINEAR

COMPARTMENTAL

163

MODELS

cients in (1 .1) and associated mean residence times . An iterative procedure on the function g(x) of (1 .6) that converges to an equilibrium point for certain compartmental models is discussed in Section 9 . 2.

BASIC INFORMATION ON x, g(x), x, AND b

Recall that if real vector y-[ y ,], then the definitions of various norms are Ilyll,=-Ely ;l, IIYII2--[Ey?]" 2 , and IIyL=-maxly,l . Bounds on g(x) may be obtained by considering either A(x)g(x) _ - b or g(x) _ - A(x) - 'b. These two equations yield the bounds

IIA(xl)II `Ilg(x)II 0, and the matrix and vector norms are compatible [27, 39] . Now, if IIA(x) - 'II can be bounded uniformly above by the constant U, and if PA(x)II c L for all appropriate x, we have hill y Ilg(x)II c

Ullbll,

giving a range for g . Note that by (1 .2) and (1 .3), the column norm IIA(x)Ili is IIA(x)ll t =max{2la ii (x)I-a„ i (x)}

for each x

E=

D.

i

Bounds on IIA(x) -I II appear in Section 4 . More explicit information may be obtained by treating the equations componentwise . Let i be any index in (1, 2, . . ., n} such that a,;(x) $ 0 for all x D and input component b ; > 0. Then for any equilibrium point i E Sl, Equation (1 .4) yields E

n

0= =a i,( i)x,+

E a,j(x)xi+b, . i=1 jti

Thus - a„(i)xi a b„ and hence x, a - b, /a,,(x) > 0 . Now let i e 11 be any equilibrium point such that component x, = 0. Then input component b ; must be zero. Let Q,= {j :j#ianda, j (x)>0forallxE11} .

1 64


Thus by (1 .4), n

S

0

aij(x)xj i

j-1

5 a,,(x)x j . jeQ I

Since a ij > 0 and x j >~ 0 in the last summation, then in addition to x ; = 0 in z, we must also have components xj = 0 for all j e Q ; . Let g(x)=-[gk(x)]k 1 . Then, the ith component equation of A(x)g(x)= -b is n

b;,

a,k(x)gk(x)

i=1,2, . . .,n .

(2 .1)

k=1

Let j3 be the sum of components b„ that is, 13 = ~ I bI113 0 . Summation over i in (2 .1), reversal of the order of summation, and use of (1 .2) gives n

F,

(2 .2)

gk(x)aOk(x) - 6 , J

k-1

where gk(x) and aok(x) are nonnegative. Some easily obtainable results from (2.2) are

$ F"-,a nk x)

(2 .3)

E ~Ig(x)lIm

and (2 .4)

R % gk(x)aok(x)

for any particular for all x e fl_ 1 ,

k.

From (2.4) we have, for any

P

a0k(X)

-i

k

such that a 0k (x) # 0 and

gk(X) i 0 .

(2 .5)

If A(x) is a single-exit compartmental matrix [1, 3] via compartment k, then (2 .2) becomes [n'`

= L gj(x)aoj(x) = gk(x)a0k(x) j=1

and the inequality (2 .5) is sharp . Define a0(x)--[a01(x), . . .,ao„(x)]T for any

165


x c- fl . When x E dl_ 1 , interpret (2 .2) as the inner product on R", (2 .6)

(g(x),ao(x)) =,B . By the Cauchy-Schwarz inequality, 13 < 1100121lao(012 .

(2 .7)

The inner product (g(x),ao(x))=/3 is especially useful at steady state

s E fl_ ,, where we have $ - (z,ao(X)) = 0 .

(2 .8)

For any x E dl, the left-hand side of Equation (2 .8) has a physical interpretation as the derivative of the total mass in the system . Let s(x) =E0. 1 x,= lxli1 be the total mass in the system, and let a(x) ds/dt . Then, premultiplying (1 .1) by [1,1, . . .,1], and using (1 .2), a(x) = R - (x, a o (x) ) .

(2 .9)

The following result is consequence of the above discussion . THEOREM 2.1

Any equilibrium point for compartmental system (1 .1) with constant input vector satisfies the equation 6 -(x,a o (x)) = 0 . If an equilibrium solution z >, 0 of model (1 .1) exists, then /3 must he in the range of the function (x,a o (x)) mapping Sl into the interval [0,o) . Example 2.1 . The condition 6 is in the range of (x,a o (x)) places bounds on the size /3=IIbII1 of the input vector b. In [11, p . 7] Cherruault presents a three-dimensional enzymatic system that can be put into compartmental form (1 .1) where 0

-k21

A(x)=

k21

0

k32

0 V k1+x2

k 32

k23 -k23

and all parameters are positive . A quick analysis of the digraph in Figure 1

166

DAVID

H.

ANDERSON

AND

TOWANNA

ROLLER

3b 2 k2,

1

10

2

3 Fie . 1 . Digraph of an enzymatic system .

associated with A(x) shows that all compartments have flow to the environment for each x 3 0 . There are no traps [16], A(x) is invertible for all x, and (1 = R ; = fl-, . The only nonzero excretion from the system is a02(X) = V/(k t +x 2 ) . The range condition of Theorem 2 .1 may be stated as

min

Vx 2

0 0 for all x e 0_ t and some d ; . Then by (2 .5), (2 .12)

dR; i, a ( j x) %gX)=x >- 0 .

Example 2.2 . Consider the chemostat model for the harvesting of bacteria 114, 381 . This model can be put into the form (1 .1) with A(x) and b given by

-D

A(x) _ 0

VX 1

k+x 1

vxt

-D- S~ k+ x2 )

b- L O

In this model, x, represents the concentration of a certain species of microorganism in the chemostat and x 2 represents the concentration of a nutrient that is essential to the growth of the microorganism but is in limited supply. The parameters v, D, S, and k are all positive . Also, if S 51, then the system is compartmental for all x e S2= R? = fl, Let 5=1, D=1, k = 1, v=2, and b 2 =10 . The digraph for the model is shown in Figure 2 . Upper bounds for z, if it exists, can be obtained through Equation (2 .12), which yields the componentwise inequalities 0 -< g(x)= x < [10,101 T• By (2 .3), IlxllL= 11g(x)IL . >- 5, but a better bound can be obtained by using (2.10); that is, 101V-2 6 Ilg(x)112 or 11x112 = Ilg(x)112 > 7 .07. At steady state, Equation (2 .8) is the most useful of all and in this case encompasses the above results, yielding Xt + x 2 = 10 . The acutal equilibrium points for this model are [0,101T and [9,11T, which are computed directly from (1 .4) . Note that since x, rep :esents the concentration of bacteria, we are not interested in a steady state for which x t is zero even though the steady

168

DAVID H. ANDERSON AND TOWANNA ROLLER

2x , a,2(x) -1

2

Ftc . 2 . Chemostat model digraph .

Fic . 3 . Chemostat model vector field .

169


state is biologically feasible . As might be expected, the steady state [0,10]7 is not stable . In contrast, [9,1]T is a stable equilibrium . See Figure 3. 1X„x 2 ]T for the chemostat is such that Suppose an equilibrium point z = x„ x 2 > 0 . The nonzero factor x t can be divided out of the first equation to get D1 v = x2 /(k + x 2) < 1 . Hence if x i * 0, then v > D must be true in the parameter list . We take this inequality as a condition for the model . This steady state has the positive components

2 Dk ) 5D v-D ' (b

Dk v-D'

The Jacobian of the model evaluated at this equilibrium point has eigenvalues - D and -Dx,0-D/v)/Sx 2 , both negative . Hence this point is locally asymptotically stable . if the first component x i in an equilibrium point is zero, then x2= b 2 /D > 0 follows immediately from the second equation of the chemostat model . The Jacobian at the point x= [0, b 2 /D] T has eigenvalues - D and - D + vb 2 /(kD + b 2 ) . The second eigenvalue is positive unless b 2 /D Dk/(v-D) . Compare this to the form of x, above . Under the condition b 2 /D > Dk /(v - D) . z = [0, b 2 /D] T is a saddle point . See Figure 3 . THEOREM 2.2

Suppose there exists an equilibrium point x for model (1 .1) in which input b t 0 . Then x and b are related by the orthogonality condition

( y ,b)=0

if A(i) T y=O .

(2 .13)

Proof. The equilibrium point x satisfies A(x)x- = - b . Thus - b is in the range space of the matrix A(x) . By linear algebra, Rge A(x) = [Null A(x)T] 1. Hence for any vector y such that A(x)Ty = 0, we must have the inner product (y,h)=0 . This is (2 .13) . Now suppose A(x) is a closed compartmental matrix so that a oi(iii)=0 for all i=1,2, . . .,n . Then we can take y --- [1,1, . . . ,1] T in A(x) Ty = 0 . The inner product (y, b> = 0 reduces to $ =b, + h, + . . . + b,, = 0, which is already known from Equation (2 .8). Now suppose A(x) represents an open compartmental system with no traps, so that A(x) is invertible . Premultiply A(x)Ty=0 by [A(i) - i]T. The result is y = 0, which is the trivial case of (2 .13). ∎

Example 2.3. This system is a two-compartment open model with the second compartment Langmuir-saturating and a trap [18, p . 186] . The

170


equation in (1 .1) has

A(x)=

-

kol - k210 - Q2x2) k,,(1 02 x 2 )

0 0

b =[b„0]T with b 1 > 0, and x in the compartmental domain dZ = (x : 0 -< x„ . The set Q- 1 is empty. Condition A(x) t y = 0 in (2 .13) 01=Y1b1=0 is satisfied for arbitrary real y 2 . However, Equation (2 .8) immediately yields the other equilibrium component : 0 = -(z,a0(x)>=b1-x,ko, . Now suppose there exists a positive constant m such that a01(x) 3

in,

for all i = 1,2, . . ., n ;

x e 11 .

(2 .14)

For instance, this would occur if domain 4 is compact, because each continuous a o ;(x) > 0 would achieve its minimum d ; > 0, and m =min d ; . From (2 .9), the total mass s(x)=llxlli of the system (1 .1) satisfies sC (3-Emx ;=/3-ms . Multiply both sides of this equation by exp(mt) to get the differential inequality 711\

s(x)-m exp(mt) c0 . j

(2 .15)

THEOREM 2 .3

If (2 .14) holds for model (1 .1), then plx(t)111-13/m]exp(mt) is a nonincreasing function for all t > 0 . Evaluating this function at t and 0 yields the total mass bound

Ilx(t)Iltcllx(0)Ilie -m.

+

m[1-e - ° ] .

(2 .16)

The right-hand side of (2 .16) is an upper bound on x ;(t) for all i . If initially 11x(0)111 c p/m, then Mx(t)Ilt c 3/m for all t > 0 .

A s& .nilar argument yields the next statement .

171

NONLINEAR COMPARTMENTAL MODELS THEOREM 2.4

Suppose there is a positive constant M such that a 0i(x)_ Ilx(0)Ilie-M, If Ilx(0)Ili

> 13/M,

then IIx(t)IIi

a [3/M

+ k[1_

e -mr]

for all t > 0 .

If each aot is constant and the same value for all i (e .g ., Example 2 .2 when S=1), then the inequalities in Theorems 2 .3 and 2 .4 become equalities. All equilibria that exist lie on the hyperplane IIxPIi = P/aur . Likewise, a lower bound on the the ith component of x(t) can be developed provided there exists a positive constant u ; such that Ia,,(x)I C u ;,

for all x e 11 .

(2 .17)

For, the ith component equation of (1 .1) yields n

b i 6h;+ Ea ii (x)xi i-t j*i

(2 .18)

= X i - a„(x) x, 5 xt + u,x, . Multiply (2.18) by exp(u,t) to get

TI

4X' ( t)-u, Jexp(ui t)S,0

(2 .19)

and the next result. THEOREM 2.5

Suppose condition (2 .17) holds for compartment i . Then [x,(t)-b,/u,] exp(u,t) is a nondecreasing function oft . Moreover, for all t > 0, x ;(t)>x,(0)e ",t+u`[l-e-"=1,

(2 .20)

and ifx ;(0)>b;/u 1 ,then x,(t)>bi/u, for alit >0 . 3.

AN EXISTENCE THEOREM BASED ON THE EXCRETION COEFFICIENTS

Equation (2 .9), o(x)=/3-, set equal to 0 is an attractive condition, since in many cases x may simply be increased or decreased until

172


equality is obtained . However, equality is a necessary rather than sufficient condition for an equilibrium point . These ideas are used by Sandberg 136, p . 274] to state a sufficient condition : If there exists a positive constant k such that a(x) < 0 for 11x111 a k, then there exists at least one equilibrium solution to the system . An important special case of Sandberg's result is now stated . THEOREM 3,1

If, for each i=1,2, . . .,n, the excretion rate coefficient a 01(x) is bounded away from zero as x, -> +o, then there exists an equilibrium point x > 0 of (1 .1). Proof. Let a0,(x) 3 d; > 0 for all x such that x ; > y; > 0 for some d ; and y, . Define d =- min d„ y =- max y„ and k =- max((3n /d, ny) . Then Ilxll, 3 k guarantees x ; > k/n for some i . For this particular i value, we have x . > (3 /d and x, > y 3 y, . Thus, a(x)=(3 - (x,a0(x)) S(3 - x,a01(x)

cR-

dd,c(3- dd=0,

and the theorem is proved .

∎

Example 3.1 . Consider the three-dimensional predator-prey model given in [8, p . 227] . In this model the predators (x 3 ) attack the adult prey (x 1 ), while the young prey (x 2 ) are fairly well protected . The model may be written in form (1 .1) where

A(x)=

-a1-hX3 e dx 3

a2 -(a,+a 2) 0

0 _C _

All parameters of the model are positive, and fi = R ; if we assume a 1 > e and h > d . We have a0i(x)=(a1-e)+(h-d)x3, a 12(x)=a 1, and a 03 (x)=c . As long as a 1 > e in the parameter list, we have that all excretion coefficients are bounded away from 0 and fl-, = Q . Thus, with any 3 x I nonnegative input vector b, Theorem 3 .1 applies to this model to guarantee the existence of at least one equilibrium point x . Theorem 2 .3 provides an upper bound for x. The real functional a(x)= (i - (x, a 0(x)) = s(x) of (2.9) divides state space into three regions : a(x) < 0, a(x) = 0, and a(x) > 0 . Level curves can be defined in state space by s(x) =-11x111=q for positive constant q . Thus, the total amount of material in the system is constant along a level curve . At points where a(x) < 0, x(t), the solution to (1 .1), crosses the level curves in a direction toward the origin as t increases . At points where a(x) > 0,


173

6 5 4 3 2

unique steady-state solution Fic . 4. The o(x) analysis for a 2-compartment model . the solution crosses the level curves going away from the origin . If u(x) = 0, the state vector is at equilibrium or else is traveling parallel to a level curve . The set {x : a(x) 6 0} is a closed set by continuity of the elements of the A(x) matrix. Much information can he gained about a particular problem if the graph of o(x) is known relative to s(x). Example 3.2.

Consider the compartmental model (1 .1) determined by

A(x)-

-(xi-9 .2x2+21 .8x,+2) 2

3x, -(3x,+4)

and b=[7,101'. Here 12 =R ;. =l-„ and this model has the unique steady state x = [0 .805,1 .809]T. For this example, s(x) = x, + x 2 and o(x) = 17x 4 +9.2x?-21 .8x? -4x 2 . Note that in Figure 4 for any point where s(x) > 6 (i .e ., any point to the right of the diagonal line x, + x 2 = 6) we have a(x) 5 0 . Thus with c = 6, we are guaranteed at least one equilibrium point to the system . 4.

EXISTENCE THEORY BASED ON THE INVERSE OF A(x)

THEOREM 4.1

Suppose Q_ 1 =S1=R ; for the given matrix A(x) of model (1 .1), and suppose there is a constant U such that IIA(x) - 'II 6 U for some matrix norm

174


and for all x in S---{xe R, : Ilxll c Ullbll} . (The vector and matrix norms are assumed to be compatible .) Then there exists a solution z e S of A(x)x+b=0. Proof. With x arbitrary but fixed in S, there is a unique solution g a. for g of the equation A(x)g = -b, because A(x) is invertible . Thus, the relation x -. g, is a function on S . Moreover, the range is in R„ because - A(x) - t _> 0 and b > 0 . Since Ig,ll=11A(x) -1bII_ d i > 0

(4 .1)

for some set (d ;}," , and for all x . Then, by (2.5), for all x c fl,

. 0 < g,(x) c ( x) _< a of d,

For any arbitrary but fixed x E 12_ 1 , bounds on IIA(x)-1II can also be constructed by the Neumann series [2; 3, p . 75] A - ' = D+DE+DE 2 +DEI+- •- C0,

(4 .2)

where D(x) -_ and E(x)=I-A(x)D(x) . Since E(x)>0 and its spectral radius p(E(x)) 0 for all x e Q_1 . Then for all x e S2_ 1,

rii(x) 6

( 4 .4)

. ao,1x)

Since T(x) =- A(x) - 1 , we obtain

Proof.

n

Tik(x) aki (x) k=t

Therefore, n

F- Tik(x)aki(x)+ 1

rii(x)aii(x)

k=l1 k*i

From (1 .2), n

rz

Tik(x)aki(x)+I =re,(x)a0t(x)+ k=1 kxi

k-1

k*i

Tli(x)aki(x) ,

with all terms nonnegative . We know that Ti,(x) 3 T,k (x) for all k =1, 2, . . ., n [3, p . 86], and thus 13 ra(x)a o ,(x) or T„ (x) E 1/a o,(x) . ∎ Anderson [3, p . 86] has shown that for each x E fl-, -ii(x) 4TU(x)

(4 .5)

for all 1 4 i, j 6 n . Note that equality can hold in Theorem 4 .2 and in (4 .5) when A(x) is single exit from compartment i [l ; 3, p . 86] . Example 4.1 .

Consider the elementary epidemic model [8, 10, 15] .C 1 =-Sx 1 x 2

-Ax,+b l

X2 = Sx1x2 -

yx2 + b2>

176


bl

S X2

1

0

2

Fu; . 5 . Digraph of the elementary epidemic model,

where x 1 and x 2 represent the susceptible and infected members of the population, respectively . The positive parameters S, A, and y are constants relating to how easily the disease is transmitted, what percentage of the population is becoming vaccinated, and what percentage of sick individuals are recovering (or dying) from the disease . This system may be written in form (1 .1), where (5X2

A(x)_

-A

0 -y

Sx,

The digraph for the model is given in Figure 5 . Since derA(x)= y(Sx 2 + A) z yA > 0 for all x 0, we have U- 1 = S1= R ; . A "particle" that begins in compartment 2 can never enter compartment 1, so T 12(x)=0 for all x 3 0 from the mean residence time interpretation . With a ,1(x)= A and a o2(x)=y, inequality (4 .5), Theorem 4 .2, and Theorem 4 .1 now guarantee the existence of at least one equilibrium point

z=-A(z) lb~

1 A

bl

0

A 62

Y

y

b2

b y

For each arbitrary but fixed x e f2 1 , let A 0(x) be the real negative eigenvalue of the n x n matrix A(x) that has maximum real part (and smallest modulus) . Then bounds on all T,j(x) can be obtained as follows.

177


For any fixed x E (1_ 1 we have r . (x)'C -1/A„(x), where T(X) is any one of the diagonal elements of T(x)=-A(x)-' .

Proof. For an arbitrary but fixed x e S2_ i , we have that T(x) 3 0. Since A„(x) denotes the eigenvalue of A(x) with smallest modulus, -1/(A„(x))p(- A(x) - ') via the Perron-Frobenius theory of nonnegative matrices . By [6, p . 271, the spectral radius

p(-A(x) - ') %p`diag(r ;;(x))~=maxTkk(x)>r, ;(x) .

∎

Theorem 4 .2 or 4 .3 can be used to get upper bounds on -A(x) - ' and IIA(x) 'II as follows . Suppose all excretion coefficients are bounded away from zero as in (4 .1). For each x E U_„ 0 < min d k < min a ak(x) < IA„(x)I,

where the last inequality uses a result in [3, p . 65] . Since A„(x) < 0 is real, then for each x E U_ t , 1 1 -1 min dk -> min aok(x) ' A„ (x) '' r,,(x) Ti(x) by Theorems 4.2 and 4.3 and Equation (4 .5) . Because [T;;(x)]=-A(x)-', the inequality in (4 .6) will provide a constant upper bound for IIA(x) - 'Jl arid all elements of -A(x) - '. THEOREM 4.4 Suppose the compartmental matrix A(x), for fixed x e d2_„ is such that a ei (x)>0 for all j = 1, 2, . . ., n . Then

1 EIIA(x) in ax aoj(x) 1

'III, 1 in ina o1 (x)

(4 .7)

I

Proof. For each such x, define the n X n matrix H(x) _ [ h 11 (x)] = A(x)T. Let z e R° . Then H(x) is a strictly (row) diagonally dominant matrix . Using

178


Blum's work [7, p . 159],

IH(x)zL3miniIh ii (x)I-

Ih i,( x)I j-t I*i

= min { -a ii (x)-

L

a ;i (x) }

IIzIIm

j=1 i' i

_ [min a oi (x)]

(4 .8)

IIzIIa .

Now let z = H(x) - 'y for any y e R" . Then (4 .8) implies

IIH(x)-'yllms

hall . (4 .9) min i a,, (x)

Hence IjA(x) '

_II[A(x) 11Tp_I~A(x)n~~

In

=IIH(xY'IImc min i aI 0i (x) by (4 .9). The other half of the desired inequality follows from the proof of Theorem 4 .3 . For, from [3, p . 65], 0 < IA"(x) < max a oi (x) . I

Then

max,aoi(x)`A"(x)-p(-A(x) ')EIIA(x)

'IIi

because A"(x) < 0 is the eigenvalue of A(x) of smallest magnitude . Example

4.2 .

∎

If Theorem 4 .4 is applied to the matrix A(x) in Example

4 .1, then 1 1 max(A,y) CIIA(x) I6E min{)L,y)

179


for any x e R 2. . Likewise, for each x 3 1> 0, the matrix in Example 3 .1 yields

I

A(x 3 )

~ h

M(x3) 0, '3 > 0, and, by Equation (2.8), (i -µ!1T'1!1 = a 3 > 0 . Hence 11x111 < b, /µ, which agrees with (4.7) and (2 .16) . Notice that by Theorem 2.4, Ilzll, >P /M, where M = µ + d and / = b, . Define the parameters

a=-

rc-(v+µ) 1 . Quantity a is positive by the epidemic threshold condition in (4 .10) . From the model equations and (2 .8), we derive the steady states (v+ 011 z 111

T'3=

re _ (d+µ)' 3 T'2= pv

b,A11301 1 d

(1-p)vi 2 T'4= µ

(4 .12)

These can be computed sequentially if Ilzll, is known in advance . This can be done as follows . On the right-hand side of the last two equations of (4 .12), put '3 and '2 in terms of Ilzllt from the second and third equations . Then add the equations of (4 .12) and solve for the only unknown, IIa!I1 (recall that fraction a/y < 1): IIXII,=

µ+ad /y -

( 4 .13)

The above AIDS model is an example (see Figure 6) of the next type of system . A cascading system is a system where the compartments are numbered from the top down and the compartmental matrix A(x) in (1 .1) is lower triangular (including the diagonal) . The following result is useful when the diagonal elements of A(x) are bounded away from zero . THEOREM 4.5

If model (1 .1) is for a cascading system in which every a ;,(x) < 0 for all x E S2, then any equilibrium point x is such that component b,+b2+ .- •+ bt I

lan (x)I

for j=1,2, . . .,n .

182


This inequality is strict if there is some index i, i < j, such that column i in A(x) is strictly diagonally dominant for all x e 0 .

Proof.

From Equation (1 .4), we find that for any k =2,3, . . .,n, k~-'` I

akklXk - bk + L a ki x i , t=,

(4 .14)

because A is lower triangular . The proof builds inductively on (4 .14), and it is probably most easily seen by looking at the following demonstration of the systematic construction process . Suppose k=5 in (4 .14) . This equation can be written as a551x5 = b5

+

1a 54 u1

Ia441x4

cb5+a5tx,+a52x2+a53x3+1'1a441xa 3 =h5+ L asiXi+h4+ i=1

3

L a4,-ti, i=1

where the inequality holds due to the column diagonal dominance of the compartmental matrix A [see (1 .3)] . The last line above follows from (4 .14) when k = 4 (and this shows the first inductive step in the proof) . We continue to repeat the process using (1 .3) and (4.14); regrouping terms in the last expression yields a551X5Eb5+b4+(a51+a41)x,+(a52+a42)x2+ a5 1 13 1 43 a331 1 3 b5+b4+(a51+a41)x1+(a52+a42)x2+1'1a331x3 = b5+b4+(a5,+a41)x1+(a52+a42)x2+b3+a31x1+a32

where the last equality uses (4 .14) when k = 3 (this completes the second step of the induction) . Now rearrange terms on x 1 ,x 2 , and again apply dominance condition (1 .3) and (4 .14) (k=2) in that order to arrive at 5 Ia551 2 5c

5

L b i + L a f1 x 1. i-2

( 4 .15)

1-2

The first component equation of (1 .4) yields x1=b1/1a111 . In (4 .15), use this expression for x1, and then diagonal dominance of the first column of

183


A to get the desired result,

5 jassfxs

E bi .

∎

i=1

5.

UNIQUENESS OF EQUILIBRIUM POINTS

It is desirable to know a priori if there exists a unique equilibrium point to model (1 .1) . Sandberg [36] has given conditions that guarantee uniqueness of an equilibrium point if it exists, but he does not fully utilize the properties (1 .3) of the system's being compartmental . Define f,(x)=k, to be the ith component function on the right-hand side of (1 .1), and use u(x) =,9-(x,ao(x)) as before in (2-9). Then, using componentwise ordering, Sandberg's conditions become (1) If u < w but u * w, then u(u) > a(w); (2) If u > w and a i = w i for any particular i, then fi(u) 3 fi (w) . If these two conditions hold for all u and w in R'+, then there is at most one equilibrium point for the system . These conditions are satisfied for compartmental models of the form described in the following theorem . Note that a ij (x) is nondecreasing on (1 if a, j (u) F a,j(v) whenever u c v in 1) . Also, a,j(x) is donor-controlled if it is a function of x j alone . THEOREM 5.1 If compartmental system (1 .1) with matrix A(x)=-[aij(x)] is such that a ;j (x) is nondecreasing on fl = R ; for all i = 0,I_ ., n and j =1, 2, . . ., a, where i # j, if each a i ,(x) is donor-controlled, and if a oj (x) is nonzero for all j and x e A, then the system has at most one equilibrium point in fl . Proof. We begin with setting up Sandberg's first condition . Assume u 5 w but u # w. Then, for all j we have u j a„ j(u) < wj a„j (w) since a

nondecreasing function of its argument . For some i we have u ; < w i since u # w, and since a oi # 0 we have u ;a oi (u) < wi a o;(w) for that particular i. Thus,

o(u)=/3-(u,ao(u)) n 1 -

L j .1

n ujanj(u) >9 -

E

j-1

wjaoj(w) - O(w) .

184


For the second condition we assume u > w and u, = w; for some particular i . Note that since a ;,(x) is donor-controlled, a, j(u)u,=a j,(u j )u j = aj,(w,)wi=a,,(w)w1 . Thus, f,(u)=

L

a ,j (u)u j +bi

j-1

n =a,,(u)u,+

a ij (u)u j +b,

E, j=1 j, i

n 3 a ii (w) w,

+ L a ij (w) wj + bi j=1 j=i

n =~aij(w)wj+bi=f;(w) .

∎

j-1 Existence and uniqueness of the equilibrium point can be proven by other means such as the contractor mapping theorem [39] . The style of construction uses the identity g(x) -g(y)={A(y) 1 -A(x) -1 }b ={A(y)-1(A(x)-A(y)]A(x)-,}b along with bounds of the type in Section 4, IIA(z) -1 E U, and a Lipschitz condition on the coefficient matrix, IIA(x)-A(y)II 0) . we know that IIx(t)I i must approach some limit as t goes to infinity . Notice that with 0 e fl, Equation (2 .9) says that the total mass s(x) acts as a Lyapunov function for model (1 .1) and the equilibrium point x = 0 because s(z) = 0, s(x) > 0 if x * z, and AX) = 0 (x) = -(x,a e (x)) 50 . Hence, x=0 is a stable equilibrium point for system (1 .1) with b = 0 . We can conclude that all solutions starting in 0 remain bounded as time passes . If Q- 1 = fl, then the system has no traps and x = 0 is globally asymptotically stable . See Example 2.1 and Figure 1 when b = 0 . All solutions starting in 0 stay bounded and approach ii = 0 as I - 119, p. 5] . If traps exist and the system is open, then a particular compartment j can still "wash out," that is, xi( t) -+ x 7 = 0 as t -> w, provided there

189


is an open path from compartment j to the environment for all x e 11 [15]. If there are internal traps, the system can be redefined by including the traps as part of the environment [22] . The "reduced" system again has the asymptotically stable equilibrium point x=0 regardless of x(0) . If the system is closed with no internal traps, then there is one steady state x, with I1xllt = IIx(0)Ih [22] . For further comments on stability when input b = 0, see Theorem 4 .1 in (151, [22], Section 4 in [31], (41, and Theorems 1 and 2 of [26] . Another approach to the question of stability is by linearization of the system . For this, we need the Jacobian of the function fix)- A(x)x + b from 11 into R" . We observe the following from the product rule for derivatives : If A(x) is then x n matrix of system (1 .1), then the i, j element of then X n Jacobian f(x) is a ri (x) 4 m-1

x m . ax ~atm(x)1'

(7 .1)

Again suppose input b=0 and z is an equilibrium point in fl- t that contains 0 . Then 0 = z = A(x)x, so z = A(x) -t0 = 0 . Thus z must be 0. Observe that f'(0)=[a11(0)]=A(0) by the above Jacobian equation. The linearized system for z=f(x) about x=0 is h(t)=f(0)h(1)=A(0)h(t) . Since x = 0 E fl_ t, all real parts of the eigenvalues of the constant compartmental matrix A(0) are negative . Hence x = 0 is locally asymptotically stable. Example 7.1 . Consider the nondimensional competition model for interacting populations [34, p . 79] of form (1 .1) with b = 0, positive parameters p, y, 3, and

A(x)

1_X1_yx 2

0

0

PO -x 2 -Sx,)

Moreover, suppose y > 1 and 3 < 1 . For A(x) to be compartmental requires 1 c x t +yxz and I c Sx t +x 2 . The compartmental domain is Sl={xeR . :1 w is his hypothesis H5 . These three hypotheses are sufficient [36, p . 276] to yield the conclusion of Theorem 7 .1 . ∎ When input b * 0, there is at least one case of stability that can be further discussed. Let Fj (x,) > 0 be a nonnegative continuous scalar function of the single variable x l > 0, for each j = 1, 2, . . . , n, and suppose that the n X n real array [a,i ] is a constant compartmental matrix . Define the donor-controlled functions

a ;,(x)

--a,I F3 (xi )/xl ,

i,j=1,2,

,n,

of x =[x„ x 2 , . . ., x„ )T on the restricted domain 12 =(x a R, : x > 0). Then the n X n matrix A(x) =[a, t(x)] satisfies properties (1 .2) and (1 .3) and thus is a nonconstant compartmental matrix. For x e Q, let F(x)=-[F1(x,), F2(x2), . . ., F„(x„)]T . The model a = A(x)x+b of (1 .1) now can be written as i=AF(x)+b, where A=[a,1 ) . There are no periodic oscillations for trajectories of this system, and the solution of the model either is unbounded or tends to the equilibrium set [31, p . 38] .


8.

191

MORE ON MEAN RESIDENCE TIMES

Some properties of the mean residence times T; 1 (x) in -A(x) -t for model (1 .1) when x e 12_, were discussed in Section 4. We continue that discussion in the present section . Note that a ;J (x) is nonincreasing on 1, if a,J (u) > a 1J(v) whenever u < v in 11 . Also, a ;J (x) is receptor-controlled if it is a function of x; alone . THEOREM 8.1

If the fractional excretion rate coefficients a 0J(x) are nondecreasing [nonincreasing] functions of x c fl, and if all fractional transfer coefficients between compartments are constant, then for u c v in Q_ 1 , we have

frd (u)c Td (v)foralli,j] .

T;J(u)3,-1(v), i,j=1,2, . . .,n

Proof. For analysis of the mean residence times we return to the Neumann series (4 .2) . With E(x) =I-A(x)D(x)=[e ;J(x)], we have e;J(x)= 0 if i=j, and e ;J(x)=-a ;J(x)/a3J(x) if i#j . Thus, from the Neumann series, the general element of T(x) =--A(x) - ' is

1 Tii(x)= aii(x)I eiJ(x)

+ k=1

eik(x)ek,(x)+

k=1(=1

eil(x)elk(x)ekJ(x)+

. ..

I

if i # j, and

Tii

lI1 eik(x) ek7(x) (x) = la,(O + kEI

eil (x) elk (x) ekl (x) +

+

kLI

.. .

1= 1

If o c v in R, and the fractional excretion coefficients are nondecreasing, then a 0J (o) c a„1(v); thus, Iai1(u)I E l a ii (v)I for all j. No other components of A have changed . Also, for l c i, j c n and i# j, 1> e,

Iaii(u)I

Iaii(v)I

la;i(v)I

eiJ(v)

Thus, r; J (u) 3 T;J(v) for all i and j by the Neumann . -ties . The proof is similar when the a 0j (x) are nonincreasing . ∎ Example 8.1 . The combat model of Example 5 .1 with a01(x)=a+ax e and a02(X) =c+yx l satisfies the hypothesis of Theorem 8 .1 since

192


the a,,(x) are nondecreasing. In particular, if ucv in 1Z_i=R2, then diag(1 / a,1/c) = - A(O)-t-> - A(u)-'> - A(v)- i for this model . The next theorem is a fairly intuitive physical result . Basically, the theorem says that if a flow rate into a compartment is increased, that compartment's mass must increase ; if a flow rate from a compartment is increased, that compartment's mass must decrease . THEOREM 8.2 Let A(x) be an n x n compartmental matrix in which all fractional excredon coefficients are constant and all fractional transfer coefficients except one between compartments are constant . Let this one variable coefficient be denoted a,j(x) if mean residence times for the receptor compartment are desired, and denoted ati(x) if mean residence times for the donor compartment are desired . With this notation, an increase (decrease) in aij(x) causes an increase (decrease) or no change in rik(x) for all k . An increase (decrease) in aii(x) causes_ a decrease (increase) or no change in rik(x) for all k . Proof. We wish to determine how the change of one fractional transfer coefficient, either from compartment i or to compartment i, affects the mean residence times of the form rik(x) where 1 c k c n . For each fixed x E= 1I_ , the element rik(x) represents the average time that a particle that begins in compartment k spends in compartment i before exiting the system . Consider a particular but arbitrary particle that begins in compartment k . We wish to determine the total amount of time that this particle spends in compartment i before it exits the system . If at some time this particle is in compartment i, leaves compartment i, travels around the system, and then returns to compartment i, there has been no effect on the total time spent in compartment i . If, however, at some time our particle leaves compartment i and is subsequently excreted from the system, the total time spent in compartment i has just been terminated . In other words, we have no concern with how long the particle has been absent from compartment i ; we are concerned only with the probability that it will return . In fact, it is simpler to use a new time scale, ?, that ceases measurement whenever the particle leaves compartment i and resumes measurement when the particle returns to compartment i . Define xi as an indicator variable such that [13, 351 1 x~t) __ 0

if a particle is in compartment i at time I otherwise

Then define I as &f(t) 2=di . 0

193


Now let At" be any small increment of time during which a given particle is in compartment i . Define (pd),,, as the probability that a particle in compartment i will go directly to compartment m from i during Af . [Actually, (pd) m is the probability that a particle in compartment i will go directly to compartment m during time At" plus the probability that it will go from compartment i directly to compartment m twice during time A? plus the probability during time A"r that this will happen three times, etc .] Define (po) m (m * i) as the probability that a particle that is in compartment m will be excreted from the system before it returns to compartment i . Finally, define (po) as the probability that a particle in compartment i will be excreted from compartment i during time interval Ai . With these definitions and the basic multiplication principle from elementary probability, we have that (pd),,, - (po)m is the probability that during the time interval Ai a particle in compartment i will enter compartment m and subsequently be lost from the system before returning to compartment i . With these given definitions, the probability p that a particle that begins in compartment i at time r" = C t will still be in compartment i at time i=i i +Ai is

P'

l ` I

1m= mxr

(pd)m(p0)m l - (po) . f1

(8 .1)

This comes from the probability rules for complements and for independent events . Since p thus defined is a probability, it may be interpreted as the average over many particles . The probability that a given particle that begins in compartment i (at time t" = 0) is in compartment i at time "t =m Ai for a positive integer m is p' by the basic multiplication principle . The total expected time spent in compartment i by a particle that began in compartment i is Tu- lim

[Ar"(p+p 2

+p 3 +

Ai-0

by the expectation formula from probability . A particle that begins in compartment k # i either enters compartment i at some time or else never enters compartment i . If the particle enters compartment i, it first enters at "r = 0, and from then on all probabilities are just as if the particle had originally begun in compartment i . With the probability that a particle that begins in compartment k enters compartment i equal to 1-(po) k , we have Tik = [ i - ( po)k1Tir

194


Let us first examine the case where ai,(x) increases . All probabilities in (8 .1) and (8 .2) are constant except for (pd)1, which has increased . Thus, -, k(x) has decreased . Similarly, if a,i (x) decreases, T, k(x) increases . Now consider the case when a ii(x) increases. All probabilities in (8.1) and (8 .2) are constant except for (po),,, . Consider all possible paths beginning with compartment m, not entering compartment i, and ending with an excretion . Probability (po)„ is the sum over all such paths of the products of the probabilities of each transfer or excretion encountered along the paths . Along any particular path, if a, i (x) increases, none of the probabilities has changed until the path enters compartment j . At this point, the probability of returning directly to i has increased and the probability of any other option has decreased . Thus, the probability of any path beginning at compartment m and ending in an excretion before returning to compartment i has decreased . Therefore (po),,, has decreased, and TWO has increased . Similarly, if ajx) is decreased . r, k (x) decreases also . ∎ Example 8.2 .

Suppose the compartmental matrix of interest is

-0 .075 6+x,/10 A(x) _-

0 .075 6+ x t /10 0

0

0

-0 .08

0 .0064

0 .08

-0 .016

The hypothesis of Theorem 8 .2 is satisfied, with the only variable transfer coefficient being a 21(x 1 )=0 .075/(6+x 1 /10), which is donor-controlled . Now a 21(2)=0 .0121 81 .3=- 11(1) . In this case, - 11(0)=80 would be the greatest lower bound for all T, 1(x,) when x 1 , 0. Theorem 8 .2 also applies to the AIDS model in Example 4 .3 . 9.

ITERATION ON g(x)

In Section 5 we introduced iteration on the function g(x) _ - A(x) - 'b of (1 .6). The context there was when g acts as a contractor on a subset of fb _ 1 and thus the sequence x(p+1)=g(x(P)), p=0,1,2, . . ., converges to an equilibrium point x of model (1 .1) . In this section we discuss the recursive process in more depth . It provides an alternative to the usual method for computing steady states of (1 .1) in the case where the nonlinear algebraic equations in (1 .4) prove difficult to solve .

195


An equivalent but better computational form for this iteration is solution of the successive linear systems A(xt")xtp + p=-b,

p=0,1,2, . . .,

(9 .1)

with a prescribed x" . The function g maps fl-, into R : and maps S2_, into itself when 11_, = fl = R'+' . If g maps a subset of f2-, into itself, then iteration using g is well defined . There are a variety of convergence results that can be proved for iteration on g(x) . The proofs build directly on information found in Section 8 . Let fl„ be a subset of (1 for a given model (1 .1) on which all excretion rate coefficients a 0 ,i( x) of ,4(x) are bounded away from zero as in (4 .1). Then 0 0 c IL_, [3, 201 . In the following, x0 -- 0 is necessary for simplicity of both the theorems and their proofs . In practice, if a good estimate of x is known, it should be used as the starting value for the process (9 .1). THEOREM 9.1

If fl„ = R' , if the excretion coefficients ao,(x) are nonincreasing functions of x, and if a!! transfer coefficients between compartments are constant, then iteration (9 .1) with xt0) 0 converges monotonically to some equilibrium point x. If x is any other steady state of the model, then x 5 x. Proof. Let u < v in R n . By Theorem 8 .1, we have - A(u) - ' < - A(v) - ', and thus g(u) c g(v). Since 12„= R°, Theorem 3 .1 guarantees the existence of at least one equilibrium point . Let that equilibrium point be denoted x* . Let xt" =-0. Then x ( o) Cx* so xto =-g(x°)cg(x*)=x* . Also, x(j) >_ x (o) , so xOt >_X" ). Thus, 0 < X(0) < xtil < xl 2) < . . . < x* . Since this sequence is monotonically nondecreasing and bounded above by x*, it converges to some equilibrium point, say x . Let x be any other equilibrium point . Then g(x)=z. Since 0 c x, x ( ' ) = g(0)< g(x)=x . Similarly, xfD) < z for all xf') in the sequence begun by x 1o) -0. Thus, z = lim X ( P ) S z . P __

∎

Now, assume that xt0 t is chosen such that it is larger (in every component) than a known upper bound for g(x) . The upper bound for g(x) can be obtained by using inequalities in either Section 2 or Section 4 . Then we have xt l t > x* in the above proof, and x (0) a xto. A nonincreasing sequence with the lower bound x* is obtained . Thus, the sequence converges . If the answers from the two starting values are the same, then the equilibrium point is unique . Otherwise, this new sequence has found the largest equilibrium point .

196


Example 9.1 . An important type of nonlinearity in compartmental models is that in which the elimination rate is capacity limited and can be modeled by Michaelis-Menten kinetics [18] . In the human body, both hepatic and renal pathways are susceptible to saturation . Consider the two-compartment model with Michaelis-Menten elimination from the first compartment and linear transfer of material between compartments . The matrix A(x) for model (1 .1) would then be

V

A(x)= [-k21 k+x 1 k 21

k12

-k 12

The excretion coefficients a 0i (x) are nonincreasing, and the transfer coefficients between compartments are constant . Thus, all the hypotheses of Theorem 9 .1 are satisfied except SZO=R'+ . This is not satisfied, because a 02(x)=0 . The hypothesis SU0=R, was needed only to guarantee the existence of an equilibrium point . With this hypothesis deleted, the conclusion to the theorem would read : If there exists an equilibrium point to the system, then the iteration converges monotonically upward from x 101 =0 to some equilibrium point z . If z is any other steady-state solution, then z < z . With k 27 =1 .2, k12 =0.7, V=20, k=2, x 101 =[0,0], and b=[0.4,0 .2]x, we obtain from iteration (9 .1) : x 01 =[0 .0600,0 .3886], xt2t=[0 .0618,0

.3916] x ,

=[[0 .0619,0 .3917] x , x (3) x (' ) =[0 .0619,0 .39181',

where x 141 is the correctly rounded equilibrium point . The iterates increase monotonically in both components . In this example, IIA(x)II1 = . The second inequality in max{21a ; 1 (x)I- a o;(x)l 612 .4 =- L for all x>0 Section 2 guarantees that the equilibrium point must obey x1+12=111111=IIg(1)11,>

![

= 0 .0484 .

If b is increased such that 0=IIb11 1 > V, the state vector x(t) will increase without bound since by (2 .9),

T F1xi = Q(x) = R - k+x1 >R - V>0. Thus, there will be no equilibrium point . (See Example 2 .1 for a related discussion .)

197


If D O =R, [or if .ft_I=R, and there exists a steady state x to system (1 .1)], if the excretion coefficients a 01(x) are nondecreasing functions of x, and if all transfer coefficients between compartments are constant, then the iteration of (9 .1) with x (0) =-0 converges to an equilibrium point i, unless there exist distinct x ) ' ) and x (' * ' ) in the interval [0,g(0)] such that g(x'" I) )= x )' ). In this case, x is between x (' ) and xt'i"L 1 .

Proof. Let u < v in R ; . Then by Theorem 8 .1, we have -A(u) -1 > - A(v) -' or g(u) > g(v) . Since (1.0 = R„ there is at least one equilibrium point z > 0 for the system by Theorem 3 .1 . Choose x (0) --- 0 . Then x 101 < x. Thus, x 1 ' 1 = g(x )0)) > g(x) = x. In general, if x ( ') > z, then x ('+' ) = g(x )' )) < i, and if x" ) < i, then x (' ") >~ i . Thus, the iterates bracket x. Since d2_ I = R„ g maps dt_ 1 into itself. Thus, xm a R, so x (' ) > 0. With x ' 2 > x (0) = 0, we have x )3) < x ( ' ), and the iterates will continue to bracket inward toward i . Thus, convergence is guaranteed unless the iterates begin to alternate between two unequal values. (This proof would work from a general nonzero x(0) if x (0) '< x ( ') _< TO a Example 9.2 . The combat model of Example 5 .1 satisfies the hypotheses of Theorem 9 .2. The transfer coefficients, a 21 and a 12 , are both zero, while the excretion coefficients, a 01(x) = a + ax e > a > 0 and a 02 (x) = c + yx 1 > c > 0, are both increasing functions of their variables and bounded away from zero . With a = 0 .9, c = 0 .8, a = 0 .7, y = 0 .65, b, = 0.3, b 2 = 0.2, and x (o) =0, we obtain from iteration (9.1) the converging sequence x )1) =[0 .333,0 .250] T ,

x' 2 =[0 .279,0 .197] T ,

=[0 .289,0 .204] T , x (3) x ( _ [0 .288,0 .203]' . s)

=[ 0 .288,0 .202]`,

Here, x(5) is the equilibrium point correctly rounded . Notice the bracketing in every position . From the discussion in Example 5 .1, we know in advance that this equilibrium point must satisfy X I +x 2 =111111=11g(i)IIt

sllA(x) - '1611b111< min (a,c)

0 .625 .

The alternation between unequal values mentioned in Theorem 9 .2 can occur. Consider the one-dimensional compartmental model (1 .1) with excretion coefficient a01(x)=-x 3 +x 2 +10x+8 and constant input b=-40.

198


Then Ax) =a,,(x)=-a01(x) and g(x) =-A(x)-tb=40/{[(x+1)x+ 10]x +8} . The (unique) steady state solution is z = 1 .4483, but g(1)= 2 while g(2) = 1 . For any starting value smaller than I or larger than 2, the iterates bracket inward toward the (1,21 interval . For any starting value inside the interval, the iterates bracket outward toward the endpoints of the interval . In a similar way the following result can be proved using Theorem 8 .2 . THEOREM 9.3

If S} o = R, and if all the transfer coefficients and excretion coefficients of AU) are constant except for one transfer coefficient a,j(x) which is receptorcontrolled and nondecreasing in x„ then the iteration x (P' t) =g(x'Pt) con Verges from xt°t to z, monotonically in the ith position . 10 .

DISCUSSION

This article examines equilibrium points x and qualitative theory of compartmental models (1 .1) on domain 11 where the transfer coefficients and excretion coefficients are state-dependent . The domain Q on which system (1 .1) is compartmental usually depends on the parameters of the model as illustrated in numerous examples . The steady states in Q-, are fixed points of the function g(x) defined in (1 .6). Information on x(t), on the range of function g, on the location of any z, and on input is given in Sections 2 and 4 . Perhaps primary among these results is the inner product form, which centers on system inputs and excretions : p - (g(x), a o (x)) =0 . An explicit equation that any equilibrium point z of the model must satisfy is o(x)=0, where a(x) =/f-(x,ao(x)) . This functional a(x) is the time derivative of the total mass s(x)=-Ilslli in the system, and s(x) is a Lyapunov function for z =0 a 11 and system (1 .1) (when there is no input) . In relating (Ax) and s(x), the desirable stability condition is to have or Do) < 0 for all sufficiently large s(x) . Existence and uniqueness of equilibrium points have been discussed in three sections . An important case is when each excretion coefficient a 0e (x) is bounded away from zero . This is not an uncommon condition to have . For instance, it occurs in the chemostat model, the predator-prey model, the epidemic model, the combat model, and the AIDS model . When each a 0i (x) is bounded away from zero, an equilibrium point exists, and constant upper bonds on IJAW - tII are easily constructed [see (4 .6)] . Moreover, since any equilibrium point z E 11 _, =11 is a fixed point of the function g(x), it is natural to investigate the application of the contractor mapping theorem to g . In view of Theorem 5 .2, this construction can always be carried out if input b is sufficiently small . If the contractor mapping theorem can be applied, then the existence and uniqueness of z is guaranteed within a


199

subset of 12- 1 . This also motivates iteration on g(x) (see Section 9), because the sequence so generated will converge to this unique x from any starting point in the subset. Convergence may be slow unless the contractor constant in (5 .1) is relatively small . Models (1 .1) that satisfy condition (6 .1) have trajectories that approach the hyperplanc (6 .3) at the rate that exp(-yt)-r0 as t-w . Condition (6.1) is a demanding hypothesis, but there are models that meet the condition, for example, any linear time-invariant model it = A x + b with no traps, or any nonlinear model where all excretion coefficients a o , are constant and the same value. The chemostat model was shown to satisfy this asymptotic result. See Figure 3 . Proposition 6 .2 presents a simple test for rejecting condition (6.1) for a particular model . Stability of the equilibrium point x =0 for model (1 .1) is now well established when 0 e 12 and input b = 0 . Possibilities for further stability research include many cases with nonzero input . This may be complicated, because with b # 0 models of type (1 .1) can exhibit both stable and unstable equilibrium points (e .g ., the chemostat model) . However, we do have the stability results of Theorem 7.1 above, Theorem 3 in [361, and the application of Lyapunov functions as in [31] (for donor-controlled compartmental systems)The elements 7,,(x) of - A(x) - ' depend on the behavior of the transfer coefficients a ii (x) in A(x). The mean residence times matrix - A(x) -1 , x e 12 1 , has nice structure and properties as was shown in Sections 4 and 8 . Moreover, these results furnished the basis for bounds on x and for the convergence theorems of Section 9 . Additional research on qualitative theory for nonlinear compartmental models can proceed along many lines, but investigation of major types of nonlinear systems that appear in applications will probably have the best chance of yielding useful results . These would include mammillary, catenary, and cascading systems, donor- and receptor-controlled systems, control systems (where transfer from one compartment to another depends on a third "control" compartment), strongly connected systems, single-exit systems, closed systems, and open systems that wash out . REFERENCES I

D . H. Anderson, Iterative inversion of single exit compartmental matrices, Comput .

Biol . Med . 9 :317-330 (1979) . D. H . Anderson, Structural properties of compartmental models, Mall,, Biosci . 58:61-81 (1982) . 3 D . H . Anderson, Compartmental Modeling and Tracer Kinetics, Lecture Notes in Biomathematics, Vol . 50, Springer, New York, 1983. 4 R . Bellman, Topics in pharmacokinetics . I : Concentration-dependent rates, Math . Biosci. 6 :13-17 (1970). 2

200


5 E . Beltrami, Mathematics for Dynamic Modeling, Academic, Boston, 1987 . 6 A. Berman and R . 3 . Plemmons, Nonnegative Matrices in the Mathematical Sciences, Academic, New York, 1979 . 7 E . K. Blum, Numerical Analysis and Computation Theory and Practice, AddisonWesley, Reading, Mass . 1972 . 8 M . Braum, C . S . Coleman, and D . A . Drew, Fds ., Differential Equation Models, Springer-Verlag, New York, 1983 . 9 E. R . Carson, C. Cobelli, and L . Finkelstein, The Mathematical Modeling of Metabolic and Endocrine Systems, Wiley, New York, 1983 . 10 V . Capasso and G . Serio, A generalization of the Kermack-McKendrick deterministic epidemic model, Math . Biosci . 42 :43-61 (1978). 11 Y . Cherruault, Mathematical Modelling in Biomedicine, Reidel, Dordrecht, 1986 . 12 S . A. Colgate, J. M . Hyman, and E . A . Stanley, Los Alamos researchers model AIDS epidemic, SIAM News 22(3), May 1989 . 13 D . G . Covell, M . Berman, and C. Delisi, Mean residence time-theoretical development, experimental determination, and practical use in tracer analysis, Math . Biosci . 72 :213-244 (1984). 14 L . Edelstein-Keshet, Mathematical Models in Biology, Random House, New York, 1988 . 15 J . Eisenfeld, On washout in nonlinear compartmental systems, Math . Biosci . 58:259-275 (1982) . 16 D. Fife, Which linear compartmental systems contain traps?, Math . Biosci . 14 :311-315 (1972). 17 M . Ghil and S . Childress, Topics in Geophysical Fluid Dynamics : Atmospheric Dynamics, Dynamo Theory and Climate Dynamics, Springer-Verlag, New York, 1987 . 18 K . Godfrey, Compartmental Models and Their Applications, Academic, New York, 1983 . 19 l . Guckenheimer and P . Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcation of Vector Fields, Springer-Verlag, New York, 1983 . 20 J . Z. Hearon, Theorems on Linear Systems, Ann . N.Y. Arced . Sci . 108 :36-68 (1963)21 J . A . Jacquez, Compartmental Analysis in Biology and Medicine, University of Michigan Press, Ann Arbor, 1985 . 22 J . A. Jacquez, and C . P. Simons, Qualitative Theory of Compartmental Systems, Preprint, Univ. of Michigan, 1989 . 23 J . A. Jacquez, C . P . Simon, J . Koopman, L . Sattenspiel, and T . Perry, Modeling and analyzing HIV transmission : the effect of contact patterns, Math . Biosci . 92 :119-199 (1988) . 24 C . R . Kennedy and R. Aris, Bifurcations of a model diffusion-reaction system, in New Approaches to Nonlinear Problems in Dynamics, P . J . Holmes, Ed ., SIAM, Philadelphia, 1980 . 25 G . S . Ladde, Cellular systems . 1 . Stability of chemical systems, Math . Biosci . 29 :309-330 (1976). 26 G . S. Ladde, Cellular systems . 11 . Stability of compartmental systems, Math . Biosci . 30 :1-21 (1976) . 27 S . J . Leon, Linear Algebra with Applications, Macmillan, New York, 1980 . 28 C . C . Lin and L . A. Segel, Mathematics Applied to Deterministic Problems in the Natural Sciences, Macmillan, New York, 1974 . 29 E . Lorenz, Deterministic nonperiodic flow, J . Atmos. Sci . 29 :130-141 (1963). 30 N . H. McClamroch, State Models of Dynamic Systems, Springer.Verlag, New York, 1980 .


201

31 H . Maeda and S . Kodama, Qualitative analysis of a class of nonlinear compartmental systems: nonoscillation and asymptotic stability, Math . Biosci . 38:35-44 (1978), 32 H . Maeda and S . Kodama, Some results on nonlinear compartmental systems, IEEE Trans. Circuits Syst . CAS-26:203-204 (1979). 33 H . Marcus-Roberts and M . Thompson, Eds ., Life Science Models, Springer-Verlag, New York, 1983 . 34 J . D . Murray, Mathematical Biology, Springer-Verlag, Berlin, 1989. 35 P. Purdue, Stochastic compartmental models: a review of the mathematical theory with ecological applications, in Compartmental Analysis of Ecosystem Models, J. H . Matis, B . C . Patten, and G . G. White, Eds ., Int . Co-op Publishing House, Brutonsville, Md ., 1979 . 36 1. W . Sandherg, On the mathematical foundations of compartmental analysis in biology, medicine, and ecology, IEEE Trans . Circuits Syst . CAS-25 :273-279 (1978). 37 1. W . Sandberg, A note on the properties of compartmental systems, IEEE Trans. Circuits Syst . CAS-25 :379-380 (1978). 38 L. A . Segel, Modeling Dynamic Phenomena in Molecular and Cellular Biology, Cambridge Univ. Press, Cambridge, 1984 . 39 G . F . Simmons, Introduction to Topology and Modern Analysis, McGraw-Hill, New York, 1963 . 40 K . W . Thornton and R . J . Mulholland, Lagrange stability and ecological systems, J. Theor . Biol. 45 :473-485 (1974).

A short parameter-fitting routine for compartmental models.

Products of Compartmental Models in Epidemiology.

Animal-to-human extrapolation using compartmental models.

Nonlinear bending models for beams and plates.

Estimating the electrotonic structure of neurons with compartmental models.

Robust global identifiability theory using potentials--Application to compartmental models.

Phase-field-crystal models and mechanical equilibrium.

Manifold learning for object tracking with multiple nonlinear models.

Methods for comparing Salmonella mutagenicity data sets using nonlinear models.

Raoult's law revisited: accurately predicting equilibrium relative humidity points for humidity control experiments.

Consistent Partial Least Squares for nonlinear structural equation models.

Comparison of linear and nonlinear models for human population dynamics.

Nonlinear age-dependent models for prediction of population growth.

Transform-both-sides nonlinear models for in vitro pharmacokinetic experiments.

Functional Nonlinear Mixed Effects Models for Longitudinal Image Data.

Disturbed nonlinear multispecies models in ecology.

Simulation of compartmental models for kinetic data from a positron emission tomograph.

Multi-compartmental biomaterial scaffolds for patterning neural tissue organoids in models of neurodevelopment and tissue regeneration.

Compensation and stability in nonlinear matrix models.

Some epidemiological models with nonlinear incidence.

Decomposition-based qualitative experiment design algorithms for a class of compartmental models.

Switching points of view in spatial mental models.

Dynamic equilibrium under periodic perturbations in simple ecosystem models.

Discovering and quantifying nontrivial fixed points in multi-field models.