MR g-ratio-weighted connectome analysis in patients with multiple sclerosis

Multiple sclerosis (MS) is a brain network disconnection syndrome. Although the brain network topology in MS has been evaluated using diffusion MRI tractography, the mechanism underlying disconnection in the disorder remains unclear. In this study, we evaluated the brain network topology in MS using connectomes with connectivity strengths based on the ratio of the inner to outer myelinated axon diameter (i.e., g-ratio), thereby providing enhanced sensitivity to demyelination compared with the conventional measures of connectivity. We mapped g-ratio-based connectomes in 14 patients with MS and compared them with those of 14 age- and sex-matched healthy controls. For comparison, probabilistic tractography was also used to map connectomes based on the number of streamlines (NOS). We found that g-ratio- and NOS-based connectomes comprised significant connectivity reductions in patients with MS, predominantly in the motor, somatosensory, visual, and limbic regions. However, only the g-ratio-based connectome enabled detection of significant increases in nodal strength in patients with MS. Finally, we found that the g-ratio-weighted nodal strength in motor, visual, and limbic regions significantly correlated with inter-individual variation in measures of disease severity. The g-ratio-based connectome can serve as a sensitive biomarker for diagnosing and monitoring disease progression.

MS group compared with that in the control group (P = 0.04; Fig. 2B, Supplementary Table S1). These included 6 motor, 3 somatosensory, 8 visual, and 13 limbic regions. Connectivity strength was averaged across all connections within each subnetwork, yielding a subnetwork-averaged NOS and g-ratio value for each individual. Inter-individual variation in terms of these subnetwork-averaged NOS and g-ratio values significantly correlated with WM lesion volumes (NOS: P = 0.007, r = −0.69; g-ratio: P = 0.006, r = −0.70) but not with disease duration or expanded disability status scale (EDSS) score.
Nodal strength. NOS-weighted nodal strength tended to reduce in the MS group relative to the control group in brain regions such as the left inferior parietal lobe, left medial orbitofrontal lobe, left insular, left thalamus, right thalamus, right caudate, right precentral lobe, and right insular (uncorrected P < 0.05; Fig. 3A, Supplementary Table S2). However, no significant differences remained after false discovery rate (FDR) correction. In contrast, for the g-ratio-based connectome, several regions were found to have significantly increased nodal strength in the MS group relative to the control group (Fig. 3B, Table 1). The most prominent increases in nodal strength in the MS group were evident in the limbic regions, including the bilateral insular, bilateral amygdala, left temporal pole, and left accumbens (Fig. 3C, Table 1).
In addition, significant positive correlations were observed between the EDSS scores and nodal strength computed using the g-ratio-based connectome in the MS group, predominantly in motor, visual, and limbic regions (Fig. 4, Table 2). For the NOS-based connectome, nodal strength did not correlate with EDSS score or disease duration.   Supplementary Table S4). Abbreviations: NOS, number of streamlines; MS, multiple sclerosis.

Discussion
In this study, we evaluated whether the tract-averaged g-ratio can provide a measure of connectivity strength, which is more sensitive to the effects of demyelination than the conventional measures of connectivity strength based on NOS, for connectome mapping studies. The tract-averaged g-ratio increased in the MS group relative to the control group along the motor, somatosensory, visual, and limbic subnetworks. By contrast, NOS decreased in the MS group relative to the control group in mostly the same subnetworks. Significant correlations were detected between WM lesion volumes and mean NOS or g-ratio across the subnetworks associated with significant between-group differences. Further, the g-ratio-based connectome showed significantly increased nodal strength in patients with MS, whereas no such differences were evident for the NOS-based connectome. Significant positive correlations were observed in patients with MS in terms of EDSS scores and g-ratio-weighted nodal strength in motor, visual, and limbic regions. Overall, our results indicate that the g-ratio-based connectome can provide greater sensitivity to the effects of myelin pathology on connectivity than connectomes mapped using conventional measures such as NOS.
Significant increases in the g-ratio within networks including motor, somatosensory, visual, and limbic regions in patients with MS could be responsible for the effects of WM demyelination. It is well-established that the g-ratio is a sensitive marker for myelination and demyelination 21,28,29,31 . In particular, an increased g-ratio has been shown to indicate demyelinating lesions in pathological studies using electron microscopy 28 and MRI 21,32 . The significant negative correlation between the tract-averaged g-ratio and WM lesion volumes (i.e., number of demyelinating lesions) in the current study supports our hypothesis, with MS known to cause local demyelinating lesions 21,32 .
The NOS-based connectomes showed alterations in MS that were largely consistent with the results of previous connectomic analyses undertaken in MS 4   www.nature.com/scientificreports www.nature.com/scientificreports/ degree of myelination 18 . Therefore, the g-ratio-based connectome offers a more straightforward method for evaluating demyelinating diseases.
The g-ratio-based connectome detected significantly increased nodal strength in the MS group, whereas the NOS-based connectome detected no differences between the two groups, suggesting that the g-ratio-based nodal strength is more sensitive at detecting abnormal network characteristics in patients with MS. Remarkably, the increased g-ratio-based nodal strength was predominantly localized to the limbic system, a region where demyelinating lesions and brain atrophy tend to occur from the early stages of MS 35,41 . Abnormalities in the WM tracts connecting structures of the limbic system have also been found in early MS with diffusion tensor imaging 42 . It was also notable that the g-ratio-based nodal strength of several brain regions, such as the motor, visual, and limbic regions, were significantly correlated with EDSS scores, implying that these regions play a key role in the development of clinical symptoms in patients with MS. EDSS is a method of quantifying disability in MS by assessing impairment in motor, sensory, visual, and other functional systems 43 . Therefore, our result suggests that the g-ratio-based nodal strength facilitates MS diagnosis from early pathological stages and could be used to monitor disease progression and treatment effects in MS.
NOS-based nodal strength reportedly reduces in MS and NOS-weighted nodal strength is useful for diagnosing and monitoring disease progression in MS 13 . In the present study, NOS-weighted nodal strength in some areas tended to decrease in the MS group, but these were not significantly different than those in the control group. The conflicting findings between the present and previous studies may be related to discrepancies in patient characteristics and inherent limitations associated with the analysis of small sample sizes. For example, a previous study 13 included patients with MS with greater disease severity (as measured by EDSS scores) and total WM lesion volume (3.5 and 10.3 mL, respectively) compared than our study (0.9 and 4.3 mL, respectively). It is likely that the patients with greater disease severity and total WM lesion volume had greater disruptions in their brain networks. Another study 44 demonstrated a correlation between clinical symptoms and anatomical location of MS lesions, suggesting that the difference in WM lesion location also contributes to the differences between the findings of the present and previous studies. However, Shu et al. (2018) did not provide details on WM lesion locations 13 ; therefore, we could not examine this issue.
Certain limitations of the present study necessitate following considerations. First, the sample size of this preliminary study was small. Large-scale, multicenter studies are required to confirm the usefulness of g-ratio-based connectome analysis for evaluating MS pathology. Second, only EDSS was evaluated as an indicator of disease severity. Despite being the most popular and widely used outcome measure for disease progression in patients with MS, EDSS has been reported to have notable limitations 45 . For example, it has been pointed out that some functional domains are not sufficiently evaluated, such as cognitive function, mood, and energy levels 46 . Therefore, in future research, it will be necessary to introduce a more detailed clinical evaluation that also considers the relationship of the functional domains with the pathology of brain networks. The associations between EDSS and nodal strength ( Fig. 4) must be cautiously interpreted because of the limited dynamic range of the EDSS scores in this study and highly skewed distribution of scores (i.e., most individuals were associated with the lowest score). Third, graph-theoretic measures other than nodal strength, such as network efficiency and modularity, should be considered in future studies to better characterize the topological effects of MS. However, caution is required when interpreting measures based on path length in the g-ratio-based connectome because the g-ratio is not a direct estimate of conduction velocity or information transfer and may thus provide no advantage relative to NOS 30 . Ideally, both g-ratio and axon diameter are needed to measure conduction velocity 47 . Lastly, fluid-attenuated inversion recovery (FLAIR) images used for plaque detection were acquired at a slice thickness of 4 mm with a 1-mm gap because of the limitation of synthetic MRI sequence. A gap of at least 1 mm is needed to reduce the cross-talk between slices. However, the current imaging guidelines for MS recommend a slice thickness of 3 mm without a gap for 2D MR acquisition 48 . Furthermore, a previous study reported a significant (8%) increase in lesion volumes when the slice thickness of MR images was reduced from 5 to 3 mm 49 . Therefore, the lesion volume detected in our study may be underestimated. This problem should be addressed in future studies.
In conclusion, the g-ratio-based connectome can complement conventional measures of connectivity strength when assessing structural connectivity in patients with MS. Our data suggest that the g-ratio-based connectome has the potential to be used as a biomarker for disease diagnosis, particularly at early disease stages, and for monitoring disease progression.

Participants.
We prospectively included 14 females with relapsing-remitting MS (mean age 42.8 ± 5.0 years) from May to November 2016. All patients were diagnosed by an experienced neurologist (KY) according to the 2010 modified McDonald criteria 50 . Disease severity was evaluated using EDSS 43 . The main demographic and clinical characteristics, including EDSS score and disease duration of patients, were also assessed by KY. For comparison, we recruited 14 age-and sex-matched healthy controls (mean age 43.2 ± 14.4 years) with no history of neurological disease, psychiatric disease, drug abuse, head trauma, or seizures and no contraindications for MRI. The clinical and demographic characteristics are shown in Table 3.
The study was approved by the Institutional Review Board of Juntendo University Hospital, Japan, and was conducted in accordance with The Code of Ethics of the World Medical Association (Declaration of Helsinki) for experiments involving humans. A written informed consent was obtained for experimentation with human subjects.
Image acquisition. MRI data were obtained using a 3.0-T system (MAGNETOM Prisma; Siemens Healthcare, Erlangen, Germany) equipped with a 64-channel head coil. Multi-shell DW-MRI and MT saturation (MTsat) images were acquired to calculate AVF and MVF, respectively. (2019) 9:13522 | https://doi.org/10.1038/s41598-019-50025-2 www.nature.com/scientificreports www.nature.com/scientificreports/ DW-MRI was obtained at b-values of 1000 and 2000 s/mm 2 along 64 directions uniformly distributed on a sphere, with a simultaneous multi-slice echo-planar imaging sequence in the anteroposterior phase-encoding direction. The same diffusion directions were used for all shells. Each DW-MRI acquisition was complemented with a non-DW volume (b = 0 s/mm 2 ). Standard and reverse phase-encoded blipped non-DW-MRI was also obtained to correct echo-planar imaging distortions 51 . The acquisition parameters of DW-MRI were as follows: repetition time, 3300 ms; echo time, 70 ms; voxel size, 1.8 × 1.8 × 1.8 mm 3 ; number of slices, 65; simultaneous multi-slice factor, 2; number of excitations, 1; and acquisition time, 6.25 min. The multi-shell DW-MRI data were visually checked in all three orthogonal views to confirm that there were no severe artifacts due to missing signals, gross geometric distortion, or bulk motion. Then, the data were corrected for eddy currents, susceptibility-induced geometric distortions, and inter-volume motion using the EDDY and TOPUP toolboxes 52 .
The dual excitation three-dimensional (3D) multi-echo fast low-angle shot sequences were performed with predominant T1, proton density (PD), and MT weighting for calculating the MTsat index. MTsat imaging was developed to improve the MT ratio (MTR) by separating MTR from the longitudinal relaxation rate (R 1 ) 53 . MTsat shows higher contrast in the brain than MTR 53 and has been shown to correlate well with quantitative MT measures 54 . It has also been shown that MTsat is less sensitive to variations among MR imaging systems than the MTR parameter and was more effective in differentiating patients with MS from healthy controls 55 . The excitation of MT-weighted images was preceded by an off-resonance Gaussian-shaped RF pulse under the following conditions: frequency offset from water resonance, 1. To enable the estimation of cortical parcels and tissue segmentation (FreeSurfer), T1-weighted images (T1WI) were also acquired using 3D magnetization-prepared rapid gradient-echo sequence. The acquisition parameters were as follows: repetition time, 15 ms; echo time, 3.54 ms; inversion time, 1100 ms; voxel size, 0.86 × 0.86 × 0.86 mm 3 ; and acquisition time, 5.14 min.
MVF calculation. First, the MTsat data were analyzed to calculate MVF using an in-house MATLAB script based on a previously described theory 53 . First, for each voxel, the apparent longitudinal relaxation rate R app 1 was calculated as follows: where S T1 and S PD denote the signal intensities of T1-and PD-weighted images, respectively; TR T1 and TR PD denote the TR of T1-and PD-weighted images, respectively; and α T1 and α PD denote the excitation flip angles of T1-and PD-weighted images, respectively. Second, apparent signal amplitude A app was calculated as follows: T1  PD T1 PD  T1 PD T1   T1  PD T1  PD  T1 PD Third, apparent MTsat δ app was calculated as follows:  where S MT , TR MT , and α MT denote signal intensity, TR, and excitation flip angle of the MT-weighted image, respectively.
MTsat is superior to conventional MTR imaging in terms of relaxation rate, inhomogeneity of RF transmit, and receive field 53,59 . Further, to further improve spatial uniformity, we corrected small residual higher-order dependencies of the MTsat on the local RF transmit field, as suggested by Weiskpof et al. 60 , as follows: where RF local is the relative local flip angle α compared with the nominal flip angle. Dual-angle method was used to calculate RFl ocal 61 . Additionally, we added two B1 maps using echo-planar imaging with flip angles of 10° and 20° acquired in approximately 10 s each. The first image was obtained after excitation with a flip angle α and second after excitation with a flip angle 2α. The first and second images had a magnitude proportional to sin(α) and sin(2α), respectively. The ratio of the two acquisitions was subsequently formed as follows: (5) from which the local flip angle α was calculated.
The MTsat used in this study required a calibration factor to estimate MVF. We determined the calibration factor based on all healthy controls involved in this study as a linearly proportional relationship between MVF and δ app , an absolute quantitative marker of myelin 54 . The calibration factor of 0.1 was subsequently used to obtain a g-ratio of 0.7 in the corpus callosum as previously suggested 62,63 . Therefore, to estimate MVF maps, we multiplied all voxels in the δ app map by 0.1. 20 . In particular, the intracellular volume fraction (Vic) and isotropic water diffusion volume fraction (Viso), indicating free water, were calculated using AMICO 64 , enabling AVF estimation as follows 21 :

AVF calculation. AVF was estimated based on parameters of the NODDI model
AVF was independently estimated for each voxel. These calculations were performed using custom MATLAB functions (MathWorks, Natick, Massachusetts). The images were linearly registered to a common space using FLIRT 65 before calculation.
MRI g-ratio calculation. A g-ratio was calculated using MVF and AVF for each voxel, generating an in vivo whole-brain g-ratio map using the following equation 21 : Before obtaining g-ratio maps, MVF maps were first affine-aligned to the corresponding non-gradient-weighted MRI images using a statistical parametric mapping software (SPM12, Wellcome Department of Cognitive Neurology, UK; http://www.fil.ion.ucl.ac.uk/spm) This was independently repeated for each individual, and spatial correspondence between voxels was ensured on the AVF and MVF maps.

Measurements of WM lesion volume and distribution.
Given that the presence of demyelinating lesions can affect connectome mapping, we evaluated the distribution and volume of WM lesions. For all patients with MS, we automatically segmented WM lesions on synthetic FLAIR images using a lesion prediction algorithm 66 , as implemented in the Lesion Segmentation Toolbox, Version 2.0.15 (http://www.applied-statistics.de/ lst.html) 67 . This algorithm comprises a binary classifier in the form of a logistic regression model trained on the data of 53 patients with MS with severe lesion patterns. As covariates for this model, we used a similar lesion belief map with the lesion growth algorithm 67 and spatial covariate that considered voxel-specific changes in lesion probability. The parameters of this model fit were used to segment lesions by providing an estimate for the lesion probability per voxel. All lesion maps were visually inspected by a neuroradiologist (AH). The total brain WM lesion volume was calculated by multiplying the lesion area by slice thickness.
Synthetic T1-weighted images of each patient were spatially normalized to Montreal Neurological Institute space, and deformation fields were saved, which were then applied to lesion maps. Because synthetic images derived from QRAPMASTER are inherently aligned 56 , no registration was required between synthetic T1-weighted images and lesion maps created on synthetic FLAIR. Normalized lesion maps were summated for all patients to create aggregate lesion maps.
Pre-processing for connectome. Figure 5 shows a schematic overview of the processing pipeline for obtaining NOS-and g-ratio-weighted connectivity matrices. DW-MRI data were preprocessed using the Functional MRI of the Brain (FMRIB) Software Library version 5.0.9 68 . First, for each subject, we used the nonlinear boundary-based registration approach in FMRIB to align 3D-T1WIs to the corresponding b0 map. Second, the Brain Extraction Tool 69 was used to remove non-brain tissue from 3D-T1WIs. Third, the tissue probability maps of WM, cortical gray matter (GM), deep GM, and CSF were obtained using the FMRIB Automated Segmentation Tool 70  www.nature.com/scientificreports www.nature.com/scientificreports/ probabilities of deep GM for all voxels within the brain were obtained using the FMRIB Integrated Registration and Segmentation Tool 71 . These tissue probability maps were processed for the multi-shell, multi-tissue constrained spherical deconvolution (MSMT-CSD) and anatomically constrained tractography framework.
Node definition. The default FreeSurfer reconstruction pipeline 72 was used to delineate 84 regions based on the Desikan-Killiany cortical atlas segmentation (Supplementary Table S5) 73 . Because of the high variability in spatial location and extent of subcortical GM parcellation mapped with FreeSurfer 72 , the subcortical GM data provided by FreeSurfer parcellation were replaced by the subcortical GM partial volume maps obtained from the FMRIB Integrated Registration and Segmentation Tool.
Edge definition. For each subject, the MSMT-CSD probabilistic tracking method 74 with multi-shell DW-MRI data (b = 0, 1000, and 2000 s/mm 2 ) was used to generate whole-brain tractograms. Anatomically constrained tractography 75 and spherical deconvolution informed filtering of tractograms 76 were applied to reduce bias in the streamline density of longer fiber pathways 75 . While streamline filtering methods can alleviate biases, it is important to note that they can also compromise the specificity with which connectome pathology can be localized 19 . The MRtrix 3.0 software package (Brain Research Institute, Melbourne, Australia, http://www.brain. org.au/software/) was used to generate tractograms.
For MSMT-CSD tracking, multiple response functions were estimated as a function of the b-value and tissue type. First, voxels were assigned to WM if the fractional anisotropy was >0.7 and tissue probability was >0.95 for WM (see section 2.7). Second, the DW-MRI profile was reoriented to ensure that the principal axis of diffusion was aligned. Third, the WM anisotropic response functions were calculated per shell by averaging the reoriented DW-MRI profiles over these voxels. Fourth, if a tissue segmentation reported the respective tissue probability to be >95% with a fractional anisotropy of <0.2, these voxels were assigned to GM and CSF. The isotropic response functions for GM and CSF were acquired by averaging DW-MRI profiles per shell. Finally, the WM fiber orientation distribution function (fODF), GM fODF, and CSF fODF were obtained using the dwi2fod command with the Figure 5. Schematic overview of the pre-processing pipeline for obtaining NOS-weighted and g-ratio-weighted connectivity matrices. NOS-weighted and g-ratio-weighted network were mapped according to the following steps. (i) MVF maps were generated from MTsat data (T1-weighted, PD-weighted, and MT-weighted images); (ii) NODDI-derived metrics, Vic and Viso, were calculated from the multi-shell DW-MRI data; (iii) AVF maps were estimated from Vic, Viso, and MVF data; (iv) g-ratio maps were calculated from AVF and MVF maps in a voxel-by-voxel manner; (v) whole-brain probabilistic tractography was performed with MSMT-CSD using multi-shell DW-MRI data and 84 connectome nodes on the Desikan-Killiany cortical atlas; (vi) NOS between each pair of nodes was enumerated; and (vii) g-ratio values were averaged across all voxels intersecting the set of streamlines interconnecting a pair of regions to yield a myelin-specific connectivity weight to contrast with NOS. This was repeated for all pairs of regions. www.nature.com/scientificreports www.nature.com/scientificreports/ msmt-csd option implemented in the MRtrix software. The maximal spherical harmonic orders (lmax) of WM, GM, and CSF were set to 8, 0, and 0, respectively.
Whole-brain MSMT-CSD tracking was performed on the WM fODFs using the second-order Integration over Fiber Orientation Distributions algorithm 77 , with the following parameters: step size, 0.9 mm; angle thresholds, 45° per step; length, 3.6-180 mm; and fiber orientation distribution threshold, 0.05. In total, 5 × 10 7 streamlines were generated through seeding from the WM fODFs. Seeding points were determined dynamically using the SIFT model 76 . For tractogram reconstructions comprising 5 × 10 7 streamlines, SIFT was also applied to filter the reconstruction from 5 × 10 7 to 5 × 10 6 streamlines.

Connectome construction.
To map the NOS-based connectome for each individual, the total number of streamlines interconnecting each pair of regions was enumerated and stored in a connectivity matrix. The g-ratio-based connectome was determined by averaging g-ratio values across all voxels intersecting the set of streamlines interconnecting a pair of regions. While averaging was performed across voxels, it is also feasible to determine a streamline-averaged g-ratio value for each streamline and subsequently average across the set of streamlines 30 . This was repeated for all pairs of regions, and the resulting tract-averaged g-ratios were stored in a connectivity matrix. Therefore, two connectivity matrices of dimension 84 × 84 were mapped for each individual, with connectivity strength estimated using the NOS-or tract-averaged g-ratio. Diagonal elements were excluded from the connectivity matrices. A NOS density threshold (T) was applied to exclude spurious links 78 . The links above the T% according to NOS were left unaltered, whereas links below the T% were set to 0. To avoid bias from adopting a single threshold, NOS-weighted nodal strengths were examined across several thresholds (10% < T < 30% in 5% steps) 79 .

Nodal strength.
For the NOS-based connectome, NOS-weighted nodal strength was calculated as where i is a given node, A ij is the NOS between nodes i and j, and N is the number of nodes in the network. Nodal strength is the simplest measure of centrality of a given node and reflects the degree of interconnectivity with other regions 80 . To estimate effect sizes for differences in nodal strength between the MS and control groups, Cohen's d was computed for the mean NOS-weighted nodal strength across a range of thresholds (10% < T < 30% in 5% steps). Cohen's d was the largest for a threshold of 30% (Supplementary Table S3); therefore, this was applied as the threshold for the NOS-weighted nodal strength. Nodal strengths were computed using relevant functions in the MATLAB Brain Connectivity Toolbox (http://www.brain-connectivity-toolbox.net/).
For the g-ratio-based connectome, nodal strength was calculated as = , where i is a given node and G ij is the average g-ratio sampled along streamlines between the nodes i and j. Thus, nodal strength was computed as the weighted average of tract-averaged g-ratio values across all regions where the weighting factor was NOS. The NOS weighting factor accentuates the contribution of tract-averaged g-ratio values associated with connections comprising many streamlines. Connections comprising relatively few streamlines may be weak or spurious; thus, this weighting moderates the contribution of these connections.
Identification of disrupted WM connectivity. The NBS 81 was used to test for differences in connectivity strength between the MS and control groups. The NBS was separately performed for the NOS-and g-ratio-based connectomes to test the null hypothesis of equality in connectivity strength between the two groups for each of the 3,321 unique pairs of regions. Further details about the NBS are described elsewhere 81 . In brief, two-sample t-test was independently performed for each connection to test the null hypothesis of equality in connectivity strength between the MS and control groups. A primary component-forming threshold (P = 0.05, t = 2.06, two-tailed t-test) was then applied to define a set of suprathreshold edges (results across different thresholds are reported in Supplementary Table S4). Next, the number of links (or size) of the remaining connected components in the network was stored. The statistical significance of the size of each connected component was evaluated with respect to an empirical estimate of the null distribution of maximal component sizes (5000 permutations), with the component size measured as the number of edges it comprised. Any connected component was reported if it remained significant at a P-value of < 0.05 after family-wise error correction.
Statistical analysis. Statistical analyses of demographic and clinical variables were performed using IBM SPSS for Windows (version 23.0; IBM Corp., Armonk, NY, USA). Demographic and clinical variables, except the EDSS score, were confirmed to be normally distributed by Kolmogorov-Smirnov test. Between-group differences were analyzed using Student's t-tests for continuous variables (e.g., age, WM lesion load, and graph-theoretic metrics) and chi-square tests for sex. Because the EDSS score was not normally distributed, Spearman's rank correlation coefficient was used to test for relationships between brain measures (e.g., WM lesion load and nodal strength) exhibiting significant between-group differences and any clinical measures (e.g., disease duration and EDSS score). FDR was used to correct for multiple comparisons using a significance threshold of P < 0.05.

Data Availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.