Research Article  Open Access
Analytical Solution of Multicompartment Solute Kinetics for Hemodialysis
Abstract
Objective. To provide an exact solution for variablevolume multicompartment kinetic models with linear volume change, and to apply this solution to a 4compartment diffusionadjusted regional blood flow model for both urea and creatinine kinetics in hemodialysis. Methods. A matrixbased 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^{5} times larger than the other eigenvalues indicating that the 4compartment model essentially reduces to the 2compartment regional blood flow model. In case of creatinine, however, the distribution of eigenvalues is more balanced (a factor of 10^{2} 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 matrixbased 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 2compartment 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 2dimensional, inhomogeneous ordinary differential equations (ODE) [4], and the closed form solution to this problem is known. Furthermore, for the variablevolume 2compartment 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 variablevolume 4compartment 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 4dimensional, 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.
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 4compartment 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: 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 we obtain The set of equations from (3) can be rewritten in matrix form 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: here is the total volume of all compartments at time , refers to the negative ultrafiltration rate, , during hemodialysis (HD), or to the positive fluid accumulation rate, , 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]: Notice that is a linear function of . Then, after having replaced elements of matrix and vector with the following values
(4) may be rewritten as
The analytical expressions for the concentrations defined by (8), for , assuming that matrix has distinct eigenvalues, are in the form where are the eigenvalues of matrix A. The matrix equation equivalent to (9) is given as 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
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
Sublemma 3. First, let one show that
Recall that the Hermite polynomial for the function of a matrix [14], also known as LagrangeSylvester formula for the case when all eigenvalues are distinct, is given as
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 righthand term of (12) into righthand side of (11), is also true:
Now, let us find, using (12), the expression for the derivative of :
To explain the result in (16) let us consider the scalar being
The derivative of the scalar function results in another function :
Utilizing again (14), we obtain that the derivative of must be
Substituting from (17) to (19) and using the result from (18) we get (16).
Comparing the righthand sides of (15) and (16), we obtain that given by (12) satisfies (11). If we substitute in (12) we obtain that also the initial condition is satisfied.
Expansion from (14) used for leads also to the identity where is the th eigenvector of matrix corresponding to the eigenvalue . Equation (20) implies that 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:
Assuming (17) and comparing (21), (22), and (12) we obtain that the analytical solution of (11) is in the form of
In the discussed dialysis models, where , we take and , and in such case the general solution (11) can be reduced to the form
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 :
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 :
Step 4. Solve another linear equation and find the set of scaling coefficients matching the initial conditions :
Step 5. Compute scaled eigenvectors to obtain coefficients :
After having computed the concentrations for all compartments, in hemodialysis modeling it is quite useful to compute the socalled equilibrated solute concentration, proportionally averaged for all compartments: 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 in the Appendix for the relevant formula.
The schematic representation of the DARBF 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.

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 a 2% radius around the central value was assumed, for 5%, for and 10%, and 50% for , , , , , , , , , and . Two conditions were added to prevent extremely irrelevant cases: (a) 0.01, (b) . 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 4compartment DARBF 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:
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.
 
: th row and th column element of ; 2: 2compartment model; 4: 4compartment model; HD: hemodialysis interval; ID: interdialytic interval; u: urea; c: creatinine; ^{†}values for 2compartment model taken from [5]. 
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^{6} 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.
 
: eigenvalue; 2: 2compartment model; 4: 4compartment model; HD: hemodialysis interval; ID: interdialytic interval; u: urea; c: creatinine; ^{†}values for 2compartment model taken from [5]; *4compartment model evaluated for . 
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 relative volume is during hemodialysis and during the interdialytic period, the term in (9) is always 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^{6}, . 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 2compartment 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)
For the 2compartment model the matrix and the eigenvalues λ for the 2compartment model have been published previously [5]. A comparison of eigenvalues from the 4 to those determined from the 2compartment 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 4compartment 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 4compartment 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)
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 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.
 
: equilibrated concentration; HD: hemodialysis interval; ID: interdialytic interval. For urea all values were below , and they remained small during ID for both solutes. 
4. Discussion
In this paper a general analytical solution for a particular class of variablevolume 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 variablevolume 2compartment model for urea kinetics as well as to the variablevolume 4compartment 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^{6} for the hemodialysis interval, the 4compartment urea kinetic model essentially is a 2compartment 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 4compartment 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 , 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 , , , , , , , , , and . 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 variablevolume 2compartment model have been presented before. In the approach presented by Grandi et al. the system of firstorder linear differential equations comparable to that given in (4) was transformed into a type of EulerCauchy secondorder differential equation and solved analytically [6]. An expansion of this method to more dimensions has not been provided.
A matrixbased strategy to solve the variablevolume 2compartment 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 diffusionlimited transfer between intra and extracellular compartments [20] comparable to that of Grandi et al. [6] or flowlimited transfer between poorly and highly perfused organ systems [5, 21].
In these variablevolume models, fluid is proportionally removed from both compartments. While this assumption is compatible with the flowlimited interpretation, the assumption is at odds with clinical and physiological understanding in the diffusionlimited interpretation of 2compartment urea kinetics. In the classic diffusionlimited 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 variablevolume 2compartment 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 flowlimited. The combination of both diffusion and flowlimited transport characteristics therefore leads to a 4compartment model presented earlier [12]. The equations can be analyzed using the matrix method.
The matrixbased 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 2compartment model this refers to solving a quadratic equation, for the 3compartment model to solving a cubic equation according to Cardano [27, 28], and for the 4compartment 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 2compartment 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 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 closedform solution for the variablevolume compartment model is presented. The solution can be applied to the variablevolume 4compartment diffusionadjusted regional blood flow model. Even though the complexity of models and mathematics considerably increases in the transition from 2 to 4compartment 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 DARBF model is provided below. See Table 1 for description of the parameters.
To simplify the notation, auxiliary parameters are introduced for the volume ratios:
The volume flows during hemodialysis (HD) and interdialytic (ID) intervals are defined as
The regional blood flows, according to [12], are
Characteristic times for extracorporeal, regional, and transmembrane flows according to [12] are and the diffusion volume flow fractions for regional and extracorporeal blood flows are given as
The diffusion volume and convective flow rates for HD and ID intervals are computed as and intercompartment and dialyzer clearances as where is the total volume at the end of HD (Table 1).
Then, the elements of (4) for the hemodialysis interval are given as while for the ID interval they are given as
The generation rate, expressed in mg/mL, in the four compartments is defined as for urea and as for creatinine, respectively.
Note that in the above formulas indices 1 through 4 refer to extracellular highflow (He), intracellular highflow (Hi), extracellular lowflow (Le), and intracellular lowflow (Li) compartments, respectively (also see Figure 1). Therefore, the accessible arterial plasma concentration for the intradialytic period is given as
For the interdialytic phase: .
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.
References
 T. A. Depner, Prescribing Hemodialysis: A Guide to Urea Modeling, Kluwer Academic Publishers, Boston, Mass, USA, 1991.
 F. A. Gotch and M. L. Keen, “Kinetic modeling in hemodialysis,” in Clinical Dialysis, A. R. Nissenson and R. N. Fine, Eds., pp. 153–202, McGrrawHill, New York, NY, USA, 4th edition, 2005. View at: Google Scholar
 J. Waniewski, M. Debowska, and B. Lindholm, “Theoretical and numerical analysis of different adequacy indices for hemodialysis and peritoneal dialysis,” Blood Purification, vol. 24, no. 4, pp. 355–366, 2006. View at: Publisher Site  Google Scholar
 P. Bugl, Differential Equations: Matrices and Models, Prentice Hall, Englewood Cliffs, NJ, USA, 1995.
 D. Schneditz and J. T. Daugirdas, “Formal analytical solution to a regional blood flow and diffusion based urea kinetic model,” ASAIO Journal, vol. 40, no. 3, pp. M667–M673, 1994. View at: Google Scholar
 F. Grandi, G. Avanzolini, and A. Cappello, “Analytic solution of the variablevolume doublepool urea kinetics model applied to parameter estimation in hemodialysis,” Computers in Biology and Medicine, vol. 25, no. 6, pp. 505–518, 1995. View at: Publisher Site  Google Scholar
 S. Eloot, A. Torremans, R. De Smet et al., “Complex compartmental behavior of small watersoluble uremic retention solutes: evaluation by direct measurements in plasma and erythrocytes,” American Journal of Kidney Diseases, vol. 50, no. 2, pp. 279–288, 2007. View at: Publisher Site  Google Scholar
 D. Schneditz, Y. Yang, G. Christopoulos, and J. Kellner, “Rate of creatinine equilibration in whole blood,” Hemodialysis International, vol. 13, no. 2, pp. 215–221, 2009. View at: Publisher Site  Google Scholar
 F. Gotch, N. W. Levin, and P. Kotanko, “Calcium balance in dialysis is best managed by adjusting dialysate calcium guided by kinetic modeling of the interrelationship between calcium intake, dose of vitamin D analogues and the dialysate calcium concentration,” Blood Purification, vol. 29, no. 2, pp. 163–176, 2010. View at: Publisher Site  Google Scholar
 V. Maheshwari, L. Samavedham, and G. P. Rangaiah, “A regional blood flow model for β2microglobulin kinetics and for simulating intradialytic exercise effect,” Annals of Biomedical Engineering, vol. 39, no. 12, pp. 2879–2890, 2011. View at: Publisher Site  Google Scholar
 D. Schneditz, M. Galach, K. Thomaseth, and J. Waniewski, “A regional blood flow model for glucose and insulin kinetics during hemodialysis,” ASAIO Journal, vol. 59, no. 6, pp. 627–635, 2013. View at: Google Scholar
 D. Schneditz, D. Platzer, and J. T. Daugirdas, “A diffusionadjusted regional blood flow model to predict solute kinetics during haemodialysis,” Nephrology Dialysis Transplantation, vol. 24, no. 7, pp. 2218–2224, 2009. View at: Publisher Site  Google Scholar
 J. R. Magnus and H. Neudecker, Matrix Differential Calculus with Applications in Statistics and Econometrics, John Wiley & Sons, Chichester, UK, 3rd edition, 2007.
 N. J. Higham, Functions of Matrices: Theory and Computation, SIAM, Philadephia, Pa, USA, 2008.
 S. L. Goldstein, J. M. Sorof, and E. D. Brewer, “Evaluation and prediction of urea rebound and equilibrated Kt/V in the pediatric hemodialysis population,” American Journal of Kidney Diseases, vol. 34, no. 1, pp. 49–54, 1999. View at: Google Scholar
 M. Debowska, B. Lindholm, and J. Waniewski, “Adequacy indices for dialysis in acute renal failure: kinetic modeling,” Artificial Organs, vol. 34, no. 5, pp. 412–419, 2010. View at: Publisher Site  Google Scholar
 A. Jung, P. Korohoda, P. Krisper, and D. Schneditz, “Relationship between kinetics of albuminbound bilirubin and watersoluble urea in extracorporeal blood purification,” Nephrology Dialysis Transplantation, vol. 27, no. 3, pp. 1200–1206, 2012. View at: Google Scholar
 J. T. Daugirdas, T. A. Depner, T. Greene, and P. Silisteanu, “Solutesolver: a webbased tool for modeling urea kinetics for a broad range of hemodialysis schedules in multiple patients,” American Journal of Kidney Diseases, vol. 54, no. 5, pp. 798–809, 2009. View at: Publisher Site  Google Scholar
 G. Lillacci and M. Khammash, “Parameter estimation and model selection in computational biology,” PLoS Computational Biology, vol. 6, no. 3, Article ID e1000696, 2010. View at: Publisher Site  Google Scholar
 D. Schneditz, B. Fariyike, R. Osheroff, and N. W. Levin, “Is intercompartmental urea clearance during hemodialysis a perfusion term? A comparison of two pool urea kinetic models,” Journal of the American Society of Nephrology, vol. 6, no. 5, pp. 1360–1370, 1995. View at: Google Scholar
 D. Schneditz, J. C. Van Stone, and J. T. Daugirdas, “A regional blood circulation alternative to inseries two compartment urea kinetic modeling,” ASAIO Journal, vol. 39, no. 3, pp. M573–M577, 1993. View at: Publisher Site  Google Scholar
 S. F. F. Santos and A. J. Peixoto, “Revisiting the dialysate sodium prescription as a tool for better blood pressure and interdialytic weight gain management in hemodialysis patients,” Clinical Journal of the American Society of Nephrology, vol. 3, no. 2, pp. 522–530, 2008. View at: Publisher Site  Google Scholar
 F. Locatelli, A. Covic, C. Chazot, K. Leunissen, J. Luño, and M. Yaqoob, “Optimal composition of the dialysate, with emphasis on its influence on blood pressure,” Nephrology Dialysis Transplantation, vol. 19, no. 4, pp. 785–796, 2004. View at: Publisher Site  Google Scholar
 J. Raimann, L. Liu, S. Tyagi, N. W. Levin, and P. Kotanko, “A fresh look at dry weight,” Hemodialysis International, vol. 12, no. 4, pp. 395–405, 2008. View at: Publisher Site  Google Scholar
 A. J. Peixoto, N. Gowda, C. R. Parikh, and S. F. F. Santos, “Longterm stability of serum sodium in hemodialysis patients,” Blood Purification, vol. 29, no. 3, pp. 264–267, 2010. View at: Publisher Site  Google Scholar
 S. W. Smye and E. J. Will, “A mathematical analysis of a twocompartment model of urea kinetics,” Physics in Medicine and Biology, vol. 40, no. 12, article 001, pp. 2005–2014, 1995. View at: Publisher Site  Google Scholar
 G. Cardano, Book Number One About the Great Art, or the Rules of Algebra, Nürnberg, Germany, 1545.
 J. P. Tignol, Galois' Theory of Algebraic Equations, World Scientific, London, UK, 2001.
 L. Ferrari, Book Number One About the Great Art, or the Rules of Algebra, G. Cardano, Ed., Nürnberg, Germany, 1545.
 http://www.galaxy.agh.edu.pl/~korohoda/applet_for_DARBF_model/MedCalApp.html.
Copyright
Copyright © 2013 Przemysław Korohoda and Daniel Schneditz. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.