IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-23, NO. 2, MARCH 1976
93
Controllability, Observability and Structural Identifiability of Multi Input and Multi Output Biological Compartmental Systems CLAUDIO COBELLI,
MEMBER, IEEE, AND
Abstract-Controllability, observability and structural identifiability of biological compartmental systems of any structure are studied in order to evaluate the significance of a multi input-multi output experiment, that is, the a priori possibility of estimating the unknown system parameters through the chosen experiment. Some basic concepts and properties of biological compartmental systems and corresponding models are rirst briefly reviewed to put compartmental models in a linear system theory context. Controllability, observability and structural identifiability are defined for compartmental systems. New theorems on compartmental systems relating controllability and observability to their structure are proved. Theorems on transfer function matrix are suitably applied to solve the problem of structural
identifiability. Solution techniques to check controliability, observability and structural identifiability of compartmental systems of any structure are given. The flow chart for a digital computer implementation of these techniques is presented. Some examples concerning tracer analysis experiments on biological compartmental systems concerning digoxin, iodine, albumin and thyroxine distribution, and calcium metabolism illustrate the application of the mathematical theory developed in the paper in order to emphasize the importance of the problem of structural identifiability in biology and medicine.
1. INTRODUCTION
IN tracer analysis of biological systems, compartmental models have been widely adopted in research, diagnosis and therapy problems. These models approximate the quantitative description of phenomena, thus allowing many parameters of biological interest to be estimated. The most widely treated problems in biomedical literature regarding compartmental models concern compartmentalization of biological systems and parameter estimation, that is, the development of models using a compartment approach and the numerical estimation of the associated parameters from
input-output tracer experiments. Actually the numerical estimation of parameters should follow a test on the a priori possibility of their estimation by the chosen input-output experiment, so that the significance of the experiment itself is tested. This test will be referred to as the solution of the problem of structural identifiability, as it is related only to the model structure; the solution determines which and/or how many parameters of the adopted model can be estimated by means of the chosen input-output experiment. Manuscript received August 2, 1974; revised April 17, 1975. C. Cobelli is with the Laboratorio per Ricerche di Dinanica dei Sistemi e di Elettronica Biomedica, Consiglio Nazionale delle Ricerche, Casella Postale 1075, 35100 Padova, Italy, and the Istituto di Elettrotecnica e d'Elettronica, Universiti di Padova, Italy. G. Romanin-Jacur is with the Laboratorio per Ricerche di Dinamica dei Sistemi e di Elettronica Biomedica, Consiglio Nazionale delle Ricerche, Padova, Italy.
GIORGIO ROMANIN-JACUR
This problem has been considered only recently [1], [2], [3]. Bellman et al. in 1970 [1] put the problem in a system theory context and stated precise general definitions. Hayek in 1972 [2] provided solution techniques for linear invariant compartmental systems open to the environment; solution is given only for the one input-one output case (not directly extendible to the multi input-multi output case) under complete controllability and complete observability conditions. In [3], the authors study structural identifiability for multi inputmulti output strongly connected compartmental systems both open and closed. Theorems relating controllability and observability to the system structure are proved and solution techniques for structural identifiability are developed; a flow chart for a digital computer implementation is presented and several application examples on biological systems are reported. The present paper studies structural identifiability for multi input-multi output compartmental linear invariant models of any structure, that is, the limiting restriction to the class of strongly connected systems [3] is removed. Briefly, in section 2 some general notes about compartmental models and their structural properties (particularly controllability, observability and structural identifiability) are given. In section 3 new theorems relating controllability and observability of compartmental systems to their structure are presented (proofs are given in the appendix). The general solution of the structural identifiability problem is discussed in section 4. In section 5 some testing criteria and the flow chart for the implementation on digital computer of the problem considered in 3 and 4 are given. Finally, in section 6 several compartmental systems of biomedical interest are analyzed from the point of view of their structural identifiability. 2. BASIC NOTES ON BIOLOGICAL COMPARTMENTAL SYSTEMS
2.1 CompartmentalSystems andModels Some basic concepts, useful for the following, will be briefly reviewed; a more comprehensive treatment can be found for instance in [4], [5]. In a biological compartmental system small deviations from its steady state are induced from the outside; under the assumptions of volume constancy for each compartment, linearity and invariance, the dynamical behavior of each compartment is described by (Fig. 1):
Xi= pki £ kixj+ ui- j*iE kiixi- koixi
O e wher where:
(1)
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, MARCH 1976
94
ai = kij,
Ui
i*j
aii = -koi-
n j=1
(4)
ki.
(5)
The eigenvalues of matrix A cannot be purely imaginary and
their real part is non positive (it can be zero when the eigenvalue is real) [6]. Matrix B is of n X rb dimension, with rb < n; if inputs are in compartments 11, i2, ib its entries can assume either value 1 or 0 under the following condition:
xoi Xi Fig. 1. Input
output material fluxes of the generic i-th compartment.
and
E bij= I Vj;
bi
(6)
,
Matrix C is of r, X n dimension with r, S n; if outputs are * h, its entries can assume taken from compartments il, either value 1 or value 0 under the following condition: 12,
E
j=l
F0 x
XO 07
xl0
ox 0
x
o x
o
LO
01
o
Fig. 2. K and Ko matrices of a compartmental system.
is the state variable associated with the i-th compartment, that is, the amount (or concentration) of material of the i-th compartment; equation (1) is also valid when xi is the amount (or specific activity) of the labeled material (in this case the system may be non linear); kii is the rate constant from the j-th compartment to the i-th compartment; koi is the rate constant from the i-th compartment to the environment; ui is the external input.
xi
The rate constants {kij} are grouped into a square matrix K of n rows, where n is the compartment number, and into a row matrix Ko of n entries. Notice that all entries of matrix K are non negative and those of the main diagonal are zero; all entries of Ko are non negative (Fig. 2). The mathematical model of a compartmental system corresponds to the usual input-state-output representation of a dynamical linear invariant system: x
=Ax + Bu
y =
Cx
(2) (3)
where x, u, y are respectively state, input and output vectors. With reference to the structure of the compartrnental system matrices A, B, C show peculiar characteristics. The order n of matrix A equals the number of compartments; its entries are related to those of the matrices K and Ko
by:
Ci1 =1 Vi;
E C.j =
0l,
pj] j1
,
.*'
C
(7)
Compartments are represented graphically either by blocks or by nodes and transferences by directed segments called branches. A path is a succession of branches such that the terminal block (node) of a branch is also the initial block (node) of the next branch. A compartmental system can be divided into subsystems made up of one or more compartments, such that every compartment in a subsystem is directly connected with at least one compartment in the same subsystem. Any subsystem S is characterized by the corresponding matrix A,. A compartmental system is said to be disjoint when it contains one or more isolated subsystems not interacting with the remaining part of the system. In the following, only non disjoint systems will be considered. A subsystem is said open (closed) when there is some exchange (no exchange) of material from the subsystem either to the environment or to another subsystem. A subsystem is strongly connected when every compartment can be reached from every other compartment in the same subsystem along at least one path. If matrix KS of the rate constants of subsystem S is considered, it can be noticed that the general entry of KS (Ks at r-th power) is given by: (8) kil I ki 12 k Ir-J5 [Kr, ij = E .
11
a2,
r-1 C
-
-
1
and then [Ks] ij > 0 implies that at least one path of length r exists from compartment i to compartment i within S. If matrix Rim) = I7. K' is considered, it can be noticed that [R(m)] iq > 0 implies that at least one path of length smaller than m, or equal to m, exists from compartment j to compartment i. If m = ns (where n. is the number of compartments in subsystem S) matrix R('s) gives the following information: a) positive entries of the j-th column denote the compartments reachable from compartment j; b) positive entries of the i-th row denote the compartments from which comnpartment i can be reached. From the previous remarks, it follows that matrix R ns) of a strongly connected subsystem has all entries positive.
95
COBELLI AND ROMANIN-JACUR: BIOLOGICAL COMPARTMENTAL SYSTEMS
In all previous definitions and remarks subsystem S may coincide with the whole compartmental system.
2.2 Controllability, Observability and Identifiability in a Structural Sense of Compartmental Systems Rigorous treatments on structural properties can be found in system theory literature (see for instance [7], [8], [9]). In the following, only tests on structural properties will be considered. In the theory of linear invariant systems as the compartmental ones considered, a system is said to be completely controllable (CC) if and only if the controllability matrix of n X (n - rb + 1) * rb dimension: (9) P = [BIABIA 2BI ... {An rbB] is of rank n, and is completely observable (CO) if and only if the observability matrix of n X (n - r, + 1) r, dimension: (10) Q = [CT IATCTI A2 TCTI * *An-rcT CT] is of rank n, where n is the order of the system (compartments number), rb and rc are the ranks of matrices B and C (inputs and outputs number, respectively). Conditions (9) and (10) correspond to the following ones: -
3. THEOREMS ON CONTROLLABILITY AND OBSERVABILITY IN A STRUCTURAL SENSE OF COMPARTMENTAL SYSTEMS In this section three new theorems relating controllability and observability of compartmental systems to their structure are presented. Proofs are reported in the Appendix. Theorem 1: The existence of at least one path reaching every uncontrolled compartment from a controlled one is a necessary condition for a compartmental system to be CC; the existence of at least one path from every unobserved compartment to an observed one is a necessary condition for a compartmental system to be CO. Note that if the system is strongly connected the necessary conditions of theorem 1 are always satisfied for any input and output; in fact, by definition, every compartment can be reached from every other one along at least one path. Theorem 2.1. A compartmental system is CC in a structural sense only if there is at least one path from a controlled compartment to every uncontrolled one. Theorem 2.2. A compartmental system is CO in a structural sense if there is at least one path from every unobservable compartment to an observed one. Theorems 2.1 and 2.2 always hold for strongly connected compartmental systems (see definition in 2.1) [3 ]. Theorems 2.1 and 2.2 allow us to state that the rank of matrices P and Q can be less than n only in a subspace of the { kie > 0} space. In fact, in the case of one input-one output systems, each of the equations det [ppT] = 0 and det [QQT] = 0 is the general expression of an hypersurface, therefore the non CC and non CO subspaces are represented by hypersurfaces. For -the multi input-multi output case, there are as many hypersurface equations as the submatrices of rank n of matrix [ppT] and [QQT], respectively. Non CC and non CO subspaces are then represented by the intersections of the hypersurfaces described by the equations relative to the considered submatrices of [PPT] and [QQT] , respectively.
(11) det[PPT] 0 (12) det [QQTI 0O. A compartmental system can be non CC and/or non CO either because of its own structure or because of a particular combination of the parameter values {k,j}. In the first case the rank of matrix P and/or Q is always less than n for whatever values of rate constants relative to the system of fixed structure (unchanged topology); it means there are no possible connections between a certain input and a certain state or between a certain state and a certain output. A compartmental system is said to be non CC (non CO) in a structural sense when it is non CC (non CO) for any values 4. STRUCTURAL IDENTIFIABILrrY OF combination of the non zero parameters {k,1}. COMPARTMENTAL SYSTEMS Input-output relations obtained from (2) and (3) are considered in order to define structural identifiability. As seen in sect. 2.2, a controllable and observable dynamical The input-output behavior of a compartmental system is linear invariant system is completely described by the transfer completely described, relatively to its CC and CO part, by the function matrix G(s) [see (13)]. The general entry [G(s)] i transfer function matrix G(s): where i = 1, r, and j = 1, rb is the transfer function between input j and output i. To find out the number of coefficients I A C * adj(sI- A) * B. G(s) = C (sI - A)-1 * B = det (sI- A) that can be estimated, the expansions of det (sI -A) and C adj(sI - A) * B are considered ([9], [10]; see also [3] for (13) details) and G(s) is given by: The knowledge of G(s) allows the estimation of a maximum number of parameters of the CC and CO part of the system [CB(sn1 +a1 sn -2 + - + an1) G(s) = (5 equal to the number of independent coefficients of the numerator and denominator polynomials of the entries of G(s). + CAB (sn 2 + a,X Sn -+ +n -2a ) A CC and CO compartmental system is said to be structur(14) ally identifiable when the chosen input-output experiment + +CAn-2B(s+a,)+CAn -B]. allows the estimation of all unknown parameters (rate [G(s)] i, has the same expression as G(s) if we take respectively constants). If the whole compartmental system is not both CC and CO, [CB] i, [CAB] ii, *- [C4n -1 B] ji. The number of parameters that can be estimated are related complete controllability and observability may hold for a part of the system; in this case the considerations on structural to the polynomials in s of the numerator and denominator of identifiability are restricted to the CC and CO part of the every [G(s)],,. Cancellations occur whenever the subsystem CC through i and CO through i is not equal to the whole comsystem. -
-
-
96
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, MARCH 1976
partmental system: more precisely, the number of simplifiable prime factors equals the number of the compartments uncontrollable through i and/or unobservable through i. Therefore the denominator associated with every [G(s)] q is a polynomial of degree np equal to the number of compartments of the associated subsystem; as the coefficient of snP is always 1, np parameters can be estimated if the output compartment is open, np - 1 if it is closed. For the numerator associated with every ij pair, the following holds: i) if [CB] y= 1, the degree is np 1; as the coefficient of snp-, is 1, np - 1 parameters can be estimated; ii) if [CA'-'BI ii = [CA'-2BIi = ... = [CB]L = 0 and [CA'B] ii # 0, the degree is np - 1 - 1; as the coefficient of snP-1-1 is a function of {ki}, the number of parameters that can be estimated are np - 1. Notice that [CA'B] ii # 0 implies [CA"' B] i # 0. Two different subsystems may have a common part; the respective denominators have q common prime factors of the first order, where q is the number of common compartments. If the common part is in cascade with the remaining parts, the numerators too may have a common factor. If [A Jfg = [A Ifg = * * * = [Am -1 ]fg = 0 and [Am Ifg > 0, where g and f are respectively the first and the last compartment of the common part (with respect to the input), then there are q - m - 1 mutually dependent coefficients in the two numerators. Therefore the number of independent coefficients of the two [G(s)] ii equals n, +n2 - 2q +m + 1 (n +n2 - 2q+ m + 2, if one common prime factor of the denominators is s) where n, (n2) is the number of coefficients obtainable from the first (the second) function transfer. Notice that the common part may coincide with one of the two subsystems. For a strongly connected system all subsystems coincide [3]. Obviously, the system is structurally identifiable when the number of parameters that can be estimated from G(s) are not smaller than the number of the unknown system parameters. The existence of redundancy, that is, the existence of more independent equations than unknown parameters, does not influence structural identifiability, but may suggest a simpler experiment and/or a more complex model.
compartment j; the positive entries of the i-th row denote compartment i and the compartments from which compartment i can be reached. Matrices R 1B and CR1 give the following information: the positive entries of each column of R IB (each row of CR1 ) denote the compartments reachable (observable) from the corresponding input (output). The system is CC (CO) if each row of R 1B (each column of CR ) has at least one positive entry. 5.2 Identifiability The number of parameters that can be estimated by means of G(s) entries may be found through the following procedure, based on the concepts developed in section 4. To recognize subsystems, matrices R1B and CRI are put in Boolean form, (R1B)b and (CR1)b, respectively, and their product T = (CR I )b X (R 1B)b is computed. Each entry [T] q is the number of compartments of the ii subsystem. Two different subsystems ij and ml have a common part if there exists at least one compartment h such that: [CR lih
-
[CR'llmh* [RlB]h
-
[RlBIhl>O
(15)
and the common part consists of all compartments {h}. The common part in cascade {w} C {h} exists if and only if kiOViE {w}, Vji {w} except possibly if, and-k1ii0V i E {w}, Vj ¢ {w} except possibly j = g. 5.3 Flow Chart and Its Use In Fig. 3 a concise flow chart of the program for a computer implementation is reported. The flow chart may be used for compartmental systems of any structure. For strongly connected compartmental systems the flow chart may be simplified [31 owing to the peculiar characteristics of such a structure (see previous sections). The main operations are:
i) loading of K, Ko, B and C; ii) checking complete controllability and observability; iii) checking structural identifiability. The program, written in PLl language, was tested on a 370/158 IBM computer through a 2741 terminal. 5. DIGITAL COMPUTER IMPLEMENTATION OF A TEST ON When the chosen experiment does not allow to identify the CONTROLLABILITY, OBSERVABILITY AND some different procedures can be adopted. When posmodel, STRUCTURAL IDENTIFIABILITY sible, either a different experiment is chosen or a different 5.1 Controllability and Observability model is developed, such that it is identifiable by the chosen Once the structure of the model has been stated, matrices K experiment; otherwise the problem can be accurately analyzed and Ko can be put in Boolean form (kie = 1 if a transfer from to look for some relations among the parameters which are to jto iexists, ki1 = 0 if it does not, i = 1, n;j 1, n). be estimated. These possibilities are to be examined always From K, RO) is obtained. taking the biological phenomenon into account; sometimes it If the entries of R(n) are all positive, the system is strongly is impossible to change the input-output experiment; in some connected: as seen in section 3, complete controllability and other cases there is no further information about the parameters, not even about the possibility of neglecting some of them. observability always hold. For systems of any structure, conditions required by Theo6. EXAMPLES rems 1, 2.1 and 2.2 (section 3) should be verified. Matrix In this section some significant non strongly connected bio+ I, where I is the identity matrix of rank n, is comRR(1) I = puted. The positive entries of the j-th column denote com- logical compartmental systems are considered, and concepts partment i and the compartments that are reachable from and methods developed in the paper are applied to analyze
COBELLI AND ROMANIN-JACUR: BIOLOGICAL COMPARTMENTAL SYSTEMS
97
Fig. 3. Flow chart for checking controllability, observability and structural identifiability of compartmental systems.
their structural identifiability. Examples of structural identifiability test on strongly connected systems can be found in [3]. The models presented in the examples are neither discussed from the biological point of view (biological references for single models are given) nor the experimental methodology, as would be needed for practical applications (design of experiment, error sources, etc.), will be examined in detail. We will emphasize the aspects concerning the application of the developed mathematical theory, particularly those related to the possible modification of the model structure and/or of the input-output experiment to achieve structural identifiability. 6.1 Digoxin Distribution A three compartment model [11 ] is reported in Fig. 4, gastrointestinal tract-I,, plasma-2 and tissues-3. The identifi'cation experiment is performed with input (injection of radioactive tracer) in compartment 1 and output in 2. [0 00 K= k2l 0 k23 Ko= [O k02 01
LO From R(3) =
k32 0
j
1= K1 it can be seen that the system is not
[1
1 2 '
3
Fig. 4. Three compartment model of digoxin distribution (gastrointestinal tract-i, plasma-2 and tissues-3).
strongly connected as not all entries are positive (see sect. 2.1). Matrices A, B, C are: 0 0 1 [-k21 A= k21 -(k32+ko2) k23 B= 0 C= [O 1 01. L 0 k32 k2 LoiJ Theorems 1, 2.1 and 2.2 of sect. 3 assure the system is CC and CO (the test is performed on R 1B and CRI, respectively, see sect. 5.1). The system has one input and one output [G(s) is scalar], therefore neither cancellations nor simplifications occur. From G(s) denominator three pararneters can be estimated (open system) and two more from the numerator as CAB * 0 (see sect. 4). The system is structurally identifiable: in fact the chosen experiment allows the estimation of four
98
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, MARCH 1976
parameters with one redundant equation. Notice that no other experiment can be performed with the model structure of Fig. 4; in fact an input other than in compartment 1 makes the system non CC, while an output either in 1 or in 3 cannot be considered. 6.2. Albumin Distribution Fig. 5. Four compartment model of albumin distribution (extravascuA four compartment model [12] is reported in Fig. 5, extralar space with slow exchange-i, plasma-2, extravascular space with rapid exchange-3 and urine-4). vascular spaces with slow exchange-i, plasma-2, extravascular may spaces with rapid exchange-3 and urine-4. An experiment be performed with input (injection of radioactive tracer) in 2 and output in 4. Following the procedure outlined in the previous example, it can be seen that the system is CC, CO and structurally identifiable with one redundant equation. If also an output in 2 is provided the system is obviously CC, Fig. 6. Three compartment model of albumin distribution. CO and structurally identifiable. In this case two subsystems must be considered, the first coinciding with the whole system and the second made up of compartments 1, 2, 3; the two subsystems have a common part in cascade, therefore simplifications occur between denominators and numerators of the respective transfer functions (see sects. 4 and 5.2). This new experiment leads to two redundant equations. A different identification experiment could be tried with input and output in 2, but in this case the model of Fig. 5 is non Fig. 7. Four compartment model of iodine distribution (thyroid-i, plasma-2, tissues-3 and urine-4). CO. However, a different model having the same dynamical strucmodel in 6. This new as Fig. be adopted behavior may ture makes the system CC, CO and structurally identifiable without redundant equations. Generally, when different experiments enable a compartmental system to be identified, as in the present example, the one to be performed is to be chosen on the basis of the technical and experimental conditions and the reliability of the performed measurements. Fig. 8. Four compartment model of thyroxine distribution (plasma-i, urine-2, extravascular reversible space-3 and extravascular non 6.3. Iodine Distribution reversible space-4). A four compartment model [13], [14] is reported in Fig. 7, thyroid-i, plasma-2, tissues-3 and urine-4. An experiment may be performed with input (injection of radioactive tracer) in 2 and output in 4. By applying the procedures of sects. 4 and 5, it can be found that the system is CC, CO and structurally identifiable with one redundant equation. If also an output in 1 is provided, the system is obviously CC, CO and structurally identifiable. Two subsystems must be Fig. 9. Three compartment model of thyroxine distribution. considered as in the previous example: after simplifications redundant equations result to be three. An experiment with input in 2 and output in 1 allows struc- ment is performed with input (injection of radioactive tracer} tural identifiability of the model, only if its structure is in 1 and output in 2, the system is CC, CO and structurally changed as in Fig. 6 of the previous example (unchanged dy- identifiable with no redundant equations. A second output in compartment 1 could be provided. Apnamical behavior). This new configuration gives no redundant plying the procedure of sects. 4 and 5, it can be seen that beequations. cause of simplifications no redundant equations are achieved 6.4 Thyroxine Distribution (only an increase of the number of observations is obtained). Considering the model of Fig. 8, one may be interested only A four compartment model [15] is reported in Fig. 8, plasma-i, urine-2, extravascular reversible spaces-3 and extra- in the exchanges of material between compartments 1 and 3: vascular non reversible spaces-4. Notice that compartment 4 the new model structure of Fig. 10 is adopted where the dycannot be observed; therefore the system is non CO with any namical relations between 1 and 3 do not change, while exexperiment. However, the different model of Fig. 9 with un- changes 1 to 4 and 1 to 2 are no longer distinguishable. Perchanged dynamical behavior may be adopted; if the experi- forming an experiment with input and output in 1, the system
COBELLI AND ROMANIN-JACUR: BIOLOGICAL COMPARTMENTAL SYSTEMS
{0
99
3
1
Fig. 10. Two compartment model of thyroxine distribution.
and structurally identifiable without redundant equations, but only the dynamical relations between 1, 2 and 3 can be determined, while exchanges from 1 to 4 and 5 cannot be distinguished. 7. CONCLUSIONS Fig. 11. Five compartment model of calcium metabolism (plasma-1, bones with rapid exchange-2, bones with slow exchange-3, urine-4 and faeces-5).
Fig. 12. Four compartment model of calcium metabolism.
Fig. 13. Four compartment model of calcium metabolism.
Fig. 14. Three compartment model of calcium metabolism.
is CC, CO and structurally identifiable without redundant equations. 6.5. Calcium Metabolism A five compartment model [16] is reported in Fig. 11, plasma-i, bones with rapid exchange-2, bones with slow exchange-3, urine4 and faeces-5. The experiment is performed with input (injection of radioactive tracer) in 1. If the outputs are taken in 1, 4 and 5 the system is CC, CO and structurally identifiable with one redundant equation. If only two compartments are observed, 1 and 4 or 1 and 5, the structure of the model is to be changed as in Figs. 12 or 13, respectively, the system becomes CC, CO and structurally identifiable with no redundant equations. If only one output in 1 is provided, the model structure must be modified as in Fig. 14. The system becomes CC, CO
The problem of structural identifiability of biological compartmental models with any structure has been analyzed. As is well known, compartmental models are widely used both in research and diagnosis problems. Input-output tracer experiments employed for their identification are often non repeatable because of induced harm, trouble or danger: therefore the problem of an a priori evaluation of the possibility of estimating all parameters of the adopted model (structural identifiability) via a given input-output experiment is of great importance. In fact, if the chosen model cannot be identified by the chosen experiment, either the experiment is to be modified or a different model (on condition it is still biologically significant) is to be adopted. Structural identifiability is studied directly connecting it to controllability and observability in a system theory context. After some review about fundamentals on biological compartmental models, new theorems relating controllability and observability to the model structure are presented. Then structural identifiability of controllable and observable compartmental systems is analyzed. Finally, techniques for checking controllability, observability and identifiability are given together with a concise flow chart for a digital computer implementation. Some examples illustrating the application of the developed mathematical theory are considered in order to emphasize the importance of the problem of structural identifiability in biology and medicine. In particular, identification of classical compartmental models via commonly used tracer experiments is analyzed and modifications of model structure and/or of input-output configuration are considered. APPENDIX
Theorem 1: Proof: by contradiction: With reference to the complete controllability of a system of n compartments, let the compartments numbered from 1 to rb be controllable, that is:
B=
0 1 .O---O 0 O 0 1 ....o ..........o ............
O O O *---
..........1 ooo......
1
.(16)
Let us suppose compartment n is not reachable along any path from any controllable compartment.
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, MARCH 1976
100
It can be seen, by inspection, that:
[A']nijO,
1= l,n - rb j= l,rb. As an example, for 1 = 1 and 1= 2, (17) is expressed tively by (18) and (19):
[A]ni =anj =knj0, as
there is
respec-
(18)
j=l,rb
path of length
no
(17)
1
partment n;
from compartment j to
com-
n
[A2]nj=
anr
ar =
r=i
+ (ann +
as both the first paths of length 1
rj r*n
knr krj
(19)
aja) knO, -
sum and kn, are zero, as there are neither nor of length 2 from compartment j to com-
partment n.
The last row of controllability matrix P is:
[B nl 1I... [B] nrb I[p4n i I -* [A Inrb [A2 Inl 1 |. [A2] 1nb * 1 [An rb Inl 1 .. I [A n-b]Inrb.- (20) From (16) and (17) it follows that all the entries of the last Q.E.D. of P are zero, therefore its rank is less than n. With reference to the complete observability, let the observed compartments be numbered from 1 to rc. Let us suppose there is no path from compartment n to any observable compartment. It can be seen, by inspection, that the last row of observability matrix Q is null, therefore its rank is less Q.E.D. than n. 7Theorem 2.1 Proof: It is necessary to prove that, under such hypotheses, matrix P cannot have a rank less than n for any combination of the parameter values {kj,} or equally that det [PPTI is not identically zero. It will be proved that [PPTI has neither its rows or columns identically null nor linearly dependent. Let the controlled compartments be numbered from I to rb. The i-th row (i > rb) of matrix P is expressed by: row
0 1 [A I il l
O *.
I
[A n-rb]il
I [A I irb I [A' I il l *l[X2 ] irbl | [An-b]irb
(21)
while for i < rb expression (21) changes only in that [PI ii = 1. Note by inspection that at least one entry of every row i > rb is non-null as there is at least one path not exceeding the length n - rb from one of the controlled compartments to every uncontrolled one. Obviously the same holds for every row i < rb The general entry on the main diagonal of [PPTI is expressed by:
[PP'l ii
n-rb
=
[PPTI ii =
E,
1=1
n-rb
E
rb
E
m =1
[A' ]m + 1,
i1
rb
i>
rb
rb
[A']mIM X
E
1=1 M=l
It is evident that the entries
on
.
the main diagonal cannot be
zero, therefore matrix ppT has neither null rows nor columns. Moreover, there is no linear dependence among the rows as every entry of a row is a sum of monomials one of which does not appear in any other entry of the same column, except for the case in which two or more rows have only one non zero element on the same column. In this last case the corresponding compartments are closed and only influenced by the same compartment; the state variables cannot vary independently, but only proportionally: however, that originates no trouble for identification, and therefore the system will be considered Q.E.D. CC all the same. Theorem 2.2 Proof: Let the observed compartments be numbered from 1 to r,. Following the same procedure used in Theorem 2.1, it can be proved that det [QQT] cannot be identically zero. ACKNOWLEDGMENT The authors wish to thank Prof. A. Lepschy, Istituto di Elettrotecnica e di Elettronica, Facolta di Ingegneria della Universita di Padova for his interest and his helpful suggestions. Thanks are also due to Prof. A. Rescigno, Department of Physiology, University of Minneapolis, Minnesota, for the careful reading of the manuscript.
REFERENCES [1] R. Bellman and K. J. Astrom, "On structural identifiability," Math. Biosci., vol. 7, pp. 329-339, 1970. 121 M. Hajek, "A contribution to the parameters estimation of a certain class of dynamical systems," Kybernetika, vol. 8, pp. 165-
172, 1972. [31 C. Cobelli and G. Romanin-Jacur, "Structural identifiability of strongly connected biological compartmental systems," Medical & Biological Engineering, vol. 13, pp. 831-838, 1975. [4] A. Rescigno and G. Segre, Drug and Tracer Kinetics. Waltham: Blaisdell Publ. Co., 1966. [5 1 J. A. Jacquez, Compartmental Analysis in Biology and Medicine. Amsterdam: Elsevier Publ. Co., 1972. [61 J. Z. Hearon, "Theorems on linear systems," Ann. N. Y. Acad. Sci., vol. 108, pp. 3648, 1963. [71 L. A. Zadeh and C. A. Desoer, Linear System Theory. New York: McGraw-Hill, 1963. [81 R. E. Kalman, "Mathematical description of linear dynamical system," J. SIAM. Control, vol. 1, pp. 152-192, 1963. [9] C. T. Chen, Introduction to Linear System Theory. New York: Holt, Rinehart and Winston, 1970. [10] F. R. Gantmacher, The Theory of Matrices. New York: Chelsea Publ. Co., 1959, vol. 2, pp. 97-98. [11] B. C. McInnis, S. El-Asfouri and T. E. Bynum, "Modeling and computer simulation to determine digoxin dosage," Proceedings of the Regulation and Control in Physiological Systems Conference, Rochester, N.Y., August 22-24, 1973, Eds. A. S. Iberall and A. C. Guyton, Pittsburgh: I.S.A., 1973, pp. 69-70. [12] T. Freeman and C. M. E. Matthews, "Analysis of the behaviour of 131I-albumin in the normal subject and nephrotic patient," in Radioactive Isotope in Klinik und Forschung, Eds. K. von Fellinger and H. Vetter, Munich: Urban & Schwarzenberg, 1958, vol. 3, p. 283. [13] J. Rotblat and R. Marcus, "Clinical uses of radioiodine," in Radioisotope Techniques, Proceedings of the Isotope Techniques Conference, Oxford, July 1951, H.M.S.O., 1953, vol. 1, p. 33. [141 G. L. Brownell, "Analysis of techniques for determination of thyroid function with radioiodine," J. Clin. Endocrin., vol. 11, p. 1095, 1951. [151 B. Bloomstedt and L. 0. Plantin, "The extrathyroid distribution of 13 1 I-thyroxine," Acta Endocr., Copnh., vol. 48, p. 5 36, 1965. [161 S. H. Cohn, S. R. Bozzo, J. Jesseph, C. Constantinides, D. R. Huene and E. A. Gusmano, "Formulation and testing of a compartmental model for calcium metabolism in man," Radiat. Res., vol. 26, p. 319, 1965.
93
Controllability, Observability and Structural Identifiability of Multi Input and Multi Output Biological Compartmental Systems CLAUDIO COBELLI,
MEMBER, IEEE, AND
Abstract-Controllability, observability and structural identifiability of biological compartmental systems of any structure are studied in order to evaluate the significance of a multi input-multi output experiment, that is, the a priori possibility of estimating the unknown system parameters through the chosen experiment. Some basic concepts and properties of biological compartmental systems and corresponding models are rirst briefly reviewed to put compartmental models in a linear system theory context. Controllability, observability and structural identifiability are defined for compartmental systems. New theorems on compartmental systems relating controllability and observability to their structure are proved. Theorems on transfer function matrix are suitably applied to solve the problem of structural
identifiability. Solution techniques to check controliability, observability and structural identifiability of compartmental systems of any structure are given. The flow chart for a digital computer implementation of these techniques is presented. Some examples concerning tracer analysis experiments on biological compartmental systems concerning digoxin, iodine, albumin and thyroxine distribution, and calcium metabolism illustrate the application of the mathematical theory developed in the paper in order to emphasize the importance of the problem of structural identifiability in biology and medicine.
1. INTRODUCTION
IN tracer analysis of biological systems, compartmental models have been widely adopted in research, diagnosis and therapy problems. These models approximate the quantitative description of phenomena, thus allowing many parameters of biological interest to be estimated. The most widely treated problems in biomedical literature regarding compartmental models concern compartmentalization of biological systems and parameter estimation, that is, the development of models using a compartment approach and the numerical estimation of the associated parameters from
input-output tracer experiments. Actually the numerical estimation of parameters should follow a test on the a priori possibility of their estimation by the chosen input-output experiment, so that the significance of the experiment itself is tested. This test will be referred to as the solution of the problem of structural identifiability, as it is related only to the model structure; the solution determines which and/or how many parameters of the adopted model can be estimated by means of the chosen input-output experiment. Manuscript received August 2, 1974; revised April 17, 1975. C. Cobelli is with the Laboratorio per Ricerche di Dinanica dei Sistemi e di Elettronica Biomedica, Consiglio Nazionale delle Ricerche, Casella Postale 1075, 35100 Padova, Italy, and the Istituto di Elettrotecnica e d'Elettronica, Universiti di Padova, Italy. G. Romanin-Jacur is with the Laboratorio per Ricerche di Dinamica dei Sistemi e di Elettronica Biomedica, Consiglio Nazionale delle Ricerche, Padova, Italy.
GIORGIO ROMANIN-JACUR
This problem has been considered only recently [1], [2], [3]. Bellman et al. in 1970 [1] put the problem in a system theory context and stated precise general definitions. Hayek in 1972 [2] provided solution techniques for linear invariant compartmental systems open to the environment; solution is given only for the one input-one output case (not directly extendible to the multi input-multi output case) under complete controllability and complete observability conditions. In [3], the authors study structural identifiability for multi inputmulti output strongly connected compartmental systems both open and closed. Theorems relating controllability and observability to the system structure are proved and solution techniques for structural identifiability are developed; a flow chart for a digital computer implementation is presented and several application examples on biological systems are reported. The present paper studies structural identifiability for multi input-multi output compartmental linear invariant models of any structure, that is, the limiting restriction to the class of strongly connected systems [3] is removed. Briefly, in section 2 some general notes about compartmental models and their structural properties (particularly controllability, observability and structural identifiability) are given. In section 3 new theorems relating controllability and observability of compartmental systems to their structure are presented (proofs are given in the appendix). The general solution of the structural identifiability problem is discussed in section 4. In section 5 some testing criteria and the flow chart for the implementation on digital computer of the problem considered in 3 and 4 are given. Finally, in section 6 several compartmental systems of biomedical interest are analyzed from the point of view of their structural identifiability. 2. BASIC NOTES ON BIOLOGICAL COMPARTMENTAL SYSTEMS
2.1 CompartmentalSystems andModels Some basic concepts, useful for the following, will be briefly reviewed; a more comprehensive treatment can be found for instance in [4], [5]. In a biological compartmental system small deviations from its steady state are induced from the outside; under the assumptions of volume constancy for each compartment, linearity and invariance, the dynamical behavior of each compartment is described by (Fig. 1):
Xi= pki £ kixj+ ui- j*iE kiixi- koixi
O e wher where:
(1)
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, MARCH 1976
94
ai = kij,
Ui
i*j
aii = -koi-
n j=1
(4)
ki.
(5)
The eigenvalues of matrix A cannot be purely imaginary and
their real part is non positive (it can be zero when the eigenvalue is real) [6]. Matrix B is of n X rb dimension, with rb < n; if inputs are in compartments 11, i2, ib its entries can assume either value 1 or 0 under the following condition:
xoi Xi Fig. 1. Input
output material fluxes of the generic i-th compartment.
and
E bij= I Vj;
bi
(6)
,
Matrix C is of r, X n dimension with r, S n; if outputs are * h, its entries can assume taken from compartments il, either value 1 or value 0 under the following condition: 12,
E
j=l
F0 x
XO 07
xl0
ox 0
x
o x
o
LO
01
o
Fig. 2. K and Ko matrices of a compartmental system.
is the state variable associated with the i-th compartment, that is, the amount (or concentration) of material of the i-th compartment; equation (1) is also valid when xi is the amount (or specific activity) of the labeled material (in this case the system may be non linear); kii is the rate constant from the j-th compartment to the i-th compartment; koi is the rate constant from the i-th compartment to the environment; ui is the external input.
xi
The rate constants {kij} are grouped into a square matrix K of n rows, where n is the compartment number, and into a row matrix Ko of n entries. Notice that all entries of matrix K are non negative and those of the main diagonal are zero; all entries of Ko are non negative (Fig. 2). The mathematical model of a compartmental system corresponds to the usual input-state-output representation of a dynamical linear invariant system: x
=Ax + Bu
y =
Cx
(2) (3)
where x, u, y are respectively state, input and output vectors. With reference to the structure of the compartrnental system matrices A, B, C show peculiar characteristics. The order n of matrix A equals the number of compartments; its entries are related to those of the matrices K and Ko
by:
Ci1 =1 Vi;
E C.j =
0l,
pj] j1
,
.*'
C
(7)
Compartments are represented graphically either by blocks or by nodes and transferences by directed segments called branches. A path is a succession of branches such that the terminal block (node) of a branch is also the initial block (node) of the next branch. A compartmental system can be divided into subsystems made up of one or more compartments, such that every compartment in a subsystem is directly connected with at least one compartment in the same subsystem. Any subsystem S is characterized by the corresponding matrix A,. A compartmental system is said to be disjoint when it contains one or more isolated subsystems not interacting with the remaining part of the system. In the following, only non disjoint systems will be considered. A subsystem is said open (closed) when there is some exchange (no exchange) of material from the subsystem either to the environment or to another subsystem. A subsystem is strongly connected when every compartment can be reached from every other compartment in the same subsystem along at least one path. If matrix KS of the rate constants of subsystem S is considered, it can be noticed that the general entry of KS (Ks at r-th power) is given by: (8) kil I ki 12 k Ir-J5 [Kr, ij = E .
11
a2,
r-1 C
-
-
1
and then [Ks] ij > 0 implies that at least one path of length r exists from compartment i to compartment i within S. If matrix Rim) = I7. K' is considered, it can be noticed that [R(m)] iq > 0 implies that at least one path of length smaller than m, or equal to m, exists from compartment j to compartment i. If m = ns (where n. is the number of compartments in subsystem S) matrix R('s) gives the following information: a) positive entries of the j-th column denote the compartments reachable from compartment j; b) positive entries of the i-th row denote the compartments from which comnpartment i can be reached. From the previous remarks, it follows that matrix R ns) of a strongly connected subsystem has all entries positive.
95
COBELLI AND ROMANIN-JACUR: BIOLOGICAL COMPARTMENTAL SYSTEMS
In all previous definitions and remarks subsystem S may coincide with the whole compartmental system.
2.2 Controllability, Observability and Identifiability in a Structural Sense of Compartmental Systems Rigorous treatments on structural properties can be found in system theory literature (see for instance [7], [8], [9]). In the following, only tests on structural properties will be considered. In the theory of linear invariant systems as the compartmental ones considered, a system is said to be completely controllable (CC) if and only if the controllability matrix of n X (n - rb + 1) * rb dimension: (9) P = [BIABIA 2BI ... {An rbB] is of rank n, and is completely observable (CO) if and only if the observability matrix of n X (n - r, + 1) r, dimension: (10) Q = [CT IATCTI A2 TCTI * *An-rcT CT] is of rank n, where n is the order of the system (compartments number), rb and rc are the ranks of matrices B and C (inputs and outputs number, respectively). Conditions (9) and (10) correspond to the following ones: -
3. THEOREMS ON CONTROLLABILITY AND OBSERVABILITY IN A STRUCTURAL SENSE OF COMPARTMENTAL SYSTEMS In this section three new theorems relating controllability and observability of compartmental systems to their structure are presented. Proofs are reported in the Appendix. Theorem 1: The existence of at least one path reaching every uncontrolled compartment from a controlled one is a necessary condition for a compartmental system to be CC; the existence of at least one path from every unobserved compartment to an observed one is a necessary condition for a compartmental system to be CO. Note that if the system is strongly connected the necessary conditions of theorem 1 are always satisfied for any input and output; in fact, by definition, every compartment can be reached from every other one along at least one path. Theorem 2.1. A compartmental system is CC in a structural sense only if there is at least one path from a controlled compartment to every uncontrolled one. Theorem 2.2. A compartmental system is CO in a structural sense if there is at least one path from every unobservable compartment to an observed one. Theorems 2.1 and 2.2 always hold for strongly connected compartmental systems (see definition in 2.1) [3 ]. Theorems 2.1 and 2.2 allow us to state that the rank of matrices P and Q can be less than n only in a subspace of the { kie > 0} space. In fact, in the case of one input-one output systems, each of the equations det [ppT] = 0 and det [QQT] = 0 is the general expression of an hypersurface, therefore the non CC and non CO subspaces are represented by hypersurfaces. For -the multi input-multi output case, there are as many hypersurface equations as the submatrices of rank n of matrix [ppT] and [QQT], respectively. Non CC and non CO subspaces are then represented by the intersections of the hypersurfaces described by the equations relative to the considered submatrices of [PPT] and [QQT] , respectively.
(11) det[PPT] 0 (12) det [QQTI 0O. A compartmental system can be non CC and/or non CO either because of its own structure or because of a particular combination of the parameter values {k,j}. In the first case the rank of matrix P and/or Q is always less than n for whatever values of rate constants relative to the system of fixed structure (unchanged topology); it means there are no possible connections between a certain input and a certain state or between a certain state and a certain output. A compartmental system is said to be non CC (non CO) in a structural sense when it is non CC (non CO) for any values 4. STRUCTURAL IDENTIFIABILrrY OF combination of the non zero parameters {k,1}. COMPARTMENTAL SYSTEMS Input-output relations obtained from (2) and (3) are considered in order to define structural identifiability. As seen in sect. 2.2, a controllable and observable dynamical The input-output behavior of a compartmental system is linear invariant system is completely described by the transfer completely described, relatively to its CC and CO part, by the function matrix G(s) [see (13)]. The general entry [G(s)] i transfer function matrix G(s): where i = 1, r, and j = 1, rb is the transfer function between input j and output i. To find out the number of coefficients I A C * adj(sI- A) * B. G(s) = C (sI - A)-1 * B = det (sI- A) that can be estimated, the expansions of det (sI -A) and C adj(sI - A) * B are considered ([9], [10]; see also [3] for (13) details) and G(s) is given by: The knowledge of G(s) allows the estimation of a maximum number of parameters of the CC and CO part of the system [CB(sn1 +a1 sn -2 + - + an1) G(s) = (5 equal to the number of independent coefficients of the numerator and denominator polynomials of the entries of G(s). + CAB (sn 2 + a,X Sn -+ +n -2a ) A CC and CO compartmental system is said to be structur(14) ally identifiable when the chosen input-output experiment + +CAn-2B(s+a,)+CAn -B]. allows the estimation of all unknown parameters (rate [G(s)] i, has the same expression as G(s) if we take respectively constants). If the whole compartmental system is not both CC and CO, [CB] i, [CAB] ii, *- [C4n -1 B] ji. The number of parameters that can be estimated are related complete controllability and observability may hold for a part of the system; in this case the considerations on structural to the polynomials in s of the numerator and denominator of identifiability are restricted to the CC and CO part of the every [G(s)],,. Cancellations occur whenever the subsystem CC through i and CO through i is not equal to the whole comsystem. -
-
-
96
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, MARCH 1976
partmental system: more precisely, the number of simplifiable prime factors equals the number of the compartments uncontrollable through i and/or unobservable through i. Therefore the denominator associated with every [G(s)] q is a polynomial of degree np equal to the number of compartments of the associated subsystem; as the coefficient of snP is always 1, np parameters can be estimated if the output compartment is open, np - 1 if it is closed. For the numerator associated with every ij pair, the following holds: i) if [CB] y= 1, the degree is np 1; as the coefficient of snp-, is 1, np - 1 parameters can be estimated; ii) if [CA'-'BI ii = [CA'-2BIi = ... = [CB]L = 0 and [CA'B] ii # 0, the degree is np - 1 - 1; as the coefficient of snP-1-1 is a function of {ki}, the number of parameters that can be estimated are np - 1. Notice that [CA'B] ii # 0 implies [CA"' B] i # 0. Two different subsystems may have a common part; the respective denominators have q common prime factors of the first order, where q is the number of common compartments. If the common part is in cascade with the remaining parts, the numerators too may have a common factor. If [A Jfg = [A Ifg = * * * = [Am -1 ]fg = 0 and [Am Ifg > 0, where g and f are respectively the first and the last compartment of the common part (with respect to the input), then there are q - m - 1 mutually dependent coefficients in the two numerators. Therefore the number of independent coefficients of the two [G(s)] ii equals n, +n2 - 2q +m + 1 (n +n2 - 2q+ m + 2, if one common prime factor of the denominators is s) where n, (n2) is the number of coefficients obtainable from the first (the second) function transfer. Notice that the common part may coincide with one of the two subsystems. For a strongly connected system all subsystems coincide [3]. Obviously, the system is structurally identifiable when the number of parameters that can be estimated from G(s) are not smaller than the number of the unknown system parameters. The existence of redundancy, that is, the existence of more independent equations than unknown parameters, does not influence structural identifiability, but may suggest a simpler experiment and/or a more complex model.
compartment j; the positive entries of the i-th row denote compartment i and the compartments from which compartment i can be reached. Matrices R 1B and CR1 give the following information: the positive entries of each column of R IB (each row of CR1 ) denote the compartments reachable (observable) from the corresponding input (output). The system is CC (CO) if each row of R 1B (each column of CR ) has at least one positive entry. 5.2 Identifiability The number of parameters that can be estimated by means of G(s) entries may be found through the following procedure, based on the concepts developed in section 4. To recognize subsystems, matrices R1B and CRI are put in Boolean form, (R1B)b and (CR1)b, respectively, and their product T = (CR I )b X (R 1B)b is computed. Each entry [T] q is the number of compartments of the ii subsystem. Two different subsystems ij and ml have a common part if there exists at least one compartment h such that: [CR lih
-
[CR'llmh* [RlB]h
-
[RlBIhl>O
(15)
and the common part consists of all compartments {h}. The common part in cascade {w} C {h} exists if and only if kiOViE {w}, Vji {w} except possibly if, and-k1ii0V i E {w}, Vj ¢ {w} except possibly j = g. 5.3 Flow Chart and Its Use In Fig. 3 a concise flow chart of the program for a computer implementation is reported. The flow chart may be used for compartmental systems of any structure. For strongly connected compartmental systems the flow chart may be simplified [31 owing to the peculiar characteristics of such a structure (see previous sections). The main operations are:
i) loading of K, Ko, B and C; ii) checking complete controllability and observability; iii) checking structural identifiability. The program, written in PLl language, was tested on a 370/158 IBM computer through a 2741 terminal. 5. DIGITAL COMPUTER IMPLEMENTATION OF A TEST ON When the chosen experiment does not allow to identify the CONTROLLABILITY, OBSERVABILITY AND some different procedures can be adopted. When posmodel, STRUCTURAL IDENTIFIABILITY sible, either a different experiment is chosen or a different 5.1 Controllability and Observability model is developed, such that it is identifiable by the chosen Once the structure of the model has been stated, matrices K experiment; otherwise the problem can be accurately analyzed and Ko can be put in Boolean form (kie = 1 if a transfer from to look for some relations among the parameters which are to jto iexists, ki1 = 0 if it does not, i = 1, n;j 1, n). be estimated. These possibilities are to be examined always From K, RO) is obtained. taking the biological phenomenon into account; sometimes it If the entries of R(n) are all positive, the system is strongly is impossible to change the input-output experiment; in some connected: as seen in section 3, complete controllability and other cases there is no further information about the parameters, not even about the possibility of neglecting some of them. observability always hold. For systems of any structure, conditions required by Theo6. EXAMPLES rems 1, 2.1 and 2.2 (section 3) should be verified. Matrix In this section some significant non strongly connected bio+ I, where I is the identity matrix of rank n, is comRR(1) I = puted. The positive entries of the j-th column denote com- logical compartmental systems are considered, and concepts partment i and the compartments that are reachable from and methods developed in the paper are applied to analyze
COBELLI AND ROMANIN-JACUR: BIOLOGICAL COMPARTMENTAL SYSTEMS
97
Fig. 3. Flow chart for checking controllability, observability and structural identifiability of compartmental systems.
their structural identifiability. Examples of structural identifiability test on strongly connected systems can be found in [3]. The models presented in the examples are neither discussed from the biological point of view (biological references for single models are given) nor the experimental methodology, as would be needed for practical applications (design of experiment, error sources, etc.), will be examined in detail. We will emphasize the aspects concerning the application of the developed mathematical theory, particularly those related to the possible modification of the model structure and/or of the input-output experiment to achieve structural identifiability. 6.1 Digoxin Distribution A three compartment model [11 ] is reported in Fig. 4, gastrointestinal tract-I,, plasma-2 and tissues-3. The identifi'cation experiment is performed with input (injection of radioactive tracer) in compartment 1 and output in 2. [0 00 K= k2l 0 k23 Ko= [O k02 01
LO From R(3) =
k32 0
j
1= K1 it can be seen that the system is not
[1
1 2 '
3
Fig. 4. Three compartment model of digoxin distribution (gastrointestinal tract-i, plasma-2 and tissues-3).
strongly connected as not all entries are positive (see sect. 2.1). Matrices A, B, C are: 0 0 1 [-k21 A= k21 -(k32+ko2) k23 B= 0 C= [O 1 01. L 0 k32 k2 LoiJ Theorems 1, 2.1 and 2.2 of sect. 3 assure the system is CC and CO (the test is performed on R 1B and CRI, respectively, see sect. 5.1). The system has one input and one output [G(s) is scalar], therefore neither cancellations nor simplifications occur. From G(s) denominator three pararneters can be estimated (open system) and two more from the numerator as CAB * 0 (see sect. 4). The system is structurally identifiable: in fact the chosen experiment allows the estimation of four
98
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, MARCH 1976
parameters with one redundant equation. Notice that no other experiment can be performed with the model structure of Fig. 4; in fact an input other than in compartment 1 makes the system non CC, while an output either in 1 or in 3 cannot be considered. 6.2. Albumin Distribution Fig. 5. Four compartment model of albumin distribution (extravascuA four compartment model [12] is reported in Fig. 5, extralar space with slow exchange-i, plasma-2, extravascular space with rapid exchange-3 and urine-4). vascular spaces with slow exchange-i, plasma-2, extravascular may spaces with rapid exchange-3 and urine-4. An experiment be performed with input (injection of radioactive tracer) in 2 and output in 4. Following the procedure outlined in the previous example, it can be seen that the system is CC, CO and structurally identifiable with one redundant equation. If also an output in 2 is provided the system is obviously CC, Fig. 6. Three compartment model of albumin distribution. CO and structurally identifiable. In this case two subsystems must be considered, the first coinciding with the whole system and the second made up of compartments 1, 2, 3; the two subsystems have a common part in cascade, therefore simplifications occur between denominators and numerators of the respective transfer functions (see sects. 4 and 5.2). This new experiment leads to two redundant equations. A different identification experiment could be tried with input and output in 2, but in this case the model of Fig. 5 is non Fig. 7. Four compartment model of iodine distribution (thyroid-i, plasma-2, tissues-3 and urine-4). CO. However, a different model having the same dynamical strucmodel in 6. This new as Fig. be adopted behavior may ture makes the system CC, CO and structurally identifiable without redundant equations. Generally, when different experiments enable a compartmental system to be identified, as in the present example, the one to be performed is to be chosen on the basis of the technical and experimental conditions and the reliability of the performed measurements. Fig. 8. Four compartment model of thyroxine distribution (plasma-i, urine-2, extravascular reversible space-3 and extravascular non 6.3. Iodine Distribution reversible space-4). A four compartment model [13], [14] is reported in Fig. 7, thyroid-i, plasma-2, tissues-3 and urine-4. An experiment may be performed with input (injection of radioactive tracer) in 2 and output in 4. By applying the procedures of sects. 4 and 5, it can be found that the system is CC, CO and structurally identifiable with one redundant equation. If also an output in 1 is provided, the system is obviously CC, CO and structurally identifiable. Two subsystems must be Fig. 9. Three compartment model of thyroxine distribution. considered as in the previous example: after simplifications redundant equations result to be three. An experiment with input in 2 and output in 1 allows struc- ment is performed with input (injection of radioactive tracer} tural identifiability of the model, only if its structure is in 1 and output in 2, the system is CC, CO and structurally changed as in Fig. 6 of the previous example (unchanged dy- identifiable with no redundant equations. A second output in compartment 1 could be provided. Apnamical behavior). This new configuration gives no redundant plying the procedure of sects. 4 and 5, it can be seen that beequations. cause of simplifications no redundant equations are achieved 6.4 Thyroxine Distribution (only an increase of the number of observations is obtained). Considering the model of Fig. 8, one may be interested only A four compartment model [15] is reported in Fig. 8, plasma-i, urine-2, extravascular reversible spaces-3 and extra- in the exchanges of material between compartments 1 and 3: vascular non reversible spaces-4. Notice that compartment 4 the new model structure of Fig. 10 is adopted where the dycannot be observed; therefore the system is non CO with any namical relations between 1 and 3 do not change, while exexperiment. However, the different model of Fig. 9 with un- changes 1 to 4 and 1 to 2 are no longer distinguishable. Perchanged dynamical behavior may be adopted; if the experi- forming an experiment with input and output in 1, the system
COBELLI AND ROMANIN-JACUR: BIOLOGICAL COMPARTMENTAL SYSTEMS
{0
99
3
1
Fig. 10. Two compartment model of thyroxine distribution.
and structurally identifiable without redundant equations, but only the dynamical relations between 1, 2 and 3 can be determined, while exchanges from 1 to 4 and 5 cannot be distinguished. 7. CONCLUSIONS Fig. 11. Five compartment model of calcium metabolism (plasma-1, bones with rapid exchange-2, bones with slow exchange-3, urine-4 and faeces-5).
Fig. 12. Four compartment model of calcium metabolism.
Fig. 13. Four compartment model of calcium metabolism.
Fig. 14. Three compartment model of calcium metabolism.
is CC, CO and structurally identifiable without redundant equations. 6.5. Calcium Metabolism A five compartment model [16] is reported in Fig. 11, plasma-i, bones with rapid exchange-2, bones with slow exchange-3, urine4 and faeces-5. The experiment is performed with input (injection of radioactive tracer) in 1. If the outputs are taken in 1, 4 and 5 the system is CC, CO and structurally identifiable with one redundant equation. If only two compartments are observed, 1 and 4 or 1 and 5, the structure of the model is to be changed as in Figs. 12 or 13, respectively, the system becomes CC, CO and structurally identifiable with no redundant equations. If only one output in 1 is provided, the model structure must be modified as in Fig. 14. The system becomes CC, CO
The problem of structural identifiability of biological compartmental models with any structure has been analyzed. As is well known, compartmental models are widely used both in research and diagnosis problems. Input-output tracer experiments employed for their identification are often non repeatable because of induced harm, trouble or danger: therefore the problem of an a priori evaluation of the possibility of estimating all parameters of the adopted model (structural identifiability) via a given input-output experiment is of great importance. In fact, if the chosen model cannot be identified by the chosen experiment, either the experiment is to be modified or a different model (on condition it is still biologically significant) is to be adopted. Structural identifiability is studied directly connecting it to controllability and observability in a system theory context. After some review about fundamentals on biological compartmental models, new theorems relating controllability and observability to the model structure are presented. Then structural identifiability of controllable and observable compartmental systems is analyzed. Finally, techniques for checking controllability, observability and identifiability are given together with a concise flow chart for a digital computer implementation. Some examples illustrating the application of the developed mathematical theory are considered in order to emphasize the importance of the problem of structural identifiability in biology and medicine. In particular, identification of classical compartmental models via commonly used tracer experiments is analyzed and modifications of model structure and/or of input-output configuration are considered. APPENDIX
Theorem 1: Proof: by contradiction: With reference to the complete controllability of a system of n compartments, let the compartments numbered from 1 to rb be controllable, that is:
B=
0 1 .O---O 0 O 0 1 ....o ..........o ............
O O O *---
..........1 ooo......
1
.(16)
Let us suppose compartment n is not reachable along any path from any controllable compartment.
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, MARCH 1976
100
It can be seen, by inspection, that:
[A']nijO,
1= l,n - rb j= l,rb. As an example, for 1 = 1 and 1= 2, (17) is expressed tively by (18) and (19):
[A]ni =anj =knj0, as
there is
respec-
(18)
j=l,rb
path of length
no
(17)
1
partment n;
from compartment j to
com-
n
[A2]nj=
anr
ar =
r=i
+ (ann +
as both the first paths of length 1
rj r*n
knr krj
(19)
aja) knO, -
sum and kn, are zero, as there are neither nor of length 2 from compartment j to com-
partment n.
The last row of controllability matrix P is:
[B nl 1I... [B] nrb I[p4n i I -* [A Inrb [A2 Inl 1 |. [A2] 1nb * 1 [An rb Inl 1 .. I [A n-b]Inrb.- (20) From (16) and (17) it follows that all the entries of the last Q.E.D. of P are zero, therefore its rank is less than n. With reference to the complete observability, let the observed compartments be numbered from 1 to rc. Let us suppose there is no path from compartment n to any observable compartment. It can be seen, by inspection, that the last row of observability matrix Q is null, therefore its rank is less Q.E.D. than n. 7Theorem 2.1 Proof: It is necessary to prove that, under such hypotheses, matrix P cannot have a rank less than n for any combination of the parameter values {kj,} or equally that det [PPTI is not identically zero. It will be proved that [PPTI has neither its rows or columns identically null nor linearly dependent. Let the controlled compartments be numbered from I to rb. The i-th row (i > rb) of matrix P is expressed by: row
0 1 [A I il l
O *.
I
[A n-rb]il
I [A I irb I [A' I il l *l[X2 ] irbl | [An-b]irb
(21)
while for i < rb expression (21) changes only in that [PI ii = 1. Note by inspection that at least one entry of every row i > rb is non-null as there is at least one path not exceeding the length n - rb from one of the controlled compartments to every uncontrolled one. Obviously the same holds for every row i < rb The general entry on the main diagonal of [PPTI is expressed by:
[PP'l ii
n-rb
=
[PPTI ii =
E,
1=1
n-rb
E
rb
E
m =1
[A' ]m + 1,
i1
rb
i>
rb
rb
[A']mIM X
E
1=1 M=l
It is evident that the entries
on
.
the main diagonal cannot be
zero, therefore matrix ppT has neither null rows nor columns. Moreover, there is no linear dependence among the rows as every entry of a row is a sum of monomials one of which does not appear in any other entry of the same column, except for the case in which two or more rows have only one non zero element on the same column. In this last case the corresponding compartments are closed and only influenced by the same compartment; the state variables cannot vary independently, but only proportionally: however, that originates no trouble for identification, and therefore the system will be considered Q.E.D. CC all the same. Theorem 2.2 Proof: Let the observed compartments be numbered from 1 to r,. Following the same procedure used in Theorem 2.1, it can be proved that det [QQT] cannot be identically zero. ACKNOWLEDGMENT The authors wish to thank Prof. A. Lepschy, Istituto di Elettrotecnica e di Elettronica, Facolta di Ingegneria della Universita di Padova for his interest and his helpful suggestions. Thanks are also due to Prof. A. Rescigno, Department of Physiology, University of Minneapolis, Minnesota, for the careful reading of the manuscript.
REFERENCES [1] R. Bellman and K. J. Astrom, "On structural identifiability," Math. Biosci., vol. 7, pp. 329-339, 1970. 121 M. Hajek, "A contribution to the parameters estimation of a certain class of dynamical systems," Kybernetika, vol. 8, pp. 165-
172, 1972. [31 C. Cobelli and G. Romanin-Jacur, "Structural identifiability of strongly connected biological compartmental systems," Medical & Biological Engineering, vol. 13, pp. 831-838, 1975. [4] A. Rescigno and G. Segre, Drug and Tracer Kinetics. Waltham: Blaisdell Publ. Co., 1966. [5 1 J. A. Jacquez, Compartmental Analysis in Biology and Medicine. Amsterdam: Elsevier Publ. Co., 1972. [61 J. Z. Hearon, "Theorems on linear systems," Ann. N. Y. Acad. Sci., vol. 108, pp. 3648, 1963. [71 L. A. Zadeh and C. A. Desoer, Linear System Theory. New York: McGraw-Hill, 1963. [81 R. E. Kalman, "Mathematical description of linear dynamical system," J. SIAM. Control, vol. 1, pp. 152-192, 1963. [9] C. T. Chen, Introduction to Linear System Theory. New York: Holt, Rinehart and Winston, 1970. [10] F. R. Gantmacher, The Theory of Matrices. New York: Chelsea Publ. Co., 1959, vol. 2, pp. 97-98. [11] B. C. McInnis, S. El-Asfouri and T. E. Bynum, "Modeling and computer simulation to determine digoxin dosage," Proceedings of the Regulation and Control in Physiological Systems Conference, Rochester, N.Y., August 22-24, 1973, Eds. A. S. Iberall and A. C. Guyton, Pittsburgh: I.S.A., 1973, pp. 69-70. [12] T. Freeman and C. M. E. Matthews, "Analysis of the behaviour of 131I-albumin in the normal subject and nephrotic patient," in Radioactive Isotope in Klinik und Forschung, Eds. K. von Fellinger and H. Vetter, Munich: Urban & Schwarzenberg, 1958, vol. 3, p. 283. [13] J. Rotblat and R. Marcus, "Clinical uses of radioiodine," in Radioisotope Techniques, Proceedings of the Isotope Techniques Conference, Oxford, July 1951, H.M.S.O., 1953, vol. 1, p. 33. [141 G. L. Brownell, "Analysis of techniques for determination of thyroid function with radioiodine," J. Clin. Endocrin., vol. 11, p. 1095, 1951. [151 B. Bloomstedt and L. 0. Plantin, "The extrathyroid distribution of 13 1 I-thyroxine," Acta Endocr., Copnh., vol. 48, p. 5 36, 1965. [161 S. H. Cohn, S. R. Bozzo, J. Jesseph, C. Constantinides, D. R. Huene and E. A. Gusmano, "Formulation and testing of a compartmental model for calcium metabolism in man," Radiat. Res., vol. 26, p. 319, 1965.
