Drug dosing for cancer therapy: A stochastic model predictive control perspective

Contents lists available at ScienceDirect Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/jtbi Drug dosing for cancer therapy: A stochastic model predictive control perspective Andrés Hernández-Rivera a,∗, Pablo Velarde b, Ascensión Zafra-Cabeza a, José M. Maestre a aDepartment of System and Automation Engineering, University of Seville, Spain bEngineering Deparment, Universidad Loyola Andalucía, Spain article i n f o Keywords: Cancer Chemotherapy Non-linear control systems Model predictive control Stochastic processes a b s t r a c t Stochastic Model Predictive Control (SMPC) is an effective decision-making method in applications where uncertainties play a significant role. This work introduces a non-linear formulation of SMPC specifically designed for cancer therapy. The proposed method considers the stochastic nature of tumor growth, non-linear dynamics, and a potential side effect of the treatment. Through one-year simulations, the results showcase the effectiveness of this strategy in controlling drug dosing. 1. Introduction Cancer is a leading cause of death worldwide and is responsible for a significant burden on health, economies, and social structures (Spanish Society of Oncology, 2020). Some forecasts estimate that between 40 % and 50 % of people could develop cancer during their lifetimes in the next century (Smittenaar et al., 2016). Cancer patients often require expensive and time-consuming treatments that reduce their quality of life. Addressing this disease requires a comprehensive and coordinated approach involving healthcare professionals, policymakers, and researchers to mitigate its impact on individuals and society. Cancer treatments depend on the type of tumor, its stage, and the patient’s condition. Among the best-known are radiotherapy (Allen et al., 2017), immunotherapy (Cappuccio et al., 2007), surgical interventions, and chemotherapy (Gustafson and Page, 2013), which is non-selective and has undesirable secondary effects that require careful consideration to minimize harm to the patient. Another common characteristic of cancer drugs is that they are usually administered following relatively static dosing patterns during treatment and may be adjusted based on side effects and the evolution of the patient. One key way to improve the efficacy of these treatments relies on developing new mathematical models of the biological systems at hand. In particular, the model of a drug’s pharmacokinetics considers its interaction with the tumor and its secondary effects; see, e.g., Malinzi et al. (2021), where different examples of mathematical models are reviewed. ∗Corresponding author at: University of Seville, Camino de los Descubrimientos, 41092, Seville, Spain. E-mail addresses: [email protected] (A. Hernández-Rivera), [email protected] (P. Velarde), [email protected] (A. Zafra-Cabeza), [email protected] (J.M. Maestre). In this regard, models provide helpful insight to reduce a drug’s toxic effects and, therefore, can be leveraged to improve its efficacy. Likewise, the use of engineering and mathematical techniques to address cancer treatment methods is rapidly developing, e.g., to predict the evolution of tumors (Ghaffari Laleh et al., 2022; Sherwin, 2019) and their average behavior in a stochastic manner (Sharpe and Dobrovolny, 2021). By using mathematical models, control algorithms may improve drug administration, as shown in works like (Liliopoulos and Stavrakakis, 2021), in which linear quadratic regulators for improving a chemotherapy-based oncological treatment are used. Simultaneously, secondary effects can be diminished if adaptive interventions are applied (Zafra-Cabeza et al., 2010). In this study, we focus on Model Predictive Control (MPC), a popular control strategy for dynamic systems, because of its capacity to deal with non-linearities, delays, and constraints on problem variables, among other reasons. MPC uses a mathematical model to predict the future evolution of the system variables and optimize a control objective subject to constraints (Camacho and Bordons, 2013). This feature makes it suitable for safety-critical applications where it is essential to guarantee that the system remains within certain operating bounds. For example, the optimal load of oncolytic virus therapy, which consists of the administration of genetically modified viruses to reduce the number of tumor cells, is optimized using MPC in Villa-Tamayo et al. (2022), Parihar et al. (2022); a non-linear MPC-based (NMPC) dosing of a combined treatment of chemotherapy and immunotherapy is presented in Chareyron and Alamir (2009); an NMPC in combination with a https://doi.org/10.1016/j.jtbi.2025.112255 Received 29 April 2024; Received in revised form 5 August 2025; Accepted 23 August 2025 Journal of Theoretical Biology 615 (2025) 112255 Available online 31 August 2025 0022-5193/© 2025 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). A. Hernández-Rivera et al. moving-horizon estimator to attain an optimal drug administration scheme is proposed in Czakó et al. (2021); and the optimal thermal dose used in the ultrasonic heating treatment of cancer is calculated using MPC in Hensley et al. (2015), Arora et al. (2002). In this context, it is essential to consider the role of uncertainties when designing and implementing the control strategy. They play a crucial role in cancer therapy because the system’s evolution is affected by non-modeled dynamics, contributing to the variability in therapy outcomes. Other sources of disturbance come from measurement errors in the data acquisition process and the dynamics of the patient, which can vary over time. Consequently, the decision-making process must incorporate stochastic variables to address these uncertainties. However, the classical formulation of MPC does not allow for the consideration of systems with uncertainties, even if some MPC schemes have been proposed to ensure stability and compliance with constraints in the presence of disturbances (Bernardini and Bemporad, 2009; Velarde et al., 2017). One way to address this problem is using the conservative min-max approach (Hans et al., 2014). A less conservative strategy is the Stochastic MPC (SMPC), which designs predictive controllers for dynamic systems subject to disturbances/uncertainties regarding the probability that a particular solution is feasible. Therefore, it is not strictly possible to speak about guaranteed feasibility in this context (Grosso et al., 2017). Scenario-based MPC is one of the options to handle uncertainties in which a single control sequence is calculated for each considered scenario, see, e.g., Schildbach and Morari (2015), Calafiore and Campi (2006). However, this approach can be conservative because extreme scenarios may affect the controller’s decisions (Velarde et al., 2017). To avoid this issue, other strategies have been considered, such as the socalled chance-constrained MPC (CC-MPC), which can also be built based on either scenarios or an explicit characterization of the uncertainty (Jurado et al., 2016; Schildbach and Morari, 2015). CC-MPC uses probabilistic modeling of system uncertainties to calculate explicit bounds on system constraint satisfaction. Moreover, CC-MPC offers advantages such as robustness, flexibility, low computational requirements, and the possibility of including the level of reliability associated with the constraints (Grosso et al., 2014; Schwarm and Nikolaou, 1999). Since CC-MPC considers the expected performance of the closed loop with probabilistic constraints instead of directly trying to assure robust constraint satisfaction, it avoids the conservativeness present in other robust MPC techniques, e.g., Scokaert and Mayne (1998), Langson et al. (2004). In this sense, the effectiveness of this stochastic formulation has been emphasized in several works in the biomedical field, such as Jurado et al. (2016), Liu et al. (2019). Within SMPC algorithms, Chance Constrained-NMPC (CC-NMPC) stands out for handling a non-linear model of the system and the probabilistic constraints involved in the optimization problem by using equivalent deterministic constraints regarding risk violation. This work applies a CC-NMPC, where variables such as the tumoral volume, the number of lymphocytes, and the final metabolites of the chemotherapy have been modeled as stochastic variables. The non-linear mathematical model of tumoral growth presented in Chen et al. (2012) has been adopted to enhance predictions of future disease behavior. The model explicitly considers tumoral growth in mice, with Tamoxifen (TM) being the administered drug. The main contribution of the work presented here lies in pursuing an enhanced treatment performance that limits the side effects of TM and achieves the highest safe reduction in tumor size through implementing a CC-NMPC controller. By employing a CC-NMPC strategy, the safety of the treatment is potentially improved, even in the presence of uncertainties. This increased robustness is accomplished by adapting the constraints to effectively handle such disturbances, enabling flexible management of the stochasticity associated with a biological system. This work presents several features that distinguish it from others, including the following highlights: •An assertive application of CC-NMPC has been made in cancer therapy. This approach is an optimal and stochastic treatment strategy that is more robust and less conservative than others. It also provides an effective way to overcome the uncertainties and complexities of cancer progression and treatment response. •The specific characteristics of non-linear systems in cancer therapy pose unique challenges. However, considering their stochastic nature, this work demonstrates how the CC-NMPC controller can be effectively customized for these systems. This adaptation is essential for successfully implementing control strategies in such complex environments and represents a significant contribution to this field. •This manuscript connects our theoretical findings and their practical applications. By demonstrating CC-NMPC’s effectiveness in a realworld context, we provide crucial insights for practitioners in the field of oncology, paving the way for more effective and personalized cancer treatments. The rest of the document is organized as follows. Section 2 describes the non-linear mathematical model of tumor behavior, the pharmacokinetics of the drug used in cancer treatment, and its side effects. Section 3 develops the complete formulation of the CC-NMPC to deal with uncertainties. To highlight the benefits of using a CC-NMPC strategy, a reliable comparison between the designed controller and other treatment schemes is carried out in Section 4. Finally, Section 5 includes some conclusions and future directions. 2. Mathematical model of the system The mathematical model used in this work considers several aspects related to the evolution and elimination of a generic tumor in a mouse: the evolution of the drug’s pharmacokinetics, the tumor’s growth, and the possible side effects the treatment might cause. As stated before, the model’s equations and parameters are derived from previous research in Chen et al. (2012), which, in turn, was derived from the works of Florian et al. (2008), de Pillis et al. (2006). The model corresponds to a saturating-rate cell-cycle model, which is key to studying the effects of cycle-specific agents like TM. It is also important to note that this model assumes a homogeneous distribution of cells within the tumor. Moreover, as stated in Chen et al. (2012), the model parameters were determined based on data taken from mice. Firstly, the evolution of the volume of cancerous cells is modeled through the use of three stages of cellular division: 𝑋𝑔 is the volume of cells in the growing state, 𝑋𝑠 corresponds to the volume of cells in the DNA synthesis state, and 𝑋𝑚 represents the volume of cells in the mitosis state. This cellular division model combines the quiescence and growth phases into the new 𝑋𝑔, while the mitotic preparation and mitosis are combined into 𝑋𝑚. These states are modeled as follows:  𝑋𝑔(𝑡) = −𝑘𝑔𝑋𝑔(𝑡) ln (Θ 𝑁(𝑡))+ 2𝑘𝑚𝑋𝑚(𝑡) ln (Θ 𝑁(𝑡)) −𝑘𝑑𝑋𝑔(𝑡)(𝑋2 𝑉+𝑐𝑋3 𝑉),(1)  𝑋𝑠(𝑡) = −𝑘𝑠𝑋𝑠(𝑡) + 𝑘𝑔𝑋𝑔(𝑡) ln (Θ 𝑁(𝑡)),(2)  𝑋𝑚(𝑡) = −𝑘𝑚𝑋𝑚(𝑡) ln (Θ 𝑁(𝑡))+𝑘𝑠𝑋𝑠(𝑡),(3) In addition, the volume of the tumor is represented by 𝑁(𝑡) = 𝑋𝑔(𝑡) + 𝑋𝑠(𝑡) + 𝑋𝑚(𝑡).(4) This tumor growth model exhibits a saturating-rate growth behavior that captures the deceleration of growth as the tumor approaches its plateau population. Additionally, using a cell-cycle model allows for the simulation of cycle-specific anticancer drugs (Chen et al., 2012). Moreover, the tumoral growth model developed in Chen et al. (2012), Florian et al. (2008), de Pillis et al. (2006) assumes that cells in the 𝑋𝑔 state (which include both quiescent and early proliferative cells Journal of Theoretical Biology 615 (2025) 112255 2 A. Hernández-Rivera et al. Table 1 Parameter for the tumor growth. Parameter Value Unit Description 𝑘𝑔0.0013 hour−1 Transfer coefficient 𝑋𝑔 to 𝑋𝑠 𝑘𝑠0.0390 hour−1 Transfer coefficient 𝑋𝑠 to 𝑋𝑚 𝑘𝑚0.0169 hour−1 Transfer coefficient 𝑋𝑚 to 𝑋𝑔 𝑘𝑑0.0062 hour−1 Cancerous cells death coefficient V8.592 ml Total blood volume c25 – Efficiency coefficient for TM 𝜃104mm3 Plateau volume of cancerous cells Table 2 Parameter for the TM pharmacokinetics. Parameter Value Unit Description 𝑘01 0.048 hour−1 Transfer coefficient from state 𝑋0 to 𝑋1 𝑘12 0.993 hour−1 Transfer coefficient from state 𝑋1 to 𝑋2 𝑘23 35.932 hour−1 Transfer coefficient from state 𝑋2 to 𝑋3 𝑘𝑟21.145 hour−1 Chemotherapy consumption coefficient in 𝑋2 𝑘𝑟339.525 hour−1 Chemotherapy consumption coefficient in 𝑋3 Table 3 Parameter for the evolution of the lymphocytes. Parameter Value Unit Description 𝛼𝐶1.21 ⋅105hour−1 Natural generation coefficient of lymphocytes 𝛽𝐶1.2⋅10−2 hour−1 Natural death coefficient of lymphocytes 𝑘𝑐0.010 ml⋅µg⋅hour−1 Chemo-induced lymphocytes death coefficient 𝑏25 – Efficiency coefficient for TM due to experimental grouping limitations) are predominantly sensitive to Tamoxifen. Although some heterogeneity in sensitivity may exist, this simplification is considered adequate for capturing the primary drug effect, as reported in the aforementioned experimental studies (Table 1). The pharmacokinetics of TM is described through the following equations:  𝑋0(𝑡) = −𝑘01𝑋0(𝑡) + 𝑢𝑐(𝑡),(5)  𝑋1(𝑡) = −𝑘12𝑋1(𝑡) + 𝑘01𝑋0(𝑡),(6)  𝑋2(𝑡) = −𝑘𝑟2𝑋2(𝑡) − 𝑘23𝑋2(𝑡) + 𝑘12𝑋1(𝑡),(7)  𝑋3(𝑡) = −𝑘𝑟3𝑋3(𝑡) + 𝑘23𝑋2(𝑡).(8) Here, 𝑢𝑐(𝑡) is the model’s input variable, and it represents the administered dose of TM, which should be below the maximum daily amount of 800 µg. The metabolization process of this drug is divided into four stages: 𝑋0, 𝑋1, 𝑋2, and 𝑋3, where the last stage denotes the concentration of 4-hidroxitamoxifen (Table 2). This mathematical model that portrays the pharmacokinetic evolution of TM is a linear time-invariant (LTI) system with first-order ODEs. In the case of a null input variable (𝑢𝑐(𝑡) = 0), the resulting system matrix has all negative eigenvalues, which ensures the stability of this linear subsystem. Moreover, when simulating an open-loop daily administration of TM, the values of each compartment (𝑋0, 𝑋1, 𝑋2, and 𝑋3) present an oscillatory dynamic that eventually converges to a stable, daily periodic orbit. If the drug administration is suddenly halted, each concentration gradually decays back to zero, which further reinforces this point. Furthermore, the model considers lymphocytes as an indicator of the chemotherapy-induced side effects of TM. In particular, they indicate the degree of degradation of the mouse’s immune system, and it can be modeled as follows:  𝐶(𝑡) = 𝛼𝐶−𝛽𝐶𝐶(𝑡) − 𝑘𝐶𝐶(𝑡)(𝑋2(𝑡) 𝑉+𝑏𝑋3(𝑡) 𝑉),(9) where 𝐶(𝑡) represents the evolution of the number of lymphocytes. It is important to note that, as a safety measure, the remaining lymphocytes should always be higher than 40 % of the initial amount (Chen et al., 2012). Therefore, the control strategy must aim to achieve the fastest tumor size reduction that allows for a safe treatment, ensuring that this lymphocyte safety threshold is never crossed (Table 3). Finally, Fig. 1 represents the proposed model and the added uncertainties that have been considered. As can be seen, the two final stages of the metabolization of the chemotherapy (𝑋2 and 𝑋3), the total tumor size (𝑁), and the number of lymphocytes (𝐶) are considered as the model’s output variables. Moreover, the input variable is the dose of TM (𝑢𝑐). Fig. 1 also shows the relationships among all system subsystems. This mathematical model has been designed to be updated hourly. These parameters’ values and detailed meanings can be found in Chen et al. (2012). 3. The control approach applied to tumor growth through chemotherapy The non-linear mathematical model presented in Section 2 has been used to simulate the evolution of the tumor-immune system in the mouse. It has been discretized using the Backward Euler method with a sample time of one hour. Moreover, a control-oriented version of the model is implemented in the controllers by adding the sources of uncertainties to the states and outputs. The system can be defined, for each time step 𝑘∈ℤ+, as 𝑥[𝑘+ 1] = 𝑓(𝑥[𝑘], 𝑢[𝑘], 𝑥[𝑘], 𝑡),(10) 𝑦[𝑘] = 𝑔(𝑥[𝑘], 𝑢[𝑘], 𝑡) + 𝑦[𝑘],(11) where 𝑢[𝑘] stands for the input variable, the amount of dosage of TM (𝑢[𝑘] = 𝑢𝑐[𝑘]). In this non-linear system, the collection of state variables is 𝑥[𝑘], and the system outputs are 𝑦[𝑘]. Furthermore, 𝑥[𝑘] represents the experimental subject-prediction model discrepancy (Chen et al., 2012), and the measurement disturbances associated with each output variable are represented by 𝑦[𝑘]. The state variables are defined as: 𝑥= [𝑋𝑔𝑋𝑠𝑋𝑚𝐶 𝑋0𝑋1𝑋2𝑋3]𝑇. The concentrations 𝑋2 and 𝑋3 as well as the total tumor volume 𝑁 and the number of lymphocytes 𝐶 are considered to be the output variables of the system (subject to the observational noise 𝑦[𝑘]). In terms of system constraints, the maximum allowable dose of TM is 800 µg according to Chen et al. (2012), that is, 0≤𝑢𝑐[𝑘]≤800 𝜇𝑔. (12) Moreover, the values of the state variables, based on Chen et al. (2012), can be constrained as: 𝑥[𝑘]≥[0004⋅1060000]𝑇(13) This article proposes using an MPC strategy to reduce the tumor size while limiting the degree of secondary effects derived from the treatment. The cost function that can address this problem is defined as a quadratic function. This function introduces a penalty that increases more severely with the magnitude of the deviation from the setpoint. This feature is particularly advantageous when it is essential to maintain proximity to the references, and more significant deviations are increasingly undesirable, i.e., 𝐽(𝑦[𝑘], 𝑢[𝑘]) = (𝑦ref −𝑦[𝑘])𝑇𝑅𝑦(𝑦ref −𝑦[𝑘]) + 𝑢[𝑘]𝑇𝑄𝑢𝑢[𝑘],(14) where 𝑦ref is the reference vector for the output variables (0 mm3 for 𝑁[𝑘], 0 𝜇𝑔∕𝑚𝑙 for both 𝑋2 and 𝑋3 and 107 for 𝐶[𝑘]). Matrices 𝑅𝑦 and 𝑄𝑢 are the weight factors for each error signal (100, 10−5, 10−5, and 1 for each output variable) and the input signal (0.1), respectively. The established weights aim to find a trade-off between removing the tumor and minimizing the adverse effects of administering the drug to the subject. The control scheme that has been implemented is depicted in Fig. 2, illustrating how it receives the output variables’ measurements and considers their references, constraints, and the risk of violation. This scheme enables the system to achieve optimal TM dosing for regulating tumoral growth. Journal of Theoretical Biology 615 (2025) 112255 3 A. Hernández-Rivera et al. This work develops two MPC-based strategies to improve treatment by considering the model’s non-linear properties. On the one hand, a standard NMPC approach is applied, while on the other, a CC-NMPC, an SMPC based on chance constraints, is formulated. They are then reliably compared. Finally, it is important to note that if a different mathematical model were used, leading to a different set of input/state/output variables, the controller’s formulation would need to be adapted accordingly. 3.1. Standard NMPC formulation The NMPC controller minimizes the stage cost function along a prediction horizon (𝑁p) at each time instant 𝑘, i.e., min 𝑢[𝑘∶𝑘+𝑁p−1] 𝑘+𝑁p−1 ∑ 𝑙=𝑘 𝐽(𝑦[𝑙], 𝑢[𝑙]),(15) subject to (10)–(13), ∀𝑙∈ [𝑘, 𝑘 +𝑁P− 1]. Moreover, this problem is solved using the 𝑓𝑚𝑖𝑛𝑐𝑜𝑛 function in MATLAB 2023b. As a result of solving the optimization problem, the vector 𝑢[𝑘∶𝑘+ 𝑁p− 1] = {𝑢[𝑘], 𝑢[𝑘+ 1],…, 𝑢[𝑘+𝑁p− 1]} is computed. However, only the first component, 𝑢[𝑘], is applied to the system, while the remaining elements are discarded. Consequently, the optimization problem (15) is repeated at the next step 𝑘+ 1 in a receding horizon fashion. As mentioned earlier, the output variables of the system are assumed to be normally distributed random variables. This implies that the 𝑖th output variable, denoted as 𝑦𝑖=(𝑦𝑖, 𝜎2 𝑦𝑖), has an average value of 𝑦𝑖 and a standard deviation of 𝜎𝑖. This disturbance results from measurement noise that the NMPC can handle. Moreover, the state variables are affected by process uncertainty, as there is a 1% normally distributed mismatch between the experimental subject and the mathematical model applied in each time step, as pointed out in Chen et al. (2012). As a consequence, the 𝑗th state variables can be expressed as 𝑥𝑗=(𝑥𝑗, 𝜎2 𝑥𝑗). All in all, the presence of both sources of uncertainties Fig. 1. Schematic of a mouse’s model of tumoral growth. Fig. 2. Closed-loop control scheme employed to regulate tumoral growth. Journal of Theoretical Biology 615 (2025) 112255 4 A. Hernández-Rivera et al. implies the realization of the stochastic variable, for both the outputs and states of the system, according to a known probability distribution for each time instant of the prediction horizon. The prediction horizon is set to one week (𝑁p= 7 days), providing a comprehensive view of the system’s behavior over a substantial period, which is crucial for anticipating and effectively managing longer-term dynamics and variations. Considering the trade-offs between robustness, practicality, and computational burden, the sampling period has been established to one day. 3.2. CC-NMPC formulation The proposed CC-NMPC combines the philosophy and benefits of an MPC controller with the probabilistic constraints necessary to deal with the process uncertainties present in this case, as evidenced by the Eqs. (10) and (11). That is, by making changes and assuming a particular risk of violation, constraints that are affected by dynamic disturbances can be rewritten in a probabilistic manner: ℙ[𝑥𝑗[𝑘]≥𝑥𝑗,min]≥1 − 𝛿𝑥, where ℙ[⋅] denotes the probability operator and (1 − 𝛿𝑥) is the given likelihood with which the constraints have to be fulfilled, i.e., 𝛿𝑥 represents the risk violation index. Remark 1. It is important to note that the parameter 𝛿𝑥 quantifies the acceptable risk of violating system constraints, ranging from 0 (no constraint violation allowed) to 1 (fully permissive). In the context of cancer therapy, a lower 𝛿𝑥 leads to a more cautious treatment strategy, which results in higher lymphocyte preservation but less tumor reduction, while a higher 𝛿𝑥 permits more aggressive tumor control at the potential cost of immune suppression. Therefore, 𝛿𝑥 offers a tunable parameter for tailoring treatment aggressiveness based on patient-specific factors. The dynamical uncertainty impacting each state variable 𝑥𝑗, where 𝑗∈ {1,2,…,8}, may be described using its cumulative distribution function (cdf), 𝜙𝑗. Individual chance constraints are applied as follows, ℙ[𝑥𝑗[𝑘]≥𝑥𝑗,min]≥1 − 𝛿𝑥⇔ ℙ[𝑓𝑗(𝑥[𝑘], 𝑢[𝑘], 𝑥[𝑘], 𝑡)≥𝑥min𝑗]≥1 − 𝛿𝑥⇔ ℙ[𝑓𝑗(𝑥[𝑘], 𝑢[𝑘], 𝑥[𝑘], 𝑡) − 𝑥𝑗[𝑘] 𝜎𝑥𝑗 < 𝑥min𝑗−𝑥𝑗[𝑘] 𝜎𝑥𝑗]< 𝛿𝑥⇔ 𝜙𝑗(𝑥min𝑗−𝑥𝑗[𝑘] 𝜎𝑥𝑗)< 𝛿𝑥⇔ 𝑥min𝑗−𝑥𝑗[𝑘] 𝜎𝑥𝑗 < 𝜙−1 𝑗(𝛿𝑥)⇔ 𝑥min𝑗−𝑥𝑗[𝑘]< 𝜎𝑥𝑗𝜙−1 𝑗(𝛿𝑥)⇔ −𝑥𝑗[𝑘]<−𝑥min𝑗+𝜎𝑥𝑗𝜙−1 𝑗(𝛿𝑥)⇔ 𝑥𝑗[𝑘]> 𝑥min𝑗−𝜎𝑥𝑗𝜙−1 𝑗(𝛿𝑥).(16) As mentioned in Section 3.1, stochastic variables 𝑥𝑗 are treated as normal distribution functions. Based on this perspective, individual chance constraints for each time step along the prediction horizon, 𝑙∈ [𝑘, 𝑘 +𝑁𝑝− 1], are written as: 𝑥𝑗[𝑙]> 𝑥min𝑗−𝜎𝑥𝑗𝜙−1 𝑗,𝑙 (𝛿𝑥),(17) where 𝜙−1 𝑗,𝑙 is the inverse cdf of the uncertainty of the state variable 𝑥𝑗 at the time instant 𝑙, which can be determined by using the prediction of the model and its statistical properties. The deterministic equivalent of the initial probabilistic constraint is derived based on the state variables presented in Section 1. The optimization problem formulated under a CC-NMPC approach can be expressed as follows: min 𝑢[𝑘∶𝑘+𝑁p−1] 𝑘+𝑁p−1 ∑ 𝑙=𝑘 𝐽(𝑦[𝑙], 𝑢[𝑙]),(18) subject to ∀𝑙∈ [𝑘, 𝑘 +𝑁p− 1], (10)–(14) and (17). Finally, only the first element of the resulting vector 𝑢[𝑘∶𝑘+𝑁p− 1] is then applied to the system. One significant advantage of this method over other stochastic MPC control strategies, such as tree-based or multiscenario MPCs, lies in its capability to reformulate the stochastic optimization problem into its deterministic equivalent through Eq. (17). Furthermore, in this case, there is a single realization of the normally distributed process uncertainty variable for each time instant along the prediction horizon. 3.3. Analysis of the statistical properties of the process uncertainty. The non-linearity of the mathematical model makes it necessary to define the behavior of its statistical properties and check how the process uncertainties are propagated along 𝑁p. As mentioned above, these dynamical uncertainties are assumed to be modeled as normal distributions. Had the mathematical model of the system been linear, this analysis would not have been necessary, as the propagation of process uncertainties would always have a normal distribution. However, due to the non-linearities presented in Section 2, the normalcy of the distribution of the propagation cannot be assured without analysis. Algorithm 1: Statistical properties analysis of the disturbances. 𝑘←1; while 𝑘≤365 do Solve NMPC and calculate the sequence of 𝑢[𝑘∶𝑘+𝑁p− 1]; while 𝑙∈ [𝑘, 𝑘 +𝑁p− 1] do for 𝑖= 1 ∶ 5000 do Apply the corresponding element of the sequence 𝑢[𝑘∶𝑘+𝑁p− 1] to day 𝑙 and collect the values of each state variable (𝑋𝑔, 𝑋𝑠, 𝑋𝑚, 𝐶, 𝑋0, 𝑋1, 𝑋2, and 𝑋3), according the control oriented model. end Obtain a suitable probability distribution of the state variables for the day 𝑙; 𝑙=𝑙+ 1; end 𝑘=𝑘+ 1; end For that, Algorithm 1 outlines the analysis of the statistical properties that can be conducted for each day of the treatment 𝑘. First, the sequence of 𝑢[𝑘∶𝑘+𝑁p− 1] is calculated using the NMPC. Afterward, for each day of the prediction horizon ∀𝑙∈ [𝑘, 𝑘 +𝑁p− 1], with 𝑁p= 7 days, the cdf of the disturbance for each state variable 𝑥𝑗 of the non-linear system is obtained. This evaluation is performed through a Monte Carlo analysis simulating 5000 cycles using the control-oriented model. As mentioned in the previous subsection, there are process uncertainties modeled as a 1 % normally distributed mismatch between the mathematical model employed for the predictions and the mouse (Chen et al., 2012). In our case and for the sake of simplicity, only the results for the state 𝐶 on day 𝑘=𝑙= 30 are represented in Fig. 3, where it is possible to observe that this state variable can be modeled as a normal distribution function. In addition to the previous analysis, two goodness-of-fit tests were conducted to verify the normal distribution of uncertainty propagation. Journal of Theoretical Biology 615 (2025) 112255 5 A. Hernández-Rivera et al. Fig. 3. Monte Carlo analysis for different days along the 𝑁𝑝. Table 4 Test results from Monte Carlo analysis along 𝑁p (𝛼= 0.01). Jarque-Bera 𝑋𝑔𝐶 𝑋2𝑋3 Days ℎ 𝑝 ℎ 𝑝 ℎ 𝑝 ℎ 𝑝 1 0 0.73534 0 0.38773 0 0.36993 0 0.42872 2 0 0.15741 0 0.69733 0 0.12957 0 0.11687 3 0 0.82486 0 0.56168 0 0.68550 0 0.79908 4 0 0.48033 0 0.55447 0 0.61687 0 0.71323 5 0 0.63429 0 0.18231 0 0.52491 0 0.46107 6 0 0.46131 0 0.63716 0 0.85877 0 0.93300 7 0 0.57450 0 0.81423 0 0.61863 0 0.62440 Shapiro-Wilk 𝑋𝑔𝐶 𝑋2𝑋3 Days ℎ 𝑝 ℎ 𝑝 ℎ 𝑝 ℎ 𝑝 1 0 0.82564 0 0.18287 0 0.91963 0 0.93092 2 0 0.19989 0 0.63980 0 0.13065 0 0.12309 3 0 0.50789 0 0.91607 0 0.44112 0 0.56271 4 0 0.48162 0 0.88210 0 0.87286 0 0.69443 5 0 0.37347 0 0.24492 0 0.71450 0 0.62246 6 0 0.59944 0 0.84169 0 0.86233 0 0.77266 7 0 0.78615 0 0.46773 0 0.59873 0 0.71535 Particularly, Jarque-Bera and Shapiro-Wilk (BenSaïda, 2014) tests were used, having the following null hypothesis: H0: The variable under study is normally distributed with unspecified mean value and standard deviation. Additionally, and for simplicity, this analysis has been performed for every state variable. However, Table 4 only displays the results for the states: 𝑋𝑔, 𝐶, 𝑋2, and 𝑋3. Analogously to Fig. 3, these tests were conducted for the entire prediction horizon, 𝑁p, starting on day 30 of the simulation. Moreover, 𝛼 denotes the significance level used for each test. The variable ℎ indicates whether the null hypothesis can be rejected based on the test results, with a zero value indicating that the null hypothesis cannot be dismissed. The p-value (𝑝) represents the probability of observing the test results, assuming that the null hypothesis is true. It becomes evident that each state variable shown consistently exhibits a normal distribution. This observation provides further support for the validity of the distributional assumptions made in the analysis. It confirms that the uncertainty propagation aligns well with a normal distribution over 𝑁p. Moreover, repeating these two tests for a CC-NMPC control sequence results in the same outcome, where the normalcy of the distribution is maintained. Finally, it is essential to note again that this analysis was necessary to check if the normally distributed process uncertainties propagated, maintaining the normality of their distribution. Furthermore, the results presented in Fig. 3 and Table 4 demonstrate that the uncertainty propagation maintains multi-variate normality. The statistical tests applied in these analyses can be further employed to confirm that the multivariate distribution continues to approximate normality throughout the prediction horizon. 4. Results and discussion To highlight the strengths and weaknesses of both NMPC-based controllers, three one-year simulation studies have been designed to compare different TM dosing strategies. The first one consists of a heuristic controller that provides the maximum dose (800 𝜇 g) when the number of lymphocytes (𝐶) is above its constraint and 0 µg in any other case. Alongside that scheme, the NMPC and the CC-NMPC were implemented with a risk of constraints violation of 𝛿𝑥= 10%. Both types of controllers use the non-linear model of the system, which provides more accurate predictions for the behavior of the tumor, given a set of initial values for each variable: 950 mm3 for 𝑋𝑔, 50 mm3 for 𝑋𝑠, 50 mm3 for 𝑋𝑚, 107 for 𝐶, and 0 µg/ml for 𝑋0, 𝑋1, 𝑋2, and 𝑋3, respectively. Furthermore, the values of the standard deviation for the measurement disturbances associated to outputs 𝑦=[𝑋2𝑋3𝑁 𝐶]𝑇, 𝜎1= 0.001 𝜇𝑔∕𝑚𝑙, 𝜎2= 0.001 𝜇𝑔∕𝑚𝑙, Journal of Theoretical Biology 615 (2025) 112255 6 A. Hernández-Rivera et al. Table 5 Comparison among the different controllers employing KPIs. Controller 𝐾𝑃 𝐼1( 𝐾𝑃 𝐼1) mm3𝐾𝑃 𝐼2( 𝐾𝑃 𝐼2) mg 𝐾𝑃 𝐼3( 𝐾𝑃 𝐼3) Heuristic 148.85 (9.95) 230.80 (1.23) 3.97 ⋅106(6 ⋅104) NMPC 149.31 (10.25) 230.25 (1.33) 3.98 ⋅106(7 ⋅104) CC-NMPC 175.62 (10.21) 214.52 (1.36) 4.18 ⋅106(7 ⋅104) Controller 𝐾𝑃 𝐼4( 𝐾𝑃 𝐼4) mm3∕mg 𝐾𝑃 𝐼5( 𝐾𝑃 𝐼5) h−1 𝐾𝑃 𝐼6( 𝐾𝑃 𝐼6) % Heuristic 3.688 (0.046) 6.95 ⋅104(4.9⋅103) 46.04 (2.32) NMPC 3.695 (0.047) 4.32 ⋅104(5.9⋅103) 44.26 (3.23) CC-NMPC 3.843 (0.050) 2.76 ⋅103(1.88 ⋅103) 5.01 (2.09) 𝜎3= 20 mm3, and 𝜎4= 105, were extracted from Chen et al. (2012), Florian et al. (2008), de Pillis et al. (2006). 4.1. Performance evaluation To compare the results and evaluate the advantages of the heuristic, NMPC, and CC-NMPC control strategies, several Key Performance Indicators (KPIs) have been introduced to statistically analyze the different controllers’ behavior and efficacy, accounting for the system’s stochasticity. As mentioned above, the mathematical model that represents the evolution of the mouse is discretized so that it updates hourly, while all of these control strategies calculate a daily dose of TM. •KPI1: Final size of the tumor. •KPI2: Total dose of TM. •KPI3: Final number of lymphocytes. •KPI4: Ratio between the volume of 𝑁(𝑡) eliminated per mg of TM. •KPI5: Number of lymphocytes below constraint per hour. •KPI6: Percentage of hours occurring in constraint violation. Table 5 provides each controller’s average values and standard deviations for each KPI obtained from the 1000 simulations performed. This statistical analysis allows for the comparison of the performances of the three controllers (heuristic, deterministic NMPC, and CC-NMPC). It is worth noting how the heuristic approach achieves the highest reduction of the tumor (KPI1), albeit with a more significant consumption of TM (KPI2) and a greater vulnerability to the uncertainties, as shown by the lower values of KPI3. The NMPC achieves a slightly safer approach than the previous strategy, as evidenced by KPI5 and KPI6, while final tumor size is minimally larger (KPI1). However, both these controllers (heuristic and NMPC) design chemotherapeutic cycles that might have severe side effects. The CCNMPC controller achieves a slightly larger final tumor size while providing a safer, more efficient approach, which is reflected in the lower consumption of the chemotherapy drug (KPI2) and larger ratio of eliminated tumor per mg of TM (KPI4). The safety level for this last treatment is evidenced by KPI3 and KPI5, as its mean value and standard deviation fall outside the constraint violation area. It is essential to highlight how, for the CC-NMPC-defined treatment, the risk of constraint violation (𝛿𝑥= 10%) is met, as the value of KPI6 for this case has an average value of 5.01 %. To further substantiate the observed performance differences between the control strategies, particularly highlighting the superior performance of CC-NMPC, two-sample t-tests were conducted. These tests aimed to assess the statistical significance of differences in drug efficiency (KPI4) and constraint violation (KPI6) between CC-NMPC and each of the other two controllers: NMPC and the heuristic approach. Firstly, a t-test was performed to evaluate the null hypothesis (𝐻0,1) that no significant difference exists in the drug efficiency. Secondly, another t-test examined the null hypothesis (𝐻0,2) that no significant difference exists in the percentage of constraint violations. The results of these tests are presented in Table 6, where the p-value indicates the likelihood of obtaining the observed results assuming the null hypothesis is valid (with an 𝛼= 0.05 for both tests). On the other Table 6 Results for t-student tests comparing the controllers for KPI4 and KPI6. Tests KPI4 (𝐻0,1) KPI6 (𝐻0,2) h p h p Heuristic vs CC-NMPC 12.51 ⋅10−8 12.92 ⋅10−12 NMPC vs CC-NMPC 16.43 ⋅10−7 12.93 ⋅10−12 hand, the h-statistic shows the test result, where h = 1 means the null hypothesis is rejected (significant result), and h = 0 means the test fails to reject it (no significant difference). The extremely low p-values across all four tests provide strong statistical evidence of a substantial difference in drug efficiency and constraint violation. These results reinforce the conclusions drawn from the data in Table 5. To illustrate the evolution of the system for one particular treatment, Fig. 4 shows the size of the tumor, 𝑁, the administered 𝑢𝑐, and the reduction in 𝐶 (log-scale Y-axis) throughout a one-year treatment for both NMPC. To further study the results for each control strategy, an analysis of the constraint violation has also been carried out in Fig. 5, where the values of the evolution of the lymphocytes along 100 one-year treatments have been represented. Each value that does not meet the system constraints is shown in red. Fig. 5a shows how the heuristic approach strategy struggles to keep the secondary effects of the TM treatment within safe levels. Furthermore, the absence of a stochastic formulation in the NMPC can induce constraint violations in the system’s closed-loop behavior. This is evidenced by how the CC-NMPC significantly reduces the constraint violations compared to the NMPC, substantially enhancing the safety of the treatment and the odds that the subject might complete it satisfactorily. Based on the gathered information and the results obtained, it can be concluded that the CC-NMPC dosing scheme outperforms the NMPC-based strategy in terms of its behavior. The former demonstrates a reliable capability to avoid constraint violations, even in the presence of uncertainties. Remark 2. The CC-NMPC’s consideration of Tamoxifen pharmacokinetics, including its decay rate, helps prevent the algorithm from producing a near-continuous therapy schedule, ensuring more controlled and effective drug delivery. This achievement is attributed to the CC-NMPC controller’s ability to reformulate probabilistic constraints over the time horizon 𝑁p, considering the statistical properties of stochastic variables and allowing for controlled risk violations through the parameter 𝛿𝑥. Finally, this controller was also implemented with two other control intervals: (a) a TM administration every two days and (b) a weekly drug treatment. Keeping the maximum TM dose at 800 µg leads to a subpar tumor reduction in case (a), 351.24 mm3, even when the maximum dose is administered every two days for the whole year. For case (b), that maximum dosage every 7 days cannot achieve a tumor reduction, leading to its rejection. Therefore, the daily administration presented in Table 5 and Fig. 4b is an acceptable trade-off between quality of control and treatment-associated burden to the patient. 4.2. Sensitivity analysis and robustness testing of the controller The performance of the CC-NMPC depends on several factors, including the choice of the weight matrix in the cost function, variations in model parameters, and the underlying system dynamics. This section analyzes the controller’s sensitivity to different weight matrix selections and evaluates its robustness across a range of model parameters through a systematic parameter sweep. Furthermore, the controller is implemented on alternative mathematical models to assess its generalizability beyond the system presented in Section 2. These analyses Journal of Theoretical Biology 615 (2025) 112255 7 A. Hernández-Rivera et al. Fig. 4. Evolution of the treatment by using SMPC controllers. provide insights into the adaptability and reliability of the controller under varying conditions. 4.2.1. Weight sensitivity analysis The results presented in Section 4.1 were calculated using a selection of weight factors in Eq. (14): 100, 10−5, 10−5, and 1 for each output variable (𝑁, 𝑋2, 𝑋3, and 𝐶), and 0.1 for the input signal (𝑢𝑐). As mentioned above, these weights were chosen as a trade-off between reducing tumor volume and minimizing the drug’s adverse effects through a trial-and-error process. One of the strengths of stochastic model predictive control is its ability to find an optimal solution within the given constraints, even in the presence of uncertainties in parameter values or unmodeled dynamics, which is a common problem in biomedical mathematical modeling. However, a sensitivity analysis using a Monte Carlo approach was conducted to further justify these values’ selection. In this analysis, the three primary cost function weights were modified one at a time, varying their orders of magnitude while keeping the remaining values of 𝑅𝑦 and 𝑄𝑢 constant. The modified weights were: the weight associated with tumor volume 𝑁 (referred to as 𝑅𝑦,1), the weight associated with the number of lymphocytes 𝐶 (𝑅𝑦,4), and the weight related to the drug dose 𝑢𝑐. The two weights associated with the other output variables (𝑋2 and 𝑋3) were kept at 10−5, as their influence on the system’s final evolution is minimal. This means that when the influence of 𝑅𝑦,1 was tested, the rest of the values of 𝑅𝑦 and 𝑄𝑢 were held constant and set to the original values presented at the beginning of this section. Fig. 6 presents the results of this analysis (represented in box plot format), from which several key insights can be drawn. First, the higher 𝑅𝑦,1 is, the drug administration cycle becomes relatively more aggressive, leading to smaller final tumor volumes, lower numbers of lymphocytes, and a more significant percentage of constraint violations. Journal of Theoretical Biology 615 (2025) 112255 8 A. Hernández-Rivera et al. Fig. 5. Constraint compliance map for each case. Second, the weights associated with 𝐶 and 𝑢𝑐(𝑡) have the opposite effect: a larger value induces a more conservative approach, resulting in larger final tumor sizes, a higher number of lymphocytes, and a lower percentage of constraint violations. These results indicate how the original weights in Eq. (14) were in the range of values that perform a good trade-off between tumor elimination and subject safety. 4.2.2. Model parameter sweep A model parameter sweep was also performed to demonstrate the robustness of the controller. In this analysis, the controller was tested on 1000 simulations, where the values of the parameters of the system’s mathematical model, presented in Section 2, were randomly selected at the beginning of each simulation. This selection was based on the standard deviations for each parameter reported in Florian et al. (2008). The test results presented in Table 7 lead to two key conclusions. Firstly, the results exhibit more significant variability than those in Table 5 for the same controller (CC-NMPC), a consequence of the diverse system dynamics introduced by the parameter sweep. These findings underscore the benefits of personalized treatment. Specifically, by adjusting the controller’s weight parameters based on expert input from healthcare professionals, treatment can be tailored to each patient’s condition, optimizing both efficacy and safety. Ultimately, automating chemotherapy cycle design aims to equip medical staff Table 7 KPIs for the parameter sweep. Controller 𝐾𝑃 𝐼1( 𝐾𝑃 𝐼1) mm3𝐾𝑃 𝐼2( 𝐾𝑃 𝐼2) mg 𝐾𝑃 𝐼3( 𝐾𝑃 𝐼3) CC-NMPC 197.89 (45.22) 218.33 (12.56) 4.11 ⋅106(7 ⋅104) Controller 𝐾𝑃 𝐼4( 𝐾𝑃 𝐼4) mm3∕mg 𝐾𝑃 𝐼5( 𝐾𝑃 𝐼5) h−1 𝐾𝑃 𝐼6( 𝐾𝑃 𝐼6) % CC-NMPC 3.757 (0.39) 7.49 ⋅104(7.25 ⋅103) 7.37 (3.15) with advanced tools that enhance decision-making and improve patient outcomes. Secondly, despite the increased variability, 𝐾𝑃 𝐼6 shows that the CCNMPC maintains the percentage of constraint violations below the 10 % threshold set during the controller’s design. This further highlights the importance of developing a control strategy that explicitly accounts for the stochastic nature of the system it regulates. 4.2.3. CC-NMPC with other mathematical models and treatments As mentioned above, the controller has been implemented on three alternative mathematical models to evaluate its generalizability. These models were selected from the literature due to their distinct approaches to studying tumor growth. Notably, two of them incorporate alternative chemotherapeutic strategies, such as combining chemotherapy with immunotherapy or antiangiogenic treatments. Journal of Theoretical Biology 615 (2025) 112255 9