Stability in a Ross epidemic model with road diffusion

1.   Introduction

Cancer is a complex disease that includes various types, such as lung, prostate, breast, and colorectal cancer. A key feature of cancer cells is their uncontrolled growth, which allows them to break through normal limits and invade nearby organs, causing the disease to spread. Cancer occurs when the body's natural control systems fail. The main aim of cancer treatment is to remove all cancer cells from the body while minimizing damage to healthy cells ^[1]. Conventional treatments like chemotherapy frequently result in drug resistance over time. Cancer cells can evolve and create mechanisms to persist, diminishing the effectiveness of these therapies as treatment continues. Many traditional therapies indiscriminately target rapidly dividing cells, causing harm to both cancerous and healthy cells ^[2]. This can lead to considerable side effects and complications, negatively impacting patients' overall quality of life. Additionally, tumors are often heterogeneous, meaning that various cells within the same tumor may respond differently to treatment. This variability complicates the effectiveness of standard therapies, as they may not effectively target all types of cancer cells present in a tumor ^[3]. Immunotherapies are becoming an essential part of comprehensive strategies for treating certain types of cancer. The goal of immunotherapy is to strengthen the body's natural defenses against cancer by improving the immune system's effectiveness. It is also important to recognize that individuals with weakened immune systems, such as those with AIDS, are at a higher risk of developing specific rare cancers ^[4]. Conventional cancer treatments often come with various side effects, including nerve damage, reduced bone marrow function, and other long-term health problems. These side effects can greatly affect patients' quality of life and may need extra management ^[5]. Although combination therapies can improve treatment effectiveness, they also present challenges. Managing the potential side effects of multiple drugs and preventing increased drug resistance are important concerns that need to be addressed. To address these gaps, it is essential to develop innovative strategies that include personalized medicine, enhanced combination therapies, and a deeper understanding of tumor biology to improve treatment outcomes for cancer patients ^[6].

To address the shortcomings of traditional cancer treatments, it is crucial to create new strategies that incorporate optimal control theory (OCT). OCT facilitates the development of mathematical models that provide insights into tumor growth and how it reacts to various therapies. Moreover, OCT allows for immediate adjustments to treatment plans based on the patient's current responses. This flexibility enhances cancer management by enabling swift modifications to treatments as tumor behavior or patient health changes. Additionally, OCT can assist in designing combination therapies that target several cancer pathways simultaneously ^[7].

Nowadays, the practical derivation of optimal control strategy can be done via different numerical methods. This includes the multiple overshooting transcription method with the interior point optimization (IPOPT) numerical solver, the state-dependent Riccati equation (SDRE), the approximate sequence Riccati equation (ASRE). Each method can help in determining optimal dosages and therefore can effectively combat the disease while preserving the integrity of the patient's immune system.

The objective of the paper is to tailor cancer treatment strategies to meet the unique needs and responses of young and old patients, ensuring the best possible outcomes while considering their varying physiological conditions. This study:

• Focuses on maximizing immune response, minimizing tumor growth, and managing side effects for young patients.

• Emphasizes preserving immune function, controlling tumor progression, and enhancing tolerability for old patients.

The paper is organized as follows: Section 2 will start by reviewing the previous attempts to fight cancer. In Section 3, the studied model is shown, along with its main parameters and terms, including the mathematical preliminaries of the OCT and the necessary conditions for optimality. This section will conclude by discussing each of the proposed techniques, along with their mathematical formulation in solving the OCP of the quadratic regulator type. Section 4 will show the proposed solutions for each case study using all proposed techniques and finish with a comparison of the results. Section 5 will provide a conclusion about the proposed technique and its solution and possible future work enhancements.

2.   Literature review

The utilization of control engineering principles in the field of cancer treatment has gained considerable significance in the literature; numerous control methods have been explored and studied. In ^[7,8,9], authors used a similar model as in ^[10]; applying optimal control treatment of chemotherapy, they used a famous nonlinear state feedback method named SDRE. However, the proposed treatment protocols of ^[7,9] were continuous (i.e., not in the form of doses), hence they cannot be applied in practice. In ^[11], the author also used an optimal control-based therapy for the model of ^[10] (with an added delay parameter in the differential equation to account for the time needed for the immune cells to be stimulated by the tumor cells), in addition to the chemotherapy drug. The results were neither satisfying nor practical. In ^[12], a fractional order version of the model in ^[13] was used, considering only immunotherapy treatment adaptive sliding mode control. The researchers in ^[14] developed a framework called control theory for therapy design (CT4TD). This framework utilizes OCT and applies it to patient-specific models of pharmacokinetics (PK) and pharmacodynamics (PD). Its application to Imatinib administration in chronic myeloid leukemia showed diversified and improved therapeutic strategies among patients compared to standard regimes. The authors in ^[15] presented optimal control in metastatic castrate-resistant prostate cancer. By adopting an optimal control theory approach, the study identified better treatment schedules that can minimize or eradicate resistant cancer cell subpopulations. Control theory, along with machine learning methods like reinforcement learning, are being applied to design cancer treatment schedules ^[16]. These models help in optimizing treatment schedules, modeling treatment responses, and incorporating patient-specific responses for optimum treatment options. The authors of ^[17] detailed how BBR's anticancer properties are regulated by several molecular pathways, including those that cause apoptosis, autophagy, and cell cycle arrest and impede invasion and metastasis. In ^[18], the application of optimal control theory in oncology was discussed, particularly in radiation and systemic therapies for cancer. It emphasized the potential of integrating patient-specific mathematical models and optimal control theory to significantly improve patient outcomes in cancer therapy.

The authors of ^[19] utilized optimal control tools and methods to examine various minimally parameterized models that illustrate the dynamics of cancer cell populations and different aspects of the tumor microenvironment under different anticancer therapies. The authors of ^[20] presented a mathematical model for breast cancer treatment, utilizing a system of ordinary differential equations to evaluate the effects of chemotherapy and a ketogenic diet. A key aspect of the research is the application of optimal control theory, which is used to determine the most effective combination of chemotherapy and a ketogenic diet in reducing cancerous cells, guided by Pontryagin's maximum principle. The study's theoretical results are supported by numerical simulations, underscoring the potential effectiveness of combining chemotherapy with a ketogenic diet in breast cancer treatment. In ^[21], four mathematical models of chemotherapy were analyzed, each based on different sets of ordinary differential equations and including treatments like chemotherapy, immunotherapy, and anti-angiogenic therapy. The study found that optimally controlled therapy can significantly alter tumor outcomes, potentially leading to eradication, unlike with standard treatments or incorrect management. Reference ^[22] introduced optimal control theory and applied it to radiation and systemic therapy, supported by literature examples. It also highlighted the potential of integrating patient-specific mathematical models with optimal control theory to significantly improve patient outcomes. In ^[23], authors described a new approach in cancer treatment that focuses on adaptive drug therapies, which have shown greater efficiency compared to traditional continuous maximum tolerated dose (MTD) methods. This adaptive approach tailors drug dosages based on the evolving state of the tumor rather than following a fixed schedule. The novelty of this method lies in its systematic optimization of adaptive policies using an evolutionary game theory model of cancer dynamics combined with dynamic programming. Specifically, the study optimized two main objectives: reducing total drug usage and shortening the time to recovery. A summary of cancer treatment studies previously discussed is presented in Table 1.

Many suggested treatment protocols involve continuous dosing, which can be hard to implement in real-life situations where discrete dosing is often preferred for better patient adherence and management of side effects. Balancing multiple goals, such as minimizing total drug use and speeding up recovery time, can be difficult and may require compromises between these goals. Some studies have employed open-loop methods, like the shooting method, which can lead to results that are not satisfactory and hard to apply in practice. The proposed treatment protocols may not easily translate to real-world situations due to the complexities of cancer progression and how patients respond to treatments.

The IPOPT can manage complex constraints and objectives, allowing for the simultaneous optimization of various treatment factors, which is essential in cancer therapy. Meanwhile, ASRE and SDRE can be tailored to optimize treatment protocols for discrete dosing instead of continuous administration. This adjustment can enhance patient adherence and help manage side effects more effectively. These methods support multi-objective optimization, enabling the simultaneous pursuit of different goals, such as minimizing drug usage and speeding up recovery times. By incorporating trade-offs into the optimization model, healthcare providers can identify a balanced approach that meets diverse treatment objectives without compromising patient care. Integrating ASRE, SDRE, and IPOPT into cancer treatment protocols can improve their practicality, refine dosing schedules, balance multiple treatment goals, and ultimately enhance patient outcomes by making therapies more responsive to individual needs.

3.   Materials and methods

3.1.   Dynamical model

While the precise understanding of the immune system's involvement in eradicating cancerous tissue remains incomplete, we can advance our comprehension of immune modulation effects by developing models that depict the interaction between tumors and the immune system based on empirical data. Various scientific publications have explored mathematical models in this context (e.g., ^{[10,14,24,25,26]}). Notably, the model proposed by de Pillis ^[10] has garnered significant attention in the literature, being the model utilized here. The model has the following main components:

1) Immune response: The model incorporates immune cells that can experience enhanced growth in response to the tumor's presence. These immune cells have the ability to eliminate tumor cells through a dynamic process. The model assumes only the $\mathrm{C}\mathrm{D}{8}^{+}$ and NK T cells.

2) Competition terms: Within the model, there is competition for available resources between NK cells and tumor cells. Additionally, $\mathrm{C}\mathrm{D}{8}^{+}$ immune cells and tumor cells engage in a predator-prey relationship, competing with each other.

3) Chemotherapy and immunotherapy: The model includes the effect of a famous chemotherapy drug (doxorubicin) on all the modeled cells and the effect of the injected tumor-activated $\mathrm{C}\mathrm{D}{8}^{+}$ T cells on the immune cells.

The model of de Pillis was further enhanced in ^[27,28,29] to include the contribution of IL-2 cytokines and other lymphocyte cells in addition to NK and $\mathrm{C}\mathrm{D}{8}^{+}$ T cells; however, due to the complexity of these models, the model in ^[10] was used instead. This is a first-order, four-dimensional ODE model that includes four states and two control variables, classified as follows:

States:

$E\left(t\right):$ Represent the $\mathrm{C}\mathrm{D}{8}^{+}$ T cells.

$T\left(t\right):$ Represent the tumor cells.

$N\left(t\right):$ Represent the natural killer (NK) cells.

$M\left(t\right):$ Depict the level or amount of the chemotherapy drug present.

Controls:

$w\left(t\right):$ Represent the injected tumor-activated $\mathrm{C}\mathrm{D}{8}^{+}$ T cells.

$v\left(t\right):$ Represent the injected doxorubicin drug.

The model is formulated in ^[10], and Table 2 gives the description of each equation. For more information about the origin of these parameters and how they were measured, please refer to ^[19,28], where the author did an astonishing job in explaining each parameter's origin in detail along with the kind of study in which it was used (i.e., either in vitro or in vivo). Table 3 gives a description of each parameter of the model ^[28,29].

3.2.   Optimal control theory

OCT is a class of modern control theory that seeks to find the best control and state trajectories that either minimize or maximize a certain criterion while satisfying the physical constraints imposed on the system. The solutions provided by OCT can be broadly classified into open-loop or closed-loop solutions. An open-loop solution represents the optimal policy (control actions) that must be followed starting from a specific state initial condition, while the closed-loop solution only gives the optimal solution as a function of the current state of the system ^[7]. Many real-life problems can be formed as a nonlinear programming problem (NLP). In other real-life problems, the decision variables are not static but rather a continuous function of time; these are called optimal control problems. The OCP has the following mathematical formulation ^[8,29]:

x: is the dependent state variable vector in ${\mathcal{R}}^{N}\times [{t}_{0}, {t}_{f}]$ ;

u: is the independent control variable vector in ${\mathcal{R}}^{M}\times [{t}_{0}, {t}_{f}]$ ;

J: ${\mathcal{R}}^{L}\times [{t}_{0}, {t}_{f}]\to \mathcal{ }\mathcal{R}$ is the optimization criterion;

L: ${\mathcal{R}}^{L}\times \left[{t}_{0}, {t}_{f}\right]\to \left[{t}_{0}, {t}_{f}\right]$ is the bath cost function;

$\varphi :\mathcal{ }\mathcal{R}\times \mathcal{R}\to \mathcal{R}$ is the final state cost function;

${g}_{1}:{\mathcal{R}}^{N}\to {\mathcal{R}}^{{M}_{1}}$ , are the state inquality constraints;

${g}_{2}:{\mathcal{R}}^{M}\to {\mathcal{R}}^{{N}_{1}}$ , are the control inquality constraints;

a: ${\mathcal{R}}^{N}\to {\mathcal{R}}^{N}$ , are the dynamical constraints of the system.

To derive the necessary conditions for optimality (NCO) ^[30,31], we will initially assume that there are no control inequality constraints. Subsequently, we can incorporate these constraints using Pontryagin's minimum principle ^[32]. A new function, known as the Hamiltonian, which defines the NCO, is presented as follows:

In the quadratic case, the function $L$ can be expressed as a sum of separate quadratic terms concerning both $x$ and $u$ . This means that the cost associated with the control inputs and state variables can be disentangled, allowing for a clear distinction in how each variable affects the overall cost.

Where $\boldsymbol{\lambda }\left(t\right)$ is a vector variable of the same size as $\boldsymbol{x}\left(t\right)$ (i.e., ${\mathcal{R}}^{N}\times [{t}_{0}, {t}_{f}]$ ), named as the costate variable or Lagrange multiplier. Their main goal is to convert the dynamic constrained OCP into an unconstrained OCP at the expense of doubling the number of state variables. For more information on the derivation of the Hamiltonian equation of (2) and the NCO, please refer to ^[30,32].

From (2), the NCO can be written as:

With the associated boundary conditions:

Where: $\varphi$ is the final state function as defined in (1) and $\delta {\boldsymbol{x}}_{f}$ , $\delta {t}_{f}$ are the differentials of the terminating state $\boldsymbol{x}\left({t}_{f}\right)$ and final time ${t}_{f}$ , respectively. If the final state $\boldsymbol{x}\left({t}_{f}\right)$ is not specified and ${t}_{f}$ is specified, which is the case in many OCP problems, then the boundary condition can be written as:

Pontryagin showed that the optimal control necessary condition can be replaced with the following more general condition ^[30]:

Equation (6) is called the Pontryagin minimum principle (PMM) and states that the Hamiltonian of the ${\boldsymbol{u}}^{\mathrm{*}}\left(t\right)$ must be less than or equal to the Hamiltonian of any other admissible control $\boldsymbol{u}\left(t\right)$ at all times ^[30,31].

Hereafter, the OCP is assumed to be of the quadratic regulation problem (QRP) form, where the objective is to drive the states $\boldsymbol{x}$ from any initial condition $\boldsymbol{x}\left(0\right) = {\boldsymbol{x}}_{0}$ to the equilibrium state $x\left({t}_{f}\right) = 0$ as ${t}_{f}\to \mathrm{\infty }$ .

It should be noted that if the desired equilibrium point of the system was not $\boldsymbol{x} = 0$ but rather $\boldsymbol{x} = {\boldsymbol{x}}_{\boldsymbol{f}}$ , then we can transfer this equilibrium point to the origin by change of variables ^[33]. The OCP problem of (1) is rewritten in the new QRP form with some changes in the constraints.

The NCO for the OCP of (7) can be written using the conditions in (3) and the PMM equating of (6) as:

Where the boundary conditions can be written from (5) as:

To solve OCPs numerically, they must first be formulated in discrete form. The OCP of (8) can be discretized by using a proper sampling time $h$ as:

Where $\boldsymbol{x}\left(t\right)$ is replaced with $\boldsymbol{x}\left(k\right)$ for the $tϵ[k\cdot h, \left(k+1\right)\cdot h]$ , the same is true for $\boldsymbol{Q}\left(k\right)$ , $\boldsymbol{R}\left(k\right)$ , $\boldsymbol{u}\left(k\right)$ ; however, the dynamic constraints are discretized using any of the famous methods of forward Euler, e.g., trapezoidal or ${4}^{th}$ order Runge-Kutta. To write the necessary conditions for optimality, we first need to write the Hamiltonian:

Then, applying the necessary conditions of (3) on Eq (11) results in

3.2.1.   Direct multiple shooting transcription method

This technique converts a continuous-time OCP into a discrete-time nonlinear programming (NLP) problem, a process known as transcription. The transcription described in Eq (10) is called direct single shooting because it only optimizes the controls $u\left(0\right), u\left(1\right), \dots u\left(k\right)$ . However, this method does not ensure that the dynamic constraints are met (i.e., $x(k+1)\ne f\left(x\left(k\right), u\left(k\right)\right)$ ). An alternative approach, called direct multiple shooting, treats the state vector $x\left(k\right)$ as a decision variable in the optimization problem. Since this problem will be solved using a computer, ${t}_{f}$ should be set to a large finite number instead of infinity. Additionally, state inequality constraints can be included directly in the multiple shooting transcription. Note that $\boldsymbol{u}\left(i\right)$ and ${\boldsymbol{u}}_{i}$ will be used interchangeably. Below is the transcription of the OCP of (6) using the direct multiple shooting method:

Note that $K$ represents the number of samples for the whole-time interval, which can be related to the continuous time interval using the relation $K = ⌈\frac{{t}_{f}-{t}_{0}}{h}⌉$ .

Where $h$ is the sampling period of the discretization process. Before getting into the details of the benchmark solver, a brief introduction to the famous mathematical framework computer algebra systems for algorithmic differentiation (CasADi) is presented. CasADi was mainly developed to solve OCPs, by first transforming them into an equivalent NLP problem using any transcription method ^[34]. Although CasADi is not used to solve NLP problems, it can be equipped with solvers that can tackle the NLP problem, which was originally formulated using CasADi without additional specific configuration for that solver. For additional information on CasADi, please refer to the original documentation ^[35]. One famous NLP solver built-in with CasADi solvers is the IPOPT ^[36].

The IPOPT method begins by transforming the constrained optimization problem into an unconstrained one by utilizing a barrier function. The complete optimization problem incorporating a barrier function can be expressed as follows:

The new task now is to find the minimum point $\boldsymbol{x}$ that minimizes (14), where $\mu$ is a manipulated parameter that determines the effect of the constraint $g\left(\boldsymbol{x}\right)$ on the optimization criterion $f\left(\boldsymbol{x}\right)$ . Ideally, $\mu$ should approach zero; however, the IPOPT problem solves the NLP problem of (14) iteratively. The algorithm finds the approximate optimal solution of each iteration by solving the corresponding first-order necessary conditions of the NLP problem. These necessary conditions are often called Karush-Kuhn-Tucker (KKT) conditions. For more details on this algorithm see ^[35].

3.2.2.   State-dependent coefficient method

The state-dependent coefficient (SDC) method utilizes linear quadratic regulator (LQR) theory to develop optimal control strategies for nonlinear systems. The integration of LQR theory with the state-dependent coefficient method provides a robust framework for optimizing control in nonlinear systems, enhancing performance and adaptability as the following farmwork:

1) In the SDC approach, the system dynamics are represented with coefficients that change according to the current state, allowing for a more precise depiction of complex, nonlinear behaviors.

2) The SDC method parallels the traditional LQR formulation by deriving a state-dependent Riccati equation. This equation includes state-dependent coefficients, facilitating the application of LQR techniques in nonlinear scenarios.

3) Solving the SDRE enables the derivation of an optimal feedback control policy that minimizes a quadratic cost function, leading to control laws that adjust dynamically in response to changes in the system state.

4) The solution process frequently employs iterative methods to improve the approximations of the Riccati equation, ensuring both computational efficiency and accuracy.

In the context of Eq (7), the system dynamics representation is as follows:

Applying the SDC parameterization yields the following system in which both the system and input matrices are explicit functions of current state variables (where the time dependences were dropped for convenience):

The dynamic system representation of Eq (16) is often called a nonlinear state-dependent affine control system. The SDC representation of Eq (16) is similar to the following linear system representation in state space ^[30]:

After defining some terminologies, and before we start with the algorithm, the LQR for LTI systems with its solution is presented shortly. The LQR solution relays on the fact that the optimal value of the constate variable vector $\boldsymbol{\lambda }$ is represented as:

Using the necessary conditions for optimality and Eq (18), it is possible to write:

Where Eq (19) is famously known as the algebraic Riccati equation (ARE), where the only unknown of this equation is the $\boldsymbol{P}$ matrix, and by solving this equation we can write the optimal control ${\boldsymbol{u}}^{\mathrm{*}}\left(t\right)$ as ^[30]:

The SDRE ^[37] technique fundamentally relies on the concept of extended linearization. Extended linearization is the action of generalizing the linear systems–based technique to be applied to the class of nonlinear systems. The SDRE technique redefines the system of (7) as:

Then, by the extended linearization concept, we can write the following set of equations, which were mimicked from the linear case. The costate-state matrix relation:

and

If the matrix $P\left({x}^{\mathrm{*}}\right)$ was obtained from Eq (23), then the control ${u}^{\mathrm{*}}\left(t\right)$ can be written in terms of the states ${x}^{\mathrm{*}}\left(t\right)$ as:

Note that since the derivation of the SDRE method rises from the extended linearization concept, the optimal control of (24) is not the optimal solution for the OCP of (21) but rather a suboptimal solution. However, the SDRE technique's suboptimal solution was proven to be very near to the optimal solution in many nonlinear systems ^[8,38]. In order for the SDRE method to be successfully applied, the ($A\left(x\right), B\left(x\right))$ must be completely controllable in the state admissible region. Note that it was assumed that the $\boldsymbol{P}\left({\boldsymbol{x}}^{\mathrm{*}}\right)$ can be solved in Eq (23). However, this is generally not always possible, since it requires very complex symbolic solvers. In practice and with the aid of a computer, the SDRE equation of (23) is often solved numerically using very fast and accurate algorithms such as Schur and Kleinman-Newton ^[38]. The numerical procedure of the SDRE starts by discretizing the system of (21) as:

(Note that the forward Euler technique was used.) Note that at each time step $k$ , the system in (8) represents a linear time invariant system, for example at $k = 0:$

The SDRE of Eq (26) is now turned into an ARE, which can be efficiently solved for $P\left({x}_{0}\right)$ , which is then used to compute the optimal control ${\boldsymbol{u}}^{\mathrm{*}}\left(0\right)$ as:

Equation (27) is then substituted into Eq (26) to obtain $\boldsymbol{x}\left(1\right)$ , and the algorithm continues, where at each time step $k$ , an ARE is solved, and the resultant $P\left(x\left(k\right)\right)$ matrix is used to compute the ${u}^{\mathrm{*}}\left(k\right)$ using (27). Hence, the SDRE algorithm can be represented with the following two equations:

The discrete SDRE equation:

The optimal control equation:

Since the control signal ${\boldsymbol{u}}^{\mathrm{*}}\left(\boldsymbol{k}\right)$ is computed independently at each time step $k$ , this allows the possibility to account for the control inequality constraints of:

The control signal ${\boldsymbol{u}}^{\mathrm{*}}\left(k\right)$ is first computed using Eq (29), then we apply the PMM principle of (7) as follows:

If the optimal control at time step $k$ was computed to be equal to ${\boldsymbol{u}}_{\boldsymbol{M}}$ (where ${\boldsymbol{u}}_{\boldsymbol{M}} > {\boldsymbol{u}}_{\boldsymbol{u}}$ in at least one component of $\boldsymbol{u}$ ), then by the PMM ^[30] principle, the Hamiltonian cannot assume a minimum value than the one achieved by setting ${\boldsymbol{u}}^{\mathrm{*}}\left(k\right) = {\boldsymbol{u}}_{\boldsymbol{u}}$ . The same idea follows for the case of ${\boldsymbol{u}}^{\mathrm{*}}\left(k\right) = {\boldsymbol{u}}_{\boldsymbol{m}} < {\boldsymbol{u}}_{\boldsymbol{l}}$ , then ${\boldsymbol{u}}^{\mathrm{*}}\left(k\right)$ is set to ${\boldsymbol{u}}_{\boldsymbol{l}}$ , we can formally write the refined optimal control equation of (29) in the following form:

Equation (31) simply means that the control single is first computed using Eq (29) and then fed to a nonlinear saturation system of lower and upper limits of ${\boldsymbol{u}}_{\boldsymbol{l}}, {\boldsymbol{u}}_{\boldsymbol{u}}$ .

3.2.3.   Approximate sequence Riccati equation

Although the SDRE offers a suboptimal robust and stable solution for the OCP, the SDC parametrization of (17) assumes that the system can be represented in the following nonlinear state-dependent affine form: $\dot{x} = A\left(x\right)\cdot x+B\left(x\right)\cdot u$ .

Nevertheless, the SDRE technique is not applicable to general nonlinear nonaffine optimal control problems (OCPs) with nonquadratic performance indices. The proposed algorithm utilizes the globally converged solution obtained from a series of approximating Riccati equations (ASRE) to explicitly design time-varying feedback controllers for the subsequent general nonlinear nonaffine state equation, as described in ^[37].

Before proceeding with the ASRE algorithm, the LQR theory for linear time-variant (LTV) systems is presented. The performance criterion and the system dynamics are represented as follows:

The SBVP, which follows from the necessary conditions of state and costate equations, can be written as:

The matrix $\boldsymbol{S}\left(t\right)$ is defined as $S\left(t\right) = B\left(t\right)\cdot {R}^{-1}\left(t\right)\cdot {B}^{T}\left(t\right)$ .

With the same boundary conditions ${\boldsymbol{x}}^{\mathrm{*}}\left(0\right) = {\boldsymbol{x}}_{0}$ of and ${\boldsymbol{\lambda }}^{\mathrm{*}}\left(\mathrm{\infty }\right) = 0$ . Following the same assumption of state-costate dependency, i.e., ${\boldsymbol{\lambda }}^{\mathrm{*}}\left(t\right) = \boldsymbol{P}\left(t\right)\cdot {\boldsymbol{x}}^{\mathrm{*}}\left(t\right)$ , the time-dependent Riccati equation is formed:

Solving for the time varying matrix $\boldsymbol{P}\left(t\right)$ results in the optimal control equation:

Although the ASRE can be applied to general nonlinear nonaffine systems with nonquadratic cost functions, the methodology derived here applies only to OCP with quadratic cost functions. The nonlinear nonquadratic optimization problem, which involves minimizing the quadratic cost of Eq (8) while adhering to the dynamics specified in Eq (32), can be converted into an equivalent linear-quadratic, time-varying problem by introducing the subsequent LTV sequence as described in Eq (38). It is important to note that the time dependency (t) has been omitted for the sake of convenience.

And the corresponding performance criterion at each iteration $i$ is:

Where the superscript $\left[i\right]$ represent the current iteration (starting from 0). The first approximation (i.e., the $i = 0$ case in (38)) does not truly represent an LTV system, since the matrices A, B are constant over the entire time $t$ and are only a function of the initial state ${\boldsymbol{x}}_{0}$ . However, for all the consecutive iterations (i.e., the $i\ge 1$ case in (38)), the system is represented by an LTV system whose matrices A, B are evaluated from the previous iteration, i.e., $i-1$ . The feedback optimal policy (control) for each iteration $i$ for the corresponding system dynamics in (38) and the performance criterion (39) takes the form:

Where ${\boldsymbol{P}}^{\left[i\right]}$ is the $N\times N$ symmetric matrix and is calculated continuously (if possible, using advanced symbolic solver tool) from the following Riccati equation:

Or numerically using appropriate sampling interval $h$ at each time step $k$ using the following ARE, which is often used in practice:

The author in ^[38] does an excellent job of describing the global convergence of the ASRE technique, which states that as the iteration index $i$ increases, the system converges to the optimal solution, where the convergence criterion is $\underset{i\to \infty }{\mathrm{lim}}‖{x}^{\left[i\right]}\left(t\right)-{x}^{[i-1]}\left(t\right)‖ = 0$ .

3.3.   Model manipulation

Before solving the OCP of the cancer therapy system, we first need to manipulate the model, so that the OCP can be formulated as a QRP and thus the derived techniques can be applied. The first step is to rewrite the model using the states $x\left(t\right)$ and controls $u\left(t\right)$ notation; the resultant reformulation is shown as ^[25]:

Note that the mapping between the original states and the new states as well as the control signals are as follows:

The next step is to find the equilibrium points of the system by setting all the dynamic equations to zero and solving for the value of the states (note that at equilibrium points, the control signals are also set to zero). Setting control signals to zero when determining equilibrium points in a cancer model enables researchers to examine the system's inherent dynamics free from external factors. This process is essential for recognizing stable states and gaining insights into how various elements affect tumor growth or reduction in a controlled setting. The solutions can be classified into the following categories ^[25]:

Tumor-free: The equilibrium point has the solution form: $\left(\frac{s}{{d}_{1}}, 0, \frac{1}{{b}_{2}}, 0\right).$

Dead: If at any equilibrium point, the normal cell population is zero, the point is considered dead. The solution form for the equilibrium point is $\left(\frac{s}{{d}_{1}}, \mathrm{0, 0}, 0\right)$ .

Next, in order to formulate the cancer therapy system as a QRP, we first need to shift the equilibrium point to the origin, and this can be easily done by defining the following new set of variables:

Where ${\boldsymbol{x}}_{\boldsymbol{f}, \boldsymbol{i}}$ is the ${i}^{th}$ component of ${\boldsymbol{x}}_{f}$ , which is given by:

Note that ${\theta }_{i}$ is just a dummy variable, and for convenience, it will be replaced by ${x}_{i}$ for $i = 1\dots 4$ . Hence, the resultant shifted system dynamic equations are:

We notice that it can be written in the nonlinear affine system format, where the matrices $A\left(x\right), B\left(x\right)$ are formulated as follows:

Where:

In order to check whether the choice of $\boldsymbol{A}\left(\boldsymbol{x}\right), \boldsymbol{B}\left(\boldsymbol{x}\right)$ pair is controllable or not, we need to compute the controllability matrix ${\boldsymbol{M}}_{\boldsymbol{c}}\left(\boldsymbol{x}\right)$ . However, computing this controllability matrix is computationally difficult. Nevertheless, it can be shown that if the states system plus the tumor-free equilibrium point is strictly positive at any time instance, then by using the controls ${u}_{1}, {u}_{2}$ , the system can always be brought to the tumor-free state ^[16,21].

3.4.   Problem formulation

We can proceed by defining the system optimization criterion and the inequality constraints. Since the cancer therapy system will be formulated as a QRP, the optimization criterion will take the following form:

Note that for implementation purposes, the upper limit of the integration is not set to infinity but rather a finite large number ${t}_{f}$ .

The state inequality constraints are of two types:

1) Non-negativity.

2) Normal cells ${x}_{3}\left(t\right)$ stay above a lower bound of ${N}_{m}$ of the total carrying capacity ^[10].

The mathematical representation of these constraints can be expressed as follows:

The control inequality constraints represent the upper and lower limits of the drug doses. The upper and lower limits are given by $1$ and $0$ , respectively. As in ^[8], these constraints can be written in mathematical terms as $0\le u\left(t\right)\le 1, tϵ[0, \infty]$ .

The $\boldsymbol{R}\left(t\right)$ matrix is simply chosen to be a constant diagonal matrix of the form: $R\left(t\right) = \left[\begin{array}{cc}{R}_{1} & 0\\ 0 & {R}_{2}\end{array}\right]$ .

The $\boldsymbol{Q}\left(t\right)$ matrix is also chosen to be a constant diagonal matrix; the resultant $\boldsymbol{Q}\left(t\right)$ matrix is:

Where ${Q}_{N}\left(t\right)$ is given by:

Since the selected Q(t) cannot be expressed as a continuous mathematical function, the QRP must be addressed in discrete time. Various numerical integration methods can be employed to convert the continuous-time OCP, including techniques like Runge-Kutta, Trapezoidal, and Simpson's rule. However, for simplicity, the commonly used forward Euler approximation was applied with an appropriate sampling time h. Therefore, the discrete QRP for the cancer chemo-immune therapy system can be formulated as follows ^[39,40]:

The QRP of the above equation is applied to two different case studies. To each case, two different therapy rules are applied. These case studies are restated again for convenience:

Case 1: In the scenario where a young patient is diagnosed with cancer but does not have any other resistant diseases, the oncologist's objective during chemotherapy is to primarily focus on decreasing the population of cancerous cells. This is because the patient's young age allows their body to compensate for the reduction in the number of normal cells and immune cells that may occur during chemotherapy. From an optimal control standpoint, minimizing the number of tumor cells is prioritized rather than maximizing the number of immune cells. This objective can be accomplished by appropriately adjusting the weighting matrices in the optimal control problem.

Case 2: In the instance where an elderly patient is diagnosed with cancer along with other refractory conditions like cardiac disease, it becomes risky to destroy normal cells (including immune cells). Consequently, the oncologist places greater emphasis on immunotherapy treatment rather than chemotherapy. The preservation of normal cells (or immune cells) takes precedence over reducing the number of cancerous cells in this situation. Similar to the previous case, the optimal control problem incorporates this prioritization by selecting appropriate weights in the weighting matrices.

In the case of younger patients, parameter weights are generally modified to focus on reducing tumor growth, as their bodies can better handle decreases in normal and immune cell counts. This results in assigning higher weights to the minimization of tumor cell populations and lower weights to the preservation of normal cells. Conversely, for older patients, the priority shifts to safeguarding healthy cells due to their greater susceptibility and possible comorbidities. Consequently, maintaining normal cell populations receives higher weights, which may result in a less aggressive strategy for tumor reduction. For patients who do not clearly fall into the young or elderly categories, parameter weights need to be tailored to their specific situations. Important factors to consider include:

• Age: Middle-aged patients may exhibit different tolerances and responses than younger or older individuals.

• Health status: Existing health conditions and overall well-being can affect the aggressiveness of the treatment approach.

• Tumor characteristics: The severity of the tumor and its responsiveness to treatments should also be taken into account.

• Patient preferences: Integrating patient values and preferences regarding treatment objectives can facilitate more personalized care.

4.   Results and discussion

We are prepared to implement the optimal control techniques outlined in the previous section to address the QRP. Before we begin, Table 4 lists the precise numerical values for the parameters of the cancer therapy OCP, while Table 5 details the parameters of the cancer model. These parameters remain constant across all examined techniques to facilitate a performance comparison with ^[41]. Each technique was applied for each case classified as continuous and dosed treatment (either Q21D or Q16D). Cancer treatment is considered continuous when it consists of ongoing therapies designed to manage the disease over a prolonged period, rather than providing a one-time or short-term intervention. This method is frequently used in situations where cancer cannot be entirely cured but can be effectively controlled or managed. In contrast, cancer treatment is categorized as discrete when therapies are administered at specific intervals, such as Q21D or Q16D, instead of continuously.

Each time, all techniques (IPOPT, SDRE, and ASRE) are applied and compared with each other and with other methods [LQR and single network adaptive critic (SNAC)] studied in ^[41], which used the same model ^[10]. Figures 1–4 contain the four states of the system (i.e., $\mathrm{C}\mathrm{D}{8}^{+}$ and NK T cells, tumor cells, and chemotherapy drug) and the two control therapies (i.e., immunotherapy and chemotherapy) and the optimization criterion of OCP for each technique.

• IPOPT: The quadratic regulation problem (QRP) was first formulated using CasADi framework; then, a problem instance was created and then solved using the IPOPT solver.

• SDRE: The solution of the ARE was computed using a very fast and efficient MATLAB built-in function called icare, which is an implicit solver for continuous-time algebraic Riccati equations.

• ASRE: Similar to the SDRE technique, the ASRE technique uses the icare MATLAB ^{[42,43,44,45,46]} built-in function to solve the ARE at each time step $k$ in each iteration $i$ of the algorithm.

Before starting to interpret and compare the results of the cancer therapy OCP, the purpose for each kind of therapy (i.e., continuous and dosed) will be elaborated. Continuous therapy is not a real type of therapy and cannot be applied in practice for obvious reasons; however, it is included in this paper for mathematical purposes.

Continuous therapy represents the actual solution for the OCP of the discrete QRP of the cancer chemo-immune therapy system and it is the one that can be used to compare the performance of the studied techniques. All optimal control techniques generated nearly the same control sequences for both 𝑢₁ and 𝑢₂. The IPOPT technique generated the best response (in terms of the value of the optimization criterion); the IPOPT solver used very sophisticated heuristic techniques in computing the optimal solution. The SDRE and the ASRE techniques came in second place as they did not find the optimal solution of the OCP due to the extended linearization theory. Looking at the same figure again, we can hardly tell the difference between the optimal control techniques studied in the optimization criterion section. This suggests that the results of these techniques are very similar.

Comparing the figures of continuous therapy solutions for all techniques for cases 1 and 2, we can notice an interesting result regarding treatment policies. Although the control weighting matrix R was identical for both cases, the response of the optimal control techniques showed a noticeable difference, where in the simulation results all these techniques agreed to use less of the chemotherapy drug and rely more heavily on the use of immunotherapy. The reason behind this difference is how the state weighting matrix Q was formulated in both cases: in case 1, there was a large weighing on the tumor population with respect to the immune cells (i.e., the $\mathrm{C}\mathrm{D}{8}^{+}$ and NK T cells) and hence the optimal control techniques produced a response that will eradicate the tumor as fast as possible. However, in case 2, the scenario changed a bit, where the static weight of the tumor cells was lowered as compared to the weight used in case 1, which can be interpreted as "we still want to eradicate the tumor cells but not as much, and we also want to keep the immune cells as high as possible". Hence, we can observe how these weighting matrices could affect the response significantly; in general, we can use time-varying matrices that focus on eradicating the tumor cells when they have a relatively large population and on increasing the immune cells when their population is relatively low. Table 6 summarizes the final values of the state variables and optimization criterion of all studied techniques for both the continuous (after 20 days) and dosed treatment (after 90 days) for case 1, while Table 7 does the same for case 2.

The SDRE and ASRE methods are considered closed-loop solutions in cancer treatment because they continually adjust treatment plans based on real-time patient monitoring. This flexibility helps ensure that therapies remain effective as the patient's condition evolves, leading to better disease management. SDRE modifies treatment inputs according to the current state of the patient and tumor dynamics. This allows for immediate changes to the treatment plan as new information is received, using feedback from how the patient responds to therapy. It dynamically optimizes dosages of chemotherapy and immunotherapy, keeping treatments effective while minimizing side effects. ASRE, on the other hand, is specifically designed for nonlinear systems, which is important in cancer treatment due to the unpredictable nature of tumor behavior. It uses a feedback mechanism to approximate optimal control based on observed patient responses. Combining open-loop and closed-loop methods for cancer treatment presents an innovative approach to enhancing therapeutic effectiveness. Open-loop methods deliver a predetermined dose without real-time feedback, which can lead to suboptimal drug levels and potential toxicity. In contrast, closed-loop systems continuously monitor patient responses and adjust treatment dynamically, improving the precision of drug delivery. Integrating these two approaches could lead to a more personalized treatment regimen. Open-loop methods could establish a baseline treatment plan based on initial assessments, while closed-loop systems could fine-tune therapy in response to ongoing patient feedback and changing conditions. This synergy could enhance the overall effectiveness of cancer therapies by ensuring that patients receive the right dose at the right time, ultimately improving outcomes and reducing adverse effects. By leveraging the strengths of both open-loop and closed-loop strategies, healthcare providers may be able to create more effective and safer cancer treatment protocols that adapt to individual patient needs throughout the course of therapy.

In contrast, closed-loop systems continuously monitor patient responses and make dynamic adjustments, enhancing the precision of drug delivery. Integrating these two approaches could lead to more personalized treatment plans. Open-loop methods can establish a baseline treatment based on initial evaluations, while closed-loop systems can refine therapy in response to ongoing patient feedback and changing conditions. This combination could improve the overall effectiveness of cancer treatments by ensuring patients receive the right dose at the right time, ultimately enhancing outcomes and reducing side effects. By utilizing both open-loop and closed-loop strategies, healthcare providers can develop safer and more effective cancer treatment protocols tailored to individual patient needs throughout their therapy.

Each method has its unique advantages and challenges in the context of cancer treatment.

• IPOPT is suitable for complex, large-scale problems but may struggle with real-time applications.

• SDRE offers adaptability and real-time optimization but requires precise modeling and can be computationally intense.

• ASRE provides a balance between efficiency and control but might sacrifice some accuracy due to its approximations.

The choice of method will depend on the specific requirements of the treatment strategy, including the need for real-time adjustments, computational resources, and the complexity of the tumor dynamics. Table 8 presents the comparison between IPOPT, SDRA, and ASRE based on different criteria.

These techniques can also be compared by computational requirements such as speed and computer memory, since in real-time application, the available memory and maximum sampling time will determine which technique should be used. Unfortunately, we cannot measure the exact memory required by MATLAB code, since MATLAB internally releases memory during code execution. Hence, only the additional required material used to execute the MATLAB code is mentioned (this includes any prewritten code that is crucial for the technique to be executed). All techniques were implemented and executed on MATLAB 2023a software run on a PC Intel(R) Core (TM) i7-8550U CPU @ 1.80GHz with 16GBs RAM and Windows 11 (64-bit) operating system. Table 9 compares the computer computational requirements for the studied optimal control techniques based on continuous therapy of case 1.

For future work, while the Pontryagin method offers a strong basis for optimal control in cancer therapies, investigating the Hamilton-Jacobi-Bellman (HJB) equation could lead to enhanced dosing strategies and improved patient outcomes. Utilizing the HJB equation might enable the development of an optimal control law that considers the entire range of patient responses instead of just immediate feedback, thereby increasing the accuracy of dosage modifications. The HJB approach is particularly useful for addressing nonlinear systems, which are prevalent in cancer treatment. It may produce different outcomes compared to the Pontryagin method, especially in complex situations with unpredictable tumor behavior. Further exploration of this approach would be valuable for a thorough understanding of optimal control in cancer therapy.

5.   Conclusions

We delved into the application of three distinct optimal control problem (OCP) techniques: IPOPT, SDRE, and ASRE. These originate from diverse mathematical areas, such as the calculus of variations and dynamic programming. Utilizing a QRP formulation, we integrated chemotherapy and immunotherapy for two cancer patients, aiming to derive the optimal therapy protocols, both continuous and dosed, for each case. A key revelation of our study is the striking similarity in results yielded by the different techniques, suggesting that each approach approximates an optimal solution closely, often classified as suboptimal. Despite the comparable outcomes, certain distinctions were noted in their applicability. For instance, open-loop methods like IPOPT are not suitable for real-time applications. Conversely, each closed-loop technique has its unique domain of efficacy.

In terms of resource utilization, the proposed method recorded a value of 580.9231 for IPOPT and 583.6575 for both SDRE and ASRE. NK T cell levels remained stable at 0.999 for C in all techniques, while D values reached 1 consistently across methods. The chemotherapy drug administration showed no active dosing in all cases, with both C and D values at 0.

Each method has its unique advantages and challenges in the context of cancer treatment.

• IPOPT is suitable for complex, large-scale problems but may struggle with real-time applications.

• SDRE offers adaptability and real-time optimization but requires precise modeling and can be computationally intense.

• ASRE provides a balance between efficiency and control but might sacrifice some accuracy due to its approximations.

It is important to utilize advanced optimization techniques to tackle the challenges and limitations of applying optimal control principles in cancer treatment. Additionally, thorough validation studies using clinical data are crucial to evaluating the effectiveness and applicability of these methods across various patient groups and types of cancer.

Author contributions

Ahmed J. Abougarair: Original draft preparation, Formal analysis, Investigation, Software; Mohsen Bakouri: Manuscript writing and review, Validation of results, Investigation; Abdulrahman Alduraywish: Conceptualization, Methodology; Omar G. Mrehel: Investigation, Formal analysis; Abdulrahman Alqahtani: Writing and review, Methodology; Tariq Alqahtani: Formal analysis, Writing and review; Yousef Alharbi: Writing and review, Visualization; Md Samsuzzaman: Formal analysis, Validation of results. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the deputyship for Research and Innovation, Ministry of Education, Saudi Arabia, grant number IFP-2022-21.

Acknowledgments

The authors extend their appreciation to the deputyship for Research and Innovation, Ministry of Education in Saudi Arabia for funding this research work through the project number IFP-2022-21.

Conflict of interest

The authors declare no conflicts of interest.