Amélie Martel

Guillaume Coiffier
email: guillaume.coiffier@chu-rennes.fr.

Aurore Bleuzen

Jean Goasguen

Michel De Bandt

Christophe Deligny

Julie Magnant

Nicole Ferreira

Elisabeth Diot

Aleth Perdriger

François Maillot

What is the best salivary gland ultrasonography scoring methods for the diagnosis of primary or secondary Sjögren's syndromes?

Keywords: Sjogren's syndrome, Salivary gland, ultrasound

Objective: To evaluate the performance of salivary gland ultrasonography for the diagnosis of primary and secondary Sjögren's syndromes (pSS and sSS). Method: Multicenter cross-sectional study on 97 patients with clinical sicca symptoms. The pSS (n = 39) and the sSS (n = 22) met the American-European Consensus Group (AECG) classification criteria. The control patients (n = 36) with sicca symptoms did not fulfill the AECG criteria. Four scores were used to evaluate the 4 major salivary gland echostructure: the Salaffi score (0-16), Jousse-Joulin score (0-4), Hocevar score (0-48) and Milic score (0-12). Results: The medians of ultrasonographic (US) scores were higher in the pSS and sSS groups than in the control group (p <0.001). The receiver-operating characteristic (ROC) curves and the positive likelihood ratio (LR+) of the four scores showed a good diagnostic performance for the US diagnosis of pSS and sSS. Respectively, for pSS and sSS, the AUC were 0.891 (95%CI 0.812-0.970) and 0.824 (95%CI 0.695-0.954) for Hocevar score, 0.885 (95%CI 0.804-0.965) and 0.808 (95%CI 0.673-0.943) for Milic score, 0.915 (95%CI 0.848-0.982) and 0.844 (95%CI 0.724-0.965) for Salaffi score, 0.897 (95%CI 0.821-0.973) and 0.851 (95%CI 0.735-0.968) for Jousse-Joulin score. This study showed an interesting inter-observer reproducibility (kappa=0.714 ± 0.131) of the US evaluation with 85,7 % agreement between reader to determine the pathological character of the salivary glands. Conclusion: Salivary gland US is a simple, non-invasive and performant imaging procedure for the diagnosis of pSS and sSS, with Salaffi, Milic and Jousse-Joulin scores.

Introduction

The primary Sjögren's Syndrome (pSS) is a chronic autoimmune disease. It is characterized by a lymphocytic infiltration and destruction of the salivary and lacrimal glands, which causes dryness of the eyes (xerophthalmia) and mouth (xerostomia) by loss of glandular functions. Asthenia, inflammatory arthralgia often complete this clinical background. Almost a third of patients have systemic involvement with pulmonary, renal, neurologic,

A C C E P T E D M

A N U S C R I P T rheumatologic, vascular or lymphatic disorders. The pSS is the systemic disease for which the incidence of non-Hodgkin's Lymphoma (NHL) is the highest [START_REF] Theander | Lymphoma and other malignancies in primary Sjögren's syndrome: a cohort study on cancer incidence and lymphoma predictors[END_REF][START_REF] Zintzaras | The risk of lymphoma development in autoimmune diseases: a meta-analysis[END_REF]. Biologically, anti-SSA/SSB antibodies are the most specific markers, but are present in only 50% to 75% of pSS diagnostics. The anti-nuclear antibodies present in 80% of pSS and the rheumatoid factors in 40% to 75% of pSS, are non specific and can be found in other systemic autoimmune diseases [START_REF] Cornec | Diagnostic accuracy of blood B-cell subset profiling and autoimmunity markers in Sjögren's syndrome[END_REF][START_REF] Fauchais | Immunological profile in primary Sjögren syndrome: clinical significance, prognosis and long-term evolution to other auto-immune disease[END_REF][START_REF] Bournia | Subgroups of Sjögren syndrome patients according to serological profiles[END_REF]. The Sjögren's Syndrome is considered secondary (sSS) when associated with other connective tissue disease.

Like other systemic diseases, there is no single and specific diagnostic test. The diagnosis of the Sjögren's syndrome is based on a range of clinical, biological, morphological and histological symptoms. The most widely used classification criteria were proposed by the American-European Consensus Group (AECG) in 2002 [START_REF] Vitali | Classification criteria for Sjögren's syndrome: a revised version of the European criteria proposed by the American-European Consensus Group[END_REF]. They are used to define and standardize the diagnosis of primary and secondary SS. Among these criteria, parotid sialography and salivary gland scintigraphy are morphological examinations that are now rarely used in clinical practice, because of the radiating character of sialography and the low specificity of scintigraphy [START_REF] Salaffi | Ultrasonography of salivary glands in primary Sjogren's syndrome: a comparison with contrast sialography and scintigraphy[END_REF]. The ultrasonography (US) appears to be an attractive imaging approach by its non-invasive, inexpensive, nonirradiating and accessible technique. It allows also morphologic evaluation. Major salivary gland US (2 parotid and 2 submandibular glands) has been developed over the past 20 years. Several studies have shown that it provides similar information to sialography while being considerably less invasive [START_REF] Salaffi | Ultrasonography of salivary glands in primary Sjogren's syndrome: a comparison with contrast sialography and scintigraphy[END_REF][START_REF] Milic | Ultrasonography of major salivary glands could be an alternative tool to sialoscintigraphy in the American-European classification criteria for primary Sjogren's syndrome[END_REF]. MRI showed a high diagnostic performance [START_REF] El Miedany | Quantitative ultrasonography and magnetic resonance imaging of the parotid gland: can they replace the histopathologic studies in patients with Sjogren's syndrome?[END_REF][START_REF] Makula | The place of magnetic resonance and ultrasonographic examinations of the parotid gland in the diagnosis and follow-up of primary Sjogren's syndrome[END_REF][START_REF] Niemelä | Ultrasonography of salivary glands in primary Sjogren's syndrome. A comparison with magnetic resonance imaging and magnetic resonance sialography of parotid glands[END_REF] with good correlation between parotid MRI and salivary gland US (kappa agreement = 0.87) [START_REF] El Miedany | Quantitative ultrasonography and magnetic resonance imaging of the parotid gland: can they replace the histopathologic studies in patients with Sjogren's syndrome?[END_REF]. However, its lack of availability and its high cost does not allow its use for diagnostic purpose, but rather for the exploration of a suspicious glandular lesion.

With US B mode, the most specific pathogenic characteristic in SS is the inhomogeneity of the glandular parenchyma. Since the first proposal of US score made in 1992 by De Vita et al. [START_REF] Vita | Salivary gland echography in primary and secondary Sjögren's syndrome[END_REF], several US scores have been proposed with different items: gland size, contour regularity, clearness of the posterior glandular border, presence of hypoechogenic zones, echogenic bands and size of the cysts. The lack of consensual score for SS diagnosis is therefore one of the first limits for the use of salivary gland US in clinical practice. In another hand, few studies have evaluated this imaging method for sSS diagnosis. Finally, the main disadvantage of US is its operator-dependent character, with currently little data on its inter-and intra-observer reproducibility.

Our objective was to compare the performance of four scoring methods for the diagnosis of pSS and sSS, and to evaluate inter and intra-observer reproducibility of salivary gland US.

Methods

This is a multicentre, cross-sectional observational study carried out in three French university hospitals (Rennes, Tours and Fort-de-France).

Patients

All included patients complained of ocular and/or oral dry syndrome. Some already had a diagnosis of pSS or sSS established according to the AECG classification criteria. Others were explored for a suspicion of SS. Patients with sicca symptoms who did not fulfill the AECG criteria were included in the control group. The inclusions have been taking place between January 2014 and November 2016, in the departments of rheumatology in Rennes, in internal medicine in Tours and Fort-de-France. All patients were informed of the objectives and procedure of the study and gave their consent. The study was approved by the ethics committee of Rennes university hospital.

Exclusion criteria were: surgery for the removal of major salivary glands (≥ 2 glands / 4), graft versus host disease, history of sarcoidosis, cervical-facial radiotherapy, hepatitis C, HTLV-1 or HIV infection. The use of drugs causing sicca symptoms was not considered as an exclusion factor because it does not interfere with the results of salivary gland US.

A C C E P T E D

M A N U S C R I P T 2.

Clinical, histopathological évaluations and laboratory investigations

A standardized questionnaire was used to evaluate the severity of the subjective sicca syndrome, rated from 1 to 6 [Appendix A, Supplementary data S1; See the supplementary material associated with this article online]. Age, sex, duration of development of the sicca symptoms (<5, between 5 and 10,> 10 years), medications providing a sicca syndrome and active smoking were notified. Schirmer's test (abnormal if ≤ 5mm/5mn), salivary flow (abnormal if ≤ 0.5g/5min) and minor labial salivary gland biopsy (abnormal if grade 3 or 4 according to Chisholm and Mason or focus score ≥1 [START_REF] Chisholm | Labial salivary gland biopsy in Sjögren's disease[END_REF]) were performed unless the patient has already had these examinations. In this case, the results were taken from the medical records.

Anti-SSA/SSB antibodies, anti-nuclear antibodies (ANA) and rheumatoid factor (RF) levels, blood protein electrophoresis assays were investigated in the medical records or performed during sicca syndrome exploration.

Ultrasonography of the major salivary glands

US examination of the parotid and submandibular salivary glands was performed on every included patient to evaluate the four scoring US methods, with a different US scanner in each center (Philips, General Electrics and Esaote). Images of the four salivary glands were obtained in transverse and longitudinal planes with a 5-18 MHz linear probe. Each of the three examiners was blinded of clinical symptoms and diagnosis: a rheumatologist in Rennes (GC), a radiologist in Tours (AB) and Fort-de-France (JG) with respectively five, fifteen and twenty years of practice of salivary gland US. Four US scores were chosen in this study for their diagnostic performance described in previous studies: Salaffi (0-16) [START_REF] Salaffi | Salivary gland ultrasonography in the evaluation of primary Sjögren's syndrome. Comparison with minor salivary gland biopsy[END_REF], Jousse-Joulin (0-4) [START_REF] Cornec | Contribution of salivary gland ultrasonography to the diagnosis of Sjogren's syndrome: toward new diagnostic criteria?[END_REF], Hocevar (0-48) [START_REF] Hocevar | Ultrasonographic changes of major salivary glands in primary Sjogren's syndrome. Diagnostic value of a novel scoring system[END_REF] and Milic (0-12) [START_REF] Milic | Major salivary gland sonography in Sjogren's syndrome: diagnostic value of a novel ultrasonography score (0-12) for parenchymal inhomogeneity[END_REF] (Figure 1 and Appendix 1, Supplementary Data S1). The four scores propose different descriptions to evaluate the glandular parenchyma with several items (hypo-echogenic areas, echogenic bands, cysts, clearness of glandular borders). The Milic score is a simplified method with only one of these items: the parenchymal inhomogeneity. The Salaffi score evaluates all these items but is less accurate than Hocevar score. Jousse-Joulin score assesses the size of cysts and echogenic bands. All the included patients benefited from a centralized reading of the US images to establish the 4 scores (GC).

The Rennes (GC) and Tours (AB) examiners performed a second reading of their US images, few months after the first reading, to evaluate the intra-observer reproducibility of the Salaffi score. The Rennes examiner (GC) read again all the images from Tours to determine the inter-observer reproducibility of the Milic and Salaffi scores. In the event of a discrepancy between Tours and Rennes results, the US score evaluated by the examiner in Rennes (GC) was used.

Statistical Analysis

The statistical analysis software used was SPSS-statistics 19.0. The sensitivity (Se), the specificity (Sp) and the positive likelihood ratio (LR+) of the salivary gland US for the SS diagnosis were calculated with the AECG criteria as gold standard. ROC curves were constructed by computing the Se and Sp of different US scores in the pSS and sSS groups. For each score, the plot determined the area under the curve (AUC). The AUC was used to evaluate the diagnostic performance. We computed also the optimal cut-off US point for each score. It corresponds to the score for which the sum of sensitivity and specificity is the highest and with the most

A C C E P T E D M A N U S C R I P T significant LR+ for SS diagnosis.
The agreement of the intra-and inter-observer US results was evaluated by Cohen's non-parametric Kappa test. Kruskall-Wallis non parameter test was used to compared continue variables. Chi-square adjusted by Bonferroni method was used to compared qualitative variables. Differences were considered significant if p<0.05.

Role of the funding source

Not applicable.

Results

Characteristics of the studied population

Five patients were excluded (3 patients were not at US appointment, 1 patient did not have the questionnaire, 1 patient died). Ultimately, 97 patients were included for the analysis. Thirty-nine of them met the AECG criteria for the diagnosis of pSS and 22 patients met the AECG criteria for the diagnosis of sSS (9 patients with systemic lupus, 9 with rheumatoid arthritis, 2 with systemic scleroderma and 2 with Sharp syndrome). There were 36 control patients with sicca symptoms who did not fulfill the AECG criteria (13 patients with associated connective tissue disease and 23 without associated systemic disease).

There was no significative difference between groups regarding the mean age, sex, tabagism, use of medications causing a sicca syndrome (Table 1).

Comparison of four US scores (Table 2)

The medians for each of the four US scores (Salaffi, Hocevar, Milic, Jousse-Joulin) was significantly higher in the pSS and sSS groups than in the control group (p<0.001). The medians scores were higher in the pSS group compared to the sSS group but with no statistically significant difference.

Comparison of performance of the four scores for diagnosis of pSS (Figure 2).

Diagnostic performance was higher for the Salaffi score with AUC of 0.915 (95%CI 0.848-0.982), followed by the Jousse-Joulin score with AUC of 0.897 (95%CI 0.821-0.973), Hocevar score of 0.891 (95%CI 0.812-0.970) and Milic score of 0.885 (95%CI 0.804-0.965). The LR+ was significant (> 11.0) for the Milic, Salaffi and Jousse-Joulin scores.

Comparison of performance of the four US scores for diagnosis of sSS (Figure 3).

The Jousse-Joulin score had the best performance with an AUC of 0.851(95%CI 0.735-0.968), followed by Salaffi with AUC of 0.844 (95%CI 0.724-0.965), Hocevar with AUC of 0.824 (95%CI 0.695-0.954) and Milic with AUC of 0.808 (95%CI 0.673-0.943). The sensitivity was 65% for the Hocevar, Milic and Salaffi scores with respective optimal cut-off points ≥ 18, 3 and 7. For the Jousse-Joulin score, the optimal cut-off point was ≥ 2 and was associated with 50% sensitivity. The specificity of the four scores was excellent, between 91.7% and 94.4%. The LR+ was significant (> 11.0) only for the Milic and Salaffi scores.

Reproducibility of US diagnosis

A C C E P T E D M A N U S C R I P T
Inter-observer reproducibility was calculated for the Salaffi and Milic scores for 28 patients. The results showed low reproducibility for the accurate count of the Salaffi and Milic scores, respectively 0.274 ± 0,083 and 0.173 ± 0,069. On the other hand, the inter-observer reproducibility was better in determining the pathological character of the US as a function of the optimal cut-off point. For the Salaffi score, the inter-observer agreement showed a kappa-value of 0.714 ± 0.131 with the optimal cut-off point ≥ 7. For the Milic score, the kappa-value was 0.721 with the optimal cut-off point ≥ 3. The agreement between reader in determining the pathological character of US was 24/28 (85.7%). The intra-observer reproducibility of the Salaffi score was evaluated on 79 patients. The results showed low intraobserver reproducibility for the accurate measurement of the Salaffi score (0.516 ± 0.098 (AB) and 0.412 ± 0.064 (GC)). However, there was a better intra-observer reproducibility to determine the pathological character of US according to the optimal cut-off score. The kappa agreement was 0.68 ± 0.127 (GC) and 0,714 ± 0.132 (AB).

Discussion

Our study confirms the high diagnostic performance of salivary gland US to discriminate pSS, sSS and symptomatic control patients with 3 US scores (Milic, Salaffi, and Jousse-Joulin). The medians of the scores were significantly higher in the pSS and sSS groups than in the control group. The Milic, Salaffi and Jousse-Joulin scores were significant and clinically relevant (LR+ > 11.0) for the diagnosis of pSS. The LR+ was also greater than 11 with the Milic and Salaffi scores for the diagnosis of sSS.

For the diagnosis of pSS, the optimal cut-off scores and the AUC were in agreement with previous studies for three scoring systems (Jousse-Joulin, Hocevar and Salaffi). Jousse-Joulin et al. [START_REF] Cornec | Contribution of salivary gland ultrasonography to the diagnosis of Sjogren's syndrome: toward new diagnostic criteria?[END_REF] also had an optimal US cutoff score ≥ 2, with a lower AUC (0.817) than in our AUC study (0.897). In the study of Hocevar et al. [START_REF] Hocevar | Ultrasonographic changes of major salivary glands in primary Sjogren's syndrome. Diagnostic value of a novel scoring system[END_REF], the optimal cut-off score was the same and the AUC was 0.860 ± 0.031, just below the AUC of our study at 0.891. Zhang et al. [START_REF] Zhang | Ultrasonographic evaluation of major salivary glands in primary Sjögren's syndrome: comparison of two scoring systems[END_REF] used the Salaffi score and found the same optimal cut-off score ≥ 7 and the AUC was 0.939, slightly higher than in our study (AUC was 0.915). On the other hand, we noted a difference of results for the simplified score system of Milic et al. [START_REF] Milic | Major salivary gland sonography in Sjogren's syndrome: diagnostic value of a novel ultrasonography score (0-12) for parenchymal inhomogeneity[END_REF]. The AUC was 0.96 ± 0.01 with an optimal cut-off score  6 in their study, while we found an optimal cut-off score 3 with an AUC of 0.885.

These four studies had also used the AECG criteria for the diagnosis of pSS, which allows a good comparison of patient groups with a pSS. Lin et al. [START_REF] Lin | Cross-sectional comparison of ultrasonography scoring systems for primary Sjogren's syndrome[END_REF] also assessed 3 US scoring systems (Hocevar, Milic and Salaffi). Hocevar score had the best LR+ (LR+ 10, AUC 0.888), whereas Milic score was recommended for its excellent balance between the diagnostic value (LR+ 6.3, AUC 0.866) and inter/intra-observer agreement value. The AUC for Salaffi score was 0.851 and LR+ to 3.8 was lower. However, this study did not use the same inclusion criteria from AECG (scintigraphy and Rose Bengal) than in our study. Futhermore, results are not comparable, because Lin et al. included sSS patients in their control group.

For the diagnosis of sSS, the AUC was evaluated between 0.808 and 0.851 according to the scoring method. Very few studies have evaluated the salivary glands US during sSS. In 1992, De Vita et al. [START_REF] Vita | Salivary gland echography in primary and secondary Sjögren's syndrome[END_REF] obtained a sensitivity of 56.8% and a specificity of 92.2%, i.e. an LR+ of 6.45 compared to the control group, but with different diagnostic criteria for sSS (Copenhagen 1975). Wernicke et al. [START_REF] Wernicke | Ultrasonography of salivary glands --a highly specific imaging procedure for diagnosis of Sjogren's syndrome[END_REF] used an US score derived from De Vita et al. with AECG criteria as gold standard. They showed a sensitivity of 63.6%, a specificity of 98% and a LR+ of 31 for the diagnosis of sSS, well beyond our results.

According to our study, the major salivary glands US appears to be efficient to discriminate a sSS from a patient with connective tissue disease or a sicca syndrome. In clinical practice, the diagnosis of sSS does not change the treatment of the associated connective tissue disease (except for local treatment). A Danish study suggested anyway that the sSS had no impact on primary disease activity or on treatment [START_REF] Haga | A study of the prevalence of sicca symptoms and secondary Sjogren's syndrome in patients with rheumatoid arthritis, and its association to disease activity and treatment profile[END_REF]. However, some studies have shown that patients with an sSS had an increased risk of developing NHL (twofold) and higher mortality compared to the group of patients with rheumatoid arthritis without sSS [START_REF] Kauppi | Elevated incidence of hematologic malignancies in patients with Sjogren's syndrome compared with patients with rheumatoid arthritis (Finland)[END_REF][START_REF] Martens | Survivorship in a population based cohort of patients with Sjogren's syndrome, 1976-1992[END_REF]. The establishment of the sSS diagnosis should therefore encourage a more regular follow-up of these patients to detect a possible NHL.

A C C E P T E D M A N U S C R I P T

The high diagnostic performance of salivary gland US to discriminate pSS and sSS was demonstrated in our study even if in our group control, 13 (36%) patients had connective tissue disease without AECG criteria. A cross-sectional study [START_REF] Luciano | Ultrasonography of major salivary glands: a highly specific tool for distinguishing primary Sjögren's syndrome from undifferentiated connective tissue diseases[END_REF] consecutively included 109 patients (55 with pSS and 54 with undifferentiated connective tissue diseases) UCTD showed a higher Milic score in pSS (mean 2.2 ± 1.8 vs. 0.2 ± 0.5). The Milic score cut-off ≥ 2 showed a LR+ of 16.25 for pSS diagnosis in this study.

Few studies have evaluated the reproducibility of salivary gland US during the pSS. Salaffi et al. [START_REF] Salaffi | Salivary gland ultrasonography in the evaluation of primary Sjögren's syndrome. Comparison with minor salivary gland biopsy[END_REF] have measured the inter-observer reproducibility of their score with two blind radiologists with a high kappa agreement of 0.88 for the submandibular and 0.93 for the parotids. Hocevar et al. [START_REF] Hocevar | Ultrasonographic changes of major salivary glands in primary Sjogren's syndrome. Evaluation of a novel scoring system[END_REF] have obtained a good inter-observer agreement (with two blind observers) with a global kappa of 0.90. A high level of agreement between the two observers was seen in evaluation of echogenicity, inhomogeneity and the presence of hypoechoic areas (k = 0.88, 0.90, and 0.88, respectively). El Miedany et al. [START_REF] El Miedany | Quantitative ultrasonography and magnetic resonance imaging of the parotid gland: can they replace the histopathologic studies in patients with Sjogren's syndrome?[END_REF] evaluated intra-observer reproducibility in only 13 patients (7 pSS and 6 controls) and using another US score (0-3) about the parotid glands. Patients were examined twice. Kappa value was 0.8. In the study of Salaffi et al. [START_REF] Salaffi | Ultrasonography of salivary glands in primary Sjogren's syndrome: a comparison with contrast sialography and scintigraphy[END_REF], intra-observer reproducibility was evaluated for only 12 patients with sicca symptoms without pSS or with pSS, for the parenchymal homogeneity (k = 0.852), echogenicity (k = 0.821), size of the glands and posterior glandular border clearness (k = 0.804). The second reading was evaluated from the recorded images. Our inter-observer reproducibility was lower on the overall score, suggesting the importance of a common learning in measuring the scores. The need of developing a consensual system is the main objective of the international US-pSS study group who reported on acquisition images, a substantial parotid glands inter-observer reliability (κ = 0.62) for echogenicity and moderate (κ = 0.52) for homogeneity. Submandibular glands interobserver reliability was moderate for homogeneity (κ = 0.46) and fair for echogenicity (κ = 0.38) [START_REF] Jousse-Joulin | Salivary gland ultrasound abnormalities in primary Sjögren's syndrome: consensual US-SG core items definition and reliability[END_REF]. We did not find an association in our study between the duration of evolution of dry syndrome in Sjogren's syndrome and a high US score of the salivary glands. Salivar glands US abnormalities assessed (Jousse-Joulin score) did not change significantly during a follow-up (nearly 2 years after an initial evaluation) for suspected pSS [START_REF] Gazeau | Time-course of ultrasound abnormalities of major salivary glands in suspected Sjögren's syndrome[END_REF].

Mossel E et al. [START_REF] Mossel | Scoring hypoechogenic areas in one parotid and one submandibular gland increases feasibility of ultrasound in primary Sjögren's syndrome[END_REF] have even proposed to score hypoechoic areas in a parotid and submandibular gland increases the feasibility of ultrasound in primary Sjögren's syndrome. The performances for the positive diagnosis are not diminishing. However, this practice in daily clinic ask question because differential diagnoses isn't taking into account such as salivary gland lymphomas, other tumors, lithiases or other inflammatory diseases of the salivary glands as sarcoidosis which, in our opinion, require evaluation of the 4 main salivary glands. The score of Jousse-Joulin seems to have a good balance between simplicity and diagnostic performance based on our results. However, the cysts were the only criterion of heterogeneity of the glandular parenchyma and not the fat involution, quite frequent in the submandibular glands. This fat involution is characterized by a difficulty in visualizing the salivary gland in its environment with disappearance of the gland borders and the presence of numerous thick hyperechoic bands. These characteristics are not taken into account in the Jousse-Joulin score resulting in a scoring of 2, whereas the scoring in the Salaffi score of the glans would be 4 (Appendix A, Supplementary data S2). Finally, the score proposed by Salaffi et al. [START_REF] Salaffi | Salivary gland ultrasonography in the evaluation of primary Sjögren's syndrome. Comparison with minor salivary gland biopsy[END_REF] has the advantage of having several items to better define the parenchymal echostructure (inhomogeneity, hypoechogenic areas, echogenic bands, cysts, visibility of the posterior glandular border), contrary to the Milic score which proposes a less detailed definition in 3 grades with a higher risk of subjectivity, and the Jousse-Joulin score which preferentially evaluates the cystic degeneration of the glandular parenchyma. The Hocevar score is very detailed but longer to carry out, with no supplementary performance compared to the scores of Salaffi and Jousse-Joulin, according to our results.

The US of the major salivary glands showed a good diagnostic performance to differentiate a primary or a secondary SS from a control patient with sicca symptoms. According to some authors, it could replace sialography and salivary scintigraphy within AECG criteria without modifying them [START_REF] Salaffi | Ultrasonography of salivary glands in primary Sjogren's syndrome: a comparison with contrast sialography and scintigraphy[END_REF][START_REF] Milic | Ultrasonography of major salivary glands could be an alternative tool to sialoscintigraphy in the American-European classification criteria for primary Sjogren's syndrome[END_REF]. Cornec et al. [START_REF] Cornec | Contribution of salivary gland ultrasonography to the diagnosis of Sjogren's syndrome: toward new diagnostic criteria?[END_REF] showed

A C C E P T E D M A N U S C R I P T
that the addition of US to the AECG criteria increased the sensitivity of the AECG criteria from 77.9% to 87% without altering the specificity (96.1% vs. 98.7%) for the pSS diagnosis. They also found that the addition of salivary gland US to the 2012 American College of Rheumatology (ACR) classification criteria increased their sensitivity from 64.4% to 84.4% without modifying the specificity (89.3% vs. 91%) [START_REF] Cornec | Salivary gland ultrasonography improves the diagnostic performance of the 2012 American College of Rheumatology classification criteria for Sjögren's syndrome[END_REF]. Takagi et al. [START_REF] Takagi | Ultrasonography as an additional item in the American College of Rheumatology classification of Sjögren's syndrome[END_REF] suggested that the US could replace all ACR criteria, but without increasing the diagnostic performance of the 2012 ACR criteria, otherwise US can improve the accuracy of the integrated ACR/EULAR system and US score (93%) by compared to ACR/EULAR criteria only (86%) [START_REF] Takagi | Combined classification system based on ACR/EULAR and ultrasonographic scores for improving the diagnosis of Sjögren's syndrome[END_REF]. In a study [START_REF] Goff | Comparison of 2002 AECG and 2016 ACR/EULAR classification criteria and added value of salivary gland ultrasonography in a patient cohort with suspected primary Sjögren's syndrome[END_REF] with 290 patients prospectively included (125 (43%) met ACR/EULAR criteria and 114 (39%) also met AECG criteria), US salivary gland was abnormal in 12% of the 165 patients fulfilling no criteria set. Including salivary gland US score among the ACR/EULAR criteria increased sensitivity from 87.4% to 91.1% when physician diagnosis was the reference standard.

In conclusion, the salivary gland US is a simple, non-invasive and efficient method for the diagnosis of pSS and sSS. However, setting up a validated consensual score, and training practitioners to evaluate this US score will still be necessary to improve its reproducibility, before using it in clinical practice within diagnostic criteria.

Figure 1 .

 1 Figure 1. Examples of Parotid Gland ultrasound in short axis. A. Normal parenchym (echogenecity comparable to the thyroid) with no inhomogeneity no hypoechogenic areas, no hyperechogenic band reflexions and regular defined borders. B. Echogenecity was comparable to the thyroid with mild inhomogeneity, a few hypoechogenic areas (< 2mm), a few hyperechogenic band reflexions and partly defined borders. (This gland was scored Hocevar = 4, Salaffi = 2, Milic = 1, and Jousse-Joulin = 2). C. Echogenecity was comparable to the thyroid with evident inhomogeneity, several hypoechogenic areas (2-6 mm), a several hyperechogenic band reflexions and partly defined borders. (This gland was scored Hocevar = 7, Salaffi = 3, Milic = 2, and Jousse-Joulin = 3). D. Echogenecity was decreased to the thyroid parenchym with grossly inhomogeneity, numerous hypoechogenic areas (2-6 mm), numerous hyperechogenic band reflexions and ill-defined borders. (This gland was scored Hocevar = 12, Salaffi = 4, Milic = 3, and Jousse-Joulin = 3).A C C E P T E D M A N U S C R I P T

Figure 2 :

 2 Figure 2: Performance for each of the four US scores for the pSS diagnosis in comparison to the control group. A: table with AUC, optimal cut-off point, Sensitivity (Se), Specificity (Sp) and Positive Likehood Ratio (LR+) for each of the four scores. B: ROC curves of the four US scores for pSS diagnosis.

Figure 3 :

 3 Figure 3: Performance of each of the four US scores for the sSS diagnosis in comparison to the control group. A. table with AUC, optimal cut-off point, Sensitivity (Se), Specificity (Sp) and Positive Likehood Ratio (LR+) for each of the four scores. B: ROC Curves of the four US scores for sSS diagnosis A.

Table 2 .

 2 Median of salivary gland US scores for each group 0-0.0] a,b [0.0-0.0] a,b [0.0-0.0] a,b pSS= primary Sjögren's syndrome; sSS= secondary Sjögren's syndrome. IQR: InterQuartile Range Kruskal-Wallis non-parametric test. a,b p< 0.001.

  table with AUC, optimal cut-off point, Sensitivity (Se), Specificity (Sp) and Positive Likehood Ratio (LR+) for each of the four scores. B: ROC curves of the four US scores for pSS diagnosis.

	A.
	US scores	AUC	Optimal	Se (%)	Sp (%)	LR+
		(CI95%)	cut off	(CI95%)	(CI95%)
			point
	Hocevar	0.891 (0.812-	 17	61.5 (45.9-	91.7 (78.2-	7.4
	Milic	0.970)	 3	75.1)	97.1)	11.8
	Salaffi	0.885 (0.804-	 7	65.8 (49.9-	94.4 (81.9-	11.3
	Jousse-	0.965)	 2	78.8)	98.5)	11.8
	Joulin	0.915 (0.848-		63.2 (48.1-	94.4 (81.9-
		0.982)		76.3)	98.5)
		0.897 (0.821-		65.8 (49.9-	94.4 (81.9-
		0.973)		78.8)	98.5)

Table 1 : Characteristics of the study population

 1

	pSS	sSS	Controls
	n (%)	n (%)	n (%)

Acknowledgments

The authors would like to thank the patients who participated in the study and the Doctor Christophe Boulnois for his proofreading and translation work in the English language.

Contributorship

All listed authors have had substantial contributions in the acquisition, analysis, and interpretation of data for the work; they were involved in revising it critically and give final approval for the version to be published. All authors agree to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Conflicts of interest

The authors do not declare any conflict of interest with this work.

A C C E P T E D


Aurélien Perera

Bernarda Lovrinčević

A comparative study of aqueous DMSO mixtures by computer simulations and integral equation theories

Several computer simulation studies of aqueous dimethylsulfoxyde with different force field models, and conducted by different authors, point out to an anomalous depressing of second and third neighbour correlations of the water-water radial distribution functions. This seemingly universal feature can be interpreted as the formation of linear water clusters. We test here the ability of liquid state integral equation theories to reproduce this feature. It is found that the incorporation of the water bridge diagram function is required to reproduce this feature. These theories are generally unable to properly reproduce atom-atom distribution functions. However, the nearideal Kirkwood-Buff integrals are relatively well reproduced. We compute the Xray scattering function and compare with available experimental results, with the particular focus to explain why this data does not reproduce the cluster pre-peak observed in the water-water structure factor. 1

Introduction

Water remains an important subject of study, because the exact nature and role of the hydrogen bond network, in relation to its peculiar properties, remains elusive [1]. This role becomes even more complicated to understand when it comes to water mixed with various types of solutes [1]. In this context, aqueous-dimethylsulfoxyde(DMSO) mixtures have attracted considerable attention, particularly the fact that this mixture has an exceptionally low freezing point depression, around T=-140 o C, and for DMSO mole fractions around x=0.33 [2]. The understanding of this property is believed to impact our current knowledge in cryoprotectancy [3], particularly in relation to the hydrogen bonding properties of water. Because of this special interest, most of the earlier studies seem to have focused in the water-DMSO clusters, as even very recent studies show the same interest [4,5,6,7,8,9]. DMSO has been modeled differently by different authors and there are numerous subsequent computer simulation studies of these particular mixtures [10,11,12,13,14,15,16,17,18,19,20]. Many of these investigations have shown that specific water-DMSO clusters occur, even at room temperature far from those related to cryogenic conditions.

In our earlier investigation of this special mixture [21], we found that, despite the various possible water-DMSO force field model combinations, all of them indicate that the water-water correlations show the same typical feature, namely a depression of the second and third neighbour correlations. This is illustrated in Fig. 1, where the main panel clearly shows that second and third neighbour correlations show oscillations around 1 with marked minima, in contrast with that of aqueous ethanol (dotted line) which is above 1, as usual with several aqueous mixtures. Albeit typical variations due to small differences between DMSO models, the universality of this feature is quite visible. While the first peak witness the strong water-water hydrogen bond, which is seen in any mixtures, the depletion of farther neighbours suggests that water has to satisfy 2 constraints: being hydrogen bonded to itself and loose the favoured tetrahedral conformation of bulk water. One solution to this problem is the formation of water string-like clusters. Indeed, such clusters are clearly visible in snapshots from simulations, and at all DMSO mole fractions x in the range 0.1 < x < 0.9, as we have been the first to point out their existence [21]. Since various types of experiments show more clearly the existence of water-DMSO bonds [8,9], one could induce that linear-water clusters must be a consequence of the first more readily visible property. Indeed, all the force field models of DMSO share the fact that the oxygen site is charged negatively while all the other sites, sulfur and the 2 methyl sites, are charged positively.

It is quite easy to visualise that the negatively charged water oxygen site would bind to these latter DMSO sites, leaving the 2 hydrogen bonds open to the bulk, hence favouring neighbouring water oxygens to attach to them. This situation will clearly disfavour tetrahedral binding of water. This is what computer simulations seem to indicate. This situation is very different from other aqueous mixtures, as for example, water-alcohol mixtures, where water is seen to form globular-like self segregated domains [22]. In Ref. [21], we have argued that it is the large size of the DMSO atom which allows water oxygen to attach to it. However, another puzzle arises: why available Xray scattering data for this mixture [23] does not show a pre-peak, which should arise from the existence of such water clusters. Indeed, the presence of such clusters induces a clear pre-peak at about k ≈ 0.75Å -1 in the water-water structure factor that we have reported in Ref. [21]. Hence, this prepeak should, in principle, contribute to the total scattered intensity. In the present paper, we explain the absence of scattering pre-peak in narrow connection to the concept of domain order introduced in Ref. [24], which leads to a compensation of the positive pre-peak in like species structure factors by the negative pre-peak of the cross species structure factors. This is always the case for many aqueous mixtures, which show the same universal local domain order. In order to break this "symmetry" and induce a scattering pre-peak, it is necessary that the interface separating the segregated water and solute domains acquires particular properties, which would enhance its contribution to the total scattering. This is the case of micelles, for example, for which water is seen to accumulate in the coronary area. We conjecture here that this is the reason why aqueous mixtures of small monools and 1-2-diols do not show scattering pre-peaks, while these tend to appear in triols, which have ester groups near the surfactant OH head group, and which attract water in the corona of the micelle.

Concerning liquid state integral equation theories (IET) approaches to realistic molecular liquids, these have been intensively tested in the early 90's, but have generally lead to disappointing results for associated liquids [25,26]. To make this point clear, both the molecular IET and the site-site IET are enable to describe the tetrahedral bonding of water through the water-water pair correlation function [26,27,28]. Furthermore, both these methods lead to a very similar incorrect shape of the radial distribution function of water [26,27,28], thus clearly demonstrating that it is more the nature of the closure approximation, rather than the use of the molecular versus site-site formalisms, which is responsible for these shortcomings. The main point we want to make here is that these theories are even less successful for mixtures involving associated liquids, precisely because association induces domain segregation, which is a more complex property to describe. In particular, we want to show that correcting the deficiencies of pure liquid description, by introducing bridge diagram for water, cannot account for the dramatic change in structure of segregated mixtures. The particular interest to test these theories for the case of aqueous DMSO, is the fact that this mixtures shows quasi-ideal Kirkwood-Buff integrals, unlike most other aqueous mixtures [29]. In that, this particular system seems to have less concentration fluctuation, and therefore should be more amenable to mean-field like approaches such as IET. In Ref. [21], we have argued that this feature of the aqueous DMSO mixtures was related to water linear aggregates forming pseudo molecular entities, leading to a weakly interacting pseudo ideal mixture of DMSO with these aggregates. Since IET are quite good for weakly interacting and fully disordered mixtures, it would be interesting to see how much of this property can be captured by these theories. Indeed, there are 2 problems here: the first is low concentration fluctuations, which is favourable for a IET description, and the second is the description of hydrogen bonded aggregates, which should be a challenge for IET.

The remainder of this paper is as follows. In the next section, we briefly remind the models used, simulation protocols and theoretical details related to the IET and their use in the present case. Details concerning the concept of domain order and calculation of the scattering functions are equally gathered in the last subsection. Section 3 details our results for the structural functions and the KBIs.

Finally, Section 4 gathers our conclusions and perspective on the present work.

2 Theoretical details 2.1 Models and computer simulations

We have essentially followed the methodology reported in Ref. [21]. Most of the data shown here is from the same simulations presented in this work. Additional structure factors have been computed in order to allow for the calculation of the scattering functions. To summarize, we have used the SPC/E water model [START_REF] Berendsen | Hermans[END_REF] and the VB DMSO model [START_REF] Vaisman | [END_REF]. We have also used in Fig. 1 the OPLS DMSO model [32]. N = 2048 particles simulations were used in the isothermal and isobaric molecular dynamics simulations with the Gromacs package [33]. Specific details can be found in Ref. [21].

Integral equation theory

We have used the site-site IET formalism, which has the same level of accuracy as the molecular version, as far as associated liquids are concerned. Although one would expect a superior accuracy of the molecular formulation, due to an exact version of the molecular Ornstein-Zernike (MOZ) equation [34,[START_REF] Hansen | Theory of Simple Liquids[END_REF], instead of the approximate site-site Ornstein-Zernike (SSOZ) version [25,[START_REF] Hansen | Theory of Simple Liquids[END_REF], it is the constraint imposed by closure relation which dominates the structural description. Indeed, the description of hydrogen bonded structures in the associated liquids is particularly constraining, requiring 3-body and higher order correlations present in the missing highly connected and irreducible bridge diagrams [START_REF] Hansen | Theory of Simple Liquids[END_REF]. This is the reason for the similarity in the poor description of both approaches for water from MOZ [27] and SSOZ [28].

The 2 equations to solve are the site-site OZ matrix equation together with a closure equation. The SSOZ equation reads

SM = I (1)
where

S = W + ρH (2)
is the generalized structure factor matrix, with elements defined as

S i α j β (k) = W i α j β (k) + ρH i α j β (k) (3)
where the index i α and j β designate atom i of molecular species α and atom j of molecular species β , respectively, and

M = W -1 -ρC (4)
The matrix W describes the intramolecular correlation, defined as rigid atomic bonds within molecular species

W i α j β (k) = sin(kd i α j α ) kd i α j α (5)
where d i α j α is the distance between atomic sites i α and j α of molecule of species α. The matrix H is related to the site-site distribution functions g i α j β (r), with

h i α j β (r) = g i α j β (r), -1, and hi α j β (k) = ´d rh i α j β (r) exp(i k. r
) the Fourier transform of h i α j β (r), and the matrix C is related to the Fourier transform of the site-site direct correlation function c i α j β (r) with the relations

H i α j β = x α x β hi α j β (k) C i α j β = x α x β ci α j β (k)
where x α is the mole fraction of species α.

For the closure equations, we have chosen both the hypernetted chain (HNC) equation and the Kovalenko-Hirata (KH) closure [START_REF] Kovalenko | [END_REF]. The HNC closure equation reads

g i α j β (r) = exp -β v i α j β (r) + h i α j β (r) -c i α j β (r) + b i α j β (r) (6)
where v i α j β (r) is the interaction between the atomic sites i α and j β (β = 1/k B T is the Boltzmann factor), and the last term b i α j β (r) represents the so-called bridge function. Setting b i α j β (r) = 0 leads to the strict HNC closure, while this function can also be tailored from the site-site correlation functions obtained through the procedure we have outlined in Ref. [37]. We will refer to this modified HNC as HNC+B.

The KH equation is a mix between the strict HNC equation above and the mean spherical approximation, based on the value of the term

Γ i α j β (r) = -β v i α j β (r)+ h i α j β (r) -c i α j β (r)
in the exponential of Eq.6:

g i α j β (r) =    exp Γ i α j β (r) when Γ i α j β (r) < 0 1 + Γ i α j β (r) when Γ i α j β (r) > 0 (7)
This ad-hoc procedure of damping the raise of the correlations in the exponential through the term Γ i α j β (r) allows to control its dramatic growth in the case of HNC, which leads to the no-solution and spinodal-like scenario often spuriously encountered in the numerical solution of these approximate closures. For this reason, this particular KH closure has found a very wide and convenient use in the community. We will discuss the occurrence of this type of scenario in the present case, below in the results sub-sections.

Finally, we have equally examined the Percus-Yevick (PY) closure, which is a linearised version of HNC

g i α j β (r) = exp -β v (SR) i α j β (r) 1 + +h i α j β (r) -c (SR) i α j β (r) (8)
where the superscript (SR) indicates that only the short range interactions and direct correlations are considered. The long range Coulomb part is removed from both functions following a standard procedure described in the literature [38]. In fact, this is also the case in the HNC and KH closures, because of the exact limiting relation:

lim r→∞ c i α j β (r) = -β v i α j β (r)
which allows to retain only the short range part of the direct correlation function in the various closure equations.

As shown in the Results section, both the HNC and the PY equations cannot be solved for the mixture case, precisely because of the behaviour of the term h i α j β (r) -c (SR) i α j β (r), which induces the growth of medium range correlations, which, in turn, promotes spurious phase separation-like behaviour. The damping procedure of the KH closure appears then as a wise method. However, we will discuss below the appropriateness of this artificial method in view of the physical phenomena in hydrogen bonding mixtures.

Both the SSOZ equation and the closure equations are solved iteratively for the unknown site-site functions h i α j β (r) and c i α j β (r). These functions are discretised on a grid of 2048 points, with r-spacing of ∆r = 0.03Å.

We find that the KH closure yields numerical solution for the aqueous DMSO mixtures for all concentrations at room temperature. Conversely, neither the PY nor the HNC closures are able to yield numerical solutions at T = 300K. While we do not expect that the linearised PY closure to be able to handle complex structures, the HNC closure deserves some comments.

The strict HNC closure is unable to yield numerical solution below very high temperatures around T = 5000K for the equimolar mixture. At such elevated temperature, HNC predicts high small-k raise of the structure factors, although no obvious long range tail seem to develop in the site-site correlations. This is a pathology of HNC closure, as noted early in the historical development of IET [39]. The mathematical origin of such behaviour is unclear at present, and it is generally admitted that this behaviour underlies a spurious spinodal [40]. However, as we have mentioned above, we do not observe that any corresponding long range tail develops in the correlation functions. At such high temperature where we can find numerical solutions, the site-site functions are nearly structureless, just like expected for any thermally disordered fluid. In fact, comparing with the KH closure results at the same temperature reveals that the first neighbour water-water correlations are overestimated in the HNC closure, leading to the k = 0 increase. In other words, it is the short range behaviour of the correlation which leads to the k = 0, and not a spinodal behaviour, contrary to the generally accepted idea about this closure.

Domain order and the scattering function

The X-ray scattering intensity I(k) can be obtained through the Debye formula [41]

I(k) =< ∑ i α , j β f i α (k) f j β (k) exp ik.(r i α -r j β ) > (9)
where the sum runs over all pairs of scattering atoms i α and j α , which are at respective spatial positions r i α and r j β , the functions f i α (k) are the atomic form factor for atom i α and depend on the type of radiation which is scattered (we use here Xray), and the symbol <...> designates an average over all possible positions of these atoms, which corresponds to a thermal average, or an ensemble average for calculational purposes. It can be shown [24,42] that this expression can be rewritten in a form which can be used in theoretical methods

I(k) = I Ideal (k) + I Intra (k) + I Inter (k) (10)
where the 3 contribution represent, respectively the ideal part

I Ideal (k) = ∑ α x α ∑ i α f i α (k) (11)
the intra-molecular part

I Intra (k) = ∑ αβ x α x β ∑ i α j β f i α (k) f j β (k)W i α j β (k)(1 -δ i α j α δ αβ ) (12)
where W i α j β (k) is given by Eq.( 5), and inter-molecular part

I Inter (k) = ∑ αβ x α x β ∑ i α j β f i α (k) f j β (k)H i α j β (k) (13)
The Kronecker delta terms in Eq.( 12) avoid the self-atomic contributions inside same species; in which case it is the ideal term in Eq.( 11) which prevails. The ideal contribution corresponds to non-interacting and non-bonded atoms. The "intra" part brings in the fact that atoms are bonded into molecules. Finally the "inter " part accounts for the correlations, hence interaction between atoms and molecules.

The Xray intensity, as written through Eqs.(10-13), can be calculated for both the cases of simulations and IET, since the correlation functions g i α j β (r) are available.

Results

Chain-like water clusters

As shown in Ref. [21], the computer simulations of different combinations of model water and model DMSO indicate that water tends to form predominantly chain-like clusters. This is illustrated in Fig. 2. This feature is not the one which is usually reported in the literature, where the accent is more on the water-DMSO dimer formation. However, these two properties are not necessarily excluding each other. What is interesting is to know which property induces the other. One view point is that, water tends to bind preferentially to water, which then dictates subsequently all other structural features. This property is reflected in the model simulation through the fact that the partial charges on the oxygen of water are quite high (q O W ≈ -0.85 for the oxygen atom, for almost all water models), and higher than for most solutes (for DMSO, the charge on the oxygen atom is q O D ≈ -0.46 and on the sulfur atom q S D ≈ 0.14, for acetone q O A ≈ -0.56, and for mono-ols q O ≈ -0.7). A logical consequence of electrostatic is that charges with highest valence will tend to dominate the interactions, leading weaker charges to adapt. Following this line of reasoning, water is more likely to bind with itself, and preserve pockets of tetrahedrally bound water molecules. This is what happens with monools [43] and acetone [44]. We would then expect something similar with DMSO as well, and at even greater scale, since the valence is smaller. This is not what is observed in Fig. 2: water seems to mix rather well with DMSO, precisely by forming chain-like clusters. Clearly, the competition between charges alters the simplistic picture presented above. High valences on the solute equally favours self solute binding, thereby enforcing the micro-segregation observed in simulation of aqueous-acetone [44] or aqueous-alcohol [24,43]. Intuitively, weaker charges on the solute should increase water segregation. This is what we observed in aqueous mixtures of "weak-water" models, where the valence of the latter weak-water was decreased to various values. Therefore, the reason for the enhanced miscibility of water in these mixtures is to be searched elsewhere.

Structure functions

We look at the site-site correlation functions and associated structure factors of pure DMSO and aqueous DMSO, as obtained from computer simulations and IET.

Pure water

The case of pure water has been reported in several previous works by many authors [25,27,28], including our own work [37]. The comparison with scattering experiments shows that classical force field such as SPC/E are quite adequate. However, structure prediction from IET cannot capture the strong tetrahedral ordering of water, and predict a structure closer to that of a simple Lennard-Jones liquid. Interestingly, and as stated above, both the molecular and the site-site produce very similar correlations, which indicates that the principal issue is not in the differences between the two methods, but the fact they both miss many body correlations. In the present work, the water-water correlations for mixtures are corrected by the pure water bridge terms, which we have computed in Ref. [37]. We show that this contribution is capital to produce the chain-like water clusters within the HNC theory.

Neat DMSO

Fig. 3 shows selected atom-atom correlation functions and corresponding structure factors for pure DMSO, as obtained from simulations (curve in green), HNC (blue), KH(red dashes) and PY(magenta). It is seen that the first neighbor correlations are generally underestimated by the theories, but not dramatically so.

The PY theory is not able to properly describe the features of the structure factors. But the principal problem with the theories is seen in small-k behaviour of the structure factors, which all show a pronounced raise, suggesting liquid-vapour pre-transitional prediction. This is most serious for the HNC closure. A close look at the long range behaviour of the correlation functions shows a general tendency to be above the asymptote 1.

Aqueous-mixtures

Selected correlations between the oxygen atoms of water and DMSO molecules, as well as associated structure factors, are displayed for typical DMSO mole fractions x=0.1, 0.5 and 0.8, in Fig. 4, Fig. 5 and Fig. 6, respectively. The color conventions are similar to Fig. 3, simulation results in green, HNC in blue and KH in red dashes. What is immediately apparent is how poor IET results are. Despite the fact that water-water correlations contain the bridge function which reproduces pure water data accurately, it is not sufficient to enforce a better agreement than the KH closure. In fact, this closure seems to reproduce the trends of the correlations for the 10% DMSO case in Fig. 4a, while the HNC+B closure is very poor. The situation is somewhat improved for the 50% and 80% DMSO cases. Perhaps the most important feature is seen in the structure factors in the small-k region, where the HNC+B closure version tends to predict the cluster pre-peak visible on the simulation data at all concentrations. The KH approximation, despite being somewhat more accurate, predicts a small-k raise at k=0 for the water-water correlations. In fact, the strict HNC closure cannot be solved for this mixture at room temperature, precisely because of spurious k=0 raise in the correlations. The partial linearisation which is made in the KH closure in Eq. (7) prevents this raise to develop in the same dramatic fashion as for the strict HNC closure. However, it is obvious that this "working recipe" does not explain how and why the inclusion of proper bridge terms would allow to obtain numerical solutions inaccessible to the HNC closure. It is only through the introduction of higher order correlation that this problem can be cured. This observation indicates how crucial such correlations are needed in order to better represent the complexity of associating mixtures.

As indicated in our previous study of aqueous-DMSO mixtures, the most important feature of the correlations is the presence of a pre-peak in the water-water structure factors, in the k-vector range k ≈ 0.7 -1Å -1 . This is due to chain-like water aggregates, as indicated by the second and third neighbour depression in the water-water correlations. These positive pre-peak contributions are accompanied with negative pre-peak contributions in the water-DMSO cross correlations. This was clearly illustrated in Fig. 7 of Ref. [21]. The presence of positive and negative pre-peaks is not specific to the present system, and are equally found in room temperature ionic liquids [45]. They correspond in fact to an underlying charge ordering between positively and negatively charged atomic groups. Despite the strong local inhomogeneity induced by such charge order, which is responsible for the micro-heterogeneity as well as the linear water aggregates found here, the global homogeneity is preserved. A signature of this global homogeneity is the perfect compensation of these positive and negative pre-peak in the scattered intensity, as we show in sub-section 4 below.

Charge order

An important property of the structure of neat DMSO is revealed in Fig. 7. It shows that the long range correlations between the oxygen and sulfur atoms with methyl atoms are out-of-phase with all the types of atom-atom the correlations. This property is typical of charge ordering [46], as recently advertised in room temperature ionic-liquids [47,48]. Charge order essentially describes the alternated disposition of positive and negatively charged atoms, most appropriately in ionic melts. But this order can equally occur between charged atomic groups in complex room temperature ionic liquids [47,48]. The hydrogen bonding in associative liquids is such a nice example of charge ordering, alternating oxygen and hydrogen atoms, which in this case enforce the tetrahedral ordering in water, or the chain-like ordering in mono-ols. Therefore, what we observe in DMSO is an alternance of oxygen and methyl. We conjecture here that it is this ordering which imposes the ordering of the water molecules, which then are not obliged to cluster into globular clusters, but rather meander between the DMSO molecules through linearly bound water chains. Indeed, the upper insert of Fig. 7 shows that charge order persists under aqueous dilution. The lower insets of Fig. 7 shows that HNC is not able to describe the strength of this charge ordering within the same spatial range, and considerably underestimate it. This could be due to high order diagrams missing. We did not try incorporating such diagrams by the procedure which we used in other cases, for the reason explained in the next sub-section.

The KH approximation, although not shown here, shows a similar behaviour.

Kirkwood-Buff integrals

Fig. 8 shows the KBI from experiments (filled dots) compared with simulations (open triangles) and IET results, HNC (full lines) and KH( dashed lines). We note that generally, the agreement between theoretical and experimental results are rather good. This is even more surprising when the IET results for structure are quite disappointing. We note that the KH results are very accurate compared to the experiments, except for G DD for small DMSO content. The HNC+B closure results are worse for G WW at high DMSO content. In fact, they seem to follow the results of the simulations rather well, which is intriguing. Indeed, the force field model we have used are exactly the same for both theoretical methods. Despite these shortcomings, this thermodynamic property related to structure is rather well reproduced. This is in sharp contrast with the KBI results for other aqueous mixtures, such as water-alcohol, for example, for which the KBI from IET are very poor. Comparing the nature of the aggregation in both systems, since aqueous alcohol mixtures tend to present quite large water-water KBI in relation with strong micro-heterogeneity and globular water aggregates, while aqueous DMSO has rather small micro-heterogeneity and linear water aggregates, it would seem that IET are better suited for weakly aggregating mixtures. This is in line with our previous work with aqueous weak-water mixture models [49].

Scattering function

Fig. 9a-b show the Xray scattering intensity from all 3 theoretical methods compared with experimental data for pure water and pure DMSO. For water, the full intensity is represented. For DMSO, another common representation is often used by experimentalists and simulators is:

∆I(k) = k. I(k) -I ideal (k) (14)
It can be seen that, despite minor discrepancies in the correlation functions, the IET manage to reproduce the total scattered intensity rather well. Fig. 10a-c show the scattering intensity ∆I(k) for mixtures, compares with the Xray data from Koga et al. [23] for different DMSO concentrations, x=0.15 (right panel), x=0.5 (middle panel) and x=0.8 (right panel). We note that the simulations are in good agreement with the experimental data, and that the IET perform quite well, specially when compared with the rather poor structure factors shown in Figs. 456. This is certainly due to favourable compensations in the sum of Eq. (13).

Perhaps the most important feature is the disappearance of the water pre-peak from the total intensity. This is due to a perfect compensation from the water-water and water-DMSO cross contributions in Eq.( 13). The positive pre-peak of the first are exactly compensated by the negative contributions from the second. This is illustrated in Figs. 11a-c, which show the total scattered intensity I(k) (in blue), together with other contributions in Eq.( 10

Discussion and Conclusion

The present work illustrates how approximate IET, which are obviously unable to account for the local ordering due to association in aqueous mixtures, can however turn out to be relatively acceptable when computing averaged quantities accessible to experiments, such as the KBI or the radiation scattering intensity. The first concerns a spatial average, and is related to various thermodynamic properties [50] , while the second is an average through summing over all sitesite contributions. IET are relatively accurate in many situations, such as simple liquids and their mixtures [51,52,53], weakly polar systems [54,55,56], ionic melts [57]. The present results indicate that their accuracy comes essentially when 2-body interactions and correlations are predominant. If clustering is present, then these theories become useless. A supplementary hint comes from the fact that IET are very unreliable to describe critical fluctuations, which means in the vicinity of phase transitions [40]. The HNC theory is not even able to predict a spinodal behaviour [40]. Since aggregation is a k-dependent concentration fluctuation, the latter which is simply the k=0 part, we do not expect these theories to be able to describe directional clustering properly. This is the principal reason of the problems reported here. Nevertheless, the fact that they can be used for weakly aggregating systems, such as aqueous-DMSO, indicate the importance of many body correlations in various other aqueous mixtures. The multiple recipe to improve IET for spherical interaction are of no use when directional clustering is present. The present work, equally indicates that accounting for pure water association is not enough to permit to describe what happens in the case of mixtures.

The no-spinodal pathology of the HNC closure indicates that recipes such as the KH closure do not bring any deep insights into the relation between directional molecular association and many body correlations. Insights into this matter are still missing.

It is equally interesting to see that radiation scattering, more specifically Xray scattering, does not show any pre-peak related to the existence the water linear chain aggregates, due to the canceling the pre-peaks in the various species-species partial contributions (positive for like species and negative for cross species), and which witness the local structure of the aqueous-DMSO mixtures, namely the linear chain aggregates. In other words, the Xray scattering is "blind" to the specific micro-structure of this particular mixture. While this chain aggregate structure is quite apparent from the water-water correlations, and seemingly independently of the DMSO models, it is not immediately deducible from spectroscopy data, which tends to enhance water-DMSO pairing, and even less from scattering data, where the pre-peak is totally absent. In a way, this result indicates that computer simulations provide a better insight into the local microscopic structure than experimental radiation scattering can provide. Perhaps neutron scattering experiments can provide additional insight into this matter, through various deuteration processes. This point remains open for subsequent investigations. Comparison between various water-water RDF reported in the literature for aqueous-DMSO mixtures. Blue line is pure SPC/E water, red line is x=0.35 from Ref. [12], green line for x=0.33 from Ref. [16],

Figure captions

jade line is for x=0.3 from Ref. [21], orange line from OPLS DMSO, and dotted gray line for aqueous ethanol x=0.3 [59]. together with the ∆I(k) representations from Eq.( 14) Experimental results in red, simulation in blue, IET results in dotted cyan. The ideal contribution from Eq.( 11) is in green and the intra-molecular from Eq.( 12) contribution in jade. For water, the experimental data for I(k) is from Ref. [58], and ∆I(k) for both water and DMSO are from from Ref. [23]. 11) is in green and the intra-molecular from Eq.( 12) contribution in jade. For water, the experimental data for I(k) is from Ref. [58], and ∆I(k) for both water and DMSO are from from Ref. [23] .

  ), namely I ideal (k) in green line and ∆I Intra (k) = I Intra (k)+I Ideal (k) in jade. The total intensities for HNC (cyan) and KH( gold) are equally shown. In addition, the partial intensities coming from the water-water (purple) I WW (k), DMSO-DMSO(brown) I WW (k) and water-DMSO (magenta) cross contributions I W D (k) are equally shown. The total intensity I(k) is then the sum of these 3 curves: I(k) = I WW (k) + I DD (k) + 2I W D (k), and one can observe that I(k) (blue) has no signature of the positive pre-peaks of I WW (k) and I DD (k) and the negative pre-peak inI W D (k) at k ≈ 1.6Å -1, which indicates that they get exactly compensated into the total intensity I(k). It is equally important to realize that the intra-molecular contribution, shown in dashed green line, is quite important, and needs to be added to the intermolecular contribution which is solely made of the site-site functions.

Fig. 1

 1 Fig.1Comparison between various water-water RDF reported in the liter-

Fig.

  Fig. Snapshot of the 50% aqueous DMSO mixture (2048 particles) evidencing linear water clusters. Water molecules are shown as red ball for oxygen atom and white for hydrogen atoms. DMSO is shown in cyan semi-transparent.

Fig. 3

 3 Fig.3 Selected atom-atom correlation functions (a) and corresponding structure factors (b) for the neat DMSO mixture, comparing simulation data (green) to HNC(blue), KH(dashed red) and PY(thin magenta) data. The various atom pairs are indicated in the panels.

Fig. 4

 4 Fig.4 Oxygen atom correlation functions (a) and associated structure factors (b) for the 10% aqueous DMSO mixture. The color conventions are as in Fig.2.

Fig. 5

 5 Fig.5 Oxygen atom correlation functions (a) and associated structure factors (b) for the 50% aqueous DMSO mixture. The color conventions are as in Fig.2.

Fig. 6

 6 Fig.6 Oxygen atom correlation functions (a) and associated structure factors (b) for the 80% aqueous DMSO mixture. The color conventions are as in Fig.2.

Fig. 7

 7 Fig.7DMSO atom-atom correlation functions demonstrating the long range charge order behaviour (see text). g OO (r) in blue, g SS (r) in purple, g SM (r) in black, g MM (r) in green, g OS (r) in red and g OM (r) in orange. The main panel is for neat DMSO, the upper two insets for 80% and 50% aqueous DMSO mixtures, and the lower inset for the HNC results.

Fig. 8 Kirkwood-

 8 Fig.8

Fig. 9

 9 Fig.9 Xray scattering intensities I(k) for neat water (a) and neat DMSO (b),

Fig. 10

 10 Fig.10 Comparison between experimental and simulated Xray scattering intensities ∆I(k) (main panel) and I(k) (inset) for aqueous mixtures, with 15% (a), 50% (b) and 80% (c) DMSO contents. Experimental results (from Ref.[23]) in red, simulation in blue, HNC results in cyan, and KH in dashed jade.

Fig. 11 Fig. 3 -Fig. 4 - 5 - 6 - 9 -

 1134569 Fig.11 Various contributions to the Xray scattering intensities I(k) (shown in thick blue line) corresponding to those shown in Fig.10, with 15% (a), 50% (b) and 80% (c) DMSO contents. The ideal and intra-molecular contributions are shown in dashed green and jade, respectively. The purple curve represents the water-water contribution I WW (k) to I(k), the brown curve the DMSO-DMSO I DD (k),and the magenta curve is for the cross water-DMSO I W D (k) contribution (see text). The vertical yellow arrows indicates the positions of the pre-peaks due to water aggregates (see text). The HNC and KH data for I(k) are shown in cyan and dashed gold lines, respectively.

11 -

 11 Various contributions to the Xray scattering intensities I(k) (shown in thick blue line), corresponding to those shown in Fig.10, with 15% (a), 50% (b) and 80% (c) DMSO contents. The ideal and intra-molecular contributions are shown in dashed green and jade, respectively. The purple curve represents the water-water contribution I WW (k) to I(k), the brown curve the DMSO-DMSO I DD (k),and the magenta curve is for the cross water-DMSO I W D (k) contribution. The vertical yellow arrows indicates the positions of the pre-peaks due to water aggregates (see text). The HNC and KH data for I(k) are shown in cyan and dashed gold lines, respectively.

Acknowledgments

We are pleased to contribute to this special issue in honor of Daan Frenkel. This work is funded in parts by the Croatian Science Foundation (Installation Research Project UIP-11-2013-4514).

.