T De Rességuier

C Roland

G Prudhomme

E Lescoute

D Loison

P Mercier

Influence of edge conditions on material ejection from periodic grooves in laser shock-loaded tin

Keywords: 62, 50, Ef -High 47, 40, Nm -Shock-wave interactions and shock effects 83, 60, Wc -Flow instabilities

In a material subjected to high dynamic compression, the breakout of a shock wave at a rough free surface can lead to the ejection of high velocity debris. Anticipating the ballistic properties of such debris is a key safety issue in many applications involving shock loading, including pyrotechnics and inertial confinement fusion experiments. In this paper, we use laser driven shocks to investigate particle ejection from calibrated grooves of micrometric dimensions and approximately sinusoidal profile in tin samples, with various boundary conditions at the groove edges, including single groove and periodic patterns. Fast transverse shadowgraphy provides ejection velocities after shock breakout. They are found to depend not only on the groove depth and wavelength, as predicted theoretically and already observed in the past, but also, unexpectedly, on the edge conditions, with a jet tip velocity significantly lower in the case of a single groove than behind a periodic pattern.

Introduction

When a shock wave propagating in a material subjected to dynamic compression breaks out at a free surface, it is reflected into a rarefaction wave and the surface undergoes a sharp acceleration. If the free surface presents geometrical defects such as scratches, pits or grooves, high velocity debris can be ejected from these defects, ahead of the main surface. Distinct defects may generate thin jets (which is sometimes referred to as microjetting) while a global roughness can lead to the expansion of a cloud of fine particles (sometimes called material ejection). Because this cloud may disrupt surface diagnostics used in shock physics (velocity interferometry, pyrometry, reflectivity) and because the impact of the ejecta can cause severe damage to nearby equipment in practical, engineering applications, this process has been widely studied both theoretically and experimentally under impact or explosive loading [1][2][3][START_REF] Zellner | Effects of shock-breakout pressure on ejection of micron-scale material from shocked tin surfaces[END_REF][START_REF] Zellner | Probing the underlying physics of ejecta production from shocked Sn samples[END_REF][START_REF] Buttler | Unstable Richtmyer-Meshkov growth of solid and liquid metals in vacuum[END_REF][START_REF] Durand | Large-scale molecular dynamics study of jet breakup and ejecta production from shock-loaded copper with a hybrid method[END_REF][START_REF] Shao | Atomistic simulations of shockinduced microjet from a grooved aluminium surface[END_REF][START_REF] Dimonte | Ejecta source model based on the nonlinear Richtmyer-Meshkov instability[END_REF][START_REF] Monfared | Experimental observations on the links between surface perturbation parameters and shock-induced mass ejection[END_REF][START_REF] Monfared | Ejected particle size measurement using Mie scattering in high explosive driven shockwave experiments[END_REF][START_REF] Durand | Mass-velocity and size-velocity distributions of ejecta cloud from shock-loaded tin surface using atomistic simulations[END_REF][START_REF] Cherne | On shock driven jetting of liquid from non-sinusoidal surfaces into a vacuum[END_REF]. In a recent paper, we used laser driven shock loading to investigate microjetting from triangular, individual grooves of micrometric dimensions in several metals, both below and above shockinduced melting [START_REF] De Rességuier | Microjetting from grooved surfaces in metallic samples subjected to laser driven shocks[END_REF]. In particular, transverse shadowgraphy allowed measuring jet tip velocities over a range of experimental conditions. Then the question arose whether such velocities might be different in the case of periodic grooves, more representative of machined or rolled surfaces of practical interest, i.e. would depend on the boundary conditions at the groove edges. Here, we report further laser shock experiments on tin samples with calibrated grooves of approximately sinusoidal shape and different edge conditions. They show unambiguously that ejection velocity does depend on these edge conditions, with clear differences between single groove and periodic pattern, and even between small and high groove number, unlike what might have been intuitively anticipated.

II. Experimental techniques

Samples were cut from (rolled) tin foils of high purity (99.99%) and 230 µm-thickness, with grain size ranging from about 40 µm to about 200 µm. Their surface was polished to mirror-like finish, then straight grooves of approximately sinusoidal profile were engraved by laser ablation. Groove depth 2h (peak to valley) was either 10 µm or 18 µm, controlled by threedimensional optical microscopy (Fig. 1). For both depths, two types of targets were prepared:

(i) regular periodic pattern of 17 parallel grooves, with an intergroove wavelength λ=58 µm (Fig. 2a), and (ii) three zone-pattern with 3 grooves, 1 groove and 20 grooves, with the same wavelength of 58 µm (Fig. 2b). A dimensionless parameter widely used to characterize such periodic geometry and to model subsequent ejection dynamics is the product kh = 2πh / λ. Here, the kh products are 0.5 and 1 for 10 µm-deep and 18 µm-deep grooves, respectively. The experiments were performed in the LULI2000 facility of the Laboratoire pour l'Utilisation des Lasers Intenses (LULI, UMR 7605, France). A high power laser pulse of 1.055 µm-wavelength, 3.4 ns-duration, about 700 J-energy, was focused on a 3.1 mm-diameter spot opposite to the grooved surface (Fig. 3). The resulting intensity was about 2.8×10 12 W/cm².

The sample was set in secondary vacuum to avoid laser breakdown in air. A thin layer of material is ablated into a plasma cloud, which expansion toward the laser source drives by reaction a short compressive pulse onto the loaded spot. The amplitude Pload, about 86 GPa, and temporal profile of this pulse were inferred at each shot from one-dimensional (1D) computations of laser-matter interaction with the ESTHER hydrocode [START_REF] Colombier | Hydrodynamic simulations of metal ablation by femtosecond laser irradiation[END_REF], using the measured profile of laser intensity as input condition, and assuming spatial uniformity of the laser energy distribution over the irradiated spot. The compression front rapidly steepens into a shock wave followed by gradual release, so that the pressure profile becomes roughly triangular (load

(A) (B) (C) 1 mm 50 µm
usually referred to as unsupported shock) and peak pressure decays with propagation distance (Fig. 4). Shock breakout pressure PSB after propagation throughout the sample thickness was inferred from the same 1D computations, since radial release waves from the periphery of the loaded spot do not affect wave propagation during this transit. 2D effects inside the grooves, over the last 10 or 18 µm before breakout, were neglected. In the shots reported here, PSB = 30.4

GPa ± 1.1 GPa. It is higher than the threshold pressure of about 20 GPa above which tin starts melting upon release [START_REF] Mabire | Shock-induced polymorphic transition and melting of tin[END_REF][START_REF] Buttler | Unstable Richtmyer-Meshkov growth of solid and liquid metals in vacuum[END_REF], so that ejection processes upon shock breakout are expected to occur in a mixed solid-liquid state.

Fast transverse shadowgraphy provides sequences of quasi-instantaneous images of material ejection behind the grooved free surface, using two high speed cameras and a wedge beamsplitter as depicted in Fig. 3. Illumination is supplied by a continuous laser of 527 nmwavelength, divided by a wedge beamsplitter, and four images are recorded by two high speed cameras with different delay times (two images each), which allows to reduce, adapt and control the interframe time at each shot. Short exposure times of 5 ns ensure minimum motion blur (typically 10 µm for a particle travelling at 2 km/s). Images are 1360×1024 pixels with equivalent size of 3.5 µm/pixel. All velocity values mentioned next are derived from initial positions and four successive positions recorded at times 250, 500, 850 and 1100 ns after the laser shot (see Appendix). Error bars in these values, essentially due to uncertainties in determining the positions of the diffuse fronts in the pictures, are estimated to about ± 3 %. All measured velocities are found to be constant within experimental uncertainty, which is consistent with ejection in vacuum, over such relatively long times after shock breakout. with a slightly lower shock breakout pressure of 26.7 GPa, where free surface velocity was 1.85 km/s and jet tip velocities were 2.99 km/s and 3.45 km/s for kh = 0.5 and kh = 1, respectively [START_REF] Buttler | Unstable Richtmyer-Meshkov growth of solid and liquid metals in vacuum[END_REF]. Such consistency is worth noting because both spatial and temporal scales are very different between both experiments. Thus, sample thickness (230 µm) and diameter of the loaded zone (3.1 mm) were about 4 mm and several cm, respectively, in Ref. [START_REF] Buttler | Unstable Richtmyer-Meshkov growth of solid and liquid metals in vacuum[END_REF]. Besides, in our laser shock experiments, the pressure pulse emerging at the grooved surface has a width at half maximum of about 10 ns (Fig. 4), whereas the duration of pressure application in explosive-driven shock loading, governed by the expansion of the detonation products, is typically of µs-order, longer by two orders of magnitude. Finally, this longer duration of Taylor shockwave loading also implies a much slower decaying rate with propagation distance than in the laser-driven shock. the same three zones as in Fig. 2b. Again, ejection velocity consistently increases with perturbation amplitude (i.e. groove depth). For a given shot, assuming homogeneous pressure loading in the central region as inferred from Fig. 5, the only difference between the three zones is the number of grooves, i.e. the boundary conditions at the edges of the zone. Unexpectedly, these edge conditions are shown to strongly affect peak ejection velocity, which is significantly lower for a single defect than for a periodic pattern, and which increases with groove number to approach the values Vj [START_REF] Monfared | Experimental observations on the links between surface perturbation parameters and shock-induced mass ejection[END_REF] and Vj [START_REF] De Rességuier | On the dynamic fragmentation of laser shock-melted tin[END_REF] reported above for the regular patterns.

IV. Discussion

While ejecta production is known to depend on the loading conditions (i.e. shock breakout pressure and temporal profile of the drive) and on the geometry (wavelength and amplitude of surface perturbations), the influence of boundary conditions reported here had never been evidenced yet. In most phenomenological pictures of material ejection, including the widely used Richtmyer-Meshkov Instability (RMI) approach [START_REF] Buttler | Unstable Richtmyer-Meshkov growth of solid and liquid metals in vacuum[END_REF][START_REF] Dimonte | Ejecta source model based on the nonlinear Richtmyer-Meshkov instability[END_REF][START_REF] Monfared | Experimental observations on the links between surface perturbation parameters and shock-induced mass ejection[END_REF], the spike raising from each concave segment (i.e. the bottom of the groove) is fed by matter flowing from the convex edges, forming bubbles on both sides. Meanwhile, rarefaction waves are induced upon breakout of the compression front at every free surface, either curved (groove edge) or planar (Fig. 7). Interactions between these waves cause tension, which will lead to fracture in a solid sample (as discussed and evidenced in Ref. [START_REF] De Rességuier | Microjetting from grooved surfaces in metallic samples subjected to laser driven shocks[END_REF]), or produce cavitation in a liquid metal [START_REF] De Rességuier | Experimental investigation of liquid spall in laser shock-loaded tin[END_REF][START_REF] De Rességuier | On the dynamic fragmentation of laser shock-melted tin[END_REF][START_REF] Loison | Skew Photonic Doppler Velocimetry to investigate the expansion of a cloud of droplets created by micro-spalling of laser shock-melted metal foils[END_REF]. Such cavitation is expected to bound the amount of material to be excavated into each jet, so that this amount would be less in a periodic pattern (Fig. 7a) than in the case of a single groove (Fig. 7b). Assuming a same incoming kinetic energy in both cases, globally transferred to the ejecta, resulting velocities would then be higher in the periodic case (lower jetted mass)

than behind a single groove, which qualitatively agrees with our observations. The apparent thickness of the central (single) jet in the pictures, about 200 µm, tends to confirm that ejecta probably include material flowing from well outside the initial groove (of much smaller width, about 50 µm). Interestingly, radiographs of ejection from explosive-loaded tin, both proton [START_REF] Buttler | Unstable Richtmyer-Meshkov growth of solid and liquid metals in vacuum[END_REF] and x-rays [START_REF] Monfared | Experimental observations on the links between surface perturbation parameters and shock-induced mass ejection[END_REF], often exhibit skew "cuts" growing from the edges of the grooved regions backward into the bulk of the sample, which look very similar to the final picture in Fig. 7b.

Throughout the pulse reflection process, additional contribution of the incident unloading wave produces further cavitation and "micro-spall" [START_REF] De Rességuier | Experimental investigation of liquid spall in laser shock-loaded tin[END_REF][START_REF] De Rességuier | On the dynamic fragmentation of laser shock-melted tin[END_REF][START_REF] Loison | Skew Photonic Doppler Velocimetry to investigate the expansion of a cloud of droplets created by micro-spalling of laser shock-melted metal foils[END_REF], not represented in Fig. 7 for clarity.

Finally, in the periodic case, neighbouring jets might interact with each other at later stages, which might limit transverse expansion and subsequent particle dispersion, and enhance axial velocity. Logically, the three-groove pattern seems to correspond to an intermediate regime between the single groove and the fully periodic case.

V. Conclusion

High power laser shock experiments on tin samples with various types of approximately sinusoidal grooves in their free surface have shown that peak ejection velocities depend not only on the groove depth and aspect ratio (which is consistent with jetting theories and former experiments) but also on the number of grooves (which was not anticipated). This dependence has been attributed to boundary effects at the groove edges, which have been qualitatively

discussed in terms of cavitation caused by wave interactions, ejected mass, kinetic energy, and possible interactions between adjacent jets. Although edge effects on material ejection were evidenced under explosive loading in the past, their influence on jet tip velocity has not been clearly observed in such conditions yet. Thus, further work is still needed for a full characterization and a quantitative interpretation of these edge effects. Finally, those laser shock experiments involving periodic, approximately sinusoidal perturbations and shock-induced melting should be directly comparable with theoretical predictions based on the RMI approach.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Three-dimensional reconstruction from optical microscopy of the sinusoidal grooves

Fig. 3 .

 3 Fig. 3. Schematic of the experimental setup. Grooves in the sample free surface are parallel to

Fig. 4 .Fig. 5

 45 Fig. 4. One-dimensional simulation of laser shock loading for a 2.79 TW/cm² laser intensity

Fig. 5 .

 5 Fig. 5. Shadowgraphs recorded behind the free surface of a 230 µm-thick tin sample subjected

Fig. 6 .

 6 Fig. 6. Shadowgraphs recorded behind laser shock-loaded tin samples with a three zone-pattern

Fig. 7 .

 7 Fig. 7. Schematic picture of jetting upon shock breakout in the case of a periodic pattern (a) or

Acknowledgements

The access to the LULI facility was allocated by the Institut Laser Plasma (ILP, FR # 2707).

We thank all the LULI staff for technical support, as well as Patrick Combis and Laurent Videau (CEA Arpajon) for providing the ESTHER code.

Appendix : analysis of the shadowgraphs

In the raw images recorded behind the sample surface at successive delay times (e.g. Fig. 5 and6), the positions of all fronts (planar surface and jet tips) were determined in pixels, converted to mm, then plotted versus time (Fig. 8), taking into account the transit time of 43 ns, inferred from the simulations, between laser shot (camera trigger) and shock breakout (beginning of surface motion). The position of the planar surface outside the shocked zone, unchanged over the time of observation, was used as a fixed reference for both cameras. The velocities were derived from a linear fit (forced to pass through the origin) of the front positions (Fig. 8).


Daniel L Silverio

Henri A Van Kalkeren

Ta-Chung Ong

Mathieu Baudin

Maxim Yulikov

Laurent Veyre

Pierrick Berruyer

Sachin Chaudhari

David Gajan

David Baudouin

Matthieu Cavaillès

Basile Vuichoud

Aurélien Bornet

Gunnar Jeschke

Geoffrey Bodenhausen

Anne Lesage

Lyndon Emsley

Sami Jannin

Chloé Thieuleux

Christophe Copéret

Tailored Polarizing Hybrid Solids with Nitroxide Radicals Localized in Mesostructured Silica Walls

Keywords: magic angle spinning dynamic nuclear polarization, mesostructured silica, dissolution dynamic nuclear polarization, hybrid organic-inorganic materials

Hyperpolarization by dynamic nuclear polarization relies on the microwave irradiation of paramagnetic radicals dispersed in molecular glasses to enhance the nuclear magnetic resonance (NMR) signals of target molecules. However, magnetic or chemical interactions between the radicals and the target molecules can lead to attenuation of the NMR signal through paramagnetic quenching and/or radical decomposition. Here we describe polarizing materials incorporating nitroxide radicals within the walls of the solids to minimize interactions between the radicals and the solute. These materials can hyperpolarize pure pyruvic acid, a particularly important substrate of clinical interest, while nitroxide radicals cannot be used, even when incorporated in the pores of silica, because of reactions between pyruvic acid and the radicals. The properties of these materials can be engineered by tuning the composition of the wall by introducing organic functionalities. out the d-DNP experiments. D.

Introduction

Dynamic nuclear polarization (DNP) provides a way to dramatically enhance NMR signals, and as a result has opened new frontiers in the characterization of a broad range of molecules (e.g., biomolecules, metabolites, surface species, complex materials, etc.). [1] In a DNP experiment, the Boltzmann polarization of electron spins of unpaired electrons (e.g., stable radicals such as nitroxides or trityls) is partially transferred to the nuclear spins of the target molecules. [2] This polarization transfer leads to NMR signal enhancement () and thereby a reduction of the acquisition time. Proton polarization can propagate throughout the entire frozen sample via spin diffusion, and cross polarization (CP) from protons can be used to obtain sensitivityenhanced NMR signals of heteronuclei (e.g., 13 C, 29 Si,etc.).

DNP is increasingly used to enhance the sensitivity of solid-state NMR spectra at cryogenic temperatures (100 K) under magic angle spinning (MAS-DNP) conditions, [2b] or for hyperpolarizing metabolites at 1.2 K that are rapidly dissolved and transferred to solution NMR spectrometers or MRI scanners (dissolution-DNP). [3] Biologically relevant hyperpolarized analytes or metabolites can be administered to cell cultures or living animals to reveal metabolic activity allowing, [3b, 3d, 4] inter alia, a reliable diagnosis of prostate cancer. [5] However, this kind of DNP experiment suffers from several limitations. Notably, the radicals must be soluble in a solvent that forms a glass at cryogenic temperature to ensure their homogeneous distribution. The radicals should not be in close contact with the analyte, as this can lead to signal quenching or, worse, reactions with the analyte. [2,6] In dissolution DNP (d-DNP) experiments, rapid removal of the radicals from the hyperpolarized solution is desirable both for the sake of biocompatibility and to avoid rapid nuclear relaxation leading to losses of signal enhancement. [7] Some attempts to alleviate these problems have recently been made by sequestering the radicals in micelles [8] or cryptands, [9] hiding the substrate in pores that are too small for the radical to enter, [10] or shielding the radicals in dendritic structures. [11] However, while they can be efficient when applicable, none of these solutions are general, and all have drawbacks.

Mesostructured silica materials [12] have emerged as an alternative strategy to alleviate these limitations. For instance, the use of mesoporous silica materials containing radicals that are regularly distributed at their pore surface obviates the use of a glassing agent to polarize an analyte in pure water because the porous structure prevents solvent crystallization and aggregation of radicals (Figure 1). [13] Another example is the use of mesoporous silica in MAS DNP to avoid segregation between radicals and colloids upon cooling. [14] However, in both systems the radical remains in close contact with the analyte and in the second case the approach requires the use of a solvent in which the radical is soluble. Recent work has also shown that radicals could be introduced within the silica matrix via a sol-gel process using a bis-silylated organic moiety that can generate radicals, thus providing a polarizing matrix for MAS DNP, though the localization of the radicals in the silica matrix and their shielding from the analyte were not investigated. [15] We reasoned that hybrid mesostructured materials containing radicals embedded in the silica walls (Figure 1 bottom and Figure 2a) [12,16] could provide a powerful approach to isolate the radical from the analyte, avoiding possible reactions and signal bleaching, while still removing the restriction to specific solvents or diluting the analyte in a co-solvent mixture. We report here the development of a HYPSO material containing radicals in the walls of the material (HYPSO-W) through a template approach and demonstrate its efficiency for MAS-DNP and d-DNP experiments. In particular, we show how engineering the wall structure and its composition by introducing organic functionalities significantly improve the polarization properties.

First, using 1, a tris-silylated precursor containing a nitroxide radical (see ESI), and various amount (X) of Si(OEt)4, the corresponding HYPSO material (HYPSO-W1 [1/X]) was prepared by a sol-gel process using Pluronic 123® (P123) and NaF as structure-directing agent and catalyst, respectively. [17] Optimal removal of P123 was achieved by washing with pyridine•HCl. [18] Using these protocols, a broad range of radical concentrations from 0.022 to 0.11 mmol/g (as determined by EPR; 49 to 58% yield) could be introduced into the silica matrix. We also synthesized materials HYPSO-W2 that contain bis(triethoxysilyl)ethane 2 (Figure 2a, bottom), as an organic dopant to help polarization transfer via 1 H-1 H spin diffusion from the radicals to the pores. [21] HYPSO-W2 [Y/1] materials were all synthesized with the same 1/200 ratio of 1 to Si, but with varying content of the organic dopant 2 as defined by

[Y/1] = [Si(OEt)4/2].
All materials HYPSO-W1 and HYPSO-W2

have large pore volumes (0.7 to 1.2 cm 3 /g), high surface areas (673 to 1066 m 2 /g), narrow pore size distributions (5 to 9 nm) and large wall thicknesses (3 to 6 nm). These properties allow them to enclose the radical 1 (molecular diameter of ca. Both trends agree with an increasing average radical-radical distance with decreasing concentration. [13a, 22] For HYPSO-W2, <T1e> and line width are essentially unaffected by varying the concentration of bissiloxyethane 2 (see ESI).

The 1 H enhancement () of the solvent peak in MAS-DNP for materials impregnated with a 4:1 mixture of D2O:H2O did not change significantly with decreasing radical concentration from 0.11 to 0.045 mmol/g, staying at ca. 50-55 (Figure 2b), in contrast to 25 for HYPSO-1 under the same conditions. Going to an even lower concentration, the enhancement drops but still remains at reasonable levels, i.e., = 34 for HYPSO-W1 [1/300] and 22 for HYPSO-W1 [1/400]. Since the enhancements observed with HYPSO-W1 [1/75] to [1/200] are approximately the same, we chose HYPSO-W1 [1/200] as the optimal material, since lower concentrations of radicals should lead to a lower quenching factors. [23] Compared with HYPSO-W1 [1/200] (Figure 2c), all HYPSO-W2 materials afforded approximately the same levels of enhancement, with HYPSO-W2 [25/1] giving

N O O O N H Si(OEt)3 3 R R = R O Si O Si O Si O Si O Si O Si O Si O Si O Si O Si O Si O Si O Si O Si
O Si O Si O Si O Si Si Si O Si O Si Si Si O Si O Si O Si O Si Radical HYPSO-W1 HYPSO-W2 (EtO)3Si 2 = Si(OEt)3
somewhat lower enhancements due to its lower concentration in radicals (See ESI for details).

Notably HYPSO-W1 [1/200] yields effective MAS-DNP enhancements in solvents in which traditional polarizing agents like TEKPol [24] and AMUPol [25] are poorly soluble, as well as solvents that are poor at forming glasses (Figure 3), properties that, as mentioned above, are generally detrimental for DNP. Examination of the enhancement afforded by HYPSO-W1 [1/200] for a mixture of D2O:H2O at different magnetic fields reveals an enhancement = 100 at 9.4 T (400 MHz proton frequency). An enhancement = 20 is obtained at 18.8 T (800 MHz proton frequency, entries 7 and 8), an expected decrease with increasing field strength for DNP processes dominated by the cross effect. [2a, 26] Overall, this material is nearly twice as effective as state-of-the-art radical-containing materials for polarization [13,15] in particular at high magnetic fields. [26] To probe how well the radicals in the HYPSO-W materials are shielded from the analytes, we first treated HYPSO-W1 [1/200] with a solution of ascorbic acid (known to reduce even sterically hindered nitroxide radicals), [27] and determined the resulting radical loading by EPR (Figure S1a). The ascorbic acid treatment resulted in only a 36% decrease of radical content.

For comparison, this treatment was applied to HYPSO-1, a material with comparable physical features but with the radical localized at the pore surface, resulting in complete destruction of the radicals according to EPR (Figure S1b). However, both ascorbic acid-treated materials gave essentially no enhancement in the DNP experiment (Figure S1a). The poor enhancement in HYPSO-W1 treated with ascorbic acid indicates that the remaining radicals are probably too deeply embedded to contribute to the polarization of the molecules present in the pores.

We therefore examined the polarization properties of HYPSO-W1 and the organic doped HYPSO-W2 materials (Figure 2d opened circles). First, the decrease in radical content upon reaction with ascorbic acid is lower for HYPSO-W2 than for HYPSO-W1 [1/200] (decrease of 22% vs. 36%). Second, and more importantly, these treated HYPSO-W2 [100/1] materials still afford a 1 H enhancement of 27 at 600 MHz, a mere 36% decrease with respect to its untreated form. When increasing the ratio of 2/Si(OEt)4 to 1/25, no loss of enhancement is observed after treatment. Overall the introduction of organic functionalities in the wall clearly helps to protect the radicals and to transfer the polarization from the radicals to the solvent molecules in the pores.

With these promising results, the HYPSO-W materials were also tested in dissolution DNP using pyruvic acid as a prototypical metabolite, [3b, 3d, 4, 28] which is known to reduce nitroxide radicals, so that previous nitroxide-containing materials failed to give rise to any polarization in d-DNP.

Dissolution-DNP experiments were carried out in a polarizer [29] operating at 1.2 K and 6.7 T equipped for cross polarization at low temperatures. [30] The optimal material is HYPSO-W2 [25/1], which affords a polarization P( 1 H) = 39.5% for neat pyruvic acid with DNP( 1 H) = 20 minutes, resulting in P( 13 C) ~17% in about 60 minutes after ascorbic acid treatment. By comparison, a trityl-containing material, compatible with neat pyruvic acid, requires a much longer polarization time of 105 minutes. [13b] The corresponding HYPSO-W1 [1/200], which lacks the ethyl fragments in its walls, affords a polarization of only P( 1 H) = 22.4% with DNP( 1 H) = 13 minutes. While ethyl groups are not mandatory to polarize the frozen solutions in the pores, their presence is very beneficial as it significantly boosts the polarization rates and improves the overall maximum polarization attained.

In conclusion, the introduction of electron polarization sources in the walls of mesostructured materials, dubbed HYPSO-W yield the highest enhancements in MAS-DNP for HYPSO materials, and are amenable to performing DNP in a wide variety of solvents. Importantly, HYPSO-W materials enable polarization of substrates in solution that would otherwise react with the nitroxides in conventional DNP formulations, e.g. neat pyruvic acid, combining the fast polarization buildup of nitroxide radicals and the stability of embedded nitroxide radicals. One critical aspect is the engineering of the structure of the material via incorporation of organic functionalities in the wall allowing more effective polarization transfer to the analyte.

  2 nm) and to polarize a large volume of analyte. The SA-XRD diffraction patterns complemented by Transmission Electron Microscopy (TEM) are consistent with an ordered hexagonal structure of the porous network for all materials except HYPSO-W2 [2/1], where the pattern is more indicative of a "wormlike" structure. The radical distribution was evaluated by EPR: decreasing the radical concentration leads to a narrowing of the EPR linewidth (13.1 G for HYPSO-W1 [1/75] to 12.3 G for HYPSO-W1 [1/400]) and an increase in the mean electron inversion recovery time (<T1e>) from 32 µs for HYPSO-W1 [1/75] to 210 µs for HYPSO-W1 [1/400].

Figure 1 .

 1 Figure 1. Comparison of materials with radicals in the pores, HYPSO-1, with materials with radicals in the walls, HYPSO-W.

Figure 2 .

 2 Figure 2. a) Synthesis of HYPSO-W1 and HYPSO-W2 materials containing radicals in the walls. MAS DNP performance (Proton enhancement (blue circles) and T1 (orange circle) of b) HYPSO-W1, c) HYPSO-W2 and d) HYPSO-W2 and HYPSO-W1 after treatment with ascorbic acid. All enhancements are determined at 600 MHz proton frequency (14.1 T) and ~100 K with 10 kHz MAS using 4:1 D2O:H2O as a solvent. T1 determined by a saturation recovery experiment. See ESI for details.

Figure 3 .

 3 Figure 3. MAS-DNP enhancement () in various solvents and at various magnetic field strengths with HYPSO-W1 [1/200].Enhancements determined at given magnetic field at ~100 K. T1 determined by saturation recovery. At 9.4 T, D2O ; H2O = 4:1 gave a lower 1 H enhancement of 71 in this case, although there was no difference in enhancement between these two solvent ratios at 14.1 T.

Acknowledgements

We acknowledge the SNF for the 600 MHz DNP spectrometer (206021_150710). This work was supported by the Swiss National Science Foundation (SNF), ETH Zürich, the Ecole Polytechnique Fédérale de Lausanne, Bruker BioSpin Switzerland AG, ERC Advanced Grant No. 320860 and EQUIPEX contract ANR-10-EQPX-47-01,and the ERC Advanced Grant "Dilute para-water". Mr. Wei-Chih Liao, Mr. Tigran Margossian, Ms. Laura Piveteau, and Dr. Martin Schwarzwälder are acknowledged for many useful discussions.


Phuoc Vy

Vincent Auffret

Miguel Castro

Pierre Badel

Michel Rochette

Pascal Haigron

Stéphane Avril

Patient-specific simulation of guidewire deformation during transcatheter aortic valve implantation

Introduction

Transcatheter Aortic Valve Implantation (TAVI) is a mini-invasive procedure to replace aortic valves in the context of Aortic Stenosis. This procedure consists in delivering a fully collapsible replacement valve (prosthesis) onto the native aortic valve through a catheter (sheath). A stiff guidewire is first inserted in the patient's aorta with access via the femoral artery. The function of the guidewire is that of a rail along which the sheath can glide. This function is essential to ensure the navigation along the aortic arch up to the heart. The first TAVI procedure was achieved in 2002 by Alain Cribier [1]. Since then, TAVI has permitted the treatment of patients who cannot benefit of an open heart surgery [2][3]. Nowadays, TAVI shows excellent 1-year outcomes for both high-risk and intermediate-risk patients [4]. Building up on these excellent results, TAVI procedures might drastically increase in the future. Despite this situation, a number of aspects remain to be improved, including the long term durability and the management of possible complications. Research efforts by medical teams are focused on identifying the mechanisms of complications, and determining the strategies to handle them [5][6]. It was shown that the risk of complications may be increased in case of inappropriate positioning of the delivery system in the valvular plane. For instance, paravalvular leaks [7][8][9], atrioventricular block [10] and coronary occlusion [11] could be related to the invalid implantation position of the prosthesis within the native aortic valve.

Patient-specific numerical simulations can help clinicians during pre-operative planning and intraoperative navigation, as shown by different reviews recently published [12][13][14]. Several studies have already proved the feasibility of simulating prosthesis deployment within patient-specific anatomy. (Bosmans et al., 2016) validated a self-expanding prosthesis deployment model by comparing it with post-operative data. They obtained an accurate prediction of the final shape of the prosthesis [15]. (Wang et al., 2014) simulated balloon-expandable valve deployment within several patient-specific aortic valve geometries [16]. They illustrated how the interaction between calcifications and Valsalva sinuses could possibly induce aortic rupture. Simulation was instrumental in the decision of the medical team for the simulated patient cases. (Auricchio et al., 2014) suggested the measurement of the coaptation area of prosthesis leaflets to estimate the risk of transvalvular leaks [17]. (Morganti et al., 2014) simulated the deformed geometry of deployed prosthesis within the aortic valve [18]. They were able to predict risks of perivalvular leaks by studying the contacts between the stent and the aortic valve leaflet.

Recently, (Morganti et al., 2015) confirmed that the implantation depth and release angle had a crucial role on determining valve anchoring, device deformation and risk of regurgitation [19]. This highlights the importance of correctly positioning the delivery system before and during deployment. However, the simulation of delivery tools (prosthesis, sheath and guidewire) before deployment has not been tackled yet despite their importance on the outcome of the procedure. Such simulation raises many difficult challenges, so we focus on a single tool as a preliminary step towards this goal. To our best knowledge, the stiff guidewire used in TAVI was never simulated, despite the valuable information that can be provided on the mechanical behavior related to tools alignment. To address this issue, we introduce a novel simulation approach of the stiff guidewire in this paper. This novel simulation approach aims at predicting guidewire deformation and positioning in the first stage of TAVI procedures. The paper is organized as follows. In section 2, the mechanical problem is first formulated from a physical analysis of the stiff guidewire insertion. Then, a numerical model solving this problem is defined. Finally, the methods for our sensitivity analysis and validation are described. The simulation results are reported in section 3 and are discussed in section 4.

Materials and methods

The purpose of the model is to predict the shape of a stiff guidewire inserted into a patient-specific aortic and ventricular geometry. The mechanical problem consists in determining the stresses and strains of a beam (the guidewire) constrained to fit within a cavity (aorta and left ventricle). The problem is solved with the Finite Element (FE) method using an implicit scheme.

Mechanical problem formulation

Assumption on motions

Respiratory and cardiac movements may distort the complete structure [20] [21] [22]. However, the proposed approach was to neglect the dynamic effects of this motion on the guidewire. In addition, guidewire motion was assumed quasi-static because its low mass produced small inertial force compared to the contact forces. Only the final position of the guidewire was predicted by the final static equilibrium of the simulation.

Mechanical properties

According to intra-operative observations, stiff guidewire insertion led to marginal deformations of the aorta. Thus, the aortic wall and the ventricle were assumed to remain rigid. Under this assumption, it was not necessary to model the influence of blood pressure, pre-stress within the vascular muscles, and the non-linear behavior of the aortic wall over large deformations. Despite this, a model with linear elastic tissues was tested for the sake of feasibility proof. Amplatz Super Stiff and Amplatz Extra Stiff guidewires, commonly recommended for the Corevalve TM and SAPIEN TM prostheses respectively, were modeled in this study. The elastic properties of the Amplatz Super stiff guidewire used in CoreValve TM implantation were characterized in prior studies [23][24]. Fig. 1.b shows the bending stiffness along the guidewire. We did not observe any permanent deformation on the stiff guidewire caused by the intervention. Excluding gothic aortic arches, we assume that the usual curvature of the aortic arch is not sufficient to cause plastic strains.

The guidewires were designed with a straight end and may perforate the ventricular wall. To minimize the perforation risk, the guidewire has to be curved at its end before the intervention. To that purpose, clinicians manually deform the guidewire to create a curved end, as shown in Fig. 1.a. The obtained deformation is permanent (plastic deformation). This pre-deformation of the guidewire tip was unknown for both patients considered in this study. Therefore, we investigated the impact of this pre-deformation through a sensitivity analysis where we varied both the angle and the radius of curvature. The straight floppy part at the end of the guidewire (Fig. 1

Contacts and boundary conditions

Interactions between the guidewire and the aortic wall result in contact forces. Friction (tangent) contact forces were neglected because the guidewire had a PTFE coating to reduce friction.

In the model, a slipping rail constrained the direction and the orientation of the proximal end (base) of the guidewire in the thoracic aorta. Translations and rotations applied by the operator were modeled as Dirichlet boundary conditions applied at the extreme distal node of the guidewire.

During actual TAVI procedures, the proximal end of the guidewire is free in the left ventricle while the distal end is moved by the operator. Usually, during its introduction through the femoral dilator, the curved end of the guidewire can take different orientations, as shown in Fig. 1.c and d. Nevertheless, it remains possible to adjust the orientation when the guidewire reaches the ventricular cavity. We defined the model in order to take into account the possibility of different orientations in the ventricular cavity.

Summary of physical parameters

The physical parameters input in the model are summarized in Tab.1. They include: the aortic geometry, the guidewire stiffness, the angle and radius of the curved tip, the axial force applied to the distal end of the guidewire and the initial orientation. A sensitivity analysis was performed on these physical parameters to determine how they affect the model predictions. This analysis is designed to provide insights on the behavior of the guidewire rather than provide exact recommendation on the parameters. Indeed, those parameters are difficult to control precisely. The geometry of the stiff guidewire is approximately created as the clinician manually deforms it without any precise tool, or selected in a restricted catalog of pre-shaped guidewire. Likewise, the orientation of the stiff guidewire can be adjusted during its insertion but it is difficult to measure.

Discretization & resolution method

Patient specific geometry

All patients who undergo TAVI at Rennes University Hospital provide informed consent to participate in the French National Transcatheter Aortic Valve Registry, FRANCE-TAVI [25], allowing anonymous collection and processing of their clinical data. These interventions are commonly performed under fluoroscopic control. The 2D visualization of the aortic root can only be achieved via the injection of a contrast agent. The images used in the present study were retrospectively collected after the procedures of 2 patients who underwent TAVI at Rennes University Hospital. For both, a CoreValve prosthesis (Medtronic TM ) was implanted from retrograde transfemoral access. These retrospective cases were chosen as both had intra-operative images of sufficient quality to permit their registration to the 3D pre-operative CT images and to be used for the comparison with numerical simulations. The intra-operative images consisted of digitally subtracted angiography (DSA) acquired during the protocol when the guidewire was inserted in the ventricle of the patient. The ascending thoracic aorta appeared clearly enough in these images thanks to the contrast injection (Fig. 2

.a).

A 3D pre-operative CT image was routinely acquired for clinical assessment before TAVI interventions. The unsynchronized CT angioscan acquisition was performed using Discovery CT750HD (General Electric TM ) (Fig. 2.b). A 3D reconstructed image of cardiac structures is shown in Fig. 2.c. In addition, the clinical protocol included a synchronized CT-scan showing exclusively the aortic valve. In our study, it was assumed that the ventricle geometry obtained in the unsynchronized CT image was fully dilated and close to diastolic state. This was verified by superimposing the synchronized and unsynchronized CT images.

3D reconstruction of the left ventricle and aorta was interactively performed using 3D-slicer [26]. The aortic section was smoothed but did not behave well at the supra-aortic trunks (Fig. 2.d). Moreover, the automatic segmentation created an irregular rough surface that is unlike the real aorta or ventricle. The surface was smoothed and imported in ANSYS using Spaceclaim, the integrated CAD module of Ansys. Only the aorta was segmented, leaving aside collateral branches and the supra-aortic trunk. The aorta was segmented from the aortic root to the distal thoracic aorta. The artifacts from the supraaortic trunks were interactively smoothed and a more realistic mesh was recovered (Fig. 2.c). This correction must be carefully performed to limit the impact on the results. We assumed that the interactions between the guidewire and the abdominal aorta, the iliac and the femoral artery had a marginal effect on the position of the guidewire at the aortic valve because of the long sheaths specifically used in TAVI procedures. This paper mainly focuses on the investigation of patient #1. The investigation of patient #2, reported in the Appendix, was aimed at the verification of the applicability of the model on different geometries.

Meshing

The geometry was meshed in ANSYS. The rigid aorta was meshed with quadrilateral and triangle contact elements (TARGE170) as shown in Fig. 1c. The default mesh had 9610 nodes for patient #1. The guidewire was meshed with quadratic 3-node beam elements based on Timoshenko beam theory (BEAM189) and line contact elements (CONTA175).

Resolution scheme

Neglecting friction and plasticity, the problem was conservative. However, due to contacts and large deflections of the guidewire, it turned to be highly nonlinear. The simulations were performed on the commercial software ANSYS 17.1 with full Newton-Raphson method and sparse matrix direct solver. Time substeps ensuring convergence of the resolution were defined automatically by ANSYS. A bisectionsubdivision of a substepwas performed when more than 26 equilibrium iterations were required to converge a substep. The minimum allowed time step was 5.10 -3 s for every simulation.

Boundary conditions and stabilization techniques

The resolution was achieved in 3 stages, as illustrated in Fig. 3. Stage I consisted in deforming the guidewire to match the shape of the centerline of the aorta. Stage II consisted in releasing the kinematic constraints applied at stage I, leading to contact between the guidewire and the aortic wall. In stage III, the end of the guidewire could be moved (axial force or rotation) to simulate the action of an operator.

At stage II, the sudden release of the kinematic constraints is often too abrupt for the Newton Raphson algorithm to face contact nonlinearities. Hence, stabilization techniques were introduced. To this end, springs were connected to the nodes of the guidewire. Longitudinal and torsional springs were compared, with torsional springs turning out to be more successful. For each stabilized node, three torsional springs prevented rotations about X, Y, and Z axis. They were defined with COMBIN14 elements. They were attached to guidewire nodes at one end and to a virtual anchoring node (created for that purpose) at the other end. Anchoring nodes had their displacement coupled with that of the respective guidewire nodes. Their rotations, however, were tied to the guidewire rotations with a delay of one loadstep. Torsional springs had high initial stiffness (11.36 N.mm.rad -1 ) so that kinematic constraints could be removed without inducing sudden large displacement. The spring stiffness was decreased at each loadstep by a decay ratio. At the end of the stage II, all torsional springs had a null stiffness.

Two options were considered for the reaction force at the distal end of the guidewire. The first option consisted in adjusting the position of the guidewire during the whole simulation to maintain the desired reaction force. This method will be referred as natural insertion. Conversely, the forced insertion consisted in maintaining the initial depth of insertion throughout stage II.

The distal end of the guidewire was arbitrarily assigned to remain tangential to the centerline of the aorta. It happened that the actual position of guidewires observed intra-operatively at the bottom of the descending thoracic aorta may noticeably deviate from the centerline and be in contact with the wall. However, this was sufficiently far from the aortic root to induce a marginal impact on the results of this study.

Contacts were activated during the stage II. Frictionless contact between patient geometry and stiff guidewire was enforced using an Augmented Lagrangian algorithm. The main parameters determining the behavior of the contact were the contact stiffness coefficient (FKN) and the penetration tolerance (FTOLN). Contact stiffness influenced the contact force applied on the bodies to prevent interpenetration. Nodes with higher penetration than the tolerance do not satisfy the contact compatibility which is also a convergence criterion. This means that the solver has to repeat the equilibrium iteration with a higher contact force. Numerous preliminary simulations showed that direct contact between the node at the tip of the guidewire and the ventricle led to divergence of the Newton-Raphson algorithm. This issue was addressed by defining a small rigid sphere at the tip of the guidewire (c.f. Fig. 1.a). In general, irregular contact zones such as kinks may lead to large variation of the stiffness matrix of the model which would lead to the failure of the resolution algorithm.

Summary of numerical parameters

Numerical parameters are summarized in Tab.1. They include: the mesh density, the number of loadsteps, the initial number of substeps, the predictor (parameter to speed up the process of solving a substep by extrapolating better initial equilibrium iteration), the depth of insertion, the decay ratio of stabilization springs, the contact stiffness and the penetration tolerance.

Sensitivity analysis

Evaluation criteria

The success of model predictions was assessed based on the following characteristics: robustness, computation time and accuracy.

Robustness

The robustness denoted the ability of the simulation to converge in a mathematically satisfying equilibrium solution regardless of the variation of the input parameters. Convergence or divergence of the solution on a specific range of parameters was a straightforward indicator of robustness. It was observed in preliminary simulations that among all physical parameters, the axial insertion force was the main restriction for robustness (as it induced compression of the guidewire that reached its critical buckling load, which is an instability that cannot be handled with implicit resolution). Thus, maximum converging insertion force was deemed to be a relevant indicator of robustness when varying other parameters and results were presented accordingly.

Computation time

The simulations on ANSYS 17.1 were performed on a HP Z800 Workstation having two processors Intel® Xeon® E5620 and 12 GB of RAM. The computation time is reported as the duration of stages I and II of the simulation.

Accuracy

We defined the accuracy of the prediction as the degree of similarity between the simulation result and intra-operative data. The primary output of the simulation was the position and angle of the simulated guidewire in the aortic valve plane. The intra-operative image was a 2D fluoroscopic image, so this comparison was performed by projecting the simulation results according to the current pose of the C-arm. An example of measurement is shown in Fig. 4.b, where the projection of the simulation is overlapped with the intra-operative view. The aortic valve plane was manually defined according to the contour of the contrast injection. The distance between the simulation and the data was measured at the aortic valve plane in pixels (pixel size for the intra-operative data used here: 0.308, 0.308 mm at the C-arm detector). The difference of angle was measured in degree.

A 3D/2D registration was necessary to bring the simulation results into the intra-operative coordinate system, as shown in Fig. 4.b. It is composed of a transformation representing the C-arm pose, a transformation positioning the preoperative CT-scan in the 3D space and a transformation corresponding to the projection performed during the generation of intraoperative images. Numerous 3D/2D registration methods can be found in the literature to obtain such transformation [27]. In this study, a rigid transform combining 3 translations and 3 rotations was iteratively defined. The C-arm pose and the intrinsic perspective projection parameters were extracted from the DICOM header of the intra-operative data in order to initialize the registration algorithm.

The aortic root of the 2D intra-operative image (Fig. 2.a) was manually contoured by an expert. The criterion for registration was the Euclidean 2D distance between the contour of the projected 3D preoperative aorta and the 2D aorta contour observed in the fluoroscopy. A distance map was computed from the 2D aortic contour with the Chamfer operator. The minimization of the distance criterion was iteratively performed using a Powell optimization algorithm until the estimate of the 3D/2D transform was satisfactory [28]. The resulting projection is shown in Fig. 4.a. Two intra-operative angiographies at different pose of the C-arm were available and considered in the registration process for patient 2.

Variations of input parameters

The simulations reported in Tab.1 are sorted according to the variation of input parameters. Sets #1-6 are related to the exploration of physical parameters. Sets #7-12 are related to numerical parameters. The geometry of patient #1 was thoroughly studied under the 12 simulation sets.

Simulation set #1 explored the curvature angle of the guidewire tip from 20° to 120° with a 20° step at a fixed radius of curvature of 30 mm. Sets #2 and 3 explored the radius of curvature from 15 to 30 mm with a 5 mm step at a fixed angle of 20° and 80° respectively. Set #4 explored the stiffness of the guidewire, which varied from 20 to 100 GPa. Set #5 explored the orientation from 0 to 180° with a fixed guidewire shape. Set #6 tested the forced insertion method to increase the initial depth of insertion. The guidewire was initially 47 mm deeper when using the forced insertion method. An insertion force of 0.1-0.2 N as indicated in Tab.1 means that the stiff guidewire was inserted to reach a reaction force of 0.20 N and then retracted. This allowed comparing the shape of the guidewire at 0.20 N and 0.10 N. Set #4 and #6 simply had the guidewire inserted as much as possible until divergence. Finally, the model with linear elastic tissues was tested for feasibility. The tissues had an elastic modulus of 2 MPa and a Poisson's ratio of 0.49. Arbitrary boundary conditions were defined to maintain the anatomy. Elastic supports (0.15 N.mm -3 ) were defined on the thoracic aorta and the apex of the ventricle. The guidewire had a curvature of 20° and a radius of 25 mm. It was inserted with a natural insertion method to reach 0.05 N and 0.1 N. The registration of the simulated deformed aorta was not performed, and accuracy evaluation was not available. Set #7 explored the contact stiffness coefficient (FKN) from 0.1 to 10. Set #8 explored the penetration tolerance (FTOLN) from 0.5 mm to 4 mm. Set #9 explored the mesh density of the guidewire (1.4 node/mm, 2.5 node/mm and 5 node/mm). The stiffness of the torsional stabilization springs was adjusted to provide the same initial level of support. Set #10 explored the decay ratio from 0.75 to 0.9. In set #11, the decay ratio was adjusted with the number of loadsteps so that the stabilization springs reached the same stiffness at their complete deactivation across all simulations of the set (10 -4 N.mm.rad -1 ). Thus, the number of loadsteps was explored while keeping the release rate of the stabilization springs consistent. Finally, set #12 explored the relationship between the number of substep and the predictor.

Results

From a qualitative point of view, the model was able to converge on a large range of parameters. Tab. 2 reports the computation times, including minimum, maximum and mean values, along with the diverged simulations. Tab. 3 reports the minimum, maximum, standard deviation and mean values of accuracy measurements, i.e., angle difference and distance rounded at the nearest pixel.

.

Sensitivity to physical parameters

Influence of the stiff guidewire properties

The sets investigating the effect of the curvature angle and radius (#1, #2, #3) converged. The maximum insertion force for convergence of set #4 is reported in Tab.4b. A linear regression analysis between the maximum insertion force and the stiffness of the guidewire was performed. The regression coefficient is 0.004 N/GPa and the coefficient of determination is 0.9984. Computation times were similar across those simulations (between 3000 s and 4000 s). The simulated guidewires (#1, #2, #3) projected on intra-operative image are shown in Fig. 5 (b, c andd). The standard deviation of the angle difference and distance were low, less than 1° of angle difference and 3px of distance. The curvature radius of the guidewire produced the highest standard deviation. The curvature angle produced the lowest standard deviation. The projections of set #4 were not reported in the figures because the simulated guidewire shapes were too similar.

Influence of the insertion conditions

Set #5 investigating the orientation converged. It required a computation time ranging from 3081 s to 3982 s. Fig. 5.e shows the guidewires of set #5. The deformation of the tip is noticeably influenced by the orientation of the guidewire. The standard deviation was 2° and 4.4 px for angle difference and distance respectively. Set #6 converged when the insertion depth was less than -18 mm or included in a range between -3.9 mm and 10.8 mm. Further range was not tested. The domain of convergence can be observed in the force versus displacement curve in Fig. 6.d. The simulated guidewires at the edge of the domain of convergence are projected in Fig. 6 (a b and c). The variation of insertion depth resulted in a large standard deviation (6.8° and 12.7 px). However, simulated guidewires from forced insertion obtained higher differences with the intra-operative guidewire. The differences produced when varying the insertion force from 0.10 N to 0.20 N (natural insertion) for set #1 are reported in Tab.4.a.

The simulated guidewires are projected in Fig. 5 (a andb). The mean values for the angle difference increased by 3.3° and distance by 1 px with the increasing force. The standard deviations remained similar for both 0.10 N and 0.20 N.

Sensitivity to numerical parameters

We observe from Tab.3 that sets #7 to #12 produce very small standard deviation (less than 0.5° and 0 px). We deduce that the variation of the numerical parameters has negligible influence on the accuracy. Instead, they show a strong influence on convergence and computation time. Excessively low contact stiffness resulted in divergence. Using predictor with a large number of substep slowed the simulation and occasionally led to divergence. The decay ratio should be adjusted so that the torsional springs are deactivated with stiffness around 10 -4 N.mm.rad -1 . The mesh density had a negligible impact on computation time. Each set is detailed in the following paragraphs.

Contact conditions

In set #7, contact stiffness coefficient of 0.1 diverged. All simulations of set #8 converged. The physical parameters were angle = 100°, radius = 30 mm, force = 0.22 N with predictor activated. The standard deviations from both contact stiffness coefficient and penetration tolerance variation were negligible, less than 0.5° for angle difference and 0 px for distance.

Model and solver parameter

Concerning set #9, the coarser mesh (1.4 node/mm) had the shortest total computation time (4056 s). The finer mesh (5 node/mm) had the longest computation time (14928 s). Maximum insertion force of the finer mesh was 0.01 N above that of the coarser mesh. A large fraction of the computation time was spent to compute the last increments of insertion force. The standard deviation was negligible for both angle difference (0.3°) and distance (0 px). The physical parameters for set #10 were angle = 20°, radius = 30 mm, force = 0.20 N with a mesh density of 0.14 node/mm and 10 substeps per loadstep. The fastest simulation ran with a decay ratio of 0.87. During the deactivation of the stabilization springs at loadstep 124, their stiffness was near 10 -4 N.mm.rad -1 . Lower decay ratio generated instability at each loadstep which resulted in an increased computation time. Higher decay ratio slowed down the removal of stabilization springs. If the spring stiffness was still too high during the loadstep of complete deactivation of stabilization springs, instabilities arose and increased the computation time. The projections of set #9 appeared undistinguishable one from another. Careful inspection showed that very few pixels were positioned differently. The standard deviation was negligible for both angle difference (0.3°) and distance (0 px). The physical parameters for set #11 were angle = 20°, radius = 25 mm, force = 0.15 N. The best computation time was 1821 s (decay ratio = 0.7929, 75 loadsteps). The projections of the simulations appeared undistinguishable except small differences in the ventricle. The standard deviation was negligible for both angle difference (0.2°) and distance (0 px). Finally, the physical parameters for set #11 were angle = 20°, radius = 30 mm, force = 0.20 N. A simulation using predictor and 10 substeps per loadstep diverged. The use of predictor and the increase of the substep number did not seem to impact the accuracy of the prediction because the standard deviations were low (0.2° for angle difference and 0 px for distance). Predictor was beneficial if the number of substeps was set to 1. Conversely, predictor severely hindered the simulation when the number of substeps was set to 10.

Discussion

Summary

This study presented a Finite Element modeling framework for predicting the shape of a stiff guidewire inserted into a patient-specific geometry of the left ventricle and aortic arch. The model was solved with an implicit resolution scheme in ANSYS. Then, the performance of the model was evaluated in terms of accuracy, robustness and speed of resolution. Nonlinearity associated with contacts proved to be a challenging problem for the implicit resolution method. Therefore, stabilization techniques were proposed to enhance the robustness of the model and overcome instabilities. The model was able to explore a wide range of parameters and provide insights on the behavior of the inserted guidewire. The results of the simulations were compared to intra-operative 2D images. Very good agreement could be obtained with respect to the intra-operative guidewire position after a thorough sensitivity analysis of the simulation problem. The guidewire properties which produced the most accurate predictions were the following: curvature of 100° and curvature radius of 30 mm, insertion force of 0.20 N with a 90° orientation.

The main conclusion of the sensitivity analysis was that the patient anatomy itself mostly determined the shape of the guidewire and had a more significant impact than the guidewire curvature properties. This can be explained by the high compliance of the guidewire in its curved part compared to the straight body of the guidewire.

Excluding the patient's geometry, the influence of the physical parameters on the prediction was intricate, and prevented from extracting general trends. This was especially the case for the shape of the guidewire in the left ventricle. Therefore, numerical simulations appeared to be useful to understand the insertion and navigation of endovascular tools in the aortic arch and left ventricle.

Contribution and comparison to bibliography

TAVI interventions heavily rely on the experience of medical teams. Thus, numerical models can be a relevant objective complementary support. Most of the published work was focused on the numerical deployment of the prosthesis regardless of its pre-deployment positioning. To the best of our knowledge, no study was performed on the position of delivery tools and prosthesis at the aortic valve. Yet, it is important to deal with the position of the prosthesis because it is a concern for the success of the deployment. Similar studies were performed by (Dumenil et al., 2012) (Gindre et al., 2015) for the insertion of guidewires during endovascular aneurysm repair [28][29]. Other studies were dedicated to the insertion of catheter during mini-invasive treatment of cerebrovascular disease [30][31][32][33] (Wang, 1997) (Lenoir, 2006) (Nowinski, 2001) (Cotin, 2015). Among them, (Lawton et al., 2000) described a computationally efficient Finite Element model based on an in-house resolution method [34]. The catheter was modeled as a thin rod and the vessels were assumed rigid like in the present study. (Schafer et al., 2009) presented a method to determine the shape of the guidewire using shortest path algorithms [35].

Each problem was unique and raised different challenges. Cerebrovascular intervention raised challenges by navigating into a complex maze of thin arteries. TAVI and endovascular aneurysm repair raised more concerns on the correct placement of a prosthesis in large vessels. During TAVI, unlike other medical contexts, the guidewire is pushed against the ventricle wall. The force exerted against the guidewire, its curvature and its properties could influence its shape. Beyond the insertion of a stiff guidewire alone, it must be kept in mind that the insertion of the prosthesis and its sheath could also deform the aortic arch. Models for cerebrovascular treatment navigation usually featured in-house implicit resolution schemes, which seemed to be successful when the vessels were rigid and little stress was generated. Conversely, models for aneurysm repair showed that explicit resolution was successful to predict the large deformations of the tools and the arteries. Likewise, simulation model for the deployment of prosthesis mainly used explicit resolution scheme so far. Despite the widespread use of explicit schemes for similar problems, we successfully used an implicit resolution scheme which is expected to provide more reliable results about stresses, strains and contact forces.

Interpretation of results

The main findings can be summarized as follows. Concerning the physical parameters, the results showed that the insertion force, orientation and curvature radius strongly impacted the guidewire deformation and position at the plane of aortic valve. The curvature angle and stiffness of the guidewire had a weak influence on the inserted shape. The tests over the numerical parameters showed that they had negligible impact on the accuracy. The base of the guidewire was assumed to have little impact on the accuracy at the aortic valve, so its position was arbitrarily chosen. As such, the base of the simulated guidewire did not match the intra-operative image. The comparison of the simulation in the ventricle was not conclusive. The intra-operative ventricle appeared radiotransparent so the quality of registration could not be assessed in this volume. However, all things considered, the model could reach very good accuracy in predicting the position and angle of the guidewire at the aortic valve.

Throughout the development of the model we observed convergence difficulties for insertion force above approximately 0.20 N using natural insertion method. A force displacement graph was plotted in Fig. 6.d and illustrated the domain where the simulation diverged. An analysis on the guidewire stiffness showed that the maximum insertion force was strongly correlated to the stiffness. However, the guidewire shape at the maximum insertion force was nearly identical whatever the stiffness. Using the forced insertion method, the simulation could overcome the instabilities and reach an insertion force between 0.3 N and 0.5 N. However, those simulations were even more different from the intraoperative 2D image, shifting the guidewire at the aortic valve more than 10 px. Beyond 0.5 N, the simulation encountered again convergence difficulties.

The insertion of a beam against a wall might be prone to buckling effects, which can create convergence difficulties for an implicit resolution method. Euler's critical buckling load is given by the formula:

(1)

Where L is the length of the beam, E is the Young's modulus of the beam (60 GPa), I is the area moment of inertia of the cross section (0.0491 mm 4 ), K is the column effective length factor wich was assumed to be equal to 2 as we considered one end highly constrained by contact force on the aortic wall (fixed) and the floppy tip of the guidewire free to move. Eq.1 is consistent with the correlation observed in set #4 between the stiffness of the guidewire and the maximum insertion force. In order to have a rough evaluation of the critical load for the guidewire, we measured that the length of the guidewire which may be subject to buckling was about 180 mm. Fig. 7 shows the major contact zones on the guidewire and the measurement (in green). We obtained a critical load around 0.22 N and the proportionality coefficient relating E to F ( equal to 0.0037 N/GPa. This theoretical critical load was in the same range order as the maximum axial forces reached before divergence of our simulations (~ 0.2 N). Coincidentally, the theoretical proportionality coefficient was also in the same range as the coefficient deduced from set #4 (0.0040 N/GPa).

Limitations

One of the main limitations of this study was related to the available data. The pre-operative 3D image was unsynchronized with heart cycle and the valve leaflets could not be extracted. The valve calcifications as seen in Fig. 2.b have been neglected and may impact the simulation. They could be added as rigid bodies. However, we observed that the simulated guidewire did not touch the segmented calcifications in those patient cases specifically. The shape of the ventricle was approximate due to heartbeat effects. The creation process of the patient mesh also increased the uncertainties on the geometry of the model. The model required information on the conditions of insertion that could not be obtained from the intra-operative data, such as the inserted length of the guidewire, the insertion force, the orientation and the shape of the curvature of the guidewire. Last, there was a possible mismatch between the model geometry and intra-operative image because of heartbeat and breathing motions. The 2D/3D registration was focused on the ascending aorta, where the contrast due to injection was highest. In the near future, the comparison of the numerical model with intra-operative images of a larger cohort of patients should provide a clearer evaluation of the reliability of the simulations and the validity of the assumptions. Finally, despite the importance of modeling the guidewire deformation, the prosthesis and the sheath should be included into the model in future studies. The Timoshenko beam formulation was not necessarily the most appropriate choice due to the occurrence of shear-locking in field-inconsistent elements. The effect of shear-locking on accuracy and robustness was mitigated by the use of 3-node quadratic element, but certainly affected the speed of the simulation. A possible direction for future works would be the improvement of the performances by including new element formulations.

Application perspectives

The simulation framework offered interesting insights on the behavior of the guidewire and provided data to improve TAVI simulation outcomes accuracy. The framework is currently not at the stage where it can provide strategies to completely correct misalignment of the prosthesis at the aortic valve but tries to contribute towards this goal. The prosthesis alignment is greatly determined by the sheath in comparison to the guidewire. The latter still contributes to both the alignment and the stability of the prosthesis. In exceptional cases, this small contribution is enough to allow experienced clinicians improving the alignment of selfexpandable prosthesis by manipulating the guidewire. The insights provided by this simulation framework could be applied to predict whether this kind of maneuver is feasible or beneficial for a given patient anatomy and a given guidewire curvature. In addition, the handling of the guidewire carries a risk of ventricle damage. The simulation framework could help studying the behavior of the contact between the guidewire and the ventricle to help choosing the guidewire curvature associated with the best stability and the safest distribution of contact forces on the ventricular wall. Currently, there is no general agreement regarding guidewire curvature. Increasingly more clinicians use pre-shaped guidewire (SAFARI TM ) and some center systematically use the smallest size for simplicity. Some clinicians report that small-sized pre-shape is prone to unwanted migration within the ventricle and provides less stability to the prosthesis deployment. The current simulation framework has only been assessed for geometrical parameters at the aortic valve. Future works could be directed towards the improvement of prediction in the ventricle. Finally, the centerline of the aortic valve is the starting point to visualize the prosthesis deployment. The simulated guidewire is a possible improvement to help the clinicians guess more accurately the alignment of the prosthesis in the pre-operative planning stages and approach more cautiously anatomies which produce severe misalignment and have higher risks of complication. As the model matures, sheath and prosthesis may be included to help the clinician making more accurate predictions. Retrograde trans-femoral access is associated with the highest survival rate, but there could be exceptional patient cases where other accesses provide safer positioning of the prosthesis. Future works could also compare the impact of vascular accesses on prosthesis positioning.

Conclusion

This paper presented a unique simulation framework to predict the angle and the position of a stiff guidewire in the aortic valve. The accuracy of the model was evaluated with patient intra-operative data and good agreement could be reached. The sensitivity analysis concluded that the shape of the guidewire was mainly determined by the patient geometry and the insertion conditions. The unique use of implicit resolution scheme in this particular context required stabilization techniques to guarantee convergence. However, divergence still occurred when the guidewire was subject to buckling. A larger number of patient cases will be simulated in the future to ascertain the reliability of this method.

Tab.1 -Variation of parameters organized as sets of simulation, each one focusing on the effect of an individual parameter. Fig. 3 -Illustration of the simulation stages. After applying initial conditions, the guidewire was constrained on the centerline. After the constraints were relaxed, the guidewire reached its final equilibrium. Then insertion force could be varied.

. . .

  .a) had a length of 40 mm. The length between the end and point C (Fig.1.a) was set to 120 mm in agreement to the usual practice in the Rennes University Hospital. This corresponded to the transition between the floppy tip and the stiff body of the guidewire (Fig.1.b).

  Fig.1 a) Model of guidewire pre-deformation; A: angle; R: radius of curvature; C: start of the curved shape; The tip has a rigid sphere to avoid convergence issues arising from contacts (c.f. section 3-d). b) Stiffness profile along the guidewire based on (Luboz et al., 2011)(Harrison et al., 2011). c) d) Definition of orientation of the tip of the guidewire.

  Fig.2 a) Intra-operative angiography of aortic root. b) Reformatted slice showing the aortic valve from pre-operative CT-scan. c) Finite Element mesh from segmentation of CT-scan d) Segmentation from CT-scan

  Fig.4 a) Projection of 3D pre-operative geometry on the intra-operative 2D view according to the registration transforms. b) Intra-operative guidewire is drawn with a yellow line; simulated guidewire appear with a cyan line. The aortic valve plane is manually drawn. The distance D and angle A of the guidewires are measured at the intersection with the valve plane.

Fig. 5 -

 5 Fig.5 -Projection of simulated guidewires. a) Set #1 exploring angle curvature at 0.10 N insertion force. b) Set #1 exploring angle curvature at 0.20 N insertion force. c) Set #2 exploring radius at 20° angle curvature. d) Set #3 exploring radius at 80° angle curvature. e) Set #5 exploring orientation

  Fig.6 a) b) c) Projections of simulated guidewires at the edges of the convergence domain reported in d).Natural insertion simulations matched closer to intra-operative data than forced insertion simulations. d) Force-displacement graph of the insertion of guidewire. Forced insertion method allows exploring isolated domain of convergence.

  Fig.7a) Contact force distribution along the simulated guidewire from set #6 (fig.6a). The green circles highlight contact zones which may produce a buckling effect. b) Measurement of the portion of the guidewire hypothetically subjected to buckling.

  Tab.3 -Accuracy evaluation of simulation sets from the projection of simulated guidewire on intraoperative image. Tab. 4 -Detailed measurements for sets #1 and #4. a) Comparison of insertion force 0.10 N and 0.20 N (set #1). b) Maximum force insertion appears linearly dependent on stiffness of the guidewire (set #4).

	Physical parameters Set 1			Simulation angle	2				radius	3	radius
	Parameter Focus s a)	Angle difference (°)	Distance (px)	1		2 Angle difference (°) 3	4 Distance (px)	5 difference (°) Distance (px) 6 Angle 7-12 orien-insertion numerical
	min Set #1		angle Insertion force 0.10N radius radius stiffness 3.6 2 -0.5 Insertion force 0.20N tation 1	depth -0.4	param.	0
	Guidewire properties max Angle difference 4.4		3 Distance		-2.3 Angle difference 6	Distance	-1.6	5
	angle (°) Angle (°) mean	(°)	3.9	20-120 2.3 (px)	20		80 -0.8 (°)	20 3.3	20 (px)	-0.7	20	20 2.3
	std	radius (mm) 20		0.3	30 15-30 0.3 3.6	15-30 1.3 3	30 2.2 0.9	25	0.9 2	25	30 2.2
		base stiffness (GPa) 40	60 3.6			60	2	60	20-100	1	60	1	60	60
	Insertion conditions Set 4 60			stiffness 3.8		5		2		orientation 6 0.9	1	depth
	Parameter	Angle 80			Distance 4		0.1-Angle	3		Distance	0.5	Angle	0
	s	force (N) difference (°) 100	0.1-0.2 (px) 3.9		0.2 0.1-0.2 difference (°) 2	(px)	max	0	0.1-0.2 difference (°) Distance (px) max 0.2 1
	orientation (°) initial depth (mm) 120 Mean min max	0 0.6	90 -40 4.4 3.9	4 5	90 -40	2 2.3	90 -40 -0.4 2.8	90 -40 13 3 0.2 0.6	0-180 -40	-0.3 1 -9.9 1.0	90 7	90 0 -40 18
	mean Standard			0.1		4.6				0.3	5.2	-5.1	9
	Numerical parameters std deviation	0.3	0.5 0.3			0.5	2	4.4 0.4	6.8 0.6	12.7
						1-6			7		8	9	10	11	12
	Focus Set b) Parameter Set #4 Contact conditions 7 Angle s difference (°) Stiffness (GPa) Angle difference (°) Distance (px) Maximum force (N) physical param. FKN FTOLN mesh decay ratio FKN 8 FTOLN 9 Distance Angle Distance Angle (px) difference (°) (px) difference (°) Distance (px) decay / loadstep substep / predictor mesh stiffness coefficient (N/mm3) 1 0.1-10 1 1 0.1 0.1 min 0 0 0.5 2 1 20 0.1 5 0.08 1 penetration tolerance 2 2 0.5-4 2 2 2 max -0.4 0 1.5 2 -0.5 1 40 0 5 0.16 2 mean -0.3 0 1 2 -0.2 1 60 -0.2 5 0.24 Model & solver parameters Mesh element size 0.4 0.7 0.7 0.4-0.7 0.7 0.7 0.7 std 0.2 0 0.5 0 0.3 0 80 0.6 4 0.32
	Decay ratio 100 Number of loadsteps Number of substeps Predictor difference (°) 10 Mean s Set Parameter Angle Standard deviation	0.87 0.87 0.1 125 125 1 1 on off difference (°) 0.87 4 125 1 on 11 0.1 4.6 Distance (px) decay Angle 0.3 0.5	0.87 0.75-0.9 0.39 125 125 1 10 on on difference (°) Distance (px) 0.7-.87 0.87 50-125 125 1 01-10 on on/off decay / (px) loadstep 12 substep / predictor Distance Angle
	min				1.4			4				1.7	4	0.7	2
	max				2.1			4				2	4	1	2
	mean				1.8			4				1.9	4	0.9	2
	std				0.3			0				0.2	0	0.2	0
								.