Computer Programs in Biomedicine 7 (1977) 21-36
© Elsevier/North-Holland Biomedical Press
A COMPUTER PROGRAM FOR THE ANALYSIS OF CONTROLLABILITY, AND STRUCTURAL
IDENTIFIABILITY
OBSERVABILITY
OF BIOLOGICAL COMPARTMENTAL
SYSTEMS
C. COBELL1, A. POLO and G. ROMANIN-JACUR Laboratorio per Rieerche di Dinamica dei Sistemi e di Elettroniea Biomediea, Consiglio Nazionale delle Rieerehe, Casella Postale 1075, 35100 Padova, Italia
A computer program to check structural identifiability of biological compartmental systems, that is the a priori possibility of estimating all unknown system parameters through a multi input-multi output tracer experiment is presented. The procedure, as based only on the adopted compartmental structure and the chosen input-output experiment, is independent of the numerical values of the parameters: therefore the program can be usefully employed before parameter estimation algorithms, to assure that all the unknown parameters evidenced in the model can be estimated from the experimental data. After a short review on compartmental models, controllability observability and structural identifiability are defined and techniques to check them are provided. The digital computer implementation of the whole procedure is discussed in detail. Some typical program runs regarding the application to biological systems are given. Compartmental Systems
Controllability
Observability
1. Introduction Compartmental models and tracer experiments are widely used for a quantitative analysis of biological systems. Given a biological system, usually, compartmental analysis involves the following steps: (i) compartmentalization of the system, i.e. on the base of the system physical structure, a block diagram characterized by compartments and intercompartmental relationships is developed; (ii) mathematical modeling of the system from the block diagram: a compartmental mathematical model consists of a set of linear differential equations describing the dynamical behavior of the system, with a set of unknown parameters (transfer rate constants) evidenced; (iii) choice of the input-output tracer experiment; (iv) parameter estimation from experimental data: the experimental samples are compared with the output values of the model whose parameters are iteratively modified through an adaptive computer procedure minimizing a suitable criterion function. Steps (i) and (iv) give difficulties for different causes: in (i) simplifying hypotheses are obviously as-
Structural ldentifiability
sumed that must be biologically significant; in step (iv) problems related to non-linear estimation are present. The whole procedure (i)-(iv) may be successful only if the adopted model is structurally identifiable, that is if all the evidenced unknown parameters can be estimated through the chosen input-output tracer experiment. An a priori identifiability analysis is necessary before step (iv). In fact, if the system is not identifiable parameter estimation is meaningless, as either it is impossible to estimate some parameters (i.e. a part of the system is not accessible) or the estimation procedure is indeterminate (i.e. some degrees of freedom are present in the solution); in such cases either a simpler model (less unknown parameters) may be adopted, though still biologically significant, or a different experiment may be planned. Therefore a check on structural identifiability of biological compartmental systems seems to be particularly useful in order to avoid repetition of tracer experiments as well as unsuccessful parameter estimation procedures. Structural identifiability of compartmental systems has recently been given an increasing attention from the theoretical point of view [ 1 - 5 ] . In this paper we present a computer program to
C. Cobelli et al., Biological compartmental systems
22
check structural identifiability of multi input-multi output biological compartmental systems, where each input enters one compartment only and each output is influenced by one compartment only. Checks on controllability and observability are also presented as they are necessary conditions for the structural identifiability. Theoretical fundamentals on this matter are widely treated in [ 3 - 5 ] and will be briefly reviewed.
/1
]=1 b i / = \ 0
(5) for uncontrolled ones
C is a constant rc × n matrix whose entries may be either 1 or 0 with the following conditions: n
]=1
cij = 1, /1
2. Basic theory and computational methods
for controlled compartments
~ ci/=\ i=1 0
i = 1, rc; for observed compartments (6) for nonobserved compartments
2.1. Compartmental models Classical biological compartmental systems are represented by [6,7]:
f¢ = A x + Bu
(1)
y = Cx
(2)
where: is an n state variable vector, where n is the number of compartments, usually representing tracer quantities in the compartments; tt is a r b ~ n vector representing input tracer quantities; is a r c ~ n vector representing measured tracer Y quantities; A is a n × n constant matrix related to the system structure matrices K and K 0, and its entries are given by: 37
ai/ = ki/
i¢]
(3) n
aii=-koi-1~1"=
kfi
(4)
where ki/is the transfer rate constant from compartment j to compartment i (compartment 0 refers to the environment). Matrices K (n × n) and K 0 (1 X n) group the rate constants between compartments and from compartments to the environment, respectively. B is a constant n × r b matrix whose entries may be either 1 or 0 with the following conditions: n
i. =~ l b O = 1,
j = 1,rb;
In compartment theory a system is called open (closed) if koi ~ 0 for at least one i (koi = 0, i = 1, n). A system is strongly connected if matrixR n = znh=l K h has no null entries.
2.2. Con trollability and observability Heuristically, a compartmental system is said to be controllable if all state variables can be modified through the inputs and is said to be observable if all outputs allow us to reconstruct all state variables. Note that controllability and observability are necessary conditions for the identifiability [4, 5]. In [4,5] it is proved that a compartmental system is completely controllable (cc) if each row o f D = (Y~= 1 Kh + 1)B has at least one positive entry and is completely observable (co) if each column o f E = C(]~= 1 Kh + 1) has at least one positive entry. If the ith row o f D (]th column of E) is identically null then compartment i (]) is not controllable (observable). Strongly connected compartmental systems are always cc and co [3]. Note that controllability and observability are considered in a structural sense, that is they are independent of the values taken by the parameters.
2.3. Structural identifiability The i n p u t - o u t p u t behavior of a cc and co compartmental system is completely described by the transfer function matrix G(s). The general entry [G(s)]i/is the transfer function of the subsystem controllable through input] (] = 1, rb) and observable
C. Cobelli et al., Biological compartmental systems through output i (i = 1, rc). The knowledge of G(s) allows the estimation of a maximum number of parameters of the cc and co part of the system equal to the number of independent coefficients of the numerator and denominator polynomials of the entries of G(s). A cc and co compartmental system is said to be structurally identifiable when the chosen i n p u t output experiment allows the estimation of all unknown parameters (rate constants). For every i/subsystem the number of estimable parameters from the denominator polynomial (see for details [4,5]) is p where p is the number of compartments of the i/subsystem; it reduces to p - 1 if the subsystem contains an integrator [4, 5]. The number of parameters which can be estimated from the numerator of an ij subsystem (see [4,5] for details) is p - r where r is such that:
[CB]i/ = [CAB]i / = ... = [CAr-IB]i/= 0
and
(7)
[CArBIi/¢ 0 except in the case r = 0 where the number is p - 1. Maybe two or more subsystems have some compartments in common: ifq is the number of common compartments then the number of estimable parameters through the denominator polynomials decreases by q (q - 1 if the common compartments include an integrator). Moreover if the common part is in cascade with the remaining parts, the number of estimable parameters from the numerator polynomials decreases by q - h - 1 where q is above and h is such that:
[A]fg = [A2]fg = ...= [Ah-l]fg = 0
[Ah ]rg ¢ 0
and
(8)
where l a n d g are respectively the first and the last compartment (with respect to the input) of the common part (see [4,5] for details). The common part between subsystems i~ and ml is obviously constituted by all compartments (d) such that: [E]i d [E]m d [D]d / [D]d I > 0 .
(9)
The common part is in cascade and h a s f a s first compartment and g as last compartment if and only if:
kdskdt>O
only f o r d = l a n d
kod kwd > 0 only for d = g
23
where s, t and v, w are compartments outside the cascade common part and belonging to the i/and ml subsystem, respectively.
2.4. Ptqncipal flow-chart A general flow-chart evidencing the whole procedure for a check of structural identifiability of compartmental systems of any structure is shown in fig. 1. In section 3 the flow-chart witl be detailed. A brief comment will be given. Input data are: (i) system structure matrices K and K 0 and (ii) input and output matrices B and C. K, K 0, B and C are given in boolean form. If the system is strongly connected the program is highly simplified (see for details [3]) as: (i) the system is always cc and co therefore no check is necessary; (ii) all subsystems coincide with the whole system, therefore all denominator polynomials of [G(s)]i/are equal and (iii) there are no common cascade parts. After a check on controllability and observability the program computes the number of the estimable parameters from denominator and numerator polynomials. As seen in section 2.3 the number of estimable parameter decreases if common parts, in cascade or not, are present among subsystems. If the number of unknown parameters is less than or equal to the number of the estimable ones, then the system is identifiable. Thus, given the system structure (K and K0) and the chosen i n p u t - o u t p u t experiment (B and C), the program gives the significance of the experiment, that is whether the parameters of the system are estimable or not. Obviously if the system is not cc or co then the system is not identifiable: in this case the same consideration as above are valid. Note that modifications in the i n p u t - o u t p u t experiment must take physical realizability into account, while modification of the system structure is to be done looking at the biological phenomenon. May be the number of estimable parameters is higher than the number of unknown ones, that is redundancy may occur. Redundancy does not affect identifiability. However, it may suggest that either a simpler experiment can be performed or a more complex model can be adopted.
C Cobelli et al., Biological compartmental systems
24
COMPUTE NUMBEROF UNKNOWN PARAMETERS CHECK STRONG CONNECTION
YES@
I
CHECK LITY CONTROLLABI AND OIBSERVABiLITY
( .... ) .o
LABLE ES
PRINT ]
( ~,o, )
NO
YES
PRINT
PRINT
STOP
STOP
Fig. 1. General flow-chart.
C. Cobelli et al., Biological compartmental systems 3. Program
description
3.1. Programming techniques The program structure strictly reflects the flowchart organization, as shown in the figs. 2 to 5. The program is segmented at some level, taking advantage of the procedure oriented characteristics of the PL/1 language, thus permitting to treat the various steps of the problem as independent blocks. Segmentation is performed according to structured programming techniques. The program is written using the "TOP-DOWN" approach; the main procedure (fig. 2) has been written first, considering the internal blocks as black boxes ("stub" procedures), then the second level procedures have been written and so on. Debugging too is facilitated, as it is restricted to one procedure at a time.
25
Program can be read from the beginning to the end without the confusing "GOTO" instructions which often cause a loss of control on the program behaviour. The only way to obtain a branch is the procedure call which always returns to the invoking point. 3.2. Main procedure description The main procedure (fig. 2) controls input data and permits several input-output experiments with the same system structure in a single program run. In the initialization section the system structure is determined; then dimensions RB and RC and values of the input output matrices, MATB and MATC, are read; a control is made on the validity of input data, as values different from 0 or 1 are not accepted; then matrices MATB and MATC are printed for possible check. At this point the problem resolution begins, according to a flag indicating whether the system is strongly connected or not; the flag is set by the initialization section according to the system structure. The cycle is repeated as long as new input data are given, i.e. as long as new input-output experiments on the same compartmental structure are to be tested. 3. 3. Initialization section
I-o l RB AND RC
MATB AND MATC
(
I
) CONN CTED INDIC, ,
NO
In this section (fig. 3) all properties of the system structure are assigned. Number of compartments NN, dimensions and values of matrices MATK and MATK0 are read with a validity check on the input data before printing. From their entries the number of unknown system parameters MM 1 is computed. Finally the strong connection is checked by means of matrix MATRN and the strong connection indicator is set: if the'system is not strongly connected MATR1 is computed for further use. 3.4. Strongly connected systems section
Fig. 2. Flow-chart of the main procedure.
In this section (fig. 4) DD is the number of degrees of freedom of the model, that is the difference between the number of the unknown parameters and the number of the mutually independent transfer function matrix coefficients. DD is initially set equal to MM 1 less the number of the denominator coefficients NN (NN - 1 if the system is closed, that is MM¢ = 0). Then the number of numerator coefficients is corn-
C. Cobelli et al., Biological compartmental systems
26
START 3
IM
I
INPUT NN ATK AND MATK®
N
I
.........
I
I
....
I
t I ,,+C~o,,TC..... I J I
NN I
NUMBER O F POSITIVE ELEMENTS OF MATC8
I °°.... ~c.~..... I
I
I
........
I
I
...........
I
i
MMI= MM+ MM~
I I
__J__ COMPUTE
I
NN
MATRN= I(MATK)I
~ I=1
O
YES
I
SET STRONGLY CONNECTED INDICATOR
IMATRI=MATRN+I I
I
1
RETURN )
Fig. 3. Flow-chart of the initialization section. Fig. 4. Flow-chart of the section for strongly-connectedsystem.
.,,,Lo.,T,
1
r
I ...........
I
C. Cobelli et al., Biological compartmental systems
27
START)
YES~
NO
I PRINT:IDENTIFIABLE SYSTEM; I IDOlREDUNDANEQUATI T ONS
PRINT:NOTIDENTIFIABLE SYSTEM ;I DDDEGREEOF SFREEDOM
i (RETURN) Fig. 5. Flow-chart of subroutine PRINT. puted by an iterative procedure and therefore subtracted. First the case of the non zero entries of MATCB is considered (see section 2.3), then the case of the zero ones for which MATCALB is to be computed. At any iteration DD is diminished by the total number CC of coefficients of numerator having degree NN - LL; the computation is stopped when the number CL of positive elements of MATCALB equals BC = RB. RC, that is MATCALB has only positive elements. Finally either the number DD of degrees of freedom if DD > 0 or the number DD of redundant equations if DD ~< 0 is printed, according to the PRINT procedure reported in fig. 5. 3.5. N o t strongly connected system section
This section starts with a check on controllability and observability and goes on with four subroutines, NUMERATOR, DENOMINATOR, SIMPLIFY and PRINT, respectively. The flow-chart of this section is reported in fig. 6. PRINT is the same as considered in 3.4 and is reported in fig. 5. The check on controllability is based on matrices MATD and MATE, respectively. As shown in section
2.2 the system is completely controllable if each row of MATD has at least one positive entry (the non controllable compartments correspond to identically null rows); the system is completely observable if each column of MATE has at least one positive entry (the non observable compartments correspond to identically null columns). If the system is not both completely controllable and completely observable, the structural identifiability check cannot be performed and the computation stops. For a controllable and observable system, subroutine NUMERATOR (fig. 7a) is called for the computation of the sum of all parameters estimable from every numerator. If MATD and MATE are reduced to binary form, then [MATD]i I = 1 means that compartment i is controllable from input l and [MATE] m/ 1 means that compartment j is observable from output m. The general ml entry of MATT = MATE × MATD gives the number of compartments of the subsystem controllable from input l and observable from m. By considering subsystem ml and the relations given in (7), the number of parameters estimable from numerator ml is given by [MATT]mr - 1 if [CB]m l = 1 and by [MATT]m / - EL if [CB]m I = [CAB]m / ...
28
C.
Cobelli et al., Biological compartmental systems
MATD~
MATR1. MATB
LIST NOT CONTROLL, MATC
• MATR1
NO
L
LIST ] NOT OBSERVABLE
0 or through a path if there exists a compartment FF such that MATRN (FF, GG) ~ 0 and MATKO (1, FF) :~ 0; a compartment FF outside the subsystem is reached if MATRN (FF, GG) > 0 where FF satisfies MATD (FF, J)" MATE (I, FF) = 0. This last check is iterated on all subsystems and DD is diminished by the total number of the denominator coefficients. In subroutine SIMPLIFY all simplifications of common factors belonging to transfer functions of different subsystems are considered. In fig. 7c the subsystem search section is represented. For every i/subsystem all possibly simplifiable ml subsystems are considered (m and l are initialized at i, respectively/', not to repeat the same operation). MATS and MATR are two n-dimension boolean row matrices which are initialized with null entries for every i~ subsystem. Their use will be explained in the following. MATS (1, FF) is put equal to 1 if compartment FF is found to be common to subsystems IJ and some ML; MATR(1, FF) is put equal to 1 whenever FF is found to belong to the cascade common part between subsystems IJand some ML. Fig. 7d shows in detail that part of subroutine SIMPLIFY which operates the search of common compartments between IJ and ML subsystems. If common compartments are found then the total number of independent coefficients is consequently reduced. Moreover common compartments possibly
29
C. Cobelli et al., Biological compartmental systems
L I..............
I
MATT=M E- MATDJ
I
¥~MAT Tli >0: I P,%,: ,,-l,c,
I
I
YES~NO
..... ~ ~' [.,.,],
[ ..............
J
I
I ......... I I
~
YES
YES
I DDcMM1 I
I ....
I
YES
I .......[
I
I
I ........ I I I I
Fig. 7b. Flow-chart of subroutine DENOMINATOR.
Fig. 7a. Flow-chart of subroutine NUMERATOR. belonging to the common cascade part between IJ and ML subsystems are searched. MATF and MATG are two n-dimension boolean matrices and are initialized with null entries for every
ML subsystem. For every FF compartment with FF = 1, NN, MATG (1, FF) is put equal to 1 if it belongs either to the IJ or the ML subsystem. If it belongs to both and MATS (1, FF) = 0 (that is the simplification relative
30
C. Cobelli et al., Biological compartmental systems
I
......
I
,I
NO
J' I
,J,, =~
I
I ::::: I M:~ [
~
[
DD= ID+T
J
v~s
J=J*l
I
MATE(M,GG)-MATD(GG,L)= e ? YES
[
YES
I ..... I
NO MATR
'ES
Fig. 7d. Flow-chart of the common compartments search of subroutine SIMPLIFY. Fig. 7c. Flow-chart of the subsystem search part of subroutine SIMPLIFY. to compartment FF has not yet been done), the cur-
rent number DD of degrees of freedom is increased by one if compartment FF neither reaches the external environment, nor is an observed compartment and
C. Cobelli et al., Biological compartmental systems
31
~VES
I
COMPUTE Vh= ~.NN
I i.=,,:=::0::o, i I
~VES
I
I ..... I
VII :
,~*T,],.,=,
L
I
COMPUTEYY:MATB(YV,J)=I CC4~PUrE Y ¥ MATy'GIIYV'I)'MATF{I'YV) z I
I ......... I
I
I
I. ,N,
HH=I
I
I
|
NO
I
I'o"S'"o~,~"T~,h
I ................ I I I ....... I 1
r>e
YES
N(;
~-7 E,,,o~,.. ~,,T~,.,=,: E"*%., o"
I............. I
IGC~4PUTEXX:MATC(I,XX)=I
I
I .......... I
t 1 .............. I
Fig. 7e. Flow-chart of the c o m m o n cascade parts search of subroutine SIMPLIFY. Fig. 7f. Flow-chart of the c o m m o n cascade parts simplification of subroutine SIMPLIFY.
6
32
C Cobelliet aL, Biological compartmental systems
reaches no compartments outside the ML subsystem. After the above considered simplification regarding denominator coefficients, if MATR (1, FF) = 0, that is FF does not belong to the common cascade part between IJ and some of the ML subsystems considered before the current one, then MATF (1, FF) = 1 and MATG (1, FF) = 0: this means that FF may possibly belong to the common cascade part between IJ and current ML subsystems. In fig. 7e the computational details for the search of the common cascade parts between numerators of different subsystems IJ and ML is presented. The first compartment YY and the last compartment XX (with respect to the input) of the common cascade part are looked for (XX and YY are initialized at zero); flag IND is initially set equal to one and is turned to zero when YY is found. Obviously the search of XX and YY is not operated if MATF = 0 as there is no common cascade part. After transposition of MATG into MATG 1 to have a column matrix, MATKG1 is computed and reduced to binary form. If the general entry of the column vector MATKG1 is non zero, then one or more compartments not belonging to the possible common cascade part, but at least to one of the two subsystems, are connected to the corresponding compartment. Integer SU1 is the number of such compartments which are directly reached from one or more compartments belonging to the possible common cascade part. SU1 = 0 means that inputs J and L coincide and the common cascade part starts in the corresponding compartment. SU1 = 1 means that only one compartment YY of the common cascade part is directly reached, then the common cascade part starts in YY. If SU1 > 1, more than one compartments of the possible common cascade part are directly reached, then they cannot belong to the actual common cascade part. In this case these compartments are no more considered and the search cycle is repeated. This iterative search ends either when compartment YY is found or when MATF is null. An analogous procedure is repeated for'the search of the last compartment XX on the common cascade part. In fig. 7f the simplification relative to the common cascade parts is operated. ZZ is the number of compartments of the common cascade part. If XX = YY either only one compartment or no compartments are included: therefore
no simplification is made. If XX :~ YY then the minimum H such that [(MATA)H]xx,Yy v~ 0 is computed, and therefore the number DD of degrees of freedom is increased by ZZ - HH - 1. Finally MATR (boolean matrix) is increased by MATF, not to repeat the simplification on the same compartments.
4. Samples of typical program runs All the examples reported in [3-5] have been tested with the program described in this paper. Here three sample runs of the structural identifiability testing procedure on three more examples of biological compartmental systems are reported. 4.1. A compartmental model o f iron kinetics
The models examined here have been proposed in [8]. In fig. 8a the compartmental structure of a first model of iron kinetics is reported, where compartment 1 represents the plasma, 2 and 3 the erythropoiesis, 4 and 5 the iron storage and 6 the hemoglobin catabolic system. The tracer experiment is performed with input and output in compartment 1. In fig. 8b the corresponding inputs (K, K 0, B, C) and outputs of the computer program are shown. The system is strongly connected, therefore only the strongly connected system section of fig. 4 is employed by the program. The model results to be identifiable with one redundant equation. A more plausible model may be adopted taking into account the red blood cells pool. The new model structure is shown in fig. 9a, where compartment 7 represents the red blood cells pool. If the intercompartmental relationship between 3 and 7 is approximated by a single transfer rate constant, the identifiability test procedure proposed here can be applied. Tracer experiments may be performed with input in 1 and outputs either in 1 or 1 and 7. In fig. 9b the corresponding inputs and outputs of the computer program are reported. The model results to be identifiable in both input-output configurations. Obviously the second experiment (outputs in 1 and 7) allows more reliable estimates of the model parameters from the experimental data (topic (iv) of section 1).
C Cobelli et al., Biological compartmental systems qATRICES
=
33
MATK AND ~IATK5
0 0 0 1 0 1 0 1 0 5 0 0 0 0 0 1 0 0 0 0 0 1 5 0 0 1 0 5 0 5 0 1 9 0 0 0 1 1 0 5 5 1 0010005
0000000
Fig. 8a. A six compartment model of iron kinetics. ~RONGLY
CONNECTED
SYSTEq
MATRICES 'IATK AND ~IATK0 ~IATRICES MATB '~ND qATC
000101 i00000 015000 150010 000100 011000
1 0 0 0 0 0 0
001000 1 0 0 0 5 0 0
~RONGLY CONNECTED SYSTEM IDENTIFIABLE ~TRICES
SYSTE'I IVITH
1 REDUNDANT E Q U A T I O N S
~4ATB AND MATC
1
MATRICES
MATB AN9 qATC
0 0 O 0 5
1 0 0 0 0 0 0
100500
1000000 0000501
IDENTIFIABLE SYSTEM t'/ITH
1 REDUNDANT EQUATIONS
BqD OF PRDGRArl
Fig. 8b. Inputs and outputs of the computer program.
Fig. 9a. A seven compartment model of iron kinetics.
IDENTIFIABLE
S Y S T E M WITH
5 REDUNDANT
EQUATIONS
END OF PROGRA~I
Fig. 9b.
Inputs and outputs of the computer program.
34
C Cobelli et al., Biological compartmental systems
4.2. A compartmental model of digitoxin kinetics A multicompartmental model for the study of the metabolic conversion of digitoxin to digoxin in man has been proposed in [9]. The compartmental structure is reported in fig. 10a, where compartment 1
represents the digitoxin in the patient's body, 2 the digoxin, 3 and 4 the urinary digitoxin and digoxin, respectively. The experiment considers an input in 1
Fig. 1 la. A four compartment model of digoxin kinetics. ~ATRICES MATK AND '4ATK0 0000 1010 01~I 0910
Fig. 10a. A four compartment model of digotoxin kinetics.
1010
NOT STRONGLY CO,"INECTED SYSTEM
VATRICES ,MATK AN') 'IATK') ~9~0 1000 I0~0 ~i00
'4ATRICES MAT~ AND MATC
1100 9010 NOT STRONGLY CONNECTED SYSTE"I ,'~DT IDE:.ITIFIABLE
SYSTEq UITH
I DEGREES OF FREEDnM
t~TP,ICES MATI3 AND 'IATC '~t",TRICES M~TB AND r.~,'%TC
901~
0001
NOT IDENTIFIABLE SYSTEM tqlTH
0190 0010
I D~GREE$ OF FREEqn'I IDENTIFIABLE SYSTEM ~'IITH
2 REDUNDANT EQUATIONS
END OF PROqRA'I END Or: PROGRAK~
Fig. lOb. Inputs and outputs of the computer program.
Fig. 1 lb. Inputs and outputs of the computer program.
c. Cobelli et al., Biological compartmental systems
and outputs from 3 and 4. Obviously both outputs in 3 and 4 are necessary to make the model observable. The model is not strongly connected and therefore the not strongly connected system section of fig. 6 is employed by the program. The model results to be not identifiable with one degree of freedom (fig. 10b). If the same model structure and input-output configuration are maintained, the only way of making the model identifiable consists in looking for external relations among the transfer rate constants. In fact in [9] the constraint k21 + k01 = given constant was assumed in order to make the system identifiable. 4.3. A tentative model o f digoxin kinetics
A tentative compartmental model of digoxin kinetics, currently under investigation by our research group, is reported in fig. 1 la, where compartment 1 represents the gastrointestinal tract, 2 the portal pool, 3 the plasma and 4 the tissues. Input-output experiments may be performed with input in 1 and outputs either in 3 or in 2 and 3. In fig. 1 lb the corresponding inputs and outputs of the computer program are reported. The model is not strongly connected and resuits to be not identifiable with the first experiment (input in 1 and output in 3). If the same model structure is maintained a more complex experiment has to be adopted in order to identify the model, that is also the portal pool is to be explored. In this case the model is identifiable with two redundant equations (note that the two subsystems have all compartments in common and compartments 1 and 2 form a common cascade part).
5. Hardware and software specifications The program runs under OS/VSI operating system on a IBM 370/158 machine. The compiler is the IBM PL/1 OPTIMIZING COMPILER. Version 1 R2.2, but no features peculiar to the optimizing version of PL/1 has been used, so the PL/1 (F) compiler should run as well. Core requirements are approximately lOOK, slightly depending on matrices dimensions. All input data are punched in free format with at least one blank between two successive values. Some explanations follow for readers unfamiliar
35
with PL/1 peculiarities. Features offered by a higher level programming language such as PL/1 have been used to minimize the number of statements. For example the clause "for all i's and/'s, if [MATT/i/> 0 then [MATT]i/ = [MATT]i/ - 1" is coded with a single statement: "MATT = MATT (MATT > 0 ) ' ; rather then with something like: "DOI=ITORC; DO J = 1 TO RB; IF MATT (I, J) > 0 THEN MATT (I, J) = MATT (I, J) - 1; END; END; (As may be noted, statements are punched one statement per card, using indentation techniques so that it is easy to see where a DO group or an IF - THEN ELSE clause begins and ends.) This reduces syntax errors at the expence of increased compilation time, but does not increase run time, the code generated being almost the same. Procedure segmentation uses internal procedures rather than external ones, to minimize arguments and parameters passing. This implies that all variables, declared in the main procedure, are known all through the program, avoiding confusing name duplication possibilities. To obtain easier reading of the source listing, variables names are chosen so that all matrices names begin with "MAT", flag variables names begin with "IND", scalar variables pertinent to the problem have two-letter names, while control variables, such as "DO" groups counters, which do not concern the problem itself, but only program control, have one letter names.
6. Mode of availability The source program listing (approximately 600 printed lines) may be obtained on request at the following address: C. COBELLI, LADSEB-CNR, Casella Postale 1075, 35100 PADOVA, Italy
36
C. Cobelli et al., Biological compartmental systems
Acknowledgements The authors would like to express their appreciation to Mr. Giampaolo Toffano for his persistence and skill in the preparations of the figures.
References [1] R. Bellman and K.J. Astrom, Math. Biosci. 7 (1970) 329. [2] M. Haiek, Kybernetika 8 (1972) 165. [3] C. Cobelli and G. Romanin-Jacur, Med. Biol. Engng. 13 (1975) 831.
[4] C. Cobelli and G. Romanin-Jacur, IEEE Trans. Biomed. Engng. 23 (1976) 93. [5] C. Cobelli, A. Lepschy and G. Romanin-Jacur, Proc. 7th IFIP Conf. on Optimization Techniques, 1975, Nice, France. Lecture Notes in Computer Science (SpringerVerlag, Berlin) 40 (1976) 88. [6] A. Rescigno and G. Segre, Drug and tracer kinetics (Blaisdell, Waltham, Mass., 1966). [71 J.A. Jacquez, Compartmental analysis in biology and medicine (Elsevier, Amsterdam, 1972). [8] G. Barosi, C. Berzuini, M. Cazzola, P. Colli Franzone, S. Morandi, M. Stefanelli, C. Viganotti and S. Perugini, J. Nucl. Biol. Med., 20 (1976) 8. [9] R.W. Jeliffe, J. Buell, R. Kalaba, R. Shridar and R. Rockwell, Math. Biosci. 6 (1970) 387.
© Elsevier/North-Holland Biomedical Press
A COMPUTER PROGRAM FOR THE ANALYSIS OF CONTROLLABILITY, AND STRUCTURAL
IDENTIFIABILITY
OBSERVABILITY
OF BIOLOGICAL COMPARTMENTAL
SYSTEMS
C. COBELL1, A. POLO and G. ROMANIN-JACUR Laboratorio per Rieerche di Dinamica dei Sistemi e di Elettroniea Biomediea, Consiglio Nazionale delle Rieerehe, Casella Postale 1075, 35100 Padova, Italia
A computer program to check structural identifiability of biological compartmental systems, that is the a priori possibility of estimating all unknown system parameters through a multi input-multi output tracer experiment is presented. The procedure, as based only on the adopted compartmental structure and the chosen input-output experiment, is independent of the numerical values of the parameters: therefore the program can be usefully employed before parameter estimation algorithms, to assure that all the unknown parameters evidenced in the model can be estimated from the experimental data. After a short review on compartmental models, controllability observability and structural identifiability are defined and techniques to check them are provided. The digital computer implementation of the whole procedure is discussed in detail. Some typical program runs regarding the application to biological systems are given. Compartmental Systems
Controllability
Observability
1. Introduction Compartmental models and tracer experiments are widely used for a quantitative analysis of biological systems. Given a biological system, usually, compartmental analysis involves the following steps: (i) compartmentalization of the system, i.e. on the base of the system physical structure, a block diagram characterized by compartments and intercompartmental relationships is developed; (ii) mathematical modeling of the system from the block diagram: a compartmental mathematical model consists of a set of linear differential equations describing the dynamical behavior of the system, with a set of unknown parameters (transfer rate constants) evidenced; (iii) choice of the input-output tracer experiment; (iv) parameter estimation from experimental data: the experimental samples are compared with the output values of the model whose parameters are iteratively modified through an adaptive computer procedure minimizing a suitable criterion function. Steps (i) and (iv) give difficulties for different causes: in (i) simplifying hypotheses are obviously as-
Structural ldentifiability
sumed that must be biologically significant; in step (iv) problems related to non-linear estimation are present. The whole procedure (i)-(iv) may be successful only if the adopted model is structurally identifiable, that is if all the evidenced unknown parameters can be estimated through the chosen input-output tracer experiment. An a priori identifiability analysis is necessary before step (iv). In fact, if the system is not identifiable parameter estimation is meaningless, as either it is impossible to estimate some parameters (i.e. a part of the system is not accessible) or the estimation procedure is indeterminate (i.e. some degrees of freedom are present in the solution); in such cases either a simpler model (less unknown parameters) may be adopted, though still biologically significant, or a different experiment may be planned. Therefore a check on structural identifiability of biological compartmental systems seems to be particularly useful in order to avoid repetition of tracer experiments as well as unsuccessful parameter estimation procedures. Structural identifiability of compartmental systems has recently been given an increasing attention from the theoretical point of view [ 1 - 5 ] . In this paper we present a computer program to
C. Cobelli et al., Biological compartmental systems
22
check structural identifiability of multi input-multi output biological compartmental systems, where each input enters one compartment only and each output is influenced by one compartment only. Checks on controllability and observability are also presented as they are necessary conditions for the structural identifiability. Theoretical fundamentals on this matter are widely treated in [ 3 - 5 ] and will be briefly reviewed.
/1
]=1 b i / = \ 0
(5) for uncontrolled ones
C is a constant rc × n matrix whose entries may be either 1 or 0 with the following conditions: n
]=1
cij = 1, /1
2. Basic theory and computational methods
for controlled compartments
~ ci/=\ i=1 0
i = 1, rc; for observed compartments (6) for nonobserved compartments
2.1. Compartmental models Classical biological compartmental systems are represented by [6,7]:
f¢ = A x + Bu
(1)
y = Cx
(2)
where: is an n state variable vector, where n is the number of compartments, usually representing tracer quantities in the compartments; tt is a r b ~ n vector representing input tracer quantities; is a r c ~ n vector representing measured tracer Y quantities; A is a n × n constant matrix related to the system structure matrices K and K 0, and its entries are given by: 37
ai/ = ki/
i¢]
(3) n
aii=-koi-1~1"=
kfi
(4)
where ki/is the transfer rate constant from compartment j to compartment i (compartment 0 refers to the environment). Matrices K (n × n) and K 0 (1 X n) group the rate constants between compartments and from compartments to the environment, respectively. B is a constant n × r b matrix whose entries may be either 1 or 0 with the following conditions: n
i. =~ l b O = 1,
j = 1,rb;
In compartment theory a system is called open (closed) if koi ~ 0 for at least one i (koi = 0, i = 1, n). A system is strongly connected if matrixR n = znh=l K h has no null entries.
2.2. Con trollability and observability Heuristically, a compartmental system is said to be controllable if all state variables can be modified through the inputs and is said to be observable if all outputs allow us to reconstruct all state variables. Note that controllability and observability are necessary conditions for the identifiability [4, 5]. In [4,5] it is proved that a compartmental system is completely controllable (cc) if each row o f D = (Y~= 1 Kh + 1)B has at least one positive entry and is completely observable (co) if each column o f E = C(]~= 1 Kh + 1) has at least one positive entry. If the ith row o f D (]th column of E) is identically null then compartment i (]) is not controllable (observable). Strongly connected compartmental systems are always cc and co [3]. Note that controllability and observability are considered in a structural sense, that is they are independent of the values taken by the parameters.
2.3. Structural identifiability The i n p u t - o u t p u t behavior of a cc and co compartmental system is completely described by the transfer function matrix G(s). The general entry [G(s)]i/is the transfer function of the subsystem controllable through input] (] = 1, rb) and observable
C. Cobelli et al., Biological compartmental systems through output i (i = 1, rc). The knowledge of G(s) allows the estimation of a maximum number of parameters of the cc and co part of the system equal to the number of independent coefficients of the numerator and denominator polynomials of the entries of G(s). A cc and co compartmental system is said to be structurally identifiable when the chosen i n p u t output experiment allows the estimation of all unknown parameters (rate constants). For every i/subsystem the number of estimable parameters from the denominator polynomial (see for details [4,5]) is p where p is the number of compartments of the i/subsystem; it reduces to p - 1 if the subsystem contains an integrator [4, 5]. The number of parameters which can be estimated from the numerator of an ij subsystem (see [4,5] for details) is p - r where r is such that:
[CB]i/ = [CAB]i / = ... = [CAr-IB]i/= 0
and
(7)
[CArBIi/¢ 0 except in the case r = 0 where the number is p - 1. Maybe two or more subsystems have some compartments in common: ifq is the number of common compartments then the number of estimable parameters through the denominator polynomials decreases by q (q - 1 if the common compartments include an integrator). Moreover if the common part is in cascade with the remaining parts, the number of estimable parameters from the numerator polynomials decreases by q - h - 1 where q is above and h is such that:
[A]fg = [A2]fg = ...= [Ah-l]fg = 0
[Ah ]rg ¢ 0
and
(8)
where l a n d g are respectively the first and the last compartment (with respect to the input) of the common part (see [4,5] for details). The common part between subsystems i~ and ml is obviously constituted by all compartments (d) such that: [E]i d [E]m d [D]d / [D]d I > 0 .
(9)
The common part is in cascade and h a s f a s first compartment and g as last compartment if and only if:
kdskdt>O
only f o r d = l a n d
kod kwd > 0 only for d = g
23
where s, t and v, w are compartments outside the cascade common part and belonging to the i/and ml subsystem, respectively.
2.4. Ptqncipal flow-chart A general flow-chart evidencing the whole procedure for a check of structural identifiability of compartmental systems of any structure is shown in fig. 1. In section 3 the flow-chart witl be detailed. A brief comment will be given. Input data are: (i) system structure matrices K and K 0 and (ii) input and output matrices B and C. K, K 0, B and C are given in boolean form. If the system is strongly connected the program is highly simplified (see for details [3]) as: (i) the system is always cc and co therefore no check is necessary; (ii) all subsystems coincide with the whole system, therefore all denominator polynomials of [G(s)]i/are equal and (iii) there are no common cascade parts. After a check on controllability and observability the program computes the number of the estimable parameters from denominator and numerator polynomials. As seen in section 2.3 the number of estimable parameter decreases if common parts, in cascade or not, are present among subsystems. If the number of unknown parameters is less than or equal to the number of the estimable ones, then the system is identifiable. Thus, given the system structure (K and K0) and the chosen i n p u t - o u t p u t experiment (B and C), the program gives the significance of the experiment, that is whether the parameters of the system are estimable or not. Obviously if the system is not cc or co then the system is not identifiable: in this case the same consideration as above are valid. Note that modifications in the i n p u t - o u t p u t experiment must take physical realizability into account, while modification of the system structure is to be done looking at the biological phenomenon. May be the number of estimable parameters is higher than the number of unknown ones, that is redundancy may occur. Redundancy does not affect identifiability. However, it may suggest that either a simpler experiment can be performed or a more complex model can be adopted.
C Cobelli et al., Biological compartmental systems
24
COMPUTE NUMBEROF UNKNOWN PARAMETERS CHECK STRONG CONNECTION
YES@
I
CHECK LITY CONTROLLABI AND OIBSERVABiLITY
( .... ) .o
LABLE ES
PRINT ]
( ~,o, )
NO
YES
STOP
STOP
Fig. 1. General flow-chart.
C. Cobelli et al., Biological compartmental systems 3. Program
description
3.1. Programming techniques The program structure strictly reflects the flowchart organization, as shown in the figs. 2 to 5. The program is segmented at some level, taking advantage of the procedure oriented characteristics of the PL/1 language, thus permitting to treat the various steps of the problem as independent blocks. Segmentation is performed according to structured programming techniques. The program is written using the "TOP-DOWN" approach; the main procedure (fig. 2) has been written first, considering the internal blocks as black boxes ("stub" procedures), then the second level procedures have been written and so on. Debugging too is facilitated, as it is restricted to one procedure at a time.
25
Program can be read from the beginning to the end without the confusing "GOTO" instructions which often cause a loss of control on the program behaviour. The only way to obtain a branch is the procedure call which always returns to the invoking point. 3.2. Main procedure description The main procedure (fig. 2) controls input data and permits several input-output experiments with the same system structure in a single program run. In the initialization section the system structure is determined; then dimensions RB and RC and values of the input output matrices, MATB and MATC, are read; a control is made on the validity of input data, as values different from 0 or 1 are not accepted; then matrices MATB and MATC are printed for possible check. At this point the problem resolution begins, according to a flag indicating whether the system is strongly connected or not; the flag is set by the initialization section according to the system structure. The cycle is repeated as long as new input data are given, i.e. as long as new input-output experiments on the same compartmental structure are to be tested. 3. 3. Initialization section
I-o l RB AND RC
MATB AND MATC
(
I
) CONN CTED INDIC, ,
NO
In this section (fig. 3) all properties of the system structure are assigned. Number of compartments NN, dimensions and values of matrices MATK and MATK0 are read with a validity check on the input data before printing. From their entries the number of unknown system parameters MM 1 is computed. Finally the strong connection is checked by means of matrix MATRN and the strong connection indicator is set: if the'system is not strongly connected MATR1 is computed for further use. 3.4. Strongly connected systems section
Fig. 2. Flow-chart of the main procedure.
In this section (fig. 4) DD is the number of degrees of freedom of the model, that is the difference between the number of the unknown parameters and the number of the mutually independent transfer function matrix coefficients. DD is initially set equal to MM 1 less the number of the denominator coefficients NN (NN - 1 if the system is closed, that is MM¢ = 0). Then the number of numerator coefficients is corn-
C. Cobelli et al., Biological compartmental systems
26
START 3
IM
I
INPUT NN ATK AND MATK®
N
I
.........
I
I
....
I
t I ,,+C~o,,TC..... I J I
NN I
NUMBER O F POSITIVE ELEMENTS OF MATC8
I °°.... ~c.~..... I
I
I
........
I
I
...........
I
i
MMI= MM+ MM~
I I
__J__ COMPUTE
I
NN
MATRN= I(MATK)I
~ I=1
O
YES
I
SET STRONGLY CONNECTED INDICATOR
IMATRI=MATRN+I I
I
1
RETURN )
Fig. 3. Flow-chart of the initialization section. Fig. 4. Flow-chart of the section for strongly-connectedsystem.
.,,,Lo.,T,
1
r
I ...........
I
C. Cobelli et al., Biological compartmental systems
27
START)
YES~
NO
I PRINT:IDENTIFIABLE SYSTEM; I IDOlREDUNDANEQUATI T ONS
PRINT:NOTIDENTIFIABLE SYSTEM ;I DDDEGREEOF SFREEDOM
i (RETURN) Fig. 5. Flow-chart of subroutine PRINT. puted by an iterative procedure and therefore subtracted. First the case of the non zero entries of MATCB is considered (see section 2.3), then the case of the zero ones for which MATCALB is to be computed. At any iteration DD is diminished by the total number CC of coefficients of numerator having degree NN - LL; the computation is stopped when the number CL of positive elements of MATCALB equals BC = RB. RC, that is MATCALB has only positive elements. Finally either the number DD of degrees of freedom if DD > 0 or the number DD of redundant equations if DD ~< 0 is printed, according to the PRINT procedure reported in fig. 5. 3.5. N o t strongly connected system section
This section starts with a check on controllability and observability and goes on with four subroutines, NUMERATOR, DENOMINATOR, SIMPLIFY and PRINT, respectively. The flow-chart of this section is reported in fig. 6. PRINT is the same as considered in 3.4 and is reported in fig. 5. The check on controllability is based on matrices MATD and MATE, respectively. As shown in section
2.2 the system is completely controllable if each row of MATD has at least one positive entry (the non controllable compartments correspond to identically null rows); the system is completely observable if each column of MATE has at least one positive entry (the non observable compartments correspond to identically null columns). If the system is not both completely controllable and completely observable, the structural identifiability check cannot be performed and the computation stops. For a controllable and observable system, subroutine NUMERATOR (fig. 7a) is called for the computation of the sum of all parameters estimable from every numerator. If MATD and MATE are reduced to binary form, then [MATD]i I = 1 means that compartment i is controllable from input l and [MATE] m/ 1 means that compartment j is observable from output m. The general ml entry of MATT = MATE × MATD gives the number of compartments of the subsystem controllable from input l and observable from m. By considering subsystem ml and the relations given in (7), the number of parameters estimable from numerator ml is given by [MATT]mr - 1 if [CB]m l = 1 and by [MATT]m / - EL if [CB]m I = [CAB]m / ...
28
C.
Cobelli et al., Biological compartmental systems
MATD~
MATR1. MATB
LIST NOT CONTROLL, MATC
• MATR1
NO
L
LIST ] NOT OBSERVABLE
0 or through a path if there exists a compartment FF such that MATRN (FF, GG) ~ 0 and MATKO (1, FF) :~ 0; a compartment FF outside the subsystem is reached if MATRN (FF, GG) > 0 where FF satisfies MATD (FF, J)" MATE (I, FF) = 0. This last check is iterated on all subsystems and DD is diminished by the total number of the denominator coefficients. In subroutine SIMPLIFY all simplifications of common factors belonging to transfer functions of different subsystems are considered. In fig. 7c the subsystem search section is represented. For every i/subsystem all possibly simplifiable ml subsystems are considered (m and l are initialized at i, respectively/', not to repeat the same operation). MATS and MATR are two n-dimension boolean row matrices which are initialized with null entries for every i~ subsystem. Their use will be explained in the following. MATS (1, FF) is put equal to 1 if compartment FF is found to be common to subsystems IJ and some ML; MATR(1, FF) is put equal to 1 whenever FF is found to belong to the cascade common part between subsystems IJand some ML. Fig. 7d shows in detail that part of subroutine SIMPLIFY which operates the search of common compartments between IJ and ML subsystems. If common compartments are found then the total number of independent coefficients is consequently reduced. Moreover common compartments possibly
29
C. Cobelli et al., Biological compartmental systems
L I..............
I
MATT=M E- MATDJ
I
¥~MAT Tli >0: I P,%,: ,,-l,c,
I
I
YES~NO
..... ~ ~' [.,.,],
[ ..............
J
I
I ......... I I
~
YES
YES
I DDcMM1 I
I ....
I
YES
I .......[
I
I
I ........ I I I I
Fig. 7b. Flow-chart of subroutine DENOMINATOR.
Fig. 7a. Flow-chart of subroutine NUMERATOR. belonging to the common cascade part between IJ and ML subsystems are searched. MATF and MATG are two n-dimension boolean matrices and are initialized with null entries for every
ML subsystem. For every FF compartment with FF = 1, NN, MATG (1, FF) is put equal to 1 if it belongs either to the IJ or the ML subsystem. If it belongs to both and MATS (1, FF) = 0 (that is the simplification relative
30
C. Cobelli et al., Biological compartmental systems
I
......
I
,I
NO
J' I
,J,, =~
I
I ::::: I M:~ [
~
[
DD= ID+T
J
v~s
J=J*l
I
MATE(M,GG)-MATD(GG,L)= e ? YES
[
YES
I ..... I
NO MATR
'ES
Fig. 7d. Flow-chart of the common compartments search of subroutine SIMPLIFY. Fig. 7c. Flow-chart of the subsystem search part of subroutine SIMPLIFY. to compartment FF has not yet been done), the cur-
rent number DD of degrees of freedom is increased by one if compartment FF neither reaches the external environment, nor is an observed compartment and
C. Cobelli et al., Biological compartmental systems
31
~VES
I
COMPUTE Vh= ~.NN
I i.=,,:=::0::o, i I
~VES
I
I ..... I
VII :
,~*T,],.,=,
L
I
COMPUTEYY:MATB(YV,J)=I CC4~PUrE Y ¥ MATy'GIIYV'I)'MATF{I'YV) z I
I ......... I
I
I
I. ,N,
HH=I
I
I
|
NO
I
I'o"S'"o~,~"T~,h
I ................ I I I ....... I 1
r>e
YES
N(;
~-7 E,,,o~,.. ~,,T~,.,=,: E"*%., o"
I............. I
IGC~4PUTEXX:MATC(I,XX)=I
I
I .......... I
t 1 .............. I
Fig. 7e. Flow-chart of the c o m m o n cascade parts search of subroutine SIMPLIFY. Fig. 7f. Flow-chart of the c o m m o n cascade parts simplification of subroutine SIMPLIFY.
6
32
C Cobelliet aL, Biological compartmental systems
reaches no compartments outside the ML subsystem. After the above considered simplification regarding denominator coefficients, if MATR (1, FF) = 0, that is FF does not belong to the common cascade part between IJ and some of the ML subsystems considered before the current one, then MATF (1, FF) = 1 and MATG (1, FF) = 0: this means that FF may possibly belong to the common cascade part between IJ and current ML subsystems. In fig. 7e the computational details for the search of the common cascade parts between numerators of different subsystems IJ and ML is presented. The first compartment YY and the last compartment XX (with respect to the input) of the common cascade part are looked for (XX and YY are initialized at zero); flag IND is initially set equal to one and is turned to zero when YY is found. Obviously the search of XX and YY is not operated if MATF = 0 as there is no common cascade part. After transposition of MATG into MATG 1 to have a column matrix, MATKG1 is computed and reduced to binary form. If the general entry of the column vector MATKG1 is non zero, then one or more compartments not belonging to the possible common cascade part, but at least to one of the two subsystems, are connected to the corresponding compartment. Integer SU1 is the number of such compartments which are directly reached from one or more compartments belonging to the possible common cascade part. SU1 = 0 means that inputs J and L coincide and the common cascade part starts in the corresponding compartment. SU1 = 1 means that only one compartment YY of the common cascade part is directly reached, then the common cascade part starts in YY. If SU1 > 1, more than one compartments of the possible common cascade part are directly reached, then they cannot belong to the actual common cascade part. In this case these compartments are no more considered and the search cycle is repeated. This iterative search ends either when compartment YY is found or when MATF is null. An analogous procedure is repeated for'the search of the last compartment XX on the common cascade part. In fig. 7f the simplification relative to the common cascade parts is operated. ZZ is the number of compartments of the common cascade part. If XX = YY either only one compartment or no compartments are included: therefore
no simplification is made. If XX :~ YY then the minimum H such that [(MATA)H]xx,Yy v~ 0 is computed, and therefore the number DD of degrees of freedom is increased by ZZ - HH - 1. Finally MATR (boolean matrix) is increased by MATF, not to repeat the simplification on the same compartments.
4. Samples of typical program runs All the examples reported in [3-5] have been tested with the program described in this paper. Here three sample runs of the structural identifiability testing procedure on three more examples of biological compartmental systems are reported. 4.1. A compartmental model o f iron kinetics
The models examined here have been proposed in [8]. In fig. 8a the compartmental structure of a first model of iron kinetics is reported, where compartment 1 represents the plasma, 2 and 3 the erythropoiesis, 4 and 5 the iron storage and 6 the hemoglobin catabolic system. The tracer experiment is performed with input and output in compartment 1. In fig. 8b the corresponding inputs (K, K 0, B, C) and outputs of the computer program are shown. The system is strongly connected, therefore only the strongly connected system section of fig. 4 is employed by the program. The model results to be identifiable with one redundant equation. A more plausible model may be adopted taking into account the red blood cells pool. The new model structure is shown in fig. 9a, where compartment 7 represents the red blood cells pool. If the intercompartmental relationship between 3 and 7 is approximated by a single transfer rate constant, the identifiability test procedure proposed here can be applied. Tracer experiments may be performed with input in 1 and outputs either in 1 or 1 and 7. In fig. 9b the corresponding inputs and outputs of the computer program are reported. The model results to be identifiable in both input-output configurations. Obviously the second experiment (outputs in 1 and 7) allows more reliable estimates of the model parameters from the experimental data (topic (iv) of section 1).
C Cobelli et al., Biological compartmental systems qATRICES
=
33
MATK AND ~IATK5
0 0 0 1 0 1 0 1 0 5 0 0 0 0 0 1 0 0 0 0 0 1 5 0 0 1 0 5 0 5 0 1 9 0 0 0 1 1 0 5 5 1 0010005
0000000
Fig. 8a. A six compartment model of iron kinetics. ~RONGLY
CONNECTED
SYSTEq
MATRICES 'IATK AND ~IATK0 ~IATRICES MATB '~ND qATC
000101 i00000 015000 150010 000100 011000
1 0 0 0 0 0 0
001000 1 0 0 0 5 0 0
~RONGLY CONNECTED SYSTEM IDENTIFIABLE ~TRICES
SYSTE'I IVITH
1 REDUNDANT E Q U A T I O N S
~4ATB AND MATC
1
MATRICES
MATB AN9 qATC
0 0 O 0 5
1 0 0 0 0 0 0
100500
1000000 0000501
IDENTIFIABLE SYSTEM t'/ITH
1 REDUNDANT EQUATIONS
BqD OF PRDGRArl
Fig. 8b. Inputs and outputs of the computer program.
Fig. 9a. A seven compartment model of iron kinetics.
IDENTIFIABLE
S Y S T E M WITH
5 REDUNDANT
EQUATIONS
END OF PROGRA~I
Fig. 9b.
Inputs and outputs of the computer program.
34
C Cobelli et al., Biological compartmental systems
4.2. A compartmental model of digitoxin kinetics A multicompartmental model for the study of the metabolic conversion of digitoxin to digoxin in man has been proposed in [9]. The compartmental structure is reported in fig. 10a, where compartment 1
represents the digitoxin in the patient's body, 2 the digoxin, 3 and 4 the urinary digitoxin and digoxin, respectively. The experiment considers an input in 1
Fig. 1 la. A four compartment model of digoxin kinetics. ~ATRICES MATK AND '4ATK0 0000 1010 01~I 0910
Fig. 10a. A four compartment model of digotoxin kinetics.
1010
NOT STRONGLY CO,"INECTED SYSTEM
VATRICES ,MATK AN') 'IATK') ~9~0 1000 I0~0 ~i00
'4ATRICES MAT~ AND MATC
1100 9010 NOT STRONGLY CONNECTED SYSTE"I ,'~DT IDE:.ITIFIABLE
SYSTEq UITH
I DEGREES OF FREEDnM
t~TP,ICES MATI3 AND 'IATC '~t",TRICES M~TB AND r.~,'%TC
901~
0001
NOT IDENTIFIABLE SYSTEM tqlTH
0190 0010
I D~GREE$ OF FREEqn'I IDENTIFIABLE SYSTEM ~'IITH
2 REDUNDANT EQUATIONS
END OF PROqRA'I END Or: PROGRAK~
Fig. lOb. Inputs and outputs of the computer program.
Fig. 1 lb. Inputs and outputs of the computer program.
c. Cobelli et al., Biological compartmental systems
and outputs from 3 and 4. Obviously both outputs in 3 and 4 are necessary to make the model observable. The model is not strongly connected and therefore the not strongly connected system section of fig. 6 is employed by the program. The model results to be not identifiable with one degree of freedom (fig. 10b). If the same model structure and input-output configuration are maintained, the only way of making the model identifiable consists in looking for external relations among the transfer rate constants. In fact in [9] the constraint k21 + k01 = given constant was assumed in order to make the system identifiable. 4.3. A tentative model o f digoxin kinetics
A tentative compartmental model of digoxin kinetics, currently under investigation by our research group, is reported in fig. 1 la, where compartment 1 represents the gastrointestinal tract, 2 the portal pool, 3 the plasma and 4 the tissues. Input-output experiments may be performed with input in 1 and outputs either in 3 or in 2 and 3. In fig. 1 lb the corresponding inputs and outputs of the computer program are reported. The model is not strongly connected and resuits to be not identifiable with the first experiment (input in 1 and output in 3). If the same model structure is maintained a more complex experiment has to be adopted in order to identify the model, that is also the portal pool is to be explored. In this case the model is identifiable with two redundant equations (note that the two subsystems have all compartments in common and compartments 1 and 2 form a common cascade part).
5. Hardware and software specifications The program runs under OS/VSI operating system on a IBM 370/158 machine. The compiler is the IBM PL/1 OPTIMIZING COMPILER. Version 1 R2.2, but no features peculiar to the optimizing version of PL/1 has been used, so the PL/1 (F) compiler should run as well. Core requirements are approximately lOOK, slightly depending on matrices dimensions. All input data are punched in free format with at least one blank between two successive values. Some explanations follow for readers unfamiliar
35
with PL/1 peculiarities. Features offered by a higher level programming language such as PL/1 have been used to minimize the number of statements. For example the clause "for all i's and/'s, if [MATT/i/> 0 then [MATT]i/ = [MATT]i/ - 1" is coded with a single statement: "MATT = MATT (MATT > 0 ) ' ; rather then with something like: "DOI=ITORC; DO J = 1 TO RB; IF MATT (I, J) > 0 THEN MATT (I, J) = MATT (I, J) - 1; END; END; (As may be noted, statements are punched one statement per card, using indentation techniques so that it is easy to see where a DO group or an IF - THEN ELSE clause begins and ends.) This reduces syntax errors at the expence of increased compilation time, but does not increase run time, the code generated being almost the same. Procedure segmentation uses internal procedures rather than external ones, to minimize arguments and parameters passing. This implies that all variables, declared in the main procedure, are known all through the program, avoiding confusing name duplication possibilities. To obtain easier reading of the source listing, variables names are chosen so that all matrices names begin with "MAT", flag variables names begin with "IND", scalar variables pertinent to the problem have two-letter names, while control variables, such as "DO" groups counters, which do not concern the problem itself, but only program control, have one letter names.
6. Mode of availability The source program listing (approximately 600 printed lines) may be obtained on request at the following address: C. COBELLI, LADSEB-CNR, Casella Postale 1075, 35100 PADOVA, Italy
36
C. Cobelli et al., Biological compartmental systems
Acknowledgements The authors would like to express their appreciation to Mr. Giampaolo Toffano for his persistence and skill in the preparations of the figures.
References [1] R. Bellman and K.J. Astrom, Math. Biosci. 7 (1970) 329. [2] M. Haiek, Kybernetika 8 (1972) 165. [3] C. Cobelli and G. Romanin-Jacur, Med. Biol. Engng. 13 (1975) 831.
[4] C. Cobelli and G. Romanin-Jacur, IEEE Trans. Biomed. Engng. 23 (1976) 93. [5] C. Cobelli, A. Lepschy and G. Romanin-Jacur, Proc. 7th IFIP Conf. on Optimization Techniques, 1975, Nice, France. Lecture Notes in Computer Science (SpringerVerlag, Berlin) 40 (1976) 88. [6] A. Rescigno and G. Segre, Drug and tracer kinetics (Blaisdell, Waltham, Mass., 1966). [71 J.A. Jacquez, Compartmental analysis in biology and medicine (Elsevier, Amsterdam, 1972). [8] G. Barosi, C. Berzuini, M. Cazzola, P. Colli Franzone, S. Morandi, M. Stefanelli, C. Viganotti and S. Perugini, J. Nucl. Biol. Med., 20 (1976) 8. [9] R.W. Jeliffe, J. Buell, R. Kalaba, R. Shridar and R. Rockwell, Math. Biosci. 6 (1970) 387.
