Classification of multiple sclerosis based on patterns of CNS regional atrophy covariance

Abstract There is evidence that multiple sclerosis (MS) pathology leads to distinct patterns of volume loss over time (VLOT) in different central nervous system (CNS) structures. We aimed to use such patterns to identify patient subgroups. MS patients of all classical disease phenotypes underwent annual clinical, blood, and MRI examinations over 6 years. Spinal, striatal, pallidal, thalamic, cortical, white matter, and T2‐weighted lesion volumes as well as serum neurofilament light chain (sNfL) were quantified. CNS VLOT patterns were identified using principal component analysis and patients were classified using hierarchical cluster analysis. 225 MS patients were classified into four distinct Groups A, B, C, and D including 14, 59, 141, and 11 patients, respectively). These groups did not differ in baseline demographics, disease duration, disease phenotype distribution, and lesion‐load expansion. Interestingly, Group A showed pronounced spinothalamic VLOT, Group B marked pallidal VLOT, Group C small between‐structure VLOT differences, and Group D myelocortical volume increase and pronounced white matter VLOT. Neurologic deficits were more severe and progressed faster in Group A that also had higher mean sNfL levels than all other groups. Group B experienced more frequent relapses than Group C. In conclusion, there are distinct patterns of VLOT across the CNS in MS patients, which do not overlap with clinical MS subtypes and are independent of disease duration and lesion‐load but are partially associated to sNfL levels, relapse rates, and clinical worsening. Our findings support the need for a more biologic classification of MS subtypes including volumetric and body‐fluid markers.

disease phenotype distribution, and lesion-load expansion. Interestingly, Group A showed pronounced spinothalamic VLOT, Group B marked pallidal VLOT, Group C small between-structure VLOT differences, and Group D myelocortical volume increase and pronounced white matter VLOT. Neurologic deficits were more severe and progressed faster in Group A that also had higher mean sNfL levels than all other groups. Group B experienced more frequent relapses than Group C. In conclusion, there are distinct patterns of VLOT across the CNS in MS patients, which do not overlap with clinical MS subtypes and are independent of disease duration and lesion-load but are partially associated to sNfL levels, relapse rates, and clinical worsening. Our findings support the need for a more biologic classification of MS subtypes including volumetric and body-fluid markers.  (Kappos et al., 2020). Numerous previous studies have shown contradictory results on volume loss differences over time (VLOT) in various CNS structures between clinical phenotypes (De Stefano et al., 2010;Tsagkas et al., 2018b;Tsagkas et al., 2020). Furthermore, a recent post-mortem study identified a separate patient group, named myelocortical MS, with demyelination restricted to the cortical gray matter and the spinal cord, which included patients from all clinical phenotypes (Trapp et al., 2018). Thus, a classification that better reflects underlying pathology is still needed and might be crucial for the effective selection of disease-modifying treatments .
In healthy adults, brain volume changes over time vary both across individuals and CNS regions (Raz et al., 2005). In previous studies, longitudinal atrophy measures have displayed a large variance between MS patients (De Stefano et al., 2010;Tsagkas et al., 2018b). Apart from these between-patient heterogeneities, several studies have shown dissociation in the progression of atrophy among different CNS compartments, suggesting an additional spatiotemporal within-subject heterogeneity Tsagkas et al., 2018b). This indicates that volume loss, progresses independently across CNS regions. Hence, we hypothesized that the anatomic patterns of atrophy progression over time could distinguish patient groups with different underlying pathomechanisms.
In this study, we aimed to reveal the presence of distinct patterns of VLOT between different CNS-structures and use these to identify patient subgroups in a large MS-cohort. We then explored the clinical, radiological, and serological features of the identified patient groups.

| Study design and participants
Clinical and MRI data of an ongoing large-scale cohort of MS patients (260 patients in total) from a single center (Multiple Sclerosis Center, University Hospital, Basel, Switzerland) were analyzed, retrospectively. All of the patients participating in our study have been previously reported Tsagkas et al., 2018a;Tsagkas et al., 2018b;Tsagkas et al., 2020). Patients were followed for up to 6 years. The MS diagnosis was made in accordance with international panel established criteria (McDonald et al., 2001). The local ethics committee approved the study (EKBB-46/04) and all patients provided written informed consent.

| Procedures
All patients underwent yearly standardized neurological examinations including the Expanded Disability Status Scale (EDSS; www. neurostatus.org) by trained and certified examiners, the dominant hand and non-dominant hand 9-hole peg test (D9HPT and ND9HPT) as well as the timed 25-ft walk test (T25fwt). The occurrence of new relapses was recorded at every visit. Serum samples were collected on the same day as the clinical visit and serum neurofilament light chain (sNfL) levels were measured by Simoa assay as previously described (Barro et al., 2018).

| MRI analysis
2.3.1 | Lesion segmentation All brain white matter lesions were segmented on T2-weighted and proton density images by trained expert observers according to the standard operating procedures used at the Medical Image Analysis Center for the analysis of clinical phase II and phase III trials .
T2-weighted lesion volume (T2LV) was calculated for the whole brain.

| Volumetry of CNS structures
All morphological analyses were performed on the T1-weighted MPRAGE brain images.
Brain white matter was computed with the fully automated tool SIENAX (version 2.6) (Smith et al., 2004). Cortical gray matter was fully automatically segmented using CIVET (version 2.1.0) (Zijdenbos, Forghani, & Evans, 2002), as also used previously in the context of several longitudinal studies (Bhagwat et al., 2018;Makowski et al., 2016;Redolfi et al., 2015;Tsagkas et al., 2020). Before segmentation with CIVET, the T1-weighted images were linearly registered to the standard stereotaxic space defined by the MNI ICBM 152 model (Mazziotta et al., 2001). The images were then corrected for intensity nonuniformity using N3 (Sled et al., 1998) and a linear registration to the model (Collins et al., 1994) was applied.
The volume of the deep gray matter nuclei including thalamus, striatum, and globus pallidus was estimated based on an established nomenclature using MAGeT as previously described (Chakravarty et al., 2013) and deployed in previous cross-sectional studies (Magon et al., 2015). Data pre-processing and detailed MAGeT application was described in a previous study .
In short, both the medical imaging netCDF (MINC) toolkit (version 2) and the advance normalization tools (ANTs) were deployed as follows: (1) bias field correction using N4-correction algorithm (Tustison et al., 2010); (2) non-local means denoising (Manjón et al., 2010); (3) affine registration using a normalized mutual information objective function (Studholme et al., 2001); and (4) brain extraction using the BEaST algorithm (Eskildsen et al., 2012). The preprocessed data were then used as input for the MAGeT algorithm. The MAGeT registration procedure is also similar to CIVET, which should at least reduce potential biases between these two algorithms that calculated the cortical, thalamic, pallidal, and striatal volumes used in this study. Quality assessment was performed on each T1-weighted image by a trained researcher to ensure segmentation correctness.
MAGeT, CIVET, and SIENAX were performed on T1-weighted images after lesion filling using the approach previously proposed by our group  in order to reduce biases related to tissue misclassification and improve the registration step (Sdika & Pelletier, 2009). The SIENAX baseline volume-correction factor was used for normalizing cortical gray matter, white matter, and deep gray matter nuclei regarding variations of head size. All analyses were performed on these corrected volumes.
Spinal cord volume analysis was performed using the established CORDIAL semi-automatic software, as described in previous methodological and clinical studies (Amann et al., 2016;Tsagkas et al., 2018b). Shortly, the segmentation was carried out over a 35 mm long spinal cord segment, starting 27 mm below the cisterna pontis, which corresponds roughly to the spinal cord volume between the Foramen magnum and the C2/C3 intervertebral disk. Segmentations were visually inspected for quality and excluded from further statistical analysis in case of segmentation errors.

| Clustering of MS patients
The mean annual loss rate of cortical gray matter, brain white matter, thalamus, striatum, globus pallidus, and the spinal cord over a maximum follow-up time of up to 6 years was determined as the average percentage of the annualized changes between all available timepoints for every patient. In order to correct for the effect of age on the mean annual atrophy rates of those regions, we calculated the residuals from a regression analysis between the atrophy rate of each region and age. Subsequently, a principal component analysis was performed on these residuals of the mean annual loss rates using the singular value decomposition method (Husson, Josse, & Pages, 2010).
Principal components represent different parts of the betweenstructure covariance and, thus, reflect between-structure patterns of CNS VLOT. Each principal component extracted using this approach highlighted the difference or similarity between volume changes across individual CNS regions.
Principal components kept in further statistical analysis had to meet the following criteria: (1) >90% of total variance retained, (2) eigenvalue >0.7, (3) "broken stick rule," (4) scree plot "elbow rule" (Jolliffe, 2002). Hierarchical clustering was then performed using Ward's-criterion on the selected principal components (Husson et al., 2010). Extraction of the resulting groups was done after visualizing the cluster dendrogram ( Figure 1). In order to validate the stability of the resulting clusters, we repeated the principal component analysis and subsequently the hierarchical clustering analysis 10,000-fold, using the same methods described above. In every repetition, a contribution weight of 0.0001 to the principal component analysis was ascribed to a random 10% of our patients and a contribution weight of 2 to another random 10%, leaving the rest 80% with a uniform contribution weight of 1. This allowed "suppression" of a smaller subject number-thus, avoiding bias against smaller clustersbut inserted significant data modifications for reliable cross-validation of our classification. After each repetition, the adjusted Rand index (Rand, 1971) and Meila's variation of information distance (Meil a, 2007) were calculated to measure the similarity between the new clustering outcome and our initial classification.

| Group comparisons
In order to approximate a normal distribution, logarithmic transformations were performed for the EDSS, D9HPT and ND9HPT, whereas an inverse transformation was performed for the T25fwt. Baseline MS severity score (Roxburgh et al., 2005), mean annual relapse number and the mean sNfL concentration were calculated. Comparisons of baseline demographic factors, clinical measurements, and number of follow-ups between groups were made using Welch's two sample t-tests and Pearson's chi-squared test with Yates's continuity correction. Between-group differences in regard to baseline MRI measures and annual rates of VLOT were performed using analyses of covariance, while correcting for age, sex, and disease duration.
Statistical analysis of the patients' MRI metrics and clinical outcomes (EDSS, D9HPT, ND9HPT, T25fwt) over 6 years was performed using separate linear mixed effect models (LMER) in order to explore between-group differences. For both model types, this was done using a random intercept and a random time-slope for each subject to allow for within-subject and between-subject variance. Independent variables were entered blockwise keeping the following sequence: first demographics, then disease duration and finally the aforementioned classification result. Each variable was tested both for its correlation to the intercept as well as to the slope over time of the dependent variable. All independent variables without statistical significance were excluded from the final model. Within-group comparisons between CNS structures were performed using paired sample ttests, whereas respective effect sizes were computed as Cohen's d for paired samples.
We performed a Cox regression analysis in a backwards stepwise fashion including sex, age, disease duration, and the classification outcome of our analysis to evaluate differences between the identified MS patient groups in terms of disease progression and time-todisease progression. Clinical disease progression was defined according to the following conventions: (1a) an increase of 1 point in the EDSS if the baseline EDSS score was ≤5.5 or (1b) an increase of 0.5 point in the EDSS if the baseline EDSS score was >5.5, and (2) no relapse in the last 12 months.

| Results of clustering analysis
After excluding MRI sessions with segmentation failures or errors (162 of 1,271 MRI) in either of the quantified CNS regions and patients with availability of only 1 MRI dataset (in total, 29 patients excluded due to missing values) as well as patients diagnosed with clinically isolated syndrome that did not convert to clinically definite MS during follow-up (5 patients in total), a total of 226 patients, and 1,080 examination datasets (including clinical, sNfL, and MRI data) were entered in the statistical analysis.
In a first step, we performed a regression analysis between the annual atrophy rates of the spinal cord, striatum, globus pallidus, thalamus, cortical gray matter, and brain white matter on one side and age on the other. Subsequently, we performed a principal component analysis using the residuals from each regression analysis. The first four principal components accounted cumulatively for 92% of the total data-variance and were kept in further statistical analysis. Table S1 displays Pearson's correlation coefficients between the residuals of the mean annual atrophy rates of the analyzed CNS regions after correcting for age. Detailed results of those four principal components are displayed in Table 1.
The dendrogram of the hierarchical clustering is shown in  Table 2. Groups did not differ in age, sex, disease duration, and clinical phenotype distribution. Patients of Group A had a higher number of untreated patients compared to Group C, whereas the distribution of treatments used in Group C was different compared to Group D.
Cross-validation of the clustering process with 10,000-fold repetition showed a clustering agreement with a mean adjusted Rand

| Within-group patterns of VLOT
Mean annual atrophy rates of the spinal cord, striatum, globus pallidus, thalamus, cortical gray matter, and brain white matter volumes as well as within-group differences between the atrophy rates of all CNS-structures and the respective size effects are reported in Table 3

| Longitudinal analysis of MRI measurements between groups
Detailed results for each analyzed CNS structure and T2LV are displayed in Table 4. Summarizing the most important results, Group A demonstrated significantly greater spinal cord and thalamic VLOT compared to all other groups. Group B showed significantly greater pallidal VLOT compared to groups C and D. Group C showed lower VLOT in all CNS structures compared to Groups A and B. Group D showed a significant cortical gray matter volume increase compared to all other groups and a significant spinal cord volume increase compared to Groups A and B; additionally, a significantly greater white matter VLOT than in Group C was shown. T2LV increase over time did not differ between groups.

| Longitudinal analysis of clinical outcomes between groups
Trends of raw values of EDSS, T25fwt, and 9HPT in Groups A, B, C, and D over 6 years are displayed in Figure 3.

| EDSS
Between-group differences of the average logEDSS were assessed using LMER analyses after correcting for sex, age, and disease duration. Group A had worse average logEDSS compared to Groups B (difference: 0.36 ± 0.10; p < .01) and C (difference: 0.48 ± 0.10; p < .001), but a similar average logEDSS compared to Group D. The average logEDSS also did not differ between Groups B, C, and D. LogEDSS change over time did not differ between all groups.

| T25fwt
Between-group differences of the average T25fwt −1 were assessed using LMER analyses after correcting for age and disease duration.

| 9HPT
Between-group differences of the average logD9HPT were assessed using LMER analyses after correcting for sex and age. Group A had worse average logD9HPT compared to Groups B (difference: 0.36 ± 0.9; p < .001) and C (difference: 0.46 ± 0.9; p < .001), but a similar average logD9HPT compared to Group D. The average logD9HPT did not differ between Groups B, C, and D. Moreover, the logD9HPT worsening over time was faster for Group A compared to Group B (difference: 0.06 ± 0.02/year; p < .05) and Group C (difference: 0.07 ± 0.02/year; p < .05), but similar compared to Group C. Changes of logD9HPT over time did not differ between Groups B, C, and D.
Between-group differences of the average logND9HPT were assessed using LMER analyses after correcting for sex, age, and disease duration. Group A had worse average logND9HPT compared to all other groups (differences vs.: Group B: 0.54 ± 0.10, p < .001, Group C: 0.65 ± 0.10, p < .001; Group D: 0.52 ± 0.14, p < .01). The average logND9HPT did not differ between Groups B, C, and D. Moreover, the logND9HPT worsening over time was faster for Group A compared to all other groups (differences vs.: Group B: 0.10 ± 0.02/year, p < .001, Group C: 0.12 ± 0.02/year, p < .001; Group D: 0.11 ± 0.03/year, p < .01). LogND9HPT worsening over time was similar between GROUPs B, C, and D.

| Disease progression
In our Cox analysis no differences between Groups A, B, C, and D were found in terms of disease progression and time-to-disease progression. T A B L E 3 Annual volume changes as well as within-group differences and effect size of within-group differences of MRI measures in patients with multiple sclerosis to higher mean sNfL compared to Group C (p = .066). All other between-group differences were not statistically significant. Raw mean serum sNfL-values between groups are shown in Figure 4.

| DISCUSSION
In this work, we aimed to reveal distinct anatomic patterns of VLOT in selected CNS structures and then use these patterns to identify patient subgroups in MS. For this purpose, spinal cord, striatal, pallidal, thalamic, cortical gray matter, and white matter volumes were measured in a large MS cohort followed-up annually over 6 years. Patterns of VLOT were extracted using a principal component analysis after removing the effect of age from the VLOT for each individual CNS structure. These patterns were then entered in a hierarchical clustering analysis for unsupervised identification of patient groups.
With this method, we identified four patient groups with distinct anatomic patterns of CNS VLOT. Moreover, we showed that these groups were associated with distinct clinical, radiological, and blood biomarker characteristics.
Our statistical approach allowed to extract four distinct patterns of VLOT. Group A was characterized by a more prominent spinothalamic VLOT compared to all other measured structures. In contrast, Group B displayed a marked pallidal VLOT. Group C was generally characterized by low volume changes and generally low betweenstructure VLOT differences. Finally, Group D showed an increase in corticospinal volume as well as a pronounced white VLOT compared to other CNS structures. Patients in these groups did not differ by sex, age, disease duration, and more importantly clinical phenotype.
All three clinical phenotypes (relapsing-remitting, secondary progressive, and primary progressive MS) were represented in each of these groups. The 10,000-fold cross-validation analysis showed good agreement between the resulting clusters, supporting the robustness of our F I G U R E 2 (a) Boxplots of mean annual volume changes of the spinal cord, striatum, globus pallidus, thalamus, cortical gray matter (GM), and brain white matter (WM) by group. Whiskers correspond to annual volume changes' 25th and 75th percentiles. (b) Effect size of within-group differences in volume loss over time between CNS-structures by group. Whiskers correspond to the effect size 95% confidence intervals. Absolute Cohen's d values are indicated with a light blue to red gradient. Shape indicates negligible, small, medium, and large effect sizes. Abbreviations: SC, spinal cord; STR, striatum; GP, globus pallidus; THLM, thalamus; cGM, cortical gray matter; WM, brain white matter results in this MS-population. In summary, our findings clearly add to the evidence currently challenging applied clinical definitions (Kappos et al., 2020;Trapp et al., 2018) and their predictive value.
One patient of our cohort was not classified in either of the four identified subgroups and was excluded from further statistical analysis. Interestingly, this patient was the only pediatric-onset MS patient of our study. Unfortunately, due to the lack of a larger pediatric-onset MS group in our cohort, we could not verify if patterns of CNS VLOT differ in this patient group compared to adult-onset MS patients. Nevertheless, this could potentially be the focus of future investigations.
Interestingly, the four identified groups-resulting from distinct patterns of VLOT across the CNS-displayed differences with regard T A B L E 4 Average CNS volumes and volume changes over time as well as between-group differences in patients with multiple sclerosis as calculated by linear mixed effect regression models over a maximum of 6 years Note: All depicted values have been calculated using linear mixed effect regression models after correcting for sex, age, and disease duration. Comparisons between groups with regard to average volumes changes and volume changes over time were performed with post hoc analysis using the Tukey test. Abbreviations: GM, gray matter; GP, globus pallidus; n, number of subjects; n.s., not significant for any pairwise comparisons between Groups A, B, and C or between CNS-structures within-groups; SC, spinal cord; SE, standard error; WM, white matter.
to their respective clinical courses during the monitoring-period of our study. Worsening of motor outcomes (walking speed and hand dexterity) was generally more severe and faster in Group A (characterized by pronounced spinothalamic atrophy) compared to all other groups.
This finding is in line with recent studies pointing to strong correlations of spinal cord and thalamic volume loss to clinical disease progression Magon et al., 2020;Tsagkas et al., 2018b). In accordance to these findings, Group A, compared to all other groups, also had a higher mean sNfL which is a specific biomarker of neuroaxonal damage and loss that also correlates with clinical progression (Barro et al., 2018). EDSS worsening over time (as a continuous variable in our linear mixed effect models) and definite EDSS progression (as binary outcome in our survival analysis) was similar between groups, which may be ascribed to the limitations of the EDSS (e.g., non-linearity, low sensitivity for clinical changes especially in advanced disease stages) combined with the small number of patients included in some of the identified patient groups. Baseline MS severity score, however, was also higher for Group A, indicating a more aggressive disease course in this group even before the monitoring-period of our study.
It is worth emphasizing that although average lesion-load was statistically higher for Group A compared to Groups B and C, T2LV which surprisingly did not correspond to demyelinated areas in histological analysis. Moreover, volumetric analysis in the MR-images of those patients revealed a significantly higher cortical gray matter parenchymal fraction and cortical thickness compared to "typical MS." Interestingly, our Group D showed a paradoxical expansion of spinal cord and cortical gray matter and faster brain white matter atrophy progression compared to other CNS structures (Figure 2). This was unique in this group compared to the other three groups in our cohort. Although increase of CNS volume is theoretically not to be expected in MS patients (especially when assessed over large periods of time, as in our study), a recent study examining large numbers of patients in a longitudinal fashion has also identified a small proportion of patients with brain volume increase (Andelova et al., 2018). Despite the fact that our study cannot provide sufficient evidence to this end, it could be speculated, that Group D corresponds to the described myelocortical MS.
The majority of current therapeutic strategies in MS have an immunomodulatory effect, which has been measured-for the greater part-in terms of a reduction of relapses and new or enlarging lesions (Tintore, Vidal-Jordana, & Sastre-Garriga, 2019). Our work displays a certain deviation from this concept with groups generally having similar characteristics in the "no evident disease activity-3" framework. In addition, a number of relapsing-remitting patients demonstrating a fast progressing CNS atrophy, and a number of secondary progressive and primary progressive MS patients sharing clinical and radiological features with relapsing-remitting patients. Since evaluation of inflammatory activity guides clinical decisions in current daily practice (Sloane, Mainero, & Kinkel, 2015), it might be important to evaluate the effect and efficacy of existing and future immunomodulatory and neuroprotective agents on these atrophy-pattern-based patient groups.
Our work demonstrated specific non-random patterns of volumetric CNS changes in MS patients. Different physiologic and pathologic mechanisms could have contributed to this. Firstly, the extent of cortical demyelination and presence of organized meningeal B-cell infiltrates, possibly contributing to subpial and cortical pathology, has been shown to vary among MS patients (Albert, Antel, Brück, & Stadelmann, 2007;Howell et al., 2011;Kooi, Geurts, van Horssen, Bø, & van der Valk, 2009). In addition, the blood-brain barrier displays clear structural and functional heterogeneity in different brain regions immune-effector mechanism in each person (Jarius et al., 2017;Metz et al., 2014). Therefore, differences in our groups may-at least partly-reflect distinct pathophysiologic mechanism. Future research may be able to classify patients using single-time point biomarkers without the need for long-term longitudinal metrics.
Some considerations for future studies should be discussed. In this proof-of-concept study, we consciously avoided to enter brain lesion-load into our principal component analysis and subsequently into our hierarchical clustering method in order to classify patients solely based on CNS VLOT patterns. This was also done in order to evaluate, whether the identified groups would differ in terms of lesion-load changes over time. Future investigations should evaluate any potential contribution of brain lesion-load or -number metrics and their localization to a potential MRI-based classification of MS patients. In addition, our dataset did not contain spinal cord MRimages (the upper cervical spinal cord volume was extracted from brain MPRAGE images), so that this metric could not be included in our work. Symptomatic and asymptomatic spinal cord lesions currently have a central role in the correct diagnosis (Brownlee, Swanton, Miszkiel, Miller, & Ciccarelli, 2016;Filippi et al., 2016;Geraldes et al., 2018;Thompson et al., 2018;Tintore et al., 2016), prognosis evaluation (Arrambide et al., 2018;Brownlee et al., 2017;Kantarci et al., 2016;Sombekke et al., 2013) and disease monitoring of MS (Zecca et al., 2016 Furthermore, with regard to the preprocessing steps of our MRI data prior to the segmentation of brain structures, the MAGeT registration procedure is similar to CIVET, which should at least reduce potential biases between these two algorithms that calculated the cortical, thalamic, pallidal, and striatal volumes used in this study. However, this step differs in SIENAX used for white matter volumes, which could have introduced bias in our volumetric assessment. Finally, despite the fact that our clustering analysis uniquely classified each patient in one group, it cannot be excluded that the identified groups represent different "states," from or to which individual patients shift during the course of the disease.
In summary, in this study, we identified patient groups with distinct patterns of regional volume change over time across the CNS in a large MS-cohort. This study suggests that these patterns of longitudinal CNS volume change at least partially correspond to different patterns of clinical progression and sNfL-levels. Our work also supports the value of MRI-assessed volume change patterns in the subclassification of MS patients and provides proof of concept for future studies.

ACKNOWLEDGMENTS
We are very grateful to all participants and staff involved in the GeneMSA cohort study, in particular Alain Thoeni. We would also like to warmly thank Prof. PhD François Husson for his help and statistical guidance in our work.

ETHICS APPROVAL STATEMENT
The local ethics committee approved the study (EKBB-46/04).

PATIENT CONSENT STATEMENT
All patients provided written informed consent.

DATA AVAILABILITY STATEMENT
Anonymized raw data are available upon reasonable request.