Niki Baccile
email: niki.baccile@upmc.fr

Guillaume Laurent

Cristina Coelho

Florence Babonneau

Li Zhao

Maria- Magdalena Titirici

Structural Insights on Nitrogen-containing Hydrothermal Carbon using Solid-State MAS 13 C and 15 N NMR

Keywords: Nitrogen Doped Carbons, Hydrothermal Carbonization, Glucose, Glycine, 13 C solidstate MAS NMR, Double quantum NMR, sugar-amino acids interactions

 C and 15 N solid state NMR is used as the main and most effective characterization technique on nitrogen-containing hydrothermal carbons obtained from glucose and glycine. This study represents a model system for other types of nitrogen-containing hydrothermal carbons, which were shown to have interesting energy-storage properties (Zhao et al. Adv. Mater., 2010, 22, 5202).

These materials are obtained either from N-containing carbohydrates or from pure carbohydrates in the presence of natural amino-containing compounds such as proteins or aminoacids. In contrast to what is generally known for this model system, high molecular weight heterogeneous polymers that are formed when sugars and amino acids combine through the Maillard reaction (e.g., melanoidins), we found an extended nitrogen-containing aromatic network which is chemically bound to a polyfuran network known to be one of the main components of the biomass derived hydrothermal carbons. In contrast to the hydrothermal carbons obtained from pure carbohydrates, these types of N-containing materials have an increased level of aromatic character already present at 180°C, after the hydrothermal treatment.

Introduction

In the past 10 years, the hydrothermal treatment of biomass has gained an increasing interest in the field of material science. Many studies focus already on both applications 1 and fundamental 2 aspects motivated by the interest in producing carbonaceous powders with tunable sizes and surface properties directly from raw, processed and even waste biomass 3 . Recent review articles can be consulted for an overview on the progress regarding this topic over the past few years 4 . The most interesting point of this approach is undoubtedly the perspective of processing waste biomass (lignins, cellulose, starch) into valuable carbonaceous materials although the main research so far has been done using pure carbohydrates normally present in biomass composition. One of the main reasons for that is the fact that the chemical reactions which transform saccharides into hydrothermal carbons are extremely complex and they had first to be understood starting from simple model systems 5 . It is widely known that saccharides dehydrolize to form furans 6 and we have already pointed out the importance of such intermediate reactants in the formation of the final material. Different techniques like FT-IR 2b , XPS and XRD 7 have been used to identify the structure of the amorphous material but, so far, the best results regarding the final structure of such carbon materials were obtained by using advanced solid state NMR techniques 2a , which attested the presence of an abundant polyfuran core instead of the expected aromatic network, previously reported by other authors 2b .

One important advantage of these carbons is undoubtedly the possibility of modifying the surface 8 or the bulk of the carbonaceous network by introducing hetero atoms such as transition metals 1i, 1q, 9 or nitrogen 10 . Hybrid metal-carbon nanostructures have been obtained with potential application in catalysis. The introduction of N-atoms within the bulk or at surface structure of these materials is extremely interesting given the multitude of applications such as, for example, pHresponsive adsorbents 11 , supercapacitors 12 , fuel cell electrodes 13 , CO2 adsorbers 14 and materials with extremely high conductivities 15 . We have recently shown the possibility of making nitrogenrich carbon nanoparticles using a mixture of glucose and ovalbumin 16 , while glucosamine and chitosan 10 can also be used as simultaneous nitrogen sources. Preliminary results show a very good and highly selective CO2 uptake suggesting that such materials could be successful candidates as CO2 sequestration agents 17 or for use in supercapacitors 18 . For this reason we want to address the problem of structural identification no matter which natural nitrogen sources, like aminosugars, aminoacids or proteins, are used. We will use both 13 C and 15 N solid state NMR experiments to identify carbon and nitrogen sites and compare the structure of a nitrogen-containing material to a nitrogen-free one, pointing out the most important structural differences. In particular, a double isotope-enriched glucose-glycine mixture will be used as model system and compared to nonenrichable ones. In addition, we will show that the hydrothermal treatment of glucose and glycine, which is the classical reaction to make melanoidin resins, actually leads to a highly complex material with an unexpected structure of very different nature than typical melanoidin resins themselves, which were recently investigated by a combination of advanced 13 C and 15 N solid state NMR experiments 19 .

Here, we will use a multinuclear NMR approach to show that the combination of longer reaction times and higher temperatures under hydrothermal conditions drastically brings the widely known Maillard reaction between glucose and glycine to its extreme condensation products characterized with a high load of aromatic nitrogen-containing groups, which were not observed before in the final raw material. Interestingly, no matter the complexity of the nitrogen source, the final material seems to have very close structural features. The better knowledge of the nitrogen-containing carbonaceous structure should help understanding the energy-storage properties of these materials and, eventually, improve and control them through synthesis.

Experimental Section

Sample HC gly-C13-N15-gly preparation. 10 mL of a deionized water solution containing 1g of fully labeled D-Glucose-13 C6 (Aldrich, CAS:110187-42-3) and 0.4 g of glycine-15 N (Aldrich, CAS: 7299-33-4) is inserted into a glass vial inside a typical PTFE-lined autoclave system for hydrothermal reaction at 180°C for 24h. After reaction, the autoclave is cooled down in a water bath at room temperature. The obtained black solid powder is then separated from the remaining aqueous solution by centrifugation (7000 rpm for 20 minutes) and dried at 80°C in an oven under vacuum during 12 hours.

Scanning Electron Microscopy (SEM) images were acquired on a LEO 1550/LEO GmbH Oberkochen provided with an Everhard Thornley secondary electron and In-lens detectors.

Elemental analysis. Chemical analysis was performed on a (C, N, O, S, H) Elementar Vario Micro Cube.

NMR experiments. 1 H and 13 C solid-state Magic Angle Spinning (MAS) NMR experiments have been acquired on a Bruker Avance 300 MHz (7 T) spectrometer using 4 mm zirconia rotors spinning at a MAS frequency of νMAS = 14 kHz. 1 H and 13 C chemical shifts were referenced relative to tetramethylsilane (TMS; δ= 0 ppm). Details on one pulse MAS, Cross Polarization (CP) MAS, Inversion Recovery Cross Polarization (IRCP), CP MAS homonuclear single-quantum doublequantum (SQ-DQ) 13 C correlation experiments have been published elsewhere 2a . 15 N one pulse experiments were performed on a Bruker Avance 500 MHz (11.7 T) spectrometer using 7 mm zirconia rotors at νMAS = 5 kHz. An anti-ring pulse programme (ARING) with a 90° pulse angle of 7.70 µs was used. Recycle delay of 30 s and number of transients of 2000 were found to be optimal conditions. 15 N{ 1 H} CP-MAS experiments were done using a recycle delay of 3 s and a contact time of 2 ms, TPPM decoupling was applied during signal acquisition and the number of transients was 7200.

Triple resonance CPMAS 13 C{ 1 H, 15 N} heteronuclear correlation was recorded on a 700 MHz (16.4 T) spectrometer with a 3.2 mm zirconia rotor (νMAS = 22 kHz) using a triple 1 H-X-Y probe. A standard double CP pulse scheme (Sup Mat Figure1) was used. Typical acquisition values were: proton 90° flip angle, t 90°= 2.75 µs,15 N{ 1 H} and 13 C{ 15 N} contact times, tc1= 5 ms and tc2= 1 ms, 64 increments in the indirect dimension, and 1200 transients. Quadrature detection was achieved using the States method 20 .

Results

As pointed out already in the introduction, given the amorphous nature of HC and their chemical heterogeneity, it is impossible to perform a clear-cut structural study using characterization techniques like XRD, infrared spectroscopy or even XPS, whose poor resolution does not provide full pieces of information. Recently, we have used solid state NMR to solve the structure of hydrothermal carbon from glucose 2a and identify nitrogen sites in glucosamine-derived HC samples 10 . Similar experiments are very helpful whenever chemical enrichment is possible, which is seldom the case, especially if raw biomass is employed in the synthesis of HC. For this reason, we concentrated our attention onto the glycine-glucose derived HC as a model system. We employed Some preliminary morphological studies of nitrogen-doped hydrothermal carbon powders obtained from glucose and glycine are very similar to the previously synthesized carbon materials 5,2a . Sup Mat Figure 2 shows SEM micrographs of typical spherical microparticles as usually obtained from the hydrothermal carbonization process. The elemental composition shows that about 8.4 ± 0.5 w% of nitrogen is incorporated inside the material (C: 63.3 ± 1.0 w%, H: 5.5 ± 1.0 w%); this value is comparable with the ones obtained in the presence of the ovalbumin protein 16 as well as from glucosamine 10 . Even if elemental analysis shows that nitrogen is successfully incorporated in the sample, nothing can be told about the nature of the N sites and whether the glycine motif is kept intact as in the case of melanoidins 19 or underwent chemical modification. 13 C and 15 N solid state NMR experiments discussed below will bring a detailed answer to these questions.

The first point to be addressed is the validity of using the glucose-glycine mixture as model compound. Figure 1 shows the 13 C CP MAS response of a set of nitrogen-containing samples obtained from the hydrothermal treatment of and glucose/albumin 16 , 15 N-glucosamine 10 (Figure 1c) the mixtures 13 C-glucose/ 15 N-glycine (Figure 1d) compared with pure glucose (Figure 1a) . We have previously observed 2a that a contact time of 3 ms is enough to homogeneously excite the whole 13 C population, allowing a qualitative but reliable comparison among the relative intensities of the 13 C resonance peaks. A close comparison among the spectra shows very similar features: the amount of C=O groups (200-210 and 175 ppm) and aliphatic carbons (10-60 ppm) is significant unless glucosamine alone is used; the aromatic region is dominated by a large intense peak between 123 and 130 ppm (shaded region in Figure 1), which is generally absent when glucose (Figure 1a) alone is processed, and two peaks at 150 and 110 ppm, assigned to, respectively, C=C-O and HC=C-O groups. The resonance in the aromatic region around 125 ppm is very intense for the HC glu-C13-N15-gly system, indicating a larger aromatic network, whose existence is hard to prove in the HC glu sample. In fact, we have previously 2 shown that hydrothermal treatment of glucose provides a carbonaceous scaffold whose amorphous structure is dominated by an extended furan-type network. For this reason and given the possibility to use 13 The structural investigation of the nitrogen-containing carbonaceous material is done according to a previously published experimental line. 2a In that occasion, different NMR techniques were combined together to identify the main constituents of the carbonaceous scaffold obtained from hydrothermal treatment of glucose. Here, a similar approach is used; in particular, one-pulse 13 C and 15 N MAS NMR experiment allows the quantitative identification of all carbon and nitrogen sites and are presented here in Figure 2a and Figure 3a; cross polarization (CP) and inversion recovery cross polarization (IRCP) experiments give an unambiguous insight on protonated carbons and nitrogen sites; in particular, IRCP helps discriminating between CH, CH2 and CH3 groups 21 .

Figure 2b and Figure 3b present 13

Carbon-proton proximities

Chemical shift analysis combined with cross-polarization experiments were used to discriminate between carbonyl, aromatic and aliphatic groups (Figure 2b). In particular, variation of contact time in CP experiments (Sup Mat Figure 3c) nicely allows to discriminate between protonated and non protonated sp2 carbons while IRCP experiments give an insight on the degree of protonation of aliphatic carbon atoms (Sup Mat Figure 3 a,b).

Carbon-nitrogen proximities

The structural analysis was completed by one-pulse nitrogen-15 and 15 N{ 1 H} cross polarization NMR experiments (Figure 3). First of all, in good agreement with elemental analysis (8.4 ± 0.5 w%) these experiments clearly show that nitrogen is well introduced inside the carbonaceous network and it is surrounded by different carbon environments. Chemical shifts at -75 and -250 ppm indicate, respectively, sp2 and sp3 nitrogen-carbon bonds.{ 1 H} 15 N CP (Figure 3b) only shows protonated nitrogen sites, whose chemical shift range (-200/-300 ppm) indicates the presence of amide or pyrrole groups 23 ; interestingly, amine groups which are initially introduced via the glycine molecule, have massively reacted since the NH2 resonance at -330 ppm is very low; additionally, the 15 N one pulse experiment (Figure 3a) shows non-protonated aromatic, pyridinic and pyrazine-like groups (0-150 ppm). This preliminary analysis shows very interesting features which make this material very different from the hydrothermal carbon obtained from pure glucose 2a ; in particular, the introduction of nitrogen induces a systematic aromatization of the carbonaceous network.

Deeper insights on the type of C-C and N-C bonds can be obtained from the 13 C{ 1 H, 15 N} triple resonance cross polarization NMR experiments shown in Figure 4. Here, an initial 15 N{ 1 H} cross polarization step at an optimum contact time (tc1= 5 ms) excites all protonated nitrogen atoms (Figure 3); then, polarization is transferred to 13 C nuclei via a second contact time (tc2= 1 ms) which was chosen long enough to detect a signal with a reasonable signal-to-noise ratio but short enough to probe the shortest possible C-N distances. The use of the proton bath to excite the 15 N spins limits the accessibility to solely N-H groups (-150/-300 ppm) with respect to the entire nitrogen population which also includes aromatic sites (0/-150 ppm) (Figure 3a). Any attempt to run a complementary 15 N{ 1 H, 13 C} experiment failed. For this reason, we stress the fact that the results obtained from this experiment only concern protonated nitrogen environments. The 2D 13 C{ 1 H, 15 N} map in Figure 4 abundant species in the material, followed by CHx-NH species. On the contrary, no specific correlation is found between (protonated) nitrogen and carbonyl groups, suggesting that amides are not relevant. This result is quite interesting as it suggests that the signal between 130.0 and 136.0 ppm mainly relates to pyrolic species.

Carbon-carbon connectivities

Two double quantum single quantum correlation experiments were recorded in order to select different types of carbon environments. In the first one (Figure 5a), a 13 C{ 1 H} contact time (tc) of 35 µs was employed before the excitation of DQ 13 C coherences while in the second one (Figure 5b) the tc was much longer (3 ms). The shorter the tc, the higher the selectivity to detect protonated carbons. In both cases the SQ-DQ 13 C excitation and reconversion time was the same (285.7 µs)

and chosen in order to observe only the direct carbon-carbon bonds 2a,22 .

At short contact time (35 µs -Figure 5a), most of the constitutive chemical units can be identified. In particular, ( 1) and ( 1a) refer to the classical furane unit connected to aliphatic groups via the α-carbon atom while (4) shows keto-aliphatic and intra-aliphatic bonds. These species were already observed on the glucose-derived hydrothermal carbon and are quite typical for these materials 2a . Additionally, more specific bonds can be identified in this nitrogen-containing sample.

First of all, the on-diagonal cross-peak at 127.5 ppm (2) now dominates the 2D map (Figure 5a,b); in particular, its on-diagonal position suggests the existence of abundant aromatic units while the off-diagonal (2a) peaks at δ= 151 and 109 ppm suggest the presence of possible phenols. The crosspeak at 127.5 ppm also indicates the existence of possible pyrazine or pyridine units that can be related to the existence of the broad peak at -75 ppm in the one-pulse 15 N spectrum in Figure 3.

Finally, specific furanic-aromatic bonds could also be proposed (2a). Pyrrole-type compounds, proposed after the 13 C{ 1 H, 15 N} HETCOR experiment, could be justified by cross-peaks (3) while NH-CH groups should be observed in (5). At longer contact times, the SQ-DQ experiment does not bring further information on the material core-network but it shows the existence of newly appearing and self-correlating cross-peaks, as outlined by the connecting black arrow. In particular, the COOH group at 176 ppm is readily connected to a cheto-bond (202 ppm) through a chemical group resonating at 106 ppm, whose exact attribution constitutes a harsh task at the moment. These cross-peaks probably identify isolated or grafted molecules whose lack of signal at short contact times rather depends on their mobility (less efficient cross polarization transfer) than on the lack of protonated environment; in fact, CP and IRCP experiments show that carbon peaks at 70.6 ppm and 106.7 ppm are indeed protonated. In Figure 5b we propose the possible molecular structure, which, according to the exact nature of R and R' groups, may be related to organic acids-derivatives (e.g., dihydroxyfumaric acid, tartaric acid or dehydroacorbic acid).

After combining the pieces of information obtained from one-pulse, CP, IRCP, HETCOR 15 N-13 C and SQ-DQ 13 C experiments presented above, the deconvolution of the one-pulse 13 C spectrum (Sup Mat Figure 4) of the HC glu-13 C- 15 N-gly sample (Figure 2a) shows not less than 25 different carbon species, whose list of resonances, chemical shift values and tentative attribution is listed in Supp Mat Table 1.

If compared with literature, the glucose/glycine system is well-known for its reactivity in hydrolytic medium undergoing Maillard reaction, which is at the origin of the non-enzymatic browning in food chemistry and in the synthesis of melanoidin resins 19,24 . Nevertheless, the conditions employed in the HC process higher temperatures (temperatures above 160°C, selfgenerated pressures -10-30 bars -and residence times between 12 and 24 hours) produce a different material: typical melanoidins keep up to 60% of glycine molecules intact and formation of pyrazines and pyridines is generally not observed while 13 C and 15 N NMR experiments on glycinederived HC's show no trace of leftover glycine and the overall amount of NH2 groups has almost quantitatively reacted further. Maillard reaction cascades are the most probable pathways that transform sugar, sugar derivatives and amines into cyclic compounds but the formation of modified pyrazines 25 , furanones, pyranones 26 and pyrroles 24b is largely documented more as degradation compounds rather than as network-forming bricks.

To go further, one should not forget that the main degradation product of glucose is hydroxymethylfurfural (HMF), which is the main reactive compound in the formation of hydrothermal carbon 2a,5 and that spectroscopic signature of furane rings is still abundant in the final nitrogen-containing HC material. The C site in α-position in furan rings and their derivatives can easily polymerize via radical, cationic or anionic reactions 27 ; nitrogen-containing derivatives of HMF can also be synthesised as precursors for linear and branched polymers 28 but the chemical processes behind are still quite unclear because the reaction between glucose, furans and glycine does not stop at the simple glycine incorporation and methods to follow the reactions in-situ are at the moment not developed yet.

Conclusions

Nitrogen-containing carbons constitute a new class of interesting energy-storage materials. The use of renewable resources in combination of hydrothermal treatment also constitutes a cheap way to introduce nitrogen within a carbonaceous scaffold. Nevetheless, the complexity of the reaction mechanism and of the final structure itself constitute a barrier to the understanding, tought the improving, of these materials.

Treating natural nitrogen-containing molecules such aminoacids or proteins with glucose or aminated saccharides such as chitosan or glucosamine under hydrothermal conditions provide the formation of nitrogenated materials up to 8-10 w%. Here, we have carried out a specific structural study on the glycine/glucose system using solid state NMR and isotopically-enriched ( 13 C and 15 N) sources. One-dimensional 13 C experiments (one-pulse and CP MAS) showed that all hydrothermal carbons have very similar structures no matter the complexity of the carbon source. Then, the presence of extended aromatic domains in coexistance with the expected furan network is documented and constituted the main feature of these materials.

Finally, the combination of one-dimensional (CP, IRCP) and two-dimensional (single quantum double quantum 13 C, double CP 13 C{ 1 H, 15 N}) experiments help drawing the carbon-to-carbon and carbon-to-nitrogen chemical arrangements. We were able to propose several possible structures which include the presence of pyrazine motifs connected to furans and we showed the quasiabsence of NH2 groups. Interestingly, the structure of such a material is very atypical for the Maillard reaction between glucose and glycine, which normally react to give melanoidin resins for which only the byproducts were shown to contain aromatic nitrogens but never the final solid.

There is no trace left of the initial glycine and glucose. Instead they completely rearranged into an extended pyrazo-furanic matrix. The exact mechanism is still ambiguous but possible reactions between the alpha-carbons of furans and the pyrazines is a possible explanation.

  C-enriched glucose and 15 N-enriched glycine, we show hereafter a full structural study of their hydrothermal derivative with our discussion based on the C-H, C-N and C-C local environments.

Figure 1 -

 1 Figure 1 -13 C CP MAS NMR experiments (tc= 3 ms) recorded on pure C13-glucose-derived carbons (HC glu) and several nitrogen-containing hydrothermal carbon materials: (a) sample HC glu; (b) HC glu-albumin; (c) HC glu-N15-glucosamine; (d) HC glu-C13-N15-gly. In (a) glucose is used alone but it can co-react with (b) ovalbumin 16 and (c) 15 N-glucosamine. 10 In (d), 13 C-glucose reacts in presence of 15 N-glycine. Note: the better signal-to-noise ratio of HC glu and HC glu-C13-N15-gly derives from 13 C enrichment.

  Figure 2b and Figure 3b present 13 C and 15 N CP MAS experiments performed at a given contact

Figure 2 -

 2 Figure 2 -Comparison between 13 C MAS NMR experiments acquired on the HC glu-C13-N15-gly system under the following conditions: a) one-pulse; b) CP performed at tct= 50 µs; c) projection from the SQ-DQ 13 C experiment at tc= 3 ms; d) projection from the HETCOR 15 N-13 C experiment.

Figure 3 -

 3 Figure 3 -15 N one-pulse and CP-MAS (with 2 ms contact time) spectra of HC glu-C13-N15-gly and 15 N projection of the 13 C{ 1 H, 15 N} HETCOR experiment.

Figure 4 -

 4 Figure 4 -2D map showing the double CP 13 C{ 1 H, 15 N} experiment on the HC glu--C13-N15-gly system and obtained with pulse sequence shown in Sup Mat Figure 1. Contact times for each cross polarization step are given in the figure. Insight on the possible chemical groups are reported in Figure 5. Note that here "R" is just a general notation to refer to C or H elements.

  provides a direct fingerprint of the nitrogen-to-carbon proximities. The 15 N-filtered 13 C projection shows two main regions, one with two peaks at 130.0 and 136.7 ppm and the other at 70 ppm while the 15 N projection shows three resonances at about -195.0, -230.0 and -250.0 ppm. As expected, these values are comparable with simple 15 N{ 1 H} CP experiments in Figure 3. A deep look at the 2D map shows that C=CH-NH groups are the most

Figure 5 -

 5 Figure 5 -HC glu-C13-N15-gly. SQ-DQ 13 C maps obtained for equal double quantum excitation and reconversion times (285.7 µs). Contact time before double quantum coherences excitation is tc= 35 µs (a) and tc= 3 ms (b).

N-glycine and13 C-glucose both as spin-1/2 nitrogen and carbon sources.


Paul Mignon
email: paul.mignon@imag.fr

Philippe Poignet
email: philippe.poignet@lirmm.fr

Jocelyne Troccaz
email: jocelyne.troccaz@imag.fr

Beveled-Tip Needle-Steering Using 3D Ultrasound, Mechanical-Based Kalman Filter and Curvilinear ROI Prediction

Keywords: Needle steering, needle detection, ultrasound, Kalman filter

This paper introduces a new robust 3D ultrasound needle detection approach integrated in a 3D needle steering system associated to a real-time path planning. The robustness of an existing algorithm is improved by limiting the needle detection to a curvilinear region of interest (ROI) using a novel mechanical-based prediction model. This linear model is also used in a Kalman filter to reduce detection noise and reject false detections. These two improvements drastically increase quality of our feedback. Finally, the 3D needle steering system is able to reach a target in phantoms with a maximal error of 0.8 mm without obstacle and 1.6 mm with obstacle.

I. INTRODUCTION

To reduce patient injuries and guarantee procedure success, needle-involved clinical procedures require very fine precision during insertion. This problem has motivated many research for robotizing needle insertion. Among them, needle steering -i.e. fine control of needle curvature for the execution of curved trajectories toward a target -has received increasing interest over the last decade.

A. Needle Steering and Path Planning

Needle steering can be performed using active or passive needles. Active needles have the ability to act on their own deformation while passive needles use external forces.

Passive needle steering can be performed by applying forces and torques on needle base or by using beveledneedle tip asymmetry to bend the needle in tissues.

In [START_REF] Dimaio | Needle Steering and Model-Based Trajectory Planning[END_REF], DiMaio introduces a 2D finite element model of tissue and needle deformations to guide a needle tip by acting on its base. A path planner based on potential fields enables reaching a specific point. Glozman [START_REF] Glozman | Flexible Needle Steering and Optimal Trajectory Planning for Percutaneous Therapies[END_REF] reduces the computation time by using a needle/tissue deformation model based on beam theory and virtual springs instead of finite element.

The needle base manipulation method has still an important drawback: needle base movement influence on needle tip tends to decrease when the needle is deeply inserted. Thus, forces applied to the needle base need to be increased with insertion depth and could generate traumas on the patient.

The second method does not suffer from this drawback and uses the deformation due to the force applied by the tissue on the beveled needle tip.

Several designs have been proposed to improve bevel-tip needle deformation or reduce patient injuries; a review of these designs can be found in [START_REF] Van De Berg | Design Choices in Needle Steering -A Review[END_REF]. A model of needle deformation is still necessary to plan needle insertion path. Webster, in 2006 [4], developed a simplified kinematic model which approximates needle trajectory by a circular path. The arc radius depends on tissue and needle mechanical properties. Duty-cycling (DC) methods, introduced in [START_REF] Engh | Flexible Needle Steering System for Percutaneous Access to Deep Zones of the Brain[END_REF], allow to modify this radius by alternating pure insertion and insertion with rotation.

Some planning algorithms based on Markov Decision Process [START_REF] Alterovitz | Steering flexible needles under Markov motion uncertainty[END_REF] or Rapidly-exploring Random Tree (RRT) [START_REF] Bernardes | Adaptive path planning for steerable needles using dutycycling[END_REF] use the previous kinematic model in 2D. Later, geometric inverse kinematic [START_REF] Duindam | Three-dimensional motion planning algorithms for steerable needles using inverse kinematics[END_REF], stochastic model [START_REF] Park | Diffusion-based motion planning for a nonholonomic flexible needle model[END_REF] or RRT [START_REF] Patil | Needle Steering in 3-D Via Rapid Replanning[END_REF] have also been used for 3D planning. Patil in [START_REF] Patil | Needle Steering in 3-D Via Rapid Replanning[END_REF] used the 3D RRT rapid-replanning strategy from [START_REF] Bernardes | Adaptive path planning for steerable needles using dutycycling[END_REF] to steer a needle while avoiding obstacle and perturbations. This fast position feedback allows to achieve closed-loop insertion and reduce errors in control and planning. For this reason, we use the same approach to guide a needle in 3D.

B. Needle Detection

Fluoroscopy [START_REF] Glozman | Image-guided robotic flexible needle steering[END_REF] or MRI [START_REF] Patel | Closed-loop asymmetric-tip needle steering under continuous intraoperative MRI guidance[END_REF] are commonly used to detect needles. However, Ultrasound (US) stays the safer, less expensive and more practical imaging modality for intraoperative purpose.

In [START_REF] Abayazid | Ultrasound-guided three-dimensional needle steering in biological tissue with curved surfaces[END_REF], Abayazid uses US imaging based on a 2D probe mounted on a motorized support to detect and steer a needle tip using the 3D RRT replanning approach. This "frame by frame" tracking method makes detection more robust but needs the needle to be roughly normal to the US plane.

In [START_REF] Adebar | Recursive estimation of needle pose for control of 3D-ultrasound-guided robotic needle steering[END_REF], Adebar extends this "frame by frame" segmentation approach to a 3D motorized probe. Color-Doppler mode is associated with an ultrasound buzzer at the needle base, to improve needle detection. An Unscented Kalman Filter (UKF) filters the detection using the non-linear unicycle model to predict needle tip displacement between two images.

We previously developed a detection method based on color-Doppler and needle rotation during insertion [START_REF] Mignon | Using Rotation for Steerable Needle Detection in 3D Color-Doppler Ultrasound Images[END_REF]. The detection is based on 3D reconstructed volume from a commercialized endorectal 3D probe. Unlike [START_REF] Abayazid | Ultrasound-guided three-dimensional needle steering in biological tissue with curved surfaces[END_REF] and [START_REF] Adebar | Recursive estimation of needle pose for control of 3D-ultrasound-guided robotic needle steering[END_REF], the needle orientation in relation to the probe is therefore not constrained by the detection algorithm.

We chose to use a Random Sample Consensus (RANSAC) algorithm inspired by [START_REF] Uhercik | Model Fitting Using RANSAC for Surgical Tool Localization in 3-D Ultrasound Images[END_REF]. This algorithm can be successfully used to extract polynomial (Bezier) curves approximating the needle shape from noisy environment.

In [START_REF] Pierre Chatelain | Real-time needle detection and tracking using a visually servoed 3D ultrasound probe[END_REF], Chatelain uses this reconstructed 3D approach and a Kalman Filter (KF) to improve needle detection given by a similar polynomial curve RANSAC. This KF uses the position continuity of each control point as the model for the needle displacement. However, as the needle is manually inserted, no prior information could be added to the filter.

Later, the limitation of the detection to a cylindrical ROI make it more robust [START_REF] Zhao | Biopsy Needles Localization and Tracking Methods in 3D Medical Ultrasound with ROI-RANSAC-KALMAN[END_REF].

C. Contributions

In this paper, we present a novel needle tracking approach based on a KF combined with a novel 3D mechanical model inspired by [START_REF] Roesthuis | Mechanics of needle-tissue interaction[END_REF]. Along with mechanical properties, this model also integrates robot information such as insertion speed and orientation of the needle to have a better estimate of the whole 3D shape deformation of a Bezier curve approximating the needle.

This mechanical model also increases the precision and robustness of our detection giving a predicted curvilinear ROI using the last needle detected in the 3D US image (see III.B).

This new approach is compared to the similar kinematic approach proposed in [START_REF] Pierre Chatelain | Real-time needle detection and tracking using a visually servoed 3D ultrasound probe[END_REF] by tracking a needle on 3D image sequences.

Finally, a 3D needle steering system merges these contributions using an existing brachytherapy robot prototype [START_REF] Hungr | A 3-D Ultrasound Robotic Prostate Brachytherapy System With Prostate Motion Tracking[END_REF] and a RRT replanning algorithm [START_REF] Bernardes | Adaptive path planning for steerable needles using dutycycling[END_REF]. Steering results obtained on phantoms with and without obstacle are summarized in section IV.C.

II. MECHANICAL-BASED KALMAN FILTER

A. Kalman Filter process

Standard needle tracking methods based on noisy images processing may lead to false detections. For that reason, we use a KF to filter the result.

The discrete KF is based on the knowledge of variables (state) evolution. This evolution is translated into a linear equation. The KF assumes that a process noise and a measurement noise, both Gaussian, corrupt respectively the prediction model and the measurement. Filtered state is computed by an optimization between the prediction and the measure. For more information about KF, see [START_REF] Kleeman | Understanding and applying Kalman filtering[END_REF].

The models described in Section II.B and II.C will be used as prediction, separately. In the mechanical KF the state represents the Bezier curve control point coordinates ⋮ with . The measurement noise, about 1 mm, corresponds to the standard deviation of each control point estimated from 100 segmentations. The process noise, around 0.1 mm, is computed from comparisons between simulations and real needle insertions.

The filtered needle segmentation corresponds to the updated estimate of the control points. The ROI is created with a prediction computed from the models (see III.B). In addition, we use a "validation gate" to reject high false detections (see [START_REF] Kleeman | Understanding and applying Kalman filtering[END_REF] for more details).

B. Mechanical-based Model

As mentioned in Section I.A, some needle deformation models use virtual springs to simulate needle-tissue interactions. For instance, in [START_REF] Roesthuis | Mechanics of needle-tissue interaction[END_REF], one of these models was developed to estimate needle deformation caused by needle asymmetry. This model considers needle as a piece-wise polynomial curve in 2D. Polynomial coefficients are computed by minimizing the energy accumulated in the system and the work produced by the beveled needle tip force.

We use a similar model and apply it on a 3D Bezier curve resulting from the detection step introduced in Section I.B. The deformation is limited by the polynomial aspect of the Bezier curve. However, simulations show that cubic or quadratic Bezier curves are sufficient to describe slight needle deformations.

Being given a needle at a time t described by a Bezier curve and its control points , every curve points could be expressed by the relation (1):

1 ⋯ ( 1 )
where the matrix depends on the order of the Bezier curve N and s is a parameter between 0 and 1. Characters in bold represent vectors and matrices.

Once the needle appears in the image, its visible base (i.e. the needle point in the image nearest to the needle base) is assumed to be fixed. It means that the first control point does not move. The visible tip moves in the direction tangent to the needle tip. Tip and base speeds are supposed to be equal. Thus we suppose there is no axial compression in the needle.

The energy of the system could be divided into the three following components:

• The energy due to needle deformation, approximated for small deformation by (2) :

. ( 2
)
where E is the Young's modulus of the needle, I is the second moment of area computed from needle section, L is the needle length approximation and represents the second derivative of in relation to s.

• The energy accumulated in the virtual springs which simulates the needle-tissue interaction. It could also be approximated for small needle deformations by (3):

‖ ‖ ( 3
)
where K is the virtual spring stiffness per length, t represents the current instant and t-1 the previous instant. is the distance traveled by the point during the insertion step and is expressed by where is the insertion speed and dt is the period between two steps. Including that term allows to consider the needle insertion movement directly in the virtual spring energy.

• The work of the force F acting on needle tip:

• 1 1 ( 4 )
The direction of the vector F is determined by last needle direction and needle roll angle given by the robot. Thus, we neglect torsion and suppose that needle base and tip roll angles are the same.

The approach consists in computing optimal control point coordinates which minimize those energies. To that purpose, we compute the system of differential equations [START_REF] Engh | Flexible Needle Steering System for Percutaneous Access to Deep Zones of the Brain[END_REF].

, , , ,

Differential equations ( 5) lead to a linear system of 3xN depending on current ( ) and previous ( ) control point coordinates.

It results in a matrix relation between current and previous control points and the force applied on the needle tip. This expression defined in [START_REF] Alterovitz | Steering flexible needles under Markov motion uncertainty[END_REF]. [START_REF] Alterovitz | Steering flexible needles under Markov motion uncertainty[END_REF] where and depends on E, I , , δs, L K and the order of the Bezier curve.

depends only on F. Using inverse or pseudo-inverse, the current control points coordinates are obtain directly from the previous ones.

The cutting force F is defined in the literature as one of the three forces applied to the needle during insertion in the tissues [START_REF] Okamura | Force modeling for needle insertion into soft tissue[END_REF] and is supposed constant at each step. F is modulated according to needle insertion DC (value between 0 and 1 corresponding to rotation time over insertion time) by the following relation 1 . corresponds to the true force applied on needle tip.

The virtual spring stiffness K reflects needle-tissue interaction. It therefore depends on the tissue and needle mechanical properties.

Since this model is linear, there is no need for a UKF as used in [START_REF] Adebar | Recursive estimation of needle pose for control of 3D-ultrasound-guided robotic needle steering[END_REF]. This model also computes the entire needle shape deformation, not only needle tip displacement, and is therefore more precise for ROI selection (see III.B).

C. Kinematic model

In order to show the added value of the presented approach and because of the similarity in the detection processes, we compare our KF to a more standard model [START_REF] Pierre Chatelain | Real-time needle detection and tracking using a visually servoed 3D ultrasound probe[END_REF]. This model is based on the control point kinematic relation between point position and speed. To have a more relevant comparison, we will add to this model information from the robot using needle tip velocity. A virtual measure of velocity vector is created from the last needle direction detected and the insertion speed given by the robot.

III. NEEDLE STEERING APPROACH

A. Control Scheme

The strategy used to steer the needle consists in two stages: the low-level control loop and high-level control loop (see Figure 2). The former concerns the robot. The latter is composed of the US acquisition, needle segmentation, the novel KF and the path planning. The path planning returns a path composed of circular arcs and is transformed into control inputs by the DC algorithm.

The insertion and rotation speeds are defined by the user. Thus, the DC algorithm takes as input the DC value and the total insertion distance. It returns several insertion steps controlling the insertion distance and rotation angle for each step. To avoid involuntary change of orientation, the rotation angle must be a multiple of 360°. Thus the rotating time has a lowest value which also limits the total insertion time step. To avoid the influence of starting and stopping motor transient phases, we also keep rotating time superior to 500 ms.

B. Selection of the Region of Interest

The mechanical or the kinematic model are used A specific thickness (called side margin), corresponding to apparent needle diameter, is added to the ROI (see Figure 4). A margin at needle tip could also be set by the user in some cases, for example when the model used to predict needle displacement does not integrate needle progression.

ROI based detection using RANSAC and KF has already been investigated [START_REF] Zhao | Biopsy Needles Localization and Tracking Methods in 3D Medical Ultrasound with ROI-RANSAC-KALMAN[END_REF] but only with a cylindrical ROI which does not evaluate needle deformations: the prediction model uses tip velocity measurement relying on speckle tracking. Here, we use robot information as well as mechanical-model and avoid potential errors from speckle tracking velocity detection.

C. Experimental Setup

PROSPER is a robot developed for brachytherapy procedures [START_REF] Hungr | A 3-D Ultrasound Robotic Prostate Brachytherapy System With Prostate Motion Tracking[END_REF]. It has a positioning module with five motors to pre-position and orient the needle tip outside the patient and an insertion/rotation module with two motors for needle steering inside tissues.

An Ultrasonix RP acquires 3D images from an endorectal US probe and sends it to a computer. A GPU-based 3D reconstruction algorithm inspired by [START_REF] Lee | Intensity-based visual servoing for nonrigid motion compensation of soft tissue structures due to physiological motion using 4D ultrasound[END_REF], using 2D Cartesian slices instead of RF signal, compensates for transducer constant velocity in acquisition. The framerate depends mainly on the 3D probe swept angle, we therefore limit it to 30°. The final 3D resolution is 0.27x0.27x1.37 mm farthest from the probe (50 mm) and the framerate is around 1 volume per second.

The user selects target position in a pre-insertion US image and can also locate obstacles corresponding to critical areas. The needle is inserted enough to make it visible in the image and is first located by the user. Then our automatic segmentation process is used to track the needle and planning is performed by RRT. These steps are repeated for each incoming image and allows to compensate for segmentation, path planning or control errors and possible robot compliance.

In order to have reusable phantoms for needle steering experiments, we selected two types of phantoms: a paraffin gel associated with a Vybar polymer (massic proportion of 25%) and PVC based phantoms. The paraffin phantoms are preferred because of the low manufacturing complexity. Moreover, the polymer assures the control on the phantom stiffness and simply heating it allows to recover its initial shape and stiffness.

IV. RESULTS

A. Model Parameters Definitions

In order to have an estimate of the cutting force in our future experiments, four different needles (2 biopsy needles, 1 brachytherapy needle and 1 Nitinol needle) have been mounted separately on a RAVEN II robot equipped with a force sensor. This sensor measures the force applied on the needle base during its insertion into PVC phantoms. Each phantom has a different apparent Young's modulus, previously measured by LASTIC [START_REF] Schiavone | LASTIC: A Light Aspiration Device for in vivo Soft TIssue Characterization[END_REF]. The norm of is extracted from a graph of the total axial force against the insertion depth.

Figure 3 represents the axial force detected on the standard brachytherapy needle (stainless steel, 17 gauges) base during the insertion in four phantoms with different elasticity. The first peak corresponds to the puncture force necessary to penetrate in the tissue. The slope after the peak corresponds to the regular friction force increase, proportional to the insertion depth. could be computed from the friction force slope and the puncture point. The same process is applied to every needle. Results are given in Table 1. The virtual spring stiffness is then estimated from a comparison between needle insertion simulations and experimental results. The four needles have been inserted in the different phantoms and the stiffness corresponding to the best fit is returned. The cutting force and spring stiffness are summarized in the Table 1.

B. Needle Tracking Using Kalman Filter and ROI

The noise in our detection comes mainly from RANSAC non-determinism and from US intrinsic noise in the image. Table 2 summarizes errors for static needle detections in 26 acquisitions in 3D. It presents shape error (mean along the needle shape), needle tip position and direction errors. Even if the standard deviation error is already very low without KF, we can see that KFs reduce the standard deviation of the errors and thus increase the precision of the detection. However, the KFs do not reduce mean errors. Thus they do not affect detection accuracy during static detections. The purpose of the following experiment is to compare tracking performance in three cases and to determine which model is the most appropriate.

In the first case, no model is applied, the ROI is then determined from last needle detection. The second case involves the kinematic model and the last case uses the mechanical model. Results have been compared to manual segmentations. A 0.5 mm nitinol needle is pre-curved at 3 mm from its tip to increase bending in the paraffin phantom. The needle is inserted at 1 mm/s. F and K for a paraffin gel phantom (80 kPa) were returned from simulations and are 1 N and 320 N/mm respectively. Side margin is set to 3.2 mm corresponding to the apparent needle thickness in the acquired volume. When no model is used, a 1 mm tip margin is added to manage needle progression. For the two cases using model, no tip margin has been added, the margin is therefore the result of the prediction.

In practice, the very low framerate (1Hz) of 3D acquisition increases the model impact for ROI selection. To emphasize further the difference between models, the needle performs a double-bend path.

Figure 6 shows tracking results of a double-bend path and the profile of the insertion speed.

The tracking without prediction model loses the tip when the needle stops for orientation change. It is caused by the presence of a constant tip margin and the hyperechogenicity of the phantom.

The kinematic model keeps tracking the needle tip during stopping time thanks to the insertion speed given as input. However, it loses the needle tip when the insertion restarts due to the change in needle orientation.

The mechanical model allows to adapt tracking to the change in orientation by using needle tip force vector as input. It allows the needle to be tracked for the rest of the insertion.

C. Needle Steering Experiments

This illustrates the efficiency of the new 3D US tracking. This section presents the 3D needle guidance toward a virtual target. For this purpose, we integrate all steps mentioned in section III.A: tracking with mechanical KF, rapid 3D RRT re-planning and DC control. Experiments are performed in two different phantoms with different elasticities (180 kPa and 80 kPa). For each phantom, needle steering is achieved with and without obstacle avoidance (see Figure 7). We use the 0.5 mm pre-curved nitinol needle and an insertion velocity of 0.5 mm/s. The error is the distance between the virtual target and the final needle tip position obtained by the automatic segmentation.

Errors found without obstacle are 0.5 and 0.8 mm. With one obstacle, errors were 1.6 and 0.7 mm. However, let us remind that the detection accuracy is around 1 mm. The working space is limited by the needle maximum curvature, which could cause some targets to be unreachable and increase final error. In this experiment only reachable targets were selected.

V. CONCLUSION

This paper presents 3D needle detection and steering using a novel mechanical-based ROI selection and KF. This approach will be tested in a very short term on more complex tissues such as heterogeneous phantoms or biological tissues. Future work also includes target and obstacles automatic tracking in order to compensate for tissue deformation.

However; preliminary results are already very encouraging: the mechanical-based model enables to predict needle change in orientation. In that way, the needle is easily tracked even when the tip performs large bending. Thus, this approach offers a more robust feedback for 3D needle steering.

VI. ACKNOWLEDGMENT

Figure 1 .

 1 Figure 1. Two successive steps of the mechanical-based needle model are visualized. ds represents the needle displacement due to insertion speed.

Figure 2 .

 2 Figure 2. Schematic of the needle steering control.

Figure 5 .

 5 Figure 5. Needle steering experimental setup. The 3D probe stay fixed during insertion and the robot is designed for brachytherapy purpose.

Figure 4 .

 4 Figure 4. ROI selection schematic: the pixels are swept along the predicted curve with a side margin of 2 pixel and a tip margin of 21pixels.

Figure 3 .

 3 Figure 3. Total axial force in relation with insertion depth for a standard brachytherapy needle in four phantoms with different elasticities.

Figure 7 .

 7 Figure 7. 3D needle steering projection on XZ and XY plans: with (b) and without (a) obstacle. Background shows the needle in 3D ultrasound slices.

Figure 6 .

 6 Figure 6. Errors on needle tip pose during three different tracking and insertion speed profile. The 180° bar corresponds to the change in direction of the double bend path.

Table 1 .

 1 Summary of the estimated cutting force F and the deduced stiffness of the virtual springs K used in the mechanical-based model

		Needle		Brachy.	Biopsy 1	Biopsy 2	Nitinol
	Diameter (mm)	1.15	0.723	0.644	0.511
	Bevel Angle (°)	36	30	30	24
		30	F(N) K (N/mm)	0.17 300	0.06 120	0.04 80	0.05 50
	Young's Modulus (kPa)	128 161 182	F(N) K (N/mm) F(N) K (N/mm) F(N) K (N/mm)	2.11 2400 2.38 2400 2.53 2600	0.78 1100 0.81 1000 1.08 1400	0.71 800 0.71 800 0.89 1000	0.30 180 0.34 230 0.44 320
		209	F(N) K (N/mm)	2.74 3000	1.40 1800	1.13 1600	0.49 400

Table 2

 2

	Errors (mean ± standard deviation) for a static needle detection
		No KF	Kinematic	Mechanical
	Shape error	1.13 ± 0.13 mm	1.14 ± 0.05 mm	1.13 ± 0.07 mm
	Tip error	1.07 ± 0.33 mm	0.91 ± 0.12 mm	1.06 ± 0.15 mm
	Direct. error	6.26 ± 0.76 °	3.68 ± 0.52 °	5.79 ± 0.37 °

This work was partly supported by the French ANR within the Investissements d'Avenir program (Labex CAMI) under reference ANR-11-LABX-0004. We also want to thank Ederson Dorileo and Abdulrahman Albakri for their help.

VII. REFERENCES