Ruoyang Zhai

Raiko Blondonnet

Ebrahim Ebrahimi

Corinne Belville

Jules Audard

Christelle Gross

Helena Choltus

Fanny Henrioux

Jean-Michel Constantin

Bruno Pereira

Loic Blanchon

Vincent Sapin

M.D Matthieu Jabaudon
email: mjabaudon@chu-clermontferrand.fr

The receptor for advanced glycation end-products enhances lung epithelial wound repair: An in vitro study

Keywords: receptor for advanced glycation end-products, lung epithelial injury, wound healing, cell migration, cell proliferation

à 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

Alteration of alveolar epithelial barrier integrity with increased epithelial permeability plays a central role in the pathophysiology of the acute respiratory distress syndrome (ARDS) [START_REF] Matthay | Intact epithelial barrier function is critical for the resolution of alveolar edema in humans[END_REF][START_REF] Ware | Alveolar fluid clearance is impaired in the majority of patients with acute lung injury and the acute respiratory distress syndrome[END_REF][START_REF] Guillot | Alveolar epithelial cells: master regulators of lung homeostasis[END_REF][START_REF] Villar | Lung Repair and Regeneration in ARDS: Role of PECAM1 and Wnt Signaling[END_REF], a severe condition associated with significant morbidity and mortality, long-term physical or psychological consequences, and high healthcare-associated costs [START_REF] Herridge | Canadian Critical Care Trials Group, Functional disability 5 years after acute respiratory distress syndrome[END_REF][START_REF] Bellani | Epidemiology, Patterns of Care, and Mortality for Patients With Acute Respiratory Distress Syndrome in Intensive Care Units in 50 Countries[END_REF][START_REF] Thompson | Acute Respiratory Distress Syndrome[END_REF]. The lung alveolar epithelial barrier is composed of alveolar type (AT)-I cells, covering 95% of the alveolar epithelial surface, and of AT-II cells. Previous research suggests that re-epithelization of the alveolar surface is a key process that is required to restore the integrity and function of the lung alveolar epithelial barrier [START_REF] Matthay | Acute respiratory distress syndrome[END_REF]. Although this repair process may involve the proliferation and migration of AT-II cells, its precise mechanisms remain underinvestigated [START_REF] Matthay | Acute respiratory distress syndrome[END_REF][START_REF] Adamson | The type 2 cell as progenitor of alveolar epithelial regeneration. A cytodynamic study in mice after exposure to oxygen[END_REF][START_REF] Mcqualter | Evidence of an epithelial stem/progenitor cell hierarchy in the adult mouse lung[END_REF].

The receptor for advanced glycation end-products (RAGE) is a membrane signal transduction receptor that can induce various intracellular pathways by binding to multiple ligands, such as S100 proteins, high-mobility group box 1 protein (HMGB1), and advanced glycation end-products (AGEs), among others [START_REF] Xie | Cellular signalling of the receptor for advanced glycation end products (RAGE)[END_REF][START_REF] Lander | Activation of the 13 receptor for advanced glycation end products triggers a p21(ras)-dependent mitogen-activated protein kinase pathway regulated by oxidant stress[END_REF][START_REF] Tafani | Hypoxia-increased RAGE and P2X7R expression regulates tumor cell invasion through phosphorylation of Erk1/2 and Akt and nuclear translocation of NF-{kappa}B[END_REF][START_REF] Huttunen | Receptor for advanced glycation end products (RAGE)mediated neurite outgrowth and activation of NF-kappaB require the cytoplasmic domain of the receptor but different downstream signaling pathways[END_REF][START_REF] Ranzato | HMGb1 promotes scratch wound closure of HaCaT keratinocytes via ERK1/2 activation[END_REF][START_REF] Hudson | Targeting RAGE Signaling in Inflammatory Disease[END_REF]. RAGE is abundantly expressed in the lung tissue, both as transmembrane full-length and soluble forms, and RAGE pathway may play a key role in lung homeostasis [START_REF] Mukherjee | Implication of receptor for advanced glycation end product (RAGE) in pulmonary health and pathophysiology[END_REF][START_REF] Queisser | Loss of RAGE in pulmonary fibrosis: molecular relations to functional changes in pulmonary cell types[END_REF]. Recent clinical and preclinical research has suggested that soluble RAGE (sRAGE) could be a biomarker of lung epithelial injury associated with severity and outcome in ARDS [START_REF] Shirasawa | Receptor for advanced glycation end-products is a marker of type I lung alveolar cells[END_REF][START_REF] Uchida | Receptor for advanced glycation end-products is a marker of type I cell injury in acute lung injury[END_REF][START_REF] Calfee | NHLBI ARDS Network, Plasma receptor for advanced glycation end products and clinical outcomes in acute lung injury[END_REF][START_REF] Mauri | Elevated plasma and alveolar levels of soluble receptor for advanced glycation endproducts are associated with severity of lung dysfunction in ARDS patients[END_REF][START_REF] Jabaudon | Soluble form of the receptor for advanced glycation end products is a marker of acute lung injury but not of severe sepsis in critically ill patients[END_REF][START_REF] Jabaudon | Soluble Receptor for Advanced Glycation End-Products Predicts Impaired Alveolar Fluid Clearance in Acute Respiratory Distress Syndrome[END_REF][START_REF] Jabaudon | Plasma sRAGE is independently associated with increased mortality in ARDS: a meta-analysis of individual patient data[END_REF]. Moreover, RAGE inhibition or AGER gene deficiency decreased lung injury in animal models of experimental ARDS from both septic or non-septic causes [START_REF] Van Zoelen | The receptor for advanced glycation end products impairs host defense in pneumococcal pneumonia[END_REF][START_REF] Ramsgaard | Lack of the receptor for advanced glycation end-products attenuates E. coli pneumonia in mice[END_REF][START_REF] Blondonnet | RAGE inhibition reduces acute lung injury in mice[END_REF][START_REF] Audard | Inhibition of the Receptor for Advanced Glycation End-Products in Acute Respiratory Distress Syndrome: A Randomised Laboratory Trial in Piglets[END_REF]. Although the RAGE pathway has a major role in initiating and perpetuating nuclear factor-kappa B (NF-κB)mediated inflammation, it also contributes, among many other potential effects [START_REF] Xie | Cellular signalling of the receptor for advanced glycation end products (RAGE)[END_REF], to repair mechanisms such as peripheral nerve regeneration through the recruitment of axonal outgrowth pathways [START_REF] Huttunen | Receptor for advanced glycation end products (RAGE)mediated neurite outgrowth and activation of NF-kappaB require the cytoplasmic domain of the receptor but different downstream signaling pathways[END_REF][START_REF] Rong | RAGE modulates peripheral nerve regeneration via recruitment of both inflammatory and axonal outgrowth pathways[END_REF] or extracellular matrix synthesis and wound repair in human bronchial epithelial cells through stimulation by HMGB1 [START_REF] Ojo | High-mobility group box 1 promotes extracellular matrix synthesis and wound repair in human bronchial epithelial cells[END_REF]. However, whether modulation of RAGE pathway might affect lung alveolar epithelial wound repair remains unknown.

In this study, we investigated whether treatment with RAGE agonists (HMGB1 and AGEs), alone or combined with a RAGE antagonist, would affect wound healing of lung alveolar epithelial cells in vitro. We further hypothesized that modulation of the RAGE pathway might affect some major processes of wound healing in these cells, such as proliferation or migration.

Materials and Methods

Reagents and antibodies

Cell culture Roswell Park Memorial Institute (RPMI 1640) medium 1640, fetal bovine serum (FBS), antibiotics (streptomycin, penicillin, amphotericin B), and phosphate buffered saline (PBS) were obtained from Fisher Scientific (Illkirch, France). Human recombinant HMGB1 was obtained from Sigma-Aldrich (Saint Quentin-Fallavier, France). AGEs (ab51995) was purchased from Abcam (Paris, France). RAGE antagonist peptide (RAP) was obtained from Calbiochem® Merck Millipore (Molsheim, France).

Cell culture

Human adenocarcinoma-derived alveolar AT-II epithelial cells (A549) were obtained from the American Tissue Culture Collection (CCL-185, LGC Standard®, Molsheim, France). After thawing, cells were cultured in a 95% humidified air atmosphere at 37°C with a CO2 concentration of 5% in a RPMI-1640 culture medium enriched with 5% FBS and supplemented with 10 mg/mL streptomycin, 10,000 U/mL penicillin, and 25 µg/mL amphotericin B (Fisher Scientific, France).

In vitro model of alveolar epithelium wound healing (scratch assay)

Cells (2x10 5 per well) were cultured for 72 h in six 9.5 cm 2 well-plates (Corning, Samois-sur-Seine, France) until 100% confluence was obtained. The tip of a 200-µL pipette was applied perpendicular to the surface of the cell monolayer to create a thin "wound". After three washes with PBS (1X), wounded cells were either left untreated (control) or treated with RAGE agonists alone or combined with RAP (12.5 µg/mL), before being incubated at 37°C with 95% humidified air and 5% CO2. The dose of RAP was chosen based on the findings from previous studies [START_REF] Audard | Inhibition of the Receptor for Advanced Glycation End-Products in Acute Respiratory Distress Syndrome: A Randomised Laboratory Trial in Piglets[END_REF][START_REF] Arumugam | S100P-derived RAGE antagonistic peptide reduces tumor growth and metastasis[END_REF][START_REF] Bongarzone | Targeting the Receptor for Advanced Glycation Endproducts (RAGE): A Medicinal Chemistry Perspective[END_REF][START_REF] Uderhardt | Resident Macrophages Cloak Tissue Microlesions to Prevent Neutrophil-Driven Inflammatory Damage[END_REF]. The wound area was observed immediately after the scratch (h0) and 12, 24, and 48 h after treatment using a Nikon polarized optical microscope with 40x magnification and a digital camera body. Images were recorded prior to their analysis using the Fiji image processing software (Wayne Rasband, National Institute of Mental Health, Bethesda, Maryland, USA) [START_REF] Schindelin | Fiji: an open-source platform for biological-image analysis[END_REF].

Measurements of the wounded area were made from edge to edge by manually delimiting the lesion surface. This experiment was repeated six times (with each condition performed in duplicate). After testing three concentrations of RAGE agonists HMGB1 (0.1, 0.3, and 0.6 µg/mL) and AGEs (10, 30, and 90 µg/mL) for their effect on wound area, subsequent experiments were performed using the minimum effective concentration of each.

Cell migration assay

Epithelial cell migration was assessed with the CytoSelect™ 96-well cell migration assay (CBA-106; 8 μm Fluorometric format; Biolabs, London, United Kingdom). At 12 and 48 h after scratch wounding, cells were suspended in a serum-free medium and placed in the upper chamber. Three hundred microliters of chemoattractant media containing 5% FBS was then added to the lower chamber. After a 24 h incubation, migratory cells that had passed through the polycarbonate membrane (8 μm pore size) were dislodged from the underside of the membrane with a detachment solution. The dislodged cells were then stained with CyQuant® GR Dye (Fisher Scientific, Hampton, USA), diluted in a lysis buffer, and quantified by fluorescence measurement at 480 nm/520 nm (Zeiss LSM 800, Carl Zeiss SAS, Marly Le Roi, France). This experiment, performed in triplicate, was independently repeated three times.

Cell proliferation assay

At 12 and 48 h after scratch wounding, A549 cells were stained with 5-Bromo-2'-deoxy-uridine (BrdU) at 10 µM using a BrdU Labeling and Detection Kit II (11299964001; Roche Diagnostics, Meylan, France). After 45 minutes, they were washed three times with PBS, and fixed with an ethanol fixative solution (50 mM glycine, pH=2, in 70% ethanol) for 20 minutes at -20°C. After three washes in PBS, the cells were incubated with an anti-BrdU antibody (1/25) for 60 minutes at 37°C. The cells were washed three times with PBS, followed by incubation with an anti-mouse-Igalkaline phosphatase (1:25 dilution) for 60 minutes at 37°C. The color substrate solution was added to the cells after three washes in PBS, and the bound anti-BrdU antibody was visualized by light microscopy (Zeiss Axio Observer, Carl Zeiss SAS, Marly Le Roi, France). Proliferation was expressed as the ratio of BrdU-positive cells to the total number of cells. This experiment, performed in triplicate, was independently repeated three times.

Statistical analysis

All analyses were performed using Prism 6 (GraphPad Software, La Jolla, CA, USA) and Stata version 15 (StataCorp, College Station, TX, USA). The tests were two-sided, with a type I error set at 5%. Continuous data were expressed as mean ± standard deviation (SD) or median [interquartile range] according to statistical distribution (assumption of normality was assessed using the Shapiro-Wilk test). Continuous parameters were compared between experimental groups using analysis of variance or the Kruskal-Wallis test when t-test assumptions were not met (normality studied using Shapiro-Wilk and assumption of homoscedasticity verified by Fisher-Snedecor test).

Random effects models were used to analyze the longitudinal evolution of variables: (i) considering between-and within-experiment variability (random subject effects: random intercept and slope) and (ii) evaluating fixed effects: group, time, and time by group interaction. The residual normality was checked for all models. Values were log-transformed for all variables to achieve normality prior to performing random effects models.

Results

Both HMGB1 and AGEs promote dose-dependent and RAGE-dependent wound healing of lung alveolar epithelial cells

We first tested different treatment concentrations of HMGB1 and AGEs for their effect on wound area after the scratch assay (Figure 1 and Table 1). Both RAGE ligands had dose-and timedependent effects towards better wound healing of epithelial cells. Compared to controls, cells treated with HMGB1 (0.1 µg/mL) had better wound healing (time by group interaction, p = 0.008) at 24 h and 48 h after injury (p = 0.04 and p <10 -4 , respectively, for post-hoc comparisons). A similar effect was found after treatment with AGEs (30 µg/mL) (time by group interaction, p <10 - 4 ), compared to untreated cells, at 24 h and 48 h after injury (p = 0.005 and p = 0.0001, respectively, for post-hoc comparisons) (p = 10 -4 ). The minimum effective doses of HMGB1 (0.1 µg/mL) and AGEs (30 µg/mL) were used in subsequent experiments.

Co-treatment with RAP (12.5 µg/mL) significantly reduced the effects of both HMGB1 (time by group interaction, p <10 -4 ) and AGEs (time by group interaction, p <10 -4 ) on wound area, with most detrimental effects on wound healing observed at 48 h for HMGB1 and at 12 h for AGEs (Figure 2 andTable 2).

HMGB1 and AGEs have distinct effects on the migration of lung alveolar epithelial cells during wound healing

After injury, cells treated with HMGB1 (0.1 µg/mL) showed higher migration at 12 h, but not at 48 h, compared to controls (Figure 3A); RAP did not influence cell migration after treatment with HMGB1 at both timepoints (Figure 3A). Cells treated with AGEs (30 µg/mL) had reduced migration at 12 h but increased migration at 48 h after injury, compared to controls. Treatment with RAP significantly counteracted these effects at both timepoints (Figure 3B and Table 3).

Both HMGB1 and AGEs increase lung alveolar epithelial cell proliferation in a RAGEdependent manner during wound healing

Treatment with HMGB1 (0.1 µg/mL) enhanced cell proliferation at both 12 h and 48 h after injury. These effects were significantly inhibited in presence of RAP at both timepoints (Figure 4A).

Treatment with AGEs (30 µg/mL) increased cell proliferation at both timepoints after injury, which was abrogated by RAP treatment at both timepoints (Figure 4B and Table 4).

Discussion

In this study, HMGB1 and AGEs promoted RAGE-dependent wound healing of lung alveolar epithelial cells after a scratch assay. In addition, both RAGE ligands increased cell proliferation in a RAGE-dependent manner. Treatment with HMGB1 increased cell migration at 12 h only, with no effect of RAP, whereas AGEs stimulated their migration in a RAGE-dependent manner 48 h after injury.

There is growing evidence supporting the involvement of RAGE pathway in mechanisms of cell and tissue repair after injury. Indeed, the AGE/RAGE axis promotes wound healing of human corneal epithelial cells [START_REF] Gross | Advanced Glycation End Products and Receptor (RAGE) Promote Wound Healing of Human Corneal Epithelial Cells[END_REF] and HMGB1/RAGE axis promotes keratinocyte healing in vitro [START_REF] Ranzato | HMGb1 promotes scratch wound closure of HaCaT keratinocytes via ERK1/2 activation[END_REF], with improved healing of diabetic human and mouse skin [START_REF] Straino | High-mobility group box 1 protein in human and murine skin: involvement in wound healing[END_REF]. In addition, HMGB1 improves scratch wound healing of human bronchial epithelial cells via toll-like receptor (TLR)-4 and RAGE-mediated increase in extracellular matrix synthesis and modulation of cell-matrix adhesion [START_REF] Ojo | High-mobility group box 1 promotes extracellular matrix synthesis and wound repair in human bronchial epithelial cells[END_REF]. However, in our dose-response experiments, higher doses of HMGB1 (>0.1 μg/mL) and AGEs (>30 μg/mL) were associated with no protective effect on wound healing, suggesting a nonlinear dose-response to increasing concentrations of alarmins [START_REF] Pittet | HMGB1 accelerates alveolar epithelial repair via an IL-1β-and αvβ6 integrin-dependent activation of TGF-β1[END_REF], such as RAGE ligands, and subsequent effects that can be either physiological or pathological to multiple cellular processes [START_REF] Xie | Cellular signalling of the receptor for advanced glycation end products (RAGE)[END_REF][START_REF] Schmidt | RAGE: a novel cellular receptor for advanced glycation end products[END_REF][START_REF] Bierhaus | Understanding RAGE, the receptor for advanced glycation end products[END_REF]. In our study, the activation of RAGE by both HMGB1 (0.1 μg/mL) and AGEs (30 μg/mL) promoted scratch wound healing whereas the use of a RAGE competitive antagonist abrogated these effects, thus supporting for the first time that previous findings of RAGE pathway as a promoter of wound healing may also hold true in the lung alveolar epithelium. Such findings should now be further explored in the clinical setting of ARDS to better understand the complex mechanisms involved in alveolar epithelium repair after injury, including the spreading, migration, and proliferation of alveolar epithelial cells [START_REF] Crosby | Balance of life and death in alveolar epithelial type II cells: proliferation, apoptosis, and the effects of cyclic stretch on wound healing[END_REF].

Here we found that both HMGB1 and AGEs promote A549 cell migration during wound healing.

Lung alveolar epithelial cells treated with HMGB1 acquired higher migration, compared to controls, at 12 h but not at 48 h. Although the significance of these discrepant findings remains uncertain, pro-migratory effects of HMGB1 may decrease with time, as previously reported by Huang et al. in another study on RAGE activation by HMGB1 [START_REF] Huang | High-mobility group box 1 impairs airway epithelial barrier function through the activation of the RAGE/ERK pathway[END_REF]. Concurrent treatment with RAP did not decrease the pro-migratory effect of HMGB1, suggesting that the activation of RAGE by HMGB1 may not be involved in epithelial cell migration during wound healing. Indeed, multiple studies support the interaction of HMGB1 with TLR-4, rather than with RAGE, as the main mechanism responsible for smooth muscle cell migration [START_REF] Yang | High mobility group box-1 induces migration of vascular smooth muscle cells via TLR4-dependent PI3K/Akt pathway activation[END_REF] or motility of human non-small lung cancer cells [START_REF] Zhu | HMGB1 induces human non-small cell lung cancer cell motility by activating integrin αvβ3/FAK through TLR4/NF-κB signaling pathway[END_REF]. In contrast, some studies suggest an important role of RAGE activation by HMGB1 in cell migration, such as in mouse endothelial progenitor cells [START_REF] Zhang | High-Mobility Group Box 1 (HMGB1) Induces Migration of Endothelial Progenitor Cell via Receptor for Advanced Glycation End-Products (RAGE)-Dependent PI3K/Akt/eNOS Signaling Pathway[END_REF], monocytes [START_REF] Vogel | Platelet-derived high-mobility group box 1 promotes recruitment and suppresses apoptosis of monocytes[END_REF], neural progenitor cells [START_REF] Meneghini | A novel role for the Receptor for Advanced Glycation End-products in neural progenitor cells derived from adult SubVentricular Zone[END_REF] or glial tumor cells [START_REF] Taguchi | Blockade of RAGE-amphoterin signalling suppresses tumour growth and metastases[END_REF]. Consistent with this hypothesis that RAGE is implicated in the process of cell migration; Huttunen et al. could inhibit transendothelial migration of melanoma cells in vivo and in vitro by blocking the HMGB1-binding site of RAGE [START_REF] Huttunen | Receptor for advanced glycation end products-binding COOH-terminal motif of amphoterin inhibits invasive migration and metastasis[END_REF], while Yang et al. reported that HMGB1 could induce migration of human dendritic cells via RAGE pathway [START_REF] Yang | High mobility group box-1 protein induces the migration and activation of human dendritic cells and acts as an alarmin[END_REF]. In addition to the RAGE-dependent beneficial effects of HMGB1 on A549 cell wound healing, we also found that treatment by AGEs (30 µg/mL) both decreased cell migration at 12 h and increased migration at 48 h, compared to controls, and that RAGE inhibition with RAP consistently counteracted these effects at both timepoints. Although future studies are needed to understand whether early pro-inflammatory effects of AGEs might explain, at least in part, decreased migration 12 h after AGEs treatment, the pro-migratory effect of AGEs at 48 h after lung alveolar epithelial cell scratch wounding is similar to previous findings in retinal pericytes [START_REF] Kim | Myricetin inhibits advanced glycation end product (AGE)-induced migration of retinal pericytes through phosphorylation of ERK1/2, FAK-1, and paxillin in vitro and in vivo[END_REF].

Cell proliferation is also an important contributor to wound healing, together with migration. Here, treatment with both HMGB1 and AGEs significantly increased cell proliferation at 12 and 48 h after scratch wounding, compared to controls, whereas simultaneous treatment with RAP counteracted these effects at both timepoints. Although this study is the first to report the importance of RAGE modulation in the proliferation of lung alveolar epithelial cells, these findings are consistent with previous in vitro studies in keratinocytes and fibroblasts [START_REF] Ranzato | HMGb1 promotes scratch wound closure of HaCaT keratinocytes via ERK1/2 activation[END_REF]. Similarly, blocking the HMGB1/RAGE axis through treatment with recombinant sRAGE, acting as a decoy receptor that prevents interaction between RAGE and its ligands, or after transfection of a truncated RAGE, was able to suppress proliferation of glial cells in mice [START_REF] Taguchi | Blockade of RAGE-amphoterin signalling suppresses tumour growth and metastases[END_REF]. We also observed that the anti-proliferative effect of RAP could be less pronounced at 12 h after treatment with HMGB1 than with AGEs (Figure 4). Indeed, a previous study found synergistic effects of RAGE and TLR activation on human bronchial epithelial wound healing, such that inhibition of both RAGE and TLR-4 was needed to abrogate the proliferative effects of HMGB1 [START_REF] Ojo | High-mobility group box 1 promotes extracellular matrix synthesis and wound repair in human bronchial epithelial cells[END_REF]. Our study is, however, unable to support that the effects of AGEs on lung alveolar epithelial cell proliferation might be more dependent on RAGE pathway than those of HMGB1. In particular, whether our use of an S100P-derived competitive antagonist for RAGE (RAP, Calbiochem® Merck Millipore, Molsheim, France) could bring different results than the use of an HMGB1-derived recombinant RAGE antagonist peptide deserves further investigation [START_REF] Lee | Production and application of HMGB1[END_REF].

This study has some limitations. First, we only used adenocarcinoma AT II-like alveolar epithelial cells as an in vitro cell model of the human lung alveolar epithelium. More translational studies of primary human epithelial cell cultures that include proinflammatory conditions and the integration of a more realistic environment (including, but not limited to, the presence of an extracellular matrix and of innate immunity cells) are now needed to further increase our understanding of the roles of RAGE pathway in mechanisms relevant to lung alveolar epithelial repair after injury, including studies of the potential roles of RAGE in cell transdifferentiation [53-55]. In the meanwhile, whether such RAGE-dependent mechanisms favoring lung epithelial wound healing could have clinical relevance in patients with acute lung injury such as ARDS also deserves further investigation. Second, we used RAP as a competitive RAGE antagonist [START_REF] Arumugam | S100P-derived RAGE antagonistic peptide reduces tumor growth and metastasis[END_REF]. However, how other strategies of RAGE inhibition, such as RAGE knockout animals, treatment with recombinant sRAGE, RNA interference or CRISPR-CAS9 methods would affect lung epithelial cell wound healing remains unknown. Finally, in this study, we did not investigate intracellular pathways (such as Mitogen-activated protein kinase ERK 1/2, JAK/STAT, JNK, Rho GTPase, PI3K/Akt, and NF-κB) through which RAGE could affect lung epithelial wound repair, and the mechanisms explaining these effects deserve further investigation [START_REF] Xie | Cellular signalling of the receptor for advanced glycation end products (RAGE)[END_REF].

In conclusion, HMGB1 and AGEs promote in vitro wound healing of human alveolar epithelial cells in a RAGE-dependent manner. These effects could be explained, at least in part, by increased cell migration and proliferation. Future studies are now warranted to better understand the precise signaling pathways through which the RAGE pathway affects lung epithelial wound repair. derived recombinant RAGE-antagonist peptide for anti-inflammatory therapy in acute lung injury, Eur. J. Pharm. Sci. 114 (2017) 275-284.

[53] H. Fehrenbach, R. Weiskirchen, M. Kasper, A.M. Gressner, Up-regulated expression of the receptor for advanced glycation end products in cultured rat hepatic stellate cells during transdifferentiation to myofibroblasts, Hepatology. 34 (2001) 943-952. and were analyzed using random effects models; when time by group interaction was significant, post-hoc comparisons between groups at each timepoint were performed using the same models. HMGB1: high-mobility group box 1 protein. AGEs: advanced glycation end-products.

. Dose-response effect of HMGB1 and AGEs on the wound healing of A549 cells. A)

Wound healing was assessed up to 48 h using a scratch assay of A549 cells treated with culture medium (Control) or with different concentrations of recombinant HMGB1 (0.1, 0.3, and 0.6 µg/mL). HMGB1 promotes wound healing, compared to the control condition within concentrations of 0.1 to 0.6 µg/mL; its effect on wound area was more pronounced with the minimum effective dose of 0.1 µg/mL, which was therefore used in subsequent experiments.

Values are reported as box plots with median values and interquartile ranges, and upper/lower limits (outliers). Measurements were performed in duplicate in six independent experiments. ****p <10 -4 . B) Wound healing was assessed up to 48 h using a scratch assay of A549 cells treated with culture medium (Control) or with different concentrations of recombinant AGEs (10, 30 and 90 µg/mL). AGEs promote wound healing, compared to the control condition at the concentration of 30 µg/mL; this minimum effective dose was therefore used in subsequent experiments. Values are reported as box plots with median values and interquartile ranges, and upper/lower limits (outliers).

[

  54] N. Demling, C. Ehrhardt, M. Kasper, M. Laue, L. Knels, E.P. Rieber, Promotion of cell adherence and spreading: a novel function of RAGE, the highly selective differentiation marker of human alveolar epithelial type I cells, Cell Tissue Res. 323 (2006) 475-488. [55] A.E. Vaughan, A.N. Brumwell, Y. Xi, J.E. Gotts, D.G. Brownfield, B. Treutlein, K. Tan, V. Tan, F.C. Liu, M.R. Looney, M.A. Matthay, J.R. Rock, H.A. Chapman, Lineage-negative progenitors mobilize to regenerate lung epithelium after major injury, Nature. 517 (2015) 621-625.

  Measurements were performed in duplicate in six independent experiments. ****p <10 -4 . HMGB1: high-mobility group box 1 protein. AGEs: advanced glycation end-products.

Figure 2 .

 2 Figure 2. The effect of RAGE modulation on lung epithelial wound healing in vitro. Wound healing was assessed in lung epithelial A549 cells treated with A) HMGB1 (0.1 µg/mL) or B) AGEs (30 µg/mL), alone or combined with RAP (12.5 µg/mL). The wound area, as measured every 12 h, is expressed as the fraction of the initial wound (100% at h0). ****p <10 -4 . RAGE: receptor for advanced glycation end-products. HMGB1: high-mobility group box 1 protein. AGEs: advanced glycation end-products. RAP: RAGE antagonist peptide.

Figure 3 .

 3 Figure 3. Effects of treatment with HMGB1 and AGEs, and of RAGE inhibition by RAP, on the migration of A549 cells during wound healing. The migration of A549 cells, as expressed as fluorescence intensity (with higher intensity corresponding to higher migration) was assessed at 12 and 48 h after a scratch assay using the CytoSelect™ 96-well cell migration assay (CBA-106; 8 μm

Acknowledgements

The authors wish to thank the staff from the Université Clermont Auvergne, CNRS, INSERM, GReD (Clermont-Ferrand, France) who helped and supported this work.

Funding

This work was supported by grants from the Auvergne Regional Council and the French Agence Nationale de la Recherche and the Direction Générale de l'Offre de Soins. RZ was funded by a PhD grant from the Ecole Doctorale des Sciences de la Vie, Santé, Agronomie, Environnement (ED SVSAE, ED 65, Université Clermont Auvergne). The funders had no influence in the study design, conduct, and analysis or in the preparation of this article.

Data availability

The datasets generated during the current study are available from the corresponding author on reasonable request.

Time postwounding (h)

Experimental conditions p value for time by group interaction p value for posthoc comparisons


Raphaël Prat

Thierry Carrard

Laurent Soulard

Olivier Durand

Raymond Namyst

Laurent Colombet

AMR-based molecular dynamics for non-uniform, highly dynamic particle simulations

Keywords: Molecular Dynamics, Micro-Jetting, AMR, MPI, OpenMP, Shockwave

à 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

Recent Molecular Dynamics (MD) simulations open up interesting perspectives in understanding the behaviour of condensed matter under shock [START_REF] Durand | Modeling from molecular dynamics simulations of ejecta production induced by shock-loaded metallic surfaces[END_REF]. This approach is complementary to and offers several advantages over the classical continuum approaches such as hydrodynamics [START_REF] Durand | Comparative simulations of microjetting using atomistic and continuous approaches in the presence of viscosity and surface tension[END_REF]. In particular, it overcomes any mesh effect while allowing any physics model but the models describing the inter-atomic interactions.

MD simulations involve N-body computations that provide the properties of a system made of particles by numerically solving the Newton's equation of motion during several timesteps: F i = m i .a i with F i the forces applied on the particle i, and m i and a i respectively the mass and the acceleration of the particle i. Forces applied on a particle are derived from the interactions with other particles modeled by a given potential F i =-∇U i , U i being the potential energy. In practice, potential is either computed by analytical formula or interpolated from tabulated values [START_REF] Wolff | Tabulated potentials in molecular dynamics simulations[END_REF].

In this paper, we mainly focus on short-range interactions, where the influence of neighboring particles is neglected beyond a given cut-off distance rcut. For simple materials such as noble gases, these interactions are modeled by pair potentials such as the Lennard-Jones potential (LJ) [START_REF] Jones | On the determination of molecular fields.II. from the equation of state of a gas[END_REF], the Morse potential [START_REF] Morse | Diatomic molecules according to the wave mechanics. II. vibrational levels[END_REF] or the Exponential-six potential [START_REF] Mason | Transport properties of gases obeying a modified buckingham (exp-six) potential[END_REF]. The Lennard-Jones potential (LJ) is usually used as a standard benchmark for MD codes due to its low computation cost. For metal materials, the Embedded Atom Model potential (EAM) or the Modified EAM (MEAM) potential [START_REF] Daw | Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals[END_REF] are more accurate, but require more computation cycles.

In theory, computing the potential for a given particle P should only involve the subset of particles included in a sphere of radius rcut centered in P . However, given that particles can move at each simulation timestep, determining such a subset for each particle P is costly. To save computing cycles, MD simulations use algorithms which approximate these subsets by considering super-sets containing more particles than those included in the rcut sphere. Two well-known techniques are used by most existing MD software: the linked cells method [START_REF] Allen | Computer simulation of liquids[END_REF], and the Verlet lists method [START_REF] Verlet | Computer "Experiments" on classical fluids. I. thermodynamical properties of lennard-jones molecules[END_REF]. The linked cells method maintains neighbor lists for each particle at each timestep with a complexity of O(n). The Verlet lists method extends the linked cells method by using a greater radius (the Verlet radius) to catch a larger neighborhood. This decreases the rate at which the lists must be refreshed.

Designing highly efficient MD simulations over parallel hardware has been a hot research topic for decades. These efforts have led to efficient and scalable software codes such as LAMMPS [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF] or Gromacs [START_REF] Berendsen | GROMACS: a message-passing parallel molecular dynamics implementation[END_REF]. Recently, these codes have evolved to cope with the evolution of supercomputers architectures towards an increasing number of cores per node, the expanding width of vector registers and the usage of non-uniform memory architectures. Current implementations of MD codes typically exploit three levels of parallelism by using a hybrid MPI+X (e.g. MPI+OpenMP) programming model together with carefully designed vector-friendly data structures. Despite these efforts, current parallelization techniques struggle with MD simulations involving large temporal variations in particle density over the simulation domain. For instance, when a material is under extreme conditions (pressure, temperature, deformation, shock, etc.) and is subject to micro-jetting or dislocation, the atom density becomes highly heterogeneous and large parts of the simulation domain suddenly become empty. In such cases, existing MD codes typically exhibit significant overhead. Moreover, their efficiency is reduced due to the sparse nature of many cells, limiting opportunities of exploiting vector capabilities of modern processors.

To compensate for such inhomogenous particle distributions, we introduce Adaptive Mesh Refinement (AMR) techniques [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF] into MD simulations. AMR is a well-known technique used in many fields to cope with such physical phenomena. It consists in dynamically adapting the mesh to better capture heterogeneity. AMR was first introduced in hydrodynamics simulations in order to reduce the computation time by refining the discretized domain (grid) where the points need to be spaced to keep their stability properties [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF]. AMR is currently used in a large amount of software [START_REF] Dubey | A survey of high level frameworks in block-structured adaptive mesh refinement packages[END_REF] for various physical problems such as astrophysics with Enzo [START_REF] Bryan | Enzo: An adaptive mesh refinement code for astrophysics[END_REF] or to solve PDE with Chombo [START_REF] Colella | Chombo software package for amr applications-design document[END_REF].

This paper introduces an efficient hybrid MPI/OpenMP parallelization taking advantage of AMR data structures to efficiently cope with highly dynamic MD simulations by reducing the imbalance between the threads and MPI processes. We leverage our previous work, based on the use of AMR techniques [START_REF] Prat | Combining task-based parallelism and adaptive mesh refinement techniques in molecular dynamics simulations[END_REF] to improve memory locality, multi-threading parallelism and we address the division of the domain simulation into subdomains assigned to the MPI processes. We demonstrate that the introduction of AMR techniques does not add overhead to communications between nodes. We also extend our task dependency graph to cope with irregular subdomain pattern. We validate our approach using an extended implementation of the ExaS-TAMP [START_REF] Cieren | ExaStamp: A parallel framework for molecular dynamics on heterogeneous clusters[END_REF] MD simulator over a cluster of 512 Intel R Xeon Phi Knight Landing (KNL). We use different input configurations and show that our approach outperforms state-of-the-art solutions while using a similar memory footprint.

The main contributions of this paper are the following:

• We introduce a new MD simulation code using AMR techniques able to efficiently exploit modern supercomputers;

• We present a task-based dependency graph strategy for AMR in order to optimize thread-based parallelization and an MPI parallelization strategy adapted to the AMR structure;

• We evaluate our implementation of the AMR using large scale productioncases and traditional benchmarks.

This paper is structured as follows. In the next Section, we briefly recall the main concepts used in MD and AMR based applications then we describe the main contribution of the presented work. We validate the relevance of our approach in Section 3 by comparing our implementation with state-ofthe-art software on several benchmarking scenarios. We present related work in Section 4, and dress concluding remarks in Section 5.

ADAPTIVE MESH REFINEMENT FOR MOLECULAR DY-NAMICS SIMULATIONS

Traditional cell-based methods can hardly cope with highly-dynamic MD simulations because they were not designed to handle huge disparities of particle density across the simulation domain. We thus propose to use AMRbased techniques to recursively subdivide the simulation domain into trees of cells, the leaves of which may be of different volumes, depending on some application-specific refinement criteria. The overall idea is to end up with small cells only where it is necessary (e.g. high density of particles), and to deal with large cells in the low density areas. We now present the most relevant features of this approach that aims at achieving high hardware efficiency as well as dynamic load balancing.

Introducing AMR in MD Simulations

As depicted in Figure 1, in our AMR-based architecture, the simulation domain is divided into subdomains using the spatial domain decomposition method [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF] and the subdomains are distributed among MPI processes. As a rule of thumb, simulations spawn one MPI process per multicore node, but other deployment schemes are possible. Unlike classic linked cells methods that split subdomains into a grid of cells, our approach splits subdomains into grids of multiple octrees called forests of octrees. Octrees are refined using an algorithm that depends on MD-specific criteria. The chosen criteria may depend on the number of particles per cell weighted by the cost of the used potential. Other criteria can be implemented according to the numerical specificities of the simulation. In 3D, each cell is recursively divided into 8 subcells of the same size, until the refinement criteria are met. The refinement must be stopped before getting the size of subcells smaller than the Verlet radius.

In each octree, the atoms' information are stored in a structure of arrays into the root cell, in order to maximize the effectiveness of compiler vectorization. Particles are sorted by their Cartesian positions (Morton index [START_REF] Morton | A computer oriented geodetic data base and a new technique in file sequencing[END_REF]) and the leaf cells correspond to contiguous sections of these arrays. During the computation of the potential energy, all the particles of an octree are processed within an OpenMP task. In other words, each octree is processed in a sequential way. Such an octree traversal strategy has proven to be very effective both on static and dynamic test cases on multicore KNL and Intel R Xeon Skylake(SKL) platforms [START_REF] Prat | Combining task-based parallelism and adaptive mesh refinement techniques in molecular dynamics simulations[END_REF].

To achieve maximal performance a tradeoff must be found between the number of octrees per MPI process (i.e. the more octrees, the more parallelism) and the size of octrees (i.e. small octrees do not exploit caches enough). Thus, at initialization time, we first evaluate the size of the small-est leaf cell, and then we try to determine a "good root cell size". This algorithm therefore searches for the optimal maximum depth of the octrees, called D max . For this purpose it uses several parameters that are mainly the Verlet radius, the minimum number of octrees over the entire domain and the maximum number of atoms per Verlet cell.

Inner-node parallelization

Inside each MPI process, we have developped two OpenMP parallelization strategies to process octrees forests over multicore hardware. The first strategy is based on a task dependency graph to limit the number of threads synchronization: octrees are processed in "waves" in such a way that neighbor octrees are never processed at the same time. This strategy is particularly well-suited when the number of octrees to process is very high with respect to the number of cores. An octree corresponds to an OpenMP task in the graph. When the number of octrees is too low, we opt for a second strategy, which uses a more synchronous approach. The overall loop that iterates on the list of octrees is parallelized using the omp for construct. Since all octrees can potentially be executed at the same time, additional thread synchronizations are needed when updating particles.

In both strategies, octrees are sequentially processed to optimize the caches reuse by a cache blocking effect and to obtain sufficient workload to be able to create tasks.

Task-based dependency graph

To take advantage of the symmetric properties of most potentials (third Newton's law: F i→j =-F j→i with F i→j the force applied by the particle i on the particle j), the forces applied between two particles are computed once for the first particle, and then reported to the second particle. This reduces the computation cost of the simulation but generates costly thread synchronizations (e.g. use of mutexes) to ensure that two particles are not updated at the same time. To address this problem, we have introduced an adaptation of the wave method [START_REF] Mangiardi | A hybrid algorithm for parallel molecular dynamics simulations[END_REF][START_REF] Prat | Combining task-based parallelism and adaptive mesh refinement techniques in molecular dynamics simulations[END_REF]. Octrees are processed in successive waves in such a way that two adjacent octrees are never processed simultaneously. This is illustrated in Figure 2, where octrees are classified on 4 waves (8 in 3D) blue, yellow, green and orange. Octrees with D max > 0 are processed wave by wave. This ensures that a particle processed by a thread is not in the neighborhood of any other particle processed by another thread. if D max = 0, there is no AMR structure and we fall back to a classical system of linked cells method. In this particular case, the number of waves may be higher, particularly because of the potential used. For example, the need to have access to neighbors' neighbours greatly increases the number of waves. It should be noted that, unlike the more traditional task blocking method [START_REF] Meyer | Efficient parallelization of short-range molecular dynamics simulations on many-core systems[END_REF], our wave method guarantees a number of waves equal to 8 if D max > 0. A straighforward approach consists in generating an OpenMP parallel loop for each wave because all tasks in the same wave can be concurrently performed without any conflict. In this case, there is always an efficiency loss at the end of each parallel loop due to the OpenMP barrier between two waves. But this method does not require to process waves in a specific order.

A solution to relax synchronizations is to overlap the execution of consecutive waves while enforcing some dependencies between neighbor octrees. This can be implemented by explicitly generating an OpenMP task graph. The construction of this graph is based on the use of dependencies between tasks with respect to their spatial location (see Figure 3). The waves are generated sequentially, but their executions can overlap. Indeed the execution of an octree belonging to the n th wave can start as soon as each adjacent octree belonging to a preceding wave has been processed.

Loop Parallelization

The task graph parallelization method exhibits high efficiency as long as the number of ready tasks is large enough to feed the cores of the underlying ), octrees classified in the wave (n) n≥1 are unlocked by specific octrees in the wave (n -1) n≥1 . In this case, a dependency between octrees 3 and 7 is added to avoid that both running at the same time.

architecture. One solution to increase the number of tasks is to artificially reduce D max . However, coping with small octrees usually has a negative impact on parallel overhead and cache usage. The second solution is to fall back to the simultaneous execution of all octrees using an OpenMP parallel loop. In this case, thread synchronizations are mandatory to preserve particle data integrity. When the number of tasks is still insufficient for both solutions, the last possibility is to iterate with OpenMP parallel loop on the leaf cells of the octrees.

Inter-node parallelization

As mentioned at the beginning of Section 2.1, the simulation domain is divided into subdomains, each subdomain being assigned to a single MPI process. There is generally a one-to-one mapping between the MPI processes and the cluster nodes. This allows our fine-grain OpenMP parallelization approach to fully tap into the potential of the underlying shared-memory architecture. Moreover, it also keeps the number of MPI processes as small as possible in order to decrease the communication cost between nodes.

Since not all neighbors of a given particle are stored locally, our approach uses the well-known ghost cells technique to grab information about missing particles. This approach is further described in the following subsection. In Section 2.3.2 we present how octrees are redistributed among MPI processes when load balancing is required. We use ghost cells around each subdomain to cache information about particles located on the border of neighbor subdomains. Such cells, which correspond to leaf cells on neighbor subdomains, do not necessarily have the same size, as illustrated on Figure 4.

Ghost Communications

Updating ghost cells requires data exchanges between nodes implemented by MPI Send/Recv communications. To speed up communications, buffers are assembled in parallel. Ghost cells are packed concurrently into send buffers (Figure 4). To avoid the usage of locks, the memory destination of each cell into a given buffer is precalculated. A similar approach is used on the receiving side to unpack ghost cells.

In contrast with the traditional approach where all the ghost cells have the same size (i.e. determined by the Verlet radius), our AMR approach transfers cells of bigger size, and thus transfers unneeded sets of particles that are actually located too far from the border of their subdomain. However, bigger-size cells do not contain large amount of particles, and transferring contiguous blocks of potentially unneeded particles has proven more effective than filtering particles to transfer only parts of some cells. To evenly distribute the workload among all MPI processes, a computational weight is evaluated for each octree based on the number of particles, their type, their depth inside the octree, the potentials used, etc. As illustrated in Figure 5, each process is in charge of a set of cells forming a contiguous space of arbitrary topological shape.

Partitioning

The possibility of using subdomains of arbitrary shape provides more flexibility when partitioning the global domain into subdomains of equivalent computational load. The set of weighted octrees is used as an input to the Zoltan [START_REF] Devine | Zoltan data management service for parallel dynamic applications[END_REF] partitioning library. The next section presents evaluations obtained using various partitioning techniques such as Shift method, Recursive Coordinate Bisection (RCB), Recursive Inertial Bisection (RIB), Space Filling Curving (SFC) and ParMetis.

EVALUATION

To evaluate the effectiveness of our approach, we present a performance study using two test cases: a stable bulk copper configuration, and a microjetting simulation. This will allow us to compare the performance of our implementation in ExaSTAMP (called AMR) with LAMMPS and with a non-AMR version of ExaSTAMP. Like LAMMPS, the non-AMR version uses a classical linked cell method combined with Verlet lists.

In this paper, we report on the performance of our MPI+OpenMP parallelization strategy, as well as on the associated load balancing mechanisms. Additional performance details about the inner-node optimizations (C++ internal structures, memory footprint, caches usage and vectorization) are presented in [START_REF] Prat | Combining task-based parallelism and adaptive mesh refinement techniques in molecular dynamics simulations[END_REF].

Experimental Framework

To validate the relevance of our approach, we evaluate the performance of our implementation against the LAMMPS package on the two selected test cases. LAMMPS is one of the most relevant classical MD simulators designed for parallel computers. It has proved to perform efficiently on benchmarks featuring several billions of atoms over tens of thousands cores. Although the standard LAMMPS version is parallelized using MPI, many packages have been developped to integrate thread parallelization such as USER-OMP or USER-INTEL. The KOKKOS package has also been integrated to optimize LAMMPS on accelerators (KNL and GPU). Load balancing can be performed by using either the shift method [START_REF] Perl | Max-min tree partitioning[END_REF] or the recursive coordinate bisection [START_REF] Berger | A partitioning strategy for nonuniform problems on multiprocessors[END_REF] (RCB) one. Despite specific optimizations for Intel KNL processors, our experiments revealed that the OpenMP version of LAMMPS remains 10 to 15% slower than the full-MPI version [START_REF] Jeffers | Intel Xeon Phi Processor High Performance Programing[END_REF]. We have notably used the following Kokkos-runtime flags recommended by LAMMPS documentation: mpirun -np nbMPI executable -k on t nbThread -sf kk -pk kokkos newton on neigh half comm no -in in.lj with newton on to use third newton law, neigh half to store a pairwise interaction in one neighbors list and comm no to use non-threaded communication.

To avoid overloading figures with uninstructive plots, we retained the hybrid version of LAMMPS in our experiments because it was the most efficient one on our target configuration.

All our experiments were run on a machine consisting of 512 KNL processors featuring 68 cores with a frequency of 1.4 Ghz. Note that 4 cores (per processor) are dedicated to for the operating system leaving 64 cores per node for user applications. Each core has a 2-wide out-of-order pipeline with support for 4 hardware threads and 512-bit vector units.

The remaining of this section is organized as follows. The first phase of our experiments consists in evaluating the performance of the two versions of ExaSTAMP against LAMMPS. We use two different test cases: the simulation of a bulk of copper (Subsection 3.2), which is a steady state configuration, and the simulation of a micro-jetting of tin (Subsection 3.3) which represents a more dynamic configuration. Both test cases use the Lennard Jones (LJ) potential. This potential is widely used to benchmark MD codes. Due to its low computational cost, the LJ potential highlights performance issues of non-computational parts of the code such as theards synchronizations, memory accesses or MPI communications.

In a second phase, we demonstrate the benefits of our approach on a more realistic test case featuring the Modified Embedded Atom Model (MEAM) potential on the micro-jetting configuration (Subsection 3.4). Although MEAM is rarely used in benchmarks due to its very high computational cost, it is one of the best potentials regarding the modelling of interactions between metallic atoms. Unfortunately, LAMMPS does not include an analytical implementation of MEAM, so only comparisons between our two versions of ExaSTAMP are reported.

Bulk of Copper

The first benchmark illustrated in Figure 6 consists in simulating 100 timesteps (∆t = 1fs) of a bulk material featuring 2.048 and 16.348 billion of atoms in thermodynamic equilibrium. This simulation has many properties such as an homogeneous density, which also remains constant throughout the whole duration of the simulation. The bulk of copper is obtained by replicating an atom pattern named lattice. A Face Centered Cubic (FCC) lattice of 3.54 Å is used with a suitable LJ potential ( = 9.340 x 10 -20 J and σ = 2.27 Å). For each atom, 39 interactions with neighboring atoms are processed due to the 5.68 Å cut-off and the symmetric property of the LJ potential. A Verlet radius of 0.3 Å is added to avoid rebuild Verlet lists. Such a homogeneous case is not a configuration able to demonstrate the advantages of AMR in MD simulations, and we merely expect ExaSTAMP AMR to exhibit decent performances thanks to its carefully optimized code design. The simulation runs on a cluster of 512 KNL (32, 768 cores) and the simulation domain is decomposed into 512 subdomains assigned to 512 MPI processes with 256 threads each (4 hyperthreads per core). Throughout this simulation, the tree depth D max is set to 2, which leads to the best results.

Overall performance numbers are reported in Most of the overall simulation time is spent inside the potential kernel. Although the number of flops computed by each version is identical, the implementation of LJ in the AMR version gives better results with 3.78 times faster than LAMMPS and 3.28 times faster than ExaSTAMP. The superiority of the AMR version on this configuration is mainly due to the avoidance of costly thread synchronizations thanks to our task-based wave execution model (see Section 2.2.1). Another reason is the optimization of the L1 cache by grouping the cells into octrees (cache blocking and Z-curve), which is only implemented in ExaSTAMP AMR. Regarding the average time consumed in updating the ghost areas, we observe that our strategy of exchanging cells of arbitrary size is not costly: the AMR and LAMMPS versions exhibit the same performance. Note that ExaSTAMP, which uses a MPI + OpenMP strategy to update ghost areas, is much slower.

Although the AMR version involves a set of techniques and algorithms specifically designed to deal with highly dynamic and heterogenous simulations, we observe that its implementation is nonetheless efficient on a steady, uniform and memory-bound test case such as Bulk of Copper.

Micro-Jetting using a Lennard Jones potential

The second test case aims at evaluating the behavior of our approach on a more dynamic and non-uniform test case: the simulation of a system generating micro-jetting.

The micro-jetting results from the reflection of a strong shock wave on a surface containing geometrical defects. The micro-jetting involves microscopic processes (break out) that are important to understand and to model for in hydrodynamics codes, for example in inertial confinement fusion simulations.

In this test case, the micro-jetting results from a shock on a tin material of 589, 228, 457 atoms, as illustrated in Figure 8, which shows the position of atoms after 1, 550, 000 timesteps of ∆t = 1fs computed using a MEAM potential. In order to make comparisons with LAMMPS, all codes use the same atom positions at timesteps 1, 550, 000 (by reading the same dump file) and use a LJ potential instead of the MEAM potential to study the impact of the non-uniform distribution of atoms. The LJ potential uses = 3.003 x 10 -20 J and σ = 2.99 Å for tin with a 10 Å cut-off. The Verlet radius is set to 1 Å for Verlet lists. The Verlet lists are rebuilt almost every 17 timesteps. In addition, the refinement algorithm is applied to each rebuilt Verlet lists.

In this test case, a high number of atoms move from one subdomain to another, increasing the volume of MPI communications. Moreover, in the AMR version, cells no longer have the same computation weight nor the same size, which challenges the underlying task scheduling system. To obtain reliable performance numbers, we run several times a simulation of 2, 000 timesteps of 1fs and the raw numbers are gathered and extracted every 100 time steps. All of the three simulations use the efficient partitioning method (Recursive Coordinate Bisection) each 500 timesteps to redistribute the octrees or cells between the subdomains.

In Figure 9, we can see that ExaSTAMP AMR performs overall 1.38 times faster than LAMMPS and 3.51 times faster than ExaSTAMP. To further understand this behavior, we report the amount of time spent in potential computation, ghost cells exchanges and Verlet lists updates. The average and maximum (scratched) time of these parts are also presented. Regarding the potential computation, we observe that AMR method clearly outperforms both the classical versions of ExaSTAMP and LAMMPS. We also note that the cost of updating octree cells in our AMR implementation is very low. As evidenced by the average time to build Verlet lists, the use of AMR techniques slightly degrades the neighbor search by introducing additional checks. The leaf cells being bigger than the cells of LAMMPS and ExaSTAMP, the number of distances checked is indeed higher.

The most expensive phase of the simulation comes from the updating ghosts cell areas. Unlike the "classic" ExaSTAMP and AMR ExaSTAMP, LAMMPS does not provide an option to desactivate the overlapping of communications, so it was necessary to find a method to determine the cost of communications in this code. To obtain a comparable communication time, we measure the time of the potential computation that we substract from the time of the global function in order to obtain the non-overlapped communication time. LAMMPS and AMR times being very close, we can conclude that the complexity of the exchanges of ghost octree cells, which heavily relies on the subdomains shape, does not degrade the overall performance of ExaS-TAMP AMR. The relative high cost of MPI communications in LAMMPS is explained by the low overlapping of these communications in the presence of a low computational potential (LJ). In addition, the AMR version is faster than ExaSTAMP because the reduction of the number of empty cells greatly reduces the time spent in the ghost cell query when building MPI messages. LAMMPS is the fastest regarding this part even if the cost of rebuilding Verlet lists has to be weighted against the total cost of the simulation, which remains low even with an inexpensive potential.

The "Other" section is the sum of the times required to compute and display logs and to implement load balancing. The AMR method has a positive impact on the time of the load balancing phase due to the low number of octrees in comparison with the number of cells.

Overall, this test case shows that ExaSTAMP AMR provides better performances than LAMMPS and confirms the relevance of introducing adaptive mesh refinement techniques in MD simulations to tackle dynamic and nonuniform configurations.

Note that our implementation could be further enhanced by allowing different D max to be used across subdomains. Currently, the number of octrees per subdomain can be highly variable, which can generate some load imbalance between MPI processes.

Micro-Jetting using a MEAM potential

We now consider using our approach on a micro-jetting production data set that requires to use the MEAM potential instead of the low cost LJ potential. The arithmetical complexity of the MEAM potential is almost 100 times greater than the complexity of the LJ potential: computing MEAM on a particle involves neighbors of its direct neighbors, and takes the possible screening of particles into account.

In our test case, the rcut is set to 4.17 Å to reduce the number of interactions. The Verlet radius remains set to 1 Å. As previously stated, LAMMPS does not provide an analytical MEAM potential and only comes with a tabulated one, whose accuracy does not meet the requirements of the simulation.

Therefore we only compare the AMR and "classic" versions of ExaSTAMP. The experimental setup is identical to the one used in previous experiments. Performances achieved by both ExaSTAMP versions are reported in Figure 10. We can observe that the overall performance of ExaSTAMP AMR is much higher: AMR is able to process 0.657x10 9 atoms per second while the "classic" version of ExaSTAMP is limited to 0.322x10 9 atom/s. The AMR version is most notably able to update ghosts more efficiently. Additionally, when the number of octree cells is too small, the AMR version is able to automatically switch to the traditional parallel strategy, which exploits inner-cell parallelism. The time spent in load balancing is lower for AMR because the algorithm has fewer objects to balance: ExaSTAMP has to deal with 81, 180, 164 cells while AMR only has to cope with 10, 087, 476 octree cells.

Note that the reconstruction of Verlet lists is more expensive in the AMR version. However, the difference with the "classic" method remains negligible and does not impact the overall performance of the code.

Finally, we have successfully used ExaSTAMP AMR to perform a complete simulation of a microjetting composed of 1.3 billion of tin atoms using MEAM (Figure 11). The simulation took almost 17 days to complete 2 millions of iterations over 512 KNLs with one post-processing output every 1000 iterations.

RELATED WORK

Many research efforts have been devoted to designing efficient MD simulations on supercomputers. LAMMPS and its associated miniMD mini application are probably the most widely used software in this domain, and have already been detailed previously. In the following, we carry on a short survey of significant MD software. dl poly 4 [START_REF] Smith | A short description of DL POLY[END_REF] brings a package of subroutines, programs and data files, written in fortran 90 and designed to make MD simulations easier and is used for production. This software is designed to perform on CPU with OpenMP and MPI, as well as on GPU. In addition there are no load balancing algorithms in dl poly 4 to compensate for a bad density distribution.

Gromacs [START_REF] Berendsen | GROMACS: a message-passing parallel molecular dynamics implementation[END_REF] is written in C with intrinsics functions for SIMD instructions and CUDA for GPU. As LAMMPS, Gromacs contains algorithms for neighbor search based on Verlet lists. Whereas it is designed for biochemical molecules like proteins, lipids and nucleic acids that have a lot of bonded interactions, an all-atom version provides short-range interactions such as Lennard Jones or Buckingham potentials. To balance the workload between MPI processes, Gromacs uses [START_REF] Hess | Gromacs 4: algorithms for highly efficient, load-balanced, and scalable molecular simulation[END_REF] the Eighth Shell Method [START_REF] Liem | Molecular dynamics simulations on distributed memory machines[END_REF].

As miniMD, several simplified versions of complex MD codes have recently been developed to serve as a testbed for various optimizations over manycore architectures. The CoMD [START_REF] Mohd-Yusof | Codesign molecular dynamics (CoMD) proxy app deep dive[END_REF] software was designed to study the dynamical properties of various liquids and solids in the context of the Ex-MatEx project. The reference implementation of CoMD uses OpenMP, MPI, OpenCL, and can use either neighbor lists or cell-lists methods.

A solution used in many other fields is to dynamically and locally adapts the mesh according to the physical phenomenon. The first adaptive mesh refinement technique [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF] has been introduced for hydrodynamics applications in order to reduce the computation time by refining the discretized domain (grid) where the points need to be spaced enough to keep their stability properties. The refinement criteria are based on the density or on a Richardson extrapolation [START_REF] Richardson | The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam[END_REF]. The AMR has been declined for N-body simulations by recursively subdividing a space into equal cells as long as the refinement criteria are fullfiled. Bodies are stored in a k-d-tree data layout, which is a special case of the binary space partitioning (BSP) trees [START_REF] Fuchs | On visible surface generation by a priori tree structures[END_REF]. It is used to find the nearest neighbors of a body or to improve the collision detection between two bodies. AMR is currently used in a large amount of software [START_REF] Dubey | A survey of high level frameworks in block-structured adaptive mesh refinement packages[END_REF] for various physical problems such as astrophysics with Enzo [START_REF] Bryan | Enzo: An adaptive mesh refinement code for astrophysics[END_REF] or to solve PDEs with BoxLib [32] and Chombo [START_REF] Colella | Chombo software package for amr applications-design document[END_REF]. To the best of our knowledge, our approach is the first to introduce AMR in molecular dynamics simulations.

CONCLUSION

Performing MD simulations over supercomputers is a major step towards a better understanding of physics in many scientific fields. In this work, we propose an efficient MPI, OpenMP and vectorization parallelization using the advantages of Adaptive Mesh Refinement structures to solve nonuniform large scale MD simulations by reducing the imbalance inside and between multicore cluster nodes. The octree-based strategy enables effective optimizations by reducing the number of empty cells, improving data reuse by cache blocking, and by using a lock-free task-based algorithm. Our experiments show that our implementation outperforms the LAMMPS state-of-art implementation by a factor of 1.38 on a micro-jetting and by a factor of 2.61 on a steady scenario.

Future works include auto-tuning to find the optimal maximum depth tree of octrees per subdomain, so as to reach the best tradeoff between the amount of parallelism and the efficiency of per-thread code sections. We are also investigating improvements to our load balancing algorithm to automatically identify the most appropriate triggering points. In order to deal with simulations with long-range interactions, one solution could be to combine our AMR structure with the octree-based architecture used by the fast mul-tipole method [START_REF] Greengard | A fast algorithm for particle simulations[END_REF]. In addition, a solution, proposed by Agullo et al [START_REF] Agullo | Task-based fmm for multicore architectures[END_REF], included a task-based parallelism with StarPU [START_REF] Augonnet | Starpu: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] runtime that matches with our task-based dependency graph. Finally, we intend to further validate our approach on new highly-dynamic configurations.

Figure 1 :

 1 Figure 1: MD software architecture adapted for AMR (2D). Particles are distributed over a number of subdomains. Each subdomain is in turn decomposed into octree cells.

Figure 2 :

 2 Figure 2: 2D representation of a subdomain of the MD simulation. The atoms are distributed into the octrees themselves included into octrees with D max = 1.

Figure 3 :

 3 Figure 3: Task dependency graph in 2D. Octrees are classified into 4 waves (8 in 3D), octrees classified in the wave (n) n≥1 are unlocked by specific octrees in the wave (n -1) n≥1 . In this case, a dependency between octrees 3 and 7 is added to avoid that both running at the same time.

Figure 4 :

 4 Figure 4: Update of a portion of the ghost octrees of domain B by domain A. In order to minimize particle transfer, only the leaf cells of the data on the border of the domain A are selected (green).

Figure 5 :

 5 Figure 5: The simulation domain is partitioned into different subdomains that can be of any topological shape.

Figure 6 :

 6 Figure 6: Bulk of coppper with 16,384,000,000 atoms of copper generated with 1600*1600*1600 FCC lattices.

Figure 7

 7 Figure 7 describes the average time spent during potential computation,

Figure 7 :

 7 Figure7: Distribution of costs between potential computation and message exchanges between sub-domains (Lower is better).

Figure 8 :

 8 Figure 8: Micro-jetting of 589, 228, 457 atoms of tin obtained with an MEAM potential after 1, 550, 000 iteration of time (∆t = 1fs).

Figure 9 :

 9 Figure 9: Time for different sections of each software during the micro-jetting simulation. For the sake of clarity, times over 200 seconds are reported on the left and the others on the right. The mean (no stripe) and maximum (stripes) time of MPI processes are given for LJ potential computation, ghost update and neighbors update. (Lower is better).

Figure 10 :

 10 Figure 10: Time for differents sections of the micro-jetting simulation with MEAM: Total (without init), MEAM potential computation, ghost udates, Neighbors update, refinement and coarsening for AMR, dynamic load balancing and IO. (Lower is better).

Figure 11 :

 11 Figure 11: Micro-jetting simulation : the shock, which propagates through the crystal, causes a phase change. This results in the ejection of a liquid tin sheet with a discontinuous edge.

Table

Table 1 :

 1 . The measurements show that ExaSTAMP AMR is 2.61 times faster than LAMMPS and 4.05/3.48 times faster than ExaSTAMP. The bulk configuration with 16.348 billion of atoms can not be generated LAMMPS due to the limitation of 32-bit integers during the velocity initialization loop, even with the BIGBIG option. Bulk of copper: number of atoms processed per second (Higher is better).

	Software	2.048 x 10 9 atoms	16.384 x 10 9 atoms
	ExaSTAMP	4.0 x 10 9	6.5 x 10 9
	ExaSTAMP AMR	16.5 x 10 9	22.6 x 10 9
	LAMMPS	6.2 x 10 9	NA

Acknowledgments

This work was funded by the French Programme d'Investissements d'Avenir (PIA) project SMICE.


Professor Pierre Labauge
email: labauge@yahoo.fr

Do disease modifying drugs (DMD) have a positive impact on the occurrence of secondary progressive multiple sclerosis? No

Keywords: Multiple Sclerosis, DMD, Long term follow, Will Rogers phenomen, Bias

RR MS evolution has changed since the beginning of the avaibility of MS disease-modifying drugs (DMD). Before concluding of an unique impact of the efficiency of DMD, careful analysis of long term studies has to be conducted. Analysis of the litterature points a few bias in the long term follow of MS patients under DMD: indication of DMD has changed since 20 years, diagnosis criteria are not the same (including the Will Rogers phenomen), and so far population are not homogeneous and comparable. Analysis criteria of the efficiency of the treatments are not the same, pending on the date of the publications.

References concerning the long term impact of DMD are in fact very limited. In addition, long term efficiency of 2 nd line treatments are not available. An other explanation of the change of MS evolution, could be the lower evolutivity of MS patients since 2 decades.

Analysis of placebo group in pivotal studies, argue to a decrease of the relapse annual rate and mean EDSS score in the more recent studies and recent MS diagnosed patients. To conclude, long term evolution of MS patients is more favorable, influence of DMD is likely, but not unique.

The main objectives of MS disease-modifying drugs (DMD) are to reduce the rate, the age of onset and the severity of the progression of MS, and so far the rate and severity of RRMS turning to SPMS. In the cohort of MS patients followed in Lyon's hospital [START_REF] Vukusic | Prognostic factors for progression of disability in the secondary progressive phase of multiple sclerosis[END_REF], the median time to start the progression of MS in RR MS patients was 19 years after the onset of the disease. French MS data (EDMUS) analysis in 2017, showed that the mean delay of onset of progression is 28 years since the onset of the disease. So clearly, progression of RR-MS has changed since 20 years. Can DMD efficiency explain the better prognosis of the MS ?.

One explanation of the delayed age of the apparition of the progression could be effectively explained by the long term efficiency of the DMD. But, analysis of long term DMD efficiency has to be studied with a careful attention. Scalfari in 2014, studied the correlation of the MS relapse annual rates [START_REF] Scalfari | Onset of secondary progressive phase and long-term evolution of multiple sclerosis[END_REF]. 730 RR MS patients were followed during 28 years. They observed that occurrence of relapses 2 years after the onset of the disease did not implicate a progression of the disease. Of the 730 patients, 389 presented 1 relapse in the 2 years after the onset, 183 : 2 relapses, 158 : 3 relapses. Of the 158 latter patients, 103 (65.2 %) turned to a SPMS, the remaining 55 did not present a SPMS. The authors concluded that influence of relapses in the beginning of the disease, was not related to a progressive form (SP).

In another study, the same goup followed retrospectively a large cohort of RR MS. They divided the patients in 3 groups : MS under DMD (868), contempory MS (> 1995 ; DMD are available, but patients did not receive treatments, n = 829), historic MS (no DMD available, < 1995, n = 959). The mean follow was respectively 5.1 ; 4 and 10. 8 years. The rate of EDSS score above of 6 was respectively 10.8 % ; 5.3 % ; 23.1 %. The authors suggest that long term evolution of MS patients with or without DMD was not so different. This study has some limits (indication of treatments was not precized, 400 patients were lost of follow, with unknown reasons, DMD patients were younger and had a more active MS).

Similary, Kinkel et al. [START_REF] Kinkel | Association between immediate initiation of intramuscular interferon beta-1a at the time of a clinically isolated syndrome and long-term outcomes: a 10-year follow-up of the Controlled High-Risk Avonex Multiple Sclerosis Prevention Study in Ongoing Neurological Surveillance[END_REF] compared the long term of follow of RR MS patients, treated by interferon after the onset of the first symptom (CIS) or 2 years after. At 10 years, the mean EDSS score was not statistically different between these two groups. Careful analysis of this study showed bias, including a low rate of patients with a long term follow (only 127 / 383), a low activity of MS. At last, EDSS score evaluation was a secondary endpoint, and not a primary endpoint.

In the same way, long term study of the BENEFIT trial [START_REF] Kappos | The 11-year long-term follow-up study from the randomized BENEFIT CIS trial[END_REF], did not find a difference in the mean EDSS score at year 5, in CIS patients treated by interferon, just after the onset or 2 years after. The number of treated patients to decrease of EDSS score of 1 was estimated to 25 patients. Ebers [START_REF] Ebers | Analysis of clinical outcomes according to original treatment groups 16 years after the pivotal IFNB-1b trial[END_REF] followed the long term clinical outcome of patients treated by Interferon (16 years), without any difference in the mean EDSS score in patients under or not DMD.

In a recent study (2015), Sormani reviewed the publications concerning the long term efficiency of the MS treatments [START_REF] Sormani | Can we measure long-term treatment effects in multiple sclerosis?[END_REF]. Each study had many bias ; different selection criteria of the patients, Will Rogers phenomen [START_REF] Sormani | Will Rogers Phenomenon in Multiple Sclerosis[END_REF], and important loss of follow. Analysis criteria of long term impact of DMD prescription, was very different in each study, ranging from mean EDSS score, apparition of progression rate, or death rate.

Another explanation of the DMD efficiency could be related to the date of DMD initiation, changing with time to the diagnosis. Capra et al. [8] , found that the delay between treatment start and diagnosis, was 10.6 years during the period of diagnosis between 1980-1990, 6.4 years during 1991-1995 period, and 0.7 years for patients with a diagnosis after 2011. Some features argue to a change of the natural history of the MS since the area of the DMD, but independently of the DMD. The comparison of the relapses rate in the placebo group of the pivotal studies showed a decrease of relapses rate, since the first studies to the more recent studies. A meta analysis of the 56 pivotal studies by Steinvorth, found that the rate of relapses decreased of 4.5 % each year [START_REF] Steinvorth | Explaining temporal trends in annualized relapse rates in placebo groups of randomised controlled trials in relapsing multiple sclerosis: systematic review and meta-regression[END_REF]. Rover et al., found the same results concerning the EDSS progression.

At last, evolution of the diagnosis criteria could induce a change of the evaluation of the disease. This phenomen is called the Will Roger phenomen. Sormani et al [START_REF] Sormani | Will Rogers Phenomenon in Multiple Sclerosis[END_REF] followed 309 CIS patients , with a mean follow of 84 months. They showed that the rate of patients with an EDSS score of 3 at 7 years, was 11 % in the patients with Poser criteria, versus 7 % using McDonald criteria, MS diagnosis: 46 % versus 27 %. The conclusions of this study are that the rate of CIS turning to MS is also pending on the used MS criteria.

One other difficulty to objective the long term evolution of the MS, is the analysis criteria of progression. The MS base group published a recent study concerning the definition of the SP MS. This large study, including 17 356 MS patients, retained three SP criteria: a progression more than 3 months, a FFS score of 2 and more, an EDSS score more than 4. Other criteria could be used, as in the primary progressive form , including EDSS progression, T25FW and 9HPT score. The association of these 3 criteria were determined as NEP, for « no evidence of progression ». Others criteria could be used as, cognitive functions, fatigue, urinary dysfunctions, eventually neur ophtalmologic evaluation, but such studies are lacking.

Interests conflicts :

Dr Labauge received fees to travel organization or board participations from Novartis, Biogen, Bayer, Merck, Teva, Genzyme/Sanofi, Roche, Alexion

  In conclusion, RRMS has changed since the beginning of the availability of DMD. Unique and direct of the DMD efficiency is likely. Neurologists have to know that studies concerning long term efficiency of DMD are important, including the different indications of DMD, heterogeneous MS populations since 25 years, (including the Will Rogers phenomen) and the change of the evolution of MS independantly of the DMD. Other questions are not solved, mainly the long term of progression under DMD. At last, no long term efficiency of the 2 nd line treatments are yet available.