Nonlinear Model Predictive Control Using Robust Fixed Point Transformation-Based Phenomena for Controlling Tumor Growth

In this paper a novel control strategy is introduced in order to create optimal dosage profiles for individualized cancer treatment. This approach uses Nonlinear Model Predictive Control to construct optimal dosage protocols in conjunction with Robust Fixed Point Transformations which hinders the negative effect of inherent model uncertainties and measurement disturbances. The results are validated by extensive simulation on the proposed control algorithm from which conclusions were drawn.


Introduction
In the past decades, physiological systems gained attention among control engineers.This particular interest covers a significant variety of topics including automated anesthesia [1], diabetes control [2] and hormonal regulation [3].Tumor growth control is no exception as cancerous diseases are responsible for 1,359,500 deaths in the European Union annually, based on [4].Because of the physiological complexity of tumor growth, the topic has been intact for decades.However in the end of the 20th century, progress in cancer research delivered mathematical models that are able to encapsulate the underlying mechanisms.An important result was Targeted Molecular Therapies (TMTs) which exploit certain cancer specific processes that corresponds to the vascluar growth of tumors, thus hampering its effect in the body of the patient.This approach led to the famous Hahnfeldt model, which was introduced by [5], and has been used thoroughly since its development among control engineers.
The primary issue with TMTs is that the developed drugs are exceedingly expensive, thus availability for the vast majority is limited.A possible solution to this issue is to apply optimal control techniques to the existing tumor growth model that incorporates the expenses of the drug while maintaining a reasonable treatment plan i.e., minimizing the costs and the treatment time.Several algorithms were proposed in the literature based on linear design methods [6] and nonlinear techniques [7].Many evidence [8,9] and practical consideration, such as unstable behaviour of the linearised system and significant discrepancies in terms of model parameters between different patients, suggest that nonlinear methods have to be employed.In spite of these facts, a novel approach is presented in this paper, which combines the optimal behaviour of the Nonlinear Model Predictive Control (NMPC) in conjunction with the error insensitive traits of the Robust Fixed Point Transformation (RFPT) based nonlinear adaptive control.In a previous work, the NMPC architecture was deployed in order to tackle the optimization issue [10].However as it was indicated in [11], the states of the system are not measurable due to large parameter uncertainties and model imperfections; therefore, an adaptive approach is suggested.The problem with the RFPT approach is that it can not guarantee an optimal dosing protocol by itself.While it is locally robust, it can not track constant references which are far from the initial value.In spite of these facts, a combined approach is presented here, which may also be applicable for other physiological systems as well that has significant parameter variations in the nominal model.
The paper first introduces the underlying model of the system and the important theoretical background of both the NMPC and the RFPT method.After the introduction, the connection is presented between the two methods with a detailed explanation of the problems that has to be solved from a control engineering aspect.In the end of the paper, numerical simulations are presented in order to verify the behaviour of the controller, on which conclusions are drawn.

System Model
Since anti-angiogenic therapy has emerged, many approaches have been considered by various authors.A comprehensive study on tumor growth models can be found in [12] from a control engineering point of view.However, advancements in cancer research have carried out new results, which indicate that the currently existing models do not incorporate important alternate vascularization methods, as pointed out in [13].In order to overcome these issues, different models were developed by [14,15] that describes the growth process accurately, which was validated by mice experiments.However, in this article, the famous Hahnfeldt model [5] is used, which does not cover every physiological aspect, but rather facilitates the comparison between the presented control algorithm and previously employed approaches.
Based on the work of [6], many simplifications can be attained on the original Hahnfeldt model, such as elimination of terms that are shown to be negligible by experiments [5], or merging variables that tend to move together [16].The modified equations of motion can be described as: where x 1 denotes the tumor volume (mm 3 ), which is the output of the system, x 2 is the volume of the tumor vasculature (mm 3 ), λ is the growth parameter of the tumor (1/day), b is the angiogenic factor (1/day), d describes the cellular blocking mechanisms of the vasculature (1/day•mm 2 ), e is the inhibition of the vasculature by the drug (kg/day•mg), and g(t) is the concentration of the administered inhibitor (mg/kg).The parameters were chosen according to [17], which are presented in Table 1.It is worth mentioning that the underlying system model describes a positive dynamical system, which entails that the control signal must remain positive for all times.Previous simulations showed that when the input signal is zero, the model attains a steady state of x 1 = x 2 = 1.734 × 10 4 (mm 3 ).This steady state yields a worst case scenario, so that the controller is designed according to this initial value in order to deal with minor initial tumor volumes as well.Because anti-angiogenic therapy does not eradicate the tumor completely, the objective is to shrink the volume to such an extent that it does not jeopardise the health of the patient.The model contains a major singularity at x 1 = x 2 = 0 (mm 3 ) that corresponds to the previously indicated behaviour of the therapy, which implies that the control objective can not be zero value.In particular the tumor is said to be in safe steady state if its volume does not exceed 10 (mm 3 ); hence from the controllers perspective it is sufficient to bring the tumor volume in this region.

Control Algorithm
In this section, the combined NMPC-RFPT approach is introduced with the necessary design steps for the simplified Hahnfeldt model as a primary example.The idea shares similar traits of the robust NMPC approach that can be seen in [18]; however, the second optimization loop is replaced by the RFPT controller.This decreases the computational burden of the dual NMPC solution and makes the design steps easier, as it does not require to construct terminal cost and sets to ensure a robust operation of the controller.This technique also overcome several obstacles concerning the RFPT method.The NMPC block provides an optimal trajectory for set point tracking, as an output of the nominal model, that can be used by the RFPT controller.This solves the set point tracking problem of the RFPT controller by providing an appropriate nominal trajectory prescription, in conjunction with the non-optimal behaviour of the RFPT technique.If the plant and the nominal model has the same parameter values, the RFPT controller can be tuned to produce minor tracking error, thus the plant trajectory does not deviate significantly from the nominal, which entails that the output is essentially optimal.If uncertainties are present in the plant, then the RFPT controller tries to reduce these deviations which renders the control loop robust.It should be also noted that the nominal trajectory, which is provided by the NMPC block, can be computed offline which makes the method computationally attractive.This means in practice that instead of solving a complex nonlinear optimization in each control cycle, one step of iteration (7) has to be computed which is faster.

The Nonlinear Model Predictive Controller
The NMPC consists of three main ingredients, namely, prediction, optimization, and execution.The controller uses the plant measurement in order to predict the behaviour of the system based on the discrete model, which is then used in a cost function.This essentially means the construction of a complicated cost function that has to be minimized by the input sequence for each prediction step, which implies the use of nonlinear optimization solvers.From the optimal control sequence, the first element is applied to the plant and the cycle repeats.The operation of the controller can be summarized in the Optimal Control Problem (OCP N n ), based on [19] as: where x 0 = x(n) is the measurement of states in the nth control cycle, N denotes the length of the optimization horizon, the value u(•) ∈ U N (x 0 ) represents a vector of predicted input values of dimension N from the feasible set U N , function is the cost function, J N is the total cost function that has to be minimized, and ) is responsible for the prediction stage with initial condition x u (0, x 0 ) = x 0 .Note that the prediction equation describes a discrete time system in the form of x + = f (x, u).This implies in particular that the continuous time system must be discretized, therefore the system model takes the form of: in which ∆t is the time step of the discretization.In (2), the left side of the prediction corresponds to the state variables in (3), namely T , and the function f (x u (k, x 0 ), u(k)) is the vector valued function on the right-hand side of the discretization equation.It is worth mentioning that different discretization methods can be employed which could lead to more accurate results, thus improving the overall robustness of the NMPC controller in exchange for computational effort.
According to [19], the cost function can be defined as where ζ and ξ are control parameters, x 1 is the tumor volume, x re f (k) is the time-varying reference signal with u re f (k) as the desired input signal at prediction time k, and • is the Euclidean norm.Equation (4) allows the use of time varying references; however, in this case a constant reference is proposed, so x re f (k) ∈ [1; 10] (because of the singularity in the model at x 1 = x 2 = 0), and u re f (k) = 0, ∀k ∈ N.For technical reasons, in every control cycle, the calculated predicted values u(•) is used in the next step by the nonlinear optimization algorithm as an initial guess.

The Robust Fixed Point Transformations Based Controller
The RFPT method is a novel nonlinear control technique which was introduced by [20] and proven to be effective in many control problems providing accurate results.Its operation can be described with the realized-response scheme.Theoretically, if one calculates the input signal for an arbitrary trajectory by inverting the dynamics of the model, then it is possible to steer the system along this trajectory by applying the calculated input signal to the plant, if there are no uncertainties between the nominal and physical systems (this will be called the desired response).However, in the vast majority of cases this is not true due to modelling imperfections, or parameter variations.This causes a realized response which differs significantly from the desired one, thus the control objective can not be achieved.The RFPT method solves this issue by deforming the desired input trajectory, so that the calculated input signal produces the desired response of the plant.A more elaborate description can be found in [20][21][22].
The method first uses a kinematic part which can be completely determined by the order of the system, that is the highest derivative of the states in the system model.The kinematic equation is given by: e int (t) = t t 0 (q n (τ) − q(τ)) dτ , where q n is the reference trajectory, q is the measured state, n is the order of the system, and Λ is a control parameter.Since the highest derivative is n = 1 and the controlled variable is x 1 , the expression takes the form of: In ( 6), x n 1 denotes the reference tumor volume that is a decreasing function with respect to time.In the combined approach this will be the output of the nominal system model of the NMPC part.After the kinematic part, the output ẋd 1 is used by the deform function which is defined as: whence the iteration is achieved by r n+1 = G r n |r d .In (7), f (r) is the measured signal, r d is the output of the kinematic block that is ẋd 1 , r d * is the solution of the fixed point problem for which the realized response of the system coincides with the desired trajectory, and A, B, K are free control parameters.Using the deformed trajectory, the output can be calculated from the inverse model, which now has to be expressed.Taking the first equation from (1) and rearranging the terms leads to: Differentiation of the previous equation yields: from which the following form can be obtained by simplification: In the simplified equation f denotes the exponential term in (9).This equation can be substituted into the second equation in the system model ( 1), together with Equation (8).By rearranging the terms, g(t) can be expressed as: which is the inverse model of the system.Applying the output of the deform function to the inverse model leads to a proper control signal g(t), which drives the system along the desired trajectory.

The Combined Approach
The last step in the design process is to determine the parameters of the controllers.A possible path is to tune the NMPC controller separately to the nominal plant by adjusting the values of ξ and ζ, so that it meets the design criterion, or using other more sophisticated techniques from the NMPC literature.If the NMPC controller is tuned, the output of the nominal system can be used as an input of the RFPT controller so that parameters A, K and B can be determined by numerical simulations on the nominal plant.Because of its novelty, the RFPT method has no exact tuning techniques yet.However, several thumb rules exist which alleviate the tuning process; B is either 1 or −1, K is usually a big number with the magnitude of 10 5 , and A = 1/(10 • K).A schematic depiction of the controller architecture can be seen in Figure 1.
Version October 9, 2018 submitted to Machines ) 2 (9) from which the following form can be obtained by simplification: In the simplified equation f denotes the exponential term in (9).This equation can be substituted into the second 140 equation in the system model ( 1), together with equation (8).By rearranging the terms, g(t) can be expressed 141 as: (11) which is the inverse model of the system.Applying the output of the deform function to the inverse model leads 143 to a proper control signal g(t), which drives the system along the desired trajectory.The last step in the design process is to determine the parameters of the controllers.A possible path is to 146 tune the NMPC controller separately to the nominal plant by adjusting the values of ξ and ζ , so that it meets 147 the design criterion, or using other more sophisticated techniques from the NMPC literature.If the NMPC 148 controller is tuned, the output of the nominal system can be used as an input of the RFPT controller so that 149 parameters A, K and B can be determined by numerical simulations on the nominal plant.Because of its novelty, 150 the RFPT method has no exact tuning techniques yet.However, several thumb rules exist which alleviates the 151 tuning process; B is either 1 or −1, K is usually a big number with the magnitude of 10 5 , and A = 1/(10 • K).

152
A schematic depiction of the controller architecture can be seen on Fig. 1.

Simulation Results
Multiple simulations were conducted in order to ensure the appropriate operation of the controller.In order to mimic real behaviour, the NMPC controller operates on the plant once every day with zero order hold.The prediction time step was 1/24 (day) with a horizon of N = 24 × 14 = 336, which predicts the evolution of the system for two weeks.Control parameters ζ = 100, ξ = 1 produced acceptable results in terms of tumor volume reduction and total inhibitor concentration, with a maximum input signal restriction of 25 (mg/kg).For this nominal trajectory, the parameters of the RFPT controller was set as K = 5.3 × 10 10 , A = 1/(10 11 ), B = −1, by empirical tuning, which produced sufficiently small tracking errors.The simulation step was ∆t = 10 −3 (day), with time interval t ∈ [0; 100].
On Figure 2 the nominal trajectory is presented which was created by the NMPC controller.After the 100th day the tumor volume was 1.083 (mm 3 ), which is close to the desired setpoint value of x 1 = 1 (mm 3 ).This result can be improved by increasing the prediction horizon; however, it affects the computational time significantly, thus it was not applied here.The tumor reduces under 10 (mm 3 ) volume in 31 (days), with an inhibitor concentration of 765.7 (mg/kg).Compared to the results in [10] this is clearly an improve in terms of concentration during the transient phase, i.e., the amount of medication given to the patient until the tumor volume reaches 10 (mm 3 ).Figure 3 illustrates the tracking error of the RFPT controller on the nominal trajectory.It can be seen that the error decreases to zero almost in ten days while its values in the transient phase does not exceed 15 (mm 3 ), which can be considered negligible.The difference between the two control signal can be scrutinised on Figure 4.It can be seen that the RFPT input signal matches its NMPC counterpart almost everywhere, thus it does not violate the bounds on the control input.It should be noted, however, that there is a jump at the beginning of the control cycle which can be attributed to numerical imprecisions.The reduction of tumor and vasculature volume can be seen on Figure 5.In order to ensure the robustness of the controller, different system parameter configurations were simulated.The parameters were chosen with Latin Hypercube Sampling from the interval [0; 1], which was then scaled to the order of magnitude for each parameter and added to their nominal values.In total fifty different parameter configuration was tested, which revealed that the system is very sensitive for parameter d, but for the other three the controller is able to manage even extreme cases as well.On Figure 6 one can see the tracking errors which are caused by the parameter uncertainties between the nominal and real models.The maximal error has increased significantly, but the error vanishes in ten days and the maxima is negligible compared to the initial value of the tumor volume.Simulations also revealed that the robustness can be increased by adjusting the value of Λ in the kinematic part; hence Λ was altered to Λ = 0.5.On Figure 7 the different administration protocols are shown.It is worth noting that their shape resembles the nominal case in Figure 4 and also their values do not violate the 25 (mg/kg) constraint prescription.

Conclusions
In this paper, a novel control method was introduced that uses the NMPC and RFPT methods in order to achieve robust optimal control of an anti-angiogenic tumor growth model.A design procedure was introduced, so that the approach can be used in different control problems that require a robust optimal control.Simulations showed that the proposed algorithm can tackle many issues posed by parameter uncertainties and modelling imperfections.The algorithm combined with previous results on discrete control with large sampling time of the Hahnfeldt model could result in a final solution to the individualized anti-angiogenic problem.Further research has to focus on auto tuning of the RFPT controller and also the possibility to include a feedback part for the NMPC part of the controller.

Figure 2 .
Figure 2. The nominal trajectory produced by the NMPC controller.

Figure 3 .
Figure 3. Tracking error of the RFPT controller.

Figure 4 .
Figure 4. Administration protocols of the RFPT and the NMPC controller.

Figure 5 .
Figure 5.The reduction of the tumor volume.