Simon D Carrière
email: simon.carriere@paca.inra.fr

Konstantinos Chalikakis

Charles Danquigny

Hendrik Davi

Naomi Mazzilli

Chloé Ollivier

Christophe Emblanch

The role of porous matrix in water flow regulation within a karst unsaturated zone: an integrated hydrogeophysical approach

Keywords: Karst, France, Hydrogeophysics, Carbonate rocks, Matrix porosity

Some portions of the porous rock matrix in the karst unsaturated zone (UZ) can contain large volumes of water and play a major role in water flow regulation. The essential results are presented of a local-scale study conducted in 2011 and 2012 above the Low Noise Underground Laboratory (LSBB -Laboratoire Souterrain à Bas Bruit) at Rustrel, southeastern France. Previous research revealed the geological structure and water-related features of the study site and illustrated the feasibility of specific hydrogeophysical measurements. In this study, the focus is on hydrodynamics at the seasonal and event time scales. Magnetic resonance sounding (MRS) measured a high water content (more than 10 %) in a large volume of rock. This large volume of water cannot be stored in fractures and conduits within the UZ. MRS was also used to measure the seasonal variation of water stored in the karst UZ. A process-based model was developed to simulate the effect of vegetation on groundwater recharge dynamics. In addition, electrical resistivity tomography (ERT) monitoring was used to assess preferential water pathways during a rain event. This study demonstrates the major influence of water flow within the porous rock matrix on the UZ hydrogeological functioning at both the local (LSBB) and regional (Fontaine de Vaucluse) scales. By taking into account the role of the porous matrix in water flow regulation, these findings may significantly improve karst groundwater hydrodynamic modelling, exploitation, and sustainable management.

Introduction

Karstified rocks cover a large portion of the world's surface [START_REF] Gunn | Encyclopedia of caves and karst science[END_REF], particularly around the Mediterranean Sea [START_REF] Bakalowicz | Ressources en eau du karst : un enjeu pour le bassin méditerranéen. [Karst water resources: an issue for the Mediterranean basin[END_REF]. Unfortunately the complexity of karst hydrosystems continues to impede sustainable water exploitation and management (e.g. [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF][START_REF] Bakalowicz | La zone d'infiltration des aquifères karstiques: méthodes d'étude, structure et fonctionnement [Infiltration zones in karst aquifers: methods of study -structure and functioning[END_REF][START_REF] Ford | Karst hydrogeology and geomorphology[END_REF][START_REF] Goldscheider | Methods in karst hydrogeology[END_REF][START_REF] White | A brief history of karst hydrogeology: contributions of the NSS[END_REF]). In the context of climate change, additional stress on water resources may require more intense exploitation of karst hydrosystem water resources. These factors and sustainable management require improved knowledge about how karst hydrosystems function. Karst is a complex medium with multi-scale heterogeneity. Water flow pathways are present throughout the entire medium, from rock matrix to fractures and karst features (Fig. 1). The hydrodynamic role of these water pathways is still poorly understood (e.g. [START_REF] Bailly-Comte | Water exchange and pressure transfer between conduits and matrix and their influence on hydrodynamics of two karst aquifers with sinking streams[END_REF] and water flow regulation within karst hydrosystems remains an important issue. On one hand, karst hydrosystems have high permeability because a large amount of water moves quickly through karst conduits and rapid pressure transfer occurs. On the other hand, karst hydrosystems provide an important buffer effect because spring discharge remains high even during long dry periods, as at the Fontaine de Vaucluse spring in France. Several hypotheses have been advanced to explain the important water flow regulation capacity of karst hydrosystems: (i) significant water storage within the saturated zone in karst conduits (e.g. [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF][START_REF] Marsaud | Structure et fonctionnement de la zone noyée des karsts à partir des résultats expérimentaux[END_REF]; (ii) water storage in the epikarst, (fractured and weathered near-surface karst) (e.g. [START_REF] Aquilina | Water storage and transfer in the epikarst of karstic systems during high flow periods[END_REF][START_REF] Ford | Karst hydrogeology and geomorphology[END_REF]; and (iii) delayed infiltration or travel time within the entire unsaturated zone (UZ) (e.g. [START_REF] Emblanch | Contribution of time tracers (Mg2+, TOC, delta C-13(TDIC), NO3-) to understand the role of the unsaturated zone: A case study -Karst aquifers in the Doubs valley, eastern France[END_REF][START_REF] Emblanch | Carbon 13 of TDIC to quantify the role of the unsaturated zone: the example of the Vaucluse karst systems (Southeastern France)[END_REF][START_REF] Mudarra | Hydrogeological functioning of a karst aquifer deduced from hydrochemical components and natural organic tracers present in spring waters. The case of Yedra spring (southern Spain)[END_REF][START_REF] Mudarra | Monitoring groundwater in the discharge area of a complex karst aquifer to assess the role of the saturated and unsaturated zones[END_REF]). However, this regulation capacity within the UZ is not related to identified geological features (e.g. fractures, matrix). Finally, few researchers consider water storage in the UZ zone to be significant enough to play a role in water flow regulation and the capacitive function of karst.

Due to these complexities, distributed hydrogeological modeling of karst systems remains difficult to implement and is not often done (e.g. [START_REF] Kiraly | Modelling karst aquifers by the combined discrete channel and continuum approach[END_REF][START_REF] Larocque | Determining karst transmissivities with inverse modeling and an equivalent porous media[END_REF][START_REF] Scanlon | Can we simulate regional groundwater flow in a karst system using equivalent porous media models? Case study, Barton Springs Edwards aquifer, USA[END_REF][START_REF] Worthington | Diagnostic hydrogeologic characteristics of a karst aquifer (Kentucky, USA)[END_REF][START_REF] Worthington | Self-Organized Permeability in Carbonate Aquifers[END_REF]. The usual approach to karst hydrosystems is the socalled "black box" model (e.g. [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF][START_REF] Marsaud | Structure et fonctionnement de la zone noyée des karsts à partir des résultats expérimentaux[END_REF]Labat et al. 2000a;Labat et al. 2000b;[START_REF] Rimmer | Modelling precipitation-streamflow processes in karst basin: The case of the Jordan River sources, Israel[END_REF][START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to Fontaine de Vaucluse (South of France)[END_REF][START_REF] Moussu | A multi-objective calibration framework for rainfall-discharge models applied to karst systems[END_REF][START_REF] Hartmann | A multi-model approach for improved simulations of future water availability at a large Eastern Mediterranean karst spring[END_REF]). However, "black box" models are site-specific. The lack of physics and geology in such models makes difficult their transfer from one karst system to another. Moreover, "black-box" modeling (e.g. [START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to Fontaine de Vaucluse (South of France)[END_REF]) and hydrochemical estimates (e.g. [START_REF] Batiot | Total Organic Carbon (TOC) and magnesium (Mg2+): two complementary tracers of residence time in karstic systems[END_REF] have indicated large residence time differences within a single karst hydrosystem.

Most techniques traditionally used in hydrogeology are of limited success in a complex and heterogeneous media such as karst [START_REF] Bakalowicz | Karst groundwater: a challenge for new resources[END_REF]. Hydrogeologists often apply techniques commonly used in surface hydrology, such as natural and artificial tracers or rainfall/runoff models.

In recent years, surface based geophysics was added to the methodological suite to improve the analysis of spatial and temporal variability of underground properties [START_REF] Berkowitz | Characterizing flow and transport in fractured geological media: A review[END_REF]. Numerous techniques are available, each with its strengths and weaknesses. [START_REF] Chalikakis | Contribution of geophysical methods to karst-system exploration: an overview[END_REF] proposed a general overview of geophysical methods for karst media and more recently [START_REF] Kaufmann | Detection and mapping of ghost-rock features in the Tournaisis area through geophysical methods -an overview[END_REF] published an overview on the use of geophysical surveys to detect ghost rock. In recent years, an increasing number of hydrogeophysical projects aimed at studies of karst hydrogeological functioning have appeared (e.g. [START_REF] Jacob | Absolute gravity monitoring of water storage variation in a karst aquifer on the larzac plateau (Southern France)[END_REF][START_REF] Worthington | Diagnostic hydrogeologic characteristics of a karst aquifer (Kentucky, USA)[END_REF]2010;[START_REF] Gondwe | Hydrogeology of the south-eastern Yucatan Peninsula: New insights from water level measurements, geochemistry, geophysics and remote sensing[END_REF][START_REF] Zhu | Challenges of using electrical resistivity method to locate karst conduits-A field case in the Inner Bluegrass Region, Kentucky[END_REF][START_REF] Deville | Caractérisation de la zone non saturée des karsts par la gravimétrie et l'hydrogéologie[END_REF][START_REF] Mazzilli | On the inclusion of ground-based gravity measurements to the calibration process of a global rainfall-discharge reservoir model: case of the Durzon karst system (Larzac, southern France)[END_REF].

This paper presents an integrated hydrogeophysical approach based on several ground-based geophysical methods combined with geological and hydrogeological techniques. This integrated approach provides additional insight into karst UZ structure and functioning. It also proposes a conceptual hydrogeological model to explain the multi-annual dynamics of the studied karst hydrosystem. This conceptual model will likely promote additional and more accurate hydrodynamic modeling. The methodological approach was developed and tested in a typical Mediterranean karst hydrosystem observatory: The Fontaine de Vaucluse -Low Noise Underground Laboratory (LSBB) karst watershed.

Experimental site

General geological and hydrogeological context

The LSBB, a nearly horizontal underground passageway originally dug for military purposes, was converted to a research laboratory in 1997 (Fig. 1 and2c). It is situated in the Fontaine de Vaucluse hydrosystem karst UZ, in the village of Rustrel. The rock cover over the passageway ranges in thickness from 0 to 519 m depending on the topography. The general saturated zone of the karst system is approximately 400 m beneath the LSBB. The 3.8 km long passageway traverses the karst medium and arbitrarily intersects faults and karst networks. As a result, the passageway also intersects some flow paths in the UZ. Since 2003 more than 61 flow points have been identified and studied within the laboratory. Three of the flow points are permanent and 58 are temporary; they are located between approximate depths of 33 m and 440 m [START_REF] Garry | Etude des processus d'écoulements de la zone non saturée pour la modelisation des acquifères karstiques. Expérimentation hydrodynamique et hydrochimique sur les sites du Laboratoire Souterrain à Bas Bruit (LSBB) de Rustrel et de Fontaine de Vaucluse[END_REF][START_REF] Blondel | Spatial and temporal tracing of groundwater in karst hydrosystems using dissolved organic matter. Experiments in and application to sites within the Underground Low Noise Laboratory (LSBB) of Rustrel[END_REF][START_REF] Barbel-Perineau | Caractérisation du fonctionnement de la zone non saturée des aquifères karstiques -Approche directe par études hydrodynamiques et hydrochimiques sur le bassin de recherche, d'expérimentation et d'observation de fontaine de Vaucluse -Laboratoire Souterrain à Bas Bruit de Rustrel -Pays d'Apt[END_REF].

One permanent flow point is located 33 m below the surface at the western extremity of the LSBB passageway (Fig. 2c). This point, called "point D," has an average discharge of approximately 130 mL/min [START_REF] Perineau | Hydrodynamic organisation of the flows in the unsatured zone of the Fontaine de Vaucluse karst system[END_REF]. The flow at point D is of great hydrogeological interest because it occurs in a featureless zone in a karst environment: no major fault, no apparent karstification (in both surface -exokarst -and at depth); furthermore, it is located several tens of meters below the so-called "epikarst" zone (according to Mangin's schema [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF]). However, it presents surprisingly smooth water dynamics compared to other flow points in the LSBB. The flow at D has sometimes increased after a rain event; however at other times flow has remained unchanged after a similar amount of rain [START_REF] Barbel-Perineau | Caractérisation du fonctionnement de la zone non saturée des aquifères karstiques -Approche directe par études hydrodynamiques et hydrochimiques sur le bassin de recherche, d'expérimentation et d'observation de fontaine de Vaucluse -Laboratoire Souterrain à Bas Bruit de Rustrel -Pays d'Apt[END_REF]. The geological context associated with this puzzling hydrogeological behavior led to a research focus on the area located above point D.

The study site is located within the Fontaine de Vaucluse karst hydrosystem in southeastern France.

The Fontaine de Vaucluse is the largest karst spring in Europe; between 1877 and 2004 it produced an average daily outlet discharge of 19 m 3 /s [START_REF] Cognard-Plancq | Historical monthly rainfall-runoff database on Fontaine de Vaucluse karst system: review and lessons[END_REF]. The catchment area is approximately 1115 km² (Fig. 2b) and the karst UZ is particularly thick; its average thickness is approximately 800 m [START_REF] Puig | Le système karstique de la Fontaine de Vaucluse[END_REF]. The large size is due to the presence of a nearly 1500 m thick massive and continuous limestone [START_REF] Masse | Contribution a l'étude de l'Urgonien (Barrémien -Bédoulien) des Monts de vaucluse et du Luberon[END_REF][START_REF] Masse | Les calcaires urgoniens de Provence ; Valanginien -Aptien inférieur ; Tome 1 : Stratigraphie -Paléontologie ; Tome 2 : Les paléoenvironnements et leur évolution[END_REF] comprising Necomanian marls to upper Aptian marls. A part of this carbonate platform is composed of reef limestone, which may reach a thickness of approximately 450 m in the region (Fig. 2d). This reef limestone contains Urgonian facies that are traditionally divided into three sub-divisions: U1, U2, and U3 [START_REF] Leenhardt | Etude géologique de la région du Mont Ventoux[END_REF]. A part of this Urgonian limestone has exceptionally high total matrix porosity for a limestone and it covers at least half of the Fontaine de Vaucluse catchment (e.g. [START_REF] Masse | Contribution a l'étude de l'Urgonien (Barrémien -Bédoulien) des Monts de vaucluse et du Luberon[END_REF][START_REF] Léonide | Diagenetic patterns and pore space distribution along a platform to outer-shelf transect (Urgonian limestone, Barremian-Aptian, SE France)[END_REF]. However, the role of this porosity in water dynamics is still unknown.

Local geological structure investigated in previous studies

The geological structure of the study site was investigated using geophysical and geological surveys [START_REF] Carrière | Combining Electrical Resistivity Tomography and Ground Penetrating Radar to study geological structuring of karst Unsaturated Zone[END_REF]) and geological modeling [START_REF] Ollivier | Improving a 3D geological model using multiple geophysical methods[END_REF]. To detect geological structure, a combination of ground penetrating radar (GPR) and electrical resistivity tomography (ERT) were used to identify geological features that may impact groundwater dynamics. The GPR results provide near surface high-resolution imaging and thus can provide relevant geological information such as stratification and fractures (Fig. 3b andc). Despite the exceptional quality of the results, GPR's investigational depth remains limited to around 12 m. ERT is able to investigate down to 40 m, but it is an integrative technique that has lower resolution than GPR. In the study area, the investigated limestone is commonly an electrically resistive formation (more than 2000 Ω.m). However, below a depth of 5-7 m, the ERT profiles reveal several zones of moderate resistivity (around 900 Ω.m). In these zones, cross-stratifications were clearly identified by GPR. The combination of both GPR and ERT results, in addition to field geological studies and geotechnical observations (recorded during underground passageway construction), lead to a well-founded geological interpretation (Fig. 3).

Cross-stratifications formed under certain sedimentation conditions (high energy, for example) may have produced a more porous limestone. However, during the first phase of the study, it was impossible to determine the exact role of the identified geological features. What are the respective roles of the detected faults, fractures, or crushed zones? Is it possible to monitor water recharge dynamics in the moderate resistivity zone? These questions were answered during the second phase of the study presented in this paper. [START_REF] Ollivier | Improving a 3D geological model using multiple geophysical methods[END_REF] performed geological modeling (Fig. 3b) using GOCAD software to combine all available geophysical, geological, and geotechnical information. This geological modeling made it possible to propose a well-founded interpretation for geological structures and to extend the analysis of the site below the depth investigated by geophysics.

An integrated hydrogeophysical approach: methods and tools

The approach used in this study combines several tools (surface based geophysics, geology and tectonics, geotechnical information, hydrodynamics, and infiltration modeling) to explain the functioning of the local hydrosystem.

Effective infiltration assessment

To assess as accurately as possible the effective infiltration entering locally into the karst hydrosystem, exchanges between soil/plant/atmosphere were modeled using CASTANEA [START_REF] Davi | Modeling carbon and water cycles in a beech forest: Part II: Validation of the main processes from organ to stand scale[END_REF][START_REF] Dufrêne | Modelling carbon and water cycles in a beech forest: Part I: Model description and uncertainty analysis on modelled NEE[END_REF]. This model calculates the effective infiltration by taking into account the following processes: canopy water interception, tree transpiration, soil and litter evaporation, dynamics of soil water content and drainage. This model is based on the Penmann-Monteith equation for calculating evapo-transpiration. The model includes the downregulation of soil drought on transpiration via stomata closure. CASTANEA is a species-specific model that considers features of the forest environment in the study area. CASTANEA was validated in a similar environment (Mediterranean shrubby forest with holm oak) at the Puéchabon site (southern France) using eddycovariance measurements [START_REF] Davi | Sensitivity of water and carbon fluxes to climate changes from 1960 to 2100 in European forest ecosystems[END_REF][START_REF] Martin | ilex L.) à la sécheresse à différentes échelles temporelles : Incidences sur la modélisation des processus[END_REF]. A requisite inventory of forest species was conducted for a representative zone measuring 120 m by 20 m to define stand-specific input parameters; it included a tree inventory to assess the standing biomass, hemispherical photographs to assess the canopy Leaf Area Index [START_REF] Davi | Effect of thinning on LAI variance in heterogeneous forests[END_REF], and soil pits to estimate soil water content.

Meteorological data were recorded at the Centre d'Information Régional Agro-Météorologique et Economique (CIRAME) station at Saint Saturnin-les-Apt (Fig. 2). Between 2004 and 2012, average daily temperatures ranged between -7.2 and 28.7°C. The average annual temperature was 13.5°C

and the average annual rainfall was 660 mm.

Flow measurement

The discharge rate at point D has been monitored since 2003 to study the response of the UZ to natural rainfall events [START_REF] Garry | Etude des processus d'écoulements de la zone non saturée pour la modelisation des acquifères karstiques. Expérimentation hydrodynamique et hydrochimique sur les sites du Laboratoire Souterrain à Bas Bruit (LSBB) de Rustrel et de Fontaine de Vaucluse[END_REF][START_REF] Blondel | Spatial and temporal tracing of groundwater in karst hydrosystems using dissolved organic matter. Experiments in and application to sites within the Underground Low Noise Laboratory (LSBB) of Rustrel[END_REF][START_REF] Barbel-Perineau | Caractérisation du fonctionnement de la zone non saturée des aquifères karstiques -Approche directe par études hydrodynamiques et hydrochimiques sur le bassin de recherche, d'expérimentation et d'observation de fontaine de Vaucluse -Laboratoire Souterrain à Bas Bruit de Rustrel -Pays d'Apt[END_REF]. Flow point D shows smooth hydrodynamic behavior even though it is located at a depth of 33 m (Fig. 4 and5).

Within the LSBB passageway, 42 other flow points have also been monitored since 2003. In this paper, we compare point D to point C; these two points show the maximum difference in terms of hydrodynamics. The measurement frequency was almost weekly (349 measurement campaigns at the end of 2012). Discharges are measured by collecting water in a graduated cylinder over a period of time.

Geophysics implementation, acquisition strategy, and field constraints

After extensive testing with a large array of geophysical methods (electromagnetic, electric, gravimetric, and nuclear magnetic resonance), [START_REF] Carrière | Etude hydrogéophysique de la structure et du fonctionnement de la zone non saturée du karst[END_REF], two efficient and accurate techniques were selected to study hydrodynamics at the local scale (few hundred meter squares):

ERT and magnetic resonance sounding (MRS). These geophysical techniques make it possible to image temporal variation in complementary physical properties of the subsurface (until 90 m). This variation may be related to variation in water content and thus may shed light on hydrodynamic behavior.

Measurements at the study site were conducted at two time scales: (i) with a seasonal time step to detect rough variations, and (ii) with a short time step during a large rain event to detect fast water flows.

At the experimental site (Fig. 4), the slope, the vegetation, and the gravel cover required extensive preparation prior to measurements. The exact position of all geophysical measurements was measured with the real time kinematic (RTK) acquisition method every 2 m using the differential global positioning system (GPS) TRIMBLE 5800.

MRS

MRS is the only ground-based geophysical method that can obtain a signal directly related to the presence of water to a depth of several meters. MRS, which is designed for a quantitative description of aquifer formations, is based on the phenomenon of hydrogen proton magnetic resonance. The MRS array generates a signal with a specific frequency to excite the hydrogen protons in water molecules, and the MRS signal received is specifically related to groundwater (e.g. [START_REF] Legchenko | Nuclear magnetic resonance as a geophysical tool for hydrogeologists[END_REF][START_REF] Vouillamoz | Localization of saturated karst aquifer with magnetic resonance sounding and resistivity imagery[END_REF]. sounding was chosen because it is located above the moderate resistivity zone detected with ERT.

The acquisition system used is a Numis Plus by IRIS instruments. The measurements were collected in summer (July 2011) and in winter (January 2012). To integrate a large area above water point D, a wide MRS array (80x80 m) was chosen for these surveys. It is important to mention that for both MRS surveys, electromagnetic noise was low compared to signal amplitude (ratio signal to noise >1.5) which ensure interpretability and comparability between measurements (e.g. [START_REF] Legchenko | Nuclear magnetic resonance as a geophysical tool for hydrogeologists[END_REF].

The details of MRS results from this site were discussed and presented by [START_REF] Mazzilli | MRS applicability for the characterisation of the unsatured zone of karst hydrosystems. First results from the LSBB site[END_REF] and [START_REF] Chalikakis | Water in karst hydrosystems unsaturated zone; MRS evidences within an integrated hydrogeophysical approach[END_REF]. These results indicated the presence of a significant amount of water in the UZ karst to a depth of 90 m (MRS sometimes indicated a water content greater than 10 %). Eight MRS were conducted at the experimental site to study the spatial variability of MRS signals. All soundings exhibit similar forms. However, significant lateral variation in water content was observed between soundings [START_REF] Mazzilli | MRS applicability for the characterisation of the unsatured zone of karst hydrosystems. First results from the LSBB site[END_REF], revealing that the karst UZ is also a heterogeneous environment in terms of water content.

ERT

The ERT technique has been widely used in karst areas because it is robust and reliable (e.g. [START_REF] Cardarelli | Electrical resistivity tomography to detect buried cavities in Rome: a case study[END_REF][START_REF] Robert | A salt tracer test monitored with surface ERT to detect preferential flow and transport paths in fractured/karstified limestones[END_REF]. The sensitivity of electrical resistivity to moisture variations and the rapidity of ERT measurements (approximately one hour per section) led to the selection of this geophysical technique to image rapid resistivity variation during heavy rainfall. During the 30 day campaign, ERT time-lapse acquisition ranged from every three hours during the rain event (17 days)

to one section a day after the rain event.

The acquisition system used is an ABEM Terrameter SAS 4000 [START_REF] Dahlin | The development of DC resistivity imaging techniques[END_REF] with four channels and 64 electrodes. Implantation of the ERT electrodes, mainly at limestone outcrops, required holes that were mechanically dug into the rock; saltwater and mud were added to ensure a good quality ground contact. The ERT section LSBB03 acquired with Gradient array is presented in this paper (Fig. 4). The section is oriented east-west. This direction is perpendicular to the general slope, is subperpendicular to one of the main fault and lineament directions, and is the most heterogeneous direction in terms of apparent resistivity spatial distribution [START_REF] Carrière | Combining Electrical Resistivity Tomography and Ground Penetrating Radar to study geological structuring of karst Unsaturated Zone[END_REF]). The Gradient array was chosen for this survey because of its robustness and its rapidity [START_REF] Dahlin | A numerical comparison of 2D resistivity imaging with 10 electrode arrays[END_REF].

A portion of the first ERT results were presented by [START_REF] Carrière | Feasibility and Limits of Electrical Resistivity Tomography to Monitor Water Infiltration Through Karst Medium During a Rainy Event[END_REF] to evaluate the effectiveness and technical limits of ERT to monitor water infiltration via previously recognized karst features under natural conditions.

Apparent (directly measured) resistivity (ρ α ) analysis is not usually considered in ERT surveys.

However, [START_REF] Carrière | Feasibility and Limits of Electrical Resistivity Tomography to Monitor Water Infiltration Through Karst Medium During a Rainy Event[END_REF] demonstrate that with the current technology, in this kind of complicated media, existing inversion schemes are not adequate. It was hoped that raw results with no artifacts related to inversion schemes could provide relevant information. If one analyzes variation in resistivity of inverted sections, it is only possible to observe general evolution of the near surface. Due to its integrative character, the inversion process smooths any fine variation even if one tries to observe only inverted resistivity variation between two consecutive time steps. In this context, apparent resistivity variations (Δρ α ) between two consecutive time steps were analyzed.

These variations were normalized by the delay (Δ T ) between consecutive measurements (ρ n and ρ n-1 ) using the following equation (Eq.1).

∆𝜌 𝛼 = ( 𝜌 𝑛-1 𝜌 𝑛 -1) 100 ∆ 𝑇 (1)

Results and interpretation

The hydrodynamics of flow points C and D as a function of recharge (rain and effective infiltration) is presented first, followed by the results of geophysical monitoring at both seasonal and event time scales.

Flow dynamics versus rain and effective infiltration

Two hydrological periods can be observed: a dry period from 2004 to 2007 and a rainy period from 2008 to 2013 (Fig. 5) (Barbel- [START_REF] Barbel-Perineau | Caractérisation du fonctionnement de la zone non saturée des aquifères karstiques -Approche directe par études hydrodynamiques et hydrochimiques sur le bassin de recherche, d'expérimentation et d'observation de fontaine de Vaucluse -Laboratoire Souterrain à Bas Bruit de Rustrel -Pays d'Apt[END_REF]. Effective infiltration modeling using CASTANEA (Fig. 5) confirms this alternation of wet and dry periods. One notes that effective infiltration is quite limited between 2004 and 2007. Rigorous modeling of effective infiltration is essential in supporting subsequent interpretations because the difference between rain and effective infiltration can be quite large and can lead to misinterpretation.

During the first dry period, the water discharge at point D generally decreases and there is no clear reaction during rain events. On the contrary, during the second rainy period, this water discharge increases and responses to rain events are more apparent. The regularity of point D is obvious when comparing variation ranges between points D and C.

Point C discharge is highly reactive to rainfall; it shows typical karstic dynamics even though it is located at a depth of 256 m (Fig. 5). Note that the reactivity of flow point D to a rain event differs from dry period (no variation) to wet period (rapid discharge increases for strong rainy events).

However, the flow rate at point C increases proportionally with effective infiltration in both dry and wet periods. Without more arguments, it is difficult to explain the difference between D and C. One can only point out that point D is located within porous Urgonian limestone and C is located within lower Barremian limestone (Fig. 2d).

Hydrogeophysical monitoring

Seasonal dynamics

MRS was used to detect the seasonal variation in storage in the karst UZ between summer 2011 and winter 2012. In processed MRS results (Fig. 3d), two layers reveal unusual water contents for karst UZ. The first interval with high water content (10 % water content according to MRS) is at a depth of approximately 15 m (above low point D); the other high water content interval (12 % water content according to MRS) is deeper than 40 m.

Although these elevated water contents are surprising for a karst hydrosystem UZ, these values agree with the high porosity estimated by other geophysical methods used at the LSBB [START_REF] Maufroy | Caractérisation et modélisation numérique de l'effet de site topographique 3D : application à la Grande Montagne de Rustrel, Vaucluse[END_REF][START_REF] Bereš | Caractérisation de l'anisotropie d'une plate-forme carbonatée karstifiée: approche méthodologique conjointe sismique et électrique[END_REF]) and total matrix porosity measured in plugs in the region (e.g. [START_REF] Masse | Contribution a l'étude de l'Urgonien (Barrémien -Bédoulien) des Monts de vaucluse et du Luberon[END_REF][START_REF] Léonide | Diagenetic patterns and pore space distribution along a platform to outer-shelf transect (Urgonian limestone, Barremian-Aptian, SE France)[END_REF]) (Fig. 2d).

Significant temporal variations in MRS water content are observed at several depth levels (Fig. 3d).

The measured variations are as high as 3 % water content. This result indicates that these storage variations could participate in seasonal water flow regulation in the karst UZ.

It is difficult to conceive that such a large volume of water could be stored in fracture and conduits porosity of karst UZ. For this reason, in this geological context where limestone can be very porous, the majority of this water content can be attributed to the porous matrix.

The high water content level detected by MRS near the depth of 15 m corresponds to location of the cross-stratifications and a moderate resistivity zone detected by GPR and ERT (Fig. 3b andc). This feature probably contains a stock of water located above point D, and supplies this water to point D.

Its hydrogeological role will be discussed in the next section.

Event dynamics

This unique large-scale ERT surface based experiment was conducted during a typical Mediterranean autumn rain event (17 days). A total of 230 mm of rain were recorded and 120 ERT time-lapse sections were measured over the same profile (LSBB03) during and after the rain event (a total of 30 days).

Apparent resistivity values have been arithmetically averaged for each section. Apparent resistivity decreased sharply during the rain event, from 1750 to 1050 Ω.m (Fig. 6b). These variations do not seem to be related to temperature variations because air temperature and water temperature (at point D) remained stable during the experiment [START_REF] Carrière | Etude hydrogéophysique de la structure et du fonctionnement de la zone non saturée du karst[END_REF]). Thus, these electrical resistivity variations can be reasonably related to variation in water content in the subsurface. Analysis of this mean apparent resistivity indicator made it possible to select the twelve critical time steps presented in Fig. 6a: before, during, and after the rain event. These results are presented (Fig. 6a) as hourly percentage changes in ρ α with a basic representation of vegetation and soil cover around the profile.

Previously recognized karst features [START_REF] Carrière | Combining Electrical Resistivity Tomography and Ground Penetrating Radar to study geological structuring of karst Unsaturated Zone[END_REF]) are also represented.

At the beginning of the rain event (Fig. 6a, section 12), ρ α decreased moderately and homogeneously along the section. This resistivity decrease could be related to moistening of near surface horizons. After the first rain event (Fig. 6a, section 23), ρ α stabilized quickly. During the following heavy rain episode (Fig. 6a, sections 45 to 78), the moistening process appeared to be quite heterogeneous and some zones resembled water pathways or preferential infiltration zones.

However, it would be unwise to link observed ρ α variations at depth with deep moistening processes.

The observed variations can also be directly influenced by near surface variations. Immediately after the rain (Fig. 6a, sections 910 and 1011), ρ α increased in some zones. This result may be related to drainage processes. These zones correlate well with zones previously identified as preferential pathways. This second observation reinforced the hypothesis of preferential pathways playing a hydraulic role that can be traced with ERT monitoring. However, it is still not possible to specify properties such as the geometry of these pathways. Other zones where such drainage processes are not identifiable may be related to zones where soil is thicker and remains moist after rainfall.

The assumed water pathways identified on the basis of apparent resistivity do not seem to be related to vegetation density or soil cover variations (Fig. 6). However, these zones appear related to fractures or fault zones previously detected by GPR [START_REF] Carrière | Combining Electrical Resistivity Tomography and Ground Penetrating Radar to study geological structuring of karst Unsaturated Zone[END_REF]. Note also that resistivity variations in fractures are different from those measured in the crushed zone located in the western part of the profile. A crushed zone may provide a larger water pathway than a fracture. It is important to emphasize that the methods used in this study do not image the pathways with certainty. Each brown point represents one electrical resistivity tomography section [START_REF] Carrière | Feasibility and Limits of Electrical Resistivity Tomography to Monitor Water Infiltration Through Karst Medium During a Rainy Event[END_REF]. PsZ is pseudo depth.

Discussion

Local scale: study area

This integrated hydrogeophysical approach at two different time scales has improved the current understanding knowledge of the hydrogeological functioning of study site. At the seasonal time step, MRS made it possible to identify water storage variation above and below flow point D. These variations are related to slow water dynamics within matrix porosity. This hypothesis is consistent with the regular flow dynamics at point D, which were observed during both dry and recharge periods (Fig. 5). Moreover, it is consistent with the large total matrix porosity within Urgonian limestone, which may reach 20 % (Fig. 2d).

At the rain event time scale, ERT monitoring revealed heterogeneity in infiltration (Fig. 6). Some fractures or faulted zones previously detected by GPR (Fig. 3) were identified as probable preferential water pathways (Fig. 7b andd).

During wet periods, when water fills these vertical structures, hydraulic head (Δh) increases above flow point D, which induces a prompt increase in flow at D as in October 2011 (Fig. 7d) during ERT monitoring (Fig. 6). During this type of rain event, this increased flow may push forward older, less mobile water according to the water schema proposed by [START_REF] Barbel-Perineau | Caractérisation du fonctionnement de la zone non saturée des aquifères karstiques -Approche directe par études hydrodynamiques et hydrochimiques sur le bassin de recherche, d'expérimentation et d'observation de fontaine de Vaucluse -Laboratoire Souterrain à Bas Bruit de Rustrel -Pays d'Apt[END_REF] who studied hydrochemistry within the LSBB.

During long dry periods, as for example, between 2004 and 2007 (Fig. 7a), flow point D did not respond to rain pulses even during extreme events (Fig. 7b). Thus, there is probably no significant hydraulic continuity between near-surface water-filled fractures and the deeper porous matrix reservoir. Thus, there is no transmission of a pressure pulse to cause increased flow at point D. This absence of connection is probably due to partial drying of the porous matrix reservoir above point D.

The evolution of the dry porous matrix to water-filled matrix may take several years. This timeframe can explain the multi-annual inertia to recharge the stock that supplies point D, as was observed between 2007 and 2010 (Fig. 5). It is particularly relevant to study the relationship between rainfall and water discharge during the period from 2004 to 2012 because recharge was highly variable during this period. A long dry period between 2004 and 2007 was followed by a more rainy period between 2008 and 2012. The annual dynamics of rainfall, effective infiltration, and water discharge presented in Fig. 8 illustrate these dry and wet periods.

Although these results are applicable to a discussion of the relationship between rainfall and output flow at LSBB flow points, the relationship with spring discharge at the Fontaine de Vaucluse should be approached with caution. However, the total annual rainfall in Saint Saturnin-les-Apt (Fig. 2b) is representative enough of inter-annual rainfall variation in the region to allow the following discussion.

These inter-annual recharge variations cause changes in discharge rates at the outlet of the karst system as a whole (Fontaine de Vaucluse spring) and at local outlets of karst sub-systems such as water points C and D in the LSBB. However, these variations are highly variable between the Fontaine de Vaucluse and water points C and D, and are summarized as follows:

-Water point C exhibits typical karst dynamics with large flow range (a factor of 10) depending on the quantity of annual effective infiltration. This relationship was clearly demonstrated in 2007 when effective infiltration was at a minimum and water flow at C was very low.

-The water point D has different dynamics because its water discharge is quite regular; the range of its annual discharge variation is only a factor of 2. The flow rate at point D decreased continuously between 2004 and 2007, after which its flow increased regularly until 2010. However, effective infiltration decreased significantly between 2008 and 2010. Thus, there is a kind of delayed dynamic. The results highlight the importance of porous matrix on water flow, even in a karst limestone aquifer system.

-The annual water discharge of the Fontaine de Vaucluse exhibits dynamics similar to those at point D (Fig. 8) although its flow variation range is slightly larger than point D, a factor of 3.5.

Similar to point D, the flow rate at the Fontaine de Vaucluse also increased between 2007 and 2010; however effective infiltration decreased significantly from 2008 to 2010. This observation suggests that water flow regulation mechanisms observed at point D are site-specific. Similar mechanisms of water flow regulation through the porous matrix may play an important role in the entire Fontaine de Vaucluse hydrosystem (Fig. 8). Considering the large volume of Urgonian porous formations in the hydrosystem (Fig. 2b), this assumption is a plausible explanation of part of the multi-annual dynamics of the Fontaine de Vaucluse spring.

Even so, at a daily time step, the Fontaine de Vaucluse spring can be highly reactive, similar to flow

Conclusions

This study was conducted at a local scale around the LSBB underground laboratory to improve the current understanding of the local karst hydrosystem. The results presented in this paper demonstrate that a considerable volume of water can be stored within matrix rock in the karst UZ. [START_REF] Masse | Contribution a l'étude de l'Urgonien (Barrémien -Bédoulien) des Monts de vaucluse et du Luberon[END_REF] measured exceptional matrix porosity for limestone (up to 20 %) on the Urgonian platform around Fontaine de Vaucluse catchment area. However, the hydrodynamic role of this matrix porosity was not yet known.

MRS surveys confirmed high free water content (5 to 10 %) within the karst UZ. This high water content indicates a large volume of water that cannot be stored within fractures or karst conduits in the UZ. This water is stored within limestone layers that were also identified by previous ERT and GPR surveys. Moreover, MRS surveys identified seasonal variations in stored water (around 3 %), which indicates that this water reserve is active. The ERT monitoring at the event scale made it possible to demonstrate that despite the presence of porous limestone, the study site has characteristics of a fractured/karstified area in which recharge processes are highly heterogeneous.

Rapid water circulation can occur through preferential water pathways and cause flow peaks at point D.

Without additional measurements, the geophysical survey cannot be used to extend the hydrodynamic properties observed at the LSBB site to the Fontaine de Vaucluse hydrosystem, either wholly or in part. That is why it is not reasonable to try to quantify the consequences of the presence of porous limestone on the Fontaine de Vaucluse spring flows. However, considering that Urgonian limestone covers more than half the surface area of the Fontaine de Vaucluse basin (Fig. 2b) and that [START_REF] Masse | Contribution a l'étude de l'Urgonien (Barrémien -Bédoulien) des Monts de vaucluse et du Luberon[END_REF] measured high matrix porosity on the entire Urgonian platform, there are necessarily repercussions in Fontaine de Vaucluse outlet flow that may explain the multi-annual dynamics of water flow (Fig. 8).

If Urgonian limestones play an important role in Fontaine de Vaucluse hydrodynamics, similar geological formations can similarly impact the hydrodynamics of other karst hydrosystems, including areas far from the Mediterranean Sea where the sedimentary conditions are similar. For example, Cenozoic carbonate rocks in the Gulf of Mexico and the Caribbean are highly porous [START_REF] Halley | High-porosity Cenozoic carbonate rocks of south Florida: progressive loss of porosity with depth[END_REF]. Therefore, on the basis of these results presented, it is essential to consider the role of the porous matrix in water flow regulation within a karst UZ. Hydrodynamic modelling should be adjusted accordingly to assess properly the groundwater resource potential and improve sustainable management and exploitation of karst hydrosystems.

Fig. 1 :

 1 Fig. 1: Karst system model showing the Fontaine de Vaucluse hydrosystem and the Low Noise Underground Laboratory of Rustrel (LSBB).

Fig. 2 :

 2 Fig. 2: (a) The Fontaine de Vaucluse basin located in France; (b) The Rustrel experimental site located in the Fontaine de Vaucluse basin; (c) Excerpt of local geological map, n°942 (after Blanc et al. 1973); (d) Total porosity of limestone outcropping on Fontaine de Vaucluse hydrosystem measured in a small sample (after Guglielmi 2010).

Fig. 3 :

 3 Fig. 3: Comparison between different types of geological and geophysical information: (a) Regional lithostratigraphic log (Masse and Fenerci-Masse 2011 (modified)); (b) 3D geological model (Ollivier et al. 2013) with 252 m long electrical resistivity tomography (ERT) section LSBB03, inverted in three iterations; (c) Zoom on ERT section LSBB03 with bedding and fractures detected by ground penetrating radar (GPR); (d) Magnetic resonance soundings (MRS) of July 2011 and January 2012 (Carrière, 2014).

Fig. 4 :

 4 Fig. 4: Experimental site of the Low Noise Underground Laboratory (LSBB): (a) location of geophysical measurements on aerial photo; (b) pseudo 3D location of geophysical measurements discussed in this paper.

Fig. 5 :

 5 Fig. 5: Precipitation and effective infiltration calculated using CASTANEA compared with flow dynamics at points C and D between 2004 and 2012.

Fig. 6 :

 6 Fig. 6: (a) Hourly change in apparent resistivity between two consecutive time steps. Positions of fractures previously detected by ground penetrating radar (GPR) and basic representation of vegetation and soil cover; (b) Evolution of mean apparent resistivity during monitoring versus rain.Each brown point represents one electrical resistivity tomography section[START_REF] Carrière | Feasibility and Limits of Electrical Resistivity Tomography to Monitor Water Infiltration Through Karst Medium During a Rainy Event[END_REF]. PsZ is pseudo depth.

  Fig. 7a-d: Synthetic functioning model of the local hydrosystem that supplies flow point D. D flow variation depends on rain conditions.

  point C. However, even if floods at the Fontaine de Vaucluse can be impressive, they represent only a limited part of the annual water discharge. The annual flow dynamics of the Fontaine de Vaucluse seem to be predominantly influenced by slow water flow within a porous matrix, as is the case at point D.

Fig. 8 :

 8 Fig. 8: Comparison of annual water discharge for Fontaine de Vaucluse (FdV) and flow points D and C with annual precipitation and effective infiltration modeled for the period of 2004 to 2012 by CASTANEA.

Acknowledgments

The authors would like to express their gratitude to CIRAME, to ORE H+ and to all the LSBB team for their technical and logistic help. This study was funded by a French ministry of education and research PhD grant. This work was performed within the framework of the FDV/LSBB observation site, part of the KARST observatory network (www.sokarst.org) initiative of INSU/CNRS, which seeks to support knowledge sharing and promote cross-disciplinary research on karst systems. We also thank the editorial board of Hydrogeology Journal, S. White, M. Saribudak and an anonymous reviewer for helping us improve this paper.


Mohamad Harastani

Oscar S Carlos

Sorzano

Dr Slavica Jonić
email: slavica.jonic@upmc.fr

Hybrid Electron Microscopy Normal Mode Analysis (HEMNMA) with Scipion

Keywords: cryo-electron microscopy, single particle analysis, continuous conformational changes, normal mode analysis, structure, dynamics, software 50-75

Hybrid Electron Microscopy Normal Mode Analysis (HEMNMA) method was introduced in 2014. HEMNMA computes normal modes of a reference model (an atomic structure or an electron microscopy map) of a molecular complex and uses this model and its normal modes to analyze single-particle images of the complex to obtain information on its continuous conformational changes, by determining the full distribution of conformational variability from the images. An advantage of HEMNMA is a simultaneous determination of all parameters of each image (particle conformation, orientation, and shift) through their iterative optimization, which allows applications of HEMNMA even when the effects of conformational changes dominate those of orientational changes. HEMNMA was first implemented in Xmipp and was using MATLAB for statistical analysis of obtained conformational distributions and for fitting of underlying trajectories of conformational changes. A HEMNMA implementation independent of MATLAB is now available as part of a plugin of Scipion V2.0 (http://scipion.i2pc.es). This plugin, named ContinuousFlex, can be installed by following the instructions at https://pypi.org/project/scipion-em-continuousflex.

In this article, we present this new HEMNMA software, which is user-friendly, totally free, and open-source.

INTRODUCTION

Cryo-electron microscopy (cryo-EM) has become comparable to X-ray crystallography with regards to the obtainable resolution of structures of biomolecular complexes, which are now increasingly determined at near-atomic resolution. [START_REF] Zhang | Near-atomic resolution using electron cryomicroscopy and single-particle reconstruction[END_REF][START_REF] Yu | 3.88 a structure of cytoplasmic polyhedrosis virus by cryoelectron microscopy[END_REF][START_REF] Liao | Structure of the trpv1 ion channel determined by electron cryo-microscopy[END_REF][START_REF] Bai | An atomic structure of human gamma-secretase[END_REF][START_REF] Bartesaghi | 2.2 a resolution cryo-em structure of beta-galactosidase in complex with a cell-permeant inhibitor[END_REF][START_REF] Nguyen | Cryo-em structure of the yeast u4/u6.U5 tri-snrnp at 3.7 a resolution[END_REF][START_REF] Natchiar | Visualization of chemical modifications in the human 80s ribosome structure[END_REF][START_REF] Zhang | Structure of phycobilisome from the red alga griffithsia pacifica[END_REF][START_REF] Liu | Near-atomic cryo-em imaging of a small protein displayed on a designed scaffolding system[END_REF][START_REF] Fan | Single particle cryo-em reconstruction of 52 kda streptavidin at 3.2 angstrom resolution[END_REF] No requirement for sample crystallization and the possibility to elucidate multiple conformations of a complex from the same sample are among the main advantages of cryo-EM. Characterizing the different conformations that can coexist is essential for understanding how the complexes function and addressing their dynamics.

To achieve near-atomic resolution of 3D reconstructions, the classical approach is to collect a large number of images of complexes (particles) at random and unknown orientations within a thin layer of vitreous ice, then, perform 2D and 3D classifications into an initially set number of classes and, finally, perform 3D reconstruction using only those particles that have the most consistent views and conformations (those that contribute to the highest-resolution class averages) while removing all other particles. [START_REF] Zhang | Structure of phycobilisome from the red alga griffithsia pacifica[END_REF][START_REF] Liu | Near-atomic cryo-em imaging of a small protein displayed on a designed scaffolding system[END_REF][START_REF] Fan | Single particle cryo-em reconstruction of 52 kda streptavidin at 3.2 angstrom resolution[END_REF][START_REF] Khatter | Structure of the human 80s ribosome[END_REF][START_REF] He | Near-atomic resolution visualization of human transcription promoter opening[END_REF][START_REF] Banerjee | 2.3 a resolution cryo-em structure of human p97 and mechanism of allosteric inhibition[END_REF][START_REF] Abeyrathne | Ensemble cryo-em uncovers inchworm-like translocation of a viral ires through the ribosome[END_REF] Such "selection" of particles may obscure information on a possibly larger conformational variability as some conformational states may be thrown away blindly instead of being elucidated. Therefore, non-classification-based methods are required to extract from images the full distribution of conformational variability (the so-called conformational space or landscape in which the images are mapped) and to assemble 3D reconstructions from identified, more or less dense regions in this space (denser regions contain more frequent conformational states and less dense regions contain less frequent states). Such methods are necessary for studying continuous conformational changes of complexes (the concept including the possibility of unequally distributed conformers, e.g. possible existence of more stable conformers) and are referred to as continuous-state methods.

Many classification-based methods (referred to as discrete-state methods) can be found in the literature. [START_REF] Van Heel | Characteristic views of e. Coli and b. Stearothermophilus 30s ribosomal subunits in the electron microscope[END_REF][START_REF] Penczek | A method of focused classification, based on the bootstrap 3d variance analysis, and its application to ef-g-dependent translocation[END_REF][START_REF] Scheres | Disentangling conformational states of macromolecules in 3d-em through likelihood optimization[END_REF][START_REF] Elad | Detection and separation of heterogeneity in molecular complexes by statistical analysis of their two-dimensional projections[END_REF][START_REF] Simonetti | Structure of the 30s translation initiation complex[END_REF][START_REF] Sorzano | A clustering approach to multireference alignment of single-particle projections in electron microscopy[END_REF][START_REF] Penczek | Identifying conformational states of macromolecules by eigen-analysis of resampled cryo-em images[END_REF][START_REF] Scheres | A bayesian view on cryo-em structure determination[END_REF][START_REF] Lyumkis | Likelihood-based classification of cryo-em images using frealign[END_REF][START_REF] Fischer | Structure of the e. Coli ribosome-ef-tu complex at <3 a resolution by cs-corrected cryo-em[END_REF][START_REF] Scheres | Chapter six -processing of structurally heterogeneous cryo-em data in relion[END_REF][START_REF] Punjani | Cryosparc: Algorithms for rapid unsupervised cryo-em structure determination[END_REF][START_REF] Haselbach | Structure and conformational dynamics of the human spliceosomal b(act) complex[END_REF][START_REF] Grant | Cistem, user-friendly software for single-particle image processing[END_REF] While the majority of the discrete-state methods require that the number of classes is set initially, setting this number in some of the methods is less arbitrary as it is based on making a balance between intraclass and interclass variances after a statistical analysis of the data, typically using an eigenvector analysis of the covariance estimated from images or from volumes reconstructed from random subsets of images. 18; 19; 21 The covariance eigenvectors estimation has also been introduced in the context of continuous conformational heterogeneity. 29; 30 The latter methods are closely related to continuous-state methods whose development started recently and is currently an active field of research. [START_REF] Jin | Iterative elastic 3d-to-2d alignment method using normal modes for studying structural dynamics of large macromolecular complexes[END_REF][START_REF] Sorzano | Hybrid electron microscopy normal mode analysis graphical interface and protocol[END_REF][START_REF] Dashti | Trajectories of the ribosome as a brownian nanomachine[END_REF][START_REF] Frank | Continuous changes in structure mapped by manifold embedding of single-particle data in cryo-em[END_REF][START_REF] Anden | Structural variability from noisy tomographic projections[END_REF][START_REF] Dashti | Ourmazd A Conformational dynamics and energy landscapes of ligand binding in ryr1[END_REF][START_REF] Lederman | Singer A Continuously heterogeneous hyper-objects in cryo-em and 3-d movies of many temporal dimensions[END_REF] We introduced the continuous-state method referred to as Hybrid Electron Microscopy Normal Mode Analysis (HEMNMA) in 2014. 31; 32 HEMNMA computes normal modes of a reference model (an atomic structure or an EM map) of a molecular complex and uses this model and its normal modes to analyze single-particle images of the complex iteratively in order to extract information on the full distribution of conformational variability from these images. In each iteration, HEMNMA aims at simultaneously solving both orientational and conformational heterogeneity of images using the reference model and its normal modes.

Although often considered as "model-free", the continuous-state methods from other groups do use a model (e.g., the model used by Dashti et al (2014) [START_REF] Dashti | Trajectories of the ribosome as a brownian nanomachine[END_REF] is EMDB-1067 density map). This model is used to first solve the orientational heterogeneity of images, by determining the particle 3D orientation and 2D shift in each image assuming conformational homogeneity of images (e.g., the particle orientation and shift in each image can be determined by the standard projection matching of the image with the model). Then, the conformational heterogeneity of images is solved (i.e., the conformations are determined) using the image orientations and shifts determined in the previous step. These methods do not refine orientations and shifts when determining conformations, assuming that the effects of orientational changes dominate those of conformational changes so that each image can be associated with the correct orientation.

An advantage of HEMNMA is a simultaneous determination of all parameters of each image (particle conformation, orientation, and shift) through their iterative optimization, which allows applications of HEMNMA even when the effects of conformational changes dominate those of orientational changes. The simultaneous determination of the parameters is done by iterative matching of the image with projections of the reference model being deformed using a linear combination of normal modes (the conformation is estimated by determining the unknown coefficients of the linear combination, i.e. by determining the displacement amplitudes along normal modes). The reference model does not need to be an atomicresolution model and can also be an EM map. This EM map can be obtained by 3D reconstruction from combined data as if the data was conformationally homogeneous (e.g., by the so-called random sample consensus approach that identifies the best EM maps from those reconstructed using random subsets of images and random image orientations, [START_REF] Vargas | Efficient initial volume determination from electron microscopy images of single particles[END_REF] maximum likelihood optimization of a number of 3D reconstructions starting from random data subsets, [START_REF] Scheres | Relion: Implementation of a bayesian approach to cryo-em structure determination[END_REF] or other methods [START_REF] Reboul | Single-particle cryo-em-improved ab initio 3d reconstruction with simple/prime[END_REF] ).

HEMNMA and other continuous-state methods are fundamentally different from classification-based (discrete-state) methods. For instance, the discrete-state method of Haselbach et al (2018) [START_REF] Haselbach | Structure and conformational dynamics of the human spliceosomal b(act) complex[END_REF] uses standard Relion 39 software for discrete classification of images into a number of 3D maps and standard principal component analysis to cluster the maps. As similar classification-based methods, it is prone to assigning images to lower-noise density maps independently of whether or not the images belong there. For further reading on the different methods, the reader can refer to some of the recent methods reviews. [START_REF] Klaholz | Structure sorting of multiple macromolecular states in heterogeneous cryo-em samples by 3d multivariate statistical analysis[END_REF][START_REF] Jonic | Computational methods for analyzing conformational variability of macromolecular complexes from cryo-electron microscopy images[END_REF][START_REF] Cossio | Likelihood-based structural analysis of electron microscopy images[END_REF][START_REF] Sorzano | Survey of the analysis of continuous conformational variability of biological macromolecules by electron microscopy[END_REF][START_REF] Serna | Hands on methods for high resolution cryo-electron microscopy structures of heterogeneous macromolecular complexes[END_REF] The first version of HEMNMA was implemented in Xmipp [START_REF] Sorzano | Xmipp: A new generation of an open-source image processing package for electron microscopy[END_REF] and was using MATLAB for statistical analysis of obtained conformational distributions and for fitting of underlying trajectories of conformational changes. 31; 32 A MATLAB-independent version of HEMNMA was recently implemented in Scipion [START_REF] De La Rosa-Trevin | Scipion: A software framework toward integration, reproducibility and validation in 3d electron microscopy[END_REF] and is currently available as part of a plugin of Scipion V2.0 (http://scipion.i2pc.es). This plugin, named ContinuousFlex, can be installed by following the instructions at https://pypi.org/project/scipion-em-continuousflex. The new HEMNMA software is totally free and open-source, and the new graphical interface is even more user-friendly than the previous one. In this article, we present this new HEMNMA software.

HEMNMA WITHIN SCIPION: STEPS, PARAMETERS, GRAPHICAL INTERFACE

The steps to follow are listed in the HEMNMA menu on the left side of the Scipion project window and are numbered from 1 to 6. When the different steps are executed, a tree-like structure of the project appears on the right side of the window, with each block corresponding to a step or a substep (Figs. 12). In this section, we describe the steps (Steps 1-6), parameters, and graphical interface of HEMNMA. Each parameter is also described in the graphical interface, via a help message that can be displayed by clicking on the question mark next to the corresponding parameter field. The majority of the parameters are set to the values that usually produce good results. These default values are visible in the graphical interface and can be modified by the user. The parameters whose values are expected to be less frequently changed are hidden by default. By selecting "Advanced" as "Expert Level" mode ("Normal" is the default mode), the values of these "advanced" parameters can be visualized (the parameter names showing up on a grey background) and modified (e.g., Fig 3A).

Step 1. Reference model: The reference model to import can be a PDB file with atomic coordinates or a density volume (e.g., an EM map or a simulated map from a PDB structure).

Both import options are provided (steps 1.a and 1.b1 in the project tree, Figs. 12). Fig. 1 and Fig. 2 show the project tree for an input atomic structure and for an input density volume, respectively. Before the next step, the input density volume must be converted into a PDBformat file. We provide a tool for the volume conversion into a set of 3D Gaussian functions (step 1.b2 in the project tree, Fig. 2), the so-called pseudoatoms, [START_REF] Jonic | Coarse-graining of volumes for modeling of structure and dynamics in electron microscopy: Algorithm to automatically control accuracy of approximation[END_REF] whose coordinates are written into a PDB file. More precisely, the conversion is performed so that the input density volume is represented with Gaussian functions of a given (input) standard deviation ("Pseudoatom radius" input parameter expressed in voxels), by minimizing the normalized volume approximation error [START_REF] Jonic | Coarse-graining of volumes for modeling of structure and dynamics in electron microscopy: Algorithm to automatically control accuracy of approximation[END_REF] to reach a given (input) value of this error ("Volume approximation error" parameter expressed in percent) (Fig. 3). Note that the default value of the "Volume approximation error" parameter is 5 %, which can be modified by first selecting "Advanced" as "Expert Level" mode to visualize the parameter (Fig. 3A). Through this conversion procedure, some noise may be removed from the input density volume, as shown elsewhere. [START_REF] Jonic | Denoising of high-resolution single-particle electron-microscopy density maps by their approximation using three-dimensional gaussian functions[END_REF] However, strong particle background noise should be removed using external tools or the masking procedures available in HEMNMA ("Mask mode" parameter in Fig. 3A allows entering the name of a binary mask file or the threshold value below which the densities should be removed). For faster processing, the volume-to-pseudoatoms conversion was parallelized and should be run using several CPU threads (sharing the same memory).

The number of threads can be specified via the "Threads" parameter (Fig. 3A). The results of the volume-to-pseudoatoms conversion can be visualized with Chimera [START_REF] Pettersen | Ucsf chimera--a visualization system for exploratory research and analysis[END_REF] (Fig. 3B) by clicking on "Analyze Results" button in the project window (Fig. 2).

After a displacement with normal modes, both an atomic reference model and a pseudoatomic reference model (from a reference density volume) are converted into a density volume whose projections are compared with images to determine the conformations in these images (Step 5). The higher the resolution of the reference model, the more accurate the normal modes and the conformational determination will be. Thus, the most accurate conformational determination is expected with input atomic models and high-resolution EM maps.

The reference atomic model can be a full-atomic model or a coarse-grained model (e.g., containing Cα atoms only). The computation of normal modes (Step 2) requires information on the shape of the entire complex. An atomic model lacking a few small regions (e.g., disordered regions) may still be used as the reference model. If the atomic model lacks large portions of a complex and an EM map can be obtained for the entire complex, this EM map should be used as the reference model. An alternative reference model could be a hybrid model of the entire complex obtained by modeling the missing parts in the atomic model or the density volume computed from this hybrid model.

Classical, classification-based approaches can be used to obtain a reference EM map for

HEMNMA to analyze the entire heterogeneous set of images or an image subset with heterogeneity reduced through classification, as shown elsewhere [START_REF] Jin | Iterative elastic 3d-to-2d alignment method using normal modes for studying structural dynamics of large macromolecular complexes[END_REF] . The use of classificationbased approaches before HEMNMA may be particularly useful in the case of combined discrete and continuous conformational heterogeneity (e.g., a mixture of flexible complexes bound and unbound with ligands).

There are no restrictions regarding the number of atoms in input PDB files or the size of input EM maps, except that the maps should have a cubic shape. There are no restrictions regarding the content of the atomic model, such as the presence of DNA, ligands, etc. However, in the case of a heterogeneous set of particles bound and unbound with ligands, it may be interesting to remove the ligand from the reference model, as shown elsewhere. 31; 51 Step 2. Normal mode analysis (NMA): Normal modes of the atomic or pseudoatomic representation of the reference model ("Input structure" field in Fig. 4A) are computed using Tirion's elastic network model [START_REF] Tirion | Large amplitude elastic motions in proteins from a single-parameter, atomic analysis[END_REF] . In this model, atoms or pseudoatoms interact if they are connected with elastic springs. The interaction cut-off distance parameter determines how many atoms or pseudoatoms will be connected with elastic springs. The following two options exist for setting this value, via "Cut-off mode" parameter (Fig. 4A): "absolute" and "relative". If "absolute" is chosen for "Cut-off mode" (unshown in Fig. 4A), the interaction cut-off distance value is an input parameter that should be set directly by the user (by setting "Cut-off distance" parameter, expressed in angstroms, which shows up only if "absolute" is chosen for "Cut-off mode"). The default value of "Cut-off distance" is 8 Å, which is the empirical recommended value for atomic structures. The "Cut-off distance" value of 8 Å implies building an elastic network in which two atoms are connected with a spring (interact with each other) only if their distance is smaller than 8 Å. If "relative" is chosen for "Cut-off mode", the interaction cut-off distance value is computed automatically based on a given (input) percentage of distances ("Cut-off percentage" parameter expressed in percent, Fig. 4A) that should be below this value. The default value of "Cut-off percentage" is 95 %, which is the empirical recommended value for pseudoatomic structures. The "Cut-off percentage" value of 95 % implies building an elastic network in which two pseudoatoms are connected with a spring only if their distance is smaller than the one below which 95 % of distances are.

The larger the interaction cut-off distance for a given complex, the more rigid the elastic network will be. A too large interaction cut-off distance produces a too rigid elastic network, which can be detected and corrected by analyzing the normal-mode animations with the provided graphical interface (NMA results viewer in Fig. 4B).

Normal modes are computed by diagonalizing a square matrix of second derivatives of the potential energy function of the structure (Hessian matrix). One dimension of the Hessian matrix is 3 times the number of atoms or pseudoatoms. For faster Hessian diagonalization with input atomic structures (large Hessians due to large numbers of atoms), the structure is split into blocks of several residues, each block having 6 degrees of freedom (3 rotations and 3 translations), which reduces the Hessian dimension to 6 times the number of blocks and is known as Rotation Translation Block (RTB) method. [START_REF] Tama | Building-block approach for determining low-frequency normal modes of macromolecules[END_REF] The number of residues per RTB block should be set by the user, and its default value in the graphical interface is 10 (Fig. 4A). The larger the number of residues per block (i.e., the smaller the number of blocks), the faster the computation of normal modes and the more rigid the elastic network will be. A too large number of residues per block will produce a too rigid elastic network, which can be detected and corrected by analyzing the normal-mode animations with the provided graphical interface (NMA results viewer, Fig. 4B).

The number of computed normal modes is 3 times the number of pseudoatoms (for a pseudoatomic structure) or 6 times the number of RTB blocks (for an atomic structure).

However, the number of modes written on the disk and animated should be set by the user ("Number of modes" parameter, Fig. 4A). It is usually enough to save animations for 20-100 modes (e.g., the most relevant modes to experimentally observed conformational changes of low-symmetry structures are generally among the first 10-20 lowest-frequency non-rigidbody modes). The computed normal modes can be inspected by visualizing their animations, collectivities [START_REF] Bruschweiler | Collective protein dynamics and nuclear spin relaxation[END_REF] and scores [START_REF] Sorzano | Hybrid electron microscopy normal mode analysis graphical interface and protocol[END_REF] , via the NMA results viewer (Fig. 4B) that can be obtained by clicking on "Analyze Results" button in the project window. The collectivity [START_REF] Bruschweiler | Collective protein dynamics and nuclear spin relaxation[END_REF] indicates how much the atoms or pseudoatoms move together with a given mode. It is computed for the entire complex without differentiating between its domains and it is normalized between 1/N (maximally localized motion) and 1 (maximally collective motion), where N is the number of atoms or pseudoatoms. Highly collective low-frequency modes have been shown to be functionally relevant [START_REF] Delarue | On the use of low-frequency normal modes to enforce collective movements in refining macromolecular structural models[END_REF][START_REF] Suhre | Norma: A tool for flexible fitting of highresolution protein structures into low-resolution electron-microscopy-derived density maps[END_REF][START_REF] Tama | Flexible multi-scale fitting of atomic structures into low-resolution electron density maps with elastic network normal mode analysis[END_REF][START_REF] Wang | Global ribosome motions revealed with elastic network model[END_REF] . Therefore, we also provide the score [START_REF] Sorzano | Hybrid electron microscopy normal mode analysis graphical interface and protocol[END_REF] that combines the collectivity and frequency criteria, by penalizing the modes with higher frequencies and lower collectivities. The score is normalized between 1/M (highest collectivity and lowest frequency motion) and 1 (lowest collectivity and highest frequency motion), where M is the number of modes.

The list of modes, with their collectivities and scores, can be obtained using "Display output Normal Modes?" in the NMA results viewer (Fig. 4B). The interface allows ordering the modes according to the increasing or decreasing values of the collectivity or score measures.

The modes most relevant to the actual conformational change are expected to be among those with the lowest scores. The number of the lowest-score modes to inspect should be larger in the case of highly symmetrical structures (e.g., icosahedral-symmetry viruses) as highly collective modes may also exist at higher frequencies in such cases. To help the user decide which modes should certainly not be selected for the image analysis (Step 5), the modes with collectivities below 0.15 (very localized motions) are unchecked in the list of modes obtained via "Display output Normal Modes?" (Fig. 4B). The value of 0.15 for this collectivity threshold can be modified by the user ("Threshold on collectivity" parameter, Fig. 4A).

The motion trajectory along a normal mode is saved in a text file by concatenating the frames (in PDB format) of the coordinate displacement along this mode. The default number of frames is 10 and the default amplitude of the displacement along normal mode is 50 (these parameters can be modified via the "Animation" tab, unopened in Fig. 4A). This trajectory file is animated with VMD, [START_REF] Humphrey | Vmd: Visual molecular dynamics[END_REF] by selecting "Display mode animation with VMD?" and specifying the mode to be animated in the "Mode number" field (Fig. 4B). VMD provides tools to save the motion trajectories in Animated GIF or MPEG movie formats that can be played with more standard movie players.

The modes can also be inspected by plotting the shifts of atoms (or pseudoatoms) along the specified mode and by plotting their maximum shifts over all modes via "Plot mode distance profile?" and "Plot max distance profile?", respectively (Fig. 4B). For instance, a too large shift of one or a very few atoms (or pseudoatoms) with respect to the shift of other atoms (or other pseudoatoms) is typically a sign that the value of the atomic (or pseudoatomic) interaction cut-off distance used for computing normal modes is not optimal. In this case, normal modes should be recomputed using a modified value of this parameter.

Step 3. Information: At this step, entitled "Stop here or continue" in the HEMNMA menu (Figs. 12), we remind the user that HEMNMA software may also be used for NMA only (the processing can be stopped after computing and analyzing normal modes at Step 2). The following steps should be performed if aiming at analyzing conformational heterogeneity in images using the normal modes computed in Step 2.

Step 4. Images: This step allows importing images that will be analyzed with normal modes.

The images should have a square shape. They should have a power of 2 pixels in each dimension when using the so-called "wavelet & splines" method for rigid-body alignment in combination with the elastic alignment (see Step 5 for more information about the alignment methods). The larger the size of images, the longer the image analysis time will be. The image size can be reduced if the speed is an issue. The image size of 128x128 pixels is usually a good compromise.

Step 5. Conformational distribution: In this step, a combined elastic and rigid-body alignment is performed between particle images and the reference model to calculate the parameters of orientation, translation, and elastic deformation (amplitudes of normal modes) of the reference model that best describe the given particle image.

Fig. 5A shows the interface ("Input" tab) for selecting the images to be analyzed ("Input particles" parameter) and the normal modes to use for the image analysis ("Modes selection" parameter) from an entire list of computed normal modes ("Normal modes" parameter). The modes to use for the image analysis can be selected as the modes with the highest collectivities or least scores, with or without using previous knowledge about possible movements. The 6 lowest-frequency modes (modes 1-6, where the index of the mode corresponds to the mode order according to the increasing frequency) describe rigid-body movements of the structure and are often referred to as rigid-body normal modes. These modes should not be selected because rigid-body movements are taken into account through the iterative combined rigid-body and elastic alignment, as explained in the next paragraph.

Usually, this alignment is performed with 1-10 selected normal modes excluding modes 1-6.

The interface for setting the parameters of the combined elastic and rigid-body alignment is shown in Fig. 5B ("Combined elastic and rigid-body alignment" tab). We use Powell's optimization (more precisely, Powell's trust region method [START_REF] Powell | Uobyqa: Unconstrained optimization by quadratic approximation[END_REF] ) to iteratively estimate the elastic deformation (normal-mode amplitudes) of the reference model by maximizing the similarity between the particle image and the best matching projection of the deformed model (objective function, in optimization terminology). In each iteration of Powell's method, an estimate of the normal-mode amplitudes is used to displace the atoms or the pseudoatoms along normal modes and the obtained, elastically deformed structure is converted into a density volume that is then rigid-body aligned with the particle image. The rigid-body alignment allows the determination of the orientation and translation parameters (rigid-body parameters) and the objective-function evaluation that is then used to better estimate the normal-mode amplitudes (elastic deformation parameters) for the next iteration of Powell's method.

The values of the parameters of Powell's method are set internally in the code and are not modifiable via the interface except for the scaling factor of the initial trust-region radius ("Elastic-alignment trust region scale" parameter in Fig. 5B). The default value of the scaling factor is 1, which generally works well. When expecting larger conformational changes, it may be interesting to increase the scaling factor (typically to a value between 1 and 2). The rigid-body alignment can be performed using one of the following two methods: 1) "projection matching" that stands for the standard reference-library-based (discrete) projection matching in real space (faster but less accurate method); and 2) "wavelets & splines" that stands for a discrete projection matching in wavelet space [START_REF] Sorzano | A multiresolution approach to orientation assignment in 3d electron microscopy of single particles[END_REF] followed by a continuous matching in 3D Fourier space based on spline interpolation [START_REF] Jonic | Spline-based image-to-volume registration for three-dimensional electron microscopy[END_REF] (slower but more accurate method; in particular, more robust to noise). The default choice is "wavelets & splines" ("Rigid-body alignment method" in Fig. 5B). The angular sampling step [START_REF] Sorzano | A multiresolution approach to orientation assignment in 3d electron microscopy of single particles[END_REF] for the reference projection library computation is 10° by default and can be modified via the interface ("Discrete angular sampling" parameter expressed in degrees, Fig. 5B).

This combined rigid-body and elastic alignment is the task that requires the most computation. Therefore, it was MPI parallelized to simultaneously process N particle images using N MPI cores on CPU computers, clusters or supercomputers. The user specifies the number of MPI cores to be used ("MPI" parameter in Fig. 5A). On an Intel Xeon CPU at 2.9

GHz, the analysis of one particle image of size 128×128 pixels with 3, 6, and 9 normal modes may respectively take 2, 5, and 10 min when using "projection matching" rigid-body alignment and 5, 13, and 25 min when using "wavelets & splines" rigid-body alignment.

Normal-mode amplitudes resulting from the image analysis with normal modes can be visualized in 1 dimension (histogram) as well as in 2 or 3 dimensions (Fig. 5C). Points in a 2D and 3D space correspond to images with assigned orientations, translations, and normalmode displacement amplitudes with respect to the reference model. Close points in this space correspond to similar conformations and vice versa.

Step 6. Dimension reduction, clusters, and trajectories:

The normal-mode displacement amplitudes obtained by image analysis in Step 5 can be projected onto a space of lower dimension (1D, 2D, or 3D) (Fig. 6), using the technique of principal component analysis (PCA) or one of several other dimension reduction techniques. 63 To perform the dimension reduction, the user specifies the conformational distribution to be analyzed (the image analysis results obtained in Step 5), the desired lower dimension, and one of the available linear or non-linear dimensionality reduction techniques via "Conformational distribution", "Reduced dimension", and "Dimensionality reduction method" fields, respectively (Fig. 6A).

The "Extra params" field (Fig. 6A) allows modifying the parameters of the dimensionality reduction techniques (the parameters and their default values are listed in the help message linked to the "Dimensionality reduction method" field). The available dimensionality reduction techniques are PCA, Kernel PCA, Probabilistic PCA, Local Tangent Space Alignment (LTSA), Linear LTSA, Diffusion Map, Linearity Preserving Projection, Laplacian Eigenmap, Hessian Locally Linear Embedding, Stochastic Proximity Embedding, and Neighborhood Preserving Embedding. 63 The dimension reduction viewer (Fig. 6B) can be obtained by clicking on "Analyze Results" button in the project window. It allows displaying the axes of the low-dimensional space (Fig. 6C) specified via the "Display normal-mode amplitudes in the low-dimensional space" field (Fig. 6B). Also, the dimension reduction viewer (Fig. 6B) allows opening the tool for making animations ("Trajectories Tool", Fig. 7A) and the tool for computing 3D reconstructions ("Clustering Tool", Fig. 8A).

While the advantage of non-linear dimensionality reduction techniques is better suitability to non-linear manifold data representations, the advantage of linear dimensionality reduction techniques (PCA, Linear LTSA, Linearity Preserving Projection, Probabilistic PCA, and Neighborhood Preserving Embedding) is that they allow mapping from the low-dimensional space back to the original space. HEMNMA uses this property to generate animated trajectories of conformational changes. More precisely, "Trajectories Tool" (Fig. 7A) allows recording and visualizing the displacement of the reference model in the low-dimensional space. The specification of the displacement trajectory to be animated requires the user's interaction. The outlier data points can be removed by providing logical (Boolean) expressions ("Expression" field in Fig. 7A). The trajectory is specified by the coordinates of 10 points in the low-dimensional space (red points in Fig. 7B). The user may select all of these 10 points (by clicking on the plot to select each point) or 8 points may be automatically placed on a line between 2 points selected by the user (the first and last points of the trajectory). The position of the initially placed points may be changed by dragging the points.

As in the case of animations of normal modes in Step 2, the trajectory is saved in a text file by concatenating the PDB-format frames of the coordinate displacement. The trajectory can be animated with VMD (Fig. 7C) or saved with VMD in Animated GIF or MPEG movie formats for playing with other movie players.

Additionally, conformational changes can be analyzed in terms of 3D reconstructions from images in the low-dimensional space. "Clustering Tool" (Fig. 8A) allows making groups of close points (Fig. 8B), corresponding to images with similar conformational states, and computing 3D reconstructions from these groups (Fig. 8C). A group of points is specified by providing logical (Boolean) expressions ("Expression" field in Fig. 8A) or by clicking on the plot and dragging to add points to the group. Each point in the selected group of points is denoted by a yellow circle (Fig. 8B). Clicking on "Create Cluster" button (Fig. 8A) saves the selected group of points and performs 3D reconstruction from this group. The saved group can be inspected by displaying the 3D reconstruction results (slices and isosurface of the reconstructed volume, Fig. 8C) and the images present in the group (the images are not shown in Fig. 8). HEMNMA uses a fast Fourier-space method for 3D reconstruction from the selected group of points. The 3D reconstruction can be performed with other reconstruction methods (available in Scipion or other software packages), using the output metadata (text) file with the rigid-body and elastic alignment parameters corresponding to the selected group of points. The 3D reconstructions represent the average states from the corresponding groups of images. One should make as homogeneous groups of points as possible with a sufficient number of points per group in order to obtain high-resolution 3D reconstructions (note that the example in Fig. 8B does not show the most optimal grouping of points).

As explained, HEMNMA provides a full conformational variability landscape before grouping images into 3D reconstructions (images with similar conformations on the landscape). As such, HEMNMA can detect the conformations and motions that are undetectable with traditional classification-based approaches, which was extensively studied elsewhere. [START_REF] Jin | Iterative elastic 3d-to-2d alignment method using normal modes for studying structural dynamics of large macromolecular complexes[END_REF]

PERSPECTIVES

A full description of conformational heterogeneity is important for both biology and drug design. Cryo-EM has been under continuous development since its beginning (early '80s), which was even accelerated with the latest instrumental developments including direct electron detector devices (DDD cameras). Allowing a routine near-atomic resolution of structures and a routine full description of conformational variability are currently two main challenges. Indeed, the recent cryo-EM advances such as DDD cameras, phase plates, and sample motion correction have reduced noise and improved contrast in images, which makes elucidation of conformations from images more accurate. Methods for a full description of continuous conformational variability are being developed. However, these methods will need to be more efficient and user-friendly to allow routine use of cryo-EM for such studies.

HEMNMA is user-friendly software that has been developed for determining the full distribution of continuous conformational variability from cryo-EM images. Still, HEMNMA needs to be made faster in order to allow an efficient high-resolution description of this variability. For instance, the analysis of 10 4 particle images based on the elastic alignment with 6 normal modes and the "wavelets & splines" rigid-body alignment may take 17 h using 128 MPI cores on 2.9 GHz Intel Xeon CPUs, which means that the analysis of 10 6 particle images would require around 70 days of use of the same 128 MPI cores. To get a highresolution description of the full conformational landscape of a complex, several millions of particle images would need to be analyzed. A combination of HEMNMA with a deep learning approach will be implemented in the future to speed up processing of such large numbers of particle images. density volume (transparent grey) and its pseudoatomic representation (spheres where the amplitudes of Gaussian functions are color-coded from white (smallest) to red (largest)). In this example, the Gaussian-function standard deviation (pseudoatom radius) and the desired volume approximation error are 2 voxels and 5%, respectively. The desired volume approximation error parameter is hidden by default ("Expert Level" is set to "Normal") and can be visualized by setting "Expert Level" to "Advanced" (panel A). The volume-topseudoatoms conversion is parallelized using threads (parallel processes sharing the same memory) and the number of threads to use can be specified in the "Threads" field (panel A).

The masking options to remove background noise from the input volume are explained in the help message under the corresponding question mark (panel A). Cluster" button (panel A) appears as a new box in the project tree (panel C), possibly after using "Refresh" button from the project window. The saved group of points can be inspected using "Analyze Results" button from the project window, which shows images in this group (not shown here) and slices of the volume reconstructed from these images (panel C). The volume isosurface can be visualized by clicking on Chimera icon in the slices display menu (panel C).

  Ultramicroscopy 103:303-317. 63. Van der Maaten L, Postma E, Van den Herik J. 2009. Dimensionality reduction: A comparative review. Tilburg centre for Creative Computing, Tilburg university. No. 2009-005.

Figure 1 :

 1 Figure 1: HEMNMA project tree for the case of analyzing images using an atomic structure

Figure 2 :

 2 Figure 2: HEMNMA project tree for the case of analyzing images using a density volume as

Figure 3 :

 3 Figure 3: Volume-to-pseudoatoms conversion. (A) Dialog box. (B) Superposition of an input

Figure 4 :

 4 Figure 4: Normal mode analysis and visualization. (A) Dialog box (note that "Cut-off mode"

Figure 7 :

 7 Figure 7: Animations after image analysis. (A) Trajectories Tool. (B) Example of trajectory

Figure 8 :

 8 Figure 8: 3D reconstructions after image analysis. (A) Clustering Tool for grouping close

ACKNOWLEDGMENTS

We acknowledge the help of José Miguel de la Rosa-Trevín, CNB-CSIC (currently at SciLifeLab, Stockholm University) with porting the HEMNMA code from Xmipp to Scipion V1.0; the support of the French National Research Agency -ANR (ANR-19-CE11-0008-01 to SJ); the access to the HPC resources of CINES and IDRIS granted by GENCI (2019-A0070710998 to SJ); and the support of French Microscopy Society -Sfµ and Sorbonne University -SU (2019 Sfµ Master internship grant and 2019 SU PhD scholarship grant to MH).

as well). (B) Normal modes viewer allowing visualizing animations of normal modes (using "Display mode animation with VMD?" for the mode specified in "Mode number") and checking their collectivities and scores (using "Display output Normal Modes?"). analysis results in 3 dimensions (amplitudes along 3 normal modes). The image analysis task is parallelized using MPI protocol and the number of MPI cores (parallel processes that do not necessarily share the same memory) can be specified in the MPI field (panel A). In panel B, the default choice for "Rigid-body alignment method" is "wavelets & splines" (the alternative method is "projection matching") and the elastic alignment is performed by Powell's trust region method. In panel C, each point corresponds to a particle image with assigned orientations, translations, and normal-mode displacement amplitudes with respect to the reference model (close points correspond to similar conformations and vice versa).