Eric Hachulla No

Le Gouellec

David Launay

Marie-H Élène Balquet H Élène Maillard

Raymond Azar

Amal Boldron

Pierre Bataille

Marc Lambert

Anne-Laure Buchdahl

Delphine Allorge

Pierre-Yves Hatron

Vincent Sobanski

Benjamin Hennart

Pierre Clerson

Sandrine Morell-Dubois

Eric Hachulla
email: eric.hachullla@chru-lille.fr

Noémie Le Gouellec

Marie-Hélène Balquet

Hélène Maillard

Esstim Investigators

Prof Eric Hachulla

Adherence to hydroxychloroquine in patients with systemic lupus: contrasting results and weak correlation between assessment tools

Keywords:

Adherence to hydroxychloroquine in patients with systemic lupus: contrasting results and weak correlation between assessment tools.

Highlights:

 The adherence rate to HCQ treatment in SLE varied between 3.2% and 32.5%  Correlations between blood HCQ-concentration and self-questionnaires are weak and agreement between methods was poor  Combining blood HCQ concentration with MASRI and MMAS-8 may help to better identify nonadherence in SLE

ABSTRACT

Objectives: Hydroxychloroquine (HCQ) is an anchor drug in the treatment of systemic lupus erythematosus (SLE). Adherence to HCQ is key for efficacy. Inaccurate evaluation of adherence could lead to non-justified switch to more expensive or less tolerated drugs.

Methods: Severe non-adherence rate to HCQ was estimated in a sample of SLE patients during a routine visit using blood HCQ concentration < 200 µg/L. Adherence was assessesd by the Medication Adherence Self-Report Inventory (MASRI) < 80/100, 8-item Morisky Medication Adherence Scale (MMAS-8) ≤ 6 /8, Health Care Provider (HCP) visual analog scale (VAS) < 80/100. Same procedures were to be repeated during a further routine visit 6 to 12 months later. We described agreement and correlations between tools and compared severely non-adherent patients and others on their characteristics.

Results:

The study involved 158 patients (86.1% females) aged 42.2±12.6 years treated with HCQ for 9.6±6.9 years. Blood HCQ concentration (mean±standard deviation) was 1046±662 µg/L at visit 1 and 855±577 µg/L at visit 2. At visit 1, the non-adherence rate varied from 3.2% (blood HCQ level < 200 µg/L) to 7.7% (MASRI), 12.4% (HCP-VAS) or 32.5% (MMAS-8). 37.8% of patients met at least one of the definitions of non-adherence. Patients' characteristics including SLE activity, damage and quality of life were similar between severely non-adherent patients and others. Correlations between blood HCQ-concentration and self-questionnaires were weak (r <0. [START_REF] Ting | Usefulness of cellular text messaging for improving adherence among adolescents and young adults with systemic lupus erythematosus[END_REF] and agreement between methods was poor.

Conclusion: Blood HCQ concentration < 200 µg/L reveals severe non-adherence. Combining blood HCQ concentration with MASRI and MMAS-8 may help to better identify nonadherence in SLE. Agreement between methods was poor and correlations with HCQ level and SLE activity were weak.

Introduction

Systemic lupus erythematosus (SLE) is a chronic inflammatory condition, notably comprising joint, skin, kidney, and sometimes, nervous system damages, with an estimated prevalence of 4.7 per 10,000 inhabitants in France [START_REF] Arnaud | Prevalence and incidence of systemic lupus erythematosus in france: A 2010 nation-wide population-based study[END_REF]. Hydroxychloroquine (HCQ), a 4-aminoquinolone antimalarial molecule, is the anchor drug of lupus treatment with a good safety profile and a favorable cost: effectiveness ratio. HCQ prevents flares [START_REF] Hydroxychloroquine Study | A randomized study of the effect of withdrawing hydroxychloroquine sulfate in systemic lupus erythematosus[END_REF], improves cutaneous manifestations [START_REF] Kuhn | Cutaneous lupus erythematosus: Update of therapeutic options part i[END_REF] and arthritis [START_REF] Williams | Comparison of hydroxychloroquine and placebo in the treatment of the arthropathy of mild systemic lupus erythematosus[END_REF], enhances survival [START_REF] Alarcon | Effect of hydroxychloroquine on the survival of patients with systemic lupus erythematosus: Data from lumina, a multiethnic us cohort (lumina l)[END_REF][START_REF] Hsu | Adherence to hydroxychloroquine improves long-term survival of patients with systemic lupus erythematosus[END_REF], and decreases organ damages [START_REF] Fessler | Systemic lupus erythematosus in three ethnic groups: Xvi. Association of hydroxychloroquine use with reduced risk of damage accrual[END_REF]. HCQ may be prescribed as a monotherapy or be combined with nonsteroidal anti-inflammatory drugs, corticosteroids, immunosuppressive drugs or biologics.

As observance to HCQ is key for efficacy, an accurate diagnosis of poor adherence may prevent inappropriate switch to more expensive or less tolerated drugs. Non-adherence to HCQ regimen increases health resource consumption [START_REF] Julian | Depression, medication adherence, and service utilization in systemic lupus erythematosus[END_REF] and the risk of subsequent acute care use [START_REF] Feldman | Medication nonadherence is associated with increased subsequent acute care utilization among medicaid beneficiaries with systemic lupus erythematosus[END_REF]. The determination of blood HCQ concentration is an objective and direct measure of the patient's adherence. A cutoff of 200 µg/L has been proposed to define non-adherence [START_REF] Costedoat-Chalumeau | Very low blood hydroxychloroquine concentration as an objective marker of poor adherence to treatment of systemic lupus erythematosus[END_REF]. Based on this threshold, the proportion of SLE patients who are poorly adherent appears to vary from 10% [START_REF] Frances | Low blood concentration of hydroxychloroquine in patients with refractory cutaneous lupus erythematosus: A french multicenter prospective study[END_REF] to 22% [START_REF] Costedoat-Chalumeau | A prospective international study on adherence to treatment in 305 patients with flaring sle: Assessment by drug levels and selfadministered questionnaires[END_REF]. Several other tools are available to assess patients' adherence. Indirect methods include medication event monitoring systems (MEMS), pharmacy records, health care provider assessment, and self-report.

Taking advantage of a regional network managing around 1500 SLE patients, we set up a study which aimed at estimating the prevalence of non-adherence using various tools including determination of HCQ blood concentration, use of well-validated questionnaires, and health care provider (HCP) assessment. We were thus able to J o u r n a l P r e -p r o o f estimate various rates of non-adherence to HCQ, to describe the characteristics of SLE patients according to their adherence to HCQ and to correlate different tools routinely used for the assessment of adherence.

Material and methods

Trial design, setting and ethics

This was a prospective, longitudinal, observational study. The study protocol was approved by a locally appointed ethics committee (Comité de Protection des Personnes -Hauts de France). Each patient was informed of the aims and constraints of the study and gave an informed consent prior to any study procedure. The study complied with the Declaration of Helsinski. The study population was composed of adult patients with a diagnosis of SLE meeting the American College of Rheumatology (ACR) criteria [START_REF] Hochberg | Updating the american college of rheumatology revised criteria for the classification of systemic lupus erythematosus[END_REF], treated with HCQ for at least 3 months at a dose  200 mg/day. Pregnant and breastfeeding women were excluded from the study. Data were collected during two routine visits at a planned interval of 6 to 12 months. SLE activity was measured using the Systemic Lupus Erythematosus Disease Activity Index (SELENA-SLEDAI) ranging from 0 to 105 [START_REF] Bombardier | Derivation of the sledai. A disease activity index for lupus patients. The committee on prognosis studies in sle[END_REF]. Damages from the diagnosis of SLE onwards were evaluated by the Systemic Lupus International Collaborating Clinics (SLICC) damage index [START_REF] Gladman | The development and initial validation of the systemic lupus international collaborating clinics/american college of rheumatology damage index for systemic lupus erythematosus[END_REF]. At both study visits, a blood sample was drawn from the patient, stored at ambient temperature and blood HCQ concentration was measured within 48 h (Toxicology Laboratory, Biology-Pathology Center, Lille) using liquid chromatography coupled to mass spectrometry detection (LC-MS). A blood HCQ concentration < 200 µg/L defined non-adherence [START_REF] Costedoat-Chalumeau | Very low blood hydroxychloroquine concentration as an objective marker of poor adherence to treatment of systemic lupus erythematosus[END_REF]. The adherence was also evaluated by the Medication Adherence Self-Report Inventory (MASRI) questionnaire [START_REF] Walsh | Responses to a 1 month self-report on adherence to antiretroviral therapy are consistent with electronic data and virological treatment outcome[END_REF], a self-administered J o u r n a l P r e -p r o o f questionnaire already used in SLE patients [START_REF] Koneru | Effectively measuring adherence to medications for systemic lupus erythematosus in a clinical setting[END_REF] with a cutoff of 80 according to the authors recommendations, and by the 8-item Morisky Medication Adherence Scale (MMAS-8) [START_REF] Morisky | Predictive validity of a medication adherence measure in an outpatient setting[END_REF][START_REF] Krousel-Wood | New medication adherence scale versus pharmacy fill rates in seniors with hypertension[END_REF][START_REF] Morisky | Improving the measurement of self-reported medication nonadherence: Response to authors[END_REF][START_REF] Berlowitz | Effect of intensive blood-pressure treatment on patient-reported outcomes[END_REF]. MASRI tries to quantify the number of missed pills and the overall adherence over the last 30 days. MMAS-8 explores the patient's behaviours regarding the treatment. Patients were considered to be poorly adherent if the MMAS-8 score was ≤ 6 [START_REF] Morisky | Predictive validity of a medication adherence measure in an outpatient setting[END_REF]. MMAS-8 allows the intentional and unintentional non-adherence to be distinguished. Intentional non-adherence refers to non-adherence that is deliberate and largely associated with patient motivation whereas unintentional non-adherence is non-adherence that is largely driven by a lack of capacity or resources to take medications (for instance, when patient cannot have access to the treatment because of its cost). In addition, the reasons underlying intentional and unintentional non-adherence are not entirely independent in that certain types of unintentional non-adherence e.g. forgetting, are logically more likely when motivation for medication is low. Patients were asked to fill in a MMAS-8 scale regarding HCQ and another one regarding other treatments. The Health Care Provider (HCP) assessed the patient adherence on a 100-mm VAS. A convenient value < 80 was used to define poor adherence as patients with compliance of at least 80% are generally considered to be compliant [START_REF] Sabaté | Improving adherence rates: Guidances for countries. In adherence to long-term therapies: Evidence for action[END_REF] . In addition, patients were asked to self-rate their adherence to the HCQ regimen and the benefit brought by the treatment on two 100-mm VAS. Quality of life was assessed with the French version of the LupusQoL scale [START_REF] Devilliers | Lupusqol-fr is valid to assess quality of life in patients with systemic lupus erythematosus[END_REF]. Anxiety and depression were explored by the HAD-Anxiety and HAD-depression scale, respectively.

Statistics

Results are expressed as mean±standard deviation or number (percentage). Percentages are calculated on observed data. The proportion of non-adherent patients was estimated with a two-sided 95% confidence interval (CI) at each visit assuming a binomial J o u r n a l P r e -p r o o f distribution. Adherent and non-adherent patients according to blood HCQ concentration at visit 1 and/or visit 2 were compared using parametric (Student's t test) and nonparametric (Mann-Whitney's U test, chi squared test, Fisher's exact test) statistical tests depending on the nature and distribution of data. The type 1 error risk was set at 5% and the tests were two-sided. The comparative analysis was conducted on patients with documented blood HCQ concentrations at both visits. Correlations between assessment tools were quantified by the Pearson's method.

Results

Population

The study included 158 patients. Blood HCQ concentration was determined in 156 patients at visit 1, in 146 patients at visit 2 and in 145 patients at both visits, spaced by 12.0±4.2 months (blood HCQ measurements were not performed in 2 patients at visit 1 and in 12 patients at visit 2 for technical reasons). Mean age was 42.2±12.6 years and 136 patients (86.1%) were females. Almost all documented patients were positive for antinuclear antibodies (152/153), and 26 (16.5%) patients had an associated antiphospholipid syndrome. Twenty-one patients (13.3%) had an associated Sjögren syndrome and 6 (3.8%) patients had a diagnosis of fibromyalgia. The SLICC damage index was 0 for 126 (81.3%) patients, 1 for 20 patients (12.9%), 2 for 8 (5.2%) patients, and 5 for one patient (0.7%). SLICC damage index was missing in 3 patients. At entry in the study, the mean SELENA-SLEDAI score was 2.7±3.7.

Patients were treated with HCQ for 9.6±6.9 years at a daily dose of 297.1±97.9 mg and 304.9±93.6 mg at visit 1 and 2, respectively. Concomitant treatments included corticosteroids (N = 86; 54.4%), immunosuppressive drugs (N = 42; 26.6%) and biologics J o u r n a l P r e -p r o o f (N = 3; 1.9%). Patients rated the benefit brought by the therapy at 70.7±22.8 and estimated their adherence to HCQ to be 88.9±17.5 on 100 mm VAS. 123 (77.9%) patients had undergone at least one blood HCQ determination in the past.

Non-adherence rates

Blood HCQ concentration (m±SD) was 1046±662 µg/L and 855±577 µg/L at visit 1 and 2, respectively (Figure 1). Two patients had undetectable blood HCQ concentration at visit 1. There were no patients with undetectable HCQ concentration at visit 2. Based on blood HCQ concentration < 200 µg/L, non-adherence rate varied over time and poor adherence was more frequent at visit 2. Five (3.2%) patients were poorly adherent at visit 1 and 14 (9.6%) at visit 2. Only one patient was poorly adherent at both visits and 18 (12.4%) patients were poorly adherent on at least one visit (Table 1). The blood HCQ concentration was not significantly different between patients who have had at least one dosage in the past (1060±654 µg/L) and those who had none (994±697 µg/L) (p=0.61).

The blood HCQ concentration was higher in patients taking more than one HCQ pill daily (1119±617 µg/L versus 917±718 µg/L) (p=0.08) and 3 out the 5 non-adherent patients identified at visit 1 were taking one pill per day.

At visit 1, the mean MMAS-8 score was 6.6±1.7 and 50 (32.5%) patients had a MMAS-8 ≤ 6. MMAS-8 score was 3.4±1.0 for intentional non-adherence and 3.3±0.9 for unintentional non-adherence. Among the 50 patients non-adherent according to the MMAS, poor adherence was intentional in 16 (10.8%) patients, unintentional in 27 (18.2%) and undetermined in 7 patients. Among the 45 patients with MAAS-8 documented for HCQ and other treatments and who were not adherent to HCQ, 34 were also poorly adherent to other medications. The mean scores for intentional and J o u r n a l P r e -p r o o f unintentional non-adherence were similar for HCQ and for other medications. Fortyseven (33.6%) patients were not adherent at visit 2 and 66 (45.5%) patients were considered as non-adherent on at least one visit. The MASRI score was 92.7±11.8 at visit 1 and 12 patients (7.7%) had a MASRI < 80. Poor adherence rate was 10.7% at visit 2 and 20 (13.8%) patients had a MASRI <80 on at least one visit. HCP-VAS was used only at visit 1 and the mean score was 87.8±12.0. Nineteen (12.4%) patients were rated less than 80.

Depending on the assessment tools and on their combinations, non-adherence rates varied from 3.2% to 37.8% at visit 1, from 9.6% to 33.4% at visit 2 and from 12.4% to 49.0% when considering both visits (table 1).

Characteristics of patients according to adherence defined by various methods

Non-adherent patients at visit 1 and/or visit 2 according to blood HCQ concentration were compared to adherent patients (Table 2). The analysis involved the 145 patients documented for blood HCQ concentration at both visits. Non-adherent patients (N=18) had no specific characteristics when compared to adherent patients. They presented more frequently an associated antiphospholipid syndrome (33.3% versus 11.8%; p=0.03) and their MASRI and MMAS-8 scores were significantly lower (p=0.01 and p=0.02, respectively). In contrast, physicians did not differentiate poor and good compliers (p=0.37). Adherence had no impact on the patients' quality of life (figure 2). Non-adherent patients defined according to MMAS-8 differed from adherent patients on some points. They were younger (p=0.0002), had a lower number of pills to take daily (p=0.03), a lower MASRI score (p<0.0001), as well as a lower score on HCP-VAS (p=0.004). Using MASRI < 80 for the definition of poor adherence, non-adherent patients were younger (p=0.04), had a lower self-assessed adherence (p=0.0006), a lower MMAS-J o u r n a l P r e -p r o o f 8 score (p<0.0001) and a lower HCP-VAS (p=0.01). Poor compliers defined by HCP-VAS < 80 had also, albeit not significant, lower MASRI and MMAS-8 scores.

Correlations and agreement between tools

Correlations between the different adherence assessment tools were weak to moderate (Table 3). MMAS-8 and MASRI were well-correlated (r=0.81 at visit 1, r=0.56 at visit 2).

Correlations between low blood HCQ concentration and other assessment tools at visit 1 did not exceed 0.20. There was no agreement between tools on the definition of nonadherent patients. Taking blood HCQ concentration< 200 µg/L as gold standard, sensitivity and specificity were 80.0% and 94.7% for MASRI < 80 at visit 1, and 38.5% and 92.1% at visit 2. Sensitivity and specificity for MMAS-8 ≤ 6 were 100.0% and 69.8% at visit 1, 42.9% and 67.5% at visit 2. MMAS-8 < 6 had a high sensitivity to predict MASRI < 80 (100.0% at visit 1 and 92.9% at visit 2) and a moderate specificity (73.0% at visit 1 and 74.8% at visit 2). Conversely, MASRI < 80 had a low sensitivity to predict MMAS-8 < 6 (24.0% at visit 1 and 29.5% at visit 2), but a high specificity (100.0% and 98.9% at visit 1 and visit 2). Blood HCQ < 200 µg/L displayed high sensitivity, but poor sensitivity, to predict MASRI < 80 or MMAS-8 < 6 independently of the visit (Table 4).

Analysis of an alternative cutoff for blood HCQ concentration

A HCQ blood cutoff of 500 µg/L has also been proposed [START_REF] Durcan | Hydroxychloroquine blood levels in systemic lupus erythematosus: Clarifying dosing controversies and improving adherence[END_REF]. Using this threshold increased the proportion of non-adherent patients at visit 1 (N=28; 19.3%), but did not improve the agreement between assessment methods: sensitivity was 14.3% for MASRI < 80 and 39.3% for MMAS-8 < 6; specificity was 93.7% and 69.0%, respectively.

Discussion

J o u r n a l P r e -p r o o f

Inaccurate assessment of adherence to HCQ of SLE patients may lead to the inappropriate prescription of more expensive or less tolerated drugs. There is a need for reliable assessment tools that can be used in practice to assess adherence. Subjective and objective measures have been proposed but are poorly correlated with each other.

Using blood HCQ determination and a cutoff of 200 g/L, we found that the nonadherence rate was around 10%. In a study involving SLE flaring patients, Costedoat-Chalumeau [START_REF] Costedoat-Chalumeau | A prospective international study on adherence to treatment in 305 patients with flaring sle: Assessment by drug levels and selfadministered questionnaires[END_REF] reported that HCQ levels were < 200 g/L in 14.5% of patients, and undetectable in 7.2%. In our study including a majority of non-flaring patients, the proportion of poor adherence was similar than those (7%) previously reported in a nonflaring SLE population [START_REF] Costedoat-Chalumeau | Very low blood hydroxychloroquine concentration as an objective marker of poor adherence to treatment of systemic lupus erythematosus[END_REF]. In contrast, and despite the use of a less stringent HCQ cutoff at 100 g/L, some authors [START_REF] Ting | Usefulness of cellular text messaging for improving adherence among adolescents and young adults with systemic lupus erythematosus[END_REF][START_REF] Iudici | Health status and concomitant prescription of immunosuppressants are risk factors for hydroxychloroquine non-adherence in systemic lupus patients with prolonged inactive disease[END_REF] reported a poor adherence rate of 29%. In a population of 70 patients with childhood-onset SLE, Ting et al found that only 29% were adherent to HCQ treatment when adherence was based on undetectable HCQ levels (< 0.1 g/L) [START_REF] Ting | Usefulness of cellular text messaging for improving adherence among adolescents and young adults with systemic lupus erythematosus[END_REF]. The educational program set up in our department for more than 10 years, individual and collective information delivered to the patients, free access to internet website of French sector on Auto-immune and Auto-inflammatory diseases (https://www.fai2r.org) may explain our apparent good results. Nevertheless, this deserves further studies.

On the 18 patients who were diagnosed as being non-adherent based on blood HCQ concentration, only one was poorly adherent at both visits. In a group-based trajectory model, among 10,406 HCQ initiators with SLE in the US Medicaid database, 17% were persistent compliers, 36% persistent non-compliers, and 47% formed two dynamic patterns of partial adherence [START_REF] Feldman | Dynamic patterns and predictors of hydroxychloroquine nonadherence among medicaid beneficiaries with systemic lupus erythematosus[END_REF]. In our study, the rate of non-adherence to HCQ varied between visit 1 and visit 2, 12 months apart, and the rate of non-adherent patients was J o u r n a l P r e -p r o o f higher at visit 2 than at visit 1. We cannot easily explain this counter-intuitive finding since visits 1 and 2 were routine visits and neither the HCQ daily dose nor the SELENA-SLEDAI score varied between visits. Moreover, patients were aware that a new determination of blood HCQ concentration would be performed at visit 2 and it could be expected that a history of recent HCQ determination would increase adherence. We did not find a better adherence in patients who have had a previous measure of HCQ blood concentration as has been reported by Durcan et al [START_REF] Durcan | Hydroxychloroquine blood levels in systemic lupus erythematosus: Clarifying dosing controversies and improving adherence[END_REF]. However a large part of our patients has had such a measure before entering in the study and comparisons may have suffered from a lack of power.

Other tools aiming to objectively measure the patient adherence have been proposed and have reported very high rates of poor adherence. The Medication Event Monitoring System (MEMS) records all openings of the jar containing the HCQ pills. Seventy-six percent of patients were shown to have taken less than 80% of the prescribed pills over a 2-year follow-up [START_REF] Marengo | Measuring therapeutic adherence in systemic lupus erythematosus with electronic monitoring[END_REF]. A similar rate was reported when analyzing pharmacy refill information of the US Medicaid data [START_REF] Feldman | Medication nonadherence is associated with increased subsequent acute care utilization among medicaid beneficiaries with systemic lupus erythematosus[END_REF]. Also based on pharmacy refill information, 51% of patients had an adherence to HCQ < 80% in a random sample of 63 SLE patients attending rheumatology clinics associated with University Medical Centers [START_REF] Koneru | Adherence to medications in systemic lupus erythematosus[END_REF].

The poor agreement between methods of assessment has already been highlighted by others, raising interesting questions. This could be due to differences in the recall periods.

The recall period is 30 days for MASRI. One MMAS-8 item refers to the day before and another one to the previous two weeks; there is no recall period for the remaining 6 items.

Blood HCQ concentration over 200 µg/L is reached after only few days of treatment [START_REF] Jallouli | Determinants of hydroxychloroquine blood concentration variations in systemic lupus erythematosus[END_REF].

Therefore the 200 µg/L threshold may be reached even in case of global poor adherence J o u r n a l P r e -p r o o f provided the patient had thoroughly taken the treatment for a few days before the visit, the so-called "white coat" compliance effect. In addition, the different tools may capture different elements of non-adherence which might explain the lack of overlap between tools. Self-questionnaires refer to perception whereas blood HCQ concentration refers to pharmacokinetics. Moreover, low blood HCQ concentration and questionnaires assess to different patterns of non-compliance. Blood HCQ concentration < 200 µg/L reflects severe poor adherence but blood HCQ concentration above this threshold does not define good adherence. Episodic omissions will not translate in low blood HCQ concentration.

Conversely MASRI and MMAS-8 may detect such behaviors. Finally, it cannot be excluded that some patients had low blood HCQ concentration for pharmacokinetic reasons despite a correct adherence to HCQ regimen (due to digestive malabsorption for example). This has been discussed for a cut-off of 500 g/L, but seems unlikely or at least very infrequent for a lower threshold like 200 g/L.

We found no specific characteristics in poorly adherent patients compared to adherent patients. This could be due to the small number of poor adherent patients and the subsequent lack of power of statistical analyses. In contrast, Costedoat-Chalumeau [START_REF] Costedoat-Chalumeau | A prospective international study on adherence to treatment in 305 patients with flaring sle: Assessment by drug levels and selfadministered questionnaires[END_REF] found that younger age at diagnosis, non-use of steroids, higher body mass index, and unemployment were independent predictors of non-adherence. These differences may be explained by the inclusion criteria. Costedoat-Chalumeau included only SLE flaring patients whereas in our study there were mostly non-flaring patients. In a cohort of SLE patients in USA, patients with sub-therapeutic blood HCQ concentration (< 500 g/L) were preferably females, less than 30 or over 60 years old, and had a lower blood concentration in vitamin D [START_REF] Durcan | Hydroxychloroquine blood levels in systemic lupus erythematosus: Clarifying dosing controversies and improving adherence[END_REF]. Using pharmacy refill information, Koneru et al [START_REF] Koneru | Adherence to medications in systemic lupus erythematosus[END_REF] showed that significant risk factors of insufficient adherence included being single, J o u r n a l P r e -p r o o f having low educational level, presence of other comorbidities, limited comprehension of physician explanations and instructions, and having to take the medication more than once daily. Using MMAS-4, the Lupus Erythematosus Long-Term Study (LuLa-study), a nationwide longitudinal study among German Caucasian patients with SLE, reported that the use of azathioprine, prednisone < 7.5 mg, being older and having an external health locus of control were predictors for high adherence. On the contrary, the general perception of medication being harmful, or addictive was detrimental [START_REF] Chehab | Medical adherence in patients with systemic lupus erythematosus in germany: Predictors and reasons for non-adherence -a cross-sectional analysis of the lulacohort[END_REF]. Using Medication possession ratio, Lee et al found that low disease activity was a predictor of poor adherence in Korean SLE patients initiating HCQ [START_REF] Lee | Compliance and persistence with hydroxychloroquine in south korean patients with systemic lupus erythematosus[END_REF]. We also found no clinical consequences of severe poor adherence to HCQ with similar SLEDAI score, SLIC damage index and quality of life, which contrasts with previous results. This could be explained by the fact that our patients were recruited during a routine visit, and/or by the low proportion of non-adherent patients and patients with active SLE whatever the adherence assessment method.

Should the blood HCQ threshold be questioned? A higher threshold such as 500 g/L has been advocated. Durcan et al reported that in a cohort of US SLE patients, 44% had subtherapeutic level (< 500 µg/L) [START_REF] Durcan | Hydroxychloroquine blood levels in systemic lupus erythematosus: Clarifying dosing controversies and improving adherence[END_REF]. In our study, using a 500 g/L threshold obviously increased the rate of non-adherent patients, but did not improve the agreement between assessment methods. Should an assessment tool be emphasized to measure HCQ adherence? Recognizing poor adherence is of paramount importance for its deleterious consequences. All assessment tools have their own limitations. The physician opinion is probably one of the worst tools to differentiate poor and good adherent patients highlighting the need for other tools.

Questionnaires offer an insight in patients' perceptions about medication intake J o u r n a l P r e -p r o o f behavior, but cannot reflect the real behavior [START_REF] Pasma | Nonadherence to disease modifying antirheumatic drugs in the first year after diagnosis: Comparing three adherence measures in early arthritis patients[END_REF]. MEMS would certainly be a good tool, since jar openings reflect the daily adherence over the whole period between visits and offer the advantage of assessing adherence over a continuum. Episodic omissions, as well as prolonged treatment discontinuations, may be detected and distinguished calling for specific interventions. Jar openings may not be synonymous of medication intake, but this method has proven to be superior to patient self-reports and pill counting for measuring adherence [START_REF] De Klerk | The compliance-questionnaire-rheumatology compared with electronic medication event monitoring: A validation study[END_REF]. Considering MEMS as a gold standard, Pasma et al found no or weak correlation with adherence questionnaire or methotrexate polyglutamate concentration in a population of patients with rheumatoid arthritis [START_REF] Pasma | Nonadherence to disease modifying antirheumatic drugs in the first year after diagnosis: Comparing three adherence measures in early arthritis patients[END_REF]. However, this cannot be routine practice for economic and practical considerations. Self-questionnaires are easy to use, but some of them are licensed; blood HCQ dosage is not expensive but is not performed in all centers.

Even if the combination of objective and subjective definitions might result in an overestimation of the non-adherence rate, over-diagnosing poor adherence would only lead to uselessly remind the importance of good adherence to patients already convinced. Since both MEMS-based and pharmacy refill-based studies have reported very high rates of non-adherence, most patients would take advantage of regular reminders. This position has also been advocated by Costedoat-Chalumeau et al [START_REF] Costedoat-Chalumeau | A prospective international study on adherence to treatment in 305 patients with flaring sle: Assessment by drug levels and selfadministered questionnaires[END_REF] in SLE flaring patients before questioning treatment efficacy and escalating the therapy.

Our study has strengths and limitations. The major strength was that patients entered the study during a routine visit. Indeed, previously informed patients could have adhered to the prescribed regimen only for a few days before the visit and been wrongly considered to be adherent, as a few days on HCQ are usually enough to increase the HCQ J o u r n a l P r e -p r o o f concentration over 200 g/L. Moreover, no patients refused to participate when informed on the aims of the study. Therefore, we are confident that our population was representative of the SLE patients routinely taken in charge by our regional network. The study also had some limitations. Blood HCQ concentration was measured on the day of the visit and measurement was untimed relative to the last dose of HCQ. However, regarding the long half-life of HCQ and the relative stability of blood HCQ concentration at steady state, a very limited impact can be expected on our results [START_REF] Al-Rawi | Steady-state pharmacokinetics of hydroxychloroquine in patients with cutaneous lupus erythematosus[END_REF]. Since there is no gold standard for the measurement of adherence to HCQ, we were unable to question the cutoff of 200 g/L for blood HCQ concentration. A further study combining MEMS as the gold standard, self-questionnaires, the calculation of the dispensed boxes by the pharmacist and blood HCQ would be of great interest to determine optimal thresholds.

Conclusion

Depending on the assessment method, the adherence rate to HCQ treatment varied between 3.2% and 32.5% with poor correlation and agreement between blood HCQ concentration, self-questionnaires on adherence and health care provider opinion. This could be explained by the different dimensions explored by each of these tools: from blood drug concentration to self-report of behaviour. Blood HCQ concentration < 200 µg/L reveals severe non-adherence. Combining blood HCQ concentration with MASRI and MMAS-8 may help to better identify non-adherence in SLE, thus guiding educational therapy.

Disclosure of interest

Acknowledgements

The authors acknowledge the other members of the ESSTIM investigators group: Thomas Quemeneur, Valenciennes; Pierre Savinel, Lens; Eric Auxenfants, Roubaix, Aurélia Lanteri, Lille. Sarahe Dehimat is acknowledged for her technical research assistance during all the protocol, Quentin Ducrocq for his MMAS-8 training and Géraldine Wojtasik for her implication in patient education. Licensing A Morisky Widget MMAS-8 License Agreement was made between the CHU de Lille -Hôpital HURIEZ (HURIEZ) 1 PLACE DE VERDUN Service de médecine interne 59037 Lille cedex ("Licensee") and MMAS Research LLC, 14725 NE 20th St. Bellevue Washington 98007.

EH, MMF, NLG, DL, HM, MHB, RA, TQ, AB, PB, ML, AL, ALB, VS, PYH, SD, GW and SMD have no conflict of interest. PC works as an independent statistician and has no conflict of interest to declare.

Grant

This work was supported by the "Conseil Régional Hauts de France", Direction de la Recherche, Enseignement Supérieur, Santé et Technologies de l'Information et de la Communication, 151 Boulevard Hoover, 59555 LILLE cedex.


Nguyen Hoang

Viet Hoang

Mai Suan Li

Junmei Wang

Philippe Derreumaux

Thi Ly Mai

Phuong H Nguyen
email: nguyen@ibpc.fr

Molecular dynamics simulation of cancer cell membrane perforated by shockwave induced bubble collapse

It has been widely accepted that cancer cells are softer than their normal counterparts. This motivates us to propose, as a proof-of-concept, a method for the efficient delivery of therapeutic agents into cancer cells while normal cells are less affected. The basic idea of this method is to use a water jet generated by collapse of the bubble under shockwave to perforate pores in the cell membrane. Given a combination of shockwave and bubble parameters, the cancer membrane is more susceptible to bending, stretching and perforating than the normal membrane because bending modulus of the cancer cell membrane is smaller than that of the normal cell membrane. Therefore, therapeutic agent delivery into cancer cells is easier than in normal cells. Adopting two well-studied models of the normal and cancer membranes, we perform shockwave induced bubble collapse molecular dynamics simulations to investigate the difference in response of two membranes over a range of shockwave impulse 15-30 mPa.s and bubble diameter 4-10 nm. The simulation shows that the presence of bubble is essential for generating water jet which is required for perforation, otherwise pores are not formed. Given a set of shockwave impulse and bubble parameters, the pore area in the cancer membrane is always larger than that in the normal membrane. But too strong shockwave and/or too large bubble results in too fast disruption of membranes, and pore areas are similar between two membrane types. The pore closure time in the cancer membrane is slower than that in the normal membrane. The implications of our results for applications in real cells are discussed in some details. Our simulation may be useful for encouraging future experimental work on novel approaches for cancer treatment.

I. INTRODUCTION

The plasma cell membrane regulates the entry and exit of substances via specific mechanisms. Small substances, such as ions, sugars and amino acids can cross the membrane via ion channels or protein pumps. Macromolecules must be carried into cells by membrane bound vesicles through the endocytosis mechanism [1]. These two mechanisms limit the entry of molecules which are not naturally needed by cells. As a consequence, although a large number of drugs is available for treating cell diseases, only a few drugs are used if they can cross the cell membrane. Therefore, a number of approaches has been proposed aimed at enhancing the transport of drugs through the cell membrane. Two well-known, vi-ral and chemical methods, have been developed but they are still limited by the low efficiency and toxicity. As a remedy, the physical approach has been recently proposed for highly efficient drug delivery with low toxicity. The basic idea of the physical approach is to create transient pores in the cell membrane, called poration, where drugs can penetrate into cell easier. Physical methods include photoporation, electroporation, magnetoporation and mechanoporation which use different kind of physical forces such as magnetic, thermal, mechanical and electrical forces to open the pore. For an excellent review of the drug delivery methods, readers are referred to recent publication [2]. Among these methods, the mechanoporation using the shockwave force is promising because it can rapidly deliver large macromolecules into cells in a local, noninvasive and cost effective manner [3]. For examples, Kodama et al. used shockwave to deliver large cytoplasmic molecules into cells. The authors suggested that the shear force generated by the shockwave temporarily af-fects the permeability of the membrane, and the impulse of the shockwave plays an important role in governing the permeability. This method can deliver large macromolecules of up to 2 × 10 6 molecular weight into cells [4]. Recently, Lopez-Marin et al. have studied the shockwave induced damage and poration in cell line HEK293 and tumor-derived cell line MCF-7, and results from scanning electron microscopy revealed transient hole-like structures after shockwave exposure. The authors also showed differences in membrane permeability of two cell lines [5]. Qi et al. showed that shockwave can trigger the release of ATP from osteosarcoma U2OS cells by increasing the membrane permeability [6]. The ability of shockwave to enhance the delivery of very large molecules suggests that shockwave should also be a promising method in gene therapy and protein delivery [7,8]. We should mention that the shockwave has been used in other medical applications. For example, high energy shockwaves have been used for more than 30 years to disintegrate urolithiasis. Extracorporeal shockwave therapy has been clinically used for many musculoskeletal conditions. It has been suggested that shockwave accelerates tissues regeneration, reduces calcification and inhibits pain receptors [9].

To further develop this sensitive method aimed at enhancing drug delivery into cancer cells, one needs a better understanding of the mechanism of shockwave interacting with the cell membrane, especially at the molecular level. Unfortunately, it is very difficult to observe experimentally the direct interaction between shockwave and membranes due to the short time scale on the order of picoseconds, and small length scale of several nanometers. Therefore, several studies have employed molecular dynamics (MD) simulations to investigate the molecular mechanism of the shockwave induced membrane damages, with and without bubbles [10][11][12][13][14][15][16][17][18][19][20][21][22][23]. Koshiyama et al. performed for the first time an all-atom nonequilibrium MD simulation under a shockwave but without bubbles and observed the penetration of waters into the hydrophobic region of the membrane, which was caused by a decrease in membrane thickness [10]. Choubey et al. performed large-scale all-atom MD simulations of lipid membranes with shockwave and nanobubbles and showed that the bubble collapse generates shear flow of water on membrane leaflets and pressure gradients across them, creating transient nanochannels through which water molecules translocate across the membrane [11]. Steinhauser and colleagues developed an advanced large and multi-scales coarse-grained simulation method using dissipative particle dynamics, and applied it to study effects of shockwave on DPPC membranes. The authors obtained a threshold shock front velocity, below which the membrane recovers from shockwave induced damage, and above that the membrane could not be recovered [12][13][14][15]. Santo et al. employed coarse-grained MD simulation to study the impact of shockwave and bubbles on the damage and recovery of lipid membranes. Interestingly, they showed that not every lipid molecule remained in the bi-layer after recovery, some lipids moved out into water and created micelles [16,17]. Adhikari et al. showed that in the absence of bubbles, high intensity shockwaves do not induce pores on membranes, but weaker impulses can lead to membrane poration in the presence of bubbles [18]. Similarly, Lu et al. performed MD simulation of the shockwave induced delivery of paclitaxel drug through a lipid membrane and showed that the paclitaxel molecule can penetrate the membrane only under the joint effect of the shockwave and nanobubble [20]. Sun et al. simulated the shockwave induced collapse of lipid-shelled nanobubbles and interestingly, they showed that compared with the cases of vacuum nanobubbles, the lipid nanobubbles could weaken the effects of shock waves [19]. Nan et al. performed coarse-grained MD simulation of the shockwave induced bubble collapse and observed not only membrane perforation but also the occurrence of nanoscale cavitation during the perforation process [21]. Hu et al. analyzed effects produced by the collapse of multiple nanobubbles in the vicinity of biomembranes in the presence of an electric field by MD simulation [22], and showed that multiple nanobubbles make it possible to create larger pores on the membrane [22]. Very recently, Wei et al. studied the impact of shockwave induced bubble collapse on the damage of cell membranes with different lipid peroxidation levels, and showed that the pore sizes increase with the peroxidation level [23]. One of important findings obtained from all these simulation studies is that the shockwave alone does not have much impact on membranes, but the presence of gas bubbles, which pre-exist or are nucleated during the shockwave propagation, plays an essential role in the membrane poration.

Our main interest is the application of the shockwave in combination with bubble in the drug delivery for cancer treatment. To date cancer treatments, including surgery, radiotherapy and chemotherapy are widely used. Extensive resection followed by adjuvant chemoradiotherapy seems to be the only treatment that modifies the survival of cancer patients. However, with the chemotherapy, the therapeutic agent affects not only cancer cells but also normal cells, causing severe side effects. Therefore, novel therapeutic methods are needed to minimise side effects on the normal cells. Our core aim is to develop such a method, based on the difference in the mechanical properties of normal and cancer membranes. Recently, a number of studies has shown that cancer cells from a large number of different organs are softer than their normal counterparts [24]. Therefore, in recent years, mechanical properties of cancer cells have been suggested as biomarkers for early cancer diagnosis, targeted for the prediction, treatment and even prevention of cancer [25]. Our research hypothesis is that if cancer cells are softer than normal cells, then we can choose shockwaves with appropriate impulses and bubbles with suitable sizes so that under the shockwave induced bubble collapse, the cancer cell membrane is maximally perforated while the normal cell membrane is minimally affected. Therefore, therapeutic agents are maximally delivered into cancer cells but minimally into normal cells.

In a previous work, we performed equilibrium MD simulations to study the elastic properties of normal cell membranes and cancer cell membranes [26]. In our simulations models, normal cell membranes have a highly asymmetric lipid composition [27], where the extracellular leaflet is mainly composed of phosphatitylcholine (PC) and sphingolipids and the intracellular leaflet is mostly composed of phosphatidylethanolamine (PE) and phosphatidylserine (PS) lipids. When normal membranes are transformed to cancer membranes the concentration of the negatively charged PS lipids is increased in the extracellular leaflet [START_REF] Connor | Proc. Natl. Acad. Sc. U.S.A[END_REF][START_REF] Utsugi | [END_REF], and the cholesterol (CHL) concentration is reduced [30]. We calculated the elastic moduli including bending, tilt and twist constants of the normal and cancer membranes. The results showed that at low cholesterol concentrations, all elastic moduli become smaller, implying that the cancer membrane is indeed softer than the normal counterpart.

In this work, we carry out nonequilibrium MD simulations of shockwave induced bubble collapse of a normal membrane and a cancer membrane whose structures are known from our previous work [26]. We vary the shockwave impulses and bubble diameters over a wide range of values to observe the difference in response of the normal and cancer membranes. The main finding is that given a combination of a shockwave and a bubble, the cancer membrane is more bent, stretched and perforated than the normal counterpart. In a previous work [26], we carried out all-atom MD simulations of five normal membrane models and five cancer membrane models. The cancer membranes were taken into account the overexpression of PS lipids in the outer leaflet and the reduction of cholesterol concentration. Results showed that at the same cholesterol concentration the bending moduli of the normal and cancer membranes are very similar. This indicates that the overexpression of the PS lipids does not affect significantly the elasticity of the normal or cancer membranes. At low cholesterol concentrations, all elastic moduli become smaller, implying that the reduction in cholesterol in cancer membranes could contribute, at least partly to the softening of cancer cells. Therefore, in this work, we consider only two membrane models among the ten membrane models used in previous work [26]: a normal membrane model and a cancer membrane model which differ in the PS lipids concentration in the outer leaflets, and in the cholesterol concentration. Each model contains four lipid types: 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC), sphingomyelin lipids (SM), 1,2-dioleoyl-sn-glycero-3-phosphoethanolamine (DOPE) and 1,2-dioleoyl-sn-glycero-3-phospho-L-serine (DOPS) lipids. For each model, we take the last equilibrium membrane structure at 2 µs obtained from previous simulations [26], and translate it along the x and y-axis to obtain a doubled size membrane with lengths (x, y) = (24,24) nm. This guarantees that the new membrane structure is in the equilibrium state, and large enough to accommodate large bubbles having diameters of ∼ 10 nm. The numbers of lipids and cholesterol of two membrane models are listed in Tab.I.

B. Shockwave induced bubble collapse simulation method

In this work, we employ the shockwave method originally developed by Koshiyama et al. [10]. The sketch of the simulation setup is shown in Fig. 1. The system is a rectangular prism box composed of two lipid membranes described above, a bubble and the whole system is solvated in water. We are interested in the lower membrane, while the upper membrane is only technically used to prevent the shockwave coming back to another side of the lower membrane due to the periodic boundary condition used in the simulation [31,32]. The distance between the membrane and the bubble is d mb , and the diameter of the bubble is d b . To generate a shockwave with an impulse I, an excess momentum M = I × A is applied to a slab of waters with a volume of A × d s . Here, A and d s are the area in the (x,y) plane and the thickness of the water slab, respectively. Initially, the distance between the slab of water and the bubble is d bs .

The all-atom CHARMM36 force field [33] and the TIP3P water model are used to model the lipids and solvent, respectively. The initial dimensions of the primary unit cell are (x, y, z) = (24, 24, 54) nm, consisting of 3.112.724 atoms in total. Starting from the initial structure of the normal membrane or cancer membrane, an equilibrium MD simulation is carried out for 100 ns in the NPT ensemble at the pressure P 0 = 1 bar and temperature T = 300 K. Then, we take the last ten equilibrium structures separated by 100 ps, and for each structure we remove water molecules in the spherical region hav- ing diameter d b near the lower membrane. This empty water space mimics a bubble. In this work, we consider four bubble sizes with d b = 4, 6, 8, 10 nm . These membrane/bubble systems are used as initial structures for the shockwave simulations. The shockwave I defined as the time integral of pressure over the shock pulse duration is varied from I = 15 mPa.s to I = 25 mPa.s. The shockwave induced bubble collapse MD simulation is performed by using the GROMACS simulation package [34] in the NVE ensemble. The equations of motion are integrated using the leapfrog algorithm with a small time step of 0.5 fs. The small time step is used to ensure the stability of the simulations. The electrostatic interactions are calculated using the particle mesh Ewald method and a cutoff of 1.2 nm [35], and the forces switched to zero from 1.0 to 1.2 nm are used for the van der Waals interactions. The nonbonded pair lists are updated every 5 fs. Each simulation is run for 50 ps and the data is saved for every 500 MD steps (every 0.25 ps) for subsequent analyses.

C. Data analysis

In our simulations the shockwave moves in the zdirection, thus to monitor the propagation, we calculate profiles of the normal pressure component, P zz (z), and temperature, T (z). To this end, we discretise the simulation box in thin slabs of width 0.5 nm which are oriented parallel to the xy-plane. By computing the average pressure and temperature in each slab, the propagation of can be monitored in details. Here, pressure is calculated by using the method developed by Ollila et al. [36], and temperature is directly calculated from velocity of atoms.

To measure the membrane pore area ∆S, we place a 2D grid with the size of each grid square of (0.1 × 0.1) nm on the (x, y) membrane surface, and count the squares on the grid having zero lipid atoms [18,23].

The order parameters of the lipid acyl chain tails are calculated as S i = 3 cos 2 θ i -1 /2, where θ i is the angle between the i-th C-H bond vector and the bilayer normal [37], and the angular brackets represent ensemble average over all lipids.

III. RESULTS

We have performed MD simulations to investigate the difference in response of the normal and cancer cell membranes to the shockwave induced bubble collapse. There are five parameters in the simulation setup including the shockwave impulse I, diameter of bubble d b , distance from bubble to membrane d mb , distance from water slap to bubble d bs , and thickness of water slap d s [Fig. 1]. Among these, two physical parameters that can be controlled experimentally are the shockwave impulse

A. Shockwave induced bubble collapse and membrane response

To describe the shockwave propagation, the collapse of the bubble and the response of the membrane in details, let us consider a representative simulation of the normal membrane using a shockwave impulse I = 21 mPa.s and a bubble with diameter d b =10 nm. We calculate various quantities including the vector field of velocity of atoms, the profiles of the pressure, kinetic energy, temperature and mass density. To provide an intuitive picture, snapshots are also visualised. These results are shown in Fig. 2 and Fig. 3.

By construction, all atoms within the water slab above the bubble [Fig. 1] are initially assigned to the same high velocity V , which is calculated as V = (I × A)/(m × N ) [10]. With the impulse I = 21 mPa.s, the area of the water slab A = 576 nm 2 , the weight of a water molecule m = 18 g/mol and the number of water molecules in that slab N = 58044 we obtain V ∼ 9 km/s with direction parallel to the z-axis and towards to the membrane. This results in an excess high temperature of ∼ 15000 K and a pressure of ∼ 44 GPa in the water slab at t = 0 ps [Figs. 2A,B]. These water atoms then move extremely fast toward the membrane as seen from the velocity vector field in Fig. 3A. This causes an increase in water density profile in the direction of the membrane with the maximum values of 1532 kg/m 3 at t ∼ 1.25 ps, z ∼ 21 nm and 1519 kg/m 3 at t ∼ 3 ps, z ∼ 15 nm, thus leaving an empty space behind them at the initial position of water slab (z ∼ 30 nm) [Fig. 2C and Fig. 3B]. This is because the initial equilibrated water molecules at 300 K above the water slab are not quick enough to occupy the space left by the water in the slab. As a consequence, the temperature and pressure in the water slab are reduced significantly from the initial values, and increased quickly in the direction of the membrane, reaching the maximum values of ∼ 3400 K and 20 GPa, respectively at t ∼ 1.25 ps and z ∼ 21 nm [Fig. 2A,B]. Along the way, a part of water in the slab enters the empty space of the bubble without any obstacle, mimicking the bubble collapse. At t = 1.25 ps, the bubble is partially collapsed [Fig. 3B] but the membrane has not been affected yet, thus the average kinetic energy of the membrane is still maintained at an equilibrium value of ∼ 4-5 kJ/mol, uniformly distributed over all lipids [Fig. 3C], and the temperature is still around the initial value of ∼ 300 K [Fig. 2A]. In contrast, the propagation of the part of waters in the slab, which do not enter the bubble, is slower due to the resistance of the water below the slap. The difference in velocity between these two parts of water, i.e, one enters the bubble and one outside the bubble, produces highspeed water jet directed towards the membrane as shown by the velocity vector field in Fig. 3A. After entering the bubble, the water jet continues to propagate, and occupies completely the bubble space at t ∼ 3 ps, mimicking the full collapse of the bubble [Fig. 3E]. Then, it hits the lipid membrane and velocity drops down to ∼ 6 km/s with the velocity vector field shown in Fig. 3D. The membrane area located just below the bubble receives directly the kinetic energy from the water jet, thus its tempera- map in Figs. 3G,I, respectively at 50 ps. The membrane is now highly stretched and pores are formed [Figs. 3H,I]. We note that the membrane is basically not translated under shockwave as seen from the position z of the membrane density profile shown in Fig. 2D. This is due to the use of two membranes technique in our simulation as described above.

Next, we calculate the average temperature and normal pressure across membrane surface as a function of time and results are shown in Fig. 4. The temperature remains at the initial value of ∼ 300 K until the arrives. Upon the arrival, the temperature increases sharply to a maximum values of ∼ 700 K at ∼ 5 ps, then it slowly decreases and reaches ∼ 560 K at 50 ps. Similarly, there is a sharp increase in pressure from the initial value ∼ 10 -4 GPa (1 bar) to the maximum value of 4.2 GPa at ∼ 5 ps, followed by a rapid decrease of pressure which reaches a minimum value of ∼ 0.35 GPa at 13 ps and then the pressure rises again toward the final value of ∼ 0.73 GPa at 50 ps. We should mention that the shockwave induced bubble collapse simulation is carried out in the NVE ensemble, and the very high initial kinetic energy of the shockwave is redistributed over all atoms, thus the temperature and pressure of the membrane do not relax to the initial values of 300 K and 1 bar, respectively. Of note, we analysed trajectories simulated with different shockwave impulses and bubble sizes, and found that the mechanism of bubble collapse and membrane response is similar for all simulations.

It is of interest to compare the targeted impulse I, which is used to generate the initial shockwave, with actual pressure impulse exerted by the shockwave on the membrane. From the pressure profile P (t) shown in Fig. 4, the pressure impulse on the membrane is calculated as I mem = t 0 P (τ )dτ , where t is the time duration of the positive phase of the shockwave [10]. The results are shown in Tab.II for the simulations using with a 10 nm bubble and different targeted impulses I. As seen, II. A comparison between the targeted impulse I and the actual pressure impulse Imem exerted on the membrane by the shockwave induced bubble collapse. v slab is the initial velocity assigned to water atoms in the water slab for a given targeted impulse. In all cases, the time duration of the positive phase of the shockwave t = 10 ps is used in the calculation. Shown are results obtained from the simulations of the normal membrane using a 10 nm bubble.

I [mPa.s] v slab [km/s]
the actual pressure impulses are very close to the targeted counterparts, indicating that the initial impulses are not dissipated much before hitting the membrane. Therefore, the membrane is indeed impacted by the targeted impulses.

B. Response of normal and cancer membranes upon bubble collapse

understood the mechanism of shockwave induced bubble collapse, we now wish to investigate the response of the normal and cancer membranes in details, focusing on the membrane bending, lipid ordering and membrane pore formation. In the previous simulation we calculated the elastic moduli of the normal and cancer membranes in equilibrium state and obtained the bending modulus K c = 10.5 × 10 -20 J for the normal membrane and K c = 8.1 × 10 -20 J for the cancer membrane [26]. This means that the cancer membrane is softer than the normal membrane. This implies that under the same water jet, the cancer membrane should be deformed more than the normal membrane. To confirm this, we carry out simulations using a 10 nm diameter bubble and a relatively weak shockwave of 15 mPa.s, such that only the membrane bending, but not pore formation, is observed. To obtain a qualitatively impression, Fig. 5 shows snapshots of the two membranes at t = 3 ps, and indeed, the cancer membrane is bent more than the normal membrane, with a bending of 2.0 nm and 1.2 nm, respectively.

To exam the response of membranes in more detail, we follow the time-evolution of the order parameters of the lipid acyl chain tails during the simulation of the normal and cancer membranes. For the DOPC, DOPE and DOPS lipids, each tail consists of 17 order parameters (i = 1 • • • 17), and for PSM lipids, each tail consists of 15 order parameters (i = 1 • • • 15). For simplicity, we sum all order parameters of each lipid type: S lipid = i S i . As an example, the time-evolution of the four S lipid in the normal and cancer membranes is shown in Fig. 6 for the simulations using a 10 nm diameter bubble and a shockwave of 21 mPa.s. As seen, both normal and cancer membranes are initially at equilibrium, and their order parameters are very similar. At t ∼ 3 ps, the bubble is completely collapsed and the water jet hits the membrane. Thus the membrane is highly disturbed and lipids become disordered. This is indicated by a largest reduction of ∼ 60% in the order parameters of all lipid types at 3 -4 ps for both membranes. Then, the intensity of the water jet is reduced because its kinetic energy is transferred to other degrees of freedom of the membrane and water. After hitting the lower membrane, the water molecules bounce back, and the membrane also tends to be pulled back to its original state as indicated by the increase in all order parameters of all lipid types from 3 ps to 18 ps. However, due to the presence of the upper membrane, the water molecules bounce back towards the lower membrane again causing the lipids to become disordered as indicated by the decrease in order parameters at t ∼ 18 -20 ps. Then, the kinetic energy of the water molecules is redistributed, and the oscillation of the water flow between the upper and lower membranes disappears, the membrane is stabilised, and the order parameters approach to horizontal values (see at 50 ps). Importantly, we observe that the cancer membrane is more disordered than the normal membrane. The order parameters of the cancer membrane are ∼ 14% lower than that of the normal membrane [Fig. 6].

To see what changes in the structure of the membrane lead to a decrease in the order parameters S i at ∼ 3 ps seen in Fig. 6, we calculate the distribution of the θ i angle between the i-th C-H bond vector and the bilayer normal for different lipid types. We found that at equilibrium these angles are mainly distributed around 0 degree, thus S i are large. Upon hit by the water jet around 3-4 ps, the distributions of θ i are broader due to the compressed leading to smaller order parameters. Therefore, the reduction in the order parameters is associated with the reduction in the membrane thickness.

The main aim of the shockwave induced bubble collapse method is to create pores on the cell membrane where drugs can be efficiently delivered into the cell. Therefore, the pore area is an important quantity to validate the efficiency of the method. We calculate the pore area ∆S (see Section II.C) in the normal or cancer membrane at equilibrium and obtain an area of ∆S ∼ 3 nm 2 . Of course, there is no well-defined pore in an equilibrium membrane, and this value is simply the sum of all small empty spaces in the (x, y) plane of a membrane at equilibrium.

First, we investigate the dependence of the pore area on the shockwave impulse, given a bubble size. Fig. 7 shows, as an example, the time evolution of the pore area induced by the collapse of a 10 nm bubble under different shockwave impulses I = 15-25 mPa.s for the normal and cancer membranes. Overall, the same feature of the time evolution of the pore area is observed for all impulses and for two membranes. At t = 0 ps when the water jet has hit the membrane yet, both membranes are at equilibrium and ∆S ∼ 3 nm 2 . Then, the shockwave induces the collapse of the bubble and the generated water jet hits the membranes. With I = 15 mPa.s, the water jet is not strong enough to perforate pores in the membrane. The membranes are slightly bent as seen from Fig. 5, and Clearly, for both membranes the pore area increases with the increase of the shockwave impulse. For instant, at a low impulse intensity of 19 mPa.s the maximum pore area is around 4.5 nm 2 but with an intensity I = 25 mPa.s the pore area is three times larger, around 13 nm 2 . We note that in all cases, largest pores are formed at t ∼ 10 -15 ps, while the maximum reduction in the lipid order parameter occurs much faster, around 3-4 ps [Fig. 6]. To explain this, we note that the membrane outer leaflet (facing the water jet) is affected immediately upon hit by the water jet, while the inner leaflet has not been affected yet [Fig. 5].

Thus, at this moment the lipid order parameters, which are calculated for lipids pertaining in both leaflets, are affected. In contrast, a pore in a bilayer are fully formed only if pores are formed in both leaflets, and this process is slow.

At weak shockwave impulses I = 15 -17 mPa.s, there is no difference in ∆S between the normal and cancer membranes because in these cases, pores are not really formed. At I = 19 mPa.s, the pore area in the cancer membrane is slightly larger than that in the normal membrane. The difference becomes more obvious with I = 21 -23 mPa.s, with the maximum pore areas in the cancer membrane are 19-23% larger than that in the normal membrane. Interestingly, with I = 25 mPa.s, the pore areas in both membranes become very similar. This is because if the shockwave is too strong then the speed of the water jet is very high, it hits and ruptures the membrane instantaneously without membrane bending and stretching, thus the results do not depend much on the membrane elasticity. Finally, we note that the difference in the average value of ∆S between the two membranes is larger than the standard deviation values, indicating that the cancer membrane is indeed more prone to perforation than normal membrane.

Next, we wish to investigate the dependence of the pore area on the bubble size given a shockwave impulse. Fig. 8 shows the time evolution of the pore area induced by the collapse of bubbles having diameters d b = 4 -10 nm under a shockwave impulse I = 23 mPa.s. This impulse is chosen because the difference in the pore area between the two membranes is most obvious [Fig. 7]. We note that in all cases, the shockwave impulse is the same thus the pressure exerted by the water jet on the membranes should be the same. However, it is clear that the larger the diameter of the bubble, the larger the diameter of the water jet and therefore the larger the area of the membrane impacted by the water jet. As a consequence, small bubbles with d b = 4 and 6 nm induce small pores with ∆S ∼ 3.5 nm 2 , while larger bubbles, d b = 8 and 10 nm, create larger pore areas. Importantly, given a bubble size the pore area in the cancer membrane is always larger than that in the normal counterpart. For example, with an 8 nm bubble, the largest pore areas are 4.5 and 6 nm 2 for the normal and cancer membranes, respectively. After 50 ps, these areas decrease to ∼ 3 and 4 nm 2 , indicating that the pore in the normal membrane is closed but that is still opened in the cancer membrane. Similar results are observed for the 10 nm bubble case [Fig. 8].

In summary, the above results show that (i) given the same bubble, a stronger shockwave impulse induces a larger pore area in the membrane, (ii) for a given shockwave impulse, a larger bubble induces a larger pore in the membrane, (iii) for a given set of the bubble and the shockwave parameters, the pore area in the cancer membrane is always larger than the pore area in the normal membrane, and (iv) if the shockwave impulse is too strong then the pore areas in the normal and cancer membranes become similar.

C. Pore closure

As shown above, the pore in a membrane induced by weak shockwave impulse or by collapse of a small bubble is quickly closed. For example, with a weak impulse of 15 mPa.s or a small 4 nm bubble, the pore area reaches the maximum value ∼ 3.5 nm 2 then quickly decreases to the initial equilibrium value of ∼ 3 nm 2 [Figs. 5,8]. However, with stronger impulses and/or bigger bubbles the pore areas are large, and the pore closure is very slow as indicated by slow relaxation of pore area to the initial equilibrium values after 50 ps.

To study the pore closure, we carry out two 100 ns MD simulations under equilibrium condition with temperature of 300 K and pressure of 1 bar for the normal and cancer membranes, starting from structures obtained at 50 ps of the shockwave simulations using shockwave impulse of 25 mPa.s and a 10 nm bubble. The initial pore areas in the normal and cancer membranes are ∼ 6.5 nm [Fig. 7]. As seen from Fig. 9, the time evolution of these pore areas under equilibrium condition is quite stable. The pore area in the normal membrane is decreased by ∼ 10%, from 6.5 to 5.8 nm, and that in the cancer membrane is decreased by only 3%, from 6.5 to 6.3 nm after 100 ns [Fig. 9]. This suggests that pore closure is rather slow, and importantly, the pore in the normal membrane tends to close more rapidly than the pore in the cancer membrane. To obtain an impression on the closing speed, we simply fit the data between 20 -100 ns shown in Fig. 9 to a linear function and obtain the closing time of ∼ 1.7 µs and 5.0 µs for the normal and cancer membranes, respectively. To explain the difference in the timescale of pore closure, we recall that the equilibrium bending modulus of the cancer membrane is smaller (K c =8.1×10 -20 J) than that of the normal (K c =10.5×10 -20 J) membrane [26]. Therefore, if two membranes are stretched by the same pressure exerted by the water jet, the normal membrane will return to its original state faster than the cancer membrane, assuming that the Hooke's law is still valid. Furthermore, we calculate the order parameters, S i , of lipids from 100 ns equilibrium trajectories after the shockwave excitation and compare results with those obtained from the equilibrium simulation without the shockwave [26]. As seen from Fig. 10, the order parameter of lipids in the normal membrane before and after shockwave excitation are quite similar. In contrast, there are differences in the order parameters, especially of the DOPC and DOPE lipids in the cancer membrane before and after shockwave excitation. This confirms that the cancer membrane relaxes to equilibrium more slowly than the normal membrane. Finally, to make contact to other works, we consider the impact of shockwave alone, i.e, without the bubble (d b = 0 nm). To this end, we perform simulations using a relatively strong shockwave with impulse of 23 mPa.s. Figs. 11(A),(B) show snapshots of the normal and cancer membranes, respectively at t = 3 ps when the shockwave just hits the membrane. As seen, the membranes are hardly affected, i.e, no bending or poration are observed. We increase the shockwave impulse up to 30 mPa.s but similar results are obtained. This indicates that the shockwave alone does not induce significantly structural changes in the membranes, in agreement with previous studies [12,18,20,38,39]. Furthermore, the response of the normal and cancer membranes to the shockwave is similar.

IV. DISCUSSION

We perform MD simulations of shockwave induced bubble collapse to investigate the difference in response of the normal and cancer cell membranes. These two membranes are different in the elastic moduli caused by the difference in the concentration of the phosphatidylserine lipids in the outer leaflet and the cholesterol in the bilayer [26].

We first discuss the impact of shockwave alone on two membranes. Under the shockwave excitation, a water flow is generated and uniformly directed towards the membrane, exerting an uniform pressure across all locations on the membrane, therefore the membrane is simply compressed but not perforated. Indeed, Berkowitz and colleagues performed shockwave induced bubble collapse simulations of a DPPC lipid membrane using the MARTINI coarse-grained lipid model and showed that a high intensity shockwave with an impulse of 18 mPa.s does not induce pores in the membrane [18]. Similarly, the simulation of Wei et al. for a single-component DOPC lipid membrane using the MARTINI force field showed slight compression of the membrane and no pores formation [23]. Lu et al. [20] and Wang et al. [39] studied DPPC lipid membranes using the all-atom GROMOS force field [START_REF] Van Gunsteren | Biomolecular Simulation: The GROMOS96 Manual and User Guide[END_REF], and showed that the membranes are not perforated under the shockwave alone. Our simulations of multi-component lipid membranes using the all-atom CHARMM36 force field show that shockwave alone does not induce pore formation, confirming previous studies using different lipid models, compositions, and force fields [18,20,23,39]. However, all these results contradict result of Espinora et al. who showed that a shockwave alone with an impulse of ∼ 0.45 mPa.s could damage the membrane to an unrecoverable state [38]. This could be due to the very small size membrane studied by Espinora et al. [38].

The presence of a bubble is essential for pore formation because it focuses the uniform water flow generated by the shockwave into a water jet which hits the membrane at a focal area. The peak pressure values at this area are very high, for example, up to 4.2 GPa in the simulation with I = 21 mPa.s and 10 nm bubble, which can perforate pores at the focal area. Using the piston method to generate the shockwave [START_REF] Wen | [END_REF], Adhikari et al. showed that the peak pressure at the membrane is ∼ 0.55 GPa in a simulation using a piston velocity of 1 km/s and a 60 nm bubble [18]. Ganzenmueler et al. performed large scale coarse-grained simulations with a piston velocity of 4.73 km/s but without a bubble and showed that the peak pressure is ∼ 3 GPa [12]. In all these studies, although the peak pressure values are very high but the durations of the shockwaves are very short, therefore the values of the pressure impulses exerted on the membrane are very small [Tab.II]. We should mention that these values are much lower than an experimentally measured value impulse of 54 Ps.s which is sufficient to deliver calcein molecules into the cell but does not cause the cell death [42].

We then study the response of the membrane as a function of the shockwave impulse and bubble size in detail. Overall, our results show that the stronger the shockwave impulse or the larger the size of the bubble, the larger the pore area. This result is similar to the earlier simulations of Wei et al. [23] and Sun et al. [19] although the authors used a different, i.e, momentum mirror method to generate the shockwave, and moreover the system sizes, lipid compositions as well as the force fields are different from our simulations.

The cancer membrane studied in this work has a lower bending modulus than the normal counterpart. As a consequence, for a given shockwave impulse and bubble combination, the cancer membrane is more bent and perforated than the normal membrane. The pore area in the cancer membrane is always larger than the pore area in the normal membrane. However, we do not observe extreme cases where the cancer membrane was perforated while the normal membrane was not. This is because the bending modulus of the cancer membrane is only ∼ 22% smaller than that of the normal membrane, thus the response of two membranes to the water jet is not much different. The largest differences of 20-23% in the pore areas between the normal and cancer membranes are observed for the cases with the shockwave impulses of 21-23 mPa.s and a bubble of 10 nm. We observe that increasing the shockwave impulse and/or bubble size to large values does not increase the difference in the pore area. This is because the pressure exerted by water jet on the membrane is too strong, lipids are immediately expelled from the membrane, and pores are formed without much bending of the membrane. This means that in this case the difference in the elasticity between the normal and cancer membranes does not really determine the difference in the pore areas. We show that the pore formed in the cancer membrane is closed more slowly than the pore in the normal membrane, implying that once deformed, the cancer membrane takes longer to relax to equilibrium than normal cell membrane. Again, this is because the cancer membrane has a lower bending modulus than the normal counterpart.

As mentioned, there have been a number of simulation studies on the interaction between the shockwave and membranes [10,11,[16][17][18][19][20][21][22][23]39]. However, to our best knowledge there is currently only one simulation of Wei et al. who studied the response of membranes with different elasticities on the shockwave induced bubble collapse. In that work, the membrane is a single-component DOPC lipid bilayer modelled by the coarse-grained MARTINI force field. The elasticity of the DOPC membrane is varied by varying the population of the peroxidized DOPC lipids. The bending modulus of the 100% peroxide lipid membrane is K c = 4.3 ×10 -20 J which is about half of that of the membrane without peroxide, K c = 8.4 ×10 -20 J. As a consequence, the authors showed that a shockwave impulse of 127.1 mPa.s is required to perforate the 100% peroxide membrane while a stronger impulse of 161.3 mPa.s is necessary for the pore formation in the 0% peroxide membrane. This result together with our result confirm that the shockwave induced bubble collapse damage threshold for low elastic membranes is lower than that required for high elastic counterparts.

It is instructive to compare the timescale of pore formation induced by shockwave induced bubble collapse with that obtained by other excitation methods. In a previous work, we simulated the effect of stable bubble cavitation on a lipid bilayer and showed that the bubble fuses with the membrane and subsequent cavitation pulls lipid molecules out of the membrane, creating pores after ∼ 30 ns [32]. In another work, we simulated the direct interaction between focused ultrasound and a lipid membrane and showed that the spatial pressure gradient between the focused and non-focused regions causes the pore formation after ∼ 225 ns [43]. If the pore opening process is slow then unwanted molecules may be able to enter the cell together with drugs during the pore opening. In this context, the shockwave method is preferable because it allows rapid drug delivery into cells.

The simulation does not show large differences in the pore areas in the normal and cancer membranes because the bending moduli of two membrane models are not much different. Is this true for real biological cells? In reality, cancer cells are much softer than normal cells [24]. For example, Lekka et al. used the scanning force microscopy to exam the elasticity of normal and cancer human bladder cells, and the results show that the Young's modulus, defined as a measure of cellular deformability, of cancer cells is one order of magnitude larger as compared to healthy cells [44]. Using high-throughput optical tweezers technique, three cell lines were compared: a nontumorigenic breast epithelial MCF10 cells, a non-motile, non-metastatic breast epithelial cancer MCF7 cells and MCF7 cells transformed with phorbol ester. The results showed a significant increase in the deformability in the transformed MCF7 cells as compared to both nonmetastatic MCF10 and non-transformed MCF7 cells [45]. Further studies provide a large database of cases showing significantly larger deformability of single cancer cells [46][47][48]. The very large difference in the elasticity between normal cells and cancer cells suggests that one can always be able to select a suitable range of parameters of the shockwave impulse and bubble such that cancer cells are maximally perforated while normal cells are minimally affected or even unaffected. Even if the normal cell is perforated, the pores will be closed more quickly, therefore if drugs are injected after the pores in the normal cells have closed, then drugs can only enter the cancer cells where pores are still open, i.e, the normal cells are safe.

We acknowledge that although this is a proof-ofconcept work, we believe that our proposed method can be realised experimentally. Indeed, shockwave has been applied to cancer therapy. Some experimental works have shown that the shockwave can suppress tumor growth and selectively kill malignant cells [15,49] . The technique that uses a laser-induced shockwave to drive a liquid microjet at a very speed has also been developed [50]. The gas bubbles have been widely used in ultrasound induced bubble cavitation aimed at enhancing drug delivery into cells. In this approach, ultrasound is used to induce the stable cavitation of microbubbles, generating microstreaming which exerts shear stresses on the cell membrane, resulting in pore formation or even disruption of the cell [51]. In this context, the combination of shockwave and bubble into an experimental method to rapidly and selectively perforate cancer cell membranes for drug delivery is quite doable.

V. CONCLUSION

We have performed MD simulations of shockwave induced bubble collapse on normal and cancer membrane models. We show that the combination of shockwave and bubble is essential for the pore formation in the membrane. That is, the perforation is due to the water nanojet generated by a collapse of bubble under shockwave. We show that given a combination of a shockwave and a bubble, the cancer membrane is more bent, and the pore area in the cancer membrane is larger than that in the normal membrane, because the cancer membrane is softer than the normal membrane. But the pore areas in the two membrane types become similar if the shockwave is too strong. Our simulation results could provide a proof-of-concept for the development of a new method that uses shockwave and bubble to rapidly deliver drugs into cancer cells while leaving normal cells less affected. For further development, it is important to obtain the pore area as a function of the bending modulus for various cancer cell membranes which have different elastic properties. This could provide a hint for optimising the shockwave impulse and bubble size in order to optimise drug delivery outcome. Finally, we note that the cell membranes considered in this work contain only lipid bilayers. It is necessary to include other components such as cytoskeleton which is responsible for providing struc-tural integrity and mechanical stability. To model such complex systems, one could resort to advanced multiscale, coarse-grained approaches [12][13][14][15]. This is our future research direction.

Figure 1 .

 1 Figure 1. Sketches of the shockwave induced bubble collapse simulation method. The rectangular prism box with lengths (x, y, z) contains two membranes and a bubble solvated in water. A slab of waters with a volume of x × y × ds is selected for the shockwave generation. A bubble with a diameter d b is placed at a distance of d mb above the membrane. The distance between the water slab and the bubble is d bs .

  I and bubble size d b . The other parameters d mb , d bs , d s are merely parameters defined only in the simulation. Thus, in this study, we only choose appropriate values, d mb = 1 nm, d bs = 2 nm and d s = 3 nm for all simulations and values of I and d b are varied to obtain different responses in the membranes. We note that two membranes are used in simulations but as explained above we are only interested in the lower membrane [Fig.1] whose results are presented below. For each simulation, ten trajectories are carried out starting from different initial structures and results are presented as an ensemble average over all trajectories.

Figure 2 .

 2 Figure 2. The profiles along the z-axis of the temperature (A), normal pressure component (B), mass density of water (C) and mass density of membrane (D) obtained at 0 ps (black), 1.25 ps (red), 3 ps (green) and 50 ps (blue). Shown are results obtained from the simulation of the normal membrane using a shockwave impulse I = 21 mPa.s and a 10 nm bubble.

Figure 3 .

 3 Figure3. The velocity field (in km/s) of atoms projected on the (x,z) plan (left panels), the snapshot of the system (middle panels) and the kinetic energy (in kJ/mol) map of atoms projected on the (x,y) plan (right panels). Shown are results at 1.25, 3 and 50 ps for the normal membrane simulated with a shockwave impulse of 21 mPa and a 10 nm in diameter bubble. In the snapshots, the water and membrane are shown in blue and red colours, respectively. The white area in the panel I indicates the pore.

Figure 4 .

 4 Figure 4. Time dependence of the average temperature (upper panel) and normal pressure (lower panel) across membrane surface along the z direction. Shown are results obtained from the simulation of the normal membrane using a shockwave impulse I = 21 mPa.s and a 10 nm bubble.

Figure 5 .

 5 Figure 5. Snapshots taken at t = 3 ps from the simulation of the normal (A) and cancer (B) membranes with a bubble having a diameter of 10 nm and a shockwave impulse of 15 mPa.s. The phosphate atoms of lipids are shown in orange and other atoms are in cyan. The bending values of the membranes are indicated.

Figure 6 .

 6 Figure 6. Time-evolution of the sum lipid order parameter S lipid of the normal (black) and cancer (red) membranes obtained from a simulation using a shockwave impulse of 21 mPa.s and a bubble with a diameter of 10 nm.

Figure 7 .

 7 Figure 7. Time evolution of the pore area ∆S induced by the collapse of a 10 nm bubble under different shockwave impulses I = 15 • • • 25 mPa.s. Show are results the normal (black lines) and cancer (red lines) membranes. The average data (over 10 trajectories) are shown by the solid thick lines and the error bars are shown by the dashed lines.

Figure 8 .

 8 Figure 8. Time evolution of the pore area induced by the collapse of bubbles with different radius d b under a shockwave impulse I = 23 mPa.s. Show are results of the normal (black lines) and cancer (red lines) membranes. The average data (over 10 trajectories) are shown by the solid thick lines and the error bars are shown by the dashed lines.

Figure 9 .

 9 Figure 9. Time evolution of the pore area in the normal (black line) and cancer (red line) membranes under equilibrium condition at 300 K and 1 bar. The initial pores at t = 0 ns were perforated by a shockwave impulse of 25 mPa.s and a 10 nm bubble.

Figure 10 .

 10 Figure 10. The order parameters averaged over the two tails sn-1 and sn-2 of four lipid types DOPC, DOPE, DOPS and PSM of the normal (left panels) and cancer (right panels) membranes. Carbon atom numbers increase in the direction of the tail termini. Shown are results before (black) and after (red) shockwave induced bubble collapse.

DFigure 11 .

 11 Figure 11. Snapshots taken at t = 3 ps from the simulation of the normal (A) and cancer (B) membranes without bubble but with a shockwave impulse of 23 mPa.s. The phosphate atoms of lipids are shown in orange and other atoms are in cyan.

VI. ACKNOWLEDGMENTS

This work has been supported by the Department of Science and Technology at Ho Chi Minh City, Vietnam (grant 13/2020/HD-PQTKHCN), the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 101034407 and the CNRS. Computational support from the IDRIS, CINES, TGCC centers (project A0110712955) are acknowledged.

DATA AVAILABILITY The data that support the findings of this study are available from the corresponding author upon reasonable request.


Benoit Rousseau

Salih Ouchtout

Jérôme Delmas

Christophe Le Bozec

Gwenaell Biotteau

Steven Le Corre

Jerôme Vicente

Ludovic Charpentier

Cindy Schick

Eric Louradour

Jerôme Delmas

Christophe Chaput

Finite Element Solver for Thermo-Hydraulic Optimization of Architected Ceramics (NExT I-Site industrial partnership program: OPTHOCERA): project overview

come

Interest of high-temperature concentrated solar heat

• Solar energy →abundant : 1% of the 3.8×10 24 joules of solar energy reaching the earth's surface every year → 40 times the electricity generated globally • But dilute, geographically non uniform, temporally intermittent, variable • Concentration 2) of solar radiation (paraboloid dishes and solar towers) →1,000°C ( with incident fluxes > 1MW/m 2 )

2) Romero, M., Steinfeld, A., 2012. Concentrating solar thermal power and thermochemical fuels. Energy Environ. Sci. 5, 9234-9245

• Development of high-temperature solar receivers with high solar-to-heat conversion efficiency and service life of 10,000 cycles under harsh conditions → thermal, mechanical, material challenges , 2021, 231, 111336 Luque et al., Solar Energy, 2018, 174, 342 Hoes et al., Energy Technol. 2019, 7, 1900484 (discrete scale model)

( ) ( ) ( ) ∫ ∫ ∫ ∞ ∞ Ω Ω - = ′ ′ ∇ 0 4 0 4 ν κ ν π κ π ν ν ν ν d d , s I d T I s q . a b a  ( ) ( ) ( ) ( ) s q s T k . s t T c cond p ′ ′ - ∇ ∇ = ∂ ∂  ρ � Ω. 𝛻𝛻𝐼𝐼 𝜈𝜈 𝑠𝑠, � Ω + 𝜅𝜅 𝑎𝑎 + 𝜅𝜅 𝑠𝑠 -𝜅𝜅 𝑠𝑠 ∫ 4𝜋𝜋 𝐼𝐼 𝜈𝜈 𝑠𝑠, � Ω 𝑃𝑃 � Ω ′ → � Ω 𝑑𝑑 � Ω ′ -
𝜀𝜀𝑛𝑛 2 𝜎𝜎 𝐵𝐵 𝑇𝑇 𝑆𝑆-𝑉𝑉 𝜕𝜕 4 -𝜀𝜀 ∫ � Ω.𝑛𝑛 𝑠𝑠 >0 𝐼𝐼 � Ω. 𝑛𝑛 𝑠𝑠 𝑑𝑑 � Ω = 0 � Ω. 𝛻𝛻𝐼𝐼 𝑠𝑠, � Ω = 0 𝛻𝛻. 16𝑛𝑛 2 𝜎𝜎 𝐵𝐵 3 𝜅𝜅 𝑎𝑎 + 𝜅𝜅 𝑠𝑠 𝑇𝑇 3 𝛻𝛻𝑇𝑇 = 𝛻𝛻. 𝑘𝑘 𝑅𝑅 𝛻𝛻𝑇𝑇

Rosseland model Homogeneous Equivalent Model (HEM)

(continous scale model)

� Ω. 𝛻𝛻𝐼𝐼 𝑠𝑠, � Ω + 𝜅𝜅 𝑎𝑎 + 𝜅𝜅 𝑠𝑠 -𝜅𝜅 𝑠𝑠 ∫ 4𝜋𝜋 𝐼𝐼 𝑠𝑠, � Ω 𝑃𝑃 � Ω ′ → � Ω 𝑑𝑑 � Ω ′ -𝜅𝜅 𝑎𝑎 1 𝜋𝜋 𝜎𝜎 𝐵𝐵 𝑛𝑛 2 𝑇𝑇 4 = 0 𝑃𝑃 � Ω ′ → � Ω = 1 4𝜋𝜋 1 -𝑔𝑔 2 1 + 𝑔𝑔 2 -� Ω. � Ω ′ � 3 2 𝜅𝜅 𝑒𝑒 = 𝜅𝜅 𝑎𝑎 + 𝜅𝜅 𝑠𝑠 = 4.8 1 -𝑝𝑝 𝑑𝑑 𝑛𝑛𝑛𝑛𝑛𝑛 𝜔𝜔 = 𝜅𝜅 𝑠𝑠 𝜅𝜅 𝑠𝑠 + 𝜅𝜅 𝑎𝑎 𝜅𝜅 𝑎𝑎 = 𝜀𝜀𝜅𝜅 𝑒𝑒

  surface • Low thermal emittance • High thermal shock resistance New concepts based on additive manufacturing coupled with radiative transfer design → Improvement of the volume propagation of concentrated solar light, reduction of thermal emission Heisel et al., Solar Energy Materials and Solar cells

  𝜅𝜅 𝑎𝑎 𝐼𝐼 𝑃𝑃𝜈𝜈 = 0 Methods for solving the steady-state RTE Avila-Marin et al., Renewable and Sustainable Energy Reviews, 2019Badri et al., Journal of Computational Physics, 2018, 360, 74-92 Badri et al. Journal of Quantitative Spectroscopy and Radiative Transfer 212 (2018): 59-74 Fast, robust, accurate parallelized vectorial FEM solver Effective radiative properties → porosity, nominal pore diameter and emissivity Local radiative properties from the bulk Detailed Simulation (DS)

𝑑𝑑Modification

  Use of the FEM solver for analyzing the coupling between conduction and radiation within a 1D slab sandwiched by a hot and a cold wall Viskanta & Grosh JHT (1962) • 1D layer of translucent conducting-radiating medium between tow parallel black walls at temperature T 1 et T 2 with θ = T/ T 1 et θ 2 = T 2 / T 1 Planck or Stark number Pure conductive regime Dominant radiative regime 𝑁𝑁 = 𝑘𝑘 𝑐𝑐𝑛𝑛𝑛𝑛𝑐𝑐 𝜅𝜅 𝑎𝑎 + 𝜅𝜅 𝑠𝑠 4𝜎𝜎 𝐵𝐵 𝑇𝑇 1 3 For a porous ceramic with absorption and scattering phenomena 𝑁𝑁 = 𝑘𝑘 𝑐𝑐𝑛𝑛𝑛𝑛𝑐𝑐 𝜅𝜅 𝑎𝑎 4𝜎𝜎 𝐵𝐵 Comparison between HEM model/exact solution from Viskanta θ 2 = 0.5 ; T 1 =2,222 K, T 2 =1,111 K, 𝜅𝜅 𝑎𝑎 = 2,625 m -1 , 𝜅𝜅 𝑠𝑠 = 0 𝑁𝑁 → 0 • Dominant radiative regime • Increase of T near the cold wall • 𝜏𝜏 • → @T=700 K , + 252 K • low 𝜅𝜅 𝑒𝑒 : for a given 𝑝𝑝, high 𝑑𝑑 𝑛𝑛𝑛𝑛𝑛𝑛 (low PPI) 𝑵𝑵 for real ceramic-based open cell foams (T=1,300 K) Interest for semi-transparent porous ceramics with high 𝑝𝑝 and high 𝑑𝑑 𝑛𝑛𝑛𝑛𝑛𝑛 331 860 tetrahedral elements for the solid phase and 1 753 680 elements for the fluid phase ; 62 processors, 128 directions 11 Comparison of the Rosseland↔HEM model and of the HEM↔DS model θ 2 = 0.5 ; T 1 =2,222 K, T 1 =1,111 K, 𝜅𝜅 𝑎𝑎 = 26,250 m -1 , 𝜅𝜅 𝑠𝑠 = 0 𝛽𝛽 𝑒𝑒𝑒𝑒𝑒𝑒 1,3 = 2𝑁𝑁 𝑃𝑃,𝛽𝛽 𝑒𝑒𝑒𝑒𝑒𝑒 2,4 → volume radiation effect more important 𝜔𝜔 𝑒𝑒𝑒𝑒𝑒𝑒 1,2 < 𝜔𝜔 𝑒𝑒𝑒𝑒𝑒𝑒 3,4 → volume radiation effect more important 𝜀𝜀 𝑒𝑒𝑒𝑒𝑒𝑒 1,2 < 𝜀𝜀 𝑒𝑒𝑒𝑒𝑒𝑒 3,4 → decreasing of the radiative losses at the front face Mora Monteros et al., Ceramics International 46 (2020) 2805without coating, case 2 8 mm, with coating, case 4 Absorption length of the incident collimated beam : effect of 𝑑𝑑 𝑛𝑛𝑛𝑛𝑛𝑛 without coating, case 2 8 mm, with coating, case 4 𝑁𝑁 𝑃𝑃,𝛽𝛽 𝑒𝑒𝑒𝑒𝑒𝑒 = 0,04422 𝑁𝑁 𝑃𝑃,𝛽𝛽 𝑒𝑒𝑒𝑒𝑒𝑒 = 0,02211 𝑥𝑥 𝑇𝑇,𝑛𝑛𝑎𝑎𝑥𝑥,𝑐𝑐1 𝑥𝑥 𝑇𝑇,𝑛𝑛𝑎𝑎𝑥𝑥,𝑐𝑐3 𝑥𝑥 𝑇𝑇,𝑛𝑛𝑎𝑎𝑥𝑥,𝑐𝑐2 𝑥𝑥 𝑇𝑇,𝑛𝑛𝑎𝑎𝑥𝑥,𝑐𝑐4 𝑥𝑥 𝐴𝐴,𝑛𝑛𝑎𝑎𝑥𝑥,𝑐𝑐3 𝑥𝑥 𝐴𝐴,𝑛𝑛𝑎𝑎𝑥𝑥,𝑐𝑐4 𝑥𝑥 𝐴𝐴,𝑛𝑛𝑎𝑎𝑥𝑥,𝑐𝑐2 𝑥𝑥 𝐴𝐴,𝑛𝑛𝑎𝑎𝑥𝑥,𝑐𝑐1 Impact 𝒅𝒅 𝒏𝒏𝒏𝒏𝒏𝒏 on the absorption of the incident radiative flux and the position of 𝑻𝑻 𝒏𝒏𝒎𝒎𝒎𝒎 Similar effect Radiation transport is stopped at the first rank of cell : after conductive (convective) heat transfer Solar-to-heat performances : OPTISOL bench Mey-Cloutier et al., Solar Energy 136 (2016) 226-235 Qualitative agreement but some experiments must still be explained 𝜂𝜂 𝑆𝑆𝑇𝑇,1 < 𝜂𝜂 𝑆𝑆𝑇𝑇,2 ↔𝑁𝑁 𝑃𝑃,𝛽𝛽 𝑒𝑒𝑒𝑒𝑒𝑒 ,1 > 𝑁𝑁 𝑃𝑃,𝛽𝛽 𝑒𝑒𝑒𝑒𝑒𝑒 ,2 𝜂𝜂 𝑆𝑆𝑇𝑇,4 < 𝜂𝜂 𝑆𝑆𝑇𝑇,3 but 𝑁𝑁 𝑃𝑃,𝛽𝛽 𝑒𝑒𝑒𝑒𝑒𝑒 ,3 > 𝑁𝑁 𝑃𝑃,𝛽𝛽 𝑒𝑒𝑒𝑒𝑒𝑒 ,4 𝜂𝜂 𝑆𝑆𝑇𝑇,4 < 𝜂𝜂 𝑆𝑆𝑇𝑇,3 ↔𝜔𝜔 𝑒𝑒𝑒𝑒𝑒𝑒,3 < 𝜔𝜔 𝑒𝑒𝑒𝑒𝑒𝑒,4

15 Numerical design approach coupled with additive manufacturing

	Numerical design of 3D digitalized porous architectures (regular, pore size gradient) Numerical design of 3D digitalized porous architectures (regular, pore size gradient)	Thermo-fluidic modelling→Tempe rature and flux fields (HEM or DS) (conduction, convection, radiation) Conductive-radiative modelling→Temper ature and flux fields (HEM and DS)	Optimization loop/textural constraints (solar-to heat efficiency, temperature profiles)
	TAMARYS network	Optimized Selected architectures → architectures elaboration with according to Planck additive number → manufacturing (stereolithography, binder jetting) elaboration by stereolithography	Performance tests : Performance tests : solar-to-heat solar-to-heat conversion conversion (OPTISOL bench, (OPTISOL bench, PROMES Odeillo, FRance). PROMES Odeillo, France).

150 mm Avila-Marin et al., Renewable and Sustainable Energy Reviews, 2019, 111 HEM : Homogeneous Equivalent Method → continuous scale approach DS : Detailed Simulation → pore scale approach

OPTHOCERA program (1 year) Radiative Transfer Equation : modelling of volume propagation of thermal and solar radiations •

  Macro-porous ceramics → semi-transparent media (absorption, scattering, emission) • Valid if the medium is homogeneous, continous and weakly spatially correlated

Chandrasekhar, Radiative Transfer, Dover, 1960

Conclusions & outlook

Thank you for your attention ! benoit.rousseau@univ-nantes.fr

Influence of the anisotropy factor

Key quantities

𝜅𝜅 𝑠𝑠

From macro-scale to meso-scale macro → meso Guévelou et al., Journal of Quantitative Spectroscopy and Radiative Transfer, 2017, 189, 329-338