Kyangwi P Malikidogo

Agnès Pallier

Frédéric Szeremeta

Célia S Bonnet
email: celia.bonnet@cnrs.fr

Olivier Sénèque
email: olivier.seneque@cea.fr

A phosphorylated

come

Introduction

Zinc is the second most abundant transition metal ion in human body after iron. 1,2 It is mostly bound to proteins that play a central role in controlling gene transcription and metalloenzyme function. Zn 2+ is a signaling ion involved in several pathways such as the immune system, 3 glutamatergic neurons, fertilization 4 and early development processes. 5,6 Zn 2+ concentration is highly regulated and disruptions in its homeostasis have been implicated in several diseases such as pancreatic, prostatic and breast cancers, 7 diabetes 5 or neurodegenerative diseases, to cite a few. Therefore the non-invasive visualization of Zn 2+ in vivo is of prime importance for the early detection of these diseases 8 and more generally to better understand the role of zinc in biomedical research.

Magnetic resonance imaging is a powerful technique due to its high resolution, which plays a central role in molecular imaging. The last years have witnessed tremendous developments of MRI contrast agents sensitive to Zn 2+ , pushed forward by its successful in vivo detection. 9 Paramagnetic Gd 3+ chelates are by far the most exploited contrast agents for Zn 2+ detection. In most of the examples, Gd 3+ is complexed by a macrocyclic derivative of DO3A (1,4,7,10 tetraazacyclododecan-1,4,7,10 tricarboxylic acid), to which a Zn 2+ sensitive moiety such as dipicolylamine, iminodiacetate or quinoline is linked. 10 The efficacy, called relaxivity, of such agents is mainly influenced by the number of water molecules directly coordinated to Gd 3+ , its water exchange rate with the bulk, and the rotational correlation time of the complex influenced by the size and rigidity of the system. Upon Zn 2+ binding, at least one of these microscopic parameters must be modified to give a responsive contrast agent. The most successful contrast agent used in vivo relies on a change of the rotational correlation time due to the formation, in the presence of Zn 2+ , of a ternary complex with Human Serum Albumin (HSA), which is the most abundant protein in the blood. Such response is maximum at medium magnetic field (20-80 MHz), where the effect of the rotation is the most important. As relaxivity is directly proportional to q, the number of water molecule coordinated to Gd 3+ , the modification of this parameter gives rise to dramatic changes upon Zn 2+ binding that are similar over the whole field range. Such strategy has been explored successfully to achieve high relaxivity increase at 7T, 11 where the resolution is better. However, the prediction of the modification of q in the presence of Zn 2+ remains very difficult and many synthetic efforts should be done to obtain one successful complex.

We have recently developed, in a completely different approach, a Zn 2+ sensor for dual MRI and luminescence detection based on a zinc finger scaffold (LZF2 Ln , Ln = Tb or Gd, see Fig. S1 of SI for its sequence). 12 The folding of the peptide upon Zn 2+ binding results in a modification of the rotational correlation time of the system that is detectable by MRI. Inspired by the work of Vázquez et al. 13 for the detection of phosphorylation using a luminescent Tb 3+ complex, we anticipated that using this versatile zinc-finger scaffold, we could develop a system based on a change in the hydration number q, taking advantage of the replacement of a Gd-bound water molecule by the phosphate group of a phosphoserine, similarly to Vazquez's probe.

We have therefore synthesized two candidate contrast agents, ZFQD Gd and ZFQE Gd . They are 27-amino acid classical  zinc finger (ZF) peptides comprising a phosphoserine and functionalized by a Gd 3+ DO3Apropylamine chelate. These two contrast agents differ only by the amino acid to which the Gd 3+ chelate is attached, which is either an aspartate (D) or a glutamate (E) (Fig. 1A). Classical  zinc fingers have no defined conformation in their Zn 2+ -free form while they adopt a compact  fold upon Zn 2+ binding to two histidines and two cysteines. [14][15][16] The DO3A propylamine ligand leaves available a coordination site for phosphoserine binding to the lanthanide ion, and we anticipated that the conformational changes upon Zn 2+ binding would modify the Gd 3+ -phosphoserine distance resulting in a change in the Gd 3+ coordination.

Results and discussion

Synthesis of Ln-ZF complexes:

The sequence of ZFQD or ZFQE is based on the consensus sequence of Zinc Fingers described by Berg et al. 14 The synthesis of ZFQE is described in Fig. 1B. The peptide was elongated on Rink Amide resin following standard Fmoc/tBu solid-phase synthesis procedures (SPPS). 18 A glutamate amino acid with its side chain protected as an allyl ester, orthogonal to other protecting groups, was insert after the Nterminal cysteine. The peptide was acetylated after the N-terminal lysine. The removal of Alloc group was performed using Pd 0 . 19 The tBu-protected DO3A-propylamine pro-ligand, which was synthesized according to literature procedures, 13,17 was then coupled to the free carboxylic acid of the peptide using PyBOP/DIEA activation. Finally, the peptide was deprotected and cleaved from the resin in a mixture of trifluoroacetic acid (TFA) and scavengers, purified by reverse phase high performance liquid chromatography (HPLC) and characterized by ESI-MS (Fig. S2 of SI). Ln 3+ complexes were obtained by reaction of ZFQE with LnCl3, xH2O in water pH 6.2 overnight. ZFQD Tb was synthesized the same way only, an Asp residue was used instead of a Glu.

Zn 2+ binding properties: First, the Zn 2+ binding properties of the ZFQE Tb and ZFQD Tb were investigated by circular dichroism (CD) and luminescence spectroscopy. CD spectra in the absence and in the presence of 1.5 equivalent of Zn 2+ were recorded in water containing 250 µM tris(2-carboxyethyl)phosphine (TCEP) at pH 7.4 (Fig. 2A and2C). The CD spectra of both Zn 2+ -free peptides are characteristic of an unstructured (random coil) classical  zinc finger peptide with an intense negative CD signal at 200 nm and a small shoulder starting at ca. 225 nm. 20 Addition of excess of Zn 2+ induces a change in the CD spectra with a minimum at 205 nm and a shoulder at 220 nm, which is characteristic of Zn-bound  zinc fingers. 20 This indicates that the presence of both Tb 3+ (DO3A-propylamide) complex and phosphoserine does not alter the behaviour and folding ability of the zinc finger peptide. Zn 2+ binding to ZFQE Tb and ZFQD Tb was investigated in HEPES buffer (10 mM, pH 7.4) containing TCEP (250 µM) by tryptophan fluorescence and Tb 3+ luminescence under 280 nm excitation (Fig. 2B and D; Fig. S3 and S4 of SI). For ZFQE Tb , tryptophan fluorescence displays almost no change upon Zn 2+ addition (Fig. S3 of SI). Contrastingly, the sensitized Tb 3+ luminescence of ZFQE Tb decreases linearly upon addition of Zn 2+ up to approximatively 0.9 equivalent and then remains constant upon further addition of Zn 2+ , thus monitoring formation of the expected 1:1 Zn/peptide complex (Fig. 2B). A 55% quenching of Tb 3+ emission is observed upon Zn 2+ binding to the peptide. For ZFQD Tb , neither tryptophan fluorescence nor Tb 3+ emission showed variation in their intensity upon Zn 2+ binding, despite Zn 2+ binding observed by CD. In order to understand these results, the hydration number, q, was determined for the Zn-free and Zn-bound forms of ZFQE Tb and ZFQD Tb by difference in Tb 3+ luminescence lifetimes in H2O and D2O (Table 1 and Fig. S5 of SI). Regarding ZFQD Tb , both the Zn-free and Zn-bound forms displayed mono-hydrated Gd chelates (q  1.0), consistent with the absence of variation in the luminescence emission of this probe. However, q values of 0.4 ± 0.2 and 1.0 ± 0.2 were obtained for ZFQE Tb and Zn-ZFQE Tb , respectively. The increase of q upon Zn 2+ binding to ZFQE Tb , is consistent with the decrease of the Tb 3+ luminescence intensity observed, indicating a change in the coordination sphere of Tb 3+ upon Zn 2+ binding. It is known that phosphate groups bind Ln 3+ and other cations in a monodentate fashion. 21,22 Therefore, we can propose that in the Zn-bound form, a water molecule is coordinated to the Tb 3+ ion in addition to the four cyclen nitrogen atoms and the four carbonyl oxygen atoms of the pendent arms, while in the Zn-free form, an equilibrium may take place between two species, one with the phosphoserine bound to Tb 3+ and the other with a water molecule coordinated instead (Fig. 3). Contrastingly, for ZFQD Tb , which has a shorter anchorage for the Tb 3+ chelate (one less methyl), the phosphoserine does not coordinate Tb 3+ in the presence or in the absence of Zn 2+ . This indicates that the length of the anchorage is crucial to observe the coordination of the phosphoserine and a significant q change. Table 1: Luminescence lifetimes and derived q values for free and Zn-loaded ZFQD and ZFQE probes in HEPES buffer (10 mM, pH 7.4) and 250 µM of TCEP. Error on lifetime is estimated ± 0.03 and error on q is estimated ± 0. Fig. 3: Principle of Zn 2+ detection: Folding of the ZFQE Ln probe and variation of q induced by Zn 2+ -binding. In the Zn 2+ free state, an equilibrium between q = 0 and q = 1 may take place.

Magnetic properties:

In light of these encouraging results, we envisioned that replacing the Tb 3+ -ion in ZFQE Tb by a Gd 3+ ion, a Zn 2+ -responsive MRI probe would be obtained. A Zn 2+ titration of ZFQE Gd in HEPES buffer pH 7.4 monitored by 1 H relaxivity at 60 MHz and 25°C was performed (Fig. 4A). The relaxivity increases up to one equivalent of Zn 2+ added and remains constant above this limit. This result is consistent with the luminescence titration of ZFQE Tb (Fig. 2B) presented above and confirms the formation of 1:1 Zn 2+ /ZFQE Gd complex. The relaxivity of the system at 60 MHz and 25°C in the absence of Zn 2+ is 10.48 mM -1 s -1 , and 25.27 mM -1 s -1 in the presence of 1 eq of Zn 2+ , which represents an increase of more than 140%. As expected, this is better than the previous zinc-finger system that we developed and based solely on a change of the rotational correlation time of the system. But it ranges also among the best system found in the literature, especially based on a change of the number of water molecule coordinated. 10,23,24 T using spin echo sequence with TE = 9.9 ms, TR = 50 ms.

In order to better understand the origin of this relaxivity response, the 1 H NMRD profiles of ZFQE Gd were recorded at two temperatures and over a wide range of magnetic fields (0.1 to 500 MHz) (Fig. 4C). First, we can see that both in the absence and in the presence of Zn 2+ , the system shows the typical hump of slowly rotating species. This was expected in the presence of zinc, as there is one water molecule coordinated to Gd 3+ , and the system is certainly endowed with a high rotational correlation time due to the folding. In the absence of zinc, this confirms that there is certainly an exchange between 0 and 1 water molecule coordinated and/or participation of second sphere water molecules to the overall relaxivity. This is consistent with the luminescence results (q = 0.4).

Conversely to what is normally obtained for systems based on a change in the hydration number, the relaxivity increase is not constant over the whole range of magnetic field. The maximum increase is obtained at medium field (40-60 MHz; ca. 150 %), while at high field (400 MHz), the relaxivity increase is 50 and 67 % at 25°C and 37°C, respectively. This shows that both the number of water molecules coordinated to Gd 3+ , but also the rotational correlation time (R), linked to the size and the rigidity of the complex plays a role in the response upon Zn 2+ binding.

This is interesting as the coordination of the phosphoserine certainly plays a role in the preorganisation of the system in the absence of Zn 2+ , but the folding in the presence of Zn 2+ still gives rise to a sizeable response. The full fitting of the data in the absence of zinc is not reasonable to obtain a reliable picture of system due to the possible dynamic of exchange. In the of Zn 2+ , it is impressive to note that the relaxivity of ZFQE Gd at 40

MHz 25°C is nearly twice that of LZF2 Gd (27.56 vs 14.64 mM -1 .s -1 ), and it is 30 % higher at 37°C (20.3 vs 15.75 mM -1 .s -1 ). 12 In the case of LZF2 Gd , 12 we had shown that the water exchange rate was limiting the relaxivity in the presence of Zn 2+ , which is not the case here as the relaxivity decreases with increasing temperature even in the presence of Zn 2+ . This was expected as the water exchange rate of this monohydrated DO3A-alkylamide system is two orders of magnitude faster than that of DO3A-monoamide. 25 Interestingly, the fitting of the data using the Solomon-Bloembergen and Morgan theory was possible without the incorporation of the Lipari-Szabo treatment, which underlines the very high rigidity of the system. The best fit parameters obtained from the analysis of the 1 H NMRD data are summarized Table 2 and Table S1. The fitted curves are presented Fig. 4C. The rotational correlation time obtained matches well with the global rotational correlation time obtained for LZF2 Gd . It shows that for ZFQE Gd in the presence of Zn 2+ , there is no loss of relaxivity due to internal flexibility of the system. This could be rationalized by the formation of a hydrogen bond between the phosphoserine and the alkylamide that would lock the Gd 3+ complex (Fig. 3). The very high relaxivities achieved here in the presence of Zn 2+ are due to a combination of a very rigid system with a high rotational correlation time and a fast water exchange rate. Moreover, the Zn 2+induced response is not influenced by the presence of physiological cations such as Ca 2+ , Mg 2+ , Na + or K + (Fig. 4B). The addition of biologically relevant concentration of Cu 2+ (0.15 eq.) does not modify the response, while the addition of 1 eq results in a 20 % relaxivity decrease (Fig. S6 of SI). Finally, phantom images were performed in vitro at 9.4 T on a MRI scanner and are presented Fig. 4D. A clear signal intensity enhancement for Zn-ZFQE Gd vs ZFQE Gd can be seen, and the T1 measurements (Fig. S7) performed in these conditions unambiguously confirm these results.

Conclusion

To conclude, we have designed a Zn-responsive contrast agent based on a zinc finger peptide. The high response achieved is based on a switch in the phosphoserine coordination to the Ln 3+ , which gives rise to a q change, and is amongst the best reported to date. The in vitro phantom images clearly show a sizeable increase of MRI efficacy at high field in the presence of Zn 2+ , due to a combined effect of an increase in the Gd 3+ hydration number and in the rigidity of the molecule. Efforts to optimise the probe response are underway. The versatility of the peptide scaffold is an undeniable asset to achieve such goal with minimal synthetic efforts as illustrated by the change of a single amino acid to alter the response. The zinc finger scaffold offers a unique opportunity to design and optimise various zinc-responsive probes.

Materials and methods

Reagents and solvent: All commercial reagents were used as received from suppliers unless otherwise indicated. NovaPEG Rink Amide resin, N-α-9-fluorenylmethoxycarbonyl (Fmoc)-protected amino acids for peptide synthesis 2-(6-Chloro-1-H-benzotriazole-1-yl)-1,1,3,3-tetramethylaminium hexafluorophosphate (HCTU) coupling reagent were obtained from Novabiochem or Iris Biotech. Cyclen was purchased from CheMatech. Other chemicals and solvents were purchased from Sigma-Aldrich.

Analyses and purification: TLC sheets ALUGRAM SILG/ ultraviolet (UV 254 ) were used for thin layer chromatography and visualized using UV light. Purification of organic compounds was performed on classic silica gel chromatography. 1 H and 13 C NMR spectra were recorded in deuterated solvents on a Bruker UltraShield Plus 400 Spectrometer operating at 400 and 100 MHz for 1 H and 13 C, respectively. All experiments were performed at 25 °C. Analytical HPLC separations were performed on an Agilent 1260 Infinity II system using XBridge BEH C18 (2.5µm, 4.6 x 7.5 mm) columns at 1 mL/min equipped with MS module. Preparative HPLC separations were performed on a VWR LaPrep system using a Waters XBridge Peptide BEH130 C18 (5 µm, 150 mm x 19 mm) column at 14 mL/min. Mobile phase consisted in a gradient of solvent A (0.1% trifluoroacetic acid (TFA) in H2O) and B (0.1% TFA in acetonitrile (MeCN)/H2O 9:1). For analytical separations, Method A consisted in 5% B during 2 min followed by a 5 to 100 % B gradient in 13 min. Eluate was monitored by electronic absorption at 214, 254 and 280. LRMS (ESI) analyses were performed on a Thermo LXQ spectrometer. Circular dichroism (CD) spectra were recorded on an Applied Photophysics Chirascan spectropolarimeter equipped with a thermo-regulated cell holder. Emission and excitation spectra and Tb 3+ luminescence lifetimes were measured on a Varian Cary Eclipse spectrometer equipped with a thermoregulated cell holder. Time-resolved Tb 3+ luminescence spectra were acquired with 100 µs time delay. All measurements were performed with aerated solutions.

Peptide synthesis: Peptides were synthetized on NovaPEG Rink Amide resin (0.46 mmol/g; 0.12 mmol scale) using standard SPPS protocols for Fmoc/tBu chemistry. The first amino acid was attached by single manual coupling, using 3-fold excess of Fmoc-Gly-OH, 2.7-fold excess of HCTU and 6-fold excess of N,Ndiisopropylethylamine (DIEA) in N,N-dimethylformamide (DMF). This was followed by a capping step (Ac2O/pyridine/DMF 1:2:7 (v/v/v), 10 mL, 5 min). Removal of the Fmoc protecting group was performed with 20 % piperidine in DMF (3×10 min). Peptide elongation was performed on an automatic peptide synthesizer (CEM Liberty1 Microwave peptide Synthesizer) using double coupling with 4-fold molar excess of Fmoc-L-amino acid, 3.8-fold molar excess of HCTU and 8-fold molar excess of DIEA during 30 min at room temperature. A capping step was performed after each coupling with Ac2O/DIEA (5 min). Fmoc removal was performed using 20% piperidine in DMF (2×10 min). The phosphoserine was introduced using Fmoc-Ser(PO(OBzl)OH)-OH. Fmoc-L-Glu(OAllyl)-OH and Fmoc-L-Asp(OAllyl)-OH were used to introduce the orthogonally-protected Glu/Asp residue used to graft the chelator on resin for ZFQE and ZFQD respectively. After the final Fmoc removal step, the N-terminus was acetylated by treatment with Ac 2 O/pyridine/DMF 1:2:7 (v/v/v, 10 mL, 5 min).

Removal of glutamate or aspartate allyl protecting group:

The removal of the allyl protecting group was performed by treating the resin twice with tetrakis(triphenylphosphine)palladium(0) (Pd(PPh3)4) (0.02 mmol, 27.6 mg) and phenylsilane (2.5 mmol, 0.36 mL) in degassed anhydrous DCM (15 mL) for 1 h in the dark. 19 The resin was then washed successively with DCM (4×3 min), DMF (2×3 min), 1% H2O in DMF (2×3 min), DMF (2×), 1% DIEA in DMF (2×3 min), DMF (2×3 min), sodium diethyldithiocarbamate in DMF (0.12 M, 2×5 min), DMF (2×3 min), DCM (4×3min).

Coupling of the (tBu)-protected macrocyclic ligand: DO3A(tBu)3-propylamine 26 (0.24 mmol, 137.5 mg) was dissolved in DMF (6 mL) and added to the resin (60 µmol) for ZFQE, then a solution of PyBOP (0.24 mmol, 71.15 mg) and DIEA (0.48 mmol, 250.3 µL) in DMF (3 mL) was added. The resin was agitated overnight at room temperature, then washed with DMF (2×), DCM (2×). The same procedure was used for ZFQD (20 µmol) and the reagents were mixed in accordance to the number of equivalents as for the ZFQE presented above (80 µmol of DO3A(tBu) 3 -propylamine derivative, 80 µmol of PyBOP and 0.16 mmol of DIEA in DMF (1 mL)).

Resin cleavage and removal of protecting groups: Removal of acid-labile protecting groups and resin cleavage were performed using TFA/H2O/triisopropylsilane/thioanisole (92.5:2.5:2.5:2.5 v/v/v/v, 20 mL) containing dithiotreitol (300 mg) for 4h. TFA was evaporated under reduced pressure and cold diethylether (Et2O) was added to precipitate the peptide, which was purified by HPLC to give ZFQE (42 mg, 18 % yield) or ZFQD (10 mg, 14 % yield) after freeze-drying. ZFQE: HPLC (anal.) tR = 6.71 min (method A); LRMS (ESI+) : average m/z = 1218.9 (3+), 914. Metallation with Ln 3+ : The zinc finger peptide (5.66 µmol, 20.7 mg) and LnCl3 (Ln = Gd or Tb; ca. 10 eq.) were dissolved in H2O and the pH was adjusted to 6.2 using NaOH. The solution was stirred overnight (after 30 min, the pH was controlled and adjusted to 6.2 if needed). Excess of Ln 3+ was removed using OASIS HLB 3cc (400 mg) LP Extraction Cartridge and the Ln complex was obtained as a white powder after freeze-drying (ca. Circular dichroism: Solutions of ZFQE Tb (~ 20 µM) and ZFQD Tb (~30 µM) without Zn 2+ and with 1.5 eq. Zn 2+ were prepared in H2O using a stock solution of the complex. TCEP (250 µM) was added as a reducing agent to ensure cysteine reduction and the pH was adjusted at 7.5 using NaOH. Cuvettes with a 0.4 cm path length were used. Spectra were recorded from 300 nm to 200 nm every 1 nm with a 2 s signal averaging for each point. Each spectrum was recorded twice, averaged and smoothed using a Stavitsky-Golay filter.

Luminescence spectroscopy: Solutions of ZFQE Tb and ZFQD Tb were prepared in a HEPES buffer (10 mM, pH 7.4) containing TCEP (250 µM) as a reducing agent. Typically, the concentration of ZFQE Tb and ZFQD Tb were 20 µM and 100 µM respectively prepared using a stock solution. Zn 2+ titrations were performed using a Zn(ClO4)2 solution (2.09 mM in water). For the determination of q, the number of water molecules coordinated to Tb 3+ , solutions of ZFQE Tb (10 µM) and ZFQD Tb (27 µM) in a HEPES buffer (10 mM, pH 7.4) containing TCEP (250 µM) were prepared in various H2O/D2O mixtures (25, 50, 75 and 100 % H2O) without and with 1.5 eq. of Zn 2+ .

Lifetimes in 100 % D2O were extrapolated from the plot of the rate constants of Tb 3+ luminescence decay (kLn =  Ln -1 ) against the fraction of H 2 O in H 2 O/D 2 O mixtures (Fig. S3). The number of water molecules coordinated to the Tb 3+ ion (q) was determined using q = 5 x (kTb/H2O -kTb/D2O -0.06). 27 Sample preparation for relaxometric studies: All samples were prepared in a HEPES buffer (100 mM, pH 7.4) containing TCEP (25 mM) as a reducing agent. The concentration of ZFQE Gd was determined using the extinction coefficients of tryptophan (ε= 5690 M -1 cm -1 at 280 nm) and extinction coefficients of tyrosine (ε= 1480 M -1 cm -1 at 280 nm) and the concentrations of Gd 3+ -containing solutions were systematically checked by MBS.

For Zn 2+ titration, 1 H relaxivity measurements were performed at 60 MHz and 25°C with [ZFQE Gd ] = 0.69 mM after addition of a ZnCl2 solution.

NMRD Profiles: All measurements were performed on a Stelar SMARTracer Fast Field Cycling relaxometer (0.01-10 MHz) and a Bruker WP80 NMR electromagnet adapted to variable field measurements (20-80 MHz) and controlled by a SMARTracer PC-NMR console. The temperature was monitored by a VTC91 temperature control unit and maintained by a gas flow. The temperature was determined by previous calibration with a Pt resistance temperature probe. The longitudinal relaxation rates (1/T1) were determined in water. At 400 MHz and 600 MHz, 1/T1 measurements were performed in a capillary tube placed in a 5 mm tube containing D2O using the inversion recovery sequence on a and a Bruker Avance III (400 MHz) equipped with a 5 mm BBO SmartProbe, and a Bruker Avance DMX (200 MHz) equipped with a 5 mm BBIZ probe, respectively. A Bruker Avance 300 MHz and a Bruker Avance III equipped with a 5 mm BBI were used for 300 MHz and 500 MHz respectively. The leastsquares fit of the 1 H NMRD data were performed using Visualiseur/Optimiseur 28,29 running on a MATLAB 8.3.0 (R2014a) platform using equations derived below. In the fitting procedure, we have fixed the rGdO distance to 2.50 Å, based on available crystal structure and ENDOR results, the rGdH distance to 3.10 Å and the closest approach of the bulk water protons to the Gd 3+ , aGdH to 3.5 Å. 30 The quadrupolar coupling constant, (1+ 2/3 ) 1/2 , has been set to the value for pure water, 7.58 MHz. 31 The diffusion constant has been fixed to 2310 -10 m 2 s -1 and its activation energy to 10 kJ mol -1 . The hydration number was 1 as determined by luminescence lifetime measurements on the Tb 3+ analogue. The water exchange rate and the water exchange enthalpy were fixed to the values obtained for similar systems (kex 298 = 11110 -6 s -1 and H ≠ = 21 kJ mol -1 ). 25 The activation energy of the correlation time for the modulation of the zero-field splitting EV has been fixed to 1 kJ.mol -1 . The fitting has been restricted to frequencies above 4.5 MHz as at low magnetic fields the SBM theory fails in describing electronic parameters and rotational dynamics of slowly rotating objects. Phantom images at 9.4 T: Phantoms of ZFQE Gd , Zn-ZFQE Gd and HEPES have been imaged together with a classical Bruker birdcage coil, 35 mm inner diameter, on a BioSpec 9.4 T spectrometer (Bruker, Wissembourg, France). A spin echo sequence was used, to acquire images of each sample with a fixed TE of 9.9 ms and a variable TR ranging from 50 ms to 3000 ms. Using ISA-Tool program, provided with Paravision 5.1, the T1 relaxation time of each phantom was calculated, using an exponential fitting, from signal mean intensity measured on the images as function of TR. Background mean intensity, referred to as noise, was also measured. Figure S5 shows the measured mean intensities as function of TR. As noise remains constant, ISA-Tool gives good T1 estimation. ZFQE Gd , Zn-ZFQE Gd and HEPES have T1 values of 655 ms, 526 ms, and 2210 ms, respectively.

Fig. 1 :

 1 Fig. 1: (A) Amino acid sequence of ZFQX Ln (X = D or E, Ln 3+ = Tb 3+ or Gd 3+ ). Zn 2+ -binding cysteines and histidines are shown in red. The amino acid bearing the Ln(DO3A-propylamide) complex is shown in brown. Phosphate group in blue was incorporated in the peptide sequence to induce a variation of q. (B) Synthetic pathway for the preparation of ZFQE Ln (or ZFQD Ln ) using Fmoc SPPS strategy. Conditions: a) Pd(PPh3)4, phenylsilane, degassed anhydrous DCM, 2x 1h, at room temperature; b) DO3A(tBu)3-propylamine 13,17 , PyBOP, DIEA, DMF, overnight at room temperature; c) TFA/H2O/triisopropylsilane/thioanisole/dithiotreitol, 4h at room temperature; d) LnCl3.xH2O, H2O, pH 6.2, overnight. ZFQD was obtained following the same procedure but using Fmoc-Asp(OAllyl)-OH instead of Fmoc-Glu(OAllyl)-OH.

Fig. 2 :

 2 Fig. 2: (A) CD spectroscopic characterization of ZFQE Tb (20 µM) without Zn 2+ (black) and in the presence of 1.5 eq. of Zn 2+ (red). (B) Evolution of time-resolved Tb 3+ luminescence emission spectra (λex= 280 nm, delay = 100 µs) upon addition of Zn 2+ in a solution of ZFQE Tb (20 µM). (C) CD spectroscopy characterization of ZFQD Tb (10 µM) without Zn 2+ (black) and in the presence of 1.5 eq. of Zn 2+ (red). (D) Evolution of time-resolved Tb 3+ luminescence emission spectra (ex = 280 nm) upon addition of Zn 2+ solution in ZFQD Tb . Solutions were prepared in unbuffered water pH 7.4 for CD and HEPES buffer (10 mM, pH 7.4) for luminescence. All samples contained TCEP (250 µM).

Fig. 4 :

 4 Fig. 4: (A) 1 H relaxivity measurements (60 MHz, 25°C) of a solution of ZFQE Gd (0.69 mM) in HEPES buffer solution (100 mM, pH 7.4, 25 mM TCEP) with increasing amount of Zn 2+ . (B) Selectivity plot obtained by measuring the relaxivity of ZFQE Gd (0.38 mM) alone (left) or in the presence of various physiological cations (1eq.) before (black) and after (grey) addition of Zn 2+ (1 eq.). (C)1 H NMRD profiles at 25°C (black symbols) and 37°C (red symbols) for Zn-free (square) and Zn-loaded (circle) ZFQE Gd (0.69 mM). The solid lines represent the least squares fit of the data points using the SBM theory. Solutions were prepared in HEPES buffer (100 mM, pH 7.4) containing TCEP (25 mM). All solutions were prepared in HEPES buffer solution (100 mM, pH 7.4, 25 mM TCEP). D) T1-weighted MR images of phantoms containing ZFQE Gd (240 µM) in the absence or in the presence of 1.5 eq of Zn2+ , in a HEPES buffer (100 mM, pH = 7.4) containing TCEP (30 mM). Images were acquired at 9.4

  3 (4+), 731.5 (5+), 609.6 (6+), 522.6 (7+) / calculated av. m/z = 1219.68 [M+3H] 3+ , 915.03 [M+4H] 4+ , 732.21 [M+5H] 5+ , 610.34 [M+6H] 6+ , 523.29 [M+7H] 7+ for M = C157H245N46O49PS2; deconvoluted mass found = 3656.0 / expected mass = 3656.05 (average isotopic composition). ZFQD: HPLC (anal.) tR = 6.71 min (method A); LRMS (ESI+): average m/z = 1214.8 (3+), 911.3 (4+), 729.3 (5+), 607.8 (6+), 521.3 (7+) / calculated av. m/z = 1214.80 [M+3H] 3+ , 911.35 [M+4H] 4+ , 729.28 [M+5H] 5+ , 607.90 [M+6H] 6+ , 521.20 [M+7H] 7+ for M = C156H243N46O49PS2; deconvoluted mass found = 3641.4 / expected mass = 3642.03 (average isotopic composition).

  65 % yield). LRMS (ESI+): average m/z = 1271.08 (3+), 953.83 (4+), 763.33 (5+), 636.33 (6+), 545.42 (7+) / calculated av. m/z = 1271.09 [M+3H] 3+ , 953.57 [M+4H] 4+ , 763.06 [M+5H] 5+ , 636.05 [M+6H] 6+ , 545.33 [M+7H] 7+ for M = C 157 H 242 N 46 O 49 PS 2 Gd; deconvoluted mass found = 3810.24 / expected mass = 3810.28 (average isotopic composition). ZFQE Tb : LRMS (ESI+): average m/z = 1271.7 (3+), 954.1 (4+), 763.6 (5+), 636.4 (6+), 535.0 (7+) / calculated av. m/z = 1271.65 [M+3H] 3+ , 953.99 [M+4H] 4+ , 763.39 [M+5H] 5+ , 636.33 [M+6H] 6+ , 545.57 [M+7H] 7+ for M = C157H242N46O49PS2Tb; deconvoluted mass found = 3811.8 / expected mass = 3811.96 (average isotopic composition).

Table 2 :

 2 Best-fit parameters obtained from the fitting of the 1 H NMRD profiles to the SBM theory. The values correspond to Er for Zn-ZFQE Gd and Eg for Zn-LZF2 Gd ; b The values correspond to  r for Zn-ZFQE Gd and  g

		Zn-ZFQE Gd Zn-LZF2 Gd
	r1 (40 MHz, 25°C)/ mM -1 .s -1 27.56	14.69
	E / kJ.mol -1 a	22 ± 4	13
	 298 / ps b	990 ± 50	1160
	E l / kJ.mol -1	-	20
	 l 298 / ps	-	145
	a for Zn-LZF2 Gd

Acknowledgements

We acknowledge the support from INBS France Life Imaging (grant ANR-11-INBS-0006 from the French ''Investissements d'Avenir'' program), Agence Nationale pour la Recherche (ANR-22-CE44-0041), DRF impulsion program of the CEA, the Labex Arcane, CBH-EUR-GS (ANR-17-EURE-0003), ITMO Cancer of Aviesan within the framework of 2021-2030 cancer control strategy, on funds administered by INSERM and La Ligue Contre le Cancer (Comités du Loiret, Loire et Cher, Eure et Loire and Morbihan, grants PM-FP/2021-407 and PM-FP/2020-244).


𝑎 Guillaume Grolleron
email: ggroller@lpnhe.in2p3.fr

Josefa Becerra González

Jonathan Biteau

Matteo Cerruti

Roger Grau

Lucas Gréaux

Talvikki Hovatta

Jean-Philippe Lenain

Elina Lindfors

Walter Max-Moerbeck

J Becerra Gonzalez

𝑏 J Biteau

𝑐 M Cerruti

𝑑 R Grau

𝑒 L Gréaux

𝑓 𝑐 T Hovatta

J.-P Lenain

𝑎 E Lindfors

𝑓 W Max-Moerbeck

𝑔 D Miceli

A Moralejo

𝑒 K Nilsson

𝑓 E Pueschel

𝑘 A Sarkar 𝑏

𝑖, 𝑗 for S Suutarinen

Variability studies of active galactic nuclei from the long-term monitoring program with the Cherenkov Telescope Array

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Blazars are active galactic nuclei (AGN) with a relativistic jet oriented toward the observer. This jet is composed of accelerated particles which can display emission over the entire electromagnetic spectrum. Spectral variability has been observed on short-and long-time scales in AGN, with a power spectral density (PSD) that can show a break at frequencies below the well-known red-noise process. This break frequency in the PSD has been observed in X-rays to scale with the accretion regime and the mass of the central black hole. It is expected that a break could also be seen in the very-high-energy gamma rays, but constraining the shape of the PSD in these wavelengths has not been possible with the current instruments. The Cherenkov Telescope Array (CTA) will be more sensitive by a factor of five to ten depending on energy than the current generation of imaging atmospheric Cherenkov telescopes, therefore it will be possible with CTA to reconstruct the PSD with a high accuracy, bringing new information about AGN variability. In this work, we focus on the AGN long-term monitoring program planned with CTA. The program is proposed to begin with early-start observing campaigns with CTA precursors. This would allow us to probe longer time scales on the AGN PSD.

38th International Cosmic Ray Conference (ICRC2023) 26 July -3 August, 2023 Nagoya, Japan

Introduction

The Cherenkov Telescope Array (CTA) will be the next generation of ground-based imaging atmospheric Cherenkov telescopes (IACT). It will improve the gamma-ray sensitivity above 20 GeV by a factor of five to ten depending on energy compared to current IACT. Therefore, CTA will bring a new outlook on the Universe in particular on its non-thermal emission.

Active Galactic Nuclei (AGN) are astrophysical sources that radiate over the entire electromagnetic spectrum. This emission comes from both thermal and non-thermal processes. One of the most important source of non-thermal emission comes from the relativistic jets, which are located near the central black hole and which can accelerate particles to the relativistic regime. In some AGN, the jet is oriented toward the observer. These objects, called blazars, show temporal variability in both their flux and their spectral properties. AGN variability can be observed on different time scales : AGN flares are responsible to the shortest scale (from minutes to hours) when the long-term variability can be observed on the year timescale [START_REF] Hovatta | Relativistic Jets of Blazars[END_REF]. The short scale AGN behavior can be described by leptonic or hadronic emission, as discussed in a study of AGN flares presented in these proceedings [START_REF] Cerruti | for the CTA Consortium, Bright blazar flares with CTA[END_REF]. Measurement of the AGN long-term behavior makes it possible to reconstruct the power spectral density (PSD) that might show a break from pink noise at low frequencies to red noise at high frequencies [START_REF] Rieger | Gamma-Ray Astrophysics in the Time Domain[END_REF]. This break frequency in the PSD has been observed in X-rays [START_REF] Edelson | A Cutoff in the X-Ray Fluctuation Power Density Spectrum of the Seyfert 1 Galaxy NGC 3516[END_REF] to scale with the ratio of the mass of the central black hole and the accretion rate [START_REF] Uttley | Measuring the broad-band power spectra of active galactic nuclei with RXTE[END_REF][START_REF] Mchardy | Active galactic nuclei as scaled-up Galactic black holes[END_REF]. Since the characterization of the PSD break at very high energies, even the brightest AGN is challenging with current IACT, CTA may be able to reconstruct this break for several of sources.

This proceeding paper focuses on the AGN long-term monitoring program of CTA, which falls within the AGN CTA Key Science Project (KSP) [START_REF] Zech | Ksp: Active galactic nuclei[END_REF]. This program consists of the monitoring of 15 selected AGN to be observed over 10 years at least once a week during their visibility period. The goal of this monitoring is to reconstruct in an unbiased way the flux distribution of AGN and the fraction of the time spent above a given flux level (the duty cycle of jetted AGN). Ultimately we also want to accurately estimate the PSD of these sources. This proceeding paper is a follow up of [START_REF] Cangemi | Probing AGN variability with the Cherenkov Telescope Array[END_REF], where we started to study the prospects of CTA through simulating AGN flares. Here the focus is on a long-term study, using simulated observations of selected AGN over 10 years. From these simulations we evaluate our ability to reconstruct the flux distribution and PSDs using CTA observations. Here we start investigating which sources are the best candidates for high-accuracy PSD reconstruction. In the following sections, we present the procedure to carry this study out. First, the injected AGN temporal and spectral models are generated, then the CTA observations are simulated, and finally a spectral model is fitted to the simulated data to reconstruct the light curve and the PSD.

Creation of AGN models from time series generation

AGN light curves can be created with a time series generator [START_REF] Emmanoulopoulos | Generating artificial light curves: revisited and updated[END_REF]. In X-rays the observed PSD of AGN can follow a pink-to red-noise PSD with a log-normal distribution of the flux [START_REF] Uttley | Non-linear X-ray variability in X-ray binaries and active galaxies[END_REF]. At this stage only a pink noise without break is simulated. We use such a model to study the capabilities of CTA to characterize PSDs. The normalization of the PSD depends on the fractional variability 𝐹 var parameter which is defined in [START_REF] Vaughan | On characterizing the variability properties of X-ray light curves from active galaxies[END_REF], the normalization is done through the equation as follows :

𝐹 var = ∫ 1 2 𝛿𝑡 1 𝑇 𝑃𝑆𝐷 (𝜈)𝑑𝜈 (1)
where 𝑇 = 10 years is the monitoring campaign duration and 𝛿𝑡 = 30 min is the sampling time used fro the generation of the time series. The cadence and the visibility presented in section 3 are considered thereafter. For this work, we use 𝐹 var = 2 to generate a log-normal time series of median equal to one, 𝜖, following the procedure described in [START_REF] Emmanoulopoulos | Generating artificial light curves: revisited and updated[END_REF]. This value has been set to produce a flux distribution in agreement with those reconstructed for Mrk 421 and Mrk 501, which are the only sources for which it is possible with current IACT [START_REF] Gréaux | STeVECat, the Spectral TeV Extragalactic Catalog[END_REF]. We model the time-dependent evolution of the AGN gamma-ray spectra Φ(𝐸, 𝑡) by assuming as spectral shape a log-parabola with exponential cut-offs, yielding:

Φ 𝑧 (𝐸, 𝑡) = Φ(𝑡) 𝐸 𝐸 0 -Γ(𝑡 ) -𝛽𝑙𝑛 𝐸 𝐸 0 𝑒 -𝐸 𝐸 cut 𝑒 -𝜏 𝛾𝛾 (𝐸,𝑧) (2)
where 𝐸 0 is the reference energy, Φ(𝑡) is the differential flux at the reference energy, Γ(𝑡) the spectral index, 𝛽 the log-parabola curvature and 𝐸 cut the cutoff energy. The last factor in Equation 2 describes the absorption of VHE photons in the extragalactic background light (EBL) with the optical depth 𝜏(𝐸, 𝑧) depending on the source's redshift, 𝑧, taken from the work of [START_REF] Domínguez | Extragalactic background light inferred from AEGIS galaxy-SED-type fractions[END_REF] and the gamma-ray energy. In this model 𝛽 and 𝐸 cut are time-independent. Their value, as well as these of the median flux Φ 0 and the spectral index Γ 0 at 𝐸 0 , are obtained for the AGN within the CTA KSP by joint modelling of archival spectra from Fermi-LAT [START_REF] Abdollahi | Incremental Fermi Large Area Telescope Fourth Source Catalog[END_REF] and IACT [START_REF] Gréaux | STeVECat, the Spectral TeV Extragalactic Catalog[END_REF]. Both the flux Φ(𝑡) and index Γ(𝑡) are time-dependent, following a harder-when-brighter behavior [START_REF] Collaboration | VHE 𝛾-ray emission of PKS 2155-304: spectral and temporal variability[END_REF][START_REF] Kapanadze | The X-ray behaviour of the high-energy peaked BL Lacertae source PKS 2155-304 in the 0.3-10 keV band[END_REF] :

Φ(𝑡) = Φ 0 [𝜖 (𝑡)] 𝜒 Γ(𝑡) = Γ 0 -𝑏 𝜒 ln 𝜖 (𝑡) (3)
The spectral variability is defined in such a way to have a constant flux at an energy 𝐸 pivot = 10 GeV [START_REF] Biteau | A window on stochastic processes and gamma-ray cosmology through spectral and temporal studies of AGN observed with H[END_REF]. To fulfill this condition, we set 𝑏 -1 = ln 𝐸 pivot 𝐸 0 and 𝜒 -1 = 1 + 𝑏 ln 𝐸 cut 𝐸 0 . We can see on Figure 1 that the spectra generated for the prototypical blazar BL Lac. The amplitude of the spectral variability, which is described by 𝐹 var and is injected at 𝐸 cut through the equation 3, is clearly visible. We can also see that the spectra are varying around the pivot at 𝐸 pivot where the amplitude of the spectral variability has been set to zero. We will use this model as an input to simulate observations with CTA.

Simulation of CTA observations and light curve reconstruction

To simulate CTA observations we are using the CtaAgnVar1 pipeline which has been designed to simulate and analyse AGN observations with CTA. This package is based on Gammapy, the analysis high-level pipeline for CTA [START_REF] Deil | Gammapy -A prototype for the CTA science tools[END_REF]. Within the CtaAgnVar pipeline, an observation sequence is created, from an input AGN time-dependent spectral model, taking into account CTA observational constraints and the visibility of the source throughout the night. The tracking of sources allows us to follow the zenith angle and to dynamically select the instrument response functions (IRF) [START_REF] Observatory | CTAO Instrument Response Functions -prod5 version v0[END_REF]. Thus, it is possible to simulate gamma-like events in a realistic manner and then fit an analytical spectral model. Furthermore, the long-term monitoring program is planned to follow the selected AGNs with a cadence of one week with 30 min of integration time. However, as most sources are not circum-polar, they are only visible between 10 to 40 weeks per year depending on their declination. The reconstruction of the light curve is done by fitting a set of spectral models to the simulation. The simplest hypothesis is a power law described by setting 𝛽 = 0 and 𝐸 𝑐𝑢𝑡 → ∞ in Equation 2. A more complex version is a log-parabola (with 𝛽 ≠ 0) or a power law with exponential cutoff (a finite 𝐸 𝑐𝑢𝑡 ), and the most complex is the combination of both (log-parabola and power law with exp. cutoff). For each time bin, a more complex model is preferred if the simplest one can be rejected at a 3𝜎 level based on a likelihood ratio test. As the log-parabola and the power law with exponential cutoff have the same complexity, the model providing the highest likelihood value is preferred. Finally, to quantify the agreement between the data and the model, a goodness-of-fit estimate has been developed based on the likelihood ratio test between the data and the model. The test statistics (TS) derived from the likelihood ratio follows a 𝜒2 distribution with TS= 0 describing a perfect fit. We defined therefore a 𝑝 𝑣𝑎𝑙𝑢𝑒 to identify possible outliers.

The simulations have been done for each source listed in the CTA AGN KSP [START_REF] Zech | Ksp: Active galactic nuclei[END_REF]. However, here we focus on three bright sources: Mrk 501 (z = 0.033), BL Lac (z = 0.069), and PKS 1510-089 (z = 0.36). The distributions of detection significance in weekly 30-min exposures are shown in Figure 2. As an example, the reconstructed light curve, for the most promising source with CTA, BL Lac, above 50 GeV is presented in Figure 3 together with the median spectrum computed over the 10 years of data. The first important result is that there is very good agreement between the injected and the reconstructed light curve with residuals around 1%.

PSD reconstruction

Having obtained the light curves, we can proceed to estimate, through periodogram, the PSDs for each source following [START_REF] Emmanoulopoulos | Generating artificial light curves: revisited and updated[END_REF]. PSDs reconstruction for PKS 1510-089 and BL Lac are presented in Figure 4. The computation with error propagation is performed following the procedure of [START_REF] Biteau | A window on stochastic processes and gamma-ray cosmology through spectral and temporal studies of AGN observed with H[END_REF]. The red line is not a fit of the PSD, but a pink noise spectrum (index of 1) with the Poisson plateau TeV -1 .s -1 .cm -2 , Γ = 2.37 ± 0.01, 𝐸 -1 cut = 1.08 ± 0.06 TeV -1 and 𝛽 statistically consistent with 0. Black points are the flux points reconstructed with current IACT obtained from [START_REF] Gréaux | STeVECat, the Spectral TeV Extragalactic Catalog[END_REF] at high frequencies. The plateau is clearly visible in the PSD of PKS 1510-089 (Figure 4, right). The normalization of the red line is done according to Equation 1. It shows that the source PSDs can be established with the CTA observation scheme, as the red line agrees well with the estimated PSDs.

CTA perspective for long-term AGN studies

In this contribution, we studied CTA's capabilities for a long-term monitoring of AGN. Figure 2 represents the significance distribution for each obervation for several interesting sources in CTA AGN KSP. The detection significance for BL Lac or PKS 1510-089 are comparable to the actual significance distribution of Mrk 501 obtained with current IACT.

In order to quantify the capabilities of CTA, we have simulated data of various AGN, applied the foreseen observation cadence, and reconstructed the PSDs. While it is possible to reconstruct the input PSD for bright sources, such as BL Lac and PKS 1510-089, the foreseen monitoring schedule (30min per source per week) might need to be optimized for fainter sources. Moreover we can clain that the duty cycle of jetted AGN might be well constrained by CTA. We have just used a simple pink-noise input spectrum here, but will present results employing more complicated spectra in a forthcoming publication. Moreover, this effort will be supported by the observing campaigns with CTA precursors. Nonetheless, the results presented here are encouraging and show that CTA will be able to probe deeply the VHE AGN long-term behavior.

Figure 1 :

 1 Figure 1: Spectra generated for BL Lac with the time generation algorithm. The color code marks the time evolution for 10 years from blue to red. Obviously, there is no variability at 𝐸 pivot as designed.

Figure 2 :

 2 Figure 2: Left: Detection significance distribution for 3 promising sources. The vertical dotted line is the median value associated to each source. Right: Comparison with the one of Mrk 501 observed by MAGIC in a low-state campaign [21]. To allow unbiased comparison between MAGIC data and CTA simulations, the pink distribution is plotted selecting observations with a flux level comparable to the flux distribution observed by MAGIC.

Figure 3 :

 3 Figure 3: Left: Reconstructed light curve and residuals computed between the simulated and reconstructed data for BL Lac for the full simulation of 10 years. Gray points are the injected values and red points are the reconstructed ones. Right: Reconstructed median spectrum for BL Lac for the 10 years of data. Spectral parameters are: Φ 0 = 7.46 ± 0.07 × 10 -10TeV -1 .s -1 .cm -2 , Γ = 2.37 ± 0.01, 𝐸 -1 cut = 1.08 ± 0.06 TeV -1 and 𝛽 statistically consistent with 0. Black points are the flux points reconstructed with current IACT obtained from[START_REF] Gréaux | STeVECat, the Spectral TeV Extragalactic Catalog[END_REF]

Figure 4 :

 4 Figure 4: PSDs estimate (blue points) of the simulated data for BL Lac (left) and PKS 1510-089 (right).The red line shows the injection PSD (module floor level at high frequencies) used to simulate the input data.

1https://gitlab.cta-observatory.org/guillaume.grolleron/ctaagnvar

Acknowledgments

This work was conducted in the context of the CTA Consortium. We gratefully acknowledge financial support from the agencies and organizations listed here: https://www.cta-observatory. org/consortium_acknowledgments/.


Vincent Stin
email: vincent.stin@espci.fr

Ramiro Godoy-Diana

Xavier Bonnet

Anthony Herrel

METHODS &

Keywords: Snake, Swimming, Anguilliform, Hydrodynamics, Volumetric PIV, Vortex

We describe a method for measuring the 3D vortical structures produced by an anguilliform swimmer using volumetric velocimetry. The wake of freely swimming dice snakes (Natrix tessellata) was quantified, revealing the creation of multiple vortices along the body of the snake due to its undulation. The 3D structure of the vortices generally consisted of paired vortex tubes, some of which were linked together to form a hairpin structure. The observations match predictions from computational fluid dynamic studies of other anguilliform swimmers. Quantitative measurements allowed us to study vortex circulation and size, and global kinetic energy of the flow, which varied with swimming speed, vortex topology, and individual characteristics. Our findings provide a baseline for comparing wake structures of snakes with different morphologies and ecologies and investigating the energetic efficiency of anguilliform swimming.

INTRODUCTION

Undulatory swimming kinematics are usually classified in four main modes involving different proportions of the body and/or caudal fin [START_REF] Lindsey | Form, function and locomotory habits in fish[END_REF][START_REF] Sfakiotakis | Review of fish swimming modes for aquatic locomotion[END_REF]). Among these, anguilliform swimming describes the motion of elongated animals, where the kinematics consist of an undulation that increases in amplitude along the body. Since Lighthill's analytical large amplitude elongated body theory [START_REF] Lighthill | Large-amplitude elongated-body theory of fish locomotion[END_REF], numerical investigations have been conducted [START_REF] Kern | Simulations of optimized anguilliform swimming[END_REF]Borazjani andSotiropoulos 2009, 2010;[START_REF] Nangia | Optimal specific wavelength for maximum thrust production in undulatory propulsion[END_REF]Battista 2020a,b;[START_REF] Khalid | Flow transitions and mapping for undulating swimmers[END_REF][START_REF] Ogunka | The ground effect in anguilliform swimming[END_REF][START_REF] Khalid | Why do anguilliform swimmers perform undulation with wavelengths shorter than their bodylengths?[END_REF]) in order to obtain a estimation of quantities such as forces, pressure fields, and swimming efficiency for an idealized swimmer. Anguilliform swimming hydrodynamics have been experimentally studied on the eels Anguilla anguilla [START_REF] Muller | How the body contributes to the wake in undulatory fish swimming: flow fields of a swimming eel (Anguilla anguilla)[END_REF]) and Anguilla rostrata [START_REF] Tytell | The hydrodynamics of eel swimming II. Effect of swimming speed[END_REF][START_REF] Tytell | The hydrodynamics of eel swimming: I. Wake structure[END_REF], the lamprey Petromyzon marinus [START_REF] Gemmell | How the bending kinematics of swimming lampreys build negative pressure fields for suction thrust[END_REF][START_REF] Du Clos | Thrust generation during steady swimming and acceleration from rest in anguilliform swimmers[END_REF][START_REF] Lehn | Volumetric analysis of lamprey hydrodynamics using synthetic aperture particle image velocimetry[END_REF], and the catfish Plotosus lineatus [START_REF] Tack | Anguilliform Locomotion across a Natural Range of Swimming Speeds[END_REF], using particle image velocimetry (PIV) in 2D.

The elongated and limbless morphology of snakes is associated with over ten different gaits used in various environments [START_REF] Jayne | What defines different modes of snake locomotion?[END_REF], with swimming kinematics first studied quantitatively by [START_REF] Taylor | Analysis of the swimming of long and narrow animals[END_REF] inspired by the work of [START_REF] Gray | Studies in animal locomotion: I. The movement of fish with special reference to the eel[END_REF]. Snakes' great diversity in size and shape makes them ideal to study anguilliform locomotion, and their numerous independent radiations into aquatic environments provide insight into the adaptation of elongated limbless tetrapods to this environment. However, apart from prey capture studies [START_REF] Van Wassenbergh | Hydrodynamic constraints on prey-capture performance in forward-striking snakes[END_REF][START_REF] Segall | Hydrodynamics of frontal striking in aquatic snakes: drag, added mass, and the possible consequences for prey capture success[END_REF], a hydrodynamic analysis based on live swimming snake data has yet to be attempted. Although 2D PIV helps to understand aquatic locomotion dynamics, 3D velocimetry is necessary to fully describe the vortex dynamics produced by a moving finite body [START_REF] Tytell | Escaping Flatland: Three-Dimensional Kinematics and Hydrodynamics of Median Fins in Fishes[END_REF], especially in out-of-plane movements as in swimming snakes. The main methods are synthetic aperture PIV (e.g., [START_REF] Lehn | Volumetric analysis of lamprey hydrodynamics using synthetic aperture particle image velocimetry[END_REF][START_REF] Mendelson | Jumping archer fish exhibit multiple modes of fin-fin interaction[END_REF], tomographic PIV [START_REF] Skipper | Characterization of hop-and-sink daphniid locomotion[END_REF], and digital defocusing PTV (e.g., [START_REF] Flammang | Volumetric imaging of fish locomotion[END_REF][START_REF] Bartol | Volumetric flow imaging reveals the importance of vortex ring formation in squid swimming tailfirst and arms-first[END_REF]. These methods instantly capture the entire 3D unsteady flow field, but flow field measurements of swimming animals are scarce. Aside from unpublished work on lampreys [START_REF] Lehn | Volumetric analysis of lamprey hydrodynamics using synthetic aperture particle image velocimetry[END_REF], the 3D wake of an anguilliform swimmer has not been characterized experimentally. We propose a method to examine the 3D flow induced by an anguilliform swimmer using a free-swimming dice snake (Natrix tessellata) as our model. Using DDPTV, we measured the wake and investigated the vortical structure. We compared results to previous computational fluid dynamics and 2D experimental PIV studies.

MATERIAL AND METHODS

Animals and swimming behaviour

Swimming kinematics and PIV data were obtained from three (one male and two females) adult captive dice snakes (Natrix tessellata). Dice snakes are semi-aquatic snakes and were chosen for their relatively small size. A first set of trials was made to characterize the swimming kinematics of Natrix tessellata as no data are available in the literature and because the kinematic parameters of swimming snakes are generally scarce. The snakes were freely swimming in a plexiglass tank (length = 400 cm, width = 40 cm, height = 30 cm, water height = 22 cm). Animals were filmed in ventral view (Phantom MIRO, 1000 fps). Kinematic parameters were extracted by digitizing the snake's midline using an custom-written MATLAB routine. We measured the following anatomical and kinematic parameters: the total length, from the tip of the snout to the tip of the tail (L), maximal body diameter (D), body mass (m), forward swimming speed (U ), tailbeat frequency (f ), peak-to-peak amplitude of the tail (A), and the body wavelength (λ). The Reynolds and Strouhal numbers were also computed as:

Re = U D ν
(1)

St = f A U (2)
where ν is the water kinematic viscosity. The general swimming kinematic parameters are summarized in Fig. 1.B and permit a broad estimation of the swimming behaviour of the test animals under laboratory conditions. Although conducted separately from the PIV experiments, these trials have the potential to be performed concurrently in a separate study, with appropriate laser filtering in place. An example of a swimming trial video is accessible in the supplementary materials.

Experimental PIV setup

In order to quantify the flow around a swimming snake we used a volumetric three-component DDPTV setup (V3V-9000-CS system, TSI ; Fig. 1.C). The experiments were conducted in a waterfilled plexiglass tank (length = 210 cm, width = 18 cm, height = 30 cm, water height = 22 cm). A narrower tank was chosen in order to increase the chances of the snake to swim through the measurement volume without over-constraining the kinematics. The water was seeded with 50 µm polyamide particles (PSP-50, Dantec Dynamics) with a concentration of around 5 × 10 -2 ppp (Cambonie and Aider 2014). The tank was illuminated from above by a 14 Hz pulsed 200 mJ dual head Nd:YAG laser (Quantel Evergreen) expanded with two -25 mm cylindrical lenses. The particles were filmed by three CCD arrays (Powerview Plus 4MP-HS) with 50 mm camera lenses (Nikon AF Nikkor, aperture f/16) mounted on a 170 mm equilateral triangular frame (V3V 9000-CS). The acquisition setup was piloted by the INSIGHT V3V-4G software and calibrated as explained by [START_REF] Troolin | Volumetric velocity measurements of vortex rings from inclined exits[END_REF]. The resulting measurement volume is the intersection between the fields of view of the cameras and the laser cone. Its dimension is approximately 18 × 18 × 12 cm and it is centered on the water tank in width (Z-axis) and height (Y-axis).

Recording procedure

During each trial the snake was placed at one side of the tank and allowed to swim. When put into the water snakes typically swam to the other side of the tank along the X direction, passing through the measurement volume while being recorded by the cameras at a frequency of 14 Hz. Preliminary trials were conducted to estimate the best ∆T between two consecutive PIV frames so that displacement of the seeding particles does not exceeded six pixels. The resulting ∆T ranged between 1 ms and 2.5 ms. As the width of the measurement volume was smaller than the width of the water tank, only the trials were the snake entirely swam through it were retained. The PIV swimming trials were performed on six different days over the course of several weeks. During the first two recording days, females were gravid which impacted their body mass and diameter. The measurements in Table 1 were made after the eggs were laid. A total of 112 sequences of free swimming Natrix tessellata with at least one tailbeat in the measurement volume were recorded. The number of sequences were respectively 31, 52 and 29 for Natrix M, Natrix F1 and Natrix F2.

Data processing

Image masking

In order to avoid artefacts induced by the presence of the snake's body on the raw images during the PIV analysis, the snake was removed from the images with the help of a custom-written masking routine in MATLAB using the Image Processing Toolbox.

PIV Processing

The images were then processed using the INSIGHT V3V-4G software in order to obtain the vector fields. The main steps of the data processing are: (1) 2D particle identification in the images of each camera using a Gaussian fitting (Levenberg-Marquardt algorithm). ( 2) 3D particle identification, the correspondence of a particle in one image to the same particle in the other two images is carried in agreement with the calibration. The fine search tolerance is set to 0.5 px and the coarse search tolerance to 1 px.

(3) The velocity is processed using the 3D displacement between two consecutive frames with a relaxation algorithm [START_REF] Pereira | Twoframe 3D particle tracking[END_REF]). ( 4) The neighbor tracking reconstruction algorithm uses the trajectories from the neighboring particles to find the probable location of a missing triplet at time t + ∆T.

(5) The randomly spaced vectors are finally interpolated using a Gaussian-weighted algorithm [START_REF] Pereira | Defocusing digital particle image velocimetry and the three-dimensional characterization of two-phase flows[END_REF], resulting in a 74 × 74 × 43 evenly spaced vector grid with a 2.5 mm resolution.

Base change

Vector fields were processed using code written in MATLAB. The swimming trials being in free swimming conditions, the movements of the snake in the measurement volume can be quite different from one sequence to another, which makes the comparison challenging. In order to tackle this problem, we changed the regular Cartesian coordinate system (XY Z) to a new rotated coordinate system (X rot Y Z rot ) for each swimming trial. For each recorded sequence the water jet created by the last tailbeat was tracked with a velocity thresholding method. The trajectory of the jet was estimated by a linear regression of the superimposed position of the jet center. A new base was then created by a rotation around the Y-axis so that the new X rot axis was co-linear to the jet's trajectory and the Z rot axis perpendicular to it.

Hydrodynamic parameters

The flow field induced by the swimming snake was monitored from its detection in the volume until it either reached the edges of the working volume or dissipated. The unsteadiness was quantified by computing the Q-criterion of the vector field. Q is defined as the second invariant of the velocity gradient tensor [START_REF] Hunt | Eddies, streams, and convergence zones in turbulent flows[END_REF][START_REF] Jeong | On the identification of a vortex[END_REF]:

Q = 1 2 (∥ Ω ∥ 2 -∥ S ∥ 2 ) ( 3
)
where Ω is the antisymmetric part and S the symmetric part of the tensor. The positive values of Q indicate the parts of the flow field where vorticity dominates over the viscous stress. The strength of the vortices was estimated by computing the circulation of the vortex cores:

Γ = C ⃗ v • d ⃗ l = S ⃗ ω • ⃗ ndS (4)
where ω is the vorticity vector and n the normal vector. The surfaces S of the vortex cores are identified by using the Γ 2 criterion [START_REF] Graftieaux | Combining PIV, POD and vortex identification algorithms for the study of unsteady turbulent swirling flows[END_REF] in the rotated X rot Y planes co-linear to the jet trajectory so that a vortex core is always present in the planes.

The vortex size d is defined by the distance between the center of the positive and negative vortex cores:

d = (x + -x -) 2 + (y + -y -) 2 (5)
where (x + , y + ) and (x -, y -) are the coordinates of the centers of the positive and negatives vortex cores in the X rot Y planes. In order to compare the results between the different individuals and sequences, the hydrodynamic parameters are non-dimensionalized by the mean forward swimming speed U of the snake during the measurement period and the snake diameter D :

Γ * = Γ U D (6) d * = d D (7)
The kinetic energy of the wake in the whole measurement volume is defined as :

E = ρ 2 V ∥ u ∥ 2 dV ( 8
)
where ρ is the water density, V the measurement volume and ∥ u ∥ the velocity magnitude. To quantify differences in the hydrodynamic parameters across individual sequences we focused on trials where at least one tail beat vortex was clearly visible and where measurements were possible (N M = 12; N F 1 = 13; N F 2 = 14).

RESULTS & DISCUSSION Vortex characterisation

Although the snakes swam with the characteristic lateral undulation of an anguilliform swimmer in the frontal plane (XZ), there were also noticeable, albeit smaller, oscillations in the sagittal plane (XY). The tail was often positioned below the body and the tailbeats were downward facing. The mean swimming speed in body lengths per second across the swimming trials was, respectively, U M = 0.51 Ls -1 , U F 1 = 0.37 Ls -1 and U F 2 = 0.34 Ls -1 (Fig. 1.D).

During a typical swimming sequence, vortices were shed after each change of direction of the swimmer's body during oscillatory movements. The displacement of the body created two vortex tubes of opposite rotation (Fig. 2.A,B) on its top and bottom, which in some cases became connected. The shed vortices were tilted at an angle of approximately 30°to 40°to the swimming direction (Fig. 2.C). Their trajectory was mostly linear with a tendency to curl. Body vortices with similar topologies and dynamics were also observed.

The circulation of the vortex cores rapidly increased in magnitude just after the tail beat and reached a maximum before decreasing (Fig. 2.D). The positive and the negative vortex cores displayed symmetrical variations during the observed time period. During its lifetime, the vortex size gradually increased (Fig. 2.E) and changed from a flattened into a more circular shape. The kinetic energy of the wake increased starting from the moment when the swimmer enters the measurement volume. After the final tailbeat of the swimmer, the perturbations dissipate and the kinetic energy decreases (Fig. 2.F).

Comparison of parameters

Across the swimming sequences, the dimensional maximal circulation ranged from around 18 to 51 cm 2 s -1 , 25 to 70 cm 2 s -1 and 22 to 88 cm 2 s -1 for Natrix M, Natrix F1, and Natrix F2, respectively. The non-dimensional circulation Γ * ranged approximately between 0.5 and 1.5. The maximum vortex size was between approximately 2.5 and 5.3 cm, 2.7 and 5.2 cm, and 2.6 and 5.4 cm for Natrix M, Natrix F1, and Natrix F2 respectively. A Spearman correlation (R=0.8, P<0.05) showed that there was an increasing monotonic trend between the maximal vortex size and the maximal circulation. As the maximal circulation is usually measured near the beginning of the vortex lifetime, this means that a high initial circulation is likely to yield a bigger vortex.

A Spearman correlation (R=0.7, P<0.05) between the maximal kinetic energy in the whole measurement volume during the entire swimming sequence and the maximal circulation of the tailbeat vortex showed an increasing monotonic trend. This means that there is a tight scaling between the local hydrodynamic parameters of the vortices (i.e., the circulation) and the global parameters (i.e., the kinetic energy of the fluid in the whole measurement field).

Discussion

As the wake is the footprint of a swimmer, it is crucial to fully capture it in order to maximise the understanding of the fluid-structure interactions between the body and the water. The 2D PIV experiments on eels by [START_REF] Muller | How the body contributes to the wake in undulatory fish swimming: flow fields of a swimming eel (Anguilla anguilla)[END_REF] in the frontal plane have led to the prediction that the 3D vortical structure of an anguilliform swimmer should resemble a vortex ring. [START_REF] Tytell | The hydrodynamics of eel swimming: I. Wake structure[END_REF] were skeptical about the true roundness of the vortex and the computational work of [START_REF] Borazjani | Numerical investigation of the hydrodynamics of anguilliform swimming in the transitional and inertial flow regimes[END_REF] suggested that the wakes were in fact vortex loops stretched in the streamwise direction or "hairpin-like" vortices. The traditional plane to measure the swimmer's wake in the existing literature is the frontal plane (XZ) without looking at the other ones (Fig. 3.A,B,C). Our experimental results in the frontal plane are in accordance with the 2D vortex behaviour in the existing anguilliform swimming studies (Fig. 3.D) with the creation of two counter-rotating vortex cores during a tailbeat cycle. The 3D velocimetry method presented here reveals that these vortex cores observable in the frontal plane are the links between the vortex tubes stretched in the swimming direction (Fig. 3.E). These links are created after each change of direction of the tail.

The observed three-dimensional vortex topologies were either hairpin-like, tubular or turbulent without any clear cohesion. We observed hairpin-like and tubular structures in most of the analysed swimming sequences. A lack of cohesion was observed essentially only when the snake kinematics involved strongly out-ofplane motions, changes of direction, or acceleration/deceleration sequences. These different vortex topologies have been predicted by CFD studies [START_REF] Kern | Simulations of optimized anguilliform swimming[END_REF][START_REF] Borazjani | Numerical investigation of the hydrodynamics of anguilliform swimming in the transitional and inertial flow regimes[END_REF][START_REF] Khalid | Why do anguilliform swimmers perform undulation with wavelengths shorter than their bodylengths?[END_REF], and associated to different kinematic parameters such as the tailbeat amplitude, frequency, and body wavelength.

As the DDPTV enables the observation of the full wake structure, our measurements revealed the different stages of the formation of the hairpin vortices. The shedding of two vortex tubes of opposite rotation immediately after a tailbeat (Fig. 2.G) creates a bridge where the vortex tubes interact (Fig. 2.H), leading to their connection and detachment before the bridge (Fig. 2.I). These steps are similar to events in coherent packets of hairpin vortices in channel flow [START_REF] Zhou | Mechanisms for generating coherent packets of hairpin vortices in channel flow[END_REF]. The rounding of the vortex structure after the bridging step is also reminiscent of the Ω-shaped vortex created by self-induction of the hairpin vortex. Despite the relatively slow frame rate, our 3D measurement allowed us to observe the bridging step. The shedding of vortices was also witnessed more anteriorly on the body. The production of body vortices has been described by [START_REF] Gemmell | How the bending kinematics of swimming lampreys build negative pressure fields for suction thrust[END_REF] for lampreys where a particular attention was given to the bending of the body. In addition to transverse motion, the body has an angular velocity in the transverse plane (YZ).

The body surface rotation helps to rotate the fluid and eases the shedding of adjacent vortices. The production of counter-rotating vortex tubes along the body has been suggested to generate additional thrust [START_REF] Fu | Transient force augmentation due to counter-rotating vortex ring pairs[END_REF] during swimming. In the case of gliding snakes, a CFD study showed that the body-induced vortices produced lift along the body [START_REF] Gong | Computational analysis of vortex dynamics and aerodynamic performance in flying-snake-like gliding flight with horizontal undulation[END_REF]). The body vortices may therefore have a double function of propelling and stabilising the swimmer. The maximal circulation measured for the three swimmers was similar to the ones measured in smaller anguilliform swimmers: ∼ 20 cm eels [START_REF] Tytell | The hydrodynamics of eel swimming: I. Wake structure[END_REF]) and ∼ 10 cm lampreys [START_REF] Gemmell | How the bending kinematics of swimming lampreys build negative pressure fields for suction thrust[END_REF]. The circulation of small hairpin vortices decreases exponentially after attaining a maximum because of the viscous interactions of the counter-rotating vortex tubes [START_REF] Sabatino | Characterizing the formation and regeneration of hairpin vortices in a laminar boundary layer[END_REF]. On the other hand, the circulation of laminar vortex rings slowly decreases after the maximum [START_REF] Rosenfeld | Circulation and formation number of laminar vortex rings[END_REF]. In our results, depending on the swimming trial, the circulation either decreases or stays relatively constant. The decreasing vortex circulation is probably due to the viscous interactions with either the counter rotating vortex tubes or surrounding vortices created by the other tail beats. The relatively constant or slowly decreasing circulation, on the other hand, behaves like a laminar vortex ring.

Studies on vortex rings show that the diameter of the vortex rapidly increases, attains a maximum, and then asymptotically decreases to a constant value [START_REF] Arakeri | Vortex ring formation at the open end of a shock tube: a particle image velocimetry study[END_REF]). The distance between the vortex cores observed here follows the same pattern as they either increase or stay relatively constant. This means that the vortices with an increasing distance are in the first part of their lifetime and the ones with a relatively constant distance are likely to be near the end of their lifetime. Our data further suggests a relationship between the circulation and vortex size.

As the experiments were done in a still-water tank, the displacement of the water was essentially caused by the swimming snakes.

The kinetic energy of the water measured is therefore the energy that was transferred by the snake minus the dissipation. Computational studies found faster swimming kinematics result in stronger vortices [START_REF] Kern | Simulations of optimized anguilliform swimming[END_REF][START_REF] Bhalla | A Forced Damped Oscillation Framework for Undulatory Swimming Provides New Insights into How Propulsion Arises in Active and Passive Swimming[END_REF]. Gravid snakes had larger body diameter and higher mass but proportionally less locomotor muscle due to eggs which may impact locomotor performance [START_REF] Seigel | Reduction in locomotor ability as a cost of reproduction in gravid snakes[END_REF]. Gravid snakes shed vortices with high circulation and size, producing the highest overall kinetic energy, suggesting they adjust swimming effort at a higher cost. Complementary results from other species may provide an explanation for the observed phenomenon of gravid snakes exhibiting less swimming activity compared to non-gravid individuals [START_REF] Aubret | Swimming and pregnancy in Tiger snakes, Notechis scutatus[END_REF].

Volumetric measurements of the wake can be used as tools for non-invasive estimations of the swimming performance [START_REF] Li | Volumetric measurements of wake impulse and kinetic energy for evaluating swimming performance[END_REF]. Coupled with 3D kinematic data, the results of this method could be used for interspecific comparisons of the swimming performance of snakes. This would lead to a better understanding of the differences of swimming efficiency in snakes with different lifestyles, shapes, and ecologies.

Fig. 1 .

 1 Fig. 1. Experimental setup and snake kinematics. (A) Superimposed images of a snake at two consecutive maxima of tail deflection. The peak-to-peak amplitude A, tailbeat frequency f and the forward velocity U are measured by finding the peaks in the lateral excursion of the tail-tip over time for at least one full tailbeat cycle. The body wave length λ is measured when two adjacent crests are clearly visible on the snake body. (B) Anatomical & kinematic parameters gathered during the swimming behaviour experiment with the asterisk meaning that the parameter was divided by the snake length. Number of swimming trials: N M = 17, N F 1 = 11, N F 2 = 11. Results are presented as means±s.d. (C) Volumetric imaging setup. a: Water tank, b: Nd:YAG laser (Quantel Evergreen), c: Measurement volume, d: Camera frame. (D) Boxplot of mean voluntary forward velocity of the swimming snakes during the PIV swimming trials (N M = 31, N F 1 = 52, N F 2 = 29).

Fig. 2 .

 2 Fig. 2. Results. (A,B,C) Trailing edge vortices in the measurement volume during a swimming trial highlighted with Q-criterion isosurfaces (Q = 80) colored by ωx. (A,B) Side view of two consecutive trailing edge vortices (∆T = 286 ms, half of a tailbeat cycle) with the corresponding snake image (Natrix F1) captured by the PIV cameras. (C) Top view of the superimposed isosurfaces shown on (A) and (B) at their respective times showcasing the angle between two consecutive vortices. (D,E,F) Measured hydrodynamic parameters during the swimming trial. (D) Positive (red) and negative (blue) vortex core non-dimensionalized circulation of the trailing edge vortex (B) from the frame before its detection (t * = 0) to the time when it reaches the edge of the measurement volume. (E) Size-adjusted vortex size evolution. Values are the mean values in three consecutive planes centered on the vortex path. Error bars are s.d. (F) Evolution of the kinetic energy in the whole measurement volume during the swimming trial. The two vertical lines correspond to the beginning and the end of t * . (G,H,I) 3D view showcasing the main steps of the behaviour of a hairpin-like vortex with the example of the trailing edge vortex in (B). (G) Two vortex tubes of opposite rotation (t * = 71 ms). (H) Bridging of the two vortex tubes (t * = 285 ms).(I) Separation of the tubes before the bridge and the shape of the structure is getting rounder (t * = 500 ms). The mean vorticity at the bridge section (ω Y = 15.08 s -1 ) is similar to the vorticity of the vortex tubes (ω X = 15.45 s -1 ). The corresponding video can be accessed in the supplementary materials.

Fig. 3 .

 3 Fig. 3. Vortex comparison. (A,B,C,D) Trailing edge vortices in the frontal plane (XZ) from different studies. The swimmers are swimming along the X axis and the visible vortex cores are the result of a full tailbeat cycle. Vector fields from (A) an eel (Tytell and Lauder 2004) (B) a lamprey (Gemmell et al. 2016) (C) an anguilliform swimmer CFD model (Kern and Koumoutsakos 2006) and (D) our results. The flow field in (D) is smoothed with a 3D Gaussian filter (σ = 1.25) to enhance the illustration of the dipolar vortex structure only (note that the quantitative computations of vortex circulation and kinetic energy were performed on the unfiltered velocimetry output). The corresponding ω Y -colored 3D plot (E) suggests that the vortex structure is more complex than it appears in a 2D plane. The grey plane is the location of the velocity field (D).

Acknowledgements

This research was supported by ANR DRAGON2 (ANR-20-CE02-0010). We thank Jake Socha and an anonymous reviewer for helpful and constructive comments on an earlier version of the manuscript.

Data availability

All data and codes are available from the authors upon request.

Competing interests

The authors declare no competing or financial interests.

Contribution

A.H., X.B. and R.G.D. conceived the study; V.S. carried out experimental work; V.S. carried out data analysis; V.S. wrote the manuscript; A.H., X.B. and R.G.D. revised the manuscript. All authors gave final approval for publication.