Application of dynamic optimisation for planning a haemodialysis process

Background The aim of the study is to show that optimization tools can be used in planning the haemodialysis process in order to obtain the most effective treatment aimed at removing both urea and phosphorus. To this end we use the IV–compartment model of phosphorus kinetics. Methods The use of the IV–compartment model of phosphorus kinetics forces us to apply new numerical tools which cope with a rebound phenomenon that can occur during haemodialysis. The proposed algorithm solves optimization problems with various constraints imposed on concentrations of urea and phosphorus. Results We show that the optimization tools are effective in planning haemodialysis processes aimed at achieving desired levels of urea and phosphorus concentrations at the end of these processes. One of the numerical experiments reported in the paper concerns patients data who experienced a rebound phenomenon during haemodialysis due to a low level of phosphorus concentration. Conclusion In order to plan haemodialysis processes one should take into account the fact that these processes, in general, are described by different equations in different regions determined by phosphorus concentrations. This follows from the fact that mechanisms modelled by IV–compartment model are activated during dialysis. Therefore, advanced numerical tools have to be used in order to simulate and optimize these processes. The paper shows that these tools can be constructed and effectively applied in planning haemodialysis processes.


Background
Hiperphosphatemia is associated with increased mortality among dialysis patients primarily due to accelerated cardiovascular disease (CVD) [1,2]. Altered calcium, phosphorus and PTH levels which accompany chronic kidney disease, as well as an excessive burden of calcium supplied by calcium-based oral phosphorus binders are responsible for exaggerated vascular calcification and lead to enhanced atherosclerosis in this population [3].
European guidelines recommend lowering the phosphorus level in the blood towards the normal range [4]. Serum phosphorus concentration in haemodialysis patients derives from phosphorus content in the diet, its elimination with residual renal function and its removal *Correspondence: Wojciech.Stecz@wat.edu.pl 1 Military University of Technology, Faculty of Cybernetics, Kaliskiego 2, 00-908 and Warsaw, Poland Full list of author information is available at the end of the article during haemodialysis (HD) sessions. Therefore, therapeutic approaches focus on two directions: the usage of oral binders which decreases phosphorus absorption from the intestinum, as well as effective phosphorus elimination by different modes of haemodialysis treatment. The problem for establishing the most efficient schedule of haemodialysis treatment is that the kinetics of phosphorus is much more complicated than other molecules such as urea or creatinine. During the first hour of an haemodialysis session, phosphorus concentration rapidly decreases, followed by a relatively stable level (plateau) until the end of haemodialysis. After haemodialysis the phosphorus level increases rapidly reaching the predialytic level within 4-8 h [5]. Therefore the most effective strategy is to extend weekly dialysis time by enhancing the length of the sessions and/or by increasing the frequency of the sessions. It is estimated that the weekly dialysis time required to avoid the usage of oral binders is 18-30 h. Other methods such as increasing blood flow rate, greater dialyzer size, increased dialysate flow rate are also used to increase phosphorus elimination however they are less effective [6].
Many mathematical models, from one to multicompartment models, have been applied to describe phosphorus kinetics ( [7][8][9][10][11]). During haemodialysis session phosphorus is removed from the plasma which is a compartment accessible for the treatment mode. The intradialytic plateau of phosphorus level suggests the existence of second phosphorus storage which is inaccessible for the dialyser and provides a continuous inflow of phosphorus resulting in the plateau of its level in the second part of haemodialysis. Therefore a two-pool model was proposed. However, it did not completely explain the constant phosphorus level during haemodialysis. Therefore Spalding et al. proposed a four-compartment model which suggested gradual phosphorus mobilization from the pools during dialysis [11]. Gotch et al proposed intracellular volume as a site of phosphorus storage [12]. Another conception of phosphorus kinetics assumes a pseudo I-compartment model where the phosphorus level during haemodialysis is determined by its clearance from distal compartment. The problem is that the rate of this clearance is dependent on the patient and clinical circumstances [13].
The aim of the study was to introduce advanced numerical tools in planning the haemodialysis process in order to obtain the most effective treatment in removing urea and phosphorus as well. The mathematical model assumes the modified IV-compartment model of phosphorus kinetics.
The organisation of the paper is as follows. "Methods" section gives the introduction to the haemodialysis problem. We concentrate on the mathematical models for urea and phosphorus kinetics. Next we present an optimisation procedure that can be used for haemodialysis. In "Results" section we present some simulation and optimisation experiments. Last section discusses our results and possible future research directions.

Models of urea kinetics
The II-compartment model depicts a patient as a twocompartment system of different volumes of solutes that are regularly cleared, i.e., intracellular fluid and extracellular fluid. Between these solutes, diffusion of substances takes place according to different substances' concentrations. II-compartment model is shown in Fig. 1.
According to [11,14] the equations describing IIcompartment model are as follows: where: UFR -ultrafiltration volume U ufr -profile of the ultrafiltration rate K ufr -ultrafiltration coefficient K urea r -urea clearance associated with residual renal function K urea D -dialyser clearance dC urea EC dt -variation of urea concentration in the extracellular fluid dC urea IC dt -variation of urea concentration in the intracellular fluid K urea IE -diffusion speed between inner and outer cellular solute through a cellular membrane -the same in both directions G urea -urea generation in the human cells according to quick variations of ions during haemodialysis -unsurveyed process, but was observed -we set G urea to zero in our simulation and optimisation tests.
In the above equations we use the Watson formula ( [15]) of intracellular to extracellular volume ratio of 2:1 (0.66 · V (0) for intracellular volume and 0.34 · V (0) for extracellular volume).

Models of phosphorus kinetics
In the paper we rely on the models described in the papers ( [7,11,16]) presenting the phosphorus kinetics. We use especially Spalding's IV-compartment model shown in Fig. 2.
As far as phosphorus kinetics is concerned we have employed the model presented in ( [7,11]): where:  Additionally, on the basis of the analysis given in [11], we model the qualitative change in equations (4)-(5) by the following algebraic equations (thereby, among other things, functions stated in [11] are continuous):  is activated only when a phosphorus level is below C PO 4 min respectively (fourth compartment).
Furthermore, we underline that a phosphorus clearance coefficient K PO 4 D is a function of two parameters: predefined by producer phosphorus clearance of the dialyser membrane K PO 4 D,base , and the clearance reduction coefficient κ for the patient (in literature, one may see that from some group of patients there are suggested values of κ [7,17], we assume value κ = 0.6 for hollow fiber dialyzers we used).
To check our model we used data of two patients such as: age, gender, body mass, height, serum urea and phosphorus levels, the dialyzers size, value of blood and dialysate flow during dialysis session, ultrafiltration volume. On the basis of the data models parameters were determined for each patient by simulating equations (1)-(8) (see Table 1). Hollow fiber dialyzers made by Allmed (Polypure series) or Gambro (Polyflux series) were used.
The study protocol was accepted by the local ethical committee. Each participant signed the informed consent.

Haemodialysis optimisation problem
The main optimisation problem considered in the paper is stated as follows. Having combined kinetic models of urea and phosphorus we look for proper concentrations of urea and phosphorus and the end of the haemodialysis process by controlling the parameters Q B , Q D , U ufr . In other words, by solving the optimal control problem we want to choose a proper dialysis membrane in order to achieve final parameters of haemodialysis.
That optimisation problem can be formulated as follows: subject to the constraints (1)-(8), the following con- and the constraints on the control variables Table 1 Parameters for simulation of haemodialysis process Input data -first patient parameters Input data -second patient parameters where L urea EC and L urea IC mean maximum admissible values for urea concentrations. L UFR min and L UFR max mean minimum and maximum admissible values for the ultrafiltration.
ufr define bound constraints on decision variables, in particular U min ufr , U max ufr limit the rate of ultrafiltration.
The stated optimisation problem is difficult to solve since it is defined by hybrid differential equations ( [18,19]). Hybrid systems are described by both discrete and continuous variables-discrete variables indicate regions in which unique systems equations are applied, continuous variables are solutions to these equations. The model of haemodialysis has three discrete states which are defined by the phosphorus concentration threshold values C PO 4 min and C PO 4 max : if C PO 4 IC has lower value than C PO 4 min then the process is in the first discrete state; if C PO 4 IC is greater, or equal to C PO 4 min but has lower value than C PO 4 max then we say that the process is in the second discrete state; eventually when C PO 4 IC assumes the values greater than C PO 4 max , the system is in the third discrete state. In each discrete state, the haemodialysis process is described by a different set of differential equations. The switch from one discrete state to another is trigged when one of the switching conditions is satisfied: The numerical procedure which we used to solve the problem (9), (1)-(8),(10)-(15) is described, to much extent, in [20] and [21]. The main features of the procedure are: a) it is based on the Radau IIa version of a Runge-Kutta method for integrating differential equations, b) it uses adjoint equations to evaluate gradients of functions defining the optimization problem.
As far as the point a) is concerned we use a Runge-Kutta method in our optimisation procedure since we are dealing with optimal control problems and for these problems these integration procedures are the most suitable (the justification of that claim is given in Chapter 6 of [22]). Our optimisation method is based on the RADAU5 procedure ( [23]) which we had to modify by incorporating into it subroutines for the location of switching points t s .
Furthermore, we had to add to the RADAU5 procedure the subroutine for the consistent evaluation of adjoint equations associated with the equations which are used by the RADAU5 procedure for the evaluation of system equations. The second modification was needed in order to realize the feature b) of our method. In our opinion the optimisation method for solving optimal ontrol problems with hybrid systems should follow the scheme in which at every optimisation procedure step system equations are solved first (due to the necessity of the switching points location) and then values and gradients of optimization functions are calculated. If we pursue the scheme we are in fact limited to the use of adjoint equations in gradient evaluations (more on that issue is in [20]).
The overview of our optimisation algorithm is as follows: whether current controls u k satisfy necessary optimality conditions. If this is the case then STOP, otherwise find new controls u k+1 (using some algorithm which refers to the gradients), increase k by one and go to Step 2.

Simulation and optimisation results
Both simulation and optimization were performed on the data of two patients. We first performed simulation of the model stated in the previous sections to verify its correctness and to analyse the relative behaviour of trajectories Haemodialysis process simulation was conducted in OpenModelica (which implements Modelica standard [24] for simulation) using integration method dassl, with tolerance 1e-6, simulation period was set to 240 time units. Exemplary simulation results of variations of an urea level and a phosphorus level in plasma for the presented equations of IV-compartment system are shown in Fig. 3. They show that while lowering phosphorus level we at the same time also decrease urea quantities. That relative behaviour of urea and phosphorus concentrations justifies our optimisation problem in which we only minimize the urea concentration at the end of the haemodialysis process (in that way we avoid solving more difficult multicriteria optimisation problem).
In the optimisation problem we used the model with the same coefficients as in the simulation experiment with the exception that control variables Q D , Q B , U ufr were not fixed but determined by the optimisation procedure outlined in the previous section.
For the first patient three optimization problems were solved for different constraints imposed on control variables -we changed the optimisation parameters such as We have observed that for the data of both patients, in order to minimise the objective function, the solver set the variable values (Q B and Q D ) to the maximum permissible values. It becomes obvious if we notice that we minimise the concentration of urea, so increasing the flow of fluids (blood and dialysate) is necessary.
However, the above conclusion is no longer valid when we impose additional constraints, for example on the value of C PO 4 IC at final time t D -in that case in addition to constraints (10)-(15) we impose the constraint (that case is stated in Table 2 as Run no. 3) In the third optimization run the upper boundary Q max D is equal to 550 [ml/min] while the obtained optimal solution of Q D is equal to 504 [ml/min]. This indicates that there is a conflict between the removal of urea and achieving the desired level of phosphorus concentration at the end of haemodialysis -using optimisation can help to resolve it. On one hand the optimization procedure tries to reach maximum level of Q D in order to remove urea as much as possible, on the other hand the maximum level Q max D does not guarantee the desired level of phosphorus concentration, eventually the optimization procedure finds the compromise value of Q D . It seems that this type of optimisation problem should be the important one, since the avoidance of low level of phosphorus concentration (and eventually the avoidance of rebound) is the important issue when planning haemodialysis.
As far as optimal values for the control U ufr are concerned, initially we assumed that For this optimization run (the data of the optimization problem are the same as the data for the first patient and run no. 1 in Table 2) values of the other optimal controls were: Q B = 300 [ml/min] and Q D = 550 [ml/min]. These results show that when we consider an optimization problem with three controls Q B , Q D and U ufr , then the optimal profile of the rate of ultrafiltration assumes the biggest possible values at the begining of a haemodialysis process and zero values afterwards. Although ultrafiltration could be modelled, we chose not to do so due to the its decreasing role in routine clinical practice. However, the simulations allows using any profile of ultrafiltration to be modelled in the haemodialysis treatment.
For the prescribed haemodialysis parameters and optimal control values, shown in Table 2, the haemodialysis results for the first patient are shown in Fig. 4. Dialysis for the first patient proceeded smoothly and the minimum concentration of phosphorus C PO 4 min in the body was not exceeded.
During dialysis, second patient experienced a rebound phenomenon due to exceeding the minimum permissible concentration of phosphorus in his body (C PO 4 min )-the optimal trajectories obtained for the second patient are shown in Fig. 5.
From these results, one can conclude that in order to achieve the optimal behaviour of haemodialysis (i.e. minimise urea concentration) one could take constant values for Q B and Q D (dependent on patient and equipment parameters) -see optimal controls in Table 2. Figure 5 presents the optimal trajectories representing parameter changes (C urea EC and C urea IC ) for the second patient. Notice that C urea EC and C urea IC steadily decrease although C PO 4 IC assumes a constant value from some time. Furthermore, Fig. 5 reaffirms our earlier observation that by forcing concentrations of C urea EC and C urea IC to lie below some low values at the end of the haemodialysis process we also guarantee that concentrations of C PO 4 IC and C PO 4 EC steadily decrease.
In the case of a phosphorus concentration for the second patient data one may spot a rebound when a minimum  In this case a phosphorus extraction mechanisms are activated to preserve a phosphorus balance. It can be seen after 120 min of haemodialysis. From that time the value of C PO 4 IC does not change and the system exhibits a 'sliding mode' during which it evolves -the details of the system behaviour in this mode are stated in [21].
A rebound process does not appear during urea removal. There is no body mechanism that preserves a minimum level of a toxic urea value. So, a proper balance of urea concentration is close to zero (theoretical value never gained because of nutrition process).

Discussion
In clinical practice, methods are required to assist nephrologists in the prescription and delivery of safe and effective haemodialysis to patients. There are a number of established physiological models to predict the effect of dialysis on certain solutes, and these allow nephrologists to anticipate the patient response to changes in haemodialysis operating parameters. However, the effects of changing these parameters on urea concentrations, phosphorus concentrations and water volumes are all dependent on each other through differential equations, and it is difficult to predict the patient response for all of these solutes simultaneously. As such, it is also difficult to predict the optimal haemodialysis settings that will result in a different solutes all simultaneously falling within pre-specified limits of safety and efficacy. The use of optimization techniques can cope with problems containing multiple and interdependent differential equations. In our paper, we illustrate the utility of a framing the issue of water and solute changes during haemodialysis as a dynamic optimization problem. We using existing physiological models of solute transport for different solutes, all imbedded within a simulation environment, in a way that results in plausible solutions to this difficult multicriteria optimization problem. Our paper illustrates the use of this technique. We presented the haemodialysis trajectories for two archetypal patients. The first represents a group that have a rather high level of phosphorus concentration at the start of dialysis with an intention to decrease levels to within, but not below, a clinically acceptable range. The second represents those with starting with a rather low concentration of phosphorus. We set Q B and Q D within acceptable ranges, and dialyser clearances to manufacturer specifications, as can be seen in Table 2. For the sake of simplicity, we set K urea r (urea clearance associated with residual renal function) to zero and urea generation during the dialysis to zero. Overall the simulations behaved appropriately and with adequate accuracy. The simulations demonstrated the expected interdependence between ultrafiltration process and the reduction of solute concentrations.
The results of our modelling show that in order to decrease the level of concentrations as much as possible we should use the maximal allowed Q B and Q D , which is not surprising. Our results also showed that the optimal strategy for solute removal is to maximize U ufr at the beginning of haemodialysis, which is also not surprising since this is when the greatest solute mass is present in the blood. However, when we need to achieve a certain reduction in urea concentrations but at the same time keep phosphorus concentrations above some prescribed level, these strategies are no longer appropriate and a different set of values for Q B and Q D are required. For such patients, the use of dynamic optimization allowed us to estimate these parameters in a reasonably precise manner to achieves a satisfactory reduction in urea concentrations, but at the same time resulting in a reduction in serum phosphate concentration that is not excessive. The use of a IV-compartment model in the simulation allowed detailed modelling of the time-concentration profile of phosphorus, and prediction of the presence, extent, and timing of rebound for given (in particular optimal) Q B and Q D . Under these conditions, our simulations showed significant rebound of phosphate between 120 and 180 min of haemodialysis.
Our experiment shows that it is feasible to develop a simple interface that uses in the background sophisticated numerical tools (such as numerical integrators for hybrid systems and optimisation procedures for optimisation problems with hybrid systems), with inputs from users that are no different from what is already undertaken in routine clinical care by medical staff while carrying out haemodialysis -routine patients' parameters, and routine haemodialysis operating parameter. The numerical procedures in OpenModelica can handle this dynamic optimization problem as a simulation with clear and comprehensible results, although it would be better if this functionality in the software was available with sliding modes (unavailable at present). Furthermore, although we used plausible and generic inputs, the model can be individualised to any patient undergoing haemodialysis, with required inputs that can be guessed from cumulative clinical experience or even measured K urea r , C PO 4 EC (0), C urea EC (0), C PO 4 min , C PO 4 max , V (0), G urea (0). Our future research pertaining to the proposed approach will include clinical trials aimed, first of all, at verifying model accuracy, in particular in the case of patients who experience phosphorus rebound during haemodialysis process. Our model contains parameters α, β related to the rebound effect, these parameters should be estimated on a group of patients by performing the standard procedure for nonlinear regression models ( [25]): on a subgroup of patients these parameters are estimated by solving the corresponding nonlinear least squares problem (see [26], [27] for the description of possible numerical procedures to carry out that task) and on the other subgroup of patients the performance of the calibrated model is examined. Note that results in ( [7]) indicate that parameters related to rebound do not vary in models used in short and conventional haemodialysis treatments so these parameters could also be estimated individually for patients after performing a trial haemodialysis and then used in the subsequent haemodialyses. Secondly, clinical trials are needed to resolve doubts about negative consequences of performing haemodialysis on the basis of solutions to optimisation problems. To this end stage patients used in clinical trials should be examined with respect to other profiles which are not directly included in optimisation models, such as pH and HCO 3 variation during dialysis ( [28]), changes in potassium, sodium, calcium levels and in uremic toxins not routinely controlled in every day practice.
The proposed approach to haemodialysis planning could be applied with other kinetics models which contain parameters that can be measured for individual patients, or could be based on estimates valid across groups of patients (so these estimates could be provided after clinical trials on these groups). Therefore, we intend to extend the presented model by including in it a submodel presenting the nature and the rate of vascular refilling during haemodialysis and ultrafiltration. In [28] such a submodel is presented that links together plasma volume, interstitial volume and the protein concentrations in both compartments. Another extension of our model can incorporate differential equations representing potassium kinetics and constraints related to a desired concentration of potassium at the end of haemodialysis (see, for example, [29]). It should be noticed that potassium concentration can exhibit rebound during haemodialysis ( [30]) so the model which includes potassium kinetics is hybrid and our optimization procedure could cope with the problem extended in this way.
All the mentioned directions of planned future research concerning the proposed approach should help to implement a coherent dialysis model suitable for use in real haemodialysis.

Conclusions
Parameters of the haemodialysis process are routinely determined on the basis of a model of the process and its simulation. Having a model of the process we can go further by employing simulations in optimization of a haemodialysis process. The paper considers several optimization problems associated with haemodialysis and shows results of solving one of them. In that way the paper indicates that applying optimization in haemodialysis is possible and that it can lead to the process improvement.