David E Cohen

Annalisa Allocca

Gilles Bogaert

Paola Puppo

Thibaut Jacqmin
email: thibaut.jacqmin@sorbonne-universite.fr

Towards optomechanical parametric instabilities prediction in ground-based gravitational wave detectors

come

I. INTRODUCTION

In 2015, the LIGO-Virgo collaboration [1][2][3][4] detected for the first time gravitational waves preceding a binary black hole coalescence [5], thus pioneering gravitationalwave astronomy. Today many other gravitational waves have been detected [6,7]. These detections have provided confirmation on the expected rate of binary black hole (BBH) mergers [8], a better understanding of BBHs population [8,9], a better limit to the mass of the graviton [10], a first direct evidence of a link between binary neutron star (BNS) mergers and short gamma-ray bursts [11], a higher precision in constraining the Hubble constant [START_REF] Abbott | Dark Energy Camera GW-EM[END_REF], and a better understanding of BNS mergers [11]. Since the first detections, improvements performed on ground based detectors yielded better detector sensitivities. Gravitational wave sources that are weaker or located further away can now be detected. Among the many improvements, increasing the light intensity in the interferometer arm cavities reduces the impact of the laser quantum phase noise, which is limiting the sensitivity in the high-frequency range. However, a laser power increase can trigger a nonlinear optomechanical effect [13,14], known as optomechanical parametric instability (OPI). This effect can jeopardize the interferometer stable operation.

During the Observing Run 1 (O1), LIGO experienced an OPI for the first time [15]: after a few seconds, the interferometer went out of lock, thus preventing further data acquisition. In this letter, we present the models that we used to compute the OPI gains for a power recycled gravitational wave detector. Compared to previous work [17], we provide a precise description of optical and mechanical modes, together with a study of the impact of losses and the thermal deformation of the mirrors. This will allow to perform parametric instabilities predictions, which is of major importance for future designs of ground-based detectors. These simulations are done for the power recycled Advanced Virgo interferometer (O3 configuration), and can be extended to any other configuration.

In sec. II, a short introduction to the model used to compute the parametric gains is given. In sec. III, we report on a detailed finite element analysis (FEA) of the mirrors, to compute precise mechanical modes frequencies and amplitudes, and estimate quality factors. An original method is then used to combine these FEA simulations with ring-down measurements performed on a subset of modes, in order to obtain accurate quality factors for all the modes. In section IV, different models for optical modes are compared: the analytical solution of the paraxial equation for purely spherical infinite size mirrors, as implemented in [17], and a brute force numerical simulation which includes finite size effects and arbitrary mirror surface shapes. In section V, we study the influence on optical modes of a thermal effect related to a local temperature increase of the mirror surface due to light absorption. Finally in section VI, we provide an example of parametric gains that are obtained in the Advanced Virgo O3 configuration, including optical losses and mechanical losses calculated with an unprecedented level of precision, and we investigate for the first time the effect of the mirror thermal deformation due to laser absorption.

II. OPTOMECHANICAL PARAMETRIC INSTABILITY

In an optomechanical cavity like one arm of a gravitational wave detector, photons from the optical zeroth order mode can be coherently scattered to a higher order transverse optical mode if a mechanical mode that sets a mirror surface into motion has its frequency ω m /2π equal to the frequency difference between the two optical modes (modulo the cavity free spectral range). This phenomenon can remove energy from the mechanical mode by annihilating phonons, and scattering photons from the zero order mode to the higher order transverse mode, thus damping the mirror motion [16]. Conversely, this phenomenon can add energy to the mechanical mode with the reverse process, thus exciting the mechanical motion. In that case, an instability can prevent the interferometer stable operation [5,13]. This instability has a threshold: it starts to grow as soon as the resonant excitation of the mechanical mode by the radiation pressure force overcomes mechanical losses.

In the following, we use the approach developed by Evans et al. [17] to simulate this effect. In this framework, the whole interaction between the three implied modes (two optical modes and a mechanical mode) is seen as a classical feedback system. This modular approach is well suited, since it can be adapted to many different interferometer configurations with the same analytical formulas. The parametric gain of the mechanical mode m is given by

R m = 8πQ m P M ω 2 m cλ ∞ n=0 [G n ]B 2 m,n (1)
where Q m is the quality factor of the mechanical mode m and ω m its frequency, P the arm-cavity optical power, λ the optical wavelength, M the mirror mass, c the velocity of light, G n is related to the scattered field optical gain of the n th optical mode and encapsulates the interferometer configuration. Finally, B m,n is the spatial overlap integral between the three involved modes. A mechanical mode is amplified if R m > 0 and damped if R m < 0. It becomes unstable if R m > 1.

FIG. 2. Geometry used for the FEA, including the ears, the anchors and the magnets attached on the mirror rear face. The suspension wires are just for sketching but not included in the simulation as they do not influence the modal frequencies

III. MECHANICAL SIMULATION

A. The spatial overlap parameter

The spatial overlap integral B 2 n,m is defined [18] as

B 2 m,n = M M eff E 00 ( r)E n ( r)µ m ⊥ ( r)d r ⊥ 2 |E 00 | 2 d r ⊥ |E n | 2 d r ⊥ , ( 2
)
where M is the mirror mass and m eff the effective mass of the mechanical mode. The integral is performed over the test mass surface (coating side). E 00 ( r) stands for the optical carrier amplitude and E n ( r) for a transverse optical For the FEA the multi-layer coating of the IM was replaced by one layer having the total thickness corresponding to the sum of the thicknesses of the high reflective and low reflective materials. The mechanical parameters used are the average values of this layer [20] mode amplitude labeled by the index n. As the interferometer is sensitive to the test mass displacement along the optical axis, only the vertical displacement µ m ⊥ ( r) is considered, where m is the mechanical mode index. The effective mass is related to the strain energy ρ e through the equation 1 2 M eff ω 2 m = ρ e r ⊥ , and effectively obtained with the formula

M eff = M < µ( r) 2 >= M 1 V µ( r) 2 d r, (3)
where µ( r) stands for the test mass displacement. The mechanical modes were computed by means of finite element analysis (FEA) developed for the actual input test mass (IM) of Advanced Virgo arm cavities [START_REF] Puppo | Evaluation of the quality factors of the mirrors in advanced virgo using the fem ansys model[END_REF]. We have used the program Ansys®Workbench TM . The IM model includes the high-reflectivity (HR) coating of the front face, the flats and the bevels. Moreover the ears and the anchors attached by silicate bonding technique are included (see figure 2). In the FEA, the multilayer optical coating is modelled as a solid 3D element having the total thickness corresponding to the sum of the thicknesses of the high reflective and low reflective materials and mechanical parameters averaged over the thicknesses of the layers. Instead of 3D shell elements, we have used 3D solid elements also for very thin materials, though more CPU time consuming, because they provide the shear deformations and energies, which are useful for getting the mechanical losses associated to the modes.

B. FEA simulations results

The flats, the ears and also the anchors play an important role. In particular, since they break the cylindrical symmetry, they lift degeneracies and increase the number of distinct mode frequencies. In this paper we will discuss the results up to 70 kHz. To estimate the accuracy of the model, we have used a set of frequencies (ν Meas ) measured on the North arm IM up to 40kHz of an IM. Fig. 3 shows relative differences (ν Meas -ν FEA )/ν FEA , versus the frequency of the FEA ν FEA . The standard deviation is 0.15%. v

We have estimated the quality factors of the mechanical modes of the IM taking into account several kinds of losses: losses of the fused silica substrate, anchors and supports of the magnets (loss angle φ FS ); coating losses (loss angle φ IMcoating ); losses of bonding layers used to attach the ears, the anchors and the magnets (loss angle φ Bonding ). The bonding layers, have a thickness of 60 nm, and are modeled as 3D solid elements. Coating losses of the IM and EM were recently measured [20]. Note that all the parameters used are given in table I. Each loss contributor is related to the energy fraction stored in the

φ Bonding • E tot = φ HCB • E bonds φ IMcoating • E tot = φ CIM • E CIM φ FS • E tot = φ Suprasil • E FS . (4)
The overall loss angle for the IM is obtained by summing up all contributors:

φ IM = φ Bonding + φ IMcoating + φ FS .
The mechanical quality factor of the IM modes then writes Q m = 1/φ IM . Fig. 4 shows the frequency dependence of the FS substrate loss and the effect of adding the optical coating and the bonding layers. The influence of the bonding term φ Bonding is strongly mode shape dependent through the deformation of the ear and anchor bulks and it is not negligible. In fact, its contribution to Q m is dominant. For this reason, from a set of Q measurements it is possible to infer the value of φ HCB by using the energy fractions calculated with the FEA. Fig. 5 shows the Q m of the IM mass computed by fitting the loss angle φ HCB by using the first set of 5 modes of the IM of the north arm and supposing that it is does not vary with the frequency. At frequencies higher than 10 kHz, the bondings have a strong damping effect, though they have a negligible effect on the thermal noise of the IM. This is a very important result for the parametric gains computation and consequently for identification of the unstable modes.

IV. TRANSVERSE OPTICAL MODES IN ARM CAVITIES

Hermite-Gauss modes (HGM) are solutions of the paraxial wave equation for infinite-sized spherical mirrors. This mode basis was used in [17] to compute the parametric gain for the LIGO interferometer. It is fast to implement as the mode shapes are provided by analytical formulas. However, it restricts the mirror model to a purely spherical shape of infinite size. In particular, it does not include the effects of the deviations from the spherical shape due to fabrication imperfections or thermal effects. Finally, it does not take account for finite size effects such as diffraction losses, which must be estimated separately.

We have computed another set of optical modes that are obtained from a numerical resolution of the paraxial equation with finite-sized mirrors [START_REF] Kogelnik | [END_REF]. This mode basis will be referred to as 'finite-sized mirror modes (FSMM)'. Contrary to HGM, FSMM are obtained directly with diffraction losses. Moreover, mirror shapes can be chosen arbitrary, which enables one to introduce any deformation of the mirrors due to thermal effects or fabrication imperfections. Note, that in this work, we did not include fabrication imperfections, which effects will be the subject of future work.

In the following, we analyze the differences in Gouy phase (or frequency), diffraction loss, and mode shape, and between the HGM and FSMM basis set.

A. Gouy phases

The Gouy phases of arm cavity modes set the optical resonant frequencies, and, thus, the OPI resonance condition ω m = δω, where δω is the difference in frequency between the zero order mode and the higher order transverse optical mode. In the case of HGM, the Gouy phase of the mode of linear index n is given by

φ Gn = O n φ G ( 5
)
where O n the order of HGM(n), and φ G is the Gouy phase of the lowest order mode HGM(1) (usually referred to as TEM 00 in the literature) given by Here, g 1 < 0 and g 2 < 0 are the g parameters of the interferometer arm cavities. For the Advanced Virgo arm cavities, φ G 2.74 rad. and can be tuned by small variations of the mirror radii of curvature. Fig. 6 shows the difference between HGM's and FSMM's Gouy phases, expressed in units of free spectral range on the left vertical axis and in units of cavity linewidth on the vertical right axis. Note that the Gouy phases have been wrapped within an interval of length π, which allows to fold all the modes within a single free spectral range. The green line splits the graph into two regions: in the above region, the deviation is more than half a cavity linewidth, and we expect the model choice to have an impact on the OPI gain, whereas in the bottom part the impact should be negligible. Thus, the critical order is 7.

φ G = arccos (- √ g 1 g 2 ). (6)

B. Diffraction losses

Diffraction losses stem from the finite size of the cavity mirrors. They are a key parameter to compute the parametric gain, since they contribute to the optical linewidth (together with material absorption losses, scattering losses, and mirror transmittance). Since low-order modes have most of the energy concentrated at the center of the mirror, their diffraction losses are small, whereas high-order modes spread over a larger surface and show higher diffraction losses. Thus, in general, high-order modes are less likely to contribute to a PI. However, note that counter-intuitively, a loss increase can some- times lead to a higher parametric gain, as explained in more details in section VI B. Diffraction losses for HGM are estimated, like in [17], by evaluating the ratio between the total light flux within the coating radius of a mirror and the total flux incident on the mirror. Figure 7 shows diffraction losses for both sets of modes. It shows that with this rough estimation method, besides the HGM(1), all HGM have their diffraction losses underestimated. However, note that the total losses (input mirror transmittance plus diffraction losses) of low-order modes are dominated by the input mirror transmittance (green line on Fig. 7). Thus, the total losses obtained with the two methods start to differ by more than 10% around mode order 5.

C. Mode amplitudes

Optical mode amplitudes are used to compute the three mode spatial overlap coefficient B mn of Eq. 1. Thus, they also directly affect the OPI gain. In order to compare the mode amplitudes of FSMM and HGM, we decompose the vectors of one basis set onto the other by using the decomposition coefficient c ij of any FSMM (index i) with any HGM (index j):

c ij = S dxdy u * i (x, y)v j (x, y), (7)
where i and j are modes integer indices, u i (resp. v j ) are the FSMM (resp. HGM) mode amplitudes, and S is the mirror coating surface. Note that the transverse profile of a FSMM is constrained on a disk (mirror coating), whereas for HGM the transverse profile is distributed over the whole plane, such that a linear superposition of FSMM will never exactly match a HGM, and a true transformation matrix between the two basis set cannot rigorously be obtained [START_REF] Siegman | Lasers[END_REF]. In Fig. 8(a), we represent |c ij | for i = 2, and j ∈ {1, 2, ..., 36}. We find that FSMM(2) is a linear combination of the two order one HGM which are the HGM(2) and the HGM(3). We find that this is true for orders below 7. Conversely, as shown in Fig. 8(b), the higher order FSMM(36) mode (shown in the inset of Fig. 8(b)) cannot be decomposed on a single order of HGM. In the presented case, it is a mixture of order 7, 9, 11, and many other higher odd orders that are not shown on the figure.

D. Conclusion

This study shows that, in the absence of mirror deformation, the HGM basis does not deviate significantly from the FSMM basis for modes of order lower than 6. For order 6 and higher, the more resource consuming FSMM basis should lead to significantly different results for OPI gains. In section VI, we compare the OPI gains obtained for the Advanced Virgo O3 configuration, with HGM and FSMM basis set.

V. THERMAL EFFECTS

The laser energy is partially absorbed both by coatings and in the bulk of mirrors. This causes a temperature gradient, which originates two effects. First, a gradient of refractive index in the bulk of input mirrors modifies the mode matching condition, but affects neither the cavity linewidth nor the mode frequencies. Second, a deformation of the mirror surface, which modifies the mode shapes and frequencies. In this part, we evaluate the impact on this second effect on the properties of cavity modes by comparing FSMM obtained for purely spherical mirrors with FSMM obtained for thermally deformed mirrors.

The deformation profile is obtained by solving the linear thermoelastic equations [START_REF] Vinet | The virgo physics book: Optics and related topics[END_REF]. Figure 9(a) shows the purely spherical and thermally deformed profiles of an Advanced Virgo input mirror, for an intracavity power of 300 kW. We fitted the central part of the deformed mirror to extract a radius of curvature. The results are given in the following table :  No thermal effect With thermal effect Input mirror 1424.6 m 1432.1 m End mirror 1695 m 1702.3 m However, note that the mirror is not spherical anymore and the result of the fit is only valid in the center. In order to evaluate the incidence of this effect on the optical cavity parameters, we computed the FSMM with and without thermal effect on the two cavity mirrors. Fig. 9(b) shows the frequency differences between the two situations. We see that optical modes acquire a significantly different Gouy phase even for very low mode orders. We checked that losses and mode amplitudes are affected only for orders higher than 7, such that the frequency shift is the main effect. In section VI, we study the impact of this phenomenon on OPI gains.

VI. PARAMETRIC GAIN COMPUTATION

A. Validation: comparison with the Finesse software

The OPI gains of all mechanical modes within the [5.7, 70.7] kHz range were computed using Eq. (1) using both the HGM and FSMM basis set. In order to validate this method, we compared our results with the one obtained with the Finesse software [START_REF] Brown | The software and source code are available at[END_REF][START_REF] Green | [END_REF]. The OPI gain obtained with the Finesse software for one mechanical mode and two arm cavity mirrors is shown in Fig. 10(a). In Fig. 10(b), we plot the relative difference with the OPI gain obtained with the Finesse software and with Eq. (1) using FSMM and HGM. We observe a difference of a few percent at maximum. Note that the slight asymmetry between the blue and red curves stems from the small parameter difference between the two arms cavities. This comparison has been performed with many other mechanical modes and showed similar results. Note that using Eq. ( 1) is much faster than using the Finesse software, and that computing the results of the following figures would not have been possible in a reasonable amount of time. Therefore, in the following we use only Eq. (1).

B. Effect of optical losses on the OPI gain

In this section we demonstrate a counter-intuitive effect of optical losses on the OPI gains. Intuitively, if optical losses increase, the OPI gains get lower since the optical linewidth also increases. Here we show that if the OPI resonance condition is not exactly fulfilled, broadening the optical mode response can increase the gain such that counter intuitively, the gain variation does not vary monotonously with the diffraction loss. This is best shown on Fig. 11(a), where the OPI gain of a mechanical mode is plotted against optical diffraction losses of the main optical contributor. In this example, the gain first increases from around 0.04 below 10 2 ppm to 0.1 at 2 × 10 4 ppm, before decreasing at higher loss values, as expected. This appears also in Fig. 11(b), where the optical gain G n of the main optical contributor to the OPI gain of the mechanical mode of Fig. 11(a) is represented as a function of the mechanical mode frequency, for two different values of diffraction losses. At low losses, the two resonance peaks are well separated, such that there is a minimum in between (black arrow). A loss increase from 0 to 21540 ppm (red bullet at maximum of Fig. 11(a)) broadens the peaks an lead to the red curve which has no minimum anymore, and which shows higher Optical gain Gn of a FSMM high order mode versus the mechanical mode frequency, which is artificially varied around the resonance condition δω = ωm (with δω = 12.5551 kHz). The blue line is for null diffraction losses, and the red one is for 21540 ppm (red bullet at the maximum of the curve in (a)). The arrow point the minimum in between the two resonance on the blue curve, where the gain increase is maximum. The gray shaded area highlights a frequency at which Gn is higher when the diffraction loss is higher.

values in a whole frequency region (gray shaded area on Fig. 11(b)). Finally, if the losses were increased further, the red curve would start lowering and the gray shaded area would vanish. In this paragraph, we study the impact of the model used to compute the optical modes. We compare the OPI gains obtained with the HGM and FSMM basis. As expected, we find that there is only a marginal difference between the two models if optical modes of order below 5 are involved. In Fig. 12, we plot the gain of a mechanical mode versus its frequency using the two optical mode basis set. This mode has been chosen for the main optical contributor to the OPI gain is an order 6 optical mode. There is a factor 3 between the two gain maxima and the two peaks are shifted by around 100 Hz, which corresponds to the optical linewidth. This is in agreement with the conclusions of Sec. IV.

D. OPI gain computation in the O3 configuration

The simulations have been performed for the Advanced Virgo configuration corresponding to that of O3. The parameters for such a configuration are shown in Table II. Measured parameters have been included rather than nominal values when they were available. Only the optical input power has been set to the nominal value of 50 W, which is the maximum value that would have been possibly used during O3 (the value effectively reached in O3 being too small to trigger any instability in the range of mechanical mode frequencies simulated in this work).

The corresponding arm-cavity power is around 300 kW. The parameters used are listed in table II.

To account for optical mode frequency uncertainties, we present OPI data in two-dimensional plots (see for instance Fig. 13(a) and (b)), where the end mirrors radii of curvature NE and WE are scanned. The color code indicates the gain value at each interferometer working point. This choice is also related to the fact that the envisioned OPI mitigation technique relies on ring heaters able to tune the end mirrors radii of curvatures [27]. In the following, we show two sets of results. Each figure is the result of the same OPI calculation but using a different set of optical modes. Fig. 13(a,b,c) shows the results for FSMM, Fig. 13(d,e,f) for FSMM including thermal effect due to coating absorption (see sec. V). NE and WE are scanned over a five meters range, which is within reach of the mirror ring heater system. In each OPI plot, the color code is chosen such that the gray scale is for gains lower than 1 (no instability), and the color scale is for R m > 1. The involved mechanical modes are indicated in the inset, and main optical modes contributors are shown below each OPI plot. Note that the result obtained with HGM are indistinguishable from that of Fig. 13(a), such that we did not include the corresponding figure. Indeed, only low order optical modes are involved here, such that HGM and FSMM give the same result. Fig. 13(b resp. d) shows the unstable mechanical modes of Fig. 13(a resp. b) on the first line. The second line shows the involved optical modes. In Fig. 13(c resp. f), we plot the number of unstable modes in the range of Fig. 13(a resp. b) versus the optical power for FSMM with) thermal effect. Modes that ring on different mirrors are counted only once. The blue curves the ratio of the area free of instability S Rm>1 to the total area S tot in Fig. 13 (a resp. b), versus the optical power. This ratio quantifies how difficult it is to escape an unstable area within the accessible radii of curvature range. The green vertical lines on Fig. 13(c andf) point the nominal power of O3 (50 W) and the power that was effectively reached (28 W).

These results show that an OPI involving mechanical modes with frequencies below 70 kHz could have been observed at the nominal power of O3, although it would be easily escaped with end-mirrors ring heaters since S Rm>1 /S tot 0.2 at 300 kW cavity power. It also shows that, at the real power of O3, it was very unlikely to observe an instability. Finally, it shows that the thermal effect caused by coating absorption has an important impact on the results at such high optical powers.

VII. CONCLUSION

In this letter, we have presented OPI gain simulations in the Advanced Virgo configuration of O3. Compared to previous work [17], we have used deeper physical modeling, including a very detailed description of mechanical modes and optical modes. The aim being that our model becomes predictive. The mechanical mode simulation includes all the mirrors details, and we implemented an original method to obtain precise quality factors value for all modes, by using the combination of FEA and ringdown measurement that are performed on a subset of modes. We have also provided a precise description of optical modes, by considering finite size effects. Our method provides directly accurate diffraction losses. Furthermore, we have shown that, counterintuitively, an optical loss increase can lead to a parametric gain increase. Our conclusion regarding optical modes is that up to order 4 (included), analytical formulas for Hermite-Gauss modes are sufficient to predict accurate OPI gains. However, if higher-order optical modes are involved, mirrors finite size effects must be accounted for. Finally, we have shown that the mirror deformation stemming from the laser absorption in mirror coatings plays an important role and must be included in the OPI simulation. These simulations pave the way towards precise optomechanical instability predictions for the current and next generations of gravitational-wave detectors.

  FIG. 1. O3 Advanced Virgo's configuration. NI and NE are respectively the input and end-mirror of the North-arm, and WI and WE are respectively the input and end-mirror of the West-arm. PR is the Power Recycling mirror, BS the interferometer beam splitter.

FIG. 3 .

 3 FIG.3. Relative differences of measured frequencies with respect to frequencies obtained with the FEA vs the FEA frequencies. The standard deviation is 0.15 • %

FSFIG. 4 .FIG. 5 .

 45 FIG. 4. Loss angles obtained from the FEA of the Input TM. The computation ha been performed up to 70 kHz.

FIG. 6 .

 6 FIG.6. Difference between HGM's and FSMM's Gouy phase, expressed in units of free spectral range on the left vertical axis, and in units of cavity linewidth on the right vertical axis. The vertical red dashed lines highlight the mode orders, which appear above the upper horizontal axis. The horizontal lower axis shows the optical mode index (modes are sorted by increasing energies). The green horizontal line has a vertical coordinate of 0.5 on the right axis.

FIG. 7 .

 7 FIG.7. Diffraction losses obtained for FSMM (cross) and estimated for HGM (circles). The vertical red dashed lines underline the mode orders, which appear in the upper horizontal axis. The horizontal lower axis shows the optical mode index (modes are sorted by increasing energies). The green line shows the input mirror transmittance.

FIG. 8 .

 8 FIG. 8. (a) |cij| terms of the decomposition of the FSMM(2) on the HGM basis: the FSMM(2) is a linear combination of the HGM(2) and the HGM(3). (b) |cij| terms of the decomposition of the FSMM(36) on the HGM basis.

FIG. 9 .

 9 FIG. 9. (a) Input mirror profiles with (dashed red line) and without thermal effect (blue solid line). (b) Frequency difference between FSMM modes with and without thermal effect, expressed in unit of FSR on the left axis, and in units of linewidth on the right axis. The green horizontal line is at vertical position 0.5 on the right axis.

FIG. 10

 10 FIG. 10. (a) Parametric gain obtained with the Finesse software for one mechanical mode (12.552 kHz) on mirror NE (blue) and WE (red). (b) Relative difference between (a) and the parametric gain Rm obtained with Eq. (1) using FSMM (solid lines) and HGM (dashed lines).

FIG. 11 .

 11 FIG. 11. (a)Parametric gain Rm of a mechanical mode of frequency 12.5551 kHz while varying artificially the diffraction losses of optical modes. (b) Optical gain Gn of a FSMM high order mode versus the mechanical mode frequency, which is artificially varied around the resonance condition δω = ωm (with δω = 12.5551 kHz). The blue line is for null diffraction losses, and the red one is for 21540 ppm (red bullet at the maximum of the curve in (a)). The arrow point the minimum in between the two resonance on the blue curve, where the gain increase is maximum. The gray shaded area highlights a frequency at which Gn is higher when the diffraction loss is higher.

FIG. 12 .

 12 FIG. 12. Parametric gain Rm of the 33.750 kHz mechanical mode (see inset) while artificially varying the mechanical frequency, using the two different optical mode basis (red: HGM, blue: FSMM). The dashed lines point the maxima of the two curves, and emphasize the height and frequency change.

FIG. 13 .

 13 FIG. 13. (a) Parametric gainRm versus radii of curvature of NE and WE, using FSMM without thermal effect, for an input power of 50 W. The gray-shaded scale highlights the parametric gains lower than 1, while the colorful scale highlights the instabilities (Rm > 1). The three unstable mechanical modes are labeled A, B, and C. Their frequencies are respectively 61.154 kHz, 61.160 kHz, and 66.784 kHz. (b) The mechanical mode shapes and associated transverse optical mode contributing the most to the aforementioned instabilities. α contributes for 93 % of A and for 77 % of B, β contributes for 6 % of A and for 15 % of B, γ contributes for 89 % of C, and δ for 10 % of C. (c) The red solid line is the number of unstable modes in the radius of curvature range of (a), with respect to the optical input power (top) or intracavity power (bottom). They become unstable from an optical input power of respectively 9 W (mode A), 26 W (mode B), and 27 W (mode C). The blue curve is the ratio of the area free of instability SRm>1 to the total area Stot in Fig.13 (a resp. b), versus the optical power. The figures (d) resp. (e) resp. (f) are the same than (a) resp. (b) resp. (c), but taking into account the thermal effect. The five unstable mechanical modes are labeled D, E, F, G, and C (same mode than without thermal effect). Their frequencies are respectively 66.888 kHz, 66.912 kHz, 61.216 kHz, 61.231 kHz, and 66.784 kHz. contributes for 84 % of D and 24 % of E, ζ contributes for 16 % of D and 76 % of E, η contributes for 100 % of F and G, and θ contributes for 97 % of C. They become unstable from an optical input power of respectively 4.6 W (mode D), 35 W (mode E), 42 W (mode F), 51 W (mode G), and 57 W (mode C). The green dashed lines show the input power reached in O3 (28 W), and the nominal power (50 W) corresponding to (a) and (b).

ACKNOWLEDGMENTS

We acknowledge Jérôme Degallaix, Ettore Majorana, François Bondu and Nelson Christensen for helping us to improve the manuscript. We thank Mikhaël Pichot Du Mezeray for providing us with the code to compute the mirror deformation due to thermal effect, and Mourad Merzougui for discussions about Ansys calculations. We also thank Slawomir Gras, Jean-Yves Vinet and Walid Chaibi for interesting discussions and advice, Andreas Freise and Anna Green for their help with the Finesse software. Finally, David Cohen is supported by the DIM ACAV+ from the Region Île-de-France.


Ghismon-De-Kasin Mayinda

Delphine Serreau

Amandine Gesbert

Fabrice Reigner

Jean- François Sutra

Anne Lespine

Guillaume Sallé
email: guillaume.salle@inrae.fr

Ivermectin treatment in lactating mares results in suboptimal ivermectin exposure in their suckling foals

Keywords: ivermectin, horse, strongyle, nematode, foal

The management of equine strongyles has become problematic over the last decade because of an increased prevalence of drug-resistant isolates worldwide. Therapeutic options are therefore limited, leaving macrocyclic lactones as the most often effective drug class. However, their lipophilic properties result in a long-lasting elimination that could favour drug resistance selection. As a result, ivermectin treatment in lactating mares could promote suboptimal exposure of their foal parasites to ivermectin, thereby selecting for more resistant worms. To test for this putative transfer, we selected two groups of six foal-mare pairs, one group of mares receiving ivermectin and the other being left untreated. We compared faecal egg count trajectories in foals from the two groups and quantified plasma ivermectin concentrations in ivermectin treated mares and their foals during seven days. Our results showed limited but sustained plasmatic exposure of foals associated with non-significant faecal egg count reduction (P = 0.69). This suggests that ivermectin treatment in lactating mares results in suboptimal exposure to the drug in their foal.

Introduction

Equine gastro-intestinal nematodes emerge as a major threat to equine health. Cyathostominae and Parascaris spp. are the most problematic species as they account for most deaths in young horses while tapeworms and Strongylus spp. seem currently under control [START_REF] Sallé | Compilation of 29 years of postmortem examinations identifies major shifts in equine parasite prevalence from 2000 onwards[END_REF].

Cyathostominae are responsible for clinical disease occurring when encysted larval stages emerge at once [START_REF] Love | Pathogenicity of cyathostome infection[END_REF]. This so-called "larval cyathostominosis" is associated with colic, protein-losing enteropathy and eventually death of the animals [START_REF] Giles | Larval cyathostomiasis (immature trichonema-induced enteropathy): a report of 15 clinical cases[END_REF]Sallé et al., 2020). Parascaris spp. mostly infect foals under the age of one year and recent observations suggested that foal infection follows a bivariate pattern with worm counts peaking at five-and nine-month of age under treatment-free conditions [START_REF] Fabiani | Dynamics of Parascaris and Strongylus spp. parasites in untreated juvenile horses[END_REF]. Complex migration of larval stages from the intestine, through the liver and lungs can result in coughing and nasal discharge [START_REF] Taylor | Veterinary Parasitology[END_REF]. Heavy infection by adult stages is associated with foal emaciation, catarrhal enteritis and it can lead to intestinal impaction or rupture [START_REF] Cribb | Acute small intestinal obstruction associated with Parascaris equorum infection in young horses: 25 cases[END_REF][START_REF] Sallé | Compilation of 29 years of postmortem examinations identifies major shifts in equine parasite prevalence from 2000 onwards[END_REF].

Control of these parasite species by anthelmintic drugs has become problematic over the last decades because of widespread emergence of drug-resistant isolates. Fenbendazole has been found inefficient against small strongyles in most countries [START_REF] Hodgkinson | The role of polymorphisms at beta tubulin isotype 1 codons 167 and 200 in benzimidazole resistance in cyathostomins[END_REF][START_REF] Lyons | Evaluation of parasiticidal activity of fenbendazole, ivermectin, oxibendazole, and pyrantel pamoate in horse foals with emphasis on ascarids (Parascaris equorum) in field studies on five farms in Central Kentucky in 2007[END_REF][START_REF] Kyvsgaard | Prevalence of strongyles and efficacy of fenbendazole and ivermectin in working horses in El Sauce, Nicaragua[END_REF][START_REF] Traversa | Efficacy of major anthelmintics against horse cyathostomins in France[END_REF][START_REF] Kumar | Benzimidazole resistance in equine cyathostomins in India[END_REF] while suboptimal ivermectin efficacy has been frequently reported for Parascaris sp. [START_REF] Lindgren | Parascaris equorum in foals and in their environment on a Swedish stud farm, with notes on treatment failure of ivermectin[END_REF][START_REF] Lyons | Evaluation of parasiticidal activity of fenbendazole, ivermectin, oxibendazole, and pyrantel pamoate in horse foals with emphasis on ascarids (Parascaris equorum) in field studies on five farms in Central Kentucky in 2007[END_REF][START_REF] Näreaho | Impaired efficacy of ivermectin against Parascaris equorum, and both ivermectin and pyrantel against strongyle infections in trotter foals in Finland[END_REF][START_REF] Laugier | Prevalence of Parascaris equorum infection in foals on French stud farms and first report of ivermectin-resistant P. equorum populations in France[END_REF][START_REF] Morris | Anthelmintic efficacy of single active and combination products against commonly occurring parasites in foals[END_REF][START_REF] Studzińska | A survey of ivermectin resistance in Parascaris species infected foals in south-eastern Poland[END_REF]. Expected efficacies of ivermectin and moxidectin are usually observed against small strongyles although failures have been reported in Brazil and western Europe [START_REF] Peregrine | Anthelmintic resistance in important parasites of horses: does it really matter?[END_REF]..

To date, risk factors favouring the selection of drug-resistant parasite populations are not fully identified. Two recent surveys in the USA [START_REF] Nielsen | Anthelmintic efficacy against equine strongyles in the United States[END_REF] or in France [START_REF] Sallé | Risk factor analysis of equine strongyle resistance to anthelmintics[END_REF] showed limited convergence in their conclusions regarding the impact of parasite management strategies, but identified pasture rotation as a beneficial strategy to limit drug resistance development.

Ivermectin and other macrocyclic lactones are slowly eliminated, thereby exposing incoming parasite infective larvae to reduced concentrations through time [START_REF] Le Jambre | Selection for anthelmintic resistance by macrocyclic lactones in Haemonchus contortus[END_REF]. These suboptimal drug concentrations promote the selection of resistant worms while eliminating susceptible individuals, in a so-called "tail-selection" process [START_REF] Le Jambre | Selection for anthelmintic resistance by macrocyclic lactones in Haemonchus contortus[END_REF]. Because of their lipophilic properties (Geary, 2005), macrocyclic lactones are found in the milk of lactating individuals including horses [START_REF] Gokbulut | Comparative plasma and milk dispositions, faecal excretion and efficacy of per os ivermectin and pour-on eprinomectin in horses[END_REF]. As a result, suckling offsprings may be significantly exposed to the drug as quantified in cattle [START_REF] Alvinerie | Plasma profile study of moxidectin in a cow and its suckling calf[END_REF]. In sheep, two independent experiments demonstrated that treatment of ewes with moxidectin could impact on the parasite population in their lambs, either reducing by 70% the establishment of susceptible Teladorsagia circumcincta [START_REF] Leathwick | Selection for anthelmintic resistant Teladorsagia circumcincta in pre-weaned lambs by treating their dams with long-acting moxidectin injection[END_REF] or inducing a 51% cut in Haemonchus contortus excretion [START_REF] Dever | Decline in faecal worm egg counts in lambs suckling ewes treated with lipophilic anthelmintics: Implications for hastening development of anthelmintic resistance[END_REF]. There is hence an opportunity that ivermectin treatment given to lactating mares as part of calendar treatments [START_REF] Sallé | Risk factor analysis of equine strongyle resistance to anthelmintics[END_REF] may affect their foal parasite communities. To the best of our knowledge, ivermectin transfer from mares to suckling foals has never been characterized.

Here, we tested the possible passage of ivermectin from the mare milk to their foals by comparing faecal egg count and plasma drug concentrations in foals whose mare were either treated or not in late summer. Our results showed non-significant faecal egg count reduction in combination with limited but sustained plasmatic exposure of foals..

Materials and methods

Experimental design

To investigate the impact of ivermectin transfer between treated mares and their foal, we used fourteen Welsh pony mares and their foals. Pregnant Welsh pony mares were allocated to two 7 ha pastures in April 2019. Pasture allocation was balanced for the stallion they were mated with, and for the foreseen date of foaling. This latter parameter was considered to obtain similar foal parasite exposure across pastures. Foaling took place from May 25 to June 23 2019. Detailed metadata are provided in supplementary Table 1.

Faecal egg count (FEC) was performed every month (from early June to late August) in the seven mare-foal pairs (supplementary Table 1) to quantify strongyle egg excretion. Foals were not sampled before they reached 14 days of age. For each pair, animals were placed in a paddock to collect fresh faecal material off the ground. In some cases, foals did not expel faecal matter directly; they were hence separated from their mother for 10 min to stimulate defecation.At the end of August (August 30, 2019), six mares with faecal egg count above 150 eggs/g [START_REF] Coles | The detection of anthelmintic resistance in nematodes of veterinary importance[END_REF] were administered ivermectin orally (Eraquell,Virbac,France,200 µg/kg bodyweight) to determine whether ivermectin treatment in the mare could affect parasite egg excretion in foals. The other 8 mares and their foals remained untreated and were used as a control group. FEC was measured 18 days later (September 17, 2019) to determine treatment efficacy in both mares and foals.

FEC measure

FEC were measured from 5 g of faecal material diluted in 70 mL of a saturated NaCl solution (specific gravity = 1.18). Filtrate was loaded on a McMaster slide and eggs were counted using light microscopy at 150X magnification (detection limit of 15 eggs/g).

Ivermectin dosage and analysis

Plasmatic ivermectin concentration was measured to i) ensure that plasmatic ivermectin exposure in mares was in line with observed FECR, ii) support ivermectin transfer to the milk and iii) evaluate plasmatic exposure in foals. For this reason, we restricted the study to the first seven days. Although ivermectin residence time in adult horses expands beyond seven days [START_REF] Pérez | Plasma Profiles of Ivermectin in Horses following Oral or Intramuscular Administration[END_REF], restricting our study to the first seven days was sufficient to address our main objectives while reducing animal use. For ivermectin dosage, milk (4 mL silica-coated tubes, BD Vacutainer ®, Plymouth, UK) and blood (4 mL heparin-coated tubes, BD Vacutainer ®, Plymouth, UK) samples were taken from six mares at 1 h, 2 h, 8 h, 24 h, 48 h, 4 days and 7 days after treatment. To measure putative ivermectin transfer in their foals (n = 6), blood samples were collected from the foals at the same times from 8 h onwards. Blood samples were centrifuged at 1500 g and 4 °C for 10 min to harvest plasma. Both milk and plasma samples were cryopreserved at -80 °C until ivermectin concentration was subsequently measured.

Ivermectin concentrations were determined by High-Performance Liquid Chromatography with fluorescence detection according to previously validated methods in plasma [START_REF] Alvinerie | Persistence of ivermectin in plasma and faeces following administration of a sustainedrelease bolus to cattle[END_REF] and in milk [START_REF] Alvinerie | Ivermectin in goat plasma and milk after subcutaneous injection[END_REF]. Data were analyzed by a non-compartmental approach using the KineticaTm software (version 4.2, innaPhase, Philadelphia, USA). Apparent and maximal concentration (C max ) and time of C max (T max ) were then determined. The partial area under the plasma concentration time curve (AUC) and the mean residence time (MRT) were calculated from 0 to 7 days by the linear trapezoidal rule. Data are expressed as mean ± standard deviation (SD) of 6 animals for mares and foals respectively.

Statistical analysis

Statistical analysis on FEC data was implemented in R v4.0.2 (R Core Team, 2016).

Treatment efficacy was assessed by FEC reduction in the treated group corrected for FEC trajectory in the control group [START_REF] Presidente | Methods for detection of resistance to anthelmintics[END_REF] :

,

where FEC 14 and FEC 0 correspond to group average FEC measured before and after treament respectively. This formula was applied to measure ivermectin efficacy in mares and in foals whose mares had been treated. Efficacy in foals was estimated from the only foals with patent strongyle infection (n = 5 in the mare treated group and 6 in the control group). Efficacy below 0% were censored to 0. Confidence intervals (95% CI) were estimated using the R-version of the eggCounts package v2.3 [START_REF] Wang | Modelling anthelmintic resistance by extending eggCounts package to allow individual efficacy[END_REF].

In addition, to determine whether FEC reduction was significantly different between experimental groups, we modeled FEC data with a marginal modeling approach using the geeM package v.0.10.1 [START_REF] Mcdaniel | Fast pure R implementation of GEE: application of the Matrix package[END_REF]. The implemented model was: FEC ij ~  0 +  1 x day ij +  2 x group ij +  3 x day ij x group ij , where FEC ij stands for the j th FEC record in individual i, with an auto-regressive correlation structure between FEC ij and FEC ij+1 ,  0 is the intercept,  1 to  3 are the regression coefficients describing the covariate effects (day, treatment group and their interaction). Because of the over-dispersed nature of FEC, a negative-binomial function was used as a link function between the linear combinations of covariates and observed FECSignificance of the interaction term (P-value < 0.05) would indicate egg count reduction between groups.

Data and R script used for this study are under the following repository: https://github.com/guiSalle/Ivermectin-in-lactating-ponies.

Results

FEC trajectory over the grazing season in mares and foals

Summary statistics for FEC measured in mares and foals are shown in Table 1 and detailed trajectories are available in supplementary Table 1 and supplementary Fig. showed increased excretion of strongyle eggs from June to mid-August, reaching 464 eggs/g (range between 0 and 2520 eggs/g) on average at that time. In foals, strongyle egg excretion peaked at the end of August (227 eggs/g on average, ranging from 0 to 570 eggs/g). Of note, one foal remained FEC negative throughout the study period.

Ivermectin treatment in mares did not significantly affect foal strongyle egg excretion

FEC reduction was measured between foals whose mare received ivermectin treatment or not over an 18-day period. Fig.Following suppressive ivermectin treatment in six mares (Table 1), reduced egg excretion was found in four foals (Fig.Table 1), although FEC increased in one foal from 30 to 45 egg/g (one foal remained FEC negative). Average FEC reduction in the foals whose mares were treated was 67% (95% CI 6.2 -85%). However, foals from the control group also demonstrated 55% reduction of their FEC (95% CI 10.4 -90.8%). FEC reduction accounting for the control group FEC trajectory was hence 0%. In agreement with this result, the marginal modeling approach found no evidence of statistical difference in FEC reduction between both groups (reduction ratio of 1.33 between foals whose mares received ivermectin and the control group; 95% CI 0.31 -5.64, P = 0.69).

Ivermectin accumulates in mare milk but not in foal plasma

To control that mare exposure to ivermectin was achieved and to explore exposure in their respective foal, plasmatic ivermectin concentration was measured in mare-foal pairs as represented on Fig. 2) and occurred 8 h after treatment (Fig.Fig. 1).

Milk ivermectin concentration followed a similar pharmacokinetic pattern than that found in plasma (Pearson's correlation r = 0.61, P< 0.001) and also peaked at 8 h post-treatment in most animals (Fig. 1). In two mares, this was slightly delayed at 24 h or 48 h post-treatment (Fig. Fig. 1). It was slightly lower than plasma concentrations, with AUC milk to AUC plasma ratios ranging between 0.27 and 1.1 (mean 0.57 ± 0.27).

On the contrary, only limited concentrations of ivermectin were measured in foal plasma (Fig. 1) with maximal concentrations ranging between 0.06 and 0.66 ng/mL (Table 2). These values showed weak correlations with their mother plasmatic (Spearman's ⍴ = 0.03, P = 0.85) or milk concentrations (Spearman's ⍴ = 0.06, P = 0.85). The total exposure of foal (AUC) represented around 3% of the mare total drug exposure, and 7 % of that found in milk (Table 2).

Discussion

Calendar drenching schemes in horses often include an ivermectin treatment given in late summer. Because of its lipophilic properties, ivermectin can be excreted into mare milk [START_REF] Gokbulut | Comparative plasma and milk dispositions, faecal excretion and efficacy of per os ivermectin and pour-on eprinomectin in horses[END_REF]. This leaves the opportunity for a suboptimal ivermectin exposure in foals. However, the transfer from mare milk to foal plasma was significant but low and under the minimal active concentration. In line with this, no FEC reduction occurred in foals whose mares were administered ivermectin relative to that with an untreated mother.

The pharmacokinetic data of ivermectin in mare plasma were in good agreement with previous studies conducted in horses that received ivermectin orally [START_REF] Gokbulut | Comparative plasma and milk dispositions, faecal excretion and efficacy of per os ivermectin and pour-on eprinomectin in horses[END_REF][START_REF] Pérez | Plasma Profiles of Ivermectin in Horses following Oral or Intramuscular Administration[END_REF]. In addition, significant ivermectin excretion in pony mare milk was found with AUC milk /AUC plasma ratio of 0.57. This is consistent with ivermectin lipophilic properties [START_REF] Prichard | Moxidectin and the avermectins: Consanguinity but not identity[END_REF]. This remains however in the common range of values (AUC milk /AUC plasma between 0.7 and 1.0) reported for cows [START_REF] Toutain | Kinetic and insecticidal properties of ivermectin residues in the milk of dairy cows[END_REF] but lower than in ewes where AUC milk /AUC plasma can be as high as 2.5 [START_REF] Imperiale | Milk excretion of ivermectin and moxidectin in dairy sheep: assessment of drug residues during cheese elaboration and ripening period[END_REF]. The lower ratio observed for our pony population when compared to other species is consistent with the lower lipid content found in equid milk, e.g. 7.4 g/L in ponies [START_REF] Pagan | Composition of milk from pony mares fed various levels of digestible energy[END_REF] and between 5 to 20 g/L in horses [START_REF] Malacarne | Protein and fat composition of mare's milk: some nutritional remarks with reference to human and cow's milk[END_REF] that contrasts with cow or ewe milk that exhibit fat content in the range of 37 g/L [START_REF] Malacarne | Protein and fat composition of mare's milk: some nutritional remarks with reference to human and cow's milk[END_REF] and 71 g/L respectively [START_REF] Wohlt | Milk yield by Dorset ewes as affected by sibling status, sex and age of lamb, and measurement[END_REF]. However, the AUC milk /AUC plasma values calculated in our study were much higher than the 0.19 value previously reported in horses [START_REF] Gokbulut | Comparative plasma and milk dispositions, faecal excretion and efficacy of per os ivermectin and pour-on eprinomectin in horses[END_REF]. Although unclear, the basis for such difference may come from different factors including variation in milk lipid content between animals. Variation in the expression level of the efflux transporter ABCG2 that is located in mammary glands during lactation and plays an important role in the milk partitioning of drugs [START_REF] Jonker | The breast cancer resistance protein BCRP (ABCG2) concentrates drugs and carcinogenic xenotoxins into milk[END_REF] including endectocides [START_REF] Jani | Ivermectin interacts with human ABCG2[END_REF] may also contribute to this discrepancy.

Ivermectin was detected in foal plasma throughout the seven-day trial. This level matched the elevated ivermectin content found in milk of ivermectin-treated mares. To our knowledge, this is the first attempt to quantify ivermectin transfer from mare to their foal. However, ivermectin concentrations estimated in foal plasma were low, with C max found below the minimal active concentration of 1 ng/mL (Geary 1983). This is compatible with the lack of FEC reduction in the group of foals whose mothers had received ivermectin. It is yet to be determined whether this residual concentration results in a suboptimal exposure of parasite populations or is of no biological relevance. Underdosing has been indeed advocated as a major factor favoring the emergence of drug resistant isolates [START_REF] Le Jambre | Selection for anthelmintic resistance by macrocyclic lactones in Haemonchus contortus[END_REF]. More importantly, it should also favor the development of resistance in Parascaris spp. populations whose seasonal prevalence peak [START_REF] Fabiani | Dynamics of Parascaris and Strongylus spp. parasites in untreated juvenile horses[END_REF][START_REF] Sallé | Compilation of 29 years of postmortem examinations identifies major shifts in equine parasite prevalence from 2000 onwards[END_REF] would match autumn treatment in the mare.

Collection of larval populations to monitor their sensitivity to ivermectin using in vitro test could validate their putative increased resistance level.

The low ivermectin plasmatic concentration measured in the foal was however contrasting with the 0.57 AUC milk /AUC plasma transfer measured in the mare. This may be underpinned by the high fractional turnover rates of body water measured in foals as a result of a low energy to water content ratio [START_REF] Oftedal | Lactation in the Horse: Milk Composition and Intake by Foals[END_REF]. This could contribute to dilute ivermectin in the foal plasmatic compartment. In contrast, mean residence time in foal plasma expanded over two weeks. The close vicinity of measured concentration with the detection limit of our method certainly contributes to inflate mean residence time. It may also reflect a slow elimination process resulting from an immature drug detoxification system in the young animal as demonstrated in humans [START_REF] Brouwer | Human Ontogeny of Drug Transporters: Review and Recommendations of the Pediatric Transporter Working Group[END_REF]. This increased half-life of the drug would promote parasite exposure to suboptimal drug concentrations and potentially aggravate drug resistance development. FECR values were all negative in foals and censored to 0%; in mares FECR values were all 100%. As a result, no confidence interval is given.

Figure captions

Table 2. Pharmacokinetic parameters of ivermectin in plasma and milk estimated for mares and foals.

Parameters were estimated from six samples collected over a seven-day trial.

PK: pharmacokinetic; AUC: area under the plasma concentration time curve; MRT: mean residence time

Fig

  Fig. 1 and supplementary Fig. 2. Maximal ivermectin concentrations in mares ranged between 18.69 and 43.33 ng/mL (Fig. 1, Table

Figure 1 .

 1 Figure 1. Kinetics of ivermectin concentration in six treated mares and their foals

Table 1 .

 1 Summary statistics of strongyle faecal egg count in foals and mareSummary statistics (mean [ minimummaximum ]) of faecal egg count measured in mares and their foals from the treated (n = 6) or control group (n = 8). They are given before or after ivermectin treatment in mares. Within-group faecal egg count reduction (FECR) is provided.

Acknowledgements

GdK Mayinda benefited from a Région Normandie (PAMIEBO project) stipend and G. Sallé is a grateful recipient of an Institut Français du Cheval et de l'Équitation grant (CYATHOMIX project).

Animal welfare statement

Experimental procedures were carried out following the guidelines of the European Council Directive (2010/63/UE) and with French laws and regulations (Articles R214-87 to R214-137 of the Rural Code and decree n°2013-118 dated February 1, 2013 published on February 7, 2013). Our project was approved by French Ministry of Higher Education and Research and the Val de Loire ethics committee under license #20513-2019043009343644 v5.

Conclusion

Our data suggests that ivermectin treatment in lactating mares leads to limited transfer into their suckling foals. This had no effect on strongyle egg excretion in their progenies. The extent of parasite exposure to ivermectin remains to be determined.


domain: spi sdv

text:
Sepaldeep Singh Dhaliwal
email: sdhaliwal@univ-lille.fr

Sarah Wilby

Sina Firouzy

Kevin Brou Boni

Martijn De Vries

Stefan Escaida Navarro

Abdelkader Belarouci

Vincent Coelen

Othman Lakhal

David Pasquier

A L Palmer

W Polak

D Jones

A Labib

J J Dobbelsteen

C Duriez

R Merzouki

CoBra robot for localized cancer treatment and diagnosis under real-time MRI

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

I. Introduction

Prostate cancer is the most common and one of the major reasons for death among men. International agency for research on cancer estimates an increase of 16.9% incidence rate and 26.2% mortality rate in prostate cancer of European population by 2030 since 2020 [START_REF]Cancer Tomorrow[END_REF]. The CoBra project is developing a semi-autonomous robotized solution aiming to perform biopsy and treatment of localized cancers under MRI-guidance, initially for prostate cancer. MRI has superior image quality compared to trans-rectal ultrasound to detect small lesions, and has been explored for two decades for localized cancer treatment under real-time (RT) intraoperative MRI.

In the literature [START_REF] Monfaredi | MRI Robots for Needle-Based Interventions: Systems and Technology[END_REF], related to the MRI-guided prostate intervention robots, most demonstrated needle placement only. Thus, demanding the manual intervention of physicians for seed delivery. In addition, the key requirement is RT intraoperative needle intervention to take an advantage of MRI-guidance and robotic systems. The strong magnetic field of MRI restricts the use of conventional actuators and sensors.

We present an integrated design of the CoBra concept with different modules as sub-component systems with specific tasks to achieve the final mission of prostate brachytherapy (BT) or biopsy (BX). The modules can be stated as: an MRI-system, a dosimetric treatment planning system, a needle trajectory planning software (NPS), and the MRconditional robot. The CoBra partners have developed an MR-conditional robot, dosimetric algorithm, curved steerable needle, MR-biopsy module, needle trajectory planning software, and proof-of-concept for active bioinspired prostate test phantom (BIP).

II. Material and methods

The MR-conditional robot has been developed to perform in-bore prostate adaptive-BT and BX with the aid of a steerable needle targeting the lesion sites under RT-MRI. The robot is actuated using ultrasonic-piezo motors (USM) in closed-loop feedback with absolute linear encoders. The patient is positioned in lithotomy with the help of Uni-Lift, maintaining the patient in-bore intraoperatively in a comfortable position throughout the procedure [START_REF] Dhaliwal | Cooperative Brachytherapy for Prostate Cancer Under MRI Guidance[END_REF]. The robot with 5 degrees-of-freedom is capable of performing oblique insertions to avoid organs-at-risk (OARs) with minimal insertion points without the aid of any templateguide. The robot is designed with a quick-lock mechanism for an easy exchange of BT or BX modules. An automated seed delivery device has been developed for the BT module actuated using USM motors.

The novelty of the concept is to deliver dose adaptively under RT target tracking intraoperatively. The curved steerable needle is able to reach a new target position when pre-defined target shifts during needle insertion, prostate motion, tissue deformation, or inflammation. Targets are planned based upon an initial dosimetry plan under RT-MRI, then tracked and updated in parallel using rapid creation of accurate synthetic-CT (sCT) datasets from the live MR images and are used to update the NPS to deliver the seed with more accuracy using curved steerable needle. Thus, the dosimetry plan integrated to NPS and robot helps to deliver dose to the new target position.

In order to attain an RT dosimetry plan using sCT datasets, the team use an augmented cycle Generative Adversarial Network (AugCGAN) algorithm [START_REF] Brou Boni | MR to CT synthesis with multicenter data in the pelvic area using a conditional generative adversarial network[END_REF]. AugCGAN is more robust with the erraticism of MR-images as compared to the standard cycleGAN. This study incorporated T2w MR and CT pelvic images of 38 patients from 5 centers. The AugCGAN was trained on 2D transverse slices of 19 patients from 3 different sites. The network was then used to generate sCT images of 19 patients coming from two other sites. Mean Absolute Errors (MAE) for each patient were evaluated between real and sCT.

The novel steerable needle has been manufactured and omnidirectional steering at the distal tip can be achieved by a jointless integrated pull-push mechanism. The degree and direction of needle steering are controlled by four piezoelectric motors at the proximal end of the robotic BT module as shown in Fig. 1. The inner needle and outer catheter with a 1.4 mm and 2.0 mm diameter, respectively. After needle placement, the inner needle can be retracted leaving behind a channel for seed delivery. Artifact formation of the needle and piezoelectric motors is evaluated in 3 T MRI. CoBra's needle path planning algorithm is designed for steerable needles based on the seed plan and prostate contours received from MRI. It clusters the seeds and creates candidate path plans to reach all seed positions, using a single insertion point. Fig. 2 is an example of Seed plan using single entry point. Acceptable results are those that allow all seed positions to be reached, while keeping a minimum, pre-set distance between the needle and the urethra. From all candidates, the one that causes the least amount of tissue damage is chosen. For CoBra robot testing, the team developed a BIP, depicting the motion, deformation, and inflammation of the prostate seen clinically during needle insertion. The BIP is connected to SOFA (Simulation open framework architecture) to estimate the changes in the prostate, thus enabling tracking of target position interactively (Fig. 3).

III. Results and discussion

CoBra robot was tested under a 3 T Philips MRI scanner at Lille Hospital, France. The absolute encoder and USM actuators resulted in normal functioning under a strong magnetic field without influencing the MR-images during the scan. Zipper artifacts were detected from an unshielded wire. Straight needle insertions were performed on Agar phantoms and BIP. Tests of the BIP and robotic needle insertion in RT-MRI showed that the BIP features could be extracted with good fidelity and the needle tip identified under T1w, T2w, and DYN-B FEE sequences. sCT datasets showed a mean MAE of 59.8 HU and 65.8 HU for first and second test sites, taking on average 8.5 s to generate a complete sCT (90 slices) on our GPU (Nvidia Quadro P6000). The in-bore testing of the BX-module and steerable needle weren't influenced by the 3 T field.

IV. Conclusions

This abstract reports the recent progress attained from the CoBra project, as it reaches the phantom testing stage. BIP helped the team to replicate the actual condition of prostate intervention, as in biological tissue, the target shifts due to multiple needle insertions and prostate motion. The CoBra concept contributes towards adaptive treatment and diagnosis, which consists of adapting the control of the robot for precise deposition of seed with a minimal perineum entry points and reducing the tissue damage. The novel curved steerable needle concept helps further to reach the targets which are inaccessible to conventional prostate brachytherapy treatment with straight needles in a grid. The CoBra concept, a multi-system model, integrates the robotic sub-systems using MRI-guidance in closed-loop feedback to implement an adaptive MR-robot control, for the treatment of localized cancers.

Figure 1 :

 1 Figure 1: Steerable needle concept under the CoBra project

Figure 2 :

 2 Figure 2: Seed plan using single entry point.

Figure 3 :

 3 Figure 3: (A) CoBra robot-guide, (B) Active prostate phantom and Phantom MRI scan.

ACKNOWLEDGMENTS

We thank our observer partners for assisting us with their components, TPS by Eckert & Ziegler BEBIG, Uni-Lift by NORAS, MR-guidewire by MaRVis, and absolute encoder by Heidenhain-Numeric Jena.

AUTHOR'S STATEMENT

This work has been co-funded by the European regional development Fund under EU Interreg 2 Seas, Grant 2S04-022.