Optimal EPO dosing in hemodialysis patients using a non-linear model predictive control approach.

Anemia management with erythropoiesis stimulating agents is a challenging task in hemodialysis patients since their response to treatment varies highly. In general, it is difficult to achieve and maintain the predefined hemoglobin (Hgb) target levels in clinical practice. The aim of this study is to develop a fully personalizable controller scheme to stabilize Hgb levels within a narrow target window while keeping drug doses low to mitigate side effects. First in-silico results of this framework are presented in this paper. Based on a model of erythropoiesis we formulate a non-linear model predictive control (NMPC) algorithm for the individualized optimization of epoetin alfa (EPO) doses. Previous to this work, model parameters were estimated for individual patients using clinical data. The optimal control problem is formulated for a continuous drug administration. This is currently a hypothetical form of drug administration for EPO as it would require a programmable EPO pump similar to insulin pumps used to treat patients with diabetes mellitus. In each step of the NMPC method the open-loop problem is solved with a projected quasi-Newton method. The controller is successfully tested in-silico on several patient parameter sets. An appropriate control is feasible in the tested patients under the assumption that the controlled quantity is measured regularly and that continuous EPO administration is adjusted on a daily, weekly or monthly basis. Further, the controller satisfactorily handles the following challenging problems in simulations: bleedings, missed administrations and dosing errors.

Introduction

According to the 2018 United States Renal Data System annual data report (USRDS 2018), approximately 2.5 million patients were treated world-wide for end-stage renal disease in 2016. In nearly all reporting countries hemodialysis (HD) was the predominant form of dialysis therapy. On December 31, 2016, there were 726,331 prevalent cases of end-stage kidney disease in the U.S., of which 63.1% were receiving HD. Due to reduced erythropoietin production in the kidneys almost all HD patients suffer from a chronically decreased number of circulating red blood cells (RBCs) and associated low hemoglobin (Hgb) levels. This condition is called anemia. Untreated anemia is associated with poor quality of life and increased morbidity and mortality. Therefore, physicians aim for a partial correction of anemia with erythropoiesis stimulating agents (ESA). In this paper we consider the treatment with epoetin alfa (EPO), a human recombinant erythropoietin produced in cell culture.

Currently, the recommended Hgb target range for anemia management is 10–12 g/dl (see Mactier et al. 2011). Albeit, the KDIGO (Kidney Disease Improving Global Outcomes) Work Group, which provides clinical practice guidelines for anemia in chronic kidney disease and guidance on diagnosis, evaluation, management and treatment of HD patients, recommends not to exceed the limit of 11.5 g/dl in general but suggests to individualize therapy for patients whose quality of life may be improved at Hgb levels above 11.5 g/dl. Usually, dialysis facilities use dosing protocols that work as follows: A starting dose gets specified and based on the resulting Hgb change and depending on whether the patient is within the Hgb target range the ESA dose gets adjusted. This one-size-fits-all approach results in about 65% of patients achieving the set Hgb target. The Dialysis Outcomes and Practice Patterns Study (DOPPS) Practice Monitor reports that in the US, since December 2014, the percentage of patients with Hgb above 12 g/dl is 14–15% while 18–20% of patients have a Hgb below 10 g/dl. Moreover, Hgb variability and cycling are well known to occur in HD patients treated with ESA (Berns et al. 2003; Fishbane and Berns 2005, 2007). According to Yang et al. (2007) a greater Hgb variability is independently associated with higher mortality. The challenge in anemia treatment is the patients’ difference in long-term Hgb response to ESA. The drug concentration in plasma influences the maturation, proliferation and apoptosis of cells in the erythroid lineage. However, these cells remain about two weeks unobservable in the bone marrow before being released into the bloodstream. Consequently, it is difficult to anticipate the resulting delayed effect of the drug administration.

The mathematical model of erythropoiesis presented in Fuertinger et al. (2013) predicts patients’ erythropoietic response. We utilize this model to design a model-based feedback controller. The successful use of control algorithms for drug dosing has been shown, for example, in Magni et al. (2007) and in Bequette (2013). Both works are about closed-loop insulin dosing (i.e., the artificial pancreas). The control strategy we use is called model predictive control (MPC), also known as moving or receding horizon control. In 2011, Brier and Gaweda have already published an MPC-based algorithm for improved anemia management that has been tested and validated in clinical studies. Unlike our approach their predictive model is based on the concept of artificial neural networks (see also Barbieri et al. 2016; Brier et al. 2010). One of the limitations of using an artificial neural network approach is the need for large training and validation data sets making it very difficult to fully tailor this approach to individual patients and to personalize target Hgb ranges for patients whose quality of life would improve from a Hgb level higher than 11.5 g/dl as suggested by the KDIGO work group.

The purpose of this work is to develop a controller scheme that is fully personalizable. Previous to the study, the used model was adapted to individual patients using clinical data on Hgb levels and EPO administration. The details of the parameter estimation procedure and its results were published in Fuertinger et al. (2018). Note that insufficient iron availability is not modelled explicitly. However, several bone marrow parameters of the erythropoiesis model that are influenced by iron availability were estimated on a per patient level. The presented feedback controller is tested on various patient parameter sets. Throughout this manuscript we assume that EPO is administered continuously. This is currently not done in clinical practice but could be achieved with an “EPO pump” similar to the programmable insulin pumps used to treat diabetes mellitus. The approach is chosen nonetheless since it yields a continuous control which provides the best situation for stabilizing a system. Having developed a functional controller scheme it can then be adapted to actual dosing schedules and the effect of reducing the dosing frequency on the achievable stability can be analyzed. The given model equations are coupled hyperbolic partial differential equations (PDEs) and the control variable enters these equations non-linearly. In the presence of a non-linear model MPC is referred to as non-linear MPC (NMPC). More details on MPC can be found, e.g., in the books by Grüne and Pannek (2011) and by Rawlings and Mayne (2009). The basic principle of MPC consists of repeatedly solving finite horizon open loop optimal control problems. In each step, an open loop problem is solved. Then, only the first component of the obtained optimal control is applied and the optimization horizon gets pushed. This allows to include measurements and to react to unforeseen disturbances or complications. Here, we assume to know the true model and that we are able to measure Hgb perfectly. We simulate (gastro-intestinal) bleedings which are a common complication in HD patients and consider a malfunction of the pump by simulating the complete failure to administer EPO for an entire day or by applying an incorrect dose. Further, we simulate different frequencies of rate change with which the pump is programmed, ranging from daily to only once a month.

The paper is organized as follows: In Sect. 2 we introduce the control variable and present the model equations. The numerical approximation of these so-called state equations is investigated in Sect. 3. Both the model and its numerical approximation are recalled from Fuertinger et al. (2013). In addition, we regularize the erythrocytes model equation to obtain differentiability required for numerically solving the optimal control problem utilizing first-order optimality conditions. In Sect. 4 we formulate the optimal control problem and the NMPC algorithm is described. In Sect. 5 we present our numerical results of the following in-silico experiments: bleedings, missed administrations or wrongly administered doses and the restriction of EPO administration rates to be constant over several weeks. We draw some conclusions in Sect. 6. Finally, all parameters used for simulations are presented in Appendix A.

The model of erythropoiesis

We start by describing the structure of the control process; compare Fig. 1. The control which is a vector of EPO administration rates can be altered to change a patient’s Hgb being the outcome. First, the EPO rates change the patient’s EPO concentration in plasma which affects the production of RBCs. The process of erythropoiesis is described by a mathematical model presented in Fuertinger et al. (2013). The different cell types during erythropoiesis are grouped into five population classes and the model allows to calculate their densities . From a control perspective these are our states. Given the density of erythrocytes we can calculate their number and finally the Hgb.

Fig. 1

Schematics of the control structure

In the following we explain the named components step by step. Note that the lemmas of this section will allow to utilize first-order optimality conditions for solving the optimal control problem numerically.

The dosing of EPO as the control variable

We first write the time-varying EPO concentration in plasma as a function of EPO administration rates. Let us consider the time interval [0, T] with large final time . We assume that EPO can be applied continuously, with a constant administration rate per day or per multiple days. The EPO rates are given in U/day, where U stands for units. The EPO concentration E(t), , in plasma is separated into a constant summand modeling the patient’s remaining endogenous erythropoietin level and a time-dependent summand resulting from the administered EPO:Note that the endogenous EPO production is estimated for each patient individually when adapting the model to patient data. For the number let the days with a constant EPO rate in , , be given asWe introduce the associated finite-dimensional control space . Only non-negative EPO rates with a given upper positive limit can be applied. Therefore, we are considering the (convex and compact) admissible setwhere ‘’ is interpreted componentwise. Throughout, vectors are denoted by boldface letters. Let be given. Then, the time- and control-dependent summand satisfies the following initial value problem (cf. Fuertinger et al. 2013, equation (12)):whereare characteristic functions of the intervals , stands for the total blood volume, is a non-negative initial condition for the exogenous EPO level and is the EPO degradation rate with half-life time . Thus, for , , the EPO concentration is given byIt follows from (4) that is a function mapping from the admissible set into C([0, T]), where C([0, T]) is the space of all continuous functions from [0, T] to . From (4) and we infer thatThus, we define the interval of possible EPO concentrationsand observe that holds for all and .

Lemma 1

The mapping E introduced in (4), has the following properties.

For every the function is continuously differentiable.

The mapping is twice continuously differentiable for any . Its gradient is given as at , , and . Further, the hessian matrix is zero.

Proof

The claims follow directly from formula (4). Since is independent of , the hessian is zero.

The PDE model of erythropoiesis

A schematic of the model is shown in Fig. 2. For the underlying assumptions and more details we refer to Fuertinger (2012), Fuertinger et al. (2013). Stem cells commit to the erythroid lineage at a constant rate . Once a stem cell has committed it passes the five shown cell classes/stages over time (if it does not die along the way): BFU-E (burst-forming unit erythroids), CFU-E (colony-forming unit erythroids), erythroblasts, marrow reticulocytes and erythrocytes (including blood reticulocytes). This means that there is a flux of cells from each population class to the subsequent one. For example, when a CFU-E cell has reached maximum age of that class it leaves the class and becomes an erythroblast with minimum maturity. The population densities depend on maturity and time. For each class an age-structured population model is given which describes the development of the respective class subject to a given EPO concentration in plasma. EPO has a direct effect on the rate of apoptosis of CFU-E cells (), the maturation velocity of marrow reticulocytes () and the mortality rate of erythrocytes (). For each class we are given an individual maturity interval , , in days. The interval boundaries are given by , , , , , , whereas the RBC lifespan is patient-dependent. Ma et al. (2017) have measured this shortened RBC lifespan in HD patients to be in the range of 37.7 to 115.8 days. Note that the first four cell classes are in bone marrow while the fifth class describes cells circulating in blood. That is why is set to 0.

Fig. 2

Schematic of the PDE model of erythropoiesis

For conciseness , we write all five state equations in one form. Suppose that for given and the EPO concentration is given as (4). Then, the equation reads with the maturity interval , the cylinder and the initial condition . The solution to (S.i) denotes the cell density of the respective cell population with maturity x at time t. The function describes the maturation velocity and is of form , where describes the profileration rate and the rate of apoptosis. Actually, the function and the sigmoid functions and depend on the bounded (patient-dependent) parameter vector with . To simplify the notation we do not indicate dependencies on . We refer to Appendix A, where all fixed and all individualized parameters, which we utilize in our numerical experiments, are listed. The individualized parameters are obtained via parameter estimation; see Fuertinger et al. (2018). The functions for the different classes read as follows:where, for , the real-valued functions and are given byand the functionstands for the erythrocytes mortality rate with an EPO threshold for neocytolysis. The (closed) non-empty interval denotes the age interval, where neocytolysis is possible.

Lemma 2

The mappings are continuously differentiable.

The claim follows directly from (6).

We denote the coupled system (S.1)–(S.5) by (S).

Total RBC population

If the erythrocytes population density is known, the total RBC population, , is given aswhere solve the state system (S).

Hgb concentration

Given a patient’s total blood volume (in ml) together with the total number of RBCs P, the Hgb concentration (in g/dl) is calculated aswhere pg denotes the mean corpuscular hemoglobin.

Regularization of the equation for the erythrocytes

In Sect. 4 we formally introduce the non-linear optimal control problem. In order to solve it numerically by utilizing first-order necessary optimality conditions we have to differentiate the state system (S) with respect to the state variable and the control variable . From Lemma 1 we already know that is continuously differentiable for every . Moreover, the mappings and are continuously differentiable by Lemma 2. However, the mapping is non-differentiable for every . Therefore, we have to regularize in order to get smooth state equations.

For that reason we introduce the Heaviside function defined asThen, the mortality rate can equivalently be written aswithFor we utilize the following regularized Heaviside function which is twice continuously differentiable. Notice that the functionis an approximation of . Thus, the function R defined in (10) can be regularized aswhich allows us to replace the non-smooth coefficient function by the smooth (with respect to E) mapping

Lemma 3

For every the mapping is continuously differentiable.

The claim follows directly from (11) because of .

In the sequel we replace in (5) by and hence by to account for the regularized fifth state equation which we denote by (). Let () be the state system (S.1)–(S.4) and ().

Numerical approximation of the state equations

The numerical solution of the age-structured population models is based on semigroup theory. We formulate the five state equations as abstract Cauchy problems which are then approximated by semigroups acting on finite dimensional subspaces. We have compared this discretization to an upwind finite difference scheme (cf. Strikwerda 2004). Due to computational speed while obtaining a similar approximation quality we have favored the semigroup based approach. We refer the reader to, e.g., Ito and Kappel (2002) and Kappel and Zhang (1993) for results on evolution operators and their approximation.

In the following we derive the discretization of the state equations by means of the general form (S.i) where we omit the dependency on i but replace by .

The state equations as abstract Cauchy problems

For we introduce the functionwhere denotes the maturity interval of the respective cell class. Hence, the sequence approximates the -distribution. Let , , be a mild solution of the non-homogeneous Cauchy problemfor given , where the linear operator , , is defined asIt can be shown that in holds for , where y is the solution to the corresponding state equation (S.1)–(S.4) or ().

Since the range for the attribute is different for the cell populations considered, it is useful to normalize these attributes such that the range of the normalized attribute is [0, 1]. In order to achieve this we set and define the mappingwith its inverse given byWe denote by the Hilbert space endowed with the weighted inner productThe induced norm is . The Hilbert spaces and are isomorphic with the isomorphism and its inverse given byApplying the operator to the sequence yieldsThen, from (12), we derive the normalized Cauchy problem in :where the operator , , is defined asand holds. Note that solves (12), where is the mild solution of the Cauchy problem (13).

Approximation of the abstract Cauchy problems

In the following we derive a discretization of the Cauchy problem (13) based on shifted Legendre polynomials. This approach was originally presented in Kappel and Zhang (1993) where they apply it to a very similar type of equation. In the considered 1D case only few basis elements are needed which yields a fast approximation. However, there are situations where one would expect diffculties due to the (oscillating) characteristics of Legendre polynomials. Higher spatial dimension of the hyperbolic equation or large maturation velocities are examples for such situations. In the following, we first recall certain characteristics and further utilized features of Legendre polynomials. For more details we refer the reader to the book by Abramowitz and Stegun (1970).

Legendre polynomials Legendre polynomials are orthogonal polynomials on withThe polynomials can be computed with Bonnet’s recursion formulawhere and Thus, the first five Legendre polynomials are of the following formTo compute the approximations for the partial differential equations we use the following formula representing the derivatives of Legendre polynomials in terms of Legendre polynomials:Furthermore, the following consequences from Bonnet’s formula (eq. (14)) are needed for the approximation:where we set

Another property of the Legendre polynomials which is used is:Approximation Let us define the basis functionsThe sequence is an orthogonal sequence in withFor we introduce the N-dimensional subspace byLet be the orthogonal projection defined asThe approximating delta distribution on is defined byFor any we define the approximating linear operator byNote that the projection needs only be applied on in case of the erythrocytes model equation due to for . Now, the discretization of (13) is given asLetwhere is a solution to (19). It follows that in for . A proof for the case can be found in Kappel and Zhang (1993, Theorem 4.3). Using (17), (18) andwe derive the Galerkin schemefor . Settingwe can express (20) as the following system of N ordinary differential equations

Theorem 1

Let given by (4) for an arbitrarily chosen . Moreover, is assumed to be piecewise continuous for every . Then, the functionis the only one that satisfies (21) at every time instance, where is continuous for every . Furthermore, holds.

We first show that (22) satisfies (21). Using it follows from Lemmas 1-3 that the mappings , and , , are continuously differentiable on [0, T] for every . Hence, the mapping is continuously differentiable for every . We obtain for any , where is continuous Further, if we choose in (22) we get . Hence, satisfies (21) at every time instance, where is continuous.

Uniqueness: Assume there exist two solutions , to (21). We set . Then, it follows that for all , where is continuous. Since and are continuous, we can extend the derivative of by for all . Furthermore, we have . From Gronwall’s inequality we infer that in [0, T], which implies in [0, T].

Remark 1

Due to Theorem 1 the solution to (19) belongs to provided is piecewise continuous for every .

To solve (21) a time integration method has to be applied. In our numerical experiments we apply the implicit Euler method.

The control-to-state operator

Let be a given control variable and be given by (4). Then, the five state variables solving () are approximated by the functionswhere the coefficient vector satisfies This is Eq. (21) with added index i. Due to Theorem 1 we define the non-linear solution operatorwhere is the generalized Cartesian product. Hence, satisfies ()–(). Then, the discretized total RBC population is given by (cf. (7))where again satisfies ()–(). Note that the equations (), can be solved consecutively. At any time the discrete state vector is of length 5N.

The optimal EPO dosing

In this section we formulate the optimal EPO dosing as an optimal control problem on the long time horizon [0, T]. For the number let the days with a constant EPO rate in , , be given as

The optimal control problem

The goal is to stabilize the Hgb around a desired value of 10.5 g/dl in order to bring and keep the Hgb into the target window of 10–12 g/dl. Values in this range are considered as safe. We choose this lower value within the range because in case of an overshoot the Hgb can not be pulled down actively but one has to wait till it decreases of itself. The optimal control problem is formulated for the number of RBCs and not for the Hgb. This means that, for each patient, prior to optimization a desired total amount of RBCs is calculated using formula (8).

In the cost functional, we measure and penalize the deviation of the total RBC population from the desired population over the whole time horizon [0, T]where the discretized total RBC population has been introduced in (24) and satisfies ()–(). Further, are regularization parameters and is a non-negative weight. Now, the optimal control problem is formulated as follows:

Remark 2

Note that () is a non-linear optimization problem. Hence, () is non-convex, so that several local minima might exist.

It can be shown that () possesses at least one optimal control in . For the sake of brevity, we do not present the proof here.

The NMPC method

Problem () can not be treated as an open-loop problem since unforeseen events and disturbances can occur. In reality, predicted and measured Hgb values will differ which has to be taken into account by means of a closed-loop controller. Moreover, patient parameters will not be constant over time which even makes readaptations of the model necessary in actual application.

Let all time intervals of constant administration rates be of same length andSuppose we are at time and consider the time horizon with , and .

Next we introduce a notation for the dayswhich belong to the current horizon ; cf. (2). At we are given the initial conditions and . Then, we define the cost functional on the time horizon with and :where the coefficient vector satisfies ()–() on the time horizon with initial conditions , . In (26) we have added a summand penalizing the deviation of the population at the final time for stability reasons. For we define the admissible setConsidering a constant EPO rate per day in our numerical experiments we set for some constant . Note that is the length of the prediction horizon. Now, the controller predicts the future evolution of the system under control over this prediction horizon, the cost functional gets minimized and we obtain an optimal control vector for the given time period. Then, (only) the first component of the optimal control vector is applied and yields new initial conditions for the next initial time point , to where the finite horizon gets pushed. In Algorithm 1 we summarize the NMPC method. The iterative computation of the control can be seen in line 5 where in each iteration the first component of the optimal open loop solution to problem () from line 4 is appended to the already given control vector. Then, this control is applied on the corresponding time interval and the initial condition for the next iteration is set to ; see line 6. If the state vectors can be measured or estimated based on measurements, the new initial condition in line 6 can be set therewith.

Numerical solution of the open-loop problem ()

It remains to discuss how the open-loop problem in line 4 of Algorithm 1 is numerically solved. In this work we apply a projected BFGS method with Armijo linesearch as described in Kelley (1999, Sect. 5.5.3). This requires to compute the gradient at . This is done using a Lagrangian-based approach as shown in Hinze et al. (2009)[Sect. 1.6.4]. The gradient is given bywith the adjoint variable , where satisfies

for , with , , andNote that the adjoint state is the Lagrange multiplier associated with the state equation.

We have presented the gradient and the adjoint equations ()–() in a very compact format. So let us explain what it actually means to calculate the gradient at :

Solve the state equations ()–() on the time horizon forward in time.

Solve the adjoint state equations ()–() on the time horizon backward in time.

Compute the gradient .

Numerical results

For our numerical experiments we have chosen the data sets from five patients which capture the main occurring characteristics. The constant endogenous erythropoietin concentration for example is once far above the threshold for neocytolysis (), twice just slightly below and twice clearly smaller. The aim was to find a general optimization setting that works for diverse data sets. The parameters can be looked up in Table 8.

Table 8

Individualized parameters for considered patients

	Patient no.
Parameter	1	2	3	4	5
c_tbv	3567	7373	4543	6102	2600
E^end	42.60	79.42	282.46	43	79
μ_2	0.016	0.009	0.008	0.011	0.010
μ_6	0.122	0.049	0.072	0.049	0.032
T_1/2	4.78	9.65	5.59	10.80	5.12
x¯_5	92.3	74.3	43.1	94.8	50.5

For computation of the population densities we scale the hyperbolic equations by which is legitimate since the equations are linear with respect to the state variables , . For discretization we use Legendre polynomials and the time step size .

We first have a look at the predicted Hgb curves when no EPO is administered. Then, we show how the NMPC algorithm is able to correct anemia under the assumption of knowing the true model and patient parameters and being able to continuously and perfectly measure Hgb.

Uncontrolled Hgb concentration

We begin by taking a look at the predicted Hgb concentrations without EPO administration. As can be seen exemplarily in Fig. 3 the Hgb levels are running in a steady state far below the target range. This state of anemia would be critical for patients since it increases cardiovascular disease and death risk (Strippoli et al. 2004).

Fig. 3

Uncontrolled Hgb levels. The dotted lines mark the target range of 10–12 g/dl

NMPC

Settings

Throughout we consider a total time period of 24 weeks, i.e. days. Excluding our last experiment, the days are given by and the length of the NMPC horizon is four weeks. This means, we set , and . Note that the length of the horizon determines the computational effort of the NMPC method: the longer the horizon, the more time it takes to calculate the numerical solution to the open-loop problem in each step. Hence, one is searching for the minimal horizon that provides stability of the NMPC closed-loop. In our simulations we have observed that four weeks serves the purpose for the named setting. The fact that the cells remain about two weeks in the bone marrow (first four classes) before released into the bloodstream explains the need for this length. In pratice, the maximum cumulative weekly EPO dose is 60,000 U with up to three administrations a week. Given our higher administration frequency we set the maximum weekly EPO dose to 7000 U which leads us to set U/day.

Given the time horizon we setwhere only the parameter is chosen in an individualized manner.

We start by testing the NMPC algorithm for different constants . This parameter penalizes control costs in the objective functional. A larger results in a stronger penalization of control costs, i.e. the controller must try to bring and keep the Hgb curve in the target range with less EPO. In Fig. 4 we present the results for patients 1 and 2. For both patients the value is the right choice for control costs penalization in the sense that it allows to accurately stabilize the Hgb around 10.5 g/dl. The corresponding EPO rates can be seen in Figs. 5 and 6. Interestingly, the smallest penalization in patient 1 results in a very specific and time sensitive administration pattern. Patient 2 requires a lot higher doses than patient 1. It is actually about 2.5 times as much; see Table 1. When looking at the uncontrolled Hgb levels in Fig. 3 patient 2 shows a higher asymptotic Hgb level than patient 1. So the higher doses are actually counterintuitive which underlines the need for a closer look into the time dynamics. For the Hgb concentration for patient 1 is still within the target range while this penalization results in a drop below the target range for patient 2. Still, for the total population for patient 2 remains as well within the target range. Using for patient 1 leads to a Hgb level below the target range. This indicates that the penalization of the control costs needs to be chosen carefully on a subgroup or even per-patient level.

Fig. 4

Optimal Hgb curves for different values for patients 1 and 2. The dotted lines mark the target range

Fig. 5

Patient 1: optimal EPO rates for different values

Fig. 6

Patient 2: optimal EPO rates for different values

Table 1

Total EPO doses for different values for patients 1 and 2

	Patient 1			Patient 2
c_γ	0.1	10	100	0.1	1	10
Total dose (U)	48,323	42,683	30,616	121,354	119,114	71,706

Our primary goal is Hgb in target. Thus, we choose such that the penalization of control costs does not block the NMPC algorithm from controlling the Hgb levels into the target range. The above test we have also done for the other patients. For patient 5 the value is fine as well while patients 3 and 4 require a weaker penalization of . Concluding, with the Hgb levels of all patients can be brought into target only that for patients 1,2 and 5 this can as well be achieved with and consequently a lower total EPO dose.

Bleeding

In the following we analyze how the NMPC algorithm handles a sudden and unforeseen (gastro-intestinal) bleeding. This is a frequent complication in HD patients. Note that the NMPC algorithm learns about the bleeding first after it has happened by an assumed measurement. In line 6 of Algorithm 1 the initial condition is overwritten based on that measurement. The results for patients 3 and 4 for a bleeding which corresponds to a drop in Hgb down to 7.5 g/dl are shown in Fig. 7. Patient 4 has a fast recovery time so that we add a later smaller bleeding to 9.0 g/dl. According to Pottgiesser et al. (2008) the mean recovery time period after a blood donation of 550 ml is 36 days (range 20–59). This average recovery time is added as a vertical arrow to Fig. 7. We run the optimization such that we keep the maximum EPO rate at 1000 U/day or after the first bleeding we increase to analyze by how much the recovery time could be decreased. In addition, we set and U/day after the first bleeding to see if this results in an overshoot.

Fig. 7

Bleeding patients 3 and 4. Up to the first bleeding equals 1000. The dotted lines mark the target range

We first analyze the results for a constant maximum EPO rate of 1000 U/day and when that rate is increased to 2000 U/day after the first bleeding; upper plots in Figs. 8 and 9. We observe that the NMPC algorithm is reacting to the bleeding by directly administering the maximum available EPO dose over a certain period of time. For U/day this higher amount is administered but for a shorter period of time. Hence, the total drug amount, see Table 2, is only 13% and 16% higher. The Hgb curves are shown in Fig. 7. The recovery time after the first bleeding for patient 4 can be significantly reduced from about 36 to about 25 days by allowing higher administration rates. For patient 3 the differences are minor. The reason is that a reduction of the recovery time for patient 3 comes along with an increase of the overshoot because the constant endogenous EPO concentration for patient 3 is around 280. For patient 4 it is only 43. Hence, by reducing the administration rate for patient 4 the controller can abruptly let drop down the EPO concentration in plasma which allows a fast recovery time without any overshoot. Even for the high maximum EPO rate of 25,000 U/day together with a minor control penalization () we do not observe an overshoot in patient 4. Note that the maximum rate is never reached; see Fig. 8c. But for patient 3 the overshoot gets more pronounced while the target range is reached some days earlier; Fig. 7b. Again, the maximum rate is not hit; Fig. 9c.

Fig. 8

Bleeding patient 4: EPO rates for different values after the first bleeding. Up to the first bleeding equals 1000. The first bleeding is at day 22 (7.5 g/d) and the second at day 80 (9.0 g/d)

Fig. 9

Bleeding patient 3: EPO rates for different values after the bleeding. Up to the bleeding equals 1000. The bleeding is at day 60 (7.5 g/d)

Table 2

Bleeding patients 3 and 4: total EPO doses for different values after the first bleeding

	Patient 4			Patient 3
u_max	1000	2000	25,000	1000	2000	25,000
Total dose (U)	79,124	91,853	187,660	83,980	94,948	199,077

Up to the first bleeding equals 1000. For patient 3 the first bleeding is at day 22 (7.5 g/d) and the second at day 80 (9.0 g/d). In case of patient 3 there is only one bleeding at day 60 (7.5 g/dl)

Missed administrations/dosing errors

In this section we consider a malfunction of the EPO pump. We simulate a complete failure to administer EPO for an entire day or that and incorrect rate is applied (not known in advance). In Fig. 10 we present the results for patient 1 with missed administrations on days 23–25, 35–39 and 110–114. In addition, on days 80–84 the maximum amount U/day gets administered by mistake. Missing an administration leads to a certain direct drop in the Hgb level. After every period of missed administrations the algorithm compensates for these by applying a higher EPO rate which is then gradually reduced in order to smoothly reincrease the Hgb concentration. After the period when the maximum rate is wrongly applied the algorithm administers nothing so that the Hgb level decreases and that an overshoot is avoided. But interestingly, it then restarts with the maximum amount before dropping on the desired level and a minor reincrease is even accepted. In doing so it achieves to balance the delayed effect of having administered no drug before and a later drastic drop-down is avoided. The total administered EPO dose is 50,853 U.

Fig. 10

Patient 1: optimal EPO rates and Hgb curve for missed treatments on days 23–25, 35–39, 110–114 and on days 80–84 the maximum EPO rate is administered by mistake

In Fig. 11 we present the results for periodically missed EPO administrations. More precisely, starting at day ten, every two weeks three days drop out. This scheme is chosen for analyzing if a periodic missing of administrations results in periodic administration rates. But as can be seen, the EPO rates after each period of missed treatments look different and the NMPC alogithm achieves to keep the Hgb level within the desired range if is set to 1500 U/day. Patient 5 demands comparatively high EPO rates so that the Hgb concentration falls below the target range if equals 1000 U/day. Note that the Hgb curve can be kept within the target range for if there are no missed treatments. The total administered EPO doses belonging to Fig. 11 are 167,848 () and 150,921 U/day ().

Fig. 11

Patient 5: optimal EPO rates and Hgb curve for missed treatments on days 10–12, 30–32, 50–52, 70–72, 90–92 and . In addition, the optimal Hgb curve for is shown in the upper plot

Constant EPO rates

Finally, we simulate different frequencies of rate change for the EPO pump. Note that less frequent adaptations could be a consequence of less frequent Hgb measurements. We investigate the effect of enlarging the period of a constant EPO rate from 1 day () to 1 week (), 2 weeks (), 3 weeks () and 4 weeks (). This requires to adjust the length of the NMPC horizon to 4 weeks (constant period 1 week), 6 weeks (constant period 2 or 3 weeks) and 8 weeks (constant period 4 weeks). In order to penalize control costs likewise for the different constant periods the parameter in Eq. (27) gets multiplied by the number of days in the respective constant period.

The results for patient 4 are exemplarily shown in Figs. 12 and 13. The Hgb curves for the different constant periods are all stabilized within the target range. The same holds true for the other patients. From a mathematical point of view a longer constant period means a restriction of degrees of freedom. Therefore, the Hgb level is best controllable for the shortest constant period. This is why the Hgb level is closest to the line of 10.5 g/dl when the controller can adjust the doses on a daily basis. And the longer the constant period is, the larger are the oscillations around this line. The total EPO doses, see Table 3, differ only little. The results for patient 4 look similar to those from patient 1 which is very interesting. Figure 5 can lead one to assume that the patient requires a very specific and time sensitive administration pattern to stabilize the Hgb level. But even for a constant period of a whole month the Hgb can be kept within the target range with negligible oscillations. For the other patients the oscillations are even smaller.

Fig. 12

Patient 4: optimal Hgb curves for different constant periods

Fig. 13

Patient 4: optimal EPO rates for different constant periods

Table 3

Patient 4: total EPO doses for different constant periods

Constant period	1 day	1 week	3 weeks	4 weeks
total dose (U)	55,598	56,964	56,721	57,695

Conclusion

The presented NMPC algorithm to correct anemia in HD patients was tested in various in-silico experiments and showed excellent performance in stabilizing simulated Hgb levels. The introduced framework uses a system of non-linear hyperbolic PDEs, which previously has been adapted to individual patients using clinically measured Hgb levels to predict individual patient response to EPO treatment (Fuertinger et al. 2018). The proposed NMPC would allow to optimize anemia treatment based on single patient data sets only. It would further allow to set individual Hgb targets for certain patient groups as proposed by the KDIGO work group. Thus, the presented work is a first step towards fully individualized anemia therapy.

The conducted in-silico experiments show that a fixed optimization setting for the controller scheme is sufficient to correct the anemia of patients with very different characteristics in the underlying prediction model. However, the penalization of the control costs needs to be tuned on a subgroup or even per-patient level to minimize the amount of administered EPO. The NMPC method requires a comparatively long horizon to account for the system’s large time delay. We have determined the required horizon length experimentally. Given this horizon length, the controller can handle the delay of the system response to treatment and achieves to stabilize Hgb levels even when presented with simulated events such as bleedings, missed administrations and EPO dosing errors. The presented controller has been tested under the assumption that treatment with EPO can be provided continuously. While this is an interesting in-silico experiment, such a therapy is currently clinically not possible as there are no “EPO pumps”, similar to the insulin pumps used for diabetes treatment, available. As a first restriction on the control we have investigated to allow changes in the EPO rate less frequently resulting in constant EPO rates over several weeks. These lead to slight oscillations of the total RBC population around the target state. However, for periods of up to four weeks the oscillations are negligible.

Different controller schemes to correct anemia in HD patients have been proposed over the last years by various groups (Barbieri et al. 2015, 2016; Brier and Gaweda 2011; Brier et al. 2010; Martínez-Martínez et al. 2014; McAllister 2017; Nichols et al. 2011). Some of them have been successfully tested in clinical trials (Brier and Gaweda 2011; Brier et al. 2010) or even been implemented in the clinical routine (Barbieri et al. 2015, 2016). Most of the proposed solutions to correct anemia are based on MPC techniques. In general, the time-varying nature of the process, long time delays in the system, high inter-individual variability in the specifics or the response of the system to treatment and the need to react to unforeseen events such as bleedings and missed doses renders the problem of correcting and stabilizing anemia in HD patients more suitable to MPC based controller approaches than PID (proportional-integral-derivative) ones. Previously presented and clinically validated MPC approaches have been based on neural networks as the underlying prediction model (Barbieri et al. 2015, 2016; Brier and Gaweda 2011; Brier et al. 2010). Such models require large training and validation sets to find the optimal weights for the neural network. The presented approach, where the underlying model is a system of coupled PDEs, although one of the more complex ones currently proposed, allows to optimize anemia treatment based on single patient data sets only. The required data from a single patient is routinely measured in clinics (as presented in Fuertinger et al. 2018). Thus, the proposed NMPC approach would provide a fully personalized anemia treatment strategy that even allows the setting of individual Hgb goals as suggested by the KDIGO guidelines. Further, estimated model parameters and their longitudinal development in individual patients potentially allow to gain further insights in the specifics of renal anemia. It might allow to better understand why some patients do not respond to treatment (e.g. short red blood cell life span versus insufficient bone marrow reaction to treatment).

In summary, the presented NMPC algorithm has the potential to bring more patients in the Hgb target range while decreasing Hgb variability and EPO utilization. However, we are still two major steps away from clinically testing the proposed NMPC approach: First, the control structure needs to be changed such that EPO is only administered during dialysis treatments (in general three times per week). With the chosen approach we will be able to analyze the effect of reducing administration times on Hgb stability. Second, the patient-model mismatch and uncertainty in parameter estimates together with measurement noise need to be addressed. In order to deal with parameter uncertainty in the underlying model, for instance, so called “robust MPC” methods would need to be incorporated into the framework. In addition, model estimates will need to be updated on a regular basis using the measured Hgb data. This is currently under investigation by the group but is beyond the scope of this publication.

Table 4

Fixed model parameters, values and units

Par.	Meaning	Value	Unit
β_1	Proliferation rate for BFU-E cells	0.6	1/day
β_2	Proliferation rate for CFU-E cells	1.2	1/day
β_3	Proliferation rate for erythroblasts	0.723	1/day
x¯_1	Maximal maturity for BFU-E cells	3	Days
x__2	Minimal maturity for CFU-E cells	3	Days
x¯_2	Maximal maturity for CFU-E cells	8	Days
x__3	Minimal maturity for erythroblasts	8	Days
x¯_3	Maximal maturity for erythroblasts	13	Days
x__4	Minimal maturity for marrow reticulocytes	13	Days
x¯_4	Maximal maturity for marrow reticulocytes	15.5	Days
α_4	Rate of ineffective erythropoiesis	0.09	1/day
α50	Intrinsic mortality rate for erythrocytes	0.002	1/day
Ω^_5	Age interval for erythrocytes, where	[14, 34]	Days
	neocytolysis is possible
ε	Regularization parameter for the mortality rate	10^-8	–
τ_E	EPO threshold for neocytolysis	80	mU/ml

Table 5

Fixed -parameters, values and units

Parameter	Meaning	Value	Unit
μ_1	Constant for the sigmoid apoptosis rate for CFU-E cells	0.5	1/day
μ_3	Constant for the sigmoid apoptosis rate for CFU-E cells	0.5	Dimension-less
μ_4,μ_5	Constants for the sigmoid maturation velocity for marrow reticulocyte	2, 0.35	Dimension-less
μ_7	Constant for the sigmoid maturation velocity for marrow reticulocytes	2.3	Dimension-less
μ_8	Constant in the mortality rate for erythrocytes	3.5·10^3	mU^3/ml^2
μ_9	Constant in the mortality rate for erythrocytes	3	mU^3/ml^2
μ_10	Constant in the mortality rate for erythrocytes	0.1	mU^3/ml^2

Table 6

Individualized model parameters and units

Parameter	Meaning	Unit
c_tbv	Total blood volume	ml
x¯_5	Maximal life span for erythrocytes	Days
S_0	Rate at which cells are committing to the erythroid lineage	1/day
E^end	Assumed constant endogenous EPO concentration in plasma	mU/ml
T_1/2	Half-life of Epoetin-α	1/day

Table 7

Individualized -parameters and units

Parameter	Meaning	Unit
μ_2	Constant for the sigmoid apoptosis rate for CFU-E cells	ml/mU
μ_6	Constant for the sigmoid maturation velocity for marrow reticulocytes	ml/mU