A class of constrained optimal control problems arising in an immunotherapy cancer remission process

1.   Introduction

Cancer is a global health issue, and treating a tumor poses numerous challenges. The selection of an appropriate treatment scheme, such as surgery and immunotherapy, generally depends on the specific characteristics of the tumor cells, including their aggressiveness and the location of their growth. Immunotherapies are employed to either enhance the body's natural immune defenses or restore the immune system functions that have been compromised. Therefore, it is significant to explore immunotherapeutic approaches to treat tumors.

In tumor immunotherapies, strategies such as an Adoptive Cellular Immunotherapy (ACI), cytokines, and other antigen-specific and non-specific agents are expolited to stimulate the immune system, eliminate abnormally proliferating cells, and enhance the immune system's ability to defend itself. ACI involves the injection of in vitro-generated immune cells that are sensitive to the tumor antigens into the host's tumorigenic environment. This process strengthens and expands the host's immune response ^[1]. Interleukin-2 (IL-2) is the primary cytokine responsible for the activation, proliferation, and differentiation of lymphocytes, and is predominantly produced by CD4$^{+}$T cells. At various stages of disease, IL-2 has been shown to enhance the activity of cytotoxic T lymphocytes (CTLs) ^[2]. ACI has been shown to reduce the tumor activity in the host and can be combined with high doses of IL-2 to halt the tumor growth until the tumor cells are eliminated ^[3].

To explore the interactions among tumor cells, immune cells, and their microenvironment, numerous researchers have conducted mathematical modeling of tumor-immune systems; see, for example, ^[4,5,6,7]. Nevertheless, the aim of this paper is to investigate tumor-immune cell interactions by applying the control parameterization strategy. Moreover, numerous researchers have investigated the applications of the optimal control theory ^[8,9,10,11] to cancer dynamics. Murray ^[12] studied a mathematical model of tumor and normal cells, thereby controlling the rate of treatment administration based on logistic and Gompertzian growth patterns. This approach aimed to minimize the tumor burden at the end of the treatment period, while the normal cell population was maintained above a certain threshold to avoid toxicity. Burden et al. ^[13] applied the optimal control theory with the goal of maximizing the effector cell and IL-2 concentrations while minimizing the number of tumor cells. Fister and Donnelly ^[14] explored the dynamics between tumor cells, immune effector cells, and IL-2 using a mathematical model. They identified the conditions under which the tumor could be eliminated and obtained the optimal control strategy by solving the optimality system. DePillis et al. ^[15] developed a mathematical model to analyze the interaction between tumor immunity and chemotherapy, thereby outlining strategies for an optimal treatment implementation. The dynamics of the model were analyzed, the optimal control strategies were established, and both quadratic and linear cases were solved. The optimal control methods related to the drug treatments were identified, and the numerical results of the optimal strategy were discussed. Based on a tumor palliation model, Chakrabarty and Banerjee ^[16] applied the classical control theory to determine how optimal external treatments with ACI and IL-2 could more effectively palliate malignant tumors, while minimizing adverse effects on the immune response. Their findings demonstrated that the combination therapy was the most effective in reducing the tumor load. Although these treatments led to a reduction in the tumor cell levels, ACI may be more beneficial than IL-2 if only one treatment option is utilized. Khajanchi and Ghosh ^[17] developed a system of nonlinear differential equations to model the interactions between cancer cells and the immune system. They formulated an optimal control problem aimed at minimizing the total tumor load and the adverse effects of drugs, while maximizing the total number of effector cells. By utilizing the deterministic optimal control theory, they evaluated the effectiveness of the optimal treatment strategies for both ACI and IL-2 therapies, thereby exploring the synergistic potential of these treatments in cancer dynamics. However, the optimal control of tumor models with constraints has not yet been considered. In particular, analyses and numerical evaluations haven't fully explored the various constraints associated with drug therapy dosages in cancer dynamics. For the optimization control methods of tumor models, most scholars provide time-varying input curves when studying, which consider the time-varying input of continuous nonlinear models. The optimal control solution obtained by this method is the ideal solution under ideal conditions, providing strong data support for practical treatment. However, due to the limitations of practical conditions, the actual input forms are all discrete, which means that in the actual treatment process, the operability is weak and it is impossible to make timely adjustments to changes in the input.

A set of valid comparison criteria was established to be used for the comparison of local collocation and control parameterization in ^[18,19]. Then, the numerical mechanisms behind both methods were examined and analyzed according to these criteria. The purpose of ^[20] was to extend the control parameterization methods to solve problems with both control-dependent and discrete delay arguments. Recently, Teo et al. ^[21] summarized the latest advancements in the control parameterization technique. They outlined methods to address the optimal control problems, both with and without constraints, and highlighted that these problems could be effectively solved by using MISER3, which is an optimal control software that utilizes control parameterization techniques. Here, the optimal control problems were approximated by the optimal parameter selection problems using control parameterization, which can be regarded as finite-dimensional optimization problems, where the control heights and switching times of the piece-wise constant function were taken as the decision variables ^[22]. To the best of our knowledge, most researchers have used piecewise continuous functions to approximate the control function in their models, while few have provided the discrete dosage forms for multiple variables. On the contrary, it is notably important to study the optimal control of unconstrained, piecewise continuous tumor models; see, for example ^{[16,17,23,24,25]}. Control parameterization discretizes the control and provides the optimal solution for the most effective control problem, leaving a certain buffer time for the next period of control changes, making the control input more accurate ^{[26,27,28,29,30]}. In addition, the optimal control problem using this method has the form of continuous state and discrete control. Compared to discrete states, this method obtains better optimal solutions. Compared to continuous control, this method has lower limitations on physical devices and stronger operability when applied in clinical treatment. Based on changes in the number of tumor cells in the body, drug control is more accurate.

In this paper, we apply the optimal control theory to derive the discrete drug dosing strategies for constrained optimal control problems, thereby focusing on two treatment modalities. Specifically, we aim to minimize the number of tumor cells and the drug's toxicity to the body, maximize the number of effector cells, and consider a variety of realistic scenarios. These include constraints on the quantity of a single drug dose and the total amount of the drug input, thus reflecting the body's tolerance to drug exposure.

The organization of this paper is as follows: in Section 2, we briefly describe the nonlinear system of ordinary differential equations that model the tumor-immune cell interactions; Section 3 presents the optimal control problem, which includes both equality and inequality constraints; in Section 4, we derive the gradient formulas and develop a numerical approach to seek the optimal control strategies; Section 5 provides numerical simulations, both with and without constraints, thereby illustrating the drug injection doses based on piece-wise constant forms; and finally, we provide the corresponding conclusions in Section 6.

2.   Tumor-immune model

In this section, we will briefly introduce the interaction model for the initial cancer dynamics as proposed by Kuznetsov et al. ^[31]. The immune response to tumors is usually dominated by cytotoxic T lymphocytes (CTLs) and natural killer cells (NKs). The interaction between effector cells ($EC$) and tumor cells ($TC$) in vitro can be described by a kinetic scheme (Figure 1). where $EC$, $TC$, $C$, $EC^*$, and $TC^*$ are the local concentrations of effector cells, tumor cells, effector cell-tumor cell conjugates, inactivated effector cells, and "lethally hit" $TC$ cells, respectively. The parameters $k_1$, $k_{-1}$, $k_2$, and $k_3$ are nonnegative kinetic constants: $k_1$ and $k_{-1}$ describe the rates of binding of $EC$ to $TC$ and detachment of $EC$ from $TC$ without damaging cells; $k_2$ is the rate at which $EC-TC$ interactions irreversibly program $TC$ for lysis; and $k_3$ is the rate at which $EC-TC$ interactions inactivate $EC$.

They established a complex nonlinear tumor model that involved a system of ordinary differential equations to describe the tumor-immune cell interactions as follows:

where $EC$, $TC$, and $C$ represent the numbers of unbound $EC$, unbound $TC$, and $EC-TC$ complexes at the tumor site, respectively; $EC^{*}$ and $TC^{*}$ indicate the number of inactivated EC and lethally hit TC at the tumor site, respectively. The total number of unhit $TC$ cells at the tumor site is denoted as $TC_{tol} = TC+C$. The parameter $s$ represents the standard rate at which mature $EC$ cells enter the $TC$ locus, which is independent of any enhancement due to the presence of $TCs$. The parameters $d _{1}$, $d _{2}$ and $d _{3}$ are positive constants that represent the rates of elimination for the $EC$, $EC^{*}$, and $TC^{*}$ cells, respectively, because of their destruction or migration from the $TC$ locus. It is assumed that the tumor is not metastatic, and therefore, there is no migration of $TCs$ or $EC-TC$ complexes from the primary site. The parameter $a$ represents the maximum growth rate of the $TC$ population, which encompasses both proliferation and death of these cells. The parameter $b$ is the maximum carrying capacity of the biological environment for the $TC$ population. The function $F(C, TC)$ characterises the rate at which the cytotoxic effector cells accumulate in the localized region of the $TC$ because of the existence of tumors.

Based on ^[16], model (2.1) can be simplified as follows:

where the function $u_1(t)$ serves as an input function for IL-2. The external administration of IL-2 is modeled by the control function $u_1(t)$, which facilitates the influx of effector cells and boosts the activity of the immune cells. The function $u_2(t)$ acts as an input function for ACI. The external administration of ACI is represented by the control function $u_2(t)$, which inhibits tumor production and directly targets the tumor cells at the tumor site. When $u_1(t) = 1$ and $u_2(t) = 0$, model (2.2) is a simplified model of the tumor without treatment.

3.   Optimal control problem

In this section, an objective function is defined for treating a tumor over a fixed treatment cycle. The therapeutic goal is to minimize the number of tumor cells and reduce the drug toxicity in terms of the side effects for the patient, while maximizing the benefit for the immune effector cells. Here, we explore the impact of controlling IL-2 ($u_1(t)$) and ACI ($u_2(t)$) treatments on the dynamics of the state equation (2.2) during tumor treatments.

Thus, the cost function for the optimal control problem is defined as follows:

where the constants $B_1$ and $B_2$ represent the costs associated with drug administration, which are also known as the weighting factors of the control. These primarily reflect the level of drug toxicity to the body. Given that drugs such as ACI and IL-2 can exhibit high toxicities, a lower toxicity is preferable. This reasoning underlies the use of a positive coefficient for the quadratic terms of the control in the objective function, as it aims to minimize the negative impact of drug toxicity.

Now, the following optimal control problem subject to system (2.2) is stated.

Problem P1: Given system (2.2), choose the optimal controls $u_1$ and $u_2$ such that the cost function (3.1) is minimized.

In the treatment process, considerations such as the patient's drug tolerance and the severity of the toxicity side effects were not taken into account. In contrast, a study by Goedegebuure et al. ^[32] involved 16 valuable patients who experienced varying levels of toxicity. This study demonstrated that different drug dosages can lead to a range of side effects, including, but not limited to, flu-like symptoms (fever, nausea, vomiting, and diarrhea), capillary leak syndrome, oliguria, transient confusion, elevated liver and kidney function, thrombocytopenia, headache, hallucinations, erythema rash, and other drug-related side effects. This indicates that even moderate drug concentrations can cause significant side effects in patients undergoing anti-cancer treatments, yet this may be considered a reasonable trade-off. Additionally, the study highlighted that the anti-tumor effects of combination drugs are significant when the dosage of each drug is kept below the maximum tolerated dose (MTD) for either drug ^[33]. This suggests that limitations on the drug input, including the amount administered per unit of time and the total amount administered during the treatment cycle, should be carefully considered in the patient treatment plans.

First, we will limit the dosages of the anticancer drugs ACI and IL-2, denoted as $\boldsymbol{u}(t): = [u_1(t), u_2(t)]$, administered per unit of time, i.e.,

where $\boldsymbol{u}^{max}: = [{u}_1^{max}, {u}_2^{max}]$ are positive constants, which represent the maximum allowable dosages for each drug in a combination therapy with each maximum value. This ensures that the dosages are within a range that patients can tolerate in terms of the side effects from a single administration of the drugs.

Second, another measure of drug toxicity considers the product of drug concentration and the duration of administration, which can be mathematically represented as the integral of the drug concentration over a designated period. When this period encompasses either the entire treatment course or a treatment cycle, the primary concern becomes the total cumulative toxicity. This approach accounts for the overlap of side effects experienced by the patient at any given time and the overall bodily mechanisms that may be adversely impacted by a prolonged exposure to side effects. When determining the total dosage of a specific drug administered, we have the following:

where $\boldsymbol{u}_{\text{sfixed}}$ are positive constants, which represent the total dosage of the drug inputs over a time period $T$.

Thus, the constrained optimal control problem subject to system (2.2) can be defined as follows.

Problem P2: Given system (2.2), choose the optimal controls $u_1$ and $u_2$ such that the cost function (3.1) is minimized subject to the constraints (3.2) and (3.3).

4.   Control parametrization

As is well-known, optimal control problems can be efficiently solved by the control parameterization technique, in which the control function is approximated by a piece-wise constant or linear functions ^[21].

The controls $u_1(t)$ and $u_2(t)$ can be approximated by the following:

where

and the new decision variables $h_{k}^{j}, k\in \{1, 2\}, j\in \{0, 1, \ldots, M\}$ satisfy the following constraint:

Let $t_{0} = 0$, $t_{M+1} = T.$ The given switching times $t_{j}, j\in \{0, 1, \ldots, M\},$ satisfy the following:

Let

where

The corresponding cost function (3.1) becomes the following:

Furthermore, system (2.2) changes to the following form:

with the initial condition

With these in mind, Problem P1 can be approximated by the following problem.

Problem PP1: Given system (4.4) with initial condition (4.5), choose an optimal $h_{k}^{j}, k\in \{1, 2\}, j\in \{0, 1, \ldots, M\}$ such that the cost function (4.3) is minimized.

5.   Gradient computation

In essence, problem (PP1) is a mathematical programming problem. It is well known that gradient-based optimization techniques are highly effective in solving such a problem. For this, the gradients of the cost function (4.3) are needed. In this section, we will present these required gradients.

Define

For $t \in\left[t_{l}, t_{l+1}\right), l\in \{j, \ldots, M\}, j\in\{0, 1, \ldots, M\}, k\in \{1, 2\},$ consider the following auxiliary dynamic systems

with the initial condition

where

Theorem 1. The gradients of the cost function (4.3) are as follows:

Proof. The proof of Theorem 1 is similar to the proof of Theorem 7.2.2 in ^[21]. □

Define

Then, Problem P2 can be approximated by the following problem.

Problem PP2: Given system (2.2), choose $h_{k}^{j}, k\in \{1, 2\}, j\in \{0, 1, \ldots, M\}$, such that the cost function (3.1) is minimized subject to the constraints (4.2) and (5.4).

In addition, the gradients of the constraints in Problem P2 are given as follows.

Theorem 2. The gradients of the constraints (5.4) are as follows:

Proof. The proof of Theorem 2 is obtained based on the chain rule. □

Generally, to maintain consistency in the time intervals before and after, we default to using the same time interval as previously mentioned when calculating the gradient of the constraints, based on the time interval of the cost function. This approach yields the gradients of the constraints, which can be solved using a standard nonlinear programming algorithm.

6.   Numerical simulation

We will now conduct numerical simulations of { Problems PP1} and PP2 to validate the analytical gradient results presented in the previous sections. To demonstrate the accuracy of our numerical simulations, we have selected parameter sets from both the tumor model and the optimal control problem. These parameters are listed in Table 1.

First, we obtain the numerical solution for system (2.2) using the fourth-order Runge-Kutta algorithm. Then, an algorithm based on the control parameterization is developed as follows:

Step1: Obtain $\psi_{k}^{j}(t)$ by solving system (5.1) with the initial constraint (5.2).

Step2: Compute the solution of systems (4.4) and (4.5) by solving system (4.4) forward in time from $t = 0$ to $t = T$ with the initial condition (4.5).

Step3: Compute the corresponding value of ${\boldsymbol{g}_{i}}(\boldsymbol{h})$ via (5.4) and compute the gradient value of the constraints (5.4) via (5.5).

Step4: Find the optimal value by using the gradient formulae from Theorems 1 and 2 in conjunction with the nonlinear programming software–Optimization Toolbox in MATLAB.

Finally, we simulate the tumor treatment process over a specified time period. To gain multiple perspectives on what is actually occurring throughout the process, we present the following simulations of the scenario.

6.1.   First scenario: no treatment

We simulate a scenario where the effector and tumor cells naturally engage in a conflict without any treatment applied. This simulation provides a geometric visualization that helps us understand the behavior of the tumor system, thus leading to insights into potential optimal control strategies. Figure 2 illustrates the dynamics of this untreated scenario. Within a fixed time horizon from 0 to 200, the effector cells oscillate, and similarly, the tumor cells fluctuate in a pattern of ebb and flow. This demonstrates that the cancer does not remit on its own without intervention.

6.2.   Second scenario: the effect of unconstrained combination therapy

We now explore the interaction between the effector and tumor cells under the influence of two combined drugs, regardless of the total dosage administered. In the treated scenario, the number of tumor cells reaches a therapeutic plateau around day 50, thus indicating a significant reduction compared to the untreated case. Figure 2 clearly shows that the combination of ACI and IL-2 markedly reduces the tumor cell and effectively suppresses their proliferation, as indicated by the yellow solid line in the first graph of Figure 3, as compared to the second graph of Figure 2. Additionally, the number of effector cells decreases due to the side effects of the drugs, as shown by the blue solid line in the first graph of Figure 3. Furthermore, the tumor cells exhibit a slight increase towards the end of the cycle noted by the upward trend at the end of the yellow solid line in the first graph of Figure 3, suggesting that the initial therapeutic effectiveness diminishes over time without controlling the total drug dosage. This reduction in drug efficacy is primarily attributed to the tumor cells developing resistance to the drugs as the treatment progresses.

6.3.   Third scenario: adding a fixed total amount of drug injection

In this section, we explore the drug toxicity by providing numerical examples from five scenarios to illustrate how the total drug input affects the number of tumor cells and immune effector cells during the treatment when two drugs are combined, thereby highlighting the accumulation of drug toxicity. We focus on the impact of the total drug amount on the tumor treatment, with the dosage of each individual drug contingent upon the constraints set by the total drug input.

Based on numerical calculations of the total drug amount, we will present a range of corresponding single drug inputs. This approach underscores that individual drug dosages also significantly affect the patient treatment outcomes.

First case: when total IL-2 and total ACI are the same, both are large

Consider the following case:

We now examine the outcomes when the total dosage of both drugs is maintained at high levels throughout the treatment cycle. The number of tumor cells was significantly reduced when high concentrations of IL-2 and ACI were administered, compared to scenarios where the total drug dosage was unrestricted. Specifically, the number of tumor cells effectively converged to zero from the first day of drug administration until about day 40, as shown by the yellow solid line in the first graph of Figure 4. Simultaneously, the number of state cells bound to tumor cells decreased, thus showing an exposure to drug toxicity, but at a relatively steady rate illustrated by the blue solid line in the first graph of Figure 4.

In this scenario, the higher single-drug dosages imply that patients with a lower tolerance for side effects can still employ this aggressive treatment strategy to quickly control the tumor progression. Initially, the ACI dosage was high, thus leading to a significant reduction in the number of tumor cell. Then, once effective, the dosage of ACI was adjusted downward, as depicted in the third graph of Figure 4. Meanwhile, the single dosage of IL-2 decreased more gradually, as shown in the second graph of Figure 4. This strategy highlights the tailored adjustment of drug dosages in response to the tumor dynamics and the patient response.

Second case: when total IL-2 and total ACI are the same, both are small

Consider the following case:

We examine a scenario where both drugs are administered at lower total dosages. Compared to the unrestricted drug treatment, administering low levels of IL-2 and ACI initially led to a fluctuation in the tumor cell count, rising and then falling from the first day of drug administration to around day 110, with an overall upward trend during this period. Subsequently, the number of tumor cells linearly decreased from day 110 to day 160 as the concentrations of IL-2 and ACI were increased. However, even at their lowest, the tumor cell counts remained above 100, 000, and a rapid rebound occurred thereafter as seen by the significant curvature at the end of the solid yellow line in the first graph of Figure 5, compared to the end of the solid yellow line in Figure 3.

In the earlier stages, when the total amount of administered drug was low, the impact on the state cells was minimal, both for the single and the total drug inputs. However, after the single dosage of ACI was increased, the number of state cells significantly declined, which is a side effect that indicates that higher drug concentrations lead to more pronounced adverse effects. Importantly, the extent of these side effects varies from person to person.

Third case: when the total IL-2 is high and the total ACI is low

Consider the following case:

In this scenario, we examine the effects when the total input of IL-2 is high and that of ACI is low. This scenario displayed a distinctive pattern in a tumor cell count increase, then decrease, then increase upon administering high IL-2 and low ACI. Initially, the number of tumor cells increased alongside the IL-2 input from the first day of drug administration up to about 50 days. After this period, despite the continued high intake of IL-2, the decline in the number of tumor cell was very gradual, thus highlighting a delayed response to the drug treatment and illustrating that the delay effects vary between different drugs.

Post 100 days, there was a rapid decrease in the number of tumor cells reaching a minimum, which was sustained over approximately 40 days, thus suggesting the combined effectiveness of the drugs, particularly as the ACI therapy was introduced, thereby enhancing the suppression of tumor growth while also resulting in a decrease in the effector cells. However, beyond 180 days, when the drug inputs were reduced to their lowest, the tumor rebounded sharply (as depicted by the distinct curvature at the end of the solid yellow line in the first graph of Figure 6), thereby surpassing the initial tumor count within 20 days. This suggests that the effectiveness of the treatment was compromised by a inadequate drug input, although this point is not the central focus of our current discussion.

Fourth case: low total IL-2 and high total ACI

Consider the following case:

In this scenario, we explored a treatment regimen characterized by a low overall intake of IL-2 and a high overall input of ACI. We observed that the effects of such a treatment were largely similar to those in the first scenario, with a notable difference being a significantly higher count of state cells in this case. Biologically, this can be mainly attributed to the higher input of IL-2 in the second scenario, which was able to boost the number of state cells within the first 40 days when low concentrations of IL-2 and high concentrations of ACI were administered, as depicted in Figure 7.

This finding suggests that clinically, to minimize drug side effects, the input of IL-2 can be appropriately reduced while using higher concentrations of ACI to achieve an enhanced therapeutic effect. Additionally, it indicates that different drugs have varying tolerance ranges, and exceeding the upper limit of the input, while yielding consistent therapeutic effects, can cause additional harm that may not justify the benefits.

6.4.   Fourth scenario: involves adjusting a specific drug concentration and executing a gradual administration of another drug

The primary goal of this scenario is to evaluate the impact of maintaining a fixed concentration of one drug while continuously adjusting the dosage of another on the tumor-immune microenvironment within the patient. The drug can effectively function by administering the correct dose to achieve a peak blood concentration. In this context, the ensuing numerical simulations not only identify the primary and secondary roles of the drugs during the treatment cycle, but also demonstrate the combined effect of both drugs more accurately, as depicted in Figures 8 and 9.

First case: maintaining the concentration of IL-2

ACI is given in increments of 10, 30, 40, 50 and 100.

The dynamics of the effector and tumor cells in Figure 8 demonstrate significant variations when the IL-2 concentration is maintained and ACI is administered in a gradient. These figures clearly show that the number of cancer cells significantly vary at different ACI concentrations. Initially, when ACI doses are limited to 30, 40, and 50, there is little variation; however, more variation becomes apparent towards the end of the cycle.

Second case: maintain the ACI concentration while administering IL-2 at doses of 10, 30, 40, 50 and 100.

Figure 9 illustrates the changing ratio of effector cells to tumor cells when the ACI concentration is maintained constant and IL-2 is administered in varying doses. Throughout the course of the treatment, there is a gradual decline in the number of effector cell, while the reduction in tumor cells is more pronounced. At this level of ACI, the drug's effectiveness remains relatively stable during the initial phases of treatment. However, with lower doses of IL-2, there is a quicker resurgence in cancer cell proliferation.

7.   Conclusions

In this study, we addressed the optimal control problem with the goal of maximizing the population of effector immune cells, minimizing the tumor cell counts, and reducing the drug toxicity. The objective of this paper was to deepen our understanding of tumor dynamics by analyzing interactions between immune effector cells, tumor cells, and their microenvironment, given a fixed total drug input. To assess the cumulative toxicity of the treatment, we introduced two external drugs (ACI and IL-2) administered over the total cycle time. The primary contribution of this work includes setting a constraint on the total drug input and establishing a discrete dosing sequence to represent the administration of anti-cancer drugs. Initially, we defined the optimal control problems, with and without constraints, based on a specified tumor model. Theoretically, we enhanced the model's descriptive power by incorporating an auxiliary system that provided gradient information for the state equations using variational and chain rule methods, thereby deriving a gradient formulation for the objective function analyzed in this paper. Subsequently, we converted the equation constraint into a standard form suitable for optimal control problems. Then, we derived the gradient formulation for the specific constraints imposed. Numerically, we explored the dynamics of the tumor system under various scenarios, including the effects of a combined ACI and IL-2 treatment in unconstrained settings and under different constraints, among other considerations.

Some research directions are outlined as follows:

1). Given the delay between drug administration and its action within the patient's specific microenvironment, we could explore this by incorporating a time-delay in the state equations.

2). Stochastic factors will be involved in the tumor-immune system to account for unpredictable scenarios that might emerge during the course of treatment.

3). Implementing secondary control measures will be considered during the drug treatment, particularly when a drug is well-tolerated or when the patient's internal environment becomes more complex. This involves selecting the optimal timing and appropriate drug for a second round of control treatment.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China under Grant 12161076 and Grant 12271307; in part by the Shandong Province Natural Science Foundation of China under Grant ZR2023MA054; in part by the Nature Science Foundation of Liaoning Province of China under Grant 2024-MS-015; in part by the Fundamental Research Funds for the Central Universities under Grant 3132024196, and Grant DUT22LAB305; in part by the China Postdoctoral Science Foundation under Grant 2019M661073; in part by the Xinghai Project of Dalian Maritime University.

Conflict of interest

The authors declare there is no conflict of interest.