State and Parameter Estimation of a Mathematical Carcinoma Model under Chemotherapeutic Treatment

: One challenging aspect of therapy optimization and application of control algorithms in the ﬁeld of tumor growth modeling is the limited number of measurable physiological signals—state variables—and the knowledge of model parameters. A possible solution to provide such information is the application of observer or state estimator. One of the most widely applied estimators for nonlinear problems is the extended Kalman ﬁlter (EKF). In this study, a moving horizon estimation (MHE)-based observer is developed and compared to an optimized EKF. The observers utilize a third-order tumor growth model. The performance of the observers is tested on measurements gathered from a laboratory mice trial using chemotherapeutic drug. The proposed MHE is designed to be suitable for closed-loop applications and yields simultaneous state and parameter estimation.


Introduction
Cancer treatment and its related fields are intensively studied subjects, since cancer is a major cause of death globally. It is estimated that there will be 21.4 million cases annually by 2030 [1,2]. The application of modern control algorithms for cancer treatment has much potential-especially in the adjustment of the dosages-however, there are difficulties as well [3]. The article investigates a facet of a larger project. The larger project aims to minimize the amount of administered chemotherapeutic drug-while retaining its effect-by applying personalized treatments. Lower amount of doses expected to reduce side effects and potentially postpone the emergence of tumor resistance. One approach is the utilization of feedback control, where the controller (which decides about the drug administration scheme) adapts to the patient. Based on mathematical models accurately describing the underlying physiological processes the professionals have the possibility to define appropriate therapies. System engineers thrive to provide a proper description of the tumor growth, however, the inter-and intrapatient variability and the effects of different drugs are cumbersome to model. Several models were introduced over the years [4][5][6][7][8][9]. In particular, we are concerned with the use of a third-order model. The model describes the dynamics of the living and dead tumor cells, as well as drug concentration [7]. The output of the system is the total tumor volume (sum of the living and the dead tumor volumes).
In tumor growth modeling there are variables that cannot be measured, or are not feasible to measure in practice on a day-to-day basis. Such variables in this particular case are the individual volume of the living and dead tumor cells and drug concentration or the model parameters characterizing the patient variability. A mathematical model gives the general behavior of the system. The general structure can be tailored to each individual by the model parameters. For instance, there is

Experimental Data
The experiments were conducted by the Membrane Protein Research Group of the Research Centre for Natural Sciences published in [14]. The mice are identified with the "PLD" acronym, indicating that the injected chemotherapeutic drug was the pegylated liposomal doxorubicin (PLD) during the experiments. In this trial tumor pieces from a Brca1-/-;p53-/-mouse tumor were transplanted orthotopically to syngenic FVB mice to reduce the time of tumor formation. The tumor volume was measured by calipers every 1-5 days. The treatment protocol was a dynamic reactive approach which was defined as follows. The mice were treated based on the decision of the professionals. The injection could be granted if the tumor volume reached 200 [mm 3 ] and at least 10 days past since the last treatment. The injected drug was the maximum tolerable dose 8 [mg/kg] every time. The experiment focused on the elimination of the primary tumor. The experiment was stopped when the tumor volume grew above the 2000 [mm 3 ] threshold.
In Figure 1, an example is given of the trial. The raw caliper measurements are represented with black circles, the treatment threshold is plotted with a dashed grey line and the solid green line is a linear interpolation between the measurements. It can be seen that the particular mouse responded well to the treatment, since the volume decreased after the injections. However, some mice develop resistance and the drug stops taking effect (meaning continuously growing volume regardless of the presence of the drug). Those who developed resistance are not taken into account in this study as the applied model does not consider resistance. As a result, data from mice PLD1, PLD2 and PLD8 from [14] are not used in this study.

The Applied Tumor Growth Model
The investigated tumor growth model is a third-order nonlinear ordinary differential equation. The model differentiates between living and dead tumor cell volumes by introducing separate state variables: x 1 [mm 3 ] and x 2 [mm 3 ], respectively. The third variable x 3 [mg/kg] describes the drug concentration in the host. The dynamics is described by the equations given in [6,7] at time t aṡ where a describes the proliferation rate of the living tumor cells, n is the tumor cell necrosis rate, w is the dead tumor cell washout rate, c is the drug depletion rate. The pharmacodynamics is affected by the b and ED 50 (median effective dose) parameters. The effect of the drug and the depletion of the drug is characeterized by equations following the dynamics of Michaelis-Menten kinetics, thus ED 50 and K B are Michelis-Menten parameters. Resistance is not modeled explicitly, but the ineffectiveness of the drug appears as a specific combination of the parameters. When the difference a − b − n > 0, the drug is not able to shrink the volume, it only decreases the growth rate. Qualitative analysis of the model has been previously done in [15]. The values of the parameters for the mice identified based on a mixed-effect model [7] are given in Table 1.

Extended Kalman Filter
For nonlinear systems one of the most widely applied observer is the EKF. First, assume a general system given in the form of:ẋ where f describes the nonlinear system dynamics, h is the output equation, w is the additive process noise or disturbance and v is the measurement noise. Equation (4) is generalized further with the augmentation of system parameters x := [x s x p ] s.t.ẋ p = 0, where x s is the original state vector and x p is the parameter vector.
Since the measurements are done only at discrete times, continuous-discrete EKF has been applied. In the prediction step the model and covariances are propagated using the forward Euler method with shorter sampling times. The equations of the prediction step in time instant t are defined aṡ where P is the error covariance matrix, F is the Jacobian matrix of the system and Q is the process noise covariance matrix. The update step in the k-th discrete step is defined by where K is the Kalman gain, H is the Jacobian of the h output equation, z is the measurement, R is the measurement noise covariance matrix and the aposthrope k indicates the a priori values of the update step. The transition from the continuous time domain to the discrete is given by: The EKF structure has been implemented using the model (1)-(3), extended with parameter estimation. The estimated parameters were selected in accordance with a previous sensitivity and identifiability analysis [16]. The results of the analysis indicate that the most sensitive parameters are a, b, n andw, thus x p := [a, b, n, w] . The identifiability analysis with the subset of parameters resulted in a structurally locally identifiable [17] system. This indicates that there is a finite number of parameter sets producing the same trajectories. For the developments of the observers the following assumption is used: w = [w s w p ] s.t. w s = 0. Practically speaking, it is assumed that the mismatch between the measurements and the model rises only from measurement noise and parameter uncertainty. With the augmented state vector the Jacobians are given as: x 7 x 6 x 1 ED 50 (ED 50 +x 3 ) 2 0 x 1

Measurement Noise Characteristics
The noise characteristics of the measurements are quantified by calculating the standard deviation of the difference between the raw and filtered measurements. The filtered values are calculated by applying a second-order low-pass Butterworth filter in zero-phase setting. The resulting standard deviation σ = 31.39 is used for the R = σ 2 measurement variance. In Figure 2 the measured and filtered pairs are shown on a logarithmic scale to better visualize the deviations at small volumes. Zero measured values are shifted to 0.1 [mm 3 ], which was chosen arbitrarily for visualizing reasons.
The tumor volumes are measured using calipers and the actual volumes are approximated. This method is less accurate and reliable when smaller volumes are being measured. Data of measurements under 50 [mm 3 ] are rather scarce too, usually from a measured zero volume at the next sampling, the volume jumps instantly to the proximity of 50 [mm 3 ].

Kalman Filter Tuning
One critical aspect of comparing the performance of different observers is the tuning. In order to extract the most performance, the elements of the Q and P[0] covariance matrix were determined during an optimization. The goal of the optimization was to minimize the root-mean-square error (RMSE) between the measured trajectories and the predicted ones: where N is the number of treated mice, z is the vector of measured values of the i-th mouse and y is the output of the system. The optimized parameters were in the main diagonal corresponding to the estimated parameters and the first element of the P[0] matrix to guarantee quicker convergence. The optimization is done by the f mincon function of the MATLAB environment from 100 initial points generated by the Latin Hypercube Sampling [18]. The optimized matrices are given as:

Moving Horizon Estimation
The MHE is an optimization-based method, where the minimization of a developed cost function is done on a moving window. A moving window of one step can coincide with the widely applied EKF. However, using a longer window and introducing some application-specific terms or constraints can be beneficial in several cases. Taking into account more measurements can result in an estimator more robust to external disturbances or delays.
The cost function can be formulated in several ways, depending on the application. General arguments and terms are summarized in Table 2, where the notations are given according to (4) and (5). The W i matrices provide weighting between the terms of the cost function, their subscripts indicates the analogy between the EKF and MHE. The weighting matrices are analogous to the inverses of the covariance matrices (8), (10). The x ol vector function is the open-loop trajectory from the initial state of the current window, the arrival cost is calculated in the k−M-th step, where M is the length of the horizon. The difference between the current parameters and the previous window is denoted by ∆p. Following the assumption that the mismatch between the model and measurements are caused by parameter uncertainty, intrapatient variability and measurement noise, the parameters p = [a, b, n, w] are optimized, similarly to the EKF design. Given the arguments, the MHE problem can be generarly formulated as follows: Compared to the general formulation, three modifications are implemented. First, as a worst case scenario the initial state of the horizon is set to the last estimated state. Secondly, the differences are normalized in order to make the weighting between the measurements and modifications easier, scaling the errors to the same order of magnitude. Division with zero is avoided: when the tumor volumes could not be exactly measured due to their small volumes, 1 [mm 3 ] is assigned to their measured value. Unit volume ensures numerical stability and is physiologically more relevant as the tumors are not completely eliminated. The charactheristics of the model is in agreement with this, as in ideal condition the volume converges to zero, but it is bound to be positive. Thirdly, a rational function (23) was introduced, which penalizes the difference less, when the measured tumor volumes are small. Justification for the second and third point is more detialed in Section 2.3.1. Considering these modifications the cost function in the k-th step is defined as where T s is the sampling time of the measurements which is irregular; it varies between 1 and 5 days. This means that the sliding window shall also be irregular, so that when the new measurement is made, the window is adjusted to the time interval which is the closest to the predefined window length. For the minimization of the cost function, the inbuilt f mincon function of the MATLAB software was applied. In each step, the initial values for the optimization were the estimates of the previous step which is often called hot start. In the first step, the initial values of the parameters are set to the nominal identified values from Table 1. Furthermore, the implementation of the cost function is extended with stopping criteria, rendering the optimization quicker or feasible at all. The intrinsic property of the model is that the tumor can grow even in cases when it is under treatment. It can occur that the optimizer selects parameter combinations which results in unrealistically large tumor volumes, thus making the solution of the differential equations unstable. To avoid such situations the propagation of such trajectories are stopped and penalized. The stopping threshold is based on clinical protocol, since the experiment has to be stopped, when the tumor volume grows larger 2000 [mm 3 ]. The cost is defined to be dependent on the volume, since otherwise assigning a constant value or infinity can make the optimizer stuck due to constant gradients of the object function. The change of parameter a is not penalized as it can be seen in (22). The reason behind this is that the tumor growth is dependent on the ax 1 (t) term. In the case of small tumor volumes, the model could not describe quickly varying tumor growth, unless more sudden changes are allowed for the a parameter. In the case of non-estimated parameters the nominal values were assigned from Table 1.

Results
The observers were evaluated on the data described in Section 2.1. As mentioned each mouse is labeled with PLD followed by an ID number. Estimation horizons of the MHE were tested in the 7-35 days range, however for better visibility only the 14 day version is shown beside the EKF. Table 3 compares the RMSE of the measured and estimated volumes between the EKF and the investigated MHEs, given for each investigated mouse. When comparing the results two points must be emphasized. First, the assumptions made during the design of the observers. Secondly, the worst case scenario is investigated, in the sense that the EKF was tuned based on an optimization, where all the mice were taken into account. On the other hand, the weights of the MHE was not optimized only a penalizing function was introduced and applied instead of other tuning process. The RMSE indicates that for horizons shorter than 20 days the MHE can achieve better approximations on average. At horizon length of 28 days the results are varying from patient to patient. Only at 35 days the EKF has a clear edge over the MHE. Increasing the horizon of the MHE deteriorated the approximations as the actual parameter set had to describe longer sections, where structural mismatch and/or parameter variation occurs due to the change in the physiology of the tumor. At this point, it has to be emphasized again that while the MHE takes into account a series of measurements, the EKF works from sample to sample, thus various horizon lengths cannot be adapted.  Figure 3 depicts the comparison of the time horizons using different measures. Three measures, namely mean absolute error (MAE), "noise" and standard deviation (SD) of each estimated parameter were selected in order to assess the performance of the observers and horizons. The MAE describes the general accuracy with respect to the measurements, the "noise" term quantifies the smoothness of the estimation by calculating the number of direction changes along the trajectory. The SD of each estimated parameter gives information about the range in which the parameter had to vary in order to provide an optimized fit. One can inspect that the MAE increases in a quasi-linear way with respect to the horizon. The smoothest trajectory is achieved by using 14-day-long horizon and the SD of parameter a is the lowest around 21 days. Concerning the penalized parameters b, w, n the lowest SD is achieved by using the 7-day-long horizon. The main difference between the EKF and MHE is that the optimization of the EKF resulted in very small gains for parameters a and n. This is anticipated as those parameters directly affect the tumor growth, with larger gains it may be possible to achieve better approximation for one patient, but it could lead to bigger divergence for a population. The divergence is avoided by the MHE as the optimization is constrained, parameter sets producing trajectories that reach 2000 [mm 3 ] are neglected. In Figure 4, an example (PLD ID:5) is given about the effect of the rational penalizing function. The function was introduced in (21). The two estimated trajectories, denoted by with ω(·) and without ω(·) are the output of the MHE with 14 days of horizon length. For visualization only the 14 day scenario is given, however, the effect is similar in the other cases, and the delay is more pronounced towards the longer horizons. It can be seen that without this weighting in the case of small volumes there is a significant delay in the response of the observer. The dashed purple lines show the value of the function with respect to time (ω(t), left subfigure) and the measured volume (ω(x), right subfigure). In the following explanatory examples we address the problem of intrapatient variability, and how the observers were able to overcome it. Figures 5-14. showcase the state and parameter estimations separately. The black circles denote the measurements, however, it is important to note that only the total tumor volume is measured, i.e., y = x 1 + x 2 . The dashed magenta lines are the results from a Stochastic Approximation of Expectation-Maximization identification of the same model [6,7]. In the case of the x 3 , the three curves overlap, because parameters corresponding to (3) are not modified, hence nominal values were used. Those parameters are neglected in accordance with the previous sensitivity analysis [16]. Compared to the fixed identified parameter set by varying key parameters of the system, the observers are able to catch additional dynamics and track the measurements more accurately. This is particularly the case where the first administration of the chemotherapeutic drug made the tumor to shrink its volume close to zero. In fact, zero volumes were measured-namely, zero volumes are represented in the data-due to such small volumes cannot be measured accurately using calipers, albeit after a certain period of time, the tumor began to grow again. This phenomenon is observed in Figures 5 and 7, where the identified sets could not model the first growth period. In Figure 9 the first growth occurs, but the discrepancy between the measurements and the model is still large. The even-numbered Figures 6-10. show how the observers tackled the problem of tracking the first growth. It is particularly noticable in the case of the MHE, which started from a high growth rate (parameter a) and reduced it gradually to be able to maintain the tumor volume close to zero. Furthermore, it can be seen that the parameter a remains close to the identified value. However, at distinct time instances, greater divergence can be experienced. These divergences usually occur when the tumor exhibits a large growth rate. It can be observed that the identified models have a better fit when the tumor didn't shrink to zero volumes after initial growth. This phenomena can be observed in Figures 11 and 13.

Discussion
The goal was to evaluate the state and parameter estimation capabilities of an MHE and an EKF which has been applied as a benchmark. In general, state and parameter estimation in the case of physiological processes have major difficulties which are presented in this study related to tumor growth estimation as well. Furthermore, the problem formulation makes it even more difficult. In the model, only the sum of the living and dead tumor volumes can be measured, often inaccurately. This limited amount of data can easily lead to observability issues mostly coming from the fact that several state and parameter combinations may result in the same output. In order to limit these issues, we applied two key components: penalizing the velocity of changing of state variables and estimated parameters, moreover, we have selected a branch of key parameters to be estimated based on sensitivity analysis. Namely, the cost function of the MHE made it possible to use additional constraints and penalize the change of parameters. The penalization was designed to be close to the physiological rate of changes to be as as realisic as possible. According to the results, we have successfully limited the set of possible state and parameter sets leading to the same output without a downturn in the estimation accuracy.
The scenario for the comparison was designed to be biased arbitrary. The results show that even with optimizing the tuning of the EKF, and not taking into account the initial state of the MHE as an additional free variable, a moving horizon estimation based observer can still outperform the former EKF. The MHE-based observers achieved lower RMSE until the horizon length was shorter than 30 days. Using horizons around 20 days, there can be cases, where the EKF outperforms the MHE, e.g., the estimation of PLD9. At 7 days the difference in RMSE is around 30 [mm 3 ], and the MHE has lower RMSE for each mouse. The inferior performance of the EKF can rise from the low gains for parameters a and n. Since the system is most sensitive to those parameters, having larger gains could result in divergence in the case of certain patients due to the high interpatient variability. The estimation accuracy of the MHE deteriorated in a quasilinear manner with respect to the window length. The increased window lengths resulted in an increased RMSE of the tumor volumes because the effects of the structural mismatch became dominant. In general, a longer window would mean that the given parameter set is valid for more data. However, the predictive capabilities are crucial in therapy optimization, thus the determination of the optimal choice for closed-loop applications is future work. The runtime of the observers are greatly different. The EKF averages 2.6 ms, while the MHE ranges from 300 ms to 413 ms depending on the horizon length. However, runtime is not a cornerstone as the injections and measurements are done once a day at a maximum rate. An important aspect of the study is that the observers are evaluated on real laboratory data, where two major difficulties are present: structural mismatch and intrapatient variability.
The results indicate the great importance of applying observers for therapy optimization algorithms. Without constant feedback about the behavior of the patient the model using only a single parameter set can introduce large discrepancies. As the tumor shifts phases, and as the heterogenity of the tumor changes its behavior changes as well. These shifts in physiology can be tracked with the observers-particularly well with the MHE-. These changes can be interpreted as trends in the parameters. Additionally, despite the optimized tuning of the EKF the MHE shows a smoother estimation of the parameters and a more accurate approximation of the measurements. These attributes facilitate the application of the MHE in closed-loop algorithms.

Conclusions
In this paper, an MHE-based observer and an EKF was designed and compared for therapy optimization and state-feedback kind of control applications. When using time horizons less than 28 days, the MHE based observer had superior performance-compared to the EKF-in tracking the measurements. Increased horizon length deteriorated the approximations as the actual parameter set had to describe longer sections, where structural mismatch and/or parameter variation occurs due to the change in the physiology of the tumor. To evaluate the performance more thoroughly advanced measurement technology is needed to distinguish live and dead tumor cells and even drug concentration. Further improvements can be done on the fine-tuning of the parameters in the cost function. Future work includes the integration of the MHE into closed-loop systems.