(38) Hwang, J.; Min, B.; Lee, J. S.; Keem, K.; Cho, K.; Sung, M. Y.; Lee, M. S.; Kim, S. Al2O3 Nanotubes Fabricated by Wet Etching of ZnO/Al2O3 Core/Shell Nanofibers. *Adv. Mater.* 2004, 16, 422-425.

(39) Yin, Y.; Lu, Y.; Sun, Y.; Xia, Y. Silver Nanowires Can Be Directly Coated with Amorphous Silica To Generate Well-Controlled Coaxial Nanocables of Silver/Silica. *Nano Lett.* 2002, 2, 427-430.

(40) Kong, X. Y.; Ding, Y.; Wang, Z. L. Metal-Semiconductor Zn-ZnO Core-Shell Nanobelts and Nanotubes. *J. Phys. Chem. B* 2004, 108, 570-574.

(41) Cao, H. Q.; Xu, Z.; Sang, H.; Sheng, D.; Tie, C. Y. Template Synthesis and Magnetic Behavior of an Array of Cobalt Nanowires Encapsulated in Polyaniline Nanotubes. *Adv. Mater.* 2001, 13, 121-123.

(42) Hu, J. Q.; Meng, X. M.; Jiang, Y.; Lee, C. S.; Lee, S. T. Fabrication of Germanium-Filled Silica Nanotubes and Aligned Silica Nanofibers. *Adv. Mater.* 2003, 15, 70-73.

(43) Moran, J. L.; Wheat, P. M.; Posner, J. D. Locomotion of electrocatalytic nanomotors due to reaction induced charge autoelectrophoresis. *Phys. Rev. Lett.* E 2010, 81, 065302.

(44) Moran, J. L.; Posner, J. D. Electrokinetic locomotion due to reaction-induced charge auto-electrophoresis. *J. Fluid Mech.* 2011, 680, 31-66.

(45) Wang, W.; Chiang, T. Y.; Velegol, D.; Mallouk, T. E. Understanding the efficiency of autonomous nano- and microscale motors. *J. Am. Chem. Soc.* 2013, 135, 10557-10565.

(46) Jang, B.; Chen, X.-Z.; Siegfried, R.; Montero Moreno, J. M.; Özkale, B.; Nielsch, K.; Nelson, B. J.; Pané, S. Silicon-supported aluminum oxide membranes with ultrahigh aspect ratio nanopores. *RSC Adv.* 2015, 5, 94283-94289.
(47) Hill, J. J.; Haller, K.; Gelfand, B.; Ziegler, K. J. Eliminating capillary coalescence of nanowire arrays with applied electric fields. *ACS Appl. Mater. Interfaces* 2010, 2, 1992-1998.

(48) Marine, N. A.; Wheat, P. M.; Ault, J.; Posner, J. D. Diffusive behaviors of circle-swimming motors. *Phys. Rev. Lett.* 2013, 100, 240701.

(49) Lee, T. C.; Alarcon-Correa, M.; Miksch, C.; Hahn, K.; Gibbs, J. G.; Fischer, P. Self-propelling nanomotors in the presence of strong Brownian forces. *Nano Lett.* 2014, 14, 2407-2412.

(50) Demirok, U. K.; Laocharoensuk, R.; Manesh, K. M.; Wang, J. Ultrafast catalytic alloy nanomotors. *Angew. Chem. Int. Ed. Engl.* 2008, 47, 9349-9351.

(51) Laocharoensuk, R.; Burdick, J.; Wang, J. Carbon-Nanotube-Induced Acceleration of Catalytic Nanomotors. *ACS Nano* 2008, 2, 1069-1075.

(52) Teeffelen, S. V.; Lowen, H. Dynamics of a Brownian circle swimmer. *Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys.* 2008, 78, 020101.

(53) Ebbens, S.; Jones, R. A.; Ryan, A. J.; Golestanian, R.; Howse, J. R. Self-assembled autonomous runners and tumblers. *Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys.* 2010, 82, 015304.

(54) Charsooghi, M. A.; Akhlaghi, E. A.; Tavaddod, S.; Khalesifard, H. R. A MATLAB program to calculate translational and rotational diffusion coefficients of a single particle. *Comput. Phys. Commun.* 2011, 182, 400-408.

(55) McKee, D. W. Catalytic decomposition of hydrogen peroxide by metals and alloys of the platinum group. *J. Catal.* 1969, 14, 355-364.
## ASSOCIATED CONTENT

### Supporting Information.

Experimental procedures, numerical modelling procedures, supplementary figures and supplementary movies are included.

This material is available free of charge via the Internet at http://pubs.acs.org.

## AUTHOR INFORMATION

### Corresponding Author

*E-mail: [bjang@ethz.ch](mailto:bjang@ethz.ch), [vidalp@ethz.ch](mailto:vidalp@ethz.ch)

### Author Contributions

S.P. and B.J.N. initiated the project. S.P. and B.J. designed the fabrication experiments. B.J. and S.W. fabricated the magnetic nanoswimmers and analyzed their motion. W.W. and B.J. conducted the simulation experiments. B.J., W.W., A.P., X.C., B.J.N., D.F. and A.F. performed the analysis of the swimmers and provided theoretical discussion. B.J., W.W., A.P., X.C., S.P., and B.J.N. supervised the work and gave critical input. All authors contributed to discussions.

### Notes

The authors declare no competing financial interests.
ACKNOWLEDGMENT

Funding support from the European Community’s Seventh Framework Programme (FP7/2007-
2013) under grant agreement 296679 (MANAQA) is acknowledged. S.P. acknowledges financial
support by the European Research Council Starting Grant “Magnetoelectric Chemonanomotorics
for Chemical and Biomedical Applications (ELECTROCHEMBOTS)”, by the ERC grant
agreement n. 336456. W.W. is grateful for the financial support from National Natural Science
Foundation of China (Grant No. 11402069) and the city government of Shenzhen (Grant No.
KQCX20140521144102503). We especially thank Prof. Mahmut Selman Sakar from the
Mechanical Engineering (School of Engineering, EPFL), and Carlos C.J. Alcantara and Ayoung
Hong from the Multi-Scale Robotics Lab (ETH Zürich) for constructive discussions.
# APPENDIX D

## MODELING OF OPTIMAL TARGETED THERAPIES USING DRUG-LOADED MAGNETIC NANOPARTICLES FOR THE LIVER CANCER (ACL12)
Modeling of Optimal Targeted Therapies using
Drug-Loaded Magnetic Nanoparticles for the Liver
Cancer

Lyès Mellal, David Folio, Karim Belharet, Antoine Ferreira

► To cite this version:

Lyès Mellal, David Folio, Karim Belharet, Antoine Ferreira. Modeling of Optimal Targeted Therapies using Drug-Loaded Magnetic Nanoparticles for the Liver Cancer. IEEE Transactions on NanoBioscience, Institute of Electrical and Electronics Engineers, 2016, pp (99), pp.10.10.1109/TNB.2016.2535380. hal-01305959

HAL Id: hal-01305959

https://hal.archives-ouvertes.fr/hal-01305959

Submitted on 22 Apr 2016

**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.
# Modeling of Optimal Targeted Therapies using Drug-Loaded Magnetic Nanoparticles for the Liver Cancer

Lyès Mellal, David Folio, Karim Belharet and Antoine Ferreira

**Abstract**—To enhance locoregional therapies for liver cancer treatment, we propose in this study a mathematical model to optimize the transcatheter arterial delivery of therapeutical agents. To maximize the effect of the treatment and minimize adverse effects on the patient, different mathematical models of the tumor growth are considered in this study to find the optimal number of the therapeutic drug-loaded magnetic nanoparticles to be administered. Three types of therapy models are considered, e.g. angiogenesis inhibition therapy, chemotherapy and radiotherapy. We use state-dependent Riccati equations (SDRE) as an optimal control methodology framework to the Hahnfeldt's tumor growth formulation. Based on this, design optimal rules are derived for each therapy to reduce the growth of a tumor through the administration of appropriate dose of anti-angiogenic, radio- and chemo-therapeutic agents. Simulation results demonstrate the validity of the proposed optimal delivery approach, leading to reduced intervention time, low drug administration rates and optimal targeted delivery.

**Index Terms**—Magnetic Resonance Navigation; Tumor Growth model; Optimal Tumor Control; Optimal drug delivery; Locoregional therapies.

## 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 [1]. Due to a lack of donors, it is even more rarely possible to perform a liver transplantation or to consider surgical procedures [2], [3]. Therefore, transarterial locoregional interventions, such as transcatheter arterial chemoembolization (TACE) or radioembolization (TARE), are procedures that provide a significant survival benefit [2]–[4]. 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 [5], [6]. 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.

<figure><img src="image_43.png" /><figcaption>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.</figcaption></figure>

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 [7], [8].

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 ⁹⁰Y have demonstrated efficient tumor targeting by using the superconducting magnet provided by MRI scanner [9]. 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 [10]–[12]. In their experiments, the authors have used FeCo magnetic particles to carry doxorubicin (cytotoxic drug) coated with a biodegradable polymer (PLGA) [11]. 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 [13]. 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 [14]–[23]. 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 [24], [25]. 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]–[22]. Therefore, the modeling of tumor treatment is realized for chemotherapy [20], [21], immunotherapy [22], anti-angiogenic therapy [15]–[17], radiotherapy [17] as well as a combination of the above [17]–[19]. For example, de Pillis and Radunskaya [22] 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 therapeutic 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 end up 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

<figure><img src="image_44.png" /><figcaption>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.</figcaption></figure>

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 *vascular 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 [14]–[23]. 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 [17]–[19]. Specifically, the underlying diffusion of tumor and endothelial cells that stimulate and inhibit angiogenesis are incorporated into a system of ordinary differential equation (ODE). Different extensions of the original Hahnfeldt’s model have been proposed [16]–[19].
<figure><img src="image_45.png" /><figcaption>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 \text{ mm}^3$ and $e_0 = 20 \text{ mm}^3$; (c) vascular step from the avascular equilibrium point $c_{\text{avasc.,}\infty} = e_{\text{avasc.,}\infty} = 625 \text{ mm}^3$; and (d) metastasis stage simulation from the vascular equilibrium point $c_{\text{vasc.,}\infty} = e_{\text{vasc.,}\infty} = 17346.5 \text{ mm}^3$.</figcaption></figure>

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:

$$ \dot{c}(t) = -\lambda_c c(t) \log \left( \frac{c(t)}{e(t)} \right) \quad (1) $$

$$ \dot{e}(t) = bc(t) - dc(t)^{2/3}e(t) \quad (2) $$

with $\lambda_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_{\infty} = e_{\infty} = \left( \frac{b}{d} \right)^{3/2} \quad (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_\infty$ 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. To understand the tumor dynamics, it is important to analyze its growth. The growth parameters set $\{\lambda_c, b, d\}$ has been identified by Hahnfeldt et al. [15] from experimental data with mice diseased with lung cancer (Lewis lung-carcinoma, LLC) for the vascular stage, and the

<table><caption>TABLE I
ANGIOGENESIS TUMOR GROWTH PARAMETERS SET.</caption><thead><tr><th>Param.</th><th>Value [15]</th><th>Units</th></tr></thead><tbody><tr><td>λ<sub>c</sub></td><td>0.192 / log(10)</td><td>day<sup>-1</sup></td></tr><tr><td>d</td><td>8.73 × 10<sup>-3</sup></td><td>day<sup>-1</sup> mm<sup>-2</sup></td></tr><tr><td>b</td><td>0.638</td><td>day<sup>-1</sup></td></tr></tbody></table><p>(a)</p><table><thead><tr><th>Stage</th><th>Avasc.</th><th>Vasc. [15]</th><th>Met.</th><th>Units</th></tr></thead><tbody><tr><td>b</td><td>0.638</td><td>5.85</td><td>8</td><td>day<sup>-1</sup></td></tr></tbody></table><p>(b)</p>

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 \text{ 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 \text{ 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_{\text{avasc.,}\infty} = e_{\text{avasc.,}\infty} = 625 \text{ 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_{\text{vasc},\infty} = e_{\text{vasc},\infty} = 17346.5 \text{ 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 $\delta(t)$, in the sampled fluid (e.g. plasma or blood). The agent concentration expression is [15]:

$$\delta(t) = \int_{0}^{t} u(s) \exp(-\lambda_{\delta}(t-s)) ds \quad (4)$$

where $u$ is the rate of administration of the therapeutic agent, and $\lambda_{\delta}$ is the elimination rate. The PD describes the effects of the agent to the body, and is related to the concentration $\delta(t)$ provided by the PK model. The classic PD model is the log-kill effect [15], [22]. 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 $\delta_a(t)$, leading to its new formulation [15]:

$$\dot{e}(t) = bc(t) - dc(t)^{2/3}e(t) - k_a\delta_a(t)e(t) \quad (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 [19]–[21]. Here, to model the effect of a cytotoxic chemotherapeutic agent of concentration $\delta_x(t)$ the angiogenesis model is used [19], and is extended as:

$$\dot{c}(t) = -\lambda_c c(t) \log \left( \frac{c(t)}{e(t)} \right) - k_{cx} c(t) \delta_x(t) \quad (6)$$

$$\dot{e}(t) = bc(t) - dc(t)^{2/3}e(t) - k_{ex} e(t) \delta_x(t) \quad (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}, \lambda_{\delta}$} characterized from experimental data are available in the literature.

3) **Radiotherapy Treatment:** Classically, the so-called tumor linear-quadratic (LQ) model [17], [18], [26] 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 [17], [18], [26]:

$$-x(t) \left( \alpha_x + \beta_x \int_{0}^{t} (u(s) \exp(-\mu(t-s)) ds) \right) u(t) \quad (8)$$

where $u(t)$ represents the fractionation schedule of radiation; $\alpha_x$ and $\beta_x$ are positive gains related to the tumor-LQ parameters; and $\mu$ 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 [18]:

$$\dot{c} = -c \left( \lambda_c \log \left( \frac{c}{e} \right) + (\alpha_c + \beta_c r) u \right) \quad (9)$$

$$\dot{e} = bc - e \left( dc^{2/3} + (\alpha_e + \beta_e r) u \right) \quad (10)$$

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 [27]. In particular, optimal control techniques have been applied for antiangiogenic therapy [16], [17], chemotherapy [21], [22], [28] or radiotherapy [17] to define the optimal treatment and drug dose. This section presents an adaptation of conventional therapies [15]–[19] 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 [29], [30]. 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 [19], [21], [23], [28]. One of these approaches consider the state-dependent Riccati equations (SDRE) framework [28], [30], [31]. Actually, SDRE has emerged as an interesting strategy for direct synthesis of nonlinear controllers [30], [31]. This optimal control methodology has been applied with the de Pillis et al. cancer model [22] in [28]. However, to the authors' knowledge, this optimal framework has not been applied to the Hahnfeldt's tumor growth formulation [15] or its extension [16]–[19].

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

$$\dot{\mathbf{x}} = \mathbf{A}(\mathbf{x})\mathbf{x} + \mathbf{B}(\mathbf{x})u \quad (11)$$

where $\mathbf{x} \in \mathbb{R}^n$ is the time-dependent state vector, with $n$ the number of states; $u \in \mathbb{R}^m$ is the (drug) control input ($m = 1$ for monotherapy); $\mathbf{A}(\mathbf{x}) \in \mathbb{R}^{n \times n}$ and $\mathbf{B}(\mathbf{x}) \in \mathbb{R}^{n \times m}$ are state-dependent matrices. In infinite horizon, the continuous-time SDRE optimization problem is then to minimize the following quadratic cost functional [29], [30]:

$$\mathbf{J} = \int (\mathbf{x}^T \mathbf{Q} \mathbf{x} + u^T \mathbf{R} u) dt \quad (12)$$

where the weights (design parameters) are state-dependent, such that $\mathbf{Q} \in \mathbb{R}^{n \times n}$ is positive-definite and $\mathbf{R} \in \mathbb{R}^{m \times m}$.
Thus, the optimal control law that minimizes this criterion is
given by:

$$
u = -\mathbf{R}^{-1}(\mathbf{x})\mathbf{B}^{\top}(\mathbf{x})\mathbf{P}(\mathbf{x})\mathbf{x} \quad (13)
$$

with $\mathbf{P}(\mathbf{x})$ is the unique, symmetric, positive-definite solution
of the algebraic state-dependent Riccati equation (SDRE) [30],
[31]:

$$
\mathbf{PA} + \mathbf{A}^{\top}\mathbf{P} - \mathbf{P}\mathbf{BR}^{-1}\mathbf{B}^{\top}\mathbf{P} + \mathbf{Q} = 0 \quad (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:

$$
\dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}) + \mathbf{g}(\mathbf{x})\mathbf{u}(t) \quad (15)
$$

To apply the above SDRE methodology, the so-called pseudo-
linearized form (11) of the nonlinear system (15) has to be
computed. More precisely, the state matrix is obtained from
$f(\mathbf{x}) = \mathbf{A}(\mathbf{x})\mathbf{x}$ by mathematical factorization. This state-
dependent parametrization of the nonlinear system is not
unique [30]. Among the alternatives, the chosen pseudo-
linearization should ensure point-wise controllability condi-
tion, where the so-called state-dependent controllability matrix
has full rank, that is:

$$
\mathcal{C} = [\mathbf{B}(\mathbf{x}) \ \mathbf{A}(\mathbf{x})\mathbf{B}(\mathbf{x}) \dots \mathbf{A}^{n-1}(\mathbf{x})\mathbf{B}(\mathbf{x})] \quad (16)
$$

$$
\mathrm{rank}(\mathcal{C}) = n \tag{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 ad-
ministration 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:

$$
\mathbf{x} = (x_1 = (c - c_{\infty}); \quad x_2 = (e - e_{\infty}); \quad x_3 = \delta_a)^{\top} \quad (18)
$$

where $c_{\infty} = e_{\infty} = x_{\infty}$ are the drug-free volume steady-state (3). A state-dependent parametrization of the angiogenesis inhibition treatment could be then written as:

$$
\mathbf{A}(\mathbf{x}) = \begin{bmatrix} a_{11} & \lambda_c \frac{x_1}{x_2} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ 0 & 0 & -\lambda_\delta \end{bmatrix}, \quad \text{and} \quad \mathbf{B}(\mathbf{x}) = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} \quad (19)
$$

with

$$
a_{11} = -\lambda_c \left( 1 + \log \left( \frac{x_1 + x_\infty}{x_2 + x_\infty} \right) \left( 1 + \frac{x_\infty}{x_1} \right) \right) \quad (20)
$$

$$
a_{13} = 0 \tag{21}
$$

$$
a_{21} = b \left( 1 + \frac{x_{\infty}}{x_1} \right) - \frac{x_{\infty}}{x_1} d (x_1 + x_{\infty})^{2/3} \quad (22)
$$

$$
a_{22} = -d(x_1 + x_\infty)^{2/3} - k_a x_3 \quad (23)
$$

$$
a_{23} = -k_a x_{\infty} \qquad (24)
$$

This pseudolinearization is used to apply the SDRE control of
the amount of anti-angiogenetic agent by considering $R = 10^3$

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

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

<figure><img src="image_48.png" /><figcaption>Fig. 4. Continuous infusion angiogenesis inhibition therapies: decreasing of (a) tumor volume and (b) endothelial over time with (c) optimal input drug.</figcaption></figure>

and $Q = \text{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_{\text{max}}$.

<table><caption>TABLE III
ENDOSTATIN PK-PD PARAMETERS [15], [32].</caption>
  <thead>
    <tr>
      <th>Param.</th>
      <th>Unit</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>k<sub>a</sub></td>
      <td>day<sup>-1</sup> conc<sup>-1</sup></td>
      <td>0.66</td>
      <td></td>
      <td></td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>&lambda;<sub>g</sub></td>
      <td>day<sup>-1</sup></td>
      <td>1.7</td>
      <td>1.37</td>
      <td>0.97</td>
      <td>0.57</td>
      <td>0.17</td>
      <td>0.017</td>
    </tr>
    <tr>
      <td>u<sub>max</sub></td>
      <td>conc day<sup>-1</sup></td>
      <td>15</td>
      <td>13</td>
      <td>11</td>
      <td>8</td>
      <td>5</td>
      <td>2</td>
    </tr>
  </tbody>
</table>

Fig. 4 illustrates optimal angiogenesis inhibition therapies
of a tumor at a vascular stage. Actually, tumors could be
<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /><figcaption>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.</figcaption></figure>

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 $\lambda_{\delta}$. Obviously, by varying $\lambda_{\delta}$ 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]–[17], endostatin antiangiogenic agents is considered here. The elimination rate of endostatin is varying from $1.7 \text{ day}^{-1}$ for a subcutaneous injection [15] to $0.017 \text{ day}^{-1}$ [32], 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], [32]. The Fig. 4a and 4b show that for a clearance of $\lambda_{\delta} = 0.017 \text{ day}^{-1}$ the tumor and the endothelial volumes vanish in about 60 days. In contrast, for $\lambda_{\delta} = 1.7 \text{ day}^{-1}$ the volumes reach a steady-state of $c_{\infty} = e_{\infty} = 476.5 \text{ mm}^3$.

However, the Fig. 4c illustrates that a continuous administration of 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 $\lambda_{\delta}$ 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 (18) 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_{\infty}) \quad (25)$$

$$a_{22} = -d(x_1 + x_{\infty})^{2/3} - k_{ex}x_3 \quad (26)$$

$$a_{23} = -k_{ex}x_{\infty} \quad (27)$$

SDRE method is then applied to control the amount of cytotoxic drug. The optimal controller parameters are defined here as: $R = 10^4$ and $Q = \text{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 [16], [17], [19] have considered experiments only. In this work the PK-PD parameters reported in Table IV are proposed based on d’Onofrio et al. analysis [19] to fit the DEB-TACE procedure pharmacokinetics [6]. 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$).

<table><caption>TABLE IV
CYTOTOXIC AGENT PK-PD PARAMETERS.</caption><thead><tr><td>Param.</td><td>Value</td><td>Units</td></tr></thead><tbody><tr><td>λ<sub>δ</sub></td><td>0.009</td><td>day<sup>-1</sup></td></tr><tr><td>k<sub>cx</sub></td><td>2</td><td>day<sup>-1</sup> conc<sup>-1</sup></td></tr><tr><td>k<sub>ex</sub></td><td>1</td><td>day<sup>-1</sup> conc<sup>-1</sup></td></tr></tbody></table>

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 \text{ conc/day}$, that is $0.0417 (\text{mg/kg})/\text{day}$ per session. Hence, after 6 sessions of drug administration the tumor vanish under a volume of $c_{\infty} = 12.67 \text{ 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 (18) and the state matrix $A(x)$ proposed in (19) could be considered, where only {$a_{13}, a_{22}, a_{23}$} parameters are modified as follows:

$$a_{13} = 0 \quad (28)$$

$$a_{22} = -d(x_1 + x_{\infty})^{2/3} \quad (29)$$

$$a_{23} = 0 \quad (30)$$
<figure><img src="image_51.png" /></figure>

<figure><img src="image_52.png" /><figcaption>Fig. 6. Cytotoxic agent direct delivery: decreasing of (a) tumor volume and endothelial over time with optimal input drug (b).</figcaption></figure>

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

<figure><img src="image_54.png" /><figcaption>Fig. 7. Radioisotope direct delivery: decreasing of (a) tumor volume and endothelial over time with optimal input drug (b).</figcaption></figure>

However, the state-dependent input matrix is now given by:

$$ \mathbf{B}(\mathbf{x}) = \begin{bmatrix} -(\alpha_c + \beta_c r)(x_1 + x_\infty) \\ -(\alpha_e + \beta_e r)(x_2 + x_\infty) \\ 1 \end{bmatrix} \quad (31) $$

The optimal amount of radioisotope is computed using the following SDRE controller parameters: $Q = \text{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 [17], [18], [26]. For the sake of simplicity, similar numerical values of the tumor-LQ parameters proposed in [17] 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_\infty = 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.

<table><caption>TABLE V
RADIOISOTOPE TUMOR-LQ MODEL PARAMETERS.</caption><thead><tr><th>Param.</th><th>Value</th><th>Units</th><th>Param.</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>α<sub>c</sub></td><td>0.7</td><td>Gy<sup>-1</sup></td><td>β<sub>c</sub></td><td>α<sub>c</sub>/10</td><td>Gy<sup>-2</sup></td></tr><tr><td>α<sub>e</sub></td><td>0.136</td><td>Gy<sup>-1</sup></td><td>β<sub>e</sub></td><td>α<sub>e</sub>/2.5</td><td>Gy<sup>-2</sup></td></tr><tr><td>μ</td><td>11.09</td><td>day<sup>-1</sup></td><td></td><td></td><td></td></tr></tbody></table>

<figure><img src="image_55.png" /><figcaption>Fig. 8. Magnetic bolus model embedding aggregates of iron-oxides microparticles. (a) General design and (b) oil-based ferrofluid bolus.</figcaption></figure>

# IV. OPTIMAL DELIVERY OF THERAPEUTIC AGENTS

## A. Structure

In 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 [13] 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 [13] for optimal size selection.

The effective drug load $V_\delta$ of a single bolus is a fraction of the total volume defined as $V_\delta = (1 - \tau_m)V$ with $\tau_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 $\tau_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 = \left[ \frac{u \cdot m}{V_{\delta} \rho_{\delta}} \right] \quad (32)$$

where $u$ is the rate of administration of the therapeutic agent; $m$ is the patient weight (considered as $m = 60$ kg hereafter); and $\rho_{\delta}$ 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 angiogenic 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 \text{ mg ml}^{-1}$. Hence, considering the upper limit $u_{\text{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 $\lambda_{\delta} = 0.017 \text{ day}^{-1}$ where the endostatin administration rate is lower ($u_{\text{max}} = 2 \text{ conc day}^{-1}$), $N_b = 4661$ magnetic boluses have to be injected continuously to administrate the daily dose.

<table><caption>TABLE VI
REQUIRED MAGNETIC BOLUSES THAT CARRY ENDOSTATIN TO TREAT A PATIENT OF AVERAGE WEIGHT $m = 60$ kg.</caption><thead><tr><th>Param.</th><th>Unit</th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>u<sub>max</sub></td><td>conc day<sup>-1</sup></td><td>15</td><td>13</td><td>11</td><td>8</td><td>5</td><td>2</td></tr><tr><td>Dose</td><td>mg day<sup>-1</sup></td><td>750</td><td>650</td><td>550</td><td>400</td><td>250</td><td>100</td></tr><tr><td>N<sub>b</sub></td><td>1 &times; 10<sup>3</sup></td><td>34.95</td><td>30.29</td><td>25.63</td><td>18.64</td><td>11.65</td><td>4.66</td></tr></tbody></table>

Similarly, if the considered scheduling of one session every week is applied (cf. Fig. 5), $N_b = 1.398 \times 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 [5], [6]. The drug eluting beads (DEB) have been easily loaded with doxorubicin to a recommended level of $25 \text{ mg ml}^{-1}$, with a maximum loading of $45 \text{ mg ml}^{-1}$ [5]. Here, we assume that the doxorubicin PK-PD model parameters is 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 [4]. $^{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 \pm 2.5 \mu 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 \times 10^{-3}$ GBq per bolus. Assuming an uniform distribution, the radiation dose intake is computed from [4] $u = \frac{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 = \left[ \frac{u \cdot m_l}{49.67 \times 3147 \times 2.5 \times 10^{-6}} \right] \quad (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-to-cell cooperative and competitive interactions [17]. 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 [16], [17], [22]. 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 [33]. 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 $\delta_{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 earlier-stage 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 high 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 $\mu$m) being selected first, followed by larger particles (500-to-700 $\mu$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 $\mu$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 [18]. To demonstrate the feasibility of MRN-assisted procedure, Bigot et al. [12] 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.

## REFERENCES

[1] R. L. Siegel, K. D. Miller, and A. Jemal, "Cancer statistics, 2015," *CA: a cancer journal for clinicians*, vol. 65, no. 1, pp. 5-29, 2015.

[2] J. M. Llovet, M. I. Real, X. Montaa, R. Planas, S. Coll, J. Aponte, C. Ayuso, M. Sala, J. Muchart, R. Sol, J. Rods, and J. Bruix, "Arterial embolisation or chemoembolisation versus symptomatic treatment in patients with unresectable hepatocellular carcinoma: a randomised controlled trial," *The Lancet*, vol. 359, no. 9319, pp. 1734-1739, 2002.

[3] A. Forner, J. M. Llovet, and J. Bruix, "Hepatocellular carcinoma," *The Lancet*, vol. 379, no. 9822, pp. 1245-1255, 2012.

[4] A. Kennedy, S. Nag, R. Salem, R. Murthy, A. J. McEwan, C. Nutting, A. Benson III, J. Espat, J. I. Bilbao, R. A. Sharma et al., "Recommendations for radioembolization of hepatic malignancies using yttrium-90 microsphere brachytherapy: a consensus panel report from the radioembolization brachytherapy oncology consortium," *International Journal of Radiation Oncology* *Biology* *Physics*, vol. 68, no. 1, pp. 13-23, 2007.

[5] A. L. Lewis, M. V. Gonzalez, A. W. Lloyd, B. Hall, Y. Tang, S. L. Willis, S. W. Leppard, L. C. Wolfenden, R. R. Palmer, and P. W. Stratford, "DC bead: in vitro characterization of a drug-delivery device for transarterial chemoembolization," *Journal of vascular and interventional radiology*, vol. 17, no. 2, pp. 335-342, 2006.

[6] M. Varela, M. I. Real, M. Burrel, A. Forner, M. Sala, M. Brunet, C. Ayuso, L. Castells, X. Montañá, J. M. Llovet et al., "Chemoembolization of hepatocellular carcinoma with drug eluting beads: efficacy and doxorubicin pharmacokinetics," *Journal of hepatology*, vol. 46, no. 3, pp. 474-481, 2007.

[7] T. Matsumoto, J. Endo, K. Hashida, H. Ichikawa, S. Kojima, S. Takashimizu, N. Watanabe, T. Yamagami, and T. Hasabe, "Balloon-occluded transarterial chemoembolization using a 1.8-french coaxial microballoon catheter for hcc: Tehnical and safety considerations," *Minimally Invasive Therapy Allied Technologies*, vol. 24, no. 2, pp. 94-100, 2015.

[8] R. Salem and K. Thurston, "Radioembolization with 90yttrium microspheres: A state-of-the-art brachytherapy treatment for primary and secondary liver malignancies- part 1: Technical and methodologic considerations," *Journal of Vascular Interventional Radiology*, vol. 17, pp. 1251-1278, 2006.

[9] U. Hafeli, S. Sweeney, B. Beresford, J. Humm, and R. Macklis, "Effective targeting of magnetic radioactive y-microspheres to tumor cells by an externally applied magnetic field. preliminary in-vitro and in-vivo resuls," *Nuclear Medecina and Biology*, vol. 22, no. 2, pp. 147-155, 1995.

[10] J.-B. Mathieu, G. Beaudoin, and S. Martel, "Method of propulsion of a ferromagnetic core in the cardiovascular system through magnetic gradients generated by an MRI system," *IEEE Trans. Biomed. Eng.*, vol. 53, no. 2, pp. 292-299, 2006.

[11] P. Pouponneau, J.-C. Leroux, G. Soulez, L. Gaboury, and S. Martel, "Co-encapsulation of magnetic nanoparticles and doxorubicin into biodegradable microcarriers for deep tissue targeting by vascular MRI navigation," *Biomaterials*, vol. 32, no. 13, pp. 3481-3486, 2011.

[12] A. Bigot, C. Tremblay, G. Soulez, and S. Martel, "Magnetic resonance navigation of a bead inside a three-bifurcation pmma phantom using an imaging gradient coil insert," *navigation (MRN)*, vol. 5, p. 7, 2014.

[13] L. Mellal, K. Belharet, D. Folio, and A. Ferreira, "Optimal structure of particles-based superparamagnetic microrobots: application to mri guided targeted drug therapy," *Journal of Nanoparticle Research*, vol. 17, no. 2, pp. 1-18, 2015.

[14] B. Ribba, N. Holford, P. Magni, I. Trocóniz, I. Gueorguieva, P. Girard, C. Sarr, M. Elishmereni, C. Kloft, and L. Friberg, "A review of mixed-effects models of tumor growth and effects of anticancer drug treatment used in population analysis," *CPT: pharmacometrics & systems pharmacology*, vol. 3, no. 5, p. e113, 2014.

[15] P. Hahnfeldt, D. Panigrahy, J. Folkman, and L. Hlatky, "Tumor development under angiogenic signaling a dynamical theory of tumor growth, treatment response, and postvascular dormancy," *Cancer research*, vol. 59, no. 19, pp. 4770-4775, 1999.

[16] A. d'Onofrio and A. Gandolfi, "Tumour eradication by antiangiogenic therapy: analysis and extensions of the model by hahnfeldt et al. (1999)," *Mathematical biosciences*, vol. 191, no. 2, pp. 159-184, 2004.

[17] A. Ergun, K. Camphausen, and L. M. Wein, "Optimal scheduling of radiotherapy and angiogenic inhibitors," *Bulletin of mathematical biology*, vol. 65, no. 3, pp. 407-424, 2003.

[18] U. Ledzewicz, A. dOnofrio, and H. Schttler, "Tumor development under combination treatments with anti-angiogenic therapies," in *Mathematical
*Methods and Models in Biomedicine*, ser. Lecture Notes on Mathematical Modelling in the Life Sciences, U. Ledzewicz, H. Schttler, A. Friedman, and E. Kashdan, Eds. Springer New York, 2013, pp. 311–337.

[19] A. d'Onofrio, U. Ledzewicz, H. Maurer, and H. Schättler, "On optimal delivery of combination therapy for tumors," *Mathematical biosciences*, vol. 222, no. 1, pp. 13–26, 2009.

[20] S. Sanga, J. P. Sinek, H. B. Frieboes, M. Ferrari, J. P. Fruehauf, and V. Cristini, "Mathematical modeling of cancer progression and response to chemotherapy," *Expert Review of Anticancer Therapy*, vol. 6, no. 10, pp. 1361–1376, 2006.

[21] L. G. de Pillis, W. Gu, K. R. Fister, T. a. Head, K. Maples, A. Murugan, T. Neal, and K. Yoshida, "Chemotherapy for tumors: An analysis of the dynamics and a study of quadratic and linear optimal controls," *Mathematical Biosciences*, vol. 209, no. 1, pp. 292–315, 2007.

[22] L. G. De Pillis and A. Radunskaya, "The dynamics of an optimally controlled tumor model: A case study," *Mathematical and Computer Modelling*, vol. 37, no. 11, pp. 1221–1244, 2003.

[23] L. Kovács, A. Szeles, J. Sápi, D. A. Drexler, I. Rudas, I. Harmati, and Z. Sápi, "Model-based angiogenic inhibition of tumor growth using modern robust control method," *Computer methods and programs in biomedicine*, vol. 114, no. 3, pp. e98–e110, 2014.

[24] A. Kansal, S. Torquato, G. Harsh, E. Chiocca, and T. Deisboeck, "Simulated brain tumor growth dynamics using a three-dimensional cellular automaton," *Journal of theoretical biology*, vol. 203, no. 4, pp. 367–382, 2000.

[25] P. Gerlee and A. R. Anderson, "An evolutionary hybrid cellular automaton model of solid tumour growth," *Journal of theoretical biology*, vol. 246, no. 4, pp. 583–603, 2007.

[26] J. F. Fowler, "The linear-quadratic formula and progress in fractionated radiotherapy," *The British journal of radiology*, vol. 62, no. 740, pp. 679–694, 1989.

[27] G. W. Swan, "Role of optimal control theory in cancer chemotherapy," *Mathematical biosciences*, vol. 101, no. 2, pp. 237–284, 1990.

[28] Y. Batmani and H. Khaloozadeh, "Optimal chemotherapy in cancer treatment: state dependent riccati equation control and extended kalman filter," *Optimal Control Applications and Methods*, vol. 34, no. 5, pp. 562–577, 2013.

[29] K. Zhou, J. C. Doyle, K. Glover et al., *Robust and optimal control*. Prentice Hall New Jersey, 1996, vol. 40.

[30] C. P. Mracek and J. R. Cloutier, "Control designs for the nonlinear benchmark problem via the state-dependent riccati equation method," *International Journal of robust and nonlinear control*, vol. 8, no. 4-5, pp. 401–433, 1998.

[31] T. Çimen, "Systematic and effective design of nonlinear feedback controllers via the state-dependent riccati equation (sdre) method," *Annual Reviews in Control*, vol. 34, no. 1, pp. 32–51, 2010.

[32] S.-U. Shin, H.-M. Cho, J. Merchan, J. Zhang, K. Kovacs, Y. Jing, S. Ramakrishnan, and J. D. Rosenblatt, "Targeted delivery of an antibody–mutant human endostatin fusion protein results in enhanced antitumor efficacy," *Molecular cancer therapeutics*, vol. 10, no. 4, pp. 603–614, 2011.

[33] P. Anand, A. B. Kunnumakara, C. Sundaram, K. B. Harikumar, S. T. Tharakan, O. S. Lai, and B. B. Aggarwal, "Cancer is a preventable disease that requires major lifestyle changes," *Pharmaceutical research*, vol. 25, no. 9, pp. 2097–2116, 2008.
# APPENDIX E

## MATHEMATICAL APPROACH FOR THE DESIGN
CONFIGURATION OF MAGNETIC SYSTEM WITH
MULTIPLE ELECTROMAGNETS (ACL15)
Mathematical approach for the design configuration of magnetic system
with multiple electromagnets

Ruipeng Chen, David Folio, Antoine Ferreira*

INSA Centre Val de Loire, Université d'Orléans, PRISME EA 4229, Bourges, France

**Abstract**

Magnetic actuation techniques and microrobots have attracted great interest since they have potential in biomedicine applications. Interventional techniques have emerged as a tool to handle a wide range of mini-mally invasive surgeries (MIS). However, current MIS procedures are constrained by the limitation of manual operation by surgeon. Thus, various microrobotic solutions including magnetic navigation systems have been proposed for MIS, which carries many potential benefits such as reduced incision, less intraoperative hemorrhaging and postoperative pain, and faster recovery time. In recent decades, many electromagnetic actuation (EMA) systems have been reported and involved to general surgery. The EMA system allows to generate efficiently magnetic source for microrobot control when its specifications are further investigated and satisfied for the desired application. To precisely manipulate the biomedical microrobot, a key issue still relies on the design of a suitable EMA platform. In this paper, we demonstrate a mathematical approach for the design configuration of magnetic system with multiple electromagnets. Especially, the required magnetic coil number has been investigated where the heading motion control, magnetic force control and their combination control are discussed respectively. The singular cases of control are pre-evaluated by a mathematical analysis of the simulated electromagnetic field. In addition, the placed positions and tilted orientations of the applied electromagnets are investigated for the optimization regarding the six typical configurations of EMA platform with 4, 6 and 8 coils. The various configurations of EMA systems have been comprehensively analyzed. Therefore, with the number of electromagnets and their optimal configuration obtained by the proposed approach, the EMA system can be initially established.

**Keywords:** Electromagnetic actuation system, Magnetic microrobot, Design methodology, Robotic magnetic platform

**1. Introduction**

Minimally invasive surgery (MIS) encompasses surgical techniques that limit the size of incisions needed and so lessen wound healing time, associated pain and risk of infection. MIS procedures have been enabled by the advance of various medical technologies. In particular, several surgical robotics systems have been developed. Such robotic platforms already play a significant role to improve patient care, though it increases surgical preparation, cost and risk of the MIS approach compared with traditional open surgery[1,

2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. For example, there are the well-known da Vinci surgical assistance robots[5, 7], developed by Intuitive Surgical1, which improves the surgeon technical skills. Moreover, unlike the dependence on using rigid instruments with dexterous distal wrists, it is commonly more attractive to use flexible or adaptive robotic tools that access internal anatomy with few skin incisions[10, 12, 13, 14, 15, 16, 17]. Especially, the mechanical parts of existing medical robotic devices are still relatively large and rigid to access and treat major inaccessible parts of the human body (e.g. in robot-assisted surgery). In parallel, the various medical robotics solutions have been developed to improve the acceptance of the use of robotics systems in clinical practices. In the meanwhile, mi-

<footnote>*Corresponding author
Email addresses: david.folioinsa-cvl.fr (David Folio), antoine.ferreira@insa-cvl.fr (Antoine Ferreira)</footnote>
crorobotics has also emerged as an attractive technology to introduce novel microsystems to further reduce trauma, create new diagnosis tools and therapeutic procedures.

Indeed, designing miniaturized and versatile microrobotic systems would allow accessing throughout the whole human body; leading to new procedures down to the cellular level; and offering localized diagnosis and treatment with greater precision and efficiency. For example, untethered microscopic devices, smaller than one millimeter, may navigate within the body for targeted therapies [18, 19, 20, 21]. Among the various actuation of microrobots, magnetic actuation is considered to be the most promising method [9, 18, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. To this aim, numerous electromagnetic actuation (EMA) systems have been proposed to control untethered magnetic microrobots for biomedical applications [18, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47]. The development of magnetic microrobotic systems circumvents the need to embed power within the microdevice [18, 19, 20]. Magnetic microrobot can be then remotely powered and actuated precisely by the use of external electromagnetic fields. This enables untethered microrobots to assist the surgeon to increase precision and dexterity of the MIS procedure [10, 18, 34, 42].

The paper aims to propose a mathematical approach for the design configuration of magnetic system using multiple electromagnets, that improves MIS procedures since the magnetic manipulation of microrobots is one of the most interesting method to assist researchers in many biomedical applications. To address this issue of designing a proper EMA system for a given application, the investigation of the appropriate design configuration carried out. The paper is then organized as follows. Section 2 presents the theoretical foundation of electromagnetism including the basic principles and necessary equations. From this basis, the electromagnetic manipulation of untethered microrobot with the magnetic torque and force control is discussed, and the necessary metrics are defined. In Section 3, the analysis of the minimum number of electromagnets for an EMA platform is investigated. In particular, singular cases are pointed out when magnetic field and gradient exhibit some linear dependencies. In Section 4, six typical configurations are simulated and evaluated for the efficient motion control of microrobot. The configuration of EMA system is optimized with respect to the considered magnetic

<figure><img src="image_56.png" /><figcaption>Figure 1: Electromagnetic manipulation: (a) illustration of the use of the magnetic force and torque on untethered magnetic microrobots; (b) composition of the magnetic field with the magnetic moment of the microrobot; and (c) the different components of the magnetic gradient.</figcaption></figure>

control objectives. In Section 5, the comparison of configurations is performed when the magnetic performances are summarized regarding the arrangement parameters of magnetic system. The paper is concluded with the Section 6.

## 2. Theoretical foundation

The aim of this section is to recall the basic principles of electromagnetism that serve as a foundation of our research works.

### 2.1. Electromagnetic Manipulation

Starting from the Maxwell's equations, and assuming that charges are either fixed or move as a steady current $J$, the governing equation of a *quasi-static magnetic field* can be described by the two following relations [48]:

$$ \nabla \cdot B = 0 \qquad (1) $$

$$ \nabla \times B = \mu_0 J \qquad (2) $$

Next, assuming that the microrobot is a magnetized body described by its magnetic dipole moment $M$, that is placed in a magnetic flux $B$, the induced magnetic force and torque are basically expressed from [49]:

$$ \mathbf{f} = (M \cdot \nabla) B \qquad (3) $$

$$ \mathbf{t} = (M \times B) \qquad (4) $$
<figure><img src="image_57.png" /><figcaption>Figure 2: Schematic representation of a single electromagnet inducing a magnetic flux $B_e(P)$ within the workspace $\Omega$ (blue square box).</figcaption></figure>

Obviously, **f** is related to the magnetic gradient $\nabla B$, whereas the magnetic torque **t** is a function of the magnetic field **B**. Moreover, to actuate the microrobot, the magnetic field must undergo either a spatial change (i.e. exhibit a spatial field gradient), or a temporal change, such as through a rotating magnetic field, oscillation and so on. Fig. 1 illustrates these basic principles. Specifically, the magnetic field and its gradient are generated by a set of n electromagnets to control the motion of undethered microrobots. The orientation of microrobot tends to be aligned with the direction of magnetic field **B**, thus the microrobots are rotated by the magnetic torque **t**. Moreover, the magnetic gradient $\nabla B$ induces the magnetic force **f** as propulsion force to move the microrobot. Hence, the magnetic object can be manipulated through the rotation and translation operations using magnetic field and gradient, respectively.

Secondly, if we assume there is no electric current flowing through the workspace occupied by the microrobot, Maxwell's equation implies that equation (2) becomes $\nabla \times B = 0$. The magnetic force (3) can be then rearranged using vector calculus into the following form:

$$ \mathbf{f} = \begin{pmatrix} \frac{\partial B}{\partial x} \\ \frac{\partial B}{\partial y} \\ \frac{\partial B}{\partial z} \end{pmatrix} \quad M = \begin{pmatrix} \frac{\partial b_x}{\partial x} & \frac{\partial b_y}{\partial x} & \frac{\partial b_z}{\partial x} \\ \frac{\partial b_x}{\partial y} & \frac{\partial b_y}{\partial y} & \frac{\partial b_z}{\partial y} \\ \frac{\partial b_x}{\partial z} & \frac{\partial b_y}{\partial z} & \frac{\partial b_z}{\partial z} \end{pmatrix} \quad M = \mathbf{G}M \quad (5) $$

where $\mathbf{G} \in \mathbb{R}^{3 \times 3}$ denotes the gradient matrix of the magnetic field $B = (b_x, b_y, b_z)$. Besides, the subscripts $x$, $y$ and $z$ explicitly refer to the basis directions of the Cartesian reference frame $\mathcal{F}_0(\mathcal{O}: x, y, z)$ linked to the workspace $\Omega$, in which all vectors are expressed (see Fig. 2).

The torque on microrobot tends to align the magnetization vector with the magnetic field. To represent vector cross products, the skew-symmetric

matrix formed of a vector can be employed, that is:

$$ \mathrm{Sk}(M) = \mathrm{Sk} \begin{pmatrix} m_x \\ m_y \\ m_z \end{pmatrix} = \begin{pmatrix} 0 & -m_z & m_y \\ m_z & 0 & -m_x \\ -m_y & m_x & 0 \end{pmatrix} \quad (6) $$

Therefore, the force (3) and torque (4) applied to magnetic microrobots can be rewritten as:

$$ \begin{pmatrix} \mathbf{f} \\ \mathbf{t} \end{pmatrix} = \begin{pmatrix} (M \cdot \nabla) \\ \mathrm{Sk}(M) \end{pmatrix} B \quad (7) $$

## 2.2. Electromagnetic control

For multiple-coil configuration, the magnetic field is induced by an EMA system consisting of several electromagnets. Hence, with a given set of n electromagnets, each of them creates a magnetic field, $B_e(P)$, at any location of the workspace: $\forall P \in \Omega$. The magnetic fields and their spatial gradients depend linearly on the currents. When $i_e$ flows through the coils $e_z$, the magnetic field can be expressed as: $B_e(P) = \tilde{B}_e(P)i_e$. The superposition principle is commonly used to compute the overall magnetic field, that is:

$$ \begin{aligned} B(P) &= \sum_{e=1}^{n} \tilde{B}_e(P) i_e = \begin{pmatrix} \tilde{b}_{x1} & \dots & \tilde{b}_{xe} & \dots & \tilde{b}_{xn} \\ \tilde{b}_{y1} & \dots & \tilde{b}_{ye} & \dots & \tilde{b}_{yn} \\ \tilde{b}_{z1} & \dots & \tilde{b}_{ze} & \dots & \tilde{b}_{zn} \end{pmatrix} \mathbf{i} \\ &= \begin{pmatrix} \mathcal{B}_x(P) \\ \mathcal{B}_y(P) \\ \mathcal{B}_z(P) \end{pmatrix} \mathbf{i} = \mathcal{B}(P)\mathbf{i} \end{aligned} \quad (8) $$

where $\mathbf{i} = (i_1 \dots i_n)^\mathrm{T}$ is the current vector, and $\mathcal{B}(P)$ is a $3 \times n$ matrix mapping the currents to the magnetic fields. Similarly, the magnetic gradient fields could be expressed as:

$$ \frac{\partial B(P)}{\partial x} = \sum_{e=1}^{n} \frac{\partial \tilde{B}_e(P)}{\partial x} i_e = \mathcal{G}_x(P)\mathbf{i} \quad (9) $$

$$ \frac{\partial B(P)}{\partial y} = \sum_{e=1}^{n} \frac{\partial \tilde{B}_e(P)}{\partial y} i_e = \mathcal{G}_y(P)\mathbf{i} \quad (10) $$

$$ \frac{\partial B(P)}{\partial z} = \sum_{e=1}^{n} \frac{\partial \tilde{B}_e(P)}{\partial z} i_e = \mathcal{G}_z(P)\mathbf{i} \quad (11) $$

where each $\mathcal{G}_{x,y,z}(P)$ is $3 \times n$ matrix mapping the current to the magnetic gradient field in the $x$, $y$, and $z$ directions, respectively.

The different mapping matrices $\mathcal{B}(P)$ and $\mathcal{G}_{x,y,z}(P)$ expressions are related to the type and
geometry of each electromagnet (e.g. length, radius, numbers of turns). For instance, the mapping for Helmholtz or Maxwell coils could be easily carried out from the Biot-Savart's law. For the other types of electromagnets, especially with a magnetic core, the mapping matrices can be calculated either from measures of the magnetic field within the EMA system's workspace or from the simulations of magnetic field using FEM models [43].

Finally, the equations of magnetic force and torque can be rearranged as follows:

$$ 
\begin{aligned}
\begin{pmatrix} \mathbf{t} \\ \mathbf{f} \end{pmatrix} &= \begin{pmatrix} \mathrm{Sk}(M)\mathcal{B} \\ M^\mathrm{T}\mathcal{G}_x(P) \\ M^\mathrm{T}\mathcal{G}_y(P) \\ M^\mathrm{T}\mathcal{G}_z(P) \end{pmatrix} \mathbf{i} = \begin{pmatrix} \mathcal{A}_t(M, P) \\ \mathcal{A}_f(M, P) \end{pmatrix} \mathbf{i} \quad (12) \\
&= \mathcal{A}(M, P)\mathbf{i}
\end{aligned}
$$

where $\mathcal{A}_t(M, P)$ and $\mathcal{A}_f(M, P)$ are the actuation matrices mapping the current vector $\mathbf{i}$ to the torque and force, respectively; and $\mathcal{A}(M, P)$ is the $6 \times n$ actuation matrix mapping to the applied magnetic wrench. These magnetic actuation matrices depend both on the position $P \in \Omega$, and on the orientation of the dipole moment $M$ of the microrobot.

According to equations (3)-(4), each column of the matrix $\mathcal{A}(M, P)$ represents the wrench on the force and torque per current unit created by each electromagnet. If there are greater than $n > 6$ electromagnets, the actuation matrix $\mathcal{A}(M, P)$ leads to a better conditioned matrix, a more isotropic workspace $\Omega$, a reduction of singularity configurations, and lower current requirements [50, 51]. In such cases, $n > 6$, the EMA system could be said "redundant" for the task. Especially, if $\mathcal{A}(M, P)$ is of full rank, for a desired force, $\mathbf{f}^*$ and torque, $\mathbf{t}^*$, the actuation currents $\mathbf{i}$ can be calculated from the pseudo-inverse:

$$ \mathbf{i} = \mathcal{A}^+(M, P) \begin{pmatrix} \mathbf{t}^* \\ \mathbf{f}^* \end{pmatrix} \quad (13) $$

If $n < 6$, the pseudo-inverse would be a least-squares approximations. Hence, for a controlled force and torque, the input current can be obtained only if the pseudo-inverse of $\mathcal{A}(M, P)$ exists. This derivation on the controlled current $\mathbf{i}$ can be similarly extended for controllers that require torque and/or force control [52].

## 2.3. Indexes for Electromagnetic System

To be able to compare such various electromagnetic coils dispositions, it is necessary to define

some metrics to get some quantitative assessments of their effectiveness. The basic idea is to set a score to a given EMA configuration to allow an effective comparative analysis. Indeed, with respect to the applications objectives, different indexes can be defined and considered.

### 2.3.1. Magnetic Field Indexes

First, the strength and homogeneity of either the magnetic field or its gradient in the workspace are the most significant metrics for medical applications. As example, a strong magnetic field strength leads to powerful torque $\mathbf{t}$ enabling reliable swimming of the helical microrobot. In the meantime, it should be homogeneous along the workspace to simplify the control strategy. Furthermore, the EMA system needs to avoid unexpected magnetic force when varying the orientation of microrobot along with the magnetic field. Similarly, a strong magnetic gradient provides an effective propulsion force $\mathbf{f}$, and its uniformity implies a homogeneous net force on microrobots. To this aim, various indexes can be used from statistical analysis of the fields distributions, such as its maximum, mean, RMS, STD, etc. To evaluate the quality of a magnetic field or its gradient, we have chosen to consider mainly the average and the uniformity indexes introduced hereafter.

**Average index.** In [40], the authors propose to characterize the strength of a sampled vector field $\varphi$ by its average (i.e. the arithmetic mean) basically defined as:

$$ \langle \varphi \rangle = \frac{1}{N} \sum_{P \in \Omega}^{N} \varphi(P) \quad (14) $$

with $N$ the number of samples, and $P$ the sampling location in the workspace $\Omega$. This arithmetic mean value reflects the strength of a global vector field $\varphi$. It can be utilized to quantify the average of the magnitudes of a magnetic field or its gradient from the corresponding samples, and then to evaluate the strength of the torque or force in the workspace.

**Uniformity index.** The isotropy index measures and yields an intuitive value for the difference between each considered vector of a vector field $\varphi$ distribution. To evaluate the isotropy, various indexes can be considered, such as the basic minimum-to-maximum ratio, or the coefficient of variation. For a vector field $\varphi$, the uniformity index $\gamma$ is usually
considered, which is defined as:

$$ \gamma(\varphi) = 1 - \frac{1}{2N\langle\varphi\rangle} \sum_{P \in \Omega}^{N} |\varphi(P) - \langle\varphi\rangle| \quad (\text{in \%}) \qquad (15) $$

where the isotropy index $\gamma$ should be commonly bounded between 0% and 100%. Then, a uniformity index close to 100% indicates that every sample of the field $\varphi$ is almost identical.

### 2.3.2. Magnetic Actuation Indexes

In the other hand, the EMA system has to actuate a magnetic microrobot in the workspace. Hence, the ability to perform or not certain motion at any location and direction is of prime importance. Basically, if there are some singularities, some motions cannot be achieved, and the DOFs of the micro-robot is reduced. For an EMA system, this analysis can be performed through the examination of its magnetic actuation matrix $A(M, P)$. As described in (12), $A(M, P)$ allows mapping the input currents $i$ to the magnetic torque and force that are applied to a microrobot possessing a magnetic moment $M$ at the location $P \in \Omega$. Hence, the actuation matrix $A(M, P)$ can be seen as a Jacobian matrix. Commonly, the Jacobian matrix of a system is used to characterize the control efficiency. Here, the columns of the actuation matrix $A(M, P)$ represents the wrench on the magnetic microrobot. As long as $A(M, P)$ is singularity-free or merely well-conditioned, the direct mapping will exist and full wrench control of the microrobot can be achieved in the workspace. Thus, the structure of $A(M, P)$ has to be analyzed to determine the "wrench-kinematic" performance of the EMA system. Therefore, the magnetic actuation matrix provides a similar information as a classic Jacobian matrix in robotics.

There are numerous studies on the kinematic performance indexes of robotic mechanisms [53, 54, 55, 56]. Most of proposed metrics were derived from the definition of manipulability index, introduced by [54].

**Manipulability index.** From classic robotic kinematics, measuring the manipulability is a well-known technique for determining the ability to maneuver in workspace. Specifically, manipulability describes the degree to which a robot can freely apply forces and torques in arbitrary directions, and quantifies the ability to perform an action quickly

<figure><img src="image_58.png" /><figcaption>Figure 3: 2D illustration of the manipulability ellipsoid: microrobots have better kinematic capacity in the direction of the major axis $\sigma_1$ of the ellipsoid, whereas the minor axis $\sigma_2$ induces worse kinematic performance.</figcaption></figure>

and skillfully [54]. To do so, the manipulability index is defined as a quality measure describing the distance to singular configurations. The approach is based on analyzing the manipulability ellipsoid that is spanned by the singular vectors of the Jacobian. Similarly, we investigate the mapping efficiency between the current input $i$ to the torque and force with this standard manipulability ellipsoid. Fig. 3 illustrates the manipulability ellipsoid, here in 2D, where the minor axis $\sigma_2$ represents the direction with the worse kinematic capacity, whereas the major axis $\sigma_1$ gives the easiest direction of force/torque transmission. Specifically, from the manipulability ellipsoid, the product of the ellipsoid's axes leads to the manipulability index which gives a measure of the maneuverability of the EMA system, and it is defined as [54]:

$$ w(\mathcal{A}) = \sqrt{\det(\mathcal{A}\mathcal{A}^\mathrm{T})} \qquad (16) $$

Commonly, the manipulability ellipsoids can be computed from SVD factorization. Let the SVD of a matrix $\mathcal{A}$ be:

$$ \mathcal{A} = U\Sigma V^\mathrm{T} \qquad (17) $$

where the $U$ and $V$ are orthogonal matrices, with $U$ a $6 \times 6$ unitary matrix and $V$ a $n \times n$ unitary matrix. In our considered system, $\Sigma$ is a $6 \times n$ singular value matrix given as follows:

$$ \Sigma = \begin{pmatrix} \sigma_1 & 0 & 0 & 0 & 0 & 0 & \cdots & 0 \\ 0 & \sigma_2 & 0 & 0 & 0 & 0 & \cdots & 0 \\ 0 & 0 & \sigma_3 & 0 & 0 & 0 & \cdots & 0 \\ 0 & 0 & 0 & \sigma_4 & 0 & 0 & \cdots & 0 \\ 0 & 0 & 0 & 0 & \sigma_5 & 0 & \cdots & 0 \\ 0 & 0 & 0 & 0 & 0 & \sigma_6 & \cdots & 0 \end{pmatrix} \qquad (18) $$

where $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_6 \ge 0$ are the singular values of $\mathcal{A}$. The singular values are significant indexes for the evaluation of system and that can be
considered as a scalar gain by which each input is multiplied to produce the corresponding output.

Since the singular values give information about the quality of the workspace [54], the product of all singular values can be used to analyze the control capability of points or configurations in certain directions. Thereby, the measure of manipulability index $w$ can be expressed as the product of the singular values:

$$w(\mathcal{A}) = \prod_{e=1}^{n} \sigma_e \qquad (19)$$

It can be noticed that $w$ is proportional to the volume of the manipulability ellipsoid. Generally, a larger measure of $w$ indicates a better conditioning of the actuation matrix $\mathcal{A}(M, P)$, that means an effective control of the microrobot. Hence, the manipulability index $w$ measures the "distance" from singularities, that represents the capability of the EMA system to control efficiently the magnetic microrobot.

Finally, as the above manipulability index (16) depends on the scale or units, the normalized manipulability can be preferred:

$$w_n = \frac{w(\mathcal{A}(M, P))}{\max_{P \in \Omega} w(\mathcal{A}(M, P))} \qquad (20)$$

*Condition number.* Another way to characterize the manipulability ellipsoid is to measure its isotropy. This is commonly achieved by computing the condition number of a matrix, which is defined as:

$$\kappa(\mathcal{A}) = \|\mathcal{A}\|\|\mathcal{A}^\dagger\| \in [1; +\infty) \qquad (21)$$

where $\mathcal{A}^\dagger$ is the Moore-Penrose pseudo-inverse of $\mathcal{A}$, and $\|\cdot\|$ denotes the Euclidean (or Frobenius) norm defined as:

$$\|\mathcal{A}\| = \sqrt{\mathrm{tr}(\mathcal{A}\mathbf{W}\mathcal{A}^\mathrm{T})} \qquad (22)$$

with $\mathbf{W}$ a weighting matrix, and especially, $\mathbf{W} = \mathbf{I}$ for the Euclidean norm.

To evaluate the control accuracy and isotropy of system, the inverse of the condition number of a matrix, also termed as *isotropic index*, was introduced by [53], and is defined as:

$$1/\kappa = \frac{\sigma_{\min}}{\sigma_{\max}} = \frac{\sigma_6}{\sigma_1} \in [0; 1] \qquad (23)$$

The isotropic index indicates how well the micro-

in some sense the dexterity of the microrobot [56]. Therefore, when $1/\kappa(\mathcal{A}(M, P))$ is close to 1, its manipulability ellipsoid is more close to be a sphere, which means the EMA system has more capacity to transmit the same force or torque in all directions (see also Fig. 3).

Let us notice that the previous indexes are only local metrics for a dedicated location $P \in \Omega$. To evaluate the system efficiency over the workspace $\Omega$, many global performance indexes have been proposed [55, 56]. Classically, to define such global performance indexes, the integral of a local performance index $\xi$ over a domain $\mathcal{W}$ is considered given by (24) or its sampled expression (25):

$$\Gamma_{\xi}(\mathcal{A}) = \frac{\int_{\mathcal{W}} \xi(\mathcal{A}) d\mathcal{W}}{\int_{\mathcal{W}} d\mathcal{W}} \qquad (24)$$

$$\Gamma_{\xi}(\mathcal{A}) = \frac{1}{N} \sum_{x \in \mathcal{W}}^{N} \xi(\mathcal{A}, x) \qquad (25)$$

When the domain is the workspace, that is $\mathcal{W} = \Omega \subset \mathbb{R}^3$ and $x = P$, and for $\xi = w_n$, $\Gamma_{w_n}(\mathcal{A})$ leads to the *global manipulability index*. Similarly, for $\xi = 1/\kappa$, $\Gamma_{1/\kappa}(\mathcal{A})$ denotes the *global isotropic index*, also termed as the *global conditioning index* [55, 56].

<figure><img src="image_59.png" /><figcaption>Figure 4: Spherical coordinates $(\rho, \theta, \varphi)$ with $\rho$ the radial distance, $\theta$ the polar angle, and $\varphi$ azimuthal angle.</figcaption></figure>

Likewise, we can evaluate the performance indexes for any direction of the magnetic moment $M$ by considering the 3D rotation group as domain, leading to: $\mathcal{W} = \mathrm{SO}(3)$, where it means the magnetic moment $M$ has 3D orientations and can be aligned to $x$-, $y$- and $z$-directions. In the following, to express the orientation of the magnetic moment, we will consider the spherical coordinate representation, as shown in Fig. 4. Obviously, when both location $P \in \Omega$ of the microrobot and the orientation of its magnetic moment $M$ are considered, the domain is defined accordingly, that is $\mathcal{W} = \Omega \times \mathrm{SO}(3)$.
In this work, the global performance metrics consisting of the average index (14), the uniformity index (15), the normalized manipulability index (20), and the (inverse) condition number (23) for EMA systems are proposed. The considered configuration of electromagnets is estimated in relation with these proposed performance metrics. Thereby, the various configurations can be quantitatively evaluated for their magnetic field distribution and activation efficiency.

## 3. DESIGN AND MODELING

### 3.1. The estimation of magnetic coil number

One of the motivation of an EMA platform is to provide the necessary DOFs for the manipulation of microrobot, mainly depending on the configuration and the number $n$ of constituting electromagnets. Intuitively, it can be shown that platforms with different arrangements of electromagnets exhibit a wide diversity of operating performances for various manipulation tasks [43, 57]. In this section, we will rigorously analyze the number of electromagnets that is required for different magnetic manipulations. Commonly, the manipulation of untethered magnetic microrobot includes the control of its orientation and position. As we mentioned in Fig. 1a, the orientation of microrobot can be simply controlled through the magnetic field. In addition, its position is reached using drilling force generated by rotating magnetic field or translating force generated by magnetic gradient. Thus, the efficient remote magnetic manipulation relies on the capability of generating proper magnetic field and gradient. Specifically, the DOFs of wireless microrobot can be evaluated through the magnetic field and its gradient that are generated by the EMA system with a given number of electromagnets.

As demonstrated in [58, 59, 60, 61], at least $n = 3$ electromagnets are required to achieve 3 DOFs pointing control. Besides, the position control can be effectively realized with at least $n = 4$ magnets in 3D workspace, but up to 5 coils are commonly used to improve the system stability [62]. As computed, 3 DOFs pointing control of the microrobot can be achieved by 2 DOFs torque manipulation. There is no torque that can be applied on the microrobot about its dipole moment main axis. Indeed, some articles have reports several methods to reach 3 DOFs torque control on the designed special microrobots [63, 64, 65]. For the minimum number of electromagnets, the non-magnetic restoring

forces, such as gravity, play also an important role, and will be considered. Moreover, some unexpected singularities involving to force control lead to some holonomic constraints, that will be investigated in the follows.

#### 3.1.1. Heading motion control

The quasi-static magnetic fields generated by electromagnets can be defined from the Maxwell's equations (1)-(2). The magnetic torque and force acting on a microrobot are described following (12). In current free space, (1) constraints the gradient matrix of the vector field B to have zero trace; and (2) constraints the gradient matrix of the vector field to be symmetric. Hence, the magnetic force (5) can be rearranged as:

$$ \mathbf{f} = \begin{pmatrix} \frac{\partial b_x}{\partial x} & \frac{\partial b_x}{\partial y} & \frac{\partial b_x}{\partial z} \\ \frac{\partial b_x}{\partial y} & \frac{\partial b_y}{\partial y} & \frac{\partial b_y}{\partial z} \\ \frac{\partial b_x}{\partial z} & \frac{\partial b_y}{\partial z} & -\left(\frac{\partial b_x}{\partial x} + \frac{\partial b_y}{\partial y}\right) \end{pmatrix} \begin{pmatrix} m_x \\ m_y \\ m_z \end{pmatrix} \quad (26) $$

This rearrangement allows simplifying the number of magnetic gradient components from nine to five. It is obvious that the magnetic force relies on the applied magnetic gradient field and the magnetic dipole moment of the microrobot. Thus, for a given magnetic dipole, the induced magnetic force only depends on the magnetic gradient that is controlled by the currents flowing through the electromagnets. The magnetic force equations can be thereby expressed as the follows to highlight the magnetic gradients as the controllable parameters:

$$ \mathbf{f} = \begin{pmatrix} m_x & m_y & m_z & 0 & 0 \\ 0 & m_x & 0 & m_y & m_z \\ -m_z & 0 & m_x & -m_z & m_y \end{pmatrix} \begin{pmatrix} \frac{\partial b_x}{\partial x} \\ \frac{\partial b_x}{\partial y} \\ \frac{\partial b_x}{\partial z} \\ \frac{\partial b_y}{\partial y} \\ \frac{\partial b_y}{\partial z} \end{pmatrix} = \mathcal{F}(M)\mathbf{g} \quad (27) $$

where $\mathcal{F}(M)$ indicates the matrix form of the dipole moment $M$ of the microrobot; and the vector field $\mathbf{g}$ includes the five magnetic gradient components.
Similarly, the applied magnetic torque is dependent on the magnetic field when the dipole moment of the microrobot is determined. Using the skew matrix (6), the Eq. (7) can be rewritten as:

$$ \begin{pmatrix} \mathbf{t} \\ \mathbf{f} \end{pmatrix} = \begin{pmatrix} \mathrm{Sk}(M) \\ \mathcal{F}(M) \end{pmatrix} \begin{pmatrix} \mathbf{B} \\ \mathbf{g} \end{pmatrix} \qquad (28) $$

Considering the current flowing through the n electromagnetic coils of the EMA system has a linearly mapping to the magnetic field (8) and gradient (9)-(11), it follows that the (12) can be then represented as:

$$ \begin{pmatrix} \mathbf{t} \\ \mathbf{f} \end{pmatrix} = \mathcal{A}(M, P)\mathbf{i} = \begin{pmatrix} \mathcal{A}_t(M, P) \\ \mathcal{A}_f(M, P) \end{pmatrix} \mathbf{i} \\ = \begin{pmatrix} \mathrm{Sk}(M) \\ \mathcal{F}(M) \end{pmatrix} \begin{pmatrix} \mathcal{B}(P) \\ \mathcal{G}(P) \end{pmatrix} \mathbf{i} \qquad (29) $$

where $\mathcal{B}$ is a $3 \times n$ matrix and $\mathcal{G}$ is a $5 \times n$ matrix, that are defined as:

$$ \mathbf{g} = \mathcal{G}(P)\mathbf{i} \\ = \begin{pmatrix} \frac{\partial \mathcal{B}_x}{\partial x} & \frac{\partial \mathcal{B}_x}{\partial y} & \frac{\partial \mathcal{B}_x}{\partial z} & \frac{\partial \mathcal{B}_y}{\partial y} & \frac{\partial \mathcal{B}_y}{\partial z} \end{pmatrix}^{\top} \mathbf{i} \qquad (30) $$

Besides, both $\mathcal{B}$ and $\mathcal{G}$ are a function of the location $P$ of the microrobot in the workspace $\Omega$. Through investigating the matrices $\mathcal{F}(M)$ and $\mathrm{Sk}(M)$, the conditioning of magnetic field and its gradient can be analyzed.

Since the orientation of microrobot tends to be aligned on the applied magnetic field, it is more convenient to linearize the system by specifying the desired field directly instead of magnetic torque. Hence, the magnetic field and inducing force can be obtained from:

$$ \begin{pmatrix} \mathbf{B} \\ \mathbf{f} \end{pmatrix} = \begin{pmatrix} \mathbb{I} & \mathbb{O} \\ \mathbb{O} & \mathcal{F}(M) \end{pmatrix} \begin{pmatrix} \mathcal{B}(P) \\ \mathcal{G}(P) \end{pmatrix} \mathbf{i} = \mathcal{A}_b(M, P) \mathbf{i} \qquad (31) $$

where $\mathbb{I}$ is an identity matrix, and $\mathbb{O}$ is a zero matrix in an appropriately size.

If the microrobot is made of permanent magnet, magnetic torque is linear with the current. If the microrobot is composed of a soft magnetic material and it does not reach the magnetic saturation, there is a linear relation between the magnetic moment and the magnetic field. Consequently, the magnetic torque becomes quadratic with the current [66]. Thus, the capability of the EMA system, such as DOFs of manipulation of the microrobot,

can be evaluated from the (31). Similarly to the (13), the input currents $\mathbf{i}$ can be also here computed using the pseudo-inverse of $\mathcal{A}_b(M, P)$ matrix when the matrix has full rank.

Furthermore, as the magnetic field and force are decoupled in (31), they can be first analyzed separately. As the magnetic field is linear with the current $\mathbf{i}$, the sole $\mathbf{B}$ field control relies on the properties of the matrix $\mathcal{B}(P)$. It can be easily shown that if the rank of $\mathcal{B}(P)$ is equal to the dimension of the workspace $\forall P \in \Omega$, that leads to a full control of $\mathbf{B}$. In other words, the number of current inputs $\mathbf{i}$ should be equal or more than the dimension of the workspace $\Omega$. In contrast, the force control requires more advanced investigations that are presented hereafter.

### 3.1.2. Magnetic force control

The matrix $\mathcal{F}(M)$ is a $3 \times 5$ non-square matrix, and it can be shown that its row rank is full. Thus, $\mathcal{F}(M)$ has a right Moore-Penrose inverse which is expressed by $\mathcal{F}(M)^{\dagger} = \mathcal{F}(M)^{\top}(\mathcal{F}(M)\mathcal{F}(M)^{\top})^{-1}$. From (26), the desired magnetic gradient vectors $\mathbf{g}^{\star}$ are thereby obtained using the right inverse when the magnetic force is required. This leads to:

$$ \begin{pmatrix} \frac{\partial b_x^*}{\partial x} \\ \frac{\partial b_x^*}{\partial y} \\ \frac{\partial b_x^*}{\partial z} \\ \frac{\partial b_y^*}{\partial y} \\ \frac{\partial b_y^*}{\partial z} \end{pmatrix} = \begin{pmatrix} \frac{(m_n)m_x}{m_m} & \frac{-m_x^2 m_y}{m_m} & \frac{-m_z}{m_{nz}} \\ \frac{m_y^3 + m_y m_z^2}{m_x^3 + m_x m_z^2} & \frac{m_x^3 + m_x m_z^2}{m_x^3 + m_x m_z^2} & 0 \\ \frac{m_m}{m_n} & \frac{-m_x m_z^2}{m_m} & \frac{m_x}{m_{nz}} \\ \frac{-m_x m_y^2}{m_m} & \frac{(m_n)m_y}{m_m} & \frac{-m_z}{m_{nz}} \\ \frac{-m_x m_z^2}{m_m} & \frac{(m_n)m_z}{m_m} & \frac{m_y}{m_{nz}} \end{pmatrix} $$

$$ \mathbf{g}^{\star} = \mathcal{F}(M)^{\dagger}\mathbf{f} \qquad (32) $$

where $m_n = m_x^2 + m_y^2 + m_z^2$, $m_m = (m_x^2 + m_y^2 + m_z^2)^2 - (m_x m_y)^2$, and $m_{nz} = m_n + m_z^2$. As one can see, the matrix $\mathcal{F}(M)$ is a function of the dipole moment $M = (m_x, m_y, m_z)^{\top}$ of the microrobot. Hence, the pseudo-inverse matrix $\mathcal{F}(M)^{\dagger}$ is also affected by the magnetic moment of microrobot.

If the magnetic moment $M$ of the microrobot is fixed (e.g. when a uniform static magnetic field $B_0$ is applied), only three independent magnetic gradients of $\mathbf{g}$ are required to produce an arbitrary force $\mathbf{f}$. However, the all five magnetic gradient components are required for continuous force control if $M$
is changing. Therefore, there are ten linear relations between magnetic gradients.

For instance, linear dependence between $\frac{\partial b_x}{\partial x}$ and $\frac{\partial b_y}{\partial y}$ is considered. Let us express this linear dependence through adding a coefficient $\alpha \in \mathbb{R}^*$, that is: $\frac{\partial b_y}{\partial y} = \alpha \frac{\partial b_x}{\partial x}$. The Eq. (26) is written with the following $\mathbf{g}$ matrix:

$$ \mathbf{g} = \left( \frac{\partial b_x}{\partial x} \quad \frac{\partial b_x}{\partial y} \quad \frac{\partial b_x}{\partial z} \quad \alpha \cdot \frac{\partial b_x}{\partial x} \quad \frac{\partial b_y}{\partial z} \right)^T \quad (33) $$

If the magnetic microrobot is aligned along the $z$-axis, the magnetic force in $z$-axis direction is expressed as:

$$ f_z = -(1 + \alpha) m_z \frac{\partial b_x}{\partial x} \quad (34) $$

Obviously, $f_z$ equals to zero when: $\alpha = -1$, that is for $\frac{\partial b_y}{\partial y} = -\frac{\partial b_x}{\partial x}$.

Let us now assume that the microrobot has a magnetic dipole moment defined by:

$$ M = \{(m, \beta m, 0)^{\top}, (m, 0, \beta m)^{\top}, (0, m, \beta m)^{\top}\} \quad (35) $$

with $\beta \in \mathbb{R}^*$ a coefficient term. Then, these values have to be taken into the (32) to investigate the existence of a solution. The necessary and sufficient condition for a matrix reversibility is that the determinant is not equal to 0. Since the determinant of a matrix is equal to the product of all of their eigenvalues, all eigenvalues must be not equal to 0 in order to enable the existence of the inverse matrix. Thus, the singularly values of the system for the considered linear dependence can be mathematically analyzed by computing whether $\det(\mathcal{F}(M)\mathcal{F}(M)^{\top}) = 0$. The determinants with regard to the magnetic dipole moments of $M$ are obtained as:

$$ \det(m, \beta m, 0) = m (\alpha^2 \beta^6 + \alpha^2 \beta^4 + 2\alpha\beta^4 + 2\alpha\beta^2 + \beta^2 + 1) \quad (36) $$

$$ \det(m, 0, \beta m) = m (\alpha^2 \beta^6 + \alpha^2 \beta^4 + 2\alpha\beta^4 + 4\alpha\beta^2 + 2\alpha\beta^2 + \beta^6 + 3\beta^4 + 3\beta^2 + 1) \quad (37) $$

$$ \det(0, m, \beta m) = m (\alpha^2 \beta^6 + 3\alpha^2 \beta^4 + 3\alpha^2 \beta^2 + \alpha^2 + 2\alpha\beta^6 + 4\alpha\beta^4 + 2\alpha\beta^2 + \beta^6 + \beta^4) \quad (38) $$

where $\det$ represents the determinant operator, such as $\det(M)$ is expressed as the determinant of $(\mathcal{F}(M)\mathcal{F}(M)^{\top})$ for the given dipole moment.

When the determinant becomes 0, the rank of matrix $\mathcal{F}(M)$ is less than 3 and the inverse matrix does not exist, thereby, the system becomes singular. Thus, the system should avoid the relationship between $\beta$ and $\alpha$ regarding the three dipole moments, respectively, meets the followings:

$$ \beta = \pm \frac{1}{\sqrt{\alpha}} \quad (39) $$

$$ \beta = \pm \frac{1}{\sqrt{-\alpha - 1}} \quad (40) $$

$$ \beta = \pm \frac{\sqrt{-\alpha(\alpha + 1)}}{\alpha + 1} \quad (41) $$

When the above equations are satisfied, one of magnetic force components becomes linear dependent to one of the others. For instance, for a dipole moment $M = (m, \beta m, 0)^{\top}$, the basic condition is: $\alpha > 0$ and not satisfy (39), and then a singularity can be avoided. The other singular cases are summarized in Table. (1).

<table><caption>Table 1: The singular cases caused by the linear dependent between $\frac{\partial b_x}{\partial x}$ and $\frac{\partial b_y}{\partial y}$</caption><thead><tr><th>$M^{\top}$</th><th>Relationship between $\beta$ and $\alpha$</th><th>Condition on $\alpha$</th></tr></thead><tbody><tr><td>$(m, \beta m, 0)$</td><td>$\beta = \pm \frac{1}{\sqrt{\alpha}}$</td><td>$\alpha > 0$</td></tr><tr><td>$(m, 0, \beta m)$</td><td>$\beta = \pm \frac{1}{\sqrt{-\alpha - 1}}$</td><td>$\alpha < -1$</td></tr><tr><td>$(0, m, \beta m)$</td><td>$\beta = \pm \frac{\sqrt{-\alpha(\alpha + 1)}}{\alpha + 1}$</td><td>$\alpha \neq -1 \land \alpha(\alpha + 1) \le 0$</td></tr></tbody></table>

Similarly, other nine linearly dependent relations: $\frac{\partial b_x}{\partial y} = \alpha \frac{\partial b_x}{\partial x}$, $\frac{\partial b_x}{\partial x} = \alpha \frac{\partial b_x}{\partial z}$, $\frac{\partial b_x}{\partial x} = \alpha \frac{\partial b_y}{\partial z}$, $\frac{\partial b_x}{\partial y} = \alpha \frac{\partial b_x}{\partial z}$, $\frac{\partial b_x}{\partial y} = \alpha \frac{\partial b_y}{\partial y}$, $\frac{\partial b_x}{\partial y} = \alpha \frac{\partial b_y}{\partial y}$, $\frac{\partial b_x}{\partial y} = \alpha \frac{\partial b_y}{\partial y}$, and $\frac{\partial b_y}{\partial y} = \alpha \frac{\partial b_y}{\partial z}$ have also been investigated. The different conditions to enable a linear relationship between the magnetic gradients terms of $\mathbf{g}$ are summarized hereafter. If these conditions are fulfilled, the matrix $\mathcal{F}(M)$ is not always full rank. Thereby, the magnetic force control becomes singular because the matrix $\mathcal{F}(M)$ is not full rank when any two magnetic gradients exist linear dependence. It means that magnetic force $\mathbf{f}$ cannot be controlled in some orientations and positions where non-full rank happens. To apply an arbitrary force $\mathbf{f}$ on the microrobot in any position $P \in \Omega$ and orientation requires the all five magnetic gradients, whose
terms must be controlled independently.

It is obvious that the magnetic gradient is a variation in the magnetic field B with respect to position. Thus, when the magnetic gradients are applied on a microrobot for providing propulsion force, it is synchronously placed in a global magnetic field. According to (4), if the magnetic field B is changing, a magnetic torque t is simultaneously generated on the microrobot, which will change its magnetic moment orientation and then affect the applied magnetic force. For an EMA system to continuously apply a desired force, it must be capable of controlling the microrobot's orientation, and updating the magnetic gradients faster than the change of dipole-moment direction. Therefore, the control of magnetic field should be considered when applying magnetic gradient to the microrobot.

### 3.1.3. Combined torque and force control

Torque control can be analyzed directly through the magnetic field properties, as expressed in (31). Commonly, the three magnetic fields components are required to generate magnetic torque to align microrobot to any directions in the workspace Ω. As discussed above, the control of the five independent magnetic gradients is the necessary and sufficient condition for a singular-free force control. Hence, the (31) can be used to analyze the combined torque and force control. In addition, from equation (31), the matrix B and matrix G are combined as a 6 × 8 matrix defined as:

$$C(M) = \begin{pmatrix} \mathbb{I} & \mathbb{O} \\ \mathbb{O} & \mathcal{F}(M) \end{pmatrix} \quad (42)$$

The first three columns are the expansion of the identity matrix I, thus there is no singular case. The last five columns are relevant to force control. The above matrix C(M) can be reduced to C(M) that is a 6 × 7 matrix by replacing a linear combination between one of the first three columns and one of the last five columns. The obtained matrix C(M) describes an EMA system where there is a linear dependence between the magnetic field and its gradient.

Moreover, the all 10 singular cases have been already investigated regarding linearly dependent relations of magnetic gradients. There is always a singular case if any one of such 10 linear dependencies exists. Besides, there are 3 linear dependencies between magnetic fields. As discussed, these 3 linear relations will lead to singular case of field

control. Thus, only the linear dependence between the first three columns and last five columns should be considered in the analysis of the matrix C. In fact, there are totally 28 linear dependencies in the matrix C(M).

Consequently, there are 10 dependent relationships between each magnetic gradients, and 3 linear dependencies between the magnetic fields. Since, there is a total of 28 linear dependencies in the matrix C, thereby, the linear combinations between fields and gradients are 15 that will be investigated hereafter.

Let us recall that an unrestrained magnetic microrobot will be aligned with the applied field direction for quasi-static manipulations. Specifically, we assume that the direction of magnetic moment changes slower compared to the rotational time of the alignment of the microrobot with the magnetic field B.

Once again, we will discuss linear dependencies on the magnetic gradient component $\frac{\partial b_x}{\partial x}$ as examples. These dependencies between $\frac{\partial b_x}{\partial x}$ and the magnetic field $B = (b_x, b_y, b_z)^T$ are defined with the following formulas: a) $\frac{\partial b_x}{\partial x} = \alpha b_x$, b) $\frac{\partial b_x}{\partial x} = \alpha b_y$, and c) $\frac{\partial b_x}{\partial x} = \alpha b_z$, with $\alpha \in \mathbb{R}^*$ a coefficient. The corresponding matrices can be expressed respectively as follows:

$$ \begin{align} C_a &= \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ \alpha m_x & 0 & 0 & m_y & m_z & 0 & 0 \\ 0 & 0 & 0 & m_x & 0 & m_y & m_z \\ -\alpha m_z & 0 & 0 & 0 & m_x & -m_z & m_y \end{pmatrix} \nonumber \\ C_b &= \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ \alpha m_x & 0 & 0 & m_y & m_z & 0 & 0 \\ 0 & 0 & 0 & m_x & 0 & m_y & m_z \\ -\alpha m_z & 0 & 0 & 0 & m_x & -m_z & m_y \end{pmatrix} \tag{43} \\ C_c &= \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \alpha m_x & m_y & m_z & 0 \\ 0 & 0 & 0 & m_x & 0 & m_y & m_z \\ 0 & 0 & -\alpha m_z & 0 & m_x & -m_z & m_y \end{pmatrix} \nonumber \end{align} $$

When $\frac{\partial b_x}{\partial x}$ and $b_x$ are linearly dependent, it is obvious that the matrix $C_a$ is not full rank. The magnetic force on the x-direction is directly related to the magnetic field when the dipole moment is $M = (m_x, 0, 0)^T$. That is the magnitude of magnetic field will be proportional to the magnetic force required in the dependent direction. However, when the force is applied to a direction, the orientation of the microrobot may be affected due to the
<figure><img src="image_60.png" /><figcaption>Figure 5: The illustration of the contradiction of microrobot orientation and magnetic field for generating magnetic force.</figcaption></figure>

linear dependence between the generated force and the applied field.

Moreover, there is still an issue on the sign of the force that can not be changed. Here, the magnetic field is used to align the microrobot magnetic moment along the same direction, that is: $m_x = \beta b_x$, with $\beta \in \mathbb{R}^*$ a coefficient. From the above analysis, the x-directed magnetic force is not linearly independent to the magnetic field, which is expressed as: $f_x = \alpha b_x m_x = \alpha \beta b_x^2$. Hence, the sign of $f_x$ can not be changed because the sign of dipole moment and the sign of field are always the same.

Therefore, such control system requires a non-magnetic restoring force in the x-direction. If a suitable nonmagnetic restoring force exists, the EMA system with the linear dependence between $\frac{\partial b_x}{\partial x}$ and $b_x$ can be used to control the microrobot in the desired direction and position, thus, at least 7 inputs are required.

Some other linearly dependent relationships between magnetic field and its gradient can cause physical impossibility. For instance, a contradiction exists between the orientation of the microrobot and the direction of magnetic field. Indeed, the orientation of microrobot should be aligned along direction of magnetic field as we assumed.

When $\frac{\partial b_x}{\partial x}$ and $b_y$ are linearly dependent, and if $M = (m_x, 0, 0)^T$, it is impossible for any EMA system to apply a magnetic force in the x-direction. In this case, the force is expressed by $f_x = \alpha b_y m_x$. Obviously, it is impossible to apply a magnetic force in the x-direction. The reason is that it is not feasible to align the microrobot along the x-axis when the magnetic field is applied in the y-direction. As shown in the Fig. 5, the desired magnetic force requires that the orientation of the microrobot is aligned in x-direction that demands the applied magnetic field is aligned to x-direction. However, the magnetic field is applied to y-direction due to

<figure><img src="image_61.png" /><figcaption>Figure 6: Schematic representation of a multiple electromagnets system inducing a magnetic flux $B(P)$ within the workspace $\Omega$ (blue square box). The length $l_e$ represents the distance between the coil center $O_e$ to the workspace center $O$.</figcaption></figure>

the assumed linear dependence. Hence, the desired magnetic force can not be produced in respect of such dependent relation. Similarly, when $\frac{\partial b_x}{\partial x}$ and $b_z$ are in linear dependence, and dipole moment is $M = (m_x, 0, 0)^T$, the magnetic force can be expressed as $f_x = \alpha b_z m_x$. This x-direction magnetic force can not be generated due to the constraint that the magnetic field aligns the microrobot to applied direction.

Other linearly dependent relations have also analyzed on $\frac{\partial b_x}{\partial y}$: d) $\frac{\partial b_x}{\partial y} = \alpha b_x$, e) $\frac{\partial b_x}{\partial y} = \alpha b_y$, and f) $\frac{\partial b_x}{\partial y} = \alpha b_z$; on $\frac{\partial b_x}{\partial z}$: g) $\frac{\partial b_x}{\partial z} = \alpha b_x$, h) $\frac{\partial b_x}{\partial z} = \alpha b_y$, and i) $\frac{\partial b_x}{\partial z} = \alpha b_z$; on $\frac{\partial b_y}{\partial y}$: j) $\frac{\partial b_y}{\partial y} = \alpha b_x$, k) $\frac{\partial b_y}{\partial y} = \alpha b_y$, and l) $\frac{\partial b_y}{\partial y} = \alpha b_z$; on $\frac{\partial b_y}{\partial z}$: m) $\frac{\partial b_y}{\partial z} = \alpha b_x$, n) $\frac{\partial b_y}{\partial z} = \alpha b_y$, and p) $\frac{\partial b_y}{\partial z} = \alpha b_z$. However, the possible linear relations only exist at following cases: $\frac{\partial b_x}{\partial x}$ and $\alpha b_x$; $\frac{\partial b_y}{\partial y}$ and $\alpha b_y$; or $\frac{\partial b_x}{\partial x} + \frac{\partial b_y}{\partial y}$ and $\alpha b_z$ (as $\frac{\partial b_x}{\partial x} + \frac{\partial b_y}{\partial y} + \frac{\partial b_z}{\partial z} = 0$). Such linear dependencies can obviously reduce the inputs for control system that requires less coils. It should be noted that reducing the number of coils only works under some certain conditions, for instance the non-magnetic restoring force is required. In addition, to the study of singularity under the linearly dependencies between input terms, with some advanced control technology application to the magnetic field control, the electromagnetic field singularity problems still can be investigated and avoided [67].

## 3.2. Magnetic field generated in the reference workspace

Considering the EMA system composed of a set of $n$ electromagnets, the generated magnetic field $\mathcal{F}_0 B(P)$ at any point $P$ in the workspace $\Omega$ can be illustrated in Fig. 6. This overall magnetic field is assumed to be the sum of the contributions of
all individual electromagnets. Recall that when the current input $i_e$ flows through the coils $e$, the corresponding magnetic field can be expressed as: $B_e(P) = \tilde{B}_e(P)i_e$. A key step is thus to compute the magnetic field produced by the electromagnets $e$. Several models, based either on numerical or analytical approaches, have been proposed in the literature. Numerical models are commonly based on maps of the magnetic field obtained either from FEM or from experimental measurement of the field $B_e$. The numerical method allows a good accuracy but costs a lot of calculation time. Analytical methods are often based on dipole approximation which offers a better computation time, or even on elliptical integrals [68, 69, 70]. Hybrid approaches, using a map of the magnetic field obtained from FEM and a fitting of an analytical model can be also considered, such as in [43]. The choice of the method commonly leads to make the best trade-off between speed and accuracy. As in this work, numerous simulations are realized, we assume that the magnetic field $B_e(P)$ induced by the electromagnet $e$ can be approximated by the magnetic point-dipole model. Specifically, the point-dipole model expresses the magnetic field $\mathcal{F}_e B_e(P)$ of the coils $e$ with respect to its own frame $\mathcal{F}_e(O_e : x_e, y_e, z_e)$, as shown in Fig. 6, and can be written as:

$$ \mathcal{F}_e B_e(P) = \frac{\mu_0}{4\pi|P|^3} \left( \frac{3(M_e \cdot P)}{|P|^2} - M_e \right) \quad (44) $$

where $M_e$ is the equivalent magnetic dipole moment related to the magnet source $e$ for a unit current input. Indeed, for the sake of clarity and simplicity, we assume that each electromagnet $e$ can be approximated by its analogous dipole moment $M_e$. Hence, the magnetic field is generated by using the point-dipole model for computation. The electromagnet is considered as a theoretical point-dipole mapping to the unit-current contribution, thus the shape or the filled core does not affect the performance analysis of the system. Especially, the manipulability indexes, conditioning indexes and existing singularities are based on the arrangements of electromagnet instead of input current.

The magnetic field $\mathcal{F}_e B_e(P)$ can be thus expressed in the reference frame $\mathcal{F}_0(O : x, y, z)$ linked to the workspace center using the homogeneous transformation:

$$ \mathcal{F}_0 B_e(P) = \mathcal{F}_0 \mathbf{T}_{\mathcal{F}_e} \times \mathcal{F}_e B_e(P) \quad (45) $$

where the homogeneous transformation matrix is

basically defined as:

$$ \mathcal{F}_0 \mathbf{T}_{\mathcal{F}_e} = \left( \frac{\mathcal{F}_0 \mathbf{R}_{\mathcal{F}_e}}{\mathbf{0}} \middle| \frac{\mathcal{F}_0 \mathbf{t}_{\mathcal{F}_e}}{1} \right) \quad (46) $$

where $\mathcal{F}_0 \mathbf{R}_{\mathcal{F}_e}$ and $\mathcal{F}_0 \mathbf{t}_{\mathcal{F}_e}$ denote the rotation and translation matrices with respect to the reference frame $\mathcal{F}_0$, respectively.

The overall magnetic field distribution in the workspace produced by a set of $n$ electromagnetic coils can be superposed, that is:

$$ \begin{aligned} \mathcal{F}_0 B(P) &= \sum_{e=1}^{n} \mathcal{F}_0 B_e(P) = \sum_{e=1}^{n} \mathcal{F}_0 \tilde{B}_e(P) i_e \\ &= \mathcal{F}_0 B(P) \mathbf{i} = (\mathcal{F}_0 \tilde{B}_1(P) \quad \dots \quad \mathcal{F}_0 \tilde{B}_n(P)) \mathbf{i} \end{aligned} \quad (47) $$

with the input currents $\mathbf{i} = (i_1, i_2, \dots, i_n)^T$.

It can be easily shown that the total magnetic field $\mathcal{F}_0 B(P)$ in the workspace can be changed not only thanks to the currents $\mathbf{i}$, but also by varying the position and/or orientation of electromagnets. If the electromagnets are dynamically moved following a control strategy, a similar current-control approach can be designed.

## 4. Configurations of magnetic systems with multiple electromagnets

The arrangement of electromagnets is also a key issue for the magnetic actuation system. With the number of electromagnets computed by the proposed approach and the derived equations, the simulations of various EMA systems can be performed mathematically. The magnetic characteristics will be investigated according to different configurations of electromagnets with the various performance metrics mentioned in section. 2.3. Thus, in this section, several typical configurations of EMA systems will be quantitatively evaluated.

In a 2D space, the flat configurations of EMA systems including $n = 4, 6$ and $8$ electromagnets are studied. The 3D-placement setups with $n = 6$ and $8$ electromagnets are simulated for 3D space. These configurations as illustrated in Fig.7 are selected to estimate the effectiveness of different arrangements of electromagnets.

As shown in Fig.7a-Fig.7c, the considered four-electromagnet system, the six-electromagnet system and the eight-electromagnet system are organized in axisymmetric configuration around the center $O$ of the workspace where their coils separated with an offset angle of 90°, 60° and 45°,
<figure><img src="image_62.png" /><figcaption>Figure 7: Representation of different basic multi-electromagnet EMA systems: (a) flat four-electromagnet system; (b) flat six-electromagnet system; (c) flat eight-electromagnet system; (d) 3D six-electromagnet system; (e) OctoMag system (eight coils) and (f) MiniMag system (eight coils).</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>Figure 8: 3D eight coils EMA systems arrangement: (a) top view on *xy*-plane and the (b) side view on *xz*-plane.</figcaption></figure>

respectively. The Fig.7d-Fig.7f illustrate the 3D-arrangement setups where the 3D six-electromagnet system, eight-electromagnet OctoMag system and eight-electromagnet MiniMag system are represented. The 3D six-electromagnet system consists of three pairs of opposing electromagnets aligned along the *x*, *y* and *z*-axis with same distance to the center *O*. The OctoMag and MiniMag systems designed and developed at the ETH Zurich [43] are both composed of *n* = 8 electromagnets. Here, the eight electromagnets have been divided into two sets of four coils referred as the upper and lower sets. As presented in Fig.8a, these two sets are equally spaced and organized around the common *z*-axis of rotation with $\alpha$ = 45°, and are pointing at the center *O* of the workspace with a fixed distance. For OctoMag system, the lower set is placed in *xy*-

<figure><img src="image_64.png" /><figcaption>Figure 9: The magnetic field *B* generated by OctoMag-like EMA system with $d_w$ = 65 mm: (a) a 3D view, (b) the *xy*-plane, (c) the *xz*-plane and (d) the *yz*-plane. The colorbar indicates the magnetic field magnitude.</figcaption></figure>

plane, with $\beta_e$ = 0 ($e$ = 1..4) and the upper are tilted $\beta_e$ = 45° ($e$ = 5..8). Whereas, the lower set of MiniMag is rotated to $\beta_e$ = 26° ($e$ = 1..4) from the *xy*-plane, and the upper set is set to $\beta_e$ = 47.5° ($e$ = 5..8).

In our study, the influence of the distance $d_w$ between the workspace center *O* and each electromagnet regarding the performance of EMA system is assessed. Also, the impact of tilted angle of the placed electromagnets is evaluated as well to optimize the system performances.

### 4.1. Optimization regarding working distance $d_w$ of electromagnets for EMA System

We consider that each coil has the equal distance $d_w$ for the considered system as represented in Fig. 8. Different arrangements of EMA system have been investigated from numerous simulations, and different metrics have been calculated to evaluate their performances. To sum up, the values of root mean square (RMS) and uniformity indexes of the magnetic field are synthesized in Table.2, and in Table.3 for its gradient.

Commonly, long distance $d_w$ from the workspace center *O* to the electromagnet leads globally to a slight more uniform magnetic field and gradient, in addition, the effective control of magnetic micro-robot is improved. Nevertheless, stronger magnetic
<table><caption>Table 2: Metrics of the magnetic field strength: its RMS (mT) and uniformity $\gamma(B)$ metric (%)</caption><thead><tr><th rowspan="2">EMA System</th><th colspan="5">Distance d<sub>w</sub> (mm)</th></tr><tr><th>60</th><th>65</th><th>70</th><th>75</th><th>80</th></tr></thead><tbody><tr><td>Flat four-coil (Fig. 7a)</td><td>13.13 (71.18)</td><td>9.39 (70.97)</td><td>6.91 (70.82)</td><td>5.21 (70.81)</td><td>4.00 (70.88)</td></tr><tr><td>Flat six-coil (Fig. 7b)</td><td>20.55 (72.55)</td><td>14.69 (72.96)</td><td>10.80 (73.03)</td><td>8.13 (72.96)</td><td>6.24 (72.85)</td></tr><tr><td>Flat eight-coil (Fig. 7c)</td><td>27.41 (73.05)</td><td>19.59 (73.26)</td><td>14.41 (73.23)</td><td>10.85 (73.09)</td><td>8.33 (72.92)</td></tr><tr><td>3D six-coil (Fig. 7d)</td><td>10.08 (51.09)</td><td>6.48 (50.00)</td><td>4.28 (49.11)</td><td>2.90 (48.39)</td><td>2.01 (47.80)</td></tr><tr><td>OctoMag (Fig. 7e)</td><td>24.86 (76.97)</td><td>18.61 (79.32)</td><td>14.49 (81.32)</td><td>11.58 (82.92)</td><td>9.44 (84.19)</td></tr><tr><td>MiniMag (Fig. 7f)</td><td>55.73 (63.86)</td><td>42.40 (67.20)</td><td>33.21 (69.89)</td><td>26.57 (72.12)</td><td>21.63 (74.02)</td></tr></tbody></table>

<table><caption>Table 3: Metrics of the magnetic field gradient strength: its RMS (mT/m) and uniformity $\gamma(B)$ metric (%)</caption><thead><tr><th rowspan="2">EMA System</th><th colspan="5">Distance d<sub>w</sub> (mm)</th></tr><tr><th>60</th><th>65</th><th>70</th><th>75</th><th>80</th></tr></thead><tbody><tr><td>Flat four-coil (Fig. 7a)</td><td>1.18 (59.89)</td><td>0.78 (64.58)</td><td>0.54 (68.95)</td><td>0.38 (72.99)</td><td>0.28 (76.78)</td></tr><tr><td>Flat six-coil (Fig. 7b)</td><td>1.48 (64.27)</td><td>0.97 (70.66)</td><td>0.68 (75.34)</td><td>0.50 (78.69)</td><td>0.38 (81.21)</td></tr><tr><td>Flat eight-coil (Fig. 7c)</td><td>1.87 (66.49)</td><td>1.27 (71.77)</td><td>0.91 (75.60)</td><td>0.67 (78.66)</td><td>0.51 (81.20)</td></tr><tr><td>3D six-coil (Fig. 7d)</td><td>1.32 (65.61)</td><td>0.84 (66.97)</td><td>0.55 (67.50)</td><td>0.37 (66.94)</td><td>0.26 (66.42)</td></tr><tr><td>OctoMag (Fig. 7e)</td><td>1.80 (36.68)</td><td>1.00 (43.02)</td><td>0.59 (49.00)</td><td>0.37 (54.40)</td><td>0.25 (59.26)</td></tr><tr><td>MiniMag (Fig. 7f)</td><td>2.17 (56.69)</td><td>1.31 (66.02)</td><td>0.88 (72.13)</td><td>0.63 (76.26)</td><td>0.47 (79.18)</td></tr></tbody></table>

field and gradient require a short length $d_w$. Thus, a compromise should be made with respect to the biomedical applications specifications. Certainly, the value of $d_w$ will affect the size of the workspace since a great length enables larger workspace dimension.

Obviously, EMA setup with $n = 8$ coils produces stronger magnetic field and gradient than configuration with fewer magnets. Furthermore, as presented in the results of tables, MiniMag arrangement generates the strongest and the most uniform magnetic field, specially in *xy*-plane. From these results, it confirms that the performances of EMA systems are not only affected by the number of coils but also by their configurations.

In order to better present the effect of the shifted distance $d_w$ of electromagnet, let us focus on the EMA system with a certain number of coils. As mentioned in section.3.1, the $n = 8$ electromagnets can be applied for the combined force and torque singular-free control with 5 DOFs (3-DOF

position and 2-DOF pointing orientation) in a 3D workspace. In the following, we only keep one variable parameter. We have chosen to set the tilted angle similar to the OctoMag EMA setup [43], while working distances of $d_w = 60$ mm, 70 mm and 80 mm are evaluated, respectively.

As presented in Fig. 9, the magnetic field $B$ distribution is mainly oriented along the $z$-axis direction. As the OctoMag setup shares the same arrangement of electromagnets in the *xy*-plane with the flat four-electromagnet configuration, the vector field in the *xy*-plane is quite homogeneous. The directions of the magnetic field vectors mainly along the $z$-axis are basically due to the upper set of electromagnets.

Fig. 10 shows the magnetic field magnitude $||B||$ in the *xy*-plane for the different working distances. The magnetic field distribution $B$ looks much less like a radial field and becomes stronger for the shorter distance $d_w$. However, if a magnetic micro-robot moves outward from the center of *xy*-plane,
<figure><img src="image_65.png" /><figcaption>Figure 10: The magnetic field magnitude $||B||$ in the xy-plane for lengths of (a) $d_w$ = 60 mm, and (b) $d_w$ = 80 mm.</figcaption></figure>

<figure><img src="image_66.png" /><figcaption>Figure 11: Magnetic field metrics over the workspace $\Omega$ for $d_w$ ranging from 60 mm to 80 mm: (a) error-bar showing the average and STD values; and (b) the uniformity index $\gamma$.</figcaption></figure>

the longer distance $d_w$ can make magnetic field more homogeneous in each direction. This means that OctoMag system is able to actuate a magnetic microrobot to any direction with less interference in long distances $d_w$. For instance, when the working distance is set to $d_w$ = 65 mm, the maximum magnetic intensity decreases from 36.99 mT with the 2D configuration to 35.86 mT with the OctoMag platform.

The magnetic field metrics are shown in Fig. 11. It appears that the average value $\langle||B||\rangle$ of the magnetic field magnitude is equivalent to the flat eight-

<figure><img src="image_67.png" /><figcaption>Figure 12: Magnetic field gradient metrics for distance $d_w$ ranging from 60 mm to 80 mm: (a) the mean and STD values; and (b) the uniformity index $\gamma$.</figcaption></figure>

<figure><img src="image_68.png" /><figcaption>Figure 13: Performance metrics of the force actuation matrix $A_f$: (a) the manipulability $w_n$ and (d) conditioning number $1/\kappa$ indexes for $P \in \Omega$ when $M$ is aligned along the x-direction; (b) $\Gamma_{w_n}(A_f)$ and (e) $\Gamma_{1/\kappa}(A_f)$ for the sampled orientation of $M$; and (c) $\Gamma_{w_n}(A_f)$ and (f) $\Gamma_{1/\kappa}(A_f)$ as function of $d_w$.</figcaption></figure>

electromagnet arrangement, however, its STD value is here improved. Furthermore, the uniformity index $\gamma(||B||)$ is more important and is increasing with the length $d_w$. Therefore, the values confirm that OctoMag-like arrangement is a promising solution to efficiently actuate a magnetic microrobot in a 3D workspace.

<figure><img src="image_69.png" /><figcaption>Figure 14: Performance metrics of the torque actuation matrix $A_t$: (a) the manipulability $w_n$ and (d) conditioning number $1/\kappa$ indexes in $P \in \Omega$ when $M$ is aligned along the x-direction; (b) $\Gamma_{w_n}(A_t)$ and (e) $\Gamma_{1/\kappa}(A_t)$ for the sampled orientation of $M$; and (c) $\Gamma_{w_n}(A_t)$ and (f) $\Gamma_{1/\kappa}(A_t)$ as function of $d_w$.</figcaption></figure>

The metrics of the magnetic field gradient are illustrated in Fig. 12. The average of magnetic field gradient behaves equivalently to flat eight-electromagnet configuration. However, as the mag-
netic field becomes more uniform, its gradient is ob-
viously lower. Next, it appears that the uniformity
of the magnetic gradient is decreasing significantly
with $d_w$, especially, for the component $\frac{\partial b_y}{\partial y}$ and $\frac{\partial b_z}{\partial z}$.

The performance metrics of the force and torque
actuation matrices in the workspace for the sampled
orientations of dipole moment *M* are presented in
Fig. 13 and 14, respectively. The manipulability
indexes *w*<sub>n</sub> of the force and torque have better per-
formance when the orientation of the microrobot
magnetic moment *M* is in the *xy*-plane (θ = 90°)
and along *z*-axis (θ = 0°). The conditioning in-
dexes 1/κ of force and torque reach a high value
when the dipole moment *M* is aligned along the *z*-
axis. From these results, it can be seen that these
global performance metrics of the both force and
torque control are increasing with *d*<sub>w</sub> in the given
OctoMag configuration. Furthermore, the impact
of changing the shifting angles of EMA system will
be estimated hereafter.

4.2. Analysis with respect to the tilted angle $\beta$ of electromagnets for EMA systems

Let us highlight the orientation angle $\beta$ of elec-
tromagnets and fix other parameters. Once again,
eight electromagnets are applied to enable reli-
able combined torque and force singular-free con-
trol with 5 DOFs. As it has been shown above that
long distance $d_w$ from the workspace center $O$ to the
electromagnet leads globally to a slight more uni-
form magnetic field and gradient, in addition, the
effectiveness control of magnetic microrobot is im-
proved. Nevertheless, stronger magnetic field and
gradient require a short length $d_w$. With a compre-
hensive consideration of strength and uniformity,
we have set the working distance to $d_w$ = 65 mm,
that is similar to the OctoMag EMA setup [43],
and a workspace of $\Omega$ = 45 mm × 45 mm × 45 mm
is considered throughout the simulations. For the
sake of simplicity, the eight electromagnetic coils of
the EMA platform are divided in two sets: four sta-
tionary electromagnets: $e$ = 1..4; and four mobile
coils: $e$ = 5..8, as illustrated in Fig. 8. These two
sets are arranged around a common axis of rotation
with an azimuth angle $\alpha_e$ = 45° ($e$ = 1..8), and are
pointing to the common center $O$ of the workspace.
The mobile sets are able to rotate their polar an-
gle $\beta_e \in [0; 90°)$ ($e$ = 5..8). Two cases are consid-
ered for the stationary coils set where the polar $\beta_e$
($e$ = 1..4) is fixed to:

1. $\beta_e = 0^\circ$, that is equivalent to the OctoMag [43] configuration;

2. $\beta_e = 26^\circ$, that is close to the MiniMag [44] arrangement.

Hence, the considered reconfigurable EMA plat-
form can vary from one to another of these config-
urations by changing the orientation of the polar
angle $\beta_e$ in their two sets of electromagnets.

In the following sections, the influence of the mo-
bile angle $\beta_e$ on the magnetic field and gradient
performance indexes, and on the magnetic actua-
tion indexes are investigated. As previously stated,
each electromagnet is computed with a point dipole
moment magnitude of: $||M_e|| = 8.178 \text{ A m}^2$, and
the microrobot is modeled from its magnetic mo-
ment $M = (1, \theta, \varphi)$, with a unit magnitude, az-
imuth $\varphi \in [0; 360^\circ)$ and polar angles $\theta \in [0; 180^\circ]$.
Let us recall that the possible collision of the coils
is not taken into account in this study.

4.2.1. Case #1: reconfigurable OctoMag-like setup

This section investigates a reconfigurable EMA
platform similar to the OctoMag system. Specifi-
cally, the mobile set of coils is able to rotate with the
same angle $\beta = \beta_e \in [0; 90^\circ)$ ($e = 5..8$). To analyze
the performance of such reconfigurable EMA setup,
different simulations are conducted, and their re-
sults are presented hereafter.

<figure><img src="image_70.png" /><figcaption>Figure 15: Performance metrics of (a)-(b) the magnetic field and (c)-(d) its gradient of reconfigurable OctoMag-like setup when the mobile coils set is rotating with $\beta \in [0; 90^\circ)$. The markers and envelopes (a)-(c) refer to the mean and standard deviation (STD) of the fields.</figcaption></figure>
<figure><img src="image_71.png" /><figcaption>Figure 16: Actuation performance indexes of the (a)-(b) force and (c)-(d) torque actuation matrices for the sampled orientations of the magnetic moment of the microrobot: (a)-(c) the global normalized manipulability index $\Gamma_{w_n}$; and (b)-(d) the global conditioning index $\Gamma_{1/\kappa}$.</figcaption></figure>

First, the magnetic field performance indexes are analyzed and shown in Fig. 15, for different moving angles $\beta \in [0; 90^\circ)$. As previously analyzed, the performance indexes are axisymmetric around the z-axis, and the magnetic field and gradient behave similarly between the $x$ and $y$ components. In particular, the average values $\langle b_x \rangle$ and $\langle b_y \rangle$ of $B$, and $\langle \frac{\partial b_x}{\partial y} \rangle$, $\langle \frac{\partial b_y}{\partial z} \rangle$ and $\langle \frac{\partial b_z}{\partial z} \rangle$ of $\nabla B$ are close to zero. Moreover, from these results, two distinguishing behaviors appear: i) for low angle $\beta \lesssim 45^\circ$, the magnetic field gradient is the strongest and more uniform; whereas ii) for higher value of $\beta$, the magnetic field becomes stronger and more uniform. Thus, the reconfigurable OctoMag-like platform exhibits a versatility to manipulate either the magnetic field or its gradient, which is only driven by simply regulating the orientation $\beta_e$ of some electromagnets.

<table><caption>Table 4: Minimum and Maximum values of the global actuation performance indexes.</caption><thead><tr><th></th><th>min</th><th>(&phi;, &theta;, &beta;)</th><th>max</th><th>(&phi;, &theta;, &beta;)</th></tr></thead><tbody><tr><td>&Gamma;<sub>w<sub>n</sub></sub>(&#x211D;<sub>f</sub>)</td><td>0.241</td><td>(0, 0, 89)</td><td>0.772</td><td>(0, 0, 22)</td></tr><tr><td>&Gamma;<sub>1/&kappa;</sub>(&#x211D;<sub>f</sub>)</td><td>0.133</td><td>(180, 180, 89)</td><td>0.739</td><td>(0, 0, 32)</td></tr><tr><td>&Gamma;<sub>w<sub>n</sub></sub>(&#x211D;<sub>t</sub>)</td><td>0.446</td><td>(90, 90, 0)</td><td>0.889</td><td>(0, 0, 0)</td></tr><tr><td>&Gamma;<sub>1/&kappa;</sub>(&#x211D;<sub>t</sub>)</td><td>0.227</td><td>(210, 90, 0)</td><td>0.973</td><td>(0, 0, 68)</td></tr></tbody></table>

Fig. 16 shows the global performance indexes for the sampled orientations of the magnetic moment $M$ with different moving $\beta$ angles. The simulation results release that if $M$ is parallel to the $xy$-plane (ie. $\theta = 90^\circ$), the global performance indexes of the torque $\mathbf{t}$ are low, whereas the force $\mathbf{f}$ is low when $M$ is aligned along the $z$-axis ($\theta = 0^\circ$ or $180^\circ$). To further investigate these results, table.4 reports the maximal and minimal values of the global performance indexes, and the corresponding angles for the sampled orientations of a unit-strength magnetic moment. It can be shown that to improve the performance of magnetic actuation, firstly, $M$ should be globally aligned along the $z$-axis. Secondly, the moving angle must be set to $\beta = 32^\circ$ in order to maximize the dexterity index $1/\kappa$ of the force, whereas a value around $\beta = 68^\circ$ is required for the torque. In contrast, the normalized manipulability of the force and torque needs a low angle below $\beta \le 22^\circ$ for improvement.

<figure><img src="image_72.png" /><figcaption>Figure 17: Statistical data of the global performance indexes of the force and torque actuation matrices for $\beta \in [0; 90^\circ)$: (a) the mean and STD, and (b) the uniformity $\gamma$ metrics.</figcaption></figure>

To exhibit the influence of the moving angle $\beta$ on the magnetic actuation performance, Fig. 17 depicts the mean and uniformity metrics of the global performance indexes $\Gamma_{w_n}$ and $\Gamma_{1/\kappa}$. The analysis of the performance indexes reveals that it is more efficient to control $\mathbf{f}$ with low angle, whereas $\mathbf{t}$ requires a higher $\beta$, especially to increase its uniformity. These simulation results demonstrate that rotating the moving angle $\beta$ could make the EMA platform more flexible to control efficiently either the force $\mathbf{f}$ or the torque $\mathbf{t}$.

Finally, the impact of the moving angle $\beta$ on the workspace size has been also investigated. Since the indexes are axisymmetric around the $z$-axis, the analysis is reported only along this direction. In particular, to evaluate the tissue penetration, the locations behind the workspace $\Omega$ could be considered. The different performance metrics are shown
<figure><img src="image_73.png" /><figcaption>Figure 18: Performance indexes along the z-axis of (a)-(b) the torque and (c)-(d) the force actuation matrices: (a)-(c) the normalized manipulability index $w_n$; and (b)-(d) the dexterity index $1/\kappa$. The black line represents the maximum value of the performance indexes.</figcaption></figure>

in Fig. 18. Fig. 18b shows that high angle $\beta$ leads to the reliable torque dexterity. Specifically, a high angle $\beta \gtrsim 50^\circ$ allows the increasing of the value of the conditioning number $1/\kappa$. The good condition- ing number of torque means that a reconfigurable EMA system is able to transmit a torque **t** along any directions more efficient. However, to enable sufficient force **f** in the high locations together with a good manipulability index $w_n$, a low angle $\beta \lesssim 30^\circ$ is required. Therefore, these results exhibit that to design a versatile EMA system with both effective force and torque control, it is necessary to be able to vary the moving angle $\beta_e$ of some coils. Such design objective seems a promising way to achieve optimal control of **f** and **t**.

### 4.2.2. Case #2: reconfigurable MiniMag-like setup

The second case corresponds to EMA platform, similar to the MiniMag system. Specifically, the mobile set of coils is able to rotate with the same angle $\beta = \beta_e \in [0; 90^\circ)$ ($e = 5..8$, as shown in Fig. 8), while the stationary set of coils is fixed to $\beta_e = 26^\circ$ ($e = 1..4$).

Fig. 19 shows the performance metrics of the magnetic field and its gradient with a moving angle $\beta \in [0; 90^\circ)$. Similarly, the high moving angle value $\beta \gtrsim 45^\circ$ provides stronger and more uniform

<figure><img src="image_74.png" /><figcaption>Figure 19: Performance metrics of (a)-(b) the magnetic field and (c)-(d) its gradient of reconfigurable MiniMag-like setup when the mobile coils set is rotating with $\beta \in [0; 90^\circ)$. The markers and envelopes (a)-(c) refer to the mean and standard deviation (STD) of the fields.</figcaption></figure>

magnetic field **B**. Especially, the reconfigurable MiniMag-like setup induces a larger magnetic field strength along the z-axis than the OctoMag-like one (cf. Fig. 15). However, the magnetic field gradi- ent requires a low moving angle $\beta$ to be stronger and more uniform. In particular, it can be seen in Fig. 19d that their uniformity metrics $\gamma$ are decreas- ing for $\beta < 26^\circ$, and increasing again for $\beta > 56^\circ$. In such case #2, it appears that it is more conve- nient to use high angles around $\beta \approx 50^\circ-60^\circ$ to get a suitable uniform magnetic field and gradient.

The performance indexes of the actuation ma- trix $\mathcal{A}(P, M)$ are evaluated for different moving $\beta$ angles, and reported in Fig. 20. From these sim- ulation results, it appears that if $M$ is parallel to the $xy$-plane (ie. $\theta = 90^\circ$), the global conditioning indexes $\Gamma_{1/\kappa}$ of both the force and torque are low, whereas the global normalized manipulability in- dexes $\Gamma_{w_n}$ are advantageous. Moreover, when $M$ is aligned along the z-axis ($\theta = 0^\circ$ mod $180^\circ$), a good $\Gamma_{w_n}$ is obtained and $\Gamma_{1/\kappa}$ has the greatest value. Next, the moving angle $\beta$ should have a low value to enable the good global performance indexes $\Gamma_{w_n}$ of the torque and the force, while their $\Gamma_{1/\kappa}$ require a value around $\beta = 60^\circ$. These trends are clearly re- vealed in Fig. 21 that shows the magnetic actuation performance indexes statistical data and uniformity
<figure><img src="image_75.png" /><figcaption>(a) Performance metrics of (b) Performance metrics of the force $\mathcal{A}_f$ the force $\mathcal{A}_f$</figcaption></figure>

<figure><img src="image_76.png" /><figcaption>(c) Performance metrics of (d) Performance metrics of the torque $\mathcal{A}_t$ the torque $\mathcal{A}_t$</figcaption></figure>

<caption>Figure 20: Actuation performance indexes of the (a)-(b) force and (c)-(d) torque actuation matrices for the sampled orientations of the magnetic moment of the microrobot: (a)-(c) the global normalized manipulability index $\Gamma w_n$; and (b)-(d) the global conditioning index $\Gamma_{1/\kappa}$.</caption>

$\gamma$.

The influence of $\beta \in [0; 90^\circ)$ along the z-axis has been evaluated and reported in Fig. 22. To maximize the performance indexes, the moving angle should be adjusted with respect to the z-depth position of the magnetic moment $M$ of the microrobot. For instance, Fig. 22d shows that there are two angle routes to provide a force with a good dexterity: one with $\beta \lesssim 35^\circ$ and the second with $\beta \gtrsim 65^\circ$. Basically, for deep location, a low angle $\beta$ enables a reliable manipulability $w_n$ and conditioning number $1/\kappa$ metrics of both the force and torque. Therefore, the tilted angle of the applied electromagnet can be determined referring the above analysis for different applications.

<figure><img src="image_77.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>(c)</figcaption></figure>

<caption>(d)</caption>

<caption>Figure 22: Performance indexes along the z-axis of (a)-(b) torque and (c)-(d) the force actuation matrices: (c)-(a) the normalized manipulability index $w_n$; and (d)-(b) the dexterity index $1/\kappa$. The black line represents the maximum value of the performance indexes.</caption>

### 4.3. Discussions

Based on the numerous simulations, it appears that MiniMag-like arrangement provides the greatest magnetic field and gradient, besides, the most effective actuation performance can also be performed among the various considered configurations. However, the magnetic field distribution uniformity index Iso of MiniMag-like arrangement is less interesting than OctoMag-like arrangement, implying some difficulties to derive the control strategy. The system is capable of providing different performances by setting the positions and orientations of electromagnets. Therefore, the choice of the basic EMA configuration can be motivated based on the applications objectives and the desired manipulation tasks.

<figure><img src="image_79.png" /><figcaption>Figure 21: Statistical data of the global performance indexes of the force and torque actuation matrices for $\beta \in [0; 90^\circ)$: (a) the mean and STD, and (b) the uniformity $\gamma$ metrics.</figcaption></figure>

## 5. Design analysis

Different configurations of electromagnetic coils can produce various magnetic field distributions. The flat multi-electromagnet EMA systems are
used for 2D manipulation. For the 3D control of microrobot, the 3D configurations are applied. In this study, the six different typical systems consisting of flat four-electromagnet, flat six-electromagnet, flat eight-electromagnet, 3D six-electromagnet, OctoMag and MiniMag configurations were simulated and compared. Based on these results, an optimal configuration can be determined for the considered application.

Obviously, the investigated EMA systems all can be reconfigurable where either the distance $d_w$ or mobile angle $\beta$ is adjusted. If either the strong magnetic field or gradient is required in 2D plane, $n = 2$ and 3 coils are necessary, respectively. The flat four-electromagnet system with short distance $d_w$ can be applied for this control with redundancy. Commonly, if 2D applications require only a weak magnetic field and its gradient, a minimum of $n = 5$ electromagnets is essential. Thus, the flat six-electromagnet or flat eight-electromagnet systems set to long distance $d_w$ can be used to improve the field distribution and the control redundancy.

Similarly, the 3D EMA systems are also formulated with different distances $d_w$ to achieve the desired performances. Especially, the efficiency of OctoMag-like and MiniMag-like systems has been further investigated with respect to the mobile angle $\beta$. When a strong magnetic field and gradient are of prime importance for the application, $d_w$ has to be set to a short distance. Moreover, to enhance the magnetic field B, the tilted angle $\beta$ has to be set to an high value. In contrast, to enable a proper magnetic field gradient $\Delta B$, different specific values of $\beta$ can be considered. Conversely, the long distance $d_w$ leads to a more uniform magnetic field and gradient, higher manipulability of torque and force, and higher dexterity of torque and force. Moreover, these performance's improvements can also be realized by implementing the corresponding tilted angles $\beta$ to superimpose their advantages with $d_w$.

From these results, various system performances can be obtained, depending on the given design parameters. All these capabilities are summarized in the chart<fnref n="1" /> given in Fig. 23 and 24. In Fig.23, it is observed that the more electromagnets lead to the stronger magnetic field strength $\langle||B||\rangle$. However, the uniformity of field does not change significant by applying more coils. If the magnetic gradient

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

<figure><img src="image_81.png" /><figcaption>Figure 23: The optimal performances for each EMA system design for (a) short and (b) long distance $d_w$.</figcaption></figure>

is prioritized for control, MiniMag, OctoMag and flat eight-electromagnet configurations are set to short working distance $d_w$. Meanwhile, MiniMag, flat eight-electromagnet and 3D six-electromagnet configurations are proposed for long distance $d_w$. Considering the uniformity of magnetic field and gradient, MiniMag system has the capability of generating both uniform field and gradient. OctoMag system can produce more uniform magnetic field as well, while its weakness lies in the uniformity of magnetic gradient. Besides, the flat configurations are also able to provide uniform magnetic gradient. The 3D six-electromagnet and flat four-electromagnet configurations carry the benefit of the good manipulability of torque. Especially in long distance $d_w$, flat four-electromagnet arrangement has significant capability of producing better manipulability of force than other configurations. With short distance $d_w$, OctoMag system can be used to provide good manipulability of the force. Whereas, the great dexterity of force is produced by using MiniMag system since the eight coils with such configuration could lead to more dexterous

<footnote n="1">¹For the sake of consistency between the different results, the magnetic field strength and gradient are here normalized.</footnote>
force control. Such improved phenomenon also occurs in the term of the dexterity of torque through applying more coils.

<figure><img src="image_82.png" /><figcaption>Figure 24: The optimal performances for the OctoMag-like and MiniMag-like designs regrading the mobile angle $\beta$ value.</figcaption></figure>

The investigation of mobile angle $\beta$ is further represented in Fig. 24. It can be observed that the magnetic field generated by MiniMag-like arrangement becomes much more strong than OctoMag-like configuration when the $\beta$ is set to a high angle value. Conversely, the lower mobile angle allows OctoMag-like system to produce a stronger and more uniform magnetic gradient $\Delta B$. Comparing in Fig. 24, MiniMag-like configuration commonly presents better performance of force and torque control. When $\beta$ is set to around $45^\circ$, MiniMag-like configuration has the significant capability to provide more manipulabilities of force and torque. When $\beta$ is rotated to higher angle, OctoMag-like systems is able to achieve high manipulability of torque control, whereas MiniMag is good at controlling force. Therefore, for the desired application and required system performance, the optimal configuration or the design inspiration could be obtained from them.

## 6. CONCLUSION

In the paper, theoretical foundation of electromagnetism has been studied, especially, the point-dipole model was analyzed for magnetic field simulation. Moreover, the electromagnetic manipulation of the untethered microrobot was investigated where the motion includes translation and rotation actuated by magnetic force and torque, respectively. Furthermore, the required minimum electromagnetic coils for EMA system were estimated

mathematically. Indeed, the linear dependencies regarding the applied field and gradient, and their resulting singular cases have been fully evaluated.

The various configurations of EMA systems have been studied and compared. The performance evaluation of the reconfigurable platform was thereby investigated. Several simulations were proceeded, and performance metrics of the actuated system were analyzed under different shifted distance $d_w$ and the tilted angle $\beta$ of electromagnets. Results show that the reconfigurable platform of electromagnets enables a variety of local magnetic field distribution. The evaluations of the force, torque, manipulability and the dexterity indexes demonstrate that the reconfigurable system provides more flexibility. Overall, the shorter distance $d_w$ can be used to generate the strong magnetic field and gradient, while their uniformities require longer distance $d_w$. Moreover, the low angle $\beta$ leads to the more efficient magnetic force control. But the magnetic torque becomes less controllable as it requires an higher angle $\beta$. Thus, when the configuration of EMA system is stationary, the system design parameters can be determined depend on the above results. If the stationary EMA platform is upgraded to a reconfigurable system by enabling the control of the distance $d_w$ and angle $\beta$, it will be resourceful to perform various tasks. Basically, for different applications, by using the mathematical approach as we demonstrated for performant electromagnetic micromanipulation platforms, the obtained EMA system could be optimized and improved concerning the specific medical application. These would help to develop more advanced navigation control strategy of biomedical magnetic microrobot for different micromanipulation tasks. In particular, it can be used for improving the convenience of minimally invasive operation.

## References

[1] F. Tendick, S. S. Sastry, R. S. Fearing, M. Cohn, Applications of micromechatronics in minimally invasive surgery, IEEE/ASME Trans. Mechatronics 3 (1) (1998) 34–42.

[2] M. J. Mack, Minimally invasive and robotic surgery, Jama 285 (5) (2001) 568–572.

[3] S. Purkayastha, T. Athanasiou, R. Casula, S. A. Darzi, Robotic surgery: a review, Hospital Medicine 65 (3) (2004) 153–159.

[4] J. Troccaz, R. Bogue, The development of medical microrobots: a review of progress, Industrial Robot: An International Journal (2008).

[5] S. Tsuda, D. Oleynikov, J. Gould, D. Azagury, B. Sandler, M. Hutter, S. Ross, E. Haas, F. Brody, R. Satava,
Sages tavac safety and effectiveness analysis: da vinci®
surgical system (intuitive surgical, sunnyvale, ca), Sur-
gical endoscopy 29 (10) (2015) 2873-2884.

[6] S.-J. Baek, S.-H. Kim, Robotics in general surgery: An evidence-based review, Asian journal of endoscopic surgery 7 (2) (2014) 117–123.

[7] L. Morelli, S. Guadagni, G. Di Franco, M. Palmeri, G. Di Candio, F. Mosca, Da vinci single site® surgical platform in clinical practice: A systematic review, The International Journal of Medical Robotics and Computer Assisted Surgery 12 (4) (2016) 724–734.

[8] J. P. Desai, R. V. Patel, A. Ferreira, S. K. Agrawal, Encyclopedia of Medical Robotics, World Scientific, 2017. 1140

[9] G.-Z. Yang, J. Bellingham, P. E. Dupont, P. Fischer, L. Floridi, R. Full, N. Jacobstein, V. Kumar, M. McNutt, R. Merrifield, et al., The grand challenges of science robotics, Science robotics 3 (14) (2018).

[10] N. Simaan, R. M. Yasin, L. Wang, Medical technologies and challenges of robot-assisted minimally invasive intervention and diagnostics, Annual Review of Control, Robotics, and Autonomous Systems 1 (2018) 465–490.

[11] C. M. Hoang, T. K. Nguyen, V. H. Le, J. Kim, C. S. Kim, Independent electromagnetic field control for practical approach to actively locomotive wireless capsule endoscope, IEEE Transactions on Systems, Man, and Cybernetics: Systems (2019) 1–13.

[12] D. Rattner, A. Kalloo, Asge/sages working group on natural orifice transluminal endoscopic surgery, Surgical Endoscopy and Other Interventional Techniques 20 (2) (2006) 329–333.

[13] J. R. Romanelli, D. B. Earle, Single-port laparoscopic surgery: an overview, Surgical endoscopy 23 (7) (2009) 1419–1427.

[14] M. M. Tiwari, J. F. Reynoso, A. C. Lehman, A. W. Tsang, S. M. Farritor, D. Oleynikov, In vivo miniature robots for natural orifice surgery: State of the art and future perspectives, World journal of gastrointestinal surgery 2 (6) (2010) 217.

[15] D. Samarasekera, J. H. Kaouk, Robotic single port surgery: Current status and future considerations, Indian journal of urology: IJU: journal of the Urological Society of India 30 (3) (2014) 326.

[16] S. U. Bae, W. K. Jeong, S. K. Baek, Current status of robotic single-port colonic surgery, The International Journal of Medical Robotics and Computer Assisted Surgery 13 (1) (2017) e1735.

[17] F. Carpi, C. Pappone, Stereotaxis niobe® magnetic navigation system for endocardial catheter ablation and gastrointestinal capsule endoscopy, Expert review of medical devices 6 (5) (2009) 487–498.

[18] B. J. Nelson, I. K. Kaliakatsos, J. J. Abbott, Micro-robots for minimally invasive medicine, Annual review of biomedical engineering 12 (2010) 55–85.

[19] M. Sitti, Voyage of the microrobots, Nature 458 (7242) (2009) 1121–1122.

[20] M. Sitti, H. Ceylan, W. Hu, J. Giltinan, M. Turan, S. Yim, E. Diller, Biomedical applications of untethered mobile milli/microrobots, Proceedings of the IEEE 103 (2) (2015) 205–224.

[21] H. Ceylan, J. Giltinan, K. Kozielski, M. Sitti, Mobile microrobots for bioengineering applications, Lab on a Chip 17 (10) (2017) 1705–1724.

[22] T. Honda, K. Arai, K. Ishiyama, Micro swimming mechanisms propelled by external magnetic fields, IEEE Transactions on Magnetics 32 (5) (1996) 5085–

5087.

[23] Y. Alapan, O. Yasa, O. Schauer, J. Giltinan, A. F. Tabak, V. Sourjik, M. Sitti, Soft erythrocyte-based bacterial microswimmers for cargo delivery, Science Robotics 3 (17) (2018) eaar4423.

[24] A. Barbot, D. Decanini, G. Hwang, The rotation of microrobot simplifies 3d control inside microchannels, Scientific reports 8 (1) (2018) 1–9.

[25] R. W. Carlsen, M. R. Edwards, J. Zhuang, C. Pacoret, M. Sitti, Magnetic steering control of multi-cellular bio-hybrid microswimmers, Lab on a Chip 14 (19) (2014) 3850–3859.

[26] F. Carpi, N. Kastelein, M. Talcott, C. Pappone, Magnetically controllable gastrointestinal steering of video capsules, IEEE Transactions on Biomedical Engineering 58 (2) (2010) 231–234.

[27] E. Diller, M. Sitti, et al., Micro-scale mobile robotics, Foundations and Trends® in Robotics 2 (3) (2013) 143–259.

[28] A. Hosney, J. Abdalla, I. S. Amin, N. Hamdi, I. S. Khalil, In vitro validation of clearing clogged vessels using microrobots, in: 2016 6th IEEE International Conference on Biomedical Robotics and Biomechatronics (BioRob), IEEE, 2016, pp. 272–277.

[29] S. Jeon, S. Kim, S. Ha, S. Lee, E. Kim, S. Y. Kim, S. H. Park, J. H. Jeon, S. W. Kim, C. Moon, et al., Magnetically actuated microrobots as a platform for stem cell transplantation, Science Robotics 4 (30) (2019) eaav4317.

[30] V. H. Le, S. Zheng, J. Han, J.-O. Park, et al., Preparation of tumor targeting cell-based microrobots carrying nir light sensitive therapeutics manipulated by electromagnetic actuating system and chemotaxis, Journal of Micro-Bio Robotics 14 (3-4) (2018) 69–77.

[31] J. Li, B. E.-F. de Ávila, W. Gao, L. Zhang, J. Wang, Micro/nanorobots for biomedicine: Delivery, surgery, sensing, and detoxification, Science Robotics 2 (4) (2017).

[32] J. Li, X. Li, T. Luo, R. Wang, C. Liu, S. Chen, D. Li, J. Yue, S.-h. Cheng, D. Sun, Development of a magnetic microrobot for carrying and delivering targeted cells, Science Robotics 3 (19) (2018) eaat8829.

[33] S. Martel, Microrobotics in the vascular network: present status and next challenges, Journal of Micro-Bio Robotics 8 (1) (2013) 41–52.

[34] J. J. Abbott, E. Diller, A. J. Petruska, Magnetic methods in robotics, Annual Review of Control, Robotics, and Autonomous Systems 3 (2015).

[35] F. Carpi, N. Kastelein, M. Talcott, C. Pappone, Magnetically controllable gastrointestinal steering of video capsules, IEEE Trans. Biomed. Eng. 58 (2) (2011) 231–234.

[36] S. Park, K. Cha, J. Park, Development of biomedical microrobot for intravascular therapy, International Journal of Advanced Robotic Systems 7 (1) (2010) 1.

[37] K. B. Yesin, K. Vollmers, B. J. Nelson, Modeling and control of untethered biomicrorobots in a fluidic environment using electromagnetic fields, The International Journal of Robotics Research 25 (5-6) (2006) 527–536.

[38] T. W. Fountain, P. V. Kailat, J. J. Abbott, Wireless control of magnetic helical microrobots using a rotating-permanent-magnet manipulator, in: 2010 IEEE International Conference on Robotics and Automation, IEEE, 2010, pp. 576–581.

[39] A. W. Mahoney, J. J. Abbott, Five-degree-of-freedom
manipulation of an untethered magnetic device in fluid using a single permanent magnet with application in stomach capsule endoscopy, The International Journal of Robotics Research 35 (1-3) (2016) 129-147.

[40] P. Ryan, E. Diller, Magnetic actuation for full dexterity microrobotic control using rotating permanent magnets, IEEE Transactions on Robotics 33 (6) (2017) 1398-1409.

[41] W. Amokrane, K. Belharet, M. Souissi, A. B. Grayeli, A. Ferreira, Macro-micromanipulation platform for inner ear drug delivery, Robotics and Autonomous Systems 107 (2018) 10-19.

[42] S. Schuerle, S. Erni, M. Flink, B. E. Kratochvil, B. J. Nelson, Three-dimensional magnetic manipulation of micro-and nanostructures for applications in life sciences, IEEE transactions on magnetics 49 (1) (2012) 321-330.

[43] M. P. Kummer, J. J. Abbott, B. E. Kratochvil, R. Borer, A. Sengul, B. J. Nelson, Octomag: An electromagnetic system for 5-dof wireless micromanipulation, IEEE Transactions on Robotics 26 (6) (2010) 1006-1017.

[44] B. E. Kratochvil, M. P. Kummer, S. Erni, R. Borer, D. R. Frutiger, S. Schürle, B. J. Nelson, Minimag: a hemispherical electromagnetic system for 5-dof wireless micromanipulation, in: Experimental Robotics, Springer, 2014, pp. 317-329.

[45] S. Jeon, G. Jang, H. Choi, S. Park, Magnetic navigation system with gradient and uniform saddle coils for the wireless manipulation of micro-robots in human blood vessels, IEEE transactions on magnetics 46 (6) (2010) 1943-1946.

[46] H. Choi, K. Cha, J. Choi, S. Jeong, S. Jeon, G. Jang, J.-o. Park, S. Park, Ema system with gradient and uniform saddle coils for 3d locomotion of microrobot, Sensors and Actuators A: Physical 163 (1) (2010) 410-417.

[47] A. Pourkand, J. J. Abbott, A critical analysis of eight-electromagnet manipulation systems: The role of electromagnet configuration on strength, isotropy, and access, IEEE Robotics and Automation Letters 3 (4) (2018) 2957-2962.

[48] D. J. Griffiths, Introduction to Electrodynamics, 4th Edition, Cambridge University Press, 2017.

[49] D. Jiles, Introduction to magnetism and magnetic materials, CRC press, 2015.

[50] E. Diller, J. Giltinan, G. Z. Lum, Z. Ye, M. Sitti, Six-degree-of-freedom magnetic actuation for wireless microrobotics, The Int. J. Robot. Res. 35 (1-3) (2016) 114-128.

[51] E. Diller, J. Giltinan, P. Jena, M. Sitti, Three dimensional independent control of multiple magnetic micro-robots, in: IEEE Int. Conf. on Intel. Robot. and Autom., 2013, pp. 2576-2581.

[52] S. Schuerle, S. Erni, M. Flink, B. E. Kratochvil, B. J. Nelson, Three-dimensional magnetic manipulation of micro-and nanostructures for applications in life sciences, IEEE Trans. Magn. 49 (1) (2013) 321-330.

[53] J. K. Salisbury, J. J. Craig, Articulated hands: Force control and kinematic issues, The International journal of Robotics research 1 (1) (1982) 4-17.

[54] T. Yoshikawa, Manipulability of robotic mechanisms, The international journal of Robotics Research 4 (2) (1985) 3-9.

[55] S. Kucuk, Z. Bingul, Comparative study of performance indices for fundamental robot manipulators, Rob. Aut.

ton. Syst. 54 (7) (2006) 567-573.

[56] J.-P. Merlet, Jacobian, manipulability, condition number, and accuracy of parallel robots, J. of Mech. Des. 128 (1) (2006) 199-206.

[57] S. Schürle, B. E. Kratochvil, S. Pané, M. A. Zeeshan, B. J. Nelson, Generating magnetic fields for controlling nanorobots in medical applications, in: Nanorobotics, Springer, 2013, pp. 275-299.

[58] K. Ishiyama, K. Arai, M. Sendoh, A. Yamazaki, Spiral-type micro-machine for medical applications, in: MHS2000. Proceedings of 2000 International Symposium on Micromechatronics and Human Science (Cat. No. 00TH8530), IEEE, 2000, pp. 65-69.

[59] R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone, J. Bibette, Microscopic artificial swimmers, Nature 437 (7060) (2005) 862-865.

[60] M. N. Faddis, W. Blume, J. Finney, A. Hall, J. Rauch, J. Sell, K. T. Bae, M. Talcott, B. Lindsay, Novel, magnetically guided catheter for endocardial mapping and radiofrequency catheter ablation, Circulation 106 (23) (2002) 2980-2985.

[61] T. Roberts, W. Hassenzahl, S. Hetts, R. Arenson, Remote control of catheter tip deflection: an opportunity for interventional mri, Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 48 (6) (2002) 1091-1095.

[62] A. J. Petruska, B. J. Nelson, Minimum bounds on the number of electromagnets required for remote magnetic manipulation, IEEE Trans. Robot. 31 (3) (2015) 714-722.

[63] C. R. Thornley, L. N. Pham, J. J. Abbott, Reconsidering six-degree-of-freedom magnetic actuation across scales, IEEE Robotics and Automation Letters 4 (3) (2019) 2325-2332.

[64] J. Giltinan, M. Sitti, Simultaneous six-degree-of-freedom control of a single-body magnetic microrobot, IEEE Robotics and Automation Letters 4 (2) (2019) 508-514.

[65] E. D. Diller, J. Giltinan, G. Z. Lum, Z. Ye, M. Sitti, Six-degrees-of-freedom remote actuation of magnetic micro-robots., in: Robotics: Science and Systems, 2014.

[66] Z. Nagy, B. J. Nelson, Lagrangian modeling of the magnetization and the magnetic torque on assembled soft-magnetic mems devices for fast computation and analysis, IEEE transactions on robotics 28 (4) (2012) 787-797.

[67] K. T. Nguyen, M. C. Hoang, G. Go, B. Kang, E. Choi, J.-O. Park, C.-S. Kim, Regularization-based independent control of an external electromagnetic actuator to avoid singularity in the spatial manipulation of a micro-robot, Control Engineering Practice 97 (2020) 104340.

[68] D. Wong, E. B. Steager, V. Kumar, Independent control of identical magnetic robots in a plane, IEEE Robotics and Automation Letters 1 (1) (2016) 554-561.

[69] M. Etiévant, A. Bolopion, S. Régnier, N. Andreff, An improved control-oriented modeling of the magnetic field, in: 2019 International Conference on Robotics and Automation (ICRA), IEEE, 2019, pp. 6178-6184.

[70] X. Jing, W. Guo, Modeling and configuration design of electromagnetic actuation coil for a magnetically controlled microrobot, Chinese Journal of Mechanical Engineering 32 (1) (2019) 63.
# INDEX OF TERMS AND NOTATIONS

## NOTATION RULES

<table><thead><tr><th>Notations</th><th>Description</th></tr></thead><tbody><tr><td><i>x</i> lower case</td><td>a scalar</td></tr><tr><td><b>v</b> bold lower case</td><td>a vector</td></tr><tr><td><b>ê</b> bold lower case</td><td>a unit vector</td></tr><tr><td><b>M</b> bold upper case</td><td>a matrix or a vector field</td></tr><tr><td><b>f</b></td><td>a force vector</td></tr><tr><td><b>t</b></td><td>a torque vector</td></tr><tr><td><b>M</b><sup><i>t</i></sup></td><td>the transpose of matrix <b>M</b></td></tr><tr><td><b>M</b><sup>-1</sup></td><td>inverse of matrix <b>M</b></td></tr><tr><td><b>M</b><sup>+</sup></td><td>pseudo-inverse of matrix <b>M</b></td></tr><tr><td><i>O</i><sub><i>n</i>×<i>m</i></sub></td><td>a zero matrix or null matrix of dimension <i>n</i> × <i>m</i></td></tr><tr><td><i>I</i><sub><i>n</i>×<i>m</i></sub></td><td>identity matrix of dimension <i>n</i> × <i>m</i></td></tr><tr><td><i>i</i></td><td>imaginary unit: <i>i</i><sup>2</sup> = −1</td></tr><tr><td><b>∇</b></td><td>gradient operator</td></tr><tr><td><b>a</b> · <b>b</b></td><td>dot product between vector <b>a</b> and <b>b</b></td></tr><tr><td><b>a</b> || <b>a</b></td><td>norm of <b>a</b></td></tr><tr><td>⟨<b>a</b>(<i>t</i>)⟩</td><td>mean value of <b>a</b>(<i>t</i>)</td></tr></tbody></table>

Unless explicitly stated otherwise, throughout this document the values are stated in SI units.
LIST OF SYMBOLS

<table><thead><tr><th>Symbol</th><th>Description</th><th>Units</th></tr></thead><tbody><tr><td>R</td><td>Microrobot radius</td><td>m</td></tr><tr><td>R<sub>v</sub></td><td>Vessel radius</td><td>m</td></tr><tr><td>V</td><td>Volume</td><td>m</td></tr><tr><td>B</td><td>Magnetic field flux density</td><td>T</td></tr><tr><td>∇B</td><td>Magnetic field flux gradient</td><td>T/m</td></tr><tr><td>E</td><td>Electric field strength</td><td>V/m</td></tr><tr><td>H</td><td>Magnetic field strength</td><td>A/m</td></tr><tr><td>j</td><td>Electric current density</td><td>A/m<sup>2</sup></td></tr><tr><td>i</td><td>Electric current</td><td>A</td></tr><tr><td>m</td><td>Magnetic moment</td><td>A m<sup>2</sup></td></tr><tr><td>M</td><td>Magnetization or magnetic polarization of a magnetic material that correspond to the volumetric density of magnetic moment</td><td>A/m</td></tr><tr><td>μ</td><td>Magnetic permeability</td><td>H/m</td></tr><tr><td>μ<sub>0</sub></td><td>Magnetic permeability of the free space: μ<sub>0</sub> = 4π × 10<sup>-7</sup></td><td>H/m</td></tr><tr><td>A</td><td>Magnetic actuation matrix</td><td>NA</td></tr><tr><td>τ<sub>m</sub></td><td>Magnetic volume ratio (τ<sub>m</sub> = V<sub>m</sub>/V)</td><td>—</td></tr><tr><td>f<sub>m</sub></td><td>Magnetic force</td><td>N</td></tr><tr><td>t<sub>m</sub></td><td>Magnetic torque</td><td>Nm</td></tr><tr><td>f<sub>d</sub></td><td>Hydrodynamic drag force</td><td>N</td></tr><tr><td>t<sub>d</sub></td><td>Hydrodynamic drag torque</td><td>Nm</td></tr><tr><td>κ</td><td>Dynamic shape factor</td><td>—</td></tr><tr><td>v</td><td>Velocity vector</td><td>m/s</td></tr><tr><td>ω</td><td>Rotational velocity vector</td><td>rad/s</td></tr><tr><td>p</td><td>pressure vector</td><td>Pa</td></tr><tr><td>R<sub>e</sub></td><td>Reynolds number that corresponds to the ratio of inertial forces to viscous forces within a fluid</td><td>—</td></tr><tr><td>η</td><td>Viscosity of a fluid</td><td>Pa · s</td></tr><tr><td>ρ</td><td>Density of a material</td><td>kg/m<sup>3</sup></td></tr></tbody></table>
## LISTS OF ACRONYMS

**2D** two-dimension

**3D** three-dimension

**CFD** computational fluid dynamics

**DOF** degrees of freedom

**EMA** electro-magnetic actuation

**ENSI** École Nationale Supérieur d'Ingénieurs

**FMM** fast marching methdod

**FEM** finite element methdod

**GPC** generalized predictive controller

**HDR** Habilitation à Diriger des Recherches

**INSA** Institut National des Sciences Appliquée

**LQI** linear quadratic integral

**LQR** linear quadratic regulator

**LTI** linear time-invariant system

**MIS** minimally invasive surgery

**MNM** micro/nano-motor

**MPC** model predictive controller

**MRN** magnetic resonance navigation

**MRI** magnetic resonance imaging

**MSRL** Multi-Scale Robotics Laboratory

**NANOMA** Nano-Actuators and Nano-Sensors for Medical Applications

**NS** Navier-Stokes equations

**PHC** Partenariats Hubert Curien

**PIANHO** Innovative Haptic Instrumental platform for 3D Nano-manipulation

**PROTEUS** Robotic Platform to facilitate transfer between Industries

**PROSIT** Robotic Platform for an Interactive Tele-echographic System

**RBC** red blood cell

**RMS** root mean square

**SD** standard deviation

**SDC** state dependent coefficients

**SDRE** state-dependent Riccati equation

**SMPC** state-space based model predictive controller

**SPIO** superparamagnetic iron oxide

**TACE** transcatheter arterial chemoembolization

**TARE** transcatheter arterial radioembolisation

**TMMC** therapeutic magnetic microcarriers

**US** (medical) ultrasound or ultrasonography

## LIST OF TERMS

**ANR** French National Research Agency, that is to say in French "Agence Nationale de la Recherche". http://www.agence-nationale-recherche.fr

**biomicromanipulation** Bionicromanipulation commonly refers to the manipulation of small biological material such as micro-organsims, cells, tissues, DNA, etc
**CNU** French National Universities Council, that is to say in French “Conseil National des Universités”, which comprises 57 sections covering different scientific disciplines. The 61st section involve IT engineering, automation and signal processing. https://www.conseil-national-des-universites.fr/cnu/

**ETH** Swiss Federal Institute of Technology in Zurich, from German “Eidgenössische Technische Hochschule”. https://ethz.ch

**hETD** Equivalent TD hours, is the reference hour to calculate the teaching duties. The rules are as follows: 1h of course = 1.5hETD, while the others, e.g. 1h of tutorial (TD) = 1h of practical work (TP) = 1hETD

**INSA CVL** The INSA Centre Val de Loire was established on 1st January 2014 following the merger of the Val de Loire ENI (National Engineering School) and Bourges ENSI (Graduate Engineering School). On the Bourges campus two engineering courses are available to students with the Industrial Risk Control (MRI) and Information Technology and Cybersecurity (STI) departments, and one apprenticeship-training course, Energy Risks and the Environment (ERE), is taught together with the Cher Chamber of Commerce and Industry (CCI) Hubert Curien CFSA (Apprentice Further Training Centre) in Bourges. https://insa-centrevaldeloire.fr

**AMiR** Division of Microrobotics and Control Engineering, that is from German “Abteilung für Mikrorobotik und Regelungstechnik”, is a department of the University of Oldenburg, Germany, headed by Prof. Dr.-Ing. Sergej Fatikov. http://www.amir.uni-oldenburg.de

**medical robotics** A medical robot is a robot used in the medical sciences to help the practitioner in diagnosis and treating patients. They include surgical robots, rehabilitation/assistive robotics, and so on.

**microfluidic** Microfluidics deals with the behaviour, precise control and manipulation of fluids that are geometrically constrained to a small scale (e.g. typically sub-millimeter) at which capillary penetration governs mass transport. It is a multidisciplinary field at the intersection of engineering, physics, chemistry, biochemistry, nanotechnology, and biotechnology, with practical applications in the design of systems in which low volumes of fluids are processed to achieve multiplexing, automation, and high-throughput screening.

**micromotor** Micromotor is a component of microrobot that can convert energy from various sources, such as magnetic fields, light, or chemical bonds, to do mechanical work in the form of directional motion.

**microrobot** Microrobot is commonly seen as a reprogrammable, microscopic machine (e.g. with dimension less than 1mm) with partly or fully self-contained capabilities entitled by on-board motion, perception, and learning.

**microrobotics** Microrobotics is commonly referred to the branch of robotics, which deals with the study and application of miniature ones, such as microscale mobile robots. The term can also be used for system capable of manipulating micrometer-sized components. While the ‘micro’ prefix has been used subjectively to mean “small”, standardizing on length scales should avoid confusion. A nanorobot would have characteristic dimensions at or below $L < 1 \ \mu\text{m}$, or manipulate components on the 1 to 1000 nm size range. Similarly, a millirobot would have characteristic length $L < 1 \ \text{cm}$; and a small-robot would have dimensions $L < 10 \ \text{cm}$. Anything above the milli-scale is, here, globally considered as a macroscale robot, where objects or phenomena are large enough to be visible with the naked eye.
**microswimmer** Microswimmer is a specific locomotion mode of a microrobot that is able to propel and do directional motion in bulk fluid.

**MAE** Department of Mechanical and Automation Engineering (MAE) of The Chinese University of Hong Kong (CUHK). http://www.mae.cuhk.edu.hk

**MSRL** Multi-Scale Robotics Lab is part of Institute of Robotics and Intelligent Systems (IRIS), ETH Zurich, Switzerland. https://msrl.ethz.ch, https://www.iris.ethz.ch

**PRISME Laboratory** Laboratoire Pluridisciplinaire de Recherche en Ingénierie des Systèmes, Mécanique, Énergétique. The PRISME Laboratory is from University of Orléans and INSA Centre Val de Loire (EPRES 4229). The PRISME laboratory seeks to carry out multidisciplinary research in the general domain of engineering sciences over a broad range of subject areas, including combustion in engines, energy engineering, aerodynamics, the mechanics of materials, image and signal processing, automatic control and robotics. The laboratory is split in 2 units: i) Fluids, Energy, Combustion, Propulsion (FECP) and ii) Images, Robotics, Automatic control and Signal (IRAUS). There are about 170 research professors, engineers, technicians and PhD. students working for this laboratory across several sites in Bourges, Orléans, Chartres, and Châteauroux. https://www.univ-orleans.fr/fr/prisme
<figure><img src="image_83.png" /></figure>
BIBLIOGRAPHY

[1] S. Hutchinson, G. D. Hager, and P. I. Corke, "A tutorial on visual servo control," *IEEE Transactions on Robotics and Automation*, vol. 12, no. 5, pp. 651-670, Oct. 1996. DOI: 10.1109/70.538972.

[2] F. Chaumette and S. Hutchinson, "Visual servo control. Part I. basic approaches," *IEEE Robotics and Automation Magazine*, vol. 13, no. 4, pp. 82-90, Nov. 2006. DOI: 10.1109/MRA.2006.250573.

[3] F. Chaumette and S. Hutchinson, "Visual servo control. Part II. advanced approaches," *IEEE Robotics and Automation Magazine*, vol. 14, no. 1, pp. 109-118, Mar. 2007. DOI: 10.1109/MRA.2007.339609.

[4] S. Fatikow and U. Rembold, *Microsystem technology and microrobotics*. Berlin Heidelberg: Springer International Publishing, 1997, ISBN: 978-3-662-03450-7. DOI: 10.1007/978-3-662-03450-7.

[5] B. J. Nelson, L. Dong, and F. Arai, "Micro/Nanorobots," in *Springer Handbook of Robotics*, B. Siciliano and O. Khatib, Eds., Berlin, Heidelberg: Springer, 2008, pp. 411-450, ISBN: 978-3-540-30301-5. DOI: 10.1007/978-3-540-30301-5_19.

[6] B. J. Nelson, I. K. Kaliakatsos, and J. J. Abbott, "Microrobots for minimally invasive medicine," *Annual Review of Biomedical Engineering*, vol. 12, no. 1, pp. 55-85, Jul. 2010. DOI: 10.1146/annurev-bioeng-010510-103409.

[7] M. Sitti, "Miniature devices: Voyage of the microrobots," *Nature*, vol. 458, no. 7242, p. 1121, 2009. DOI: 10.1038/4581121a.

[8] M. Sitti, *Mobile Microrobotics*, ser. Intelligent Robotics and Autonomous Agents. MIT Press, Jun. 2017, ISBN: 9780262036436.

[9] J. J. Abbott, Z. Nagy, F. Beyeler, and B. J. Nelson, "Robotics in the Small, Part I: Microbotics," *IEEE Robotics and Automation Magazine*, vol. 14, no. 2, pp. 92-103, Jun. 2007. DOI: 10.1109/MRA.2007.380641.

[10] M. J. Mack, "Minimally invasive and robotic surgery," *Jama*, vol. 285, no. 5, pp. 568-572, 2001. DOI: 10.1001/jama.285.5.568.

[11] P. Gomes, Ed., *Medical Robotics: Minimally Invasive Surgery*, ser. Woodhead Publishing Series in Biomaterials 51. Cambridge: Elsevier Science, Oct. 2012, ISBN: 9780857097392.

[12] J. Troccaz, Ed., *Medical Robotics*, ser. ISTE. Wiley Online Library, Mar. 2013, ISBN: 9781118563335. DOI: 10.1002/9781118562147.
[13] M. Kroh and S. Chalikonda, Eds., *Essentials of robotic surgery*. Springer International Publishing, Oct. 27, 2014, 218 pp., ISBN: 978-3-319-09564-6.

[14] H. Ryan and S. Tsuda, "History of and current systems in robotic surgery," in *Essentials of Robotic Surgery*, M. Kroh and S. Chalikonda, Eds., Springer International Publishing, 2015, pp. 1–12, ISBN: 978-3-319-09564-6. DOI: 10.1007/978-3-319-09564-6_1.

[15] G. M. Grimsby, M. E. Dwyer, M. A. Jacobs, M. C. Ost, F. X. Schneck, G. M. Cannon, and P. C. Gargollo, "Multi-institutional review of outcomes of robot-assisted laparoscopic extravesical ureteral reimplantation," *The Journal of urology*, vol. 193, no. 5, pp. 1791–1795, 2015. DOI: 10.1016/j.juro.2014.07.128.

[16] R. H. Taylor, A. Menciassi, G. Fichtinger, P. Fiorini, and P. Dario, "Medical robotics and computer-integrated surgery," in *Springer handbook of robotics*, B. Siciliano and O. Khatib, Eds., Cham: Springer International Publishing, 2016, pp. 1657–1684, ISBN: 978-3-319-32552-1. DOI: 10.1007/978-3-319-32552-1_63.

[17] J. Li, B. E.-F. de Ávila, W. Gao, L. Zhang, and J. Wang, "Micro/nanorobots for biomedicine: Delivery, surgery, sensing, and detoxification," *Science Robotics*, vol. 2, no. 4, 2017.

[18] H. Alemzadeh, J. Raman, N. Leveson, Z. Kalbarczyk, and R. K. Iyer, "Adverse events in robotic surgery: A retrospective study of 14 years of fda data," *PloS one*, vol. 11, no. 4, e0151470, 2016. DOI: 10.1371/journal.pone.0151470.

[19] N. Simaan, R. M. Yasin, and L. Wang, "Medical technologies and challenges of robot-assisted minimally invasive intervention and diagnostics," *Annual Review of Control, Robotics, and Autonomous Systems*, vol. 1, pp. 465–490, 2018.

[20] J. Hwang, J.-y. Kim, and H. Choi, "A review of magnetic actuation systems and magnetically actuated guidewire- and catheter-based microrobots for vascular interventions," *Intelligent Service Robotics*, vol. 13, no. 1, pp. 1–14, Jan. 1, 2020. DOI: 10.1007/s11370-020-00311-0.

[21] F. Tendick, S. S. Sastry, R. S. Fearing, and M. Cohn, "Applications of micromechanics in minimally invasive surgery," *IEEE/ASME Transactions on Mechatronics*, vol. 3, no. 1, pp. 34–42, 1998.

[22] M. C. Carrozza, P. Dario, and L. P. S. Jay, "Micromechatronics in surgery," *Transactions of the Institute of Measurement and Control*, vol. 25, no. 4, pp. 309–327, 2003.

[23] S.-J. Baek and S.-H. Kim, "Robotics in general surgery: An evidence-based review," *Asian journal of endoscopic surgery*, vol. 7, no. 2, pp. 117–123, 2014.

[24] F. Courreges, P. Vieyres, and G. Poisson, "Robotized tele-echography," in *Teleradiology*, Berlin, Heidelberg: Springer International Publishing, 2008, pp. 139–154, ISBN: 978-3-540-78871-3. DOI: 10.1007/978-3-540-78871-3_13.

[25] S. Kumar and E. Krupinski, *Teleradiology*. Springer International Publishing, Jul. 25, 2008, ISBN: 9783540788713.

[26] S. Avgousti, E. G. Christoforou, A. S. Panayides, S. Voskarides, C. Novales, L. Nouaille, C. S. Pattichis, and P. Vieyres, "Medical telerobotic systems: Current status and future trends," *BioMedical Engineering OnLine*, vol. 15, no. 1, p. 96, Aug. 12, 2016. DOI: 10.1186/s12938-016-0217-7.

[27] G.-Z. Yang, J. Bellingham, P. E. Dupont, P. Fischer, L. Floridi, R. Full, N. Jacobstein, V. Kumar, M. McNutt, R. Merrifield, et al., "The grand challenges of science robotics," *Science Robotics*, vol. 3, no. 14, eaar7650, 2018. DOI: 10.1126/scirobotics.aar7650.
[28] K. Ishiyama, M. Sendoh, and K. I. Arai, "Magnetic micromachines for medical applications," *Journal of Magnetism and Magnetic Materials*, Proceedings of the Joint European Magnetic Symposia (JEMS'01), vol. 242-245, pp. 41-46, Apr. 1, 2002. DOI: 10.1016/S0304-8853(01)01181-7.

[29] H. Ceylan, J. Giltinan, K. Kozielski, and M. Sitti, "Mobile microrobots for bioengineering applications," *Lab on a Chip*, vol. 17, no. 10, pp. 1705-1724, 2017. DOI: 10.1039/C7LC00064B.

[30] I. Kassim, L. Phee, W. S. Ng, F. Gong, P. Dario, and C. A. Mosse, "Locomotion techniques for robotic colonoscopy," *IEEE Engineering in Medicine and Biology Magazine*, vol. 25, no. 3, pp. 49-56, May 2006. DOI: 10.1109/MEMB.2006.1636351.

[31] G. Ciuti, A. Menciassi, and P. Dario, "Capsule endoscopy: From current achievements to open challenges," *IEEE Reviews in Biomedical Engineering*, vol. 4, pp. 59-72, 2011. DOI: 10.1109/RBME.2011.2171182.

[32] W. Gao, R. Dong, S. Thamphiwatana, J. Li, W. Gao, L. Zhang, and J. Wang, "Artificial micromotors in the mouse's stomach: A step toward in vivo use of synthetic motors," *ACS nano*, vol. 9, no. 1, pp. 117-123, 2015. DOI: 10.1021/nn507097k.

[33] A. Servant, F. Qiu, M. Mazza, K. Kostarelos, and B. J. Nelson, "Controlled in vivo swimming of a swarm of bacteria-like microrobotic flagella," *Advanced Materials*, vol. 27, no. 19, pp. 2981-2988, 2015. DOI: 10.1002/adma.201404444.

[34] O. Felfoul, M. Mohammadi, S. Taherkhani, D. de Lanauze, Y. Zhong Xu, D. Loghin, S. Essa, S. Jancik, D. Houle, M. Lafleur, L. Gaboury, M. Tabrizian, N. Kaou, M. Atkin, T. Vuong, G. Batist, N. Beauchemin, D. Radzioch, and S. Martel, "Magneto-aerotactic bacteria deliver drug-containing nanoliposomes to tumour hypoxic regions," *Nature Nanotechnology*, vol. 11, no. 11, pp. 941-947, 11 Nov. 2016. DOI: 10.1038/nnano.2016.137.

[35] S. Jeon, S. Kim, S. Ha, S. Lee, E. Kim, S. Y. Kim, S. H. Park, J. H. Jeon, S. W. Kim, C. Moon, B. J. Nelson, J.-y. Kim, S.-W. Yu, and H. Choi, "Magnetically actuated microrobots as a platform for stem cell transplantation," *Science Robotics*, vol. 4, no. 30, 2019. DOI: 10.1126/scirobotics.aav4317.

[36] J. Li and D. J. Mooney, "Designing hydrogels for controlled drug delivery," *Nature Reviews Materials*, vol. 1, no. 12, p. 16071, 2016. DOI: 10.1038/natrevmats.2016.71.

[37] E. M. Purcell, "Life at low Reynolds number," *American Journal of Physics*, vol. 45, no. 1, pp. 3-11, 1977. DOI: 10.1119/1.10903.

[38] M. Setti, "Microscale and nanoscale robotics systems: Characteristics, state of the art, and grand challenges," *IEEE Robotics and Automation Magazine*, vol. 14, no. 1, pp. 53-60, Mar. 2007. DOI: 10.1109/MRA.2007.339606.

[39] A. A. Solovev, Y. Mei, E. Bermúdez Ureña, G. Huang, and O. G. Schmidt, "Catalytic microtubular jet engines self-propelled by accumulated gas bubbles," *Small*, vol. 5, no. 14, pp. 1688-1692, 2009. DOI: 10.1002/smll.200900021.

[40] A. Ghosh and P. Fischer, "Controlled propulsion of artificial magnetic nanostructured propellers," *Nano letters*, vol. 9, no. 6, pp. 2243-2245, 2009. DOI: 10.1021/nl900186w.

[41] V. Garcia-Gradilla, J. Orozco, S. Sattayasamitsathit, F. Soto, F. Kuralay, A. Pourazary, A. Katzenberg, W. Gao, Y. Shen, and J. Wang, "Functionalized ultrasound-propelled magnetically guided nanomotors: Toward practical biomedical applications," *ACS nano*, vol. 7, no. 10, pp. 9232-9240, 2013. DOI: 10.1021/nn403851v.
[42] V. Magdanz, M. Medina-Sánchez, Y. Chen, M. Guix, and O. G. Schmidt, "How to improve spermbot performance," *Advanced Functional Materials*, vol. 25, no. 18, pp. 2763-2770, 2015. DOI: 10.1002/adfm.201500015.

[43] P. Dario, R. Valleggi, M. Carrozza, M. Montesi, and M. Cocco, "Microactuators for microrobots: A critical survey," *Journal of Micromechanics and microengineering*, vol. 2, no. 3, p. 141, 1992.

[44] J. J. Abbott, K. E. Peyer, L. Dong, and B. J. Nelson, "How should microrobots swim?" In *Springer Tracts in Advanced Robotics*, vol. 66, Springer Berlin Heidelberg, 2010, pp. 157-167. DOI: 10.1007/978-3-642-14743-2_14.

[45] T. Honda, K. Arai, and K. Ishiyama, "Micro swimming mechanisms propelled by external magnetic fields," *IEEE Transactions on Magnetics*, vol. 32, no. 5, pp. 5085-5087, Sep. 1996. DOI: 10.1109/20.539498.

[46] H. Zhou, C. C. Mayorga-Martinez, S. Pané, L. Zhang, and M. Pumera, "Magnetically Driven Micro and Nanorobots," *Chemical Reviews*, Mar. 31, 2021. DOI: 10.1021/acs.chemrev.0c01234.

[47] K. B. Yesin, K. Vollmers, and B. J. Nelson, "Modeling and Control of Untethered Biomicrorobots in a Fluidic Environment Using Electromagnetic Fields," *The International Journal of Robotics Research*, vol. 25, no. 5-6, pp. 527-536, May 2006. DOI: 10.1177/0278364906065389.

[48] S. Martel, J.-B. Mathieu, O. Felfoul, A. Chanu, E. Aboussouan, S. Tamaz, P. Pouponneau, L. Yahia, G. Beaudoin, G. Soulez, and M. Mankiewicz, "Automatic navigation of an untethered device in the artery of a living animal using a conventional clinical magnetic resonance imaging system," *Applied Physics Letters*, vol. 90, no. 11, pp. 11-13, 2007. DOI: 10.1063/1.2713229.

[49] T. W. Fountain, P. V. Kailat, and J. J. Abbott, "Wireless control of magnetic helical microrobots using a rotating-permanent-magnet manipulator," in *IEEE International Conference on Robotics and Automation (ICRA)*, IEEE, 2010, pp. 576-581.

[50] B. E. Kratochvil, M. P. Kummer, S. Erni, R. Borer, D. R. Frutiger, S. Schürle, and B. J. Nelson, "MiniMag: A hemispherical electromagnetic system for 5-DOF wireless micromanipulation," *Springer Tracts in Advanced Robotics*, vol. 79, pp. 317-329, 2014. DOI: 10.1007/978-3-642-28572-1_22.

[51] M. P. Kummer, J. J. Abbott, B. E. Kratochvil, R. Borer, A. Sengul, and B. J. Nelson, "Octomag: An electromagnetic system for 5-dof wireless micromanipulation," *IEEE Transactions on Robotics*, vol. 26, no. 6, pp. 1006-1017, 2010. DOI: 10.1109/TRO.2010.2073030.

[52] H. Choi, J.-o. Park, and S. Park, "Electromagnetic Actuated Micro and Nano Robots," in *The Encyclopedia of Medical Robotics*, vol. 2, World Scientific, 2018, pp. 1-32. DOI: 10.1142/9789813232280_0006.

[53] W. Amokrane, K. Belharet, M. Souissi, A. B. Grayeli, and A. Ferreira, "Macro-micromanipulation platform for inner ear drug delivery," *Robotics and Autonomous Systems*, vol. 107, pp. 10-19, 2018.

[54] D. Li, F. Niu, J. Li, X. Li, and D. Sun, "Gradient-Enhanced Electromagnetic Actuation System With a New Core Shape Design for Microrobot Manipulation," *IEEE Transactions on Industrial Electronics*, vol. 67, no. 6, pp. 4700-4710, Jun. 2020. DOI: 10.1109/TIE.2019.2928283.
[55] Z. Yang and L. Zhang, "Magnetic Actuation Systems for Miniature Robots: A Review," *Advanced Intelligent Systems*, vol. 2, no. 9, p. 2000082, 2020. DOI: 10.1002/aisy.202000082.

[56] O. Erin, M. Boyvat, M. E. Tiryaki, M. Phelan, and M. Sitti, "Magnetic Resonance Imaging System-Driven Medical Robotics," *Advanced Intelligent Systems*, vol. 2, no. 2, p. 1900110, 2020. DOI: 10.1002/aisy.201900110.

[57] J. Wang, *Nanomachines: fundamentals and applications*. Wiley Online Library, 2013.

[58] S. Pané, J. Puigmartí-Luis, C. Bergeles, X. Z. Chen, E. Pellicer, J. Sort, V. Počepcová, A. Ferreira, and B. J. Nelson, "Imaging Technologies for Biomedical Micro- and Nanoswimmers," *Advanced Materials Technologies*, vol. 1800575, pp. 1-16, 2018. DOI: 10.1002/admt.201800575.

[59] C. R. Reid, H. MacDonald, R. P. Mann, J. A. Marshall, T. Latty, and S. Garnier, "Decision-making without a brain: How an amoeboid organism solves the two-armed bandit," *Journal of The Royal Society Interface*, vol. 13, no. 119, p. 20160030, 2016. DOI: 10.1098/rsif.2016.0030.

[60] R. W. Carlsen and M. Sitti, "Bio-hybrid cell-based actuators for microsystems," *Small*, vol. 10, no. 19, pp. 3831-3851, 2014. DOI: 10.1002/smll.201400384.

[61] J. A. Sethian, *Level Set Methods and Fast Marching Methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science*. Cambridge university press Cambridge, Jun. 1, 1999, 400 pp., ISBN: 0521645573.

[62] J.-B. Mathieu, G. Beaudoin, and S. Martel, "Method of propulsion of a ferromagnetic core in the cardiovascular system through magnetic gradients generated by an MRI system," *IEEE Transactions on Biomedical Engineering*, vol. 53, no. 2, pp. 292-299, 2006. DOI: 10.1109/TBME.2005.862570.

[63] F. Martini, *Fundamentals of Anatomy and Physiology*, 11th ed. New York: Pearson Education, Limited, 2018, 1272 pp., ISBN: 9780134396026.

[64] Y. C. Fung, *Biomechanics: Mechanical Properties of Living Tissues*, 2nd ed. Springer Science & Business Media, 1993, 585 pp., ISBN: 978-1-4757-2257-4.

[65] Y. C. Fung, *Biomechanics: Circulation*, 3rd ed. Springer Science & Business Media, 1997, 587 pp., ISBN: 978-1-4757-2696-1.

[66] W. F. Boron and E. L. Boulpaep, *Medical Physiology*. Elsevier Health Sciences, Dec. 3, 2008, 4220 pp., ISBN: 978-1-4377-2017-4.

[67] E. Gauger and H. Stark, "Numerical study of a microscopic artificial swimmer," *Physical Review E*, vol. 74, no. 2, p. 021907, Aug. 2006. DOI: 10.1103/PhysRevE.74.021907.

[68] S. Tamaz, R. Gourdeau, A. Chanu, J.-B. Mathieu, and S. Martel, "Real-time MRI-based control of a ferromagnetic core for endovascular navigation," *IEEE Transactions on Biomedical Engineering*, vol. 55, no. 7, pp. 1854-1863, 2008. DOI: 10.1109/TBME.2008.919720.

[69] L. Arcèse, M. Fruchard, and A. Ferreira, "Endovascular Magnetically-Guided Robots: Navigation Modeling and Optimization," *IEEE Transactions on Biomedical Engineering*, vol. 54, no. 4, pp. 977-987, 2012. DOI: 10.1109/TBME.2011.2181508>.

[70] R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone, and J. Bibette, "Microscopic artificial swimmers," *Nature*, vol. 437, no. 7060, pp. 862-865, 7060 Oct. 2005. DOI: 10.1038/nature04090.

[71] A. Yamazaki, M. Sendoh, K. Ishiyama, K. Ichi Arai, R. Kato, M. Nakano, and H. Fukunaga, "Wireless micro swimming machine with magnetic thin film," *Journal of*
*Magnetism and Magnetic Materials*, Proceedings of the International Conference on Magnetism (ICM 2003), vol. 272-276, E1741-E1742, May 1, 2004. DOI: 10.1016/j.jmmm.2003.12.337.

[72] B. Behkam and M. Sitti, "Design Methodology for Biomimetic Propulsion of Miniature Swimming Robots," *Journal of Dynamic Systems, Measurement, and Control*, vol. 128, no. 1, pp. 36-43, Mar. 1, 2006. DOI: 10.1115/1.2171439.

[73] K. E. Peyer, L. Zhang, and B. J. Nelson, "Bio-inspired magnetic swimming microrobots for biomedical applications," *Nanoscale*, vol. 5, no. 4, pp. 1259-1272, 2013. DOI: 10.1039/C2NR32554C.

[74] A. F. Tabak and S. Yesilyurt, "Improved kinematic models for two-link helical micro/nanoswimmers," *IEEE Transactions on Robotics*, vol. 30, no. 1, pp. 14-25, 2014. DOI: 10.1109/TRO.2013.2281551.

[75] J. Nam, W. Lee, J. Kim, and G. Jang, "Magnetic Helical Robot for Targeted Drug-Delivery in Tubular Environments," *IEEE/ASME Transactions on Mechatronics*, vol. 22, no. 6, pp. 2461-2468, Dec. 2017. DOI: 10.1109/TMECH.2017.2761786.

[76] S. Palagi and P. Fischer, "Bioinspired microrobots," *Nature Reviews Materials*, vol. 3, no. 6, pp. 113-124, 2018. DOI: 10.1038/s41578-018-0016-9.

[77] L. Zhang, J. J. Abbott, L. Dong, K. E. Peyer, B. E. Kratochvil, H. Zhang, C. Bergeles, and B. J. Nelson, "Characterizing the swimming properties of artificial bacterial flagella.," *Nano letters*, vol. 9, no. 10, pp. 3663-7, 2009. DOI: 10.1021/nl901869j.

[78] S. Martel, M. Mohammadi, O. Felfoul, Zhao Lu, and P. Pouponneau, "Flagellated magnetotactic bacteria as controlled MRI-trackable propulsion and steering systems for medical nanorobots operating in the human microvasculature," *The International Journal of Robotics Research*, vol. 28, no. 4, pp. 571-582, 2009. DOI: 10.1177/0278364908100924.

[79] I. S. Khalil, V. Magdanz, S. Sánchez, O. G. Schmidt, L. Abelmann, and S. Misra, "Magnetic Control of Potential Microrobotic Drug Delivery Systems: Nanoparticles, Magnetotactic Bacteria and Self-Propelled Microjets," in *Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)*, 2013, pp. 5299-5302. DOI: 10.1109/EMBC.2013.6610745.

[80] I. S. M. Khalil, V. Magdanz, S. Sanchez, O. G. Schmidt, and S. Misra, "Precise Localization and Control of Catalytic Janus Micromotors Using Weak Magnetic Fields," *International Journal of Advanced Robotic Systems*, vol. 12, no. 1, p. 2, Jan. 2015. DOI: 10.5772/58873.

[81] A. X. Lu, Y. Liu, H. Oh, A. Gargava, E. Kendall, Z. Nie, D. L. DeVoe, and S. R. Raghavan, "Catalytic Propulsion and Magnetic Steering of Soft, Patchy Microcapsules: Ability to Pick-Up and Drop-Off Microscale Cargo," *ACS Applied Materials & Interfaces*, vol. 8, no. 24, pp. 15676-15683, Jun. 2016. DOI: 10.1021/acsami.6b01245.

[82] D. H. Kim, P. S. S. Kim, A. A. Julius, and Min Jun Kim, "Three-dimensional control of engineered motile cellular microrobots," in *IEEE International Conference on Robotics and Automation*, May 2012, pp. 721-726. DOI: 10.1109/ICRA.2012.6225031.

[83] Y. Alapan, O. Yasa, O. Schauer, J. Giltinan, A. F. Tabak, V. Sourjik, and M. Sitti, "Soft erythrocyte-based bacterial microswimmers for cargo delivery," *Science Robotics*, vol. 3, no. 17, eaar4423, Apr. 25, 2018. DOI: 10.1126/scirobotics.aar4423.

[84] M. Aubry, W.-A. Wang, Y. Guyodo, E. Delacou, J.-M. Guigner, O. Espeli, A. Lebreton, F. Guyot, and Z. Gueroui, "Engineering E. coli for Magnetic Control and the Spatial
Localization of Functions," *ACS Synthetic Biology*, vol. 9, no. 11, pp. 3030-3041, Nov. 20, 2020. doi: 10.1021/acssynbio.0c00286.

[85] A. Cavalcanti, B. Shirinzadeh, T. Fukuda, and S. Ikeda, "Nanorobot for Brain Aneurysm," *The International Journal of Robotics Research*, vol. 28, no. 4, pp. 558-570, Apr. 1, 2009. doi: 10.1177/0278364908097586.

[86] L. P. Chua, J. Zhang, and T. Zhou, "Numerical study of a complete anastomosis model for the coronary artery bypass," *International Communications in Heat and Mass Transfer*, vol. 32, no. 3, pp. 473-482, Feb. 1, 2005. doi: 10.1016/j.icheatmasstransfer.2004.03.016.

[87] A. K. Politis, G. P. Stavropoulos, M. N. Christolis, F. G. Panagopoulos, N. S. Vlachos, and N. C. Markatos, "Numerical modeling of simulated blood flow in idealized composite arterial coronary grafts: Steady state simulations," *Journal of Biomechanics*, vol. 40, no. 5, pp. 1125-1136, Jan. 1, 2007. doi: 10.1016/j.jbiomech.2006.05.008.

[88] C. Vlachopoulos, M. O'Rourke, and W. W. Nichols, *McDonald's Blood Flow in Arteries: Theoretical, Experimental and Clinical Principles*. CRC Press, Jul. 29, 2011, 758 pp., ISBN: 978-1-4441-2878-9.

[89] M. Zamir, *The Physics of Pulsatile Flow*. Springer Science & Business Media, Dec. 6, 2012, 232 pp., ISBN: 978-1-4612-1282-9.

[90] A. Afshari, M. Akbari, D. Toghraie, and M. E. Yazdi, "Experimental investigation of rheological behavior of the hybrid nanofluid of MWCNT-alumina/water (80%)-ethylene-glycol (20%)," *Journal of Thermal Analysis and Calorimetry*, vol. 132, no. 2, pp. 1001-1015, May 1, 2018. doi: 10.1007/s10973-018-7009-1.

[91] U. Dinnar, *Cardiovascular Fluid Dynamics*. CRC Press, Aug. 15, 2019, 241 pp., ISBN: 978-1-00-069688-2.

[92] D. Quemada, "Rheology of concentrated disperse systems II. A model for non-newtonian shear viscosity in steady flows," *Rheologica Acta*, vol. 17, no. 6, pp. 632-642, Nov. 1, 1978. doi: 10.1007/BF01522036.

[93] A. R. Pries and T. W. Secomb, "Blood Flow in Microvascular Networks," in *Comprehensive Physiology*, American Cancer Society, 2011, pp. 3-36, ISBN: 978-0-470-65071-4. doi: 10.1002/cphy.cp020401.

[94] A. R. Pries, D. Neuhaus, and P. Gaehtgens, "Blood viscosity in tube flow: Dependence on diameter and hematocrit," *American Journal of Physiology-Heart and Circulatory Physiology*, vol. 263, no. 6, H1770-H1778, Dec. 1, 1992. doi: 10.1152/ajpheart.1992.263.6.H1770.

[95] A. S. Popel and P. C. Johnson, "Microcirculation and Hemorheology," *Annual Review of Fluid Mechanics*, vol. 37, pp. 43-69, Jan. 1, 2005. doi: 10.1146/annurev.fluid.37.042604.133933.

[96] O. K. Baskurt and H. J. Meiselman, "Blood rheology and hemodynamics," *Seminars in Thrombosis and Hemostasis*, vol. 29, no. 5, pp. 435-450, Oct. 2003. doi: 10.1055/s-2003-44551.

[97] J. Happel and H. Brenner, *Low Reynolds Number Hydrodynamics*, ser. Mechanics of Fluids and Transport Processes. Dordrecht: Springer Netherlands, 1983, vol. 1, ISBN: 978-90-247-2877-0. doi: 10.1007/978-94-009-8352-6.

[98] E. Guyon, J. P. Hulin, L. Petit, and C. D. Mitescu, *Physical Hydrodynamics*. Oxford University Press, 2015, 547 pp., ISBN: 978-0-19-870245-0.
[99] J. R. Womersley, "Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known," *The Journal of Physiology*, vol. 127, no. 3, pp. 553-563, Mar. 28, 1955. DOI: 10.1113/jphysiol.1955.sp005276.

[100] D. W. Holdsworth, C. J. D. Norley, R. Frayne, D. A. Steinman, and B. K. Rutt, "Characterization of common carotid artery blood-flow waveforms in normal human subjects," *Physiological Measurement*, vol. 20, no. 3, pp. 219-240, Jan. 1999. DOI: 10.1088/0967-3334/20/3/301.

[101] T. S. Anirudhan and S. Sandeep, "Synthesis, characterization, cellular uptake and cytotoxicity of a multi-functional magnetic nanocomposite for the targeted delivery and controlled release of doxorubicin to cancer cells," *Journal of Materials Chemistry*, vol. 22, no. 25, pp. 12888-12899, Jun. 6, 2012. DOI: 10.1039/C2JM31794J.

[102] S. Kim, F. Qiu, S. Kim, A. Ghanbari, C. Moon, L. Zhang, B. J. Nelson, and H. Choi, "Fabrication and Characterization of Magnetic Microrobots for Three-Dimensional Cell Culture and Targeted Transportation," *Advanced Materials*, vol. 25, no. 41, pp. 5863-5868, 2013. DOI: 10.1002/adma.201301484.

[103] X. Dong and M. Sitti, "Controlling two-dimensional collective formation and cooperative behavior of magnetic microrobot swarms," *The International Journal of Robotics Research*, vol. 39, no. 5, pp. 617-638, Apr. 1, 2020. DOI: 10.1177/0278364920903107.

[104] J. J. Abbott, E. Diller, and A. J. Petruska, "Magnetic Methods in Robotics," *Annual Review of Control, Robotics, and Autonomous Systems*, vol. 3, no. 1, pp. 57-90, 2020. DOI: 10.1146/annurev-control-081219-082713.

[105] D. J. Griffiths, *Introduction to Electrodynamics*, 4th ed. Cambridge University Press, Jun. 29, 2017, 624 pp., ISBN: 978-1-108-35714-2.

[106] D. Jiles, *Introduction to Magnetism and Magnetic Materials*. CRC Press, Sep. 18, 2015, 616 pp., ISBN: 978-1-4822-3888-4.

[107] A. Sommerfeld, *Electrodynamics*, ser. Lectures on Theoretical Physics. Academic Press, 1964.

[108] E. P. Furlani, *Permanent magnet and electromechanical devices: materials, analysis, and applications*. Academic press, 2001, 540 pp., ISBN: 978-0-12-269951-1. DOI: 10.1016/B978-0-12-269951-1.X5000-1.

[109] L. D. Landau, J. Bell, M. Kearsley, L. Pitaevskii, E. Lifshitz, and J. Sykes, *Electrodynamics of continuous media*, ser. Course of Theoretical Physics. Elsevier Science, 1984, vol. 8.

[110] F. M. White, *Viscous Fluid Flow*, 3rd ed. New York: McGraw-Hill, Jan. 5, 2005.

[111] J. W. Maccoll, "Aerodynamics of a spinning sphere," *The Aeronautical Journal*, vol. 32, no. 213, pp. 777-798, Sep. 1928. DOI: 10.1017/S0368393100136260.

[112] W. L. Haberman and R. M. Sayre, "Motion of rigid and fluid spheres in stationary and moving liquids inside cylindrical tubes," 1143, 1958. [Online]. Available: http://hdl.handle.net/1721.3/48988.

[113] R. Chhabra, S. Agarwal, and K. Chaudhary, "A note on wall effect on the terminal falling velocity of a sphere in quiescent Newtonian media in cylindrical tubes," *Powder Technology*, vol. 129, no. 1-3, pp. 53-58, Jan. 8, 2003. DOI: 10.1016/S0032-5910(02)00164-X.

[114] G. Kasper, T. Niida, and M. Yang, "Measurements of viscous drag on cylinders and chains of spheres with aspect ratios between 2 and 50," *Journal of Aerosol Science*, vol. 16, no. 6, pp. 535-556, Jan. 1985. DOI: 10.1016/0021-8502(85)90006-0.