Magnetic concentric tube robots: Introduction and analysis

In this paper, we propose a new type of continuum robot, referred to as a magnetic concentric tube robot (M-CTR), for performing minimally invasive surgery in narrow and difficult-to-access areas. The robot combines concentric tubes and magnetic actuation to benefit from the ‘follow the leader’ behaviour, the dexterity and stability of existing robots, while targeting millimetre-sized external diameters. These three kinematic properties are assessed through numerical and experimental studies performed on a prototype of a M-CTR. They are performed with general forward and inverse kineto-static models of the robot, continuation and bifurcation analysis, and a specific experimental setup. The prototype presents unique capabilities in terms of deployment and active stability management, while its dexterity in terms of tip orientability is also among the best reported for other robots at its scale.


Continuum robots for minimally invasive surgery
Continuum robots are frequently considered for minimally invasive surgery (MIS) due to their kinematic advantages (Burgner-Kahrs et al., 2015). They usually consist of a slender elastic backbone that is deformed by actuation torques or forces, in order to control their shape and tip pose (see Figure 1). In particular, it is possible to deploy continuum robots with a 'follow the leader' (FTL) behaviour. The robot backbone can be actuated so that its shape corresponds to the path followed by the tip without relying on contacts with its environment (Garriga-Casanovas and Rodriguez y Baena, 2018;Palmer et al., 2014). As a result, continuum robots are well-suited for deployment of medical tools through a patient's anatomy while avoiding sensitive areas, followed by precise manipulation at the operation site. They have successfully been considered in a number of interventions including vascular surgery, cardiac surgery (Gosline et al., 2012), neurology (Burgner et al., 2013), otolaryngology (Simaan et al., 2009), foetal (Dwyer et al., 2017) and abdominal interventions (Burdette et al., 2010).
MIS is constantly evolving towards interventions in smaller and more difficult to access operation sites. Two examples of interventions showing this trend are olfactory cells inspection (Girerd et al., 2018), and middle-ear interventions (Dahroug et al., 2018;Fichera et al., 2017). They involve manipulating medical tools in millimetrescale cavities along several tens of millimetre, inside the nose and the ears, for tissue manipulation and imaging. In particular, such tools must be deployed through tortuous and narrow orifices, while ensuring the integrity of sensitive wall tissues. For example, the olfactory cleft is covered with fragile neural receptors (Kavoi and Jameela, 2011) which would be damaged in case of contact with the tool, resulting in a degraded olfactory sense. In the case of middle-ear surgery, the tympanic membrane and facial nerves can be damaged, resulting in hearing losses and facial distortion (Olszewska et al., 2004). The tools must be manipulated at the operation site while complying with the constraints imposed by the patient's anatomy, which requires dexterity. Being able to rotate tools with a large angular displacement around the operation site is an important feature in tissue manipulation (Simaan et al., 2009) and imaging (Goldman et al., 2013). Limited dexterity may imply several complications such as partial resection and thus regrowth of unhealthy tissues, such as cholesteatoma (Olszewska et al., 2004), and wrong diagnostics due to partial imaging (Mowatt et al., 2011).
The evolution of MIS involves four kinds of requirements for the design of continuum robots. First, there is a need for robots with millimetre diameters, in order to navigate through narrow environments. Second, FTL behaviour is required to access the operation site and to manage the safety throughout the tool's deployment. Third, high dexterity is required in term of tip orientation capability, which has been referred as 'orientability' in the literature (Li et al., 2017;Wu et al., 2017). Fourth, performing FTL deployments and tip orientation control requires the robot to be stable during its motion, in particular in the vicinity of sensitive tissues. Many other requirements can be relevant for the design of continuum robots for minimally invasive surgery. Load bearing capabilities and stiffness properties are of importance in applications involving tissue resection and manipulation for example. Medical devices must also comply to requirements in terms of sterilization, usability and integration to the surgical workflow. In this paper, we focus on the four requirements detailed above, demonstrating an interesting concept of continuum robot for minimally invasive surgery rather than providing a readyto-use and specific medical tool.

Magnetic concentric tube robot
Continuum robots have been classified (Burgner-Kahrs et al., 2015) according to the actuation strategy used to deform the backbone. Each actuation strategy has its own advantages and drawbacks, and it confers to the robot different performances and features as shown in Table 1. Three criteria are used according to the requirements we introduced. The first criterion is the capability of achieving FTL behaviour during deployment. The behaviour can be perfect (P); that is, the robot backbone follows the tip path exactly, or approximated (A), depending on the intrinsic kinematic properties of the robot. The second criterion is the highest angular displacement that can be achieved by the robot tip, which is directly related to orientability (Li et al., 2017). The third criterion is the lowest external diameter achieved among the existing prototypes.
Amongst the intrinsic and extrinsic actuation strategies, only the concentric tube (CT) and tendon driven (TD) actuation have been considered to perform FTL deployments, the last one not being capable of perfect deployments. The angular displacements seem limited for the most part to ±120°, except for robots considering TD and magnetic (M) actuation. External diameters below 1 mm were only achieved with robots using magnetic and electroactive polymer (EAP) actuation. Consequently, there is no design of continuum robot that can achieve simultaneously perfect FTL deployments, over 120°angular displacement and have a millimetre diameter with a single actuation strategy. Considering hybrid strategies seems promising to obtain new robot design with improved performances. In particular, the hybrid design that combines tendon-driven actuation and concentric tubes seems to provide interesting compromises. However, the prototypes presented in Amanov et al. (2017) and Swaney et al. (2016) cannot satisfy the three requirements simultaneously. In particular, performing perfect FTL behaviour requires to assemble two robots concentrically while integrating tendons, making its miniaturization challenging. In this work, we introduce and investigate the use of concentric tubes with magnetic actuation in order to benefit from the combined advantages of both concentric tube robots (CTR) and magnetic continuum robots (M-CR).
A CTR consists of a telescopic assembly of pre-curved elastic tubes (Webster et al., 2006;Dupont et al., 2010) that are usually made out of Nitinol (NiTi). When assembled, the tube interactions create forces and torques along the backbone, which deform the robot. The tubes are rotated and translated at their base using actuators to control the robot shape and tip pose. CTRs are interesting because they can be deployed in a perfect FTL manner for specific shapes and configurations of the tubes  Garriga-Casanovas and Rodriguez y Baena, 2018) and because they can be fabricated with millimetre external diameters due to the simple assembly of the tubes. Thanks to these two features, they have been proven to be useful for deployment through narrow and sensitive areas such as the olfactory cleft (Girerd et al., 2018) and the brain (Comber et al., 2016). However, CTRs also exhibit complex kinematic behaviour that can make them cumbersome to control. Most notably, they suffer from elastic instabilities due to the torsion of the tubes (Dupont et al., 2010;Gilbert et al., 2016;Ha et al., 2016). For given rotation and translation of the tubes, the robot can have multiple configurations and may snap from one configuration to another. These unstable phenomena are usually avoided in order to ensure patient safety. It is realized through careful tube design, which implies limiting the tube pre-curvature (Kim et al., 2014;Bergeles et al., 2015;Ha et al., 2017) and thus the tip angular displacement. These phenomena have recently been considered as useful for tasks requiring a minimal amount of force such as suturing (Riojas et al., 2018). The high dynamics induced by the snapping behaviour are then used to generate tip forces that could not be produced in a quasistatic fashion due to the robot's flexibility.
M-CRs, also called magnetic catheters in the literature, consist in a flexible backbone along which magnetic elements are fixed Tunay, 2004). An external magnetic field is generated using an electromagnetic navigation system (eMNS), in order to generate torques and forces at specific locations along the backbone. Controlling the magnitude, the orientation, and the gradient of the magnetic field alters the robot's shape and the tip's pose. Due to low bending stiffness and the application of torques directly on the backbone, M-CRs can achieve high angular displacements during magnetic actuation (Chautems et al., 2018). M-CRs can also be fabricated at sub-millimetre scales when passive elements are used (Charreyron et al., 2019). Achieving FTL behaviour with M-CRs has never been considered to the best of our knowledge.
The continuum robot with hybrid actuation we introduce and analyze in this paper is designated as a magnetic concentric tube robot (M-CTR). It is designed to combine the properties of CTRs and M-CRs. As shown on Figure 1, M-CTRs consist of a telescopic assembly of tubes on which one or several passive magnetic elements are fixed.
The tubes can be pre-curved or initially straight and are composed of materials with different stiffness properties, such as NiTi or silicone. M-CTRs are actuated by rotating and translating the tubes and by varying the magnetic field applied on the robot. This concept aims at obtaining robots with millimetre diameter and with FTL deployment capabilities and orientability that could not be achieved using CTR and M-CR technologies alone. In addition, the M-CTR may possess unique abilities in its interaction with the environment. Indeed, magnetic actuation has an effect on concentric tube elastic stability. Applying external forces and torques on a CTR can stabilize or destabilize the robot as shown in Ha et al. (2016). As a result, active stability management may be possible, to ensure patient safety while generating higher tip forces.

Contributions
In this paper, the concept of M-CTRs was investigated throughout three studies, each leading to different contributions. First of all, FTL deployment capability of the hybrid continuum robot was assessed. As it requires the use of magnetic actuation to achieve FTL capability for the first time, a deployment strategy is proposed and validated numerically. Second, the robot orientability was evaluated. We especially consider a scenario inspired from tissue manipulation and imaging, where the robot tip must be rotated at the operation site about axes fixed with respect to the target. This lead us to evaluate orientability numerically while considering physical phenomena such as tube torsion and deformation due to gravity and magnetic actuation for The denomination of the classes is inspired from Burgner-Kahrs et al. (2015). A and P indicates approximate and perfect FTL deployment, respectively. The values correspond to the best results reported in the literature and may come as a result from different papers. Indices 1, 2 and 3 are used to indicate which paper the values are from.
the first time. Open-loop control was also proposed for rotating the tip around the target which was validated experimentally. Third, the M-CTR stability was evaluated numerically and experimentally, as well as active control of the robot's stability. A strategy for identifying adequate magnetic fields is provided and validated. Length of the subsection k of section j α i Rotation angle of tube i B Magnetic field G Vector built with the independent components of the magnetic field gradient R 0 ,R b ,R i World, backbone and tube i frames Cartesian position of the backbone Q Quaternion vector representing the backbone orientation u Darboux vector representing the backbone curvature b u i Darboux vector representing the pre-curvature of tube i θ i Torsion angle of tube i k bi , k ti Bending and torsional stiffness of tube i τ m ,f m Magnetic torques and forces m Dipole moment f g Gravity forces ρ i Linear mass of tube i g Gravity acceleration λ p ,λ q Lagrange multipliers x State space vector as a function of s q Actuation inputs vector N Number of nodes used to discretize the robot backbone k index of node x k State space vector evaluated at node k p N ,Q N Position and orientation of the robot tip e 3 Unitary vector e 3 ¼ ½ 0 0 1 T f Relative angle between the magnet direction and the magnetic field in the (y 0 , z 0 ) plane f M Maximum angle f for which the robot stays stable Operators a x , a y , a z Components of vector a, a ¼ ½ a x a y a z T a=fa 1 ,a 2 g Vector a without the components a 1 and a 2 jaj Euclidean norm of vector a a ×b Cross product of a and b ab Matrix product of a and b =(a) Gradient of a The paper is organized as follows. The methods and material used to realize the three studies are described in section 2. The results concerning the FTL deployment capabilities, orientability and stability of the M-CTR are described in sections 3, 4 and 5, respectively. According to these results, the interest of the proposed concept of continuum robot is discussed in section 6.

Notations
The abbreviations and notations used in the rest of the paper are gathered in Table 2.

Methods and materials
FTL deployment capabilities, orientability and stability are assessed using a combination of numerical and experimental approaches. The numerical evaluations are based on a general kineto-static model of M-CTR.

Kineto-static model
The kineto-static model relates the M-CTR configuration to its actuation inputs and external forces. It is elaborated for arbitrary M-CTR designs following the steps introduced in Peyron et al. (2018).
2.1.1. Description and assumptions. M-CTR kineto-statics rely on the interaction of tubes with different stiffness and magnetic properties. We introduce a specific parametrization of the robot which can take into account these different properties. The M-CTR is composed of n tubes which are indexed from the outermost to the innermost. When assembled, the tubes form the backbone of the robot, which is then composed of n sections with a different number of tubes in interaction. In addition, due to the magnetic elements along the tubes, the sections can be decomposed into several subsections with different magnetic and mechanical properties. This robot segmentation is represented for the design under investigation on Figure 2. The 3 tubes form 3 sections, and section 3 is composed of two subsections due to the presence of the magnet at the tip of the inner tube. The sections are indexed from the proximal to the distal end. The subsections are indexed from the proximal to the distal end of the corresponding section, the subsection l of section j being labelled jl. The length of subsection jl is denoted ΔL jl . The first subsection length ΔL j1 is determined by the relative translation of the tubes, while the others are fixed and depend on the length of the magnetic elements and their location along the tubes.
The actuation inputs are the proximal rotation and translation of each tube, and the magnetic field generated on the magnetic elements. The proximal rotation of tube i is denoted α i . The tube proximal translations are parametrized by ΔL j1 , j = 1…n. We consider that the magnetic field can be controlled at one location in the workspace, which is generally the case for eMNSs, and that it is not homogeneous. As a result, it varies across the workspace, including along the robot backbone. The magnetic field generated on the magnetic elements situated at an arc-length s from the robot base is denoted B(s). These actuation inputs are related to the robot configuration under the following assumptions: · The compression and shear deformations are negligible, which is true for thin and slender continuum robots Kratchman et al., 2017). · The materials composing the tubes have a linear and isotropic behaviour . · The tubes are guided in the transmission section so that they remain straight, which is generally the case for CTRs (Dupont et al., 2010). · The robot evolves in free-space. Gravity is the only external force applied on the robot backbone.
2.1.2. M-CTR configuration. The M-CTR configuration comprises the robot backbone position, orientation and curvature, and the torsion angle of each tube. These geometrical entities vary along the backbone, which is parametrized by the arc-length s. Dependency in s is not mentioned in the following for sake of compactness. The backbone position is denoted p and is expressed in the eMNS reference frame denoted R 0 ðO,x 0 ,y 0 ,z 0 Þ. The robot orientation is parameterized with a Bishop frame R b ðp,x b ,y b ,z b Þ where z b is tangent to the backbone. The rotation between R 0 and R B is represented by the unit quaternion vector Q which satisfies Compared to the Euler angles used in Peyron et al. (2018), this representation is singularity-free, induces quadratic equations more efficiently solved with numerical tools and is classically considered for M-CRs and elastic rods (Tunay, 2013;Edelmann et al., 2017). The backbone curvature is represented by a Darboux vector expressed in R b and denoted as b u. Given the construction of R b , the coordinate along z b of b u always equals 0, so that b u ¼ ½ u x u y 0 . The torsion angle θ i of tube i is defined by attaching a frame R i ðp,x i ,y i ,z i Þ at each cross section of the tube so that z i is tangent to the tube centre-line. The pre- curvature of tube i expressed in this frame is denoted by i b u i . Frame R i is then obtained by applying a rotation (θ i , z b ) on R b . In the following, we use left superscript 0, b and i to indicate in which frames vectors are expressed.
Position, orientation and curvature vectors constitute a redundant representation of the robot configuration. As demonstrated in Nikravesh et al. (1985), orientation and curvature are linked by the following relation where prime denotes the derivative with respect to s and L q is a matrix depending on Q. Its expression is given in Appendix B.
The backbone position and orientation are also linked. Since shear and extension are neglected, the axis z b remains tangent to the backbone. The relation between p and Q is then where 0 R b is the rotation matrix from R b to R 0 and e 3 ¼ ½ 0 0 1 T .
2.1.3. Mechanical equilibrium. The equilibrium conditions relate the robot stiffness to the forces induced by the magnetic field and gravity. We first express them along subsection jl, where n j tubes interact. The robot stiffness depends on the stiffness of each tube present in the subsection, which is expressed as a stiffness matrix K i for tube i where k bi and k ti are, respectively, the bending and torsional stiffness of tube i. The magnetic field and the gravity exert magnetic forces f m , torques τ m and gravity forces f g which are distributed along the robot's sections. With m as the dipole moment of the magnetic element present on subsection jl, ρ i the linear mass of tube i and g the gravity acceleration, the forces and torques are written where n j is the number of tubes interacting along subsection j, Fð 0 R b b mÞ is a 3 × 5 matrix formed with the components T contains the 5 independent components of the magnetic field gradient as described in Petruska and Nelson (2015). The equations describing mechanical equilibrium are obtained by first computing the total potential energy due to the robot deformation, and the forces and torques applied on the backbone. The Euler-Lagrange formula is then applied to this energy formulation considering the geometrical constraints (1), (2) and (3). The equations consist in a set of differential equations which are written as The first equation describes the equilibrium in torsion of tube i. Since all tubes can rotate freely with respect to each other, it is written for each tube present on the subsection. The torsional magnetic torque τ mz applies on tube i if a magnet is fixed on this tube, and equals ð b R 0 0 τ m Þ T e 3 . For the other tubes, τ mz = 0. The other equations describe the equilibrium in bending of the robot backbone. The two matrices λ q and λ p contain the Lagrange multipliers due to the geometrical constraints. Multipliers λ q are related to constraints C q1 , C q2 and λ p to constraint C p . They correspond, respectively, to the internal moments and forces due to the gravity and the magnetic field. The intermediate steps of the model derivation are provided in Appendix B.
The resulting equilibrium equations are valid on the subsection jl. Following the approach in Ha et al. (2016) and Peyron et al. (2019a), we extend them to the whole length of the robot by introducing virtual tubes when needed. This means that each tube is considered to cover the entire length of the robot, but with different mechanical and magnetic properties according to the presence of the tubes and magnetic elements. When a tube i is not present along section j, its bending stiffness and linear density are considered as zero and its torsional stiffness as constant. This implies that the stiffness, linear density and dipole moment are piece-wise constant functions of the arc-length s. Equations (7) are then valid for every section of the M-CTR.
2.1.4. Boundary conditions. The differential equation (7) and the geometrical constraints (2) and (3) are subject to proximal (at s = 0) and distal (at s = L) boundary conditions. First, the backbone position and orientation at s = 0, denoted by p 0 and Q 0 , respectively, appear in the following proximal boundary condition The rotation of tube i at its proximal extremity also induces a proximal condition on the tube torsion angle θ i . In particular, the tube accumulates torsional deformations in the transmission section, resulting in a torsion angle at s = 0 that is lower than the actuation angle α i . Since the tubes are considered as straight along the transmission section, their torsional curvature is constant along this section (Dupont et al., 2010). The torsion angle of tube i at s = 0 is then written where the transmission length β i depends on the tube and section length as follows Since no force or torque is applied at the robot tip, the tubes are not deformed in torsion at their distal extremity, leading to the distal boundary condition This assumption also implies a distal condition for the internal moments and forces λ q and λ p The equations describing the robot geometry (1-3) and mechanical equilibrium (7), and the boundary conditions (8-12), constitute the kineto-static model of M-CTR. This model allows for computing the 16 + n robot states, denoted by

Numerical framework
CTRs and M-CRs are known to present multiple configurations for a single set of actuation inputs. Loss of stability can also be encountered during robot deployments. Therefore, numerical analysis is performed with the framework for cardinality and stability analysis of continuum robots developed in Peyron et al. (2019a). It has been used successfully to evaluate the FTL deployment, orientability and stability of CTRs (Peyron et al., 2019a(Peyron et al., , 2019b and to study the magneto-elastic phenomena of M-CRs (Peyron et al., 2018).

Discretization.
The kineto-static model of a M-CTR is a set of non-linear differential equations with proximal and distal boundary conditions. In order to solve it, the model is discretized using finite differences. First, the robot is discretized as a finite number N of nodes along its backbone. In particular, each section is discretized with 20 nodes, which proved to result in a reasonable accuracy in previous work. Each node is indexed from the proximal to the distal extremities of the M-CTR. The states composing the robot configuration are then evaluated at each node. As a result, the robot configuration is represented with a finite number of states that we gather in a matrix X ¼ ½ x 1 … x N , where x 1 is the state vector evaluated at node 1. Finally, the differential equations are evaluated at each node by replacing the derivatives by finite differences. A central second order finite difference scheme is used for the second order derivative of θ i , as proposed in Peyron et al. (2019a). Backward and forward first order finite difference schemes are used for (p, Q) and (λ q , λ p ), respectively, similarly to Peyron et al. (2018), in order to take into account respective proximal and distal boundary conditions. After discretization, the kineto-static model consists in a set of non-linear equations of the general form is the pose of the robot distal tip. We then obtain an implicit kineto-static model, which is used to solve forward and inverse kinematics with the same set of equations. Forward kinematics are obtained by considering q as the input and the robot configuration X as the output. Inverse kinematics are obtained by considering the tip pose (Q N , p N ) as the input and Y and q as outputs.
2.2.2. Kinematic and stability analysis. FTL deployment, orientability and stability analysis were performed on Matlab R2015a (Mathwork Inc.) using the Matcont Toolbox, which implements a continuation method, a step size control algorithm and bifurcation analysis. FTL deployment capabilities were assessed by simulating deployments using the continuation method. During the simulation, the section lengths were increased sequentially and the forward kinematics were solved at each step with a prediction and correction process. As a result, robot configurations are computed while being robust to non-linear behaviours as demonstrated in the case of CTR in Peyron et al. (2019a).
Orientability was evaluated by using the same continuation process as for the inverse kinematics. The tip position p N was then fixed, and the tip orientation Q N was varied while the corresponding robot configurations were computed. During the simulation, the orientation variation applied at each step was automatically controlled using a dedicated algorithm, in order to ensure convergence of the model solution. In particular, the variation decreased automatically and ultimately vanished when approaching orientation limits, beyond which the inverse kinematics do not have solutions. As a result, achievable tip orientations were obtained.
Stability changes due to magnetic actuation were evaluated by detecting the appearance of multiple robot configurations for the same actuation inputs. As with CTRs, unstable phenomena observed in M-CRs are characterized by cardinality changes as shown in Peyron et al. (2018). These changes appear as bifurcations during the continuation process. Following Peyron et al. (2018;2019a), actuation inputs were varied, the forward kinematics were solved and the appearance or vanishing of bifurcations was assessed with bifurcation analysis.

Design of M-CTRs
In order to simplify the analysis, a M-CTR with 3 tubes and controlled with a homogeneous magnetic field was considered. The general structure of the considered robot is depicted in Figure 2. It is composed of two pre-curved tubes and one straight inner tube with a magnetic element fixed at the tip. This results in a concentric assembly of a 2-tube CTR and a M-CR. The two pre-curved tubes were actuated in rotation and translation, and the inner tube in translation only. Since the magnetic field is considered as homogeneous, its gradient is null and does not generate forces on the magnetic element. It can generate torques to bend the robot with 2 degrees of freedom (d.o.f.). As a consequence, the M-CTR has a total of 7 d.o.f. The robot motion due to these d.o.f. is presented in the video of Extension 1 (see Appendix A, Table 6).
Two objectives were considered for the determination of the tube geometrical, mechanical and magnetic properties. The first objective is to allow the investigation of FTL deployment capabilities, orientability and stability. The second objective is to keep consistency with the targeted application, which requires a robot with millimetre diameter and capable to perform in-vivo imaging. The initial shape of the pre-curved tubes was chosen in order to deploy the robot with a perfect FTL behaviour, following the conditions given in Gilbert et al. (2015). The tubes pre-curvature was then optimized in order to observe stability changes under actuation, so that the influence of the magnetic field on the robot stability could be observed. The optimization was conducted while limiting the tube deformation to ε = 0.5% for which the Nitinol is linear (Iasnii and Junga, 2018). The expressions of tube deformation and elastic stability limit in Dupont et al. (2010) were used as objective functions. A Nelder-Mead numerical method, implemented in the fsolve function in Matlab, was used for the optimization. Finally, the inner tube stiffness was chosen so that large deformations of the distal subsection 31 were obtained under magnetic actuation, maximizing tip angular displacement of the M-CTR and therefore its orientability.
The resulting robot is shown in Figure 5. The two precurved tubes consist of NiTi tubes with planar and constant pre-curvatures. They were obtained from straight tubes 3.89 × 10 À4 2.00 × 10 À3 1.20 × 10 À6 4.88 × 10 À1 k ti (N.m 2 ) 2.88 × 10 À4 1.50 × 10 À4 8.94 × 10 À7 3.84 × 10 À1 ρ i (kg.m À3 ) 6.45 × 10 3 6.45 × 10 3 1.14 × 10 3 7.46 × 10 3 m z (N.m À2 .T À1 ) 0 0 0 1.17 × 10 4 Maximum section lengths Section jl 11 21 31 32 ΔL jl (m) 50.00 × 10 À3 50.00 × 10 À3 13.00 × 10 À3 3.00 × 10 À3 Initial position and orientation (Euroflex Inc.) which were pre-shaped with a heating process. The tube pre-curvatures after fabrication were measured with a FARO arm and the values are reported in Table 3. The external diameter of the outer tube is 1.07 mm. An optical fibre of 250 μm diameter (SMF28, Corning Inc.) and a permanent ring magnet were assembled on the M-CTR. The use of optical fibres is consistent with in-vivo imaging scenarios (Girerd et al., 2018), and the desired level of flexibility and linearity of material behaviour (Antunes et al., 2008). The magnet is chosen with a conventional axial magnetization, which allows to demonstrate the targeted properties while making the results simple to interpret. A 4 mm diameter magnet is used, which allowed to generate enough magnetic torques considering the specifications of the eMNS. All design parameters are gathered in Table 3.

Experimental setup
The experimental analysis was conducted using the setup presented in Figure 3 It is composed of a M-CTR prototype and its actuation unit (A2), the CardioMag eMNS  for magnetic field generation (A1) and a vision-based measurement system composed of a stereo-vision camera system (A3) and a front camera (A4). The M-CTR actuation unit for tube rotation and translation was specifically designed to be compatible with the Car-dioMag. The different components are interfaced with a master computer using the Robot Operating System (ROS). They are presented in more detail in the following.
2.4.1. Actuation unit. The actuation unit is inspired from Hendrick et al. (2015b) but is designed to be compatible with the eMNS. The actuation unit is presented in Figure 4. The tubes are mounted on two translation stages (B2) actuated in rotation and translation using electro-mechanical motors (B1). The entirety of the structure is composed of a non-magnetic material. The structural parts were 3D printed using an ABS polymer. Glass and polymer bearings were used to guide the different axes. The screws used for translation were made out of non-magnetic steel (Inox 1.430 1), and the axes for rotation out of aluminium. The tubes were fixed to the translation stages using a 3-finger clamp (B3) machined from a block of ABS. The CTR actuation consists of four stepper motors (17HM15-0904S, Gotronic) located at the unit base, which are not sensitive to the magnetic fields used in this work. The motors are controlled using an Arduino Mega 2560 board. The actuation unit is equipped with sensors in order to calibrate the tube translation and rotation. Mechanical switches (B4) are used to calibrate the tube translations. In order to calibrate the tube rotations, a mechanical coding system comprising two pins was developed. A first planar pin (B7) is glued to the tubes on a surface plate, so that the pin plane corresponds to the tube plane. A second cylindrical pin (B8) allows for fixing the orientation of the clamp, and thus of the tube, with respect to the gear (B6). An optical sensor (B5) is used to detect a coloured marker placed on one face of the gear (not represented here), which initializes the tube rotation.
2.4.2. Electromagnetic navigation system. The CardioMag is composed of eight fixed coils supplied with independently controlled currents. It is able to generate any magnetic field or field gradient at one point of its workspace, with a maximum magnetic field amplitude of 80 mT . The magnetic field varies from that point across the workspace and was modelled in . This model was used to compute the current reference for each coil, according to the desired magnetic field and application point. It was also used to predict the magnetic field along the robot B(s), and its gradient for given currents in the electromagnets.
2.4.3. Vision-based measurement system. A vision-based strategy is adopted to measure the M-CTR pose without disturbing the magnetic field. A stereo tracking system was  implemented using two Basler A602fc cameras (656 × 490 pixels, 15 Hz) fixed on the ceiling (A3). Coloured hollow beads of 6 mm diameter were fixed to the robot (C1). The centres of the beads are detected using standard image processing functions, providing a position accuracy of 1 mm. The positions of the beads are expressed in frame R 0 , represented on Figure 5. The transformation between the cameras and R 0 is measured using an April tag (C2), and the AprilTag library (Olson, 2011).
The tip orientation is measured when needed using an April tag fixed on the magnet (C3). Orientation estimation of the marker using the stereo system is not sufficient. A third Basler A602fc camera, focused on the robot tip (A4), provided the missing orientation measurement. The camera was positioned manually in order to maximize the April tag detection during its motion. The camera pose with respect to the robot frame is not known. It is thus used to measure tip position and angular displacements instead of absolute tip position and orientation. The measurement errors were determined to be less than 1 mm and 1°, respectively. They correspond to less than 2% of the robot length and 5% of the tip angular displacement, respectively, which is reasonable for assessing experimentally the orientability and stability of M-CTR.

FTL deployment
In this section, a strategy for deploying M-CTR in a FTL manner using magnetic actuation is developed, which is analyzed and validated numerically.

Perfect FTL deployment of a M-CTR
When a path is composed of sections with constant curvatures, FTL deployment with a continuum robot is possible. The robot must be composed of sections with constant curvature and variable length. CTRs can fulfil these two requirements by considering tubes with constant precurvature, planar or helical, and by placing them so that their curvatures are aligned or opposed . However, since the tube curvature is fixed, the robot can follow a limited number of paths following the tube alignments.
Magnetic actuation can be used to overcome this limitation with the considered design of M-CTR. Aligning or opposing the pre-curved tubes implies that the robot is contained in a plane. Without loss of generality, we here choose to work in the (y 0 ,z 0 ) plane as shown on Figure 6(a). Applying a magnetic field in this plane generates a constant bending moment along the backbone (Chautems et al., 2018). As a result, neglecting the effect of gravity, the curvature of the robot backbone stays constant along the sections, but varies according to the magnetic field. Considering that the robot is in the (y 0 , z 0 ) plane and the permanent magnet has an axial magnetization, the expression of the backbone curvature in equation (7) Controlling the magnetic field so that the magnetic torque τ mx stays constant during the deployment leads to a constant value of b u, and the condition for perfect FTL behaviour is fulfilled. The torque depends on the magnetic field magnitude and orientation with respect to the magnet. We chose to keep the magnetic field magnitude constant during the deployment. Therefore, the magnetic field direction is computed so that its relative orientation with the magnet f (see Figure 6(a)) stays constant with the relation where Q x (f) is the quaternion representation of the rotation about x 0 of angle f, and R(Q) is the rotation matrix associated to the quaternion Q. The two quaternion vectors are multiplied using the Hamilton product operator. The effectiveness of this strategy was analyzed by simulating the deployment for several values of f. The pre-curved tubes were placed so that their curvatures were aligned, that is, α 1 = α 2 = 0. The magnetic field was assumed to be controlled at the magnet centre and its magnitude was fixed at jBj ¼ 80 mT. The gravity was set to 0. The section lengths were increased sequentially from the proximal to the distal section until their maximum value was reached (see Table 3). The deployment consists of 3 steps during which Δ L 11 , Δ L 21 and Δ L 31 are varied, respectively. The section length variations and the resulting robot shape are computed using the continuation method. The error of the FTL deployment is evaluated as the mean distance between the desired backbone node position along the path and the actual node position for each value of section length.
The results of the simulations are presented on Figure 6, where robot configurations during the deployment and the followed paths are represented. A video of the complete deployment is provided in Extension 2 (see Appendix A, Table 6). The M-CTR exhibits a perfect FTL behaviour as expected. It is able to follow the desired paths with a mean FTL error during the deployments below 0.54 mm, which represents 0.47% of the robot total length. Such an error is negligible with respect to the values reported in the literature (Garriga-Casanovas and Rodriguez y Baena, 2018). Moreover, the M-CTR distal subsection 31 can be deployed along a curved path although the inner tube is initially straight, and different curvatures of this path can be obtained for different values of f. This is an improvement over the deployment capabilities of CTRs, for which this subsection would remain straight.

Influence of the magnetic field orientation
The range of paths that can be followed with perfect FTL behaviour depends on the achievable range of f. Following equation (14), the magnetic bending moment is maximum when the magnetic field is orthogonal to the magnet axis, that is, when jfj ¼ π=2. In practice however, the maximum value of f is limited by the magneto-elastic instabilities inherent to the behaviour of M-CR (Peyron et al., 2018;Tunay, 2004). When the magnetic field is oriented in the robot plane, there is a critical magnetic field orientation above which the robot becomes unstable and can have multiple configurations for the same inputs. This critical orientation then determines the maximum value of f, denoted f M that can be achieved while conserving stability of the M-CTR.
The value of f M is computed by considering the worst case scenario in terms of magneto-elastic stability. Tunay (2004) proved that the longer the robot, the higher the risks of becoming unstable. The M-CTR is then considered fully deployed for the evaluation of f M . We use then the same process as described in Peyron et al. (2018). The angle f is varied with the continuation process until a stability change is detected. The value of f M corresponds to the value of f at the stability change, which equals f M = 6.31°with the parameters of the M-CTR prototype. The corresponding range of paths is computed with equations (2), (3) and (14) and is presented in Figure 7. The paths exhibit large variations in curvature for the distal section and small variations for the proximal ones. This is due to the low bending stiffness of the inner tube with respect to the NiTi tubes and the magnetic field magnitude. For low values of f, the bending moment induced by the magnetic field is too small to affect the shape of the NiTi tubes, but is high enough to deform the inner tube significantly.
As a consequence, the curvature of the robot distal section can be varied continuously by changing the relative orientation between the magnetic field and the magnet. Moreover, significant variations of curvature were observed. The minimum radius of curvature was 5.4 mm, which is interesting considering the scale and the geometry of environments in which such type of robot could be deployed. This implies that tip angular displacements of ±160°can be  achieved, which is comparable to the best displacements presented in the literature according to Table 1.

Influence of the inner tube stiffness
The inner tube stiffness is an important quantity for the deployment capabilities of the M-CTR. A low inner tube stiffness with respect to the NiTi tube sets the global shape of FTL paths to have constant curvatures along subsections 11 and 21, and a variable curvature for subsection 31. The stiffness also has a significant impact on the magneto-elastic phenomena and then on the value of f M . The evolution of f M according to k b3 is presented in Figure 8. The values of stiffness were chosen in the range [1 × 10 À8 , 1 × 10 À4 ] N.m 2 , with a logarithmic distribution, to consider materials such as Nylon, NiTi, and stainless steel. For the three highest values of k b3 , the robot does not experience any stability changes under magnetic field rotations. Therefore, we report on the graph the value of f which produces the maximum magnetic torque, which is equal to 90°, as explained previously. As the stiffness decreases, f M decreases rapidly and stays below 3°for k b3 < 5 × 10 À7 N.m 2 .
This variation of f M implies a variation of the bending moment applied on the robot, which impacts the deformation of the pre-curved tubes due to magnetic actuation.
The paths which can be followed in a FTL manner are computed for each value of stiffness and for the extreme values of f denoted f M and À f M . A subset of these paths is presented in Figure 9. Large values of k b3 imply large values of f M and consequently large values of the magnetic bending moment. As a result, the NiTi tubes are deformed significantly by the magnetic field. Such values also imply that variations of subsection 11 curvature are smaller. Changing the value of f globally modifies the shape of the path, while the curvature experiences small variations from one section to another. As the stiffness decreases, so does f M and the deformation of the NiTi tubes. For k b3 ≤ 5.62 × 10 À7 N.m 2 , their deformation is negligible and the curvature of section 2 is high, leading to paths with similar shapes than the paths obtained with the M-CTR prototype. Note that in practice this curvature is also limited by self-collisions between the robot tip and its backbone, which are not accounted for here.
Changing the inner tube stiffness allows for deployment of the robot with a perfect FTL behaviour along paths with different shapes. In particular, increasing k b3 allows for the robot to withstand a higher magnetic torque while staying stable, high enough to deform the NiTi tubes. This results in a global modification of the path, with significant variations of the path position according to the magnetic field orientation.

Orientability
In this section, the orientability of M-CTR is analyzed numerically and demonstrated experimentally in a tip orientation scenario using open-loop control.

Orientability definition and investigation strategy
In the literature, the orientability of a robot at one point of its workspace p d is defined as the number of orientations achievable by the robot tip while maintaining its position at p d . In most related work, continuum robots are considered to have 6 d.o.f., which allows for full control of their tip pose (Chikhaoui et al., 2016;Li et al., 2017;Wu et al., 2017). In these works, the evolution of orientability across the robot workspace and the impact of actuation redundancy were studied. Orientability evaluation is based on solving simple kinematic models of continuum robot based on the constant curvature assumption (Chikhaoui et al., 2016;Li et al., 2017;Wu et al., 2017). As a result, phenomena such as torsion of the backbone and deformation due to external forces are neglected. The tube torsion has been considered in a preliminary work in Peyron et al. (2019b) for CTRs. Orientability evaluation is performed here for the first time accounting for the effect of external forces produced by the magnetic field and the gravity.
The study of M-CTR orientability is performed in three phases. First, an orientability map across the robot workspace is generated. A simplified forward static model is used that neglects torsion deformations and gravity effects. Second, orientability is evaluated numerically at several locations in the workspace by solving the inverse kinetostatic problem in (13) with our numerical framework, accounting for the first time for torsion deformations and gravity effects. Third, these results are used to achieve tip orientation control of the M-CTR prototype and to demonstrate its orientability experimentally.
In order to simplify the study of M-CTR orientability, a minimal set of actuation inputs among the 7 available variables was used. Since the tubes can be rotated at their base, the control of the axial tip orientation is trivial. Consequently, we are not interested in this rotation and only 5 d.o.f. are considered. Two actuation inputs are considered as constant during the tip rotation, which are α 3 and ΔL 31 . In particular, subsection 31 is completely deployed in order to maximize the angular displacement under magnetic actuation and α 3 is set to 0. The minimal set of actuation inputs is then [α 2 , ΔL 11 , ΔL 21 , B].

Orientability evaluation
As an initial step, the orientability map is generated following the method proposed in Wu et al. (2017). The forward kinematics of the M-CTR are derived assuming constant curvature, and neglecting gravity and deformations due to tube torsion. We also assume that the magnetic field only deforms subsection 31. Subsection 31 has a constant curvature u mag and it deforms in a plane of osculating angle f mag . In that case, the curvature along the backbone is The robot configuration is obtained by integrating (2, 3) with this expression of curvature.
The actuation space was discretized into 40 × 10 6 sets of actuation inputs using the Monte Carlo method, and accounting for the ranges of variation summarized in Table 4. The model was solved for each set of actuation inputs, leading to a set of robot tip positions and orientations. The workspace was then discretized into 4 × 4 × 4 mm cubic regions, of which the centres correspond to the desired position p d . The robot configurations which tip positions lie in a same cubic region were then collected, producing a set of tip orientations.
The set of tip orientations was finally used to evaluate the orientability at the corresponding p d . The tip orientations are represented on a service sphere of area A S located at the tip position. This sphere is discretized into a finite number of patches, representing the different possible tip orientations. The area A R of the service region is determined from the number of patches intersected by the backbone tangent at the tip. The orientability D(p d ) is then computed as The best robot orientability is obtained when D(p d ) gets close to 1.
The M-CTR orientability was evaluated for each cubic region, leading to the heat map presented on Figure 10. The workspace of the M-CTR is symmetric with respect to the (y 0 , z 0 ) plane, where the orientability is maximum. Therefore, we represent orientability in the (y 0 , z 0 ) plane only. The maximum orientability is equal to 0.98. It indicates that almost any tip orientation can be achieved with the robot tip. Three locations, labelled 1, 2 and 3, were chosen to evaluate M-CTR tip orientability. Location 1 was chosen at the centre of the workspace to facilitate experimental evaluation while maximizing the orientability. Location 2 and 3 were obtained by setting (α 2 , α 3 ) = (0, 0) to limit stability issues and by moving progressively towards the workspace boundaries. The service spheres at the 3 Table 4. Ranges of actuation inputs considered during the actuation space discretization.
Subsection lengths and curvature are, respectively, in m and m À1 , respectively. Figure 10. Orientability of the m-CTR across its workspace. The three locations considered for orientability assessment are represented with red squares. Configurations leading to these locations according to the constant curvature model are also represented.
locations are represented on Figure 11, and the orientability equals 0.93, 0.83 and 0.38, respectively.
Beyond the generation of the orientability map, the continuation method is used with the developed kinetostatic model to evaluate M-CTR orientability, by considering 1D rotations of the tip starting from an initial configuration of the robot. The initial configurations of each location are selected from the set computed previously, so that the robot tip is in the corresponding cubic region and the robot is in the (y 0 , z 0 ) plane. They are represented on Figure 10. We consider then the tip to be oriented around x 0 , y 0 and z 0 . The corresponding rotations are represented with quaternions Q x , Q y and Q z , respectively, so that where a x , a y and a z are the rotation angles. Moreover, the target around which the tip is rotated is considered at a distance jdj ¼ 3 mm from the robot tip, which corresponds to the thickness of the April tag used to measure the orientation. To perform the 1D rotation, the inverse kinetostatic model is solved from (13) by considering the desired tip orientation Q N and position p N where Q x|y|z is to be replaced with the desired rotation quaternion in (18), Q N,0 is the initial tip orientation and Q θ ¼ ½ cosðθ 1,N =2Þ 0 0 sinðθ 1,N =2Þ T represents the tip axial rotation of angle θ 1,N . The initial robot configuration computed with the constant curvature model is used as an initial guess for the inverse kinematics solution. The numerical framework is finally used to vary the rotation angles (a x , a y , a z ) and to compute the corresponding robot configurations. Robot configurations obtained at location 1 are presented in Figure 12. The initial M-CTR configuration is labelled 1 on the figure. Configurations (2,3), (4,5) and (6,7) are obtained when rotating the tip around x 0 , y 0 and z 0 , respectively. We obtain robot configurations that lead to the same tip position with different tip orientations. Since the initial tip tangent is almost aligned with y 0 , rotating about   this axis leads only to small backbone motions. The angular displacements obtained for the 3 locations are presented in Table 5. The results show significant differences of angular displacement between the axes of R 0 . The M-CTR tip orientability is not isotropic in the workspace. Interestingly, relatively small angular displacement is obtained about x 0 at the location 1, where the orientability is expected to be maximum. This is probably due to the impact of gravity and tube torsion, which could not be anticipated with simpler constant curvature forward kinematics. The major part of the obtained values is in the range [90°, 360°], which is higher than the angular displacements reported in the literature (Goldman et al., 2013;Simaan et al., 2009). According to these simulations, the M-CTR can be used to rotate a tool fixed at a given position about a specific axis with large amplitudes.

Orientation control
A demonstration of the M-CTR tip rotation capabilities was performed experimentally at location 1 in the workspace. We focus especially on tip rotation around x 0 and z 0 , since variations of the robot shape are the most important. We perform the tip rotations by following a pre-computed path in actuation space with quasi-static open-loop control. The block diagram of the control strategy is presented on Figure 13. The path was first planned in the task space, that is, for desired tip position p d and tip orientation angle a x , a z . Starting from the initial configuration, a z (resp. a x ) was increased with a constant step size towards its maximum value, moving from configuration 1 to 6 (resp. 1 to 3) on Figure 12. The maximum value of a x is chosen as 25°, in order to leave some margin with the maximum angular displacement evaluated in simulation. The maximum value of a z is limited to 90°to ensure the April tag detection. Then angle a z (resp. a x ) was then decreased towards its minimum value (configuration 6 to 7, minimum values of À90°, resp. 0°) and increased again to return to the initial tip angle (configuration 7 to 1). The step size was computed so that the entire rotation cycle is achieved in 40 steps. The path in task space was then expressed in the actuation space using the M-CTR inverse kinematic model. The actuation inputs computed during the simulation are interpolated according to the desired path. The desired magnetic field B d was then generated at p d with the CardioMag, and the actuation inputs α 2 , α 3 , ΔL 11 and ΔL 21 were sent to the M-CTR actuation unit. The steps were performed with a period of 4s, so that the robot reaches its steady state at the end of each step. Each cycle was repeated 5 times. Videos of the robot motion during the rotation cycles are presented in Extension 3 (see Appendix, Table 6). In the videos, the visual markers are removed and the open-loop trajectory is played continuously instead of step by step. During the rotation cycle, the robot tip pose was estimated using the camera A5 (see Figure 3). Details concerning the computation of the tip pose are provided in appendix C.
The mean values of the measured tip orientations over the 5 cycles are projected on service spheres shown in Figure 14. We observe that the points are globally contained in a plane, of which the normal vector corresponds to the tip rotation axis. It indicates that the tip is rotated around a fixed axis, as expected. The experimental angular displacement and the rotation errors were computed following the equations provided in appendix. We obtained angular displacements of Δa = 22.24 ± 0.01°and Δa = 136.36 ± 1.2°for the rotations around x 0 and z 0 , respectively. They represent 11.04% and 24.24% of the desired angular displacement around the two axes. Rotation errors for the two axes were ξ = 0.18°and ξ = 7.76°, respectively, which represent 0.81% and 4.97% of the angular displacement. Higher angular displacement and rotation errors for the rotation around z 0 can be attributed to friction which is higher since the middle tube is rotated.
The evolution of the tip displacement error ε is presented as a function of the desired tip rotation angle in Figure 15. Each point corresponds to the mean error value obtained for the corresponding tip angle over the 5 cycles. The mean deviation from these mean values is 0.22 mm. The mean position error for tip rotations about x 0 and z 0 are 0.8 mm and 6.3 mm, respectively, which represents 0.93% and 7.33% of the robot length. Considering that the standard errors of kineto-static models of continuum robots are up to 5% of the robot length, the errors obtained here seem acceptable for a first prototype.
As a result, the proposed open-loop control can be used to perform M-CTR tip rotation about the axes of R 0 while staying at the same position in the workspace. As predicted by the model, angular displacements surpassed 90°for the rotation around z 0 , showing the key advantage of M-CTRs in terms of orientability.

Stability
In this section, the stability modulation capabilities of M-CTR are studied by developing a modulation strategy and validating it both numerically and experimentally.

Magnetic field for stability modulation
Unstable phenomena of CTRs are similar to beam buckling, as demonstrated in Gilbert et al. (2016). When the tubes are placed in opposition, increasing the tube interaction length is analogous to increasing compression forces on a straight beam. This increase eventually leads to a critical interaction length (resp. a critical compression force) for which cardinality and stability changes are observed. The configuration where the tubes are placed in opposition is particularly important because it is the worst case scenario from a stability point of view. When tubes are rotated with respect to one another and their interaction length is gradually increased, this is the first configuration to become unstable, as shown in Hendrick et al. (2015a).
From the definition of stability, the M-CTR is stable when the tubes are placed in opposition, if it returns to the same configuration when a perturbation is applied. The robot can be ensured to keep this configuration by applying forces and torques on the backbone, which has been shown for a buckled beam with tendons (Berlin and Sussman, 1994), shape memory alloy (Baz et al., 1992) and piezoelectric actuators (Schaeffner et al., 2016). In the case of the M-CTR, the magnetic field is used to generate these restoring torques. For the analysis, the inner tube is considered as fully retracted, so that the magnetic torques apply directly on the pre-curved tubes. For the analysis, the inner tube is considered as fully retracted, i.e. its tip aligns with the tip of the middle tube and ΔL 31 ¼ 0mm. As a consequence, the magnetic torques apply directly on the pre-curved tubes. As the NiTi tubes are significantly stiffer than the inner tube, the robot does not experience the magneto-elastic instabilities mentioned in section 3. The robot is considered as contained in the (y 0 , z 0 ) plane when the pre-curved tubes are in opposition. Our strategy to generate the restoring torques is to consider a magnetic field in this plane as well, of which the relative angle with the magnet axis is f, as introduced in section 3, Figure 6. The buckling of the precurved tubes implies out-of-plane motion of the robot and of the magnet. This motion is represented with a rotation of angle γ about y b . Assuming that the magnetic field is homogeneous and following equation (5), the magnetic torque applied on the magnet is written in R b where c γ and s γ (resp. c f and s f ) denote the cosine and sine of angle γ (resp. f). We are especially interested in the component about y b , which induces out-of-plane motion. When f 2 [ À π/2, π/2], τ m .y b is of opposite sign to γ. The magnetic torque induces tip motion opposed to the out-of-plane perturbation and thus restores the initial configuration of the M-CTR. On the contrary, for f 2 [π/2, 3π/2], τ m .y b and γ have the same sign. The magnetic torque favours out-of-plane motion. Choosing a magnetic field aligned with the magnet when the tubes are in opposition (f = 0) leads therefore to the largest restoring torque and improves the M-CTR's stability. Choosing a magnetic field opposed to the magnet (f = π) deteriorates the robot's stability. These two magnetic fields are denoted B 1 and B 2 in the following.

Stability analysis
The impact of the fields B 1 and B 2 was determined by performing a stability analysis with the numerical method. In the literature, the stability of CTRs was studied by detecting stability and cardinality changes during tube rotation, or during tube deployment when they are assembled in opposition. Only the first scenario has been used to study stability experimentally. The second scenario consists in placing the robot in its buckling configuration, which was proven to be sensitive to small defects in experimental setups for M-CR (Singh et al., 2013). Therefore, we consider the first scenario.
First, the kineto-static model of M-CTR is solved with the tubes being placed in the buckling configuration to  obtain the magnet orientation. The tubes are placed in opposition ((α 2 , α 3 ) = (180, 0)°) and deployed until the critical interaction length is reached. In that way, the precurved tubes are at a stability limit, and the impact of the magnetic field on this limit can be directly observed. This critical length is computed with the criterion in Hendrick et al. (2015a) and equals 43.2 mm. The tube 2 is also fully deployed. The magnetic fields B 1 and B 2 are computed considering the magnet orientation and a magnitude jBj ¼ 80 mT with the relation (15). They are then applied and considered as constant during the tube rotation. Starting from the configuration where (α 2 , α 3 ) = (0, 0)°, the tube rotation α 2 is varied with continuation while α 3 is held constant.
The resulting diagrams for the case without magnetic fields and with the fields B 1 and B 2 are presented in Figure 16 in terms of the tip coordinate along x 0 . The case without magnetic fields reveals a S-shaped diagram with a vertical tangent at α 2 = 180°. As expected, the robot is at the stability limit for this tube rotation. In presence of B 1 , the diagram is similar but with a lower tangent steep angle at α 2 = 180°, indicating that the robot is stable during the full tube rotation. In presence of B 2 , two limit points (LP) appear during the tube rotation, the configurations between these points being unstable.
We can draw two conclusions from these simulations. First, the magnetic field has a significant impact on the robot stability, even if it does not deform the NiTi tubes significantly. This is due to the fact that in the buckling configuration, the robot behaviour is sensitive to small perturbations, produced here by the magnetic torque. Second, the strategy for selecting magnetic fields is efficient. The fields B 1 and B 2 stabilize and destabilize the robot, respectively.

Experimental analysis
We reproduced the scenario with the experimental setup. The tubes were initially aligned and deployed until the critical interaction length was reached. The critical interaction length of the prototype differs from the predicted    value. It was evaluated using a trial and error process, by slowly translating the tubes and performing full rotation of tube 2, until unstable phenomena were observed. We determined a value of 15 mm. The significant difference can be attributed to the friction between the tubes, which is not taken into account in the criterion of Hendrick et al. (2015a). The magnetic fields B 1 and B 2 were computed with this value of critical interaction length. They were generated with the CardioMag at the location corresponding to the tip position when the tubes are in opposition. The tube 2 was then rotated from 0°to 360°and in the inverse direction with 80 equally-spaced angular steps. Since the magnetic field is constant during the tube rotation and the stepper motors have a higher bandwidth than the CardioMag, the robot converges faster towards its steady state. Therefore, the steps were sent to the actuation unit with a period of 2 s. The rotation cycle was repeated 5 times. During the tube rotation, the tip position was measured using the stereoscopic vision system.
The resulting experimental diagrams are represented in Figure 17, and videos of the robot motion during tube rotation are provided in Extension 4 (see Appendix A, Table 6). In the videos, the tube rotation was performed continuously. Each point corresponds to the mean value of the tip coordinate over the 5 cycles. The mean deviation is 0.35 mm. When no magnetic field is applied, the robot becomes unstable near α 2 = 180°. When applying B 1 , the robot stays stable during the full tube rotation in the two directions. Hysteresis is visible on the corresponding diagram, which is due to friction between the tubes. When applying B 2 , the robot snaps with higher amplitudes and dynamics.
In conclusion, the numerical results were confirmed experimentally. The magnetic field has a significant impact on the robot stability, and it can be used to provoke or to prevent snapping. The strategy for computing the magnetic fields also showed its benefits experimentally, despite the presence of unmodelled phenomena such as friction.

Conclusion and perspectives
In this paper, M-CTRs are introduced and analyzed as candidates for performing MIS at millimetre-scales and in difficult-to-access areas. The interest of the M-CTR is assessed on a representative design of robot through numerical and experimental studies. These studies required the establishment of a generic kineto-static model, the application of a numerical framework and the implementation of a specific experimental setup. They led to a number of results concerning three important performance criteria in the targeted medical application, namely the FTL deployment capabilities, the orientability and the stability.
The robot can be deployed along planar paths with a perfect FTL behaviour by adjusting the magnetic field orientation during the deployment. A continuous set of paths with different curvatures can be followed by using magnetic actuation. In particular, the M-CTR can achieve distal radii of curvature down to 5.4 mm and angular displacements of ±160°all with a backbone which external diameter does not exceed 1 mm. This angular displacement is on par with the best performing continuum robots at that size. In addition, the global shape of the paths can be modulated by changing the inner tube stiffness.
The M-CTR can be used to orient a tool at its tip around a target. Numerical studies show that the tip orientability depends on the tip position in the workspace, and the axis around which the tip is rotated. For three locations in the workspace, angular displacements in the range [90°, 360°] are obtained which are comparable to and even greater than the angular variations reported in the literature. M-CTR orientability was also demonstrated experimentally by implementing open-loop control of the tip pose. Angular displacements of 136.4°were obtained, which exceed the state of the art. Moreover, open-loop control is promising as a control strategy for M-CTRs, given its simplicity.
Magnetic fields can be used to alter the stability of the M-CTR. In particular, the robot can be stabilized or destabilized by proper selection of the magnetic fields, which was demonstrated numerically and experimentally.
This work opens several perspectives that we intend to develop in the future. The M-CTR will be designed for assessment in medical applications such as olfactory cells inspection and middle-ear surgery. We will in particular focus on the optimization of the number of tubes, the tube properties and the development of accurate control for teleoperated or autonomous procedures. The robot will be composed of a navigation section, with tubes dedicated to the FTL deployment, and an exploration section dedicated to tip pose control as proposed for CTR in Mitros et al. (2018). The number, pre-curvatures and stiffness of the tubes in the navigation section will be optimized to conform with the specific geometries of nasal cavities and ear canals. The exploration section will be designed to fit the dexterity requirements. Also, the dimensions of the magnet will be optimized in order to fit the required 1~mm external diameter. This will be a challenge as the magnet volume should be large enough to produce a sufficient magnetic torque on the concentric tubes, which in turn depends on the tube properties and the application requirements. Potential solutions to this problem are to include the constraint on the magnet size in the tube optimization, and to consider other magnetic elements such as tubes composed of flexible ferro-magnetic material (kim et al., 2019). Concerning control, the proposed open-loop control shows promising results but its accuracy could still be improved. Closed-loop control will be developed by considering sensor integration and in-vivo constraints. Also, we will consider the use of inhomogeneous magnetic field and the field gradient as an actuation input, in order to increase the number of degrees of freedom and the achievable robot shapes. Finally, the fundamental properties of M-CTRs will be further investigated. With the current robot structure, the FTL deployment strategy will be validated and analyzed experimentally. The magnet at the tip will be replaced by magnets with different magnetization directions, in order to expand the deployment, orientation and stability modulation capabilities. On the long term, other architectures of M-CTRs will be studied. They will be generated by varying the number of tubes, the number and the location of the magnets along the backbone, and the mechanical, geometrical and magnetic properties of each element.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Agence Nationale de la Recherche within the Biomedical Innovation program (NEMRO ANR-14-CE17-0013), the Investissements d'Avenir (Robotex ANR-10-EQPX-44, Labex CAMI ANR-11-LABX-0004) and EUR EIPHI (Contract No. ANR-17-EURE-0002), the Swiss National Science Foundation through grant number 200020B_185039 and the ERC Advanced Grant 743 217 Soft Micro Robotics (SOMBOT).

Supplemental Material
Supplemental material for this article is available online.

ORCID iDs
Quentin Peyron  https://orcid.org/0000-0003-2628-1464 Quentin Boehler  https://orcid.org/0000-0001-7573-3135 for each robot state coordinate x k in x ¼ Â θ 1 … θ n u x u y Q T p T Ã ∂ ∂s where C q ¼ ½ C q1 C q2 and ha,bi denotes the scalar product between vectors a and b. This leads to the following equilibrium equation The relation between τ q and τ m is obtained by comparing the equilibrium equations with the static model of CTR developed in Lock et al. (2010). In that work, the special Cosserat rod equations are used, leading to the following expression of the backbone curvature where m and n are, respectively, the internal torques and forces due to distributed torques b τ and forces 0f, and v ¼ ½ 0 0 1 T . We conclude from the first and last line of equation (30) that the Lagrange multipliers λ q and λ p correspond, respectively, to the internal moment À b m and to the internal force 0 n . As a consequence, multiplying both side of the second line of equation (30) by 2L T q , and considering distributed magnetic torques b τ ¼ b R 0 0 τ m , leads to By comparing equations (31) and (29), we deduce that τ q relates to the magnetic torque 0τ m by the relation As a side note, by injecting equations (21) and (22) into equation (31), we can also verify the expressions of S q and S p C Appendix: Experimental measure of tip orientation ad angular displacement C.1 Computation of the tip pose. The front camera A5 (see Figure 3) gives a relative estimation of the robot tip pose with respect to its initial configuration. In order to compute the tip pose, the measured position and angular displacements are thus reported on the initial tip pose predicted by the kineto-static model. As a consequence, the position around which the tip is rotated and the axis of rotation differs from the measurements by a constant and unknown translation and rotation. We assume in the following that the mean position of the tip during the experiment, denoted p corresponds to p d . To enforce this, a constant translation t is applied at each measured tip position We also assume that the tip is rotated about the desired axis. The rotation axis measured by the front camera is evaluated by constructing the service sphere, finding the orientation plane which best fits the experimental points, and computing the normal vector z p to this plane. The plane is computed using the leastsquares method. A constant rotation can then be defined between the desired rotation axis, x 0 or z 0 , and z p . This constant rotation is applied to each measured orientation.
The camera pose is chosen in order to maximize the visibility of April tags during the rotation cycle. It was not possible to find a pose for which the April tag was always visible during the rotation about z 0 , due to the large angular displacement of the M-CTR tip. As a consequence, the camera pose was chosen to measure half of the rotation cycle, that is, for a z ≥ 0. The robot showed a symmetric behaviour during the other half of the cycle. The symmetric results for a z ≤ 0 are constructed, then presented in the following figure.
C.2 Angular displacement and rotation error. The angular displacement achieved during tip orientation and the rotation error are defined according to the points and vectors introduced in Figure 18. The angular displacement is computed by projecting the experimental points p in the orientation plane. The projected points are denoted p + . The points marking orientation limits, denoted p + 1 and p + 2 , are then extracted as presented on Figure 19. The angular displacement Δa is finally computed using the relation The rotation error, denoted ξ, is defined as the angle between the points on the service sphere before and after projection. It is defined by