- Email: [email protected]
Normal-mode identifiability analysis of linear compartmental systems in linear stages
Normal-Mode linear
Identifiability
Compartmental
Analysis
of
Systems in linear
Stages
J. P. NORTON Department of Electronic and Electrical Engineering, Uniuersi& of Birmingham, P. 0. Box 363, Birmingham B15 2TT, England. Received 28 Februav
1979; reoised 29 October 1979
ABSTRACT The identifiability of linear, time-invariant compartmental models from specified experiments with noise-free observations is analyzed in terms of normal modes. The formulation of the analysis reduces the problem to that of checking if a solution by a succession of linear stages is possible. This allows the information available from a proposed experiment and from prior knowledge of model parameters, primarily those taken as zero, to be checked easily for redundancy, and unique (global) identifiability established or disproved. The analysis covers any combination of impulse and step forcing, and observation of any set of state variables. Examples are given, and the new method is compared with those based on the transfer-function matrix and the Markov parameter matrix. The features giving rise to nonunique but distinct models (local identifiability) are also explained in terms of normal modes, and illustrated by examples.
1.
INTRODUCTION
The conditions for compartmental models of a system with linear, time-invariant dynamics to be identifiable from specified experiments have been examined in many recent papers. In particular, the problem of identifiability from error-free observations has received a great deal of attention since its discussion by Bellman and Astrom [ 11. Some justification for the introduction of yet another technique is therefore needed. It will be provided towards the end of this introduction, when the problem has been described and existing methods for investigating identifiability have been reviewed. The system is described by a state equation a(t)=Ax(t)+Bu(t),
the deterministic
x(0) = 0,
aspects of
(1)
relating state variables x(t), one per compartment, to forcing &r(t), normally an impulse or a step input into one compartment. An algebraic MATHEMATICAL
BIOSCIENCES
50:95-115
OElsevier North Holland, Inc., 1980 52 Vanderbilt Ave., New York, NY 10017
95
(1980)
0025-5564/80/0500!35
+ 2 1$02.25
96
observation
J. P. NORTON
equation Y(l)
=
Cx(t>
specifies the observed outputs y(t), virtually always separately observed state variables, and their measurement gains. With the contents of compartment i used as state variable x,(t), the elements a0 of the state matrix A are related to the transfer rate constants k, from compartment j to compartment i and k, from compartmentj to the environment by [2] a,=k,, au = - k, -
i=1,2 5
1=l
k,,
,..., p, j=
j=1,2
,..., p,
1,2 )...) p,
i#j,
(3)
(4)
i+j
so identification of the rate constants is equivalent to identification of the state matrix. The structure of the model is defined by the number of compartments, p; by the transfers between compartments and to the environment assumed to occur, as indicated by the rate constants not taken as zero; and by the input and output connections specified by B and C. Structural identifiability has usually been defined so as to avoid the separate problem of choosing the input waveforms. Bellman and Astrom [l] point out the fact that if a linear, time-invariant system is identifiable at all, it is identifiable from its impulse-response matrix. Identifiability analyses therefore tend either to address the problem of identification once the impulse-response or transfer-function matrix is known [3-51, or achieve the same end by defining identifiability with reference to output behavior over a whole class of input forcing functions [6-71. By contrast, the technique introduced in this paper aims to test specifically the identifiability of a model with given input waveforms. To keep the exposition simple while still covering cases of practical importance, the forcing will be assumed to be by impulses, step functions, or both. A number of analyses of structural identifiability have related the unknown rate constants either to the coefficients of the transfer function or the Laplace transforms of the observed outputs [ 1, 3, 81 or directly to the amplitudes and exponents of the normal modes appearing in the observations [9- 111. In both approaches the equations in the rate constants are nonlinear. Because of their nonlinearity, it is sometimes hard to see whether redundancies exist or nonunique solutions arise. The importance of checking for redundancy is pointed out by Delforge [12]. Grewal and Glover [6] define a model as identifiable at a point in parameter space where the corresponding output, for some admissible initial
IDENTIFIABILITY
ANALYSIS
IN LINEAR
STAGES
97
state and forcing, is distinguishable from that corresponding to any other point. The resulting identifiability test is a rank test on the Jacobian matrix of the Markov parameter matrix
While the method is conceptually attractive, the Jacobian is large for models of modest size. For instance, to test the identifiability of five free parameters in a three-compartment model with two observed states and one compartment forced directly, one must compute and test a 12 X 5 Jacobian. Walter et al. [lo] propose an alternative with a smaller matrix to test, but the algebra involved in obtaining it is quite complicated even for a threecompartment model. To test whether a model is globally identifiable by any of these methods, an algebraic trial solution must be carried out, either explicitly or implicitly in the course of checking whether fixing the parameters gives a unique response. Even for a parametrization as simple as that by the rate constants, (3) and (4) the relations between the parameters and the observations are not simple and consistent enough from case to case to make a globalidentifiability test a routine matter. An alternative is to derive from the observations a set of equations for rows of the modal matrix (the matrix whose columns are the eigenvectors of the system) and columns of its inverse. Together with the eigenvalues, which can be observed directly, these two matrices determine the state matrix and thence all the rate constants. Introducing the modal matrix and its inverse increases the number of unknowns, but as shown later it simplifies the form of the equations to the point where the presence of redundancy and the effects of prior knowledge of some rate constants can be seen quite readily. To be specific, the analysis reduces to the examination of sets of linear equations for solubility. Because of this linearity, the analysis quickly confirms or refutes the ability of a given experiment to identify a given model uniquely (globally) and makes clear the effects of modifications to the model or experiment if the original proposal is inadequate. Furthermore, the reasons for local but not global identifiability-that is, isolated but not unique parameter values-become obvious during the analysis. Examples will be used to illustrate these features.
98
2.
J. P. NORTON
RELATIONS BETWEEN AND EIGENVALUES
RATE
CONSTANTS,
EIGENVECTORS
For the linear, time-invariant system given by (1) the eigenvalues X, and, up to an arbitrary scaling factor, the eigenvectors m, are defined by Am, =himi, so, collecting we have
the eigenvectors
i=1,2 >...>p>
together
as columns
(6)
of the modal matrix M,
(7)
AM=MA
where A is the diagonal matrix with ith principal diagonal element hi. Here the eigenvalues are taken to be distinct. Smith and Mohler [13] discuss conditions for distinctness, but the assumption is not in practice restrictive. For instance, systems with two or more zero eigenvalues, corresponding to pure sinks, are excluded. They may nevertheless be handled by lumping the flows to the sinks with those to the environment and testing the identifiability of the rest of the model, then noting if enough sink states are observed for their inflow rate constants to be identified. Nonzero repeated eigenvalues can be treated as differing by very small known amounts. The eigenvalues are also real and nonpositive because the rate constants are nonnegative and the system is stable. Denoting by N the inverse of M (which exists, since the eigenvalues are distinct), the rate constants between compartments are given by ku = element (ij)
of MAN = rrRnj, i=1,2 )...) p,
j=1,2
)..., p,
i#j,
where r: and nj are row i of M and column j of N. The rate constants environment are k,=
-
5 ao= i=l
-
5 r’Anj,
(8) to the
j= 1,2 ,...,p.
i-1
These expressions are used both to find the unknown rate constants once the rows of M and columns of N have been found, and to incorporate any prior knowledge of rate constants. By definition of N, (8) and (9) are accompanied by r,Tnj= S,, where 8, is the Kronecker
i= 1,2 ,...,p,
j= 42, . . . .p,
delta, zero except for i=j
when it is one.
(10)
IDENTIFIABILITY
ANALYSIS
IN LINEAR
99
STAGES
The response of each state variable to any combination of impulse and step-function forcing can, as in the next section, be written readily in terms of A, M, and N, making it clear what information is provided by the experimental choice of inputs and observations. The fact that the observations allow eigenvectors and eigenvalues of the state matrix and hence, in some cases, the matrix itself to be found has been noted by Walter et al. [5], but they do not point out the possibility of redundancy or the need for systematic examination of prior knowledge to determine when it renders the model identifiable. Because (8) (9) and (10) are bilinear in ri and n/, it is simple to detect redundancy arising when they are to be solved to give rows of M and columns of N one at a time, each stage consisting of the solution of p equations linear in the unknown vector. As seen below, A is found immediately from the observations. 3.
RESPONSE
IN TERMS
OF NORMAL
MODES
Consider the response of all the state variables to unit-impulse forcing of one compartment j. Responses of the observed variables to this or more complicated forcing are then found by scaling, superposition, and (if necessary) integration. In general x(t) = e”‘x(O)+
JteA(‘-‘)Bu(7)d7,
(1’)
0
so for zero initial conditions and forcing v,6(t), where vj is zero but for 1 in row j and S(t) is a unit impulse, the response is x(t) = eA’vj = MeA’Nvj = Me%,,
t >o,
(‘2)
so P
xi(t) = rTeAfnj =
2 q--l
miqntie~*.
(‘3)
Hence A is obtainable from the exponents of any observed variable containing all the normal modes. Failure of some modes to appear is discussed later. However many state variables are observed separately, each gives p normal mode coefficients in terms of nj and one ri. If nj or any one ri is scaled at will, all the others follow, providing the measurement gains are known. More generally, impulse forcing @(t) of compartments j~Ji and step forcing 4 S&-s(t) dt of compartments j EJ~ gives a response in compartment i
100 The constant normal-mode
J. P. NORTON
term in xi gives no information, coefficients E,q
as we have by summing
the
(‘5) On the right-hand side, the first term is known and the second is identical to the constant term. Redundancy occurs in thep equations for Eiq either if all the state-matrix elements uij, ~EJ,, (perhaps none) are known in advance-for, using (7) and (lo),
-or if there are no step inputs, so that the second term on the right-hand side of (15) vanishes. It is important to note that, as a consequence of these possible redundancies, the identifiability of a particular model may be altered by the choice of input signals. This fact is obscured by the assumption, in the usual definition of identifiability, that the transfer-function or impulse-response matrix is available. Together (8) (9) (lo), and (14) specify all the information available after an identification experiment. The test for unique identifiability involves first noting the rows of M and columns of N found directly from the normalmode coefficients, then counting the equations linear in the remaining rows and columns, deieting any not usable in successive linear-equation solutions for the unknown rows and columns because of redundancy as explained later. If there remain at least as many such linear and independent equations as unknowns, the model is uniquely identifiable. If not, any redundant set of the equations gives rise to a further equation obtained by setting the determinant of their coefficients to zero. Such an equation may be linear in an unknown row or column and allow the solution to proceed. When a solution in successive linear stages is not possible, it is usually straightforward to see what change to the model or experiment would help. The following two-compartment example illustrates the basic ideas. Example 1. Compartment 1 of a two-compartment system is perturbed by a unit-impulse input, and its response is observed. The transfer rate constant k,, from compartment 1 to the environment is known to be zero. The measured response of compartment 1 is
c,x,(t) where ci is the measurement
= c,r~e”‘n,,
gain. The observed
time constants
are xi and
IDENTIFIABILITY
ANALYSIS
IN LINEAR
101
STAGES
h,, the diagonal elements of A. Because the scaling of each eigenvector column of M is arbitrary, rT can be taken for convenience as [ 1 l] provided it contains no zeros. This is so, as both the normal modes appear in the observed state variable. If c, is known, n, follows from the measured amplitudes of the exponential normal-mode components of x,(t). The equations in the remaining unknowns r2 and nz are rTn2 = 1
r:n, = 0,
rrn2 = 0,
from (10); and from prior knowledge
of k,,, using (9), we have
(rl + r,)=An,
= 0.
With r, and n, known already, r, can be found uniquely r2‘[ nl
An,]=[O
from
-rTAn,].
At this stage, M is completely known, and the analysis has shown that the model is globally identifiable. If numerical values of the rate constants were being computed, n, would follow from
[I [I 6 6
and the unknown
rate constants
O
n2=
1 ’
from
k,2 = rrhn,, k,, = r,Thn,,
k,,=
-(r,
+r,)*An,.
Comparisons between methods of analysis on such a small example are perhaps pointiess, but it is worth noting that use of the transfer function requires algebraic formation and then partial-fraction expansion of C[sZA]-‘B, looking out for pole-zero cancellation while doing so, and then examining the resulting equations to see if they are soluble for the rate constants. The Markov-parameter-matrix approach requires formation and differentiation with respect to each unknown of CA ‘B for values of i up to 2p - 1, and then examination of the Jacobian matrix to see the conditions for full rank. The modal solution above is not complicated compared with these. If c, had been unknown in Example 1, the proportions of the elements of II, would have been given by the ratio of the normal-mode amplitudes, and
102
J. P. NORTON
the scaling by r:n, = 1 from (10). This latter equation is redundant when c, is known beforehand. Since the effects of unknown measurement gains are straightforward to assess, they will be taken for brevity as known and unity in what follows. Had k,
rather than k,, been zero, the prior-information
equation would
have been (r,
+r2)TAn2=0,
and instead of solving linear equations
immediately
to hand for r2,
would have been used to give
I
6
I
=o.
(rl +r21TA
since n, is not entirely zero. This is a linear equation in r2, of the form
f Tr2= g, where f and g are known, depending
only on r, and A, and g is nonzero.
Hence r2 is found from
[ 32=[
3.
The rest of the solution is as before. Two types of special case exist. The first is when not all normal modes appear in the observed state variables. It causes no great difficulty, as seen in Sec. 4. The second, and the harder to analyze, is when definite but nonunique solutions are obtainable. One such example is given by Skinner et al. [14] and several by Jones and Godfrey [15], using the nonlinear equations in the rate constants given by the coefficients in the rational polynomial Laplace transforms of the observed state variables. The mechanism by which nonunique solutions arise is not made clear by the Laplace-transform analysis. It is discussed in terms of eigenvectors and eigenvalues
in Sec. 7.
IDENTIFIABILITY
4.
ANALYSIS
OBSERVABILITY
IN LINEAR
STAGES
103
AND CONTROLLABILITY
The object of this section is just to examine the implications of incomplete observability and controllability on the workability of the normalmode analysis, not to give a general review of their relevance to identifiability. Pohjanpalo [16] gives a simple check for the excitation of a given mode in a particular state variable. Writing the matrix exponentials in (8) as h4e”‘N and Meh(‘-7)N, it is seen that if some mi9 is zero, normal mode q does not appear in xi(t), and if some nsi is zero, neither direct forcing of compartmentj nor an initial value for xi excites normal mode q; the system is partially unobservable from x,(t) and uncontrollable from inputj respectively. The occurrence of zeros in M or N poses no problem in testing identifiability by the normal-mode method. Failure of normal mode q to appear in xi as a result of forcing j does not indicate on its own which of miq and nti is zero. However, unless another observed state variable does include that mode, X, and consequently the system as a whole are unidentifiable. If the mode appears in any other observed state x,,, nd is not zero, so mi4 is zero and either mk or nti may be chosen freely to scale eigenvector q. Zeros in rows of M or columns of N that are not seen directly in the normal-mode coefficients but found later on in the calculations are as tractable as nonzero values in principle, and will generally occur as small nonzero values in practice, because of experimental and modeling error. With distinct eigenvalues, failure to excite mode q is indicated by row q of NB being zero, and failure of mode q to appear in the observed outputs implies that column q of CM is zero. These conditions are the modal corollaries of the rank criteria of Kalman for complete controllability and observability, to which a number of authors refer [13, 171. 5.
REDUNDANCY
DUE TO ZEROS IN THE STATE MATRIX
Returning to the general case, any extra equations needed to make the model identifiable are given by prior knowledge of rate constants, particularly knowledge that some are zero. Zero rate constants corresponding to zero entries in the state matrix are also the principal cause of redundancy in the bilinear equations to be solved for the rows of M and colunms of N. If a,* is known a priori,
(17) can be used with other equations either knowing n* to find r, or knowing rr to find n,. Taking the first possibility and noting that An, = col. s of AN =col. s of NA = 2 q-1
a,,n,,
(18)
104
J. P. NORTON
the coefficient of the unknown r, in (17) is seen to be a linear combination of columns of N. Apart from (17) there are p equations
Cn, = arq,
q=1,2
>...,p,
(19)
in r,, and perhaps others like (17) but with s different. To avoid redundancy in the set of p linear equations to be solved for r,, any one of (19) replaced by (17) must have q such that uqs is nonzero for at least one of the equations (17) included in the set. If not, none of the coefficients of r, is a function of nq, so the p equations are linearly dependent, being functions of at most p - 1 columns of N. Among other things this implies that u,~= 0 cannot be used in place of r;n, = 1. Similarly, in finding n, any equation from
rqTn= s,,, replaced equation
s=1,2 >...,p>
by (17) must have q such that a,, is nonzero (17) in the set to be solved for n,.
(20) for some such
Example 2. A 3-compartment model has k,, known, so aI3 is known. The responses of x, and xa to an impulse input into compartment 1 are observed, with known measurement gains. On choosing the scaling of rI or r3, n, is found immediately from the normal-mode coefficients mlqnql or msqnql in x, or x,; then r3 or r, follows from the normal-mode coefficients in the other observed state. It is proposed next to find n3 from
rT 6
I [I n3=
0 1
a13
:
rrh
However, rrR=a,,rr+
a&+
a,,$,
so if a,* is zero, the third equation in n3 is implied by the other two. Here the redundancy is not affected by whether a,3 is zero or not, but if x2 rather than x3 had been observed, giving r, instead of r3, a zero value for aI3 would have made the set
redundant, whereas a nonzero been immaterial.
value would not. In that case, a,2 would have
IDENTIFIABILITY
Knowledge
ANALYSIS
IN LINEAR
of a rate constant
to the environment,
5 r:An,
105
STAGES
k,, gives
= - k,.
(2’)
q=l
This is potentially useful to help find any one of the rows of M. Redundancy arises just as before except that rTns = 1 can usually be replaced, since a,, is nonzero unless compartment s is a pure sink. The sole complication in checking the usefulness of (21) is that the check extends over all q. The other potential use for (21) to help find 9, is significant only if
i q=l
r,’
is known before all the rows of M are known individually, since otherwise n, can be found using (10) without recourse to (21). For this to be so, it must be possible to replace at least one of the equations E rFns= 1,
s= 1,2 )...) p,
(22)
q=l
[which follows from (lo)] with an equation (21) not dundancy. The analysis proceeds as before, but with
replacing any one of the rows of M, and the corresponding replaced by (22). 6.
TESTS FOR UNIQUE
introducing
re-
equations
(10)
IDENTIFIABILITY
All the experimental and prior information concerning the transfer rate constants appears as bilinear equations in rows of M and columns of N. They are linear in any one row or column, and if required more such equations may be obtainable, as in Example 1 with k,, zero, from the determinants of any consistent but redundant subsets of the bilinear equations. A unique solution is guaranteed if the equations can be solved in stages, each involving only collectingp equations in which all but one row or column have already been found. The selection of such a subset can be viewed as replacing any of the p equations (10) in the current unknown which also contains another unknown row or column by one from (17) or (21) which does not. The rows and columns given directly by the choice of eigenvector scaling are known at the start of this process.
106
J. P. NORTON
To test the adequacy of the information available in a given experiment, each equation (17) or (21) giving a transfer rate constant a priori can be checked for redundancy in each of its potential uses as a replacement for an equation (10). By this means one experiment is tested, but it is worth noting that for a realistic number of compartments (up to about three if the exponential components of each observed response are to be resolved with acceptable accuracy), it is quite feasible to set out exhaustively all completely linear solutions using every adequate collection of prior information for a given combination of forcing and observations. In other words, all the uniquely identifiable model structures compatible with the experimental conditions can be found. The factors which may disqualify from use an equation contributed by prior knowledge of a rate constant will now be listed. Consider first (17): r:An,= a,* = k,s.
(‘7)
r:n, = S,,
(‘9)
It may not be used to replace
in forming
a set of equations
for rr if
(i) x, is observed, so that r, is found immediately from the coefficients of the normal modes in x, or chosen to fix the scaling of the eigenvectors, and (17) is not needed; (ii) the only forcing input is into compartment q, so that n4 is found immediately from the normal-mode coefficients of an observed state variable or chosen to fix the eigenvector scaling (only if n4 were unknown would there be any point in removing r:n, = S,, from the set to be solved for r,; this disqualification does not apply if there are two or more forcing inputs); (iii) q = s, so that if n, is known and use of (17) is possible, nq is also known and it is not necessary to replace (19); (iv) Us=is zero, so that the replacement would introduce redundancy as explained in Sec. 5; or (v) q= t and u,~ is zero, so that redundancy would arise as explained in Sec. 5. In the same way, (17) may not replace rTns = S,,
(20)
IDENTIFIABILITY
ANALYSIS
IN LINEAR
STAGES
107
in finding n, if (i) the sole forcing is of compartment (ii) x4 is observed; (iii) q= t; (iv) a,, is zero; or (v) q = s and u,~ is zero. The small modifications when k, rather earlier. The possible redundancy explained in all a0 are known, with xi observed and Commonly there are no step inputs, so and the sum
s;
than k, is known have been noted Sec. 3 may be checked by noting if impulse forcing of compartment j. ri for each observed state variable
over the impulse inputs are found directly from the observed normal-mode coefficients. Consequently, if we rewrite (10) by replacing any one nj,j EJ,, with
then the rewritten
equations
(23) are redundant and are deleted. It is simplest to carry out all the checking for redundancy in terms of n; in such cases. The following examples show how global identifiability is tested by examining the feasibility of a unique solution for the modal matrix and its inverse. Example 3. A model used in pharmacokinetic experiments on the hepatobiliary system has 3 compartments: the plasma, liver, and biliary processes. An impulse input (rapid intravenous injection) forces compartment 1 directly; xi and xg are observed; and k,,, k,,, k,,, and k,, are taken as zero, defining the structure of the model. The forcing amplitude is 1, and the observation gains c, and cs are known. Choose r, to be a vector of ones; then n, follows from the coefficients of the exponentials in x,, and hence r, from those in xj. Check each known-
108 rate-constant
J. P. NORTON
equation
for usefulness:
(i) k,, gives (ri + r, + r,)Thn, = 0, which cannot be used as r;An, = known quantities to accompany r;n, = 0 and rln,= 1, because a3, is zero and so An, is linearly dependent on n, and n2. Its only potential use is to accompany rln, = 0 and r:n, = 0 in finding r2. (ii) k,, =0 gives r;An, =O, which is useful only to replace rln,=O, i.e. to accompany rFn3 = 0 and rTn3 = 1 in finding n3. (iii) k,,=0 gives r:An, = 0, which cannot replace rcn, =0 or rcn, = 1 in finding r,, because doing so would give a redundant set of equations, and cannot replace an equation with unknown coefficients of n3 in finding n3, since both rl and r3 are already known, and replacing r:n3 = 0 is pointless. This item of prior information is therefore useless. (iv) k,,=0 gives rTAn, =O, completely redundant, as r3, A, and n, are already known. Seven equations (10) remain once the redundant rTn, = 1 and r:n, = 0 are removed, so with the two useful items of information k,, =0 and k,, =0, there are 9 equations in the 9 unknowns in rz, n2, and n3. Solution by linear stages gives n3 from
rT
I I [I rTh
0
0 ,
n,=
1
rT
then r, from
rT[4
An,]=[O
n3
0
-(r1+r3)TAn,].
By now M is known uniquely, so the model is uniquely identifiable. What is more, it is clear which items of prior information are essential and which are redundant. It is interesting to compare the analysis above with the transfer-function and Markov-parameter-matrix investigation of the same example. The transfer-function approach gives
Y(s)=
[ :,
=
;
1(s+k,,)(s +k,, +k,, +4112 element
[
=-[
q,,,,1-I[i] (1,l)
of [sl-A]-’
element (3,1) of [sl-A]-’
A
b,kz
IDENTIFIABILITY
ANALYSIS
IN LINEAR
= (s + k,,)(s - a)(s - P),
109
STAGES
say
The observed eigenvalues are - k,,, o, and fl. Here the attribution of -kc,, is unambiguous because it appears in x3(t) but not, because of pole-zero cancellation, in x,(t). The amplitudes of the observed exponential components are the partial-fraction coefficients of X,(S) and X,(s), i.e. the residues (Y, and /3. For X,(s) they are at - ka,, e II=
o+ke2+k,Z+k32 o--P
and
e _ P+koz+ki,+k,, 12-
P-a
’
and for X,(S) they are
Examination of the equations, with (Y, p, and k,, known, shows that the remaining rate constants can be found uniquely in the order k,,, k3*, k,, in the and k,,. It also becomes apparent that there are two redundancies coefficients e3,, e32, and e33, and one in e,, and e,2. In spite of the simplification due to - k,, being an eigenvalue, the algebra is not negligible, and some features do not emerge clearly. For instance, prior knowledge of k,, and k,, is shown by the modal analysis to be unnecessary, but a reworking (with very complicated algebra, as the reader is invited to confirmj is needed to verify this by the transfer-function method, and, more to the point, the original transfer-function analysis does not bring these redundancies to one’s attention. The Markov-parameter-matrix method requires the formation of the derivatives of terms up to CA’B with respect to the rate constants. Merely writing out CA3B for this example will convince the reader that the analysis is lengthy and its structure is obscured by algebraic detail. The next example is intended to show how a modal solution by linear stages may require new linear equations to be derived from singular subsets of the original equations in rows of M and columns of its inverse. Example 4. A 3-compartment model has k,, = k,, = 0, so that the peripheral compartments 2 and 3 interact with the central compartment 1 but not with each other. Compartment 1 does not excrete to the environment, so k,, = 0. A unit impulse forces compartment 1 and xi and x2 are observed. Eigenvector scaling and the amplitudes of the exponentials in xi and x2 determine r,, r,, and n,. There are 7 nonredundant equations (IO), 3
110
J. P. NORTON
prior-rate-information equations, and 9 unknowns in r3, n2, and n3. However, a linear solution is not immediately possible for any of these, as the only substitution which would produce a set of 3 equations in one unknown, r;An, = 0 for rrn,= 1, gives a redundant set
6 6 :
r;h
I [I 0 0 . 0
n3=
Instead, all redundant sets of 3 equations linear in the same row or column are examined. They consist entirely of equations with right-hand side zero, and give rise to
In1
n3
An,\=0
from equations
linear in rz,
(244
In1
n2
An,l=O
from equations
linear in r3,
(23
=o
from equations
linear in n,,
PC)
=o
from equations
linear in n,,
PD)
=o
from equations
linear in n3,
(29
6 6 (r,+r,+r,)=A
rT rT r:A rT rT r;A
Substitution of rf’n2=0 and r$n,= 1 into Eq. (24B) gives a quadratic in any element of 9, so there are two solutions for n,. One of them is easily seen to be a scaled version of n3, which appears in Eq. (24A) just as n, appears in Eq. (24B), and satisfies rrn3 =0, but has rTn,=O. Consequently there is in fact only one solution for n,. The nonredundant set r3‘[ nl
n2
An,]=[
0
0
-(r1+r2)TA]
then gives r,, so A4 is complete and the model is globally identifiable. that r3 could not have been found from the set
6b which is redundant.
n2
An2]=[o
o 01,
Notice
IDENTIFIABILITY
ANALYSIS
IN LINEAR
STAGES
111
A similar and quicker alternative is to use Eq. (24D) with r:n, = 0 and r;An, = - (r, + r,)TAn, to find r,; again two solutions arise, one of them being a scaled version of r, because Eq. (24C) is identical in form to Eq. (24D) and r:n, = 0 to rrn, = 0. With r3 known, M is complete. 7.
LOCAL IDENTIFIABILITY
When a particular model has been found not to be uniquely identifiable from a given experiment, it may still be locally identifiable, that is, there may exist a number of distinct solutions for some or all of the rate constants. If so, background knowledge or other experimental work may suggest which is the correct one. Nonunique solutions arise in two ways, through singular sets of equations or through incomplete coupling in the state matrix. The first of these is best understood by reference to an example. Example 5. A 3-compartment model has k,,, and kz3 zero. A unit impulse is injected into compartment 1, and the responses of compartments 1 and 3 are observed with known measurement gains. From these observations A, r,, r,, and n, can be found immediately, leaving 9 unknowns making up r,, n2, and n3. The usual checks show that knowledge of ko, allows substitution of r,TAn, = -(r, +r,)TAn, for r,Tn,= 1 or rTn3=0, depending hand, the equation
(25)
on the order of solution.
On the other
r:An, = k,, = 0 cannot be used immediately, as r, and II, are both unknown, and it is redundant when r, is unknown and A and n3 known, being implied by rTn, = 0 and rTn3= 0 from (lo), because An, = a,+, + a,+, + a33n3= a+,
+ a3,n,.
The three equations in r, could therefore not be solved even if n3 had previously been found, but the linear relation between An,, n,, and n3 gives an equation in ng:
b
n3
An,]=O.
With A and n, known, this is used together with rrn3 = 0 and rrn3 = 1 to give a quadratic in each element of n,. Either of the two resulting solutions for n, can be used in rTn3= 0 with (25) and r;n, = 0 to give r,; then n, follows from (10). As there are unique solutions for r,, r,, n,, and A, but two sets of
112
J. P. NORTON
consistent solutions for r2, n2, and n3, the equations (8) give k3, uniquely but give two solutions for each of the other rate constants. An alternative order of solution uses k,, = 0 and (lo), or equivalently
to obtain
which with (25) and r:n, =0 gives two solutions from (10).
for r,; then n, and n3 follow
A singular set of equations in one unknown row of M or column of N may also be obtained from a set of p + 1 consistent equations. If they are all in the same row and in only two columns, one known and the other unknown, or ali in the same column and in only two rows, one of them unknown, then a singular set of equations can be produced by scaling and subtraction. Example 6. In a 3-compartment model, k,, and k,, are both zero, and as the input is injected into compartment 1, n, is found immediately. The two zero rate constants give
and
and two more equations namely
in r, + r, + r,, nl, and n3 can be obtained
(r, +r,+r,)*n,=
1
(r,+rz+r3)Tn3=
1.
and
Together
these four equations
imply that
(rl+b+rdT[An,
An3 n,-n,]=[o
o 01,
from (lo),
IDENTIFIABILITY
ANALYSIS
IN LINEAR
STAGES
An,
n, -n3 1=O.
113
so that (An,
The other mechanism which may lead to a finite number of distinct solutions is the existence of a normal mode which appears in only one state variable, or a state variable which shows only one mode. Mode j appears only in xi(t) if and only if columnj of the state matrix A is zero but for uii, i.e., compartmentj loses material only to the environment. Correspondingly only element j is nonzero in eigenvectorj and I+ and (4) with either (6) or (18) gives A, = aii = - k,.
(26)
The rate constant identified thus is nonunique in cases where without observing any other state variable one cannot decide which of the observed eigenvalues is 4. Similarly if row j of A has only ui, nonzero, i.e., compartment j receives material only from outside the system, then rj and rowj of N have only elementj nonzero, and (26) is again true. Here also 4. may not be identifiable uniquely, for unless 3(t) is observed, one may not know which normal mode corresponds to h, and if xj(t) alone is observed, the other modes are not seen. system with k,, zero has compartment Example 7. A two-compartment 1 forced by an impulse and compartment 2 observed. If the forcing amplitude bi and observation gain cz are both known, free eigenvector scaling and the amplitudes of the observed exponentials give r, and n,, and the eigenvalues, both appearing in x2, give A. Because A has only one nonzero element in row 1 (and column 2), r, has only its first (and n, its second) element nonzero. Accordingly rrnl= 1 gives r,,= l/n,,. However, this value depends on which mode was labeled mode 1 in calculating n, 1, so there are two distinct solutions for r, (and nJ, and the model is only locally identifiable from this experiment. If also it is known that k,, is zero, giving r;An, = - r,TAn,
(27)
the positions of the zeros in A, M, and N are not altered, and r, is calculated as above. Now, however, only one solution for r,, is consistent with (27), which may be written r 11=
r22A2n2, -rzl-
7’ 1
11
J. P. NORTON
114 so this model is uniquely identifiable from the experiment, of the mode attribution being resolvable.
8.
as a consequence
CONCLUSIONS
By the use of normal modes and introducing the rows of the modal matrix and columns of its inverse as unknowns, both the response of a compartmental system to impulse and/or step forcing and any prior information on parameter values have been written as equations which are bilinear in the unknowns. In consequence it is relatively straightforward to check for redundancy and the feasibility of a stage-by-stage linear solution, and hence to determine whether the model is uniquely identifiable. The alternatives, use of Laplace-transform or time-domain expressions for the responses written directly in terms of the parameters (rate constants) to be identified, yield in general nonlinear equations ill suited to checking for redundancy or multiple solutions. The use of normal modes also makes it relatively easy to see what modifications, if any, are needed to bring about identifiability. The origins of nonuniqueness of identification are clarified by normalmode analysis. Because the procedure has been reduced to an algebraically simple and systematic one, computer checking of a wide variety of experiment model combinations for global and local identifiability may be useful. This work was carried out during a Visiting Fellowship at the University of Warwick. Comments by Mr. R. P. Jones and Dr. K. R. Godfrey of the University of Warwick and particularly Dr. R. F. Brown of the University of New South Wales are gratefully acknowledged. The comments of a reviewer have been most helpful in improving the presentation. REFERENCES R. Bellman and K. J. Astrom, On structural identifiability, Math. Biosci. 7:329-339 (1970). J. A. Jacquez, Comparfmental Analysis in Biology and Medicine, Elsevier, Amsterdam, 1972. J. J. Di Stefano III and F. Mot-i, Parameter identifiability and experiment design: thyroid hormone metabolism parameters, Amer. J. Physiol. 233:R134-R144 (1977). K. Glover and J. C. Willems, Parametrizations of linear dynamical systems: canonical forms and identifiability, IEEE Tram. Automatic Control 19:64(X645 (1974). E. Walter, G. Le Cardinal, and P. Bertrand, On the identifiability of linear state systems, Math. Biosci. 31:131-141 (1976). M. S. Grewal and K. Glover, Identifiability of linear and nonlinear dynamical systems, IEEE Trans. Automafic Control 21:833-837 (1976). C. Cobelli, A. Lepschy, and G. Romanin-Jacur, Identifiability of compartmental systems and related structural properties, Math. Biosci. 44: I- 18 (1979).
IDENTIFIABILITY
8 C. Cobelli identifiability 9 10
ANALYSIS and
G.
Romanin-Jacur,
of multi-input
12
the formulation of models, J. Delforge, The problem
15
16 17
Controllability, multi-output
Biosci. 11:203-247 (1971). E. Walter, P. Bertrand, and G. Le Cardinal, IEEE and nonlinear dynamical systems,” (1979). M. Berman
14
and
and R. Schoenfeld,
115
STAGES
IEEE Trans. Biomed Eng. 23:93-100 (1976). S. I. Rubinow and A. Winzer, Compartment
11
13
IN LINEAR
Invariants
observability
biological analysis:
and
structural
compartmental an inverse
problem,
systems, Math.
Comments on “Identifiability of linear Trans. Auromatic Control 24~329-331
in experimental
J. A&. Phys. 27: 1361-1370. of structural identifiability
data on linear kinetics of a linear
and
compartmental
system: solved or not?, Math. Biosci. 36: 119-125 (1977). W. D. Smith and R. R. Mohler, Necessary and sufficient conditions in the tracer determination of compartmental system order, J. Theoref. Biol. 57: l-21 (1976). S. M. Skinner, R. E. Clark, N. Baker, and R. A. Shipley, Complete solution of the 3-compartment model in steady state after single injection of radioactive tracer, Amer. J. PhysioI. 1961238-244 (1959). R. P. Jones and K. R. Godfrey, Mathematical of linear compartmental models, University Report 74, 1978. H. Pohjanpalo and B. Wablstrom,
analysis of the structural identifiability of Warwick Control Theory Centre
The uniqueness
of linear compartmental
Inrernat. J. Sysfewu Sci. 8:619-632 (1977). M. Milanese and G. P. Molino, Structural identifiability and pathophysiological information from the kinetics 26:175-190
(1975).
systems,
of compartmental models of drugs, Math. Biosci.
Identifiability
Compartmental
Analysis
of
Systems in linear
Stages
J. P. NORTON Department of Electronic and Electrical Engineering, Uniuersi& of Birmingham, P. 0. Box 363, Birmingham B15 2TT, England. Received 28 Februav
1979; reoised 29 October 1979
ABSTRACT The identifiability of linear, time-invariant compartmental models from specified experiments with noise-free observations is analyzed in terms of normal modes. The formulation of the analysis reduces the problem to that of checking if a solution by a succession of linear stages is possible. This allows the information available from a proposed experiment and from prior knowledge of model parameters, primarily those taken as zero, to be checked easily for redundancy, and unique (global) identifiability established or disproved. The analysis covers any combination of impulse and step forcing, and observation of any set of state variables. Examples are given, and the new method is compared with those based on the transfer-function matrix and the Markov parameter matrix. The features giving rise to nonunique but distinct models (local identifiability) are also explained in terms of normal modes, and illustrated by examples.
1.
INTRODUCTION
The conditions for compartmental models of a system with linear, time-invariant dynamics to be identifiable from specified experiments have been examined in many recent papers. In particular, the problem of identifiability from error-free observations has received a great deal of attention since its discussion by Bellman and Astrom [ 11. Some justification for the introduction of yet another technique is therefore needed. It will be provided towards the end of this introduction, when the problem has been described and existing methods for investigating identifiability have been reviewed. The system is described by a state equation a(t)=Ax(t)+Bu(t),
the deterministic
x(0) = 0,
aspects of
(1)
relating state variables x(t), one per compartment, to forcing &r(t), normally an impulse or a step input into one compartment. An algebraic MATHEMATICAL
BIOSCIENCES
50:95-115
OElsevier North Holland, Inc., 1980 52 Vanderbilt Ave., New York, NY 10017
95
(1980)
0025-5564/80/0500!35
+ 2 1$02.25
96
observation
J. P. NORTON
equation Y(l)
=
Cx(t>
specifies the observed outputs y(t), virtually always separately observed state variables, and their measurement gains. With the contents of compartment i used as state variable x,(t), the elements a0 of the state matrix A are related to the transfer rate constants k, from compartment j to compartment i and k, from compartmentj to the environment by [2] a,=k,, au = - k, -
i=1,2 5
1=l
k,,
,..., p, j=
j=1,2
,..., p,
1,2 )...) p,
i#j,
(3)
(4)
i+j
so identification of the rate constants is equivalent to identification of the state matrix. The structure of the model is defined by the number of compartments, p; by the transfers between compartments and to the environment assumed to occur, as indicated by the rate constants not taken as zero; and by the input and output connections specified by B and C. Structural identifiability has usually been defined so as to avoid the separate problem of choosing the input waveforms. Bellman and Astrom [l] point out the fact that if a linear, time-invariant system is identifiable at all, it is identifiable from its impulse-response matrix. Identifiability analyses therefore tend either to address the problem of identification once the impulse-response or transfer-function matrix is known [3-51, or achieve the same end by defining identifiability with reference to output behavior over a whole class of input forcing functions [6-71. By contrast, the technique introduced in this paper aims to test specifically the identifiability of a model with given input waveforms. To keep the exposition simple while still covering cases of practical importance, the forcing will be assumed to be by impulses, step functions, or both. A number of analyses of structural identifiability have related the unknown rate constants either to the coefficients of the transfer function or the Laplace transforms of the observed outputs [ 1, 3, 81 or directly to the amplitudes and exponents of the normal modes appearing in the observations [9- 111. In both approaches the equations in the rate constants are nonlinear. Because of their nonlinearity, it is sometimes hard to see whether redundancies exist or nonunique solutions arise. The importance of checking for redundancy is pointed out by Delforge [12]. Grewal and Glover [6] define a model as identifiable at a point in parameter space where the corresponding output, for some admissible initial
IDENTIFIABILITY
ANALYSIS
IN LINEAR
STAGES
97
state and forcing, is distinguishable from that corresponding to any other point. The resulting identifiability test is a rank test on the Jacobian matrix of the Markov parameter matrix
While the method is conceptually attractive, the Jacobian is large for models of modest size. For instance, to test the identifiability of five free parameters in a three-compartment model with two observed states and one compartment forced directly, one must compute and test a 12 X 5 Jacobian. Walter et al. [lo] propose an alternative with a smaller matrix to test, but the algebra involved in obtaining it is quite complicated even for a threecompartment model. To test whether a model is globally identifiable by any of these methods, an algebraic trial solution must be carried out, either explicitly or implicitly in the course of checking whether fixing the parameters gives a unique response. Even for a parametrization as simple as that by the rate constants, (3) and (4) the relations between the parameters and the observations are not simple and consistent enough from case to case to make a globalidentifiability test a routine matter. An alternative is to derive from the observations a set of equations for rows of the modal matrix (the matrix whose columns are the eigenvectors of the system) and columns of its inverse. Together with the eigenvalues, which can be observed directly, these two matrices determine the state matrix and thence all the rate constants. Introducing the modal matrix and its inverse increases the number of unknowns, but as shown later it simplifies the form of the equations to the point where the presence of redundancy and the effects of prior knowledge of some rate constants can be seen quite readily. To be specific, the analysis reduces to the examination of sets of linear equations for solubility. Because of this linearity, the analysis quickly confirms or refutes the ability of a given experiment to identify a given model uniquely (globally) and makes clear the effects of modifications to the model or experiment if the original proposal is inadequate. Furthermore, the reasons for local but not global identifiability-that is, isolated but not unique parameter values-become obvious during the analysis. Examples will be used to illustrate these features.
98
2.
J. P. NORTON
RELATIONS BETWEEN AND EIGENVALUES
RATE
CONSTANTS,
EIGENVECTORS
For the linear, time-invariant system given by (1) the eigenvalues X, and, up to an arbitrary scaling factor, the eigenvectors m, are defined by Am, =himi, so, collecting we have
the eigenvectors
i=1,2 >...>p>
together
as columns
(6)
of the modal matrix M,
(7)
AM=MA
where A is the diagonal matrix with ith principal diagonal element hi. Here the eigenvalues are taken to be distinct. Smith and Mohler [13] discuss conditions for distinctness, but the assumption is not in practice restrictive. For instance, systems with two or more zero eigenvalues, corresponding to pure sinks, are excluded. They may nevertheless be handled by lumping the flows to the sinks with those to the environment and testing the identifiability of the rest of the model, then noting if enough sink states are observed for their inflow rate constants to be identified. Nonzero repeated eigenvalues can be treated as differing by very small known amounts. The eigenvalues are also real and nonpositive because the rate constants are nonnegative and the system is stable. Denoting by N the inverse of M (which exists, since the eigenvalues are distinct), the rate constants between compartments are given by ku = element (ij)
of MAN = rrRnj, i=1,2 )...) p,
j=1,2
)..., p,
i#j,
where r: and nj are row i of M and column j of N. The rate constants environment are k,=
-
5 ao= i=l
-
5 r’Anj,
(8) to the
j= 1,2 ,...,p.
i-1
These expressions are used both to find the unknown rate constants once the rows of M and columns of N have been found, and to incorporate any prior knowledge of rate constants. By definition of N, (8) and (9) are accompanied by r,Tnj= S,, where 8, is the Kronecker
i= 1,2 ,...,p,
j= 42, . . . .p,
delta, zero except for i=j
when it is one.
(10)
IDENTIFIABILITY
ANALYSIS
IN LINEAR
99
STAGES
The response of each state variable to any combination of impulse and step-function forcing can, as in the next section, be written readily in terms of A, M, and N, making it clear what information is provided by the experimental choice of inputs and observations. The fact that the observations allow eigenvectors and eigenvalues of the state matrix and hence, in some cases, the matrix itself to be found has been noted by Walter et al. [5], but they do not point out the possibility of redundancy or the need for systematic examination of prior knowledge to determine when it renders the model identifiable. Because (8) (9) and (10) are bilinear in ri and n/, it is simple to detect redundancy arising when they are to be solved to give rows of M and columns of N one at a time, each stage consisting of the solution of p equations linear in the unknown vector. As seen below, A is found immediately from the observations. 3.
RESPONSE
IN TERMS
OF NORMAL
MODES
Consider the response of all the state variables to unit-impulse forcing of one compartment j. Responses of the observed variables to this or more complicated forcing are then found by scaling, superposition, and (if necessary) integration. In general x(t) = e”‘x(O)+
JteA(‘-‘)Bu(7)d7,
(1’)
0
so for zero initial conditions and forcing v,6(t), where vj is zero but for 1 in row j and S(t) is a unit impulse, the response is x(t) = eA’vj = MeA’Nvj = Me%,,
t >o,
(‘2)
so P
xi(t) = rTeAfnj =
2 q--l
miqntie~*.
(‘3)
Hence A is obtainable from the exponents of any observed variable containing all the normal modes. Failure of some modes to appear is discussed later. However many state variables are observed separately, each gives p normal mode coefficients in terms of nj and one ri. If nj or any one ri is scaled at will, all the others follow, providing the measurement gains are known. More generally, impulse forcing @(t) of compartments j~Ji and step forcing 4 S&-s(t) dt of compartments j EJ~ gives a response in compartment i
100 The constant normal-mode
J. P. NORTON
term in xi gives no information, coefficients E,q
as we have by summing
the
(‘5) On the right-hand side, the first term is known and the second is identical to the constant term. Redundancy occurs in thep equations for Eiq either if all the state-matrix elements uij, ~EJ,, (perhaps none) are known in advance-for, using (7) and (lo),
-or if there are no step inputs, so that the second term on the right-hand side of (15) vanishes. It is important to note that, as a consequence of these possible redundancies, the identifiability of a particular model may be altered by the choice of input signals. This fact is obscured by the assumption, in the usual definition of identifiability, that the transfer-function or impulse-response matrix is available. Together (8) (9) (lo), and (14) specify all the information available after an identification experiment. The test for unique identifiability involves first noting the rows of M and columns of N found directly from the normalmode coefficients, then counting the equations linear in the remaining rows and columns, deieting any not usable in successive linear-equation solutions for the unknown rows and columns because of redundancy as explained later. If there remain at least as many such linear and independent equations as unknowns, the model is uniquely identifiable. If not, any redundant set of the equations gives rise to a further equation obtained by setting the determinant of their coefficients to zero. Such an equation may be linear in an unknown row or column and allow the solution to proceed. When a solution in successive linear stages is not possible, it is usually straightforward to see what change to the model or experiment would help. The following two-compartment example illustrates the basic ideas. Example 1. Compartment 1 of a two-compartment system is perturbed by a unit-impulse input, and its response is observed. The transfer rate constant k,, from compartment 1 to the environment is known to be zero. The measured response of compartment 1 is
c,x,(t) where ci is the measurement
= c,r~e”‘n,,
gain. The observed
time constants
are xi and
IDENTIFIABILITY
ANALYSIS
IN LINEAR
101
STAGES
h,, the diagonal elements of A. Because the scaling of each eigenvector column of M is arbitrary, rT can be taken for convenience as [ 1 l] provided it contains no zeros. This is so, as both the normal modes appear in the observed state variable. If c, is known, n, follows from the measured amplitudes of the exponential normal-mode components of x,(t). The equations in the remaining unknowns r2 and nz are rTn2 = 1
r:n, = 0,
rrn2 = 0,
from (10); and from prior knowledge
of k,,, using (9), we have
(rl + r,)=An,
= 0.
With r, and n, known already, r, can be found uniquely r2‘[ nl
An,]=[O
from
-rTAn,].
At this stage, M is completely known, and the analysis has shown that the model is globally identifiable. If numerical values of the rate constants were being computed, n, would follow from
[I [I 6 6
and the unknown
rate constants
O
n2=
1 ’
from
k,2 = rrhn,, k,, = r,Thn,,
k,,=
-(r,
+r,)*An,.
Comparisons between methods of analysis on such a small example are perhaps pointiess, but it is worth noting that use of the transfer function requires algebraic formation and then partial-fraction expansion of C[sZA]-‘B, looking out for pole-zero cancellation while doing so, and then examining the resulting equations to see if they are soluble for the rate constants. The Markov-parameter-matrix approach requires formation and differentiation with respect to each unknown of CA ‘B for values of i up to 2p - 1, and then examination of the Jacobian matrix to see the conditions for full rank. The modal solution above is not complicated compared with these. If c, had been unknown in Example 1, the proportions of the elements of II, would have been given by the ratio of the normal-mode amplitudes, and
102
J. P. NORTON
the scaling by r:n, = 1 from (10). This latter equation is redundant when c, is known beforehand. Since the effects of unknown measurement gains are straightforward to assess, they will be taken for brevity as known and unity in what follows. Had k,
rather than k,, been zero, the prior-information
equation would
have been (r,
+r2)TAn2=0,
and instead of solving linear equations
immediately
to hand for r2,
would have been used to give
I
6
I
=o.
(rl +r21TA
since n, is not entirely zero. This is a linear equation in r2, of the form
f Tr2= g, where f and g are known, depending
only on r, and A, and g is nonzero.
Hence r2 is found from
[ 32=[
3.
The rest of the solution is as before. Two types of special case exist. The first is when not all normal modes appear in the observed state variables. It causes no great difficulty, as seen in Sec. 4. The second, and the harder to analyze, is when definite but nonunique solutions are obtainable. One such example is given by Skinner et al. [14] and several by Jones and Godfrey [15], using the nonlinear equations in the rate constants given by the coefficients in the rational polynomial Laplace transforms of the observed state variables. The mechanism by which nonunique solutions arise is not made clear by the Laplace-transform analysis. It is discussed in terms of eigenvectors and eigenvalues
in Sec. 7.
IDENTIFIABILITY
4.
ANALYSIS
OBSERVABILITY
IN LINEAR
STAGES
103
AND CONTROLLABILITY
The object of this section is just to examine the implications of incomplete observability and controllability on the workability of the normalmode analysis, not to give a general review of their relevance to identifiability. Pohjanpalo [16] gives a simple check for the excitation of a given mode in a particular state variable. Writing the matrix exponentials in (8) as h4e”‘N and Meh(‘-7)N, it is seen that if some mi9 is zero, normal mode q does not appear in xi(t), and if some nsi is zero, neither direct forcing of compartmentj nor an initial value for xi excites normal mode q; the system is partially unobservable from x,(t) and uncontrollable from inputj respectively. The occurrence of zeros in M or N poses no problem in testing identifiability by the normal-mode method. Failure of normal mode q to appear in xi as a result of forcing j does not indicate on its own which of miq and nti is zero. However, unless another observed state variable does include that mode, X, and consequently the system as a whole are unidentifiable. If the mode appears in any other observed state x,,, nd is not zero, so mi4 is zero and either mk or nti may be chosen freely to scale eigenvector q. Zeros in rows of M or columns of N that are not seen directly in the normal-mode coefficients but found later on in the calculations are as tractable as nonzero values in principle, and will generally occur as small nonzero values in practice, because of experimental and modeling error. With distinct eigenvalues, failure to excite mode q is indicated by row q of NB being zero, and failure of mode q to appear in the observed outputs implies that column q of CM is zero. These conditions are the modal corollaries of the rank criteria of Kalman for complete controllability and observability, to which a number of authors refer [13, 171. 5.
REDUNDANCY
DUE TO ZEROS IN THE STATE MATRIX
Returning to the general case, any extra equations needed to make the model identifiable are given by prior knowledge of rate constants, particularly knowledge that some are zero. Zero rate constants corresponding to zero entries in the state matrix are also the principal cause of redundancy in the bilinear equations to be solved for the rows of M and colunms of N. If a,* is known a priori,
(17) can be used with other equations either knowing n* to find r, or knowing rr to find n,. Taking the first possibility and noting that An, = col. s of AN =col. s of NA = 2 q-1
a,,n,,
(18)
104
J. P. NORTON
the coefficient of the unknown r, in (17) is seen to be a linear combination of columns of N. Apart from (17) there are p equations
Cn, = arq,
q=1,2
>...,p,
(19)
in r,, and perhaps others like (17) but with s different. To avoid redundancy in the set of p linear equations to be solved for r,, any one of (19) replaced by (17) must have q such that uqs is nonzero for at least one of the equations (17) included in the set. If not, none of the coefficients of r, is a function of nq, so the p equations are linearly dependent, being functions of at most p - 1 columns of N. Among other things this implies that u,~= 0 cannot be used in place of r;n, = 1. Similarly, in finding n, any equation from
rqTn= s,,, replaced equation
s=1,2 >...,p>
by (17) must have q such that a,, is nonzero (17) in the set to be solved for n,.
(20) for some such
Example 2. A 3-compartment model has k,, known, so aI3 is known. The responses of x, and xa to an impulse input into compartment 1 are observed, with known measurement gains. On choosing the scaling of rI or r3, n, is found immediately from the normal-mode coefficients mlqnql or msqnql in x, or x,; then r3 or r, follows from the normal-mode coefficients in the other observed state. It is proposed next to find n3 from
rT 6
I [I n3=
0 1
a13
:
rrh
However, rrR=a,,rr+
a&+
a,,$,
so if a,* is zero, the third equation in n3 is implied by the other two. Here the redundancy is not affected by whether a,3 is zero or not, but if x2 rather than x3 had been observed, giving r, instead of r3, a zero value for aI3 would have made the set
redundant, whereas a nonzero been immaterial.
value would not. In that case, a,2 would have
IDENTIFIABILITY
Knowledge
ANALYSIS
IN LINEAR
of a rate constant
to the environment,
5 r:An,
105
STAGES
k,, gives
= - k,.
(2’)
q=l
This is potentially useful to help find any one of the rows of M. Redundancy arises just as before except that rTns = 1 can usually be replaced, since a,, is nonzero unless compartment s is a pure sink. The sole complication in checking the usefulness of (21) is that the check extends over all q. The other potential use for (21) to help find 9, is significant only if
i q=l
r,’
is known before all the rows of M are known individually, since otherwise n, can be found using (10) without recourse to (21). For this to be so, it must be possible to replace at least one of the equations E rFns= 1,
s= 1,2 )...) p,
(22)
q=l
[which follows from (lo)] with an equation (21) not dundancy. The analysis proceeds as before, but with
replacing any one of the rows of M, and the corresponding replaced by (22). 6.
TESTS FOR UNIQUE
introducing
re-
equations
(10)
IDENTIFIABILITY
All the experimental and prior information concerning the transfer rate constants appears as bilinear equations in rows of M and columns of N. They are linear in any one row or column, and if required more such equations may be obtainable, as in Example 1 with k,, zero, from the determinants of any consistent but redundant subsets of the bilinear equations. A unique solution is guaranteed if the equations can be solved in stages, each involving only collectingp equations in which all but one row or column have already been found. The selection of such a subset can be viewed as replacing any of the p equations (10) in the current unknown which also contains another unknown row or column by one from (17) or (21) which does not. The rows and columns given directly by the choice of eigenvector scaling are known at the start of this process.
106
J. P. NORTON
To test the adequacy of the information available in a given experiment, each equation (17) or (21) giving a transfer rate constant a priori can be checked for redundancy in each of its potential uses as a replacement for an equation (10). By this means one experiment is tested, but it is worth noting that for a realistic number of compartments (up to about three if the exponential components of each observed response are to be resolved with acceptable accuracy), it is quite feasible to set out exhaustively all completely linear solutions using every adequate collection of prior information for a given combination of forcing and observations. In other words, all the uniquely identifiable model structures compatible with the experimental conditions can be found. The factors which may disqualify from use an equation contributed by prior knowledge of a rate constant will now be listed. Consider first (17): r:An,= a,* = k,s.
(‘7)
r:n, = S,,
(‘9)
It may not be used to replace
in forming
a set of equations
for rr if
(i) x, is observed, so that r, is found immediately from the coefficients of the normal modes in x, or chosen to fix the scaling of the eigenvectors, and (17) is not needed; (ii) the only forcing input is into compartment q, so that n4 is found immediately from the normal-mode coefficients of an observed state variable or chosen to fix the eigenvector scaling (only if n4 were unknown would there be any point in removing r:n, = S,, from the set to be solved for r,; this disqualification does not apply if there are two or more forcing inputs); (iii) q = s, so that if n, is known and use of (17) is possible, nq is also known and it is not necessary to replace (19); (iv) Us=is zero, so that the replacement would introduce redundancy as explained in Sec. 5; or (v) q= t and u,~ is zero, so that redundancy would arise as explained in Sec. 5. In the same way, (17) may not replace rTns = S,,
(20)
IDENTIFIABILITY
ANALYSIS
IN LINEAR
STAGES
107
in finding n, if (i) the sole forcing is of compartment (ii) x4 is observed; (iii) q= t; (iv) a,, is zero; or (v) q = s and u,~ is zero. The small modifications when k, rather earlier. The possible redundancy explained in all a0 are known, with xi observed and Commonly there are no step inputs, so and the sum
s;
than k, is known have been noted Sec. 3 may be checked by noting if impulse forcing of compartment j. ri for each observed state variable
over the impulse inputs are found directly from the observed normal-mode coefficients. Consequently, if we rewrite (10) by replacing any one nj,j EJ,, with
then the rewritten
equations
(23) are redundant and are deleted. It is simplest to carry out all the checking for redundancy in terms of n; in such cases. The following examples show how global identifiability is tested by examining the feasibility of a unique solution for the modal matrix and its inverse. Example 3. A model used in pharmacokinetic experiments on the hepatobiliary system has 3 compartments: the plasma, liver, and biliary processes. An impulse input (rapid intravenous injection) forces compartment 1 directly; xi and xg are observed; and k,,, k,,, k,,, and k,, are taken as zero, defining the structure of the model. The forcing amplitude is 1, and the observation gains c, and cs are known. Choose r, to be a vector of ones; then n, follows from the coefficients of the exponentials in x,, and hence r, from those in xj. Check each known-
108 rate-constant
J. P. NORTON
equation
for usefulness:
(i) k,, gives (ri + r, + r,)Thn, = 0, which cannot be used as r;An, = known quantities to accompany r;n, = 0 and rln,= 1, because a3, is zero and so An, is linearly dependent on n, and n2. Its only potential use is to accompany rln, = 0 and r:n, = 0 in finding r2. (ii) k,, =0 gives r;An, =O, which is useful only to replace rln,=O, i.e. to accompany rFn3 = 0 and rTn3 = 1 in finding n3. (iii) k,,=0 gives r:An, = 0, which cannot replace rcn, =0 or rcn, = 1 in finding r,, because doing so would give a redundant set of equations, and cannot replace an equation with unknown coefficients of n3 in finding n3, since both rl and r3 are already known, and replacing r:n3 = 0 is pointless. This item of prior information is therefore useless. (iv) k,,=0 gives rTAn, =O, completely redundant, as r3, A, and n, are already known. Seven equations (10) remain once the redundant rTn, = 1 and r:n, = 0 are removed, so with the two useful items of information k,, =0 and k,, =0, there are 9 equations in the 9 unknowns in rz, n2, and n3. Solution by linear stages gives n3 from
rT
I I [I rTh
0
0 ,
n,=
1
rT
then r, from
rT[4
An,]=[O
n3
0
-(r1+r3)TAn,].
By now M is known uniquely, so the model is uniquely identifiable. What is more, it is clear which items of prior information are essential and which are redundant. It is interesting to compare the analysis above with the transfer-function and Markov-parameter-matrix investigation of the same example. The transfer-function approach gives
Y(s)=
[ :,
=
;
1(s+k,,)(s +k,, +k,, +4112 element
[
=-[
q,,,,1-I[i] (1,l)
of [sl-A]-’
element (3,1) of [sl-A]-’
A
b,kz
IDENTIFIABILITY
ANALYSIS
IN LINEAR
= (s + k,,)(s - a)(s - P),
109
STAGES
say
The observed eigenvalues are - k,,, o, and fl. Here the attribution of -kc,, is unambiguous because it appears in x3(t) but not, because of pole-zero cancellation, in x,(t). The amplitudes of the observed exponential components are the partial-fraction coefficients of X,(S) and X,(s), i.e. the residues (Y, and /3. For X,(s) they are at - ka,, e II=
o+ke2+k,Z+k32 o--P
and
e _ P+koz+ki,+k,, 12-
P-a
’
and for X,(S) they are
Examination of the equations, with (Y, p, and k,, known, shows that the remaining rate constants can be found uniquely in the order k,,, k3*, k,, in the and k,,. It also becomes apparent that there are two redundancies coefficients e3,, e32, and e33, and one in e,, and e,2. In spite of the simplification due to - k,, being an eigenvalue, the algebra is not negligible, and some features do not emerge clearly. For instance, prior knowledge of k,, and k,, is shown by the modal analysis to be unnecessary, but a reworking (with very complicated algebra, as the reader is invited to confirmj is needed to verify this by the transfer-function method, and, more to the point, the original transfer-function analysis does not bring these redundancies to one’s attention. The Markov-parameter-matrix method requires the formation of the derivatives of terms up to CA’B with respect to the rate constants. Merely writing out CA3B for this example will convince the reader that the analysis is lengthy and its structure is obscured by algebraic detail. The next example is intended to show how a modal solution by linear stages may require new linear equations to be derived from singular subsets of the original equations in rows of M and columns of its inverse. Example 4. A 3-compartment model has k,, = k,, = 0, so that the peripheral compartments 2 and 3 interact with the central compartment 1 but not with each other. Compartment 1 does not excrete to the environment, so k,, = 0. A unit impulse forces compartment 1 and xi and x2 are observed. Eigenvector scaling and the amplitudes of the exponentials in xi and x2 determine r,, r,, and n,. There are 7 nonredundant equations (IO), 3
110
J. P. NORTON
prior-rate-information equations, and 9 unknowns in r3, n2, and n3. However, a linear solution is not immediately possible for any of these, as the only substitution which would produce a set of 3 equations in one unknown, r;An, = 0 for rrn,= 1, gives a redundant set
6 6 :
r;h
I [I 0 0 . 0
n3=
Instead, all redundant sets of 3 equations linear in the same row or column are examined. They consist entirely of equations with right-hand side zero, and give rise to
In1
n3
An,\=0
from equations
linear in rz,
(244
In1
n2
An,l=O
from equations
linear in r3,
(23
=o
from equations
linear in n,,
PC)
=o
from equations
linear in n,,
PD)
=o
from equations
linear in n3,
(29
6 6 (r,+r,+r,)=A
rT rT r:A rT rT r;A
Substitution of rf’n2=0 and r$n,= 1 into Eq. (24B) gives a quadratic in any element of 9, so there are two solutions for n,. One of them is easily seen to be a scaled version of n3, which appears in Eq. (24A) just as n, appears in Eq. (24B), and satisfies rrn3 =0, but has rTn,=O. Consequently there is in fact only one solution for n,. The nonredundant set r3‘[ nl
n2
An,]=[
0
0
-(r1+r2)TA]
then gives r,, so A4 is complete and the model is globally identifiable. that r3 could not have been found from the set
6b which is redundant.
n2
An2]=[o
o 01,
Notice
IDENTIFIABILITY
ANALYSIS
IN LINEAR
STAGES
111
A similar and quicker alternative is to use Eq. (24D) with r:n, = 0 and r;An, = - (r, + r,)TAn, to find r,; again two solutions arise, one of them being a scaled version of r, because Eq. (24C) is identical in form to Eq. (24D) and r:n, = 0 to rrn, = 0. With r3 known, M is complete. 7.
LOCAL IDENTIFIABILITY
When a particular model has been found not to be uniquely identifiable from a given experiment, it may still be locally identifiable, that is, there may exist a number of distinct solutions for some or all of the rate constants. If so, background knowledge or other experimental work may suggest which is the correct one. Nonunique solutions arise in two ways, through singular sets of equations or through incomplete coupling in the state matrix. The first of these is best understood by reference to an example. Example 5. A 3-compartment model has k,,, and kz3 zero. A unit impulse is injected into compartment 1, and the responses of compartments 1 and 3 are observed with known measurement gains. From these observations A, r,, r,, and n, can be found immediately, leaving 9 unknowns making up r,, n2, and n3. The usual checks show that knowledge of ko, allows substitution of r,TAn, = -(r, +r,)TAn, for r,Tn,= 1 or rTn3=0, depending hand, the equation
(25)
on the order of solution.
On the other
r:An, = k,, = 0 cannot be used immediately, as r, and II, are both unknown, and it is redundant when r, is unknown and A and n3 known, being implied by rTn, = 0 and rTn3= 0 from (lo), because An, = a,+, + a,+, + a33n3= a+,
+ a3,n,.
The three equations in r, could therefore not be solved even if n3 had previously been found, but the linear relation between An,, n,, and n3 gives an equation in ng:
b
n3
An,]=O.
With A and n, known, this is used together with rrn3 = 0 and rrn3 = 1 to give a quadratic in each element of n,. Either of the two resulting solutions for n, can be used in rTn3= 0 with (25) and r;n, = 0 to give r,; then n, follows from (10). As there are unique solutions for r,, r,, n,, and A, but two sets of
112
J. P. NORTON
consistent solutions for r2, n2, and n3, the equations (8) give k3, uniquely but give two solutions for each of the other rate constants. An alternative order of solution uses k,, = 0 and (lo), or equivalently
to obtain
which with (25) and r:n, =0 gives two solutions from (10).
for r,; then n, and n3 follow
A singular set of equations in one unknown row of M or column of N may also be obtained from a set of p + 1 consistent equations. If they are all in the same row and in only two columns, one known and the other unknown, or ali in the same column and in only two rows, one of them unknown, then a singular set of equations can be produced by scaling and subtraction. Example 6. In a 3-compartment model, k,, and k,, are both zero, and as the input is injected into compartment 1, n, is found immediately. The two zero rate constants give
and
and two more equations namely
in r, + r, + r,, nl, and n3 can be obtained
(r, +r,+r,)*n,=
1
(r,+rz+r3)Tn3=
1.
and
Together
these four equations
imply that
(rl+b+rdT[An,
An3 n,-n,]=[o
o 01,
from (lo),
IDENTIFIABILITY
ANALYSIS
IN LINEAR
STAGES
An,
n, -n3 1=O.
113
so that (An,
The other mechanism which may lead to a finite number of distinct solutions is the existence of a normal mode which appears in only one state variable, or a state variable which shows only one mode. Mode j appears only in xi(t) if and only if columnj of the state matrix A is zero but for uii, i.e., compartmentj loses material only to the environment. Correspondingly only element j is nonzero in eigenvectorj and I+ and (4) with either (6) or (18) gives A, = aii = - k,.
(26)
The rate constant identified thus is nonunique in cases where without observing any other state variable one cannot decide which of the observed eigenvalues is 4. Similarly if row j of A has only ui, nonzero, i.e., compartment j receives material only from outside the system, then rj and rowj of N have only elementj nonzero, and (26) is again true. Here also 4. may not be identifiable uniquely, for unless 3(t) is observed, one may not know which normal mode corresponds to h, and if xj(t) alone is observed, the other modes are not seen. system with k,, zero has compartment Example 7. A two-compartment 1 forced by an impulse and compartment 2 observed. If the forcing amplitude bi and observation gain cz are both known, free eigenvector scaling and the amplitudes of the observed exponentials give r, and n,, and the eigenvalues, both appearing in x2, give A. Because A has only one nonzero element in row 1 (and column 2), r, has only its first (and n, its second) element nonzero. Accordingly rrnl= 1 gives r,,= l/n,,. However, this value depends on which mode was labeled mode 1 in calculating n, 1, so there are two distinct solutions for r, (and nJ, and the model is only locally identifiable from this experiment. If also it is known that k,, is zero, giving r;An, = - r,TAn,
(27)
the positions of the zeros in A, M, and N are not altered, and r, is calculated as above. Now, however, only one solution for r,, is consistent with (27), which may be written r 11=
r22A2n2, -rzl-
7’ 1
11
J. P. NORTON
114 so this model is uniquely identifiable from the experiment, of the mode attribution being resolvable.
8.
as a consequence
CONCLUSIONS
By the use of normal modes and introducing the rows of the modal matrix and columns of its inverse as unknowns, both the response of a compartmental system to impulse and/or step forcing and any prior information on parameter values have been written as equations which are bilinear in the unknowns. In consequence it is relatively straightforward to check for redundancy and the feasibility of a stage-by-stage linear solution, and hence to determine whether the model is uniquely identifiable. The alternatives, use of Laplace-transform or time-domain expressions for the responses written directly in terms of the parameters (rate constants) to be identified, yield in general nonlinear equations ill suited to checking for redundancy or multiple solutions. The use of normal modes also makes it relatively easy to see what modifications, if any, are needed to bring about identifiability. The origins of nonuniqueness of identification are clarified by normalmode analysis. Because the procedure has been reduced to an algebraically simple and systematic one, computer checking of a wide variety of experiment model combinations for global and local identifiability may be useful. This work was carried out during a Visiting Fellowship at the University of Warwick. Comments by Mr. R. P. Jones and Dr. K. R. Godfrey of the University of Warwick and particularly Dr. R. F. Brown of the University of New South Wales are gratefully acknowledged. The comments of a reviewer have been most helpful in improving the presentation. REFERENCES R. Bellman and K. J. Astrom, On structural identifiability, Math. Biosci. 7:329-339 (1970). J. A. Jacquez, Comparfmental Analysis in Biology and Medicine, Elsevier, Amsterdam, 1972. J. J. Di Stefano III and F. Mot-i, Parameter identifiability and experiment design: thyroid hormone metabolism parameters, Amer. J. Physiol. 233:R134-R144 (1977). K. Glover and J. C. Willems, Parametrizations of linear dynamical systems: canonical forms and identifiability, IEEE Tram. Automatic Control 19:64(X645 (1974). E. Walter, G. Le Cardinal, and P. Bertrand, On the identifiability of linear state systems, Math. Biosci. 31:131-141 (1976). M. S. Grewal and K. Glover, Identifiability of linear and nonlinear dynamical systems, IEEE Trans. Automafic Control 21:833-837 (1976). C. Cobelli, A. Lepschy, and G. Romanin-Jacur, Identifiability of compartmental systems and related structural properties, Math. Biosci. 44: I- 18 (1979).
IDENTIFIABILITY
8 C. Cobelli identifiability 9 10
ANALYSIS and
G.
Romanin-Jacur,
of multi-input
12
the formulation of models, J. Delforge, The problem
15
16 17
Controllability, multi-output
Biosci. 11:203-247 (1971). E. Walter, P. Bertrand, and G. Le Cardinal, IEEE and nonlinear dynamical systems,” (1979). M. Berman
14
and
and R. Schoenfeld,
115
STAGES
IEEE Trans. Biomed Eng. 23:93-100 (1976). S. I. Rubinow and A. Winzer, Compartment
11
13
IN LINEAR
Invariants
observability
biological analysis:
and
structural
compartmental an inverse
problem,
systems, Math.
Comments on “Identifiability of linear Trans. Auromatic Control 24~329-331
in experimental
J. A&. Phys. 27: 1361-1370. of structural identifiability
data on linear kinetics of a linear
and
compartmental
system: solved or not?, Math. Biosci. 36: 119-125 (1977). W. D. Smith and R. R. Mohler, Necessary and sufficient conditions in the tracer determination of compartmental system order, J. Theoref. Biol. 57: l-21 (1976). S. M. Skinner, R. E. Clark, N. Baker, and R. A. Shipley, Complete solution of the 3-compartment model in steady state after single injection of radioactive tracer, Amer. J. PhysioI. 1961238-244 (1959). R. P. Jones and K. R. Godfrey, Mathematical of linear compartmental models, University Report 74, 1978. H. Pohjanpalo and B. Wablstrom,
analysis of the structural identifiability of Warwick Control Theory Centre
The uniqueness
of linear compartmental
Inrernat. J. Sysfewu Sci. 8:619-632 (1977). M. Milanese and G. P. Molino, Structural identifiability and pathophysiological information from the kinetics 26:175-190
(1975).
systems,
of compartmental models of drugs, Math. Biosci.