Bias Introduced by Multiple Head Coils in MRI Research: An 8 Channel and 32 Channel Coil Comparison

Neuroimaging MRI data in scientific research is increasingly pooled, but the reliability of such studies may be hampered by the use of different hardware elements. This might introduce bias, for example when cross-sectional studies pool data acquired with different head coils, or when longitudinal clinical studies change head coils halfway. In the present study, we aimed to estimate this possible bias introduced by using different head coils to create awareness and to avoid misinterpretation of results. We acquired, with both an 8 channel and 32 channel head coil, T1-weighted, diffusion tensor imaging and resting state fMRI images at 3T MRI (Philips Achieva) with stable acquisition parameters in a large group of cognitively healthy participants (n = 77). Standard analysis methods, i.e., voxel-based morphometry, tract-based spatial statistics and resting state functional network analyses, were used in a within-subject design to compare 8 and 32 channel head coil data. Signal-to-noise ratios (SNR) for both head coils showed similar ranges, although the 32 channel SNR profile was more homogeneous. Our data demonstrates specific patterns of gray and white matter volume differences between head coils (relative volume change of 6 to 9%), related to altered image contrast and therefore, altered tissue segmentation. White matter connectivity (fractional anisotropy and diffusivity measures) showed hemispherical dependent differences between head coils (relative connectivity change of 4 to 6%), and functional connectivity in resting state networks was higher using the 32 channel head coil in posterior cortical areas (relative change up to 27.5%). This study shows that, even when acquisition protocols are harmonized, the results of standardized analysis models can be severely affected by the use of different head coils. Researchers should be aware of this when combining multiple neuroimaging MRI datasets, to prevent coil-related bias and avoid misinterpretation of their findings.


INTRODUCTION
Large multicenter data samples are increasingly used to establish and reproduce MRI neuroimaging findings. Although pooling MR imaging data contributes to increased study power, the reliability and results from such studies may be compromised by the use of different hardware elements (Cannon et al., 2014;Smith and Nichols, 2018). For example, changing head coils during a longitudinal study, or combining cross-sectional data acquired with different head coils may introduce a coil-related bias (Focke et al., 2011;Pardoe et al., 2016). Studies that depict the effect of using multiple head coils are currently limited to analysis of T1 weighted imaging data (Focke et al., 2011;Pardoe et al., 2016), demonstrating a difference in gray matter volume (Focke et al., 2011), and cortical thickness (Pardoe et al., 2016). Still, the effects on other quantitative MRI variables, for instance obtained with diffusion tensor imaging or resting state functional MRI, are unknown, but highly important for multicenter or longitudinal studies using different types of hardware. Identification of the brain regions affected by coilrelated bias is essential, not only to increase awareness, but more importantly to avoid misinterpretation of results from studies using multiple MRI hardware elements (Li et al., 2018). In the present study, we aimed to estimate and depict the impact of using different receive-only phased array head coils (8 channel head coil and 32 channel head coil) on T1 weighted, DTI and resting state functional MRI data using a within-subject design in a large cohort of cognitively healthy subjects.

Study Procedure and Participants
For the present study, we included 77 cognitively healthy participants who underwent MRI of the brain on a 3Tesla Philips Achieva scanner (Philips Medical Systems, Best, Netherlands) at the Leiden University Medical Center, Leiden, Netherlands. The MRI protocol contained T1 weighted, DTI and resting state functional MR images, acquired with both an 8 channel SENSE head coil (8CH) and an 32 channel SENSE head coil (32CH) within one MRI session (coil geometry is displayed in Supplementary Figure S1). During acquisition, optimal image quality was obtained by using the incorporated "constant level appearance" (CLEAR) inhomogeneity correction algorithm on the scanner. We display one raw dataset for all sequences from both coils from a representative healthy participant in Supplementary Figure S2.
Cognitively healthy participants were included in the context of the prospective longitudinal frontotemporal dementia risk cohort (FTD-RisC) in which families with autosomal dominant inherited FTD gene mutations are followed using standardized assessment protocols including an MRI of the brain every year, as described previously (Dopper et al., 2014;Papma et al., 2017). To confirm cognitively healthy status of all participants, Mini Mental State Examination [MMSE (Folstein et al., 1975)] and the Frontal Assessment Battery [FAB (Dubois et al., 2000)] are reported as cognitive screening measures and the Neuropsychiatric Inventory [NPI-Q (Cummings et al., 1994)] and Frontotemporal Dementia Rating Scale [FRS (Mioshi et al., 2010)] are reported as behavioral screening questionnaires.

Signal to Noise Assessment
For signal to noise ratio's (SNR) assessment, we acquired proton density weighted single-slice images with one noise-only image, i.e., without radiofrequency pulses, in one healthy volunteer, following the procedures described by Wiggins et al. (2006). For both coils, we acquired the images in the axial, coronal and sagittal direction, with the following parameters: Repetition time (TR) = 200 ms, echo time (TE) = 3.1 ms, field of view (FOV) = 220 × 220 × 3 mm, flip angle = 20 • (for noise scan 0 • ), slice thickness = 3 mm, voxel size 0.85 × 0.85 × 3.0 mm, number of averages = 10, acquisition time = 9 min and 24 s.

Diffusion Imaging
Diffusion imaging was performed in 60 non-collinear gradient directions using single shot echo planar imaging. The phase encoding direction was anterior to posterior for both coils. The following parameters were used for the 8CH coil: 60 b = 1000 s/mm 2 , TR = 8250 ms, TE = 80 ms, FOV = 256 × 208 × 140 mm, flip angle = 90 • , slices = 70, voxel size = 2 × 2 × 2 mm, SENSE = 2.0, one b = 0 s/mm2 acquisition, scan time = 8 min and 48 s, 2 signal averages. For the 32CH coil, we increased the number of slices to 80, to have sufficient coverage to include the cerebellum in the imaging volume for our longitudinal FTD-RisC study (Papma et al., 2017). As a result, the TR for the 32CH coil was 9245 ms, FOV was 256 × 232 × 160 mm, and acquisition time increased to 9 min and 52 s with a SENSE factor of 2.3. Other parameters were identical between both head coils.

MRI Processing
Before image preprocessing and analysis, we checked the scans thoroughly for image quality and the presence of artifacts. Data processing and statistical analyses were carried out using Functional Magnetic Resonance Imaging of the Brain Software Library (FSL) version 5.0.8. (Jenkinson et al., 2012).

Signal to Noise Ratio
For both the 8CH coil and 32CH coil, we isolated an average signal image and one noise-only image for each orientation plane. We subdivided the signal and noise images into nonoverlapping regions of interest (ROI) of 16 by 16 voxels. Next, we calculated the mean signal of the ROI using the averaged signal image, and the standard deviation of the noise of the ROI using the noise-only image. Since the noise images were amplitudereconstructed, the measured standard deviation was corrected for the Rician noise distribution (Haacke et al., 1999). Ultimately, for each ROI, SNR was calculated according to: Third, we translated the SNR ROI matrices into color-coded maps, in order to visualize the SNR distribution throughout the brain for both coils.

Structural Imaging
To assess the influence of the head coil on gray and white matter volume measurements, we applied the standard voxelbased morphometry (VBM) pipeline as implemented in FSL.
Preprocessing of the T1 weighted images included brain extraction followed by radiofrequency (RF) inhomogeneity correction, tissue segmentation and realignment to Montreal Neurological Institute (MNI) standard space using non-linear registration. We performed quality control to ensure good brain extraction, that was not different between both head coils. Next, FMRIB's Automated Segmentation Tool (FAST) was used for correction for spatial intensity variations, also known as bias field or RF inhomogeneity, and segmentation of the T1 weighted images (Zhang et al., 2001). The corrected, segmented gray matter images were re-registered non-linearly to a study-specific template with a balanced set of 8CH and 32CH coil images. The registered partial volume images were divided by the Jacobian of the warp field to correct for any local expansion or contraction. An isotropic Gaussian kernel with a sigma of 3 mm, which corresponds to a full width at half maximum kernel (FWHM) of approximately 7 mm, was used to smooth the gray matter segmentations. We also applied the VBM processing pipeline to the white matter segmentations, resulting in registered, corrected and smoothed white matter images for voxel-wise analyses.

Diffusion Imaging
Diffusion scans were corrected for motion artifacts and eddy currents by alignment to the b = 0 image using the FMRBIB Diffusion Toolbox. The tensor was fitted each voxel to create fractional anisotropy (FA) and mean diffusivity (MD) images. Subsequently, we applied standard tract-based spatial statistics (TBSS) as implemented in FSL (Smith et al., 2006). FA images were aligned to standard space using non-linear registration and averaged into a mean FA image. To create a study-specific FA mask, we thresholded the mean FA image with a minimum value of FA ≥ 0.2. This binarised FA mask was applied to voxel-wise comparisons of FA and MD between coils.

Resting State Functional MRI
For preprocessing of the resting state fMRI scans, we applied the fMRI Expert Analysis Tool (FEAT) as implemented in FSL, consisting of motion correction with MCFLIRT and spatial smoothing with a kernel of 6 mm FWHM. The datadriven independent component analysis (ICA) based Automatic Removal of Motion Artifacts (ICA-AROMA) approach was used to identify and remove noise components from the resting state fMRI data (Pruim et al., 2015). After denoising, high pass temporal filtering was performed with a cut-off frequency of 0.01 Hz. The functional resting state images were registered to the corresponding T1 weighted images using Boundary-Based Registration and were subsequently registered to the 2 mm isotropic MNI standard space using non-linear registration with a warp resolution of 10 mm. Voxel-based functional connectivity was studied in a standardized manner using the eight standard Beckmann resting-state functional networks of interest (Beckmann et al., 2005), i.e., the medial and lateral visual system network, the primary auditory network -also known as the salience network -, the sensory motor network, the default mode network, the executive control network and the left and right dorsal visual processing stream networks. To further account for noise, white matter and CSF templates were included in the analyses as regressors. Functional connectivity of each network of interest was calculated using dual regression, as previously described (Hafkemeijer et al., 2017). In short, the eight standard resting state networks (Beckmann et al., 2005) were used as a reference. Voxel-based resting state functional connectivity was determined in terms of similarity of the BOLD fluctuations in the brain in relation to characteristic fluctuations in the standard resting state networks. With dual regression, individual time series were first extracted for each template, using the resting state networks, and the two additional white matter and cerebrospinal fluid maps, in a spatial regression against the individual fMRI data set (regression 1). The resulting matrices described temporal dynamics for each template and individual. Next, the temporal regressors were used to fit a linear model to the individual fMRI data set (regression 2), to estimate the spatial maps for each individual. This results in 3D images for each individual, with voxel-wise z-scores representing the functional connectivity to each of the predefined standard networks.

Statistical Analysis
For all analyses, we designed within-subject paired sample t-tests with each subject's mean effect to analyze head coil differences in gray and white matter volume, FA, MD, and resting state functional connectivity. We performed voxel-based non-parametric permutation testing (Nichols and Holmes, 2002) with 5000 permutations using FSL-randomize. The statistical threshold was set at p < 0.05, using threshold-free cluster enhancement (TFCE) technique and family-wise error (FWE) correction to correct for multiple comparisons across voxels. We quantified the severity of head coil differences by calculating effect sizes and percentage of change.

Voxel-Specific Scaling Factors
For our own longitudinal clinical study (Dopper et al., 2014;Papma et al., 2017), we aimed to create voxel-based scaling factors to correct for the use of different head coils. For the T1weighted, DTI and resting state fMRI images, we calculated and validated voxel-specific scaling images. Procedures are described in detail in the Supplementary Material. In short, we separated our sample into a template dataset (n = 39) and validation dataset (n = 38), matched for age and sex. For the template dataset, we merged and averaged the images into a mean 8CH coil image and 32CH coil image. The averaged 32CH coil image was divided by the averaged 8CH coil image, resulting in a voxel-based scaling factor. We reduced noise by applying a median filter with a kernel of 5 mm. Next, the 8CH coil images from the validation set were multiplied with the voxel-based scaling factor, equalizing the 8CH coil images to the signal intensity of the 32CH coil images. We repeated previous described statistical analysis on the validation set to complete verification of the scaling factor.

Participants
In total 77 participants were included in this study (Table 1).
Cognitive and behavioral screening tests confirmed a cognitively healthy status of all participants.

Signal to Noise Ratio
Visualization of the SNR ROI matrices revealed a quite homogeneous distribution of SNR throughout the brain using the 32CH coil, with the highest SNR in posterior cortical areas of the brain. SNR of the 8CH coil was highest in the frontal lobe, but dropped in central and medial areas (Figure 1).

Structural Imaging
Quality control showed no differences in brain extraction of the T1 weighted images between the 8CH and 32CH coil. Gray matter volumes obtained with the 32CH coil were larger than obtained with the 8CH coil (p FWE < 0.05, effect size = 2.096, increase = 6.2%), particularly in the middle and inferior frontal lobe, the superior and middle temporal lobe, the anterior insular cortex, the temporo-parietal junction, the paracingulate, and the cuneus (yellow areas in Figure 2A). Gray matter volumes appeared smaller using the 32CH coil (p FWE < 0.05, effect size = 2.571, decrease = 8.9%) in frontal and deeper cerebral areas, such as the medial temporal lobe, medial frontal lobe, basal ganglia, posterior insular cortex, anterior cingulate, superior frontal cortex, occipital lobe, and the cerebellum (blue areas in Figure 2A). In the white matter, we found larger white matter volumes in subcortical and posterior cortical regions using the 32CH coil compared with the 8CH coil (p FWE < 0.05, effect size = 1.951, increase = 8.5%; see yellow areas in Figure 2B). White matter volumes were smaller using the 32CH coil in frontotemporal regions (p FWE < 0.05, effect size = 1.637, decrease = 6.1%; see blue areas in Figure 2B).

Diffusion Tensor Imaging
We found higher FA values (p FWE < 0.05, effect size = 2.197, increase = 5.7%) with the 32CH coil compared with the 8CH coil in all tracts of the right hemisphere and frontal tracts of the left hemisphere, such as the forceps minor, anterior parts of the uncinate fasciculus, anterior thalamic radiation, inferior fronto-occipital fasciculus, and superior longitudinal fasciculus (yellow areas in Figure 3A). On the contrary, we found lower FA values using the 32CH coil (p FWE < 0.05, effect size = 2.038, decrease = 5.0%) in part of the left-sided posterior tracts, such as the forceps major, the corticospinal tract, the inferior longitudinal fasciculus, and the central and posterior parts of the anterior thalamic radiation (blue areas in Figure 3A). MD was lower using the 32CH coil compared to 8CH coil (p FWE < 0.05, effect size = 1.871, decrease = 4.6%) in the entire right hemisphere, and some tracts of the left hemisphere located in the prefrontal and the occipital lobe (red-yellow areas in Figure 3B). Using the 32CH coil, MD was higher (p FWE < 0.05, effect size = 1.952, increase = 4.7%) in all tracts of the left hemisphere, except for the prefrontal and occipital lobe (blue areas in Figure 3B).

Resting State Functional MRI
Resting state functional connectivity was predominantly higher when using the 32CH coil compared with the 8CH coil (yellow areas in Figure 4), between the medial visual network and the lateral occipital cortex, calcarine cortex, and lingual gyrus (p FWE < 0.05, effect size = 0.392, increase = 13.5%; Figure 4A), between the lateral visual network, the lateral occipital cortex and the occipital pole (p FWE < 0.05, effect size = 0.637, increase = 27.5%; Figure 4B), between the default mode network and the lateral occipital cortex (p FWE < 0.05, effect size = 0.505, increase = 9.4%; Figure 4E), and between the dorsal visual stream networks and regions of the lateral occipital cortex (right: p FWE < 0.05, effect size = 0.368, increase = 13.6%; Figure 4G and left: p FWE < 0.05, effect size = 0.391, increase = 12.04%; Figure 4H). Functional connectivity was lower with the 32CH coil in the executive control network and a small area in the frontal pole (p FWE < 0.05, effect size = 0.608, decrease = 23.1%; blue area in Figure 4F). No differences in functional connectivity between both coils were found neither in Frontiers in Neuroscience | www.frontiersin.org FIGURE 1 | Color-coded SNR brain maps in the sagittal, axial and coronal direction for 8 channel head coil (8CH) and 32 channel head coil (32CH). Colorbar represents SNR values, ranging 0-85. the auditory, or salience network (Figure 4C), nor the sensorymotor network (Figure 4D).

Voxel-Specific Scaling Factors
Results are described in detail in the Supplementary Material. For T1-weighted imaging, after applying the scaling factor, VBM analyses showed almost complete removal of the head coil differences throughout the brain for both GM and WM (Supplementary Figure S3). Validation of the DTI scaling factors showed a successful harmonization of the FA images, removing all significant coil differences. For MD, head coil differences were reduced, except for some small areas at the forceps major and right thalamus (Supplementary Figure S4). For resting state analyses, validation of the scaling factors showed extensive reduction of coil differences in the medial and lateral visual networks, and complete removal of all significant differences in the default mode network, executive control network (salience network), and the left and right dorsal visual stream networks (Supplementary Figure S5).

DISCUSSION
In this study, we demonstrated that quantitative results of standard processing pipelines for T1-weighted MRI, diffusion tensor imaging and resting state fMRI will be severely affected by the use of different head coils. Paired-wise group analyses between coils revealed different patterns of gray and white matter volume, white matter connectivity and functional connectivity.  Voxel-based morphometric analysis of T1-weighted imaging data, revealed smaller apparent gray matter volume for the 32CH coil in the outer frontal, temporal and parietal cortex, the inner cerebellum, the precuneus, and posterior cingulate cortex compared to the 8CH coil. Gray matter volume was larger using the 32CH coil in the occipital lobe, the central layer of the frontal, temporal, parietal cortex, the peripheral layer of the cerebellum, and subcortical areas. A less extensive but similar pattern of gray matter volume differences was previously found using two identical scanners with, respectively, an 8-channel and a 12channel head coil (Focke et al., 2011). Compared to gray matter VBM, we found an opposite pattern of head coil differences for the white matter VBM, meaning that in areas where gray matter volume was larger, white matter volume was smaller and vice versa. The visual overlap of gray and white matter differences may be a result of the spatial smoothing (7 mm). When the contrast between gray and white matter is unclear, a higher level of smoothing is necessary to account for the uncertainties in partial volume estimation. Our results indicate that differences in gray/white matter contrast of the images leads to differences in tissue classification during segmentation (Zhang et al., 2001;Tohka, 2014). Three factors could have led to altered image contrast: (1) despite using inhomogeneity correction prior to segmentation, the SNR of the 8CH and 32CH coil has affected the probabilities of gray and white matter (Wiggins et al., 2006;Reiss-Zimmermann et al., 2013). For example, in our case, higher SNR in posterior cortical areas for the 32CH coil increased the probability of a voxel being white matter, and the same principle may explain the results in the frontal areas of the 8CH coil, (2) both coils have differences in signal distribution due to the coil geometry (Blamire, 2008) and the CLEAR algorithm does not fully correct for these (Yun et al., 2007), or (3) there has been a difference in the effective b1 distribution (Marques et al., 2010). As we demonstrate here, differences in image contrast in structural images could pose serious problems for studies combining MRI hardware elements and need to be equalized before tissue segmentation and partial volume estimation to prevent methodological errors.
We demonstrated that for DTI, measured FA was higher for nearly the entire white matter when comparing the 32CH coil with the 8CH coil, except for parts of the left temporal and parietal lobe. For MD, the results were hemispherical dependent, i.e., the 32CH coil showed higher MD compared to the 8CH coil in the white matter of the left hemisphere but lower MDs in the right hemisphere. The asymmetry in especially MD metrics was unexpected and inexplicable. Previous research demonstrated that MD may have more variance and be less reproducible than FA, even within sites (Fox et al., 2012;Jovicich et al., 2014;Helmer et al., 2016). Interestingly, the unexpected pattern in diffusion metrics was different from the SNR profiles of both coils and volumetric results. Therefore, we emphasize that SNR profiles and pattern of changes that occur in T1-weighted data cannot be translated to DTI data. Our results underline previous studies that already demonstrated the possible pitfalls of pooling diffusion data (Pagani et al., 2010;Vollmar et al., 2010;Takao et al., 2012;Venkatraman et al., 2015). Harmonization methods have been investigated in a number of studies, in attempt to overcome the problems with pooled DTI data, but sufficient harmonization has proven to be difficult (Mirzaalian et al., 2015;Venkatraman et al., 2015;Pohl et al., 2016;Fortin et al., 2017). Note that we slightly increased our FOV to allow coverage of the cerebellum for the 32CH coil data, which also increased TR, and could have influenced our results. This is, however, expected to be a very minor effect, since both TRs are significantly longer than five times the T1 of the white matter.
Resting state fMRI analyses showed increased functional connectivity between multiple posterior located networks and brains areas using the 32CH coil, corresponding to the SNR profile. Decreased functional connectivity using the 32CH coil was found between a small frontal area and the executive control network. Previous studies demonstrated around 10 percent variability between scanner hardware (Costafreda et al., 2007;Friedman et al., 2008;Kaza et al., 2011), but others could not detect significant differences in resting state networks between the 8CH and 32CH Philips head coils (Paolini et al., 2015). As our current sample size (n = 77) is larger than previous studies (n = 26), we assume that the increased power in our study allowed for detection of coil differences in resting state functional connectivity. Concordant with the results from our study, past research demonstrated that the 32CH coil has particularly increased SNR in posterior cortical areas compared to coils with less channels, caused by a difference in coil geometry (Wiggins et al., 2006;Reiss-Zimmermann et al., 2013). Resting state functional MRI studies may benefit from increased SNR in the 32CH coil, especially when posterior cortical areas are of interest.
Since changing head coils during longitudinal clinical research might introduce bias, we aimed to create voxel-based scaling images (Supplementary Material) for our own longitudinal clinical study (Dopper et al., 2014;Papma et al., 2017). The unique within-subject design in a large cohort (n = 77) with stable acquisition parameters allowed for the use of voxel-specific information for VBM, DTI, and resting state fMRI scaling. We validated the use of the scaling factors on independent data and demonstrated that coil differences can be substantially reduced when using voxel-based scaling in all modalities (see Supplementary Material). Our voxel-specific scaling factors may be an interesting harmonization method for within-subject variation. In previous literature, many harmonization methods have aimed to equalize neuroimaging data, all with their own advantages and shortcomings (Fennema-Notestine et al., 2007;Keihaninejad et al., 2010;Takao et al., 2011;Chen et al., 2014;Griffanti et al., 2014;Salimi-Khorshidi et al., 2014;Feis et al., 2015;Mirzaalian et al., 2015;Venkatraman et al., 2015;Pardoe et al., 2016;Pohl et al., 2016;Fortin et al., 2017;Li et al., 2018). For example, additional complexity in statistical models may cause decreased sensitivity for the actual outcome of interest (Li et al., 2018). Neuroimaging studies combining multiple MRI head coils and other hardware elements should be aware of confounding factors and be committed to use robust, sensitive and validated methods to deal with these factors (Jack et al., 2008;Smith and Nichols, 2018;Boeve and Rosen, 2019). We aim to deepen our research into neuroimaging harmonization methods and our scaling factors in the near future.
Major strengths of this study are the large sample size and the within-subject study design of two protocols within one scanning session using similar acquisition parameters. Despite our best efforts, some confounding factors are essentially inevitable, such as scanner drift, re-positioning of the participants heads inside the different coils and use of head cushions (Littmann et al., 2006). Other factors that could have influenced the results in our study may be habituation of the subject and scanner warm-up, especially for the resting state fMRI, and the anisotropic voxel size of the T1 weighted sequence. We are aware that extrapolating the results from our study to other head coils or other vendors may be difficult. Instead, we emphasize that our study may be appreciated as increasing awareness for the variability and possible bias in quantitative MRI metrics in data originating from different hardware elements. This is especially important for clinical research, where data acquired with different hardware elements is increasingly pooled.
In conclusion, this study provides evidence that the results of standard analysis models are severely compromised when data from different head coils is combined, or when head coils are changed during longitudinal clinical studies, even though acquisition protocols are completely harmonized. Studies combining neuroimaging MRI data with multiple head coils or other MRI hardware elements should be aware that measurements of gray and white matter volume, white matter connectivity and functional connectivity will differ between head coils and should handle these confounding factors with caution.

DATA AVAILABILITY
The datasets for this manuscript are not publicly available yet as they are part of the ongoing longitudinal FTD-RisC study and the international GENFI cohort. Requests to access the raw MRI datasets should be directed to JS (j.c.vanswieten@erasmusmc.nl). Voxel-specific scaling factors from this study are available upon reasonable request to AH (a.hafkemeijer@lumc.nl).

ETHICS STATEMENT
The study has been carried out in accordance with the Declaration of Helsinki and has been approved by the Medical and Ethical Review Committees of the Erasmus MC University Medical Center, Rotterdam, Netherlands and the Leiden University Medical Center, Leiden, Netherlands. Written informed consent has been obtained from all participants.

AUTHOR CONTRIBUTIONS
JLP designed the study, collected and analyzed the data, and was involved in data interpretation and writing of the manuscript. YT, MO, JG, and JMPa contributed to data analyses, interpretation and writing of the manuscript. EE, JMPo, LJ, LM, ED, and JS were involved in data collection and writing of the manuscript. MB and SR contributed to data interpretation and writing of the manuscript. AH was involved in study design, analyses and interpretation of the data, and writing of the manuscript.