Christophe Botella

Maximilien Servajean

Pierre Bonnet

Alexis Joly

Botella Christophe

Servajean Maximilien

Bonnet Pierre

Joly Alexis

Overview of GeoLifeCLEF 2019: plant species prediction using environment and animal occurrences

Keywords: LifeCLEF, biodiversity, environmental data, species recommendation, evaluation, benchmark, Species Distribution Models, methods comparison, presence-only data, model performance, prediction, predictive power

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Overview of GeoLifeCLEF 2019: plant species prediction using environment and animal occurrences 1 Introduction

The automatic prediction of the species most likely to be observed at a given location is an important issue for many areas such as biodiversity conservation, land management or environmental education. First, it could improve species identification processes and tools by reducing the list of candidate species observable at a given site (whether automated, semi-automatic or based on traditional field guides or flora). More generally, it could facilitate biodiversity inventories and compliance with regulatory obligations for the environmental integration of development projects. Finally, it could be used for educational purposes through biodiversity discovery applications offering functionalities such as contextualized educational pathways.

In the context of LifeCLEF evaluation campaign 2019 [START_REF] Joly | Lifeclef 2019: Biodiversity identification and prediction challenges[END_REF], the objective of the GeoLifeCLEF challenge is to evaluate the state of the art of species prediction methods over the long term and with a view to reproducibility. To achieve this, the challenge freely provides researchers with large-scale, documented and accessible data sets over the long term. Concretely, the aim of the challenge is to predict the list of species that are the most likely to be observed at a given location. Therefore, we provide a large training set of species occurrences and a set of environmental rasters that characterize the environment in a quantitative and qualitative way at any position in the territory. Indeed, it is usually not possible to learn a species distribution models directly from spatial positions because of the limited number of occurrences and the sampling bias. What is usually done in ecology is to predict the distribution of species based on a representation in environmental space, typically a characteristic vector composed of climatic variables (mean temperature at that location, precipitation, etc.) and other variables such as soil type, land cover, distance to water, etc. GeoLife-CLEF's originality is to encourage the extension of this approach to learning a more complex representation space that takes into account various input data such as environmental descriptors, their spatial structure and the known biotic context. Therefore, we provide tools to facilitate the extraction of environmental tensors that can be easily used as input data to models such as convolutional neural networks.

In 2019, the provided data was significantly enriched and several methodological improvements have been made. In more details, the new features introduced are as follows:

1. Pl@ntNet occurrences: to increase the amount of plant occurrences in the training set, we completed the publicly available data from the GBIF6 with user-generated observations of the Pl@ntNet mobile application [START_REF] Affouard | Pl@ ntnet app in the era of deep learning[END_REF]. These data are clearly noisier and more biased than conventional occurrence data but they can be filtered by the confidence level of the taxonomic automatic classifier used in the app and they have the advantage of being produced in huge quantities. 2. Occurrences of other kingdoms: to investigate how knowledge of the presence of non-plants organisms can help predict the presence of plants species, we provided a large training set of occurrences from other kingdoms coming from the GBIF platform. 3. A better quality test set: to ensure the reliability of our evaluation, the occurrence data of the test set were restricted to expert data with the highest species identification certainty and high geographical accuracy (lower than 50 m). Last but not least, the test occurrences were sampled in order to avoid, as much as possible, biases of spatial coverage and in the species representation. By this way, it contributes to give relatively more importance to rare species and scarce areas.

In the following sections, we describe in more details the data produced and the evaluation methodology used. We then present the results of the evaluation and the analysis of these results.

Dataset

Train occurrences

Pl@ntNet raw data. (PL complete) This data is directly pulled from [START_REF] Botella | Pl@ntnet queries 2017-2018 in france[END_REF]. Pl@ntNet7 is a smartphone app using machine learning to identify plant species from pictures submitted by a broad public of users. For each submission, also called a query, the Pl@ntNet algorithm answers a distribution of probability values across the targeted taxonomic referential. If the users allows it, the query's geolocation is also stored. In the provided training data, we used all accurately geolocated queries (with maximum 30 meters uncertainty) in France from the beginning of 2017 to the end of October 2018. Each geolocated occurrence is labelled with the species of higher identification probability. This dataset is thus very heterogeneous in species identification quality, due to the high variability of the image quality submitted by users. The confidence score is provided to Geo-LifeCLEF participants as specific field in this dataset, who can use it to account for identification uncertainty in their models. This data set contains 2,377,610 occurrences covering 3,906 plant species.

Pl@ntNet filtered data. (PL filtered) We proposed a filtered version of the previous dataset based on species identification quality. We only kept the occurrences for which the first species probability value was above 0.98. This score has been determined by expert to give a reasonable degree of identification confidence. This set of 237,087 occurrences covers 1,364 species.

GeoLifeClef 2018. (GBIF) Train and test occurrences datasets from the previous year edition [START_REF] Botella | Overview of geolifeclef 2018: location-based species recommendation[END_REF] were merged to feed the current challenge. Those plants occurrences were extracted from the Global Biodiversity Information Facility8 . This set of occurrences is around ten times smaller than the Pl@ntNet dataset, as shown in Figure 1. Within this dataset, occurrences are often aggregated on a same geographic point, which denotes uncertain or degraded geolocation. However, the geolocation certainty field is often missing. It contains 281, 952 occurrences covering 3, 231 plant species.

Occurrences of other kingdoms. (GBIF) This data source is made of species that are not plants, but may interact somehow with plants (e.g. trophic, pollination, symbiosis, use of plant as habitat or shelter), and are thus likely to carry interesting correlations with plant species presences. None of those species are in the list of species to predict in the test set (which are only plant species). Those occurrences have also been extracted from the GBIF; based on the following filters: { Basis of record: Human, Location : include coordinates, Country or area : France }. We extracted occurrences from 7 non-plant taxonomic groups:

-Chordata/ Aves (8,000,000).

-Chordata/ Mammalia (1,300,000) -Chordata/ Amphibia (300,000) -Chordata/ Reptilia (200,000) -Arthropoda/ Insecta (3,250,000) -Arthropoda/ Arachnida (70,000) -Fungi/ Basidiomycota (50,000) It contains 10,618,839 occurrences in total covering 23,893 taxa.

Taxonomic and geographic filters applied to all datasets. Because scientists do not name species by the same way in all regions of the world, many official lists of species names, called referentials, co-exist. There are no exact matching between them (in particular because of the new scientific knowledge acquired during the period between the creation of two separate lists) except those suggested by the scientific latin names themselves. In our case, the distinct data sources don't use the same referentials. Furthermore, distinct species names might be considered as redundant (synonyms) in some referentials. GBIF uses its own referential made from several taxonomic referentials, and GBIF occurrences may not be at the species taxonomic level, but at sub-species, or genus, etc. Pl@ntNet data includes occurrences from several plants taxonomic referentials (like The Plant List9 , GRIN10 , the French National plant list, etc.). Thus, for attributing species identifiers in GeoLifeCLEF, it was important to first match all occurrences names to a single taxonomic referential adapted for the French Flora. We chose to use Taxref v1211 referential. We only kept names matching Taxref v12 according to an exact matching algorithm (R script provided on Github12 ). Some true species might have been lost due to distinct spelling between the GBIF taxonomy and Taxref. We only kept points falling inside the French territory (Polygon from GADM 13 ) or inside a 30 meters buffer zone, to account for geolocation uncertainty. Finally, occurrences were randomly shuffled to avoid any bias introduced by their order of use.

Environmental data

Geographic rasters. The geographic and environmental data proposed to participants are a compilation of geographic rasters. The variables represented are often used for the purpose of species distribution modelling, especially for plants.

The nature of values stored in the rasters are quantitative (bioclimatic, topological, hydrographical and evapo-transpiration variables), ordinal (pedological variables) or categorical (land cover). The rasters are extracted from the data repository of Botella [START_REF] Botella | A compilation of environmental geographic rasters for sdm covering france (version 1)[END_REF], where readers can find a detailed description. Tensors extraction. To facilitate the learning of representations taking into account the spatial structure of the environment, we provided a Python toolbox14 allowing to extract local environmental tensors from any position in the rasters. By default, it extracts for each raster a 64x64 pixels patch centered on the target position and aggregate the patches from all rasters in the form of a tensor of size nx64x64 where n is the number rasters.

Test data

We have chosen an independent and unpublished source dataset of occurrences for the test set. It is extracted from the SILENE database maintained by the Conservatoire Botanique Mediterranéen15 . Those observations come from various providers including the conservatory himself, but also national parks, botanical associations or impact study consultants. We removed species (i) that were not present in the train set, (ii) vulnerable species according to the SINP referential "espèces sensibles"16 , (iii) and species that are at least vulnerable according to the IUCN red list 17 . This dataset has a high degree of identification certainty because only botanical experts contribute to it. Its geolocation certainty is under 50 meters. We used random weighted selection scheme to draw 25,000 test occurrences among the 700,000 of the initial set noted S. We compute, for each occurrence s i in S a weight w i :

w i = 1/(n i × r i )
Where r i is the number of species in the neighborhood of s i defined by a circle of radius d. n i is the total number of occurrences in the neighborhood. We define the spatial scale d = 2 kilometers. With these weights and the following algorithm, we guaranty that (i) test occurrences are uniformly distributed in the geographic space at scale 2d, (ii) there is as many occurrences of each present species on neighborhoods of radius 2d. We then draw the test occurrences from S without replacement, through the following algorithm:

-Initialize the bag of test occurrences S := S and the test set T = ∅.

-Randomly draw an occurrence in S , say i.

-Draw a scalar z ∼ U (0, max(w 1 , ..., w | S|)).

-If z < w i , remove i from S and add it to T , otherwise leave it in S .

-Stop if |T | = 25000, otherwise we go back to step (1).

Task description

For every occurrence of the test set, the evaluated systems must return a list of 50 species maximum, ranked without ex-aequo. The main evaluation metric used is the top 30 accuracy (TOP30). We provide its expression hereafter:

TOP30 : 1 Q Q q=1 1 rankq≤30
where Q is the total number of query occurrences x q in the test set and rank q is the rank of the correct species y(x q ) in the ranked list of species predicted by the evaluated method for the occurrence x q . A secondary metric is the Mean Reciprocal Rank (MRR), a statistic measure for evaluating any process that produces a list of possible responses to a sample of queries ordered by probability. The reciprocal rank of a query response is the multiplicative inverse of the rank of the correct answer. We provide its expression hereafter :

MRR : 1 Q Q q=1 1 rank q
The MRR was used as main metric during last year edition. We compute it this year, in order to enable comparisons between two campaigns.

Participants and methods

61 participants registered to the challenge through the online platform, among which 5 participants managed to submit runs in times. A total of 44 runs were submitted. All participants runs methods are characterized by their types of model architecture, the occurrences and input data they used in table 6. In the following paragraph, we describe in more details the methodology of each team.

LIRMM, Inria, Univ. Paul Valery, Univ. Montpellier, France, 4 runs, [START_REF] Negri | Plant prediction from cnn model trained with other kingdom species (geolifeclef 2019: Lirmm team)[END_REF] : This team used a single deep convolutional neural network architecture derived in four models. All models take as input the default environmental tensors extracted by the provided python toolbox (see section 2.1), with a one-hot encoding transformation for each category of the land cover variables (clc), inducing 77 layers images in the input of the model. The chosen architecture was an Inception V3 ( [START_REF] Szegedy | Rethinking the inception architecture for computer vision[END_REF]). Models were trained as classifiers, using a softmax output and a cross-entropy loss (also known as multinomial logistic regression). Model of run 27006 was trained on all occurrences of PL complete and glc18 datasets, while models 27004 used PL complete with identification score ≥ 0.7, and 27005 used PL complete with identification score ≥ 0.98 (filtered dataset). Furthermore, runs 27004 and 27005 were only trained on a subset of the occurrences: a sample of around 30K occurrences was drawn according to the same selection procedure as for the test set. Thus, all those models predicted only plant species. On the contrary, model 27007 was trained on all occurrences datasets including PL complete, glc18 and also noPlants. This one was trained to predict plant species and many animal species.

SaraSi, EcoSols, UMR 1222 INRA -Montpellier SupAgro, France, 5 runs, [START_REF] Si-Moussi | Species recommendation using environment and biotic associations[END_REF] : This team used mainly two types of models: a convolutional neural network (CNN) based on the environmental tensors in the same spirit as LIRMM (27086, 27087, 27088) with a customized architecture, and a deep neural network using only a vector of co-occurrences of non-plants taxa as input (27089, 27082). The CNN model architecture separates the feature extraction depending on the type of variables that is deal with. Indeed, it apply distinct convolutional layers to the three categories of environmental patches (continuous, ordinal and categorical). The extracted features are concatenated and used as input in a series of fully-connected layers. A noticeable technique of "categories embedding" was used for the categorical and ordinal patches. It transforms the one-hot encoded patches in a lower number of continuous valued matrices. Also, they addressed the class imbalance of the training set by optimizing a weighted cross-entropy loss so that occurrences of more abundant species were less numerous. They trained this model on the P L complete dataset (27086) and on a reduced version of this dataset to test set species (27088). the run 27087 was like 27086 but trained longer. For the other approach they implemented a customized version of the Continuous Bag of Words model [START_REF] Mikolov | Efficient estimation of word representations in vector space[END_REF]. The input is a set of identifiers of the non-plant "super-taxa" occurring in the neighborhood. An embedding vector associated to the set of "super-taxa" is learned. A "super-taxa" is an aggregation of many species assumed to share a same type of interaction with plants. They were determined through experts knowledge.

SSN CSE, SSN College of Engineering of Chennai, and VIT University of Vellore, India, 12 run, [START_REF] Krishna | Species recommendation using machine learning -geolifeclef 2019[END_REF] : This team tackles the challenge with classical machine learning techniques. They relied on three datasets : (i) spatial position of the occurrences only, (ii) spatial position and punctual environmental vector at the position of the occurrence, (iii) spatial position and vector of the average value of the environmental variables within a 16x16 pixels square centered on the occurrence. As a baseline, the authors first propose a probabilistic model where the probability of a species depends on its frequency in the whole training set (Const. prior). In addition, the authors relied on three categories of models. They first used random forest with spatial coordinates only as input (27102), and boosted trees (XGBoost: 26997, 26996, 27O13, 27012, 26988) and artificial neural network (27069, 27070, 27064, 27067) for using either spatial positions, environmental vectors or both. For one neural network, the authors split the features in 5 groups and trained a neural network per group for which predictions are then combined to form a single model.

Atodiresein, Faculty of Computer Science, "Alexandru Ioan Cuza" University, Romania, 20 runs [START_REF] Atodiresei | Location-based species recommendationgeolifeclef 2019 challenge[END_REF] : This team based their runs on standard machine learning algorithms: nearest neighbors (K-NN), random forests (Rand. For.), boosted trees (XGBoost) and deep neural networks (ANN). Those algorithms were applied to either the P L complete or P L trusted datasets. They used either the spatial coordinates or the environmental punctual values of a selection of 29 environmental variables, or the concatenation of coordinates and variables. All combinations of algorithms, occurrences data and input data were evaluated on a validation set and the best of them were submitted. They also carried ensemble predictions from those models (runs 26969, 26970, 26958, 27062, 26960, 26971, 26961, 26964, 26968). A partial explanation of the low performances of their runs is that they only answered a short list of species (maximum 5) for each test occurrences, which lowers down performances a lot, especially for the top30 metric.

Lot of Lof, Inra, France, 3 runs, [START_REF] Monestiez | Location-based species recommendation -geolifeclef 2019 challenge[END_REF] : This team used occurrences density estimation based on log-linear spatial in-homogeneous Poisson point processes (PPP). They used a restricted set of environmental variables to model the distribution of occurrences based on expert knowledge: etp, alti, chbio 5, chbio 12, awc top, bs top, slope and aggregated clc in 5 land covers categories. They built their models with the 141 test species having the most occurrences in the P L trusted dataset. Run 27124 is the standard PPP, while runs 27123 and 27063 apply different corrections for spatial sampling bias.

Results and discussion

The TOP30 and MRR evaluation scores achieved by all submitted runs are provided in Figures 3 and4 (numerical values of the TOP30 are also replicated in the third column of Table 6). As a complementary analysis, Figure 5 displays the average TOP30 accuracy obtained for each species in the test set as a function of the number of occurrences of this species in the test set.

These results contributes to drive the following findings:

The occurrences of the other kingdoms significantly improve plants prediction. This can be observed from the comparison of run 27007 and run 27006 of the LIRMM team which are all things equal except the use of the occurrences of other kingdoms. The TOP30 increases from 0.136 to 0.177, which represents an improvement of 30%. The use of the occurrences of the other kingdoms is therefore the main cause of the best performances obtained by this team with regard to the SaraSi team. From the ecological point of view, this suggests that the biotic interactions (competition, predation, facilitation) between plant species and other biological groups play a very important role in determining the distribution of the species. From a deep learning point of view, it means that the convolutional neural network is able to transfer a consistent knowledge from the domain of the other kingdoms to the plant domain. An architecture that aim at predicting so many species through mutual neurons (as run 27007) might be a more efficient design for learning those relationships than using the co-occurrences as input data (as did runs 27089, 27082). It would be interesting to investigate this by comparing the latter strategy with a model taking both environmental patches and co-occurrences as input. Weighting the loss by species is better for predicting rare species. The CNN models learnt by the SaraSi team were based on a weighted crossentropy loss penalizing the classes with more samples as a way to compensate class imbalance. Interestingly, it can be seen in Figure 5 that this significantly increased the ability of the model 27086 to predict the species having few occurrences compared to the winner CNN (run 27007) from LIRMM. Run 27086 is better than 27007 for more than 80% of the species. LIRMM team gave equal weights to all occurrences in the loss for training model 27007. It also shows how the most represented species hide the performances on the majority of species, which rarely occur. Giving more balanced weights across species is certainly important to achieve more robust predictions because the observation preferences across species vary a lot from one biodiversity dataset to another, as it is the case here between Pl@ntNet, the GBIF and SILENE.

The more complex the model, the better the prediction. The analysis of the column "model" of Table 6 suggests that, at least models using environmental inputs, can be ranked according to their performance as: (i) Convolutional Neural Network (CNN), (ii) Boosted trees (XGboost), (iii) Deep Neural Network (ANN), (iv) Poisson point processes, (v) K-Nearest Neighbors. This clearly shows a gradient from the models that integrate the most complex input data (CNN having the most complex with many channels of environmental images) and the most flexible architectures (CNN, XGBoost and ANN can fit very complex functions of their input data), to the models that are the most constrained by their input data (environmental vectors only) and with simple architectures (log-linear model of PPP, no optimized parameters for K-NN). This shows that the size of the available datasets and the complexity of the problem give a real advantage to complex statistical learning methods. More specifically, once again CNN results far exceeded those of the other methods which reinforces the results obtained in the last edition of the challenge. The CNN are likely to extract complex features of spatio-environmental patterns in their highest level neurons which are more suited to describe species habitats than environmental variables designed by experts. They may also captures spatial configurations of habitats that favor certain dispersion mechanisms, e.g. source-think colonization, or detect signatures of particular trophic assemblages.

The training of CNN can fail. Although the best models were based on CNNs, not all CNNs obtained so good results. Indeed, some runs based on CNNs were even worst than the prior ranking of species according to their global abundance (see 27004 ≤ 26821). Furthermore, non-submitted CNN models mentioned in a participant working note did perform less in validation than simpler approaches (see [START_REF] Krishna | Species recommendation using machine learning -geolifeclef 2019[END_REF] 3.4). Model design (architecture, selection of environmental channels, management of categorical variables), regularization (optimization algorithm, use of dropout, learning rate and stopping rule policy), training data (especially size, see runs 27004 and 27005) and occurrence weighting scheme de-termine jointly the implementation success. Results of the MRR show that performances were globally lower than last year. Indeed, last year average MRR of the ten best runs was 0.039 while it is 0.024 this year. This large global performance gap is probably due to the difficulty of the test set, given that last year dataset was included in the training data. We note that the test set was not identically distributed, firstly because it was located on the Mediterranean region only, but also because the occurrences were sampled to avoid spatial and species biases. We know that all models predict less well rare species and under-sampled areas. Thus, this drop in overall performance supports the idea that the new test set has succeeded in giving greater importance to rare species and sub-sampled areas. In absolute terms, the best run gives the good answer 20% of the times in its top-30. Thus, roughly speaking, even the best model gives generally a large majority of wrong species in its top-30 list. To give an order of comparison, the database Sophy [START_REF] Ruffray | sophy", une banque de données phytosociologiques; son intérêt pour la conservation de la nature. Actes du colloque "Plantes sauvages et menacées de France: bilan et protection[END_REF] contains more than 35,000 exhaustive plant species inventories on plots generally not exceeding 400m 2 , and covers a wide range of environments in France. According to it, the species diversity in such plots is 25 in average and rarely exceeds 70. There is thus large room for improvement in automated predictions.

Conclusion and perspectives

We now come back on the main outcomes of this task and discuss its perspectives. LIRMM best CNN successfully integrated many non-plants species occurrences in their models predictions to better extract spatio-environmental patterns that more robustly predict plants species. It suggests that the global biotic assemblage highly determine the plant assemblage through underlying species interactions, and the multi-species prediction proved again to be a good deep learning strategy to account for it. This is the main new outcome of this year's edition. However, there should be significant room for improvement in the implementation of this approach. Indeed, LIRMM indicated that the winning model training couldn't be finished for time constraints reasons. Furthermore, light and customized models architectures accounting for the different variables natures seem more adapted to the problem than heavily parameterized state-of-the-art image classification architectures. Indeed, SaraSi customized CNN architecture has performed better than the related LIRMM Inception V3 CNN with the same output. Merging the strengths of both strategies promises good improvements in the future. A rich source of information that remains unexploited for this task is the high resolution satellite images data. For example, today, 50 cm resolution satellite images are freely available for research all over the french territory through the National Institute of Geography (IGN) 18 . Including such images as input in the current models would inform them about very local land cover type and thus give much finer resolution prediction, if one can efficiently handle the size of this data. The philosophy of the evaluation was to favor models that are more robust to biases in the training data, especially the imbalance of species representation and the heterogeneous spatial coverage, both consequences of the reporting process heterogeneity. We can say that it is a success concerning species imbalance representation. Indeed, SaraSi achieved remarkably stable performances even for rare species through a per class weighting scheme in the cost function. A next step would be to account for spatial sampling heterogeneity, as we have seen that all methods still struggle a lot with scarcely reported areas. Regarding the evaluation process on this problem globally, we put an effort this year in the quality of the occurrences identification, and corrected for the species imbalance bias and heterogeneous spatial coverage (due to the reporting heterogeneity). Our new evaluation strategy was quite discriminant across the methods, and lowered globally the computed results. In absolute terms, we have also seen that even the best model tends to rank a lot of relevant species (i.e. probably absent from the surroundings) before the good one. The problem of spatial prediction of plant species lists is objectively far from being solved. Still, with the new areas of improvements that the task results pointed out, we are optimistic about the future methodological advances on the problem of location based species prediction.

Table 1. Results and summarized methodology description of all runs submitted to Ge-oLifeCLEF 2019. Symbols and abbreviations: A + B means that variables/data B was added to A. A\B means that variables/data B where removed from A. complete ∩ test means that only test species occurrences from the complete dataset were used. Products (×) and exponent notations in column "model archi." decompose an ensemble methods with its different models. Occurrences: complete=P L complete,filtered=P L f iltered, all plants=P L complete + P L f iltered + glc18, all=P L complete + P L f iltered + glc18 + nonP lants. Covariates in model input: "enviro. tensors"=environmental tensors with spatial neighborhood", "enviro."=punctual values of environmental variables, "coord."= spatial coordinates.

Fig. 1 .

 1 Fig. 1. Number of occurrences per training dataset. Trusted occurrences were determined from Pl@ntNet species identification engine certainty score.

Fig. 2 .

 2 Fig. 2. Patch extracted at the city of Brest, France.

Fig. 3 .

 3 Fig. 3. Average Top30 accuracy per run and participant. It was computed over the 25,000 test occurrences. This was the official ranking metric for the task.

Fig. 4 .

 4 Fig. 4. Mean Reciprocal Rank per run and participant. It was computed over the 25,000 test occurrences.

Fig. 5 .

 5 Fig. 5. Top30 accuracy averaged per species abundance class for the two best CNN models. Species were ranked by decreasing number of occurrences in the test set and then aggregated in 14 classes of abundances. For run 27086, each occurrence is weighted inversely proportional to the abundance of its species in the loss function.

https://www.gbif.org/

https://plantnet.org

https://www.gbif.org/

http://www.theplantlist.org/

https://www.ars-grin.gov/

https://inpn.mnhn.fr/programme/referentiel-taxonomique-taxref?lg=en

https://github.com/maximiliense/GLC19/blob/master/GITHUB_taxonomic_and_ spatial_filtering.R

https://gadm.org/

https://github.com/maximiliense/GLC19

http://flore.silene.eu/index.php?cont=accueil

http://www.naturefrance.fr/languedoc-roussillon/ referentiel-des-donnees-sensibles

https://geoservices.ign.fr/documentation/geoservices/


Denis Efimov

Andrey Polyakov

Alexander Aleksandrov

Discretization of Homogeneous Systems Using Euler Method with a State-Dependent Step

Keywords: Homogeneous systems, Discretization, Euler method

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

In many applications non-asymptotic convergence rates are preferable, and the theory of homogeneous dynamical systems are frequently used to guarantee these fast rates [START_REF] Bernuau | On homogeneity and its application in sliding mode[END_REF]. Examples of such an application of homogeneity can be found in the sliding mode control domain [START_REF] Levant | Homogeneity approach to high-order sliding mode design[END_REF], in the designs seeking for control or estimation in finite-time [START_REF] Grüne | Homogeneous state feedback stabilization of homogeneous systems[END_REF][START_REF] Hong | H∞ control, stabilization, and input-output stability of nonlinear systems with homogeneous properties[END_REF][START_REF] Bhat | Geometric homogeneity with applications to finite-time stability[END_REF][START_REF] Moulay | Finite time stability of differential inclusions[END_REF][START_REF] Li | Global finite-time stabilization by dynamic output feedback for a class of continuous nonlinear systems[END_REF][START_REF] Shen | Semi-global finite-time observers for nonlinear systems[END_REF][START_REF] Yiguang | Finitetime input-to-state stability and applications to finite-time control design[END_REF] or in fixed-time [START_REF] Andrieu | Homogeneous Approximation, Recursive Observer Design, and Output Feedback[END_REF][START_REF] Cruz-Zavala | Uniform robust exact differentiator[END_REF][START_REF] Polyakov | Nonlinear feedback design for fixed-time stabilization of linear control systems[END_REF][START_REF] Polyakov | Finite-time and fixed-time stabilization: Implicit lyapunov function approach[END_REF]. Another reason of popularity of homogeneous systems is due to that they take an intermediate place between linear and nonlinear dynamics [START_REF] Bacciotti | Lyapunov Functions and Stability in Control Theory[END_REF] (the linear dynamical systems are also homogeneous). The common advantages of homogeneous systems include scalability of trajectories, robustness to measurement noises and exogenous disturbances [START_REF] Yu | On the asymptotical stability of solutions of nonstationary differential equation systems with homogeneous right hand sides[END_REF][START_REF] Peuteman | Averaging results and the study of uniform asymptotic stability of homogeneous differential equations that are not fast time-varying[END_REF][START_REF] Bernuau | Verification of ISS, iISS and IOSS properties applying weighted homogeneity[END_REF], e.g. as the linear dynamics. From another side, these systems are This work was partially supported by ANR 18 CE40-0008 (Project DIGITSLID), the Government of Russian Federation (Grant 08-08) and Saint Petersburg State University (Project Id 37569826).

Email addresses: denis.efimov@inria.fr (Denis Efimov), andrey.polyakov@inria.fr (Andrey Polyakov), a.u.aleksandrov@spbu.ru (Alexander Aleksandrov).

described by essentially nonlinear differential equations, where some interesting nonlinear phenomena can be observed, namely finite-or fixed-time convergence rates for negative and positive degrees of homogeneity, respectively, and robustness with respect to delays [START_REF] Efimov | Weighted homogeneity for time-delay systems: Finite-time and independent of delay stability[END_REF][START_REF] Zimenko | A note on delay robustness for homogeneous systems with negative degree[END_REF].

The models of industrial plants are usually described by continuous-time ordinary differential equations, then frequently the control and estimation algorithms are developed in continuous-time setting. Hence, design and analysis of stability or performance of closed-loop systems are performed in continuous time, while for a verification or implementation, the solutions of closed-loop systems have to be calculated in a computer or in a digital controller, which are usually operating in a discrete time. Different numerical discrete-time approximation methods and time discretization schemes are proposed for continuous-time differential equations in order to realize such an implementation [START_REF] Allen | Numerical Analysis for Applied Science[END_REF][START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF]. Their variety is originated by the fact that these approaches guarantee the preservation of performances obtained in continuous time for discrete-time approximations only under some particular hypothesis, then their application is case-dependent. For the control or estimation algorithms designed using the theory of homogeneous sy-stems, which possess an increased rate of convergence with respect to linear systems, the problem of implementation and derivation of solutions becomes important being in the focus of many works, especially in the sliding mode control community [START_REF] Levant | On fixed and finite time stability in sliding mode control[END_REF][START_REF] Efimov | Realization and discretization of asymptotically stable homogeneous systems[END_REF].

The conventional Euler method is a first-order numerical routine for solving ordinary differential equations with a given initial value and a fixed time step, which represents the most basic explicit/implicit method of numerical integration and the simplest Runge-Kutta method, which is extremely popular in the control community. In particular, in the recent papers [START_REF] Efimov | Realization and discretization of asymptotically stable homogeneous systems[END_REF][START_REF] Levant | Stability and robustness of homogeneous differential inclusions[END_REF] it has been shown that application of the explicit Euler method (with a fixed time step) for the global approximation of solutions of homogeneous systems with non-zero degree is problematic (see also [START_REF] Levant | On fixed and finite time stability in sliding mode control[END_REF]), and the implicit Euler scheme has a better perspective (as it has been also indicated in [START_REF] Acary | Implicit euler numerical simulations of sliding mode systems[END_REF][START_REF] Brogliato | Globally stable implicit euler time-discretization of a nonlinear single-input sliding-mode control system[END_REF][START_REF] Huber | Lyapunov stability and performance analysis of the implicit discrete sliding mode control[END_REF][START_REF] Miranda-Villatoro | Multivalued robust tracking control of lagrange systems: Continuous and discrete-time algorithms[END_REF]). In other papers, the conditions of convergence and stability of the explicit and implicit Euler methods have been studied for linear systems (the notion of A-stability) [START_REF] Germund | A special stability problem for linear multistep methods[END_REF][START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF], or for particular classes of nonlinear systems [START_REF] Gonzalez | Backward euler discretization of fully nonlinear parabolic problems[END_REF][START_REF] Merlet | Convergence to equilibrium for the backward euler scheme and applications[END_REF][START_REF] Yu | Preservation of stability under discretization of systems of ordinary differential equations[END_REF].

Following the results obtained in [START_REF] Efimov | Realization and discretization of asymptotically stable homogeneous systems[END_REF], the present work is devoted to application of explicit and implicit Euler discretization schemes for approximation of solutions of homogeneous stable dynamical systems. The case of a state-dependent scaling of the time discretization step is considered, and different conditions for the existence (in the implicit case) and convergence to zero of approximations for the explicit and implicit Euler integration methods are derived. It is shown that introduction of a state-dependent re-scaling of the time discretization step allows the finite-or fixed-time rates of convergence to be recovered by the discrete-time approximations of solutions, but for an infinite number of steps (in a finite number of steps the convergence to a vicinity of the origin is obtained). Absolute and relative errors (closeness of the approximations to real solutions) for the explicit and implicit Euler integration schemes are investigated using the homogeneity theory, and it is shown that the relative errors are globally bounded and by decreasing the initial discretization step it is possible to make them arbitrary small, hence, the proposed modification of the Euler method can be indeed used for calculation of solutions of homogeneous dynamics. The preliminary formulation of these results (without proofs) appears in [START_REF] Efimov | On statedependent discretization of stable homogeneous systems[END_REF]. In the present work, the proofs are given and, in addition, since stable homogeneous systems admit certain robustness with respect to external inputs [START_REF] Bernuau | Verification of ISS, iISS and IOSS properties applying weighted homogeneity[END_REF], it is shown that the discretized approximations obtained with the proposed methods keep this property.

The outline of this note is as follows. The problem statement and some preliminary results are introduced in Section 2. Some basic properties and relations between solution approximations are studied in Section 3. The convergence conditions are established in Section 4. Pro-perties of relative and absolute errors of approximation of solutions of homogeneous systems by the explicit and implicit Euler methods are investigated in Section 5. Analysis of stability robustness of approximations with respect to external disturbances is presented in Section 6. Simple illustrating example is considered in Section 7.

Notation

• R + = {x ∈ R : x ≥ 0},

Preliminaries

In this work the following nonlinear system is considered:

ẋ(t) = f (x(t)), t ≥ 0, (1)
where x(t) ∈ R n is the state, f : R n → R n ensures forward existence and uniqueness of the system solutions at least locally, f (0) = 0. For an initial condition x 0 ∈ R n define the corresponding solution by X(t, x 0 ) for any t ≥ 0 for which the solution exists. If f is discontinuous, then the solutions are understood in the Filippov's sense [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF] (it is assumed that f (x) ∈ Φ(x), where Φ is the corresponding Filippov's set-valued extension).

Following [START_REF] Roxin | On finite stability in control systems[END_REF][START_REF] Khalil | Nonlinear Systems. NJ 07458[END_REF][START_REF] Polyakov | Nonlinear feedback design for fixed-time stabilization of linear control systems[END_REF], let Ω ⊂ R n be an open set, 0 ∈ Ω.

Definition 1 At the steady state x = 0 the system (1) is said to be (a) stable if for any x 0 ∈ Ω the solution X(t, x 0 ) is defined for all t ≥ 0, and for any > 0 there is δ > 0 such that for any x 0 ∈ Ω, if x 0 ≤ δ then X(t, x 0 ) ≤ for all t ≥ 0;

(b) asymptotically stable if it is stable and for any κ > 0 and > 0 there exists T (κ, ) ≥ 0 such that for any x 0 ∈ Ω, if x 0 ≤ κ then X(t, x 0 ) ≤ for all t ≥ T (κ, );

(c) finite-time stable if it is stable and finite-time converging from Ω, i.e. for any x 0 ∈ Ω there exists 0 ≤ T < +∞ such that X(t, x 0 ) = 0 for all t ≥ T . The function T 0 (x 0 ) = inf{T ≥ 0 : X(t, x 0 ) = 0 ∀t ≥ T } is called the settling time of the system (1);

(d) fixed-time stable if it is finite-time stable and sup x0∈Ω T 0 (x 0 ) < +∞.

The set Ω is called a domain of stability/attraction.

If Ω = R n , then the corresponding properties are called global stability/asymptotic stability/finite-time/fixedtime stability of (1) at x = 0.

Similarly, the stability notions can be defined with respect to a compact invariant set, by replacing the distance to the origin in Definition 1 with the distance to a set.

Weighted homogeneity

Following [START_REF] Zubov | On systems of ordinary differential equations with generalized homogenous right-hand sides[END_REF][START_REF] Bacciotti | Lyapunov Functions and Stability in Control Theory[END_REF][START_REF] Levant | Homogeneity approach to high-order sliding mode design[END_REF], for strictly positive numbers r i , i = 1, n called weights and λ > 0, define:

• the vector of weights r = (r 1 , . . . , r n ) , r max = max 1≤j≤n r j and r min = min 1≤j≤n r j ;

• the dilation matrix function Λ r (λ) = diag{λ ri } n i=1 , note that ∀x ∈ R n and ∀λ > 0 we have Λ r (λ)x = (λ r1 x 1 , . . . , λ rn x n ) ; • the r-homogeneous norm x r = ( n i=1 |x i | ρ r i ) 1 ρ
for any x ∈ R n and ρ ≥ r max , it is not a norm in the standard sense, since the triangle inequality is not satisfied for • r , however there exist σ, σ ∈ K ∞ such that σ( x r ) ≤ x ≤ σ( x r ) ∀x ∈ R n ;

• the sphere and the ball in the homogeneous norm

S r ( ) = {x ∈ R n : x r = } and B r ( ) = {x ∈ R n : x r ≤ } for ≥ 0. Definition 2 A function g : R n → R is r-homogeneous with degree µ ∈ R if ∀x ∈ R n and ∀λ > 0 we have: λ -µ g(Λ r (λ)x) = g(x).
A vector field f : R n → R n is r-homogeneous with degree ν ∈ R, ν ≥ -r min if ∀x ∈ R n and ∀λ > 0 we have:

λ -ν Λ -1 r (λ)f (Λ r (λ)x) = f (x),
which is equivalent for i-th component of f being a r-homogeneous function of degree r i + ν.

System (1) is r-homogeneous of degree ν if the vector field f is r-homogeneous of the degree ν.

The homogeneous norm • r is an example of r-homogeneous function with degree 1.

Theorem 3 [START_REF] Zubov | On systems of ordinary differential equations with generalized homogenous right-hand sides[END_REF][START_REF] Rosier | Homogeneous Lyapunov function for homogeneous continuous vector field[END_REF] For the system (1) with r-homogeneous and continuous function f the following properties are equivalent:

• the system (1) is asymptotically stable;

• there exists a C 2 r-homogeneous Lyapunov function V : R n → R + such that ∀x ∈ R n and ∀λ > 0

α 1 ( x ) ≤ V (x) ≤ α 2 ( x ), L f V (x) = ∂V ∂x (x)f (x) ≤ -α( x ), λ -µ V (Λ r (λ)x) = V (x), µ > r max , for some α 1 , α 2 ∈ K ∞ and α ∈ K.
The requirement on continuity of the function f has been relaxed in [START_REF] Bernuau | On extension of homogeneity notion for differential inclusions[END_REF].

Euler schemes

Let us consider a class of systems in (1) with sup ξ∈Sr( 1)

f (ξ) < +∞
and admitting the following hypothesis:

Assumption 1 Let (1) be r-homogeneous with a degree ν ∈ R and asymptotically stable.

In this work, in order to obtain a discretization of the solutions of (1), we will propose modifications of implicit and explicit Euler methods [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] with state-dependent scaling of the time discretization step. To this end, select a basic discretization step h > 0, define a sequence of time instants t i for i = 0, 1, . . . such that t 0 = 0 and t i+1 -t i > 0, and denote by x i an approximation of the solution X(t i , x 0 ) at the corresponding time instant (i.e.

x i X(t i , x 0 ) and x 0 = X(0, x 0 )), then the approximation x i+1 calculated in accordance with the explicit Euler method is given by:

x i+1 = x i + h x i ν r f (x i ), (2)
t i+1 = t i + h x i ν r
for i = 0, 1, . . . , while the approximation calculated by the implicit Euler method comes from:

x i+1 = x i + h x i+1 ν r f (x i+1 ), (3)
t i+1 = t i + h x i+1 ν
r for i = 0, 1, . . . In the algorithms (2) and (3) it is assumed that x i = 0 or x i+1 = 0 since in the opposite case the discretization stops at the equilibrium due to f (0) = 0.

Remark 4

In order to clarify the relations of the results obtained in this work and in [START_REF] Efimov | Realization and discretization of asymptotically stable homogeneous systems[END_REF], and also to motivate the selected design, for the system (1) under Assumption 1 let us consider a function

ψ x0 (t) = t 0 X(s, x 0 ) ν r ds
for any x 0 ∈ R n , which is well-defined and invertible when the trajectory stays out of the origin, i.e. for all t ∈ [0, T 0 (x 0 )) where T 0 : R n → R + ∪ {+∞} is possibly infinite time of convergence to the origin (as in Definition 1). Denote Y (τ,

x 0 ) = X(ψ -1 x0 (τ ), x 0 ) for τ ∈ [0, ψ x0 (T 0 (x 0 ))), then d dτ Y (τ, x 0 ) = d ds X(s, x 0 ) s=ψ -1 x 0 (τ ) 1 X(ψ -1 x0 (τ ), x 0 ) ν r = f (X(ψ -1 x0 (τ ), x 0 )) X(ψ -1 x0 (τ ), x 0 ) ν r = f (Y (τ, x 0 )) Y (τ, x 0 ) ν r . (4)
Obviously, for any x 0 ∈ R n the asymptotic stability properties of the corresponding solutions X(t, x 0 ) of the system (1) for t ∈ [0, T 0 (x 0 )) and Y (τ, x 0 ) of ( 4) for τ ∈ [0, ψ x0 (T 0 (x 0 ))) are interrelated, and the system (4) is r-homogeneous with degree 0:

f (Λ r (λ)y) Λ r (λ)y ν r = Λ r (λ) f (y) y ν r
for any y ∈ R n and λ > 0. As it has been shown in [START_REF] Efimov | Realization and discretization of asymptotically stable homogeneous systems[END_REF], for the case ν = 0 and a properly established discretization step h > 0 the explicit and implicit Euler methods provide an admissible approximation of the solution y i Y (τ i , x 0 ) with τ i = ih for i = 0, 1, . . . :

y i+1 = y i + h y i ν r f (y i ), y i+1 = y i + h y i+1 ν r f (y i+1 ).
Hence,

τ i+1 -τ i = h = ti+1 ti X(s, x 0 ) ν r ds,
where t i are the corresponding instants of discretization in the original time t, and under an assumption that X(s, x 0 ) ν r const for s ∈ [t i , t i+1 ], i.e. the discretization interval is sufficiently small relative to the velocity of the system, we obtain 2) and (3), respectively. In this work such a scheme is used directly for the system (1) to approximate the solution X(t i , x 0 ), which is a significant difference from [START_REF] Efimov | Realization and discretization of asymptotically stable homogeneous systems[END_REF]. This intuition also highlights our motivation to reduce an Euler computation scheme to work with a system of degree zero.

t i+1 -t i h X(ti,x0) ν r or t i+1 -t i h X(ti+1,x0) ν r for (
Therefore, it is a well-known fact that for ν = 0 with h → 0 both methods approach the real solution [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF], i.e. x i → X(t i , x 0 ) in ( 2) and (3) with h → 0 over any compact time interval. In [START_REF] Efimov | Realization and discretization of asymptotically stable homogeneous systems[END_REF] for the case of the conventional explicit and implicit Euler methods with a constant discretization step h and ν = 0 it has been shown that outside of a vicinity of the origin for ν < 0 or a vicinity of infinity for ν > 0 the same results can be obtained.

In the sequel, the problem of global convergence to zero of the approximations {x i } ∞ i=0 derived in ( 2) and ( 3) is studied for systems in (1) satisfying Assumption 1 with ν = 0 (the case ν = 0, or without scaling of the discretization step, has been analyzed in [START_REF] Efimov | Realization and discretization of asymptotically stable homogeneous systems[END_REF]). For this goal we need to establish some auxiliary properties of approximations {x i } ∞ i=0 obtained in ( 2) and (3).

Relations between approximations obtained for different initial conditions

One of the main features of ( 2) and ( 3) is as follows:

Proposition 5 Let the system (1) be r-homogeneous with a degree ν ∈ R. If {x i } ∞ i=0 is a sequence generated by ( 2) or (3) for the time instants {t i } ∞ i=0 with the step h and the initial state x 0 , then for any λ > 0, y i = Λ r (λ)x i is a sequence obtained by ( 2) or (3), respectively, for the instants λ -ν {t i } ∞ i=0 with the step h and the initial state

y 0 = Λ r (λ)x 0 .
PROOF. Fixing λ > 0 and multiplying both sides of (2) by the dilation matrix Λ r (λ) we obtain:

Λ r (λ)x i+1 = Λ r (λ)x i + h λ ν x i ν r λ ν Λ r (λ)f (x i ) = Λ r (λ)x i + h Λ r (λ)x i ν r f (Λ r (λ)x i ),
or, equivalently, y i+1 = y i + h yi ν r f (y i ) that gives the desired result. The proof for the scheme (3) is the same.

Note that y i is an approximation of X(λ -ν t i , y 0 ) for scaled instants of time.

Corollary 6 Let the system (1) be r-homogeneous with a degree ν ∈ R. Let for all x 0 ∈ S r (1) there exist sequences {x i } ∞ i=0 obtained by ( 2) or (3) with the step h > 0 and the initial state x 0 possessing one of the properties:

sup i≥0 x i < +∞;
(5)

lim i→+∞ x i = 0. ( 6
)
Then for any y 0 ∈ R n there exist sequences {y i } ∞ i=0 generated by ( 2) or (3) with the step h and the initial state y 0 possessing the same property.

PROOF. For any y 0 ∈ R n there exists λ > 0 such that y 0 = Λ r (λ)x 0 for some x 0 ∈ S r (1), next the result is a direct consequence of Proposition 5 since y i = Λ r (λ)x i .

Norm deviations

Another property of the proposed approximation algorithms ( 2) and (3) consists in their finite increment: Lemma 7 Let the system (1) be r-homogeneous with a degree ν ∈ R. Then there are h 0 > 0 and 0 < c < c < +∞ such that for any h ∈ (0, h 0 ]: 2) or (3) (provided that a solution x i+1 exists).

c x i r ≤ x i+1 r ≤ c x i r for all x i ∈ R n in (
PROOF. Starting with (2), let x i ∈ R n and denote x i = Λ r ( x i r )ξ i for some ξ i ∈ S r (1), then the equation ( 2) can be rewritten as follows:

x i+1 = Λ r ( x i r ) (ξ i + hf (ξ i )) .
Therefore, x i+1 r = x i r ξ i + hf (ξ i ) r and there exists h 0 > 0 such that 0 < c < c < +∞ for For (3), let x i ∈ R n and denote x i+1 = Λ r ( x i+1 r )ξ i+1 for some ξ i+1 ∈ S r (1), then the equation (3) can be rewritten:

x i = Λ r ( x i+1 r ) (ξ i+1 -hf (ξ i+1 )) ,
and the statement follows repeating the steps as for [START_REF] Yu | On the asymptotical stability of solutions of nonstationary differential equation systems with homogeneous right hand sides[END_REF].

Corollary 8 Let the system (1) be r-homogeneous with a degree ν ∈ R. Select h 0 > 0, then there exists 0 < c < +∞ such that for any h ∈ (0, h 0 ] and for all x i ∈ R n , x i+1 r ≤ c x i r in (2) or c x i r ≤ x i+1 r in (3) (provided that a solution x i+1 exists).

PROOF. This result follows repeating the arguments of the proof of Lemma 7 by taking into account that for a generic f it might be inf h∈(0,h0] inf ξ∈Sr(1) ξ + hf (ξ) r = inf h∈(0,h0] inf ξ∈Sr(1) ξ -hf (ξ) r = 0 without a restriction on the amplitude of h 0 .

4

Convergence of sequences {x i } ∞ i=0 generated by Euler methods

In this section we will study the stability and convergence features of {x i } ∞ i=0 only. The quality of the corresponding approximations (the closeness of {x i } ∞ i=0 to the continuous-time solutions X(t i , x 0 )) will be considered in the next section.

According to Theorem 3, under Assumption 1 for the system (1) there is a twice continuously differentiable and r-homogeneous Lyapunov function V : R n → R + of the degree µ > -ν such that a = -sup ξ∈Sr(1)

L f V (ξ) > 0, 0 < b = sup ξ∈Br(1) ∂V (ξ) ∂ξ < +∞, (7)
0 < c 1 = inf ξ∈Sr(1) V (ξ), 0 < c 2 = sup ξ∈Sr (1)
V (ξ),

c 1 x µ r ≤ V (x) ≤ c 2 x µ r ∀x ∈ R n .

Convergence of the explicit Euler scheme (2)

The principal statement of this subsection is as follows:

Theorem 9 Let Assumption 1 be satisfied, then there exists h 0 > 0 such that for any discretization step h ∈ (0, h 0 ] the sequences {x i } ∞ i=0 obtained by (2) for any initial state x 0 ∈ R n and the step h possess the following properties:

(a) sup i=0,1,... x i r < γ x 0 r for some γ ∈ (0, +∞); (b) lim i→+∞ x i r = 0;

(c) for ν = 0 the sequence {x i } ∞ i=0 has an exponential convergence rate, for ν < 0 for any x 0 ∈ R n the time of convergence to the origin t x0 +∞ = lim i→+∞ t i is finite, and for ν > 0 the time of convergence from any x 0 ∈ R n to B r (1) is also finite independently of x 0 .

PROOF. The proof is given in [START_REF] Efimov | Proofs for "Discretization of Homogeneous Systems Using Euler Method with a State-Dependent Step[END_REF].

The property obtained for ν > 0 is related with the fixed-time rate of convergence (to the unit ball).

Remark 10

The maximal admissible discretization step h 0 can be evaluated as

h 0 = min{1, a v , 2 c 2 a },
where

a = inf y∈Sr(1)
∂V (y) ∂y f (y), c 2 = sup y∈Sr(1)

V (y), v = sup y∈Sr(1)
sup θ∈(0,1)

f (y) ∂ 2 V (ξ) ∂ξ 2 ξ=y+θf (y) f (y)
are calculated on the unit sphere for known f and V . Such an estimate is a lower bound on h 0 .

Existence of approximations in (3)

Existence of some x i+1 ∈ R n for any x i ∈ R n in the explicit case (2) is straightforward, but it is not the case of (3). According to the result of Corollary 6, it is enough to find the conditions of existence of x i+1 for all x i ∈ S r (1) in [START_REF] Yu | Preservation of stability under discretization of systems of ordinary differential equations[END_REF]. Note that in a general case, it is difficult to provide some simple conditions for existence and uniqueness of x i+1 in the equation (3) for any x 0 ∈ S r (1), but homogeneity may further simplify the solution under additional mild restrictions on f .

Proposition 11 Let f be r-homogeneous of degree ν ∈ R and continuous on S r (1). Then there is h 0 > 0 such that for all h ∈ (0, h 0 ) the equation (3) for any

x i ∈ R n has a solution x i+1 ∈ R n .
PROOF. The proof is given in [START_REF] Efimov | Proofs for "Discretization of Homogeneous Systems Using Euler Method with a State-Dependent Step[END_REF].

Thus, for the proposed implicit Euler method with statedependent discretization step (3), the global existence of solutions may be guaranteed by selecting h sufficiently small provided that f is homogeneous and continuous.

Convergence of the implicit Euler scheme (3)

Similar to (2) results can be obtained for (3) and even with some improvement under additional restrictions. Note that in this subsection we assume that the sequences {x i } ∞ i=0 exist in (3) for all x 0 ∈ R n , e.g. the conditions of Proposition 11 are verified.

Theorem 12 Let Assumption 1 be satisfied, then there exists h 0 > 0 such that for any discretization step h ∈ (0, h 0 ] the sequences {x i } ∞ i=0 obtained by (3) for any initial state x 0 ∈ R n and the step h possess:

(a) sup i=0,1,... x i r < γ x 0 r for some γ ∈ (0, +∞); (b) lim i→+∞ x i r = 0;

(c) for ν = 0 the sequence {x i } ∞ i=0 has an exponential convergence rate, for ν < 0 for any x 0 ∈ R n the time of convergence to the origin t x0 +∞ = lim i→+∞ t i is finite, and for ν > 0 the time of convergence from any x 0 ∈ R n to B r (1) is also finite independently of x 0 .

If the matrix ∂ 2 V (ξ) ∂ξ 2
≥ 0 for all ξ ∈ R n , where V : R n → R + is a r-homogeneous Lyapunov function satisfying [START_REF] Bernuau | On homogeneity and its application in sliding mode[END_REF], then h 0 > 0 can be selected arbitrary.

PROOF. The proof is given in [START_REF] Efimov | Proofs for "Discretization of Homogeneous Systems Using Euler Method with a State-Dependent Step[END_REF].

The requirement on nonnegative definiteness of the second derivative of V is related with the condition of convexity level set of V imposed in [START_REF] Efimov | Realization and discretization of asymptotically stable homogeneous systems[END_REF].

Absolute and Relative Errors of Discretized Homogeneous Systems

Conventional characteristics of the discretization accuracy are studied in this section for the proposed algorithms ( 2) or ( 3) and for homogeneous systems. To this end denote Ξ h (x 0 ) = X(h x 0 -ν r , x 0 ) as the value of the solution of (1) with the initial condition x 0 ∈ R n evaluated in ( 2) or (3) after one iteration with the discretization step h > 0 (at t 1 = h x 0 -ν r with t 0 = 0). Denote by Ξh (x 0 ) the estimated value derived by ( 2) or (3) for the same x 0 and h > 0 (note that Ξh (x 0 ) = x 0 + h x0 ν r f (x 0 ) in the case of (2)), then [START_REF] Dahlquist | Numerical methods in scientific computing[END_REF]:

• the absolute error is the magnitude of the difference between the exact value and its approximation:

∆ h (x) = Ξ h (x) -Ξh (x) , ∆ h r (x) = Ξ h (x) -Ξh (x) r ;
• the relative error expresses how large the absolute error is compared with the exact value:

δ h (x) = ∆ h (x) Ξ h (x) , δ h r (x) = ∆ h r (x) Ξ h (x) r .
The errors are given for two different norms, the conventional one • and the homogeneous norm • r , the former one is used habitually for evaluation of a discretization method precision, while the latter norm suits better for analysis of homogeneous systems. Equivalence of these norms has been shown in [START_REF] Efimov | Realization and discretization of asymptotically stable homogeneous systems[END_REF], that is why in this note we will study only ∆ h r (x i ) and δ h r (x i ).

These error functions admit the following useful properties:

Theorem 13 Let the system (1) be r-homogeneous of degree ν ∈ R and Ξh (x) be calculated by the explicit (2) or implicit (3) Euler scheme for x ∈ R n and h > 0.

Then the functions ∆ h r (x) and δ h r (x) are r-homogeneous of degree 1 and 0, respectively. PROOF. First of all note that for a homogeneous system (1) [START_REF] Bacciotti | Lyapunov Functions and Stability in Control Theory[END_REF]:

X(t, Λ r (λ)x) = Λ r (λ)X(λ ν t, x) ∀t ≥ 0, ∀x ∈ R n for any λ > 0, then Ξ h (Λ r (λ)x) = X(h Λ r (λ)x -ν r , Λ r (λ)x) = X(λ -ν h x -ν r , Λ r (λ)x) = Λ r (λ)X(h x -ν r , x) = Λ r (λ)Ξ h (x),
hence Ξ h is r-homogeneous vector field of degree 0. According to Proposition (5) the vector fields Ξh (x) are also r-homogeneous of degree 0 in ( 2) or (3).

Consider the case of the explicit Euler scheme (2), then

∆ h r (x) = Ξ h (x) -Ξh (x) r = Ξ h (x) -x -h x -ν r f (x) r
and for any λ > 0 we obtain:

∆ h r (Λ r (λ)x) = Ξ h (Λ r (λ)x) -Λ r (λ)x -h Λ r (λ)x -ν r f (Λ r (λ)x) r = λ Ξ h (x) -x -h x -ν r f (x) r = λ∆ h r (x),
which implies that the error function ∆ h r (x) is r-homogeneous function of degree 1. Similarly,

δ h r (Λ r (λ)x) = ∆ h r (Λ r (λ)x) Ξ h (Λ r (λ)x) r = ∆ h r (x) Ξ h (x) r = δ h r (x)
for any λ > 0 and the error function δ h r (x) is r-homogeneous function of degree 0.

For the case of the implicit Euler method (3) we have:

∆ h r (x) = Ξ h (x) -Ξh (x) r = Ξ h Ξh (x) -h Ξh (x) -ν r f ( Ξh (x)) -Ξh (x) r
and repeating the same arguments (as it has been shown above, Ξ h (x) and Ξh (x) are r-homogeneous of degree 0) the following facts can be substantiated:

∆ h r (Λ r (λ)x) = λ∆ h r (x), δ h r (Λ r (λ)x) = δ h r (x).
Any homogeneous function of degree 0 is globally bounded (it may be discontinuous) if its maximal amplitude is finite being evaluated on S r (1). Therefore, if for any initial conditions x ∈ S r (1) the error δ h r (x) stays sufficiently small for a reasonable selection of h (i.e. the one step error of usual Euler discretization approaches is small on the sphere), then the explicit (2) and the implicit (3) Euler schemes provide a uniformly bounded relative error δ r globally. Boundedness of δ r implies that the difference between Ξ h (x) and Ξh (x) stays of the order Ξ h (x) (roughly speaking proportional to x). Indeed, assume that f is continuously differentiable on S r (1). In this scenario, the second derivative of the solution X(t, x 0 ) exists and continuous, and for all x 0 ∈ S r (1):

Ξ h (x 0 ) = X(h x 0 -ν r , x 0 ) = x 0 + h Ẋ(0, x 0 ) + h 2 2 Ẍ(θ, x 0 ) = x 0 + hf (x 0 ) + h 2 2 Ẍ(θ, x 0 ),
where θ ∈ [0, h] and the Lagrange reminder of Taylor series expansion was used. Then for x 0 ∈ S r (1) and the explicit discretization algorithm (2):

∆ h r (x 0 ) = Ξ h (x 0 ) -Ξh (x 0 ) r = x 0 + hf (x 0 ) + h 2 2 Ẍ(θ, x 0 ) -x 0 -hf (x 0 ) r = h 2 2 Ẍ(θ, x 0 ) r
for h ≤ 1 and r max = 1 (always can be assumed without loosing generality). By differentiability of f , there exists a constant ς > 0 such that sup

θ∈[0,1], x0∈Sr (1)
Ẍ(θ, x 0 ) r ≤ ς, hence, ∆ h r (x 0 ) ≤ ς 2 h 2
for all x 0 ∈ S r (1) and h ∈ (0, 1]. The quadratic in h convergence of the absolute error ∆ h r (x 0 ) (the same can be shown for the relative error δ h r (x 0 )) implies approaching of the real solution X(h, x 0 ) by its approximation given in (2) [START_REF] Dahlquist | Numerical methods in scientific computing[END_REF]. Remark 14 Note that the time step in (2) and (3) is state-and degree-dependent. In particular, if x i r 1 (it is sufficiently big) and ν < 0, then t i+1 -

t i = x i -ν r h h in (2)
and the time step can be too large, which is also related with the obtained accuracy estimates in Theorem 13. It is worth stressing that anyway the convergence of these algorithms is not influenced, and this observation deals only with the sampling of discretization. Therefore, for big amplitudes of x i in the case of ν < 0 it is admissible to use the conventional Euler methods (without a scaling of the time discretization step), which may provide a reasonable accuracy of approximation of the solutions under a more regular sampling [START_REF] Efimov | Realization and discretization of asymptotically stable homogeneous systems[END_REF][START_REF] Levant | Stability and robustness of homogeneous differential inclusions[END_REF].

Robustness with respect to external inputs

Let us consider a version of the system (1) extended by external inputs:

ẋ(t) = F (x(t), d(t)), t ≥ 0, (8)
with F : R n+m → R n and F (x, 0) = f (x), where d : R + → R m is a measurable and essentially bounded function of time. For x 0 ∈ R n and an input d : R + → R m denote a corresponding solution of (8) as X(t, x 0 , d). The following hypothesis will be imposed on (8):

Assumption 2 There exist rj > 0, j = 1, m and ν ∈ R such that

F (Λ r (λ)x, Λ r (λ)d) = λ ν Λ r (λ)F (x, d)
for all x ∈ R n , d ∈ R m and λ > 0, and there exists σ ∈ K ∞ such that sup y∈Sr( 1)

F (y, d) -F (y, 0) ≤ σ( d r ) for all d ∈ R m .
If Assumption 1 is also satisfied, since F (x, 0) = f (x), the degrees ν in these assumptions coincide. Moreover, if assumptions 1 and 2 are both verified, then the system (8) possesses the input-to-state stability (ISS) property with respect to d [START_REF] Bernuau | Verification of ISS, iISS and IOSS properties applying weighted homogeneity[END_REF], [START_REF] Bernuau | Robustness of Homogeneous and Homogeneizable Differential Inclusions[END_REF], and the zone of asymptotic convergence (the asymptotic gain) can be evaluated as

x r ≤ c 2 c 1 1 µ d r,∞ σ -1 ( a 2b )
, where constants µ, a, b, c 1 and c 2 are introduced in (7),

d r,∞ = ess sup t≥0 d(t) r .
Let us show that (2) and (3) preserve the same stability performance for [START_REF] Bernuau | On extension of homogeneity notion for differential inclusions[END_REF]. Note that being applied to (8) these algorithms can be formulated as follows:

x i+1 = x i + h x i ν r F (x i , d i ), (9)
t i+1 = t i + h x i ν r or x i+1 = x i + h x i+1 ν r F (x i+1 , d i+1 ), (10)
t i+1 = t i + h x i+1 ν r for i = 0, 1, . . . , respectively, where as before x i is an estimate of X(t i , x 0 , d) and d i = d(t i ).

Theorem 15 Let assumptions 1 and 2 be satisfied, then there exists h 0 > 0 such that for any discretization step h ∈ (0, h 0 ] the discrete-time systems [START_REF] Bernuau | Verification of ISS, iISS and IOSS properties applying weighted homogeneity[END_REF] or (10) are ISS 1 , and for any initial state x 0 ∈ R n and any bounded input d : R + → R m the corresponding sequences {x i } ∞ i=0 enter in the set where

x i r ≤ c 2 c 1 1 µ d r,∞ σ -1 ( a 2b )
.

PROOF. If Assumption 1 is verified, then there exists a Lyapunov function V (x) for the system (1) that satisfies the conditions [START_REF] Bernuau | On homogeneity and its application in sliding mode[END_REF].

Let us check the increment of this Lyapunov function V (x) on the sequences produced by (9) (i.e. let us verify that V (x) is an ISS Lyapunov function for (9) [START_REF] Jiang | Input-to-state stability for discrete-time nonlinear systems[END_REF]):

V (x i+1 ) -V (x i ) = V (x i + h x i -ν r F (x i , d i )) -V (x i ) = x i µ r [V (y i + hF (y i , δ i )) -V (y i )],
where y i ∈ S r (1) such that x i = Λ r ( x i r )y i , d i = Λ r ( x i r )δ i for some δ i ∈ R m , and the first property from Assumption 2 has been used. Next, similarly to the proof of Theorem 9, for some θ ∈ [0, 1] we obtain:

V (y i + hF (y i , δ i )) = V (y i ) + ∂V (ξ) ∂ξ ξ=yi hF (y i , δ i ) + h 2 2 F (y i , δ i ) ∂ 2 V (ξ) ∂ξ 2 ξ=yi+θhF (yi,δi) F (y i , δ i ) = V (y i ) + ∂V (ξ) ∂ξ ξ=yi hf (y i ) + ∂V (ξ) ∂ξ ξ=yi h[F (y i , δ i ) -F (y i , 0)] + h 2 2 F (y i , δ i ) ∂ 2 V (ξ) ∂ξ 2 ξ=yi+θhF (yi,δi) F (y i , δ i ) ≤ V (y i ) -ha + hbσ( δ i r ) + h 2 2 F (y i , δ i ) ∂ 2 V (ξ) ∂ξ 2 ξ=yi+θhF (yi,δi) F (y i , δ i ),
1 See [START_REF] Jiang | Input-to-state stability for discrete-time nonlinear systems[END_REF] for the definition of this property and also for its equivalent Lyapunov characterizations.

where the estimates from [START_REF] Bernuau | On homogeneity and its application in sliding mode[END_REF] and the second property from Assumption 2 are used on the last step. Denote

w = sup y∈Sr(1) sup δ r ≤σ -1 ( a 2b ) sup ι∈[0,1] F (y, δ) ∂ 2 V (ξ) ∂ξ 2
ξ=y+ιF (y,δ)

F (y, δ), obviously w ∈ (0, +∞) (recall that sup y∈Sr(1) f (y) < +∞ and F (y, δ) ≤ F (y, 0) + σ( δ r ) for all y ∈ S r and δ ∈ R m by Assumption 2), and assume that

δ i r ≤ σ -1 ( a 2b ), (11)
then

V (y i + hF (y i , δ i )) -V (y i ) ≤ -ha + hbσ( δ i r ) + h 2 2 w ≤ -h a 2 + h 2 2 w provided that h ≤ 1. In addition, if h ∈ (0, h 0 ] with h 0 = min{1, a 2w }, then V (x i+1 ) -V (x i ) ≤ -h a 4 x i µ r
under the restriction [START_REF] Bhat | Geometric homogeneity with applications to finite-time stability[END_REF], which can be equivalently rewritten as

x i r ≥ d r,∞ σ -1 ( a 2b ) or V (x i ) ≥ c 2 d r,∞ σ -1 ( a 2b ) µ .
Therefore, the system (9) is ISS [START_REF] Jiang | Input-to-state stability for discrete-time nonlinear systems[END_REF], and any sequence {x i } ∞ i=0 generated by ( 9) reaches for the zone

c 1 x i µ r ≤ V (x i ) ≤ c 2 d r,∞ σ -1 ( a 2b ) µ ,
which is desired.

For the system (3) the analysis is similar:

V (x i+1 ) -V (x i ) = V (x i+1 ) -V (x i+1 -hF (x i+1 , d i+1 )) = x i+1 µ r [V (y i+1 ) -V (y i+1 -hF (y i+1 , δ i+1 ))],
where x i+1 = Λ r ( x i+1 r )y i+1 for y i+1 ∈ S r (1) and d i+1 = Λ r ( x i+1 r )δ i+1 for some δ i+1 ∈ R m . Next, all the rest steps coincide.

Therefore, the schemes (2) and (3) (or ( 9) and ( 10)) for approximation of solutions of a homogeneous system (8) keep unchanged the convergence rates and stability margins in the presence of perturbations.

Remark 16 Note that for an input d(t) = 0 the corresponding trajectory may visit the origin, but does not stay there due to disturbance presence. In such a case the algorithms (9) and (10) will decrease to zero the amplitude of the time step for the systems with negative degree or augment it till infinity for the systems with positive degree following their construction. In order to avoid this issue, the upper and lower limits 0 < h < h < +∞ on

h x i -ν r or h x i+1 -ν
r must be imposed in (9) and (10), respectively.

Example

Consider the super-twisting benchmark system:

ẋ1 = -k 1 |x 1 | 0.5 sign(x 1 ) + x 2 + d, ẋ2 = -k 2 sign(x 1 ) + f, where x = [x 1 x 2 ] ∈ R 2 is
the state vector, k 1 > 0 and k 2 > 0 are the system parameters; d ∈ R is the external disturbance and f ∈ R is the matched perturbation that has a known upper bound with k 2 > sup t≥0 |f (t)|; ν = -0.5 is the homogeneity degree for r = [1 0.5] . The signal f can be embedded in the system equations by considering instead a differential inclusion, and it is easy to check that then the super-twisting algorithm is an example of the system (8), which verifies assumptions 1 and 2 (for r = 0.5 and σ(s) = s). Since it is a discontinuous system, then we are going to compare for this model the explicit method (9) and the same approach without state-dependent scaling:

x i+1 = x i + hF (x i , d i ), t i = ih (12)
for i = 0, 1, . . . , which was previously studied in many works and it is recommended for the use in the sliding mode control community [START_REF] Levant | Homogeneity approach to high-order sliding mode design[END_REF][START_REF] Efimov | Realization and discretization of asymptotically stable homogeneous systems[END_REF][START_REF] Levant | Stability and robustness of homogeneous differential inclusions[END_REF][START_REF] Levant | On fixed and finite time stability in sliding mode control[END_REF] (a more sophisticated solution is given in [START_REF] Livne | Proper discretization of homogeneous differentiators[END_REF]).

For simulations k 1 = 8 and k 2 = 5 are selected with

f (t) = k 2 3 (1 + sin(2t 2 )),
and the run of the algorithms (9) (for i = 0, . . . , N with the maximum number of iterations N = 10 4 ) and ( 12) have been performed for the same initial condition x 0 = [10 10] with the discretization step h = 10 -2 and h = h -1 = 10 10 from Remark 16. The results of simulation for the case d(t) = 0 are presented in Fig. 1, where x(t) r is plotted in a logarithmic scale, the result of (9) corresponds to the blue line, the result of ( 12) is shown by the red one. As we can conclude from these simulations, the method (9) outperforms [START_REF] Brogliato | Globally stable implicit euler time-discretization of a nonlinear single-input sliding-mode control system[END_REF] since the former is finite-time converging, while the latter one cannot converge to the origin. It is important to highlight that these discretization algorithms have a similar computational complexity (comparing with [START_REF] Brogliato | Globally stable implicit euler time-discretization of a nonlinear single-input sliding-mode control system[END_REF], in [START_REF] Bernuau | Verification of ISS, iISS and IOSS properties applying weighted homogeneity[END_REF] it is necessary to additionally compute the value of the norm x i ν r , which only needs the power and the addition operations). The results of simulation for d(t) = 1 + cos(10t) are plotted in Fig. 2, where (9) is approaching the origin much closer, while [START_REF] Brogliato | Globally stable implicit euler time-discretization of a nonlinear single-input sliding-mode control system[END_REF] generates bounded trajectories loosing the precision in a vicinity of zero. All these results confirm the theoretical finding of the paper.

Conclusions

In this work a series of results has been established devoted to application of the explicit and implicit Euler methods with state-dependent scaling of the discretization step for discrete-time approximation of solutions of homogeneous systems. The main contributions can be summarized as follows:

• The obtained approximations can be made globally converging for asymptotically stable homogeneous systems provided that the discretization step is selected sufficiently small, in addition they preserve finite-time or fixed-time rates of convergence in infinite number of iterations (theorems 9 and 12). The maximal admissible discretization step can be evaluated on the sphere (Remark 10). • For any value of the discretization step, the implicit and explicit Euler methods provide a good approximation of the system solutions (Theorem 13). • Both approaches, (2) and (3) (( 9) and ( 10)), preserve the input-to-state stability of the homogeneous dyn- [START_REF] Bernuau | On extension of homogeneity notion for differential inclusions[END_REF] with the same asymptotic gain (Theorem 15).

ξ

  + hf (ξ) r since f (ξ) is bounded for any ξ ∈ S r (1).

Fig. 1 .Fig. 2 .

 12 Fig. 1. The results of simulation for d(t) = 0

  where R is the set of real number.• |•| denotes the absolute value in R, for a vector x ∈ R n the symbol x denotes the Euclidean norm, and the corresponding induced matrix norm A for a matrix A ∈ R n×n . • A continuous function α : R + → R + belongs to the class K if α(0) = 0 and the function is strictly increasing. The function α : R + → R + belongs to the class K ∞ if α ∈ K and it is increasing to infinity. • The identity matrix of dimension n × n is denoted as I n , and diag{r i } n i=1 is a diagonal matrix with the elements on the main diagonal equal r i .• A sequence of integers 1, 2, ..., n is denoted by 1, n.


Abbas Abou Hamdan

Sébastien Dementin

Pierre-Pol Liebgott

Oscar Gutiérrez-Sanz

Pierre Richaud

Antonio L De Lacey

Marc Rousset

Patrick Bertrand

Laurent Cournac

Christophe Léger

Understanding and Tuning the Catalytic Bias of Hydrogenase

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Understanding and tuning the catalytic bias of hydrogenase

Abbas Abou Hamdan, a Sébastien Dementin, a Pierre-Pol Liebgott, a Oscar Gutierrez-Sanz, b Pierre Richaud, c ABSTRACT: When enzymes are optimized for biotechnological purposes, the goal often is to increase stability or catalytic efficiency. However, many enzymes reversibly convert their substrate and product, and if one is interested in catalysis in only one direction, it may be necessary to prevent the reverse reaction. In other cases, reversibility may be advantageous because only an enzyme that can operate in both directions can turn-over at a high rate even under conditions of low thermodynamic driving force. Therefore, understanding the basic mechanisms of reversibility in complex enzymes should help the rational engineering of these proteins. Here, we focus on NiFe hydrogenase, an enzyme that catalyses H2 oxidation and production, and we elucidate the mechanism that governs the catalytic bias (the ratio of maximal rates in the two directions). Unexpectedly, we found that this bias is not mainly determined by redox properties of the active site, but rather by steps which occur on sites of the proteins that are remote from the active site. We evidence a novel strategy for tuning the catalytic bias of an oxidoreductase, which consists in modulating the rate of a step that is limiting only in one direction of the reaction, without modifying the properties of the active site.

The four Michaelis parameters (two maximal rates and two values of Km) which characterize the forward and reverse reactions of a one-substrate one-product enzyme are related to each other and to the equilibrium constant of the reaction by the Haldane equation. 1 The forward and reverse maximal rates sometimes differ so much that certain enzymes were designated as "one-way enzymes". 2 The origin of such kinetic asymmetry, referred to as "catalytic bias", has rarely been investigated. Jencks proposed that directionality may result from the destabilization of the enzyme-substrate complex, which would decrease the energy required to reach the transition state in the forward direction. 2 This ""Circe effect"" is controversial 3 and has found no echo in the case of oxidoreductases, whose directionality is always discussed by comparing the potential of the substrate/product redox couple with the potential of either the active site or the redox centers of the electron transfer (ET) chain, when there is one (see examples in supplementary information).

In trying to explain the catalytic bias from a single property of the enzyme (the potential of either the active site or the ET chain), one implicitly assumes that a single redox step, ET either between the substrate and the active site or to/from an electron relay, determines both maximal rates. However, the catalytic cycle of oxidoreductases involves various steps (substrate binding, product release, proton and electron transfers, active-site chemistry) and it may occur that the rate limiting step (rls) is not the same when the enzyme works forward or backward. Two different steps may define the two maximal rates, and their ratio. Demonstrating that this can occur requires that the rls be defined in both directions in a series of variants that exhibit different catalytic preferences. Hereafter, we do so by characterizing a series of Desulfovibrio fructosovorans (Df) NiFe hydrogenase mutants. Previously, we have shown that the WT enzyme catalyzes H2 production and oxidation at similar maximal rates, and that narrowing the substrate channel using site directed mutagenesis has no effect on the maximal rate of H2 oxidation. 4,5 Here we show that these mutations slow H2 production up to 100-fold. In redox titrations, the active site of the mutants that have little reductive activity behaves like in the WT enzyme. We use a novel method based on the isotope-exchange assay to determine the rates of H2 release from the active site to the solvent, and we conclude that this step limits H2 production whereas H2 entry does not determine the maximal rate of H2 oxidation. Conversely, the previous observation 5 that modifying the electron transfer chain selectively slows H2 oxidation shows that electron transfer limits the rate of H2 oxidation but not H2 production. 6 This is the first demonstration, on a specific example, that slowing a step that is rate limiting only when the enzyme works in one direction is a general mechanism for biasing the enzyme in the other direction, independently of the redox properties of the cofactors.

We prepared the enzyme samples as described previously. 5 H2 oxidation rates were measured using a spectrophotometric assay at pH 8, 30°C, with 50mM oxidized methyl viologen (MV), under one atm. of H2. 5 The Km for H2 were measured electrochemically. 5 The voltammograms in fig 2 were obtained with the enzyme bound to a rotating graphite electrode. 7 H2 production was followed using mass spectrometry (MS), in a solution initially saturated with Ar, with a saturating concentration of re-duced MV (2mM) at pH 7.2, 30°C; since H2 inhibits H2production by NiFe hydrogenases, 8 we extrapolated the initial rates using a plot of 1/rate against 1/[H2]. We performed the H + /D + exchange experiments at pH 7.2, 30°C, 4 and FTIR titrations at pH 8. 9 Previous investigations using crystallography and molecular dynamics have predicted that a network of hydrophobic channels guides H2 to/from the buried active site of Df NiFe hydrogenase. 10 Substituting Leu122 and Val74, which shape a bottleneck in the channel (fig. 1a), strongly decreases the rates of intramolecular transport of H2, CO and O2 in both directions. 4,5 The decrease of the H2-binding bimolecular rate constant has no effect on the maximal rate of H2 oxidation obtained by extrapolation to infinite concentration of H2. 5 Furthermore, the mutations considered below moderately increase the Michaelis constant for H2 (5 to 30-fold for L122F-V74I and V74M, respectively), but the Michaelis constant is small in the WT (about 10 matm of H2) and these mutations decrease less than 2-fold the H2 oxidation rate under one atm. of H2. 5 Hence, substrate binding does not limit H2 oxidation under 1 atm. of H2 in the WT and in these mutants. Figure 2 compares electrochemical signals obtained with the WT enzyme and two mutants (L122F-V74I and V74M, pdb 3CUS and 3H3X) attached to rotating-disc graphite electrodes. 7,8,11 Hydrogen oxidation and production are detected as positive and negative currents in a single experiment where the electrode potential is swept across a wide range. The absolute magnitude of the current is meaningless because it is proportional to turnover rate times the unknown electroactive coverage, but figure 2 clearly shows that the mutations decrease the ratio of oxidation over reduction currents: unlike the WT, the mutants are biased towards H2 oxidation. This illustrates that the term "bias" does not refer to "the direction of the reaction", which is imposed by thermodynamics, but to the ratio of rates (here, currents) measured for the same reaction proceeding in opposite directions. Since the conditions chosen for the two measurements are different, the ratio of rates is not an equilibrium constant. 12 In electrochemical experiments, thermodynamics only forces the "open circuit potential" (OCP) to equate the reduction potential of the H + /H2 couple given by the Nernst equation; it is indeed independent of the enzyme (fig 2). 8,11 The results of solution assays with oxidized or reduced methyl viologen (MV) quantify the bias of the mutants. Figure 3a shows the maximal rates of H2 against the maximal rates of H2 production. The mutants have the same oxidation activity as the WT (only V74C is slightly more active 9 ) whereas their maximal H2 production rates vary by two orders of magnitude. The ratio of maximal rates (oxidation over production) ranges from 2.5 for the WT to 200 for the V74M mutant. Note that the exact val-ue of this bias is not an intrinsic property of the enzyme: it depends on which reaction is catalyzed (it is not the same for hydrogenase exchanging electrons with MV, an electrode, or its physiological partner), and on the conditions that are chosen to measure the two rates.

Table 1. Redox properties of the WT enzyme and two mutants.

Ni-A, Ni-B and Ni-SU are inactive forms of the enzyme, whereas Ni-SI, Ni-C and Ni-R are catalytic intermediates. 13 All reduction potentials are in mV vs SHE. Typical errors are ± 15mV. To examine whether the decrease in H2 production activity was due to a mutation-induced increase in the reduction potential of the active site, we compared the redox properties of the active site in the WT enzyme and the V74C and V74M mutants. Redox titrations of the active site states of wild type hydrogenase and the mutants V74M and V74C were monitored by FTIR in a spectroelectrochemical cell as described before. 9 The FTIR signatures of these mutants are similar to those of the WT enzyme. Supplementary fig. S1 shows the redox titration of the CO band of the Ni-C state for the three hydrogenases, from which the reduction potentials of the active site catalytic intermediates could be determined. Table 1 shows the formal reduction potentials of all the active site states. All results are similar, which suggests that the low reductive activities of the mutants are not the consequence of the potential of the active site being shifted upward compared to the WT value.

Enzyme

In search of the rate limiting step of H2 production, we determined the rate constants of H2 release, under the same conditions as H2 production, using the isotope exchange assay, where the enzyme transforms D2 into HD and eventually H2 using protons from the solvent. 4 The reaction is followed using mass spectrometry (fig 4). Assuming that all diatomic molecules diffuse to the active site with a bimolecular rate constant kin and from the active site with a 1st-order rate constant kout, and that H + /D + exchange at the active site proceeds with a 1storder rate constant k, we predicted 3 that both the concentration of D2 and the isotope content T=[D2]+[HD]/2 should decrease exponentially with time, as observed in figure 4, with rates kD and kT, respectively, which depend on kin, kout, k, e0 (the enzyme concentration) and c0 (the initial concentration of D2). The analysis that we previously reported allowed the determination of only the ratio kout/k (ref 4). We now show in supplementary information that under certain conditions, which apply to the case of Df NiFe hydrogenase, the rate of H2 release kout can be deduced from the data using eq. 1:

When kout/k is small, little HD is released, kD≈kT (fig 4b), and kout ≈kDc0/e0. We used this method with all above mentioned mutants to determine the rates of release of H2 (kout), which we plotted against the maximal rate of H2 production in fig 3b (note the log-log scale). All data points are close to the y=x line, demonstrating that H2 production is limited by H2 release from the active site: the mutants have little production activity because the channel is obstructed, and H2-diffusion also determines the rate of H2production in the WT enzyme.

Conversely, we previously made the same enzyme a better catalyst of H2 production than H2 oxidation by modifying its ET chain. 6 The enzyme has a series of three FeS clusters, which "wires" the active site to the enzyme's redox partner (fig 1). The surface-exposed (distal from the active site) [4Fe4S] cluster is coordinated by His184; replacement with Gly184 or Cys184 slows H2 oxidation 33 and 60-fold, respectively, whereas the rates of H2 production decrease only 1.3 and 2.1-fold (Table 1 in ref 6). These mutants are biased towards H2 production: both mutations decrease the ratios of maximal rates (oxidation over production) by a factor of about 27. The H184G mutant has an open coordination site on one Fe ion of the distal cluster; when exogenous imidazole binds to this cluster, the H2-oxidation rate increases from 3% to 33% of that of the WT, whereas the H2-production rate is not affected 6 : binding of imidazole to the distal cluster of H184G partially repairs the ET chain and restores the native catalytic bias. It has long been suggested that ET limits the rate of H2-oxidation in WT NiFe-hydrogenase 14 and this was recently supported by measurements of ET rates in the enzyme from D. fructosovorans. 15 That a modification of the distal cluster biases the enzyme in the direction of H2 production shows that ET is rate limiting only for H2-oxidation.

Regarding the WT and mutant hydrogenases discussed here, our data show that H2-diffusion in the gas channel is the rate limiting step of H2-production but not H2oxidation, whereas the latter is limited by electron transfer. Altering one of these steps selectively affects one of the two reactions. We therefore evidence a novel mechanism for tuning the catalytic bias of an oxidoreductase which is independent of the properties of the active site, but requires that the rls be different under the two different sets of conditions that are used to drive catalysis in one direction and the reverse. This contrasts with all previous explanations of the catalytic bias of oxidoreductases (see examples in supplementary information) which proposed that it may be determined by the redox properties of the cofactor(s).

The above results, and previous evidences that enzymes" turnover rates may be limited by steps other than active site chemistry (e.g. proton transfer, electron transfer, substrate release or lid-opening 16 ), emphasize the need to study all steps of the reaction rather than only active site chemistry. Identifying the rls in a catalytic reaction, keeping in mind that active site chemistry may be fast, is a difficult task. Yet this is a prerequisite if one aims at understanding global kinetic properties, such as rate enhancement, proficiency, or bias.

Hydrogenases could be used as H2 oxidation or production catalysts in biotechnological devices 17 if the enzymes that can be produced in large amounts were not inhibited by O2. This has motivated research on the inhibition mechanism, and Df NiFe hydrogenase mutants that proved more resistant to O2 than the WT enzyme have been characterized; this includes the V74C and V74M mutants studied here. 9,18 That mutations which increase O2 tolerance may also change the catalytic bias of the enzyme will have to be considered in studies which aim at optimizing this biological catalyst.

ASSOCIATED CONTENT

Complete references 4-6,9,15,18 ; examples of discussions in the literature of what may determine the bias of oxidoreductases ; demonstration of eq 1 ; redox titration curves of the WT and the V74C and V74M mutants. This material is available free of charge via the Internet at http://pubs.acs.org.

  Abbas Abou Hamdan, a Sébastien Dementin, a Pierre-Pol Liebgott, a Oscar Gutierrez-Sanz, b Pierre Richaud, c Antonio L. De Lacey, b Marc Rousset, a Patrick Bertrand, a Laurent Cournac, c * and Christophe Léger. a * (a) Laboratoire de Bioénergétique et Ingénierie des Protéines, Institut de Microbiologie de la Méditerranée, CNRS and Aix-Marseille Université. 31 chemin Joseph Aiguier, 13402 Marseille Cedex 20, France. (b) Instituto de Catalisis y Petroleoquimica (CSIC), 28049 Madrid, Spain (c) CEA, DSV, IBEB, Laboratoire de Bioénergétique et Biotechnologie des Bactéries & Microalgues, Saint Paul Lez Durance / CNRS, UMR Biologie Végétale & Microbiologie Environnementales, Saint Paul lez Durance, France / Aix-Marseille Université, Saint Paul lez Durance, France.

Figure 1 .

 1 Figure 1. The structure of Df NiFe hydrogenase (pdb accession code: 1YQW). Close-ups show the structural elements that are involved in the catalytic cycle: the gas channel (a), the active site (b), a chain of aminoacids that are putatively involved in proton transfer (c), and the FeS clusters that wire the active site to the redox partner (d).

Figure 2 .

 2 Figure 2. Bidirectional electrocatalysis by the WT (red), L122F-V74I (blue) and V74M (green) forms of Df NiFe hydrogenase. The conditions (10% H2, pH 5.5, 40°C, scan rate 10mV/s, electrode rotation rate 3krpm) ensure that H2 oxidation (at high electrode potential) and production (at low potential) are detected in the same experiment. At the open circuit potential (OCP), the rates of oxidation and reduction exactly cancel each other. The blank was recorded with no adsorbed enzyme.

Figure 3 .

 3 Figure 3. Kinetic properties of the hydrogenase channel mutants. Panel a shows the maximal rate of H2 oxidation plotted against the maximal rate of H2 production for 6 forms of the enzyme, and demonstrates that the mutants are biased towards H2 oxidation. Panel b shows the rates of H2 release (measured using eq 1 from experiments such as those in fig 4) against the maximal rates of H2 production, and demonstrates that the former limits the latter.

Ni

Figure 4 .

 4 Figure 4. Isotope-exchange assay of the WT (Panel a) and L122F-V74I mutant (Panel b). The changes in concentrations were used to determine kout using eq. 1

Funding Sources

Our work is funded by the CNRS, CEA, Aix-Marseille Université, ANR, City of Marseilles and Région Provence Alpes Cote d'Azur (PACA) and the Spanish Ministerio de Ciencia (Project CTQ2009-12649). We acknowledge support from the "Pole de compétitivité Capénergie". AAH thanks the CNRS and Region PACA for funding his PhD.

AUTHOR INFORMATION

Corresponding Author christophe.leger@imm.cnrs-mrs.fr, laurent.cournac@ird.fr

Present Addresses

Laurent Cournac: IRD, UMR Eco&Sols, 2 Place Viala, 34060 Montpellier Cedex 02, France.

Author Contributions

AAH, SD and PPL purified the proteins. AAH carried out the electrochemical experiments. AAH, SD, PR and LC carried out all solution assays and analyzed the data. OG and AdL carried out the spectro-electrochemical titrations. PB, LC, MR and CL initiated research. PB, SD and CL wrote the paper. All authors have given approval to the final version of the manuscript.