Mathematical Techniques for Solving Analytically Large Compartmental Systems

Paper MATHEMATICAL TECHNIQUES FOR SOLVING ANALYTICALLY LARGE COMPARTMENTAL SYSTEMS Guillermo Sanchez* and Jesus Lopez-Fidalgo† (output). Nevertheless, any compartment will be considered as a potential initial and final compartment. The intake could be a fixed or a random quantity, depending on time or even other features of the substance and the context. In this paper a constant rate of transfer from one compartment to another is assumed. In order to solve a general model, the network will be decomposed in subsystems. When possible, subsystems will be decomposed in catenary unidirectional chains. A general introduction to this theory can be found in Anderson (1983) and Jazquez (1985). Assume a general compartmental model with compartments denoted by numbers, i ⫽ 1, 2, 3, . . ., n ⫹ 1. A dummy compartment, say n ⫹ 1, is considered as the environment, receiving flow but not giving input again to the system. This compartment will include the flow corresponding to the disintegrating rate of each compartment. Let k ij be the rate of transfer from compartment i to compartment j. Thus, the total rate of outputs from n⫹1 compartment i will be k i ⫽ 兺 j⫽1, j⫽i k i, j , i ⫽ 1, 2, 3, . . ., n ⫹ 1. If k i ⫽ 0, then compartment i is a trap compartment. This is always the case of the dummy compartment n ⫹ 1, k n⫹1 ⫽ 0, since k n⫹1, j ⫽ 0, j ⫽ 1, 2, . . ., n. Let q i (t) be the retention in compartment i at time t. Let b i (t) be the input coming from the environment to compartment i at time t. The derivative of q i (t) represents the variation of the retention in compartment i at time t: Abstract—The flow of radioactive particles inside the body from internally deposited radioisotopes in people exposed to inhalation, ingestion, injection, or other ways is usually evaluated using compartmental models. The biokinetic models included in the documents of the International Commission on Radiological Protection such as International Commission on Radiological Protection 66 and 78 involve many compartments. Usually numeric methods are applied. Very often analytical solutions are not possible. New computer programs that include symbolic capability can be used to solve compartmental systems. In this paper some techniques are developed in order to make feasible a computer program that gives not only faster and more accurate solutions, but also analytic solutions for these kind of models. The main idea is to make a partition of subsystems and solve them sequentially. The concept of pseudotrap compartments in a subsystem is crucial at this point. Impulse (acute), constant, and continuous (such as exponential) intakes are considered. This technique has been applied to develop a computer code called Humorap to solve the International Commission on Radiological Protection 66 and 78 models. Health Phys. 85(2):184 –193; 2003 Key words: modeling, biological factors; biokinetics; radiation protection; dose, internal INTRODUCTION A compartmental model is a network where the nodes are compartments connected by arrows designing the flow of some substance from one to another. In particular, the flow and retention of some kind of “substance” will be considered. There are initial compartments where the intake (input) of the substance takes place and final compartments from which the substance is eliminated dq 1 共 t 兲 ⫽ dt 冘k r,1 q r 共 t 兲 r ⫺ 冘k 1, j q 1 共 t 兲 ⫹ b i共 t 兲 j ⫽ 冘k r,1 q r 共 t 兲 ⫺ k 1q 1共 t 兲 ⫹ b 1共 t 兲 (1) r * ENUSA, Fca Juzbado, Apdo 328, 37080-Salamanca, Spain; Department of Statistics, Plaza de los Caidos s/n, 37008 — Salamanca, Spain. For correspondence or reprints contact G. Sanchez, ENUSA, Fca Juzbado, Apdo 328, 37080-Salamanca, Spain, or email at guillerm@usal.es. (Manuscript received 3 December 2001; revised manuscript received 10 February 2003; accepted 5 May 2003) ··· † dq n 共 t 兲 ⫽ dt 0017-9078/03/0 Copyright © 2003 Health Physics Society 冘k r r,n q r 共 t 兲 ··· ⫺ ··· ··· 冘 k q 共t兲 ⫹ b 共t兲 ⫽ 冘 k q 共t兲 ⫺ k q 共t兲 ⫹ b 共t兲 n, j n i j r,n n r 184 ··· n n i Analytical solutions for compartmental systems ● G. SANCHEZ and the variation of quantity eliminated at time t for the whole system is dq n⫹1 共 t 兲 /dt ⫽ 冘k r r,n⫹1 q r 共 t 兲 . This differential equation can be written as dq 共t兲 ⫽ A q共t兲 ⫹ b共t兲, dt (2) where the matrix A will have non-positive elements in the diagonal and non-negative elements elsewhere: A⫽ 冢 ⫺k1 k12 ··· k1,n k1,n⫹1 k2,1 ⫺k2 ··· k2,n k2,n⫹1 ··· ··· ··· ··· ··· kn,1 kn,2 ··· ⫺kn kn,n⫹1 冣 0 0 ··· , 0 0 冢冣 冢 冣 Every column sums up to zero, and the last column is zero since there is no output from the dummy compartment. The input from the environment to compartment i, b i (t), is considered separately in the equation and could be any kind of function. Usually, this function is an impulse (acute) input at time t ⫽ 0, or constant (chronic). The general solution of this differential equation, given an initial value q(0) ⫽ q0, is well known: 冕 t e共 t⫺ ␶ 兲 Ab共␶兲d␶. J. LOPEZ-FIDALGO 185 ␭ jt q i (t) ⫽ 兺 j a i, j e ⫹ ␤ i (t), where ␤ i (t) depends on the expression of b(t), which is assumed to be known. If b(t) is constant (chronic), or multiexponential, or impulse input then ␤ i (t) will also be multiexponential. Frequently, many of these terms are negligible leading to much simpler models. The convolution theorem, q(t) ⫽ 兰 t0 q0 (t ⫺ ␶ ) b ( ␶ ) d ␶ , provides a general solution, q(t), with timedependent input, b(t), when the solution q0 (t) for a unit impulse input is known. CATENARY UNIDIRECTIONAL SYSTEMS b1 共t兲 q1 共t兲 b2 共t兲 q2 共 t兲 b共t兲 ⫽ · · · , q共t兲 ⫽ · · · . (3) bn 共t兲 qn 共 t兲 0 qn⫹1 共t兲 q 共t兲 ⫽ etAq0 ⫹ AND (4) A catenary unidirectional system is the most simple compartmental model. It is a sequence of compartments in such a way that each one receives flow from the previous one, with rate k i , and it gives flow to the next one, with rate k i⫹1 . The first one receives flow from the environment, b 1 (t), and only the last one gives flow to the environment, with rate k n . Disintegrating rates, k i,0 , i ⫽ 1, 2, . . ., n, for each compartment will be considered as well. Fig. 1 shows this situation. In fact, if it is possible large systems are decomposed in catenary chains. Then, where it has been said flow to the environment it would be said to flow to the environment or other compartments outside the chain. The differential equations for this model are now much simpler: dq 1 共 t 兲 ⫽ b 1共 t 兲 ⫺ K 1q 1共 t 兲 , dt (6) dq r 共 t 兲 ⫽ k r⫺1 q r⫺1 共 t 兲 ⫺ K r q r 共 t 兲 , r ⫽ 2, . . . , n dt 0 This is an abstract formula to work with. A simple assumption on A is made now in order to show what the solution looks like in this case. Let A be diagonalizable with real eigenvalues. This happens when rank (A ⫺ ␭ i I) ⫽ n ⫺ m i , where m i is the multiplicity of the eigenvalue ␭ i . For instance, if all the eigenvalues are different, this is automatically satisfied. Let ⌳ ⫽ diag (␭1, . . ., ␭n , 0) be the diagonal matrix of the eigenvalues and S be the matrix with the eigenvectors e1 , . . ., en , en⫹1 as columns, then e tA ⫽ S e t⌳ S⫺1. The solution of the n⫹1 homogeneous equation is q(t) ⫽ 兺 i⫽1 ␥ i ei e ␭ it , where ⫺1 ␥ ⫽ S q0 , and the general solution is 冘␥ee n q 共 t兲 ⫽ i i i⫽1 ␭ it 冕 dq n 共 t 兲 ⫽ k nq n共 t 兲 , dt where K i ⫽ k i ⫹ k i0 , i ⫽ 1, 2, . . ., n and K 0 ⫽ 0. The vector of inputs is now bT (t) ⫽ [b 1 (t), 0, . . ., 0]. If the rates of transfer, K i , in a unidirectional catenary system are all different, then the matrix A for this system is diagonalizable. As a matter of fact, an t ⫹ S e共 t⫺ ␶ 兲 ⌳ S⫺1 b共␶兲d␶. (5) 0 Notice that ␥ n⫹1 ⫽ 0 if 兩An⫹1,n⫹1 兩 ⫽ 0. Summarizing, the retention in compartment i at time t can be modeled as Fig. 1. Catenary unidirectional system with n compartments. The retention in each compartment q 1 (t), q 2 (t), . . ., q n (t), . . . is computed in the usual way. 186 Health Physics explicit solution for the retention is derived straightforward for the case of impulse input b 1 (Skrable et al. 1974): 写k 兲冘 i⫺1 q i共 t 兲 ⫽ b 1共 i p p⫽1 j⫽0 冢 e ⫺K jt 写 i p⫽0,p⫽j 共 K p ⫺ K j兲 冣 i ⫽ 1, 2, . . . , n. (7) The advantage of a catenary unidirectional model is that the differential equations can be solved sequentially as each one involves only the retention of the previous compartment. Thus, it would be very convenient to decompose a system in catenary unidirectional chains, if possible. From a theoretical point of view it can be said that if A in the general model given by eqn (2) is diagonalizable and the inputs are multiexponential, then the system is equivalent to a unidirectional catenary one. INTERCONNECTED SUBSYSTEMS Frequently, in a family of compartmental systems there are subsystems that are common to all of them. In order to solve these models for different possible values of the rates of transfer, using computer programs, it is convenient to make a partition of the whole network in subsystems and take advantage of the solution of the common subsystems for several models. This makes the program faster and it can be applied to a variety of networks. A general review of mathematical procedures using matrix methods has already been made (Polig 2001). It is stated that mathematical software, as is the case of Mathematica (Wolfram Research, Inc., Champagne, IL) allows symbolic calculation, making feasible obtaining analytical solutions. So far, the usual methods have been numerical. Thus, it is needed to define properly what a subsystem is and how to solve the problem “independently” in the different subsystems. This has to be done sequentially. The main idea to define a subsystem is the concept of pseudotrap compartments. They have to satisfy the following conditions: 1. They are the only compartments in the subsystem which have flow to another subsystem; 2. There is not flow from them to other compartments in their own subsystem; 3. Each of them has flow to one or more compartments in just one other subsystem; 4. Recycling involving compartments of different subsystems is not allowed. August 2003, Volume 85, Number 2 An example is given in Fig. 2, where the pseudotrap compartments of subsystem P satisfying these conditions are i, i⬘, and i⬙. There is recycling when the flow returns to the original compartment throughout others. Thus, a set of compartments ᐀ ⫽ {␬i,i 1, ␬ i 1, . . ., ␬ i p,i } forms a cycle if 0 ⰻ ᐀. If none of the subscripts is repeated more than once, it is a proper cycle. Some authors (Godfrey 1983) call them unilateral circulation models. For the sake of simplicity, it will be denoted ᐀ ⫽ {i, i 1 , i 2 , . . ., i p }. A pseudotrap compartment in a subsystem, P, is always associated with another specific subsystem, B, but flow to this compartment from a third subsystem is also allowed. For a better understanding of what follows have a look at Figs. 2 and 3. A subsystem is solved independently considering the retention in the pseudotrap compartments as if they did not flow to any other compartments. Once a subsystem is solved, then the pseudotrap compartment becomes a compartment of the second subsystem, say B⬘. The external input for this new compartment will be the corresponding retention considered in the first subsystem, as well as some possible input coming directly from the environment or other subsystems. Let P and B be subsystems with the conditions above with transference from a pseudotrap compartment Fig. 2. Subsystem process, step 1. The retention in compartments q i (t), q i⬘ (t), and q i⬙ (t) are computed in different ways. Analytical solutions for compartmental systems ● G. SANCHEZ AND J. LOPEZ-FIDALGO 187 of a subsystem will be multiexponential. This means the input for the next subsystem will be multiexponential. This is the reason which makes important these kind of input functions. RETENTION IN A COMPARTMENT WITH KNOWN INPUTS Fig. 3. Subsystem process, step 2. i in P to some others in B. Once the subsystem P has been solved, q Pi (t) will denote the retention in i assuming there is not flow to any compartment in B. Subsystem B⬘ includes i and other possible pseudotrap compartments from P or other subsystems. A precise statement can be made here to clarify the process. Flow from a third subsystem C to compartment i is possible. If subsystem C has been solved previously and the transference from it has already been considered in the solution of P, then subsystem C will not be considered any more in order to solve B⬘. If, otherwise, subsystem C has not been solved before, then the corresponding transference from C to i would have to be considered at the moment of solving B⬘ as well as other possible inputs from the environment. A possible transference from a compartment in B to i will be considered as usual in the solution of B⬘. From this point of the process on the compartment i will be considered as a compartment of subsystem B⬘ and the problem will be solved as usual considering the input to i as follows. From the general theory for impulse, constant, or multiexponential inputs, the retention in a trap compartment, comes from the solution of the corresponding differential equation, say q Pi (t) ⫽ d ⫺ 兺 s d s e ⫺c st , where d, d s , and c s are the coefficients of a particular solution. Considering i as a compartment of B, the retention in this compartment can be computed from the differential equation, dq i (t)/dt ⫽ b i (t) ⫹ [dq Pi (t)/dt] ⫺ k i q i (t). The rate k i is computed considering the compartments in subsystem B having flow from i. The input is now dq Pi (t)/dt ⫽ 兺 s d s c s e ⫺c st and, following the precision made above, any other possible inputs, b i (t). The solution to this differential equation is the same as that of the general problem with variable inputs. In Figs. 2 and 3 the whole process is shown in a very general case assuming the properties given before. If the input from the environment is constant or impulse, then the retention in a pseudotrap compartment In this section the problem of solving a subsystem in an efficient way is considered. Let j be a compartment receiving radioactivity from i. The retention in j considering only the input coming from i will be denoted by q (i) j (t). Assuming that q i (t) is known then the variation of the retention in j due to i is dq (i) j (t)/dt ⫽ k i, j q i (t) ⫺ k j q (i) (t) ⫽ k q (t). Using b̃ i, j i, j i (t) as the input quantity, j the usual differential equation is obtained: dq j共 i 兲 共 t 兲 ⫽ b̃ i, j 共 t 兲 ⫺ k j q 共j i 兲 共 t 兲 , dt (8) which is solved as usual, q j共 i 兲 共 t 兲 ⫽ q j共 i 兲 共 0 兲 e ⫺k i, jt ⫹ 冕 t k i, j q i 共 t 兲 e ⫺k i, j共 t⫺ ␶ 兲 d ␶ . 0 (9) The final retention in compartment j will be q j (t) ⫽ 兺 { j:k i, j⫽0} q (i) j (t). Therefore, assuming all the inputs to j are known, the retention in j can be easily computed. If there is no recycling, this suggests a method to compute the whole system sequentially: 1. Choose a compartment, say i ⫽ 1, which receive known input only from the environment. Then the variation of retention in it is dq 1 (t)/dt ⫽ b 1 (t) ⫺ k 1 q 1 (t), that is a form of eqn (8) and the solution will be as in eqn (4). This has to be done with all the compartments of this type; 2. Choose a compartment, say i ⫽ 2, receiving flow only from the environment and/or compartments of type 1. Now the variation on the retention in this compartment will be dq 2 (t)/dt ⫽ b 2 (t) ⫹ 兺 {i:k 1,2⫽0} k i,2 q i (t) ⫺ k 2 q 2 (t), which is again a form of eqn (8). All the compartments of this type should be solved in this step; and 3. Choose a compartment, say i ⫽ 3, receiving flow only from the environment, and/or compartments of type 1, and/or compartments of type 2. Then dq 3 (t)/dt ⫽ b 3 (t) ⫹ 兺 {i:k i,3⫽0} k i,3 q i (t) ⫺ k 3 q 3 (t), which is again a form of eqn (8). All the compartments of this type should be solved in this step. This process is possible in each step and ends if, and only if, there is not recycling in the system. In this case, 188 Health Physics the type of inputs considered above leads to multiexponential solutions. In the next section the case of recycling is considered. SYSTEMS WITH RECYCLING It can be assumed without loss of generality that ᐀ ⫽ { 1, 2, 3, . . ., r ⫺ 1, r} is a proper cycle. The part of the matrix corresponding to this cycle is then A᐀ ⫽ 冢 ⫺k1 k1,2 ··· 0 0 ⫺k2 ··· 0 0 0 ··· 0 ··· ··· ··· ··· 0 0 ··· ⫺ kr⫺1,r 冣 kn,1 0 ··· . kr It will be assumed that all the inputs to the cycle are known. Let b̃ i (t) be the total input to compartment i. Then dq ᐀ (t)/dt ⫽ b̃ ᐀ (t) ⫹ A ᐀ q ᐀ (t), where q T᐀ (t) ⫽ (q 1 , . . ., q r ) and b̃ T᐀ (t) ⫽ (b̃ 1 , . . ., b̃ r ). The solution is obtained as usual, q ᐀ (t) ⫽ e tA ᐀ q ᐀ (0) ⫹ 兰 t0 e (t⫺ ␶ ) A ᐀ b̃( ␶ )d ␶ . An explicit solution of a cycle is not easily derived in general. This means, anyway, a significant reduction in computation if a decomposition in small cycles is possible without recycling between them. In this case the procedure given in the last section can be applied here changing single compartments for cycles. A necessary and sufficient condition in the general matrix can be given for such decomposition. Without loss of generality two cycles can be assumed, ᐀1 ⫽ {1, 2, 3, . . ., r ⫺ 1, r 1 } and ᐀2 ⫽ {r 1 ⫹ 1, r 1 ⫹ 2, . . ., r 1 ⫹ r 2 }. Then A ᐀ 1 and A ᐀ 2 have to be of the form eqn (10) to be proper cycles. Thus, the partitioned matrix corresponding to this two cycles is A ᐀ 1 A 2,1 A 1,2 A ᐀ 2 冊 . a) The matrix-exponential e At has been computed using e tA ⫽ S e t⌳ S⫺1 ; b) When b(t) ⫽ b (constant) has been applied eqn (12), in the place of the well known 兰 t0 e (t⫺ ␶ )Ab( ␶ )d ␶ ⫽ A⫺1 (e tA ⫺ I) b. It has been like that because to compute A⫺1 requires no trap compartments (Anderson 1983): 冕 t (11) There is no recycling between them if, and only if, A 2,1 ⫽ 0 or A 1,2 ⫽ 0. Decomposition of a system verifying these conditions could be not unique. A convenient decomposition in order to get faster computation will look for as many cycles as possible, as well as low order for them. HUMORAP: A PROGRAM FOR SOLVING ICRP 66 AND 78 MODELS The previous ideas have been applied to develop a program, called Humorap (Sanchez 2002), which solves the International Commission on Radiological Protection (ICRP) 66 and 78 biokinetic models analytically. Humorap is a package of Mathematica. It has been programmed using the symbolic capability of Mathematica. The solution of the models has been obtained applying eqn (4) using the mathematic criteria which it follows: 冕 t e 共 t⫺ ␶ 兲 Ab 共␶兲d␶ ⫽ b e␶ Ad␶; 0 (10) 冉 August 2003, Volume 85, Number 2 (12) 0 c) When the terms of b(t) are sum of exponentials or 0, then each term of e (t⫺ ␶ ) A b(␶) is exponential or zero and they can be solved using eqn (13). It has been programmed as a rule to be applied whenever this pattern appears: a 冕 t e bt⫹ 共 c⫺b 兲 ␶ d ␶ ⫽ a 共 e bt ⫺ e ct 兲 b⫺c if b ⫽ c; 0 (13) d) When b(t) are not impulsive inputs, constants or exponentials q(t) is obtained applying the convolution theorem. To describe the specific method applied for Humorap to solve ICRP 66 and 78 it is convenient to summarize them: ● ICRP 66 describes the compartmental model of the human respiratory tract (RT) applied to the intake of radioactive aerosols by inhalation (Fig. 4). In this model the material is deposited in the respiratory tract: in compartments labeled “Particles in Initial State” (PIS) and in ET1. From each PIS compartment the material is transferred into the body fluids at an absorption rate s p . It is also simultaneously transferred from PIS (at a rate s pt ) to a corresponding compartment labeled “Particles in Transformed State” (PTS). The flow goes from 1, in PIS, to 1, in PTS, from 2, in PIS, to 2, in PTS, and so on. We can consider that each compartment 1-13 in PIS has a mirror compartment 1-13 in PTS. In each compartment in PTS the isotope is dissolved at a constant rate s p into the body fluids (usually the blood). For instance: the total transfer rate for AI2 in PIS will be K AI2 ⫽ k 2,4 ⫹ s pt ⫹ s p , and for AI2, in PTS, K⬘AI2 ⫽ k 2,4 ⫹ s t . This general model of RT is common to any element. The absorption rates {s pt , s p , s t } are related with the chemical form of the element. ICRP 66 establishes three types of materials according to its absorption behavior: fast, F, moderate, M, and slow, S. Analytical solutions for compartmental systems ● G. SANCHEZ AND J. LOPEZ-FIDALGO 189 Fig. 4. The dashed arrow from subsystem PIS to subsystem PTS means that the flow goes from each compartment in PIS to the compartment with the same number in PTS. The hollow arrow, f, means a flow from each compartment in subsystem PIS or PTS to the “Body fluids.” A simple arrow, 3, means flow from a single compartment to other. ● ICRP 78 expands ICRP 66 including all the compartments that take part in the metabolic process starting in the intake and ending in the fecal and urine excretion. The intake can be by inhalation, ingestion, or injection (Fig. 5). ICRP 78 applies for gastrointestinal tract (GI) the same model as in ICRP 30. The models are specific for groups of elements. ICRP 78 establishes three generic groups: i) hydrogen, cobalt, ruthenium, cesium, and californium, ii) strontium, radium, and uranium (Fig. 6), 190 Health Physics August 2003, Volume 85, Number 2 considered trap compartments. Subsystem P is divided in to two subsubsystems: P1) PIS is solved as a subsubsystem, with compartmental matrix A1. The compartments are numbered as it is shown in Fig. 4. Compartment B is number 14 and ST (first compartment of GI) is number 15. The retention in each compartment for an impulse intake I0 in t ⫽ 0 is given by q A1 共t兲 ⫽ I0 IDFetA1 , Fig. 5. ICRP 78 Biokinetic Model. A dashed arrow means that the flow is possible for some isotopes and not for others. A hollow arrow means flow from a subsystem while a simple one means flow from a single compartment. (14) where IDF ⫽ (IDF 1 , . . ., IDF 9 , 0, IFD 11 , IDF 12 , 0, 0, 0) are the initial deposition factors. They are functions of the activity median aerodynamic diameter (AMAD), physiological factors of the subject as well as various conditions of exposure. IDF i ⫽ 0 in compartments i ⫽ {10, 13, 14, 15}, because they are not receiving flow from outside the system. IDF can be calculated following the procedure described in ICRP 66 or obtained from annex F of ICRP 66. P2) PTS is solved as a subsubsystem, with compartmental matrix A2. It is taken into account that from each compartment there is a flow to PTS from PIS with a rate of transfer s pt , which is the same for all compartments. That is, the input into i⬘ from i is b i ⬘ (t) ⫽ s pt q (t), so the retention in PIS compartments is given by 冕 t q A2 共t兲 ⫽ spt e共 t⫺ ␶ 兲 A2 qA1 共␶兲d␶. (15) 0 Fig. 6. Biokinetic model for uranium, lead and strontium (based on figure 8 of ICRP 78). and, iii) thorium, neptunium, plutonium, americium, and curium. For the elements of each group the same model is applied. However, some parameters are specific of the element. The criteria applied in the development of the program consist of dividing the metabolic model in subsystems that are solved as follows: 1. Subsystem ET1. It is solved as a separated compartment; 2. Subsystem P. It represents the RT (except ET1) which has already been described (Fig. 4). The ending compartments of this subsystem are B and ST (which is the first compartment of GI). Both of them are It can be solved applying the rule given by eqn (13). The division in a compartment i in PIS and other compartment i⬘ in PTS is only for mathematical purposes. Each couple i ⫺ i⬘ should be considered as an only compartment to give the content in a specific compartment. Therefore, the retention in compartments of subsystem P will be qA (t) ⫽ qA1 (t) ⫹ qA2 (t). It is important to realize that the retention in compartments B and ST have the pattern q Ac B (t) ⫽ d B ⫺ 兺 s Ac d s e ⫺c st and q ST (t) ⫽ d ST ⫺ 兺 r d r e ⫺c rt (notice that in a pseudo-trap compartment the decay constant ␭R ⫽ 0); 3. Subsystem B. It represents the rest of the model (Fig. 5). In this case, because B and ST in subsystem P were assumed pseudo-trap compartments, the subsystem P (RT) is replaced with an input flow to ST and to B given by b ST (t) ⫽ ⭸q RT,ST (t)/⭸t ⫽ 兺 r c r d r e ⫺c rt and b B (t) ⫽ ⭸q RT,B (t)/⭸t ⫽ 兺 s c s d s e ⫺c st . Therefore, numbering B as compartment 1 and ST as n, in eqn (4), b( ␶ ) ⫽ (兺 r c r d r e ⫺c r␶ , 0, . . ., 0, 兺 r c r d r e ⫺c r␶ ) and developing e (t⫺ ␶ )A (兺 r c r d r e ⫺c r␶ , 0, . . ., 0, 兺 r c r d r e ⫺c r␶ ) is obtained so that all the terms are a sum of exponentials or 0. Then, eqn (4) can be solved applying the rule given in eqn (13). Analytical solutions for compartmental systems ● G. SANCHEZ So far the radioactive decay has not been taken into account, but it will be included in the final solution by using qR(t) ⫽ q(t) exp(⫺␭Rt) where ␭R is the decay constant of the isotope. By default Humorap applies biokinetic parameters for reference workers that are given in ICRP 78, but these parameters can be modified by users. The general model represents the inhalation, but the injection and ingestion intakes can be evaluated as well. Humorap can be downloaded from http://web.usal.es/⬃guillerm/biokmod. htm. The most important features of Humorap are as follows: 1. It gives analytic solutions, apart from numeric (Fig. 7); 2. All parameters (deposition factors, rates of transfer, rate of absorption, . . .) can be modified by users. It can be useful for sensitivity analysis and for fitting experimental data; 3. It can be used for any kind of continuous inputs (exponential, periodic, etc.), even for random inputs and not only for impulse and chronic inputs as most AND J. LOPEZ-FIDALGO 191 software do (For instance, qContinuous [b(t), inputimpulse, t, t 1 ] gives the content—activity, mass, etc.—in one compartment or region for a continuous intake b(t), e.g., a ⫹ cos (b t), during a time t 1 , input-impulse is the retention function for a single intake 1 in t ⫽ 0); 4. The user can build special compartmental models (e.g., a new model for an element not included in ICRP 78) in a very easy way; and 5. The package runs in the spreadsheet Excel (Microsoft Corporation, Seattle, WA), using MathLink for Excel (MathLink for Excel by Wolfram Research, Inc.) (Fig. 8). APPLICATION FOR URANIUM Here is an illustrative example that consists of solving the ICRP 66/78 uranium model with Humorap. For RT the ICRP 66 model (Fig. 4) is applicable. The content in RT compartments as a function of the time t can be obtained with the Humorap function qRT[I 0 , Fig. 7. The example corresponds to the thoracic lung retention and the daily urine and fecal excretion for worker exposed to inhalation of uranium class S, and AMAD 5 ␮m. (The negligible terms have been eliminated, it has also been assumed ␭R ⫽ 0). 192 Health Physics August 2003, Volume 85, Number 2 Fig. 8. Humorap using Excel as interface. The example corresponds to the thoracic lung retention for worker exposed to inhalation of uranium (S, AMAD 5 ␮m). The input given by user corresponds to the cells in grey background. IDF, FRA, t, ␭ R , options] where I 0 is the input at time t ⫽ 0, IDF are the initial deposition factors, FRA are the fractional rates of absorption and ␭ R is decay constant. With options the default rate transfer can be modified. So qRT[3, AMAD5, S, 30, 0] gives the content in each RT compartment for a reference worker 30 d after having intaken 3 Bq of aerosols of UO2 (class S) and AMAD ⫽ 5 ␮m, ␭ R is taken “0” because for 238U, 235U, and 234U decay constant ␭ R ⬇0. To include all compartments, apart from RT, the uranium model (UM), that is shown in Fig. 6, should be used. The RT and UM are connected through ST and B as usual. The compartmental matrix for UM, using ICRP 78 rate of transfer values, is included in Humorap. However, in some situations the user can be interested in defining the compartmental matrix (e.g., adding new compartments or modifying the rate of transfer default values). In these cases it can be used the Humorap function CompartMatrix[n, {{i, j, k ij }}], where n is the number of compartments and {i, j, k ij } means flow from i to j with a rate of transfer k ij . Using this function the compartmental matrix for UM will be given in the Mathematica style as CompartMatrix[22, {{1, 3, 10.5}, {1, 2, 0.245}, {1, 16, 15.43}, {1, 8, 2.94}, {1, 9, 0.0122}, {1, 18, 0.122}, {1, 6, 0.367}, {1, 4, 1.63}, {1, 5, 0.0735}, {1, 10, 2.04}, {1, 13, 1.63}, {3, 1, 8.32}, {2, 1, 0.347}, {8, 16, 0.099}, {9, 1, 0.00038}, {6, 1, 0.092}, {6, 7, 0.00693}, {4, 1, 0.0347}, {5, 1, 0.000019}, {13, 1, 0.0693}, {13, 14, 0.0693}, {7, 1, 0.00019}, {11, 12, 0.00578}, {12, 1, 0.000493}, {15, 1, 0.0000821}, {14, 13, 0.0173}, {14, 15, 0.00578}, {10, 1, 0.0693}, {11, 10, 0.0173}, {10, 11, 0.0693}, {16, 17, 12}, {18, 19, kULI}, {19, 20, kLLI}, {21, 1, f1 kSI/(1 ⫺ f1)}, {21, 18, kSI}, {22, 21, kST}}]. where the number of compartments is given in Fig. 6 and the rate of transfer, in d⫺1, in Table A.10.1 of ICRP 78. The solution given by Humorap for an impulse intake “1” in t ⫽ 0 for workers exposure to inhalation of Analytical solutions for compartmental systems ● G. SANCHEZ Table 1. Predicted values (Bq per Bq intake) for inhalation by workers of 234U, 235U, 238U (AMAD 5 ␮m and class S). t 1 2 3 4 5 6 7 8 9 10 20 30 40 50 60 70 80 90 100 200 300 400 500 600 700 800 900 1,000 2,000 Lungs Daily urinary ⫺2 6.43 ⫻ 10 6.27 ⫻ 10⫺2 6.19 ⫻ 10⫺2 6.13 ⫻ 10⫺2 6.07 ⫻ 10⫺2 6.01 ⫻ 10⫺2 5.95 ⫻ 10⫺2 5.90 ⫻ 10⫺2 5.84 ⫻ 10⫺2 5.79 ⫻ 10⫺2 5.32 ⫻ 10⫺2 4.94 ⫻ 10⫺2 4.63 ⫻ 10⫺2 4.39 ⫻ 10⫺2 4.18 ⫻ 10⫺2 4.02 ⫻ 10⫺2 3.88 ⫻ 10⫺2 3.76 ⫻ 10⫺2 3.67 ⫻ 10⫺2 3.12 ⫻ 10⫺2 2.82 ⫻ 10⫺2 2.57 ⫻ 10⫺2 2.34 ⫻ 10⫺2 2.14 ⫻ 10⫺2 1.96 ⫻ 10⫺2 1.80 ⫻ 10⫺2 1.65 ⫻ 10⫺2 1.52 ⫻ 10⫺2 7.28 ⫻ 10⫺3 ⫺4 7.04 ⫻ 10 4.41 ⫻ 10⫺5 2.60 ⫻ 10⫺5 2.37 ⫻ 10⫺5 2.20 ⫻ 10⫺5 2.05 ⫻ 10⫺5 1.92 ⫻ 10⫺5 1.80 ⫻ 10⫺5 1.70 ⫻ 10⫺5 1.60 ⫻ 10⫺5 1.02 ⫻ 10⫺5 7.72 ⫻ 10⫺6 6.44 ⫻ 10⫺6 5.69 ⫻ 10⫺6 5.18 ⫻ 10⫺6 4.81 ⫻ 10⫺6 4.52 ⫻ 10⫺6 4.28 ⫻ 10⫺6 4.09 ⫻ 10⫺6 3.18 ⫻ 10⫺6 2.82 ⫻ 10⫺6 2.55 ⫻ 10⫺6 2.33 ⫻ 10⫺6 2.13 ⫻ 10⫺6 1.96 ⫻ 10⫺6 1.80 ⫻ 10⫺6 1.65 ⫻ 10⫺6 1.52 ⫻ 10⫺6 7.52 ⫻ 10⫺7 Daily fecal 1.14 ⫻ 10⫺1 1.63 ⫻ 10⫺1 8.39 ⫻ 10⫺2 3.52 ⫻ 10⫺2 1.40 ⫻ 10⫺2 5.65 ⫻ 10⫺3 2.47 ⫻ 10⫺3 1.27 ⫻ 10⫺3 8.21 ⫻ 10⫺4 6.47 ⫻ 10⫺4 4.39 ⫻ 10⫺4 3.50 ⫻ 10⫺4 2.81 ⫻ 10⫺4 2.28 ⫻ 10⫺4 1.86 ⫻ 10⫺4 1.53 ⫻ 10⫺4 1.27 ⫻ 10⫺4 1.07 ⫻ 10⫺4 9.06 ⫻ 10⫺5 3.30 ⫻ 10⫺5 2.43 ⫻ 10⫺5 2.12 ⫻ 10⫺5 1.89 ⫻ 10⫺5 1.70 ⫻ 10⫺5 1.52 ⫻ 10⫺5 1.37 ⫻ 10⫺5 1.23 ⫻ 10⫺5 1.11 ⫻ 10⫺5 3.88 ⫻ 10⫺6 uranium aerosols UO2 (class S) and AMAD 5 ␮m is shown in Fig. 7. These solutions can be applicable to 238 U, 235U, and 234U because the decay constants for these isotopes is ␭⬇0. The numeric values for different days after an impulse intake “1” are shown in Table 1. For the lungs, they are also shown in Fig. 8, which includes the values for acute (1 Bq) and chronic intake (1 Bq d⫺1). DISCUSSION In this paper the rates of transfer are known. Then the analytical solution for the retention is given explicitly. If the inputs are impulse, the retention will be AND J. LOPEZ-FIDALGO 193 m j⫽1 ␭ it multiexponential. Let’s say q i (t) ⫽ 兺 a i, j e , where ␭ j is a known linear combination of the rates of transfer. They could be statistically estimated if it is possible to measure q i (t) empirically for different values of t. The best choice for these times is a matter of optimal experimental design (Lopez-Fidalgo and Rodriguez-Diaz 1998). After that, least squared estimates could be suitable. This is usually called the inverse problem. Having an explicit solution of a model is held to estimate the corresponding parameters and the rates of transfer for a specific problem. Acknowledgments—This research is supported by the grant from Junta de Castilla y Leon SA004/01. REFERENCES Anderson DH. Compartmental analysis and tracer kinetics. Berlin: Springer-Verlag; 1983. Godfrey K. Compartmental models and their application. London: Academic Press; 1983. International Commission on Radiological Protection. Human respiratory tract model for radiological protection. Oxford: Pergamon Press; ICRP Publication 66; 1994. International Commission on Radiological Protection. Individual monitoring for internal exposure of workers. Oxford: Pergamon Press; ICRP Publication 78; 1997. Jazquez JA. Compartmental analysis in biology and medicine. Ann Arbor: The University of Michigan Press; 1985. Lopez-Fidalgo J, Rodriguez-Diaz JM. Characteristic polynomial criteria in optimal experimental design. In: Advances in model oriented data analysis. Proceedings of the MODA5. New York: Springer-Verlag; 1998:31–38. Polig E. Modeling the distribution and dosimetry of internal emitters: A review of mathematical procedures using Matrix Methods. Health Phys 81:492–501; 2001. Sanchez G. Humorap: A program for solving the compartmental models of ICRP 66 and 78. In: Radioprotección Extra2002. Proceedings of the IX Meeting of Spanish Society of Radioprotection. Madrid: Senda; 2002:365–368 (in Spanish). Skrable KW, French C, Chabot G, Major A. A general equation for the kinetics of linear first order phenomena and suggested applications. Health Phys 27:155–157; 1974. f f