- Email: [email protected]
Transfer function matrix of a compartmental model
Computer Programs in Biomedicine 12 (1980) 141-155
141
Elsevier/North-HollandBiomedicalPress
TRANSFER FUNCTION MATRIX OF A COMPARTMENTAL MODEL A p r o g r a m for the rapid writing o f its s y m b o l i c e x p r e s s i o n Anna.lisa BOSSI, Claudio COBELLI *, Livio COLUSSI and Giorgio ROMANIN-JACUR * Computing Center, University of Padova and * LADSEB-CNR, Padova, Italy A method for the fast writing of the symbolic expression of the transfer function matrix of a compartmental rrfodel has been implemented on a FORTRAN G program. Special care has been devoted to the reduction of computations and memory by means of original subroutines. The program input is very simple and the output very clear to general users. The program may be employed for the study of all input-output relations of the considered compartmental model. Details on the computational methods adopted for its implementation are given. Application examples are reported. A priori identifiability
Compartmental model
Transfer function matrix
1. Introduction
Linear time invariant compartmental models are widely used for quantitative analysis of biological systems, especially in conjunction with radioactive tracer experiments [1]. The transfer function matrix, by providing a synthetic expression of all relations among model inputs and outputs, is usually considered as one of the most powerful analytical tools for the study of the dynamics of compartmental models [2]. As an example, an important property which can be studied on the basis of transfer function matrix is a priori or structural identifiability [3], a necessary condition for model parameter estimation from the experimental data of a given input-output experiment. The writing of the transfer function matrix is generally a difficult task, as it normally requires the inversion of a matrix in symbol form. If high-level computer languages, specifically built up for symbolic manipulation (e.g., LISP) are employed, a very large amount of computer memory and computation time is needed even for small models. Moreover, not all computers may work in LISP or similar languages. Therefore the usefulness of a fast method, to be
implemented on any computer or mini-computer, is rather obvious. In the present paper we give a computer program, written in FORTRAN language, which writes the symbolic expression of the transfer function matrix of a linear time invariant compartmental model, without resorting to matrix inversion, and requires no special features to be run, even for models of large size. In section 2 some theoretical fundamentals are reviewed; in section 3 computational methods adopted for the program implementation are presented, while a full program description is given in section 4; some application examples are reported in section 5 and hardware and software specifications are given in section 6; some concluding remarks are reported in section 7. 2. Basic theory 2.1. Generalities on compartmental dynamics
Compartmental system dynamics is described by mass balance equations [1]. In particular an ideal or first order input-output perturbation experirnent on a constant steady state (e.g., tracer experiments)leads to the linear time invariant equations [1]: nB
Address correspondence and reprint requests to: Dr. C. Cobelli, LADSEB-CNR, Corso Staff Uniti 4, 35100 Padova, Italy.
n
JCi= ~ kijxj + ~ bilUl- ~ kjixi j=l l=1 1=0 jd:i j~i
0010-468x/80/0000-0000/$02.25 © Elsevier/North-HollandBiomedicalPress
i = l,n
(1)
142 nc
Ym = ~
i=1
m = l, n c
CmiXi
where xi is the deviation of the i-th compartmental quantity from its steady state value, ut is the l-th perturbation input (of which only fraction bil is allowed to enter compartment i), kij is the transport rate constant from compartment i into i, Ym is the m-th output (accessible port for measurement) and riB, nc are the numbers of inputs and outputs, respectively. Equivalent to eq. (1), a compartmental model can be described by a diagram, e.g, see fig. 1 for a three compartment model of bilirubin metabolism in man. Equations (1) correspond to a particular parametrization of the input-output state equations of a linear time-invariant dynamical system (see [2] for details): ~= Ax+Bu
y
= Cx
(2)
The information provided by an input-output experiment may be conveniently represented through the n c X n B transfer function matrix: G(s) = C(sI - A)-I
B-
1
d e t ( s l - A)
C(sI - A)
B
(3)
i.e., the generic (m,/) entry of G(s) reflects an experiment performed on the system between the l-th input and the m-th output (where I = 1,2 ..... n B and m = 1,2 ..... nc). It is to be noted that the computation of the symbolic expression of G(s) for a generic linear system is usually a hard task, as it requires the inversion of a matrix of dimension equal to the order of the system (i.e., the number of compartments). 2.2. Some useful properties Two properties of compartmental models, input and output connectability will be considered here. We can define them as follows: With reference to the compartmental graph (e.g., fig. 1) a compartment (and its associated state variable) is input (output) connectable if there exists a path from an input to the considered compartment (from the considered compartment to an output). A compartmental model is input (output) connectable if all its compartments are input (output) connectable. If a model is not both
kl 3
k21 -
k3~
~
k02
k12
Fig. 1. Compartmental diagram of the model of bilirubin metabolism in man; Q) = blood; (~) = liver; (~) = tissues.
Tracer test input and measurement take place in compartment (~.
input and output connectable, then the non-input connectable compartments and the non-output connectable ones may be neglected, as they do not influence the input-output relation (3); in fact, as shown in [2], the transfer function matrix depends only on the input and output connectable compartments. In the same way, for every single experiment between input l and output m considered separately, it is possible to determine the corresponding (m, l) submodel, formed by all compartments which are both input. connectable with respect to the l-th input only and output connectable with respect to the m-th output only. 2.3. The transfer function matrix o f an input-output connectable (IOC) compartmental model Besides the compartmental graph, the flow graph corresponding to description of type (2) may be considered (fig. 2). Taking this flow graph into consideration, Mason's formula [4] may be employed to write the transfer function matrix (3). Because of the particular expression of entries ai/for a compartmental model (as pointed out also in fig. 2) the transfer function matrix may be written in a comparatively simple way (see [4] for all theoretical details). Given a single input-single output, IOC compartmental model, the corresponding transfer function is characterized as follows: (i) Both denominator and numerator are polynomials in s whose coefficients are sums of monomials multilinear in the kii's, bil'S , Crni'S with positive sign and numerical coefficient equal to one.
143 The generic coefficient ofs n-g is the sum of all products of: an output-coefficient emi, an input coefficient bit,~and a monomial of degree g in the kii's, such that every considered product contains:
k12~~k21
(1) The expression of an input-output path (output coefficient times input coefficient times all kifs corresponding to the arcs of the path); (2) A denominator monomial not containing any arcs outgoing from a compartment of the considered path.
,3,//
\,,3
Fig. 2. Signal flow graph corresponding to the compartmental diagram of fig. 1 :
2=Ax+Bu;y=Cx where ::
-k21 - k31 A=
k21
k13
-k12 - k02
k31
B =
kz2
0
C = [0
0 -k13
0
c]
(ii) The denominator is a polynomial in s of degree n (where n is the number of compartments) with positive coefficients for all powers of s; the coefficient of sn is always 1 ; the coefficient of the generic sn-g is the sum of all monomials of degree g in the kifs such that every monomial contains neither pairs of the type kde kfe (corresponding to arcs outgoing from the same compartment) nor r-tuples of the type kill2 , ki2i3 .... kiril (corresponding to closed paths or 'cycles'). (iii) The numerator is a polynomial in s of degree n d - 1 (where d is the length, i.e., the number of intercompartmental arcs, of the shortest path between input and output).
If a multiqnput-multi-output IOC compartmental model is considered, then for every (m, l) submodel the corresponding transfer function may be written as above, and constitutes the (m, l) entry of the transfer function matrix G(s). All denominators of G(s) entries are factors of the common denominator det(sl-A), as it may be realized from (3). The symbolic expression of the common denominator, which is often useful, coincides with the denominator of the • whole multi-input-multi-output IOC compartmental model. 2.4. General outline o f the proeedure to write the transfer function matrix From a conceptual point of view, the whole procedure may be summarized in the following steps. Step 1: Input Matrices K (which contains all kij's), Ko (all koj'S), B and C are read by the program. Step 2: Single denominators and numerators For every input-output pair the corresponding input-output connectable submodel is determined as explained in subsection 2.2 and the corresponding denominator and numerator are written as explained in subsection 2.3. Step 3 : Input-output connectability Input-output connectability of the whole model is checked; if the model is not IOC then those compartments which are either non-IC or non-OC and all arcs outgoing from them are no more considered
144 Step 4: Common denominator The common denominator det(sl-A) of the transfer function matrix is written as explained in subsection 2.3.
3. Computational methods
As seen in section 2, the starting point to build up denominator and numerator coefficients is the construction of a list, L1, of simple cycles, and of a list, L2, of paths of the compartmental graph. Afterwards list L1 must be enlarged by adding all pairs of parameters of the type kdekfe related to arcs outgoing from the same compartment. Many algodthms for finding simple cycles and paths are available in the literature [4,5]. As denominator coefficients of different powers of s are characterized by different degree, it is possible to put all of them in the same list D. In order to build up list D all r-tuples k i l] l , k i 2J 2 , ..., kirJr containing no elements of list L1 as factors are to be found. A rude and systematic approach to build up list D would enumerate all 2 t - 1 r-tuples (where t is the number of non-zero transport rate parameters) and check all of them. We can observe that if a generic r-tuple contains an element of list L1 (and therefore cannot be included in D) the same happens for all s-tuples (s > r) obtained from the former by adding some klj. values. In order to speed up the procedure we designed an algorithm which enumerates and checks only those r-tuples which cannot be directly eliminated as above. The algorithm we have adopted starts by choosing an arbitrary order for the kii-values: k i l ] l ... kitJt. The first r-tuple in the enumeration is kiljl, The generic i + 1-th r-tuple in the enumeration is obtained from the i-th one according to the following rules: (i) - If the i-th r-tuple has been included in D and the r-th parameter does not coincide with kitJt we add its successor in the chosen order, thus obtaining the i + 1-th s-tuple (s = r + 1);
choosen order, thus obtaining the i + l-th s-tuple (s ---r); (iii) - If the r-th parameter in the i-th r-tuple is kitJt , cancel it; - if the resulting r - 1-tuple is not empty we replace the r - 1-th parameter by its successor in the choosen order to obtain the i + 1-th s-tuple (s = r - 1); - if the resulting r - 1-tuple is empty the algorithrn stops. For the same reasons given for the denominator also the numerator monomial may be put in a single list N. We build list N as follows. For every path in L2 we form the r-tuple R of all kij-values corresponding to arcs outgoing from the compartments crossed by the path and s-tuple S of all parameters corresponding to arcs of the path. Then we look for all tuples of list D without any ki/in common with R. All such tuples are enlarged by adding s-tuple S and by associating the input and output coefficients (b/t and Cmi, respectively) related to the path under consideration. To familiarize with the writing of denominator and numerator lists, all steps described above related to the compartmental graph of fig. 1 are reported in table 1, where tuples enumerated and not included in D, respectively in N, are cancelled by a slash. For memorization of lists D and N, by taking the chosen order for the ki/-values into account, we represent every tuple as a binary sequence of fixed length t, where the ones correspond to the kij-values which are in the tuple. Each sequence is memorized in a single word as can be seen in table 1 (in second column). This memorization permits us to save space and time. Indeed we can use boolean operations (AND, OR) on words in order to manipulate tuples. The OR operation applied to a couple of words representing different tuples gives us a word representing their union. The AND operation applied to two words representing tuples gives us a word representing their intersection.
Example (ii) - If the i-th r-tuple has not been included in D and the r-th parameter does not coincide with kitlt we replace it by its successor in the
With the chosen order for the
kij-values
of table 1 :
k12, k13, k21, k31, k02 (with reference to fig. 1,2) we have:
145
Tuple
Representation
Table 1 The building up of denominator and numerator monomials of the transfer function G(s) of the model in fig. 1
k12klako2 klak21ko2
11001 01101
- Chosen order for the ki] 's: k12 k13 k21 k31 ko2
Intersection
AND
kl2klako2 Nklak21ko2 =kl3k02
AND(11001, 01101) =0100 1
List L 1 :
-
Union
OR
kl2klak02 uklak21k02 =kl2klak21k02
OR(ll001, O1101) = 11101
k12 kl3 kl2 k21
10100 01010 10001 00110
k21 k31 k02 k31
- Denominator tuples k12 k12
kl3
step (i) step (i) step (ii) step (ii) step (iii) step (ii) step (i) step (iii) step (iii) step (i) step (i) step (ii) step (iii) step (ii) step (iii) step (i) step (iii) step (i) step (iii)
grr--kl~-0~31 k12
kl3
k12
k31
k~3 k13
k21
4.1. Main program
k13
k21
The main program SISCOMP performs the following actions (fig. 3): The data taken from the compartmental graph are entered and stored via subroutine LEGGI. The simple cycles are found by a call to subroutine CICLI. The list L1 is created in compact form. For every input-output pairs three calls are performed. The first one to subroutine CONTRL in order to rule out the compartments not connected to the particular input-output pair; the second one to subroutine SDENOM in order to compute and to print the single denominator; and the third one to the subroutine SNUMER in order to compute and to print the numerator. At last the common denominator is computed and printed by a call to subroutine SDENOM.
k13 k21 k21 ka~ kax ko2
ko2
We notice also that the monomials' degree can be easily computed by counting the number of ones in the above compact memorization. 4.
Program
description
4.2. ,Subroutine LEGGI
Subroutine LEGGI reads (and writes down for a further check) the following data: N = Number of compartments
gm~-4q-1
ko2
ko2
ko2 ko2
Numerator tuples: - List L2: contains only one zero-tuple -
k21
R
- S
empty
List N : ck12
c-
k12
ccc-
k12
c-
~-s----k~
cc cccce --
/ ~ kl3 ko2 kcr
kl3 kl3
k13
ko2
k02
k31
146
READ INPUT DATA
< CICLI
I
> >
FIND SIMPLE CYCLES
I NUMERATE THE ARCS (AND ASSOCIATED PARAMETERS) AND PRINT THE NUMERATION
I COMPUTE ALL THE PAIRS OF ARCS OUTGOING FROMTHE SAME COMPARTMENT
I BUILD LIST L1 IN COMPACT FORM
I ZlNP = 1
ZOUT
-
I
= I
-1 DIVIDUALIZE ZINP-ZOUT> \~ONNECTABLE SUBMODEL
/
I ARTIAL DENOMN I ATOR
[ UMERATOR
I
i
J
ZOUT = ZOUT + 1
YES~
?
ZINP = ZINP + i
Y E S ~
PRINTTHERULEDOUT PARAMETERSIFANY
I
PUTECOMMONDENOMINATOR MONOMIALSANDPRINTCOMMON
Fig. 3. Flo wch art of the main program SISCOMP.
NINP = Number of inputs N O U T = Number of outputs K, K ~ = Binary matrices of non-zero transport rate parameters, respectively, among compartments and towards the environment. K, KO are stored in a single array A ( N × (N + 1)) INP = Binary matrix of inputs (N × NINP) OUT = Binary matrix of outputs (NOUT ×N) 4.3. Subroutine CICLI
Subroutine CICLI computes all simple cycles of the compartmental graph using Weinblatt's [5] algorithm for the simple cycles of a digraph (fig. 4). This subroutine uses directly or indirectly subroutines CODA, COPIA, CONFR, CONCAT. According to Weinblatt's suggestion we use linked memorization for cycles and paths. For example, the cycle (ala2) (a2a3) (a3a4) (a4al) is memorized as shown in fig. 5. In our implementation the maximum number of nodes and arcs has been stated equal to 60, so that they can be numbered by binary numbers of at most 6 bits; 12 more bits are sufficient to store pointers, as the size of information dedicated to dynamic storage is 4000. In such a way we can store both the compartment number and the pointer in 18 bits. To facilitate handling of this memorization we used three functions: UNP1, UNP2, IMP: UNP1 (x) unpacks the pointer from the x-th word of the dynamic storage. UNP2 (x) unpacks the compartment number from the x-th word. IMP (a, p) packs in a single word the compartment number and the pointer p. The tops of cycles (see fig. 5) and paths are stored in arrays, the top of free dynamic storage is stored in variable IPLIB. A subroutine ADDLIB is used to return to the free dynamic storage those lists which are no more necessary. Subroutine CODA computes the function tail and the subroutine COPIA makes a copy of a list; the function CONFR tests for a merge of a list into another list. Wienblatt's algorithm in the original version needs recursive calls to the subroutine CONCAT (fig. 6). To realize the algorithm in Fortran we developed an iterative version of the algorithm by means of some
147
INITIALIZE DYNAMICSTORAGE~Erl INITIALIZE: SV( :-Z :~,,,N i
1 COUNTEROF CYCLES I --
SV(IV)=O
I IPTT=IPLIB IPLIB~MEH(IPLIB) MEr,I(IPLIB)=IMP(IV,O)
I [ FNDTT=IV--I
i
I
SV(END "FD=I
L :MEM(IPLIB) MEM(IPLIB)=IMP(ENDA,IPIT) IPTT:IPLIB IPLIB=LL
MARK ARC (ENDILLNDA) ENDTI=ENDA IV=ENDA
I
RECUR=O
SV(IV)=O CODA FIND HEAD OF CYCLE
L=IPLIB IPLIB=rIEM(L) MEM[L)=IMP(IV,L)
I
<
COPlA > COPYTHE CYCLEFOUND; SET [PI:TOP OF CYCLE
[
I NICpPc;N(CN~;~=IP1
l
I
>
I
IPP(Z)=L
I CONCAT
<
-
>
[-
L=UNPI(IPTT) MEM(IPTT)=IPLIB IPLIB=IPTT IPTT=L
i L
ENDTT=UNP2(IPTT)
I
Fig. 4. Flowchart of subroutine CICLI.
NO ~
YES
148
al
[
Pl
[
a2
]
P2
I
%
I
o
sents the compact memorization of the set of submodel's kii's. It cumulates in the variable VRIM the kii's of all the submodels. So at the end VRIM will represent the input-output connectability of the whole model (see fig. 9).
ISUC The function ISUC (I) gives the first parameter in the set ARES following the parameter I. Fig. 5. Example of cycle memorisation.
COUNT
stacks to store partial results at the various levels of recursive calls.
The function COUNT gives the number of ones in a computer word.
4.3. Subroutines for the construction of the denominator and the numerator polynomials
4.4. The output subroutines: SCRIVI, SCRIVD, SCRIVN
SDENOM Subroutine SDENOM generates the r-tuples of ki/values in the chosen order as shown in section 3.1 and builds the. fist D of denominator monomials. At the same time it computes and stores the monomial degrees. At last it prints the denominator by means of a call to subroutine SCRIVD (see fig. 7).
SNUMER The subroutine SNUMER computes the list N of numerator monomials and prints them by means of a call to the subroutine SCRIVN (see fig. 8).
We took particular care in printing output results. The two output subroutines SCRIVD and SCRIVN, starting from degree g = 1 up to the maximal monomials degree, look at the list D or N for all the monomials of degree g, and for every of them call the subroutine SCRIVI. The subroutine SCRIVI converts the monomial from the internal representation to the external representation as alphanumeric string by adding suitable parenthesis and operator symbols, and packs this string in the output line. The output line is successively written, and emptied once full.
5. Application examples
5.1. Example 1 CONTRL The subroutine CONTRL determines for an i n p u t output pair the related input-output connectable submodel and computes variable ARES which repre-
A simple example, related to the model considered in section 2, fig. 1, is first considered, in order to show how the program works and also how the symbolical expression of G(s) can be usefully employed. As seen from the computer results in fig. 10, the transfer function (1 × 1 matrix) is given by:
Cl 1[s2 + (k02 + k,2 + k13)s + (ko2k13 + kx2kx3)] G(S) = s3 +(k21 + k31 + k02 +/612 + kl3)S 2 + (k2xk02 + k21k13 + k31k02 + k31k12 + ko2k13 + k12k13)s +(k2xko2k13) (4)
149
ICH=I LEVEL OF RECURSIVE CALL
NNCT(1)=O NNCP(1)=O --1 -I
LL=IPP(ICH) LL=ToP OF PATH TO BE EXAMIN~AT THIS LEVEL
I
I
I:NNCP(ICH) I--
I=I+1
J
NO I
ADDLIB
:;NCT~~ :~:EPLAI~HEXAMINED/ <
I
ADDLIB
YES
~
No
CANCEL THE CYCLE-TAILS STORED AT LEVEL
I
ICH=ICH-Z
ICH
]
J
~
BUILD UP THE TAIL OF THE I-TH CYCLE
YES
NNCT(ICH)= NNCT(ICH)+ I
I
~
ADD THIS LAST TAIL IN THE ICH-TH LIST OF CYCLE-TAILS
ENDCP=LASTNODEOF [-TN CYCLE
Fig. 6 (continued on page 150).
ES
]
150
I RECUR=I
I NNCP(ICH)=I
I I I
ICH=ICH+I
Y E S ~
BO
[
BUILD UP NEW CYCLE
I YES
]
I
NO ALREABEEN DYFOUND?
NNCT(ICH)=O
I NNCP(ICH)=O
I
I
l ~ L ~ p ~ T ~ W o ~ A ~ PATH
IPP(ICH)=KPI
ADD NEH CYCLE TO THE CYCLE-LIST
CANCEL THE CYCLE FOUND
[
Fig. 6. Flowchart of subroutineCONCAT.
as b i x = 1 due to physical considerations. Once the considered experiment (input and output in compartment 1) has been performed, then the transfer function can be written from the analysis of input and.output behaviour: ~'(s: + (hs + # o ) G(s) = s3 + a2s2 + a l s + ao
(5)
ferent solutions in kla and (k12 + ko2): k13 =61,2 - ~ 1
+~1
4/30)
(7.1)
k12 + k02 = 62,1 = l(fl 1 -T-V~I - 4/30) By difference we have from eq. (6.4) and (6.2): k21 + k31 = or2 -/31
where t~o,t~l, o~2,#o,/31, 3' are real numbers. By equating (4) and (5) a system of six non-linear equations in the ki/'s and Cxx can be written: cll =7
(6.1)
ko2 + k12 + k13 = fll
(6.2)
ko2k13 + k12k13 =/30
(6.3)
k2a + k31 + ko2 + k12 + k l a = or2
(6.4)
k21k02 + k2akla + k3lk02 + kaak12 + ko2k13 + kx2kxa = oq
(6.5)
k21ko2kaa = ao
(6.6)
We can note that cl a is immediately given by eq. (6.1); eq. (6.2) and (6.3) permit us to obtain two dif-
(7.2)
(8)
From eq. (7.1) and (6.6):
k21ko2 = ao/61,2
(9)
By difference from eq. (6.5), (6.3) and (9): o~o k21ka3 + k~l(ko2 +k12) =o~1 - / 3 o - - ~1,2
(10)
By substituting eq. (7.1) and (7.2) in (10): 61,2k21 + t~2,1k31 =°~1 - / ~ 0 -
oto
= ~1,2
(11)
From the system of equations made up of (8) and (11) both k21 and k31 can be obtained. Finally, ko2 is obtained from eq. (9) and k12 from (7.2) by substitution. In conclusion it is possible to know, from the analysis of the symbolic expression of G(s), that two solu-
151
I z=°
I
I
I ~ol I
I
I
[
J : ISUC(O)
SET THE J-TH BIT IN PILA(1)
I IND(1)=J -i
I I
z=z+l
I
I
I
DENOM(Z,L)= PILA(1)
I
I
DENOM(Z,2)=COUNT(PILA(I ))I
I
I:I+Z
I
J = ISUC(IND(I-1))
I I
I I l, YES
~ ~
I SET IN PILA(1):PILA(I-I) ADDING THE J-TH BIT
I IND(1)=J
I
I= I - 1
Jo
I
YF$
<~ SCRIVID - - I NO
ISUC(IND(1))
YES
SET IN PILA(1) THE J-TH BIT AND THE I N ~ ( D - T H BIT
I [
Fig. 7. F l o w c h a r t
of subroutine
SDENOM.
I
IND(1)=J
152
0
COMPUTE ALL THE INPUT-OUTPUT PATHS AND STORE THEM IN THE COLUMNS OF THE MATRIX CAMrl MOREOVER STORE IN THE TWO VEC-I TORS IBB, ICC THE RELATED INPUT-OUTPUT PARAMETER TO BE I,mFD IN THE SUBROUTINE SCRIVINI
IF I zs AN INPUT COMPARTMENT
I
THEN SET MARK:(1) = i
I
ELSE SET MARK:(1) = 0
I
I
IF I IS AN OUTPUT COMPARTMENT THEN SET MARK2(1) = 1 ELSE SET MARK2(1) ~ 0
F
NUM= 0
"-I
I
IND =
I=I
"-I
FOR EVERY PAIR
SET CAM WITH THE COMPACT REPRESENTATION I-TH PATH
0
[
I (I,J)
OF COMPART-
MENTS, IF ( ! , J ) Is AN ARC THEN: I - IF MARKI{I)=I AND MARKI(J)=0 I THEN SET MARKI{J)=I AND IND=I [ - IF flARK2(J)=I AND MARK2{I)=0 I THEN SET HARK2(1)=: AND IND=I
OF THE
SET SUM WITH THE COMPACT REPRESENTATION OF THE ARCS OUTGOING FROMTHE COMPARTMENTS ON THE PATH
YES~
1 NUM = NUM + 1
I NUMER(NUM,I) = CAM NUMER(NUM,2)=COUNT(CAM)
FOR EVERY I SET I'IARKI(I) = 'ARK1( ) + NARK2(1)
I
r
]
FOR EVERY DENOMINATOR MONOMIAL DEN NOT I SHARING PARAMETERS WITH SUM SETI
NUM=NUM+I NUMER(NUM,I)=0R{CAM,DEN) NUMER(NUM,2)=COUNT(NUMER(NUM,I)
[
I=I+l
I
II
I
I A~ESoo I
FOR EVERY ARC (I,J) WITH ASSOCIAITE NUMBER LL IF MARKI(1) = 2 THENSET ARES: ARES U LL
]
[
I VRIM= OR (ARES, VRIM)
I
[ I
Fig. 9. Flowchart o f subroutine CONTRL. SCRIVIN
Fig. 8. Flowchart o f subroutine SNUMER.
tions for the k/i-values and one solution for cl 1 are obtainable from the experiment. Therefore the considered model is system identifiable, but not parameter-identifiable [3]. 5.2. Example 2 The second application example refers to a five compartment model proposed for the study of
thyroid hormone distribution and metabolism [6], reported in fig. 11, where: O = Amount of oral dose in undissolved form; @ = Amount of oral dose in dissolved form; ® = Hormone in plasma; (~) = Hormone in tissues with fast exchange; ® = Hormone in tissues with slow exchange. Two injections, the first (oral) in O and the second (intravenous) in Q are considered, while only (~) is accessible for measurement. The performing of the considered experiments will provide the ratio between the amount of orally injected hormone and the amount of absorbed hormone, which may be
153 MATHIX
1
1
0
0
0
0
I
0
MATPI.~
,~AT~IA
K0
k
MATRIX C
0
(}
~FRCb NUMEhAT I OI'J I 2 3 4
5
I~'JPUT I
( 2, ( 3', ( 0, ( I, ( l,
i) I) 2} 2) 3)
OuTwUI
directly obtained from the knowledge of transport rate parameters k21, k32, ko2 [6]. As may be easily seen from fig. 12, if G(s) were obtained according to expression (3), it would require the inversion of a symbolic matrix of order 5, containing 12 parameters (ki/'s, ci's and bil'S ). To perform such an inversion the symbolical computation of 25 determinants of 4 × 4 matrices are needed and this requires obviously a considerable amount of work even by a computer. On the contrary, by employing the proposed method without mittrix inversion, the transfer function matrix (which has a very long expression, as seen in fig. 12) can be provided quickly by the computer (2 s CPU on a CDC Cyber 76). By considering that bll and b32 in fig. 12 are equal to one due to physical considerations, it is possible to rewrite some of the coefficients as follows:
]
t)ENOMINA I OW
Input 1, Output 1
I
Denominator 1 XI*X3*XI*XS÷X(*~3*X2*X4*X~*Xf*X4*X5
XI*X~*X5
k21 +k32 + ko2 +P1
(12)
k21 (k32 + ko2) + (k21 + k32 + ko2) el + P2
(13)
Numerator C(1,1)*~(1,])
c13k32k21
(14)
C(l,l)*~(1,l)~(x3÷x,÷xb)
Input 2, Output 1 C(I,1)*~(I,1)*(~3*xS÷X4*Xb) F ~ . 10. C o m p u t e r o u t p u t r e h t e d to the c o m p ~ t m e n t model o f fig. 1.
Denominator 1
e,
(15)
P2
(16)
Numerator c13
Fig. 11. A five c o m p a r t m e n t model o f thyroid h o r m o n e distrib u t i o n and metabolism.
(17)
From eq. (15), (16), (12) and (13) two different solutions for k21 and (k32 + ko2) are possible; therefore, as k32k31 is obtainable from eq. (17) and (14),
154 4
4"
4 ~ x t
A j. <
rn -< 4.
.< +
X II,
*.-~®
O
÷
x
÷ 3. ~
x
+
c.3 ,<
÷
I~XX
e~
÷
K
~<
odx~
÷
x c 0.) +
uh
~-
x
~÷
2
~ x ~ ~ X X
."ix>( ÷ ~-
*
÷
~.1~÷
212
E
I
a3 +
.<
x
÷ +
<
;
÷
.<
÷
*
.<
..<
,It x
÷
g ,-n
..3
~ d '(
N
~
°
m x ~ ~
r~
r.3 .o
*
+
r~.. x : ~
÷ • >(,
+
z O--
~
*
°
+
~
E
m
~
x
~
~
8
×
2
z O Z
z
,
i* ~r.-
÷
~
+
..~
o. m~
I.~ x +
~ 0 ~ 0 0 ~ 0 ~ 0 0
~
+~r~-
c~
÷ o3x
x x
~ x x 4 •~ .It
< r x ~ x
. ~ x x x
~.~ x
~
155 two solutions for ka2 and k02 are possible. Therefore the experiment permit to estimate uniquely c 1a but for k21, k32 and ko2, two solutions are possible. From the practical point of view, one of the two solutions was chosen on the basis of previous knowledge and physiological considerations. It may be easily seen from above (see also [6]) that the two i n p u t - o u t p u t experiments are necessary to estimate c13, k21, k32 and ko2, while all the other model parameters cannot be estimated (they were not o f interest for the particular purpose of the study, anyway). From this second example it should be evident that the symbolical expression of G(s) is useful for planning an i n p u t - o u t p u t experiment as it permit us to know in advance whether the interesting model parameters can in principle be estimated by means of a chosen experiment.
All above limits are easily removable by means of a DIMENSION statement. The number of arcs is bounded by the variable IeL (word length), which must be preset to the size o f the computer word by a DATA statement. The program requires about 23 000 words of computer memory.
7. Conclusion A computer program for the rapid writing of the symbolic expression o f the transfer function matrix o f a compartmental model has been described. It is particularly practical and may be easily employed by a user as it may be run on many computers without requirement of outstanding features. The transfer function matrix in symbolic form may be usefully employed for any i n p u t - o u t p u t study of the model, in particular to analyze its a priori (structural) identifiability.
6. Hardware and software specifications The program is currently running on CDC CYBER 76 and IBM 370/158 computers. The program is written in FORTRAN G and contains about 400 statements. The following upper limits have been set for program variables: 40 10 10 4000 200 20 300 - 1000
Compartments Inputs Outputs Words of dynamic storage Cycles Recursive calls Elements in list L1 Denominator monomials for every denominator - 1000 Numerator monomials for every numerator - 100 Paths for every i n p u t - o u t p u t pair -
-
-
R e f e r e n c e s
[1] C. CobeUi, A. Lepschy and G. Romanin Jacur, Identifiability in tracer experiments, Fed. Proc. 39 (1980) 91. [2] C. CobeUi, A. Lepschy and G. Romanin Jacur, Identifiability of compartmental system and related structural properties, Math. Biosci. 44 (1979) 1. [3] C. Cobelli and J.J. Di Stefano III, Parameter and structural identif'mbility concepts and ambiguities: a critical review and analysis, Am. J. Physiol. 239 (1980) R7. [4] A. Bossi, C. Cobelli, L. Colussi and G. Romanin Jacur, A method of writing symbolically the transfer matrix of a compartmental model, Math. Biosci. 43 (1979) 187. [5] H. Weinblatt, A new search algorithm for finding the simple cycles of a finite directed graph, J. Assoc. Comput. Mach. 19 (1972)43. [6] J.J. Di Stefano III and P.H. Mak, On model and data requirements for determining the bioavailability of oral therapeutic agents: Application to gut absorption of thyroid hormones, Am. J. Physiol. 236 (1979) R137. [7] P.D. Berk, R.B. Howe, J.R. Bloomer and N,J. Berlin, Studies of bilirubin kinetics in normal adults, J. Clin. Invest. 48 (1969) 2176.
141
Elsevier/North-HollandBiomedicalPress
TRANSFER FUNCTION MATRIX OF A COMPARTMENTAL MODEL A p r o g r a m for the rapid writing o f its s y m b o l i c e x p r e s s i o n Anna.lisa BOSSI, Claudio COBELLI *, Livio COLUSSI and Giorgio ROMANIN-JACUR * Computing Center, University of Padova and * LADSEB-CNR, Padova, Italy A method for the fast writing of the symbolic expression of the transfer function matrix of a compartmental rrfodel has been implemented on a FORTRAN G program. Special care has been devoted to the reduction of computations and memory by means of original subroutines. The program input is very simple and the output very clear to general users. The program may be employed for the study of all input-output relations of the considered compartmental model. Details on the computational methods adopted for its implementation are given. Application examples are reported. A priori identifiability
Compartmental model
Transfer function matrix
1. Introduction
Linear time invariant compartmental models are widely used for quantitative analysis of biological systems, especially in conjunction with radioactive tracer experiments [1]. The transfer function matrix, by providing a synthetic expression of all relations among model inputs and outputs, is usually considered as one of the most powerful analytical tools for the study of the dynamics of compartmental models [2]. As an example, an important property which can be studied on the basis of transfer function matrix is a priori or structural identifiability [3], a necessary condition for model parameter estimation from the experimental data of a given input-output experiment. The writing of the transfer function matrix is generally a difficult task, as it normally requires the inversion of a matrix in symbol form. If high-level computer languages, specifically built up for symbolic manipulation (e.g., LISP) are employed, a very large amount of computer memory and computation time is needed even for small models. Moreover, not all computers may work in LISP or similar languages. Therefore the usefulness of a fast method, to be
implemented on any computer or mini-computer, is rather obvious. In the present paper we give a computer program, written in FORTRAN language, which writes the symbolic expression of the transfer function matrix of a linear time invariant compartmental model, without resorting to matrix inversion, and requires no special features to be run, even for models of large size. In section 2 some theoretical fundamentals are reviewed; in section 3 computational methods adopted for the program implementation are presented, while a full program description is given in section 4; some application examples are reported in section 5 and hardware and software specifications are given in section 6; some concluding remarks are reported in section 7. 2. Basic theory 2.1. Generalities on compartmental dynamics
Compartmental system dynamics is described by mass balance equations [1]. In particular an ideal or first order input-output perturbation experirnent on a constant steady state (e.g., tracer experiments)leads to the linear time invariant equations [1]: nB
Address correspondence and reprint requests to: Dr. C. Cobelli, LADSEB-CNR, Corso Staff Uniti 4, 35100 Padova, Italy.
n
JCi= ~ kijxj + ~ bilUl- ~ kjixi j=l l=1 1=0 jd:i j~i
0010-468x/80/0000-0000/$02.25 © Elsevier/North-HollandBiomedicalPress
i = l,n
(1)
142 nc
Ym = ~
i=1
m = l, n c
CmiXi
where xi is the deviation of the i-th compartmental quantity from its steady state value, ut is the l-th perturbation input (of which only fraction bil is allowed to enter compartment i), kij is the transport rate constant from compartment i into i, Ym is the m-th output (accessible port for measurement) and riB, nc are the numbers of inputs and outputs, respectively. Equivalent to eq. (1), a compartmental model can be described by a diagram, e.g, see fig. 1 for a three compartment model of bilirubin metabolism in man. Equations (1) correspond to a particular parametrization of the input-output state equations of a linear time-invariant dynamical system (see [2] for details): ~= Ax+Bu
y
= Cx
(2)
The information provided by an input-output experiment may be conveniently represented through the n c X n B transfer function matrix: G(s) = C(sI - A)-I
B-
1
d e t ( s l - A)
C(sI - A)
B
(3)
i.e., the generic (m,/) entry of G(s) reflects an experiment performed on the system between the l-th input and the m-th output (where I = 1,2 ..... n B and m = 1,2 ..... nc). It is to be noted that the computation of the symbolic expression of G(s) for a generic linear system is usually a hard task, as it requires the inversion of a matrix of dimension equal to the order of the system (i.e., the number of compartments). 2.2. Some useful properties Two properties of compartmental models, input and output connectability will be considered here. We can define them as follows: With reference to the compartmental graph (e.g., fig. 1) a compartment (and its associated state variable) is input (output) connectable if there exists a path from an input to the considered compartment (from the considered compartment to an output). A compartmental model is input (output) connectable if all its compartments are input (output) connectable. If a model is not both
kl 3
k21 -
k3~
~
k02
k12
Fig. 1. Compartmental diagram of the model of bilirubin metabolism in man; Q) = blood; (~) = liver; (~) = tissues.
Tracer test input and measurement take place in compartment (~.
input and output connectable, then the non-input connectable compartments and the non-output connectable ones may be neglected, as they do not influence the input-output relation (3); in fact, as shown in [2], the transfer function matrix depends only on the input and output connectable compartments. In the same way, for every single experiment between input l and output m considered separately, it is possible to determine the corresponding (m, l) submodel, formed by all compartments which are both input. connectable with respect to the l-th input only and output connectable with respect to the m-th output only. 2.3. The transfer function matrix o f an input-output connectable (IOC) compartmental model Besides the compartmental graph, the flow graph corresponding to description of type (2) may be considered (fig. 2). Taking this flow graph into consideration, Mason's formula [4] may be employed to write the transfer function matrix (3). Because of the particular expression of entries ai/for a compartmental model (as pointed out also in fig. 2) the transfer function matrix may be written in a comparatively simple way (see [4] for all theoretical details). Given a single input-single output, IOC compartmental model, the corresponding transfer function is characterized as follows: (i) Both denominator and numerator are polynomials in s whose coefficients are sums of monomials multilinear in the kii's, bil'S , Crni'S with positive sign and numerical coefficient equal to one.
143 The generic coefficient ofs n-g is the sum of all products of: an output-coefficient emi, an input coefficient bit,~and a monomial of degree g in the kii's, such that every considered product contains:
k12~~k21
(1) The expression of an input-output path (output coefficient times input coefficient times all kifs corresponding to the arcs of the path); (2) A denominator monomial not containing any arcs outgoing from a compartment of the considered path.
,3,//
\,,3
Fig. 2. Signal flow graph corresponding to the compartmental diagram of fig. 1 :
2=Ax+Bu;y=Cx where ::
-k21 - k31 A=
k21
k13
-k12 - k02
k31
B =
kz2
0
C = [0
0 -k13
0
c]
(ii) The denominator is a polynomial in s of degree n (where n is the number of compartments) with positive coefficients for all powers of s; the coefficient of sn is always 1 ; the coefficient of the generic sn-g is the sum of all monomials of degree g in the kifs such that every monomial contains neither pairs of the type kde kfe (corresponding to arcs outgoing from the same compartment) nor r-tuples of the type kill2 , ki2i3 .... kiril (corresponding to closed paths or 'cycles'). (iii) The numerator is a polynomial in s of degree n d - 1 (where d is the length, i.e., the number of intercompartmental arcs, of the shortest path between input and output).
If a multiqnput-multi-output IOC compartmental model is considered, then for every (m, l) submodel the corresponding transfer function may be written as above, and constitutes the (m, l) entry of the transfer function matrix G(s). All denominators of G(s) entries are factors of the common denominator det(sl-A), as it may be realized from (3). The symbolic expression of the common denominator, which is often useful, coincides with the denominator of the • whole multi-input-multi-output IOC compartmental model. 2.4. General outline o f the proeedure to write the transfer function matrix From a conceptual point of view, the whole procedure may be summarized in the following steps. Step 1: Input Matrices K (which contains all kij's), Ko (all koj'S), B and C are read by the program. Step 2: Single denominators and numerators For every input-output pair the corresponding input-output connectable submodel is determined as explained in subsection 2.2 and the corresponding denominator and numerator are written as explained in subsection 2.3. Step 3 : Input-output connectability Input-output connectability of the whole model is checked; if the model is not IOC then those compartments which are either non-IC or non-OC and all arcs outgoing from them are no more considered
144 Step 4: Common denominator The common denominator det(sl-A) of the transfer function matrix is written as explained in subsection 2.3.
3. Computational methods
As seen in section 2, the starting point to build up denominator and numerator coefficients is the construction of a list, L1, of simple cycles, and of a list, L2, of paths of the compartmental graph. Afterwards list L1 must be enlarged by adding all pairs of parameters of the type kdekfe related to arcs outgoing from the same compartment. Many algodthms for finding simple cycles and paths are available in the literature [4,5]. As denominator coefficients of different powers of s are characterized by different degree, it is possible to put all of them in the same list D. In order to build up list D all r-tuples k i l] l , k i 2J 2 , ..., kirJr containing no elements of list L1 as factors are to be found. A rude and systematic approach to build up list D would enumerate all 2 t - 1 r-tuples (where t is the number of non-zero transport rate parameters) and check all of them. We can observe that if a generic r-tuple contains an element of list L1 (and therefore cannot be included in D) the same happens for all s-tuples (s > r) obtained from the former by adding some klj. values. In order to speed up the procedure we designed an algorithm which enumerates and checks only those r-tuples which cannot be directly eliminated as above. The algorithm we have adopted starts by choosing an arbitrary order for the kii-values: k i l ] l ... kitJt. The first r-tuple in the enumeration is kiljl, The generic i + 1-th r-tuple in the enumeration is obtained from the i-th one according to the following rules: (i) - If the i-th r-tuple has been included in D and the r-th parameter does not coincide with kitJt we add its successor in the chosen order, thus obtaining the i + 1-th s-tuple (s = r + 1);
choosen order, thus obtaining the i + l-th s-tuple (s ---r); (iii) - If the r-th parameter in the i-th r-tuple is kitJt , cancel it; - if the resulting r - 1-tuple is not empty we replace the r - 1-th parameter by its successor in the choosen order to obtain the i + 1-th s-tuple (s = r - 1); - if the resulting r - 1-tuple is empty the algorithrn stops. For the same reasons given for the denominator also the numerator monomial may be put in a single list N. We build list N as follows. For every path in L2 we form the r-tuple R of all kij-values corresponding to arcs outgoing from the compartments crossed by the path and s-tuple S of all parameters corresponding to arcs of the path. Then we look for all tuples of list D without any ki/in common with R. All such tuples are enlarged by adding s-tuple S and by associating the input and output coefficients (b/t and Cmi, respectively) related to the path under consideration. To familiarize with the writing of denominator and numerator lists, all steps described above related to the compartmental graph of fig. 1 are reported in table 1, where tuples enumerated and not included in D, respectively in N, are cancelled by a slash. For memorization of lists D and N, by taking the chosen order for the ki/-values into account, we represent every tuple as a binary sequence of fixed length t, where the ones correspond to the kij-values which are in the tuple. Each sequence is memorized in a single word as can be seen in table 1 (in second column). This memorization permits us to save space and time. Indeed we can use boolean operations (AND, OR) on words in order to manipulate tuples. The OR operation applied to a couple of words representing different tuples gives us a word representing their union. The AND operation applied to two words representing tuples gives us a word representing their intersection.
Example (ii) - If the i-th r-tuple has not been included in D and the r-th parameter does not coincide with kitlt we replace it by its successor in the
With the chosen order for the
kij-values
of table 1 :
k12, k13, k21, k31, k02 (with reference to fig. 1,2) we have:
145
Tuple
Representation
Table 1 The building up of denominator and numerator monomials of the transfer function G(s) of the model in fig. 1
k12klako2 klak21ko2
11001 01101
- Chosen order for the ki] 's: k12 k13 k21 k31 ko2
Intersection
AND
kl2klako2 Nklak21ko2 =kl3k02
AND(11001, 01101) =0100 1
List L 1 :
-
Union
OR
kl2klak02 uklak21k02 =kl2klak21k02
OR(ll001, O1101) = 11101
k12 kl3 kl2 k21
10100 01010 10001 00110
k21 k31 k02 k31
- Denominator tuples k12 k12
kl3
step (i) step (i) step (ii) step (ii) step (iii) step (ii) step (i) step (iii) step (iii) step (i) step (i) step (ii) step (iii) step (ii) step (iii) step (i) step (iii) step (i) step (iii)
grr--kl~-0~31 k12
kl3
k12
k31
k~3 k13
k21
4.1. Main program
k13
k21
The main program SISCOMP performs the following actions (fig. 3): The data taken from the compartmental graph are entered and stored via subroutine LEGGI. The simple cycles are found by a call to subroutine CICLI. The list L1 is created in compact form. For every input-output pairs three calls are performed. The first one to subroutine CONTRL in order to rule out the compartments not connected to the particular input-output pair; the second one to subroutine SDENOM in order to compute and to print the single denominator; and the third one to the subroutine SNUMER in order to compute and to print the numerator. At last the common denominator is computed and printed by a call to subroutine SDENOM.
k13 k21 k21 ka~ kax ko2
ko2
We notice also that the monomials' degree can be easily computed by counting the number of ones in the above compact memorization. 4.
Program
description
4.2. ,Subroutine LEGGI
Subroutine LEGGI reads (and writes down for a further check) the following data: N = Number of compartments
gm~-4q-1
ko2
ko2
ko2 ko2
Numerator tuples: - List L2: contains only one zero-tuple -
k21
R
- S
empty
List N : ck12
c-
k12
ccc-
k12
c-
~-s----k~
cc cccce --
/ ~ kl3 ko2 kcr
kl3 kl3
k13
ko2
k02
k31
146
READ INPUT DATA
< CICLI
I
> >
FIND SIMPLE CYCLES
I NUMERATE THE ARCS (AND ASSOCIATED PARAMETERS) AND PRINT THE NUMERATION
I COMPUTE ALL THE PAIRS OF ARCS OUTGOING FROMTHE SAME COMPARTMENT
I BUILD LIST L1 IN COMPACT FORM
I ZlNP = 1
ZOUT
-
I
= I
-1 DIVIDUALIZE ZINP-ZOUT> \~ONNECTABLE SUBMODEL
/
I ARTIAL DENOMN I ATOR
[ UMERATOR
I
i
J
ZOUT = ZOUT + 1
YES~
?
ZINP = ZINP + i
Y E S ~
PRINTTHERULEDOUT PARAMETERSIFANY
I
PUTECOMMONDENOMINATOR MONOMIALSANDPRINTCOMMON
Fig. 3. Flo wch art of the main program SISCOMP.
NINP = Number of inputs N O U T = Number of outputs K, K ~ = Binary matrices of non-zero transport rate parameters, respectively, among compartments and towards the environment. K, KO are stored in a single array A ( N × (N + 1)) INP = Binary matrix of inputs (N × NINP) OUT = Binary matrix of outputs (NOUT ×N) 4.3. Subroutine CICLI
Subroutine CICLI computes all simple cycles of the compartmental graph using Weinblatt's [5] algorithm for the simple cycles of a digraph (fig. 4). This subroutine uses directly or indirectly subroutines CODA, COPIA, CONFR, CONCAT. According to Weinblatt's suggestion we use linked memorization for cycles and paths. For example, the cycle (ala2) (a2a3) (a3a4) (a4al) is memorized as shown in fig. 5. In our implementation the maximum number of nodes and arcs has been stated equal to 60, so that they can be numbered by binary numbers of at most 6 bits; 12 more bits are sufficient to store pointers, as the size of information dedicated to dynamic storage is 4000. In such a way we can store both the compartment number and the pointer in 18 bits. To facilitate handling of this memorization we used three functions: UNP1, UNP2, IMP: UNP1 (x) unpacks the pointer from the x-th word of the dynamic storage. UNP2 (x) unpacks the compartment number from the x-th word. IMP (a, p) packs in a single word the compartment number and the pointer p. The tops of cycles (see fig. 5) and paths are stored in arrays, the top of free dynamic storage is stored in variable IPLIB. A subroutine ADDLIB is used to return to the free dynamic storage those lists which are no more necessary. Subroutine CODA computes the function tail and the subroutine COPIA makes a copy of a list; the function CONFR tests for a merge of a list into another list. Wienblatt's algorithm in the original version needs recursive calls to the subroutine CONCAT (fig. 6). To realize the algorithm in Fortran we developed an iterative version of the algorithm by means of some
147
INITIALIZE DYNAMICSTORAGE~Erl INITIALIZE: SV( :-Z :~,,,N i
1 COUNTEROF CYCLES I --
SV(IV)=O
I IPTT=IPLIB IPLIB~MEH(IPLIB) MEr,I(IPLIB)=IMP(IV,O)
I [ FNDTT=IV--I
i
I
SV(END "FD=I
L :MEM(IPLIB) MEM(IPLIB)=IMP(ENDA,IPIT) IPTT:IPLIB IPLIB=LL
MARK ARC (ENDILLNDA) ENDTI=ENDA IV=ENDA
I
RECUR=O
SV(IV)=O CODA FIND HEAD OF CYCLE
L=IPLIB IPLIB=rIEM(L) MEM[L)=IMP(IV,L)
I
<
COPlA > COPYTHE CYCLEFOUND; SET [PI:TOP OF CYCLE
[
I NICpPc;N(CN~;~=IP1
l
I
>
I
IPP(Z)=L
I CONCAT
<
-
>
[-
L=UNPI(IPTT) MEM(IPTT)=IPLIB IPLIB=IPTT IPTT=L
i L
ENDTT=UNP2(IPTT)
I
Fig. 4. Flowchart of subroutine CICLI.
NO ~
YES
148
al
[
Pl
[
a2
]
P2
I
%
I
o
sents the compact memorization of the set of submodel's kii's. It cumulates in the variable VRIM the kii's of all the submodels. So at the end VRIM will represent the input-output connectability of the whole model (see fig. 9).
ISUC The function ISUC (I) gives the first parameter in the set ARES following the parameter I. Fig. 5. Example of cycle memorisation.
COUNT
stacks to store partial results at the various levels of recursive calls.
The function COUNT gives the number of ones in a computer word.
4.3. Subroutines for the construction of the denominator and the numerator polynomials
4.4. The output subroutines: SCRIVI, SCRIVD, SCRIVN
SDENOM Subroutine SDENOM generates the r-tuples of ki/values in the chosen order as shown in section 3.1 and builds the. fist D of denominator monomials. At the same time it computes and stores the monomial degrees. At last it prints the denominator by means of a call to subroutine SCRIVD (see fig. 7).
SNUMER The subroutine SNUMER computes the list N of numerator monomials and prints them by means of a call to the subroutine SCRIVN (see fig. 8).
We took particular care in printing output results. The two output subroutines SCRIVD and SCRIVN, starting from degree g = 1 up to the maximal monomials degree, look at the list D or N for all the monomials of degree g, and for every of them call the subroutine SCRIVI. The subroutine SCRIVI converts the monomial from the internal representation to the external representation as alphanumeric string by adding suitable parenthesis and operator symbols, and packs this string in the output line. The output line is successively written, and emptied once full.
5. Application examples
5.1. Example 1 CONTRL The subroutine CONTRL determines for an i n p u t output pair the related input-output connectable submodel and computes variable ARES which repre-
A simple example, related to the model considered in section 2, fig. 1, is first considered, in order to show how the program works and also how the symbolical expression of G(s) can be usefully employed. As seen from the computer results in fig. 10, the transfer function (1 × 1 matrix) is given by:
Cl 1[s2 + (k02 + k,2 + k13)s + (ko2k13 + kx2kx3)] G(S) = s3 +(k21 + k31 + k02 +/612 + kl3)S 2 + (k2xk02 + k21k13 + k31k02 + k31k12 + ko2k13 + k12k13)s +(k2xko2k13) (4)
149
ICH=I LEVEL OF RECURSIVE CALL
NNCT(1)=O NNCP(1)=O --1 -I
LL=IPP(ICH) LL=ToP OF PATH TO BE EXAMIN~AT THIS LEVEL
I
I
I:NNCP(ICH) I--
I=I+1
J
NO I
ADDLIB
:;NCT~~ :~:EPLAI~HEXAMINED/ <
I
ADDLIB
YES
~
No
CANCEL THE CYCLE-TAILS STORED AT LEVEL
I
ICH=ICH-Z
ICH
]
J
~
BUILD UP THE TAIL OF THE I-TH CYCLE
YES
NNCT(ICH)= NNCT(ICH)+ I
I
~
ADD THIS LAST TAIL IN THE ICH-TH LIST OF CYCLE-TAILS
ENDCP=LASTNODEOF [-TN CYCLE
Fig. 6 (continued on page 150).
ES
]
150
I RECUR=I
I NNCP(ICH)=I
I I I
ICH=ICH+I
Y E S ~
BO
[
BUILD UP NEW CYCLE
I YES
]
I
NO ALREABEEN DYFOUND?
NNCT(ICH)=O
I NNCP(ICH)=O
I
I
l ~ L ~ p ~ T ~ W o ~ A ~ PATH
IPP(ICH)=KPI
ADD NEH CYCLE TO THE CYCLE-LIST
CANCEL THE CYCLE FOUND
[
Fig. 6. Flowchart of subroutineCONCAT.
as b i x = 1 due to physical considerations. Once the considered experiment (input and output in compartment 1) has been performed, then the transfer function can be written from the analysis of input and.output behaviour: ~'(s: + (hs + # o ) G(s) = s3 + a2s2 + a l s + ao
(5)
ferent solutions in kla and (k12 + ko2): k13 =61,2 - ~ 1
+~1
4/30)
(7.1)
k12 + k02 = 62,1 = l(fl 1 -T-V~I - 4/30) By difference we have from eq. (6.4) and (6.2): k21 + k31 = or2 -/31
where t~o,t~l, o~2,#o,/31, 3' are real numbers. By equating (4) and (5) a system of six non-linear equations in the ki/'s and Cxx can be written: cll =7
(6.1)
ko2 + k12 + k13 = fll
(6.2)
ko2k13 + k12k13 =/30
(6.3)
k2a + k31 + ko2 + k12 + k l a = or2
(6.4)
k21k02 + k2akla + k3lk02 + kaak12 + ko2k13 + kx2kxa = oq
(6.5)
k21ko2kaa = ao
(6.6)
We can note that cl a is immediately given by eq. (6.1); eq. (6.2) and (6.3) permit us to obtain two dif-
(7.2)
(8)
From eq. (7.1) and (6.6):
k21ko2 = ao/61,2
(9)
By difference from eq. (6.5), (6.3) and (9): o~o k21ka3 + k~l(ko2 +k12) =o~1 - / 3 o - - ~1,2
(10)
By substituting eq. (7.1) and (7.2) in (10): 61,2k21 + t~2,1k31 =°~1 - / ~ 0 -
oto
= ~1,2
(11)
From the system of equations made up of (8) and (11) both k21 and k31 can be obtained. Finally, ko2 is obtained from eq. (9) and k12 from (7.2) by substitution. In conclusion it is possible to know, from the analysis of the symbolic expression of G(s), that two solu-
151
I z=°
I
I
I ~ol I
I
I
[
J : ISUC(O)
SET THE J-TH BIT IN PILA(1)
I IND(1)=J -i
I I
z=z+l
I
I
I
DENOM(Z,L)= PILA(1)
I
I
DENOM(Z,2)=COUNT(PILA(I ))I
I
I:I+Z
I
J = ISUC(IND(I-1))
I I
I I l, YES
~ ~
I SET IN PILA(1):PILA(I-I) ADDING THE J-TH BIT
I IND(1)=J
I
I= I - 1
Jo
I
YF$
<~ SCRIVID - - I NO
ISUC(IND(1))
YES
SET IN PILA(1) THE J-TH BIT AND THE I N ~ ( D - T H BIT
I [
Fig. 7. F l o w c h a r t
of subroutine
SDENOM.
I
IND(1)=J
152
0
COMPUTE ALL THE INPUT-OUTPUT PATHS AND STORE THEM IN THE COLUMNS OF THE MATRIX CAMrl MOREOVER STORE IN THE TWO VEC-I TORS IBB, ICC THE RELATED INPUT-OUTPUT PARAMETER TO BE I,mFD IN THE SUBROUTINE SCRIVINI
IF I zs AN INPUT COMPARTMENT
I
THEN SET MARK:(1) = i
I
ELSE SET MARK:(1) = 0
I
I
IF I IS AN OUTPUT COMPARTMENT THEN SET MARK2(1) = 1 ELSE SET MARK2(1) ~ 0
F
NUM= 0
"-I
I
IND =
I=I
"-I
FOR EVERY PAIR
SET CAM WITH THE COMPACT REPRESENTATION I-TH PATH
0
[
I (I,J)
OF COMPART-
MENTS, IF ( ! , J ) Is AN ARC THEN: I - IF MARKI{I)=I AND MARKI(J)=0 I THEN SET MARKI{J)=I AND IND=I [ - IF flARK2(J)=I AND MARK2{I)=0 I THEN SET HARK2(1)=: AND IND=I
OF THE
SET SUM WITH THE COMPACT REPRESENTATION OF THE ARCS OUTGOING FROMTHE COMPARTMENTS ON THE PATH
YES~
1 NUM = NUM + 1
I NUMER(NUM,I) = CAM NUMER(NUM,2)=COUNT(CAM)
FOR EVERY I SET I'IARKI(I) = 'ARK1( ) + NARK2(1)
I
r
]
FOR EVERY DENOMINATOR MONOMIAL DEN NOT I SHARING PARAMETERS WITH SUM SETI
NUM=NUM+I NUMER(NUM,I)=0R{CAM,DEN) NUMER(NUM,2)=COUNT(NUMER(NUM,I)
[
I=I+l
I
II
I
I A~ESoo I
FOR EVERY ARC (I,J) WITH ASSOCIAITE NUMBER LL IF MARKI(1) = 2 THENSET ARES: ARES U LL
]
[
I VRIM= OR (ARES, VRIM)
I
[ I
Fig. 9. Flowchart o f subroutine CONTRL. SCRIVIN
Fig. 8. Flowchart o f subroutine SNUMER.
tions for the k/i-values and one solution for cl 1 are obtainable from the experiment. Therefore the considered model is system identifiable, but not parameter-identifiable [3]. 5.2. Example 2 The second application example refers to a five compartment model proposed for the study of
thyroid hormone distribution and metabolism [6], reported in fig. 11, where: O = Amount of oral dose in undissolved form; @ = Amount of oral dose in dissolved form; ® = Hormone in plasma; (~) = Hormone in tissues with fast exchange; ® = Hormone in tissues with slow exchange. Two injections, the first (oral) in O and the second (intravenous) in Q are considered, while only (~) is accessible for measurement. The performing of the considered experiments will provide the ratio between the amount of orally injected hormone and the amount of absorbed hormone, which may be
153 MATHIX
1
1
0
0
0
0
I
0
MATPI.~
,~AT~IA
K0
k
MATRIX C
0
(}
~FRCb NUMEhAT I OI'J I 2 3 4
5
I~'JPUT I
( 2, ( 3', ( 0, ( I, ( l,
i) I) 2} 2) 3)
OuTwUI
directly obtained from the knowledge of transport rate parameters k21, k32, ko2 [6]. As may be easily seen from fig. 12, if G(s) were obtained according to expression (3), it would require the inversion of a symbolic matrix of order 5, containing 12 parameters (ki/'s, ci's and bil'S ). To perform such an inversion the symbolical computation of 25 determinants of 4 × 4 matrices are needed and this requires obviously a considerable amount of work even by a computer. On the contrary, by employing the proposed method without mittrix inversion, the transfer function matrix (which has a very long expression, as seen in fig. 12) can be provided quickly by the computer (2 s CPU on a CDC Cyber 76). By considering that bll and b32 in fig. 12 are equal to one due to physical considerations, it is possible to rewrite some of the coefficients as follows:
]
t)ENOMINA I OW
Input 1, Output 1
I
Denominator 1 XI*X3*XI*XS÷X(*~3*X2*X4*X~*Xf*X4*X5
XI*X~*X5
k21 +k32 + ko2 +P1
(12)
k21 (k32 + ko2) + (k21 + k32 + ko2) el + P2
(13)
Numerator C(1,1)*~(1,])
c13k32k21
(14)
C(l,l)*~(1,l)~(x3÷x,÷xb)
Input 2, Output 1 C(I,1)*~(I,1)*(~3*xS÷X4*Xb) F ~ . 10. C o m p u t e r o u t p u t r e h t e d to the c o m p ~ t m e n t model o f fig. 1.
Denominator 1
e,
(15)
P2
(16)
Numerator c13
Fig. 11. A five c o m p a r t m e n t model o f thyroid h o r m o n e distrib u t i o n and metabolism.
(17)
From eq. (15), (16), (12) and (13) two different solutions for k21 and (k32 + ko2) are possible; therefore, as k32k31 is obtainable from eq. (17) and (14),
154 4
4"
4 ~ x t
A j. <
rn -< 4.
.< +
X II,
*.-~®
O
÷
x
÷ 3. ~
x
+
c.3 ,<
÷
I~XX
e~
÷
K
~<
odx~
÷
x c 0.) +
uh
~-
x
~÷
2
~ x ~ ~ X X
."ix>( ÷ ~-
*
÷
~.1~÷
212
E
I
a3 +
.<
x
÷ +
<
;
÷
.<
÷
*
.<
..<
,It x
÷
g ,-n
..3
~ d '(
N
~
°
m x ~ ~
r~
r.3 .o
*
+
r~.. x : ~
÷ • >(,
+
z O--
~
*
°
+
~
E
m
~
x
~
~
8
×
2
z O Z
z
,
i* ~r.-
÷
~
+
..~
o. m~
I.~ x +
~ 0 ~ 0 0 ~ 0 ~ 0 0
~
+~r~-
c~
÷ o3x
x x
~ x x 4 •~ .It
< r x ~ x
. ~ x x x
~.~ x
~
155 two solutions for ka2 and k02 are possible. Therefore the experiment permit to estimate uniquely c 1a but for k21, k32 and ko2, two solutions are possible. From the practical point of view, one of the two solutions was chosen on the basis of previous knowledge and physiological considerations. It may be easily seen from above (see also [6]) that the two i n p u t - o u t p u t experiments are necessary to estimate c13, k21, k32 and ko2, while all the other model parameters cannot be estimated (they were not o f interest for the particular purpose of the study, anyway). From this second example it should be evident that the symbolical expression of G(s) is useful for planning an i n p u t - o u t p u t experiment as it permit us to know in advance whether the interesting model parameters can in principle be estimated by means of a chosen experiment.
All above limits are easily removable by means of a DIMENSION statement. The number of arcs is bounded by the variable IeL (word length), which must be preset to the size o f the computer word by a DATA statement. The program requires about 23 000 words of computer memory.
7. Conclusion A computer program for the rapid writing of the symbolic expression o f the transfer function matrix o f a compartmental model has been described. It is particularly practical and may be easily employed by a user as it may be run on many computers without requirement of outstanding features. The transfer function matrix in symbolic form may be usefully employed for any i n p u t - o u t p u t study of the model, in particular to analyze its a priori (structural) identifiability.
6. Hardware and software specifications The program is currently running on CDC CYBER 76 and IBM 370/158 computers. The program is written in FORTRAN G and contains about 400 statements. The following upper limits have been set for program variables: 40 10 10 4000 200 20 300 - 1000
Compartments Inputs Outputs Words of dynamic storage Cycles Recursive calls Elements in list L1 Denominator monomials for every denominator - 1000 Numerator monomials for every numerator - 100 Paths for every i n p u t - o u t p u t pair -
-
-
R e f e r e n c e s
[1] C. CobeUi, A. Lepschy and G. Romanin Jacur, Identifiability in tracer experiments, Fed. Proc. 39 (1980) 91. [2] C. CobeUi, A. Lepschy and G. Romanin Jacur, Identifiability of compartmental system and related structural properties, Math. Biosci. 44 (1979) 1. [3] C. Cobelli and J.J. Di Stefano III, Parameter and structural identif'mbility concepts and ambiguities: a critical review and analysis, Am. J. Physiol. 239 (1980) R7. [4] A. Bossi, C. Cobelli, L. Colussi and G. Romanin Jacur, A method of writing symbolically the transfer matrix of a compartmental model, Math. Biosci. 43 (1979) 187. [5] H. Weinblatt, A new search algorithm for finding the simple cycles of a finite directed graph, J. Assoc. Comput. Mach. 19 (1972)43. [6] J.J. Di Stefano III and P.H. Mak, On model and data requirements for determining the bioavailability of oral therapeutic agents: Application to gut absorption of thyroid hormones, Am. J. Physiol. 236 (1979) R137. [7] P.D. Berk, R.B. Howe, J.R. Bloomer and N,J. Berlin, Studies of bilirubin kinetics in normal adults, J. Clin. Invest. 48 (1969) 2176.