# Probabilistically Certified Region of Attraction of a Tumor Growth Model with Combined Chemo-and Immunotherapy

Kaouther Moussa, Mirko Fiacchini, Mazen Alamir

► To cite this version:

Kaouther Moussa, Mirko Fiacchini, Mazen Alamir. Probabilistically Certified Region of Attraction of a Tumor Growth Model with Combined Chemo-and Immunotherapy. International Journal of Robust and Nonlinear Control, 2022, 32 (11), pp.6539-6555. 10.1002/rnc.6158. hal-03651829

HAL Id: hal-03651829

https://hal.science/hal-03651829

Submitted on 26 Apr 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Probabilistically Certified Region of Attraction of a
Tumor Growth Model with Combined Chemo- and
Immunotherapy

Kaouther Moussa¹, Mirko Fiacchini², and Mazen Alamir²

¹Univ. Polytechnique Hauts-de-France, LAMIH, CNRS, UMR 8201, F-59313 Valenciennes, France, kaouther.moussa@uphf.fr

²Univ. Grenoble Alpes, CNRS, Grenoble INP, GIPSA-lab, 38000 Grenoble, France, {mirko.fiacchini, mazen.alamir}@grenoble-inp.fr

April 13, 2022

Abstract

The recent progress in immunology lead to a considerable interest in modeling cancer dynamics in order to better understand and analyze such complex systems. Many works have been carried out in order to design cancer treatment protocols using mathematical models. One of the main complexities of such models is the presence of different types of uncertainties, which remains less considered in the literature. This paper deals with the estimation of regions of attraction (RoAs) under parametric uncertainties for a cancer growth model with combined therapies. We propose a framework of probabilistic certification, based on the randomized methods, in order to derive probabilistically certified RoAs of a cancer growth model. The model considered in this paper describes the interaction between a tumor and the immune system in presence of a combined chemo- and immunotherapy treatment, with considerations on pharmacokinetics and pharmacodynamics (PK/PD) of both treatments.

The last decades witnessed a considerable progress in experimental and clinical im-
munology as well as in modeling the immune system dynamics [13]. The progress in
cancer dynamics modeling motivated researchers to apply control approaches in order to
schedule cancer treatments using optimal control strategies. We can find in the literature
many works regarding the application of optimal control approaches on cancer treatment
problems. For instance [10], where optimal protocols for anti-angiogenic therapy were
investigated, or [6] where linear controls were designed for a tumor-immune interactions
model with chemotherapy delivery. However, only few works addressed the problem of
handling parametric uncertainties. One can cite for example, [1] where a robust feedback
scheme is proposed to schedule anti-angiogenic treatment combined with chemotherapy,
[18] where an H∞-based robust control was applied to the same model and [2] where a
general framework for probabilistic certification of cancer therapies was proposed.

The estimation of the region of attraction for cancer models is an interesting problem since it provides a set of possible initial conditions (tumor volume and immune density for example) that can be driven to a desired target set (benign region). This problem becomes complex when dealing with nonlinear systems and even more challenging for uncertain systems. There are some works which dealt with the problem of estimating the RoA for cancer models, we cite for example [8] and [29], where the authors proposed different Lyapunov functions based approaches, to estimate the domain of attraction of the tumor free equilibrium point corresponding to autonomous cancer growth models, where no therapies are considered, see also [22] and references therein. However, only few works considered model uncertainties, in particular in [24], an iterative method to estimate the robust RoA was presented. Furthermore, the work in [9] dealt with the estimation of regions of attraction for cancer models, using rational Lyapunov functions and switching control strategies. The latter work has been further investigated in [7] in order to take into account model disturbances, based on input to state (ISS) Lyapunov functions. Moreover, the model proposed in [9] has been used for data fitting in [14] for real-life experiments. Finally, it is important to highlight the fact that the estimation of robust RoAs is based on the worst-case scenario analysis leading to potentially pessimistic design, this because the worst-case is considered no matter how small its probability of occurrence is.

The presence of parametric uncertainties can drastically affect the efficiency of a nom-
inal controller as well as the size of the estimated RoAs. In this work, we propose a
framework to probabilistically certify the existence of a control structure that drives the
states corresponding to tumor cells and immune density from an initial state set to a
certified target set. This probabilistic certification framework is based on the randomized
methods proposed in [3] and [4], which, unlike the robust classical design, avoids focusing
on few unlikely extremely bad scenarios allowing to overcome the conservatism of the
robust RoA design. The methodology that we propose in this paper consists mainly of
two steps. Firstly, we derive an ordered sequence of sets and their control strategy such
that the states can be driven from a set to the previous one with a certain probabilistic
guarantee. The appropriate choice of the first set allows to insure that the union of the
sets is a probabilistically certified approximation of the RoA. The second step consists in
providing a global certification on the probability of convergence to the initial certified
target set. The randomized methods have been already used in [15] in order to determine
the stability region for nonlinear deterministic systems, without taking into consideration
parametric uncertainties.

The model that we investigate here is a modified version of the classical Stepanova
one [27] that has been extensively used in the literature, we cite for example [11], [19]
and [21] where optimal control methodologies were proposed to schedule chemo- and
immunotherapy administration profiles. Furthermore, [25] proposed a model predictive
control scheme to design chemo- and immunotherapy administration schedules. In [26],
the authors proposed a robust model predictive control scheme, in order to consider direct drug targeting pharmacokinetic uncertainties as well as system model mismatches. In this paper, we model the concentration of the chemotherapy agent in the plasma and the tumor site via a pharmacokinetics compartmental model, we also model the pharmacokinetics of immunotherapy, as well as the pharmacodynamics of both drugs. Although the classical Stepanova model has been widely used, it has never been investigated in the literature to estimate the controlled region of attraction of its corresponding tumor free equilibrium. Therefore, we aim at pointing out the importance of uncertainties considerations in the RoA estimation for such models.

This paper is organized as follows: In Section 1, the dynamical cancer model and the problem of RoA probabilistic certification are introduced. Section 2 recalls the randomized algorithms approach for probabilistic certification. In Section 3, a framework for RoA probabilistic certification is proposed, based on the randomized methods presented in [3] and [4]. In Section 4, the proposed RoA probabilistic certification framework is applied to the considered cancer model. Finally, Section 5 summarizes the contribution that we present in this paper.

# 1 Problem statement

The following nonlinear dynamical system describes the interaction between a tumor and the immune system in presence of chemotherapy and immunotherapy treatments:

$$ 
\begin{aligned}
\dot{x}_1 &= \mu_C \left(1 - \frac{x_1}{x_\infty}\right) x_1 - \gamma_X x_1 x_2 - \kappa_X \left(\frac{x_4^{\gamma_c}}{x_4^{\gamma_c} + C_{50c}^{\gamma_c}}\right) x_1 \\
\dot{x}_2 &= \mu_I (1 - \beta x_1) x_1 x_2 + \alpha_Y - \delta_Y x_2 - \eta_Y x_2 x_3 + \kappa_Y \left(\frac{x_5^{\gamma_i}}{x_5^{\gamma_i} + C_{50i}^{\gamma_i}}\right) x_2 \\
\dot{x}_3 &= -(k_1 + k_2) x_3 + s_1 \frac{u_1(t)}{V_1} \\
\dot{x}_4 &= k_{12} \frac{V_1}{V_2} x_3 - k_2 x_4 \\
\dot{x}_5 &= -c_i x_5 + s_2 u_2(t), \\
x(0) &= (x_1(0), x_2(0), x_3(0), x_4(0), x_5(0)) = x_0,
\end{aligned}
\quad (1) $$

where the variables are defined as follows:

$x_1$ The number of tumor cells ($10^6$ cells)

$x_2$ The density of effector immune cells (ECs) (dimensionless)

$x_3$ The concentration of chemotherapy in the plasma ($\mu g \cdot mL^{-1}$)

$x_4$ The concentration of chemotherapy in the tumor effect site ($\mu g \cdot mL^{-1}$)

$x_5$ The concentration of immunotherapy in the immune cells site ($\text{mg} \cdot \text{mL}^{-1}$)

$u_1$ The dosage of the chemotherapeutic agent (Etoposide) ($\mu g \cdot \text{day}^{-1}$)

$u_2$ The dosage of the immuno-stimulator agent (Nivolumab) ($\text{mg} \cdot \text{mL}^{-1} \cdot \text{day}^{-1}$)
<figure><img src="image_2.png" /><figcaption>Figure 1: Temporal open-loop control structure for each cycle, in black and yellow, respectively, the immunotherapy and the chemotherapy profiles.</figcaption></figure>

This model is an extension of the model presented in [11] that has the advantage of being a low dimensional system that nevertheless includes the main aspects of cancer-immune interactions, and it has been widely used in the literature for cancer drug scheduling. In many models it is assumed that the drug concentration is equal to its dosage which is an oversimplification [20]. Therefore, we revisited the model proposed in [11] by adding a pharmacokinetic (PK) compartment (involving $x_3$ and $x_4$) that allows to model the concentration of chemotherapy in the plasma and the tumor effect site. Furthermore, we incorporated to this model the pharmacodynamics of chemotherapy using a Hill function, the equations and the parameters values have been taken from [5] and [16] for the Etopside drug. Other types of chemotherapy can be considered such as Pegylated Lipo-somal Doxorubicin, whose main corresponding medical parameters have been presented in [12]. Moreover, we extended the model with the pharmacokinetics of immunotherapy which contains only one compartment represented by $x_5$ which is the concentration of immunotherapy in the immune cells site, since we consider that the immunotherapy does not have a direct inhibition effect on the tumor. The pharmacodynamics effects have also been incorporated using a Hill function. The parameters values related to the PK/PD of immunotherapy have been taken from [17].

Table 1 summarizes the definitions of the model parameters and their nominal values. Furthermore, the parameters $s_1$ and $s_2$ allow to scale the drugs effects. Note that we consider examples of possible treatment protocols, nevertheless, it is worth emphasizing that in this paper, we focus on the assessment of a methodology that remains applicable for different parameters values, treatment strategies and equations.

Let's denote by $x = (x_1, x_2, x_3, x_4, x_5)$ and $u = (u_1, u_2)$ respectively, the state and the control input vectors. In this paper, we consider a cycle-based treatment, where
the drugs are injected following $N_C$ therapeutic cycles. Each cycle having two phases, a hospitalization period lasting 5 days, where the patient receives one infusion per day, and a rest period where the patient recovers. Figure 1 shows a typical temporal combined control structure, the different notations in this figure are defined as follows:

<table><caption>Table 1: Definitions and nominal values of the parameters used in model (1).</caption><thead><tr><th>Parameter</th><th>Definition</th><th>Numerical value</th><th>Unit</th><th>Reference</th></tr></thead><tbody><tr><td>μ<sub>C</sub></td><td>Tumor growth rate</td><td>1.0078</td><td>10<sup>7</sup>·day<sup>-1</sup></td><td>[23]</td></tr><tr><td>μ<sub>I</sub></td><td>Tumor stimulated<br>proliferation rate</td><td>0.0029</td><td>day<sup>-1</sup></td><td>[23]</td></tr><tr><td>α<sub>Y</sub></td><td>Rate of immune<br>cells influx</td><td>0.0827</td><td>day<sup>-1</sup></td><td>[23]</td></tr><tr><td>β<sub>Y</sub></td><td>Inverse threshold</td><td>0.0040</td><td>(-)</td><td>[23]</td></tr><tr><td>γ<sub>X</sub></td><td>Interaction rate</td><td>1</td><td>10<sup>7</sup>·day<sup>-1</sup></td><td>[11]</td></tr><tr><td>δ<sub>Y</sub></td><td>Death rate</td><td>0.1873</td><td>day<sup>-1</sup></td><td>[23]</td></tr><tr><td>κ<sub>X</sub></td><td>Chemotherapeutic<br>killing parameter</td><td>1</td><td>10<sup>7</sup>·day<sup>-1</sup></td><td>[11]</td></tr><tr><td>κ<sub>Y</sub></td><td>Immunotherapy<br>administration parameter</td><td>1</td><td>10<sup>7</sup>·day<sup>-1</sup></td><td>[11]</td></tr><tr><td>x<sub>∞</sub></td><td>Fixed carrying capacity</td><td>780</td><td>10<sup>6</sup>·cells</td><td>[11]</td></tr><tr><td>η<sub>Y</sub></td><td>Chemo-induced loss<br>on immune cells</td><td>1</td><td>mL/(μg·day)</td><td>[23]</td></tr><tr><td>γ<sub>c</sub></td><td>Patient response/ resistance<br>to Etopside</td><td>2.5</td><td>(-)</td><td>[5]</td></tr><tr><td>C<sub>50c</sub></td><td>Half-effect concentration<br>of Etopside</td><td>10</td><td>μg·mL<sup>-1</sup></td><td>[5]</td></tr><tr><td>k<sub>1</sub></td><td>Chemotherapy clearance<br>rate from plasma</td><td>1.6</td><td>day<sup>-1</sup></td><td>[5], [16]</td></tr><tr><td>k<sub>2</sub></td><td>Chemotherapy clearance<br>rate from tumor</td><td>0.8</td><td>day<sup>-1</sup></td><td>[5], [16]</td></tr><tr><td>k<sub>12</sub></td><td>Link process between<br>plasma and the tumor</td><td>0.4</td><td>day<sup>-1</sup></td><td>[5], [16]</td></tr><tr><td>V<sub>1</sub></td><td>Plasma volume</td><td>25</td><td>L</td><td>[5], [16]</td></tr><tr><td>V<sub>2</sub></td><td>Effect site volume</td><td>15</td><td>L</td><td>[5], [16]</td></tr><tr><td>γ<sub>i</sub></td><td>Patient response/ resistance<br>to Nivolumab</td><td>2.5</td><td>(-)</td><td>[17]</td></tr><tr><td>C<sub>50i</sub></td><td>Half-effect concentration<br>of Nivolumab</td><td>32·10<sup>-6</sup></td><td>mg·mL<sup>-1</sup></td><td>[17]</td></tr><tr><td>c<sub>i</sub></td><td>Clearance rate of Nicolumab</td><td>11.6/24</td><td>day<sup>-1</sup></td><td>[17]</td></tr></tbody></table>
<table><caption>Table 2: Definitions of the different notations in Figure 1.</caption><thead><tr><th>Parameter</th><th>Definition</th></tr></thead><tbody><tr><td>σ<sub>I</sub></td><td>Duration of immunotherapy infusion</td></tr><tr><td>σ<sub>C</sub></td><td>Duration of chemotherapy infusion</td></tr><tr><td>d<sub>I</sub></td><td>Concentration of immunotherapy infusion</td></tr><tr><td>d<sub>C</sub></td><td>Concentration of chemotherapy infusion</td></tr><tr><td>T = 5</td><td>Hospitalization duration (days)</td></tr><tr><td>T<sub>c</sub> = 15</td><td>Cycle duration (days)</td></tr></tbody></table>

For a given treatment cycle, the therapeutic profile considered in this paper is completely defined by the following control parametrization θ:

$$ \theta = (\nu_C, \sigma_C, d_C, \sigma_I, d_I). \qquad (2) $$

In cancer treatment design, we usually have some constraints to satisfy, they can be defined either on the states or on the control inputs. These constraints allow to prevent from drug toxicity and excessive immune weakening. In this paper, we consider the following constraints for all $t \in [0, T]$, with $T \in \mathbb{R}_+$:

$$ x_2(t) \geq c, \text{ with } c \in \mathbb{R}_+, \qquad (3) $$

$$ 0 \leq u_1(t) \leq 1, \qquad (4) $$

$$ 0 \leq u_2(t) \leq 1, \qquad (5) $$

where (3) is a health constraint on the minimal density of immune cells. The constraints on $u_1(t)$ and $u_2(t)$ for all $t$, see (4) and (5), are normalized drug toxicity constraints, they can be satisfied by properly choosing the parametrization $\theta$ (namely, $d_C$ and $d_I$) of the control input $u$. Therefore, we will consider only the constraint (3), since the satisfaction of the other constraints can be monitored by a proper choice of $\theta$.

The uncontrolled model (1) (for $u = (0, 0)$) has two locally asymptotically stable equilibriums points. The macroscopic malignant equilibrium is $x_m = (766.4, 0.018, 0, 0, 0)$ and the benign one is $x_b = (41.45, 0.954, 0, 0, 0)$. In general, the objective of the treatment is to drive the state initial conditions to the region of attraction of the benign equilibrium (safe region), without constraints violation. We are specifically interested in characterizing the set of initial conditions (tumor volume and immune density) from which the trajectories of (1) can be driven to the safe region under parametric uncertainties.

In [23], we proposed a methodology to characterize the controlled region of attraction of model (1) with bang-bang controls (without pharmacokinetics). Then, we used this approach to derive and estimate the robust region of attraction. In this paper, we propose to derive a probabilistically certified RoA for model (1), that is based on chance-constrained problems, tolerating some constraints violations provided that their corresponding probability is small enough.
**Definition 1** We denote by $\Omega_0 \subseteq \mathbb{R}_+^5$ a set of initial conditions $x_0$, such that the state trajectories $(x_1, x_2)$ corresponding to model (1), when no control is applied (i.e. $u = (0,0)$), converge to their corresponding benign equilibrium in spite of all uncertainties realizations, with a given confidence probability. This set can be seen as a targeted safe set for each treatment. Note that since $u = (0,0)$ and the initial conditions $x_3(0), x_4(0), x_5(0)$ are all equal to zero, the set $\Omega_0$ is technically determined in $\mathbb{R}_+^2$.

The region of uncertainties is defined such that the parameters are normally distributed in the interval $[0.9p_{nom}, 1.1p_{nom}]$, where $p_{nom}$ stands for the vector of nominal parameters that will be properly defined in Section 4.

**Problem 1 (Estimation of a probabilistically certified RoA)** We aim at computing a sequence of sets $\{\Omega_k\}_{k=1}^{N_C}$, for $N_C$ therapeutic cycles. Those sets are determined in the space of $(x_1, x_2)$, representing respectively, the cancer burden (defined by the number of cancer cells) and the ECs density, such that, in the family of control parametrizations that we consider, there exists a therapeutic protocol that drives, with a desired probability, the states from $\Omega_{k+1}$ to $\bigcup_{j=0}^k \Omega_j$ without safety constraints violations.

We denote by $\hat{\Omega}_0^{p_{nom}}$ an estimation of the region of attraction of the benign equilibrium for $u = (0,0)$, when nominal model parameters (in Table 1) are considered.

## 2 Overview on probabilistic certification

The randomized algorithms were presented in [28], [3] and [4] in order to solve optimization problems with probabilistic constraints satisfaction. In contrast to the standard robust control design, which is based on the worst-case scenario analysis leading hence to potentially pessimistic design, the randomized methods provide the possibility to avoid focusing on the worst scenarios if their probability of occurrence is small. Therefore, this framework is very interesting from the cancer treatment point of view, since the latter involves many uncertainties that have to be considered.

This section aims at briefly recalling the main key-points of the randomized methods that are important for the assessment of the approach that we propose in this paper, in which we present a framework of estimation of probabilistically certified region of attraction for a cancer therapies dynamical model.

Consider the following optimization problem :

$$ 
\begin{array}{ll}
\min_{\theta \in \Theta} & J(\theta) \\
\text{s.t.} & \forall p \quad g_c(\theta, p) = 0,
\end{array}
\qquad (6)
$$
where $\theta \in \Theta \subset \mathbb{R}^{n_{\theta}}$ is the decision variable (which can be a parametrization of a control law) and $p$ is the uncertainties vector following the probability measure $\mathcal{P}$ defined in the set $\mathbb{P}$ (the vector $p$ can contain for example model parameters that are considered to be uncertain), $J$ is the cost to be minimized. In terms of control design for dynamical systems, the cost $J$ can involve the states, the input variables, their respective integrals with respect to time or any combination of these indicators. Finally, $g_c$ is an indicator function on the violation of some given constraints and is defined as follows:

$$g_c(\theta, p) := \begin{cases} 0 & \text{if all the constraints are satisfied} \\ 1 & \text{otherwise} \end{cases}$$

The randomized method consists in replacing the original hard problem in (6) by the following relaxed problem:

$$\begin{array}{ll} \min_{\theta \in \Theta} & J(\theta) \\ \text{s.t.} & \mathrm{Pr}_{\mathcal{P}}\{g_c(\theta, p) = 1\} \le \eta, \end{array} \qquad (7)$$

where the constraint is on the probability of constraints violation, giving therefore a soft constraint in the sense that we can accept a value of $\theta$ which minimizes the cost $J$, even if the constraints are violated for some realizations of $p$, provided that the probability of these realizations is less than or equal to $\eta$ (small enough). Even though the constraint in (7) simplifies the previous constraint in (6), the computation of the violation probability remains expensive. Authors in [3] and [4] proposed a simplification which consists in replacing the probability by the mean value of $g_c$ over $N_p$ drawn independent identically distributed (i.i.d.) samples of $p$ in $\mathbb{P}$ according to the probability distribution $\mathcal{P}$. Therefore, the simplified optimization problem is the following:

$$\begin{array}{ll} \min_{\theta \in \Theta} & J(\theta) \\ \text{s.t.} & \displaystyle\frac{\sum_{i=1}^{N_p} g_c(\theta, p^{(i)})}{N_p} \le \frac{m}{N_p}, \end{array} \qquad (8)$$

where $m$ is the number of allowed constraints violations. In [3] and [4], several bounds on $N_p$ are given such that the fulfillment of the constraint in (8) implies that the probability condition in (7) is satisfied with a confidence probability greater than or equal to $1 - \delta$. The bounds that are derived on $N_p$ involve the precision $\eta$ and the confidence of fulfillment $\delta$.

In this paper, we are interested in specific control structures, since cancer treatment schedules are often defined by cycles with a hospitalization period where the patient receives several drug infusions and a rest period for recovery. Therefore, it is more adequate in this case to consider that the controls are parametrized by a discrete variable $\theta$ with cardinality $n_{\theta} \in \mathbb{N}$. This choice of $\theta$ simplifies the optimization problem (7), since it can be solved by a simple enumeration. In this case, the following proposition from [4] holds:

**Proposition 1** Let $m \in \mathbb{N}$ be any integer representing the number of accepted failures.
Let $\delta \in (0, 1)$ be a targeted confidence parameter. Take $N_p$ satisfying

$$N_p \geq \frac{1}{\eta} \left( m + \ln \left( \frac{n_{\theta}}{\delta} \right) + \left( 2m \ln \left( \frac{n_{\theta}}{\delta} \right) \right)^{\frac{1}{2}} \right) \quad (9)$$

then any solution of (8) in which $\{p^{(j)}\}_{j=1}^N$ are i.i.d. following the probability distribution $P$ satisfies the constraint in (7) with a probability greater than or equal to $1 - \delta$

The inequality (9) is mathematically based on the binomial distribution. In this section, we presented a concise overview of the basic theoretical aspects of this methodology, the readers interested in further details might refer to [4].

It is interesting to notice that the bound on $N_p$ provided by Proposition 1 does not depend on the dimension of $p$ which is useful when having many uncertain parameters in the certification problem. Furthermore, as we can see in Table 3, since the confidence parameter $\delta$ affects the bound logarithmically, we can have a highly confident certification with a tractable number of random samples.

<table><caption>Table 3: The evolution of the number of samples $N_p$ required to achieve the certification, with respect to the confidence design parameter $\delta$ and the number of control parametrizations $n_{\theta}$, for $\eta = 10^{-2}$ and $m = 1$.</caption><thead><tr><th>n<sub>θ</sub></th><th>δ = 0.1</th><th>δ = 0.01</th><th>δ = 0.001</th></tr></thead><tbody><tr><td>10</td><td>864</td><td>1162</td><td>1451</td></tr><tr><td>100</td><td>1162</td><td>1451</td><td>1732</td></tr><tr><td>1000</td><td>1451</td><td>1732</td><td>2008</td></tr><tr><td>10000</td><td>1732</td><td>2008</td><td>2280</td></tr></tbody></table>

Furthermore, for a specific desired confidence parameter $\delta = 10^{-3}$, Table 4 provides an idea on the evolution of the number of trials $N_p$ that should be performed for each possible control law $\theta$, with respect to the precision parameter $\eta$ and the number of control parametrizations $n_{\theta}$. Therefore, the total number of simulations is $N_{sim} = N_p \cdot n_{\theta}$.

<table><caption>Table 4: The evolution of the number of samples $N_p$ required to achieve the certification, with respect to the precision design parameter $\eta$ and the number of control parametrizations $n_{\theta}$, for $\delta = 10^{-3}$ and $m = 1$.</caption><thead><tr><th>n<sub>θ</sub></th><th>η = 0.1</th><th>η = 0.01</th><th>η = 0.001</th></tr></thead><tbody><tr><td>10</td><td>146</td><td>1451</td><td>14503</td></tr><tr><td>100</td><td>174</td><td>1732</td><td>17312</td></tr><tr><td>1000</td><td>201</td><td>2008</td><td>20073</td></tr><tr><td>10000</td><td>228</td><td>2280</td><td>22796</td></tr></tbody></table>
This approach provides a powerful pragmatic tool allowing to certify control strategies. In [2], a randomized methods based framework for probabilistic certification of feedback control strategies has been proposed for a combined cancer therapy model.

# 3 Probabilistic certification of ROA

In this section, we will establish a framework of RoA probabilistic certification, based on the randomized methods presented in the previous section. We propose to use this general framework in order to probabilistically certify the existence of a control structure which allows to drive initial states from a given set to a target set under parametric uncertainties.

Let us rewrite system (1) into the following form:

$$ \dot{x} = F(x, u, p), \qquad (10) $$

where $p$ is the vector of parameters that model (1) involves. Furthermore, we consider that the variables of system (10) are subject to the following constraints:

$$ x \in \mathbb{X}, \quad x(T) \in \Omega, \quad u \in \mathbb{U}, \qquad (11) $$

where the sets $\mathbb{X}$ and $\mathbb{U}$ are defined as follows:

$$ \mathbb{X} = \{x \in \mathbb{R}_+^5 | x_2 \ge c\}, \qquad (12) $$

$$ \mathbb{U} = \{(u_1, u_2) \in \mathbb{R}_+^2 | 0 \le u_1, u_2 \le 1\}, \qquad (13) $$

and $\Omega \subseteq \mathbb{R}_+^5$ represents the target set that will certified to be safe for each treatment cycle.

**Remark 1** Note that since the states $x_3, x_4$ and $x_5$ vanish after the treatment period, when $u = (0,0)$, the certified safe set $\Omega$ will be determined technically in the space $(x_1, x_2)$ which is a subset of $\mathbb{R}_+^2$.

As previously mentioned, we consider that the control inputs are parametrized by a vector $\theta$ which lies in a discrete set $\Theta$ with cardinality $n_{\Theta} \in \mathbb{N}$. This choice of $\theta$ fits particularly to the case of cancer therapy design, since some of the parameters involved in the treatment scheduling are naturally quantified.

Suppose that the parameters vector $p$ is a random variable following the probability distribution $P$ that we denote $p \sim P$. Given a set $\Gamma \subseteq \mathbb{R}_+^5$, containing initial conditions $x_0$ (Remark (1) holds for $\Gamma$ since the initial conditions $x_3(0), x_4(0), x_5(0)$ are null), and a parameterization of the input $\theta \in \Theta$, let's consider the following optimization problem:

$$ \begin{array}{ll} \min_{\theta \in \Theta} & J(\theta) \\ \text{s.t.} & \forall (x_0, p) \in (\Gamma \times \mathbb{P}) \quad g_c(\theta, x_0, p, \Omega) = 0, \end{array} \qquad (14) $$
where $J(\theta)$ is a cost function to be minimized. In terms of cancer treatment design, this function can be a combination of many objectives that one seeks to achieve, for example reducing the quantity of injected drugs, to prevent from toxicity, or reducing the duty cycle in order to reduce the hospitalization duration. $g_c$ is the failure indicator function, defined on the state trajectories of (10). Note that $g_c$ depends on $\Omega$ in order to emphasize the fact that the terminal constraint, expressing that the health indicators belong to the certified target set, is a part of the failure indicator that will be properly defined in the sequel. The function $g_c$ is deterministic such that, for a given initial state $x_0$, an input parametrization $\theta$ and a model parameters vector $p \in \mathbb{P}$, it is equal to one if the constraints (11) are violated, zero otherwise. Problem (14) then aims at selecting the optimal control strategy such that no constraints violation occurs.

As previously explained, the randomized method consists in replacing the original problem in (14) by the following chance-constrained problem tolerating some violations:

$$ \begin{array}{ll} \min_{\theta \in \Theta} & J(\theta) \\ \text{s.t.} & \mathrm{Pr}_{\mathcal{X}_0(\Gamma) \times \mathcal{P}}\{g_c(\theta, x_0, p, \Omega) = 1\} \le \eta, \end{array} \quad (15) $$

where the constraint is on the probability of violation, with respect to the distribution of $x_0$ on $\Gamma$, that we denote $\mathcal{X}_0(\Gamma)$, and $p \sim \mathcal{P}$. This problem gives therefore a chance-constrained formulation in the sense that we can accept a vector $\theta$ which minimizes the cost $J$, even if the constraints are violated for some realizations of $(x_0, p)$, provided that the probability of these violations is lower than $\eta$, hence small enough.

Since problem (15) is hard to solve, it can be simplified into the following problem, employing the empirical mean instead of the probability of the constraints violation, Given $\Gamma \subseteq \mathbb{R}^5$:

$$ \begin{array}{ll} \min_{\theta \in \Theta} & J(\theta) \\ \text{s.t.} & \displaystyle\sum_{i=1}^{N} g_c\left(\theta, x_0^{(i)}, p^{(i)}, \Omega\right) \le m, \\ & (x_0, p)^{(i)} \sim (\mathcal{X}_0(\Gamma) \times \mathcal{P}), \forall i = 1, \dots, N, \end{array} \quad (16) $$

where $m$ is the maximum number of allowed constraints violation. If $N$ satisfies the condition in (9), then the solution of (16) satisfies the constraint in problem (15) with a probability higher than $1 - \delta$.

As previously explained, problem (16) can be solved by simple enumeration which means that, for each possible treatment protocol defined by $\theta$, we simulate the model (1) for $N$ samples of the uncertain parameter vector $p$, in order to compute

$$ \sum_{i=1}^{N} g_c\left(\theta, x_0^{(i)}, p^{(i)}, \Omega\right) $$

and thereby to select the treatment protocols satisfying the latter constraint. The solution is the treatment protocol defined by $\theta$ corresponding to the minimum cost $J(\theta)$. In the
next section, we will explain how the iterative resolution of problems of the type (16)
allows one to generate a sequence of sets {$\Omega_k$}$_{k=1}^{N_C}$ such that the constraints violation on
passing from $\Omega_{k+1}$ to $\bigcup_{j=0}^k \Omega_j$ is smaller then $\eta$ with a certain desired confidence probability
$1-\delta$.

**3.1 Algorithm for RoA estimation**

Given a target set $\Omega \subseteq \mathbb{R}_+^5$, our objective is to certify that the set $\Gamma$ is such that there exists a control parametrization $\theta$, for which at least $100 \cdot (1 - \eta)\%$ of the trajectories of (10), generated by the distributions of the initial states $x_0 \in \Gamma$ and the uncertain parameters $p$, converge to $\Omega$ at time $T$, while satisfying constraints (11), with a confidence higher than $1 - \delta$. Any solution of (16) defines a local control strategy that satisfies the constraints while minimizing the cost $J(\theta)$.

$\Gamma$ generator

We suppose that we have a generator of sets $\Gamma \subseteq \mathbb{R}_+^5$ with a parametrized geometry providing a family of nested potential sets $\Gamma$, then we can compute the biggest one that is probabilistically certified through (16). In the case under study, we consider that the sets $\Gamma$ have a polytopic form.

Therefore, starting from $\Omega_0$ which is in the certified region of attraction of many be-
nign equilibriums without therapies, an iterative procedure can be designed to generate
the sequence $\{\Omega_k\}_{k=0}^{N_C}$ such that the trajectories starting in $\Omega_{k+1}$ end in $\bigcup_{j=0}^k \Omega_j$ with the
desired probability and without violating the constraints. In particular, we will consider
sequences of sets such that $\Omega_k \cap \Omega_{k+1} = \emptyset$. Then, we keep doing this certification process
until given $\Omega_{k-1}$, the set $\Omega_k$ is empty. Once the RoA probabilistic certification algorithm
terminates, the candidate to be a probabilisitically certified RoA is the set $\Omega_C = \bigcup_{i=1}^{N_C} \Omega_i$.

Note that, if $x_0 \in \Omega_k$, for $k = 1, \dots, N_C$, this means that the trajectory of length $T$ will end in $\bigcup_{j=0}^{k-1} \Omega_j$ without violating the constraint with a certain probability, but no direct probabilistic guarantee is given regarding the convergence to the set $\Omega_0$.

It is not straightforward to derive a probabilistic bound on driving the states directly
from the last set of the sequence $\Omega_{N_C}$ to $\Omega_0$. This is because the latter probability in-
volves the accuracy and confidence parameters, $\eta$ and $\delta$. Another reason is that there is
no guarantee that, given the initial state distribution $\mathcal{X}_0(\Omega_k)$, the distribution of the state
at the end of the k-th therapeutic cycle is $\mathcal{X}_0(\Omega_{k-1})$, for which the probabilistic validation is performed. However, after deriving the sequence of certified sets, we can approximate the probability of driving the states from $\Omega_{N_C}$ to $\Omega_0$, with the corresponding certified control strategy, using Monte-Carlo simulations.

<caption>**Algorithm 1 Sequence of probabilistically certified sets**</caption>

**Input:** $\Omega_0$

$k \leftarrow 0$

**while** $\Omega_k \neq \emptyset$ **do**

$\Omega \leftarrow \bigcup_{j=0}^{k} \Omega_j$

**repeat**

Generate non-empty $\Gamma$ such that $\Omega \cap \Gamma = \emptyset$

**until** (16) is unfeasible for $\Gamma$

**if** $\exists \Gamma$ such that (16) is feasible **then**

$\Omega_k \leftarrow \Gamma$

**else**

$\Omega_k \leftarrow \emptyset$

**end if**

$k \leftarrow k + 1$

**end while**

$N_C \leftarrow k - 1$

**Output:** $\Omega_C \leftarrow \bigcup_{i=0}^{N_C} \Omega_i$

Finally, by using Algorithm 1, we can obtain a sequence of certified sets, such that the output is the candidate to be a probabilistically certified RoA $\Omega_C$.

Note that since we focus on the estimation of the RoA for a specified control parametrization, the use of a cost function is not relevant although it could have been used in the case where some of the parameters defining the control are kept free.

# 4 Probabilistically certified RoA for a cancer model

As previously explained, considering $N_C$ treatment cycles, our objective consists in estimating the probabilistically certified RoA of model (1) that we denote $\Omega_C$. To this end, we certify a sequence of successive disjoint sets such that their union is the candidate to be a probabilistically certified RoA.

Moreover, the temporal control profiles that we consider correspond only to the hospitalization period (see Figure 1), meaning that the rest period is not included in the decision variable $\theta$ defined in Section 1, since we assume that this parameter can be esti-
mated afterwards depending on the health conditions of the patient.

Therefore, we propose a feedback control strategy that can be seen in an implicit way, such that at the end of each therapy period, we measure the states (patient health and tumor volume) and depending on the certified set $\Omega_k$ where this measure lies, we can estimate the maximal possible recovery time ($T_c - T$) that the patient can take. At the end of the rest period, the certified therapy corresponding to this set is then applied, we keep doing this process until we reach the safe region $\Omega_0$.

The initial condition $x_0$ is assumed to be uniformly distributed in the set $\Gamma$ while the following parameters vector is assumed to be uncertain:

$$p = (\mu_C, \mu_I, \gamma_X, \kappa_X, \kappa_Y, \delta_Y, \alpha_Y, \beta_Y, \eta_Y) \quad (17)$$

and normally distributed in the following interval:

$$[0.9p_{nom}, 1.1p_{nom}], \quad (18)$$

where $p_{nom}$ is the nominal value of each parameter, previously presented in Table 1.

The failure indicator function, which determines whether the constraints (3)–(5) are satisfied or not, is defined on $x(t|x_0, p, \theta)$ which is the state trajectory of (1) for a given control parametrization $\theta$ and a random sample of $x_0$ and $p$. We denote by $x(T|x_0, p, \theta)$ the state trajectory evaluated at the end of the hospitalization period. Therefore, the failure indicator is defined as:

$$g_c(\theta, x_0, p, \Omega) := \begin{cases} 0 & \text{if } x_2(t|x_0, p, \theta) \ge c \ \forall t \text{ and } x(T|x_0, p, \theta) \in \Omega \\ 1 & \text{otherwise} \end{cases} \quad (19)$$

where $\Omega$ is a probabilistically certified target set which can be seen as the safe region to reach at the end of the cycle.

Using Algorithm 1, we can derive a sequence of probabilistically certified sets providing the probabilistically certified RoA. Firstly, we need to derive an initial target set $\Omega_0$, in order to initialize the certification algorithm.

## 4.1 Probabilistically certified initial target set $\Omega_0$

**Definition 2** Given $p \in \mathbb{P}$ (drawn according to the probability distribution $\mathcal{P}$), $x_0$ following a uniform distribution on $\Omega_{start} \subseteq \mathbb{R}_+^5$, and $\Omega_{end} \subseteq \mathbb{R}_+^5$ a certified target set in a neighborhood of benign equilibrium points of (1), generated by the realizations of $p$ according to the probability distribution $\mathcal{P}$. We define the certified safe reachability by:

$$\Pr_{\mathcal{U}(\Omega_{start}) \times \mathcal{P}} \{x_2(t|x_0, p) \ge c, \forall t > 0 \text{ and } x(T|x_0, p) \in \Omega_{end} \} > 1 - \eta. \quad (20)$$
This means that the state trajectories having initial conditions in $\Omega_{start}$ converge with a given probability to the set $\Omega_{end}$, after a time horizon $T$, in spite of all parametric uncertainties and without constraints violations. Note that $\Omega_{start}$ can be equal to $\Omega_{end}$, in this case $\Omega_{start}$ is slightly different than a probabilistically certified invariant set, since we don't require that the trajectories starting in $\Omega_{start}$ stay in it, but rather to converge to it after some time $T$.

Firstly, we certify $\Omega_{eq}$ such that:

$$ \mathrm{Pr}_{\mathcal{U}(\Omega_{eq}) \times \mathcal{P}}\{x_2(t|x_0, p) \geq c, \forall t > 0 \text{ and } x(T|x_0, p) \in \Omega_{eq}\} > 1 - \eta. \quad (21) $$

Then, given $p$ and $x_0$ following a uniform distribution on $\Omega_0$, that we denote $\mathcal{U}(\Omega_0)$, $\Omega_0$ is determined such that:

$$ \mathrm{Pr}_{\mathcal{U}(\Omega_0) \times \mathcal{P}}\{x_2(t|x_0, p) \geq c, \forall t > 0 \text{ and } x(T|x_0, p) \in \Omega_{eq}\} > 1 - \eta, \quad (22) $$

for a given time $T$. Note that the set $\Omega_{eq}$ is derived to be used as a target set for the determination of $\Omega_0$. In order to provide an estimation of $\Omega_{eq}$, we draw the distribution of the benign equilibriums of model (1) for many parameters vector samples (selected according to the probability distribution $\mathcal{P}$). Then, we choose a geometry for $\Omega_{eq}$ surrounding the benign equilibriums of the sample shown in Figure 2. Finally, we expand this set until (21) is not satisfied.

<figure><img src="image_3.png" /><figcaption>Figure 2: Probabilistically certified sets $\Omega_0$ for different horizons $T$ (days). The red circles are the benign equilibrium points of model (1) for 1000 parameters vector samples.</figcaption></figure>

After finding a proper geometry for the set $\Omega_{eq}$ such that it satisfies (21), we use Algorithm 1 in order to provide an estimation of the certified set $\Omega_0$. Note that in this
case $\mathcal{X}_0(\Gamma)$ corresponds to $\mathcal{U}(\Omega_0)$ since we assume that $x_0$ is uniformly distributed on $\Omega_0$, and the target set for the states at time $T$ denoted $\Omega$ in the definition of $g_c$ corresponds to $\Omega_{eq}$. Furthermore, since we deal with an uncontrolled system, (7) turns out to be a feasibility problem, where we need only to guarantee the probability condition in (22) by using the empirical mean over $g_c$ for $N$ i.i.d. samples of $(x_0, p)$ mentioned in (16), with $\theta = 0$, $n_\Theta = 1$, and with the bound $N$ given by (9), for $m = 1$, $\eta = 10^{-2}$ and $\delta = 10^{-3}$.

<figure><img src="image_4.png" /><figcaption>Figure 3: Phase portrait of (1), the colored trajectories represent the states $(x_1, x_2)$ for different initial conditions, the blue dots correspond to the 3 equilibrium points of the system, estimated nominal uncontrolled RoA $\hat{\Omega}_0^{p_{nom}}$ in dashed cyan and the estimated certified initial target set $\Omega_0$ for $T = 60$ in blue.</figcaption></figure>

We assume that the set $\Omega_0$ to be certified has the same geometry as the estimated nominal uncontrolled region of attraction $\hat{\Omega}_0^{p_{nom}}$ (derived in [23]) that we shrink until (22) is not satisfied given the confidence probability $1 - \delta$. There is clearly no guarantee that the set $\Omega_0$ that we obtain is the biggest possible certified set, however, in this case, proving the existence of a set $\Omega_0$ satisfying (22) is enough, since $\Omega_0$ is only used as a target set for the Algorithm 1 allowing therefore to compute the sequence of certified sets.

Figure 2 shows the probabilistically certified RoA of the benign equilibriums $\Omega_{eq}$, the estimated uncontrolled nominal region of attraction $\hat{\Omega}_0^{p_{nom}}$ and the initial probabilistically certified target set $\Omega_0$ for different $T$. Figure 3 shows the phase portrait of (1) with both
the estimated nominal RoA $\hat{\Omega}_0^{p_{nom}}$ without control, and the certified initial target set $\Omega_0$ for $T = 60$. We can see that the $\Omega_0$ is smaller than $\hat{\Omega}_0^{p_{nom}}$ which shows the effects of parametric uncertainties consideration.

### 4.1.1 Validation of the estimation of $\Omega_0$

In order to validate the estimation of the target set $\Omega_0$, we carry out 5000 Monte-Carlo simulations by randomly selecting the initial states as well as the model parameters according to their respective probability distributions. We can notice that in Figure 4 there are only 11 trajectories that converge to the malignant equilibrium, violating thereby the specified constraints. This corresponds to 99.78% of successful trajectories, validating therefore the imposed probabilistic bound on $\eta$.

<figure><img src="image_5.png" /><figcaption>Figure 4: Monte-Carlo simulations to validate the certified target set $\Omega_0$.</figcaption></figure>

## 4.2 Probabilistically certified region of attraction $\Omega_C$

We denote by $\Omega_C$ the probabilistically certified region of attraction of system (1). We initialize Algorithm 1 with $\Omega_0$ in order to derive the sequence of probabilistically certified sets providing the certified RoA for model (1).

We consider that the decision variable $\theta$ is defined by the following variables:

$$ \left\{ \begin{array}{l} \sigma_I \in \{0, 0.16, 0.32, 0.48, 0.64, 0.8\}, \\ \sigma_C = 0.5, \quad \nu_C = 0.2, \\ d_I \in \{0, 0.25, 0.5, 0.75, 1\}, \\ d_C \in \{0, 0.11, 0.22, 0.33, 0.44, 0.56, 0.67, 0.78, 0.89, 1\}. \end{array} \right. $$

### Complexity analysis and computation time

The cardinality of $\Theta$ is $n_{\Theta} = 300$ giving the bound $N \ge 1863$ according to (9), for $m = 1$, $\eta = 10^{-2}$ and $\delta = 10^{-3}$. The number of simulations to be performed for each set certi-
fication is $N_{sim} = N \cdot n_{\Theta} = 558900$. The required computational time to perform $N_{sim}$ simulations is less than 10 min using Matlab coder toolbox. Therefore, 1 simulation requires around 1.1 ms on a computer with an Intel(R) Core(TM) i5-10310U and a 2.21 GHz CPU. However, it is important to highlight two main points, the first one is the fact that the number of simulations in our case is due to the choice of solving the optimization problem by enumeration, since the cardinality of the input set $\Theta$ is relatively low. It is definitely possible to solve such problems using iterative algorithms (such as gradient descent) in order to reduce the computation time. The second point is that since the problem is completely scalable, it is possible to solve it using parallel computing which considerably reduces the computational time.

Figure 5 shows the 3 certified cycles for $T = 5$ days obtained using Algorithm 1, nominal and robust RoAs that have been estimated using the method presented in [23], where bang-bang control strategies were considered. We can see that, as the number of cycles increases, the certified RoA gets closer to the robust controlled one denoted $\Omega_R$. Even if the certified RoA remains smaller, it is important to recall that both the nominal and the robust RoAs in [23] have been estimated using continuous infusions of drugs, and do not consider the PK/PD effects. Note also that changing the type of the sets can change the size of the certified RoA.

<figure><img src="image_6.png" /><figcaption>Figure 5: Probabilistically certified RoAs for 3 administration cycles, the estimated nominal RoA $\Omega_u^{p_{nom}}$ and the estimated robust RoA $\Omega_R$ using continuous drugs infusions.</figcaption></figure>
4.2.1 Validation of the estimation of Ω$_{C}$

We approximated the probability of driving the states from Ω₃ to Ω₀ using 5000 Monte-Carlo simulations. We obtained that 99.4% of the trajectories of (1) having initial conditions in Ω₃ converge to Ω₀ using the probabilistic certified control strategies that we derived. Figure 6 shows the phase portrait of the 5000 Monte-Carlo trajectories. We can notice that only a small part of these trajectories violate the minimal constraint on immune cells density. The trajectories violating this constraints are presented in Figure 7.

<figure><img src="image_7.png" /><figcaption>Figure 6: Monte-Carlo simulations to validate the certified sequence of controls with their respective sets, the blue polytope is the set Ω₃ where the initial states were selected.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Monte-Carlo simulations, the trajectories that violate the minimal constraint on immune cells density, the blue polytope is the set Ω₃ where the initial states were selected.</figcaption></figure>
# 5 Conclusion

In this paper, we presented a framework of probabilistic certification for regions of attraction which is based on the randomized methods, allowing to overcome the conservatism of worst-case robust approaches by proposing a tractable problem with probabilistic constraints. This framework has been used to derive a certified region of attraction for a cancer growth model. Furthermore, we provided a validation on the probability of driving the states to the certified safe target set with its corresponding control strategy.

The main advantages of this framework is that it is less conservative than the worst-case approach, since it is more tolerant to constraints violations in the presence of uncertainties. In addition to this, the methodology that we presented in this paper provides the control strategy corresponding to each certified initial states set, which allowed to validate the estimations using Monte-Carlo simulations. Furthermore, this approach is flexible in terms of computational complexity, since we can considerably reduce the computational time by solving the optimization problems iteratively instead of using enumeration or by using parallel computing since the problem is completely scalable.

The probabilistic certification of regions of attraction can be seen as a tool to tune the several parameters of the treatment protocols by properly choosing the model parameters and their distributions, the geometry of the regions of attraction to be certified and the control parametrization. All these choices impact the size of the certified region of attraction.

One of the future works in terms of RoA certification methodology is to enrich the set generator to consider different geometries, in order to approximate the biggest certifiable RoA. From a medical point of view, a future perspective would be to consider the synergy between the different injected drugs and to model their combined effects on the different compartments, in order to solve more challenging and seemingly realistic problems.

# Data Availability Statement

Data sharing is not applicable to this article as no new data were created in this study.

# References

[1] ALAMIR, M. Robust feedback design for combined therapy of cancer. *Optimal Control Applications and Methods* 35, 1 (2014), 77–88.

[2] ALAMIR, M. On probabilistic certification of combined cancer therapies using strongly uncertain models. *Journal of Theoretical Biology* 384 (2015), 59–69.
[3] ALAMO, T., TEMPO, R., AND CAMACHO, E. F. Randomized strategies for probabilistic solutions of uncertain feasibility and optimization problems. *IEEE Transactions on Automatic Control* 54, 11 (2009), 2545–2559.

[4] ALAMO, T., TEMPO, R., LUQUE, A., AND RAMIREZ, D. R. Randomized methods for design of uncertain systems: Sample complexity and sequential algorithms. *Automatica* 52 (2015), 160–172.

[5] BUKKURI, A. Optimal control analysis of combined chemotherapy-immunotherapy treatment regimens in a PKPD cancer evolution model. *Biomath* 9 (2020).

[6] DE PILLIS, L. G., GU, W., FISTER, K. R., HEAD, T., MAPLES, K., MURUGAN, A., NEAL, T., AND YOSHIDA, K. Chemotherapy for tumors: An analysis of the dynamics and a study of quadratic and linear optimal controls. *Mathematical Biosciences* 209, 1 (2007), 292–315.

[7] DOBAN, A. I. *Stability domains computation and stabilization of nonlinear systems: implications for biological systems*. PhD thesis, Eindhoven University of Technology, 2016.

[8] DOBAN, A. I., AND LAZAR, M. Domain of attraction computation for tumor dynamics. In *53rd IEEE Conference on Decision and Control* (2014), pp. 6987–6992.

[9] DOBAN, A. I., AND LAZAR, M. A switching control law approach for cancer immunotherapy of an evolutionary tumor growth model. *Mathematical Biosciences* 284 (2017), 40–50.

[10] D'ONOFRIO, A., AND CERRAI, P. A bi-parametric model for the tumour angiogenesis and antiangiogenesis therapy. *Mathematical and Computer Modelling* 49, 5–6 (2009), 1156–1163.

[11] D'ONOFRIO, A., LEDZEWICZ, U., AND SCHÄTTLER, H. *On the Dynamics of Tumor-Immune System Interactions and Combined Chemo- and Immunotherapy*. Springer Milan, 2012, pp. 249–266.

[12] DREXLER, D. A., FERENCI, T., LOVRICS, A., AND KOVÁCS, L. Modeling of tumor growth incorporating the effect of pegylated liposomal doxorubicin. In *INES 2019 : IEEE 23rd International Conference on Intelligent Engineering Systems* (2019).

[13] EFTIMIE, R., GILLARD, J. J., AND CANTRELL, D. A. Mathematical models for immunology: current state of the art and future research directions. *Bulletin of Mathematical Biology* 78, 10 (2016), 2091–2134.

[14] GUILLEREY, C., NAKAMURA, K., PICHLER, A., BARKAUSKAS, D., KRUMEICH, S., STANNARD, K., MILES, K., HARJUNPÄÄ, H., YU, Y., CASEY, M., DOBAN, A., LAZAR, M., HARTEL, G., SMITH, D., VUCKOVIC, S., TENG, M., BERGSAGEL, P., CHESI, M., HILL, G., MARTINET, L., AND SMYTH, M. Chemotherapy followed by anti-CD137 mAb immunotherapy improves disease control in a mouse myeloma model. *JCI Insight* 5, 14 (2019).
[15] HOKAYEM, P. F., MASTELLONE, S., AND ABDALLAH, C. T. *Nonlinear systems stability via random and quasi-random methods*. Springer London, London, 2006, pp. 365–379.

[16] ILIADIS, A., AND BARBOLOSI, D. Optimizing drug regimens in cancer chemotherapy by an efficacy-toxicity mathematical model. *Computers and Biomedical Research 33* (2000), 211–226.

[17] IONESCU, C. M., GHITA, M., COPOT, D., DEROM, E., AND VERELLEN, D. A minimal PKPD interaction model for evaluating synergy effects of combined NSCLC therapies. *Journal of Clinical Medicine 9* (2020), 1–17.

[18] KOVACS, L., SZELES, A., SAPI, J., DREXLER, D. A., RUDAS, I., HARMATI, I., AND SAPI, Z. Model-based angiogenic inhibition of tumor growth using modern robust control method. *Computer Methods and Programs in Biomedicine 114*, 3 (2014), e98–e110.

[19] LEDZEWICZ, U., AND FARAJI, M. On optimal protocols for combinations of chemo-and immunotherapy. In *51st IEEE Conference on Decision and Control* (Hawaii, USA, December 2015), pp. 7492–7497.

[20] LEDZEWICZ, U., AND SCHÄTTLER, H. Optimal controls for a model with pharmacokinetics maximizing bone marrow in cancer chemotherapy. *Mathematical Biosciences 206*, 2 (2007), 320–342.

[21] LEDZEWICZ, U., AND SCHÄTTLER, H. On the role of the objective in the optimization of compartmental models for biomedical therapies. *Journal of Optimization Theory and Applications 187* (2020), 305–335.

[22] MEROLA, A., COSENTINO, C., AND AMATO, F. An insight into tumor dormancy equilibrium via the analysis of its domain of attraction. *Biomedical Signal Processing and Control 3*, 3 (2008), 212–219.

[23] MOUSSA, K., FIACCHINI, M., AND ALAMIR, M. Robust domain of attraction estimation for a tumor growth model. *Applied Mathematics and Computation 410* (2021), 126482.

[24] RIAH, R., FIACCHINI, M., AND ALAMIR, M. Iterative method for estimating the robust domains of attraction of non-linear systems: Application to cancer chemotherapy model with parametric uncertainties. *European Journal of Control 47* (2019), 64–73.

[25] SHARIFI, N., OZGOLI, S., AND RAMEZANI, A. Multiple model predictive control for optimal drug administration of mixed immunotherapy and chemotherapy of tumours. *Computer Methods and Programs in Biomedicine 144* (2017), 13–19.

[26] SHARIFI, N., ZHOU, Y., AND HOLMES, G., C. Y. Overcoming channel uncertainties in touchable molecular communication for direct drug targeting assisted immunotherapy. *IEEE Transactions on NanoBioscience 19*, 2 (2020), 249–258.
[27] STEPANOVA, N. Course of the immune reaction during the development of a malignant tumour. *Biophysics* (1980).

[28] TEMPO, R., CALAFIORE, G., AND DABBENE, F. *Randomized algorithms for analysis and control of uncertain systems, with applications (2nd ed.)*. London: Springer-Verlag, 2013.

[29] ZAREI, M., JAVADI, K., AND KALHOR, A. Perturbed tumor immunotherapy domain of attraction estimation via the arc-length function. In *2018 25th National and 3rd International Iranian Conference on Biomedical Engineering (ICBME)* (2018), pp. 1-6.