Pierre Develter
email: pierre.develter@onera.fr

Jonathan Bosse
email: jonathan.bosse@onera.fr

Olivier Rabaste
email: olivier.rabaste@onera.fr

Philippe Forster
email: pforster@parisnanterre.fr

Jean-Philippe Ovarlez
email: jean-philippe.ovarlez@onera.fr

Sur

Sur la probabilité de fausse alarme du Filtre Adapté Normalisé pour la détection de cibles hors-grille

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Classiquement, la détection de signaux à paramètres inconnus est réalisée grâce à un test du Rapport de Vraisemblance Généralisé (GLRT) : les paramètres inconnus sont remplacés par leurs estimateurs au sens du maximum de vraisemblance (MV) dans le test du rapport de vraisemblance. Quand des expressions analytiques du MV ne sont pas disponibles pour les paramètres d'intérêt du signal, la plupart des stratégies de détection simplifient le problème en discrétisant l'espace des paramètres en une grille. Cependant, les vrais paramètres n'ont aucune raison de se trouver précisément sur la grille, puisqu'ils sont distribués continûment. Cela crée une désadaptation entre les paramètres sous test et les vrais paramètres de la cible, et détériore les performances en détection de la plupart des tests de l'état de l'art supposant les cibles sur-grille. Dans cet article, nous illustrons ce problème en contexte radar où les paramètres inconnus peuvent inclure le Doppler, le retard ou la direction d'arrivée.

L'impact du hors-grille est particulièrement spectaculaire pour des détecteurs comme le Filtre Adapté Normalisé (NMF) [START_REF] Scharf | Signal detection in Gaussian noise of unknown level : an invariance application[END_REF]. Ce test est utilisé sous hypothèse bruit Gaussien de niveau inconnu. Il est aussi beaucoup utilisé en contexte radar adaptatif sur signaux non Gaussiens, par exemple quand le bruit est distribué selon une distribution symétrique elliptique complexe (CES) [START_REF] Pascal | On a SIRV-CFAR detector with radar experimentations in impulsive noise[END_REF]. Dans certains cas, la probabilité de détection peut descendre à 0 même pour des hauts Rapports Signal sur Bruit (RSB) [START_REF] Rabaste | Off-grid target detection with Normalized Matched Subspace Filter[END_REF], en particulier à basses P F A , usuelles en contexte radar.

Pour résoudre ce problème, la solution la plus évidente consiste à tester tout le support des paramètres continus, et pas juste la grille : c'est le vrai GLRT hors-grille. Cependant, à notre connaissance, l'expression analytique de la statistique de l'hypothèse nulle et sa P F A est inconnue dans la littérature pour ce GLRT : puisqu'il s'agit du maximum d'un ensemble infini de v.a. non indépendantes, le résultat n'est pas trivial.

Quelques travaux proposent des solutions analytiques optimales pour d'autres types de désadaptation comme [START_REF] Besson | Detection of a signal in linear subspace with bounded mismatch[END_REF][START_REF] De Maio | Fractional QCQP with applications in ML steering direction estimation for radar detection[END_REF], mais ils ne sont pas optimaux pour la désadaptation que nous considérons. D'autres approches incluent l'approximation du GLRT hors-grille, avec une approche par sous-espace [START_REF] Rabaste | Off-grid target detection with Normalized Matched Subspace Filter[END_REF] ou une méthode inspirée par l'écartométrie [START_REF] Develter | Off-grid radar target detection with the normalized matched filter : A monopulse-based detection scheme[END_REF]

H 0 : r = n , H 1 : r = α d(θ) + n , (1)
d(θ) = 1 √ N 1, e 2iπθ , . . . , e 2iπ(N -1)θ T . (2)
Ce modèle de vecteur cible est souvent rencontré en radar en problème de détection d'une cible après intégration du retard où plus généralement quand le problème consiste en l'estimation d'une sinusoïde complexe noyée dans du bruit.

On montre que pour σ supposé inconnu sous H 0 et α, σ supposés inconnus sous H 1 , avec θ supposé connu, le test GLRT correspondant est le filtre adapté normalisé (NMF) suivant [START_REF] Scharf | Signal detection in Gaussian noise of unknown level : an invariance application[END_REF] :

d(θ) H Γ -1 r 2 d(θ) H Γ -1 d(θ) r H Γ -1 r H1 ≷ H0 w 2 .
(

) 3
Ce test est aussi utilisé en radar adaptatif en contexte non Gaussien [START_REF] Pascal | On a SIRV-CFAR detector with radar experimentations in impulsive noise[END_REF][START_REF] Conte | Asymptotically optimum radar detection in compound-Gaussian clutter[END_REF], par exemple quand le bruit est distribué selon une distribution symétrique elliptique complexe (CES) . Sa statistique, dans ce cas, est la même que dans le cas Gaussien. De manière équivalente, l'équation (3) peut être réécrite avec les vecteurs blanchis normalisés :

s(θ) H u Γ -1/2 d(θ) et u = Γ -1/2 r Γ -1/2 r .
La relation P F A -seuil correspondante est connue :

P F A = 1 -w 2 N -1 . (5)
Quand le NMF est testé en un point θ différent du vrai paramètre θ 0 de la cible, on dit que la cible est hors-grille. Cela induit une désadaptation (θ ̸ = θ 0 ) entre le vrai vecteur cible de la cible s(θ) et le vecteur cible s(θ 0 ) sous test. Malheureusement, il a été montré dans [START_REF] Rabaste | Geometrical design of radar detectors in moderately impulsive noise[END_REF] que le NMF est très sensible à une désadaptation du vecteur testé, ce qui peut potentiellement amener à une grave détérioration de la performance en détection : en particulier, pour une désadaptation plus grande que le seuil de détection, la probabilité de détection asymptotique tend vers 0 pour des RSB élevés : ce phénomène apparaît pour des P F A aussi élevées que 10 -3 dans la cellule de résolution de largeur 1/N , car la réponse du NMF tombe en dessous du seuil en bord de cellule [START_REF] Rabaste | Off-grid target detection with Normalized Matched Subspace Filter[END_REF]. Pour corriger ce problème, on doit donc estimer θ.

Avec θ inconnu, la procédure naturelle du GLRT mène à

GLRT (u, D) = max θ∈D s(θ) H u 2 H1 ≷ H0 w 2 , (6)

Une relation P F A -seuil analytique avec une interprétation géométrique

Hotelling [START_REF] Hotelling | Tubes and spheres in n-spaces, and a class of statistical problems[END_REF] a établi une méthodologie pour l'étude de tests statistiques sur la sphère réelle. Elle peut être utilisée pour évaluer la P F A du GLRT dans le cas réel : cette approche est présentée en section 3.1. Cette approche ne peut pas être transposée directement dans le cas complexe de (4), et nous étudions une extension dans la section 3.2.

L'approche géométrique de Hotelling dans le cas réel

Dans cette section, nous étudions tout d'abord le cas réel. Notons que l'expression du NMF (4) a une interprétation géométrique simple. C'est en effet le cosinus carré de l'angle entre le vecteur cible s (θ) et le signal reçu u. Le seuil peut être vu comme le cosinus carré de l'angle cos -1 w. Quand l'angle entre les vecteurs est sous cet angle limite, une cible est détectée. Par ailleurs, puisque u a été blanchi, il est uniformément θ comme l'union des deux calottes précédentes. On a une fausse alarme dès que u ∈ SC θ . En conséquence, la P F A est le rapport de la surface des calottes sphériques sur la surface de la sphère unité.

Considérons le GLRT dans l'équation [START_REF] Develter | Off-grid radar target detection with the normalized matched filter : A monopulse-based detection scheme[END_REF]. Une fausse alarme survient quand u ∈ T = θ∈D SC θ . Alors la P F A est le rapport entre la surface de T et la surface de la sphère unité. Un exemple dans R 3 est représenté en figure 1 sur S 2 . T est l'union de deux tubes symétriques plongés dans la sphère, de rayon géodésique cos -1 w autour des deux courbes s(θ) et -s(θ) avec des demi-calottes sphériques ajoutées aux extrémités.

Calculer la P F A revient à trouver la surface de l'ensemble des points dans S N -1 qui sont à la distance Euclidienne donnée par 2(1 -w) des courbes.

Dans [START_REF] Hotelling | Tubes and spheres in n-spaces, and a class of statistical problems[END_REF], Hotelling établit une formule pour la surface des tubes sur la sphère réelle à n dimensions S n-1 .

Extension de l'approche d'Hotelling au cas complexe

Malheureusement le résultat d'Hotelling n'est pas directement applicable au GLRT considéré. Comme nous allons le voir, quand on transpose notre problème complexe dans R 2N la variété résultante est alors à deux dimensions, car la valeur absolue dans le terme droit du GLRT donné en équation (6) implique alors une recherche du maximum de la phase.

En effet, pour n'importe quel scalaire réel α ∈ [0, 2π], remarquons que Re(s (θ)

H ue -iα ) ≤ s (θ) H u , ces deux quan-tités étant égales pour α = s (θ) H u. On trouve alors, en décomposant s (θ) = s r (θ) + is i (θ) et u = u r + iu i en parties réelles et complexes :

Re(s (θ) H ue -iα ) = s r (θ) T u r + s i (θ) T u i cos α = γ 1 (θ) T u cos α + γ 2 (θ) T u sin α , où γ 1 (θ) = s r (θ) s i (θ) , γ 2 (θ) = -s i (θ) s r (θ) et u = u r u i est un
vecteur de bruit réel de dimension 2N tiré uniformément sur S 2N -1 sous H 0 . Le GLRT (6) peut alors être écrit comme :

max α,θ γ 1 (θ) T u cos α + γ 2 (θ) T u sin α H1 ≷ H0 w . ( 7
)
La nature bi-dimensionnelle de la variété réelle à inspecter est directement apparente avec cette formulation : en notant γ(α, θ) = γ 1 (θ) cos α + γ 2 (θ) sin α, la région d'acceptation dans le cas complexe est un nouveau tube T autour de la variété à deux dimensions γ(θ, α).

Le résultat d'Hotelling ne couvre pas ce type de variété puisqu'il ne donne la surface d'un tube qu'autour d'une courbe. Cependant, dans [START_REF] Johnstone | On Hotelling's formula for the volume of tubes and Naiman's inequality[END_REF], ce résultat est étendu au cas spécial de variété bi-dimensionnelles d'intérêt pour nous.

Théorème 3.1 [10] Pour i ∈ [1, 2], soit γ i : [0, t 0 ] → S n-1 des courbes régulières. Supposons γ 1 (t) T γ 2 (t) = 0 pour tout t. Soit Z(t) = γ 1 (t) T u 2 + γ 2 (t) T u 2 1/2
où u est distribué uniformément sur S n-1 . Alors pour 0 < w < 1 :

P max 0≤t≤t0 Z (t) > w ≤ (1 -w 2 ) (n-2)/2 + Γ n 2 w (1 -w 2 ) (n-3)/2 2π 3/2 Γ n -1 2 × t0 0 2π 0 ∥ γ1 (t) cos ω + γ2 (t) sin ω∥ 2 - γ1 (t) T γ 2 (t) 2 1/2 dω dt , (8)
où γi (t) est la dérivée de γ i (t) par rapport à t. Quand il n'y a pas de recouvrement, cette inégalité devient une égalité (par recouvrement, on entend une intersection non nulle des sections droites de la variété dans les régions de décisions de l'hypothèse H 1 sur S n-1 ).

Il apparaît que nous pouvons reformuler le problème de manière à remplir les conditions du théorème. Avec nos notations, on peut vérifier que

s (θ) H u 2 = γ 1 (θ) T u 2 + γ 2 (θ) T u 2 , (9)
de sorte que le théorème nous donne la P F A désirée (quand il y a égalité). On obtient alors notre résultat principal : -Sous l'hypothèse d'un bruit blanc (Γ = σ 2 I), cette équation se simplifie en :

P F A = (1 -w 2 ) N -1 + (10) π 3 Γ(N ) w (1 -w 2 ) N -3 2 Γ N -1 2 N 2 -1 1 2 (θ 2 -θ 1 ) .
-Sous l'hypothèse d'un bruit coloré (Γ ̸ = σ 2 I), (8) peut être évaluée numériquement.

Le premier terme dans [START_REF] Johnstone | On Hotelling's formula for the volume of tubes and Naiman's inequality[END_REF] représente la surface des deux calottes sphériques aux extrémités du tube. Ainsi, il est égal à la P F A connue du NMF exprimée en équation [START_REF] De Maio | Fractional QCQP with applications in ML steering direction estimation for radar detection[END_REF]. Le second terme montre l'influence de la variété induite par la nature horsgrille du problème.

Résultats numériques

Vérifions la validité de la formule [START_REF] Johnstone | On Hotelling's formula for the volume of tubes and Naiman's inequality[END_REF]. La figure 2 présente la relation P F A -seuil donnée dans l'équation [START_REF] Johnstone | On Hotelling's formula for the volume of tubes and Naiman's inequality[END_REF] et les seuils calculés numériquement en utilisant 10 8 échantillons de bruit blanc Gaussien pour un vecteur cible de taille N = 10. La recherche continue sur le domaine D est remplacée par une recherche discrète en suréchantillonnant la case d'un facteur 30, où D = [0, 1/N ] est une des cellules de résolution usuelle de Fourier.

La formule correspond parfaitement quand la P F A est assez faible (c'est à dire quand le seuil est assez haut). La formule n'est pas valide quand la P F A est proche de 1 à cause de phénomènes de recouvrement de surface (sa valeur dépasse même 1). Cependant, des P F A aussi élevées n'ont pas d'intérêt pratiques pour les applications les plus communes en radar. Une analyse du phénomène de recouvrement sera donnée dans un prochain article.

Conclusion

Dans cet article, nous avons répondu au problème de la détection hors-grille en utilisant le NMF-GLRT en élaborant une relation P F A -seuil analytique, valide pour la plupart des applications usuelles notamment en radar. Dans de futurs travaux, nous chercherons mathématiquement les conditions garantissant qu'il n'y a pas de recouvrement pour une P F A donnée pour fournir un domaine et des conditions de validité précises.

Références

FIGURE 1 -

 1 FIGURE 1 -T plongé sur la sphère unité S 2 dans R 3 composé d'un tube (en violet) autour d'une variété à une dimension θ∈D s(θ), deux calottes demi-sphériques (en vert) et leurs opposés dont le contour est en pointillés. SC θ est dessiné en bleu.

Corollaire 3 . 1 . 1

 311 En l'absence de recouvrement (régimes à faibles P F A ), la P F A du GLRT (6) pour un domaine de recherche D = [θ 1 , θ 2 ] avec le vecteur cible d (θ) défini par (2) est donnée par l'équation[START_REF] Rabaste | Geometrical design of radar detectors in moderately impulsive noise[END_REF].

3 FIGURE 2 -

 32 FIGURE 2 -Comparaison entre la relation P F A -seuil théorique de l'équation (10) et la relation P F A -seuil empirique pour N = 10. La relation hors-grille (5) est aussi tracée à des fins de comparaison. Tests pour θ ∈ [0, 1/N ].

  Notations : Les matrices sont en gras et en majuscules, les vecteurs en gras. Pour toute matrice A ou vecteur, A T est la transposée de A et A H est la transposée Hermitienne de A. I est la matrice identité et CN (µ, Γ) est la distribution normale complexe circulaire de moyenne µ et de matrice de covariance Γ. S n-1 est la sphère unité dans R n . L'opérateur partie réelle d'un nombre complexe est dénoté par Re(.). L'opérateur u est l'angle du nombre complexe u. La mission principale d'un radar consiste à détecter la présence d'une cible d'amplitude complexe d ∈ C N perturbée par du bruit additif n. Ce problème peut être formalisé par un test d'hypothèse binaire :

	grille et la formulation du vrai GLRT. La section 3 présente le
	calcul de la relation P F A -seuil. Dans la section 4, la validité du
	résultat présenté est vérifiée sur simulation numérique.
	2 Formulation du problème
	, par exemple. Ces solu-
	tions sont sous-optimales : elles sont conçues pour des situa-
	tions pratiques où l'utilisation d'une approximation précise du
	GLRT est irréaliste car son manque de forme analytique peut
	mener à un coût de calcul élevé.
	Dans cet article, nous obtenons une relation P F A -seuil pour
	le NMF hors-grille grâce à une analyse géométrique du problè-
	me. Cette relation est valide pour de faibles P F A d'intérêt dans
	des applications courantes en radar.
	La section 2 présente le modèle du signal, le problème hors-

H1 ≷ H0 w 2 , (4) où s(θ) = Γ -1/2 d(θ)


Priscilla Galinié

Camille Eyssartier

Christophe Sauret

Mickael Tordjman

Marie-Line Pissonier

Robert Carlier

Wafa Skalli

Claudio Vergari
email: c.vergari@gmail.com

In-vivo characterization of the lumbar annulus fibrosus in adults with ultrasonography and shear wave elastography

Keywords:

à la diffusion de documents scientifiques 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.

Introduction

Low back pain (LBP) is the leading cause of disability in the world [START_REF] Hartvigsen | What low back pain is and why we need to pay attention[END_REF]. While pain is the main symptom, and the cause is often difficult to diagnose, the intervertebral disc (IVD) and its degeneration often play a role. This is not surprising since the IVD is a key component of spine's biomechanics, and at the same time it is avascular and therefore it has low endogenous repair capability. The degeneration usually occurs progressively with aging, and is characterized by a dehydration of the tissue, which can become more fibrotic and tend to lose its microarchitecture [START_REF] Galbusera | Ageing and degenerative changes of the intervertebral disc and their impact on spinal flexibility[END_REF][START_REF] Cao | 3D characterization of morphological changes in the intervertebral disc and endplate during aging: a propagation phase contrast synchrotron microtomography study[END_REF]. This also results in a stiffening of the spinal functional unit, which is due to a combination of disc's height loss and changes in mechanical properties.

Magnetic resonance imaging (MRI) is the gold standard to assess the disc. It can provide information on disc morphology and hydration through standard imaging protocols, while more advanced methods are being developed to estimate mechanical properties through shear wave elastography [START_REF] Co | Magnetic resonance elastography of intervertebral discs: spin-echo echo-planar imaging sequence validation[END_REF][START_REF] Cortes | Magnetic resonance elastography of intervertebral disc -a new biomarker for disc degeneration[END_REF]. However, MRI is difficult of access in clinical routine, and it might not be adapted for early detection of disc alteration or long-term follow-up of disc pathology.

Ultrasound shear wave elastography (SWE) has shown promising results for the non-invasive assessment of disc in vivo. It is based on the measurement of shear wave speed (SWS) in the annulus fibrosus (AF), which is directly related to the tissue's mechanical properties [START_REF] De Reuver | Ultrasound shear wave elastography of the intervertebral disc and idiopathic scoliosis: a systematic review[END_REF][START_REF] Gennisson | Ultrasound elastography: principles and techniques[END_REF], and it can be coupled with conventional ultrasonography which can provide information on the AF microstructure [START_REF] Kakitsubata | Sonographic characterization of the lumbar intervertebral disk with anatomic correlation and histopathologic findings[END_REF][START_REF] Langlais | Microstructural characterization of annulus fibrosus by ultrasonography: a feasibility study with an in vivo and in vitro approach[END_REF][START_REF] Marshburn | New Heights in Ultrasound: first Report of Spinal Ultrasound from the International Space Station[END_REF]. SWE has been applied to List of abbreviations: AF, annulus fibrosus; BMI, body mass index; IVD, Intervertebral disc; MRI, Magnetic resonance imaging; NP, Nucleus pulposus; SWE, Shear wave elastography; SWS, Shear wave speed. lumbar discs in children, demonstrating a stiffening of the AF in adolescent idiopathic scoliosis and recovery after arthrodesis [START_REF] Vergari | Shear wave elastography of lumbar annulus fibrosus in adolescent idiopathic scoliosis before and after surgical intervention[END_REF].

In adults, ultrasound measurements were limited to the cervical discs [START_REF] Vergari | Non-invasive biomechanical characterization of intervertebral discs by shear wave ultrasound elastography: a feasibility study[END_REF], while lumbar discs were neglected because of the technical challenges posed by their depth within the abdominal cavity and the nature of the abdominal tissues that ultrasound must traverse. However, lumbar IVDs in adults present a specific clinical interest because they are often subject to degeneration with aging. The present work aimed to assess the feasibility and reliability of SWE and conventional ultrasonography on adult lumbar discs.

Material and methods

Participants

Thirty asymptomatic adults were included (14 women and 16 men, age range: 22 to 67 years old). Subjects had an average body mass index (BMI) of 24 ± 4 kg/m 2 , with values ranging between 17 (considered underweight) and 31 kg/m 2 (considered obese). The exclusion criteria were the following: any musculoskeletal pathology, including back pain, previous spine surgery, or pregnancy at the time of measurement. All the participants signed an informed consent, as approved by the ethical committee CPP Ile de France VI (authorizations N. 2018-A01926-49 and N. 6001), and the study was conducted in accordance with the principles embodied in the Declaration of Helsinki and in accordance with local statutory requirements.

IVD imaging and protocol

A Supersonic Mach 30 (Supersonic Imagine, France) was used for the simultaneous conventional ultrasonography and SWE acquisitions. The device provides a direct measurement of SWS in the tissue, in the form of a colored chart (Fig. 2). SWS is related to the tissue shear modulus (µ) by the following relationship: SWS = ̅̅̅̅̅̅̅ ̅ μ/ρ √ , where ρ is tissue mass density (~1000 kg/m 3 ) [START_REF] Gennisson | Ultrasound elastography: principles and techniques[END_REF]. It was previously shown in vitro that SWS of the annulus is related to disc's mechanical properties [START_REF] Vergari | Intervertebral disc characterization by shear wave elastography: an in vitro preliminary study[END_REF].

Data were acquired by one operator in L3-L4 intervertebral disc, using the "general" acquisition mode. Subjects were in supine position, arms along their trunk (Fig. 1), and they were asked to maintain a steady breathing during the acquisition. The aortic bifurcation into iliac arteries was used as an anatomical landmark to locate the L4 level. The probe was tilted to be aligned with the IVD plane, and it was considered aligned when lamellae were visible. At least 3 clips of approximately 10 s were acquired for each subject, for a total of approximately 30 frames containing the elastography and ultrasonography panels (Fig. 2). A L18-5 probe was first used in all patients and, if the signal was too attenuated (i.e., bad quality of the elastographic chart), a lower-frequency L10-2 probe was used. All SWS measurements were made in "penetration mode", while depth, brightness, size and location of the region of interest were adapted to the subject.

Shear wave elastography data processing

SWE data were analyzed using a previously described protocol and software [START_REF] Vergari | Lumbar annulus fibrosus biomechanical characterization in healthy children by ultrasound shear wave elastography[END_REF]. In brief, the acquired clips were exported and transformed into single frames and custom software was used to calculate the average SWS in the AF in each frame. A region of interest was selected in the elastographic chart and it was semi-automatically tracked in the sequence of images, in order to ensure that the same region was tracked in all frames. The region of interest was manually placed in the AF, where the lamellae were visible. The average SWS was calculated by converting the colour red-green-blue values of each pixel in the region of interest to speed, using the colour chart available in each clip (Fig. 2). Then, an average was calculated for each movie (~10 images), and finally an overall average SWS per subject. Images or regions showing signal saturation (dark red or dark blue pixels) or with missing signal were excluded.

Ultrasonography data processing

Lamellar thickness was measured using a previously described technique [START_REF] Langlais | Microstructural characterization of annulus fibrosus by ultrasonography: a feasibility study with an in vivo and in vitro approach[END_REF]. Every image was processed with a histogram equalization filter to improve the contrast. Then, an operator traced three segments over the annulus in a radial direction, starting from the AF outer border to the AF inner border. Gray values along each segment were plot, and local maxima and minima were detected automatically. Midpoints between each consecutive local maximums and minimums were computed, and the lamellar thickness was calculated as the distance in pixels between two consecutive midpoints. Pixels were transformed in micrometers using the scale available on the side of each image Fig. 2).

Reliability and statistics

Two additional operators, one with a 7-year experience and the other Data distributions of lamellar thickness and SWS were not normal according to Shapiro-Wilk normality test (p < 0.05), and therefore, correlations between SWS, age and BMI were assessed using Spearman's test. Statistical significance was set at α = 0.05. Results are reports as average ± standard deviation. All processing was carried out in Matlab (The MathWorks Inc., Natick, MA, USA).

Results

Shear wave elastography

One subject was excluded from further analysis because of a previously undiscovered scoliosis, while four subjects were excluded because their disc could not be imaged properly. Demographic data and the main results of the remaining 25 subjects was reported in Table 1. The subjects excluded due to imaging issues had a BMI between 24 and 30 kg/m 2 , with discs at depths between 3 and 6 cm. Average disc depth was 3.7 ± 1 cm, with depths ranging from 2.5 to 6.1 cm. Disc depth was correlated with BMI (R = 0.7, p = 0.0007).

Average SWS in AF at the L3-L4 level was 4.0 ± 0.9 m/s. SWS was significantly correlated with BMI (R = 0.4, p = 0.04) but not with age (p = 0.6). However, as shown in Fig. 3, both correlations were significant when pooled with previously published data on children [START_REF] Vergari | Shear wave elastography of lumbar annulus fibrosus in adolescent idiopathic scoliosis before and after surgical intervention[END_REF]. Gender did not have a significant effect on SWS (p = 0.7). Inter-operator reproducibility of SWS was 8.7% and intra-operator repeatability was 5.6%.

Ultrasonography

Eight subjects were excluded from this analysis because it was not possible to clearly image the lamellae. On the remaining twenty-two subjects, average lamellar thickness was 255 ± 27 µm, and it was correlated with SWS (R = 0.77, p < 0.0001). Inter-operator reproducibility of lamellar thickness measurement was 24.7 µm (9.6%) and intraoperator repeatability was 19.6 µm (7.6%). Average pixel size of the ultrasound images was 71 µm.

No significant effect of age or BMI were observed (p > 0.05), but both correlations were significant when values were pooled with previous data on children (Fig. 4).

Discussion

In this study, ultrasonography and SWE were applied for the first time to adult lumbar disc to assess the structural and mechanical properties of AF. Ultrasonography gave access to the lamellar thickness while SWE measured the SWS, which is related to tissue mechanical properties. While the theoretical relationship between SWS and me-

chanical properties is straightforward (SWS = ̅̅̅̅̅̅̅ ̅ μ/ρ √
), it relies on several hypotheses which are not necessarily valid in the annulus fibrosus (homogeneity, infinite medium, etc.). For this reason, in this work only SWS was provided and not its corresponding shear modulus or, under the hypothesis of incompressibility, its Young's modulus E (E ≈ 3µ) [START_REF] Gennisson | Ultrasound elastography: principles and techniques[END_REF]. Therefore, SWS should be considered as a biomarker related to the intervertebral disc mechanical properties.

Reproducibility of SWE was 8.7%, with negligible operator effect, which is only slightly higher than uncertainties reported for applications in children (6.2%, [START_REF] Vergari | Lumbar annulus fibrosus biomechanical characterization in healthy children by ultrasound shear wave elastography[END_REF]). SWE measurement was not possible in 4 subjects (13% of the cohort), which is not negligible. It is our opinion that the challenge in SWE of the disc in adults only partially depends on disc depth in the abdomen, but tissue characteristics and ultrasound absorption also play a role. These subjects had similar BMI and disc depths to the rest of the cohort. In other words, it appears that measurement is not feasible in one out of four subjects with BMI > 24 kg/m 2 . Age does not seem to play a factor in feasibility, since these subjects' ages ranged between 27 and 66 years, nor does disc depth. It is possible that lower frequency probes (1-6 MHz) could allow SWE assessment of the annulus fibrosus, but resolutions obtained at this frequency might not be sufficient to characterize the lamellae.

Measurement of average lamellar thickness was not feasible in six subjects (20% of the cohort). Three of them coincided with the patients for which SWS measurement was not feasible, suggesting a limitation of ultrasound penetration. However, the remaining three subjects had low to normal BMI (between 17 and 26 kg/m 2 ). While the lack of visibility of the lamellae could be still due an insufficient ultrasound penetration, it could also be due to a loss of architecture of the AF due to some stage of disc degeneration, and therefore a loss of lamellar organization. This would be consistent with the fact that the IVD tends to slowly and progressively degenerate with age starting during childhood [START_REF] Wu | Intervertebral Disc Aging, Degeneration, and Associated Potential Molecular Mechanisms[END_REF], but this hypothesis should be confirmed with a study comparing ultrasound and MRI measurements in healthy and degenerated discs. Previous works also found that average SWS for children and adolescents was 3.0 ± 0.3 m/s [START_REF] Langlais | Shear-wave elastography can evaluate annulus fibrosus alteration in adolescent scoliosis[END_REF] whereas an average of 4.0 ± 0.8 m/s was found in the present study including adults from 22 to 67 years old. Fig. 3 confirms that this difference is due to the different age of the two cohorts. It also shows that aging affects SWS in the annulus fibrosus, corresponding to an increase of tissue stiffness with age. This is consistent with the known natural stiffening of the spine occurring with aging [START_REF] Galbusera | Ageing and degenerative changes of the intervertebral disc and their impact on spinal flexibility[END_REF][START_REF] Wu | Intervertebral Disc Aging, Degeneration, and Associated Potential Molecular Mechanisms[END_REF]. However, this aging effect observed with SWS could also include, or hide, a component of progressive but asymptomatic degeneration of the disc. For instance, previous studies proved that pathologic discs in adolescent idiopathic scoliosis are stiffer than healthy ones [START_REF] Langlais | Shear-wave elastography can evaluate annulus fibrosus alteration in adolescent scoliosis[END_REF]. This hypothesis is also corroborated by the increase in variability that can be observed with age, suggesting that older discs are more affected by their history. However, this is still speculative at this stage, and it is impossible to state if some of the discs were pathological, and in particular the ones showing extreme SWS values observed, for instance higher than 5 m/s. Therefore, further studies should clarify the relevance of SWS as an early biomarker of disc degeneration by comparison with other imaging methods.

A previous study found and average lamellar thickness of 229.7 ± 91.5 µm for subjects aged from 10 to 16 years old [START_REF] Langlais | Microstructural characterization of annulus fibrosus by ultrasonography: a feasibility study with an in vivo and in vitro approach[END_REF], which is similar to the value obtained in this study (257.4 ± 27 µm). However, Fig. 4 suggests that ages do play a role in the lamellar organization of AF, as previously suggested by Marchand and Ahmed [START_REF] Marchand | Investigation of the laminate structure of lumbar disc anulus fibrosus[END_REF]. Lamellar thickness was also correlated with SWS (Fig. 4); this could be expected, as it confirms that AF microstructure plays a role in its macroscopic mechanical properties. It is possible that thicker lamellae correspond to more densely packed collagen fibers, and hence a stiffer tissue which would show higher SWS.

The main limitation of the proposed technique is that measurement was not possible in some subjects. This should be further investigated, to ascertain if only disc depth and the properties of the overlying tissue limit ultrasound penetration, or if the actual IVDs presented properties that impeded ultrasound penetration. In the latter case, these discs might be the ones which would benefit more from further clinical assessment, and the impossibility of acquiring SWE images could represent a biomarker itself.

Another limitation is the uncertainty concerning the measured spine level, since the aortic bifurcation corresponds to the L4 level in only 64% of the population. Spine radiographs can in some cases help in determining the location of the discs, for instance relative to the iliac crests.

Conclusions

This study showed that biomechanical and structural characterization of adult AF with shear wave elastography and ultrasonography is feasible and reliable. SWS and lamellar thickness in AF increased with age, which might be due to disc aging or slow progression of disc degeneration; this aspect requires further investigation.

SWE and ultrasonography are non-invasive and radiation-free methods, easy to access in clinical routine, and therefore they show great potential to obtain biomarkers of disc health.

Fig. 1 .

 1 Fig. 1. -Setup for measurements in the lumbar intervertebral disc of a subject.

Fig. 2 .

 2 Fig. 2. -Example of data obtained with the shear wave elastography (left). The lamellae are visible in the region of interest in the annulus fibrosus (AF). The nucleus pulposus (NP) appears hypoechogenic.

Fig. 3 .

 3 Fig. 3. -Correlation between shear wave speed at the L3-L4 level and subjects' age and body mass index. Data of asymptomatic children from (Vergari et al., 2020) were included.

Fig. 4 .

 4 Fig. 4. -Correlation between lamellar thickness at L3-L4 level and subjects' body mass index, age and disc shear wave speed.

Table 1

 1 Demographic data.

	Sex	14 women and 16 men
	Age (years)	36 ± 16.5
	Weight (m)	68.8 ± 13.2
	Height (kg)	1.7 ± 0.1
	Body mass index (kg/m 2 )	23.2 ± 3.4
	Disc depth (cm)	3.7 ± 0.9
	Annulus shear wave speed (m/s)	4.0 ± 0.9
	Lamellar thickness (µm)	255.01 ± 26.9

Funding

The authors are grateful to the BiomecAM chair program on subjectspecific musculoskeletal modelling (with the support of ParisTech and Yves Cotrel Foundations, Société Générale, Covea and Proteor). The funding body had no involvement in study design, in the collection, analysis and interpretation of data; in the writing of the report; nor in the decision to submit the article for publication.

Ethical approval

All the participants signed an informed consent, as approved by the ethical committee CPP Ile de France VI (authorizations N. 2018-A01926-49 and N. 6001), and the study was conducted in accordance with the principles embodied in the Declaration of Helsinki and in accordance with local statutory requirements.

CRediT authorship contribution statement

Declaration of Competing Interest

No conflict of interest to disclose.


M Benchekroun
email: mouna.benchekroun02@gmail.com

B Chevallier
email: baptiste.chevallier7@gmail.com

V Zalc

D Istrate

D Lenne

N Vera

The Impact of Missing Data on Heart Rate Variability Features: A Comparative Study of Interpolation Methods for Ambulatory Health Monitoring

Keywords: Heart Rate Variability, HRV analysis, real time, Inter beat intervals, IBI, RR intervals, wearables, e-health

Wearable devices have facilitated the remote measurement of heart rate variability (HRV), a promising indicator of various physiological and psychological states such as stress, sleep and other conditions. However, errors during the transmission or acquisition can lead to missing data, which can affect HRV features and cause false medical diagnosis. Interpolation is a useful technique for handling missing data, but the choice of interpolation method must be carefully considered. Therefore, the objective of this study is to investigate the impact of four interpolation methods (Nearest Neighbour, Linear, Shape-preserving piecewise cubic Hermite, and cubic spline) on HRV features when data is deliberately deleted. It is an expansion of a previously published work on HRV data imputation. The study utilizes a real-time approach to data interpolation and HRV analysis. The results indicate that the choice of interpolation method significantly affects HRV features, with varying effects depending on the percentage of missing data. Additionally, the study proposes to adapt the interpolation method based on both the percentage of missing values and the targeted HRV feature for maximum performance.

Introduction

With the rise of telemedicine and healthcare wearables, scientists are eager to collect every trackable parameter from the human body throughout different physiological signals. One widely used signal is the Heart rate variability (HRV), now used as an indicator of different physiological states and pathologies [START_REF] Acharya | Heart rate variability: a review[END_REF]. Its time and frequency domain analysis can give insights into autonomic nervous function. They provide information about the sympathetic-parasympathetic balance and cardio-vascular health [START_REF] Cardiology | Heart rate variability: standards of measurement, physiological interpretation and clinical use[END_REF].

HRV measures the variation in the time interval between two consecutive heartbeats, known as inter beat intervals (IBI) or RR intervals. They correspond to the time elapsed between two successive R-waves of the QRS complex, characterizing ventricular depolarization, on an ECG signal.

In an ideal situation, HRV analysis is performed with RR interval time series including only pure sinus beats, normally recorded by a 12 lead ECG. However, RR intervals are now usually measured thanks to wearable ECGs or photoplethysmographs (PPG) as a substitute of the gold standard ECG used in hospitals.

Thanks to such wearables, it is now possible to passively record heart activity continuously, opening the way to easier remote health monitoring during user's daily life.

However, for a reliable HRV analysis, these RR time-series should be carefully edited to identify gaps and abnormal heart beats beforehand.

In this paper, we investigate the impacts of data imputation using different interpolation methods on HRV features. We remove an increasing amount of data from an originally perfect HRV signal. The deleted values are then handled by four interpolation methods (Nearest Neighbour, Linear, Shape-preserving piece-wise cubic Hermite and cubic spline). Finally, an estimation error was computed to compare HRV features from reconstructed signals against those computed from the original signal. The goal is to identify the best interpolation method, that yields the lowest error in both time and frequency domains, based on the signal's quality and percentage of missing values. Ultimately, the best approach may be to choose the interpolation method according to the percentage of missing data in each HRV window analysis independently.

Context

The main downside to HRV assessment through wearables is the data quality that is often cor-rupted. Errors occur during the acquisition, the transmission or the storage, thus leading to an important data loss and unintended changes to the original HRV signal. Ectopic beats also introduce a bias into HRV features. When they are not caused by a physiological phenomenon such as premature ventricular contractions (PVC) or premature atrial contractions (PAC), they can occur due to a false QRS detection on the ECG signal or a missed beat.

Such artifacts represent a significant problem in the interpretation of HRV features making it sometimes even impossible. Therefore, they need to be addressed beforehand for a reliable HRV analysis [START_REF] Kim | Effect of missing RR-interval data on heart rate variability analysis in the time domain[END_REF].

Previous studies on the subject suggested different pre-processing methods for RR time series including filtering, deletion and interpolation. Each of these solutions however has its own disadvantages.

The main issue with the deletion approach is the signal depletion since the ectopic beats are removed without being replaced. The remaining RRintervals are just merged together which increases the abrupt changes in the beat to beat variability and the disruptions in the natural fluctuation [START_REF] Peltola | Role of editing of rr intervals in the analysis of heart rate variability[END_REF].

Deletion may be enough for HRV analysis in the time domain but is not sufficient for frequency domain analysis. Re-sampling, which is essential for analyzing HRV in the frequency domain, may produce outliers if the RRI time series contain missing values. [START_REF] Eguchi | RR interval outlier processing for heart rate variability analysis using wearable ECG devices[END_REF] Interpolation on the other hand roughly preserves the overall recording duration and the number of beats, but the beat manipulation does introduce changes that affect HRV analysis. The most used methods are Spline and Linear interpolation. Although they do not significantly affect the power spectral density (PSD), they may produce RRI outliers due to oscillation of the interpolation function, especially when using the spline function [START_REF] Eguchi | RR interval outlier processing for heart rate variability analysis using wearable ECG devices[END_REF].

Besides, authors in [START_REF] Clifford | Quantifying errors in spectral estimates of hrv due to beat replacement and resampling[END_REF], as well as many others, found that interpolation introduces low frequency components (LF) and reduces high-frequency com-Figure 1: Processing approach for HRV analysis ponents (HF) power, thus altering frequency domain HRV features.

Conventional HRV analysis is usually performed in four steps [START_REF] Acharya | Heart rate variability: a review[END_REF], including both deletion of ectopic beats and interpolation of missing values. If HRV is derived from an ECG sensor, the first step is R wave extraction from ECG signal and RR interval computation. Secondly, RR interval evaluation. The purpose of this step is to exclude RR intervals that do not meet some physiological criteria. The deleted RRIs are then replaced by interpolation. Time domain features can be computed from the output of the first two steps, whereas frequency domain features require data re-sampling for spectrum analysis. Another commonly used method is to compute time domain features without interpolation as done in [START_REF] Eguchi | RR interval outlier processing for heart rate variability analysis using wearable ECG devices[END_REF].

There are other studies that propose more advanced methods for imputing missing HRV data [START_REF] Benchekroun | Preprocessing methods for ambulatory hrv analysis based on hrv distribution, variability and characteristics (dvc)[END_REF][START_REF] Stankus | A novel adaptive noise elimination algorithm in long rr interval sequences for heart rate variability analysis[END_REF][START_REF] Zhang | A novel machine learning-enabled framework for instantaneous heart rate monitoring from motionartifact-corrupted electrocardiogram signals[END_REF], however, these methods tend to be quite complex and computationally expensive. They could include sophisticated mathematical algorithms, machine learning models, and other complex techniques that require significant computational resources and expertise to implement. Despite their added complexity, these methods have the potential to provide more accurate and reliable results in comparison to simpler methods. Nevertheless, the use of these methods may not be practical in certain real-world scenarios, such as in resource-limited settings or real-time analysis.

Paper contribution. The particularity of the present paper is the real-time approach to data deletion and interpolation. Both steps are performed iteratively in order to simulate a real-time HRV data acquisition and processing scenario.

This study is one of the first to investigate the effects of missing data on real-time HRV analysis. Data acquisition with missing values is simulated, and the missing values are replaced in real-time using various interpolation techniques before HRV analysis. The HRV features derived from the reconstructed signal are then compared to those from the original RR time-series.

The main purpose is to identify the best approach for processing the RR time-series in real time, based on the percentage of missing data in each HRV window. The real-time aspect is actually vital for continuous health monitoring.

Besides, to the best of our knowledge, this would be one of the first papers to investigate the effect of a very large amount of missing data (up to 70%) on HRV analysis. Recent developments in wearable devices have heightened the need for such studies since wearables produce a huge number of abnormal beats due to motion artifacts as well as missing data due to connectivity problems.

Methods

HRV signals used in this study are considered ideal thanks to the automatic R peak detection, on ECG signals, which was manually corrected by a specialist. These signals did not contain missing nor ectopic peaks. The first step was to delete values from these signals in order to simulate a reallife, ambulatory low quality data acquisition. Since HRV features are usually computed from windows of the signal, data was deleted independently and iteratively from each HRV analysis window. The exact procedure is explained below and depicted in figure 3. Next, the deleted values are replaced using different interpolation methods. In the same way as for data deletion, interpolation is also done iteratively for each window analysis. After the signal has been degraded and reconstructed, time and frequency domain features are computed from both original and reconstructed signals; and compared using the Mean Absolute Percentage Error (MAPE). All the steps were implemented on Python. Figure 2 summarizes the overall process of the study and each step is further detailed in the following sections.

Dataset

The dataset used is from the MIT-BIH Normal Sinus Rhythm RR Interval Database (nsr2db) available on PhysioNet [START_REF] Goldberger | Physiobank, physiotoolkit, and physionet: Components of a new research resource for complex physiologic signals[END_REF].

The database includes beat annotation files for long-term ECG recordings of 54 subjects in normal sinus rhythm (30 men, aged 28.5 to 76, and 24 women, aged 58 to 73). The original ECG recordings were digitized at 128Hz, and the beat annotations were obtained by automated analysis with manual review and correction [START_REF] Goldberger | Physiobank, physiotoolkit, and physionet: Components of a new research resource for complex physiologic signals[END_REF]. In this paper, RR segments including only normal beats between 0.3s and 1.3s were used (45-200bpm).

Missing values simulation

The objective was to delete the same percentage of data from all analysis windows. By doing so, we can directly evaluate the effect of each percentage of missing data on HRV features.

As depicted in figure 3, all windows had the same percentage of missing values but the distribution of deleted RRIs is completely randomised. Details are provided on the analysis window as well as the deletion procedure below.

HRV window :

In order to compute time domain and frequency domain HRV features, the RR timeseries were split into 5min segments, with a 1min sliding window (4min overlap). The choice of a sliding window is to address the discontinuities observed at the edges of each window. It also means a new set of HRV features is available every minute, which is closer to a real-time HRV analysis for continuous health monitoring.

Deletion procedure :

Since the goal is to evaluate the effect of interpolation on HRV features, the same percentage of missing data was removed from each window used to compute HRV features. The steps for the deletion procedure are explained in the pseudo code below.

Algorithm 1 RR deletion procedure

1: Randomly delete P % of the data in the first 5min window

2:
for Each new window i do 3:

Compute N , total number of values to be deleted N = W indowLength×P 100 4:

Determine N overlap number of deleted data in the 4min overlap.

5:

Compute the number of values still to be deleted from the sliding window :

N sliding = N -N overlap 6:
Randomly remove N sliding from the last minute of the window 7: end for In the first iteration, N values are randomly deleted from the first window of the signal. For each window after the first one, the initial step is to compute the total number (N) of data that should be deleted in order to reach the deletion percentage (P). The number of missing values in the 4 min overlap, deleted in the previous iteration, is then computed (N overlap ), and serves to determine the number of data to randomly remove from the last minute of the window N sliding = N -N overlap . At the end of the loop, all the analysis windows had the same percentage of randomly deleted values. The beats were removed away from the window's edges in order to avoid extrapolation problems. Other than this, there was no condition on the number of consecutive beats to be deleted, nor on their positions. The deletion procedure is completely random. It is however obvious that the higher the percentage of deleted data, the larger (and more numerous) the gaps with successive missing beats.

Interpolation methods

The missing RR intervals deleted in the last step were then replaced using four different interpolation methods. Interpolation was also performed independently and iteratively for each HRV window. In each iteration, only the values in the last minute of the window are interpolated since the values in the first 4 minutes (overlap) were filled in the previous iteration. This is in order to simulate a real-time data acquisition and processing. Interpolation methods used in this study are listed below:

• Nearest Neighbour (NN): Zero-order interpolation method that assigns the value of the nearest existing RR interval to the missing beat.

• Linear: First order interpolation method. Derives a straight line connecting the adjacent RR intervals and calculates the missing beats based on the line.

• Shape-preserving piecewise cubic Hermite interpolating polynomial (PCHIP):

A piecewise cubic polynomial determined by the given data and their specified derivatives at the interpolation points [START_REF] Moler | Splines and pchips[END_REF].

P (x k ) = y k , P (x k+1 ) = y k+1 (2) P ′ (x k ) = d k , P ′ (x k+1 ) = d k+1 (3)
The main idea is to determine the slopes d k so that the function values do not overshoot the data values [START_REF] Moler | Splines and pchips[END_REF]. One of the potential ways to determine d k , used in this paper, is briefly explained below.

If δ k and δ k-1 have opposite signs or if either of them is zero, then x k is a discrete local minimum or maximum, so d k is set to be equal to zero. In (figure 4a), the green curved line is the shape-preserving interpolant, formed from two different cubics.

The two cubics interpolate the center value and their derivatives are both zero there [START_REF] Moler | Splines and pchips[END_REF].

On the other hand, if δ k and δ k-1 have the same sign, then d k is a weighted harmonic mean, with weights determined by the lengths of the two intervals around x k .

w 1 + w 2 d k = w 1 δ k-1 + w 2 δ k (4)
where

w 1 = 2h k + h k-1 , w 2 = h k + 2h k-1 .
(h k denotes the length of the k th subinterval:

h k = x k+1 -x k ) and h k-1 the length of the (k -1) th interval.
At the breakpoint, the reciprocal slope of the Hermite interpolant is the weighted average of the reciprocal slopes of the piecewise linear interpolant on either side (figure 4b). The shape-preserving interpolant is formed from the 2 cubics that interpolate the center value and that have slope equal to d k there [START_REF] Moler | Splines and pchips[END_REF].

• Cubic Spline: One popular third degree in-Figure 4: Slopes for PCHIP. δ k and δ k-1 are the two slopes of the piecewise linear interpolant on either side of the breakpoint terpolation method is the cubic spline interpolation, where data points are estimated by fitting a third degree polynomial. A spline is also a piecewise cubic Hermite that is exceptionally smooth, in the sense that the first and second derivatives of consecutive polynomials are equal and thus continuous, ensuring smoothness of the resulting curve. This avoids the problem of the straight polynomial interpolation that tends to induce distortions on the edges of the polynomials [START_REF] Moler | Splines and pchips[END_REF].

The Pchip and the spline methods both perform piecewise cubic Hermite interpolation. They only differ in how the slopes of the interpolant are computed, thus leading to different behaviors when the underlying data has flat areas or undulations.

After the interpolation step, HRV features were estimated on the reconstructed data and compared to the original HRV set from the original signal.

The error was then estimated through the mean absolute percentage error (Mape) in order to identify the best interpolation approach.

HRV analysis

To choose the best imputation approach for HRV signals, the impact of interpolation on multiple HRV features was evaluated through an estimation error to compare features from reconstructed signals to those from the original signal. Features mostly used in literature were selected, they can be separated into two categories, time domain and frequency domain features. (In this study, ectopic beats created by interpolation are not filtered before HRV analysis. SDRR will be referred to as SDNN since the formula is the same.)

SDNN is mostly computed over 24H periods, however, researchers have found significantly shorter periods of analysis to be relevant [START_REF] Baek | Reliability of ultra-short-term analysis as a surrogate of standard 5-min analysis of heart rate variability[END_REF]. In this case 300 seconds (5min) periods were used. Considered as gold standard in quantification of the cardiac risk [START_REF] Cardiology | Heart rate variability: standards of measurement, physiological interpretation and clinical use[END_REF], reflection of both sympathetic nervous system (SNS) and parasympathetic nervous system (PNS) activity can be measured on SDNN which makes it one of the most useful features of HRV analysis.

SDN N = N i=1 (RR i -RR) 2 N -1 (5)
Where :

RR = 1 N N i=1 (RR i ) (6)
RMSSD is the root mean square of successive differences between normal heartbeats. Like SDNN, it takes only normal IBI as an input. This feature reflects more PNS activation than SDNN does.

RM SSD

= N -1 i=1 (RR i -RR i+1 ) 2 N -1 (7)
where N is the number of RR intervals in the signal.

Frequency domain :

Several methods can be used for frequency domain analysis such as Fast Fourier Transform (FFT), auto regressive modeling (AR) or wavelet transform. In this study, we tested both FFT and AR and compared the results to Kubios, the reference software for HRV analysis for validation purposes. Beside the simplicity of implementation, we opted for FFT since the results from Python were the closest to Kubios.

The goal of frequency domain analysis using any of the methods cited above is always to separate HRV signal spectrum into four components which are Ultra Low Frequency (≤ 0.003Hz), Very Low Frequency (0.003 -0.04Hz), Low Frequency (0.04 -0.15Hz), and High Frequency (0.15 -0.4Hz) [START_REF] Cardiology | Heart rate variability: standards of measurement, physiological interpretation and clinical use[END_REF], (respectively ULF, VLF, LF and HF).

Since ULF and VLF generally require long periods of recording not suitable for real-time analysis, they will not be included in this study. Also, their physiological correlates are still unknown which makes them less relevant for e-health applications.

HF and LF on the other hand can be assessed on 1 to 2 min windows respectively [START_REF] Cardiology | Heart rate variability: standards of measurement, physiological interpretation and clinical use[END_REF]. Their ability to reflect the overall cardiac health and the state of the autonomic nervous system (ANS) has been proven by many studies [START_REF] Shaffer | An overview of heart rate variability metrics and norms[END_REF][START_REF] Acharya | Heart rate variability: a review[END_REF], in different contexts including stress [START_REF] Salahuddin | Dependence of heart rate variability on stress factors of stress response inventory[END_REF][START_REF] Benchekroun | Comparison of stress detection through ECG and PPG signals using a random forestbased algorithm[END_REF] and sleep [START_REF] Michail | Eeg and hrv markers of sleepiness and loss of control during car driving[END_REF].

Evaluation metrics

The difference between HRV features from the reconstructed data and those from the original signal was assessed by the Mean Absolute Percentage Error (MAPE) [START_REF] Stankus | A novel adaptive noise elimination algorithm in long rr interval sequences for heart rate variability analysis[END_REF]. Other metrics can also be used but the idea behind choosing the MAPE is to avoid mutual cancellation of the positive and negative errors. Moreover, since each HRV parameter has a wide range [START_REF] Kim | The effect of missing rr-interval data on heart rate variability analysis in the frequency domain[END_REF], normalization by the actual value allows the comparison of differently scaled time-series data.

M ape = 1 n n t=1 A t -E t A t (8)
where : n = number of times the summation iteration happens, which corresponds to the number of HRV windows. A t = Actual value, from the original RR timeseries. E t = Estimated value, from reconstructed signal.

Another interesting parameter to look at is the number of ectopic beats created by the interpolation. As explained before, non physiological beats should be filtered and, eventually, replaced before HRV analysis. The replacement method (ie : interpolation) should not be creating more ectopic beats. We assessed the percentage of abnormal RR intervals (P ectopic ) in the reconstructed signals as follows :

P ectopic = Number of ectopic beats Signal Length (9)

Results and discussion

In this paper, 24 RR time-series of 50min duration were analysed for a total of 1104 HRV windows of 5min duration. To investigate the effect of missing data on HRV features, the same percentage of RRintervals was removed from each window starting from 10% up to 70% of missing values with a 10% step.

The deleted beats were then replaced by four different interpolation methods explained in section 3. An example of data interpolation is shown in Figure 5.

The cubic spline interpolation overshoots the data at some points as can be seen in figure 5. This is due to the requirement for equal second order derivatives at every point. By eliminating this condition, it is possible to prevent, or at least reduce, the overshooting as done by the Pchip method.

Time domain features

According to the results in table 1, SDNN seems to be less sensitive to interpolation. It was the least affected with an estimation error not greater than 5% even with a huge number (70%) of missing data. The same conclusion was found by authors in [START_REF] Salo | Ectopic beats in heart rate variability analysis: effects of editing on time and frequency domain measures[END_REF]. RMSSD on the other hand is much more sensitive to interpolation. The estimation error increases almost linearly with the percentage of missing data in the original signal. As can be seen in table 1, the Nearest Neighbour interpolation yields the least error for both SDNN and RMSSD compared to other interpolation methods.

Mape (%)

Missing Since SDNN is the standard deviation of each RR interval from the mean RR duration, it reflects the LF component in some way whereas the RMSSD correlates with the HF since it uses the difference between successive beats. This may explain why SDNN is much less sensitive to interpolation than RMSSD. In fact, NN interpolation acts as a low-pass filter since it produces flat-like shapes [START_REF] Salo | Ectopic beats in heart rate variability analysis: effects of editing on time and frequency domain measures[END_REF]. In situations where the heart rate is relatively stable and does not vary abruptly, the NN interpolation is most likely to preserve the heart rate variability.

When the percentage of missing data exceeds 50% however, it has been found that the best results for RMSSD estimation are achieved without editing the RR tachograms, i.e without replacing the missing data by any of the interpolation methods used in the study.

[19] also concluded that RMSSD does not require any interpolation to obtain reliable estimations, but they found the threshold to be at 30% instead.

Table 2 summarises RMSSD estimation errors by nearest neighbour approach against no interpolation. Not editing RR time-series yields better RMSSD estimation than editing more than half the data. This however should be verified when the acquisition includes different contexts that may cause the heart rate to vary a lot. The decrease of the Mape when the percentage of missing data increases may be due to the lower number of compared windows. When the missing values are not replaced by any interpolation, remaining RR intervals are just merged. This makes the RR signal much shorter and thus reduces the number of HRV windows.

Frequency domain features are clearly much more sensitive to interpolation as can be seen from table 3. Linear and Pchip interpolation perform almost equally and yield the least estimation error for LF, HF and LF/HF . They are thus considered to be the best interpolation methods for frequency domain features.

Generally speaking, physiological variables such as the Autonomic cardiovascular regulation operates at sufficiently low frequencies [START_REF] Saul | Beat-to-beat variations of heart rate reflect modulation of cardiac autonomic outflow[END_REF] that nothing would be lost using a linear or a Pchip approach. Unless there is a physiological reason to suppose a non-linear trend, linear seems to assume less than the other methods.

Contrary to the time domain analysis, the cubic spline interpolation gives the worst results with an error almost two times greater than all the other interpolation methods for frequency domain features. This can be explained by the fact that cubic splines are prone to severe oscillation and they overshoot at intermediate points. The overshooting introduces many ectopic beats thus increasing the HF components. It has been found in [START_REF] Berntson | Ecg artifacts and heart period variability: don't miss a beat![END_REF] that the presence of only one ectopic beat in a 2 min ECG recording introduces an increase in the HF power of around 10%.

[22] however compared linear, spline backward and forward interpolation. They opted for spline as the best interpolation method. One potential explanation for this is the small percentage of missing data (10%) simulated in the signal and the deletion approach that was not completely random. A maximum of four successive missing values was set. with regards to RMSSD estimation (Table 2). The latter method outperforms the first one when the percentage crosses the 50% threshold. Generally speaking, the Pchip interpolation seems to do well in most cases. It preserves the linear trend of the data while adding very light waves. As explained in [START_REF] Acharya | Heart rate variability: a review[END_REF], the structure generating the RR signal is not only simply linear, but also involves nonlinear contributions. The Pchip interpolation thus seems to better mimic the RR timeseries trend.

Conclusion

In the time domain, nearest neighbour interpolation gives the best results for up to 50% of edited data. Beyond 50%, the best estimation was achieved when the deleted data was not replaced. It seems better not to use any interpolation for RMSSD beyond this threshold. In the frequency domain however, the lowest errors of HRV feature estimation are obtained using linear or Pchip interpolation.

If only one approach had to be chosen for a good overall estimation, the Pchip would be privi-leged because it preserves the linear trend and the slightly non linear contributions in the RR timeseries.

Since HRV features are used for preventive health and users' well-being, it is fundamental to know the effect of missing data on these parameters. The findings of this study, namely the best interpolation methods based on the percentage of missing beats could be used for a data-driven decision-making strategy to decide whether reliable conclusions can be drawn from the signal. This preprocessing step, including filtering and interpolation, is fundamental before any HRV analysis can be performed. It enables continuous passive monitoring of users' cardiovascular activity in a non-obtrusive way despite a relatively poor data quality.

Limits and Perspectives

It is worth bearing in mind that interpolation remains at a mathematical level. Physiological implications and interpretations could further be explored but are outside the scope of this paper. The need and efficacy of interpolation in general should be assessed against the end-goal of HRV analysis. Moreover, in real-life acquisitions, the exact number of missing data in a time gap is unknown.

On the other hand, many additional aspects could be investigated in a future work. Based on the findings described above and table 4, a potential good approach may be using a combination of different interpolation methods chosen based on the HRV feature and the percentage of missing data in each HRV segment. It would be interesting to measure the estimation error of such an approach including different interpolation methods based on the percentage of missing data in each window. The effect of interpolation on other HRV features such as the total spectral power, and Non linear features could also be investigated.

Additionally, it would be very interesting to identify an upper limit for missing beats, in each HRV window, beyond which any interpolation would be pointless. This upper limit would depend once again on the context and the purpose of HRV analysis in the first place. It would help decide whether an HRV segment can be used for a reliable diagnosis or should be discarded.

Figure 2 :

 2 Figure 2: General flowchart of the study including the three major steps of data deletion, interpolation, followed by feature extraction and comparison using the MAPE.

Figure 3 :

 3 Figure3: Example of the deletion procedure on an HRV signal (green signal). We used 5 minutes windows with 4 minutes overlap, which is the recommended and mostly used window size in literature[START_REF] Cardiology | Heart rate variability: standards of measurement, physiological interpretation and clinical use[END_REF]. In the first iteration, a percentage of RRIs is deleted from the first red window, depicted as red and black asterisks. In the next iteration, RRIs already deleted in the 4min overlap (black asterisks) are maintained and only RRIs located in the last minute of the blue window are deleted in order to reach the same percentage of missing data in all windows. Deleted RRI in the second iteration are depicted as blue asterisks. Since both windows have the same deleted data in the 4min overlap segment. The first minute (red arrow) of the i window and the last minute (blue arrow) of the i + 1 window have the same percentage of missing data.

  Time domain :Two of the most known indices were chosen, which are SDNN and RMSSD, for the time domain analysis. SDNN stands for Standard Deviation of Normal to Normal beats. Normal to normal means that ectopic and other abnormal beats have to be removed beforehand. Variations of SDNN such as Standard deviation of RR intervals (SDRR) are sometimes used. The formula is the same, the only difference is that RR time series-for SDRR-include abnormal or false beats.

Figure 5 :

 5 Figure 5: Example of interpolation for 50% missing data. The red arrow indicates the ectopic beats created by the Cubic Spline interpolation

Table 1 :

 1 Mean absolute percentage error of estimated time domain HRV features from 10 to 70% missing data

		% HRV Features	NN	Linear	Pchip	Spline
		RMSSD	3.84 ± 1.15	6.69 ± 0.96	6.56 ±0.99	5.13 ±1.03
	10%	SDNN	0.87 ± 0.31	0.91 ± 0.40	0.86 ± 0.38	1.03 ±0.44
		% ectopic	0	0	0	0
		RMSSD	6.98 ± 2.88	12.76 ± 2.59	12.49 ±2.65	9.51 ±2.38
	20%	SDNN	1.36 ± 0.59	1.43 ± 0.55	1.3 ± 0.48	1.89 ±1.01
		% ectopic	0	0	0	0
		RMSSD	10.17 ± 3.94	19.89 ± 3.23	19.39 ± 3.29	14.84 ± 3.17
	30%	SDNN	1.70 ± 0.52	2.28 ± 1.03	1.96 ± 0.87	2.93 ± 1.42
		% ectopic	0	0	0	0.5
		RMSSD	13.99 ± 4.45	27.63 ± 3.91	26.92 ±4.05	26.11 ± 26.4
	40%	SDNN	2.08 ± 0.52	3.18 ± 1.04	2.56 ± 0.83	7.42 ± 15.25
		% ectopic	0	0	0	0.7
		RMSSD	17.3 ± 6.34	33.66 ± 7.66	32.54 ± 7.48	25.83 ± 6.94
	50%	SDNN	2.63 ± 0.91	4.21 ± 1.57	3.45 ± 1.32	6.60 ± 3.01
		% ectopic	0	0	0	1.3
		RMSSD	20.7 ± 7.87	41.63 ± 8.98	40.87 ± 9.07	31.9 ± 7.37
	60%	SDNN	3.47 ± 1.39	5.31 ± 2.09	4.35 ± 1.83	9.59 ± 6.58
		% ectopic	0	0	0	1.5
		RMSSD	25.57 ± 8.25	50.58 ± 10.6	49.4 ± 10.4	40 ± 8.4
	70%	SDNN	2.63 ± 0.91	4.21 ± 1.57	3.45 ± 1.32	6.60 ± 3.01
		% ectopic	0	0	0	1.5

Table 2 :

 2 Mean absolute percentage error of estimated RMSSD for 50%, 60% and 70% missing data

Table 4

 4

	summarises the best interpolation ap-
	proach for some HRV feature at a specific range
	of missing data. At exactly 50% of missing beats,
	NN and no interpolation approach perform equally

Table 3 :

 3 Mean absolute percentage error of estimated frequency domain HRV features from 10% to

	70% missing data

Table 4 :

 4 Best interpolation approach for HRV features based on the percentage of missing data.

Acknowledgment :

The study was carried out in the frame of a collaboration between two PhD researches. The first one is funded by a grant from the French Ministry of Education The second PhD is funded by the industrial sponsor, Core for Tech and the National Association of Research and Technology in France (ANRT).

Conflict of interest :

None


Aurélie Antunes

Denis Servent

Rémy Béroud

Michel De Waard

Evelyne Benoit

From identification to multi-scale pharmacological evaluation of venom peptides against pain

à la diffusion de documents scientifiques 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.

From identification to multi-scale pharmacological evaluation of venom peptides against pain

Aurélie Antunes 1,2 , Denis Servent 1 , Rémy Béroud 2 , Michel De Waard 3 , Evelyne Benoit 1 Pain care is a global public health priority and current pain relievers, although effective in many cases, are too often associated with side effects. A promising approach to manage pain and avoid central side effects is to target the dorsal root ganglia (DRG) neurons, located in the periphery. These neurons are wellknown to express various proteins involved in pain transmission, including the 1.7 subtype of voltage-gated sodium channels (Na V 1.7), a major and genetically validated target for chronic pain treatment. Thanks to their high affinity and selectivity, venom peptides of animal origin have been identified as an original source of antinociceptive molecules.

The objective of my PhD is to identify and characterize, by multi-scale studies, a large number of novel antinociceptive venom peptides possessing high selectivity and potency for Na V 1.7. For this purpose, an electrophysiological screening on cells overexpressing the human Na V 1.7 allowed to identify novel hits from a unique and large venom collection of Smartox Biotechnology. The different studies, currently undertaken on some selected peptides, aim :

(1)to chemically synthesize these peptides (2) to characterize their selectivity profile and safety issue, notably on cell lines overexpressing Na V subtypes (from 1.1 to 1.8) and on native systems such as the mouse cardiomyocytes and neuromuscular system, (3) to specify their efficiency and mechanism of action on mouse DRG neurons, (4) to evaluate in vivo their analgesic properties in vigilant mice (5) to determine their pharmacophore and optimize their functional properties, and (6) to provide evidence for the clinical potential of the best candidates in correcting de novo identified Na V channelopathies.

Abstract

Based on a collection of numerous and various venoms, this research strategy makes possible to identify and characterize different toxins of therapeutic interest in the context of their development as antinociceptive agents.

The evaluation of functional properties of cyriotoxin-1a reveals that the peptide exhibits marked antinociceptive properties, as evidenced in DRG neurons by its effects on TTX-sensitive Na V channels (Na V 1,7) and on tactile and heat sensitivity of mice in vivo This allows us to consider a clinical development of the studied toxin.

Conclusion In vitro

In vivo

1

  Service d'Ingénierie Moléculaire pour la Santé, ERL CNRS/CEA 9004, CEA de Saclay, Université Paris-Saclay, Gif sur Yvette; 2 Smartox Biotechnology, Saint-Egrève; 3 Institut du Thorax, Inserm UMR 1087 / CNRS UMR 6291, Nantes, France