Valéry Ozenne
email: valery.ozenne@ihu-liryc.fr

Pierre Bour
email: pierre.bour@ihu-liryc.fr

Baudouin Denis De Senneville

Solenn Toupin
email: solenn.toupin@siemens-healthineers.com

Alexis Vaussy
email: alexis.vaussy@siemens-healthineers.com

Matthieu Lepetit-Coiffé
email: matthieu.lepetit-coiffe@siemens-healthineers.com

Pierre Jaïs
email: pierre.jais@ihu-liryc.fr

Hubert Cochet
email: hubert.cochet@chu-bordeaux.fr

Bruno Quesson
email: bruno.quesson@u-bordeaux.fr

Matthieu Lepetit

Assessment of left ventricle magnetic resonance temperature stability in patients in the presence of arrhythmias

Keywords: Magnetic Resonance Thermometry, Proton Resonance Frequency Shift, Temperature Mapping, Arrhythmia, Motion Correction, Patient, Left Ventricle

à 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

Catheter-based radiofrequency ablation (RFA) is a widely accepted method for the clinical treatment of cardiac arrhythmia, with 200 000 procedures performed in Europe in 2014 [START_REF] Raatikainen | Statistics on the use of cardiac electronic devices and electrophysiological procedures in the European Society of Cardiology countries: 2014 report from the European Heart Rhythm Association[END_REF]. However, a significant recurrence rate of arrhythmia is observed: up to 50% of the patients with atrial fibrillation (AF) [START_REF] Ganesan | Long-term outcomes of catheter ablation of atrial fibrillation: a systematic review and meta-analysis[END_REF][START_REF] Heart Rhythm | EHRA/HRS expert consensus statement on cardiac resynchronization therapy in heart failure: implant and followup recommendations and management[END_REF], and 37% with ventricular tachycardia (VT) experience arrhythmia recurrence [START_REF] Tanner | Catheter ablation of recurrent scar-related ventricular tachycardia using electroanatomical mapping and irrigated ablation technology: results of the prospective multicenter Euro-VT-study[END_REF]. For those patients, a redo procedure is required. Arrhythmia recurrence is mainly caused by incomplete ablation/isolation of pathological tissues [START_REF] Berte | VT Recurrence After Ablation: Incomplete Ablation or Disease Progression? A Multicentric European Study[END_REF] or by the recovery/reconnection of the targeted area [START_REF] Darby | Recurrent Atrial Fibrillation After Catheter Ablation: Considerations For Repeat Ablation And Strategies To Optimize Success[END_REF][START_REF] Lin | Electrophysiologic findings and long-term outcomes in patients undergoing third or more catheter ablation procedures for atrial fibrillation[END_REF]. At the present time, RFA are conducted under X-ray guidance, providing very limited inputs on cardiac anatomy and poor information on the created lesion. As of now, electro-anatomic mapping is the key monitoring tool for 3D localization of the catheter and identification of electrical dysfunctions. However, it does not provide anatomical information such as myocardial wall thickness or substrate properties (e.g. presence of scars).

As illustrated in a recent case study by Berte and colleagues [START_REF] Berte | Impact of septal radiofrequency ventricular tachycardia ablation: insights from magnetic resonance imaging[END_REF], transmural ablation in the septum can be difficult to achieve. The ability to create transmural and permanent lesions is, however, of high importance to ensure efficient treatment of the arrhythmia [START_REF] Wittkampf | RF catheter ablation: Lessons on lesions[END_REF]. In such a situation, MRI-guided electrophysiology (EP) procedures embedding real-time assessment of the lesion formation might enable the physician to precisely target the remaining viable tissues to be ablated. Real-time thermometry and thermal dose computation could therefore provide a therapy end-point, ensuring transmurality and contiguity of the lesions. Once validated, cardiovascular magnetic resonance (CMR) thermometry could benefit from RF-power retro control loop, in order to maximize the heat transfer through the myocardial wall and to stop RF delivery once transmurality is achieved.

Recently, MRI-guided EP has been reported as the only alternative to make catheter procedures safer and more efficient [START_REF] Koopmann | Why hesitate introducing real-time magnetic resonance imaging into the electrophysiological labs?[END_REF]. Paetsch et al. [START_REF] Paetsch | Cardiovascular Magnetic Resonance-Guided Electrophysiological Interventions: Radiofrequency Ablation of Typical Atrial Flutter[END_REF] and Chubb et al. [START_REF] Chubb | Development, Preclinical Validation, and Clinical Translation of a Cardiac Magnetic Resonance -Electrophysiology System With Active Catheter Tracking for Ablation of Cardiac Arrhythmia[END_REF] have performed successful ablations in respectively five and nine patients using a fully MR-compatible mapping and ablation catheter at 1.5 T. However, no real-time assessment of the lesion formation was performed. Lesion formation monitoring can be achieved in real-time using MR-thermometry which may improve both safety and efficacy of catheter ablation procedures [START_REF] Kolandaivelu | Noninvasive assessment of tissue heating during cardiac radiofrequency ablation using MRI thermography[END_REF]. The ideal RF lesion has to be transmural while preserving the surrounding healthy tissue. The mean and maximum left ventricular myocardial wall thickness is about 6 mm and 8 mm respectively [START_REF] Kawel | Normal left ventricular myocardial thickness for middle-aged and older subjects with steady-state free precession cardiac magnetic resonance: the multi-ethnic study of atherosclerosis[END_REF]. Therefore, to obtain an accurate temperature estimate, CMR thermometry must provide a spatial resolution in the range of the millimeter, while ensuring fast acquisition and compensate for cardiac motions, respiratory motions, and associated magnetic susceptibility artifacts.

Recent CMR thermometry studies [START_REF] Ozenne | Improved cardiac magnetic resonance thermometry and dosimetry for monitoring lesion formation during catheter ablation[END_REF][START_REF] Toupin | Feasibility of real-time MR thermal dose mapping for predicting radiofrequency ablation outcome in the myocardium in vivo[END_REF][START_REF] De Senneville | Feasibility of fast MR-thermometry during cardiac radiofrequency ablation[END_REF] have used an echo planar imaging (EPI) acquisition sequence with ECG triggering and respiratory compensation, using dedicated MR-compatible instrumentation. In Toupin et al. [START_REF] Toupin | Feasibility of real-time MR thermal dose mapping for predicting radiofrequency ablation outcome in the myocardium in vivo[END_REF], a correlation coefficient of 0.87 indicating good agreement between cumulative thermal dose images computed from temperature maps and actual lesion sizes was reported. In this study, imaging was performed with an interpolated in-plane spatial resolution of 0.8 mm and a 3 mm slice thickness. CMR thermometry was also evaluated on healthy volunteers and showed temporal standard deviation of temperature around 2°C in the left ventricle. This variability (also referred as uncertainty) was considered sufficient in perspective of temperature monitoring during cardiac arrhythmia treatment with RFA. The current thermometry pipeline uses data acquired during a learning phase to correct motion-related artifacts on temperature images during the ablation. Thus, unobserved situations (like arrhythmia and irregular breathing events) during the learning step may result in a decreased temperature precision. In the literature, CMR thermometry in humans has only been validated on volunteers without heating either at 3T [START_REF] Grissom | Hybrid referenceless and multibaseline subtraction MR thermometry for monitoring thermal therapies in moving organs[END_REF][START_REF] Hey | Towards optimized MR thermometry of the human heart at 3T[END_REF] or at 1.5T [START_REF] Ozenne | Improved cardiac magnetic resonance thermometry and dosimetry for monitoring lesion formation during catheter ablation[END_REF][START_REF] Toupin | Feasibility of real-time MR thermal dose mapping for predicting radiofrequency ablation outcome in the myocardium in vivo[END_REF][START_REF] De Senneville | Feasibility of fast MR-thermometry during cardiac radiofrequency ablation[END_REF][START_REF] Toupin | Combination of principal component analysis and optical-flow motion compensation for improved cardiac MR thermometry[END_REF]. In Grissom et al. [START_REF] Grissom | Hybrid referenceless and multibaseline subtraction MR thermometry for monitoring thermal therapies in moving organs[END_REF], a hybrid image model for PRF-shift thermometry has been validated in healthy volunteers in the left ventricle using a spiral acquisition. In de Senneville et al. [START_REF] De Senneville | Feasibility of fast MR-thermometry during cardiac radiofrequency ablation[END_REF] and in Hey et al. [START_REF] Hey | Towards optimized MR thermometry of the human heart at 3T[END_REF], temperature mapping was performed in eight and ten healthy volunteers respectively using an EPI acquisition. Susceptibility effects were corrected with a multi baseline model approach. In all mentioned studies, neither the presence of arrhythmia nor the presence of irregular breathing have been observed and measured.

In the present study, we assess the reliability of the CMR thermometry by measuring temperature changes in patients during free-breathing acquisition on a clinical scanner. The temperature precision was assessed by measuring on the left ventricle the temporal standard deviation and temporal mean of consecutive temperature map during 3 min. In presence of arrhythmia, since RR intervals may vary over the duration of the scan, we exploit a new data processing method to characterize the stability of the cardiac cycles from ECG recordings and sort RR intervals into distinct arrhythmic events. Finally, to account for irregular respiratory and cardiac motion that may impact the precision of the temperature estimate, we propose to improve the precision of the CMR thermometry using two new algorithms that automatically detect and reject corrupted dynamic acquisitions. These strategies were tested and combined in one patient.

Methods

Patient population

The study group consisted of fourteen patients referred for CMR assessment of known or suspected heart disease. Patients were aged 59 ± 17 years (mean ± standard deviation) ranging from 15 to 76 years (ten men). Each patient was informed about the protocol approved by ethic committee and consented to participate to the study. The inclusion criteria were patients requiring a clinical MRI examination, affiliated to healthcare insurance, male and female, adults of 18 years old or higher, patients between 15 years old and 18 year old after parent authorization, given consent after reading information notice of the protocol.

The exclusion criteria were as follows: patient under 15 years old, or unable to give consent, or patient already enrolled in another clinical study, or patient deprived of liberty by judicial or administrative decision, or adult protected by law, or patient not eligible to MRI examination.

MR-thermometry acquisition

All imaging was performed on a 1.5 T MRI scanner (Avanto, Siemens Healthineers, Erlangen, Germany) using two 16-channel cardiac coils surrounding the thorax. A fat-saturated, single-shot EPI sequence was combined with GRAPPA parallel acquisition technique (acceleration factor of 2) to achieve a 1.6x1.6x3 mm 3 spatial resolution (except for the patients #1 and #2 in which the spatial resolution was 2x2x3 mm 3 ), TE/TR/FA = 18-20 ms/110 ms/60°, bandwidth = 1568 Hz/px, ETL = 32.5 ms, echo spacing =0.77 ms, FOV=180x180 mm 2 , and 6/8 partial Fourier. The acquisition was synchronized on the ECG and the acquisition was repeated over 250 consecutive heart beats during free breathing.

Four to five slices (spacing between slices = 0.3 mm) were acquired sequentially at each heart beat in coronal orientation with phase-encoding set to head feet direction. The number of slices was adjusted to ensure that all slices could be fit in a single heartbeat. Thus, the total acquisition time of one multi-slice acquisition range between 440 ms to 550 ms. Saturation slabs surrounded the FOV in the phase encoding direction to avoid aliasing and two additional saturation bands were set parallel to the imaging slices to reduce the signal of blood inside the cavities, as suggested in [START_REF] Ozenne | Improved cardiac magnetic resonance thermometry and dosimetry for monitoring lesion formation during catheter ablation[END_REF]. A standard shimming box centered on the heart was used.

Image reconstruction and temperature mapping

Image reconstruction and additional processing were implemented into the Gadgetron framework [START_REF] Hansen | Gadgetron: an open source framework for medical image reconstruction[END_REF], as previously described in [START_REF] Ozenne | Improved cardiac magnetic resonance thermometry and dosimetry for monitoring lesion formation during catheter ablation[END_REF]. Online data streaming, reconstruction and processing were seamless from the user perspective. Images were sent back on-the fly to the scanner immediately after reconstruction without any user interaction and simultaneously to a remote host computer for display on a dedicated thermometry software (Thermoguide™, IGT, Pessac, France).

A brief reminder of the thermometry image processing method is provided. All the following steps below were implemented into the Gadgetron framework and processed in real-time:

1) In-plane zero-filling at 0.8x0.8x3 mm 3 voxel size was added to the pipeline to reduce partial-volume artifacts [START_REF] Du | Reduction of partial-volume artifacts with zerofilled interpolation in three-dimensional MR angiography[END_REF][START_REF] Todd | The effects of spatial sampling choices on MR temperature measurements[END_REF], allowing better visualization of the myocardium.

2) The first ten acquisitions (from repetition numbers 1 to 10) were used to determine the reference image located at the middle course of the respiratory cycle. Non-rigid registration [START_REF] Black | A framework for the robust estimation of optical flow[END_REF] was performed on magnitude images to correct in-plane displacement resulting from respiratory motions

3) The following thirty dynamic acquisitions (learning steps, from repetition numbers 11 to 40) were registered to the reference image and motion fields were stored in a library. A parameterized description of the phase variations was performed using a principal component analysis of the motion fields [START_REF] De Senneville | MR-guided thermotherapy of abdominal organs using a robust PCA-based motion descriptor[END_REF]. MR-thermometry was then computed by subtracting the current phase image from a synthetic reference phase image, built by the model. Temperature data were available from repetition numbers 41 to 250. 4) Compensation of spatial-temporal drift of the phase and low pass temporal filtering (first order low pass Butterworth filter and cutoff 0.14 Hz frequency) were implemented [START_REF] Ozenne | Improved cardiac magnetic resonance thermometry and dosimetry for monitoring lesion formation during catheter ablation[END_REF].

ECG Recording and Synchronization

ECG was recorded using standard 3-lead ECG acquisition provided by the MR manufacturer. A homemade algorithm [START_REF] Pan | A real-time QRS detection algorithm[END_REF] was designed to identify the QRS peaks with an adaptive threshold and to synchronize the ECG signal with the temperature map acquisition. Processing of the ECG and synchronization were done retrospectively, after the acquisition. ECG was available for all patients but could not be recorded and post-processed for the patients #1, #2 and #3 due to technical issues. Since the acquisition time stamp of each image is included in the image header, the RR duration intervals were estimated in these cases from the time delays between successive acquisitions of the same slice.

Categorization of beats

As proposed by Contijoch et al. [START_REF] Contijoch | Quantification of Left Ventricular Function With Premature Ventricular Complexes Reveals Variable Hemodynamics[END_REF], categorization of beats based on the duration of two consecutive RR intervals were used to quantify the prevalence of arrhythmia for each patient (see Figure 1A). The RR duration of the loading beat RRi-1 was plotted versus the RR duration of the current beat RRi allowing the identification of the presence of distinct beat types in clusters. Clusters were defined using a kmeans algorithm with 1 to 5 groups function implemented in MATLAB (The MathWorks, Inc., Natick, Massachusetts, United States). In sinus rhythm, a single cluster was observed whereas multiple clusters indicate the presence of distinct beat types. To quantify the prevalence of arrhythmia during the acquisition, mean and standard deviation of the RR duration was calculated. Then, the prevalence of irregular beat type was calculated as the percentage of beats identified in that beat type relative to the total number of observed beats along the acquisition sequence. In order to investigate the impact of arrhythmia on the temperature measurements, data were categorized in three groups of patients: patients in sinus rhythm, patients with fewer than 15% prevalence and patients with more than 15% prevalence.

The prevalence of arrhythmia is provided for each patient on the last column of Table 1.

Rejection of corrupted images due to irregular respiratory motions

Irregular respiratory motions may lead to out of plane motion and result in erroneous images in the time series. The proposed thermometry method assumes that slices are positioned along the respiratory motion. Violation of this assumption leads to erroneous temperature information and corresponding images need to be removed [START_REF] Hey | Towards optimized MR thermometry of the human heart at 3T[END_REF]. The following automatic detection algorithm is proposed (see Figure 1B), based on the analysis of phase images:

-A Sobel operator (3x3 convolution kernel) was applied to each registered phase image to compute its gradient image in both directions (Gx, Gy). The phase image after applying the Sobel operator (Psobel) was then calculated (Figure1B). Using this method, phase images showing a large number of phase wraps were weighted with a higher value.

-A binary mask computed on registered magnitude image (with a threshold value equal to 10% of the maximal value) was applied on the Psobel image to exclude voxels showing a low SNR (imprecise phase definition). The resulting image is named the "gradient change of the phase" (Gphase) image. A typical mask is shown in Figure1B.

-The percentage of the remaining voxels in the Gphase image was used as a surrogate for phase image quality. The Gphase metric is defined as the sum of all voxels in the Gphase image normalized by the total number of voxels (excluding the "zero" voxels outside the binary mask).

A rejection strategy was finally implemented based on the temporal variation of the Gphase metric. The mean and standard deviation of the Gphase metric was calculated at the end of the learning step. Then the Gphase metric was subtracted by the mean value over the learning step to obtain the baseline value. The threshold of the rejection criteria was empirically chosen as six times the standard deviation of the Gphase metric during the learning step. During the following dynamic acquisitions, all phase acquisitions with a metric superior to the threshold of rejection criteria were removed from the temperature calculation.

Rejection of corrupted images during arrhythmia-related beats

Using offline processing in MATLAB, we investigated in one patient (#14) the possibility to reject dynamic acquisitions based on the timing of acquisition window relative to the beat duration. If a premature contraction occurs, the RR interval is shortened, resulting in incorrect acquisition timing.

First, the clustering of the ECG has been used to identify images acquired during an ectopic beat.

Categorization of beats were computed to identify the presence of clusters. Then, as an example of strategy, the phase images corresponding to the arrhythmia-related beats have been removed from the temperature pipeline.

Data analysis

In absence of heating, the heart temperature is constant, therefore precision of the method can be assessed by measuring on the left ventricle the temporal standard deviation (also referred as temperature uncertainty) and the temporal mean of consecutive temperature maps [START_REF] Grissom | Hybrid referenceless and multibaseline subtraction MR thermometry for monitoring thermal therapies in moving organs[END_REF][START_REF] Weidensteiner | Stability of real-time MR temperature mapping in healthy and diseased human liver[END_REF][START_REF] Weidensteiner | Real-time MR temperature mapping of rabbit liver in vivo during thermal ablation[END_REF][START_REF] Rieke | Referenceless PRF shift thermometry[END_REF][START_REF] Kickhefel | Accuracy of real-time MR temperature mapping in the brain: a comparison of fast sequences[END_REF]. Temperature changes can remain within an interval of one or two degrees as temperature increase required for catheter ablation is generally higher than 20°C [START_REF] Haines | The biophysics of radiofrequency catheter ablation in the heart: the importance of temperature monitoring[END_REF][START_REF] Haines | Tissue heating during radiofrequency catheter ablation: a thermodynamic model and observations in isolated perfused and superfused canine right ventricular free wall[END_REF]. To evaluate temperature precision without heating, temporal mean (µT) and standard deviation (σT) of temperature in each voxel were computed over the 210 dynamic acquisitions following the learning step. Theoretically the temporal mean (µT) should be null and its value reflects a potential offset due to phase drift for instance. The temporal standard deviation (σT) reflects the stability of the method. Region of interests (ROIs) were manually drawn on the left ventricle in each slice. After image reconstruction and processing, residual voxels displaying local phase wraps may remain and lead to high apparent temperature uncertainty on temperature images. To account for that, voxels with σT > 10 °C were removed from the data analysis. The number and percentage of these voxels per slice were reported. Table 2 (see Supporting Information) summarizes information about temporal standard deviation (σT) and temporal mean (µT) for each patient and for each slice. The number of voxels in the ROI on the LV is also indicated as well as the percentage of myocardial voxels displaying σT below 2°C and higher than 10°C (the latter being outliers excluded from the data analysis).

Unless stated otherwise, the results from the data analysis have been computed with all dynamic acquisitions (without the proposed rejection methods). Finally, the σT and µT distribution were computed using all slices for each patient, the results were plotted in Figure 2A and 2B with box and whisker graphical representation using five quantiles corresponding to: 10%(T10), 25% (T25), 50% (T50), 75% (T75) and 90% (T90) of the distribution of the σT and µT values, respectively. Finally, a regional analysis was performed as follow: for each slice on each patient, we segmented the left myocardium into twelve angular sectors of equal aperture, centered on the barycenter of the ROI (see Figure S1A). In each sector, the mean of the temporal standard deviation of the temperature was computed and averaged over slices and then over patients. The percentage of excluded voxels was also computed using the same methodology.

Results

Table 1 summarizes information about each patient, including actual scan duration, the RR duration intervals and the number of arrhythmia-related beats. The mean ± standard deviation (SD) heart rate across all patients was 68 ± 8 bpm, ranging from 54 to 90 bpm. The mean ± SD of the total scan duration was 3 min 50 s ± 30 s. The corresponding mean ± SD of the total time where temperature mapping was available was 3 min 13 s ± 25 s. Five patients were in stable sinus rhythm without important variations of the RR duration during MRI scanning. Nine patients showed ectopic beats on a sinus rhythm background: six had sporadic premature ventricular or atrial contractions while three had high frequency of premature ventricular or atrial contractions. For each group, a representative case study is detailed followed by the data evaluation on all patients.

Thermometry in patient #5 under sinus rhythm

Patient #5 was in sinus rhythm with RR duration intervals of 1031 ± 37 ms. ECG recording (Figure 3A) and categorization of beats (Figure 3B) resulted in a single beat cluster. Averaged registered magnitude images (over 240 dynamic acquisitions) and phase images measured in coronal view are shown in Figures 3C and3D. Registered magnitude images did not show image inconsistency and had welldelimited edges between the heart and the adjacent organs, suggesting that the motion correction was sufficiently accurate to compensate the respiratory motion. On each slice, a phase wrap was visible close to the heart-liver-lung interface (indicated by an orange arrow on Figure 3D). Figure 3E displays the temperature uncertainties (σT) overlaid on an averaged registered magnitude image (Figure 3C) with a temperature uncertainty cut off of 10 °C. Temperature precision appears spatially consistent with most of the voxels displaying a σT of approximately 2°C over the whole myocardium. Temperature uncertainty was higher than 10 °C in areas with remaining phase wraps, indicated by orange arrows in slices 4 and 5. The fraction of myocardial voxels displaying σT higher than 10°C is less than 2% while the fraction of myocardial voxels displaying σT lower than 2°C reaches 89% (see Table 2). The σT mean ± SD was 1.4 ± 0.8 °C, 1.8 ± 0.8°C, 1.4 ± 0.5°C, 1.4 ± 0.5 °C and 1.3 ± 0.5 °C, for each of the five slices, respectively (and 1.4 ± 0.7 °C over all slices). The average number of voxels in the ROI over all slices included 3654 ± 554 voxels (mean ± SD) and the mean temperature over time (µT) were 0.1 ± 1.4 °C (mean ± SD). Typical examples of temperature data over time are displayed in Figure S3: 150 voxels located in the ROI surrounding the LV were randomly selected and displayed using a color code corresponding to the temperature uncertainty.

Thermometry in patient #7 with less than 15% of arrhythmia events

Patient #7 had multiple events of arrhythmia during the acquisition sequence (Figure 4A) which resulted in 4 beat clusters (Figure 4B). The prevalence of arrhythmia was low with 20 beats over 255 beats for a total acquisition time of 3 min 48 s. For this patient, the acquisition window was adjusted between midsystole to mid-diastole (trigger delay set to 260 ms) to reduce the impact of arrhythmia on thermometry and four slices were acquired at each heartbeat (Figure 4C). In each slice, the average registered image showed well-delimited edges between the myocardium and adjacent organs which are indicative of the registration quality. Phase images show good quality during the cardiac contraction (Figure 4D) resulting in stable temperature all over the myocardium (Figure 4E). We note that blood signal was not attenuated in this configuration except for the last slice (Figure 4C). The σT mean ± SD in each slice was 1.4 ± 0.8 °C, 1.4 ± 0.8°C, 1.9 ± 0.8°C and 1.9 ± 0.9 °C (1.6 ± 0.8 °C over all slices), respectively. ROIs over all slices included 6398 ± 511 voxels with a mean temperature (µT) of 0.3 ± 1.8 °C over time. The number of voxels was increased due to the largest extent of the myocardium during systole.

Thermometry in patients #12 and #6 with more than 15% of arrhythmia events

Patient #12 had multiple events of arrhythmia (Figure 5A) which resulted in 4 beat clusters (Figure 5B).

In total, 89 beats over 257 were associated with arrhythmia over a total acquisition time of 3 min 34 s.

Blurring on the last two magnitude images (slices 4 and 5) might be due to cardiac cycle shortening and image registration issues. Unstable temperature estimation near the apex was observed and may be related to a thinner myocardium wall and the presence of phase wraps indicated by a red and orange arrows, respectively. However, phase image quality was found acceptable (Figure 5D) but not sufficient to compute accurate temperature maps over the whole myocardium (Figure 5E). Most of the voxels covering the apex showed a temperature uncertainty higher than 5 °C. The σT mean ± SD in each slice was 3.3 ± 1.7 °C, 2.7 ± 1.4 °C, 3.0 ± 1.5 °C, 3.3 ± 1.6 °C and 3.6 ± 1.5 °C (3.2 ± 1.6 °C over all slices). The mean ± SD size of the ROIs all slices was 2814 ± 165 voxels, the temperature mean (µT) was 0.2 ± 2.8 °C. The number of voxels was lower in this patient due to a smaller myocardium. Typical examples of temperature data over time are displayed in Figure S4 for patient #12.

Patient #6 had permanent ongoing arrhythmia (Figure 6A) which resulted in 5 beat clusters (FIG 6B).

Four slices were acquired at each heart beat due to shorter RR duration. Total acquisition time was 3 min 35 s. Blurring on magnitude image was observed with not well-defined edges which denotes insufficient or inaccurate registration (Figure 6C) although phase images show acceptable quality. In this regard, as arrhythmia influence the cardiac contraction, the images processed by the registration algorithm acquired during ectopic beat can differ from the images acquired without arrhythmia. We note the presence of phase wraps indicated by orange arrows affecting the temperature mapping. Most of the voxels remained below 2°C uncertainty but blurring on magnitude images foresees potential voxels misregistration on the myocardium, particularly in slices 1 and 2. The σT mean ± SD in each slice was 1.6 ± 0.8 °C, 1.7 ± 0.8°C, 2.5 ± 1.3 °C and 2.0 ± 0.6 °C (1.9 ± 1.0 °C over all slices). The mean ± SD size of the ROIs over all slices was 2981 ± 107 voxels, the temperature mean over time (µT) was 0.2 ± 1.7 °C.

Data analysis in each group of patients

The histogram distribution of temporal standard deviation (σT) is displayed for each acquisition and for each slice in Figure S2 using a log scale along the x-axis for each patient and each slice. Voxels excluded from the data analysis (σT >10°C) are located on the right side of the red line. Figure 2 displays the data analysis of σT and µT over all patients. Over all patients under sinus rhythm, the mean standard deviation of temperature σT was 1.6 ± 0.3 °C. 63 ± 94 voxels per patients for each slice were excluded due to residual phase unwraps, over a total of 4284 ± 1752 voxels (excluded voxels represent 1% and are indicated in dark red). Over all patients with fewer than 15% arrhythmic events, the temperature uncertainty was found similar to values reported in patients under sinus rhythm. The mean standard deviation of temperature σT over these patients was 1.7 ± 0.4 °C. 33 ± 20 voxels for each slice were excluded over a total of 4512 ± 1340 voxels (excluded voxels represent less than 1%). Not surprisingly, temperature standard deviation over the three patients showing more than 15% of arrhythmic events was found higher than for patients in sinus rhythm or patients with fewer than 15% of arrhythmic events.

The mean standard deviation of temperature σT over these three patients was 2.5 ± 0.6 °C. Moreover, 144 ± 113 voxels per patients for each slice were excluded over a total of 3286 ± 678 voxels (excluded voxels represent less than 4 %). In all patients, temperature uncertainty remained below 2°C in more than 50% of the voxels of the left ventricle in twelve over fourteen patients. In six over fourteen patients (Figure 2C) it represent 75% of the voxels. To perform a regional analysis of the temperature uncertainty, Figure S1B displays a color-coded representation of the mean value of the temporal standard deviation of the temperature in each of the twelve segments of the left myocardium. In patients under sinus rhythm and with fewer than 15% of arrhythmic events, values remained around 2°C in each segment, with a maximal value of 2.3 °C. In patients with high prevalence of arrhythmia, higher values are reported, with the segment #5 displaying temperature uncertainties of 4.2°C, corresponding of the liver-lung-heart interface. The second row in Figure S1C displays the percentage of excluded voxels.

For most of the segments, less than 5% of the voxels were removed. In the segment #5, the percentage of removed voxels was higher than 25%, with a maximal value of 60% for patient #12 due to the presence of a phase wrap and a low SNR in this area (see the orange arrow in Figure 5C).

Automatic rejection of corrupted phase images due to irregular respiratory motions

To illustrate the impact of large magnetic field variations on temperature mapping, the magnitude and phase images with the corresponding ECG, and Gphase metric of patient #14 were plotted in Figure 7.

The magnitude images in the dynamic acquisitions #221 and #222 showed larger respiratory amplitudes and lower signal to noise ratios. The corresponding phase images were characterized by the presence of phase wraps. A small variation on the ECG-V2 caused by the sudden movement can also be observed.

The calculation of the Gphase metric allows the automatic identification and exclusion of this data in the thermometry pipeline. It is important to underline that the presence of ectopic beats did not modify substantially the phase images and so the Gphase metric. Specificity of the method was necessary in order to keep the maximum number of dynamic acquisitions to follow temperature elevation.

This method has been applied on all datasets. In most of the datasets (11/14), we denoted flat or regular Gphase metrics as illustrated in Figure S5. Regular variation due to breathing is well visible in patient #2, #3, #8 as well as sudden variation in patient #10 and # 14. In three datasets, peak strikes corresponding to the presence of phase wraps were observed in the Gphase metric, resulting in removal of four images in the time series at maximum. Figure 8 summarizes the data improvement for the three patients in which dynamic acquisitions have been rejected. In all three patients, the temporal standard deviation of the temperature was slightly decreased (Figure 8A) without modification of the mean temporal value of temperature (Figure 8B). However, the percentage of voxels with a temporal standard deviation higher than 10°C was reduced (Figure 8C) for each patient. The temporal standard deviation of the temperature before and after data removal using the Gphase rejection method were found statistically different (Student's t-test, P<0.001).

Potential of combining rejection methods of arrhythmia and irregular respiratory events

Using offline processing in MATLAB, we investigated in patient #14 the possibility to reject dynamic acquisitions based both on the timing of acquisition window relative to the beat duration and the presence of irregular breathing. In Figure 8D, timing of the data collection is plotted in red dashed box, showing that the last slices were acquired during systole instead of diastole for each point located left of the vertical dashed line. The clustering of the ECG of this patient resulted in five clusters composed of a main cluster and four minor clusters. As an example of strategy, the phase images corresponding to the four minor clusters have been removed from the temperature pipeline (41 acquisitions in total). The resulting standard deviation of temperature is improved (Figure 8A) and the mean value (Figure 8B) is nearly identical. The percentage of voxels displaying a temperature standard deviation below 2°C was increased by approximately 18 % and reached 74 % (Figure 8C, last vertical bar). Figure 8E left displays the temperature standard deviation map of the last slice acquired on patient #14 without rejection of any image, whereas the central image displays the resulting map after rejection of images using Gphase metric (irregular breathing) and image on the right is the result after both Gphase metric and keeping only images of the main cluster (displayed in Figure 8D). The temporal standard deviation of the temperature is clearly improved on left ventricle and left atrium. The temporal standard deviation of the temperature distributions before and after data removal using the clustering rejection method were found statistically different (Student's t-test, P<0.001). Although the temporal resolution of the thermometry was impacted due to the rejected data, the mean and maximal delay between successive acquisitions were increased from 1.1 s and 2.1 s to 1.3 s and 8.3 s respectively for the case study #14.

Discussion

Quantitative monitoring of temperature evolution, by providing valuable information of the lesion formation, is a key element for the development of future MR-guided RF ablation procedures. In this paper, we present the first human non-heating baseline thermometry results using conventional MR equipment. The objective is to characterize fast changes in temperature measurements caused by noise and physiological motion. We previously developed the processing pipeline and validated on a large animal model [START_REF] Ozenne | Improved cardiac magnetic resonance thermometry and dosimetry for monitoring lesion formation during catheter ablation[END_REF][START_REF] Toupin | Feasibility of real-time MR thermal dose mapping for predicting radiofrequency ablation outcome in the myocardium in vivo[END_REF]. Here, we transferred it to a clinical environment and applied to patients in complement addition to a routine clinical MR scan. The raw data were streamed in real-time to a separate workstation and the images were reconstructed and processed online using the Gadgetron framework.

The resulting images were streamed online back to the MR database and to a separate computer for display. Images were collected in patients showing some degree of arrhythmia and under free breathing conditions.

The results show that cardiac MR-thermometry is feasible in real-time with acceptable uncertainty.

Image quality, and temperature uncertainty were precisely described (Figure 3, 4, 5, 6) in four representative patients (with stable sinus rhythm, with sporadic and high frequency of premature ventricular, with an infarct). For each patient, one lead ECG was plotted during 80 seconds to depict the variety of cardiac rhythm. Additionally, an evaluation of the temperature measurement precision (Figure 2, Figure S1, Table 2) was performed for each slice and patient. No significant difference was found when comparing the distribution of σT over the slices for each patient. In all patients, the histogram distributions (Figure S2) show a unimodal narrow distribution with the exception of patient #12.

Although the latter presents events of arrhythmia, additional phenomenon could explain the larger distribution observed and related decrease in temperature precision. Multiple factors affected the phase signal quality: the myocardial wall was thinner in this patient (Figure 5D), the presence of persistent wrapped phase as well as a low SNR were noticed in the lower right part of the myocardium (Figure 5C). Although the cohort of patients was limited to fourteen, we were able to map temperature in the LV myocardium in each patient with acceptable uncertainty (smaller than 2°C in the majority of LV voxels). These patients had different heart rate (from 54 to 90 bpm) and some of them had ectopic beats together with irregular breathing, both conditions being expected to degrade thermometry. The previous studies of CMR thermometry in humans have been validated in healthy volunteers in free breathing without irregular cardiac contraction or sudden respiratory motion and on large animals under assisted ventilation and full anesthesia. In this study, we showed that thermometry is feasible and accurate in free breathing awake patients. However, the resulting standard deviation of temperature was found higher in patients than in the preclinical studies where physiological conditions were better controlled.

In this study we have investigated this MR-thermometry sequence where the acquisition is triggered to the diastole. During this period of time, the heart is in resting state and shows limited movements. An alternative strategy would be to trigger the acquisition in systole, thereby reducing the impact of arrhythmia on temperature uncertainty (as shown in patient #7). However, this could expose the images to motion artifacts during heart contraction.

To account for irregular respiratory and cardiac motion, we propose new algorithms that are compatible with online data processing. Removal of corrupted dynamic acquisition is a necessary process to ensure reliable temperature measurement and thermal dose computation [START_REF] Hey | Towards optimized MR thermometry of the human heart at 3T[END_REF]. In Hey et al. [START_REF] Hey | Towards optimized MR thermometry of the human heart at 3T[END_REF], the Pearson intercorrelation coefficient computed on module images was used to remove images containing artefacts from erroneous slice tracking. In the present study, phase wraps in presence of large amplitude motion were observed. Since temperature can be strongly affected by phase wraps that are difficult to compensate in mobile organ (like in the heart), we proposed here a new metric (the Gphase metric) that identifies irregular respiratory events using the phase images of the thermometry sequence. We have shown that discarding a limited number of images in the time series (4/250 at maximum) improves the precision of temperature estimate (Figure 7). The sensitivity in the Gphase metric could be optimized by changing the rejection threshold value. In this first implementation, the rejection threshold value was set to remove unambiguous erroneous images as shown in Figure 7 and S4. Thus, the number of rejected images was very low and the temporal resolution of the thermometry was preserved. To identify irregular cardiac events, RR intervals of ECG tracings can be divided into clusters using an already published method in a different context (binning of data to create cine images for each cluster). We showed from four different examples that ECG of patients can be categorized from stable rhythm (sinus) to different degrees of arrhythmia through a simple 2D scatter plot that can be populated dynamically.

In our current implementation, such a plot was computed offline. However, the analysis of ECG may be performed prior to MR-thermometry in order to optimize the acquisition parameters (e.g. delay between QRS detection and data collection, number of slices to acquire for each heartbeat, etc.). Thus, the influence of arrhythmia on the uncertainty of the temperature estimate may also be reduced using an optimized acquisition protocol for each patient. In our current implementation, such a pre-analysis could not be performed since raw ECG were not yet available in the ISMRMRD format and thus could not be exploited online. Nevertheless, this method is conceptually compatible with real-time analysis and allows to decide whether data can be kept or rejected for temperature calculation. Therefore, we propose a basic method to reject data that do not belong to the main cluster of RR intervals to improve thermometry precision, at the cost of a reduced number of available images in the time series. The resulting mean frame rate was 1.3 s, which appears acceptable in the context of monitoring RF ablations.

However, the maximal delay between consecutive acquisitions was 8.3 s. It could be reduced by rejecting only the cluster of beats where the RR interval is shorter than the acquisition time. Although removal of data reduces the monitoring performances, no temperature monitoring is available during RFA in current clinical practice. Therefore, even in presence of highly arrhythmic situation, online display of remaining temperature images would be valuable. Pacing of the heart during RFA may also help in improving thermometry performance, as already proposed on preclinical study [START_REF] Toupin | Feasibility of real-time MR thermal dose mapping for predicting radiofrequency ablation outcome in the myocardium in vivo[END_REF]. Both retrospective strategies operate on different data stream and information. The Gphase metric uses the phase image quality while the arrhythmia rejection uses the RR interval duration. For instance, in presence of arrhythmia as shown in patient #14 in Figure 7, the Gphase metric kept all images with good phase quality even those that correspond to ectopic beats. Thus, this implementation offers flexibility and independence to find an appropriate tuning for each rejection strategy.

Limitations

The central objective of this paper is to assess the stability of phase (and therefore temperature estimate) in patients in absence of heat-based therapy, but potentially exhibiting arrhythmia. We show that temperature stability under these conditions appears acceptable, in view of the achieved local temperature during clinical treatment (> 65°C) using catheter RFA. However, other sources of artifacts may corrupt image quality such as susceptibility artifacts of the catheter itself. It was shown in a previous study using similar thermometry pipeline that reliable temperature and thermal dose mapping can be retrieved during RF ablation in vivo on sheep [START_REF] Ozenne | Improved cardiac magnetic resonance thermometry and dosimetry for monitoring lesion formation during catheter ablation[END_REF][START_REF] Toupin | Feasibility of real-time MR thermal dose mapping for predicting radiofrequency ablation outcome in the myocardium in vivo[END_REF]. A second potential issue is related to potential the local contrast changes on magnitude images associated with temperature increase. This can be attributed to change in T1 and/or reduction of signal intensity due to intra-voxel dephasing (in case of non-uniform temperature distribution within a voxel). As a result, motion compensation relying on algorithms based on the conservation of luminance on magnitude images may be erroneous. This effect was addressed in a previous work [START_REF] Toupin | Combination of principal component analysis and optical-flow motion compensation for improved cardiac MR thermometry[END_REF] using a dedicated algorithms providing robust thermometry against local intensity variations on magnitude images. A third potential issue is out-of-plane motion during breathing. However, in patients under resting condition, the risks of out-of-plane motion during regular breathing appears limited. In the present study, no catheter was inserted into the heart. However, in presence of a catheter, the slice orientation may be refined using local contrast changes provided by the catheter or using catheter-tracking methods.

Accounting for region variability is important in the context of cardiac ablation. Our results indicate some regional variability in particular in the free wall (heart-lung interface, Figure 6E) or at the apex (heart-liver-lung interface, Figure 5E andS1). In these areas, phase wraps, characterized by fast fluctuations of the phase, challenge motion compensation, phase modeling and phase subtraction and impact both accuracy and precision of the method. In particular, the unwrapping process which performs 2π addition or subtraction is impacted by both the presence of fast temporal and spatial fluctuation variation and imperfect motion compensation in presence of sudden physiological changes. Fortunately, these areas are limited and additional strategies could be set up to improve phase quality (modification of the orientation plane or of the shim box). A further study with a larger cohort should be performed to quantify precisely the regional variability using for instance a bullseye representation of the left ventricle. Images were only acquired in coronal orientation due to the initial long duration cardiovascular magnetic resonance exam. Alternative planes could have been chosen as long as the major component of the respiratory motion is set into the imaging plane for future motion compensation.

Alternatively, respiratory navigated acquisitions may be performed to compensate for irregular breathing, as suggested in Toupin et al [START_REF] Toupin | Feasibility of real-time MR thermal dose mapping for predicting radiofrequency ablation outcome in the myocardium in vivo[END_REF] allowing the visualization of natural cardiac perspective (SAX, LAX). Types of arrhythmia are numerous and repercussions on image acquisition and arrhythmia rejection tools should be reported for a larger population to access its effectiveness. Moreover, more sophisticated image processing strategies may also be envisioned, such as the creation of reference data sets for each cluster in order to reduce the effects of susceptibility variations associated with cardiac motion in case or irregular cardiac rhythm (e.g. see results on patient #14 on Figure 8D andE). Such a promising strategy will require further developments.

Finally, none of the patients had a sustained type of arrhythmia such as atrial fibrillation. And we did not account for the impact of irregular rhythms that may arise when ablating the left ventricle. The proposed post processing that select images for rejection is predicated on a series of ECG patterns during the training phase of imaging, it is possible that duration of RR intervals will be significantly different during catheter ablation making the RR-duration-based cluster analysis less effective. Again, further studies will be needed to test this thermometry sequence during such AF or during catheter ablation and this was beyond the scope of the current study. It is also important to keep in mind that the small temporal fluctuation voxel is not critical regarding the temperature increase during ablation (>20°C). Moreover, it must be noticed that in practice, temperature evolution during heating or cooling is spatially and temporally continuous, thus, the presence of sporadic fluctuations can be tolerated.

Conclusion

Real-time CMR thermometry provides real-time visualization of temperature distribution within the myocardium. This study indicates that this technique is feasible in patients with a good precision even in presence of some degree of arrhythmia. Two strategies were implemented to improve the robustness of the pipeline to the presence of ectopic beats and sudden respiratory motion. In perspective of future MR-guided radiofrequency ablation, such a method may provide an online display of lesion formation during the procedure with the objective of creating a transmural lesion of the arrhythmogenic substrate by coagulation necrosis. It is expected to help in optimizing radiofrequency delivery by preventing insufficient (no lesion) or excessive heating (myocardium injuries) thus improving treatment efficiency and patient safety.

Figure and table legends

Table 1: Patient information. The total scan duration and the total time where temperature mapping is available during acquisition is indicated. The RR interval are reported in milliseconds and beats per minute and the corresponding number of beats extracted from ECG processing and finally the arrhythmia prevalence in beat numbers and percentages. Prevalence of arrhythmia is high and the total acquisition time is 3 min 35 s. This patient #6 had RR duration intervals of 853 ± 286 ms. We note the presence of phase wraps (indicated by the orange arrows) affecting the temperature mapping. In (E), the temporal standard deviation of the temperature is overlaid on the averaged registered magnitude images from (C).Reference phase image is shown in (D). Blurring on magnitude images inside the cavity with not well-defined edges may be due to potential voxels misregistration on the myocardium. The largest cluster in blue has been used for the calculation. Temporal standard deviation of temperature was calculated without rejection of the temperature maps (E left), after rejection using the Gphase metric (E middle) and after rejection combining the Gphase metric and clustering (E right). Without rejection, the presence of phase wraps affects locally the temporal standard deviation of temperature in the liver and in the left ventricle. After removing four dynamic acquisitions using the Gphase metric, the temporal standard deviation of temperature is improved. Combining both rejection methods allows the decrease of the temperature uncertainty in the left ventricle and also offers some improvement in thinner structure like the aorta or the left atria (indicated by white arrows).

Supporting information legend:

Table 2: Temperature measurement precision. The temporal standard deviation (σT) and temporal mean (µT) for each patient and for each slice are reported. The number of voxels in the ROI surrounding the LV is reported as well as the fraction of myocardial voxels displaying σT below 2°C and higher than 10°C (removed outliers). 150 voxels located in the ROI on the LV were randomly selected using the rand function from MATLAB. Voxels temperature uncertainty (σT) <1, <2, <5, <10 °C and >10°C were plotted in dark blue, pale salmon pink, orange and red, respectively (using the same color code as Figure 6C). The red curves (σT >10°C) were not used in the data analysis. The spatial location of the voxels is indicated using green dots (σT <10°C) and red dots (σT >10°C) in the average magnitude images for each slice. Patient #1 to #5 are shown as well as the worst-case scenario in patient #12.

Figure 1 :

 1 Figure 1:A. Categorization of beats. The ECG was recorded using standard 3-lead acquisition. A home-made algorithm was designed to identify the QRS peaks with an adaptive threshold, in order to associate ECG signal with acquisition of MR thermometry data. A categorization of beats based on the duration of two successive RR intervals was used to quantify the prevalence of arrhythmia for each patient. The RR duration of the loading beat RRi-1 was plotted versus the RR duration of the current beat RRi to allow clustering of the different beat types. B1. Calculation of the gradient change of the phase (Gphase)image. A Sobel filter for edge detection is applied in the registered phase image (a): two 3×3 convolution kernels (to calculate spatial derivatives in x and y direction) are applied on registered phase image. The resulting convolutions named Gx and Gy (b) were combined together to find the absolute magnitude (c) of the gradient images defined as Psobel. A binary mask (e) based on the averaged registered magnitude image (d) is then multiplied to the Psobel image in order to remove the low SNR areas. The resulting image is named the gradient change of the phase Gphase image (f). B2. Comparison of the Gphase metric at two different dynamic acquisitions. (g) Phase images indicate the presence of a large motion of the heart between both acquisitions. While the phase image is flat in dynamic acquisition 216, the phase image in dynamic acquisition 221 shows numerous phase wraps in the myocardium and in the liver. (h)The registered magnitude and phase images were processed to calculate the Gphase images and finally the Gphase metric was computed. We note a huge increase of the Gphase metric in the presence of phase wraps.

Figure 2 :

 2 Figure 2: Data analysis of temperature mapping in sixteen patients. Box-and-whisker plots of the distribution of temperature σT (A) and offset µT (B) for all patients. Temperature data were taken on the left ventricle using all dynamic images excluding voxels where σT was higher than 10°C. Percentage of voxels with temperature uncertainty below 1, 2, 5, 10 °C is plotted in a diagram representation (C) and excluded voxels (σT>10°C) are indicated in red.

Figure 3 :

 3 Figure 3: Patient in sinus rhythm. A) 80 seconds of the ECG trace (red) along with the detected QRS peaks (blue triangles). The EPI gradient artifacts on the ECG are visible during the acquisition (indicated by the black arrows). B) Two-dimensional RR duration plot with one beat cluster. The delay between the QRS and the beginning of acquisition is indicated by a green boxe and for the time of acquisition of each slice is indicated by the red boxes. The patient #5 was in sinus rhythm with RR duration intervals of 1031 ± 37 ms. In (E), the temporal standard deviation of the temperature is overlaid on the averaged registered magnitude images from (C). Reference phase image is shown in (D). A total of 210 dynamic acquisitions were used for calculating temperature uncertainty σT and 240 dynamic acquisitions for averaging the registered magnitude images. Dotted lines show contours of the myocardium used for data analysis of temperature data.

Figure 4 :

 4 Figure 4: Patient with arrhythmia. A) 80 seconds of the ECG trace (red) with the detected QRS peaks (blue triangles). Arrhythmic events are visible at t = 33 s and t = 36 s. Acquisition window has been adjusted from mid-diastole to mid-systole in order to limit the impact of arrhythmia on MRthermometry. B) Two-dimensional RR duration plot with four beat clusters. Prevalence of arrhythmia is low with 20 beats over 255 beats for a total acquisition time of 3 min 48 s. This patient #7 had RR duration intervals of 904 ± 94 ms. In (E), the temporal standard deviation of the temperature is overlaid on the averaged registered magnitude images from (C). Reference phase image is shown in (D). Dotted lines show contours of the myocardium used for data analysis of temperature data.

Figure 5 :

 5 Figure 5: Patient with arrhythmia. (A) 80 seconds of the ECG trace (red) with the detected QRS peaks (blue triangles). (B) Two-dimensional RR duration plot with four beat clusters. Prevalence of arrhythmia is high with 89 beats over 257 beats for a total acquisition time of 3 min 34 s. This patient #12 had RR duration intervals of 842 ± 102 ms. In (E), the temporal standard deviation of the temperature is overlaid on the averaged registered magnitude images from (C). Reference phase image is shown in (D). The red arrow shows the myocardial infarction and the orange arrow show a wrapped phase and a low SNR in the lower right part of the myocardium. Blurring on magnitude images may be due to potential voxels misregistration on the myocardium in slices 4 and 5.

Figure 6 :

 6 Figure 6: Patient with permanent arrhythmia: (A) 80 seconds of the ECG trace (red) with the detected QRS peaks (blue triangle). (B) Two-dimensional RR duration plot with the five beat clusters.Prevalence of arrhythmia is high and the total acquisition time is 3 min 35 s. This patient #6 had RR duration intervals of 853 ± 286 ms. We note the presence of phase wraps (indicated by the orange arrows) affecting the temperature mapping. In (E), the temporal standard deviation of the temperature is overlaid on the averaged registered magnitude images from (C).Reference phase image is shown in (D). Blurring on magnitude images inside the cavity with not well-defined edges may be due to potential voxels misregistration on the myocardium.

Figure 7 :

 7 Figure 7: Detection of spontaneous motion inducing phase wraps. The magnitude and phase images are plotted with the corresponding ECG and Gphase metric of patient #14. Dynamic acquisitions are indicated in small boxes. Two events of interest are presented. The first event (in green) is related to arrhythmia and did not impact the magnitude and phase images. A second event (in purple) related to a spontaneous respiratory motion inducing modification of the surrounding magnetic field. Magnitude images in dynamic acquisitions 221 and 222 showed larger respiratory amplitudes. The corresponding phase images are characterized by the presence of numerous phase wraps. The calculation of the Gphase metric allows the identification of the phase wraps and the exclusion of the dynamic acquisitions 221 and 222.

Figure 8 :

 8 Figure 8: Impact of rejection of dynamic acquisition using the Gphase metric on temporal standard deviation of temperature. Statistics with and without rejection of dynamic acquisition are shown. Boxand-whisker plots of the distribution of temperature σT (A) and offset µT (B) over all patients. The data were retrieved from the same ROIs on the left ventricle. Percentages of voxels temperature uncertainty below 1, 2, 5, 10 °C are plotted in a diagram representation (C) and percentages of excluded voxels are indicated in red. Impact of rejection using the Gphase metric and the clustering method on patient #14 on slice 5: (D) clustering of RR duration intervals in five clusters (indicated in purple, blue, orange, yellow and red).The largest cluster in blue has been used for the calculation. Temporal standard deviation of temperature was calculated without rejection of the temperature maps (E left), after rejection using the Gphase metric (E middle) and after rejection combining the Gphase metric and clustering (E right). Without rejection, the presence of phase wraps affects locally the temporal standard deviation of temperature in the liver and in the left ventricle. After removing four dynamic acquisitions using the Gphase metric, the temporal standard deviation of temperature is improved. Combining both rejection methods allows the decrease of the temperature uncertainty in the left ventricle and also offers some improvement in thinner structure like the aorta or the left atria (indicated by white arrows).

Figure S1 :

 S1 Figure S1: Regional analysis of the temperature uncertainty: (A) segmentation of the left myocardium into twelve angular sectors of equal aperture, centered on the barycenter of the ROI. (B) For each group of patient, the mean value of the temporal standard deviation of the temperature is displayed in each of the twelve segments of left myocardium. (C) The percentage of removed voxels is indicated using the same graphical representation.

Figure S2 :

 S2 Figure S2: Histogram distribution of temporal standard deviation (σT): the histogram distribution of σT is displayed for each patient and for each slice. A log scale ranging from 0 to 100 is used along the x axis. Two vertical lines corresponding σT =1°C (black line) and σT =10°C (red line) are indicated. Points located on the right side of the red line (σT >10°C) correspond to the voxels that have been removed from the data analysis. The temperature uncertainties (σT) mean ± SD over the ROI on the LV is indicated after outliers' removal.

Figure

  FigureS3and S4: Temperature as function of time. 150 voxels located in the ROI on the LV were randomly selected using the rand function from MATLAB. Voxels temperature uncertainty (σT) <1, <2, <5, <10 °C and >10°C were plotted in dark blue, pale salmon pink, orange and red, respectively (using the same color code as Figure6C). The red curves (σT >10°C) were not used in the data analysis. The spatial location of the voxels is indicated using green dots (σT <10°C) and red dots (σT >10°C) in the average magnitude images for each slice. Patient #1 to #5 are shown as well as the worst-case scenario in patient #12.

Figure S5 :

 S5 Figure S5: Variation of the Gphase metric as function of the dynamic acquisition: for each acquisition the Gphase metric is plotted (black line). A rejection strategy was implemented based on the temporal variation of the Gphase metric. The threshold of the rejection criteria was defined as six times the standard deviation of the Gphase metric during the learning step. During the following image acquisition, all phase acquisitions with a metric superior to the threshold (black dotted lines) were removed from the temperature calculation (red dots).

Table 1

 1 presence of arrhythmia but quantification NA due to a technical problem during ECG processing

			Scan	T. mapping RR interval		Total	Arrhythmia
								beats nb
	Patient Age Sex Duration	duration	ms	bpm	beat nb	(%)
	1	76	M 3 min 36	3 min 02	NA	NA	NA	*
	2	50	M 3 min 37	3 min 02	NA	NA	NA	-
	3	44	F 4 min 44	3 min 58	NA	NA	NA	*
	4	62	M 3 min 34	3 min 00	856 ± 53	70.1 ± 4.4 253	5 (2%)
	5	75	M 4 min 16	3 min 36	1031 ± 37 58.2 ± 2.1 250	-
	6	42	M 3 min 35	3 min 00	853 ± 286 70.4 ± 23.6 273	273 (100%)
	7	71	M 3 min 48	3 min 10	904 ± 94	66.4 ± 7.3 255	20 (8%)

*

Acknowledgements

We would like to thank Lorena Sanchez Blanco, Vincent Dejarnac and Xavier Pillois for helping with the ethics approval. We also would like to thank Remi Dubois for helping with the ECG processing.

Funding

This work has been implemented through the project P3-Stroke, an EIT Health innovation project. EIT Health is supported by EIT, a body of the European Union. This work received the financial support from the French National Investments for the Future Programs ANR-10-IAHU-04 (IHU Liryc), Laboratory of Excellence ANR-10-LABX-57 (TRAIL) and the research programs ANR-11-TecSan-003-01 (TACIT) and Equipex ANR-11-EQPX-0030 (MUSIC) and ANR-17-CE19-000 (CARTLOVE).

Availability of data and material

Not applicable.

Consent for publication

The protocol and CCPRB approval include permission to use de-identified patient images in scientific communications.

Competing interests

ST, AV and MLC are employees of Siemens Healthcare. No other authors have competing interests to declare.

Authors' contributions

VO conducted the study, analyzed the data and did the primary writing of the manuscript. VO, PB, ST, and MLC participated to the design and to the optimization of the MR-thermometry pulse sequence.

VO, ST and BDS contributed to the design of the image processing pipeline including compensation of motion and associated artifacts. VO, AV, PB, and PJ contributed to the ECG analysis. HC and PJ were responsible for the human study and assisted in the interpretation of the results. BQ conceived of the study, and participated in its design and coordination. VO, PB, BDS, ST, AV, MLC, PJ, HC and BQ drafted the manuscript. All authors read and approved the final manuscript.


Sylvain Boust

Houssein El Dirani

François Duport

Laurène Youssef

Yannick Robert

Larrue Alexandre

Camille Petit-Etienne

Eric Vinet

Sébastien Kerdiles

Erwine Pargon

Compact optical frequency comb source based on a DFB butt-coupled to a silicon nitride microring

Keywords: ²Optical frequency comb, photonic integrated circuits, silicon nitride, semiconductor lasers /e²

Q Q i Q e A.


Laura Piqué-Fandiño

Miquel Hernández

Fernando V Ramirez

BIRTH SEASONALITY IN THE BAKA PYGMIES

RESULTS

There is a significant birth pattern in the Baka Pygmy

  The Baka Pygmies in Moange le Bosquet are a semi-nomadic population living in equatorial Africa with an economy based on hunts and gathers. During the year they develop different economical activities. Pygmies share a common ancestor and splitted from Bantu-speaking populations at around 60,000 years BP [5] and present a particular growth, slower rate of growth during infancy for Baka, which explains their pygmy phenotype. [6-8] Patterns in births have been observed in all kind of human populations [1-3], but less studies have been focused on smaller hunter-gatherer populations. The patterns differ across populations, and are usually explained by environmental variables, in combination with social and economic factors [4]. The two main weather variables that have been related with fluctuations in births are temperature (behavioural and biological) and rainfall (farming seasons, scarcity of food). Also cultural factors such as marriages, labour migration, religion festivities or holidays, among others, have been related to birth patterns in western populations. To assess the presence of a specific pattern of birth seasonality in the Baka Pygmies from Moange le Bosquet, Cameroon. To determine the factors that influence the pattern of birth in this population.The slopes graph represents the increase or decrease of births from one month to the following. We assigned 1 if the number of births increased, -1 if it decreased and 0 if it did not the line is in the negative area, the number of births decrease, and it increases in the positive area. If it crosses the 0 going from negative to positive, it's a minimum, and from positive to negative, it's a maximum. The method was designed to minimize the effect of long term change in birth number over the years fluctuations, by looking only whether birth numbers increase or decrease between months.

  population of Moange le Bosquet (P<0,05) Months of rainfall peaks (May and October), correspond with peaks in birth (Fig. 1). The mechanisms through which rainfall affects births are unknown, but are not likely to be related to the farming cycle or food scarce because of the Baka's lifestyle. Year-round oscillation in temperature is too small (2,7º C) to affect conceptions or births. At least through the mechanisms described in previous literature. Cultural factors like marriages or holidays are ruled out as factors causing the birth seasonality because of the characteristics of the population. The economical activities developed by the Baka, such as hunting seasons and collecting wild mango, could affect season of birth. 1. Cogwill, U. M. (1996). The Season of Birth in Man. Man, 1 (2), 232-240 2. Lam, D. A. & Miron, J. J. A. (1991). Seasonality of Births in Human Populations. Biodemography and Social Biology, 38 (1-2), 57-78 3. Leslie, P. W. & Fry, P. H. (1989). Extreme Seasonality of Births among Nomadic Turkana Pastoralists. American Journal of Physical Anthropology, 79 (1), 103-115. 4. Bailey, R. C., Jenike, M. R., Elliston, P. T., Bantley, G. R., Harrigan, M. A., & Peacock, N. R. (1992). The Ecology of Birth Seasonality among Agriculturalists in Central Africa. Journal of Biosocial Science, 24 (3), 393-412. 5. Verdu, P. et al. (2009). Origins and Genetic Diversity of Pygmy Hunter-Gatherers from Western Central Africa. Current Biology, 19, 312-318. 6. Ramirez Rozzi, F. V., Koudou, Y., Froment, A., Le Bouc, Y. & Botton, J. (2015). Growth pattern from birth to adulthood in African pygmies of known age. Nature Communications, 6:7772. 7. Ramirez Rozzi, F. V. (2016). Diversity in tooth eruption and life history in humans: illustration from a Pygmy population. Scientific reports, 6, 27405. 8. Ramirez Rozzi, F. V. (2018). Reproduction in the Baka pygmies and drop in their fertility with the arrival of alcohol. Proceedings of the National Academy of Sciences of the United States of America, 115 (27), 1-9.

Figure 1 .

 1 Figure 1. Monthly slope showing the increase /decrease of number of births whithin months. Monthly precipitation (mm) and temperature (ºC) means from 1961-1990 from Batouri weather station. Seasons of different economical activities carried out by the Baka.

  Baka Pygmy individual, taken by Dr. Ramirez Rozzi with informed consent from the individual.

Acknowledgements

We thank Dr. Samuel Pavard for his great contribution to the study, as well as his advice and recommendations.