Viet Hoang

Mai Suan Li

Philippe Derreumaux

Junmei Wang

Toan T Nguyen
email: nguyen@ibpc.fr

S Nangia

Phuong H Nguyen

Molecular mechanism of ultrasound interaction with a blood brain barrier model

The brain is strictly protected by the blood brain barrier preventing the crossing of therapeutics to treat brain diseases. The high and low intensity focused ultrasound methods have been used to open temporarily the blood brain barrier, facilitating the transport of drugs. The methods are very promising because the opening is transient, localized and noninvasive. However, the molecular mechanism of the opening is unknown, and this limits the development and application of these methods. With this in mind, we carry out a molecular dynamics simulation study to understand the interaction of ultrasound with the cell membrane and the tight junction. Our minimal blood brain barrier model is composed of two lipid bilayers, mimicking two portions of neighbouring cells, connected together by a tight junction formed by a pair of two cis-dimers of the claudin-5 protein. Using an experimental ultrasound frequency of 50 MHz, simulations show that at low intensities, ultrasound does not impact the structure of the cell membranes and tight junction, implying that the direct interaction of ultrasound with the blood brain barrier is not responsible for the experimentally observed opening. At high intensities, the ultrasound pulls the monolayers of individual cell membrane lipid bilayers apart, creating air compartments inside the bilayers. This reduces the free energy barrier for the translocation of drugs across the lipid bilayer and enhances drug permeability. At very high intensities, the two monolayers are largely separated, resulting in cell damage and implying that the blood brain barrier is primarily opened at the experimentally observed damaged areas.

I. INTRODUCTION

The blood brain barrier (BBB) acts as the permeability barrier between the blood capillary and the brain [START_REF] Daneman | The blood-brain barrier[END_REF]. It is composed of endothelial cells connected together by junctional complexes consisting of tight junctions and adherens junctions. As a general rule, only small lipid soluble molecules with a molecular mass ≤ 500 Da may cross the BBB. Otherwise, different molecules may gain access to the brain only via certain endogenous transport systems within the BBB [START_REF] Janigro | Are you in or out? Leukocyte, ion, and neurotransmitter permeability across the epileptic blood-brain barrier[END_REF][START_REF] Vykhodtseva | Progress and problems in the application of focused ultrasound for blood-brain barrier disruption[END_REF]. As a consequence, only about 5% of ∼ 7000 known potential drugs for the brain treatment can cross the BBB [START_REF] Pardridge | Molecular Trojan horses for bloodbrain barrier drug delivery[END_REF].

As a remedy, several methods have been developed aimed at enhancing the permeability of drugs through the BBB [START_REF] Gross | Intrastriatal transplantation of microcarrier-bound human retinal pigment epithelial cells versus sham surgery in patients with advanced Parkinson's disease: a double-blind, randomised, controlled trial[END_REF][START_REF] Marks | Safety and tolerability of intraputaminal delivery of CERE-120 (adeno-associated virus serotype 2-neurturin) to patients with idiopathic Parkinson's disease: an openlabel, phase I trial[END_REF][START_REF] Pardridge | Reengineering biopharmaceuticals for targeted delivery across the blood-brain barrier[END_REF][START_REF] Matsukado | Intracarotid low dose bradykinin infusion selectively increases tumor permeability through activation of bradykinin B2 receptors in malignant gliomas[END_REF][START_REF] Rapoport | Functional brain imaging to identify affected subjects genetically at risk for Alzheimer's disease[END_REF]. Among these methods, the use of focused ultrasound is very promising as it is truly transient, localized and noninvasive for temporarily opening the BBB [START_REF] Hynynen | Noninvasive MR imaging-guided focal opening of the blood-brain barrier in rabbits[END_REF][START_REF] Mesiwala | High-intensity focused ultrasound selectively disrupts the blood-brain barrier in vivo[END_REF][START_REF] Baseri | Multi-Modality Safety Assessment of Blood-Brain Barrier Opening Using Focused Ultrasound and Definity Microbubbles: A Short-Term Study[END_REF][START_REF] Mcdannol | Temporary disruption of the blood-brain barrier by use of ultrasound and microbubbles: safety and efficacy evaluation in rhesus macaques[END_REF][START_REF] Horodyckid | Safe long-term repeated disruption of the blood-brain barrier using an implantable ultrasound device: a multiparametric study in a primate model[END_REF]. The very first study was carried out in 1955 by Barnard et al. to study the BBB permeabil-gether, after many decades, the induced-HIFU molecular mechanism of BBB opening is still unclear.

One of the reasons that prevents further development of the HIFU method is probably due to the fact that a safer method was latter developed by Hynynen and colleagues in 2001. In this method, a low intensity focused ultrasound (LIFU) is used in combination with microbubbles to open the BBB. Microbubbles concentrate ultrasound energy to desirable areas, thus low intensity ultrasound is sufficient to open the BBB, and cell lesions could therefore be minimised [START_REF] Hynynen | Noninvasive MR imaging-guided focal opening of the blood-brain barrier in rabbits[END_REF]. Since then, a large number of studies have been carried both in in vitro and in vivo, focusing on different aspects. One of the important aspects is the fundamental interaction mechanisms of LIFU with microbubbles, and effects of bubble cavitation on the BBB [START_REF] Sheikov | Cellular mechanisms of the blood-brain barrier opening induced by ultrasound in presence of microbubbles[END_REF]. Another aspect is the influence of ultrasound parameters such as pressure amplitude, frequency, burst length, pulse repetition frequency, sonication duration, and parameters of microbubbles such as dose, diameter and type, on the BBB opening mechanism [START_REF] Choi | Optimization of Blood-Brain Barrier Opening in Mice using Focused Ultrasound[END_REF][START_REF] Hynynen | Local and reversible blood-brain barrier disruption by noninvasive focused ultrasound at frequencies suitable for trans-skull sonications[END_REF][START_REF] O'reilly | Impact of standing wave effects on transcranial focused ultrasound disruption of the blood-brain barrier in a rat model[END_REF][START_REF] Sheikov | Effect of focused ultrasound applied with an ultrasound contrast agent on the tight junctional integrity of the brain microvascular endothelium[END_REF]. Also, many works have focused on the quantification of the permeability of the disrupted brain tissues [START_REF] Vlachos | Permeability assessment of the focused ultrasound-induced blood-brain barrier opening using dynamic contrast-enhanced MRI[END_REF][START_REF] Vlachos | Permeability dependence study of the focused ultrasound-induced bloodbrain barrier opening at distinct pressures and microbubble diameters using DCE-MRI[END_REF], test the delivery of new drugs, evaluate and monitor the response of treatments [START_REF] Chen | Novel magnetic/ultrasound focusing system enhances nanoparticle drug delivery for glioma treatment[END_REF][START_REF] Liu | Blood-brain barrier disruption with focused ultrasound enhances delivery of chemotherapeutic drugs for glioblastoma treatment[END_REF][START_REF] Jordao | Antibodies Targeted to the Brain with Image-Guided Focused Ultrasound Reduces Amyloid-β Plaque Load in the TgCRND8 Mouse Model of Alzheimer's Disease[END_REF][START_REF] Raymond | Ultrasound Enhanced Delivery of Molecular Imaging and Therapeutic Agents in Alzheimer's Disease Mouse Models[END_REF]. Nevertheless, the exact molecular mechanism of the LIFU induced BBB opening is still unknown.

Compared to experiments, theoretical studies of the ultrasound induced BBB opening are still at their infancy. Several mathematical models have been developed, which are essentially composed of coupled partial differential equations and solved numerically [START_REF] Nhan | Modeling localized delivery of Doxorubicin to the brain following focused ultrasound enhanced blood-brain barrier permeability[END_REF][START_REF] Wiedemair | On Ultrasound-Induced Microbubble Ossillation in a Capillary Blood Vessel and its Implications for the Blood-Brain Barrier[END_REF][START_REF] Hosseinkhah | Microbubbles and Blood-Brain Barrier Opening: A Numerical Study on Acoustic Emissions and Wall Stress Predictions[END_REF][START_REF] Hosseinkhah | A Three-Dimensional Model of an Ultrasound Contrast Agent Gas Bubble and its Mechanical Effects on Microvessels[END_REF][START_REF] Hosseinkhah | Mechanisms of Microbubble-Vessel Interactions and Induced Stresses: A Numerical Study[END_REF][START_REF] Qin | Acoustic response of compliable microvessels containing ultrasound contrast agents[END_REF][START_REF] Wiedemair | The breakup of intravascular microbubbles and its impact on the endothelium[END_REF]. These studies provide valuable information on the effect of bubbles cavitation and ultrasound on the vessel wall. However, relying on continuum mechanics, these models cannot capture detailed conformational changes of the BBB at the molecular level. At the molecular level, a number of simulations employing coarse-grained and atomistic molecular models have been carried out to study the effects of the bubble inertial cavitation and shock waves on the cell membranes. These studies have confirmed that inertial cavitation can induce the formation of membrane pores [START_REF] Ganzenmuller | ShockWave Induced Damage in Lipid Bilayers: A Dissipative Particle Dynamics Simulation Study[END_REF][START_REF] Choubey | Poration of Lipid Bilayers by Shock-Induced Nanobubble Collapse[END_REF][START_REF] Schanz | Molecular Dynamics Simulations of Cavitation Bubble Collapse and Sonoluminescence[END_REF][START_REF] Nomura | Nanobubble Collapse on a Silica Surface in Water: Billion-Atom Reactive Molecular Dynamics Simulations[END_REF][START_REF] Koshiyama | Structural Change in Lipid Bilayers and Water Penetration Induced by Shock Waves: Molecular Dynamics Simulations[END_REF][START_REF] Santo | Shock Wave Induced Collapse of Arrays of Nanobubbles Located Next to a Lipid Membrane: Coarse Grained Computer Simulations[END_REF][START_REF] Santo | with a Phospholipid Membrane: Coarse-Grained Computer Simulations[END_REF][START_REF] Fu | Sonoporation at Small and Large Length Scales: Effect of Cavitation Bubble Collapse on Membranes[END_REF][START_REF] Adhikari | Mechanism of membrane poration by shock wave induced nanobubble collapse: a molecular study[END_REF][START_REF] Agrawal | 145 A pressure-transferable coarse-grained potential for modeling the shock Hugoniot of polyethylene[END_REF]. Recently, a molecular dynamics (MD) simulation study based on a coarse-grained model explored the possibility of opening a simple tight junction due to shock wave induced bubble collapse [START_REF] Goliaei | Opening of the Blood-Brain Barrier Tight Junction Due to Shock Wave Induced Bubble Collapse: A Molecular Dynamics Simulation StudyJ[END_REF].

The molecular mechanisms of the BBB opening under HIFU and LIFU irradiations are very complex, because they results from interaction of ultrasound and bubble interaction with the cellular membranes and the tight junctions. To simulate these interactions, in principle we must have a bubble model, an ultrasound simulation method and a molecular BBB model. With this in mind, we have recently developed a bubble model and a focused ultrasound simulation method [START_REF] Viet | Nonequilibrium All-Atom Molecular Dynamics Simulation of the Ultrasound Induced Bubble Cavitation and Application to Dissociate Amyloid Fibril[END_REF][START_REF] Viet | Rayleigh-Plesset Equation of the Bubble Stable Cavitation in Water: A Nonequilibrium All-Atom Molecular Dynamics Simulation Study[END_REF][START_REF] Viet | Molecular mechanism of the cell membrane pore formation induced by bubble stable cavitation[END_REF]. This allows us to study the molecular mechanism of the interaction between ultrasound with lipid cell membrane, in the presence of bubble [START_REF] Viet | Molecular mechanism of the cell membrane pore formation induced by bubble stable cavitation[END_REF] or without bubble [START_REF] Viet | Interaction mechanism between the focused ultrasound and lipid membrane at the molecular level[END_REF]. For the molecular BBB model, Nangia and coworkers have recently obtained the structure of the claudin-5 protein via homology modelling of the experimental structure of the claudin-15 protein. Their MD simulations have shown that claudin-5 monomers embedded in seven membranes with different lipid compositions readily form stable cisdimers, which subsequently form contiguous strands and high-order assemblies [START_REF] Irudayanathan | Molecular Architecture of the Blood Brain Barrier Tight Junction Proteins-A Synergistic Computational and In Vitro Approach[END_REF]. These cis-dimers are used to construct the trans interfaces, which are formed when claudin-5 proteins on adjacent cells interact head-on to form the tight junction assembly, by using docking methods in combination with symmetric refinement [START_REF] Irudayanathan | Architecture of the paracellular channels formed by claudins of the blood-brain barrier tight junctions[END_REF]. Fundamental pore features, such as diameter and length of the pore were characterized in details [START_REF] Rajagopal | Computational Nanoscopy of Tight Junctions at the Blood-Brain Barrier Interface[END_REF].

As a first attempt to provide more insights into the ultrasound induced BBB opening molecular mechanism due to HIFU or LIFU, in this work, based on the above mentioned achievements, we carry out MD simulations aimed at understanding the impact of ultrasound alone on the cell membranes and tight junction using the minimal BBB model. In particular, we wish to understand which part of the BBB model: membrane or tight junction resists more to the ultrasound. The impacts of the bubble cavitation, and the construction of more sophisticated BBB models will be carried out in next studies.

II. METHODS

A. The minimal BBB model

The family of claudin proteins forms tight junctions in endothelial and epithelial cells [START_REF] Piontek | Formation of tight junction: determinants of homophilic interaction between classic claudins[END_REF]. Among 27 known members of the claudin family, the claudin-5 is the key tight junction protein with an expression level that is much higher than other claudin proteins [START_REF] Nitta | Size-selective loosening of the blood-brain barrier in claudin-5-deficient mice[END_REF][START_REF] Rossa | Overexpression of claudin-5 but not claudin-3 induces formation of trans-interaction-dependent multilamellar bodies[END_REF][START_REF] Cording | In tight junctions, claudins regulate the interactions between occludin, tricellulin and marvelD3, which, inversely, modulate claudin oligomerization[END_REF]. Thus, we use only claudin-5 protein to construct the tight junction. Our minimal BBB model consists of two lipid bilayer membranes connected together by a paracellular pore formed by two claudin-5 cis-dimers. The molecular structure is shown in Fig. 1 and briefly described here. First, the structure of the claudin-5 monomer is obtained via homology modelling of the experimental structure of the claudin-15 protein (PDB ID: 4P79). The claudin-5 protein consists of 218 residues, forming four transmembrane helix bundles (TM1-4), two extracellular loops (ECL1, ECL2) and an intracellular loop (ICL). The structure of a claudin-5 cis-dimer in the phosphatidylcholine (POPC) membrane was obtained from previous MD simulations (model B) [START_REF] Irudayanathan | Molecular Architecture of the Blood Brain Barrier Tight Junction Proteins-A Synergistic Computational and In Vitro Approach[END_REF]. In this structure, two claudin-5 monomers are mediated by the interaction between leucine residues 83, 90 of TM2 of one monomer and residues 124,131 of TM3 of another monomer, and flanking aromatic π-π interactions of Phe127 and Trp138 of two TM3. Then, a pair of two cis-dimers is used to construct the paracellular pore by using the docking methods in combination with symmetric refinement. This trans pair is stablized by hydrogen-bonds and salt bridges between ECL contacts. This minimal BBB model captures This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. the fundamental interactions responsible for the BBB tight junction assembly [START_REF] Irudayanathan | Architecture of the paracellular channels formed by claudins of the blood-brain barrier tight junctions[END_REF][START_REF] Rajagopal | Computational Nanoscopy of Tight Junctions at the Blood-Brain Barrier Interface[END_REF], and serves as a starting point for simulating effect of ultrasound on the BBB. We note that our BBB model is quite similar to that from Ref. [START_REF] Goliaei | Opening of the Blood-Brain Barrier Tight Junction Due to Shock Wave Induced Bubble Collapse: A Molecular Dynamics Simulation StudyJ[END_REF], which is composed of a pair of two claudin-15 dimers at the TJ We use the all-atom CHARMM36 force field [START_REF] Mackerell | All-atom empirical potential for molecular modeling and dynamics studies of proteins[END_REF] to model the protein and POPC lipid, and the TIP3P water model [START_REF] Jorgensen | Comparison of simple potential functions for simulating liquid water[END_REF] to describe solvent. The whole model is solvated in a water box consisting of ∼ 180 000 waters. The initial dimensions of the unit cell are (L x , L y , L z ) = [START_REF] Hynynen | Noninvasive MR imaging-guided focal opening of the blood-brain barrier in rabbits[END_REF][START_REF] Hynynen | Noninvasive MR imaging-guided focal opening of the blood-brain barrier in rabbits[END_REF][START_REF] Ballantine | Progress and problems in the neurological applications of focused ultrasound[END_REF] nm. Starting from this structure, an equilibrium MD simulation is carried out for 100 ns in the NPT ensemble with the pressure P 0 = 1 bar and temperature T = 300 K, employing the GROMACS simulation package [START_REF] Lindahl | GROMACS 3.0: A Package for Molecular Simulation and Trajectory Analysis[END_REF]. The last structure is used as initial structure for ultrasound simulations.

PLEASE CITE THIS ARTICLE AS

B. The ultrasound simulation method

In a conventional simulation, the pressure of the system is maintained at a desired value P 0 by using a barostat. Currently, the Berendsen barostat is widely used as it is simple in both theory and implementation. Basically, it consists of an external bath so that the global pressure of the system is weakly coupled to this pressure bath by rescaling periodically the lengths of the system [START_REF] Berendsen | Molecular-Dynamics With Coupling to an External Bath[END_REF]

r i → µr i (i = 1 • • • N ) and L → µL, (1)
where N , r i and L are, respectively, the number of atoms in the system, the coordinate of the i-th atom and the length of the system box. The scale factor is given by

µ = 1 - β∆t τ p (P (t) -P 0 ) 1/3 , (2)
where P (t), β, ∆t, and τ p are the instantaneous pressure in the system, the isothermal compressibility, the integration time step and the temperature coupling constant, respectively. After scaling lengths through Eq. 1, the system volume becomes V → [1 -β∆t τp (P (t) -P 0 )]V . This ways, the instantaneous pressure P (t) always fluctuates around the desired value P 0 .

In an ultrasound simulation, the instantaneous pressure of the system will be increased and decreased around the reference value P 0 following the compression and rarefaction of the sound wave. We take into account this effect of the ultrasound wave, which has the form

p(t) = A sin(2πωt), (3)
by using a modified scale factor µ = 1 -β∆t τ p (P (t) -P 0 -A sin(2πωt))

1/3 , (4)
where ω and A are the frequency and amplitude of the ultrasound. The lengths of the system are then scaled as usual (Eq.1). This guarantees that the pressure of the system is always equal to the pressure of the ultrasound.

The GROMACS simulation package [START_REF] Lindahl | GROMACS 3.0: A Package for Molecular Simulation and Trajectory Analysis[END_REF] coupled to our code of ultrasound is used for all the simulations. In all simulations, the ultrasound with frequency ω = 20 MHz is used, and the amplitude is varied from 100 to 500 bar. The reference pressure, P 0 = 1 bar, the pressure coupling constant, τ p = 1 ps and the isothermal compressibility β = 4.5 × 10 -5 bar -1 are used. To ensure that the damage to the membrane is not due to heat generated by work done by ultrasound, we couple both membrane and water to the heat bath at 300 K employing the Berendsen coupling method [START_REF] Berendsen | Molecular-Dynamics With Coupling to an External Bath[END_REF] with a temperature coupling constant of 0.1 ps. The equations of motion are integrated using the leapfrog algorithm with a small time step of 2 fs. The electrostatic interactions are calculated using the particle mesh Ewald method and a cutoff of 1.4 nm [START_REF] Darden | Particle Mesh Ewald: An N•log(N) Method for Ewald Sums in Large Systems[END_REF]. A cutoff of 1.4 nm is used for the van der Waals interactions. The nonbonded pair lists are updated every 5 fs. The data is saved for every 25 ps for subsequent analyses. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS

III. RESULTS

A. Response of the system to the ultrasound

In all simulations, the ultrasound period τ = 20 ns (50 MHz) is used, and the amplitude is varied in the range A = 100 -500 bar. To obtain a first impression, Fig. 2(A) shows, as an example, the time evolution of the ultrasound p(t) with amplitude A = 100 bar. The negative and positive periods are called the rarefaction and compression phases, respectively. For each ultrasound amplitude, the simulation is run for 120 ns, i.e. the system is excited by five ultrasound pulses, and the response of the system to the ultrasound is shown in Fig. 2. Let us first present the results for the case A = 400 bar. During the first rarefaction phase, i.e. p(t) ≤ 0 (t ≤ τ /2), the system is expanded with an increase in the volume from the initial value V ∼ 1650 nm 3 to V ∼ 1700 nm 3 at 5 ns [Fig. 2(B)]. As a consequence, the pressure inside the system becomes negative, and decreases from the initial equilibrium value of 1 bar to -400 bar [Fig. 2(C)]. As a result, atoms are pulled apart, thus the attractive interactions, including the long-ranged Coulombic and van der Waals interactions become dominant. Thus the potential energy of the system is decreased from -1.30 ×10 6 to -1.29 ×10 6 kJ/mol [Fig. 2(D)]. Because the system is coupled to the heat bath, its temperature is always fluctuating around the reference value of 300 K [Fig. 2(E)]. Within the time interval τ /4 ≤ t ≤ τ /2, the ultrasound intensity is reduced, but still in the rarefaction phase, thus the system is relaxed, as indicated by the decrease in the volume, pressure and potential energy, and achieved the equilibrium state at t = τ /2 = 10 ns. Next, the ultrasound is in the compression phase, i.e. p(t) ≥ 0 (τ /2 ≤ t ≤ τ ), the system is compressed as shown by the decrease and increase in the volume and pressure, respectively. The atoms are pushed closely, thus the short-ranged Coulombic and van der Waals interactions are dominated, therefore the potential energy of the system is increased from -1.30 ×10 6 to -1.31 ×10 6 kJ/mol [Fig. 2(D)]. After t ≥ 3τ /4, the ultrasound intensity is reduced, but still in the compression phase, the system volume is increased, the pressure is reduced and the potential energy is increased. At t = τ , the ultrasound is vanished and the system reaches the equilibrium state, finishing one ultrasound irradiation period. The response of the system is repeated for the next four ultrasound periods. As seen from Fig. 2, the response of the system to the weaker ultrasounds with A = 100 or 200 bar, is similar to that of A = 400 bar, but less intense. However, by increasing the ultrasound intensity to A = 450 bar, the system is able to oscillate with the ultrasound within the first period (t ≤ 20 ns), but then exploded during the second ultrasound period. That is, at t = 25 ns, when the ultrasound is fully expanded, the system volume suddenly increases to a large value, the pressure is reduced immediately to zero, and the potential energy becomes very large because atoms are largely separated apart. Although the system is coupled to the heat bath, as the expansion takes place very fast, the velocity of atoms are not rapidly rescaled and the temperature of the system is significantly increased. Technically, the simulation is crashed.

B. Response of the claudin-5 tight junction to the ultrasound

To understand the effect of ultrasound on the claudin-5 tight junction, we calculate the root-mean-square displacement (RMSD) of two cis-dimers (each anchored to each membrane) with respect to the initial structure for different ultrasound intensities. As seen from Fig. 3(A), the RMSDs undergo small increases of ∼ 0.25 nm compared to the initial structure, despite large variations in the ultrasound intensities from 100 to 450 bar. Interestingly, although the system is exploded with A = 450 bar as shown above, but the RMSD does not change much, indicating that the overall structure of the tight junction is well-maintained. Fig. 3(B) shows the total number of the intermolecular hydrogen bonds (Hbond) between two dimers. We observe that the number of Hbond fluctuates in time, but it does not follow the oscillation of the ultrasound pressure. After 120 ns, two dimers are slightly separated, as indicated by a small reduction in the number of Hbond from 5 to 3 for A = 400 bar. Even at 450 bar where the system is exploded, there is no large change in the number of Hbond, indicating that the tight junction is not disrupted.

Next, we wish to understand the impact of the ultra-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. sound on the secondary structures of the claudin-5 proteins. To this end, we calculate the population of the secondary structures of two dimers by using the STRIDE program [START_REF] Frishman | Knowledge-Based Protein Secondary Structure Assignment Proteins[END_REF] and the results are shown in Fig. 4. As seen, the extracellular loops of claudin-5 proteins, which are mainly in the β structure, contribute initially to ≈ 15% of the secondary structure population. The transmembrane helix bundles contribute to ≈ 65%, and ≈ 20% are attributed to the turn and coil structures of the linkers and of the intracellular loops. We observe that the sec- ondary structures hardly undergo any changes, irrespective large variations in the ultrasound intensities, even at 450 bar where the system is exploded. All these results show that both 3D and 2D structures of individual claudin-5 proteins, and of the whole tight junction are largely maintained, and the tight junction is not disrupted by the ultrasound irradiation.

PLEASE CITE THIS ARTICLE AS

C. Response of the lipid bilayer membranes to ultrasound

To investigate the response of the lipid bilayers to the ultrasound, we first calculate the carbon-hydrogen order parameters of the lipid tails S CH = 3 cos 2 θ-1 /2, where θ is the angle between a C-H bond vector and the bilayer normal. The angular brackets represent molecular and temporal ensemble averages. For a POPC lipid, the number of C-H vectors in the first, s n-1 , and second, s n-2 , tails are 15 and 17, respectively. To obtain an overall picture on the order of each tail, we simply take the average of the order parameters of all C-H vectors pertaining to that tail. The results are shown in Figs.5 for two tails. As seen, for the ultrasound intensities A = 100, 200 and 400 bar, the order parameters of both tails do not undergo any oscillation with ultrasounds, indicating that lipids are quite ordered. However, with A = 450 bar, the order parameters drop to zeros at t = 25 ns, a moment when the system is exploded, indicating significantly structural changes in the membranes. To understand this, we calculate the shortest intra-distance between two monolayers of each bilayer, and the shortest inter-distance between two bilayer membranes. The latter measures the gap This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. of the tight junction between two cells. The results are shown in Fig. 6. As seen, two monolayers of individual lipid bilayer membranes are always separated by intradistances of ≈ 0.18 nm [Figs.6(A), (B)], and two lipid bilayer membranes of two adjunct cells are stay intact at ≈ 0.9 nm for ultrasounds with A = 100, 200 and 400 bar. A visual inspection of a snapshot at 120 ns of the trajectory with A = 400 bar shown in Fig. 7(A) confirms these results, that is the membranes and the claudin-5 tight junction structures are well-maintained, and similar to the initial structures shown in Fig. 1(B). However, with A = 450 bar, we observe that the intra-distances between two monolayers of individual bilayers increase and become suddenly large at t = 25 ns when the system is exploded. At the same time, the inter-distance between two bilayer membranes is reduced from 0.9 to 0.5 nm. To understand this, we show in Fig. 7 two snapshots just before collapsing at t = 24.5 and 24.7 ns, and one snapshot at the collapsed moment at 25 ns. We observe that when the system is about to explode, the two monolayers of individual membranes are separated apart, and this pushes two lipid bilayer membranes of two cells to each others [Figs. 7(B),(C)] . When the system is fully exploded, two monolayers are highly separated, two cells get closer, and lipids become disorder. This explains the reduction in the order parameters shown in Fig. 5. Nevertheless, the claudin-5 tight junction is hardly affected, and still maintained by two monolayers of two adjunct cells [Fig. 7(D)].

PLEASE CITE THIS ARTICLE AS

In the HIFU experiments the direct interaction of ultrasound with the BBB is the cause of the BBB disruption. For the LIFU experiments, it is widely believed that ultrasound induces the stable and/or inertial cavitation of bubble, which in turn exerts acoustic shear forces and/or shockwave on the BBB, leading to the BBB opening. However, in principle, ultrasound also interacts directly with the cells and tight junction. This results in complicated BBB opening mechanisms, especially at the molecular level. Thus, our strategy is to split this complex problem into simpler tasks. In this work, our primary aim is to understand how the lipid bilayer cell membrane and the tight junction respond to the ultrasound and whether the cell membrane and/or the claudin-5 junction are affected by ultrasound. Before discussing the results, the limitations of our study should be pointed out. First, our BBB model is very simple, composed of only two lipid bilayers, mimicking two portions of neighbouring cells. The tight junction is composed of a paracellular pore formed by two cis-dimers of claudin-5 proteins. Nevertheless, given the fact that currently there is no experimental molecular structure of the BBB, this computational model still can serve as a starting point for MD simulations, which can provide some insights into the BBB opening mechanism. Second, the experimental frequency of the ultrasound is usually tuned between ω = 0.5-50 MHz, the ultrasound pressure is in the range 10 -50 bar for LIFU experiments [START_REF] Aryal | Ultrasound-mediated blood-brain barrier disruption for targeted drug delivery in the central nervous system[END_REF][START_REF] Cheng | Theranostic Strategy of Focused Ultrasound Induced Blood-Brain Barrier Opening for CNS Disease Treatment[END_REF], and 100 -1000 bar for HIFU experiments [START_REF] Canney | Acoustic characterization of high intensity focused ultrasound fields: A combined measurement and modeling approach[END_REF]. The ultrasound irradiation duration is usually in the second to minute timescales [START_REF] Aryal | Ultrasound-mediated blood-brain barrier disruption for targeted drug delivery in the central nervous system[END_REF][START_REF] Cheng | Theranostic Strategy of Focused Ultrasound Induced Blood-Brain Barrier Opening for CNS Disease Treatment[END_REF]. To be close to experiments, we use an ultrasound frequency of 50 MHz for our simulations. To be computationally feasible with current computer technology, the simulation time can only be in the nanosecond scale. This means that our system is only irradiated by several ultrasound cycles. We use the typical HIFU intensity in the range 100-500 bar for our simulations, which are strong enough to allow us to obtain response of the BBB within the simulation timescales. For low intensities used in LIFU experiments, the timescale of simulations would be much longer to obtain any significant responses of the BBB, and this is beyond our current capability.

Our simulations show that with weak ultrasound intensities A ≤ 400 bar, the structures of both lipid cell membranes and claudin-5 tight junction are well-maintained. This implies that for the HIFU experiments, the direct interaction of ultrasound with BBB does not yield to the opening neither at the cells nor tight junction. In the context of the LIFU experiments, where the ultrasound intensity is even weaker, this interaction can be safely excluded in the interpretation of the BBB opening mechanism, confirming that the interaction of the bubble cavitation with BBB is the main cause as hypothesised. With stronger ultrasound intensities, say A = 450 bar, This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. we observe that two monolayers within individual lipid bilayers are detached following the expansion of the simulation box during the ultrasound rarefaction phase, and this could result in the damage of the cell membrane. In contrast, the tight junction between two cells is hardly affected. This suggests that the BBB disruption likely takes place at the cell surface instead of at the cell junction.

PLEASE CITE THIS ARTICLE AS

To explain this, we calculate the interaction potential energy acting on two monolayers M1 and M2 of a bilayer shown in Fig. 8. As seen, the monolayer M1 is stabilised by potential energy V(M1) which is contributed by four interactions V(M1) = V(M1-C1) + V(M1-W1) + V(M1-C2) + V(M1-M2), and the monolayer M2 is stabilised by three interactions V(M2)=V(M2-C1)+ V(M2-W1) + V(M2-M1), where V(M1-C1) is the interaction potential energy between the monolayer M1 and the claudin cisdimers C1 embed in M1, V(M1-W1) is the interaction energy between M1 and waters between two cells, and so on. Fig. 8 shows the time evolution of V(M1) and V(M2) for an example case A = 400 bar. As seen, these potential energies oscillate with the ultrasound, i.e, their values increase and decrease using the ultrasound rarefaction and compression phases, respectively. However, on average, V(M1) is ∼ -2000 kJ/mol lower than V(M2), suggesting that the monolayer M1 is more stable than M2. Therefore, while the M1 layer is anchored by the tight junction, the M2 layer moves following the expansion of the system during the ultrasound rarefaction phase, resulting in the detachment of the layer M2 from the layer M1.

In the real BBB, the tight junction is much more complicated, composed of many claudin-5 proteins and other proteins. Of course, this will increase the tightness of the junction but also increase the interaction between cell membranes with transmembrane proteins of the tight junction. The question is then whether the BBB is still disrupted at the cell membrane position? To this end, we construct an empirical BBB model whose tight junction is composed of four pairs of two cis-dimers, each pair is located on the corner of a square with edge lengths of 10 nm [Fig. 9]. This is obtained by translating a pair of two cis-dimers of the above model [Fig. 1] in the membrane surface for 10 nm along the x and y-axis. We then carry out simulations using the same parameters as used for the previous model. Fig. 9 shows snapshots after six ultrasound irradiation periods, i.e. at 120 ns for A = 400 and 450 bar. Consistent with the results of the single cis-dimer pair model, the structure of the four cis-dimer pairs model is well-maintained with A = 400 bar, but with A = 450 bar, two monolayers of individual bilayer membranes is detached, suggesting that the BBB is also damaged at the cell plasma membrane rather than at the tight junction.

Let us discuss how the separation of monolayers of the cell membrane results in the enhancement of the drug permeability. The separation of two monolayers in a bilayer results in large spaces at the middle of the lipid bilayer (Figs. observed experimentally, in the context of drug delivery studies using liposomes. That is, under LIFU irradiation, gas bubble nuclei may be formed in the hydrophobic region of the lipid bilayer of the liposomes. These nuclei grow until they permeate the membrane, forming a transient pore through which the drug is released [START_REF] Huang | Acoustically active liposomes for drug encapsulation and ultrasound-triggered release[END_REF][START_REF] Schroeder | Ultrasound, liposomes, and drug delivery: Principles for using ultrasound to control the release of drugs from liposomes[END_REF]. At the molecular level, the presence of the large spaces alternates the diffusive process of drugs across the membrane. Indeed, a number of simulations have shown that the free energy profile for translocation across the lipid bilayer of drugs usually exhibits a high barrier at the middle of the two monolayers [START_REF] Shinoda | Permeability across lipid membranes[END_REF][START_REF] Orsi | Permeability of Small Molecules through a Lipid Bilayer: A Multiscale Simulation Study[END_REF][START_REF] Xiang | Molecular dissociation processes in lipid bilayers: a molecular dynamics simulation[END_REF][START_REF] Wang | Expoloring gas permeability of cellular membranes and membrane channels with molecular dynamics[END_REF]. If two monolayers are separated, then this leaves a large space in the middle of the bilayer, resulting in high population of drugs in this area, i.e. the free energy barrier is reduced, and therefore the drug permeability is increased. We find that the molecular mechanism of the ultrasound induced permeability of BBB is similar to that observed in liposomes, but the origin of the large spaces formation may be different. It is expected that considerable heterogeneity exists in the preparations of liposomes, and it is unlikely that all li-posomes are unilamellar. Thus, large spaces are likely formed at the flexible areas of the membrane of liposomes. For the BBB, because the outer layer of the cell membrane is anchored by the tight junction, but the inner layer is contracted during the rarefaction phase of the ultrasound, therefore this separation results in large spaces inside the cell membrane bilayer [Fig. 8].

As reviewed earlier in this paper, under HIFU irradiation, several experiments observed the BBB disruption within the damaged areas [START_REF] Bakay | Ultrasonically produced changes in the blood brain barrier[END_REF][START_REF] Patrick | Ultrasound and the blood-brain barrier[END_REF], but other results showed the BBB disruption without damage of the brain parenchyma [START_REF] Ballantine | Progress and problems in the neurological applications of focused ultrasound[END_REF][START_REF] Vykhodtseva | Histologic effects of high intensity pulsed ultrasound exposure with subharmonic emission in rabit brain in vivo[END_REF]. In both cases, it was unclear whether the BBB disruption is due the cell damage or opening of the tight junction [START_REF] Mesiwala | High-intensity focused ultrasound selectively disrupts the blood-brain barrier in vivo[END_REF]. Our simulations could suggest that the BBB may be disrupted due to the damaged of the cells, but the tight junction is unlikely happened.

A large body of literature confirms that LIFU-induced BBB opening involves at least four mechanisms: transcythosis, transendothelial opening, i.e. fenestration and pore formation, tight junction opening at low ultrasound intensities, and passage through damaged cell at high ultrasound intensities [START_REF] Sheikov | Cellular mechanisms of the blood-brain barrier opening induced by ultrasound in presence of microbubbles[END_REF][START_REF] Sheikov | Brain arterioles show more active vesicular transport of blood-borne tracer molecules than capillaries and venules after focused ultrasound-evoked opening of the blood-brain barrier[END_REF]. Our simulation study suggests that with low ultrasound intensities used in LIFU experiments, the impact of the direct interactions of ultrasound with cells and with tight junctions can be negligible. Thus, the interaction between bubble and the BBB plays an essential role. Indeed, it has been suggested that the acoustic microstreaming generated by stable bubble cavitation stimulates the transcythosis, transendothelial opening, and the inertial cavitation may induce cell membrane pore formation. The tight junction opening may due to the volume expansion of microbubbles in the vessel, which in turn exerts direct mechanical forces on the BBB. However, at high ultrasound intensity, the direct interaction of ultrasound with the BBB becomes significant, leading to the cell damage as observed by experiments [START_REF] Sheikov | Cellular mechanisms of the blood-brain barrier opening induced by ultrasound in presence of microbubbles[END_REF][START_REF] Sheikov | Brain arterioles show more active vesicular transport of blood-borne tracer molecules than capillaries and venules after focused ultrasound-evoked opening of the blood-brain barrier[END_REF], and confirmed by our simulations.

In conclusions, we have performed a molecular dynamics simulation study to understand the molecular impacts of HIFU and LIFU on the BBB opening. The study shows that even with a relatively high ultrasound intensity of 400 bar, the structures of the cell membranes and tight junction were well-maintained, implying that the direct interaction between LIFU with the BBB is not responsible for the BBB opening, that has been observed experimentally. At higher intensities (A ≥ 450 bar), the rarefaction of ultrasound pulls two monolayers of individual cell membrane lipid bilayers apart, creating air compartments inside the bilayer. This reduces the free energy barrier for translocation across the lipid bilayer of drugs, thus enhances the drug permeability. At very high intensities, two monolayers are largely separated, resulting in cell damage, implying that the BBB is primarily disrupted at the damaged areas as observed by experiment. The MD simulations of the full system composed of ultrasound, bubbles and BBB will be carried out to understand impacts of bubble cavitation on the This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. BBB model. Also, more detailed BBB models will be constructed. These works may open the door of drugs delivery for Alzheimer's diseases [START_REF] Nasica-Labouze | Amyloid beta-protein and Alzheimer's disease: when computer simulations complement experimental studies[END_REF][START_REF] Doig | Inhibition of protein aggre-gation and amyloid formation by small molecules[END_REF]. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS

PLEASE CITE THIS ARTICLE AS

Figure 1 :

 1 Figure 1: (A) The molecular structure, viewed from two angles, of two claudin-5 cis-dimers (green and orange) forming a paracellular pore. (B) The molecular structure of our minimal BBB model where cis-dimers are embed into two lipid membranes of two neighbouring cells. The water is filled between and outside two lipid bilayers.

Figure 2 :

 2 Figure 2: Time evolution of various quantities, including the ultrasound (A), the system volume (B), the pressure (C), the potential energy (D) and the temperature (E). Shown are results obtained with different ultrasound intensities A = 100, 200, 400 and 450 bar. The period of the ultrasound is τ = 20 ns (50 MHz).

Figure 3 :

 3 Figure3: The time evolution of the RMSD of the claudin-5 tight junction with respect to the initial structure (A), and the total number of intermolecular Hbonds between two dimers (B). Shown are results obtained from simulations using ultrasound with period of 20 ns and various intensities.

Figure 4 :

 4 Figure 4: The time evolution of the population of various secondary structures. For each structure, shown is the total populations of four claudin-5 proteins pertaining to the junction. Shown are results obtained from simulations using ultrasound with period of 20 ns and various intensities.

2 Figure 5 :

 25 Figure5: Time evolution of the carbon-hydrogen order parameters of the first sn-1 (upper) and second sn-2 lipid tails, respectively (lowver). These results are averaged over 15 and 17 order parameters of C-H vectors of the first and second lipid tails, respectively, and obtained from simulations using ultrasound with period of 20 ns and various intensities.

Figure 6 :

 6 Figure 6: Time evolution of the shortest intra-distances between two monolayers of the first (A) and second (B) cell membrane lipid bilayers. The shortest inter-distance between two lipid bilayer cell membranes is shown in the (C). Shown are results obtained from simulations using ultrasound with period of 20 ns and various intensities.

Figure 7 :

 7 Figure 7: Snapshots of the system at 120 ns obtained by using ultrasound intensity of 400 bar (A), and at 24.5, 24.7 and 25 ns obtained by using ultrasound intensity of 450 bar (B, C, D respectively). The ultrasound period of 20 ns is used in all simulations. While the bilayers are disrupted, the tight junction is well-maintained.

Figure 8 :

 8 Figure 8: (Upper) A structure of the BBB model where components are labeled as M1 and M2 for two monolayers of a cell membrane bilayer; W1 and W2 for waters between two bilayers and above the membrane M2, respectively; C1 and C2 for two claudin-5 dimers embed in two cell membrane lipid bilayers. (Lower) Time evolution of the potential energy acting on M1 (black) and M2 (red). For clarity, only phosphate groups of the lipids are shown in beads and waters are not shown.

  DOI:10.1063/5.0010667

Figure 9 :

 9 Figure 9: Snapshots of the system whose tight junction is composed of four pairs of claudin-5 cis-dimers at 120 ns obtained by using ultrasound intensity of 400 bar (left), and at 25 ns obtained by using ultrasound intensity of 450 bar (right). The ultrasound period is 20 ns in all simulations. Waters are not shown for clarity. The separation of two monolayers is clearly seen from the snapshot at 450 bar. For clarity, only phosphate groups of the lipids are shown in beads and waters are not shown.

  DOI:10.1063/5.0010667This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI:10.1063/5This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI:10.1063/5This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI:10.1063/5the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI:10.1063/5This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI:10.1063/5This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI:10.1063/5

Acknowledgments This work has been supported by the Department of Science and Technology at Ho Chi Minh City, Vietnam (grant 10/2018/HD-KHCNTT), CNRS, the World Bank and the Vietnam Ministry of Science and Technology FIRST project (grant 13/FIRST/1.a/VNU1), the Polish Narodowym Centrum Nauki (NCN, grant 2019/35/B/ST4/02086), the National Science Foundation (NSF, grant SI2-SEE-1534941), the National Institute of Health (NIH-R01GM118508) and the CINES center for providing computer facilities (project A0080711440).

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


Romain Courault
email: romain.courault@ird.fr

Marianne Cohen
email: marianne.cohen@sorbonne-universite.fr

Mathilde Pottier
email: mathilde.pottierulr@gmail.com

Disentangling the influence of geographical laws and sampling biais to model distribution of Birch tree from Open-acess Biodiversity Dataset (OBDs) in Swedish Lapland

Keywords: Biogeography, Geographical Laws in Spatial Ecology, Biodiversity Databases and Datasets, GBIF, Birch tree, Swedish Lapland, Habitat Niche Modelling, Climate Change

The current biodiversity crisis, combined with climate change are major issues requesting specific monitoring of plants communities' responses in terms of geographical distribution. Nowadays, large open-access biodiversity datasets (OBDs) such as the Global Biodiversity Information Facility (GBIF) are commonly used to describe, explain and predict fauna and flora geographical distribution. They constitute new opportunities, but stay related to major uncertainties about sampling biases, driven by the concentration of various biodiversity data records and associated data providers. Taking the example of a widely studied tree (Betula pubsescens Ehrh.), in a scientifically well-funded region (Swedish Lapland), we discuss those peculiar issues in the frame of geographical laws (e.g. respectively spatial autocorrelation, heterogeneity and similarity) at macro-and micro-regional scales. After spatial and temporal filtering on georeferenced records and discussion on sampling strategies heterogeneity, tests of spatial autocorrelation (Moran's I index) has been conducted on Birch tree records provided by major institutions, comparatively. Pearson Khi-2 (χ 2 ) test is thus applied on the generated grid to confront number of Birch tree records with accessibility factors (e.g. artificial land cover, roads, protected natural areas) Thus, a micro-regional analysis is conducted to quantify Birch tree records in vegetation classes where this tree species is supposed to be dominant. At the macro-regional scale, results show the high spatial variability of sub-datasets according to institution providing records from the studied GBIF OBDs (with higher autocorrelation results for large contributors). This spatial variability, and high spatial autocorrelation effects appears to be partly explained at macro-, micro-and local scale by the distribution of human accessibility and facility factors (e.g. roads, cities etc). In this study case, exploring OBDs for an extensively addressed tree species, in a significantly funded region/country was particularly useful, demonstrating the relevance of spatial similarity law to differentiate adequately sampling biases efforts from "natural" spatial autocorrelation.

I. INTRODUCTION

A. Geographical distribution of species, geographical distribution of sampling biases : deciphering tracability of Open-access Biodiversity Datasets records with general principles in geography

At a time when there is an urgent need to model and simulate geographic space, and at the same time, a large amount of data is available online to feed these models, the geographer has to ask himself the following question: how to disentangle the influence of geographic laws and sampling bias related to these large open-source data and avoid misinterpretation? This is particularly true for the modelling of vegetal species distribution, particularly for food species bringing important ecosystem services for the society. The distribution range of these species is expected to move with climate change, with dramatic changes for socio-ecosystems, interacting with biodiversity crisis [START_REF] Pearson | Shifts in Arctic vegetation and associated feedbacks under climate change[END_REF], particularly in arctic and subarctic regions [START_REF] Pfeifer | 1.5, 2, and 3 global warming: visualizing European regions affected by multiple changes[END_REF], IPCC 2013, IBPES 2019). The current global environmental crisis constitutes an urgent need to model and simulate geographic space to get prepared and seeking for resiliency at all geographical levels. Meanwhile, publicly available geographical data about environmental conditions and biodiversity (also known as Primary Biodiversity Data, hereafter named Open-access Biodiversity Datasets, OBDs), Anderson, Araùjo et al., 2020) are currently easily and freely accessible, coming from various data aggregators such as GBIF (2016), DataDryad, IUCN, Neotoma, NatureServe, Conservation of Arctic Flora and Fauna, CHELSA-Project, WorldClim, etc [START_REF] Sillero | Common mistakes in ecological niche models[END_REF].

Under these conditions, it is necessary for geographers and other disciplines interacting with space management (biology, ecology, engineering, economics; political and law sciences etc) to question and to put in perspective the reliability of such easy-to-process datasets and associated software [START_REF] Sillero | Common mistakes in ecological niche models[END_REF]Zizka, Antonelli et al., 2020). Without mentioning the influence of spatial correlation and covariation between species locations, and interpolated abiotic factors/datasets used to model patterns of distribution or migration, it is thus particularly crucial to sum on to the influence of geographic laws in sampling bias related to OBDs and having hindsight on those datasets, quantitatively important in georeferenced records (~10 9 individual species records, Andersen, Araùjo et al., 2020) but qualitatively peculiarly questionable when modeling those individuals, taken without any ponderation when processed [START_REF] Sillero | Common mistakes in ecological niche models[END_REF]. The other side is also true, three mixed groups contributing to nourish, aggregate and use OBDs (data providers, data aggregators, data users), with little respect (or no time) to carefully providing guidance on first-hand biodiversity datasets, which will be later used for a numerous of usages (Andersen, Araùjo et al., 2020).

In an absolute way, OBDs mainly consists on georeferenced records of specie, which are in a first sight one location of one individual belonging to one taxa (fauna, flora, fungus, etc). Lack of information regarding metadata (e.g. year, data acquisition conditions, sampling efforts, and research questionings) constitutes gaps and biases, which are common for OBDs. They can be summarized in five interlaced types [START_REF] El -Gabbas | Improved species-occurrence predictions in data-poor regions: using largescale data and bias correction with down-weighted Poisson regression and Maxent[END_REF]Zizka & Antonelli, 2020):

-Biological biases (taxonomic bias represented by biologist/ecologist selectivity about species of interests, taxonomic misidentification, quality monitoring not being provided)

-Spatial biases (remote areas and proximity effects; errors when geo-referencing/inaccuracies and approximations on field and hereafter; georeferenced collections reflecting collection/museum location instead of where those occurrences have been found; "road-side-bias" effect and more widely physical accessibility and)

-Environmental biases (harsh environments, "avoided" studied areas because of their difficult conditions for researchers) -Temporal biases (absence of temporal information; inaccuracy, seasonal biases when sampling)

-Socio-economic and political biases (sampling biases linked to research history and funds, conflicts, human rights)

A striking example of multi-factorial sampling bias at a global scale is the more important number of GBIF records in well-financed countries, compared to species-rich countries [START_REF] El -Gabbas | Improved species-occurrence predictions in data-poor regions: using largescale data and bias correction with down-weighted Poisson regression and Maxent[END_REF]. Models and their cross-validation, and more deeply firsthand datasets as input of models, such as available in OPDs constitutes important opportunities which should be used in the framework of spatial distribution laws and documentation on their use in such framework disseminated (Cotello et al., 2013in Anderson, Araujo et al. 2020). Another important gaps of knowledge found in the literature is represented by the "few attempts [which] has been made to compare the geographic sampling bias among datasets […] [to quantify] the effect size of specific bias factors and compare it among them" (Zizka, Antonelli et al., 2020).

B. Birch tree (Betula pubescens Ehrh. 1791), a key-stone species of subalpine woody tundras

Biogeographically, the Birch tree (Betula pubescens Ehrh. 1791) is widespread across Eurasia in mountainous areas of the boreal zone (Tutin, Burges et al., 1993). According to [START_REF] Andersson | The last European landscape to be colonised: a case study of land-use change in the far north of Sweden 1850-1930[END_REF], Betula pubescens is present in Swedish Lapland with two subspecies: Betula pubescens ssp.

pubescens, named Birch tree and Betula pubescens (Ehrh) ssp. tortuosa (Ledeb), named mountain Birch. Previously they were considered as two different species named Betula pubescens Ehrh. and Betula czerepanovii n.i. Orlova). The first one is significantly present in the boreal forest, in monospecific stands or as subdominant species in conifer forests and more generally in Middle and Northern Eurasia. The second one is present on the North-Western Atlantic fringe of the distribution map of Betula ssp., at the foothills of the Fennoscandian chain. Unfortunately, in the OBDs used in this paper, the information about subspecies is available for only 8% of the records. Subspecies cannot be distinguished at the regional scale of Swedish Lapland, but only Birch tree is represented at the micro-regional and local scales of our study (section II.A).

The high genetic variability of mountain Birch renders it particularly flexible, although the species is hygrophilous and prefers clay, explaining its affinity with peats, moorlands to be later colonized by it. Patchworks of wood Birch tree woods are a-zonal subalpine tundra ranging from low to middleelevation in northern Europe and Siberia. Indeed, the tree species and associated floristic communities form an ecotone between northern boreal forests and scattered, perhaps vascular-free tundra at high altitude and/or high latitude (Rydin, Snoeijs et al. (eds), 1999). Along the northern boreal zone in Eurasia, Betula genus are often considered as the last trees to survive to such cold conditions, and form then treelines. Nowadays, observational and experimental studies in palynology, genetic and ecology have demonstrated the high variability of the Birch tree-line, and its high potential to invade higher altitudes, in partly due to the direct effect of warming on recruitment [START_REF] Kullman | Tree limit dynamics of Betula pubescens ssp. tortuosa in relation to climate variability: evidence from central Sweden[END_REF][START_REF] Truong | Recent invasion of the mountain birch Betula pubescens ssp. tortuosa above the treeline due to climate change: genetic and ecological study in northern Sweden[END_REF][START_REF] Bryn | 20th century Betula pubescens subsp. czerepanovii tree-and forest lines in Norway[END_REF]Kullman, 2017). This present study, focused on Betula pubescens, including its two subspecies, aims to better assess the present and future distribution of the species, being the interlaced links between tree biomass and palatable vegetable for reindeer and for reindeer herding systems. Indeed, even if young leaves of mountain Birch might be consumed by subarctic and arctic herbivorous, and particularly semidomesticated reindeer [START_REF] Forbes | The ecological role and geography of reindeer (Rangifer tarandus) in northern Eurasia[END_REF], impacts of invasion and densification of Birch tree woods on biotopes and other biological communities might be important (Courault, 2018b;[START_REF] Agnan | Distribution of trace and major elements in subarctic ecosystem soils: Sources and influence of vegetation[END_REF]. Although direct inter-specific interactions are complex between Mountain Birch and reindeer, tree-line rising, biomass densification and phenological aspects of the tree are hypothetically susceptible to interfere on reindeer biological rhythm and migration in spring and summer (Kumpula, Stark et al., 2010). Beyond that, it is important to survey and assess the species distribution in the past, present and future, among which habitat modeling niche techniques could be useful to infer and predict on its future distribution at different ecological and geographical scales.

C. Geographical laws in the frame of Openaccess spatial ecology : questioning and objectives of the study

Spatial autocorrelation is a property of the First Law of Geography, Tobler's Law 'Everything is related to everything else, but near things are more related than distant things' (Tobler, 1970[START_REF] Sillero | Common mistakes in ecological niche models[END_REF]. Auto-correlation is necessary for modeling descriptive and explanatory species distribution areas [START_REF] Segurado | Consequences of spatial autocorrelation for nichebased models[END_REF], Dormann 2007, Dormann et al. 2007, De Marco et al., 2008[START_REF] Sillero | Common mistakes in ecological niche models[END_REF], but misleading because of biases and quality heterogeneity of OBDs. Nevertheless, for ecologists and biogeographers, spatial autocorrelation might appear to a vague concept, often confused with spatial clustering of records and partly diluted with filtering records [START_REF] Sillero | Common mistakes in ecological niche models[END_REF]. Hence, spatial clustering might appear to be either a result of 1) data aggregation of various sources 2) sampling efforts coming from different data providers 3) species' distribution in its natural habitat.

The second point linking biodiversity research, such as habitat niche modeling, and geographical laws is the essential bridge between biotic and abiotic relationships in space. Strong environmental gradients are theoretically related to more realistic models (e.g. increased explanatory or predictive power of the model (Seoane et al., 2005[START_REF] Sillero | Common mistakes in ecological niche models[END_REF], based on the geographic similarity principle, corresponding to the "contagious model" opposed to the "random model" (spatial heterogeneity) of species distribution by [START_REF] Gounot | Méthodes quantitatives d'étude de la végétation[END_REF]. Indeed, Varela et al. (2014) showed that filtering by environmental criteria provides better results [START_REF] Sillero | Common mistakes in ecological niche models[END_REF]. This justify the study area as well as the studied species, Betula pubescens.

Into this conceptual framework, the case study is based on Birch tree (Betula pubescens Ehrh.), whose biological/ecological characteristics are relatively well studied and geographically enquired at the scale of its current distribution. Thus, our hypothesis is that spatial biases coming from OPDs are still detectable and won't provide sufficient and reliable geographical information to firmly describe and predict its present and future distribution using habitat niche modeling techniques. In particular, accessibility and facility factor would be likely to be quantified at various geographical scale, such as macro-, micro-and local levels in Swedish Lapland.

The objective of this paper is to test to what extent and under what conditions OPDs' sampling biases can be confused with geographical laws, and to propose methods to separate their respective influence. We will present the case study, the material and the methods mobilized (section II), then the results obtained at 3 scales, macro-regional, microregional and local (section III) and will conclude with a discussion and the methodological perspectives opened for geographers by this study.

II. STUDY AREA, MATERIAL AND METHODS

A. Macro-and micro-regional study areas

The macro-regional study area, represented by the Swedish Lapland (e.g. Sameby in northern Sàmi) is complex ecologically, historically and geographically [START_REF] Forbes | The ecological role and geography of reindeer (Rangifer tarandus) in northern Eurasia[END_REF]Ojala, 2009, Löf, 2013). Administratively speaking, the macro-regional studied area stretches from Svealand to Norrland old Swedish provinces, and comprise as well some Norwegian territories (= 51 reindeer herding communities, e.g. samebyar). Indeed, the studied area corresponds to the Swedish part of the Lapland whose ecosystems and landscapes were economically, culturally and socially constructed and valued by Sàmi people, and fiercely negotiated with Scandinavian authorities along the last centuries [START_REF] Manker | The nomadism of the mountain lapps: the siidas and their migratory routes in 1945[END_REF][START_REF] Svonni | Samisk markanvändning och MKB[END_REF]. Having Finno-Ugrian origins, linguistic heritage and culture, cultural landscapes have been (and are still) mainly devoted to migratory reindeer herding, with a specialization of pastoral system according to clans and families' belongings.

Landscapes of Swedish Lapland, as synoptically represented in the Figure 1, present a diversity according to their location within 4 main gradients/ecoclines: latitudinal, longitudinal, continental, altitudinal [START_REF] Courault | La fragmentation des paysages de l'élevage des rennes. Une étude de cas en Laponie suédoise[END_REF][START_REF] Courault | Les Paysages Culturels de L'élevage de Rennes en Scandinavie Face au Changement Global, Une Approche Multi-Scalaire (Laponie Suédoise[END_REF]. More precisely, vegetation zonation of Swedish Lapland ecosystems expands from the Southern Boreal zone, to Middle-and Northern Boreal subzones, to Subalpine and Alpine belts compose landscapes of the Fennoscandian mountains, devoted to calving areas and summer pastures for reindeers. Zonal as well as a-zonal vegetal landscapes are highly dependent on the climate (colder above the Arctic Circle, winder and milder in glacial valleys influenced by the Gulf Stream), mineral soils (and associated bedrocks), hydrology and (micro-) topography (Rydin, Snoeijs et al. (eds), 1999). At this level, describing the current and future distribution of birch tree (Betula pubescens Ehrh.) is particularly interesting in both fundamental and applied research in biogeography, being the last "tree" to remain at the Subalpine belt and the frame of current climate change.

The micro-regional study area, as displayed in the Figure 1, is located in the northern edge of the Swedish Lapland. Four samebyar (reindeer herders' communities) are represented within this area, and are characterized by mountainous reindeer herding pastoral systems (e.g. wider vital areas of semi-domesticated reindeers, wider activities of owners, and wider contemporary disturbances ;[START_REF] Manker | The nomadism of the mountain lapps: the siidas and their migratory routes in 1945[END_REF][START_REF] Courault | La fragmentation des paysages de l'élevage des rennes. Une étude de cas en Laponie suédoise[END_REF]. In terms of vegetal landscapes, this area is particularly interesting because of its transitional situation between the Boreal zone, where Scots pines are highly exploited, as well as the alpine and subalpine belts, where winter and summer tourism might constitute troubles and tensions for reindeer pastoral systems as well as for local biotic communities (Skarin, Danell et al., 2008). This microregional study area has been selected for multiple reasons, because of the naturalist entanglement between "traditional" activities, modern ones, juridical issues in terms of environmental protection and wildlife watching by scientists (INTERACT-network of natural science stations). Furthermore, the area has been chosen because where field work carried out by [START_REF] Courault | Les Paysages Culturels de L'élevage de Rennes en Scandinavie Face au Changement Global, Une Approche Multi-Scalaire (Laponie Suédoise[END_REF] provides particular expertise on OBD quality and Betula genus is represented by Mountain Birch (Betula pubescens ssp tortuosa) Fig. 1. Sudied areas in Sweden: macro-regional scale (Swedish Lapland, Sameby), and micro-regional scale of analysis of mountain Birch records from GBIF dataset (Sources: Swedish Saami Parliement, RenGIS 2.0 sametinget.se; Lantmäteriet, SLU; realization R.Courault 2020)

B. Material 1) Acquisition of georeferenced Birch tree records, semi-automatized filtering and study area scaling

Using the Rgbif package (R Studio), georeferenced occurrences of Birch tree (Betula pubescens Ehrh.) were downloaded at the scale of Eurasia (e.g. 100 000 records; GBIF). For spatial biases which are analyzed here, cleaning records (e.g. removing removing and keeping occurrences of one or several species) using distances and uncertainties of their locations appears to be one solution "to dilute the effect of uneven sampling effort across the study area" (Anderson and Raza, 2010, Boria et al., 2014in El-Gabbas & Doorman, 2018). With the Clean Coordinate package (R Studio) spatial tests were applied to flag and to remove:  Record without any coordinates  Records belonging to museum' and institutional' collections, herbarium and revealing their location instead of the exact place where records were sampled (non-human in-situ observations)

 Non accurately georeferenced records, be it for national capital, province/administrative region centroids  Geographical locations with more than one record  Spatial uncertainty of coordinates being above the threshold of 1000 meters (Figure 2).

Once spatial tests conducted using R packages, georeferenced occurrences Birch tree has been imported in Qgis 3.14.16 to intersect Eurasian occurrences at the extent of Swedish Lapland, considered as the extent of all Swedish reindeer herding communities' borders (source: sametinget.se; RenGIS 2.0). Thus, after dataset visual control (Figure 2), an attribute query deleted georeferenced Birch tree being recorded before 1995, until today (2020). This temporal filtering has mainly ecological/biogeographical goals in terms of habitat niche modeling. The preprocessed dataset gives then 12 265 records (Table 1). Fig. 2. Example of data visualization and balance between sampling effort (on field) of Birch tree GBIF dataset, its geographical uncertainties on the one hand and statistical needs for records and niche species modelling on the other hand. Note the systematic and squared sampling effort/strategy, in the southwestern part of the study area, presenting a significant spatial uncertainty (e.g. 2.5km , institution: ArtDataBanken). Nevertheless, this systematic sampling minimize the bias due to presence only data set, the species absence is noticeable in two areas. Even if overall autocorrelation (scale of the study area) might be over-represented, geographical information missed by an arbitrary cleaning constitutes a serious compromise to take into account in ecological modeling. udied areas in Sweden: macro-regional scale (Swedish Lapland, Sameby), and micro-regional scale of analysis of Birch tree records from GBIF dataset (Sources: GBIF, Swedish Saami Parliement, RenGIS 2.0 sametinget.se; Lantmäteriet, SLU; realization R.Courault 2020)

2) Acquisition of geographical data on land cover, accessibility and facility factors

The portal Copernicus -Land Monitoring Service (https://land.copernicus.eu/pan-european/corine-land-cover) has been used to acquire the Corine Land Cover map of 2018 (CLC 2018), in its shapefile format. When pasting the CLC 18 at the macro-regional extent of the study area, an attribute query is made to select land cover categories related to artificial zones and human facilities, excluding other land covers from the analysis: For roads and railways, divagis.org has been used to get Norwegian and Swedish transportation networks in shapefile (https://www.diva-gis.org/gdata).

The last two datasets (Natural protected areas and roads and railways) have been intersected with the extent of the Swedish Lapland, and merged in particular for transportation networks.

Figure 3 shows that human observations of Betula pubescens are concentrated in the boreal forest, but extend to the South-western Botnia Gulf (Birch tree) and to the Northwestern foothills and valleys of the Fenno-Scandian chain (Mountain Birch). Geographical bias effects are visible at the scale of a territory. For those recorded human observations of Birch tree, concentration effects are noticeable, e.g. Norrbotten county and around Östersund, respectively in the northern and southern parts of the macro-regional study area.

C. Methods 1) Spatial and statistical analysis at the macroregional scale

For the macro-regional scale and analysis conducted on Birch tree records and tested geographical features (autocorrelation, accessibility and facility factors), grids computing has been chosen for several reasons. Firstly, grid cells avoid statistical noise which could be due to pre-or postprocessing, such as interpolation methods [START_REF] Ruete | Displaying bias in sampling effort of data accessed from biodiversity databases using ignorance maps[END_REF]. Secondly, it technically simplifies counting number of occurrences of both studied phenomenon (here Birch tree) and explanatory factors (its distribution among institution providers, as well as accessibility/facility factors which could potentially explain Birch tree distribution as recorded by GBIF). a) Grids generations and absence/presence countings Using the tool "Create grid" from Gdal on QGis 3.14.16, grids were generated at the spatial extent of each dataset corresponding to institutions records of the Birch tree, at the wider scale of the study area (e.g. Swedish Lapland). As represented in figure 4, the wider grid used to count records of Birch tree and associated factors is displayed by the Artdatabanken providing institution. Other generated grids correspond to the extent of each providing institution. To include the modal spatial uncertainty given by records in the GBIF dataset, spatial resolution of grids has been set up to 1*1 km for every grid (see Table II for number of columns and rows; projection system of the GIS: SWEREF 99 TM).

For such macro-regional level analysis (autocorrelation test and influence of accessibility/facility factors) those counting include as well other geographical factors as listed in the Material section: artificial land cover (CLC 18), Natural protected areas, roads and railways.

b) Moran's I spatial autocorrelation test per providing institution

Morans'I spatial autocorrelation tests have been conducted by differentiating records of Birch tree Betula pubscens Ehrh. from various institutions sources. As shown in the table I, number of occurrences greatly vary accordingly. In this way, only institutions having recorded more than 20 occurrences of Birch tree are compared in Moran's I spatial autocorrelation test. Those are listed below and displayed in figure 4 I, figure 4)

The tool "Spatial Autocorrelation (Moran's I)" of ArcMap 10.4 is used to quantify the degree of spatial autocorrelation / dependence between records of the same dataset (e.g records emanating from the same institution) and to compare results between institutions having recorded Birch tree at the scale of the Swedish Lapland.

Because of the polygonal geometry induced by grids (each cell composing grids having a number of records), spatial relation conceptualization parameter has been set to 'Contiguity Edges Corners'; these to better assess spatial relationships between contiguous cells (sides, or cells).

We chose to standardize processed grids and associated counts, to better allow comparisons between Moran's I spatial autocorrelation tests applied for the different institutions. Row standardization is set because generated grids, even though having the same granularity (e.g. 1000 * 1000 m) differ in spatial extents, depending on records distribution/institutional sub-dataset.

c) Chi-2 (χ2) test between number of Birch tree records and accessibility/facility factors

Always at the macro-regional scale (e.g. the entire Swedish Lapland), we were willing to test the degree of dependence between all records (without any distinction of data provider/institution) and factors reflecting scientific/human accessibility and facilities to fieldwork and naturalist observations (e.g. transportation networks, artificial areas such as housings, natural protected areas as an attractive factor, see subpart Introduction). We used the grid generated at the level of all Birch tree records (e.g. grid having the most important extent; "Grids generations subpart"), which cover the entire Swedish Lapland (e.g. 489 291 cells of 1km*1km in SWEREF-99 TM) to check the number of co-occurrences of Birch records and transportation, protected areas and artificial zones. Once counting done, the attributary table has been exported as a CSV file, and then processed with XLStat 2020 to produce a contingency table on those ordinal dependent / explicative variables. A Pearson Chi-2 (χ2) is then calculated on the contingency table, giving as outputs χ2 by cell, observed and theoretical frequencies, observed, critical and p-values of the overall χ2.

2) Spatial and statistical analysis at the micro-regional and local scales a) Chi-2 (χ 2 ) test between number of Birch tree records and vegetal communities categories

At the micro-regional scale (figure 1), we wanted to test the independence hypothesis between records of Birch tree and type of ecological habitat as described by the Swedish ecological map (see subpart Material). Thus, number of occurrences by polygon and vegetation communities' categories have been enumerated and summarized using both Qgis 3.14 and XlStats 2020. A contingency table has been produced for the 44 different habitats. Artificial areas were excluded from the analysis. χ 2 test is then conducted on those categorical variables.

b) Study case of Caluna AB data provider at a local scale with photointerpretation

A study case has been conducted at a local scale in the right bank of the Torneträsk lake, near from the Research Station in Environmental Sciences and the eponymous Natural Park (Abisko naturvetenskapliga station, INTERACT network and Swedish Polar Research Secretariat, Kiruna municipality, studied area ~1200 ha). GBIF records of Birch tree come from Calluna AB (Table I), aligned along the railway track, which nearly follows the E10 highway. Contours of the tree formations has been digitized (~970 polygons), apart from one cloudy area standing at the center of the study area (source of the imagery: Google Earth, exploited by ArcGIS Pro 2.6). Two subareas, located respectively to the west and to the east, of comparable surface area (~600 ha) are then analyzed by 1) calculating number of Birch tree and their distance in meters to the railway line; frequencies of records by category of photo-interpreted land covers (e.g. Birch forests, water bodies, snow patches) and associated percentages of records by hectares of forests.

III. RESULTS

A. Spatial and statistical analysis at the macro-regional scale a) Moran's I spatial autocorrelation test per providing institution

Table II summarize Morans' I spatial autocorrelation tests by institutional providers having quantitatively contributed to the GBIF dataset of Birch tree at the studied macro-regional scale. The highest Morans' I index is found for Artabanken records (i=+0,15; z-score = 212; p-value <0.0001). Given the geographical distribution of Birch tree records (figure 3 and4), two main concentration effects can be detectable, the first almost corresponding to the borders of the Norrbotten county, the second cluster covering the administrative borders of the Jämtland county. The last one condense a lot of records, as confirmed by the second highest Morans' I index for Länsstyrelsen Jämtland (records provided by this institution being in blue, figure 4). Thus, in the southern part of Swedish Lapland, effects of the accessibility on sampling biases appear to be particularly strong, records being intensively clustered around the administrative capital of the county, Östersund. More broadly, other providing institutions might have clustering effects quantified by Morans' I test, although being too weak or even non-significant (for Morans' I test, z-score, or p-value, see Table II).

b) Chi-2 (χ 2 ) test between number of Birch tree records and accessibility/facility factors

As already mentioned, biodiversity records as available in OBD such as GBIF, are often the reflect of sampling strategies and efforts, whose might be biased by accessibility and facility factors such as transportation, housing, Natural protected areas [START_REF] El -Gabbas | Improved species-occurrence predictions in data-poor regions: using largescale data and bias correction with down-weighted Poisson regression and Maxent[END_REF], Zizka, Antonelli et al. 2020). Figure 5 plots Birch tree records counted at the macroregional extent, cells of the biggest macro-regional grid summarizing number of records of Birches tree, Artificial areas, Road and Railways, Natural protected areas. Quite obviously, number of cells having one or more presence of roads/railways are related to the number of artificial areas as the scattered plot of the Figure 5 shows (upper-right top and lower-left bottom of the Figure 5). For relationship between number of presence of Birch trees and external factors (second column of the Figure 5), presence of artificial areas, roads and Natural protected areas might be classed as follow, given number of presences: 1) presence/number of roads 2) presence/number of artificial areas 3) presence/number of Natural protected areas. This is confirmed by respective standard deviations. This is as well confirmed by χ 2 tests applied on this macroregional dataset, which are highly significant (Table III). Birches trees are significantly over-observed in the grid cells containing a road, an artificial zone or a protected area, in descending order of observed χ 2 value.

The over-representation of birch observations near roads (93% of the observed χ 2 value) can be considered as a sampling bias, as the road facilitates the human observations. However, roads are built in flat valleys on moderately wet terrain to facilitate construction and maintenance and to increase driver safety on straight and level roads. In Swedish Lapland, such glacial valleys correspond to the most frequently common habitat and geographical distribution of the tree.

Birch trees presence in cells containing one or more artificial area (97% of the χ 2 observed value) are overrepresented as well (observed χ 2 = 485, critical χ 2 = 3.8, p-value <0.0001) and can be interpreted in the same way, both as the result of sampling bias, since observations are more easily made close to an inhabited area, but also because the inhabited areas are located in birch habitat.

On other hand, the χ 2 value is less important for birch observations in grid cells where there is a protected area (81% of the observed χ 2 value, e.g. 236, critical χ 2 =3.8, p-value <0.0001). The effect of sampling bias, with an observer preference for protected areas, is not reinforced by a preferential birch habitat effect, unlike the two previous cases.

) test between number of mountain Birch records and vegetal communities categories

Table IV displays main results of χ 2 tested conducted between ecological habitat of the Birch trees records and the Swedish vegetation map, at the micro-regional scale of our studied area. For 18 ecological habitats, mainly translated and summarized as peaty soils, covered with conifers, heathland, deciduous forests and/or mosses, Birch trees records are not present. In 21 habitats, Birch tree is observed in less than 4% of the polygons, corresponding to high altitude habitats, be it peaty soils, heathlands and conifer forests.

The three ecological habitats in which Birch trees has been more frequently recorded correspond to "Moss rich broadleaf forest" (60 observations in 6% of the polygons), and in a lesser extent "Deciduous forests on moist to hydromorphic soils" and "Deciduous forest on lichen-covered ground" (5 to 6 observations, in 2.3 to 2.8% of the total sum of polygons).

The χ 2 test is very significant (observed value 383.1, critical value 5.99, p<0.0001, Table IV), the overrepresentation of Birch trees in Broadleaved forests matching up with 92% of the observed χ 2 value.

b) Study case of Caluna AB data provider at a local scale with photointerpretation

As shown in the Figure 6, Birch records are concentrated on a narrow strip inferior to 100 m from the railway line (93.4 m on average), while photo-interpreted Birch stands are much more widespread. To the eastern part of the studied local site, the south shore of the line is very sloping, explaining the very low number of occurrences (N=2). The 11 other observations are made to the north, between the railway line and the E10 road, i.e. in an easily accessible area. Birch stands extend over more than one kilometer and a half (up to 1.7 km), but the numerous water bodies and wetlands (21% of the surface) probably constituted obstacles for observers. In the west, observations are made on the south side along a secondary road between a train station and the E10 highway, then observations are made on the north side when the south slope is rocky and partially snowy, without pedestrian access. The number of mountain Birch observations is not proportional to the area occupied by forest formations. Their density is 0.12 in the east and 0.29 in the west. This result strongly suggests that observations and records were more numerous while the presence of the species was more remarkable because of its scarcity within the landscape (e.g. observer-bias).

IV. DISCUSSION AND PERSPECTIVES a) OpenAccess biodiversity databases (OBDs) spatial variation in regard of geographical laws

By previous results, a discussion about the relevance of the general laws in geography, especially in biogeography and spatial ecology has been raised. Indeed, the principle of geographical similarity (as defined as similar configurations of two locations having similar features) has been well illustrated when looking at records distribution of Birch tree even before habitat niche modeling. Mainly, this study was dedicated to carefully analyze metadata and their ecological meaning applied to habitat niche modeling (e.g. abiotic factors describing and explaining records distribution). Geographical similarity has been evoked all along the study, in particular at the micro-regional and local scale of spatial analysis. Such similarity was geographicaldependent, mainly calling for altitudinal and topography factors, where mountain Birch woodlands oftenly are found in the shallows of glacial valleys, all along the shore of lakes, open waters and watercourses simultaneously with human facilities such as roads, railways and housings. Finally, law of similarity appears to be very useful to distinguish autocorrelation effects, either be natural or human-induced through sampling biases. Geographical laws allow to disentangle the role of sampling bias, expressed by autocorrelation law, and the role of environmental conditions, corresponding to the geographical similarity principle.

b) Applicability of the study for habitat niche modeling and perspectives

The main descriptor being latitude and longitude of mountain Birch records, analyzing links between institution providers, given by the GBIF dataset has been highly relevant. Spatial autocorrelation having compared different providing institutions has been promising in the frame of the aim: at the macro-regional scale, the distribution of Birch tree records for the Norrbotten county (Artdatabanken as providing institution) appears to be more suitable for habitat niche modeling. To support those first conclusions, and considering the dominant administrative and institutional (public) source of records, administrative boundaries of Swedish Lapland (from province, to counties until municipalities) would be interesting to select a dataset where records are not spatially clustered, at least not because of human factors (e.g. accessibility, facility).

c) Learning for Geography

From this experience, we can preconize a reasoned use of Open Access databases for modeling in geography, regardless of the type of data. First of all, it is important to keep in mind geographic laws as a theoretical framework, even if they are not understood or known by the producers of the data and the scientists using them. It is sometimes the case for OBDs, only the law of spatial auto-correlation has been identified as an obstacle (Silero and Barbosa, 2020), whereas the principle of geographic similarity was of great importance and formed the basis for reliable modeling. Procedures for correcting sampling bias must likewise be reexamined in the light of an analysis that balances the difference between real bias and a location corresponding to a certain geographical reality, according to the principle of geographical similarity. For instance, caution should be exercised if it is necessary to select or remove self-correlated observations based on geographic objects, when these geographic objects meet the environmental conditions for the presence of the observations.

  Sport and leisure facilities Thanks to the SLU and Lantmäteriet, Vegetation fjällkedjan vektor (shapefile of mountainous vegetal communities) and Vegetation Norrbotten vektor (shapefile of boreal vegetal communities of the Norrbotten province) were downloaded (https://www.lantmateriet.se/) and merged as geoprocessing/ Concerning Natural protected areas, Nationally Designated Areas (CDDA) have been downloaded from the European Environmental Agency Open Access data portal (https://dd.eionet.europa.eu/datasets/latest/CDDA) to get conservation areas of Swedish Lapland.

Fig. 3 .

 3 Fig. 3. Cartographic representation of filtered records of Birch tree and tested accessibility/facility factors (Sources: GBIF, Corine Land Cover 2018 Copernicus, European Environmental Agency and DivaGIS; realization R.Courault 2020)

Fig. 4 .

 4 Fig. 4. Distribution of Birch tree records according to the most contributive institutions of the GBIF dataset (Sources: GBIF; realization R.Courault 2020) Others data providers (e.g. institutions) presenting a low number of records are mainly represented by consult companies, NGO's and one municipality (e.g. ECOSIA, MFU, NA, NTNU -VM, Piteå kommun; Table I, figure 4)

Fig. 5 .

 5 Fig. 5. Plots of Birch trees records counted by cell of 1 km*1 km at the macro-regional scale and number of records of presences of accessibility/facility factors ("Sum_artif_CLC18" : Number of presences of artificial land covers and facilities; "Sum_bet_pub" : Number of presences of mountain Birch; "Sum_nat_areas" : Number of presences of natural protected areas; "Sum_pres_roads" : Number of presences of roads and railways ; realization R.Courault 2020

Fig. 6 .

 6 Fig. 6. Local photointerpretation displaying relationship between Birch records, transportation and Birch stands (Sources: GBIF; realization M.Cohen 2020)

TABLE I

 I

	.	NUMBER OF RECORDS BEFORE SEMI-AUTOMATIZED
	SPATIAL, BIOLOGICAL AND TEMPORAL FILTERING OF BIRCH TREES, GIVEN BY
	INSTITUTION PROVIDER AT THE MACRO-REGIONAL SCALE (SWEDISH
	LAPLAND)
	Institution source/provider	Σ records before semi-automatized spatial and biological filtering (1995-2020) Σ records %
	Artdatabanken		9831	80.1
	Calluna AB		50	0.4
	Ecofact		3	0.02
	Gävleborgs
	Botaniska		986	8
	Sällskap
	Jämtlands
	Botanika		25	0.2
	Sällskap
	Länsstyrelsen Jämtland		404	3.3
	Länsstyrelsen Norrbotten		15	0.12
	Länsstyrelsen Västerbotten		102	0.83
	Länsstyrelsen Västernorrland		81	0.7
	MFU		2	0.02
	NA		2	0.02
	NTNU-VM		7	0.06
	NV		743	6.1
	O		11	0.09
	Piteå kommun		3	0.02
	TOTAL		12265

TABLE II

 II

	.	MORANS' I SPATIAL AUTOCORRELATION TESTS FOR
	INSTITUTIONAL PROVIDERS SUB-SAMPLES OF BIRCH TREES AT THE MACRO-
	REGIONAL SCALE (SWEDISH LAPLAND)
			Geo-statistical parameters and results
	Providing institution		Σ columns	Σ rows	Σ records of mountain Birch	Moran's I index	z-score	p-value
	Artdatabanken	565	866	6504	+0.15	212 <0.0001
	Gävleborgs Botaniska sällskap	136	92	32	+0.03	7.1 <0.0001
	Jämtlands Botaniska sällskap	111	76	25	-0.0002 -0.04 0.96
	Länsstyrelsen Jämtland		138	127	392	+0.09	23.9 <0.0001
	Länsstyrelsen Västerbotten		329	290	102	+0.03	35 <0.0001
	Länsstyrelsen Västernorrland	170	192	81	+0.05	21.9 <0.0001

TABLE III

 III ABSENCE/PRESENCE OF BIRCH TREES RECORDS BY CELL AND PRESENCE/ABSENCE OF ACCESSIBILITY FACTOR AT THE MACRO-REGIONAL SCALE (SWEDISH LAPLAND)

	.	PEARSON	CHI-SQUARED	TABLE	(Χ 2 )
	Mountain
	Birch
	absence/

BETWEEN

pre sence Accessibility factor absence/presence No roa d Roa d No protec ted area Protec ted area No artificial ized land cover Artificial ized land cover No Birch

		0.4 16	6.1	0.381	1.734	0.068	4.3
	Presence of Birch	46. 02	670. 23	42.3	192.2	7.502	473.15
	p-value	<0.0001	<0.0001	<0.0001
	Observed χ 2	721.06	236		485
	Critical χ 2	3.8		3.8		3.8
	B. Spatial and statistical analysis at the micro-regional
	scale
	a) Chi-2 (χ 2

TABLE IV

 IV

	.	CHI-SQUARED VALUES (Χ 2 ) BETWEEN ECOLOGICAL
	HABITATS CATEGORIES AND MOUNTAIN BIRCH RECORDS COUNTED AMONG
	THOSE AT THE MICRO-REGIONAL SCALE (NORTHERN LAPLAND)

Ecological habitat Χ2 Mountain Birch Χ2 Σ absences Σ presences Total Habitats without Mountain Birch tree

		0.1	14	14.06
	Habitats where Mountain Birch is lowly observed	0.12	16.7	16.73
	Deciduous forests habitat	2.55	349.77	352.31
	Total	2.77	380.33	383.1

V. ACKNOWLEDGMENTS

We want to acknowledge deeply the TRAC3S project (Sorbonne-Université, Institut de la Transition Environnementale) and its carriers (Marie Alexis, Maryse Rouelle-Castrec) for their financial and scientific support in initiating this study, and we hope it will be useful for future projects. We also want to warmly thank Cecilia Svantesson for her very accurate translation of the Swedish Vegetation Maps categories from Swedish to French. Also, the GDR AREES 2012 (Groupe de Recherche Arctique, Environnement et Sociétés, CNRS Meudon, LGP) for its financial support in the translation of the Swedish Map of Vegetation. All our acknowledgements as well for the Abisko Research Station and the SLU (INTERACT network) without fieldworks and GIS datasets won't be available. Finally, we want to greet Saami people and the Sametinget/Swedish Saami Parliement for the time spent and sharing.