Analytical Solution of Multicompartment Solute Kinetics for Hemodialysis

Abstract

Objective. To provide an exact solution for variable-volume multicompartment kinetic models with linear volume change, and to apply this solution to a 4-compartment diffusion-adjusted regional blood flow model for both urea and creatinine kinetics in hemodialysis. Methods. A matrix-based approach applicable to linear models encompassing any number of compartments is presented. The procedure requires the inversion of a square matrix and the computation of its eigenvalues λ, assuming they are all distinct. This novel approach bypasses the evaluation of the definite integral to solve the inhomogeneous ordinary differential equation. Results. For urea two out of four eigenvalues describing the changes of concentrations in time are about 10⁵ times larger than the other eigenvalues indicating that the 4-compartment model essentially reduces to the 2-compartment regional blood flow model. In case of creatinine, however, the distribution of eigenvalues is more balanced (a factor of 10² between the largest and the smallest eigenvalue) indicating that all four compartments contribute to creatinine kinetics in hemodialysis. Interpretation. Apart from providing an exact analytic solution for practical applications such as the identification of relevant model and treatment parameters, the matrix-based approach reveals characteristic details on model symmetry and complexity for different solutes.

1. Introduction

Compartment models are popular in pharmacokinetics and, as a special application, in hemodialysis, where such models serve to quantify treatment dose [1–3]. Most of that kinetic analysis has been done for urea usually described by 2-compartment models. However, unlike most pharmacokinetic models, the volume of compartments cannot be assumed as constant because of ultrafiltration of excess volume within and accumulation of volume between hemodialysis treatments. The effects on solute concentrations and solute balance caused by volume changes are not negligible. Still, the problem can be expressed as 2-dimensional, inhomogeneous ordinary differential equations (ODE) [4], and the closed form solution to this problem is known. Furthermore, for the variable-volume 2-compartment model urea concentrations have been presented as explicit functions of time and model parameters so that the concentrations in the two compartments at any time can be computed in a single step [5, 6]. That approach was based on the variation of the constants method.

While urea, a solute of little toxicity, is a useful marker of uremia, solutes with limited transfer between compartments such as phosphate, creatinine, 𝛽₂-microglobulin, or glucose are of greater clinical interest [7–11]. Recently, one of us presented a variable-volume 4-compartment model to describe the kinetics of both urea and creatinine using physiological principles of solute transport and only one parameter (membrane permeability) to distinguish between the two solutes [12] (Figure 1). Again, the model can be presented as a 4-dimensional, inhomogeneous ODE, for which the general solution involves integration and which is not defined in the closed form. Thus, the concentrations of the solutes in the four compartments are not obtained without extensive computations.

Figure 1 

Schematic representation of the four-compartment DA-RBF model [12]; note that the blood flow is not exposed; for detailed description see formulas in the Appendix.

It is the purpose of this paper to provide an exact analytical solution for multicompartment kinetic problems, assuming linear changes in volume, and to apply this solution to the recently described 4-compartment model for both urea and creatinine kinetics in hemodialysis.

2. Methods

The 𝑁-compartment model with time dependent volume values, 𝑣_𝑛(𝑡), and concentration of the solute, 𝑐_𝑛(𝑡), can be described by the following set of differential equations:𝑑𝑣𝑑𝑡^𝑛(𝑡)𝑐_𝑛(𝑡)=𝑘_𝑛1𝑐₁(𝑡)+𝑘_𝑛2𝑐₂(𝑡)+⋯+𝑘_𝑛𝑁𝑐_𝑁(𝑡)+𝑔_𝑛,𝑛=1,2,…,𝑁, (1) where coefficients 𝑘_𝑛𝑚 describe the solute exchange or removal in terms of solute clearance, typically in mL/min, and 𝑔_𝑛 models solute mass input to the 𝑛th compartment in terms of generation rate, typically in mg/min. After substitution using 𝑏_𝑛𝑚=𝑘_𝑛𝑚𝑘for𝑚≠𝑛_𝑛𝑛−𝑑𝑣𝑑𝑡^𝑛(𝑡)for𝑚=𝑛, (2) we obtain 𝑣_𝑛(𝑡)𝑑𝑐𝑑𝑡^𝑛(𝑡)=𝑏_𝑛1𝑐₁(𝑡)+⋯+𝑏_𝑛𝑛𝑐_𝑛(𝑡)+⋯+𝑏_𝑛𝑁𝑐_𝑁(𝑡)+𝑔_𝑛,𝑛=1,2,…,𝑁. (3) The set of equations from (3) can be rewritten in matrix form 𝐯(𝑡)∘𝑑𝐜(𝑡)(𝑡)𝑑𝑡=𝐁𝐜+𝐠, (4) where capital symbols are used for matrices, small symbols in bold represent vectors, and ∘ denotes the Hadamard product [13].

Assume that the volume changes are described by linear functions of time: 𝐯(𝑡)=𝑉₀+𝑄𝑡𝐫∶^𝑁_𝑛=1𝑟_𝑛=1; (5) here 𝑉₀ is the total volume of all compartments at time 𝑡=0, 𝑄 refers to the negative ultrafiltration rate, 𝑄<0, during hemodialysis (HD), or to the positive fluid accumulation rate, 𝑄>0, during the interdialytic interval (ID), and 𝑟_𝑛 refers to the fractional volume of each compartment.

Next, let us normalize the volume to the initial value [5]: 𝑤(𝑡)=𝑉₀+𝑄𝑡𝑉₀. (6) Notice that 𝑤(𝑡) is a linear function of 𝑡. Then, after having replaced elements of matrix 𝐁 and vector 𝐠 with the following values 𝑎_𝑛𝑚=𝑏_𝑛𝑚𝑟_𝑛𝑄,𝑏_𝑛=𝑔_𝑛𝑟_𝑛𝑄, (7)

(4) may be rewritten as 𝑉₀𝑄𝑤(𝑡)𝑑𝐜(𝑡)(𝑡)𝑑𝑡=𝐀𝐜+𝐛. (8)

The analytical expressions for the concentrations defined by (8), for 𝑛=1,2,…,𝑁, assuming that matrix 𝐀 has 𝑁 distinct eigenvalues, are in the form 𝑐_𝑛(𝑡)=𝑥_𝑛1𝑤^𝜆1(𝑡)+𝑥_𝑛2𝑤^𝜆2(𝑡)+⋯+𝑥_𝑛𝑁𝑤^𝜆𝑁(𝑡)+𝑑_𝑛, (9) where 𝜆_𝑛 are the eigenvalues of matrix A. The matrix equation equivalent to (9) is given as 𝐜(𝑡)=𝐗𝑤^𝝀(𝑡)+𝐝, (10) where 𝑤^𝝀(𝑡) is a column vector consisting of elements 𝑤^𝜆𝑛(𝑡).

2.1. The Proof

Theorem 1. To prove that for a system described by (8), the solution is provided by (9) or (10).

Lemma 2. Consider the general matrix form of the vector differential equation to be solved as 𝑑𝑝+𝑞𝑡𝐜(𝑡)(𝑡)𝑑𝑡=𝐀𝐜+𝐛,𝑝+𝑞𝑡≠0. (11)
For simplification assume that all N eigenvalues of matrix 𝐀 are distinct.
One will show that, under assumption that 𝐛 belongs to the image of a transform defined by matrix 𝐀 which is true for the considered model, the solution of (11) is in the form of 𝐜(𝑡)=𝑝+𝑞𝑡𝑝^(1/𝑞)𝐀𝐜₀−𝐝+𝐝,where𝐝satisfiescondition∶𝐛=−𝐀𝐝. (12)

Sublemma 3. First, let one show that 𝐀𝑢^𝐀=𝑢^𝐀𝐀. (13)
Recall that the Hermite polynomial for the function of a matrix 𝑓(𝐀) [14], also known as Lagrange-Sylvester formula for the case when all 𝑁 eigenvalues are distinct, is given as 𝑓(𝐀)=^𝑁_𝑛=1𝑓𝜆_𝑛^𝑁_{𝑘=1𝑘≠𝑛}𝐀−𝜆_𝑘𝐈𝜆_𝑛−𝜆_𝑘, (14) where 𝜆_𝑛 is the nth eigenvalue of matrix 𝐀, so that 𝑓(𝜆_𝑛) is the scalar function of a scalar variable. Assuming that 𝑢 is scalar, and taking 𝑓(𝐀)=𝑢^𝐀, one obtains that 𝑢^𝐀 is a sum of scalars multiplied by integer powers of 𝐀, showing that (13) must be true. This observation ends the proof of the sublemma.

With (13) we obtain that the following identity, which results from substituting the right-hand term of (12) into right-hand side of (11), is also true: 𝐀𝐜(𝑡)+𝐛=𝐀𝑝+𝑞𝑡𝑝^(1/𝑞)𝐀𝐜₀+𝑝+𝑞𝑡𝑝^(1/𝑞)𝐀𝐛. (15)

Now, let us find, using (12), the expression for the derivative of 𝐜(𝑡): 𝑑𝐜(𝑡)=1𝑑𝑡𝐀𝑝+𝑞𝑡𝑝+𝑞𝑡𝑝^(1/𝑞)𝐀𝐜₀−𝐝. (16)

To explain the result in (16) let us consider the scalar 𝑢 being 𝑢=𝑢(𝑡)=𝑝+𝑞𝑡𝑝^1/𝑞. (17)

The derivative of the scalar function 𝑢^𝜆(𝑡) results in another function 𝑓(𝜆): 𝑓(𝜆)=𝑑𝑢𝑑𝑡^𝜆(𝑡)=𝜆𝑢(𝑡)𝑢^𝜆(𝑡)𝑑𝑢(𝑡),𝑑𝑑𝑡wherefor𝑝+𝑞𝑡≠0∶𝑢(𝑡)=𝑑𝑑𝑡𝑑𝑡𝑝+𝑞𝑡𝑝^1/𝑞=1𝑞𝑝𝑝+𝑞𝑡𝑝+𝑞𝑡𝑝^1/𝑞𝑑𝑑𝑡𝑝+𝑞𝑡𝑝=1𝑞𝑝𝑢(𝑡)𝑞𝑝+𝑞𝑡𝑝=1𝑢(𝑡).𝑝+𝑞𝑡 (18)

Utilizing again (14), we obtain that the derivative of 𝑢^𝐀(𝑡) must be 𝑓(𝐀)=𝑑𝑢𝑑𝑡^𝐀(𝑡)=1𝑢(𝑡)𝐀𝑢^𝐀(𝑡)𝑑𝑢(𝑡).𝑑𝑡 (19)

Substituting 𝑢(𝑡) from (17) to (19) and using the result from (18) we get (16).

Comparing the right-hand sides of (15) and (16), we obtain that 𝐜(𝑡) given by (12) satisfies (11). If we substitute 𝑡=0 in (12) we obtain that also the initial condition is satisfied.

Expansion from (14) used for 𝑓(𝐀)=𝑢^𝐀 leads also to the identity 𝑢^𝐀𝐱_𝑛=𝐱_𝑛𝑢^𝜆𝑛, (20) where 𝐱_𝑛 is the 𝑛th eigenvector of matrix 𝐀 corresponding to the eigenvalue 𝜆_𝑛. Equation (20) implies that 𝑢^𝐀𝑁_𝑛=1𝐲_𝑛𝑠_𝑛=^𝑁_𝑛=1𝐲_𝑛𝑠_𝑛𝑢^𝜆𝑛, (21) where 𝐲_𝑛 is the normalized 𝑛th eigenvector and 𝑠_𝑛 is the scaling coefficient relating 𝐲_𝑛 to 𝐱_𝑛. As it is always possible to represent the 𝑁-dimensional vector in the relevant basis, we can write the following: 𝐜₀−𝐝=^𝑁_𝑛=1𝐲_𝑛𝑠_𝑛. (22)

Assuming (17) and comparing (21), (22), and (12) we obtain that the analytical solution of (11) is in the form of 𝐜(𝑡)=𝐝+^𝑁_𝑛=1𝐲_𝑛𝑠_𝑛𝑢^𝜆𝑛(𝑡)=𝐲₁𝑠₁𝑢^𝜆1(𝑡)+𝐲₂𝑠₂𝑢^𝜆2(𝑡)+⋯+𝐲_𝑁𝑠_𝑁𝑢^𝜆𝑁(𝑡)+𝐝. (23)

In the discussed dialysis models, where 𝑝+𝑞𝑡=(𝑉₀/𝑄)𝑤(𝑡)=(𝑉₀+𝑄𝑡)/𝑄, we take 𝑞=1 and 𝑝=𝑉₀/𝑄, and in such case the general solution (11) can be reduced to the form 𝐜(𝑡)=𝑤^𝐀(𝑡)𝐜₀−𝐝+𝐝. (24)

Thus, in (23) 𝑢(𝑡) should be replaced with 𝑤(𝑡), which ends the proof of the theorem.

2.2. The Computational Recipe

The practical computational algorithm to compute the coefficients in (9) or in (10) is presented below.

Step 1. Solve the linear equation and find 𝐝: 𝐀𝐝=−𝐛⟹𝐝=𝐀⁻¹(−𝐛). (25)

Step 2. Find the eigenvalues 𝜆_𝑛 of matrix 𝐀 and check if they are distinct. If they are, proceed with Step 3; if not, which should be an extremely rare case, return to the model parameters and introduce small changes in their values; then start from formulating the set of equations.

Step 3. Find the corresponding eigenvectors 𝐲_𝑛, and form the following square matrix 𝐘: 𝐘=𝐲₁𝐲₂⋯𝐲_𝑁. (26)

Step 4. Solve another linear equation and find the set of 𝑁 scaling coefficients 𝑠_𝑛 matching the initial conditions 𝐜₀=𝐜(𝑡=0): 𝐘𝐬=𝐜₀−𝐝⟹𝐬=𝐘⁻¹𝐜₀−𝐝. (27)

Step 5. Compute scaled eigenvectors 𝐱_𝑛 to obtain coefficients 𝑥_𝑛𝑚: 𝐱_𝑛=𝑠_𝑛𝐲_𝑛⟹𝐗=𝐱₁𝐱₂⋯𝐱_𝑁. (28)

After having computed the concentrations for all compartments, in hemodialysis modeling it is quite useful to compute the so-called equilibrated solute concentration, proportionally averaged for all compartments: 𝑐_eq(𝑡)=^𝑁_𝑛=1𝑟_𝑛𝑐_𝑛(𝑡), (29) which represents the overall state of the patient in terms of solute concentration and which is directly related to treatment dose [15–17]. Another concentration of interest is that of the solute in the arterial plasma representing the concentration accessible to direct experimental measurements. See (A.12) in the Appendix for the relevant formula.

The schematic representation of the DA-RBF model, for which the described algorithm was developed, is presented in Figure 1.

The basic computations were performed for two sets of reference model parameters, for urea and creatinine, contained in Table 1.

Table 1 

Initial conditions, times, and model parameters of the DA-RBF model adapted from [12].

To verify whether the assumptions required in the designed algorithm are reasonable and to check whether the observed properties are more general, a comparative study was performed with model parameters randomly varied around the reference values within the physiologically justified range. Pseudo random values were generated for assumed intervals according to uniform distribution. Thus, for 𝑓_𝑝 and 𝑓_pw a 2% radius around the central value was assumed, for 𝑓_ew 5%, for 𝑓_QH and 𝐸_𝑥 10%, and 50% for 𝑄_𝑐, 𝑄_𝑓, 𝑄_𝑥, 𝑉_𝑥, 𝐻, 𝑓_VH, 𝑓_ecv, 𝑉_𝑡, 𝐶₀, and 𝐺₂₄. Two conditions were added to prevent extremely irrelevant cases: (a) 0.01𝑄_𝑥<𝑉_𝑡, (b) 𝑉_uf<0.125𝑉_𝑡. For example, 𝑄_𝑐 ranged from 2.9 L/min to 8.7 L/min. In such a way two sets of 100,000 models were obtained, separately for urea and creatinine.

3. Results

For the purpose of this study it was assumed that the whole treatment cycle was balanced with regard to volume; that is, the total water volume uptake between treatments was the same as the volume removed during ultrafiltration [12].

For the 4-compartment DA-RBF model the set of explicit expressions for matrix 𝐁 is provided in Appendix. Then, the matrix describing the set of differential equations, (8), for the HD and ID periods, respectively, consists only of the constant values and takes the following form: 𝐀_HD=𝑎₁₁𝑎₁₂𝑎₁₃0𝑎₂₁−𝑎₂₁𝑎00₃₁0𝑎₃₃𝑎₁₂00𝑎₂₁−𝑎₂₁,𝐀_ID=𝑎₁₁𝑎₁₂𝑎₁₃0𝑎₂₁−𝑎₂₁𝑎−100₃₁0𝑎₃₃𝑎₁₂00𝑎₂₁−𝑎₂₁.−1 (30)

The resulting values of relevant 𝑎₁₁ to 𝑎₄₄ elements as functions of the model parameters, taken from Table 1, for both hemodialysis and interdialytic intervals and for urea as well as for creatinine are given in Table 2.

Table 2 

Elements 𝑎_nm of A.

The corresponding eigenvalues are summarized in Table 3. Notice that all eigenvalues are negative for the interdialytic phase and positive for the hemodialysis phase. Also notice the range in eigenvalues for different solutes (up to the range of 1 × 10⁶ for urea) and phases of a complete treatment cycle. In the considered cases all the eigenvalues proved to be distinct, which resulted in four different eigenvectors.

Table 3 

Eigenvalues λ of A.

The above observations were confirmed for all comparative, randomly generated models. In particular, in all modeled cases matrix 𝐀 proved to be diagonalizable, and the eigenvalues and eigenvectors were always real (not complex).

Since for 𝑡>0 relative volume 𝑤 is 𝑤<1 during hemodialysis and 𝑤>1 during the interdialytic period, the term 𝑤^𝜆 in (9) is always 0<𝑤^𝜆<1 during both hemodialysis and interdialytic periods because of positive or negative eigenvalues, respectively. However, with very large values of λ such as with urea where 𝜆₃ and 𝜆₄ are in the range of 1 × 10⁶, 𝑤^𝜆≈0. This indicates that the split in intra- and extracellular spaces is ineffective in the case of urea (i.e., the intercompartment clearance is very high) and that urea follows 2-compartment kinetics. Notice the close correspondence of intra- and extracellular urea concentrations (Figure 2(a)). In case of creatinine, however, the ratio of the largest to the smallest eigenvalue is ≈100 during hemodialysis. This indicates that all compartments contribute to overall creatinine kinetics during hemodialysis and somewhat less during the interdialytic phase. Notice the separation of concentrations in all four compartments throughout hemodialysis (Figure 2(b)).

(a)

(b)

Figure 2 

Time course of solute concentrations (urea (a); creatinine (b)) in the four compartments during hemodialysis and during a 100 min postdialytic phase using model parameters from Table 1. Notice that extra- and intracellular concentrations are superimposed in case of urea (a) and separated in case of creatinine (b). 𝑐Li: low-flow intracellular; 𝑐Le: low-flow extracellular; 𝑐Hi: high-flow intracellular; 𝑐He: high-flow extracellular; 𝑐eq: equilibrated over throughout all compartments.

For the 2-compartment model the matrix 𝐀 and the eigenvalues λ for the 2-compartment model have been published previously [5]. A comparison of eigenvalues from the 4- to those determined from the 2-compartment model shows close correspondence for urea. Eigenvalues are in the same order of magnitude such as 𝜆₁ for the intradialytic phase with values of 15.4 and 12.3, respectively (Table 3). The small difference originates from differences in parameters assumed in published 2- and 4-compartment models which were obtained in different studies [5, 12]. When the fraction of extracellular volume was assumed as close to 1, equivalent to eliminating the effect of intracellular sequestration in the 4-compartment model, the eigenvalues for the urea model remained essentially unchanged.

Figure 3 presents the sensitivities, expressed as absolute values, of two crucial eigenvalues and related coefficients that may have impact on the shapes of the concentration profiles when the values presented in Table 1 are varied by a small fraction. Notice that changes in 𝐸_𝑥 directly related to dialysis also affect the runs during the ID phase because the concentrations at the end of the HD phase serve as starting points for the ID stage.

(a)

(b)

Figure 3 

Sensitivity of computed coefficients 𝑥, 𝑑, and λ, with regard to model parameters for urea (upper bars in light red) and creatinine (lower bars in dark blue) models. Normalized sensitivities are indicated on the horizontal axis.

The equilibrated concentration (29) representing the weighted average of compartmental concentrations shows the dominating impact of the first two eigenvalues. Table 4 presents the relative error of 𝑐_eq(𝑡) computed at the beginning of HD for the set of 100,000 simulations using only two eigenvalues and corresponding coefficients compared to that computed with the complete set of eigenvalues. Notice the negative error with omission of positive values.

Table 4 

Relative error for 𝑐_eq, in %, for the creatinine model during the HD phase, when 𝑐_eq was computed using only two eigenvalues and corresponding coefficients using parameters from 100,000 simulations.

4. Discussion

In this paper a general analytical solution for a particular class of variable-volume multicompartment solute kinetic models is presented. The solution is based on a matrix approach applicable to linear models encompassing any number of compartments assuming that all eigenvalues of the matrix are distinct. The presented solution is based on a finite volume change typical for hemodialysis. The solution in absence of a volume change requires a different approach the discussion of which is beyond the scope of this paper. The detailed procedure is also provided.

One purpose of mathematical modeling is to characterize the system by its structure such as the number of compartments and their interaction and by its parameters such as the distribution volumes and rate constants [18]. Many parameters are usually inaccessible to direct measurement and must be obtained by parameter identification, that is, by fitting an appropriate model output to observable experimental data. This procedure involves recurrent numerical solution of the ODE for different sets of parameter values until a chosen set of parameters provides the best fit. While parameter identification in more complex models is a problem of its own [19], the numerical solution of a given ODE for each parameter set to be evaluated in the process of parameter identification makes this task very laborious and time consuming, even for powerful personal computers. It therefore pays to replace numerical for exact analytical procedures wherever possible.

The solution developed in this study was applied to the variable-volume 2-compartment model for urea kinetics as well as to the variable-volume 4-compartment model for urea and creatinine kinetics presented elsewhere [5, 12]. The latter has also been used to describe the kinetics of 𝛽₂-microglobulin [10].

Apart from its mathematical use, the qualitative examination of matrices and eigenvalues provides a good means to compare models and to judge the effective number of compartments. Based on the eigenvalues ranging from 12.3 to 1.1 × 10⁶ for the hemodialysis interval, the 4-compartment urea kinetic model essentially is a 2-compartment urea kinetic model. In the case of creatinine, however, the eigenvalues range from 7.5 to 725.9 for the hemodialysis interval and are much closer to each other so that the 4-compartment structure is justified (Figure 2). The models with varied parameters indicated comparable properties.

Results presented in Table 4 confirm that computations can considerably be simplified for experimental applications. Reduction of parameters to compute 𝑐_eq, important from the diagnostic and planning point of view, leading to just two eigenvalues with corresponding coefficients, should provide acceptable accuracy.

The sensitivity study shows that in clinical practice, where collection of the complete set of personalized model parameters is not always possible, the most attention should be paid to an accurate estimation of 𝑓_QH, 𝑄_𝑐, 𝑉_𝑡, 𝑉_uf, 𝑄_𝑥, 𝐸_𝑥, 𝐶_𝑜, 𝐺, 𝑓_pw, and  𝑓_VH. However, limited sensitivity to other parameters indicates that uncertainty in such parameters does not considerably change the modeling and treatment outcome.

Analytical solutions for the variable-volume 2-compartment model have been presented before. In the approach presented by Grandi et al. the system of first-order linear differential equations comparable to that given in (4) was transformed into a type of Euler-Cauchy second-order differential equation and solved analytically [6]. An expansion of this method to more dimensions has not been provided.

A matrix-based strategy to solve the variable-volume 2-compartment model has been presented by one of the authors [5]. This approach has been applied to two physiologically distinct interpretations of urea kinetics, either assuming diffusion-limited transfer between intra- and extracellular compartments [20] comparable to that of Grandi et al. [6] or flow-limited transfer between poorly and highly perfused organ systems [5, 21].

In these variable-volume models, fluid is proportionally removed from both compartments. While this assumption is compatible with the flow-limited interpretation, the assumption is at odds with clinical and physiological understanding in the diffusion-limited interpretation of 2-compartment urea kinetics. In the classic diffusion-limited model the two compartments refer to intra- and extracellular volumes, and since excess fluid is best removed under close to isotonic conditions during hemodialysis [22–24] and accumulated under isotonic conditions by matched ingestion of salt and water between treatments [25], fluid is more or less exclusively removed from and added to the extracellular compartment, leaving the volume of the intracellular compartment unchanged. The classic variable-volume 2-compartment urea kinetic model therefore assumes a constant intracellular volume during hemodialysis and between dialysis treatments [1]. To account for this effect Smye and Will allocated 95% of the volume changes to the extracellular compartment and provided an approximated solution for this model [26].

For solutes actually sequestered in the intracellular space because of limited membrane permeability, transport is both diffusion- and flow-limited. The combination of both diffusion- and flow-limited transport characteristics therefore leads to a 4-compartment model presented earlier [12]. The equations can be analyzed using the matrix method.

The matrix-based strategy presented in this paper is akin to that presented in [5] albeit distinct from the previous method regarding an important technical point. As in the previous approach, the solution requires the computation of the eigenvalues of 𝐀. For the 2-compartment model this refers to solving a quadratic equation, for the 3-compartment model to solving a cubic equation according to Cardano [27, 28], and for the 4-compartment model to solving a quartic equation according to Tignol and Ferrari [28, 29]. For models involving more than four compartments the eigenvalues have to be found by numeric approximation using dedicated software such as Matlab (The MathWorks Inc., Natick, Massachusetts, USA). As in the previous approach, the solution requires the computation of 𝐀⁻¹, the inverse of matrix 𝐀. The manual inversion of a 2 by 2 matrix is easily done for the 2-compartment model but prohibitively tedious for matrices of higher dimensions. Dedicated software such as Matlab or internet resources can be used for that problem. The computational tool, to compute the coefficients and concentrations, developed in Java, may be accessed at [30]. An important difference to the previous approach, however, is the absence of computing the definite integral (the integration of the inverse transformation matrix from the start of the observation phase to the point of interest, equation C15 in [5]) over the time course. If the integral is not solved analytically as in [5] the integration has to be carried out numerically, defeating the purpose of an exact solution and consuming considerable computational time.

5. Conclusions

In conclusion, a closed-form solution for the variable-volume 𝑁-compartment model is presented. The solution can be applied to the variable-volume 4-compartment diffusion-adjusted regional blood flow model. Even though the complexity of models and mathematics considerably increases in the transition from 2- to 4-compartment models, the availability of an exact solution should help in practical applications such as the identification of relevant model and treatment parameters in hemodialysis.

Appendix

The complete set of formulas for the DA-RBF model is provided below. See Table 1 for description of the parameters.

To simplify the notation, auxiliary parameters are introduced for the volume ratios: 𝑟₁=𝑓_VH𝑓_ecv,𝑟₂=𝑓_VH1−𝑓_ecv,𝑟₃=1−𝑓_VH𝑓_ecv,𝑟₄=1−𝑓_VH1−𝑓_ecv. (A.1)

The volume flows during hemodialysis (HD) and interdialytic (ID) intervals are defined as 𝑄_uf=𝑉_uf𝑡_𝑑,𝑄_in=𝑉_uf𝑡_di−𝑡_𝑑,then𝑄=−𝑄_uf𝑄forHD,_inforID. (A.2)

The regional blood flows, according to [12], are 𝑄_𝑠=𝑄_𝑐−𝑄_𝑓,𝑄_𝐻=𝑄_𝑠𝑓_QH,𝑄_𝐿=𝑄_𝑠1−𝑓_QH. (A.3)

Characteristic times for extracorporeal, regional, and transmembrane flows according to [12] are 𝜏_𝑥=𝑉_𝑥𝑄_𝑥,𝜏_𝐻=𝑉_𝑏𝑓_VH𝑄_𝐻,𝜏_𝐿=𝑉_𝑏1−𝑓_𝑉𝐻𝑄_𝐿,𝑘𝜅=^𝑠𝑓_ew/𝑓_pw(𝐻/(1−𝐻))+1𝑓_ew, (A.4) and the diffusion volume flow fractions for regional and extracorporeal blood flows are given as 𝑓_𝑗𝑥=1−𝐻𝑓_𝑝𝑓_pw+𝐻𝑓_𝑝1−exp−𝜅𝜏_𝑥𝑓_ew,𝑓_𝑗𝐻=(1−𝐻)𝑓_pw+𝐻1−exp−𝜅𝜏_𝐻𝑓_ew,𝑓_𝑗𝐿=(1−𝐻)𝑓_pw+𝐻1−exp−𝜅𝜏_𝐿𝑓_ew,𝑓_𝑗𝑣=𝑓_QH𝑓_𝑗𝐻+1−𝑓_QH𝑓_𝑗𝐿,𝑓_𝑗𝑎=𝑓_𝑗𝑣𝑄_𝑠+𝑓_𝑗𝑥𝑄_𝑥𝑄_𝑠+𝑄_𝑥,𝑓_𝑗𝑓=𝑓_𝑗𝑥𝑄_𝑥𝑄_𝑓+𝑓_𝑗𝑎𝑄1−^𝑥𝑄_𝑓. (A.5)

The diffusion volume and convective flow rates for HD and ID intervals are computed as 𝐹_𝐻=𝑄_𝐻𝑓_𝑗𝐻,𝐹_𝐿=𝑄_𝐿𝑓_𝑗𝐿,𝐹_𝑐=𝑄_𝑐𝑓_𝑗𝑎,𝐹_𝑓=𝑄_𝑓𝑓_𝑗𝑓,𝐹_uf𝐻=𝑄_uf𝑓_VH,𝐹_uf𝐿=𝑄_uf1−𝑓_VH,𝐹_ufHi=𝑄_uf𝑓_VH1−𝑓_ecv,𝐹_ufLi=𝑄_uf1−𝑓_VH1−𝑓_ecv,𝐹_in𝐻=−𝑄_in𝑓_VH,𝐹_in𝐿=−𝑄_in1−𝑓_VH,𝐹_inHe=−𝑄_in𝑓_VH𝑓_ecv,𝐹_inLe=−𝑄_in1−𝑓_VH𝑓_ecv,𝐹_inHi=−𝑄_in𝑓_VH1−𝑓_ecv,𝐹_inLi=−𝑄_in1−𝑓_VH1−𝑓_ecv (A.6) and intercompartment and dialyzer clearances as 𝐾_𝑐=𝑘_𝑠1−𝑓_ecv𝐾𝑉,_𝐻=𝐾_𝑐𝑓_VH,𝐾_𝐿=𝐾_𝑐1−𝑓_VH,𝐾_𝑑=𝑄_𝑥𝑓_𝑗𝑥𝐸_𝑥+𝑄_uf1−𝐸_𝑥, (A.7) where 𝑉 is the total volume at the end of HD (Table 1).

Then, the elements of 𝐁 (4) for the hemodialysis interval are given as 𝑏₁₁=−𝐹_𝐻−𝐾_𝐻+𝐹_𝐻𝐹_𝐻+𝐹_uf𝐻𝐹_𝑐−𝐹_𝑓+𝐾_𝑑,𝑏₁₂=𝐾_𝐻,𝑏₁₃=𝐹_𝐻𝐹_𝐿+𝐹_uf𝐿𝐹_𝑐−𝐹_𝑓+𝐾_𝑑,𝑏₂₁=𝐾_𝐻−𝐹_ufHi,𝑏₂₂=−𝐾_𝐻+𝐹_ufHi,𝑏₃₁=𝐹_𝐿𝐹_𝐻+𝐹_uf𝐻𝐹_𝑐−𝐹_𝑓+𝐾_𝑑,𝑏₃₃=−𝐹_𝐿−𝐾_𝐿+𝐹_𝐿𝐹_𝐿+𝐹_uf𝐿𝐹_𝑐−𝐹_𝑓+𝐾_𝑑,𝑏₃₄=𝐾_𝐿,𝑏₄₃=𝐾_𝐿−𝐹_ufLi,𝑏₄₄=−𝐾_𝐿+𝐹_ufLi, (A.8) while for the ID interval they are given as 𝑏₁₁=−𝐹_𝐻−𝐾_𝐻+𝐹_𝐻𝐹_𝐻𝐹_𝑐−𝐹_𝑓+𝐹_inHe,𝑏₁₂=𝐾_𝐻,𝑏₁₃=𝐹_𝐻𝐹_𝐿𝐹_𝑐−𝐹_𝑓,𝑏₂₁=𝐾_𝐻,𝑏₂₂=−𝐾_𝐻+𝐹_inHi,𝑏₃₁=𝐹_𝐿𝐹_𝐻𝐹_𝑐−𝐹_𝑓,𝑏₃₃=−𝐹_𝐿−𝐾_𝐿+𝐹_𝐿𝐹_𝐿𝐹_𝑐−𝐹_𝑓+𝐹_inLe,𝑏₃₄=𝐾_𝐿,𝑏₄₃=𝐾_𝐿,𝑏₄₄=−𝐾_𝐿+𝐹_inLi. (A.9)

The generation rate, expressed in mg/mL, in the four compartments is defined as 𝑔₁𝑔=0,₂=𝐺_urea,𝑔₃𝑔=0,₄=0 (A.10) for urea and as 𝑔₁𝑔=0,₂𝑔=0,₃𝑔=0,₄=𝐺_creatinine (A.11) for creatinine, respectively.

Note that in the above formulas indices 1 through 4 refer to extracellular high-flow (He), intracellular high-flow (Hi), extracellular low-flow (Le), and intracellular low-flow (Li) compartments, respectively (also see Figure 1). Therefore, the accessible arterial plasma concentration for the intradialytic period is given as 𝑐_𝑝,art(𝑡)=𝑐_He(𝑡)𝐹_𝐻+𝐹_uf𝐻+𝑐_Le(𝑡)𝐹_𝐿+𝐹_uf𝐿𝐹_𝑐−𝐹_𝑓+𝐾_𝑑. (A.12)

For the interdialytic phase: 𝐾_𝑑=𝐹_uf𝐻=𝐹_uf𝐿=0.

Conflict of Interests

The authors declare that they have no conflict of interests.

Acknowledgment

The help of Dawid Worek in developing the stand alone Java calculator [30] is gratefully acknowledged.