## 2 Materials and Methods
### Personalized Chemotherapy
#### 2.6.3 Optimization Problem  
The aim of cancer treatment is to minimize the tumor cell population \(T(t)\), maintain a healthy population of natural killer cells \(N(t)\) while minimizing the total amount of chemotherapy administered during the treatment to reduce its adverse effect on healthy cells. The objective of designing the chemotherapy dosing protocol is to ensure the transition of the state variables from their initial values \(N(0),\ L(0),\ T(0)\) to the coexisting state \(N_{0},L_{0},T_{0}\), which represents cancer remission [59]. The therapy is implemented as a finite sequence of \(n_{d}\) individual (discrete) doses \(D_{i}\), delivered periodically as short, uniformly spaced impulses (boluses) with a period of \(T_{D}\) days. Each dose is represented as a rectangular pulse of the administration rate signal \(v(t)\), lasting \(0.1\times T_{D}\) and having an amplitude of \(D_{i}\times\frac{10}{T_{D}}\), ensuring that the pulse's area under the curve equals \(D_{i}\). The aim of the design is to find vector \(\mathbf{D}\in\mathbb{R}^{n_{D}}\) of \(n_{D}\) positive doses  
\[\mathbf{D}= \begin{pmatrix}D_{1}&D_{2}&D_{3}&\cdots&D_{n_{D}}\end{pmatrix} ^{\text{T}}, \tag{34}\]  
minimizing the criterion  
\[ \begin{array}{c}J_{D}(\mathbf{D})=\sum_{k=1}^{n_{D}\times n}[w_{N}\big{(}N\big{(} kT_{s},\mathbf{D},\mathbf{\hat{\theta}}\big{)}-N_{0}\big{)}^{2}+w_{L}\big{(}L\big{(}kT_{s}, \mathbf{D},\mathbf{\hat{\theta}}\big{)}-L_{0}\big{)}^{2}+\\ w_{T}\big{(}T\big{(}kT_{s},\mathbf{D},\mathbf{\hat{\theta}}\big{)}-T_{0}\big{)}^{2}+w_{ u}u\big{(}kT_{s},\mathbf{D},\mathbf{\hat{\theta}}\big{)}^{2}],\end{array}  \tag{35}\]  
where \(w_{N}>0\), \(w_{L}>0\), \(w_{T}>0\), \(w_{u}>0\) are the weighting coefficients. Since \(\mathbf{D}\) represents positive quantities of chemotherapy to be administered, the additional constraints are imposed as \(\mathbf{D}>0\). Because the responses \(N\big{(}kT_{s},\mathbf{D},\mathbf{\hat{\theta}}\big{)}\), \(L\big{(}kT_{s},\mathbf{D},\mathbf{\hat{\theta}}\big{)}\), \(T\big{(}kT_{s},\mathbf{D},\mathbf{\hat{\theta}}\big{)}\) are not linear with respect to the decision vector \(\mathbf{D}\), the traditional quadratic programming problem cannot be formulated. Therefore, given that our problem is finite-dimensional, we opt for a numerical optimization while \(N\big{(}kT_{s},\mathbf{D},\mathbf{\hat{\theta}}\big{)}\), \(L\big{(}kT_{s},\mathbf{D},\mathbf{\hat{\theta}}\big{)}\), \(T\big{(}kT_{s},\mathbf{D},\mathbf{\hat{\theta}}\big{)}\), \(u\big{(}kT_{s},\mathbf{D},\mathbf{\hat{\theta}}\big{)}\) are obtained by the numerical solution of differential equations given by (1) using the 4 th-order Runge-Kutta method. To determine the vector \(\mathbf{D}\) that minimizes (35), numerical optimization was performed using the fmincon function available in MATLAB Optimization Toolbox (R2023b, [[https://www.mathworks.cn/products/optimization.html](https://www.mathworks.cn/products/optimization.html)]([https://www.mathworks.cn/products/optimization.html](https://www.mathworks.cn/products/optimization.html))) [60]. This function implements the interior-point method to solve constrained numerical optimization problems [61; 62]. Interior-point methods represent a class of optimization algorithms that focus on exploring the interior of the feasible region. They use a barrier function (typically logarithmic) to penalize points near the boundary, ensuring iterates remain within the feasible region defined by inequality constraints. The barrier parameter is gradually reduced, allowing the solution to approach the boundary as it converges to the optimum. These methods leverage a central path, a trajectory defined by the solutions of the barrier-augmented optimization problem as the barrier parameter is reduced, to approach the optimal solution satisfying the primal and dual conditions (Karush-Kuhn-Tucker conditions).