Clifford W Fong
email: cwfong@internode.on.net

Mechanism, structure activity relationships of lipid peroxidation of cell membranes and brain protection for cerebral ischemia by Edaravone and Edaravone analogs -a quantum mechanical study

Keywords: Edaravone, Edaravone analogs, lipid peroxidation, structural activity relationships, brain protection, oxidative stress, free radicals, quantum mechanics

Mechanism, structure activity relationships of lipid peroxidation of cell membranes and brain protection for cerebral ischemia by Edaravone and Edaravone analogs -a quantum mechanical study

Clifford W Fong

Introduction

Oxidative stress which results from an imbalance between the biological effects caused by reactive oxygen species (ROS) and the body's ability to detoxify the ROS or repair the resulting damage. Changes in the normal redox state of cells can cause toxic effects through the production of peroxides and free radicals that damage all components of the cell, including proteins, lipids, and DNA. Oxidative stress is thought to be involved in the development of cancer, Parkinson's disease, Alzheimer's disease, atherosclerosis, heart failure, myocardial infarction, attention deficit hyperactivity disorder, diabetes mellitus, hypertension, and other human diseases. [START_REF] Yoshikawa | What Is Oxidative Stress?[END_REF][START_REF] Noda | Cancer and oxidative stress[END_REF][START_REF] Hegedűs | Redox control of cancer cell destruction[END_REF] Lipid peroxidation is a biomarker of oxidative stress. [START_REF] Ito | Measurement and Clinical Significance of Lipid Peroxidation as a Biomarker of Oxidative Stress: Oxidative Stress in Diabetes, Atherosclerosis, and Chronic Inflammation[END_REF] Lipid peroxidation is well known to be a chain reaction generally initiated by the formation of a free radical particularly from the allylic position of poly-unsaturated fatty acids (PUFA) in the lipid bilayer. [START_REF] Ha Hassan | Lipid Peroxidation End-Products as a Key of Oxidative Stress: Effect of Antioxidant on Their Production and Transfer of Free Radicals[END_REF][START_REF] Repetto | Lipid Peroxidation: Chemical Mechanism, Biological Implications and Analytical Determination[END_REF][START_REF] Niki | Lipid peroxidation: Mechanisms, inhibition, and biological effects[END_REF] The allylic (and particularly the bisallylic) free radical is a particularly stable species owing to resonance electron delocalization. The fatty acids in phospho-and glycolipids of cell membrane lipid bilayers usually contain about 16 to 18 carbon atoms and can be saturated or poly-unsaturated PUFAs. The degree of unsaturation has a major influence on membrane fluidity. PUFAs affect numerous physiological processes related to cardiovascular, neurological and immune function as well as cancer [START_REF] Benatti | Polyunsaturated fatty acids: biochemical, nutritional and epigenetic properties[END_REF]9]. The human brain contains about 25-30% PUFAs, and PUFA content in all brain structures is strongly dependent on the dietary PUFAs. [9a] PUFAs can be involved in the modulation of the physical properties of biological membranes via alteration of membrane lipid composition. [9b] Lipid peroxidation is a key mediator of many pathologies including inflammation, cancer, neurodegenerative disease, as well as ocular and kidney degeneration. [START_REF] Gaschler | Lipid peroxidation in cell death[END_REF][START_REF] Feng | Unsolved mysteries: How does lipid peroxidation cause ferroptosis?[END_REF][START_REF] Agmon | Modeling the effects of lipid peroxidation during ferroptosis on membrane properties[END_REF][START_REF] Ov Taso | Lipid peroxidation products and their role in neurodegenerative diseases[END_REF][START_REF] Da Butterfield | Evidence that amyloid beta-peptideinduced lipid peroxidation and its sequelae in Alzheimer's disease brain contribute to neuronal death[END_REF] Lipid peroxidation is also the key downstream feature of ferroptosis, a form of regulated iron dependent nonapoptotic cell death.

Antioxidants or free radical scavengers can be used as chemotherapies to ameliorate diseases caused by lipid peroxidation. Edaravone (Edv) is a free radical scavenger FDA approved in 2017 for the treatment of amyotrophic lateral sclerosis (ALS). It has been previously investigated for the treatment of ischemic stroke, reperfusion injury, and myocardial infarction as it possesses antioxidant and anti-apoptotic properties. It was approved in Japan in 2001 as a drug to treat acute phase cerebral infarction, and in 2015 it was approved for ALS. It can cross the bloodbrain barrier BBB to exert nootropic and neuroprotective effects. [START_REF] Watanabe | How is edaravone effective against acute ischemic stroke and amyotrophic lateral sclerosis?[END_REF][START_REF] Petrov | ALS Clinical Trials Review: 20 Years of Failure. Are We Any Closer to Registering a New Treatment[END_REF][17][START_REF] Fong | Improved Edaravone delivery to the brain and crossing the blood brain barrier: using quantum mechanics[END_REF] Edv is also effective in ameliorating lipid peroxidation directly. [START_REF] Fong | Improved Edaravone delivery to the brain and crossing the blood brain barrier: using quantum mechanics[END_REF] The mechanisms by which Edv is effective against lipid peroxidation is unknown, direct scavenging of the hydroxyl radical and or H 2 O 2 involvement has been postulated. [START_REF] Ha Hassan | Lipid Peroxidation End-Products as a Key of Oxidative Stress: Effect of Antioxidant on Their Production and Transfer of Free Radicals[END_REF][START_REF] Repetto | Lipid Peroxidation: Chemical Mechanism, Biological Implications and Analytical Determination[END_REF][START_REF] Niki | Lipid peroxidation: Mechanisms, inhibition, and biological effects[END_REF][START_REF] Watanabe | How is edaravone effective against acute ischemic stroke and amyotrophic lateral sclerosis?[END_REF][START_REF] Ahmadinejad | Molecular Mechanisms behind Free Radical Scavengers Function against Oxidative Stress[END_REF] In one study the cytoprotective effects of Edv were investigated in SHSY5Y dopaminergic neurons, and the modification of various proteins suggested that Edv scavenges H 2 O 2 and protects cells against oxidative stress via direct induction of Peroxiredoxin-2 and inhibition of apoptosis. [START_REF] Ms Jami | Edaravone leads to proteome changes indicative of neuronal cell protection in response to oxidative stress[END_REF] The antioxidative effect of Edv in cerebral infarction has been documented in human clinical trials, and oxidized low-density lipoprotein (OxLDL) has been shown to be a plasma biomarker of the therapeutic efficacy of Edv patients with cerebral infarction. The antioxidant effect of Edv in ALS patients has also been demonstrated. [START_REF] Watanabe | How is edaravone effective against acute ischemic stroke and amyotrophic lateral sclerosis?[END_REF] Study objectives: Using a previously documented linear free energy relationship method, evaluate the structure activity of a wide series of Edaravone analogs and derive a mechanistic predictive model for the in vitro lipid peroxidation of cell membranes. Investigate whether the same methodology is applicable to in vivo brain protection.

Results

Watanabe et al [START_REF] Watanabe | Structure-activity relationship of 3-methyl-1-phenyl-2-pyrazolin-5-one (edaravone)[END_REF] have previously reported the free radical scavenging effects of Edaravone and a wide range of analogs (over 50) and their antioxidant effects on lipid peroxidation. The lipid peroxidation-inhibitory activity was measured by the inhibition of auto-oxidation of the Wistar rat brain homogenate according to the thiobarbituric acid test. A more limited in vivo study was also carried out.

This study uses the data from Watanabe's study, and applies our previously developed linear free energy relationship method to the in vitro inhibition of lipid peroxidation and the in vivo study.

We have previously described a model that has been shown to apply to a wide range of drug transport, binding, metabolic and cytotoxicity properties of cells and tumours (Equation 1).

The model is based on establishing linear free energy relationships between the four drug properties and various biological processes. Equation 1 has been previously applied to passive and facilitated diffusion of a wide range of drugs crossing the blood brain barrier, [START_REF] Fong | Permeability of the Blood-Brain Barrier: Molecular Mechanism of Transport of Drugs and Physiologically Important Compounds[END_REF] the active competitive transport of tyrosine kinase inhibitors by the hOCT3, OATP1A2 and OCT1 transporters, and cyclin-dependent kinase inhibitors and HIV-1 protease inhibitors. The model also applies to PARP inhibitors, the anti-bacterial and anti-malarial properties of fluoroquinolones, and active organic anion transporter drug membrane transport, and some competitive statin-CYP enzyme binding processes. There is strong independent evidence from the literature that ΔG desolvation , ΔG lipophilicity , the dipole moment and molecular volume are good inherent indicators of the transport or binding ability of drugs. [START_REF] Fong | Permeability of the Blood-Brain Barrier: Molecular Mechanism of Transport of Drugs and Physiologically Important Compounds[END_REF][START_REF] Fong | The extravascular penetration of tirapazamine into tumours: a predictive model of the transport and efficacy of hypoxia specific cytotoxic analogues and the potential use of cucurbiturils to facilitate delivery[END_REF][START_REF] Fong | Screening anti-colorectal cancer drugs: free radical chemotherapy[END_REF][START_REF] Fong | Free radical anticancer drugs and oxidative stress: ORAC and CellROXcolorectal cancer cells by quantum chemical determinations[END_REF][START_REF] Fong | The role of free radicals in the effectiveness of anti-cancer chemotherapy in hypoxic ovarian cells and tumours[END_REF][START_REF] Fong | Free radicals in chemotherapy induced cytotoxicity and oxidative stress in triple negative breast and ovarian cancers under hypoxic and normoxic conditions[END_REF][START_REF] Fong | Role of stable free radicals in conjugated antioxidant and cytotoxicity treatment of triple negative breast cancer[END_REF][START_REF] Fong | Toxicology of platinum anticancer drugs: oxidative stress and antioxidant effect of stable free radical Pt-nitroxides[END_REF][START_REF] Fong | The effect of desolvation on the binding of inhibitors to HIV-1 protease and cyclin-dependent kinases: Causes of resistance[END_REF][START_REF] Fong | Physiology of ionophore transport of potassium and sodium ions across cell membranes: Valinomycin and 18-Crown-6 Ether[END_REF][START_REF] Fong | Statins in therapy: Understanding their hydrophilicity, lipophilicity, binding to 3hydroxy-3-methylglutaryl-CoA reductase, ability to cross the blood brain barrier and metabolic stability based on electrostatic molecular orbital studies[END_REF][START_REF] Fong | Predicting PARP inhibitory activity -A novel quantum mechanical based model[END_REF][START_REF] Fong | A novel predictive model for the anti-bacterial, anti-malarial and hERG cardiac QT prolongation properties of fluoroquinolones[END_REF][START_REF] Fong | Statins in therapy: Cellular transport, side effects, drug-drug interactions and cytotoxicity -the unrecognized role of lactones[END_REF][START_REF] Fong | Drug discovery model using molecular orbital computations: tyrosine kinase inhibitors[END_REF] Equation 1:

Transport or Binding or Cytotoxicity = ΔG desolv,CDS + ΔG lipo,CDS + Dipole Moment + Molecular Volume

Eq 1 uses the free energy of water desolvation (ΔG desolv,CDS ) and the lipophilicity free energy (ΔG lipo,CDS ) where CDS represents the non-electrostatic first solvation shell solvent properties. CDS may be a better approximation of the cybotactic environment around the drug approaching or within the protein receptor pocket, or the cell membrane surface or the surface of a drug transporter, than the bulk water environment outside the receptor pocket or cell membrane surface. The CDS includes dispersion, cavitation, and covalent components of hydrogen bonding, hydrophobic effects. Desolvation of water from the drug (ΔG desolv,CDS ) before binding in the receptor pocket is required, and hydrophobic interactions between the drug and protein (ΔG lipo,CDS ) is a positive contribution to binding. The lipophilicity ΔG lipo,CDS is calculated from the solvation energy in n-octane or n-octanol. In some biological processes, where biological reduction may be occurring, and the influence of molecular volume is small, the reduction potential (electron affinity) has been included in place of the molecular volume. In other processes, the influence of some of the independent variables is small and can be eliminated to focus on the major determinants of biological activity.

We have recently used this model to develop a predictive model of the transport and efficacy of hypoxia specific cytotoxic analogues of tirapazamine and the effect on the extravascular penetration of tirapazamine into tumours. [START_REF] Fong | The extravascular penetration of tirapazamine into tumours: a predictive model of the transport and efficacy of hypoxia specific cytotoxic analogues and the potential use of cucurbiturils to facilitate delivery[END_REF] It was found that the multiparameter model of the diffusion, antiproliferative assays IC 50 and aerobic and hypoxic clonogenic assays for a wide range of neutral and radical anion forms of tirapazamine (TPZ) analogues showed: (a) extravascular diffusion is governed by the desolvation, lipophilicity, dipole moment and molecular volume, similar to passive and facilitated permeation through the blood brain barrier and other cellular membranes, (b) hypoxic assay properties of the TPZ analogues showed dependencies on the electron affinity, as well as lipophilicity and dipole moment and desolvation, similar to other biological processes involving permeation of cellular membranes, including nuclear membranes, (c) aerobic assay properties were dependent on the almost exclusively on the electron affinity, consistent with electron transfer involving free radicals being the dominant species.

The model (eq 2) has also been recently applied to triple negative breast and ovarian cancers where transient and stable free radicals are involved in the cytotoxic oxidative stress processes. The electron affinity of the various drugs, along with the water desolvation, lipophilicity and dipole moment, has been shown to be an important predictor of cytotoxic efficacy. [START_REF] Fong | The role of free radicals in the effectiveness of anti-cancer chemotherapy in hypoxic ovarian cells and tumours[END_REF][START_REF] Fong | Free radicals in chemotherapy induced cytotoxicity and oxidative stress in triple negative breast and ovarian cancers under hypoxic and normoxic conditions[END_REF][START_REF] Fong | Role of stable free radicals in conjugated antioxidant and cytotoxicity treatment of triple negative breast cancer[END_REF][START_REF] Fong | Toxicology of platinum anticancer drugs: oxidative stress and antioxidant effect of stable free radical Pt-nitroxides[END_REF] In our recent study of ORAC and CellROX free radical anticancer drugs and oxidative stress in colorectal cancer cells [START_REF] Fong | Screening anti-colorectal cancer drugs: free radical chemotherapy[END_REF][START_REF] Fong | Free radical anticancer drugs and oxidative stress: ORAC and CellROXcolorectal cancer cells by quantum chemical determinations[END_REF] we found that eq 2 was also applicable.

The antioxidant capacity (or chemical reactivity) of the drugs can be assessed by the HOMO-LUMO energy gap, since we have recently shown that this gap is linearly related to the LD 50 toxicity of various drugs which are involved in oxidative stress processes in some cancers. [START_REF] Fong | Screening anti-colorectal cancer drugs: free radical chemotherapy[END_REF][START_REF] Fong | Role of stable free radicals in conjugated antioxidant and cytotoxicity treatment of triple negative breast cancer[END_REF][START_REF] Fong | Toxicology of platinum anticancer drugs: oxidative stress and antioxidant effect of stable free radical Pt-nitroxides[END_REF] Eq 2 can then be represented by substituting the electron affinity by the HOMO-LUMO energy gap for analogs where no radical species is formed. Equation 2.

Oxidative Stress or Oxidative Properties = ΔG desolv,CDS + ΔG lipo,CDS + Dipole Moment + Electron Affinity (or HOMO-LUMO gap) Evaluation methodology

Given the successful transport, binding or cytotoxicity model as described by eq 1 and 2 for wide range of drugs, the method used was to calculate the molecular parameters for the lipid peroxidation-inhibitory activity (IC 50 μM) of Edv and Edv analogs taken from Watanabe [START_REF] Watanabe | Structure-activity relationship of 3-methyl-1-phenyl-2-pyrazolin-5-one (edaravone)[END_REF]:

(1) the free energy of water desolvation (ΔG desolv,CDS ), (2) the lipophilicity free energy (ΔG lipo,CDS ) in n-octane, (3) the dipole moment in water, (4) the molecular volume in water, and (5) the adiabatic electron affinity, AEA, (or reduction potential) in water [START_REF] Repetto | Lipid Peroxidation: Chemical Mechanism, Biological Implications and Analytical Determination[END_REF] or the HOMO-LUMO energy gap in water. These values are shown in Table 1. In vivo results for the brainprotective activity (recovery of the electroencephalogram and survival time) in cerebral ischemia in rats for the 9 Edv analogs which showed brain protective activity are also shown in Table 2 .

Stepwise multiple regression was then applied to seek out which drug molecular properties had the largest effect on the in vitro IC 50 values for Edv and Edv analogs shown in Table 1 and for the in vivo EEG or survival times in Table 2.

Stepwise multiple regression analysis showed that eq 3 was the best description when tested against ΔG desolv,CDS, ΔG lipo,CDS, Dipole Moment, Molecular Volume, Adiabatic Electron Affinity or the HOMO-LUMO energy gap. It was also clear that a number of analogs were statistical outliers, particularly where the phenyl ring (X) was ortho substituted (methyl, methoxy, chloro, hydroxyl, compounds 22,37,29 and 33), and where R4 was CH 2 -imadazoyl (compound 56). For compound 47, R2 = thiophenyl, it was found that the C-S bond was very labile upon electron attachment, indicating unusual reducing or antioxidant reactivity. Compound 18, R1 = 4Pyridyl, was a clear outlier, and was tested as the protonated species, but this species was not the cause of the outlying behavior. Compound 37, X= 4Hydroxymethyl was also a large outlier. It can be concluded that some steric effect is involved in lipid peroxidation causing the majority of outliers. It has been previously noted that Edv is an excellent scavenger of most studied radical species (eg hydroxyl, alkoxyl, peroxyl etc) except the phenyl radical, which has been attributed to steric hindrance. [START_REF] Watanabe | How is edaravone effective against acute ischemic stroke and amyotrophic lateral sclerosis?[END_REF] (a) In vitro results The corresponding where AEA is substituted for the HOMO-LUMO gap was shown to be a less accurate relationship, owing to a weaker relationship between IC 50 and AEA compared to IC 50 and HOMO-LUMO energy gap. Molecular volume was poorly related to IC 50 for this data set.

Eq 3 Inhibition of lipid peroxidation IC 50 for Edv and 35 analogs in

(b) In vivo results

(i) Stepwise regression analysis tested against

Discussion (a) In vitro results:

The mechanism by which antioxidant drugs such as Edv and the Edv analogs govern lipid peroxidation depends primarily on how such drugs enter the cell lipid membrane from an aqueous bulk environment which has a dielectric of 80. Phospholipid bilayers typically have a dielectric gradient, that typically vary from 1-2 (at the hydrophobic core of the bilayer) to about 4-30 at the headgroup/ester region.

[37] For passive diffusion to occur, we have previously shown that desolvation of the drug is required to allow entry into the phospholipid bilayer. The degree to which such drugs can then further permeate deeper into the hydrophobic core of the bilayer depends on the lipophilicity of the drugs as well as the dielectric gradient and the molecular volume of the drugs. [START_REF] Fong | Permeability of the Blood-Brain Barrier: Molecular Mechanism of Transport of Drugs and Physiologically Important Compounds[END_REF][START_REF] Fong | The extravascular penetration of tirapazamine into tumours: a predictive model of the transport and efficacy of hypoxia specific cytotoxic analogues and the potential use of cucurbiturils to facilitate delivery[END_REF][START_REF] Fong | Screening anti-colorectal cancer drugs: free radical chemotherapy[END_REF][START_REF] Fong | Free radical anticancer drugs and oxidative stress: ORAC and CellROXcolorectal cancer cells by quantum chemical determinations[END_REF][START_REF] Fong | The role of free radicals in the effectiveness of anti-cancer chemotherapy in hypoxic ovarian cells and tumours[END_REF][START_REF] Fong | Free radicals in chemotherapy induced cytotoxicity and oxidative stress in triple negative breast and ovarian cancers under hypoxic and normoxic conditions[END_REF][START_REF] Fong | Role of stable free radicals in conjugated antioxidant and cytotoxicity treatment of triple negative breast cancer[END_REF][START_REF] Fong | Toxicology of platinum anticancer drugs: oxidative stress and antioxidant effect of stable free radical Pt-nitroxides[END_REF][START_REF] Fong | The effect of desolvation on the binding of inhibitors to HIV-1 protease and cyclin-dependent kinases: Causes of resistance[END_REF][START_REF] Fong | Physiology of ionophore transport of potassium and sodium ions across cell membranes: Valinomycin and 18-Crown-6 Ether[END_REF][START_REF] Fong | Statins in therapy: Understanding their hydrophilicity, lipophilicity, binding to 3hydroxy-3-methylglutaryl-CoA reductase, ability to cross the blood brain barrier and metabolic stability based on electrostatic molecular orbital studies[END_REF][START_REF] Fong | Predicting PARP inhibitory activity -A novel quantum mechanical based model[END_REF][START_REF] Fong | A novel predictive model for the anti-bacterial, anti-malarial and hERG cardiac QT prolongation properties of fluoroquinolones[END_REF][START_REF] Fong | Statins in therapy: Cellular transport, side effects, drug-drug interactions and cytotoxicity -the unrecognized role of lactones[END_REF][START_REF] Fong | Drug discovery model using molecular orbital computations: tyrosine kinase inhibitors[END_REF] These factors are represented by eq 1 and shown schematically in Figure 1.

Compared to ions, water molecules actually have a relatively large permeability through the bilayer, as evidenced by osmotic swelling. The anomalously large permeability of water through bilayers is still not completely understood and continues to be controversial. Lipids with an unsaturated tail disrupt the packing of those with only saturated tails in the bilayer. The resulting bilayer has more free space and is, as a consequence, more permeable to water and other small molecules. It is well known that the diffusion of water through phospholipid bilayers is fastest just outside the membrane surface (reduced by a factor of ca 2 compared to the bulk water environment), and decreases rapidly further into the more hydrophobic regions of the bilayer by an order of magnitude. About six water molecules are tightly bound to lipid molecules. In PUFA molecules in lipid bilayers, water molecules overlap with the unsaturated regions. [START_REF] Hansen | A Miskowiec Diffusion of water and selected atoms in DMPC lipid bilayer membranes[END_REF][START_REF] Saiz | Structural properties of a highly polyunsaturated lipid bilayer from molecular dynamics simulations[END_REF][START_REF] Huster | Water permeability of polyunsaturated lipid membranes measured by 17 O NMR[END_REF] Inhibition of lipid peroxidation involves the antioxidant drug which has penetrated the bilayer (as represented by eq 1) then interacting with an allylic PUFA free radical site produccd (presumably) by a hydroxyl free radical (as shown in Figure 2). In the present study, Watanabe's experimental results [START_REF] Watanabe | Structure-activity relationship of 3-methyl-1-phenyl-2-pyrazolin-5-one (edaravone)[END_REF] were from autooxidation of rat brain homogenate. The reactivity between an antioxidant Edv or Edv analog and an allylic free radical PUFA can be described by the AEA or HOMO-LUMO gap of the antioxidants, both of which are measures of the ease of reaction with a PUFA free radical. Eq 2 can then be used to describe the inhibition process, particularly where the influence of molecular volume of the drugs is small.

It has been found that after the incorporation of PUFAs into human heptoma cell line, HepG2, cellular lipid peroxidation correlated with the PUFA conformation in the phospholipid membrane. [START_REF] Araseki | Lipid Peroxidation of a Human Hepatoma Cell Line (HepG2) after Incorporation of Linoleic Acid, Arachidonic Acid, and Docosahexaenoic Acid[END_REF] For example, NMR analysis of docosahexaenoic acid -phosphatidylcholine liposomes, has shown the possibility of a wide variety of specific conformation of the PUFA docosahexaenoic acid, including back-bended (hairpin-like), helical, and angle-iron conformations. [START_REF] Everts | 1 H and 13 C NMR of multilamellar dispersions of polyunsaturated phospholipids[END_REF][START_REF] Se Feller | Polyunsaturated fatty acids in lipid bilayers intrinsic and environmental contributions to their unique physical properties[END_REF] Eq 3 then quantifies the relationship between the inhibition of lipid peroxidation, IC 50 , and drug desolvation, lipophilicity, the dipole moment (properties representing how the drug penetrates the bilayer) and the reactivity of the drug towards a presumably allylic free radical PUFA (as represented by the HOMO-LUMO gap). It was noted that steric factors contributed to several clear outliers from eq 3, namely where the phenyl ring (X) was ortho substituted (methyl, methoxy, chloro, hydroxyl, compounds 22,37,29 and 33), and where R4 was CH 2 -imadazoyl (compound 56). These outliers suggest that the origin of these outliers is due to some conformational factors of the allylic free radical PUFA which prevents these outlying antioxidant drugs from closely interacting with the reaction site of the allylic free radical PUFA (which can occur for the drugs described by eq 3).

Further insights into the mechanism can be gleaned from antioxidant properties of EdvC18 (see Figure 3, where R 2 =linear C 18 H 37 , R 3 =R 4 =H), which was studied by Minelli. [START_REF] Minnelli | Characterization and Antioxidant Properties of a New Lipophilic Derivative of Edaravone[END_REF] Table 1 shows that EdvC18 has a larger desolvation and lipophilicity (and molecular volume as expected) than comparable drugs in Table 1 (eg compound 44, R 2 =isobutyl). Minelli showed that EdvC18 showed greater localization in the bilayer with respect to Edv, resulting in enhanced lipid peroxidation-inhibitory activity when the oxidation is induced in the lipid phase. It is noted that compound 43 (R 2 =Methyl) has a higher IC 50 =20.6 than compound 44 which has the longer isobutyl side chain, IC 50 =5.7. These data are consistent with Minelli's results. As we have shown there are steric interactions between the allylic free radical PUFA and the antioxidant inhibitor for R1 (ortho-phenyl substituents) and R4, it can be concluded that bulky substituents on R2 or R3 may not exert significant steric hindrance to the interaction of the inhibitor with the allylic PUFA. Watanabe [START_REF] Watanabe | How is edaravone effective against acute ischemic stroke and amyotrophic lateral sclerosis?[END_REF] has discussed the importance of lipophilicity to Edv activity, particularly the R1 (phenyl), R2 and R3 substituents, but this study shows that R4 substituents may create steric activity problems.

Inhibition of lipid peroxidation involves the antioxidant drug which has penetrated the bilayer (as represented by eq 1) then interacting with an allylic PUFA free radical site produccd (presumably) by a hydroxyl free radical. Watanabe's experimental results [START_REF] Watanabe | Structure-activity relationship of 3-methyl-1-phenyl-2-pyrazolin-5-one (edaravone)[END_REF] were for the autooxidation of rat brain homogenate by Edv and the various Edv analogs in the keto form. A critical issue is the form of the Edv and its analogs within the lipid bilayer, since it has been demonstrated that Edv exerts it strong antioxidative effect in aqueous media mainly by the anionic form derived from the emol tautomer. [START_REF] Watanabe | How is edaravone effective against acute ischemic stroke and amyotrophic lateral sclerosis?[END_REF][START_REF] Pérez-González | On the Outstanding Antioxidant Capacity of Edaravone Derivatives through Single Electron Transfer Reactions[END_REF] However the keto tautomer of Edv is more hydrophobic than the enol tautomer, and would be expected to permeate better than the enol form within the lipid bilaye. Formation of the anionic form of Edv from any enol tautomer of Edv that did penetrate the lipid layer would be minimal in that largely non-aqueous environment. [START_REF] Fong | Improved Edaravone delivery to the brain and crossing the blood brain barrier: using quantum mechanics[END_REF] The anionic form of Edv would not permeate the lipid bilayer owing to a large desolvation requirement that prevents charged species from permeating the blood brain barrier. [START_REF] Fong | Permeability of the Blood-Brain Barrier: Molecular Mechanism of Transport of Drugs and Physiologically Important Compounds[END_REF] It is also known that Edv is a much poorer HOO• scavenger in lipid media than it is in aqueous media.

[46]

(b) In vivo results:

The in vivo data is complicated by the fact that Edv and its analogs are known to exist in various tautomeric form, particularly the enol and keto forms, with the former known to easily deprotonate to form an enolate anion at physiological pH levels. Edv can exist as the neutral and anionic forms with their relative concentrations at pH 7.4 are 28.5% and 71.5% respectively in water. [START_REF] Watanabe | How is edaravone effective against acute ischemic stroke and amyotrophic lateral sclerosis?[END_REF] The in vivo brain protection from the keto Edv tautomers for cerebral ischemia in rats is described in eqs 4 and 5 for survival times and EEG recovery times, with the corresponding enol tautomers described in eq 6 and 7. Eqs 6 and 7 are stronger statistically than eqs 4 and 5, possibly indicating that the enolate tautomers are the likely active species providing brain protection in rats. Brain protection is better indicated by EEG recovery times (which directly measure brain function) rather than survival times (which may include toxicity in other organs as well as the brain).

Eqs 5 and 7 for EEG recovery are comparable, both showing a negative dependency on desolvation (or a positive dependency on drug solvation, which may indicate a dependency on bioavailablity). Eq 5 shows a reversed dependency on lipophilicity for the keto tautomers versus the enol tautomers, which seems anomalous as a greater lipophilicity would be expected to enhance permeability of the blood brain barrier (BBB). Eq 7 has the expected dependency on lipophilicity of the enol tautomers. It is noted that eq 3 for the in vitro results for the keto forms of Edv and its analogs does show a positive dependency on lipophilicity. A substantial difference may be expected between the in vitro conditions where the keto tautomeric form of Edv and its analogs were used as a starting point, and the longer equilibration times involved in the in vivo conditions where substantial formation of the enol tautomeric form can occur in the rat blood and brain. This situation offer a possible explanation for the results in eq 7, ie that the in vivo EEG data is due to the enol form of Edv and its analogs being the dominant species rather than the keto forms. Keto and enol tautomerism of Edv and its analogs has been well documented, and the lipophilicity of the enol form is comparable to the keto form, and so would be expected to have similar permeability through lipd bilayer membranes. [START_REF] Watanabe | How is edaravone effective against acute ischemic stroke and amyotrophic lateral sclerosis?[END_REF][START_REF] Fong | Improved Edaravone delivery to the brain and crossing the blood brain barrier: using quantum mechanics[END_REF] Eqs 4 and 6 apply to whole animal toxic responses (or survival times) to the drugs, and eq 6 is the stronger correlation, probably indicating that HOMO-LUMO gap of the enol form of the Edv analogs is the dominant determinant driving animal survival time. We have recently shown the HOMO-LUMO energy gap is linearly related to the LD 50 toxicity of various drugs which are involved in oxidative stress processes in some cancers. [START_REF] Fong | Screening anti-colorectal cancer drugs: free radical chemotherapy[END_REF][START_REF] Fong | Role of stable free radicals in conjugated antioxidant and cytotoxicity treatment of triple negative breast cancer[END_REF][START_REF] Fong | Toxicology of platinum anticancer drugs: oxidative stress and antioxidant effect of stable free radical Pt-nitroxides[END_REF]

Conclusions

The in vitro inhibition of lipid peroxidation by Edaravone and a wide range of Edaravone analogs has been found to be linearly dependent upon the desolvation (ΔG desolv,CDS ), lipophilicity (ΔG lipo,CDS ), the dipole moment (properties representing how the drug penetrates the bilayer) and the reactivity of the drugs towards a presumably allylic free radical PUFA (as represented by the HOMO-LUMO gap). Steric factors where the phenyl ring (X) was ortho substituted and where R4 was CH 2 -imadazoyl have been identified. These steric effects are consistent with some conformational factors of the intermediate allylic free radical PUFA which prevents these antioxidant drugs from closely interacting with the reaction centre of the allylic free radical PUFA.

The in vivo brain protection for cerebral ischemia in rats exhibited by Edaravone and Edaravone analogs appears to best attributed to the enol tautomeric form rather than the keto tautomeric form. The EEG recovery time of rats is positively correlated with lipophilicity in the blood brain barrier and desolvation of the enol tautomeric form of Edaravone and its analogs. The animal survival time is strongly related to the HOMO-LUMO gap of the enol tautomeric form of Edaravone and its analogs after treatment with these drugs.

This in vitro and in vivo study is consistent with a mechanism that is dominated by (a) the penetration of the drug into the cell membrane which requires initially desolvation of the drug as it moves from the bulk aqueous environment into the lipid bilayer, then the inherent lipophilicity of the drug as penetrates deeper into the more hydrophobic centre of the bilayer, and (b) the reactivity of the drug in the oxidized regions of the lipid bilayer, which can be assessed by the HOMO-LUMO energy gap. It is step (a) that is slow and rate determining. Once the drug penetrates into the lipid bilayer, reaction with allylic radical PUFAs will be rapid.

Experimental Methods

All calculations were carried out using the Gaussian 09 package. Energy optimizations were at the DFT/B3LYP/6-31G(d,p) (6d, 7f) level of theory for all atoms in water. Selected optimizations at the DFT/B3LYP/6-311 + G(d,p) (6d, 7f) level of theory gave very similar results to those at the lower level. Optimized structures were checked to ensure energy minima were located, with no negative frequencies. Energy calculations were conducted at the DFT/B3LYP/6-31G(d,p) (6d, 7f) for neutral compounds using the IEFPCM/SMD solvent model. With the 6-31G* basis set, the SMD model achieves mean unsigned errors of 0.6 -1.0 kcal/mol in the solvation free energies of tested neutrals and mean unsigned errors of 4 kcal/mol on average for ions. [START_REF] Av Marenich | Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions[END_REF] The 6-31G** basis set has been used to calculate absolute free energies of solvation and compare these data with experimental results for more than 500 neutral and charged compounds. The calculated values were in good agreement with experimental results across a wide range of compounds. [START_REF] Rayne | Accuracy of computational solvation free energies for neutral and ionic compounds: Dependence on level of theory and solvent model[END_REF][START_REF] Rizzo | Estimation of Absolute Free Energies of Hydration Using Continuum Methods: Accuracy of Partial Charge Models and Optimization of Nonpolar Contributions[END_REF] Adding diffuse functions to the 6-31G* basis set (ie 6-31 + G**) had no significant effect on the solvation energies with a difference of less than 1% observed in solvents, which is within the literature error range for the IEFPCM/SMD solvent model. HOMO and LUMO calculations included both delocalized and localized orbitals (NBO).

Electron affinities (EA) in eV in water were calculated by the SCF difference between the optimised/relaxed neutral and optimised radical species method as previously described. [START_REF] Fong | Permeability of the Blood-Brain Barrier: Molecular Mechanism of Transport of Drugs and Physiologically Important Compounds[END_REF][START_REF] Fong | The extravascular penetration of tirapazamine into tumours: a predictive model of the transport and efficacy of hypoxia specific cytotoxic analogues and the potential use of cucurbiturils to facilitate delivery[END_REF][START_REF] Fong | Screening anti-colorectal cancer drugs: free radical chemotherapy[END_REF][START_REF] Fong | Free radical anticancer drugs and oxidative stress: ORAC and CellROXcolorectal cancer cells by quantum chemical determinations[END_REF][START_REF] Fong | The role of free radicals in the effectiveness of anti-cancer chemotherapy in hypoxic ovarian cells and tumours[END_REF][START_REF] Fong | Free radicals in chemotherapy induced cytotoxicity and oxidative stress in triple negative breast and ovarian cancers under hypoxic and normoxic conditions[END_REF][START_REF] Fong | Role of stable free radicals in conjugated antioxidant and cytotoxicity treatment of triple negative breast cancer[END_REF][START_REF] Fong | Toxicology of platinum anticancer drugs: oxidative stress and antioxidant effect of stable free radical Pt-nitroxides[END_REF][START_REF] Fong | The effect of desolvation on the binding of inhibitors to HIV-1 protease and cyclin-dependent kinases: Causes of resistance[END_REF] It has been shown that the B3LYP functional gives accurate electron affinities when tested against a large range of molecules, atoms, ions and radicals with an absolute maximum error of 0.2 eV. [START_REF] Jc Rienstra-Kiracofe | Atomic and molecular electron affinities: Photoelectron experiments and theoretical computations[END_REF] It is noted that high computational accuracy for each species in different environments is not the focus of this study, but comparative differences between various species is the aim of the study. The literature values for drug efficacy and toxicity used in the stepwise multiple regressions have much higher experimental uncertainties than the calculated molecular properties. The statistical analyses include the multiple correlation coefficient R 2 , the F test of significance, standards errors for the estimates (SEE) and each of the variables SE(ΔG desolCDS ), SE(ΔG lipoCDS ), SE(Dipole Moment), SE (EA) as calculated from "t" distribution statistics. Residual analysis was used to identify outliers. 1 and2 and the keto and enol tautomers referred to in Table 2 (5.15)

1.17

Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. Schematic representation of drug permeation into and through the phospholipid bilayer of a cell membrane showing dielectric gradient of bilayer

Table 1 was: IC 50 = 4.19ΔG desolv,CDS + 2.20ΔG lipo,CDS -2.39Dipole Moment -2.95(HOMO-LUMO) + 72.16

 1

	Where R 2 = 0.401, SEE = 7.73, SE(ΔG desolvCDS ) = 1.24, SE(ΔG lipoCDS ) = 1.53, SE(Dipole Moment) =0.92, SE(HOMO-LUMO) =
	2.50, F=5.20, Significance=0.0025

-34.74ΔG desolv,CDS + 22.79ΔG lipo,CD -47.81

  ΔG desolv,CDS, ΔG lipo,CDS, Dipole Moment, Molecular Volume, Adiabatic Electron Affinity or the HOMO-LUMO energy gap of the in vivo brain protection for cerebral ischemia in rats with molecular properties of Edaravone and Edaravone analogs in the keto tautomeric form resulted in eqs 4 and 5: Eq 4 in vivo brain protection for cerebral ischemia in rats, keto tautomers:Where R 2 = 0.520, SEE = 29.92, SE(ΔG desolvCDS ) = 16.65, SE(ΔG lipoCDS ) = 14.71, F=3.26, Significance=0.110 Eq 5 in vivo brain protection for cerebral ischemia in rats, keto tautomers:

	Survival Time =

EEG Recovery Time = -10.22ΔG desolv,CDS -10.44ΔG lipo,CDS -88.85

  Stepwise regression analysis tested against ΔG desolv,CDS, ΔG lipo,CDS, Dipole Moment, Molecular Volume or the HOMO-LUMO energy gap of the in vivo brain protection for cerebral ischemia in rats with molecular properties of Edaravone and Edaravone analogs in the enol tautomeric form resulted in eqs 6 and 7: Eq 6 in vivo brain protection for cerebral ischemia in rats, enol tautomers:

	Where R 2 = 0.604, SEE = 9.44, SE(ΔG desolvCDS ) = 5.25, SE(ΔG lipoCDS ) =4.64, F=4.58, Significance=0.061
	(ii) Survival Time = -17.70 HOMO-LUMO gap + 111.81
	Where R 2 = 0.562, SEE = 8.74, SE(HOMO-LUMO) = 5.52, F=10.27, Significance=0.012

Eq 7 in vivo brain protection for cerebral ischemia in rats, enol tautomers: EEG Recovery Time = -36.44ΔG desolv,CDS +40.05ΔG lipo,CDS + 116.00

Where R 2 = 0.611, SEE = 27.14, SE(ΔG desolvCDS ) = 12.60, SE(ΔG lipoCDS ) =18.62, F=5.50, Significance=0.036

Table 1 . Lipid peroxidation inhibition of Edaravone (Edv) and Edv analogs (keto tautomers) with molecular properties of Edv and Edv analogs (see Figure 3 for structures) No* IC 50 µM ΔG desolv,CDS kcal/mol ΔG lipo,CDS kcal/mol DM D Molec Vol cm 3 /mol AEA eV HOMO- LUMO

 1 .

	eV

Footnotes: No* refers to ref [START_REF] Watanabe | Structure-activity relationship of 3-methyl-1-phenyl-2-pyrazolin-5-one (edaravone)[END_REF] and R2, R4 and X refers to substituents in structures in Figure 3. Values in parentheses are for the enol tautomers, and values without parentheses are for the keto tautomers (see Figure 3)


Alexandre Paris
email: alexandre.paris@cea.fr

Philippe Heinrich

Raphael Paris

Stéphane Abadie

Raphaël Paris

Anak Krakatau, Indonesia, Landslide and Tsunami: Preliminary Modeling Results

Keywords: Tsunami, Landslide, Anak Krakatau, Simulation

published or not. The documents may come

Introduction

Landslide tsunamis are not very well known to the general public, although they are quite common and devastating. This phenomenon can happen wherever weakened volumes of rocks or sediments lie, i.e in deltas, coastal cliffs, rivers, fjords and lakes [START_REF] Masson | Submarine landslides: processes, triggers and hazard prediction[END_REF]. As the velocity of the terrain deformation due to a landslide can reach very high values (up to 100 m.s -1 according to [START_REF] Satake | Three-dimensional reconstruction and tsunami model of the Nuuanu and Wailau giant landslides, Hawaii[END_REF]), the tsunami energy may be initially of the same order as that generated by a major earthquake [START_REF] Okal | A theoretical comparison of tsunamis from dislocations and landslides[END_REF], but the generated amplitudes generally decrease quickly and the tsunami effects are local.

One of the largest landslide tsunamis may probably be the Storegga slide, off Norway, with its volume between 2400 and 3200 km 3 , 8200 years ago [START_REF] Haflidason | The Storegga slide: architecture, geometry and slide development[END_REF][START_REF] Bondevik | The Storegga slide tsunami -comparing field observations with numerical simulations[END_REF]. Another famous case is the 1958 Lituya Bay landslide-generated mega-tsunami [START_REF] Miller | Giant waves in Lituya Bay, Alaska[END_REF]. After a 8.3magnitude earthquake, a volume of 30.6 million m 3 collapsed in Gilbert Inlet, generating a tsunami which run-up height reached 524 m on the opposite side [START_REF] Fritz | Lituya Bay landslide impact generated mega-tsunami 50th anniversary[END_REF].

2018 Tsunami in Sunda Strait, Indonesia 3 Until now, one of the deadliest recent landslide tsunamis is the Papua New Guinea one, 1998, during which over 2100 people died after a 4 km 3 landslide collapsed at a depth of 550 m, generating a tsunami with run-up heights up to 15 m (Heinrich et al., 2001a;[START_REF] Synolakis | The slump origin of the 1998 Papua New Guinea tsunami[END_REF]. More recently, the 2014 Tangjiaxi, China (0.16 million m 3 , 3 deaths, 9 people missing and 11 injured [START_REF] Huang | Analysis of the Tangjiaxi landslide-generated waves in the Zhexi Reservoir, china, by a granular flow coupling model[END_REF]), the 2017 Karrat Fjord, Greenland (50 million m 3 , 4 people killed [START_REF] Paris | Numerical modeling of the June 17, 2017 landslide and tsunami events in Karrat Fjord, west Greenland[END_REF]) events and now the 2018 Anak Krakatau, Indonesia collapse remind us of the threat landslide tsunamis represent, even if their volume appear relatively small compared to historical events such as Ritter Island in 1888 (5 km 3 [START_REF] Cooke | Eruptive history of the volcano at Ritter Island[END_REF][START_REF] Johnson | Large-scale volcanic cone collapse; the 1888 slope failure of Ritter volcano[END_REF]) and Oshima-Oshima in 1741 (2.4 km 3 [START_REF] Satake | The 1741 Oshima-Oshima eruption: extent and volume of submarine debris avalanche[END_REF].

A possible way to better understand a physical phenomenon is to use numerical modeling. A review of landslide-generated tsunami models can be found in [START_REF] Heidarzadeh | The State-of-the-Art Numerical Tools for Modeling Landslide Tsunamis: A Short Review[END_REF] and Yavari-Ramshe and Ataie-Ashtiani (2016). Two approaches are used to simulate landslide-generated tsunamis: simulate the tsunami propagation considering the bottom deformation due to the landslide as a boundary condition for the water surface elevation or simulate both the landslide and the tsunami in a single model. In the first category, the landslide movement is reproduced as a bathymetry deformation in time (see [START_REF] Harbitz | Model simulations of tsunamis generated by the Storegga slides[END_REF], [START_REF] Grilli | Tsunami generation by submarine mass failure. I: Modeling, experimental validation, and sensitivity analyses[END_REF] or [START_REF] Tappin | The Papua New Guinea tsunami of 17 July 1998: anatomy of a catastrophic event[END_REF]). In the second category, the landslide is modeled together with the free surface using a rheology law. Among them we can find Newtonian fluid approaches as in the numerical experiments of [START_REF] Fine | Numerical Modeling of Tsunami Generation by Submarine and Subaerial Landslides[END_REF], or in the modeling of the 1979 Nice events [START_REF] Assier-Rzadkiewicz | Numerical modelling of a landslide-generated tsunami: The 1979 Nice event[END_REF], or in the Franz et al.'s (2015) study of the 2006 Nicolet, Québec, Canada, landslide. Non-Newtonian models such as the Bingham model of [START_REF] Skvortsov | Numerical simulation of the landslidegenerated tsunami in Kitimat Arm, British Columbia, Canada, 27[END_REF] or the BING model (e.g. the modeling of the 1888 Brattora, Norway, landslide tsunami [START_REF] L'heureux | The 1888 shoreline landslide and tsunami in Trondheimsfjorden, central Norway[END_REF] or the 2014 submarine landslide at Statland, Norway [START_REF] Glimsdal | The 29th January 2014 submarine landslide at Statland, Norway-landslide dynamics, tsunami generation, and run-up[END_REF]) can also be used. Finally, landslides rheologies can be modeled as granular flows, as in the study of Reunion Island landslide-tsunamis by [START_REF] Kelfoun | Landslide-generated tsunamis at Réunion Island[END_REF] or the Güímar debris avalanche simulation of [START_REF] Giachetti | Numerical modelling of the tsunami triggered by the Güìmar debris avalanche, Tenerife (Canary Islands): Comparison with field-based data[END_REF].

Landslides are mostly simulated using 2D depth-averaged (2DH) models, although more complex and time-consuming 3D computations may be envisaged [START_REF] Løvholt | Oceanic propagation of a potential tsunami from the La Palma Island[END_REF][START_REF] Abadie | Numerical modeling of tsunami waves generated by the flank collapse of the Cumbre Vieja Volcano (La Palma, Canary Islands): Tsunami source and near fiedl effects[END_REF][START_REF] Horrillo | A simplified 3-D Navier-Stokes numerical model for landslide-tsunami: Application to the Gulf of Mexico[END_REF]. Nevertheless, they are generally spatially restricted to the landslide and wave generation area. Tsunami propagation can be realized in 2DH using shallow water equations [START_REF] Jiang | The coupling of a submarine slide and the surface waves which it generates[END_REF][START_REF] Harbitz | Numerical simulations of large water waves due to landslides[END_REF] or Boussinesq equations [START_REF] Tappin | Did a submarine landslide contribute to the 2011 Tohoku tsunami?[END_REF][START_REF] Harbitz | Rockslide tsunamis in complex fjords: From an unstable rock slope at Åkerneset to tsunami risk in western Norway[END_REF][START_REF] Grilli | Modeling of SMF tsunami hazard along the upper US East Coast: detailed impact around Ocean City, MD[END_REF]. Here, we use the 2D depth-integrated model AVALANCHE [START_REF] Heinrich | Near-field modeling of the July 17, 1998 tsunami in Papua New Guinea[END_REF][START_REF] Hébert | Numerical modeling of the September 13, 1999 landslide and tsunami on Fatu Hiva Island (French Polynesia)[END_REF][START_REF] Friant | Numerical simulation of the last flank-collapse event of Montagne Pelée, Martinique, Lesser Antilles[END_REF] that has been successfully employed to simulate subaerial or submarine landslides [START_REF] Rodriguez | Owen Ridge deep-water submarine landslides: Implications for tsunami hazard along the Oman coast[END_REF][START_REF] Poupardin | The 1979 Submarine Landslide-Generated Tsunami in Mururoa, French Polynesia[END_REF][START_REF] Paris | Numerical modeling of the June 17, 2017 landslide and tsunami events in Karrat Fjord, west Greenland[END_REF], considering the landslide as a granular flow following a Coulomb frictional law and using shallow water or Boussinesq equations.

On December 22, 2018 at 13:50 UTC (20:50 local time) the southwestern flank of Anak Krakatau volcano (Sunda Strait, Indonesia) collapsed to the sea and generated a tsunami that devastated the coasts of Java and Sumatra, killing more than 430 people and damaging thousands of houses and boats, as reported by BNPB (Badan Nasional Penanggulangan Bencana1 ). Human casualties and material losses were recorded all around the Sunda Strait, on the islands of Java and Sumatra, with run-up heights reaching up to 14 m [START_REF] Tdmrc | The latest update from post-Sunda Strait tsunami survey[END_REF]. A summary of the available observations data is provided in Section 3. This Anak Krakatau event was recently studied by [START_REF] Grilli | Modelling of the tsunami from the December 22, 2018 lateral collapse of Anak Krakatau volcano in the Sunda Straits, Indonesia[END_REF], using the 3D model NHWAVE [START_REF] Ma | Shock-capturing non-hydrostatic model for fully dispersive surface wave processes[END_REF][START_REF] Ma | A two-layer granular landslide model for tsunami wave generation: Theory and computation[END_REF][START_REF] Kirby | The 27 April 1975 Kitimat, British Columbia, submarine landslide tsunami: a comparison of modeling approaches[END_REF] for the landslide simulation and the tsunami generation, then the 2D model FUNWAVE-TVD [START_REF] Shi | A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF] for the tsunami propagation. It was also simulated by [START_REF] Heidarzadeh | Numerical modeling of the subaerial landslide source of the 22 December 2018 Anak Krakatoa volcanic tsunami, Indonesia[END_REF] using the COMCOT model (Cornell Multi-grid Coupled Tsunami Model [START_REF] Liu | Computer Programs for Tsunami Propagation and Inundation[END_REF][START_REF] Wang | An analysis of 2004 Sumatra earthquake fault plane mechanisms and Indian Ocean tsunami[END_REF]) and an initial sea surface elevation as the landslide source. We present a comparison between water heights of the two latter mentioned studies and our simulation.

In this paper, the landslide-generated water waves are calculated by AVALANCHE and compared with available tsunami observations around Sunda Strait (tide gauge records, wave heights, and flow depths, inundation heights and run-up heights reported from three field surveys). The influence of the rheology on water wave heights is discussed using different basal friction angle values and finally the dispersive effects are assessed through numerical tests comparing a shallow water to a Boussinesq model.

Geological setting

The Krakatau volcanic complex is a group of islands that lies in the middle of the Sunda Strait, between Java and Sumatra. The tectonic setting of the strait is characterized by an extensional regime that formed different grabens at the southeastern end of the Sumatra fault zone [START_REF] Harjono | Seismicity of the Sunda strait: evidence for crustal extension and volcanological implications[END_REF][START_REF] Susilohadi | Structures and sedimentary deposition in the Sunda Strait, Indonesia[END_REF]. With a water depth not exceeding 200 m, the eastern part of the strait is relatively shallow, compared with the 1000 m deep Semangko graben to the west (see Figure 2 in [START_REF] Susilohadi | Structures and sedimentary deposition in the Sunda Strait, Indonesia[END_REF]).

Krakatau is the only active volcano of a south-southwest to north-northeast volcanic line that extends across the strait from Ujung Kulon in western Java to Rajabasa in eastern Sumatra [START_REF] Nishimura | Neotectonics of the Straits of Sunda, Indonesia[END_REF]. The 1883 calderaforming eruption completely reshaped the morphology of the Krakatau volcanic complex (e.g. [START_REF] Simkin | Krakatau 1883: The Volcanic Eruption and its Effects[END_REF]), and the present-day active edifice Anak Krakatau ('child of Krakatau' in Indonesian) was built on the steep northeast wall of the submarine caldera that was formed during the 1883 eruption [START_REF] Deplus | Inner structure of the Krakatoa volcanic complex (Indonesia) from gravity and bathymetry data[END_REF]. The 5×4 km 2 large rectangular caldera is characterized by a flat bottom of 200-240 m deep. The location of Anak Krakatau on the northeast rim of this steep-sided submarine basin led several authors to question its stability [START_REF] Camus | Petrologic evolution of Krakatau (Indonesia): implications for a future activity[END_REF][START_REF] Deplus | Inner structure of the Krakatoa volcanic complex (Indonesia) from gravity and bathymetry data[END_REF][START_REF] Giachetti | Tsunami hazard related to a flank collapse of Anak Krakatau volcano, Sunda Strait, Indonesia[END_REF]. During this 1883 eruption, a tsunami was generated, reaching 15 up to 40 m run-up heights in the Sunda Strait [START_REF] Nomanbhoy | Generation mechanism of tsunamis from the 1883 Krakatau eruption[END_REF][START_REF] Choi | Simulation of the transoceanic tsunami propagation due to the 1883 Krakatau volcanic eruption[END_REF] and killing more than 35000 people [START_REF] Sigurdsson | Pyroclastic flows of the 1883 Krakatau eruption[END_REF].

Anak Krakatau first emerged from the sea in January 1928 [START_REF] Stehn | The geology and volcanism of the Krakatau Group[END_REF].

From 1928 to 1959 phreatomagmatic activity progressively formed a 1.7 km large, 150 m high tuff-ring [START_REF] Sudradjat | The morphological development of Anak Krakatau Volcano, Sunda Strait[END_REF]. In 1960, eruptive activity shifted to a Strombolian style, thus building a new cone on the SW flank of the initial tuff-ring. During the last 60 years, recurrent lava flows and pyroclastic fall deposits contributed to the growth of Anak Krakatau Island, and the cone was 350 m high before its collapse in December 2018. Photographs captured from the coast of Java2 confirmed that the cone had grew by almost 100 m between July 2016 and December 2018. There are many coastal resorts, harbors and cities within a radius of 60 km around the volcano (e.g. Kalianda in Sumatra, Anjer, Labuhan and Sumur in Java (Figure 1)).

Anatomy of the event

Before any satellite revealed the effects of the on-going volcanic eruption and collapse, photographs captured by Øystein Lund Andersen in Anjer-Kidu (western coast of Java, 50 km east of Krakatau volcano) show a strong strombolian activity (Figure 2a), followed by a phase of ash emission that blurred the entire volcano (Figure 2b). The collapse was recorded on a seismic signal at 20:55 local time, i.e. 13:55 UTC [START_REF] Walter | Precursors and processes culminating in the Anak Krakatau December 2018 sector collapse and tsunami[END_REF]. The flank collapse which removed 93.8 million m 3 of subaerial volcanic rocks from the western side of the volcano (Gouhier and Paris, 2019) followed a period of 6 months of volcanic activity and rapid growth of the volcanic cone [START_REF] Gouhier | SO2 and tephra emissions during the December 22, 2018 Anak Krakatau flank-collapse eruption[END_REF][START_REF] Walter | Precursors and processes culminating in the Anak Krakatau December 2018 sector collapse and tsunami[END_REF].

Local effects of the tsunami generated by the flank collapse are visible on Rakata Island (Figures 3a, b andc) and Sertung Island (Figure 3d). Taking into account the trimline, the run-up height is estimated to 25 up to 30 m on Rakata and Sertung islands, located respectively 4 km south and 3 km north-northwest of Anak Krakatau Island. About 33 min after the collapse, the tsunami was recorded at Marina Jambu, Banten, east of the volcano on Java Island (Figure 1), with a first water wave height of 60 cm and 6 min later a 20 cm high wave was recorded at Ciwandan, north of Marina Jambu (Figure 4). At the same time, on Sumatra Island, a tide gauge located at Kota Agung, NW of the Anak Krakatau, recorded a wave height of 35 cm. Finally, 58 min after the collapse, the tsunami was recorded in the Bandar Lampung Bay at Panjang, north of the volcano, with a wave height of 37 cm. Note that for the two locations of Marina Jambu and Ciwandan, on Java Island, the first wave is not the highest one. In both cases, the second wave height is twice or more the first one and reached 140 cm at Marina Jambu, and 40 cm at Ciwandan. Thanks to the work of the Copernicus program3 , preliminary summary maps about the consequences of the tsunami were made available. These maps give an idea about inundation distances and damages on houses. On Java Island, the most affected areas are Labuhan, with up to 300 m of inundation distance and Tanjung Lesung (see Figure 1 Locations, inundation distances and recorded water heights are summarized in Figure 1 and Table 1.

Finally, the field surveys conducted by the KKP (Kementerian Kelautan dan Perikanan, Ministry of Marine Affairs and Fisheries) south of Tanjung Lesung and between Labuhan and Anjer-Kidu [START_REF] Muhari | The December 2018 Anak Krakatau Volcano Tsunami as Inferred from Post-Tsunami Field Surveys and Spectral Analysis[END_REF], the BMKG (Badan Meteorologi, Klimatologi, dan Geofisika) at Tanjung Lesung and Carita, and the Tsunami and Disaster Mitigation Research Center of Syiah Kuala University confirmed that the areas of Tanjung Lesung, Labuhan, Carita and Anjer-Kidu on Java Island and Kalianda on Sumatra Island were among the most impacted, with destroyed brick walls and boulder(s) displaced [START_REF] Tdmrc | The latest update from post-Sunda Strait tsunami survey[END_REF]. These different field surveys confirmed that the waves penetrated up to 300 m inland at Tanjung Lesung (run-up height of 5 m asl), 170 m at Carita (run-up height of 6.2 m asl) and 330 m at Sakarame, north Carita (run-up height of 4.6 m asl). A local run-up height of 13.5 m asl for an inundation distance of 125 m was measured in a location between Tanjung Lesung and Sumur. Another remarkable run-up height of 12.5 m was measured at Cipenyu Beach by [START_REF] Takabatake | Field survey and evacuation behaviour during the 2018 Sunda Strait tsunami[END_REF]. The locations of the places studied by these field surveys are listed in Tables 2 and3.

Methods

Landslide model

In the code AVALANCHE [START_REF] Heinrich | Numerical simulation of the December 1997 Debris Avalanche in Montserrat, Lesser Antilles[END_REF][START_REF] Paris | Numerical modeling of the June 17, 2017 landslide and tsunami events in Karrat Fjord, west Greenland[END_REF] the flank collapse is modeled using the one-phase grain-flow model of [START_REF] Savage | The motion of a finite mass of granular material down a rough incline[END_REF] and a Coulomb basal friction. The equations of conservation of mass and momentum are solved in a (x , y ) coordinate system linked to the topography (Figure 5):

∂h s ∂t + ∂ ∂x (h s u s ) + ∂ ∂y (h s v s ) = 0; (1) ∂ ∂t (h s u s ) + ∂ ∂x (h s u 2 s ) + ∂ ∂y (h s u s v s ) = - 1 2 κ ∂ ∂x (gh 2 s cos θ) + κgh s sin θ x + F x ; (2) ∂ ∂t (h s v s ) + ∂ ∂x (h s v s u s ) + ∂ ∂y (h s v 2 s ) = - 1 2 κ ∂ ∂y (gh 2 s cos θ) + κgh s sin θ y + F y (3)
where h s is the slide's thickness in a direction perpendicular to the slope, u = (u s , v s ) the depth-averaged velocity vector parallel to the slope, κ = 1 -ρ w /ρ s where ρ w and ρ s are the water and rock densities with a ratio ρ s /ρ w = 1.5 (for the subaerial part of the slide, κ is equal to 1), θ(x, y) the local steepest slope angle, θ x and θ y the slope angles along the x and y axes respectively, and F = -κgh s cos(θ) tan(φ)u/ u the friction forces, where φ is the friction angle of the sliding materials. Curvature terms representing the effects of coordinate transformations [START_REF] Savage | The dynamics of avalanches of granular materials from initiation to runout. Part I: Analysis[END_REF] are considered as second-order terms in this paper. In this set of equations, water acts on the slide only through a buoyancy term and any drag contribution is neglected.

The basal friction angle is adjusted through a sensitivity study to fit with the observed water waves at the four tide gauges that recorded the tsunami.

Values of 1, 2, 5 and 10°were tested and results were analyzed in the near-field (Gauge 1, Figure 6a) and in the far-field (Gauge 3, Figure 6b and Gauge 6, Figure 6c). Gauges locations are presented in Figure 9 and listed in Table 1.

Tsunami model

As in [START_REF] Paris | Numerical modeling of the June 17, 2017 landslide and tsunami events in Karrat Fjord, west Greenland[END_REF], tsunami generation is modeled by solving shallow water equations (Equations ( 4), ( 5) and ( 6)) during 80 seconds after the landslide triggering. Beyond 80 seconds, Boussinesq equations (Equations ( 7), ( 8) and ( 9)) are solved in order to take into account any possible dispersive effects during the propagation. The shallow water equations solved by AVALANCHE are written as:

∂η ∂t + ∂(hu) ∂x + ∂(hv) ∂y = - ∂d ∂t ; (4) ∂u ∂t + u ∂u ∂x + v ∂u ∂y = -g ∂η ∂x + F x ;
(5)

∂v ∂t + u ∂v ∂x + v ∂v ∂y = -g ∂η ∂y + F y ( 6
)
where η is the surface elevation, h=η+d the water column height, d is the depth, u and v the depth-averaged velocities along the x and y axes respec-tively and F x and F y the friction and Coriolis forces, which are assumed to be negligible in our modeling.

Following [START_REF] Løvholt | Oceanic propagation of a potential tsunami from the La Palma Island[END_REF], the Boussinesq equations read: where u t and v t the time-derivatives of u and v.

∂η ∂t + ∂(hu) ∂x + ∂(hv) ∂y = - ∂d ∂t ; (7) ∂u ∂t + u ∂u ∂x + v ∂u ∂y = -g ∂η ∂x + d 3 3 ∂ 2 u t ∂x 2 + ∂ 2 v t ∂x∂y + d 2 u t ∂ 2 d ∂x 2 + u t ∂ 2 d
Both the landslide and shallow water equations, which are very similar, are solved by Godunov's finite-volume scheme, extended to second order by a Van Leer scheme (Heinrich et al., 2001a;[START_REF] Labbé | Refined numerical modeling of the 1979 tsunami in Nice (French Riviera): Comparison with coastal data[END_REF]. This numerical scheme was validated with an analytical solution in [START_REF] Mangeney | Analytical Solution for Testing Debris Avalanche Numerical Models[END_REF].

Boussinesq equations are solved using a finite-difference scheme for spatial derivatives together with a Crank-Nicolson scheme for the temporal discretization. This latter scheme is based on an iterative procedure that uses centered differences for linear terms and forward differences for advection terms. The implicit momentum equations are solved by alternating implicit sweeps in the

x and y components using an ADI method (Alternating Direction Implicit).

For a given direction, the dispersion terms in the other direction are discretized explicitly. For each direction (x and y), a tridiagonal matrix is then solved at each iteration, following [START_REF] Pedersen | Documentation of a global Boussinesq solver[END_REF].

The sea-bottom deformation due to the landslide, ∂d/∂t in Equation ( 7), is computed as a forcing term: E, respectively, see Figure 9).

∂d ∂t = 1 cos θ ∂h s ∂t (

Results

Landslide simulation

The sensitivity study realized on the basal friction angle shows that its influence on the generated wave heights decreases with distance. The basal friction angle plays a major role in the near-field (Gauge 1, Figure 6a) but its influence is minor in the far-field (Gauges 3 and 6, Figure 6b andc). There is no significant difference between 1 and 2°, at least in terms of computed wave height at the coasts. Using a friction angle of 2°, calculated water heights are of 75 cm at Marina Jambu and 33 cm at Kota Agung, whereas observed ones are of 60 and 35 cm respectively. A friction angle of 10°results in wave heights of 56 and 16 cm respectively. Based on all the tide gauges, the best fit is obtained with a friction angle of 2°. This low value is consistent with the one used in [START_REF] Giachetti | Tsunami hazard related to a flank collapse of Anak Krakatau volcano, Sunda Strait, Indonesia[END_REF] and other studies about landslides on volcanoes slopes (e.g. [START_REF] Friant | Numerical simulation of the last flank-collapse event of Montagne Pelée, Martinique, Lesser Antilles[END_REF] with 7°for the flank collapse of Montagne Pelée (Martinique, Lesser Antilles), [START_REF] Kelfoun | Landslide-generated tsunamis at Réunion Island[END_REF] with values between 3 and 5°for different landslides scenarios envisaged at La Réunion Island or [START_REF] Giachetti | Numerical modelling of the tsunami triggered by the Güìmar debris avalanche, Tenerife (Canary Islands): Comparison with field-based data[END_REF] with values between 1.3 and 3.9°for reproducing the Güìmar debris avalanche (Tenerife, Canary Islands)).

A large part of the simulated landslide collapses in about 40 seconds (Fig- ure 8c), the volcano summit decreasing from about 350 m of altitude to about 120 m (Figures 7 and8), which is concordant with the topographic reconstruc-tion of [START_REF] Gouhier | SO2 and tephra emissions during the December 22, 2018 Anak Krakatau flank-collapse eruption[END_REF]. Beyond 40 seconds, the calculated ground displacement is negligible and no longer has any effect on the water surface deformation. The landslide covers a horizontal distance of about 4000 m from the volcano to the 1883 submarine caldera, with an average velocity of 35 m.s -1 .

Figure 7 shows both the numerical initial and final states of the volcano, with a cross-section illustrating the major topographic change induced by the flank collapse.

As shown by the temporal evolution of the landslide tip and the associated tsunami front (Figure 8), velocities are quite similar, which can signify efficient transfer process. But further investigation may be needed as the part of the landslide with the largest thickness situated in smaller depth may be more important for energy transfer than the slide tip, and the wave is a bit too fast to stay in contact with this area.

Tsunami simulation

Computed results are analyzed through a maximum surface elevations map (Figure 9) covering the entire domain. On the tsunami path, synthetic gauges are located at places where tide gauges recorded the tsunami, i.e. Marina Jambu and Ciwandan for Java, Kota Agung and Panjang for Sumatra (Gauges 3, 4, 6 and 7 respectively, see Figure 9 and Table 1) and results are compared with detided data (the tide filter is based on the MATLAB package T_TIDE [START_REF] Pawlowicz | Classical tidal harmonic analysis including error estimates in MATLAB using T_TIDE[END_REF]) in Figure 10. Other synthetic gauges are used, near the volcano, at Tanjung Lesung and Kalianda (Gauges 1, 2 and 5 respectively, see Figure 9 and Table 1) in order to compare the results to the field surveys data in these areas.

During the first 20 seconds of collapse the tsunami wave is being generated until it reaches a maximum positive elevation of 80 m (Figure 8). Sertung

Island is the first place to be hit after 80 seconds of wave propagation. It is located 3 km from the volcano which results in a wave celerity of about 38 m.s -1 . Rakata Island is hit in about the same time, with a speed of 50 m.s -1 .

Then the wave travels around the Sunda Strait and reaches Gauge 1 (5 min after the collapse, with 5.6 m of wave height (Figure 11a)), Gauge 2 (Tanjung Lesung, 29 min, 1.7 m (Figure 11b)), Gauge 3 (Marina Jambu, 34 min, 75 cm (Figure 10a)), Gauge 4 (Ciwandan, 42 min, 22 cm (Figure 10b)), Gauge 5 (Kalianda, 49 min, 1.2 m (Figure 11c)), Gauge 6 (Kota Agung, 45 min, 33 cm (Figure 10c)) and Gauge 7 (Panjang, 1:06 hours, 20 cm (Figure 10d)).

For the first wave, simulated water heights on the synthetic gauges are quite consistent with the tide gauges records (75 cm computed vs 60 cm recorded at Marina Jambu, 22 cm vs 20 cm at Ciwandan, 33 cm vs 35 cm at Kota Agung), excepted at Panjang (20 cm vs 37 cm). In addition, according to the concert video, the height of the first water wave is estimated to be 1 to 2 m, which is consistent with the computed water height of 1.7 m at Tanjung Lesung (Figure 11b, Gauge 2). However, computed arrival times are late compared to the recorded ones (1 min at Marina Jambu, 3 min at Ciwandan, 5 min at Kota

Agung and 8 min at Panjang). The general wave pattern is well reproduced at Marina Jambu but poorly at Ciwandan, Kota Agung and Panjang.

Dispersion assessment

In order to assess and highlight the potential dispersive effects in the Sunda Strait, shallow water simulation was performed with AVALANCHE and com-pared with the present Boussinesq simulation (Figure 12). Time series calculated by the two models are compared at three synthetic gauges (Gauges 1, 8 and 9, Figure 9) located in deep water at depths of 100 (Figure 12a), 1140 (Figure 12b) and 1960 m (Figure 12c). Close to the source, both simulations produce approximately the same first water wave as depicted in Figure 12a with water wave heights of about 5.5 m. However, dispersion effects are noticed in the far-field and in deep ocean. Compared to shallow water simulations, the period of the first wave increases during the propagation, ranging from about 1 min (Gauge 1, Figure 12a) to about 1 min 30 s (Gauge 8, Figure 12b) and about 2 min (Gauge 9, Figure 12c).

Discussion

This study presents preliminary results of the 22 December 2018 Anak Krakatau collapse and tsunami simulations, comparing the computed numerical results with several observed data such as flow depths or recorded surface elevations. The main goal of this study was to better understand the landslide that occurred and assess its volume, despite the lack of submarine information.

With the collapse of 150 million m 3 modeled following a Coulomb frictional law, the obtained water heights are quite consistent with the observed ones, all over the Sunda Strait, and the computed amplitudes of the first wave match the tide gauges-extracted amplitudes (with errors ranging from 6% at Kota

Agung to 25% at Marina Jambu), excepted at Panjang (error of 46%).

In comparison with the study of [START_REF] Grilli | Modelling of the tsunami from the December 22, 2018 lateral collapse of Anak Krakatau volcano in the Sunda Straits, Indonesia[END_REF], our results (water heights and time delays of the first wave) are very similar at the four tide gauges, with water heights differences of 50% at Marina Jambu, 12% at Ciwandan, 10% at Kota Agung and 0% at Panjang. In the near-field, our results Table 2 and Table 3 show computed water heights at places studied by the field surveys of Takabatake et al. ( 2019) and the BMKG, and the KKP respectively. Although the comparison between inundation or run-up heights with water wave heights may be hard to realise, we see for instance that at Cipenyu Beach, where highest run-up values were measured, we obtain the highest water heights. The computed wave heights are a little bit shorter than the observed ones by the BMKG survey, with 3 to 4 m vs 5 to 6 m. At the 5 locations close to the shoreline, i.e. distance lower than 30 m (Carita Lagon for BMKG survey (Table 2); Karang Suraga, Bulakan, Tanjung Jaya and Banyuasih for the KKP survey (Table 3)), the measured flow depths are on the same order that the computed tsunami height at the coast. Note that the bathymetric resolution (180 m interpolated to 25 m for the child grids) used in this study does not allow to compare precise inundation or run-up heights.

Although the results are encouraging, some limitations can be pointed out. First of all, the landslide reconstruction may be subject to discussion. It has been realized based on satellite images of the 23rd of December for the subaerial part and a quasi constant slope hypothesis for the submarine part.

Unfortunately, the intense volcanic activity of the days following the tsunami has completely changed the shape of the volcano and it has made precise reconstruction impossible. Nevertheless, tsunami simulation associated to this landslide is rather consistent with observed water heights. A sensitivity study on the basal friction coefficient φ (Figure 6) suggests that these water heights may vary quasi linearly with φ in the near-field, the smaller the friction angle, the higher the water wave heights. In the far-field along the coasts of Java or Sumatra, differences of water heights between φ=1°or φ=2°are not significant.

Second, the time delays of wave arrivals at the four tide gauges (Gauges 3, 4, 6 and 7, Figure 9) could indicate that the landslide volume and its dynamics are unknown. However, sensitivity tests (not shown in this study) on the landslide volume produce similar results in terms of wave arrivals. In addition, our travel times are similar to those calculated by [START_REF] Giachetti | Tsunami hazard related to a flank collapse of Anak Krakatau volcano, Sunda Strait, Indonesia[END_REF] and [START_REF] Grilli | Modelling of the tsunami from the December 22, 2018 lateral collapse of Anak Krakatau volcano in the Sunda Straits, Indonesia[END_REF] for landslide volumes of 280 million m 3 and 270 million m 3 respectively.

A possible explanation would be inaccuracies of the bathymetric data in this area. Travel times should thus be interpreted carefully. More information on the landslide will be available thanks to forthcoming surveys in the caldera.

Third, this study is a first attempt of simulating the collapse of A. K. volcano using 2D depth-averaged models both for the landslide and the tsunami.

Several complex phenomena are not taken into account such as mixing of the slide with the surrounding water, soil erosion or dissipation of water waves due to wave breaking or friction.

Finally, according to [START_REF] Glimsdal | Dispersion of tsunamis: doest it really matter?[END_REF], the dispersive effects can be estimated using the dispersion parameter:

τ = 4h 0 2 L λ 3 (11)
where h 0 is the depth at the source, L the distance to the coast of interest and λ the source width or in other terms the wavelength.

The dispersion parameter is ∼1.4 at Kalianda up to ∼3.4 at Kota Agung, with h 0 = 250 m, L between 45 and 110 km (see Table 1) and λ = 2000 m (see Figure 8b), and suggests that the propagation is highly dispersive. Results of the comparison between the shallow water and Boussinesq simulations (Figure 12) confirm that the first wave is subjected to dispersion, losing highfrequency components and being stretched by dispersive effects. After 25 min of propagation in deep ocean (Figure 12c), the period is approximately two times longer and the amplitude decreases by about 50% compared to the one calculated by the shallow water model. The use of the Boussinesq model in this study is therefore relevant.

Conclusion

The 22 December, 2018 Anak Krakatau southwestern flank collapse is modeled by a 2D depth-integrated code with a Coulomb frictional law and a basal friction angle of 2°. The generated tsunami is propagated using a Boussinesq model. This study highlighted some points and revealed some issues:

-The computed water heights from the present modeling fit well with those recorded at tide gauges of Marina Jambu, Ciwandan and Kota Agung (errors of 25%, 10% and 6% respectively), with the video recorded at Tanjung

Lesung and with the different observed data from the KKP, BMKG and Takabatake et al. ( 2019) field surveys.

-Results presented in this study are very similar to those obtained by [START_REF] Grilli | Modelling of the tsunami from the December 22, 2018 lateral collapse of Anak Krakatau volcano in the Sunda Straits, Indonesia[END_REF] both in the near-field and in the far-field. However, their model is 3D, suggesting proper representation of the governing processes in both studies. Heights and wavelengths of the first wave are also very close to those obtained by [START_REF] Heidarzadeh | Numerical modeling of the subaerial landslide source of the 22 December 2018 Anak Krakatoa volcanic tsunami, Indonesia[END_REF] in the near-field.

-Delays between simulated and observed travel times may be attributed to inaccuracies in the bathymetric data, particularly in the Bandar Lampung Bay.

-Despite the lack of submarine information at the bottom of the volcano, the collapse of the reconstructed volume of 150 million m 3 produces water waves that fit well with the observed results.

-The sensitivity study on the basal friction angle shows that its influence on the generated wave heights decreases with distance and that there are no significant differences of water heights between basal friction coefficients of 1 or 2°.

-Finally, a comparison between a Boussinesq model and a shallow water propagation highlights dispersive effects in the Sunda Strait that appear to be significant for the first waves. In deep waters, their periods progressively increase whereas their amplitudes progressively decrease. (green) and 9' (magenta) used by [START_REF] Grilli | Modelling of the tsunami from the December 22, 2018 lateral collapse of Anak Krakatau volcano in the Sunda Straits, Indonesia[END_REF] and presented in Figure 9b. Results are very similar to those presented in Figure 5e of [START_REF] Grilli | Modelling of the tsunami from the December 22, 2018 lateral collapse of Anak Krakatau volcano in the Sunda Straits, Indonesia[END_REF].

  for locations), where the water entered up to 430 m inland. South of Tanjung Lesung, a Youtube video 4 shows a water wave, that can be estimated between 1 and 2 m high, destroying a concert stage. Elsewhere, the inundation distances reached approximately 35 m north of Anjer-Kidu, 50 m at Teluk Lada, 170 m at Carita, 200 m south of Anjer-Kidu, and 250 m in Sumur. On Sumatra Island, the water penetrated the land up to 150 m at Lampung, 250 m at Kalianda and 280 m at Taman Agung.

  landslide and tsunami in the near-field (i.e. Anak Krakatau volcano and the other islands of the Krakatau archipelago), a pre-collapse topography and bathymetry of[START_REF] Gouhier | SO2 and tephra emissions during the December 22, 2018 Anak Krakatau flank-collapse eruption[END_REF] was used. The precollapse topography of Anak Island was derived from the DEMNAS (national digital elevation model of Indonesia, spatial resolution of 0.27 arc-second using the vertical datum EGM2008, provided by the Indonesian Geospatial Agency, and available at http://tides.big.go.id/DEMNAS/index.html). The original DEMNAS was slightly modified in order to include the latest growth of the edifice, as seen on photographs taken in August and November 2018 5 , and satellite images (e.g. Sentinel-2 image captured on 30 September 2018, and PlanetScope image captured on 17 December 2018). Pre-collapse bathymetry is from[START_REF] Deplus | Inner structure of the Krakatoa volcanic complex (Indonesia) from gravity and bathymetry data[END_REF]. As explained in[START_REF] Gouhier | SO2 and tephra emissions during the December 22, 2018 Anak Krakatau flank-collapse eruption[END_REF], the contour of the collapse scar was inferred from a Sentinel-1A image captured ∼8:30 hours after the collapse (22/12/2018 at 22:33:44 UTC) and photographs taken by Susi Air flight crew the day after (23/12/2018).Note that the subaerial volume estimated by[START_REF] Gouhier | SO2 and tephra emissions during the December 22, 2018 Anak Krakatau flank-collapse eruption[END_REF] likely corresponds to a minimum value because there is no data available on post-collapse bathymetry and the submarine extent of the collapse scar. Considering a slope between 5 and 8°in the continuity of the subaerial landslide part, the total volume approaches 150 million m 3 . This volume is two times smaller than the volume of 270 million m 3 modeled by[START_REF] Grilli | Modelling of the tsunami from the December 22, 2018 lateral collapse of Anak Krakatau volcano in the Sunda Straits, Indonesia[END_REF]. Pre and post-collapse topography and bathymetry around the Anak Krakatau are presented in Figures7c and drespectively.The bathymetric grids are built from the BATNAS one, with a spatial resolution of 180 m (available at http://tides.big.go.id/DEMNAS/Batnas.php).The model uses a system of multiple grids (coarse grid over deep water regions and fine grids over coastal regions) to model local effects of bathymetry. The parent nested grid covers the Sunda Strait, from eastern of Sumatra to western of Java Island (see Figures 1 or 9 for the grid footprint). Five child grids are used with a resolution of approximately 25 m: the first domain covers the Kota Agung Bay, the second the Bandar Lampung Bay, the third the Kalianda area from northern of Taman Agung to southern of Lampung, the fourth the Teluk Lada area from western of Tanjung Lesung to southern of Anjer-Kidu, and the fifth the Anjer-Kidu area from Marina Jambu to eastern of Ciwandan (grids A, B, C, D and

  are also very similar to the time series calculated at five additional gauges used by[START_REF] Grilli | Modelling of the tsunami from the December 22, 2018 lateral collapse of Anak Krakatau volcano in the Sunda Straits, Indonesia[END_REF] (see Figure9bfor gauges locations). The amplitudes of the generated water wave calculated by[START_REF] Grilli | Modelling of the tsunami from the December 22, 2018 lateral collapse of Anak Krakatau volcano in the Sunda Straits, Indonesia[END_REF] are obviously larger in the near-field since the authors consider a landslide volume of 270 million m 3 . Nevertheless, we obtain the same arrival times, periods and to a lesser extent the same wave behaviour (see Figure13in our paper and Figure5in[START_REF] Grilli | Modelling of the tsunami from the December 22, 2018 lateral collapse of Anak Krakatau volcano in the Sunda Straits, Indonesia[END_REF]). Our results are also similar in the near-field to those of[START_REF] Heidarzadeh | Numerical modeling of the subaerial landslide source of the 22 December 2018 Anak Krakatoa volcanic tsunami, Indonesia[END_REF]. Considering a forward modeling trial-error approach, the latter ones obtain a wave height of about 100 to 150 m with a wavelength of 1.5 to 2 km. Their results are close to our first simulated wave characterized by a maximum amplitude of 80 m and a wavelength of about 2 km (see Figure8in this study and Figure12in[START_REF] Heidarzadeh | Numerical modeling of the subaerial landslide source of the 22 December 2018 Anak Krakatoa volcanic tsunami, Indonesia[END_REF]).

Fig. 1

 1 Fig. 1 Location map of Anak Krakatau volcano in the Sunda Strait (brown volcano icon) and of the different observations of the event: inundation distance ID (less than 100 m in yellow, between 100 and 300 m in orange and more than 300 m in red), recorded height H at tide gauges (blue icons) and photographs (grey icons).

Fig. 2

 2 Fig. 2 Photographs taken in Anjer-Kidu, 50 km east of the volcano, a before the landslide at 19:05 local time and b after the landslide at 21:16 local time (https://www.oysteinlundandersen.com/krakatau-volcano-witnessing-the-eruptiontsunami-22december2018/).

Fig. 3 Fig. 4

 34 Fig. 3 Photographs of inundation effects in the near-field on a, b, c Rakata Island and dSertung Island (Figure1), close to the Anak Krakatau volcano (James Reynolds @EarthUn-cutTV on Twitter: https://twitter.com/earthuncuttc/status/1083305942228160513). The base of trees is estimated to be now 25 to 30 m above sea level.

Fig. 6

 6 Fig. 6 Time series calculated a at Gauge 1 near the Anak Krakatau volcano, b at Marina Jambu tide gauge (Gauge 3) and c at Kota Agung tide gauge (Gauge 6), using friction angles φ of 1°(black line), 2°(red line), 5°(green line) and 10°(blue line).

Fig. 8

 8 Fig. 8 Snapshots along a W-E cross-section of the simulated water wave (blue to red scale) and landslide (white to dark blue scale) at a t=0s, b t=20s, c t=40s and d t=60s.The post-collapse topo-bathymetry is represented in black solid line. There is a ×5 vertical exaggeration.

Fig. 10 Fig. 11

 1011 Fig. 10 Comparison between recorded data (in red) and simulated water surface elevations (in black) at a Marina Jambu, b Ciwandan, c Kota Agung and d Panjang (Gauges 3, 4, 6and 7, respectively, see Figure9).

Fig. 12 Fig. 13

 1213 Fig. 12 Comparison of time series calculated by the shallow water (in black) and the Boussinesq (in red) models at Gauges a 1, b 8 and c 9 (Figure 9). The water depths of these gauges are: a 100 m, b 1140 m and c 1960 m.

Table 1

 1 Summary

	Gauges	Coordinates	Depth	Travel time (min)		Wave Height (cm)
	location	(Longitude , Latitude)	(m)	Giachetti et al. (2012)	Comp.	Rec. Comp.	Rec.
	1-Near the volcano	105.31376 , -6.14967	100		5		560
	2-Tanjung Lesung	105.64894 , -6.47980	1.5	28	29		170	100-200
	3-Marina Jambu	105.84263 , -6.18953	1.8	33	34	33	75	60
	4-Ciwandan	105.95513 , -6.01441	1.1	42	42	39	22	20
	6-Kota Agung	104.62044 , -5.50120	1.0		45	39	33	35
	5-Kalianda	105.55701 , -5.70134	2.9	47	49		120
	7-Panjang	105.29036 , -5.44579	1.3	68	66	58	20	37

of results sorted by arrival times. Computed (Comp.) results are compared with recorded (Rec.) data and Giachetti et al.'s (2012) results. Gauges number refer to Figure 9.

https://bnpb.go.id/volume-tubuh-gunung-anak-krakatau-berkurang-jumlah-korbantsunami-bertambah, last accessed 10

september 2019

https://www.oysteinlundandersen.com/krakatau-volcano-witnessing-the-eruptiontsunami-22december2018/ 2018 Tsunami in Sunda Strait, Indonesia

https://emergency.copernicus.eu/mapping/list-of-components/EMSR335

https://www.youtube.com/watch?v=2ERXCR86GU4

https://www.oysteinlundandersen.com/krakatau-volcano/krakatau-eruption-seenfrom-anyer-west-java-17th-november-2018/

Acknowledgements We thank Abdul Muhari, from the Directorate of Sustainable Utilization of Coastal Zone and Small Islands, Ministry of Marine Affairs and Fisheries, Indonesia, for providing us with the tide gauges data and the KKP field survey data, and Gegar Prasetya, for providing us the BMKG field survey data.

Alexandre Paris et al.

Table 2 Computed water heights (W.H.) (m) at the shoreline compared with flow depths (F.D.), inundation heights (I.H.) and run-up heights (R.H.) identified by the [START_REF] Takabatake | Field survey and evacuation behaviour during the 2018 Sunda Strait tsunami[END_REF] and BMKG surveys. Values in parenthesis are the distances from the shoreline.

Values in brackets are the maximum inundation distances measured by the different team surveys.

Places


Enrique Garcia Moreno-Esteva
email: enrique.garciamoreno-esteva@helsinki.fi

Jessica F A Salminen-Saari

Miika Toivanen

Markku S Hannula

A picture is worth a thousand words: visualizing collaboration through gaze synchrony graphs

Keywords: Math problem-solving in class, collaboration, mobile gaze-tracking, gaze synchrony, analytic methods

Introduction and research questions

In our case study, in which we attempt to understand the reasons why collaborative behavior arises, we turned to eye tracking during a live lesson involving a mathematics problem solving session in a classroom, and measured shared attention via gaze synchrony. Gaze synchrony is a measure of how much the locus of gaze of the different participants in the study overlaps. Together with a qualitative analysis of the episode, we tried to determine what produces gaze synchrony and why. Our data collection took place in a classroom and not in a lab, making our results ecologically valid. Our resulting picture, an annotated graph, is worth a thousand words when supplemented with specific knowledge about the key events it identifies. With the graph we are able to locate moments when (gaze) focused communication and interaction takes place, and identify what brings about that communication and interaction.

Background

Collaboration is a process whereby two or more people work together towards achieving a particular goal. Intellectual work has become increasingly collaborative. There is remarkable consensus among educational policy makers, that the future labor market needs primarily nonroutine analytic and non-routine interpersonal skills, including problem solving, interpersonal communication, and self-management (OECD, 2013). Although the nature and learning of such 21st century skills (see Binkley et al. 2010) has been studied both in psychological and educational research there are some areas that have received little attention, if any. In our study we address the question of what brings about collaborative behavior during non-routine mathematics problem solving in class.

Methodology

In order to study patterns of collaboration during mathematics problem solving sessions in the classroom, student teachers at the University of Helsinki were asked to solve a non-routine problem during a lesson while wearing mobile eye-tracking devices. The problem consists of finding the optimal way to connect four cities located at the vertices of a square with cable, using the shortest amount of it (the optimal solution is given in Figure 1b). Several kinds of data were collected in addition to the eye-tracking data, such as notebook scribble, screen content when laptops or tablets were used, video, voice, and a stimulated recall interview with the subject wearing the eye-tracking devices during which they could see the eye-tracking recording as a form of stimulus. The students work first alone, then in pairs, and then in groups of four. The role of the teacher is to encourage students to find the best solution possible without providing information about it. Four students were wearing the gaze tracking glasses. The gaze tracking data was processed to measure gaze synchrony and produce the graph in Figure 1a. The red line represents a statistically generated mean of the measure of gaze synchrony. The blue line represents the amount of overlap as given by the gaze tracking data relative to this statistically derived mean, in standard deviation units. The green lines are three standard deviations away from the red mean curve.

Results and discussion

The graph above and the careful analysis of other audio and video data allow us to construct a storyline that describes the problem-solving episode and which highlights collaborative behavior.

Gaze synchrony occurs at the peaks. Leading to these moments is the need of one or more of the students to verify their belief that they have identified a plausible solution (something we call goal proximity). When this belief wanes, students return to work on their own in hope of finding a better solution and gaze synchrony is lost, as indicated by the troughs in the figure. For example, in the picture, towards the right, there is a long trough followed by a high plateau. During the trough, the students are working alone, but then, one of the students thinks he has found the best solution and asks everyone to look at the board to verify this, and the verification process ensues. The graph helped us identify this pattern of collaboration in which discovery and verification alternate, with discovery taking place while students work on their own, and verification signaled by increased gaze synchrony, occurring while the group works together and the visual attention of the group members is more focused on the same targets. We believe that with this tool we might be able to identify, with more data, other collaboration patterns occurring in problem solving sessions in the mathematics classroom and provide teachers with valuable insight informing her interventions.

Figure 1 :

 1 Figure 1: (a) 3-way gaze synchrony graph (b) the optimal solution to the problem