Marine Fauny
email: m.fauny@chru-nancy.fr

Caroline Morizot
email: c.morizot@chru-nancy.fr

Edem Allado
email: e.allado@chru-nancy.fr

Frank Verhoeven
email: fverhoeven@chu-besancon.fr

Eliane Albuisson
email: e.albuisson@chru-nancy.fr

Marie Semaan
email: marie.semaan@balamand.edu.lb

Astrid Pinzano
email: astrid.pinzano@univ-lorraine.fr

Isabelle Chary-Valckenaere

Damien Loeuille
email: d.loeuille@chru-nancy.fr

Isabelle Chary- Valckenaere

Consequences of spinal ankylosis on bone trabecular fragility assessed on CT scans in patients with ankylosing spondylitis. A retrospective study

Keywords: Ankylosing spondylitis, osteoporosis, vertebral fracture, bone mineral density, CT scan, bone attenuation coefficient, ankylosis, mSASSS

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Ankylosing spondylitis (AS) is an inflammatory rheumatic disease that mostly affects young adults. The symptomatology often concerns the axial skeleton and, less frequently, the peripheral skeleton (hands, feet) and the enthesis, with night pain and a prolonged morning stiffness. The pathophysiology of bone formation (syndesmophytes, enthesophytes, bone bridges) is complex [START_REF] González-Chávez | Molecular mechanisms of bone formation in spondyloarthritis[END_REF][START_REF] Wendling | New bone formation in axial spondyloarthritis[END_REF] and coexists with bone resorption (bone erosion, osteoporosis). Due to this bone resorption, the prevalence of vertebral fracture (VF) and of its neurological complications is increased in the AS population (odds ratio 3.26) [START_REF] Maghraoui | Extra-articular manifestations of ankylosing spondylitis: prevalence, characteristics and therapeutic implications[END_REF]. AS patients have a higher risk of osteoporosis, with a prevalence varying from 11.7 to 34.4% (14 to 27% in the spine and 4 to 14% in the hip) [START_REF] Ramírez | Prevalence and risk factors for osteoporosis and fractures in axial spondyloarthritis: A systematic review and meta-analysis[END_REF][START_REF] Van Der Weijden | High prevalence of low bone mineral density in patients within 10 years of onset of ankylosing spondylitis: a systematic review[END_REF][START_REF] Ghozlani | Prevalence and risk factors of osteoporosis and vertebral fractures in patients with ankylosing spondylitis[END_REF][START_REF] Briot | Bone oedema on MRI is highly associated with low bone mineral density in patients with early inflammatory back pain: results from the DESIR cohort[END_REF][START_REF] Schett | Joint remodelling in inflammatory disease[END_REF][START_REF] Ardizzone | Ankylosing spondylitis and osteoporosis[END_REF][START_REF] Klingberg | Osteoporosis in ankylosing spondylitis -prevalence, risk factors and methods of assessment[END_REF]. A significant bone loss was reported in 51% of patients with early spondyloarthritis at either lumbar spine or hip [START_REF] Fechtenbaum | Baseline MRI inflammation is not a determinant of 5-year bone mineral density loss in patients with early spondyloarthritis[END_REF]. Osteoporosis could be explained by less general mechanical stress secondary to stiffness or ankylosis, coupled with reduced local mobility at the ankylosed disco-vertebral units [START_REF] Klingberg | Osteoporosis in ankylosing spondylitis -prevalence, risk factors and methods of assessment[END_REF][START_REF] Klingberg | Bone microarchitecture in ankylosing spondylitis and the association with bone mineral density, fractures, and syndesmophytes[END_REF][START_REF] Pray | Bone Mineral Density and Fracture Risk in Ankylosing Spondylitis: A Meta-Analysis[END_REF][START_REF] Klingberg | Vertebral fractures in ankylosing spondylitis are associated with lower bone mineral density in both central and peripheral skeleton[END_REF][START_REF] Ghozlani | Prevalence and risk factors of osteoporosis and vertebral fractures in patients with ankylosing spondylitis[END_REF][START_REF] Geusens | The prevalence of vertebral fractures in spondyloarthritis: relation to disease characteristics, bone mineral density, syndesmophytes and history of back pain and trauma[END_REF][START_REF] Deminger | Which measuring site in ankylosing spondylitis is best to detect bone loss and what predicts the decline: results from a 5-year prospective study[END_REF].

Dual X-ray absorptiometry (DXA) is considered the gold standard for the diagnosis of osteoporosis. In AS patients, bone mineral density on spine DXA can be overestimated due to the presence of new bone formation (syndesmophytes, bone bridge), especially with respect to anteroposterior incidence [START_REF] Geusens | The prevalence of vertebral fractures in spondyloarthritis: relation to disease characteristics, bone mineral density, syndesmophytes and history of back pain and trauma[END_REF][START_REF] Deminger | Which measuring site in ankylosing spondylitis is best to detect bone loss and what predicts the decline: results from a 5-year prospective study[END_REF][START_REF] Malochet-Guinamand | Prevalence and risk factors of low bone mineral density in spondyloarthritis and prevalence of vertebral fractures[END_REF][START_REF] Wildberger | Impact of lumbar syndesmophyte on bone health as assessed by bone density (BMD) and bone texture (TBS) in men with axial spondyloarthritis[END_REF][START_REF] Magrey | Utility of DXA scanning and risk factors for osteoporosis in ankylosing spondylitis-A prospective study[END_REF].

Thoracic or thoraco-abdomino-pelvic computed tomography scans (CT scans), are often performed in AS patient follow-up to explore extra-articular comorbidities, therapeutic adverse events or before immunosuppressive drug initiation. Spine CT scans can also be performed to evaluate acute back pain in the ankylosed spine. Trabecular bone density evaluation is available on CT scans through the scanographic bone attenuation coefficient (SBAC) of lumbar vertebrae [START_REF] Pickhardt | Opportunistic screening for osteoporosis using abdominal computed tomography scans obtained for other indications[END_REF]. The SBAC permits the estimation of trabecular bone density, avoiding cortical bone, excessive bone formation (syndesmophytes, bone bridge) and vascular calcifications. The same diagnostic performance for osteoporosis screening has been seen for all the lumbar vertebrae, but L1 is the easiest identifiable [START_REF] Pickhardt | Opportunistic screening for osteoporosis using abdominal computed tomography scans obtained for other indications[END_REF]. Recent studies reported that an SBAC-L1 ≤ 145 Hounsfield units (HU) was more reliable than a T-score ≤ -2.5 SD for identifying patients at risk of VF in rheumatoid arthritis (RA) and the general population [START_REF] Pickhardt | Opportunistic screening for osteoporosis using abdominal computed tomography scans obtained for other indications[END_REF][START_REF] Perrier-Cornet | Opportunistic screening for osteoporosis using thoraco-abdomino-pelvic CT-scan assessing the vertebral density in rheumatoid arthritis patients[END_REF]. Furthermore, SBAC-L1 allowed us to identify 50% of sclerodermic patients with bone fragility and was statistically significantly associated with severity criteria of the disease [START_REF] Fauny | Vertebral fracture prevalence and measurement of the scanographic bone attenuation coefficient on CT-scan in patients with systemic sclerosis[END_REF]. CT scans offer an opportunity to explore the trabecular bone density in the AS population and to study the effect of vertebral ankylosis on trabecular bone density.

Our primary objective was to determine the impact of lumbar vertebral ankylosis on trabecular bone fragility through the SBAC measure for each lumbar vertebra according to their ankylosis status: ankylosed versus non-ankylosed. The secondary objective was to determine whether the presence of vertebral ankylosis on mSASSS was associated with a higher risk of reaching an SBAC under the fracture threshold (≤145 HU).

Patients and methods

Population

This descriptive and retrospective study was performed on AS patients followed at the Nancy University Hospital, according to the ASAS 2009 or New York-modified criteria [START_REF] Van Der Linden | Evaluation of diagnostic criteria for ankylosing spondylitis. A proposal for modification of the New York criteria[END_REF][START_REF] The | Evaluation of diagnostic criteria for ankylosing spondylitis: a comparison of the Rome, New York and modified New York criteria in patients with a positive clinical history screening test for ankylosing spondylitis[END_REF]. To be included in this study, patients should have performed radiographs of the pelvis and cervicodorso-lumbar spines as well as a thoracic, TAP or abdomino-pelvic CT scan within a delay of two years between the two exams. The records were selected from the pre-biologic evaluation between 2009 and March 2017. Demographic characteristics (age, gender, smoker status, etc.), clinical data (disease duration, activity score, etc.), biological data (C-reactive protein (CRP)) and treatments were collected from the complete medical record.

AS radiographic assessment

The radiographs were anonymized and read on Osirix software (v6.5.1-64 bits). The frontal pelvic radiographs or radiographs from the lumbo-sacral junction were analyzed for sacroiliitis by two independent senior readers (MF, FV), with adjudication by a third rheumatologist (DL) in case of discordance, according to the modified New York classification [START_REF] Van Der Linden | Evaluation of diagnostic criteria for ankylosing spondylitis. A proposal for modification of the New York criteria[END_REF][START_REF] The | Evaluation of diagnostic criteria for ankylosing spondylitis: a comparison of the Rome, New York and modified New York criteria in patients with a positive clinical history screening test for ankylosing spondylitis[END_REF]. Sagittal spinal radiographs (cervical and lumbar) were evaluated for the presence of squaring, syndesmophytes or ankylosis (bone bridge) for anterior vertebral corners, from C2 inferior to T1 superior and from T12 inferior to S1 superior (score 1 for squaring, 2 for syndesmophytes, 3 for ankylosis), according to mSASSS score, ranging from 0 to 72 [START_REF] Creemers | Assessment of outcome in ankylosing spondylitis: an extended radiographic scoring system[END_REF][START_REF] Wanders | What is the most appropriate radiologic scoring method for ankylosing spondylitis? A comparison of the available methods based on the Outcome Measures in Rheumatology Clinical Trials filter[END_REF]. Structural spinal involvement was defined by an mSASSS ≥ 2, corresponding to the presence of at least one syndesmophyte. The mSASSS was performed by three rheumatologists (MF, FV, DL). The diagnosis of syndesmophyte was considered if two of the three readers agreed. The mSASSS was calculated with the adjudication for each corner, without considering grade 1, in view of its poor reproducibility in previous studies [START_REF] Baraliakos | Progression of radiographic damage in patients with ankylosing spondylitis: defining the central role of syndesmophytes[END_REF][START_REF] Baraliakos | Development of a radiographic scoring tool for ankylosing spondylitis only based on bone formation: Addition of the thoracic spine improves sensitivity to change[END_REF][START_REF] Kim | Do we really need to evaluate entire cervical spines for squaring score in modified stoke ankylosing spondylitis spinal score?[END_REF]. Ankylosed lumbar vertebra was defined by the presence of bone bridges with its two adjacent vertebrae. Partial ankylosed vertebra was defined by the presence of a bone bridge with only one of its adjacent vertebrae. The level and number of ankylosed lumbar vertebrae were then evaluated.

Reproducibility of the measures (mSASSS and sacroiliitis)

Intrareader reproducibility was performed on 30 exams and interreader reproducibility on all exams between the three readers for sacroiliitis and mSASSS.

Densitometric scanographic bone evaluation (scanographic bone attenuation coefficient of L1: SBAC-L1) (Figure 1)

All CT scans were performed at the Nancy University Hospital and were read on OSIRIX software (v6.5.1-64 bits). A total of 73 CT scans were available: 67 TAP CT scans, 3 thoracic CT scans and 3 abdominopelvic CT scans.

The SBAC was measured on axial sections through the pedicles in the bone window for each lumbar vertebra. The largest elliptical region of interest (ROI) was drawn in the trabecular bone and provided the average bone mineral density (in HU). This evaluation was performed by only a reader (MF) because we have previously shown the excellent intra-and interreader reliability of the SBAC-L1 (kappa > 0.9) [START_REF] Fauny | Vertebral fracture prevalence and measurement of the scanographic bone attenuation coefficient on CT-scan in patients with systemic sclerosis[END_REF]. Pickardt et al. showed that CT scans performed similarly from T12 to L5 for the SBAC measure [START_REF] Pickhardt | Opportunistic screening for osteoporosis using abdominal computed tomography scans obtained for other indications[END_REF].

A threshold of 145 HU was used because it identified 96.6% of the patients with VF in a general population, whereas the DXA (with T-score ≤ -2.5 SD) identified only 39% of the patients with VF in the same population [START_REF] Pickhardt | Opportunistic screening for osteoporosis using abdominal computed tomography scans obtained for other indications[END_REF]. This threshold of 145 HU was used to maintain an acceptable balance between sensitivity and specificity in this population, in which a higher risk of bone fragility is suspected.

Ethics approval

All of the data used were obtained from the medical records. No supplementary examination was performed for patients to meet the inclusion criteria. This study is registered to the Information Technology and Freedoms Commission for the University Hospital of Nancy (file number: 2019PI007) and was designed in accordance with the general ethical principles outlined in the Declaration of Helsinki. The protocol of this study was approved by the Information Technology and Freedoms Commission. All patients gave their consent for the use of their medical data during the time period they received medical care at the University Hospital.

Statistical analysis

Both descriptive and comparative analyses were conducted by accounting for the nature and distribution of the variables. Qualitative variables were described as frequencies and percentages; quantitative variables were evaluated with the mean ± standard deviation (SD) or with the median and interquartile range (IQR). For quantitative variables, Student's t test or the Mann-Whitney U test were used. For qualitative variables, the chi-squared test was used, along with Fisher's exact test if necessary. The association between the trabecular bone density of each lumbar vertebra (L1 to L5), measured by the SBAC, and structural spinal involvement was analyzed by logistic binary regression with the odds ratio (OR) and its 95% confidence interval (95% CI). The Kruskal-Wallis test was performed to study the relationship between the SBAC of the lumbar vertebra and the ankylosed status.

To analyze the intrareader and interreader reliabilities, we used Cohen's kappa method. The risk α was established as 0.05. IBM™ SPSS Statistics v23 was used for the data analysis.

Results

Population

Of 1503 spondyloarthritis patients followed between 2009 and 2017 and screened for biologic treatment initiation, 73 patients fulfilled the inclusion criteria. The mean age was 60.3 years (+/-10.7), with a large predominance of men (89%). The disease duration was 24.6 years (+/-24), and the mean BASDAI, BASFI and ASDAS were 7.2 (+/-10.5), 43.5 (+/-23.7) and 3.4 (+/-1.2), respectively. Patients were mainly treated with TNF-blockers (58.9% of the patients), and NSAIDs (67.1% of the patients). A total of 9.6% took corticosteroids (Table 1).

AS radiographic assessment (Figure 2)

All patients presented sacroiliitis on radiographs. The kappa coefficients were 0.71 (95% CI: 0.49-0.89) and 0.82 (95% CI: 0.7-0.93) for the right sacroiliac joint and 0.54 (95% CI: 0.30-0.75) and 0.69 (95% CI: 0.55-0.82) for the left sacroiliac joint, for inter-and intrareader reproducibility, respectively. Concerning spinal involvement, the cervical spine assessment was technically limited by spine deformity and ankylosis. Indeed, an optimal assessment of this segment was obtained for 79.5% of the patients. Fourteen patients (19.2%) had no spinal lesions. The mean mSASSS was 20.7 (±21.3), and 60 patients (82.2%) presented an mSASSS ≥ 2. Thirty-seven patients (50.7%) had at least one bone bridge on the spine, and six patients (8.2%) had a bamboo spine. For the lumbar segment, 17 (23.3%) and 13 (17.8%) patients had complete or partial vertebral ankylosis, respectively. The intra-and interreader reproducibility for spinal structural involvement was excellent. For the mSASSS, the kappa coefficients were 0.97 (95% CI: 0.93-0.98) and 0.97 (95% CI: 0.95-0.98) for intra-and interreader, respectively. Adjudication resulted in the reclassification of three patients in the mSASSS-group.

Densitometric scanographic bone evaluation of lumbar vertebrae (SBAC)

Two patients had no available evaluation at L2, L3, L4 and L5 levels, one at L5 vertebra and one at L1 vertebra, due to severe spine deformation. All four patients were in the mSASSS+ group. The mean SBAC was 141.1 HU (±45) for L1, 130.4 HU (±52.9) for L2, 127.9 HU (±52.2) for L3, 124.3 HU (±58.2) for L4 and 135 HU (±56.7) for L5. There was not significant difference between the mean SBAC for each lumbar vertebra between patients with or without vertebral fractures (p non-significant); 7 of the 9 patients with VF have a SBAC for most lumbar vertebrae under the fracture threshold. We found no statistically significant association between SBAC from L1 to L5 and mSASSS status. However, there was a significant difference in SBAC values between patients with bone bridges versus patients without bone bridges (p < 0.01).

Compared to patients without bone bridges, patients with bone bridges of the lumbar vertebrae more frequently presented an SBAC under the fracture threshold of 145 HU (for each lumbar vertebra) (p < 0.05). Thirty-four patients (91.9%) with bone bridges on mSASSS presented an SBAC under the fracture threshold for one of the five lumbar vertebrae (Table 2). We did not find any correlation between the bone mineral density or T-score and the SBAC for each lumbar vertebra.

Association between SBAC of the lumbar vertebra and ankylosed status

The mean lumbar SBAC progressively decreased according to the degree of lumbar vertebra ankylosis from 150.5 HU (±44.9) in non-ankylosed vertebrae to 137.9 HU (±44.

Discussion

This study is the first to show that severe spinal involvement, such as bone bridges in cervical or lumbar spine segments on mSASSS is associated with lower lumbar SBAC values. More precisely, local partial and complete ankylosis of any lumbar vertebra were associated with lower values of lumbar SBAC (L1 to L5), suggesting deterioration of the trabecular bone density in AS. In our AS patients, nearly 60% presented bone fragility with an SBAC≤145 HU.

The population of this study was similar to other AS populations, such as the OASIS or GESPIC cohorts [START_REF] Protopopov | Relevance of structural damage in the sacroiliac joints for the functional status and spinal mobility in patients with axial spondyloarthritis: results from the German Spondyloarthritis Inception Cohort[END_REF]. Our patients were more frequently treated with TNF-blockers (58.9%), but NSAID intake was similar to that of other cohorts (approximately 67%). In the OASIS cohort, 19% of the patients were mSASSS- [START_REF] Ramiro | Higher disease activity leads to more structural damage in the spine in ankylosing spondylitis: 12-year longitudinal data from the OASIS cohort[END_REF], which was similar to our population (with 18.7% of mSASSS-). However, the mean mSASSS was higher (20.7) in our population than in the OASIS cohort (10.8), probably because our study included older patients with a longer disease duration [START_REF] Ramiro | Higher disease activity leads to more structural damage in the spine in ankylosing spondylitis: 12-year longitudinal data from the OASIS cohort[END_REF]. We have also shown that osteoporosis screening is not automatically performed in routine care, despite the recommendations of French Society for Rheumatology (SFR) about the screening and management of comorbidities, including osteoporosis in spondyloarthritis patients [START_REF] Wendling | update of French Society for Rheumatology (SFR) recommendations about the everyday management of patients with spondyloarthritis[END_REF]. In our patients, only 22% had undergone DXA, whereas more than 72% of them had at least one clinical or biological risk factor for osteoporosis. One of the limits of this study was the low number of DXA performed in our population. We did not find any correlation between the bone mineral density or T-score and the SBAC for each lumbar vertebra. Due to the low number of DXA, we are not able to conclude about this correlation. All of the patients with osteoporosis or low BMD (osteopenia) on DXA have a SBAC for each lumbar vertebra under the fracture threshold of 145 HU. Some studies showed an interest in the use of trabecular bone score (TBS), not influenced by syndesmophytes and correlated to lumbar BMD [START_REF] Wildberger | Impact of lumbar syndesmophyte on bone health as assessed by bone density (BMD) and bone texture (TBS) in men with axial spondyloarthritis[END_REF][START_REF] Boussoualim | Evaluation of bone quality with trabecular bone score in active spondyloarthritis[END_REF].

The mSASSS evaluates spinal damage in AS from the antero-inferior corner of C2 to the anterosuperior corner of T1 and from the antero-inferior corner of T12 to the antero-superior corner of S1 [START_REF] Creemers | Assessment of outcome in ankylosing spondylitis: an extended radiographic scoring system[END_REF][START_REF] Wanders | What is the most appropriate radiologic scoring method for ankylosing spondylitis? A comparison of the available methods based on the Outcome Measures in Rheumatology Clinical Trials filter[END_REF]. In our study, the cervical evaluation was technically limited due to severe ankylosis and deformation; however, 79.5% of the patients had an optimal assessment of the cervical spine. As previously reported, we obtained very good reproducibility for mSASSS with a kappa coefficient of 0.97 for intra-and interreaders [START_REF] Claudepierre | Reliability of mSASSS scoring in everyday practice in DESIR-cohort study centres: cross-sectional study of agreement with trained readers[END_REF][START_REF] Biagioni | Reliability of radiographic scoring methods in axial psoriatic arthritis[END_REF][START_REF] Salaffi | Radiological scoring methods for ankylosing spondylitis: a comparison between the Bath Ankylosing Spondylitis Radiology Index and the modified Stoke Ankylosing Spondylitis Spine Score[END_REF]. Only anterior vertebral corners were evaluated with the mSASSS. The posterior segment of the spine was not evaluated; however, various studies have shown that 20 to 50% of AS patients had facet joint ankylosis in the cervical and/or lumbar spine associated with anterior spinal involvement [START_REF] Maas | Incorporating assessment of the cervical facet joints in the modified Stoke ankylosing spondylitis spine score is of additional value in the evaluation of spinal radiographic outcome in ankylosing spondylitis[END_REF][START_REF] Maas | Radiographic damage and progression of the cervical spine in ankylosing spondylitis patients treated with TNF-α inhibitors: Facet joints vs. vertebral bodies[END_REF]. On radiographs, posterior ankylosis assessment is difficult due to the two-dimensional projection of the lumbar facet joints. In our study, we were not able to compare patients with isolated cervical or lumbar bone bridges due to the small number in each group (only 7 patients with isolated cervical bone bridges and 9 with isolated lumbar bone bridges).

We did not show any association between SBAC and mSASSS. This could be explained by the fact that mSASSS is a composite score assessing different structural damages: squaring, erosion, syndesmophyte and bone bridge, reflecting different levels of structural severity. Squaring and erosion were excluded from our analysis due to poor intra-and interreader reliabilities, as previously mentioned in the literature [START_REF] Baraliakos | Progression of radiographic damage in patients with ankylosing spondylitis: defining the central role of syndesmophytes[END_REF][START_REF] Baraliakos | Development of a radiographic scoring tool for ankylosing spondylitis only based on bone formation: Addition of the thoracic spine improves sensitivity to change[END_REF][START_REF] Kim | Do we really need to evaluate entire cervical spines for squaring score in modified stoke ankylosing spondylitis spinal score?[END_REF]. In this study, we showed that SBAC significantly decreased when vertebrae presented at least one bone bridge in the cervical or lumbar spine. Moreover, when we studied the local interaction between the degree of lumbar vertebra mobility and the lumbar SBAC (L1-L5) values, AS patients with at least one bone bridge presented a statistically significant risk of having an SBAC under the fracture threshold compared to AS patients with no ankylosis or only syndesmophytes without bone bridges. Our study did not respond to a potential deleterious effect of syndesmophytes on SBAC since the number of patients without spinal lesions was too low (n=14).

The use of SBAC is interesting in AS patients because this technique permits direct evaluation of the trabecular bone density, avoiding cortical artifacts, such as osteophytes and syndesmophytes, as well as other artifacts (vascular calcifications) frequently observed in elderly AS patients. Furthermore, the SBAC measure is very reproducible [START_REF] Fauny | Vertebral fracture prevalence and measurement of the scanographic bone attenuation coefficient on CT-scan in patients with systemic sclerosis[END_REF][START_REF] Fauny | Study of vertebral fracture Scanographic Bone Attenuation Coefficient in rheumatoid arthritis and ankylosing spondylitis vs. controls[END_REF] and easy to perform. This method could be used in daily practice to consider initial opportunistic bone screening on CT scans already performed during AS patients' follow-up (treatment side effects, comorbidity evaluation, emergency, etc.). The SBAC was not available for four patients. Among these four patients, SBAC of L2, L3, L4 and L5 levels were not available for two patients, SBAC of L1 for one patient and SBAC of L5 for another patient because spinal deformation made the measurement impossible in the region of interest. For these patients, the SBAC evaluation was performed on the available vertebrae. All of these patients were in the mSASSS+ group.

This study is the first to evaluate the consequence of spinal ankylosis on trabecular bone density using opportunistic CT scans in AS patients. Complete ankylosis of the lumbar vertebra and partial ankylosis of L3 and L4 were associated with lower values of SBAC, suggesting deterioration of the trabecular bone structure in ankylosed vertebrae. The deterioration of trabecular bone attenuation could be explained by a loss of mechanical strain related to ankylosis (partial or total). Indeed, we know that mechanical strain on the bone increases its density, as in physical activity [START_REF] Moreira | Physical exercise and osteoporosis: effects of different types of exercises on bone and physical function of postmenopausal women[END_REF][START_REF] Birge | The role of exercise in preventing osteoporosis[END_REF][START_REF] Kannus | Physical loading, exercise, and bone[END_REF][START_REF] Pfeifer | Musculoskeletal rehabilitation in osteoporosis: a review[END_REF]. Moreover, bone ankylosis may also reflect patients with a more severe disease and more biological inflammation [START_REF] Aouad | Structural progression in axial spondyloarthritis[END_REF], resulting in less physical activity and mobility. Local or systemic inflammation could also increase spinal bone loss. This hypothesis may be suggested in our population since patients with bone bridges had more severe disease activity measured by the ASDAS. This study raises questions about the sensitivity to changes in SBAC over time and the effect of biologic treatments and antiosteoporotic drugs on trabecular bone mineralization. Further longitudinal studies are needed to answer these questions.

In conclusion, we showed that CT scans is an alternative exam for screening bone fragility (through the SBAC of each lumbar vertebrae) in AS patients. The SBAC measure is highly reproducible and can be easily performed in daily practice on opportunistic CT scans (performed for any other indication). We showed that patients with bone bridges on mSASSS (cervical or lumbar segments) and patients with lumbar ankylosis (total or partial) had a higher risk of presenting an SBAC under the fracture threshold. Patients with spinal ankylosis should be screened for osteoporosis, regardless of other osteoporosis risk factors, using opportunistic CT scans and DXA (the gold standard). The bone loss kinetics and structural progression should also be investigated in further longitudinal studies.

Figure 1 :Figure 2 :

 12 Figure 1: SBAC-L1 measurement on normal spines (A), spines with syndesmophytes (B) and ankylosed spines (C) (sagittal lumbar spine radiographs and CT scan in the bone window and axial CT scan section of L1 with SBAC-L1 measured on the ROI)

  The mean duration between radiographs and CT scans was 75.6 days (±166.4). Concerning osteoporosis screening, 53 patients (72.6%) had at least one clinical or biological risk factor for osteoporosis. Sixteen (21.9%) underwent spine DXA during their follow-up. Among them, five patients had osteoporosis on DXA when the 3 sites (spine, total hip and femoral neck) were taken into account. VFs were depicted in 9 patients on radiographs. DXA images were available for three of them. One patient presented osteoporosis (T-score: -2.7 SD),

another patient was classified with osteopenia (T-score: -1.9 SD), and the last one had a normal T-score (+0.6 SD). Among the 73 patients, 11 (15.1%), 16 (21.9%) and 12 (16.4%) took calcium tablets, vitamin D supplementation and specific anti-osteoporotic drugs, respectively. The anti-osteoporotic treatments were bisphosphonates (10/12: zoledronic acid, alendronate, risedronate, pamidronate), denosumab (1/12) and strontium renalate (1/12).

Table 1 :

 1 : ankylosing spondylitis, ASDAS: Ankylosing Spondylitis Disease Activity Score, BASDAI: Bath ankylosing spondylitis disease activity index, BASFI: Bath ankylosing spondylitis functional index, CRP: C-reactive protein, mSASSS: modified Stoke Ankylosing Spondylitis Spinal Score, NS: non-significant, NSAIDs: Non-steroidal anti-inflammatory drugs, TNF: tumor necrosis factor Age in the mSASSS+ group and disease duration are the median (IQR), whereas age in the general population, BASDAI, BASFI, and ASDAS are the mean (± SD). In the comparisons, Student's t-test or the Mann-Whitney U test were used. For qualitative variables, the chi-square test with, if necessary, the exact calculation of Fisher was used. The risk α was established as 0.05 Demographical and clinical characteristics of the 73 AS patients according to the mSASSS (mSASSS +/mSASSS -) and bone bridge status

				Total population	mSASSS + mSASSS -	p	Bone bridge +	Bone bridge -	p
				n=73	n=60	n=13		n=37	n=23
	Demographic and clinical	Age (n=73)	60.3 (60)	62 (61)	52.7 (57)	NS	62 (57.5-73) 59 (53-67)	NS
			Men (n=73)	65 (89)	53 (88.3)	12 (92.3)	NS	34 (91.9)	19 (82.6)	NS
			Smoker (n=48)	31 (42.5)	25 (41.7)	6 (46.2)	NS	17 (45.9)	8 (34.8)	NS
		Alcohol intake (n=28)	4 (5.5)	4 (6.7)	0 (0)	NS	4 (10.8)	0 (0)	NS
	AS characteristics	Disease duration (n=71) 24 (12-34)	25 (13-35)	15.5 (12-25.3)	NS	27 (13.8-35.8)	24 (11-35)	NS
			BASFI (n=58) 43.5 (23.7)	44.5 (22.2)	41.2 (30.8)	NS	46.7 (22.2)	42 (22.4)	NS
			BASDAI (n=62) 7.2 (10.5)	7.1 (10.5)	8.4 (11.3)	NS	7.2 (12.7)	6.9 (7.3)	NS
			ASDAS (n=21)	3.4 (1.2)	3.4 (0.9)	3 (2.6)	NS	3.7 (1)	2.81 (0.6)	<0.05
			HLA B27 (n=69)	46 (63)	39 (65)	7 (53.8)	NS	26 (70.3)	13 (56.5)	NS
			CRP (n=69)	31 (42.5)	28 (46.7)	3 (23.1)	NS	18 (48.6)	10 (43.5)	NS
			mSASSS (n=73) 20.7 (21.2)	25.1 (20.8)	0 (0)	0.0001	37.2 (17.8)	5.9 (4.2)	0.0001
	Mean grade of left sacroiliitis (n=73)	3.3 (0.8)
	Mean grade of right sacroiliitis (n=73)	3.4 (0.6)
	Treatments	Corticosteroids (n=49)	7 (9.6)	5 (8.3)	2 (15.4)	NS	3 (8.1)	2 (8.7)	NS
			NSAIDs (n=51)	49 (67.1)	44 (73.3)	5 (38.5)	NS	30 (81.1)	14 (60.9)	0.05
	PPI Proton Pump Inhibitor (n=15)	15 (20.5)	14 (23.3)	1 (7.7)	/	9 (24.3)	5 (21.7)	/
			TNF-blockers (n=43)	43 (58.9)	36 (60)	7 (53.8)	/	23 (62.2)	13 (56.5)	/
	Osteoporosis		≥1 risk factor(n=73)	53 (72.6)	46 (76.7)	7 (53.8)	NS	29 (78.4)	17 (73.9)	NS
	Vitamin D and/or calcium (n=73)	18 (24.6)	14 (23.3)	4 (30.8)	NS	9 (24.3)	5 (21.7)	NS
	Anti-osteoporotic drug (n=72)	12 (16.4)	9 (15)	3 (23.1)	NS	8 (21.6)	1 (4.3)	NS

AS

Table 2 :

 2 Scanographic bone evaluation of the lumbar spine according to the mSASSS and bone bridge status. : Hounsfield unit; mSASSS: modified Stoke Ankylosing Spondylitis Spinal Score; NS: non-significant; SBAC: scanographic bone attenuation coefficient SBAC values expressed as the mean (SD), threshold ≤ 145 HU in n (%)

		Total	mSASSS +	mSASSS -		Bone bridge	Bone bridge
		population	N=60	N=13	P	+ N=37	-	P
		N=73					N=23
	SBAC-L1	141.1 (45)	138.1 (44.8) 154.8 (44.9)	NS	123.9 (41.1) 160.4 (41.9)	< 0.01
	≤ 145 HU	42 (57.5)	36 (60)	6 (46.2)	NS	27 (73)	9 (39.1)	< 0.01
	SBAC-L2	130.4 (52.9) 127.2 (55.1) 144.5 (41.1)	NS	111.5 (48.8) 149.8 (50.5)	< 0.01
	≤ 145 HU	42 (57.5)	35 (58.3)	7 (53.8)	NS	26 (70.3)	9 (39.1)	< 0.05
	SBAC-L3	127.9 (52.2) 124.9 (54.4) 141.3 (39.9)	NS	105.2 (53.1) 151.3 (40.1)	< 0.001
	≤ 145 HU	47 (64.4)	39 (65)	8 (61.5)	NS	29 (78.4)	10 (43.5)	< 0.01
	SBAC-L4	124.3 (58.2) 119.7 (61.2)	145 (38.8)	NS	99.7 (61.2)	152.5 (41.1)	< 0.001
	≤ 145 HU	45 (61.6)	37 (61.7)	8 (61.5)	NS	27 (73)	10 (43.5)	< 0.05
	SBAC-L5	135 (56.7)	132.1 (59.4) 147.8 (42.3)	NS	106.6 (51.7) 163.4 (46.8)	< 0.001
	≤ 145 HU	40 (54.8)	34 (56.7)	6 (46.2)	NS	29 (78.4)	5 (21.7)	< 0.001

HU

Table 3 :

 3 Number and percentage of lumbar vertebrae with an SBAC≤145 HU according to structural involvement (syndesmophytes, partial or total ankylosis) Hounsfield units. L1/L2/L3/L4/L5: first/second/third/fourth/fifth lumbar vertebra. NS: non-significant. SBAC: scanographic bone attenuation coefficient. SD: standard deviation

		SBAC L1		SBAC L2		SBAC L3		SBAC L4		SBAC L5
		≤145 HU	p	≤145 HU	p	≤145 HU	p	≤145 HU	p	≤145 HU	p
		n (%)		n (%)		n (%)		n (%)		n (%)
	Without structural impairment	20 (47.6)		21 (47.7)		20 (54.1)		19 (47.5)		20 (43.5)
	Syndesmophyte(s)	5 (71.4)	NS	5 (62.5)	NS	8 (53.3)	NS	9 (69.2)	NS	5 (50)	0.05
	Partial ankylosis	8 (66.7)	NS	7 (63.6)	NS	9 (90)	<0.01	8 (80)	<0.01	7 (87.5)	NS
	Total ankylosis	9 (75)	<0.05	10 (100)	<0.01 10 (90.9) <0.001	9 (90)	<0.001	8 (88.9)	<0.01
	HU:

Acknowledgments

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Funding

None

Availability of data and material

Data are available upon reasonable request.

Declarations Ethics approval and consent to participate

All of the data used were obtained from medical records. No examinations were performed for patients to meet the inclusion criteria. This study is registered to the Information Technology and Freedoms Commission for the University Hospital of Nancy (file number: 2019PI007) and was designed in accordance with the general ethical principles outlined in the Declaration of Helsinki. The protocol of this study was approved by the Information Technology and Freedoms Commission for the University Hospital of Nancy. All patients gave consent for the use of their medical data during the time period they received medical care at the University Hospital.

Competing interest

All of the authors declare that they have no conflicts of interest.

Authors' contributions

All authors meet the journal's criteria for authorship and have made substantial contributions to all of the following: (1) conception and design of the study, acquisition of data, or analysis and interpretation of data; (2) drafting the article or critical revision for important intellectual content; and (3) final approval of the version to be submitted.

Submission declaration

This work has not been published previously, it is not under consideration for publication elsewhere, this publication is approved by all authors and tacitly or explicitly by the responsible authorities where the work was carried out, and, if accepted, it will not be published elsewhere in the same form, in English or in other language, including electronically without the written consent of the copyright-holder.


Pierre-Olivier Brissaud
email: pierre-olivier.brissaud@thalesgroup.com

Jérôme François
email: jerome.francois@loria.fr

Isabelle Chrisment
email: isabelle.chrisment@loria.fr

Thibault Cholez
email: thibault.cholez@loria.fr

Olivier Bettan
email: olivier.bettan@thalesgroup.com

Encrypted HTTP/2 Traffic Monitoring: Standing the Test of Time and Space

Keywords: Traffic Analysis, HTTP2, HTTPS monitoring, Privacy, TLS, Encrypted Traffic

Encrypted HTTP/2 (h2) has been worldwide adopted since its official release in 2015. The major services over Internet use it to protect the user privacy against traffic interception. However, under the guise of privacy, one can hide the abnormal or even illegal use of a service. It has been demonstrated that machine learning algorithms combined with a proper set of features are still able to identify the incriminated traffic even when it is encrypted with h2. However, it can also be used to track normal service use and so endanger privacy of Internet users. Independently of the final objective, it is extremely important for a security practitioner to understand the efficiency of such a technique and its limit. No existing research has been achieved to assess how generic is it to be directly applicable to any service or website and how long an acceptable accuracy can be maintained.

This paper addresses these challenges by defining an experimental methodology applied on more than 3000 different websites and also over four months continuously. The results highlight that an off-the-shelf machine-learning method to classify h2 traffic is applicable to many websites but a weekly training may be needed to keep the model accurate.

I. INTRODUCTION

HTTPS is nowadays the default protocol for 59.4% 1 of websites in the Internet mostly due to the users' expectation regarding privacy and security after the scandals revealed over the last decade and the increase of data leakage and identity theft impacting individuals.

Among the websites served by HTTPS, according to W3Techs 2 , 41.4% run over HTTP/2 and this number has increased by more than 30% over one year. Moreover, 77 websites of the top 100 websites 3 use HTTP/2, which makes HTTP/2 the future protocol for encrypted web traffic and that will replace HTTP/1.1 over TLS. In this paper, we thus consider the analysis of HTTP2 traffic over TLS also called h2.

Thanks to large deployment of HTTPS and under the guise of privacy, malicious users can also unfortunately hide abnormal or even illegal use of a service. In our prior work [START_REF] Brissaud | Transparent and service-agnostic monitoring of encrypted web traffic[END_REF], we have demonstrated that machine learning algorithms can 1 https://w3techs.com/technologies/details/ce-httpsdefault 2 https://w3techs.com/technologies/details/ce-http2/all/all 3 results from our own tests help fight against this kind of malicious behavior. We have proposed a classification approach which, combined with a proper set of features, is able to identify the incriminated HTTP/2 traffic even when encrypted. It can be used to detect specific forbidden user actions such as a keyword search performed on some services. Quite rightly, we cannot exclude that our technique can be exploited by malevolent users to track user activities and so it indirectly increase the privacy risks. In all cases, the scope of the applicability of the technique and its limit are important to understand. On one hand, a security practitioner may rely on it to trigger alerts and need to know how accurate the technique is in different scenarios. On the other hand, a user expecting privacy would need to know how much is compromised, that may possibly lead to new protocol refinements.

The robustness of a classification model is always questionable and is mainly impacted in our case by two main factors. First, the traffic generated for a same request can change over time (how and what data are retrieved) and thus impacting the validity of models learned. Second, the generalization of a model can only be validated thanks to a large number of different instances to be tested. We refer to these two problems as test-of-time and test-of-space respectively. A practical limitation to address these issues is the availability of extensive datasets (in time and space).

This paper thus extends our previous work [START_REF] Brissaud | Transparent and service-agnostic monitoring of encrypted web traffic[END_REF] by investigating the two questions mentioned above. Besides, as we consider the de facto standard for web nowadays, our reported results aim at evaluating the viability of our technique on-theshelf, i.e. with no specific parameter customization, and at large scale. Therefore it assesses the viability of monitoring encrypted HTTP2.0 traffic.

To summarize our main contributions are three-fold:

• Define a test-of-time and test-of-space methodology for HTTPS classification including the specification of a crawling campaign to collect relevant datasets. The collected dataset can be accessible on-demand (its size is too large to be publicly accessible on our servers). They consists of the pcap and HTML files but also screenshots of each accessed website. • Evaluate our H2 classifier over around four months (testof-time) and the impact of a regular re-training.

• Evaluate our H2 classifier on more than 3000 websites (test-of-space). The paper is organized as follows. Section II presents related works and background, in particular our classification technique omitting in-depth details available in [START_REF] Brissaud | Transparent and service-agnostic monitoring of encrypted web traffic[END_REF]. Our test methodology is given in Section III along with dataset descriptions. Section IV and Section V respectively give the evaluation over time and over different services. Section VI concludes the paper.

II. RELATED WORK AND BACKGROUND

A. Related Work

Privacy is not guaranteed by encryption especially with HTTPS. Already in 1996, Wagner et al. [START_REF] Wagner | Analysis of the ssl 3.0 protocol[END_REF] found vulnerabilities for SSL protocol against traffic analysis. The encrypted traffic analysis topic is a wide research area and has been approached from different angles.

Velan et al. [START_REF] Velan | A survey of methods for encrypted traffic classification and analysis[END_REF] present in their survey classification methods which allow detecting applications or protocols inside a encrypted flow. On their side, Shbair et al. [START_REF] Shbair | A multilevel framework to identify https services[END_REF] proposed a framework for detecting web services over HTTPS.

Pescape et al. [START_REF] Pescape | Anonymity services tor, i2p, jondonym: Classifying in the dark (web)[END_REF] showed they could distinguish the anonymized networks (Tor, i2P, JonDonym), from each other and to classify types of applications. Website fingerprinting on anonymous networks such as Tor has been extensively studied in literature like in [START_REF] Panchenko | Website fingerprinting in onion routing based anonymization networks[END_REF]- [START_REF] Hayes | k-fingerprinting: A robust scalable website fingerprinting technique[END_REF]. These authors used different clustering methods and features, mostly based on the packet size distribution, for their fingerprints.

Another part of the research activities has also concerned mobile traffic analysis for identifying the services installed on a mobile device. AppScanner [START_REF] Taylor | Robust smartphone app identification via encrypted network traffic analysis[END_REF] profiles 110 mobile applications with an accuracy between 73% and 96%. In this area, particular activities or behaviors might be detected within the flow [START_REF] Conti | Analyzing android encrypted network traffic to identify user actions[END_REF]- [START_REF] Saltaformaggio | Eavesdropping on fine-grained user activities within smartphone apps over encrypted network traffic[END_REF]. On top of that, specific analysis for some applications have been put forward, e.g. for Netflix [START_REF] Reed | Identifying https-protected netflix videos in real-time[END_REF] or Skype [START_REF] Wright | Spot me if you can: Uncovering spoken phrases in encrypted voip conversations[END_REF], [START_REF] Wright | Language identification of encrypted voip traffic: Alejandra y roberto or alice and bob?[END_REF]. These works are built on the awareness of the encoding techniques for video and audio and their impact on the encrypted domain.

Finally we designed an HTTP/2 traffic classifier (H2CLassifier) [START_REF] Brissaud | Transparent and service-agnostic monitoring of encrypted web traffic[END_REF] in order to detect specific behaviors inside a predefined service over HTTP/2. To the best of our knowledge, H2Classifier is the first solution able to identify user actions on a service protected by HTTPS (with HTTP/2). In this paper we evaluate the H2Classifier over time and space.

For the case of the space it is hard to compare with the other works because we don't detect the same concept. But for the time evaluation it can be tested for different works. However, this question is not usually discuss in the paper presented before. In the work of Juarez [START_REF] Juarez | A critical evaluation of website fingerprinting attacks[END_REF] which discuss about some limitations of the traffic analysis the topic of efficiency over the time is handle. Especially they show that in their case the accuracy drop by 50% in under one month for their classification. For this reason we deal with this question in this paper to evaluate the impact and find solutions.

B. Traffic Classification: H2Classifier

H2Classifier's goal is to passively and transparently infer actions of encrypted web services, more especially HTTPS using HTTP2. In this work, user actions are actually requests performed by users such as making use of search fields.

Our method leverages the Random Forest (RF) algorithm [START_REF] Breiman | Random forests[END_REF]. For each individual service s in the set S of all services to be monitored, our method will automatically learn a classification model related to particular actions defined a priori. As an example S could be defined as {Google Images, Amazon, Instagram}.

In a nutshell, our technique relies on collected traffic associated with particular actions on a given service s. This traffic contains multiple instances of the same action in order to apply supervised learning with Random Forest and then be able to classify a new unknown instance of collected traffic. Therefore, it is necessary to establish the list of monitored user actions in order to generate and collect the associated traffic, which also prevent us to detect any other (private) actions but the monitored ones.

The main issue of all machine learning-based techniques is the definition of features specific to the addressed problem. In our case, we have to identify the user request and so the content of a request that is contained in TLS records. Due to encryption, only meta-data features can be extracted and helpful [START_REF] Brissaud | Transparent and service-agnostic monitoring of encrypted web traffic[END_REF]:

• Connection statistics (3 features): Number of incoming and outgoing TLS records, total number of bytes exchanged at the TLS layer. • Burst statistics (10 features): the minimum, maximum, standard deviation, mean and median of (1) number of bytes sent by the server between two client packets and

(2) number of records send by the server between 20 records.

• Number of distinct sizes of incoming and outgoing TLS records (two features). • Top 20 most representative sizes of TLS records (2×20):

the 20 TLS record sizes with the lowest frequencies, from incoming and outgoing packets. • Size distribution: all TLS record size frequencies for both the request and response. It is a theoretically highly dimensional vector (up to 2×18,432 by definition) but in practice, most values are never observed and have no impact on the result and are so discarded from the feature list.

The rational behind this feature set is to capture indirectly some estimates about the size of loaded objects after a web page is requested, as for example images. These object sizes are actually very discriminating due to their high variance in size and so helpful to fingerprint the content of a page, especially if it is composed of many objects. As an example, thumbnails corresponding to products in an online store are usually related to the user search among other criteria such as her own history and profile. Finally, the feature set is independent of TLS version (1.2 or 1.3) making H2Classifier compatible with last standards.

III. METHODOLOGY AND DATASETS

A. Methodology

In our previous paper [START_REF] Brissaud | Transparent and service-agnostic monitoring of encrypted web traffic[END_REF], the results demonstrate a high accuracy of our technique (> 94%) for five major services in the Internet namely Google Search, Google Images, Google Maps, Amazon and Instagram. However, the robustness of our technique might be questionable, especially to verify how much accurate our technique applies to other web services (test-of-space) and how long the learned RF models can be valid (test-of-time). Indeed, contents in the Internet can be frequently updated; a database of meaningful models (for each service to be monitored) must thus be maintained. This would require to collect traffic representative of the new contents and retraining the models. As an additional question, knowing the root causes of the classifier errors is important in order to assess if it is due to a content change or a modification of the technical environment (protocol version, server software, application libraries used, etc.).

To address these objectives without a bias, the configuration of parameters strictly follows recommendations made in our previous work. For RF classification hyper-parameters, the number of trees and the maximum depth of trees are set to 400 and 50 respectively. Only the 300 most representative features are used based on the Gini importance (identify the features which mostly reduce the impurity based on the Gini criterion).

For evaluating the robustness of the classifier over time and space, two datasets have been constituted and are detailed in Sections III-C and III-D. We make them available on demand 4to allow reproducible research and benefit to other works. Each dataset includes all the packets of every captured web page as well as the related screenshot, the HTML export and some additional information.

B. Crawling tool

Our data crawler is implemented in Python based on Selenium 5 coupled with geckodriver6 in order to control a Firefox web browser. This setup allows the automated loading of HTTPS pages while the tcpdump 7 program is used in parallel to capture all the packets. In the following of the paper we refer to packets captured from the loading of a page on a service as a trace.

Our tool is able to automaticaly detect and fill-out a search bar and thus to automate the search of keywords on random websites. During the first loading of a web page, it extracts the HTML of the page and looks for an input tag of type "text", "search" or "url". Most of the time, a search bar has a special

C. Temporal Dataset

For the temporal evaluation, we actively crawled 500 keywords of four major services (Amazon, Instagram, Google and Google Images). Four individual traces are collected for each couple of keyword and service per day as summarized in Table I. In total, 8000 daily traces are stored during a 121days period (about 4 months) except a few of maintenance days operated in our lab infrastructure (June 21th to 26th, July 21th and 22th).

For Amazon, Google and Google Images, the list of keywords is identical and was built from the dataset of J. McAuley [START_REF] Mcauley | Image-based recommendations on styles and substitutes[END_REF]. It contains a list of exact Amazon product names. From this list we extract the most frequently 3-grams of words and use them as keywords. The benefit of using those keywords is that we get a lot of related results. Regarding Instagram, the top 500 account names sorted by followers have been used. The complete list of the exact keywords used for the different services is not detailed here but is available on the web page of our dataset.

D. Service-wide Dataset 1) Service selection:

The service wide-dataset is built by requesting many popular websites. Starting from an initial list of the most visited US websites 10 , those with no h2 (HTTPS with HTTP/2) or search bars have been excluded. In addition, a few websites cannot be accessed or generated errors (e.g. timeout) in the capture. As a result, our dataset is composed of 3096 websites crawled between July 15th and August 14th 2019. The crawling campaign has been spread over time only for scalability reason by opposition to the previous dataset where time is carefully taken into consideration. Similarly, parallelization was performed thanks to multiple IP addresses. For ethical research reason and avoid flooding the targeted websites, queries have been alternated within a window of 20 websites. Hence, a trace was collected every 15 seconds and a website was queried every 5 minutes in average.

2) Keyword selection: For each visited website, we use 20 keywords and capture 20 traces for each of them (so a total of 400 traces per website) as stated in Table II. Keywords must be carefully selected due to the variety in website content. Irrelevant keywords will lead to meaningless webpages, as for example with no search results whereas our goal is to identify pages associated with particular user actions and so representative of realistic requests. Thus, we defined a distinct keyword list for each type or category of websites based on the FortiGuard classification 11 .

As highlighted in Figure 1, the distribution of the categories of the 3096 tested websites is not balanced. We focused so on defining specific keywords by hand for the 12 most represented categories whereas random words from Oxford English dictionary have been considered in other cases.

IV. TEST-OF-TIME

The goal of this experimentation is to evaluate the efficiency of H2Classifier over time and if its accuracy can be maintained by regularly updating the model based on the dataset described in Section III-C.

A. Fixed learned models

The traces collected during the 15 first days are reserved for the training stage of the RF models (one for each service).

Results are presented in Figure 2. The grey areas correspond to maintenance days where no data has been collected. As a recall, 500 keywords are tested four times a day for each service. The value thus represents the mean accuracy, i.e. the proportion of traces that has been correctly classified. The global trend is a slow decrease of the accuracy over the time. This observation acknowledged our intuition as the content of the different pages is updated over time and, obviously, the more time spent, the more content might change.

Amazon and Google Images services have a smooth moderate decrease over time. The definition of an acceptable 11 https://fortiguard.com/webfilter 0 5 -2 6 -2 0 1 9 0 6 -0 2 -2 0 1 9 0 6 -0 9 -2 0 1 9 0 6 -1 6 -2 0 1 9 0 6 -2 3 -2 0 1 9 0 6 -3 0 -2 0 1 9 0 7 -0 7 -2 0 1 9 0 7 -1 4 -2 0 1 9 0 7 -2 1 -2 0 1 9 0 7 -2 8 -2 0 1 9 0 8 -0 4 -2 0 1 9 0 8 -1 1 -2 0 1 9 0 8 -1 8 -2 0 1 9 0 8 -2 5 -2 0 1 9 0 9 -0 1 -2 0 1 9 0 9 -0 8 -2 0 1 9 Dates accuracy depends on the type of application. For example, if the objective is to block some users accessing illicit contents, high accuracy is desired. In case the objective is to extract a general trend of user activities, lower accuracy can be acceptable. Regarding Google search engine, we observe a fast decrease in the first days till a stable value around 0.2. Therefore, in addition to a higher difficulty of monitoring this service as highlighted in [START_REF] Brissaud | Transparent and service-agnostic monitoring of encrypted web traffic[END_REF], the content seems to be more dynamic and invalidates the test of time for this particular case.

Instagram follows the same decreasing trend profile as Amazon and Google Images except three noticeable drops. We assume here two hypotheses: (1) the content of pages related to a monitored keyword has been modified or (2) the manner to send content has changed, for example using new library/software including updates at server-side. The first hypothesis was manually invalidated based on stored screenshots. No clear content update can be observed. Therefore, Instagram service might have been updated during our experimentation. Because accuracy was restored a few weeks after the first drop (mid June 2019), a plausible explanation is a restoration to a past configuration of the service. Amazon exhibits a similar drop toward the end of the period with the same plausible explanation of software changes altering the features.

B. Dynamic models

The previous experiment shows a decrease over a relatively low number of days. Although re-training the RF models can be done offline, we evaluate in this section an appropriate update frequency.

In this experiment, the models are updated every n d days assuming the last 15 days of traces for learning (sliding window). In Figure 3, n d is set to 5 and 10 days for comparison purpose. The vertical doted lines indicate the days when the models are updated. For readability, exact dates are not shown but are the same as in Figure 2.

In both cases, n d = 5 or n d = 10, retraining the models annihilates the impact of time on the classifier. As expecting, with a lower number of days n d between training, the accuracy In mostly all cases, re-learning the models regularly (about every week) with new captures can maintain the accuracy over the time.

C. Discussion

In general, accuracy in our current study is lower than our previous research work [START_REF] Brissaud | Transparent and service-agnostic monitoring of encrypted web traffic[END_REF]. It is thus natural to identify and discuss the reasons that can explain this situation. We noticed two major differences between the two evaluations:

• The division of traces for training and testing in [START_REF] Brissaud | Transparent and service-agnostic monitoring of encrypted web traffic[END_REF] are random but, for a given keyword, they are all captured successively in a few minutes. In this paper, a more realistic scenario is considered. The traces for learning and testing are collected over completely disjointed period of time. Actually, even a single day difference has an impact. • The selected keywords in this paper have been voluntary defined based on trends and topicality in order to have more volatile content and evaluate our technique in a worst case scenario. It leads to different accuracy per single keyword as shown in Figure 4 focusing on the 10% worst and best keywords (keywords with the worst/best mean accuracy for the whole experiment). Although a large span in the accuracy can be observed, it also clearly shows that for a part of keywords (top 10%) our technique is very efficient even in the case of Google.

We expect that these observations are also helpful for future research in traffic analysis in order to better customize the scenarios, according to predefined criteria (worst case analysis, realistic case) and for researchers interested in our datasets.

V. TEST-OF-SPACE

This section presents the evaluation of the classification for more than 3000 other services (test-of-space) based on the dataset presented in Section III-D.

As explained, a RF model is built for each service knowing that the built dataset provides 20 keywords per services. Because there are 20 traces by keyword, the learning stage is based on 16 traces while the remaining four are kept for testing.

The cumulative mean accuracy is reported in Figure 5 in comparison to a random classifier with an accuracy of 0.05 (i.e. select a random class out of 20). This figure shows that 50% of the websites have an accuracy higher or equals to 0.9. On the remaining 50%, 40% have an accuracy higher than 0.5.

Less than 2% of tested websites leads to an accuracy bellow the random cut-off. Through a manual analysis of the captured screenshots, two potential causes of misclassification have been identified:

• Whatever the specified keyword, the website always returns the same webpage. It is due to a non relevant keyword used or to an internal problem at the server side (e.g. maintenance in progress). • The returned page is different for each keyword but the page is constituted of a very few objects although H2Classifier has been designed assuming complex webpages composed of numerous objects as presented in section II-B. For example, the news website www.metro.us search results are displayed as a simple text list without images. For the very same reason, accuracy highly differs between Google (search engine) and Google Images as shown in section IV.

Reminding there was no particular selection of tested websites and that for 50% of them, the accuracy is higher than 0.9, the evaluated classification technique is quite generic and can be used off-the-shelf for a large amount of services. It further confirms its efficiency which was solely evaluated on five services in [START_REF] Brissaud | Transparent and service-agnostic monitoring of encrypted web traffic[END_REF]. In this paper, we evaluate the practical application of H2Classifier assuming its deployment in real scenarios where HTTPS traffic must be monitored over long time period and on many web services. Although H2Classifier has been specifically designed to monitor user activities in h2-based services (HTTP2 over TLS), the designed methodology and outcomes (results and discussions) are helpful for other encrypted traffic classification or fingerprinting techniques leveraging machinelearning. To summarize, the application of this type of technique is relevant and efficient for many different websites, assuming obviously the method has not been fitted to a particular one by design. However, the volatility of content in Internet as well as likely server software updates require regular re-learning, around every week in our case. As a future work, our plan is to evaluate if classification errors can be exploited to automatically detect falsified websites.

Fig. 1 :

 1 Fig. 1: Categories of the 3096 websites composing our URL list

Fig. 2 :

 2 Fig. 2: Mean accuracy over time for each day

Fig. 3 :

 3 Fig. 3: Mean accuracy every day with new learning every n d days

Fig. 4 :Fig. 5 :

 45 Fig. 4: Mean accuracy every day with new learning every 5 days for the 10% worst (top) and best (lower) keywords

TABLE I :

 I Temporal dataset

	Test-of-time
	Services	Amazon, Instagram, Google, Google Images
	Number of keywords per service	500
	# traces per keyword per day	4
	Number of hosts for capture	1
	Total size of dataset	1.8 TB
	attribute name='q' 89 but we also consider the attribute names
	search, research or searchbar.

TABLE II :

 II Services dataset

	Test-of-space
	Number of Services	3096 different
	Number of keywords per service	20
	Number of traces per keyword	20
	Number of hosts for capture	6
	Total size of dataset	6 TB
	N e w s a n d M e d ia B u s in e s s S h o p p in g In f o r m a t io n T e c h n o lo g y E n t e r t a in m e n t E d u c a t io n S p o r t s H e a lt h a n d W e ll n e s s S o c ie t y a n d R e s t a u r a n t a n d D in in g L if e s t y le

http://betternet.lhs.inria.fr/datasets/

https://www.seleniumhq.org/

https://github.com/mozilla/geckodriver

https://www.tcpdump.org/

https://www.w3schools.com/tags/att_input_type_search.asp

https://developer.mozilla.org/en-US/docs/Web/HTML/Element/input/ search\#Basic_example

https://www.quantcast.com/top-sites/US

ACKNOWLEDGMENT

This work has been partially supported by the NATO Science for Peace and Security Programme under grant G5319 ThreatPredict, the European Union's Horizon 2020 research and innovation programme under grant agreement No 830927 and from the French National Research Agency (ANR), project number ANR-19-CE39-0011-01.


Floriane Poignant

Hela Charfi

Chen-Hui Chan

Elise Dumont

David Loffreda

Benoit Gervais

Michaël Beuve
email: michael.beuve@univ-lyon1.fr

Monte Carlo simulation of free radical production under keV photon irradiation of gold nanoparticle aqueous solution. Part II: Local primary chemical boost

Keywords: local primary chemical boost gold nanoparticles, photon irradiation, water radiolysis, free radical production, Monte Carlo simulation

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Over the past two decades, the use of high-atomic-number (Z) nanoparticles (NPs) has been of high interest due to their radiosensitizing properties. Several pre-clinical studies have shown the efficiency of gold NPs (GNPs) to enhance the effect of radiation therapy, in particular when using low energy (keV) X-rays 5 [START_REF] Hainfeld | The use of gold nanoparticles 275 to enhance radiotherapy in mice[END_REF][START_REF] Her | Gold nanoparticles for applications in cancer radiotherapy: Mechanisms and recent advancements[END_REF]. The origin of NP enhanced radiotherapy may originate from complex physical, chemical and biological mechanisms. The relative contribution of each stage remains under investigation [START_REF] Brun | Actual questions raised by nanoparticle radiosensitization[END_REF]. Theoretical approaches, such as Monte Carlo (MC) simulations, may help to better understand how early physical and chemical stages could impact a biological system and lead, for instance, to cell 10 death increase.

MC studies have shown that, for keV X-rays, a significant increase of dose deposition may be reached for sufficiently large gold concentration [START_REF] Cho | Estimation of tumour dose enhancement due to gold nanoparticles during typical radiation treatments: a preliminary Monte Carlo study[END_REF][START_REF] Cho | The dosimetric feasibility of gold nanoparticle-aided radiation therapy (GNRT) via brachytherapy using lowenergy gamma-/x-ray sources[END_REF]. At such energies, the photo-electric effect dominates for high-Z materials and the photon interaction probability with GNP is much larger than the interaction 15 probability with biological tissues. Other MC studies have suggested that the heterogeneity of dose deposition in the vicinity of NPs may further contribute to the enhancement, due to an hypothetical release of low-energy electrons.

Following a photo-ionization of a high-Z NP, desexcitation processes result in the emission of a cascade of secondary electrons called Auger cascade. It might 20 be responsible for a high energy deposition in the vicinity of the GNP, and has often been associated with GNP efficiency, by presumably inducing an increased biological effectiveness similar to that observed in hadrontherapy [START_REF] Cho | Estimation of tumour dose enhancement due to gold nanoparticles during typical radiation treatments: a preliminary Monte Carlo study[END_REF][START_REF] Cho | The dosimetric feasibility of gold nanoparticle-aided radiation therapy (GNRT) via brachytherapy using lowenergy gamma-/x-ray sources[END_REF][START_REF] Mcmahon | Nanodosimetric effects of gold nanoparticles in megavoltage radiation therapy[END_REF][START_REF] Mcmahon | Biological consequences of nanoscale energy deposition near irradiated heavy atom nanoparticles[END_REF][START_REF] Lechtman | A Monte Carlo-based model of gold nanoparticle radiosensitization accounting for increased radiobiological effectiveness[END_REF].

While there have been many MC studies of the physical stage and dose deposition, simulation studies of the chemical stage are scarce [START_REF] Xie | Simulation on the molecular radiosensitization effect of gold nanoparticles in cells irradiated by x-rays[END_REF][START_REF] Tran | Geant4 Monte Carlo simulation of absorbed dose and radiolysis yields enhancement from a gold nanoparticle under MeV proton irradiation[END_REF][START_REF] Rudek | Radio-enhancement by gold nanoparticles and their impact on water radiolysis for x-ray, proton and carbon-ion beams[END_REF]. In particular,

25

GNP impact on water radiolysis has hardly been investigated, although it was suggested to be a crucial step that may connect the physical effect to the biological consequences [START_REF] Brun | Actual questions raised by nanoparticle radiosensitization[END_REF]. In a previous work [START_REF] Poignant | Monte Carlo simulation of free radical production under keV photon irradiation of gold nanoparticle aqueous solution[END_REF], we studied, on a macroscopic scale, the impact of GNPs on free radical yields for 20-90 keV photon irradiation.

We showed that for a GNP concentration of 1 mg • mL -1 , the increase of yields 30 is well correlated to the macroscopic dose enhancement. Yields were found to mostly depend on the photon energy and, to a lesser extent, on the NP radius for a given mass concentration of GNPs. As previous MC studies at nanometric scale pointed out the importance of dose heterogeneity and, in particular, the high energy deposition in the vicinity of the GNP [START_REF] Mcmahon | Energy dependence of gold nanoparticle radiosensitization in plasmid DNA[END_REF][START_REF] Lechtman | A Monte Carlo-based model of gold nanoparticle radiosensitization accounting for increased radiobiological effectiveness[END_REF][START_REF] Zygmanski | Dependence of Monte Carlo microdosimetric computations on the 325 simulation geometry of gold nanoparticles[END_REF][START_REF] Lin | Gold nanoparticle induced vasculature damage in radiotherapy: comparing protons, megavoltage photons, and kilovoltage photons[END_REF][START_REF] Martinov | Heterogeneous multiscale Monte Carlo simulations for gold nanoparticle radiosensitization[END_REF][START_REF] Zygmanski | Nanoscale radiation transport and clinical beam 335 modeling for gold nanoparticle dose enhanced radiotherapy (GNPT) using X-rays[END_REF][START_REF] Sakata | Geant4-DNA track-structure simulations for gold nanoparticles: The importance 340 of electron discrete models in nanometer volumes[END_REF],

35
we wished to focus, in this study, on the production of chemical species at the nanoscale. In particular we searched for a possible boost of radical production near the GNP, which might lead to local enhancement of radical recombination.

The goal of this work was thus to investigate, in a systematic way, the production of free radicals for keV photon irradiation in water around one GNP, following geometry are presented. Section 3 presents the radial distribution of chemical species around a NP. Section 4 discusses and concludes on the results. 2. Material and methods

Geometry

The geometry considered for the calculation is described in Fig. 1. We defined a GNP placed at the center of a volume V filled with liquid water. The half length of V was chosen to be larger than the maximum range of the most energetic 55 photo-electron in water. The GNP was ionized by a photon whose energy varied from 20 to 90 keV. The half length of the volume was hence larger than 160 µm.

The systematic study of the NP-size impact on chemical species yields were performed with four NP radii (5 nm, 12.5 nm, 25 nm and 50 nm). To analyze the specificity of gold material, the calculation was also performed for a spherical 60 NP made of water with the same radius.

Monte Carlo simulation for water radiolysis

Our calculation was based on the Monte Carlo simulation MDM, and the full description of the model is available elsewhere [START_REF] Poignant | Monte Carlo simulation of free radical production under keV photon irradiation of gold nanoparticle aqueous solution[END_REF][START_REF] Poignant | Theoretical derivation and benchmarking of cross sections for low-energy electron transport in gold[END_REF][START_REF] Gervais | Numerical simulation of multiple ionization and high LET effects in liquid water radiolysis[END_REF]. In short, 3 consecutive stages were performed to calculate the yield of chemical species. The 65 physical stage consisted of an event-by-event tracking of electrons in the different media (water or gold) down to thermalization energy at 300 K (37 meV).

Photons were not explicitly tracked, and we only considered the photoelectric effect and the Compton scattering as they dominate at the investigated photon energies. The photon ionization impact was randomly distributed within the 70 NP. At the end of the physical stage, the surrounding water molecule may be

found either excited (H 2 O * ) or ionized (H 2 O + , H 2 O 2+ and H 2 O -)
, and the medium contained thermalized electrons. The electrons in GNP for which the energy was too low to escape the GNP surface were no longer tracked, as they did not lead to additional water molecule excitation or ionization.

75

This stage was followed by a physico-chemical stage during which the medium relaxes, leading to the production of primary chemical species. The different branching ratios of each molecular rearrangement are given elsewhere [START_REF] Gervais | Numerical simulation of multiple ionization and high LET effects in liquid water radiolysis[END_REF]. At the end of this stage, the outcome consisted of spatial distribution of chemical species in water, with a majority of e - aq , H 3 O + and • OH.

80

Finally, during the chemical stage, the primary chemical species diffused and interacted with each other. More than 50 chemical reactions were considered.

In our simulation, this phase was simulated up to 10 -6 s. The method used to calculate the yields and to optimize the computing time, as well as the list of standard chemical reactions of water radiolysis, are available elsewhere and will 85 not be discussed in this paper [START_REF] Gervais | Numerical simulation of multiple ionization and high LET effects in liquid water radiolysis[END_REF][START_REF] Colliaux | O2 and glu-350 tathione effects on water radiolysis: a simulation study[END_REF][START_REF] Colliaux | Simulation of ion-induced water radiolysis in different conditions of oxygenation[END_REF][START_REF] Poignant | Monte Carlo simulation of free radical production under keV photon irradiation of gold nanoparticle aqueous solution[END_REF]. In this study, we did not consider any possible specific reactions at the surface of the GNP.

Calculated quantities

We calculated 4 different quantities to characterize the production and distribution of chemical species around one NP. We calculated (1) the mean number 90 of chemical species X following one single ionization event in the NP R 1 X ; (2) the radial number of chemical species, which is the differential of this number in the radial distance

∂R 1 X ∂r ;
(3) the radial concentration of the chemical species X, c X (r), as the concentration of X averaged over the spherical shell of radius r and of thickness δr and centered around the GNP, as represented in dashed 95 green and blue in Fig. 1. This quantity is related to the previous one through the relation,

c X (r) = 1 4/3π((r 3 + δr) -r 3 ) r+δr r ∂R 1 X ∂r dr . (1)
Finally, we calculated (4) the spherical concentration defined as concentration per ionized NP, averaged over sphere centered on the nanoparticle and limited by the nanoparticle surface on one side, and the radial distance r on the other 100 side, as represented in blue in Fig. 1,

C X (r) = 1 4 3 π(r 3 -R 3 NP ) r RNP ∂R 1 X ∂r dr . (2)
To have a common origin for all NP radii, the results are presented as a function of the distance from NP surface, x = r-R NP . We shall refer to c X (r) and C X (r) as c X (x) and C X (x).

Results

105

Average number of chemical species per nanoparticle ionization

Tab. 1 gives the average number of chemical species (R 1

• OH and R 1 H2O2 ) for various E photon , two NP radii (R NP = 5 and 50 nm), calculated both at 10 -12 and 10 -6 s following the photon absorption.

For the smallest GNPs, the number of • OH produced at 10 -12 s varied from 721 110 (E photon = 20 keV), up to 2912 (E photon = 70 keV), due to the increase of the photo-electron energy. It then dropped at E photon = 90 keV, as there is a high fluorescence emission in this case. The mean free path of the emitted photons is very large compared to the volume of the simulated system and the probability of interaction is low. This latter was neglected in our simulations. If the GNP 115 was embedded in a volume of water larger than the photon mean free path, the photon would deposit its energy in the volume, but most of the time far away (several mm at least) from the ionized GNP. The same variation of R 1

• OH and R 1

H2O2 was obtained with regard to E photon , regardless of other parameters (e.g., radius, time ...). The number of chemical species decreased with increasing GNP 120 radius due to the energy loss in the GNP. This was particularly visible at 20 keV: at 10 -12 s, the number of • OH decreased from 721 to 571 when the radius of the GNP increased from 5 nm to 50 nm. The same trends were observed with

H 2 O 2 .
For WNP, the number of chemical species was hardly impacted by the NP size.

125

The number of chemical species, at fixed time, first decreased from 20 keV to 50 keV, before increasing again. The first drop was due to the fact that, at 20 keV, the photo-electric effect remained non-negligible, in which case all the photon energy was deposited. As the Compton effect became dominant, only part of the photon energy was deposited. This was compensated by the fact 130 that the average Compton electron energy increased with the photon energy, so that the number of species per photo-ionization increased above 50 keV.

The number of chemical species generated per photo-ionization was systematically higher for GNPs than WNPs, with two exceptions. For E photon = 20 keV and R NP = 50 nm, the number of chemical species for GNPs was below that of 135 WNP, due to a particularly high energy absorption within the GNP. At 10 -6 s, E photon = 20keV, R NP = 5 nm, and regardless of the time, the number of • OH for GNP was also below that of WNP, while it was not the case of H 2 O 2 . This denotes the fact that in this specific case, the recombination of • OH was particularly important, leading to a high number of H 2 O 2 . Finally, as expected, 140 the largest difference between GNPs and WNPs was obtained at 70 keV. In this case, the photo-electron emitted from the GNP had the largest kinetic energy, while the Compton electron in water still had a limited kinetic energy. The radial concentration was maximum near the NP surface and dropped severely for higher values of r. This decrease was proportional to 1/r 2 . It is due to the dependency in r 2 of spherical shell volume when r increases, and to a lesser extent to the differential number of chemical species. Such tendency was similar 150 regardless of the photon energy, both at 10 -12 s (not shown here) and 10 -6 s.

The global shape of the radial distribution was not specific to the material, as WNP showed the same tendency.

Due to • OH recombination and diffusion, the radial concentration decreased with increasing t over the whole r range for • OH. For H 2 O 2 (not shown here), 155 c H2O2 (x) dropped for x < 200 nm over the time period 10 -12 to 10 -6 . This drop was due to outward diffusion, which is of the order of 100 nm for 1 µs. escape the NP and less species are produced close to the surface.

The spatial distribution of chemical species also depended on photon energy.

To better analyze this dependency, we calculated the radial number of chemical species, instead of the radial concentration as the latter was strongly impacted by the r 2 normalization. Fig. 3 displays the radial number of chemical species 170

• OH,

∂R 1
• OH ∂r , as a function of the distance to the NP surface x.

For GNP, some structures are visible around x = 100 nm and 1000 nm, regardless of the photon energy. The first structure came from LLS (S > N) Auger electrons emitted with an energy lower than 2 keV, which had ranges lower than ∼ 160 nm in water. The second structure came from LMS (S > M) Auger elec-175 trons emitted with an energy between 8-12 keV, which corresponded to a range in water around 1000-2000 nm. Another structure is visible for large x, which corresponded to the end of photo-electron path. As shown in Tab. 2, these photo-electrons have different energies depending on the photon energy and the shell from which they are ejected. This phenomenon was particularly visible for 180 E photon = 90 keV. As the probability to ionize the K-shell was overwhelming in this case, many photo-electrons were ejected with an energy around 9 keV.

These electrons therefore contributed to the boost of the number of chemical species observed around 1000 nm. For WNP, the progressive decrease of ∂R 1

• OH ∂r with

x was related to the energy 185 distribution of Compton electrons. An exception was noticeable for 20 keV: as the photo-electric effect remained non-negligible, the radial number of chemical species was a combination of the radial number following a photo-electric event producing a ∼ 19.5 keV photo-electron, and the radial number following a Compton event producing a broad energy distribution below 1400 eV.

190

When comparing GNP and WNP, the first striking fact was that ∂R 1

• OH ∂r remained higher at longer distances (x > ∼ 10 000 nm) for GNP. This is explained by the fact that, on average, the Compton electron had an energy lower than a photo-electron energy. At very small distance (x < 100 nm), the number of chemical species was higher for GNP compared to WNP when R NP = 5 nm. in keV) and range of the photo-electron in water (R. in µm) following an ionization event in gold with a photon of energy E photon .

To better quantify the relative difference between GNP and WNP, we calculated the difference of spherical concentration of GNP and WNP, to that of WNP,

C X,GNP (x) -C X,WNP (x) C X,WNP (x) , (3)
for the two extrema sizes of NPs investigated here. Results are shown in Tab. 3 for • OH. The results confirmed that, for R NP = 50 nm, the spherical concen-200 tration of WNPs for small volumes (10 and 100 nm thick) were higher than that of GNP, except when the photon energy was equal to 90 keV. This was due to non negligible energy loss in the GNP. When the considered volume increased (x = 1000 nm), the concentration for GNPs was higher than that of WNP. For even larger volumes (x = 10 µm), the concentrations were higher for GNPs,

205

with one exception (E photon = 20 keV). This translated the fact that, on average, GNPs generated more chemical species than WNP following one photon ionization event, but far away from the GNP (i.e., beyond 100 nm).

Discussion and conclusion

In this work, we investigated the production and spatial distribution of radical 210 species following keV-photon ionization of a GNP or a WNP embedded in a volume of water. We focus on the local production of radicals, meaning that molecules following an ionization event in GNPs induced a minor boost in the radial number. This boost appeared around 100 nm and 1000 nm.

(3) Both • OH and H 2 O 2 had very similar distributions. In particular, H 2 O 2 230 production was not particularly higher due to Auger electrons. The Auger electrons are indeed significant for GNP ionization, but their effect is rather comparable to that of Compton electrons produced by WNP ionization.

(4) The concentration of chemical species was not necessarily higher for GNP compared to WNP near the ionized NP (within a few hundreds of nm), despite 235 the Auger electron cascade when the GNP was large. The absorption of lowenergy electrons within the GNP produced a drop of the number of chemical species near the surface with respect to the case of smaller GNPs.

Attempting an extension of these results to biological systems, one may expect that, on average, a cellular target containing or located near a large ionized 240 GNP (x < 200 nm) does not suffer a higher concentration of chemical species, compared to the ionization of a water molecule. However, it may experience a higher occurrence of this event, as the cross section for gold-photon interaction is higher than that of water-photon interaction at these energies. In other words, the major effect of gold is not to induce a high concentration boost of chemical 245 species compared to water, but rather to increase the chances of having an ionization, due to a high probability of photon-gold interaction compared to photon-water interaction.

In interpreting our results from a biological perspective, it is important to keep in mind that biological media contain many scavenging elements. Besides, 250 cells might catalyse complex reactions inducing an additional boost of chemical species production. Our calculation may therefore only characterize the primary chemical events at very short time generated by one photo-ionization event. In vitro data showed that GNPs often appear clustered together in vesicles [START_REF] Chithrani | Determining the size and shape dependence of gold nanoparticle uptake into mammalian cells[END_REF][START_REF] Chithrani | Gold nanoparticles as radiation sensitizers in cancer therapy[END_REF][START_REF] Bobyk | Photoactiva-tion of gold nanoparticles for glioma treatment[END_REF][START_REF] Jeynes | Investigation of gold nanoparticle radiosensitization mechanisms using a free radical scavenger and protons of different energies[END_REF][START_REF] Retif | Monte Carlo simulations guided by imaging to predict the in vitro ranking of radiosensitizing nanoparticles[END_REF], which might roughly be represented by large GNPs.

255

Our calculated spherical concentrations may be used to estimate the primary chemical boost experienced by cellular components when they internalize GNPs.

Typically, a size of 10 nm may represent the dimension of a membrane; a size of 100 nm may represent the dimensions of an organelle (mitochondria, lysosome); a size of 1000 nm may represent the dimensions of a large organelle (large lyso-260 some); a size of 10 000 nm may represent a nucleus. As our calculation showed, the number of chemical species following an ionization event decreases for large GNPs, especially in its vicinity. At first sight, it thus appears less advantageous to have aggregates rather than isolated GNPs. However, if GNPs accumulate in a specific zone of the cell, an accumulation of damages may occur near the 265 GNPs, due to a high probability of photon-gold interaction. Such accumulation might be deleterious for the cell. The authors thank Jean-Michel Moreau for his precious help.

Acknowledgment

40Figure 1 :

 1 Figure 1: Schematic view of the geometry. Orange dots represent excited or ionized water molecules or chemical species produced following the NP ionization. On the right, the dashed green and blue area represents the scoring volume used to calculate radical species concentration in spherical concentric shells, c X (x) and C X (x).

  50

Figure 2 :

 2 Figure 2: Radial concentration c• OH (x) as a function of x (E photon = 50 keV and t = 10 -6 ), for two R NP . Solid line shows results for a GNP and dashed line for a WNP.

Fig. 2

 2 Fig. 2 shows the radial concentration c• OH (x) of • OH following an ionization 145

  -𝑅 NP (nm) 𝑥 = 𝑟 -𝑅 NP (nm)

Figure 3 : 1 •

 31 Figure 3: Radial number of chemical species ∂R 1 • OH ∂r as a function of x. In each panel, results are given for GNP (solid line) and WNP (dashed line). Results are displayed for E photon = 20, 50, 70 and 90 keV. R NP = 5 nm, and t = 10 -6 s.

  This work was supported by the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-270 IDEX-0007) operated by the French National Research Agency (ANR). This work has been supported by the Fondation ARC pour la recherche sur le cancer.

Table 1 :

 1

			10 -12 s			10 -6 s
		• OH	H 2 O 2		• OH	H 2 O 2
	Energy GNP WNP GNP WNP GNP WNP GNP WNP
				R NP = 5 nm
	20 keV	721	673	33	28	256	288	135	118
	50 keV 2070	460	81	20	934	198	350	81
	70 keV 2912	501	112	21	1385	211	482	88
	90 keV 1560	618	63	25	671	264	270	108
				R NP = 50 nm
	20 keV	571	674	25	28	212	288	107	118
	50 keV 1960	455	76	19	906	196	329	80
	70 keV 2820	500	107	21	1365	211	462	88
	90 keV 1401	617	56	26	625	263	239	108

Average number of chemical species produced per GNP/WNP ionization for different times, photon energies and GNP radii. The statistical uncertainty was below 1 %.

  On the contrary, H 2 O 2 increased for larger radii. Increasing the NP radius R NP decreased c• OH (x) below 200 nm, regardless of photon energy, time t, chemical species or material. There are two reasons for 160 this. First, by definition of the radial concentration, the number of species is normalized to the volume of the shell. For a given value of x, this volume increases with increasing R NP . Second, the part of the absorbed energy inside the NP, which increases with increasing R NP , corresponds to the stopping of low-energy electrons. When increasing the NP radius, less low-energy electrons

	165

Table 2 :

 2 Probability (P. in %) of ionizing the shell K,L or above, photo-electron energy (E.

	h
	Atomic shell		K			L			≥ M
	E photon	P. E. R. P.	E.	R.	P.	E.	R.
	20 keV	0	-	-	77	6-8	0.9-1.5 23 > 17 5.3-7.4
	50 keV	0	-	-	77 17-18	5.3-6	23 > 47	35-42
	70 keV	0	-	-	77 56-58	49-53	23 > 67	68-74
	90 keV	77 9 1.8 16 76-78	88-92	7 > 87 > 111

On the contrary, it was lower for GNP compared to WNP when R NP = 50 nm.

x 10 nm 100 nm 1 000 nm 10 000 nm 100 000 nm t = 10 -6 s -R NP = 5 nm [START_REF] Gervais | Numerical simulation of multiple ionization and high LET effects in liquid water radiolysis[END_REF] (2) Chemical species produced by the interaction of Auger-electrons with water