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

Computer Methods and Programs in Biomedicine, 37 (1992) 205-214

205

© 1992 Elsevier Science Publishers B.V. All rights reserved 0169-2607/92/$05.00 COMMET 01271 Section II. Systems and programs

Simulation of compartmental models for kinetic data from a positron emission tomograph P a m e l a G. C o x s o n a, E . M . S a l m e r o n

2 R.H.

H u e s m a n 1 a n d B.M. M a z o y e r 3

1 Center for Functional Imaging, Lawrence Berkeley Lab, Berkeley, CA 94720, USA, 2Arcueil, France and 3 Service Hospitalier Fr~ddric Joliot, Orsay, France

Linear compartmental models are used to describe the disposition of radio-labelled compounds in regions of interest in the mammalian body, based on a time sequence of measurements from a positron emission tomograph (PET). In this paper we show how closed form solutions for the model equations have been incorporated into a computer program for simulation and parameter estimation. A typical PET data example is included to illustrate the implementation and compare the closed form method with a numerical ode solution method. Compartmental models; Tracer kinetics; Positron emission tomography;Matrix exponential; Convolution

I. Introduction We have developed a package of Fortran-77 code to simulate compartmental models and fit their parameters based on kinetic data from a positron emission tomograph. By a simulation we shall mean a calculation of the perfect (noise-free) data that would result from a particular model and a given set of parameters. For linear compartmental models, our simulations are computed either by numerical integration of the differential equation or by direct evaluation of a closed form solution to the differential equation. Here we describe the latter approach. It is an alternative which requires significantly less computational effort than numerical integration.

Correspondence: Pamela G. Coxson, Staff Scientist, Center for

Functional Imaging, Mail Stop 55-121, Lawrence Berkeley Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA.

Closed form solutions for systems with two or three compartments are often calculated explicitly. They provide an efficient means to simulate low order systems. A general approach applicable to higher order systems with distinct eigenvalues and constant or impulsive inputs was reported in [1]. Our approach follows along these same lines, extending the analysis to systems with non-distinct eigenvalues and without constraints on the type of input. The removal of these restrictions makes it possible to apply the results to a broader class of models, including those which arise in kinetic analysis of data from a positron emission tomograph (PET). It is this example of PET data analysis which has motivated the work presented here. P E T input functions consist of blood concentrations of the injected radioisotope measured virtually continuously over the course of an experiment. Each PET measurement of blood or tissue concentration is actually an integral over a short interval of time, and this adds to the program-

206

ruing complexity. Mazoyer et al. [2] showed that the integral nature of the data should not be ignored. In our work, the problems of solving the differential equation and then integrating the solution over the P E T time intervals are not separated. We present a numerically tractable solution of the combined problem, which is in closed form up to numerical determination of the eigenstructure of the kinetic rate matrix. Section 2 describes the P E T kinetic model considered here, along with its closed form solution. In Section 3, an example is given of data from a typical dynamic P E T brain study. Sections 4 through 7 describe how we incorporated the closed form solution into our parameter estimation and simulation program. Section 8 discusses our implementation in the context of the example from Section 3.

2. System equations and closed form solution Compartmental systems are routinely used to model the uptake, distribution and elimination of radio-labelled compounds in particular regions of interest in the mammalian body. When the transfer rates between compartments are proportional to concentration, the system is said to have first order kinetics and the resulting differential equation is linear. The general equation of the linear model is:

The closed form solution of equation (1) is:

Q_.( t ) = foteA~IYu( t -- z )d'r + eAtQ( O) = ( e A t F ) * u ( t ) +eAtQ_.(O)

(2)

where * denotes a convolution. The first factor of the convolution, eAtl~ is the impulse response vector for the compartments; that is, the solution of equation (1) when ~ff(0)= 0 and u(t) is the impulse 3(t). A P E T scanner accumulates counts of radioactive emissions over short intervals of time (seconds). Tomographic reconstruction yields images which represent the regional distribution of radiotracer, and adding the values from a contiguous set of picture elements gives a measurement which is proportional to the amount of compound in the region of interest during the data acquisition interval. Thus the measurement function for the system is given by:

y , = [ 1 , 1. . . .

,1]ftf'+'~ff(t)dt

(3)

where [t t, tt+ 1] is the data acquisition interval, and the scalar product sums the elements of ~ff. The simulated measurement is then yt:[1,

1.....

1]ft'+'[(eAtF-~)*u(t)+ eAt~ff(O)]dl tl

d(ff(t) dt

AQ_.(t) + Fu(t)

(1)

where:

• u(t)

is the input function varying with time, • ~ff(t) is the vector containing the amount of labelled compound in each compartment of the model, • A is the system matrix of transfer rates between compartments, and • /~ is a vector of transfer rates from the input to each compartment.

(4) There are two substantial numerical hurdles involved in evaluating this closed form solution for y~: evaluation of the matrix exponential terms (eA'/Tand en'~ff(0)) and evaluation of the integral of the convolution. In Section 4, we discuss the first issue and describe the triangular and Jordan form transformations used here. In Sections 5 and 6, we show that both of these transformations allow us to e x p r e s s eAtt~ and eAt~ff(O) as Sl~(t) and S(°)E(t), respectively, where S and S (°) are n × n constant matrices and /~(t) is an n-di-

207

mensional vector of scalar functions. Section 7 addresses the integration problem.

Region of Interest

3. Example of a kinetic PET study The blood (u(t)) and P E T measurement (YJ(ti+ l -ti)) for a human brain study are plotted in Fig. 1. Each P E T data point represents the concentration of counts found to originate in the right frontal cortex of the P E T measurement slice. The patient was injected with [18F]fluorodeoxyglucose (FDG) at time 0. P E T measurements were acquired over 40 min according to the following schedule: twelve 5-s files, six 10-s files, four 30-s files, six 60-s files, and six 300-s files. During the same 40-min period, 274 blood sampies were measured. The plot shows only 20 min and a subset of the data points, so that the pattern of activity in the first few minutes can be seen. This experiment is one of many F D G studies carried out in our laboratory (see, for example, Jagust et al. [3]). A model for glucose utilization developed by Sokoloff et al, [4] was used to fit this data. The

k2

Fig. 2. Compartmental model for [lSF]fluorodeoxyglucose. A compartmental representation of the Sokoloff model with four rate constants: K I = uptake rate of F D G from the blood pool; k 2 = washout rate of non-phosphorylated FDG; k 3 = rate of phosphorylation; and k 4 = rate of dephosphorylation.

compartmental structure of the model is shown in Fig. 2. There are four model parameters: K 1 (uptake of glucose from the blood), k 2 (washout of non-phosphorylated glucose), k 3 (rate of phosphorylation by the enzyme hexokinase), and k 4 (rate of dephosphorylation). In terms of the model equations in Section 2,

k3

and

FDG Time Activity Curves Blood Levels + 10 (f'dled circles) PET Data (open circles) Fitted Model (solid line)

80~

-k 4 '

~ff(O) = 0.

0

'

(5)

An iterative procedure which invokes the simulation routine many times in each iteration is used to determine the rates {k i} that best fit the data.

7060o~

40-

8

4. The matrix exponential function

50,-o

3020100

"6'loeo•

•

•

•

•

•

•

•

•

I

I

I

I

I

4

8

12

16

20

Time (minutes) Fig. 1. M e a s u r e m e n t s of the 18 F-labelled c o m p o u n d F D G are plotted over t i m e . P E T data were acquired with the D o n n e r 600 crystal tomograph at the Lawrence Berkeley Laboratory. The 274 blood data points have been interpolated to match t h e P E T m e a s u r e m e n t times.

Matrix exponentiation, unlike the scalar exponential function, is not a built-in feature of computer languages. The subtleties involves in numerical evaluation of e TM w e r e discussed at length in [5]. A major point of that paper is that the best method for computing a matrix exponential depends strongly on the application. When multiple evaluations with different times t are required, as in this application, efficiency is increased by taking advantage of similarity transformations which make the time dependence explicit. A substantial computational effort is required to find the scalar

208 functions of time which produce each element of the matrix e A t , but that effort pays off when the formulas are used repeatedly• The matrix A is similar to H if they are related by

A = P H P -1.

(6)

Similar matrices have the same eigenvalues and the exponentials of similar matrices are similar under the same transformation P:

e At = P e m P -1

(7)

We use two transformations for which e m has a simple expression as a function of the scalar e x p o n e n t i a l s e at of the eigenvalues, At, of At. The first is a relative of the Schur decomposition in which H is upper triangular: A1

hi2

0

A2

hln

"'"

I

H= h n -- l,n

0

...

0

A~

J

(8)

where {Ai} is the set of distinct eigenvalues of A. In the second case, where some of the eigenvalues are equal, we use the Jordan Canonical Form:

The transformation P and the matrix H may be complex even if A is real. The eigenvalues are either real numbers or pairs of complex conjugates• Given the transformed matrix H and equation (7), we see that each entry of the impulse response vector, eAtl~ is a linear combination, determined by P and/Y, of the entries of e m. Thus, the impulse response can be expressed as a matrix product,

eAtlF= Sl~( t ),

w h e r e / ~ ( t ) is the vector of distinct entries of em. Similarly, the initial term has an expression,

eA'(ff( O) = S(°)E( t ).

0

0

J2

"'"

5. Solution by triangularization for systems with distinct eigenvalues

Su

where H is a block diagonal matrix, and each block Ji of H has the form hi

1.

0

0

Ai

1.

(13)

(9) 0

o

t)

0

H= ...

(12)

The length of the vector /~(t) is at most n 2, the number of entries of em. Since many of the entries of H are zero for both the triangular and Jordan forms, it is clear that a more compact expression is possible• Below, for both the triangular and Jordan forms, we derive n × n matrices S and S {°) which satisfy equations (11) and (12), where /~(t) is an n-dimensional vector of functions tke a' of time, t, and the eigenvalues, A, of A. As a result, equation (2) becomes

6_ ( t ) = s E ( t ) • u( t ) + J1

(11)

"'"

J/=

0

0 0

Ai

1.

• "•

0

,~i,

(10)

Let H be an upper triangular matrix such that

A = PHP-1. The elements of H are denoted hij where h i i = )ti, the i th eigenvalue of A. We assume Ai m Aj for i # j , and defer to Section 6 for consideration of systems with non-distinct eigenvalues• The Schur decomposition (see, for example, [6]) guarantees that such an H can be found by unitary transformation. Although P could be taken to be unitary, our analysis will not use this fact.

2O9 If we let Z(t) denote P-l~ff(t), then the differential equation for Z is obtained by applying P-1 to equation (1) is

The convolution of two distinct exponential functions is a linear combination of the two functions. Direct computation yields the identity:

OZ

e*" * e*, ' * u( t )

dt (t)

=P-1A~_(t)+e-lFu(t)

1

- --[ea"*u(t) A i - hj

= P - IAPZ( t ) + P - IFu( t ) = HZ(t) +Ru(t)

(14)

with initial condition Z(0)=P-l~ff(0) and /~= P - I / ~ H and P can be computed numerically, and since H is triangular, we determine a closed form solution for equation (14) by convolution and back substitution. The individual differential equations are (from the bottom up):

-eait*u(t)]

(17)

Using this identity and induction on [k, k = n . . . . . 1), it is easy to show that Zk(t) can be expressed as a linear combination of e xit and e xit * tt(t), with j > k:

zk(t ) = ~ qtkje x / * u ( t ) + ~ qSkjeaJt, j=k

(lS)

j=k

dz n dt (t) =AnZn(t ) +rnu(t)

where the coefficients qJkj are determined recursively from the relations

dzn_ 1

gJkj = 0 if k > j

(19)

r.

(20)

dt

(t) = A . _ , z . _ , ( t ) +hn_lnZn(t ) 6nn =

+ rn_ lU( t ) 1

~kj dz k dt ( t ) =A~zk(t ) +

~ hkjZj(t ) +rku(t ) j=k+l

dt (t) = a l Z , ( t ) + ~ h l j z j ( t ) + q u ( t ) j=2 where r i denotes the i th element of the vector/~. The n 'h equation is solved immediately:

n

(21)

E IlJjl if j < n l=j+ 1

(22)

&kj = 0 if k > j

(23)

~b,,,, = z,,(O)

(24)

1 ¢DkJ

- - -

J

~

h j -- h k / = k + l

dpjj = zj(O) -

~

hklcktj if k < j

chit if j < n

(25)

(26)

l=j+ l

E hkieXk'*zj( t ) + r k e x k ' * u ( t ) j=k+l +eXk'z~(0)

hktt~tj if k < j

n

tlljj = rj --

(15)

The solution of the k 'h equation (k < n ) is a convolution of e ~kt with the input function u(t) and with z/s from lower in the triangle:

zk(t)=

Aj - Ak t=k+l

and 4'kj are determined from the closely related formulas:

dz 1

z . ( t ) = r.e x.' * u ( t ) + eX-tz.(0).

J

(16)

The recursion formula above was obtained directly by back-solving the triangular system of differential equations. However, the existence of

210

such a recursive formula can be inferred from a result of Parlett [7]. In vector notation, equation (18) can be expressed compactly by

Z(t) = ~E(t)*u(t)

+~E(t)

(27)

where ~ and @ are the n × n matrices with ~. and &k~,respectively, in the k th row and Jk~ column. E ( t ) is the vector with jth element E j ( t ) • " = e Ajt . The original system vector Q-~( t ) is related to Z-(t) by Q~(t) = P Z ( t ) , and therefore

= 1 . . . . . r/b, k = 0 . . . . . n i - 1, where n b is the number of blocks and n i is the size of the i th block. If we let N~ denote the first index of the i th Jordan block (that is, N 1 = 1 and N i = 1 + Eij-=~ni) then the entries of /~ are given explicitly by E j ( t ) = tke ~,t, where i satisfies N i _
1

Ni+l-(k+l)

n

E c=N i

PiL, E

Srj = ~

Pt,- 1+k,wfw,

(31)

w=l

(i,j) entry of P-~. T h e coefficient of t k e xit in the r th entry of the initial term is

w h e r e pi~ 1 denotes the

Q_.( t ) = SIE( t ) • u( t ) + S'°)lE( t ) ,

(28)

where S = p q t and S ~°) = Pq).

1

N'+l-(k+l)

n

L'=N,

w=l

Pv +k,wqw(0) •

-r)

6. Solution based on the Jordan Canonical Form for systems with repeated eigenvalues The Jordan matrix Ht is block diagonal and its exponential is the block diagonal matrix of exponentials e J~t. The exponential of a Jordan block is computed directly from the identity

(32)

As in the case of the triangular system, we again obtain the system solution as a function of an n-dimensional v e c t o r / ~ ( t ) of simple functions of the eigenvalues of A, in this case tkeAjt:

Q.(t) = S I E ( t ) * u ( t )

+S'°)l~(t).

(33)

(29)

7. A recursive formula for calculation of the measurement function

where N, is the n~ × n i matrix Ji - Aft, which has l's on the super diagonal and O's elsewhere. When expanded this yields:

The net result of our transformations has been to express the system solution in terms of n 'basis' functions of the form e At or t~e ~t, where the As are the distinct eigenvalues of A. Now we can formulate the m e a s u r e m e n t function Yt in equation (4) in terms of integrals of convolutions of the basis functions with the input function:

k=o

eJi t : eAi t

k!

1

t

0

1

'

t 2

tki -1

2

( k i - 1)! tki-2

t

• . •

Yt = [1, 1 . . . . . 11 ( k i -

0 ,0

.-.

1 0

2) !

t 1 (30)

From here, it is clear that /~(t) can be taken to be the n-dimensional vector with entries t k e ;tit, i

=t1,1

. . . . .

I.

It

0

211

If the convolution /~*u is available in closed form, then Yt can be readily computed. However, in most cases of practical interest, this is not the case. PET input functions are generally given as discrete data points which are fit piecewise. The computational effort is significantly reduced by changing the order of integration and using the closed from integral of /~ to obtained a single, rather than double, integral:

The second integral of equation (37) looks new, but can actually be constructed from the previously computed value of (G*u)(tt). To see this, we first observe that for any t and 7,/~(t + 7) can be factored into C(t)t~(7), where C(t) is a block diagonal matrix, with blocks conforming in size to the Jordan blocks. (The size of all of the Jordan blocks is 1 × 1 for the case of distinct eigenvalues.) C(t) is derived below:

Ej( t + 7) = eAi¢'+~)(t + 7) k

fti'+'(E~* u)(t)dt

k

tl tl+l~ =/oft; (t-7)dtu(7)d7 k

tl+l

=eait E (k)tk-r(7re'~'r)

tl+l

+ f t, f~ E(t-r)dtu(7)d7 = (G* u)(tt+l) - (G*

u)(tt)

(38)

r=0

(35)

From the above, it follows that the C(t) is:

i th

Ci(t)

where

G( t)=foE( S)ds

(36)

= eait Because of the convolution of ~ff with u, each piece of u is involved in all subsequent measurements yl. Additional simplification is desired to avoid repeatedly integrating the same pieces of u. Note that for any partition [0 ~ t~, t 2 , . . . , t N = T] of the time interval [0, T], (G * u)(0) = 0, and therefore, the first simulated measurement has only one convolution, ~ff* u(t2). For each subsequent interval, ( G , u)(t t) will be available from the previous interval, so that only one convolution, G , u ( t l + l ) , is computed. This remaining convolution integral can be broken into two pieces: a weighted integral of u(t) on the current interval [tl, tt+l] , and a second integral which involves only previous values of u(t).

2

o

o

...

1

0

'-.

2t

1

block of

!/

. . . . .

(39) Now G(t + s ) can be decomposed as G ( t ) + as follows:

C(t)G(s), 6(t+s)

((~r*u)(tt+,)= fti"+'(ff(tt+,-7)u(r)d7 tl

+ fo G(t,+,-r)u(r)d7

(37)

= 6(t) +

(4o)

212

Thus, as we claimed earlier, the second integral in (37) can be simplified:

TABLE 1 C o m p u t a t i o n a l effort Closed f o r m solutions

fot'G(tl+l -- ~')u(~')d, 11

=fo(G(tl+l-tl) + C( tt+ 1 - tt)G( t t - "r) )u( ~")d~" = 6(t,+,

+

-

- t , ) ( 6 , u)(t,)

(41)

which is a linear transformation of the previously evaluated convolution ( G * u)(t t) plus a multiple of the integrated input function. In summary, we see that when the measurements {Yt, ! = 1, 2 , . . . , N} are computed in sequence, the only integrations which are needed to compute the measurement on a particular interval are the product

6(t,+,

(42)

from equation (41) and the first integral of equation (37):

ft,+, G(tt+l-~

- r)u(r)dr

(43)

tl

Evaluation of equation (43) depends on the specific form of the input function u(t). In the appendix we show how this double integral was computed for piecewise linear input functions.

8. Implementation We use ACM Algorithm 560 written by Bo Kagstrom and Axel Ruhe ([8] and [9]) to obtain the triangular and Jordan forms of the rate ma-

N u m e r i c a l O D E solver

Initial values

Fitted values

K 1=0.2

K I=0.12570

k 2 = 0.3

k 2 = 0.19454

k 3 = 0.1

k 3 = 0.08580

k 4=0.1

k 4=0.00644

K l=l.0

K I=0.12570

k 2 = 1.0

k 2 = 0.19454

k 3=0.0 k 4 = 0.0

k 3=0.08580 k 4 = 0.00644

K 1=0.3 k 2=0.3

K 1=0.12570

k 3=0.0

k 3=0.08580

k 3 = 0.08580

k 4 = 0.3

k 4 = 0.00644

k 4 = 0.00644

k 2=0.19454

CPU Fitted values t i m e (s) 1

5.2

K I = 0.12570 k 2 = 0.19454

CPU t i m e (s) 1

10.3

k 3 = 0.08580 k 4 = 0.00644

8.9

K t = 0.12570 k 2=0.19454

18.9

k 3 = 0.08580 k 4 = 0.00644

11.6

K t = 0.12570 k 2 = 0.19454

23.0

I On a S U N Sparcstation IPX Model 4/50.

trix. This code was modified slightly to make it compatible with the rest of our code. The P E T measurement double integral was implemented in house using the equations derived above. Closed form solutions are the default option in our simulation and parameter estimation package. These programs are invoked daily within our research group for the analysis of data from as many as 20 or 30 regions of interest in each P E T slice. The closed form solution method takes about half as much computational effort when compared with the numerical differential equation solver. Results of several fits to the data from Section 3 are shown in Table 1. The two columns compare the closed form solution techniques described here with a Runge-Kutta type numerical differential equation solver. The separate rows of the table show the results obtained with three different sets of initial conditions. The third set of initial conditions produces an initial system matrix with non-trivial Jordan form having two eigenvalues at 0.3. A region of interest in a P E T data slice generally includes a small fraction of blood data due to vasculature in the tissue. This is usually evidenced in the P E T data set by a peak coinciding with the time of the peak in blood data. In the

213 brain, this fraction is usually 0.03 to 0.06. The fits in Table 1 have assumed a vascular fraction of 0.05. Figure 1 shows the PET counts predicted by the model with the parameters in Table 1. We have made our programs available to any interested party in a package called RFIT. Currently it is used by PET groups at the University of Utah, the University of Pennsylvania, the VA Hospital in Palo Alto, California, and Hopital d'Orsay in Orsay, France.

~p and Eq are related by: 4(l~,

U)

p--1

uq+r

- 4 + A A ' ") + r=0 E ( q + r ) !Ap-r (47) The integral of t% At has a simple expression in terms of the eless functions:

abt ke~tdt 9. Appendix: calculations for piecewise linear input

i:0

b-a)

(48)

Below, we describe the formula used to evaluate equation (43) for piecewise linear input functions u(t). The desired integral is

9.2. Solution in terms of eless functions ft,+,~(tt+l _ r)u(r)dr

(44)

II

We let {%, m = 1, 2 . . . . . M t} denote a partition of [t t, tt+ l] such that u(t) is linear on each piece. That is, u(t)=am'r+b m on ['rm,'rm+l]. Then I is a sum of integrals I m over each piece.

or, from the definition of ~ff,

ft,tl+l f ll+l E ( t - ' r ) d t u ( z ) d z .

(45) I =

Evaluation of this integral entails parallel evaluations for the entries e~'('-')(t- z) k of / ~ ( ( t - ~'), as shown below.

=

f"+lf"+'eA'('-"( t

EMt

-

dtu( "r)dz

f'm+lft/+lehi(t_r)(

t __ T ) k

m = 1 7"m

Xdt( am'C+ bin)dr

9.1. The eless functions

mt

For numerical stability, it is desirable to obtain an expression for the integral of equation (43) in terms of what we have called "eless" functions, Wq(A, u), defined below, d~l(1, u) is the built-in function expml available in UNIX.

E ~

(49)

m=l

where

' Wq(A, u) - ~-~ e ~ " -

ri

(46)

t _ z)k dt( amz + bm)d~" (50)

214

We obtained the following general formula for I m in terms of the eless functions. I m = eai(tt+l-r,,,+O( _ 1) k k[(7 m - tl+|)k-s+2(S X

-- 1)

(k-s+2)!

am(S - 2) (rm - tl+l)} X {am"r I + b m + 2(k - s + 3)

X <(.Ai, Tin+ 1 -- 'I'm) ) am(k + 2)! -F

Blood Institute of the U.S. Department of Health and Human Services under Grant No. HL25840. Dr Salmeron was supported by La Fondation pour la Recherche Medicale and l'Institute National de Recherche en Informatique et Automatique (INRIA).

2

O'Ok+3(~i, Tm+ 1 -- Tm)

--k!OXgk+l(/~i, 7m+ 1 --

× Li: ( a m r +

tt+l) (51)

A step by step annotated derivation of this result is available upon request

Acknowledgements

This work was supported by the Director, Office of Energy Research, Office of Health and Environmental Research of the U.S. Department of Energy under Contract No. DE-AC0376SF00098 and by the National Heart, Lung and

References [1] D.Z. D'Argenio, A. Schumitzky and W. Wolf, Simulation of linear compartment models with application to nuclear medicine kinetic modeling, Comput. Methods Programs Biomed. 27, (1988) 47-54. [2] B.M. Mazoyer, R.H. Huesman, T.F. Budinger, et al. Dynamic PET data analysis, J. Comput. Assist. Tomogr. 10 (1986) 645-653. [3] W.J. Jagust, J.P. Seab, R.H. Huesman, et al. Diminished glucose transport in Alzheimer's disease; dynamic PET studies, J. Cerebral Blood Flow Metab. 11 (1990) 323-330. [4] L. Sokoloff, M. Reivich, C. Kennedy, et al. The [t4C]deoxyglucose method for the measurement of local cerebral glucose utilization: theory, procedure, and normal values in the conscious and anesthetized albino rat, J. Neurochem. 28 (1977) 897-916. [5] C. Moler and C. Van Loan, Nineteen dubious ways to compute the exponential of a matrix, SIAM Rev. 20 (1978) 801-836. [6] G.H. Golub and C.F. Van Loan, Matrix Computations, The Johns Hopkins University Press Baltimore, 1983. [7] B.N. Parlett, A recurrence among the elements of functions of triangular matrices, Linear Algebra Appl. 14 (1976) 117-121. [8] B. Kagstrom and A. Ruhe, An algorithm for numerical computation of the Jordan normal form of a complex matrix, ACM Trans. Math. Software 6 (1980) 398-419. [9] B. Kagstrom and A. Ruhe, Algorithm 560: JNF, an algorithm for numerical computation of the Jordan normal form of a complex matrix IF2], ACM Trans. Math. Software 6 (1980) 437-443.