Application of dynamic optimisation for planning a haemodialysis process

Abstract

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 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 multi–compartment models, have been applied to describe phosphorus kinetics ([7–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.

Methods

Models of urea kinetics

The II–compartment model depicts a patient as a two–compartment 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.

Fig. 1

II-compartment kinetic model of toxins and substances

According to [11, 14] the equations describing II–compartment model are as follows:

$$ {\begin{aligned} \frac{{dC}_{EC}^{urea}}{dt} & = \frac{K_{IE}^{urea} \cdot \left(C_{IC}^{urea}-C_{EC}^{urea}\right) - C_{EC}^{urea} \cdot \left(K_{D}^{urea}+K_{r}^{urea}+K^{ufr}\right) }{0.34 \cdot V(0)-UFR} \end{aligned}} $$

(1)

$$ {\begin{aligned} \frac{{dC}_{IC}^{urea}}{dt} & = \frac{K_{IE}^{urea}\cdot\left(C_{EC}^{urea}-C_{IC}^{urea}\right)+G^{urea}}{0.66 \cdot V(0)} \end{aligned}} $$

(2)

$$ {\begin{aligned} \frac{dUFR}{dt} & = U_{ufr} \end{aligned}} $$

(3)

where:

UFR - ultrafiltration volume

U_ufr - profile of the ultrafiltration rate

K^ufr - ultrafiltration coefficient

$K_{r}^{urea}$ - urea clearance associated with residual renal function

$K_{D}^{urea}$ - dialyser clearance

$\frac {{dC}_{EC}^{urea}}{dt}$ - variation of urea concentration in the extracellular fluid

$\frac {{dC}_{IC}^{urea}}{dt}$ - variation of urea concentration in the intracellular fluid

$K_{IE}^{urea}$ - 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.

Fig. 2

IV-compartment kinetic model of toxins and substances

As far as phosphorus kinetics is concerned we have employed the model presented in ([7, 11]):

$$\begin{array}{*{20}l} \frac{{dC}_{EC}^{{PO}_{4}}}{dt} & = & \frac{K_{IE}^{{PO}_{4}} \cdot (C_{IC}^{{PO}_{4}}-C_{EC}^{{PO}_{4}}) - K_{D}^{{PO}_{4}} \cdot C_{EC}^{{PO}_{4}}}{0.34 \cdot V(0)-UFR} + \end{array} $$

(4)

$$\begin{array}{*{20}l} && K_{3}^{{PO}_{4}} + K_{4}^{{PO}_{4}} \\ \frac{{dC}_{IC}^{{PO}_{4}}}{dt} & = & \frac{K_{IE}^{{PO}_{4}} \cdot (C_{EC}^{{PO}_{4}}-C_{IC}^{{PO}_{4}}) }{0.66 \cdot V(0)} \end{array} $$

(5)

where:

$\frac {{dC}_{EC}^{{PO}_{4}}}{dt}$ - variation of phosphorus concentration in the extracellular fluid

$\frac {{dC}_{IC}^{{PO}_{4}}}{dt}$ - variation of phosphorus concentration in the intracellular fluid

$K_{IE}^{{PO}_{4}}$ - diffusion speed between inner and outer cellular solute through a cellular membrane - the same in both directions

$K_{D}^{{PO}_{4}}$ - phosphorus clearance coefficient in dialyzer

$K_{3}^{{PO}_{4}}$ - phosphorus generation coefficient - used after exceeding the upper boundary of phosphorus concentration for a healthy person

$K_{4}^{{PO}_{4}}$ - phosphorus generation coefficient - used after exceeding the lower boundary of phosphorus concentration for a healthy person

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):

$$\begin{array}{*{20}l} K_{4}^{{PO}_{4}} &=& \alpha \cdot \max \left (C_{min}^{{PO}_{4}}-C_{IC}^{{PO}_{4}},0\right) \end{array} $$

(6)

$$\begin{array}{*{20}l} K_{3}^{{PO}_{4}} &=& \beta \cdot \max \left (C_{max}^{{PO}_{4}}-C_{IC}^{{PO}_{4}},0\right) \end{array} $$

(7)

where:

α,β - coefficients that describe the speed of returning the phosphorus concentrations to the initial levels. $K_{3}^{{PO}_{4}}$ is activated only when a phosphorus level is below $C_{max}^{{PO}_{4}}$ value (third compartment). And $K_{4}^{{PO}_{4}}$ is activated only when a phosphorus level is below $C_{min}^{{PO}_{4}}$ respectively (fourth compartment).

Furthermore, we underline that a phosphorus clearance coefficient $K_{D}^{{PO}_{4}}$ is a function of two parameters: predefined by producer phosphorus clearance of the dialyser membrane $K_{D,base}^{{PO}_{4}}$, 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).

$$ K_{D}^{{PO}_{4}} = \kappa \cdot K_{D,base}^{{PO}_{4}}. $$

(8)

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.

Table 1 Parameters for simulation of haemodialysis process

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:

$$ \min_{Q_{B},Q_{D},U_{ufr}}C_{EC}^{urea}(t_{D}) $$

(9)

subject to the constraints (1)–(8), the following constraints at final time t_D

$$\begin{array}{*{20}l} & C_{EC}^{urea}(t_{D}) &\leq L^{urea}_{EC} \end{array} $$

(10)

$$\begin{array}{*{20}l} & C_{IC}^{urea}(t_{D}) &\leq L^{urea}_{IC} \end{array} $$

(11)

$$\begin{array}{*{20}l} L_{min}^{UFR} \leq & UFR(t_{D}) &\leq L_{max}^{UFR}, \end{array} $$

(12)

and the constraints on the control variables

$$\begin{array}{*{20}l} Q_{B}^{min} \leq & Q_{B}(t) &\leq Q_{B}^{max} \end{array} $$

(13)

$$\begin{array}{*{20}l} Q_{D}^{min} \leq & Q_{D}(t) &\leq Q_{D}^{max} \end{array} $$

(14)

$$\begin{array}{*{20}l} U_{ufr}^{min} \leq & U_{ufr}(t) &\leq U_{ufr}^{max},\ t\in [0,t_{D}], \end{array} $$

(15)

where $L^{urea}_{EC}$ and $L^{urea}_{IC}$ mean maximum admissible values for urea concentrations. $L_{min}^{UFR}$ and $L_{max}^{UFR}$ mean minimum and maximum admissible values for the ultrafiltration. $Q_{B}^{min}$, $Q_{B}^{max}$, $Q_{D}^{min}$, $Q_{D}^{max}$, $U_{ufr}^{min}$, $U_{ufr}^{max}$ define bound constraints on decision variables, in particular $U_{ufr}^{min}$, $U_{ufr}^{max}$ 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_{min}^{{PO}_{4}}$ and $C_{max}^{{PO}_{4}}$: if $C_{IC}^{{PO}_{4}}$ has lower value than $C_{min}^{{PO}_{4}}$ then the process is in the first discrete state; if $C_{IC}^{{PO}_{4}}$ is greater, or equal to $C_{min}^{{PO}_{4}}$ but has lower value than $C_{max}^{{PO}_{4}}$ then we say that the process is in the second discrete state; eventually when $C_{IC}^{{PO}_{4}}$ assumes the values greater than $C_{max}^{{PO}_{4}}$, 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: $C_{IC}^{{PO}_{4}}(t_{s})=C_{min}^{{PO}_{4}}$, or $C_{IC}^{{PO}_{4}}(t_{s})=C_{max}^{{PO}_{4}}$, for some switching time t_s.

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:

• 1. Set initial values of the controls u₀ (u=(Q_B,Q_D,U_ufr)) and k=0.

• 2. Integrate hybrid system equations for given initial conditions and the current controls u_k, while integrating system equations locate the switching points $\{{t_{s}^{l}}\}$. Calculate values of all functionals which define the optimisation problem.

• 3. Evaluate adjoint equations for each optimisation problem’s functional. Having trajectories of adjoint equations evaluate gradients of the problem functionals.

• 4. On the basis of calculated gradients determine 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.

Results

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 $C^{urea}_{EC}$, $C^{urea}_{IC}$, $C^{{PO}_{4}}_{EC}$ and $C^{{PO}_{4}}_{IC}$.

The coefficients used in simulation are given in Table 1. Notice that we assumed the following constant values for the control variables: Q_D=500[ml/min] and Q_B=280[ml/min]. A constant ultrafiltration rate (U_ufr=4.08[ml/min] for a first patient, U_ufr=11.25[ml/min] for a second patient) was assumed throughout haemodialysis. These values, together with the other haemodialysis parameters, may be considered as typical.

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).

Fig. 3

Simulation results for IV-compartment model (first patient). a IV-compartment model presenting variations of urea concentration during haemodialysis. b IV-compartment model presenting variations of phosphorus concentration during haemodialysis

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 $Q_{D}^{min}$, $Q_{D}^{max}$, $Q_{B}^{min}$ and $Q_{B}^{max}$. In the first run we set $Q_{D}^{max}=550$ [ml/min] and $Q_{B}^{max}=300$ [ml/min], in the second run $Q_{D}^{max}=650$ [ml/min] and $Q_{B}^{max}=400$ [ml/min]. In the second run we increased the control parameters $Q_{D}^{max}$ and $Q_{B}^{max}$ to achieve assumed final concentrations of urea.

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_{IC}^{{PO}_{4}}$ 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)

$$\begin{array}{@{}rcl@{}} L_{IC}^{PO} & \leq & C_{IC}^{{PO}_{4}}(t_{D}). \end{array} $$

(16)

Table 2 Parameters for optimisation of haemodialysis process

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_{D}^{max}$ 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 $U_{ufr}^{min}= 0$ [ml/min] and $U_{ufr}^{max}= 20$ [ml/min] and the optimization procedure found the following optimal trajectory for U_ufr: U_ufr=20.00 [ml/min] on the subinterval [0,48), U_ufr=4.91 [ml/min] on the subinterval [48,62) and U_ufr=0 [ml/min] on the subinterval (62,240] (the time values 48 and 62 follows from the fact that we used piecewise constant approximations for controls and in the approximation the number of subintervals was equal to 10). 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_{min}^{{PO}_{4}}$ in the body was not exceeded.

Fig. 4

Optimisation results for the first patient. a State trajectories for the urea $C_{EC}^{urea}$ and $C_{IC}^{urea}$ for the first run. b State trajectories for the phosphorus $C_{IC}^{{PO}_{4}}$ and $C_{EC}^{{PO}_{4}}$ for the first run (c) State trajectories for the urea $C_{EC}^{urea}$ and $C_{IC}^{urea}$ for the second run. d State trajectories for the phosphorus $C_{IC}^{{PO}_{4}}$ and $C_{EC}^{{PO}_{4}}$ for the second run

During dialysis, second patient experienced a rebound phenomenon due to exceeding the minimum permissible concentration of phosphorus in his body ($C_{min}^{{PO}_{4}}$)— the optimal trajectories obtained for the second patient are shown in Fig. 5.

Fig. 5

Optimisation results for the second patient. a State trajectories for the urea $C_{EC}^{urea}$ and $C_{IC}^{urea}$. b State trajectories for the phosphorus $C_{IC}^{{PO}_{4}}$ and $C_{EC}^{{PO}_{4}}$ (rebound of $C_{EC}^{{PO}_{4}}$ can be seen in 120 min of haemodialysis)

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_{EC}^{urea}$ and $C_{IC}^{urea}$) for the second patient. Notice that $C_{EC}^{urea}$ and $C_{IC}^{urea}$ steadily decrease although $C_{IC}^{{PO}_{4}}$ assumes a constant value from some time. Furthermore, Fig. 5 reaffirms our earlier observation that by forcing concentrations of $C_{EC}^{urea}$ and $C_{IC}^{urea}$ to lie below some low values at the end of the haemodialysis process we also guarantee that concentrations of $C_{IC}^{{PO}_{4}}$ and $C_{EC}^{{PO}_{4}}$ steadily decrease.

In the case of a phosphorus concentration for the second patient data one may spot a rebound when a minimum phosphorus level is reached (at the level 0.70 [mmol/L]). 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_{IC}^{{PO}_{4}}$ 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_{r}^{urea}$ (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_{r}^{urea}$, $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³ 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.