Lyès Mellal

David Folio

Karim Belharet

Antoine Ferreira

Modeling of Optimal Targeted Therapies using Drug-Loaded Magnetic Nanoparticles for the Liver Cancer

Keywords: Magnetic Resonance Navigation, Tumor Growth model, Optimal Tumor Control, Optimal drug delivery, Locoregional therapies

HAL is

I. INTRODUCTION

Cancer is known as one of the major causes of morbidity and death worldwide. Especially, the liver cancer continues to be a major cause of mortality, and its incidence is still increasing [START_REF] Siegel | Cancer statistics[END_REF]. Due to a lack of donors, it is even more rarely possible to perform a liver transplantation or to consider surgical procedures [START_REF] Llovet | Arterial embolisation or chemoembolisation versus symptomatic treatment in patients with unresectable hepatocellular carcinoma: a randomised controlled trial[END_REF], [START_REF] Forner | Hepatocellular carcinoma[END_REF]. Therefore, transarterial locoregional interventions, such as transcatheter arterial chemoembolization (TACE) or radioembolization (TARE), are procedures that provide a significant survival benefit [START_REF] Llovet | Arterial embolisation or chemoembolisation versus symptomatic treatment in patients with unresectable hepatocellular carcinoma: a randomised controlled trial[END_REF]- [START_REF] Kennedy | Recommendations for radioembolization of hepatic malignancies using yttrium-90 microsphere brachytherapy: a consensus panel report from the radioembolization brachytherapy oncology consortium[END_REF]. These minimally invasive procedures aim to restrict a tumor's blood supply (embolization) together with a drug delivery (radioembolization or chemoembolization) into an artery supplying a tumor. In particular, the use of drug eluting beads (microspheres) that carry the chemotherapeutic agent ensures the delivery a lower side-effect than systemic chemotherapy [START_REF] Lewis | DC bead: in vitro characterization of a drug-delivery device for transarterial chemoembolization[END_REF], [START_REF] Varela | Chemoembolization of hepatocellular carcinoma with drug eluting beads: efficacy and doxorubicin pharmacokinetics[END_REF]. The success of both TACE or TARE procedures requires a critical mass of drug eluting beads or radioactive Yttrium-90 microspheres to implant in vessels at the tumor periphery, as illustrated in Fig. 1. Due to anatomical size constraints the delivery catheter is limited in terms of accessibility, improper catheter placement and insufficient extrahepatic artery embolization. As a This paper extends the work presented at the IEEE International Conference on Intelligent Robotics and Systems, Hamburg, Germany, 2015. L. Mellal, D. Folio and A. Ferreira is with the INSA Centre Val de Loire, Univ. Orléans, K. Belharet is with the HEI Centre, PRISME Laboratory EA 4229, Bourges, F-18020, France. E-mail: antoine.ferreira@insa-cvl.fr.

Microvascular networks

Size: 100-300µm

Fig. 1. Schematic representation of the principle locoregional therapy for the liver cancer. Access is gained from the femoral artery, and the catheter is passing through the abdominal aorta to the hepatic (liver) artery for delivering the therapeutic agents (drug eluting beads) close to the tumor. Then, the beads are cutting off the blood supply (embolization), as well as delivering the drugs to the tumor (chemo-or radio-therapy). Without proper control some therapeutic agents could end up in the systemic flow and cause injuries.

consequence, the current procedures demonstrated insufficient deposition results close to the tumor vessels promoting tumor growth. Furthermore, some therapeutic agents end up in the systemic flow, and destroy healthy liver tissue via ischemia, radiation damage and/or drug inflammation [START_REF] Matsumoto | Balloonoccluded transarterial chemoembolization using a 1.8-french coaxial microballoon catheter for hcc: Tehnical and safety considerations[END_REF], [START_REF] Salem | Radioembolization with 90yttrium microspheres: A state-of-the-art brachyterapy treatment for primary and secondary liver malignancies-part 1: Technical and methodologic considerations[END_REF] .

Recently, magnetic navigation using a clinical magnetic resonance imaging (MRI) scanner has been introduced as a novel technique to improve the above mentioned loco regional procedures. It consists to serially steer single magnetic agents (delivered through the catheter tip) to the disease site to reach a given drug dose or radioactive isotope concentration (Fig. 1). As example, radioactive magnetic biodegradable nanoparticles that incorporate both magnetite and the β-emitter 90 Y have demonstrated efficient tumor targeting by using the superconducting magnet provided by MRI scanner [START_REF] Hafeli | Effective targeting of magnetic radioactive y-microspheres to tumor cells by an externally applied magnetic field. preliminary in-vitro and in-vivo resuls[END_REF]. In the same way, Martel et al. have proposed the use of an upgraded clinical MRI scanner to enable magnetic resonance navigation (MRN) of therapeutic agents to realize direct tumor targeting [START_REF] Mathieu | Method of propulsion of a ferromagnetic core in the cardiovascular system through magnetic gradients generated by an MRI system[END_REF]- [START_REF] Bigot | Magnetic resonance navigation of a bead inside a three-bifurcation pmma phantom using an imaging gradient coil insert[END_REF]. In their experiments, the authors have used FeCo magnetic particles to carry doxorubicin (cytotoxic drug) coated with a biodegradable polymer (PLGA) [START_REF] Pouponneau | Coencapsulation of magnetic nanoparticles and doxorubicin into biodegradable microcarriers for deep tissue targeting by vascular MRI navigation[END_REF]. These preliminary trials point out the need to know precisely the number, the size, the shape and the steering properties of the therapeutic agents to be injected with respect to the developmental stage of the tumor [START_REF] Mellal | Optimal structure of particles-based superparamagnetic microrobots: application to mri guided targeted drug therapy[END_REF]. In this study, to maximize the effect of the treatment and to minimize adverse effects on the patient, mathematical models of the tumor growth are considered to find the number of magnetic agents to be administered. Indeed, an ultimate goal in the clinical setting is to use models to help design optimal therapeutic strategies. Actually, mathematical modeling offers interesting tools that could give insights into a better understanding and control of these open clinical problems. In particular, modeling of cancer behavior is an active research field for biologists, mathematicians and engineers. Different approaches are used in the mathematical modeling of cancer and its control [START_REF] Ribba | A review of mixedeffects models of tumor growth and effects of anticancer drug treatment used in population analysis[END_REF]- [START_REF] Kovács | Model-based angiogenic inhibition of tumor growth using modern robust control method[END_REF]. For instance, some researchers investigated the tumor growth model by using cellular automata which can include very specific characteristics of the tumor, patient and drug effectively in the model [START_REF] Kansal | Simulated brain tumor growth dynamics using a three-dimensional cellular automaton[END_REF], [START_REF] Gerlee | An evolutionary hybrid cellular automaton model of solid tumour growth[END_REF]. Most of mathematical formulations are mainly constructed using ordinary differential equations (ODE) to exhibit the cancer dynamics and their response to the therapeutic agents [15]- [START_REF] De Pillis | The dynamics of an optimally controlled tumor model: A case study[END_REF]. Therefore, the modeling of tumor treatment is realized for chemotherapy [START_REF] Sanga | Mathematical modeling of cancer progression and response to chemotherapy[END_REF], [START_REF] De Pillis | Chemotherapy for tumors: An analysis of the dynamics and a study of quadratic and linear optimal controls[END_REF], immunotherapy [START_REF] De Pillis | The dynamics of an optimally controlled tumor model: A case study[END_REF], anti-angiogenic therapy [15]- [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF], radiotherapy [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF] as well as a combination of the above [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF]- [START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF]. For example, de Pillis and Radunskaya [START_REF] De Pillis | The dynamics of an optimally controlled tumor model: A case study[END_REF] set out an ODE system to depict the dynamics of tumor growth by means of the populations of tumor, normal and immune cells. To the best of the author's knowledge, no studies have investigated the mathematical modeling of locoregional procedures. This is mainly due to modeling difficulties. First, as the tumor behavior (avascular, vascular, and metastatic stages) and medication level are highly patient-dependent, different therapeutical materials may be considered as navigable agents, e.g. drug eluting microspheres, biologically active agents, chemical mediators of cell function, viral vectors or genetic material. Second, a large dispersion of the therapeutic agents occurs along the feeding vessels to the tumor(s) that endup in the systemic flow. It renders difficult the estimation of drug concentration deposited within the tumor. Third, the tolerability of conventional locoregional therapies seems to be affected by the type of regimen and the frequency of the treatment. The model should take into account constraints on drug delivery rates and state of the patient by solving an optimization problem at regular time intervals. This work aims to fill this gap by extending conventional therapy models for the transcatheter tumor therapy using MRN procedure for hepatocellular carcinoma (HCC). This paper is organized as follows: first, Section 2 presents the mathematical model framework for cancer therapy. Then, Section 3 proposes an optimal direct targeting therapy optimization problem. Section 4 demonstrates the applicability of the design framework in the context of MRN-assisted procedures. Finally, Section 5 discusses the main limitations when translating the theoretical results to practical experiments.

II. MATHEMATICAL MODELING OF CANCER A. Background

Cancer is a group of diseases involving basically abnormal cell growth (neoplasia). With over 100 different known forms that could affect all human organs, cancer is not a unique pathology. Hence, cancers are classified by the type of malignant cell, such as blastoma, carcinoma, sarcoma, lymphoma and leukemia. Therefore, many treatments against cancer exist, including surgery, chemotherapy, radiation therapy, or palliative care. The chosen treatment(s) depends on the type, location, and grade of the cancer as well as the person's states.

Thus, the understanding of the cancer evolution is a main issue in its fight. To this aim, it is mandatory to investigate the neoplasia process. Actually, neoplasia commonly forms a mass that is referred as tumor. Basically, a tumor evolution could be described in three stages, as illustrated in Fig. 2. First, a tumor appears and starts to grow to obtain its nutriments from its immediate environment. At this step, known as avascular growth, it is usually difficult to detect it through medical screening. The tumor could then continue to expand by seeking additional resources that allows it entering to the second stage: the vascular growth. At this step, the body provides new nutriments, and endothelial cells migrate to the tumor to form new microvascular networks. This process, referred as angiogenesis, allows the tumor to continue its development. In this vascular stage, symptoms may appear and the tumor becomes detectable on medical images. Finally, the tumor could evolve to the third step: the metastasis stage. In this step, the tumor cells are spread to other organs, and it becomes more difficult to treat them.

B. Tumor Growth Models

Different cancer models have been proposed [START_REF] Ribba | A review of mixedeffects models of tumor growth and effects of anticancer drug treatment used in population analysis[END_REF]- [START_REF] Kovács | Model-based angiogenic inhibition of tumor growth using modern robust control method[END_REF].

One key aspect in the cancer evolution is the angiogenesis process. Commonly, in this context, the Hahnfeldt's tumor growth model [15] is considered, as it could be applied for different cancer therapies [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF]- [START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF]. Specifically, the underlying diffusion of tumor and endothelial cells that stimulate and inhibit angionesis are incorporated into a system of ordinary differential equation (ODE). Different extensions of the original Hahnfeldt's model have been proposed [START_REF] Onofrio | Tumour eradication by antiangiogenic therapy: analysis and extensions of the model by hahnfeldt et al[END_REF]- [START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF]. Let c denotes the volume of cancer cells, and e the volume of endothelial cells that supplies the tumor with nutriments. The tumor angiogenesis evolution could be then described by the following ODE system:

ċ(t) = -λ c c(t) log c(t) e(t) (1)
ė(t) = bc(t) -dc(t) 2/3 e(t) (2)
with λ c the tumor growth rate; b the vascular endothelial cells birth rate; and d the vascular endothelial cells death rate. Especially, the tumor follows a Gompertzian law [15]: its growth saturates at a maximal volume, and we get the following steady-state:

c ∞ = e ∞ = b d 3/2 (3)
The endothelial cell birth (b) and death (d) rate depend mainly on the type of tumor and the patient. Obviously, the system model states that the tumor cannot increase over c ∞ volume, and then does not take into account the evolution to the metastasis stage. As the system model ( 1)-( 2) is mainly devoted to the vascular stage, to take into account the other stages the endothelial birth rate b could be considered as piece-wise continuous. corresponding values are reported in Table I. Fig. 3a shows the evolution of the tumor and endothelial volume growth for each stage. First, we consider that a malignant neoplasm appears with an initial volume c 0 = 20 mm 3 . We assume that nutriments sought by the tumor from its immediate environment could be equivalent to an initial volume of endothelial cells of e 0 = 20 mm 3 . Hence, in the avascular stage the endothelial state simulates the resources absorbed by the tumor (see Fig. 3b). The tumor volume increases with a Gompertzian curve, and after 100 days reaches the equilibrium value: c avasc.,∞ = e avasc.,∞ = 625 mm 3 , predicted by (3). After 200 days, endothelial cell migrates to the tumor to form a microvasculature networks around the tumor. The avascular equilibrium values is used as the initial condition of the vascular stage (cf. Fig. 3c). Once again, an equilibrium is reached: c vasc.,∞ = e vasc.,∞ = 17 346.5 mm 3 . If no treatment is administered or it fails to block the growth of the tumor, the metastasis stage could occur (cf. Fig. 3d).

C. Cancer Therapy Modeling

This section addresses the mathematical modeling of targeted therapy for optimal administration of drug-loaded magnetic nanoparticles, namely therapeutic vectors. When a therapeutic vector is delivered to a host, two different types of processes are involved: pharmacokinetics (PK) and pharmacodynamics (PD). The PK model characterizes what the body does with the agent, and is commonly modeled through the concentration of the therapeutic agents δ(t), in the sampled fluid (e.g. plasma or blood). The agent concentration expression is [15]:

δ(t) = t 0 u(s) exp (-λ δ (t -s)) ds ( 4
)
where u is the rate of administration of the therapeutic agent, and λ δ is the elimination rate. The PD describes the effects of the agent to the body, and is related to the concentration δ(t) provided by the PK model. The classic PD model is the logkill effect [15], [START_REF] De Pillis | The dynamics of an optimally controlled tumor model: A case study[END_REF]. The following presents how the PK and PD are related to the tumor growth, according to considered cancer therapy. 1) Angiogenesis Inhibition Treatment: Basically, for such treatment only the endothelial growth (2) will be antagonized by an anti-angiogenic agent, such as endostatin, whose plasma concentration is δ a (t), leading to its new formulation [15]:

ė(t) = bc(t) -dc(t) 2/3 e(t) -k a δ a (t)e(t) (5)
with k a the angiogenic drug killing parameter. Table III shows the PK parameters of the endostatin.

2) Chemotherapy Treatment: The literature provides different mathematical formulations of the tumor functioning and response to chemotherapy [START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF]- [START_REF] De Pillis | Chemotherapy for tumors: An analysis of the dynamics and a study of quadratic and linear optimal controls[END_REF]. Here, to model the effect of a cytotoxic chemotherapeutic agent of concentration δ x (t) the angiogenesis model is used [START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF], and is extended as:

ċ(t) = -λ c c(t) log c(t) e(t) -k cx c(t)δ x (t) (6)
ė(t) = bc(t) -dc(t) 2/3 e(t) -k ex e(t)δ x (t) (7)
with k cx and k ex the PD log-kill parameters. The case k ex = 0 means that the cytotoxic drug does not have any effects on the endothelial cells. Contrarily to angiogenesis inhibition model, no PK-PD parameters {k cx , k ex , λ δ } characterized from experimental data are available in the literature.

3) Radiotherapy Treatment: Classically, the so-called tumor linear-quadratic (LQ) model [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF], [START_REF] Ledzewicz | Tumor development under combination treatments with anti-angiogenic therapies[END_REF], [START_REF] Fowler | The linear-quadratic formula and progress in fractionated radiotherapy[END_REF] is used to characterize the effects of ionizing radiation on the tumor and endothelial cells. Thus, the damage of radiation on a cell x is modeled in the following form [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF], [START_REF] Ledzewicz | Tumor development under combination treatments with anti-angiogenic therapies[END_REF], [START_REF] Fowler | The linear-quadratic formula and progress in fractionated radiotherapy[END_REF]:

-x(t) α x + β x t 0 (u(s) exp (-µ(t -s)) ds) u(t) (8)
where u(t) represents the fractionation schedule of radiation; α x and β x are positive gains related to the tumor-LQ parameters; and µ is the cell repair rate. As the radiation damage the tumor and the endothelial cells, the model for radiotherapy could be expressed with the following ODE set [START_REF] Ledzewicz | Tumor development under combination treatments with anti-angiogenic therapies[END_REF]:

ċ = -c λ c log c e + (α c + β c r) u (9)
ė = bc -e dc 2/3 + (α e + β e r) u [START_REF] Mathieu | Method of propulsion of a ferromagnetic core in the cardiovascular system through magnetic gradients generated by an MRI system[END_REF] with r(t) the concentration of the radiation.

III. OPTIMAL DIRECT TARGETING THERAPY

In traditional treatment regimes the amount of administered therapeutic agents is very important for patient's survival. Actually, the therapy does not treat only the tumor, it often kills some healthy tissues or causes them serious damage. Hence, the dosage of the therapy must be carefully adjusted to minimize side-effects, while maximizing the capability to destroy the tumor. To address this issue, control theory is basically employed [START_REF] Swan | Role of optimal control theory in cancer chemotherapy[END_REF]. In particular, optimal control techniques have been applied for antiangiogenic therapy [START_REF] Onofrio | Tumour eradication by antiangiogenic therapy: analysis and extensions of the model by hahnfeldt et al[END_REF], [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF], chemotherapy [START_REF] De Pillis | Chemotherapy for tumors: An analysis of the dynamics and a study of quadratic and linear optimal controls[END_REF], [START_REF] De Pillis | The dynamics of an optimally controlled tumor model: A case study[END_REF], [START_REF] Batmani | Optimal chemotherapy in cancer treatment: state dependent riccati equation control and extended kalman filter[END_REF] or radiotherapy [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF] to define the optimal treatment and drug dose. This section presents an adaptation of conventional therapies [15]- [START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF] to the case of magnetically controlled drug-loaded nanoparticles.

A. Optimal Control Design

Optimal control deals with the problem of finding a regulator such that a certain optimality criterion is achieved [START_REF] Zhou | Robust and optimal control[END_REF], [START_REF] Mracek | Control designs for the nonlinear benchmark problem via the state-dependent riccati equation method[END_REF]. Given cancer models presented in Section II, the issue is to find an optimal control input u to decrease the tumor size while minimizing total drug administered. However, the different therapy models exhibit a highly nonlinear and complex nature. Thus, different optimization techniques have been used for stabilizing tumor growth while minimizing the administered therapeutic agent [START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF], [START_REF] De Pillis | Chemotherapy for tumors: An analysis of the dynamics and a study of quadratic and linear optimal controls[END_REF], [START_REF] Kovács | Model-based angiogenic inhibition of tumor growth using modern robust control method[END_REF], [START_REF] Batmani | Optimal chemotherapy in cancer treatment: state dependent riccati equation control and extended kalman filter[END_REF]. One of these approaches consider the state-dependent Riccati equations (SDRE) framework [START_REF] Batmani | Optimal chemotherapy in cancer treatment: state dependent riccati equation control and extended kalman filter[END_REF], [START_REF] Mracek | Control designs for the nonlinear benchmark problem via the state-dependent riccati equation method[END_REF], [START_REF] ¸imen | Systematic and effective design of nonlinear feedback controllers via the state-dependent riccati equation (sdre) method[END_REF]. Actually, SDRE has emerged as an interesting strategy for direct synthesis of nonlinear controllers [START_REF] Mracek | Control designs for the nonlinear benchmark problem via the state-dependent riccati equation method[END_REF], [START_REF] ¸imen | Systematic and effective design of nonlinear feedback controllers via the state-dependent riccati equation (sdre) method[END_REF]. This optimal control methodology has been applied with the de Pillis et al. cancer model [START_REF] De Pillis | The dynamics of an optimally controlled tumor model: A case study[END_REF] in [START_REF] Batmani | Optimal chemotherapy in cancer treatment: state dependent riccati equation control and extended kalman filter[END_REF]. However, to the authors' knowledge, this optimal framework has not been applied to the Hahnfeldt's tumor growth formulation [15] or its extension [START_REF] Onofrio | Tumour eradication by antiangiogenic therapy: analysis and extensions of the model by hahnfeldt et al[END_REF]- [START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF].

1) SDRE Optimal Regulation: Let us consider a system represented by the following pseudo-linear equations in statedepend coefficient (SDC) form:

ẋ = A(x)x + B(x)u (11)
where x ∈ R n is the time-dependent state vector, with n the number of states; u ∈ R m is the (drug) control input (m = 1 for monotherapy); A(x) ∈ R n×n and B(x) ∈ R n×m are statedependent matrices. In infinite horizon, the continuous-time SDRE optimization problem is then to minimize the following quadratic cost functional [START_REF] Zhou | Robust and optimal control[END_REF], [START_REF] Mracek | Control designs for the nonlinear benchmark problem via the state-dependent riccati equation method[END_REF]:

J = x T Qx + u T Ru dt (12)
where the weights (design parameters) are state-dependent, such that Q ∈ R n×n is positive-definite and R ∈ R m×m .

Thus, the optimal control law that minimizes this criterion is given by: u = -R -1 (x)B (x)P(x) x

with P(x) is the unique, symmetric, positive-definite solution of the algebraic state-dependent Riccati equation (SDRE) [START_REF] Mracek | Control designs for the nonlinear benchmark problem via the state-dependent riccati equation method[END_REF], [START_REF] ¸imen | Systematic and effective design of nonlinear feedback controllers via the state-dependent riccati equation (sdre) method[END_REF]:

PA + A P -PBR -1 B P + Q = 0 (14)
Let us notice that the SDRE control design is similar to the well known linear quadratic regulator (LQR), where each of the above matrices are state-dependent.

2) Pseudolinearization: The nonlinear ordinary differential equation (ODE) of the considered cancer therapy could be expressed in the following general nonlinear dynamical system:

ẋ(t) = f (x) + g(x)u(t) (15)
To apply the above SDRE methodology, the so-called pseudolinearized form (11) of the nonlinear system (15) has to be computed. More precisely, the state matrix is obtained from f (x) = A(x)x by mathematical factorization. This statedependent parametrization of the nonlinear system is not unique [START_REF] Mracek | Control designs for the nonlinear benchmark problem via the state-dependent riccati equation method[END_REF]. Among the alternatives, the chosen pseudolinearization should ensure point-wise controllability condition, where the so-called state-dependent controllability matrix has full rank, that is:

C = B(x) A(x)B(x) . . . A n-1 (x)B(x) (16) rank(C) = n (17)

B. Control of Angiogenesis Inhibitor

The objective here is to apply the SDRE optimal control strategy to reduce the growth of a tumor through the administration of appropriate dose of anti-angiogenic agents. As mentioned, the drug-free equilibrium is given by (3). Thus, the state vector x is designed to shift the equilibrium point to the origin, that is:

x = (x 1 = (c -c ∞ ); x 2 = (e -e ∞ ); x 3 = δ a ) (18)
where c ∞ = e ∞ = x ∞ are the drug-free volume steadystate (3). A state-dependent parametrization of the angiogenesis inhibition treatment could be then written as:

A(x) =   a 11 λ c x1 x2 a 13 a 21 a 22 a 23 0 0 -λ δ   , and B(x) =   0 0 1   (19)
with

a 11 = -λ c 1 + log x 1 + x ∞ x 2 + x ∞ 1 + x ∞ x 1 (20)
a 13 = 0 (21)
a 21 = b 1 + x ∞ x 1 - x ∞ x 1 d (x 1 + x ∞ ) 2/3 (22)
a 22 = -d (x 1 + x ∞ ) 2/3 -k a x 3 (23) a 23 = -k a x ∞ ( 24
)
This pseudolinearization is used to apply the SDRE control of the amount of anti-angiogenetic agent by considering R = 10 3 and Q = diag (50, 0, 0). This Q and R values are chosen to take into account the volume of the tumor, while preventing high doses. Especially, the drug intake u has to be saturated to an upper limit u max . detected mainly from this step, and the treatment is applied at the beginning of this stage (Fig. 3). Let us recall that our main motivation is to propose a direct targeted drug delivery using magnetic boluses. Therefore, our drug intake approach differs significantly than the original subcutaneous injection technique investigated by Hahnfeldt et al. [15]. The direct targeting could be modeled through the drug clearance rate λ δ . Obviously, by varying λ δ the amount of drug that is not filtered out of the blood and delivered to the tumor can be modeled. Similarly to the works achieved in [15]- [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF], endostatin antiangiogenic agents is considered here. The elimination rate of endostatin is varying from 1.7 day -1 for a subcutaneous injection [15] to 0.017 day -1 [START_REF] Shin | Targeted delivery of an antibodymutant human endostatin fusion protein results in enhanced antitumor efficacy[END_REF], which corresponds to the efficiency of the drug in targeted treatment. The different values of PK-PD models parameters of endostatin are summarized in the Table III based on the results given in [15], [START_REF] Shin | Targeted delivery of an antibodymutant human endostatin fusion protein results in enhanced antitumor efficacy[END_REF]. The Fig. 4a and4b show that for a clearance of λ δ = 0.017 day -1 the tumor and the endothelial volumes vanish in about 60 days. In contrast, for λ δ = 1.7 day -1 the volumes reach a steadystate of c ∞ = e ∞ = 476.5 mm 3 . However, the Fig. 4c illustrates that a continuous adminisof endostatin has to be realized. For instance, a reliable therapy procedure would be a session of 30 min every 7 days. Nevertheless, due to the endostatin limit u max , it is not possible to treat the tumor with such treatment sequence. Indeed, Fig. 5 shows the evolution of the tumor and the corresponding endostatin intake (without the upper bound). As one can see, only a low elimination rate λ δ enables a noticeable decrease of the tumor volume. These results exhibit that the angiogenesis inhibition could be significantly improved thanks to direct targeting.

C. Chemotherapeutic Agents Direct Delivery

The same reasoning as above is used here to design optimal chemotherapy. Especially, the chemotherapy treatment model is similar to the angiogenesis inhibition model. Hence, the same state vector [START_REF] Ledzewicz | Tumor development under combination treatments with anti-angiogenic therapies[END_REF] and state-dependent parametrization ( 19) could be considered, where only {a 13 , a 22 , a 23 } parameters are modified as follows:

a 13 = -k cx (x 1 + x ∞ ) ( 25
)
a 22 = -d (x 1 + x ∞ ) 2/3 -k ex x 3 (26) a 23 = -k ex x ∞ ( 27
)
SDRE method is then applied to control the amount of drug. The optimal controller parameters are defined here as: R = 10 4 and Q = diag (1, 0, 0).

Contrarily to the anti-angiogenetic model, there are not PK-PD parameters based on experimental data for the model presented in paragraph II-C2. Actually, most previous works [START_REF] Onofrio | Tumour eradication by antiangiogenic therapy: analysis and extensions of the model by hahnfeldt et al[END_REF], [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF], [START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF] have considered experiments only. In this work the PK-PD parameters reported in Table IV are proposed based on d'Onofrio et al. analysis [START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF] to fit the DEB-TACE procedure pharmacokinetics [START_REF] Varela | Chemoembolization of hepatocellular carcinoma with drug eluting beads: efficacy and doxorubicin pharmacokinetics[END_REF]. This choice allows to address direct targeting (low clearance) with a highly cytotoxic drug (e.g. doxorubicin) combined with embolization of the microvascular network feeding in nutriments the tumor (k ex > 0).

Param. Value

Units

λ δ 0.009 day -1 kcx 2 day -1 conc -1 kex 1 day -1 conc -1
A tumor at vascular stage is then treated thanks to the cytotoxic agent (e.g. doxorubicin), and the results are depicted in Fig. 6. Here, the sampling time is settled to 1 h. In particular, the therapy schedule is a session of 30 min every 7 days. To limit the amount of cytotoxic drug, the intake is saturated to u max = 2.5 conc/day, that is 0.0417 (mg/kg)/day per session. Hence, after 6 sessions of drug administration the tumor vanish under a volume of c ∞ = 12.67 mm 3 . demonstrating the efficiency of cytotoxic agents.

D. Direct Delivery of Radioisotope

Based on the same modeling, the optimal control of radiotherapy is designed on the same previous formulation. Hence, the same state vector [START_REF] Ledzewicz | Tumor development under combination treatments with anti-angiogenic therapies[END_REF] and the state matrix A(x) proposed in [START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF] could be considered, where only {a 13 , a 22 , a 23 } parameters are modified as follows: However, the state-dependent input matrix is now given by:

a 13 = 0 (28)
a 22 = -d (x 1 + x ∞ ) 2/3 (29
B(x) =   -(α c + β c r) (x 1 + x ∞ ) -(α e + β e r) (x 2 + x ∞ ) 1   (31)
The optimal amount of radioisotope is computed using the following SDRE controller parameters: Q = diag (1, 0, 1) and R = 10 3 . The tumor-LQ model is used to characterize the effects of radiation on the tumor and endothelial cells [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF], [START_REF] Ledzewicz | Tumor development under combination treatments with anti-angiogenic therapies[END_REF], [START_REF] Fowler | The linear-quadratic formula and progress in fractionated radiotherapy[END_REF]. For the sake of simplicity, similar numerical values of the tumor-LQ parameters proposed in [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF] are considered, and summarized in Table V. As previously, we consider a cancer at vascular stage. Fig. 6 illustrates the corresponding radiotherapy treatment. To limit the amount of radiation, the radioisotope dose is saturated to u max = 10 Gy. Once again, the radiotherapy schedule is a session of 30 min every 7 days. Thus, after 10 sessions the tumor vanish under a volume of c ∞ = 0.195 mm 3 . As expected, radiotherapy is able to cure efficiently the cancer but the required number of sessions is higher compared to chemotherapy.

IV. OPTIMAL DELIVERY OF THERAPEUTIC AGENTS

A. Structure standard transcatheter tumor delivery using MRN, the therapeutic agents are constituted by magnetic boluses of microagglomerations of iron-oxide nanoparticles and drug loadings (shown in Fig. 8). The optimal size and shape of the magnetic magnetic microcarriers have been determined in [START_REF] Mellal | Optimal structure of particles-based superparamagnetic microrobots: application to mri guided targeted drug therapy[END_REF] through simulations and experiments. Spherical boluses of 500 µm of diameter composed of SPIO particles (BioMag BM547, Bang Laboratories, Inc.) of 8 µm in size were chosen. To be controllable, the bolus has its optimal nondimensional magnetophoretic number C mt ratio greater than 1. The reader may refer to [START_REF] Mellal | Optimal structure of particles-based superparamagnetic microrobots: application to mri guided targeted drug therapy[END_REF] for optimal size selection.

The effective drug load V δ of a single bolus is a fraction of the total volume defined as V δ = (1 -τ m )V with τ m the so-called magnetization rate. When injected in the vascular phantom placed in the tunnel of the MRI scanner, the magnetic bolus will take the form of a prolate sphere due to the strong dipole-to-dipole interactions between the nanoparticles. Due to the B 0 field of the scanner, the equivalent magnetization M s of the aggregates is aligned with the B 0 field. We consider here that the magnetic material (Vm) is saturated to its maximal value M s = 56 emu/g under a magnetic field B 0 = 3.0 Tesla. The injected magnetic volume τ m V is computed as a tradeoff between the encapsulated drug volume and the necessary steering force. Finally, the number of magnetic boluses is estimated from

N b = u m V δ ρ δ ( 32
)
where u is the rate of administration of the therapeutic agent; m is the patient weight (considered as m = 60 kg hereafter); and ρ δ is the agent load density.

B. Endostatin Angiogenesis Inhibitor

The objective is to determine the amount N b of required boluses to reduce the tumor size by administrating an antiangiogenic agent. In paragraph III-B endostatin agent is used as angiogenesis inhibitor. Endostatin is a naturally occurring 20 kDa polypeptide, identified as a C-terminal fragment of collagen (type 18). Commonly, endostatin agent is available in solution concentration of 1 mg ml -1 . Hence, considering the upper limit u max given in Table III, the Table VI presents the maximum of required number of magnetic boluses that have been to be conveyed continuously. Even in the low clearance case with λ δ = 0.017 day -1 where the endostatin administration rate is lower (u max = 2 conc day -1 ), N b = 4661 magnetic boluses have to be injected continuously to administrate the daily dose. Similarly, if the considered scheduling of one session every week is applied (cf. Fig. 5), N b = 1.398 × 10 4 boluses have to be administered at the first session.

C. Cytotoxic Drug (doxorubicin)

Chemotherapy is achieved by direct targeting of highly cytotoxic drug, such as doxorubicin. Doxorubicin is an anthracycline tumor antibiotic that works by intercalating DNA. This drug has already shown its efficiency in DEB-TACE procedure since a decade [START_REF] Lewis | DC bead: in vitro characterization of a drug-delivery device for transarterial chemoembolization[END_REF], [START_REF] Varela | Chemoembolization of hepatocellular carcinoma with drug eluting beads: efficacy and doxorubicin pharmacokinetics[END_REF]. The drug eluting beads (DEB) have been easily loaded with doxorubicin to a recommended level of 25 mg ml -1 , with a maximum loading of 45 mg ml -1 [START_REF] Lewis | DC bead: in vitro characterization of a drug-delivery device for transarterial chemoembolization[END_REF]. Here, we assume that the doxorubicin PK-PD model paramis close to the cytotoxic agent given in Table IV. For a patient of average weight m = 60 kg a maximum dose of 120 mg has to be conveyed per session (cf. Fig. 6). Therefore, a maximum of N b = 1036 magnetic boluses have to be administered.

D. Radioisotope (Yttrium)

The radiotherapy consists to administrate microspheres containing a radioisotope. Usually, yttrium-90 ( 90 Y) microspheres are administered in directed cancer radiotherapy [START_REF] Kennedy | Recommendations for radioembolization of hepatic malignancies using yttrium-90 microsphere brachytherapy: a consensus panel report from the radioembolization brachytherapy oncology consortium[END_REF]. 90 Y is a high-energy beta-emitting isotope with no primary gamma emission and a half-life of 2.67 days. For instance, TheraSphere (MDS Nordion, Canada) are biocompatible and nonabsorbable glass microspheres loaded with 90 Y. The radius of TheraSphere beads is about r p = 12.5 ± 2.5µm, and then one magnetic bolus could embed up to 3147 radioisotope microspheres. Each 90 Y bead contains an activity of 2500 Bq, leading to 7.8675 × 10 -3 GBq per bolus. Assuming an uniform distribution, the radiation dose intake is computed from [START_REF] Kennedy | Recommendations for radioembolization of hepatic malignancies using yttrium-90 microsphere brachytherapy: a consensus panel report from the radioembolization brachytherapy oncology consortium[END_REF] u = 49.67A m l with A the radioisotope activity (in GBq), and m l the mass of the liver (in kg). Therefore, the number of magnetic boluses that carries 90 Y beads is obtained from:

N b = u • m l 49.67 × 3147 × 2.5 × 10 -6 (33)
We consider that the radioisotope tumor-LQ model parameters given in Table V corresponds to TheraSphere. For a typical patient with a liver mass of m l = 2 kg, a maximum of N b = 52 boluses have to be administered in one session to carry the upper limit of 10 Gy (cf. Fig. 7).

V. DISCUSSION

Locoregional therapies are considered as the best treatments in patients with unresectable hepatocellular carcinoma (HCC). The various mathematical models of cancer therapies for HCC that have been evaluated so far describe some interactions between tumor cells and blood vessels. However, the internal cell behavior is not explicitly expressed, whereas further sources of complexity arise from its internal cell-tocell cooperative and competitive interactions [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF]. Especially, due to the use of ODE system the model also does not include the spatial characteristic of the tumor growth along the surrounding tissue. Particularly, this work considers medically reasonable parameters set given in Table I-to-V taken from the literature. Actually, only few data are available to estimate the parameter sets, especially in the case of human models [START_REF] Onofrio | Tumour eradication by antiangiogenic therapy: analysis and extensions of the model by hahnfeldt et al[END_REF], [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF], [START_REF] De Pillis | The dynamics of an optimally controlled tumor model: A case study[END_REF]. Mathematically modeling the dynamics of a tumor allows us to better understand tumor activity by considering critical biological influences, time delays, human intervention and tumor adaptation. Our model parameters are computed on the basis of experimental data given in [15]. The environmental factors such as habits, smoking, alcohol consumption, and infections, have a profound influence on cancer development. As example, it is possible to regulate the nuclear transcription factor, NF-kB, by controlling the vascular endothelial cells birth rate, b, expressed in the model [START_REF] Anand | Cancer is a preventable disease that requires major lifestyle changes[END_REF]. Furthermore, the control of the drug input u offers the possibility to model the drug overdose and/or toxicity.

In the proposed direct targeting therapy, the drug input is saturated to u max values which should be related to the upper bound of the drug concentration u max . For instance, the most commonly used dosage of doxorubicin is to not exceed δ max = 2 mg/kg to have a very low risk of overdose. Finally, the theoretical simulation results point out some benefits when considering MRN-assisted locoregional therapy. At an earlierstage disease, TACE and TARE treatments are commonly used. Actual TACE treatments with doxorubicin at baseline (cTACE) or doxorubicin-loaded eluting beads (DEB-TACE) deliver a dose of drug over a longer period of time (3 months) into the tumor. The schedule of treatment is settled to 3 months, 6 months, and then every 6 months thereafter. The size of DEB is chosen according to the particular study, usually with smaller particles (100-to-300 µm) being selected first, followed by larger particles (500-to-700 µm). When performing MRN-assisted procedure using embolic magnetic boluses releasing doxorubicin, similar benefits (drug rate and frequency of treatment) are found but the technique may prevent bile duct injury by avoiding overembolization that will relentlessly lead to nontarget embolization. In the case of embolic radioembolization in TARE treatment, glass or resine microspheres of 20-to-30 µm in size are used. A high number of aggregated spheres (glass: 1.2 million microspheres/3 GBq; resin: 40-to-60 millions of microspheres/3 GBq) are required leading to deposition in the liver as a large number of discrete clusters. The use of MRN could significantly improve the technique by reducing the number and the timeless of sequences, and most importantly, by providing heterogeneous radiation dose distribution patterns. At an advanced-stage disease, Sorafenib, a multikinase inhibitor with antiangiogenic properties of Sorafenib, in combination with locoregional therapy if of increased interest to specialists who care for patients with unresectable hepatocellular carcinoma [START_REF] Ledzewicz | Tumor development under combination treatments with anti-angiogenic therapies[END_REF]. To demonstrate the feasibility of MRN-assisted procedure, Bigot et al. [START_REF] Bigot | Magnetic resonance navigation of a bead inside a three-bifurcation pmma phantom using an imaging gradient coil insert[END_REF] have demonstrated that the propulsion of N b = 1530 magnetic boluses inside a multi-bifurcation phantom required 33 min of duration. These experiments are in agreement with our simulation results concerning TACE and TARE treatments since only 1036 and 52 magnetic boluses are required.

VI. CONCLUSION

The paper described an optimal design strategy for innovative transcatheter arterial magnetic bolus delivery using MRN procedure. To maximize the effect of the treatment and minimize adverse effects on the patient, a mathematical model has been developed to find the number of magnetic boluses with respect to the growth of a tumor and the carried therapeutic agent. The application of such mathematical design framework is the optimization of anti-cancer treatment modalities by its optimal use with locoregional tumor therapies. Indeed, the preliminary results demonstrate that MRN may improve the conventional TACE and TARE treatments by increasing the local drug concentration, reducing the number and duration of treatment sequences, and decreasing the side-effects related to toxicity in healthy tissues. Undergoing work considers the locoregional magnetic bolus delivery of various sizes of aggregates in small animals.

  i) The tip of the catheter positioned in the upstream of the tumor ii) the embolizations block the tumor feed iii) the eluting beads release the therapeutic agent Tumor Liver iv) some beads end up in the systemic flow

Fig. 2 .

 2 Fig. 2. Representation of the liver cancer stages: (i)avascular stage: the tumor is small and found in one part of the liver without connection with the blood vessels or lymph nodes; (ii) vascular stage: tumor has grown into the blood vessels of the liver; (iii) metastases stage: tumor has spread into other organs.

Fig. 3 .

 3 Fig. 3. (a) Tumor growth rate changes at different stages of the liver cancer. Evolution of the tumor and endothelial volumes without drug administered: (b) vascular growth from c 0 = 20 mm 3 and e 0 = 20 mm 3 ; (c) vascular step from the avascular equilibrium point cavasc.,∞ = eavasc.,∞ = 625 mm 3 ; and (d) metastasis stage simulation from the vascular equilibrium point cvasc.,∞ = evasc.,∞ = 17 346.5 mm 3 .

Fig. 4 .

 4 Fig. 4. Continuous infusion angiogenesis inhibition therapies: decreasing of (a) tumor volume and (b) endothelial over time with (c) optimal input drug.

  This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNB.2016.2535380

Fig. 5 .

 5 Fig. 5. Angiogenesis inhibition therapies with a session of 30 min every 7 days: (a) evolution of the tumor volume with (b) an optimal input drug.

) a 23 Fig. 6 .

 236 Fig. 6. Cytotoxic agent direct delivery: decreasing of (a) tumor volume and endothelial over time with optimal input drug (b).

Fig. 7 .

 7 Fig. 7. Radioisotope direct delivery: decreasing of (a) tumor volume and endothelial over time with optimal input drug (b).

Fig. 8 .

 8 Fig. 8. Magnetic bolus model embedding aggregates of iron-oxides microparticles. (a) General design and (b) oil-based ferrofluid bolus.