Ludovic Chaput

Edithe Selwa

Eddy Elisée

Bogdan I Iorga
email: bogdan.iorga@cnrs.fr

‡ Eddy Elisée

Blinded evaluation of cathepsin S inhibitors from the D3RGC3 dataset using molecular docking and free energy calculations

Keywords: docking, ranking, scoring function, free energy calculations, Gold, Vina, Glide, Gromacs, PMX, cathepsin S inhibitors, D3R, Drug

INTRODUCTION

Drug Design Data Resource (D3R, https://drugdesigndata.org/) organizes, on a regular basis, blinded prediction challenges with the aim to evaluate the performance of tools and protocols that are used in real-life computer-aided drug discovery projects. To achieve this, datasets presenting different levels of difficulty are presented to the community, which is asked to predict, in "blind" conditions, the binding modes and the relative affinities of compounds.

The D3R Grand Challenge 3, which was held in 2017, was focused on a single protein, cathepsin S (CatS, Figure 1), a target with multiple potential applications that has received much attention during the recent years [START_REF] Thurmond | Cathepsin S inhibitors as novel immunomodulators[END_REF][START_REF] Link | Advances in cathepsin S inhibitor design[END_REF][START_REF] Wiener | Recent advances in the design of cathepsin S inhibitors[END_REF][START_REF] Lee-Dutra | Cathepsin S inhibitors: 2004-2010[END_REF][START_REF] Wilkinson | Cathepsin S: therapeutic, diagnostic, and prognostic potential[END_REF]. In Phase 1A the participants were asked to predict the crystallographic poses of 24 CatS ligands, the affinity ranking for 136 CatS ligands and the relative binding affinities for a designated free energy subset of 33 compounds. In phase 1B, the organizers released the corresponding receptor structures for the 24 CatS ligands, and asked to predict the crystallographic poses using the released receptor structures. In Phase 2 the participants were required to provide the same affinity predictions for the subset of 136 ligands and relative free energy for the subset of 33 compounds as in Phase 1, but taking into account the additional structural data (24 new protein-ligand complexes) released at the end of Phase 1B. The CatS dataset containing the 136 ligands used for ranking prediction can be found in the Electronic Supplementary Material (Figure S3).

Figure 2 shows the chemical structures of the 24 compounds from CatS dataset for which the pose predictions were required. Additionally, the participants were asked to predict the relative affinities for one homogeneous subset of compounds that is suited for free energy calculations.

The structures of the 33 CatS inhibitors included in the free energy subset are presented in Figure 3.

METHODS

Protein structures. We found 38 crystal structures of CatS available in the Protein Data Bank (PDB) [START_REF] Berman | The Protein Data Bank[END_REF] from which 31 were human CatS and 7 were mouse CatS (see Table S1 in the Electronic Supplementary Material for the complete list). These structures were superposed, showing no conformational variability for the backbone and for most sidechains. The only exceptions were the sidechains of residues Lys64, Arg141 and Asn161, which are very flexible, adopting multiple conformations, and of residues Phe146 and Phe211, which adopt two major orientations each. The residues Lys64 and Arg141 are relatively far from the binding site and therefore their flexibility was not taken into account for the docking process. In these conditions, we generated four representative conformations of CatS (named CatS_conf1 to CatS_conf4) starting from the structure 1GLO (apo) and adjusting the orientations of Phe146 and Phe211. These conformers, with the ions and solvent molecules manually removed, were superimposed on the reference structure CatS (containing the SO4 ion) provided by the D3R Grand Challenge 3 organizers, in order to conserve the same coordinate system through the whole process.

Hydrogen atoms were added using Hermes, the graphical interface of Gold v5.2.2 software [START_REF] Verdonk | Improved proteinligand docking using GOLD[END_REF],

or with Maestro, the graphical interface of Schrödinger software (http://www.schrodinger.com), prior to docking.

Ligands. A training set of CatS ligands was constituted by gathering from ChEMBL, PubChem

Assay and BindingDB databases the compounds with known activity (IC50) for CatS (1591 unique compounds). The structures of these compounds were retrieved in SMILES format and were converted into three-dimensional MOL2 files using Ligprep v37014 (Schrödinger, http://www.schrodinger.com). The same conversion from SMILES into MOL2 formats was applied to the D3RGC3 dataset of CatS inhibitors provided by organizers. In all cases, the protonation state for the compounds was adjusted at physiological pH using LigPrep (Schrödinger, http://www.schrodinger.com).

Preliminary docking. In the preliminary analysis step, several docking software and scoring functions have been tested for their ability to predict relative affinity ranking and reproduce the protein-ligand complexes: Gold [START_REF] Verdonk | Improved proteinligand docking using GOLD[END_REF] with the GoldScore, ChemScore, ChemPLP and ASP scoring functions, Vina [START_REF] Trott | AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading[END_REF] and Glide (Schrödinger, http://www.schrodinger.com). Default parameters were used in all cases for docking, except with Gold, where a search efficiency of 200% was used in order to better explore the conformational space. For Glide, the standard precision (SP) mode was used. The binding sites were considered with Gold as spheres with a 20 Å radius around a virtual point with coordinates (-7,12,-8). With Vina, the binding sites were defined as a 25 x 25 x 25 Å 3 cube centered on the same virtual point. With Glide, the binding site was defined with a inner box of 13 x 12 x 13 Å 3 and an outer box of 30 x 29 x 30 Å 3 . Glide attempts to position the ligand center within the inner box volume, while the outer box is defined to contain all the ligand atoms. The docking calculations were carried out with the four CatS conformers described above and with the two structures (containing a DMSO molecule and a SO4 ion, respectively) provided by the organizers. All structures were rigid during the docking process.

Ten docking poses were generated for each ligand from the training set and from the CatS D3RGC3 dataset.

Phase 1A docking. From the preliminary analysis, Gold with the ASP scoring function were identified as the most adapted for the given system, and therefore were used for Phase 1 predictions. The conformation CatS_conf4 gave the best results in the preliminary docking step (data not shown) and was used for the docking of most compounds from the CatS D3RGC3 dataset. Notable exceptions are: i) for compounds CatS_4 and CatS_6, which contain a different scaffold, the structure 3MPF was used; ii) for compound CatS_14, the structure provided by organizers containing the DMSO molecule was used; iii) for compounds CatS_2, CatS_17, CatS_20, CatS_22, CatS_23 and CatS_24, the structure provided by organizers containing the SO4 ion was used. The two flip/flop forms of the amide side chain from Asn161 were considered, and then best pose (in terms of score or RMSD) selected. The same Gold docking parameters as described above were used, except that 100 docking poses were generated for each ligand. In the post-docking processing procedure, docking complexes were used as input structures for 10 ns molecular dynamics simulations in explicit solvent using Gromacs and AMBER 99SB force field. During the simulation, constraints were applied on the protein backbone and on the ligand core fragment that was also constrained during docking. The first 1ns was discarded and the remaining snapshots clusterized (using the g_cluster tool available in Gromacs) in order to provide a single cluster. The representative structure of this cluster was selected for submission (1 structure per compound).

Phase 1B docking. After the release of the 24 crystal structures corresponding to the compounds CatS_1 to CatS_24 (but without the coordinates of the ligands), we carried out docking calculations for each of these structures with its native ligand using the same conditions as in Phase 1A (Gold with the ASP scoring function, 100 docking poses generated for each ligand).

Phase 2 docking. Our experience from Phase 1A using Gold/ASP with scaffold match constraints showed that the atoms with constraints applied were not always well aligned with the reference structure, even when very strong constraints are used (increase of the weight parameter, that determines how closely ligand atoms fit onto the scaffold, up to 2000). Therefore, in Phase 2 we used the Glide/SP protocol described above, with the exception of scaffold constraints that were applied on the common substructure of ligands with either CatS_15, CatS_20 or CatS_24. The best Glide/SP docking scores were selected for Phase 2 ranking prediction submission (1 per compound). The related docking poses for the 33 CatS compounds from the free energy dataset were used as the starting points for following free energy calculations.

Free energy calculations. The complexes selected in the previous step represented the input for relative alchemical free energies calculations, using the protocol described in a previous study [START_REF] Selwa | Blinded evaluation of farnesoid X receptor (FXR) ligands binding using molecular docking and free energy calculations[END_REF]. Hybrid structures and topologies were built using a modified version of the MOL2FF package developed in our team. CatS_67 was used as reference structure for the free energy calculations. Briefly, equilibrium 10 ns MD simulations were performed for the two states (corresponding to lambda 0 and 1), using Gromacs [START_REF] Pronk | GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit[END_REF] and OPLS-AA force field [START_REF] Kaminski | Evaluation and Reparametrization of the OPLS-AA Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides[END_REF][START_REF] Robertson | Improved Peptide and Protein Torsional Energetics with the OPLSAA Force Field[END_REF].

Snapshots from the equilibrium runs were extracted to spawn 100 simulations of 50 ps each to alchemically morph between the two states of the system. The work values over every nonequilibrium transition were extracted using some scripts from the PMX software [START_REF] Gapsys | pmx: Automated protein structure and topology generation for alchemical perturbations[END_REF][START_REF] Gapsys | Calculation of binding free energies[END_REF][START_REF] Gapsys | Accurate and Rigorous Prediction of the Changes in Protein Free Energies in a Large-Scale Mutation Scan[END_REF] in order to estimate the free energy differences relying on the Crooks Fluctuation Theorem and utilizing Crooks-Gaussian Intersection as estimator. As described previously [START_REF] Selwa | Blinded evaluation of farnesoid X receptor (FXR) ligands binding using molecular docking and free energy calculations[END_REF], when the two ligands considered for the alchemical transformation had the same charge, separate calculations in explicit solvent were carried out for the transformation of the protein-ligand complex and of the ligand alone. The relative free energy of binding is then calculated as the difference between the corresponding work for these two transformations. When the two ligands had different charges, a single box containing the protein-ligand complex and the ligand alone, separated by 30 Å, was used. The ligand from the complex was converted in one direction, whereas the ligand alone was converted simultaneously in the opposite direction. Thus, the overall charge of the system is conserved during the whole simulation, and the global work for this system represents the relative free energy of binding between ligands. Submission preparation. Protein structures were converted into PDB format for submission using UCSF Chimera [START_REF] Pettersen | UCSF Chimera --a visualization system for exploratory research and analysis[END_REF], and the docking poses were converted into MOL format using CORINA v3.60 (http://www.molecular-networks.com/), as the MOL format corresponds to the SDF output format in CORINA.

Graphics. Chemical structures were depicted using CACTVS Chemoinformatics Toolkit v3.409 (Xemistry, http://www.xemistry.com/), images for protein structures were generated using PyMol 1.8.1 (Schrödinger, http://www.pymol.org/).

Statistics. Statistics were computed using the R-package (http://www.r-project.org).

Chemoinformatics. Tanimoto similarities were computed using Canvas (Schrödinger).

RESULTS AND DISCUSSION

During the last years we developed a protocol for docking and virtual screening that proved to be relatively successful in our participation to different blind prediction challenges: SAMPL3 (2011) [START_REF] Surpateanu | Evaluation of docking performance in a blinded virtual screening of fragment-like trypsin inhibitors[END_REF], SAMPL4 (2013) [START_REF] Colas | Virtual screening of the SAMPL4 blinded HIV integrase inhibitors dataset[END_REF], CSAR (2014) [START_REF] Martiny | Blind pose prediction, scoring, and affinity ranking of the CSAR 2014 dataset[END_REF], D3R Grand Challenge (2015) [START_REF] Selwa | Molecular docking performance evaluated on the D3R Grand Challenge 2015 drug-like ligand datasets[END_REF] and D3R Grand Challenge 2 (2016) [START_REF] Selwa | Blinded evaluation of farnesoid X receptor (FXR) ligands binding using molecular docking and free energy calculations[END_REF]. This protocol involves two steps, the first one representing a preliminary analysis of publicly available structural and biochemical data in order to identify the most appropriate docking software and scoring function for the system of interest. In the second step, we use these parameters for the prediction of binding modes (pose prediction) and of relative affinities of ligands (scoring). Enhanced genetic algorithm parameters for docking (a search efficiency of 200%) proved to be beneficial in our previous studies [START_REF] Surpateanu | Evaluation of docking performance in a blinded virtual screening of fragment-like trypsin inhibitors[END_REF][START_REF] Colas | Virtual screening of the SAMPL4 blinded HIV integrase inhibitors dataset[END_REF][START_REF] Martiny | Blind pose prediction, scoring, and affinity ranking of the CSAR 2014 dataset[END_REF][START_REF] Selwa | Molecular docking performance evaluated on the D3R Grand Challenge 2015 drug-like ligand datasets[END_REF][START_REF] Selwa | Blinded evaluation of farnesoid X receptor (FXR) ligands binding using molecular docking and free energy calculations[END_REF], especially for big and flexible ligands. Therefore, the same parameters were used in this work to explore as exhaustively as possible the docking conformations of ligands from the CatS D3RGC3 dataset.

Preliminary analysis

The superposition of the 38 crystal structures of CatS identified in the PDB (Table S1) showed very little conformational variability on backbone and on most side chains, with the exception of a few residues shown in Figure 4. Among these structures, 19 contain a ligand in the active site, and 4 of these ligands (599, 23Z, N2A and N2D from the structures 3IEJ, 3KWN, 3MPE and 3MPF, respectively) share a common scaffold with some compounds from the CatS D3RGC3 dataset. The three-dimensional structure of protein in these structures is well conserved, with the exception of two side chains of Phe211 and Phe146 that adopted alternative orientations, as found in 3 representative structures (PDB codes 2F1G, 3N4C, 3MPF). Therefore, we used an apo structure of CatS (PDB code 1GLO) to build 4 conformers, representing the four combinations of representative conformations for Phe211 and Phe146 side chains. In addition, 2 structures (containing a DMSO molecule and a SO4 ion) were provided by the D3RGC3 organizers to take into account specific ligand-DMSO or ligand-SO4 interactions. These six structures were used in the preliminary docking calculations. Gold with the GoldScore, ChemScore, ChemPLP and ASP scoring functions, Vina and Glide.

The analysis of the docking results showed that no docking program was able to provide a good prediction for the relative affinity ranking. The Spearman's and Kendall's correlation coefficients between docking scores and pIC50 values of ligands from the training set are reported in Table S2.

Surprisingly, the Spearman's and Kendall's correlation coefficients (Table S2) seem to indicate in some cases an anti-correlation between docking scores and pIC50. Indeed, the scatter plot Score versus pIC50 (Figure S2) seems to indicate the existence of two clusters of compounds (two distinct clouds of points). To investigate this observation, we calculated the distance matrix between the compounds from the training set, based on the molprint 2D descriptor, the Tanimoto distance (Canvas, Schrödinger), and a hierarchical cluster analysis (R).

These results highlighted the existence of a large cluster of compounds (colored in blue) that tend to have favorable docking scores despite low to moderate pIC50 (Figures S1 andS2, red/blue dendrogram and scatterplot). This blue cluster can explain in part some anti-correlations observed between the docking scores and pIC50 values (Table S2). Further calculation of Tanimoto distances between the compounds from the CatS D3RGC3 and the training set identified 292 compounds with a Tanimoto score >0.5 compared with at least one of the 141 CatS D3RGC3 compounds, and all these compounds belong to the blue cluster (Figures S1 andS2). The Spearman's and Kendall's correlation coefficients for the 292 compounds from the blue cluster (Table S2) show that Vina provided the best results, with correlation coefficients between 0.29 and 0.47 (Spearman's) and between 0.21 and 0.34 (Kendall's), depending on the CatS conformation used for docking. Other docking programs and scoring functions did not provide good correlations for the 292 compounds from the blue cluster.

The ability of the docking programs to reproduce crystallographic poses was also evaluated.

The 141 compounds from the D3RGC3 CatS dataset were docked on the four representative CatS conformations using different docking software and scoring function combinations and the RMSD were calculated between the common substructure of these compounds with the ligand 599 from the structure 3IEJ (Figure 5). The mean RMSD values for the 141 compounds from the D3RGC3 CatS dataset are reported in the Table 1, showing that Gold with the ASP scoring function clearly outperformed the other docking programs and scoring functions for all the 4 CatS conformers, with mean RMSD values between 2 and 3 Å. In these conditions, Gold with the ASP scoring function and the CatS_conf4 structure were selected for the next step. Table 1. Analysis of preliminary docking poses for the 141 compounds from the CatS dataset using four representative conformations of CatS. For all 10 docking poses of each ligand, RMSD values (Å) were calculated using the common substructure between the ligand and compound 599 (see Figure 5). The best RMSD value for each ligand was retained and the overall mean values are shown below.

CatS_conf1

CatS_conf2 CatS_conf3 CatS_conf4

Gold/Goldscore

Pose prediction and scoring (Phase 1A)

The 141 CatS ligands were docked with Gold and the ASP scoring function, generating 100 poses for each ligand. The CatS_conf4 conformer was used for docking the vast majority of compounds, with some exceptions (see the Methods section for details).

Three versions of the protocol were tested: i) one without any constraints; ii) one with constraints on the common substructure between the ligand and one of PDB ligands 599, 23Z and N2A and iii) one with the same constraints as above, followed by a 10 ns molecular

Pose prediction (Phase 1B)

In Phase 1B, we carried out docking calculations using the 24 crystal structures corresponding to the compounds CatS_1 to CatS_24 that were released at the end of Phase 1A (but without the coordinates of the ligands) with their native ligands. The results from Table 2 (entries 9 and 10)

show that the re-docking of ligands into their native structures provided only a slight improvement in the pose prediction for the "best RMSD" submission (from 4.36 Å to 3.95 Å) and a very bad pose prediction for the "best score" submission (8.49 Å instead of 5.19 Å).

Scoring and free energy calculations (Phase 2)

In phase 2, we only predicted the relative ranking for the 33 compounds of the free energy set (Table 2, entries 11 and 12). We used the Glide/SP protocol with scaffold constraints applied on the common substructure of ligands with either CatS_15, CatS_20 or CatS_24. The poses with the best docking score (1 pose per compound) were submitted for Phase 2 ranking prediction.

The same CatS complexes were also used for relative free energy calculations, which were carried out using a protocol described in a previous paper [START_REF] Selwa | Blinded evaluation of farnesoid X receptor (FXR) ligands binding using molecular docking and free energy calculations[END_REF]. The compound CatS_67 was used as reference, and the calculations were designed in order to minimize the transformations involving a change in the global charge of the ligand (Figure 7). Our previous study [START_REF] Selwa | Blinded evaluation of farnesoid X receptor (FXR) ligands binding using molecular docking and free energy calculations[END_REF] has highlighted some potential OPLS-AA force field issues involving non-optimal dihedrals associated with the sulfonamide group that are now fixed (manuscript in preparation), the simulations from this study using optimal parameters. The correlation coefficients obtained in Phase 2 are not as good as those obtained in Phase 1A, indicating that Glide is not performing well in ranking the poses from this free energy subset (Table 2, entry 11). Similarly, the free energy calculations do not improve the quality of the ranking prediction compared with the direct ranking prediction from the scoring function (Table 2, entry 12). However, much better results are obtained when considering only the cocrystallized ligands (RMSEc of 0.32±0.24 kcal/mol, Kendall's τ of 1.0±0.47, Spearman's ρ of 1.0±0.49, Pearson's r of 0.83±0.49), indicating a possible bad positioning of the peripheral groups that may affect the correct ranking of the compounds.

CONCLUSIONS

In this challenge we used a protocol involving an initial assessment of the most appropriate docking software and parameters using structural and functional information available in public databases followed by the docking of a CatS inhibitors dataset to provide a prediction for the binding modes and ranking of ligands. The size and conformational flexibility of ligands were important factors to take into account, and the docking calculations afforded reasonable overall pose predictions, which are however dependent on the specific nature of each ligand. In general, the rigid core of the ligands was well positioned, but the conformation of the flexible substituents was more difficult to predict, in spite of the high value of search efficiency used in the docking process. As we also observed during the previous D3R challenges, the correct ranking of docking poses is still a problematic issue. Post-processing of docking poses with molecular dynamics simulations in explicit solvent can significantly improve the prediction. No significant improvement in the ranking prediction, compared with the direct ranking obtained from the scoring function, was provided by free energy calculations.

Electronic Supplementary Material. The Electronic Supplementary Material contains the list of CatS crystal structures from the PDB, the chemical structures of the scoring CatS D3RGC3 dataset, and the plots showing the performance of our submissions. for the CatS D3R Grand Challenge 3 dataset. Our submission is highlighted in green. See Table 2 in the main text for details.

Figure 1 .

 1 Figure 1. Solid surface representation of a representative crystal structure of cathepsin S (PDB

Figure 2 .

 2 Figure 2. Chemical structures of the 24 CatS ligands included in Phase 1 for pose prediction.

Figure 3 .

 3 Figure 3. Chemical structures of the 33 CatS ligands included in the free energy set.

Figure 3 .

 3 Figure 3. Chemical structures of the 33 CatS ligands included in the free energy set (continued).

Figure 4 .

 4 Figure 4. All 38 CatS structures available in the PDB superimposed (cartoon representation).

Figure 5 .

 5 Figure 5. a) Ligand 599 in the active site of cathepsin S (crystal structure, PDB code 3IEJ); b)

Figure 6 .

 6 Figure 6. Comparison of our Phase 1A best ranking poses without constraints (submission

Figure 7 .

 7 Figure 7. Diagram showing the pathways used for relative free energy calculations. The

Figure S6 .

 S6 Figure S6. Performance of Phase 1B pose prediction submissions for the CatS D3R Grand Challenge 3 dataset (RMSD for pose 1, mean value over all compounds). Our submissions are colored in green. See Table2in the main text for details.

Figure S7 .

 S7 Figure S7. Performance of Phase 1A ranking prediction submissions (Kendall Tau, mean value over all compounds) for the CatS D3R Grand Challenge 3 dataset. Our submissions are highlighted in green. See Table2in the main text for details.

Figure S8 .

 S8 Figure S8. Performance of Phase 2 ranking prediction submissions (Kendall Tau, mean value over all compounds, incomplete submissions) for the CatS D3R Grand Challenge 3 dataset. Our submission is highlighted in green. See Table2in the main text for details.

Figure S9 .

 S9 Figure S9. Performance of Phase 2 free energy submissions (Kendall Tau, mean value over all compounds)for the CatS D3R Grand Challenge 3 dataset. Our submission is highlighted in green. See Table2in the main text for details.

Table S1 .

 S1 Cathepsin S structures available in the PDB

	PDB code	Resolution (Å)	Ligand	Organism
	1GLO	2.2	-	Homo sapiens
	1MS6	1.9	BLN	Homo sapiens
	1NPZ	2.0	C1P	Homo sapiens
	1NQC	1.8	C4P	Homo sapiens
	2C0Y	2.1	-	Homo sapiens
	2F1G	1.9	GNF	Homo sapiens
	2FQ9	1.65	-	Homo sapiens
	2FRA	1.9	-	Homo sapiens
	2FRQ	1.6	-	Homo sapiens
	2FT2	1.7	-	Homo sapiens
	2FUD	1.95	-	Homo sapiens
	2FYE	2.2	-	Homo sapiens
	2G6D	2.5	-	Homo sapiens
	2G7Y	2.0	-	Homo sapiens
	2H7J	1.5	-	Homo sapiens
	2HH5	1.8	GNQ	Homo sapiens
	2HHN	1.55	GNQ	Homo sapiens
	2HXZ	1.9	-	Homo sapiens
	2OP3	1.6	TF5	Homo sapiens
	2R9M	1.97	Y11	Homo sapiens
	2R9N	2.0	-	Homo sapiens
	2R9O	2.0	-	Homo sapiens
	3IEJ	2.18	599	Homo sapiens
	3KWN/5QC4	2.1	23Z	Homo sapiens
	3MPE/5QBY	2.25	N2A	Homo sapiens
	3MPF/5QBV	1.8	N2D	Homo sapiens
	3N3G	1.6	93N	Homo sapiens
	3N4C	1.9	EF3	Homo sapiens
	3OVX	1.49	O64	Homo sapiens
	4P6E	1.8	2FC	Homo sapiens
	4P6G	1.58	2FZ	Homo sapiens
	4BPV	2.0	-	Mus musculus
	4BQV	1.7	-	Mus musculus
	4BS5	1.25	MG2	Mus musculus
	4BS6	1.2	-	Mus musculus
	4BSQ	1.96	QQV	Mus musculus
	4MZO	1.47	-	Mus musculus
	4MZS	1.85	-	Mus musculus

Table S2 .

 S2 Kendall's τ and Spearman's ρ correlations for the whole training set (black) and of a subset of compounds having a Tanimoto similarity score of >0.5 with any compound from the CatS D3RGC3 dataset (blue). See also FigureS1.

		CatS_conf1		CatS_conf2		CatS_conf3		CatS_conf4
		Kendall	Spearman	Kendall	Spearman	Kendall	Spearman	Kendall	Spearman
	Gold	-0.05	-0.07	-0.07	-0.10	-0.10	-0.14	-0.06	-0.09
	(Goldscore)	0.04	-0.18	-0.13	-0.21	0.04	-0.23	0.07	-0.09
	Gold	-0.20	-0.29	-0.20	-0.28	-0.22	-0.33	-0.18	-0.26
	(Chemscore)	-0.12	0.07	-0.14	-0.20	-0.16	0.06	-0.07	0.10
	Gold (ASP)	-0.07 0.04	-0.10 0.07	-0.09 -0.05	-0.13 -0.08	-0.07 0.00	-0.10 0.00	-0.08 0.07	-0.12 0.10
	Gold	-0.09	-0.14	-0.08	-0.12	-0.12	-0.18	-0.10	-0.14
	(ChemPLP)	-0.02	-0.02	-0.11	-0.16	-0.08	-0.11	0.06	0.09
	Glide (SP)	0.11 0.12	0.17 0.17	-0.03 -0.13	-0.04 -0.19	-0.03 -0.06	-0.05 -0.10	0.07 0.02	0.10 0.03
	Vina	0.04 0.25	0.06 0.35	0.00 0.21	0.00 0.29	0.03 0.29	0.04 0.42	0.04 0.34	0.06 0.47

Acknowledgments. We thank Prof. Bert de Groot for helpful discussions. This work was supported by the Laboratory of Excellence in Research on Medication and Innovative Therapeutics (LERMIT) [grant number ANR-10-LABX-33], by the JPIAMR transnational project DesInMBL [grant number ANR-14-JAMR-0002] and by the Région Ile-de-France (DIM Malinf).

dynamics simulation in explicit solvent. The details of all submissions are presented in Table 2.

The "best RMSD" submission without constraints (entry 3) and the "best score" submission with constraints (entry 5) provided similar performances (RMSD of 4.36 Å and 4.34 Å, respectively).

The post-processing of docking poses with molecular dynamics (MD) simulations (entry 7)

clearly improved the quality of the prediction (RMSD of 2.96 Å). Thus, it seems that such MD simulations are beneficial for a better positioning of flexible substituents. The performance of pose prediction is variable within the dataset, depending on the nature of the ligand, as shown in Figure 6. When only half of the ligand structure is correctly positioned

Electronic Supplementary Material

Blinded evaluation of cathepsin S inhibitors from the D3RGC3 dataset using molecular docking and free energy calculations Ludovic Chaput, 1,2, ‡ Edithe Selwa, 1, ‡ Eddy Elisée, 1, ‡ Bogdan I. Iorga 1,* 2 in the main text for details.


Saoussen Oueslati

Pascal Retailleau

Ludovic Marchini

Camille Berthault

Laurent Dortet

Rémy A Bonnin

Bogdan I Iorga

Thierry Naas
email: thierry.naas@aphp.fr

Role of the Arginine 214 in the substrate specificity of OXA-48

Keywords: oxacillinase, carbapenemase, OXA-232, antibiotic resistance, beta-lactamase

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

Introduction

The intensive use of antibiotics to treat infections led to the emergence of multidrug resistant pathogens especially in Gram-negative bacteria (GNB). Among -lactams, carbapenems are considered as last resort antibiotics to treat severe infections caused by GNB [START_REF] Nordmann | Global Spread of Carbapenemase-producing Enterobacteriaceae[END_REF]. The major mechanisms of carbapenem resistance in Enterobacterales are (i) the association of a ß-lactamase that very weakly hydrolyze carbapenems (ESBL or a cephalosporinase) with decreased outer-membrane permeability or (ii) inactivation by specific carbapenem-hydrolyzing -lactamases, carbapenemases [START_REF] Nordmann | Global Spread of Carbapenemase-producing Enterobacteriaceae[END_REF]. Genes encoding for these enzymes are mostly harbored by plasmids explaining their rapid spread. -Lactamases are classified into 4 groups (Amber's classes A to D) based upon sequence homology [START_REF] Nordmann | Global Spread of Carbapenemase-producing Enterobacteriaceae[END_REF][START_REF] Bush | The ABCD's of β-lactamase nomenclature[END_REF]. The clinically-relevant carbapenemases in Enterobacterales belong to class A (KPC-type), B (NDM, VIM, IMP) and D (OXA-48-like). OXA-48, initially identified in Turkey (3), has since rapidly spread in the Mediterranean area, middle East, Europe, India and is now turning into a major global threat [START_REF] Poirel | OXA-48-like carbapenemases: the phantom menace[END_REF]. OXA-48 hydrolyzes penicillins including temocillin, narrowspectrum cephalosporins, and also carbapenems at low rate, but spares expanded-spectrum cephalosporins (ESC) e.g. ceftazidime and cefepime [START_REF] Docquier | Crystal structure of the OXA-48 beta-lactamase reveals mechanistic diversity among class D carbapenemases[END_REF]. Along with the spread of OXA-48, several variants have been reported that differ from OXA-48 by amino-acid (AA) substitutions or deletions (http://bldb.eu/BLDB.php?class=D#OXA) mostly located in the 5-6 loop [START_REF] Naas | Beta-lactamase database (BLDB)structure and function[END_REF]. OXA-181 that differs from OXA-48 by 4 amino-acid substitutions (T103A, N110D, E169Q and S171A) is the second most prevalent OXA-48 variant. OXA-232 that differs from OXA-181 by an additional substitution R214S in the 5-6 loop, is particularly interesting, as its carbapenem-hydrolyzing activity was significantly impaired, as compared to OXA-181 or OXA-48 [START_REF] Oueslati | Heterogeneous hydrolytic features for OXA-48-like β-lactamases[END_REF]. These results suggested a pivotal role of R214 in the hydrolysis of carbapenems. To further investigate the role of R214 in the hydrolytic profile of OXA-48-like carbapenemases, the X-ray structure of OXA-232, and the steady-state kinetic parameters of in-vitro generated OXA-232 mutants at AA position 214 were determined.

Results

Susceptibility testing. To determine the effect of the substitutions at the AA position 214 on the hydrolysis of imipenem and temocillin, MIC values of E. coli expressing OXA-48 and OXA-232 and their respective point mutant derivatives were determined.

The substitution R214S in OXA-48 led to a phenotype similar to that of OXA-232, e.g. reduced MICs for temocillin and imipenem (Table 1). Conversely, OXA-232-S214R restored MICs for temocillin and imipenem similarly to OXA-48/181. When AA214 is substituted with an uncharged AA, such as G and L, the MIC values were similar to those of the OXA-232. The most interesting results were obtained with substitutions by negatively charged AA at pH 7.0 such as aspartic acid and glutamic acid. Indeed, MIC values for imipenem and temocillin were remarkably affected, similar to E. coli Top10 control isolate.

In a second step, we focused on the analysis of mutants of OXA-232: OXA-232-S214K; OXA-232-S214E and OXA-232-S214G, which are representative for each amino acid group: polar positively charged (K), polar negatively charged (E) and uncharged (G).

MIC values for other -lactams were performed and compared to those of OXA-232, OXA-181 (OXA-232-S214R) and OXA-48 (Table 2). Overall, MIC values for benzylpenicillin and cephalothin were not affected and MIC values of the mutants for the other carbapenems (meropenem and ertapenem) vary in the same way as for imipenem.

Biochemical properties determination.

To further characterize the impact of the nature of the residue at position 214 on the hydrolytic profile, steady-state kinetic parameters of the 3 OXA-232 mutants (OXA-232-S214G; OXA-232-S214E and OXA-232-S214K) were determined and compared to those of OXA-48, OXA-181 and OXA-232 [START_REF] Docquier | Crystal structure of the OXA-48 beta-lactamase reveals mechanistic diversity among class D carbapenemases[END_REF][START_REF] Oueslati | Heterogeneous hydrolytic features for OXA-48-like β-lactamases[END_REF]. Overall, the kinetic studies revealed three patterns (Table 3). OXA-232-S214G exhibited hydrolytic activity towards all tested -lactams similarly to those of OXA-232. The mutant OXA-232-S214K possessed a higher catalytic efficiency (k cat / K m = 83 mM -1 .s -1 ) for imipenem of 4fold as compared to OXA-232 (k cat / K m = 23 mM -1 .s -1 ), due to a weak increase of the turnover of the enzyme (k cat ). The most interesting result was observed with the OXA-232-S214E mutant. Indeed, the substitution by a negatively charged amino acid led to a drastic increase of the K m (>2000 µM) of at least 200-folds as compared to OXA-232 (K m = 9 µM), thus decreasing the catalytic efficiency. Moreover, this mutant totally lost its hydrolytic properties for temocillin. All these results were concordant with the observed MIC values. Taken together, our results confirm that the hydrolysis of imipenem depends on the nature of the residue at position 214. It appears in light of these results that positively charged amino-acids favor hydrolysis of imipenem.

Crystal Structure of OXA-232. The crystals of OXA-232 were obtained at basic pH in 30% w/v PEG3000 and in the presence of 0.2 M Li 2 SO 4 (Figure 2 A). Glycerol was added to the mother liquor to cryoprotect the crystals before flash-annealing for data collection.

They belonged to hexagonal space group, P 6 2 , with 43.8 % solvent as calculated from the Matthews coefficient of 2.19 Å 3 Da -1 [START_REF] Matthews | Solvent content of protein crystals[END_REF][START_REF] Docquier | Crystal structure of the OXA-48 beta-lactamase reveals mechanistic diversity among class D carbapenemases[END_REF]. At that time, the best available OXA-48 model in the PDB was 3HBR (1.9Å) (5) obtained from crystals grown at lower pH (7.5 vs 8.5), in PEG 4000 then dipped into cryoprotectant ethyleneglycol, and which belonged to the monoclinic space group P 2 1 with 48.3% solvent (Matthews coefficient, 2.38 Å 3 Da -1 ). After deletion of 22 N-ter amino acids, one half of R214S-3HBR coordinates were used as quasi-homologous model for phasing the hexagonal crystal to 2.2 Å resolution. The refined electron density appeared neat throughout the entire backbone and most of side chains of the biological dimer per asymmetric unit of the crystal. Alternate positions for Q41 and Q53 were observed as well for D230, whose one position permitted a salt bridge with R107 at the dimer interface. The average Real-space correlation factor calculated by Sfcheck software (10) is 0.932. Compared to 3HBR, up to three additional C-ter amino acids could be displayed in chain B of OXA-232 (only one in chain A) but not the 9-His-tag chain, disordered in the surrounding solvent.

Additionally, some minor discontinuities observed in some parts of 3HBR main chains could be fixed in 5HFO (N50D/K51D, T99C, S150A h-bonding D148A). The rmsd values with the two 3HBR dimers range between 0.3452 and 0.4197 Å (Superpose, [START_REF] Krissinel | Secondary-structure matching (SSM), a new tool for fast protein structure alignment in three dimensions[END_REF] 196-199; β5: 204-212; β6: 219-227; β7: 232-240). 98.2% of all residues were inside the favored regions of the Ramachandran plot, and 1.8% in the allowed regions. The quaternary structure of OXA-232 is dimeric as observed for OXA-48 around a non-crystallographic twofold axis, in compliance with our data of size-exclusion chromatography. It was already shown that less monomer surface (200Å 2 ) was buried in that dimer formation compared to OXA-10 and it was built upon an intermolecular β-sheet involving β4 from each subunit linked by reciprocal H-bonds between the A199 NH and the A199 carbonyl O of one chain with the corresponding one related by the ncs two-fold axis. Several other H-bonds and salt bridges also participated in the dimeric interface stabilization, and among these interactions, an anionic binding site lying on the ncs two-fold axis is tweezed by two facing R206 from ncs related β5 sheets. This contrasts with the cation-binding site described at the interface of OXA-10, involving other residues replaced in OXA-48. In the original 3HBR structure, a water molecule was placed but on the faith of spherical electron density shape, low B-factor and recurrent presence of halogens in other class D beta-lactamase structures (but not upon anomalous signal), it was substituted by a residual chloride coming from previous purification steps. The same misinterpretation was also done in in the structure 4S2K [START_REF] King | Molecular Mechanism of Avibactam-Mediated β-Lactamase Inhibition[END_REF], where chlorine was also reported in the crystallizing medium composition. If sulfate ions have also been seen making that salt bridge in class D structures grown in the presence of large amount of sulfate salts (e.g. structures 6NLW and 5FDH), here six such anions have been spotted essentially at the molecule surface as well as five opportune glycerol molecules. Carbamylated lysines (KCX 73) can be described in a similar environment in both structures. Regarding the mutation of R214S (the R214 residue was shown in water-mediated interaction with the avibactam drug in the active site in the structure 4S2K [START_REF] King | Molecular Mechanism of Avibactam-Mediated β-Lactamase Inhibition[END_REF], this has not structurally disturbed the β5-β6 hairpin and its direct environment, just adding few solvent molecules in the created void (2-3 waters, and one glycerol molecule in chain A) and modifying the L158B and I215B side chain conformers. Deeper into both active sites, respective carbamylated lysine (KCX 73) can be described in a similar environment in both structures OXA-232 and OXA-48. We note however the presence of a sulfate ion in interaction with the same residues that recognize the sulfate moiety of the avibactam (Figure 2B). More recently the structure of another variant, OXA-181; which differs from OXA-232 by five substitutions T103A, N110D, E169Q, S171A and R214S, has been released (PDB code 5OE0 (13)) at a similar resolution (2.05 Å). The rmsd with OXA-232 is even lower than the former overlay (0.224 and 0.252Å) in relation to the fact that structure crystallized in the same unit cell dimensions and space group P6 2 than OXA-232. However, five residue mutations are on the protein surface that could potentially modulate the contacts leading to the crystal growth and crystallization conditions were slightly different with a neutral pH and PEG-mme 5000 instead of PEG 3000 or PEG 4000 for 3HBR. Like OXA-232, sulfate has been used in low quantity in the crystallization drop and one anion was reported in the final 5HFO model bound to the amide N of R186A as for OXA-232. Variability in β-lactamase crystal forms remains tenuous to explain or even predict, since crystals of OXA-48 in complex with inhibitors grew in the same unit cell dimensions (the choice of P32 space group instead of P62 could be interpreted by the author's intention to refine two more monomers and active sites independently) in different crystallizing solution conditions (MPD was the major component) and in absence of sulfate.

Molecular modelling.

An in-silico study was performed to identify the structural determinants that could explain the experimentally determined differences between the hydrolytic profiles of OXA-232 and its variants in comparison with OXA-48. The OXA-232 structure was used as a starting point and the mutations S214G, S214K and S214E were modeled based on predicted low energy conformations [START_REF] Dunbrack | Rotamer libraries in the 21st century[END_REF]. The resulting models showed that lateral chains of serine and lysine mutants are positioned in the same axis and direction as the lateral chain of arginine, without any clashes. The glutamate mutant would be positioned differently, oriented towards the active site, establishing hydrogen bonds with the backbone nitrogen and with the sidechain hydroxyl of T213. This is the only conformation that avoids clashes of this glutamate mutant with the neighboring residues I215 and D159, the latter being involved in OXA-48 in a salt bridge with R214 that maintains a closed conformation of the active site, which would stabilize the substrates (including carbapenems) in a conformation compatible with an efficient hydrolysis, thus influencing the carbapenem turnover [START_REF] Docquier | Crystal structure of the OXA-48 beta-lactamase reveals mechanistic diversity among class D carbapenemases[END_REF]. For OXA-232 and its variants the salt bridge R214-D159 cannot be formed, and this may explain the decrease of k cat values compared to OXA-48.

Molecular docking calculations of imipenem and temocillin on the R214E mutant of OXA-48 provided explanations for the different hydrolytic parameters that were observed experimentally. Imipenem showed two alternative binding modes, one that was known, with a single ionic interaction between the imipenem carboxylate and the sidechain of R250, compatible with a nucleophilic attack by S70 (Figure 2C), and a second one, presumably more stable, with two ionic interactions, one between the imipenem carboxylate and the sidechain of R250 as in the previous binding mode, and a second between the positively charged R2 substituent of imipenem and the mutated residue E214 (Figure 2C). This latter binding mode is not productive, as it is not compatible with a nucleophilic attack by S70, thus increasing the apparent K m value. To confirm this hypothesis, we have determined the IC 50 for imipenem with OXA-232 and with the mutant OXA-232 S214E, the only variant with an increased K m value. Interestingly, the values obtained (2 µM and 0.059 µM, respectively)

showed that the mutant OXA-232 S214E has an IC 50 ~34-fold lower as compared to that of OXA-232. On the other hand, temocillin is the only β-lactam tested possessing a negative charge on the R1 substituent. Our docking calculations showed that this negative charge would establish strong unfavorable ionic interactions with the mutated residue E214, thus precluding the binding in a conformation compatible with a nucleophilic attack by S70 (Figure 2D).

Discussion

OXA-48-producing Enterobacterales are now endemic in many countries, such as Turkey, Middle East, North Africa, India and have widely spread across Europe. Since the first description of OXA-48, several variants have been described [START_REF] Nordmann | Global Spread of Carbapenemase-producing Enterobacteriaceae[END_REF][START_REF] Naas | Beta-lactamase database (BLDB)structure and function[END_REF]. These variants can be classified into 3 groups according to their hydrolysis profile. Most of them, including OXA-181 or OXA-162 [START_REF] Oueslati | Heterogeneous hydrolytic features for OXA-48-like β-lactamases[END_REF], have an enzymatic activity similar to OXA-48 [START_REF] Docquier | Crystal structure of the OXA-48 beta-lactamase reveals mechanistic diversity among class D carbapenemases[END_REF]. The second group, represented by OXA-163 [START_REF] Oueslati | Heterogeneous hydrolytic features for OXA-48-like β-lactamases[END_REF], OXA-247 [START_REF] Gomez | Intrapatient emergence of OXA-247: a novel carbapenemase found in a patient previously infected with OXA-163-producing Klebsiella pneumoniae[END_REF] and OXA-405 [START_REF] Dortet | Genetic and Biochemical Characterization of OXA-405, an OXA-48-Type Extended-Spectrum β-Lactamase without Significant Carbapenemase Activity[END_REF], have no carbapenemase activity, but instead a marked hydrolytic activity against ESC, similarly to OXA-ESBLs [START_REF] Poirel | Diversity, epidemiology, and genetics of class D beta-lactamases[END_REF].

Finally, the third group represented by OXA-244 [START_REF] Potron | Characterisation of OXA-244, a chromosomally-encoded OXA-48-like β-lactamase from Escherichia coli[END_REF] and OXA-232 [START_REF] Oueslati | Heterogeneous hydrolytic features for OXA-48-like β-lactamases[END_REF] exhibit an overall reduced activity towards all ß-lactams including carbapenems as compared to OXA-48. The amino acid sequence comparison of OXA-48-variants suggested a link between the primary structure and the function of these enzymes. Indeed, all OXA-48-variants with an OXA-ESBLs phenotype (loss of carbapenem hydrolysis and the gain of activity towards ESC) have amino-acid deletions in the 5-6 loop (Figure 1). This observation suggests that this loop plays a role in the substrate specificity. The phenotypic study of OXA-244 (which differs from OXA-48 by only one substitution R214G) [START_REF] Potron | Characterisation of OXA-244, a chromosomally-encoded OXA-48-like β-lactamase from Escherichia coli[END_REF] and the enzymatic study of OXA-232 (which differs from OXA-181 by a single substitution R214S) underline that the residue 214 is crucial for the carbapenem and temocillin hydrolysis. To confirm this hypothesis, we generated mutants of OXA-48 and OXA-232 at the position 214, and analyzed their hydrolytic profiles including OXA-232-S214R and OXA-48-R214S. Thus, the substitution of S214R for OXA-232 re-stored hydrolytic properties toward temocillin and imipenem.

Similarly, the substitution of R214S for OXA-48, led to a drastic decrease of imipenem and temocillin hydrolysis. In order to better understand this phenomenon, 5 mutants at position 214 of OXA-48 and OXA-232 were generated. These substitutions corresponded to amino acid that were: polar positively charged (214K), polar negatively charged (214E and 214D), non-polar (214L) and glycine (214G). Overall, the substitutions with D, E, L and G led to a decreased MICs of imipenem and temocillin. The substitution S214K in OXA-232 led to increased MIC values, but still lower than those of OXA-48. The biochemical study of OXA-232-S214G and OXA-232-S214K showed that these 2 substitutions did not affect the apparent affinity toward -lactams (K m values were similar to those of OXA-232 and OXA-48) but have an effect on the acylation or deacylation steps of the substrate catalysis (k cat values were similar to those of OXA-232 but smaller than those of OXA-48). The catalytic efficiencies (k cat /K m ) of these two mutants were in agreement with the MIC values. The most significant differences were observed with OXA-232-S214E, with a drastic decrease in the apparent affinity for imipenem (K m was at least 200 fold higher as compared to OXA-232 and OXA-48) and a loss of temocillin hydrolysis. The glutamate substitution seems to have a direct effect on the apparent affinity of imipenem. Analysis of the 3D structure of OXA-48 showed that R214 interacts with D159 via a salt bridge [START_REF] Docquier | Crystal structure of the OXA-48 beta-lactamase reveals mechanistic diversity among class D carbapenemases[END_REF], which maintains the shape and the network of water molecules within the binding site. Our molecular modeling study revealed that with the G214 or S214 mutations, this interaction with D159 is lacking, which presumably increases the flexibility of this part of the binding site. In contrast, K is a positively charged polar amino-acid, which appears to maintain the interaction with D159, even with a lateral chain shorter than R. In the mutant E214, we evidenced an alternative binding conformation of imipenem that is not compatible with a nucleophilic attack by S70.

In this case, imipenem acts as an inhibitor with an IC50 value of 0.059µM, which explains the significantly higher K m values determined experimentally for this substrate. We also showed that, in the same E214 mutant, the unfavorable interaction between the negatively charged R1 substituent of temocillin with the E214 residue precludes the binding of this antibiotic, thus explaining the structural basis responsible for the lack of temocillin hydrolysis observed experimentally.

Overall, we demonstrated that the AA position 214 in OXA-48-like β-lactamases is critical for the carbapenemase activity but also for the temocillin hydrolysis. This point is of outmost clinical importance and explains why detection of OXA-244 or OXA-232-producers remains a challenge for clinical microbiology laboratories [START_REF] Potron | Characterisation of OXA-244, a chromosomally-encoded OXA-48-like β-lactamase from Escherichia coli[END_REF][START_REF] Hoyos-Mallecot | OXA-244-Producing Escherichia coli Isolates, a Challenge for Clinical Microbiology Laboratories[END_REF][START_REF] Potron | Genetic and biochemical characterisation of OXA-232, a carbapenemhydrolysing class D β-lactamase from Enterobacteriaceae[END_REF]. Indeed, these isolates do not grow on ChromID Carba Smart (bioMérieux, Marcy L'Etoile, France), one of the most used type of medium for the screening of CPEs [START_REF] Hoyos-Mallecot | OXA-244-Producing Escherichia coli Isolates, a Challenge for Clinical Microbiology Laboratories[END_REF]. The ChromID Carba Smart is a biplate containing on one side a carbapenem and on the other temocillin, two substrates that are only weakly hydrolyzed by OXA-244 and OXA-232 [START_REF] Hoyos-Mallecot | OXA-244-Producing Escherichia coli Isolates, a Challenge for Clinical Microbiology Laboratories[END_REF].

Materials and methods

Bacterial strains. The clinical strain Escherichia coli LIEU (20) expressing the OXA-232 lactamase was used to clone the bla OXA-232 gene. E. coli TOP10 (Invitrogen, Saint-Aubin, France) was used for cloning and mutagenesis experiments and E. coli BL21 DE3 was used for overexpression experiments (Novagen, Fontenay-sous-Bois, France).

Antimicrobial agents, susceptibility testing and microbiological techniques.

Antimicrobial susceptibilities were determined by the disc diffusion technique on Mueller-Hinton agar (Bio-Rad, Marnes-La-Coquette, France) and interpreted according to the EUCAST breakpoints, updated in May 2018 (http://www.eucast.org). MICs were determined using the Etest technique (bioMérieux). Antibiotics were purchased from Sigma (Saint-Quentin-Fallavier, France) except temocillin (Eumedica, Brussels, Belgium) PCR, cloning, site-directed mutagenesis and DNA sequencing. The recombinant plasmids pTOPO-bla OXA-232 and pTOPO-bla OXA-48 , obtained from a previous study [START_REF] Oueslati | Heterogeneous hydrolytic features for OXA-48-like β-lactamases[END_REF], were used as a template for the site-directed mutagenesis assays and specific primers were designed for the different mutations, using the program QuickChange Primer Design (Agilent Technologies).

QuikChange II Site-Directed Mutagenesis Kit (Agilent Technologies) was used, following the manufacturer's recommendations, in order to substitute the AA214 into a Glycine (G), a Lysine (K), a Leucine (L), an Aspartic acid (D), a Glutamic acid (E), a Serine (S) and an Arginine (R). Mutagenesis reaction products were transformed in E. coli TOP10 (Invitrogen, Saint-Aubin, France) and selection was performed on TSA plate containing Kanamycin (50g/ml). For the production of OXA-232, genes were amplified by PCR using the forward primer OXA 23-256NdeI (5'-aaaaaCATATGaaggaatggcaagaaaacaaa-3') and the reverse primer OXA XhoI-stop (5'-aaaaaCTCGAGggggaataatttttcctgtttgag-'3) to increase the purification yield.

For the mutants OXA-232 S214E, OXA-232 S214G and OXA-232 S214K the forward primers OXA NdeI (5′-aaaaaCATATGTTGGTGGCATCGATTATCGG-3′) and reverse primer OXA XhoI-stop

(5-aaaaaCTCGAGGAGCACTTCTTTTGTGATGGC-3') were used.

Recombinant plasmids were cloned into pET41b (for OXA-232 wt) allowing the expression of the enzyme with an His-tag and pET9a (for OXA-232 mutants) vector (Invitrogen, Life

Technologies, Cergy-Pontoise, France), then transformed into E. coli BL21 DE3 (Novagen).

All the recombinant plasmids were extracted using the Qiagen miniprep Kit and sequenced using a T7 promoter and M13 reverse primers or T7 terminator (depending on the plasmid) with an automated sequencer (ABI Prism3100; Applied Biosystems). The nucleotide sequences were analyzed using software available at the National Center for Biotechnology Information website (http://www.ncbi.nlm.nih.gov).

Protein purification. Overnight culture of E. coli BL21 DE3 harboring recombinant plasmids pET41b-OXA-232 or pET9a-OXA-232-mut were used to inoculate 2 L of LB medium broth containing 50 mg/L kanamycin. Bacteria were cultured at 37°C until an OD of 0.6 at 600 nm was reached. The -lactamase was induced overnight with 0.2 mM IPTG as inducer and cultures were centrifuged at 6000 g for 15min. The pellets were resuspended with the binding buffer 25mM phosphate sodium pH 7.4, 300 mM K 2 SO 4 , 10 mM imidazole for OXA-232 wt and with the buffer 20 mM Bis-Tris H 2 SO 4 (pH 7.2) for mutants of OXA-232.

Bacterial cells were disrupted by sonication and the bacterial pellet was removed by two consecutive centrifugation steps at 10 000 g for 1h at 4°C; the supernatant was then centrifuged at 96 000 g for 1h at 4°C. OXA-232 wt was purified with a NTA-Nickel column (GE Healthcare, Freiburg, Germany) by using the elution buffer 25 mM phosphate sodium pH 7.4, 300 mM K2SO4, 500 mM imidazole. Mutants of OXA-232 were purified by using 2 anion-exchange chromatography HiTrap TM QHP GE Healthcare (20 mM Bis-Tris H 2 SO 4 pH 7.2, then 20 mM piperazine H 2 SO 4 pH 9.5) [START_REF] Oueslati | Heterogeneous hydrolytic features for OXA-48-like β-lactamases[END_REF]. Finally, a gel filtration step was performed for the purifications of all -lactamases with 100mM sodium phosphate buffer pH 7 and 150 mM NaCl with a Superdex 75 column (GE Healthcare, Freiburg, Germany). The protein purity was estimated by SDS-PAGE and the pooled fractions were dialyzed against 10 mM Tris-HCl pH 7 for the mutants and against 0.1 M Hepes (pH 7.5) for OXA-232 wt and concentrated using vivaspin 20 (10 000 MWCOPES Sartorius) columns. Protein concentration was determined according to the Bradford method using the Bio-Rad protein assay standard II kit (Bio-Rad, Marnes-la-Coquette, France) with bovine serum albumin (BSA) as a standard.

Kinetics assays. Steady-state kinetic parameters were determined using a spectrophotometer ULTROSPEC 2000 (Amersham Pharmacia Biotech) and performed at 30°C in 100 mM phosphate buffer (pH 7) for 10 minutes [START_REF] Oueslati | Heterogeneous hydrolytic features for OXA-48-like β-lactamases[END_REF]. The disappearance of substrate (-lactam antibiotics) was monitored at the specific wavelengths and converted to initial velocities using the specific extinction coefficients. The k cat and K m values were determined by analyzing lactam hydrolysis under initial-rate conditions using the Eadie-Hoffstee linearization of the Michaelis-Menten equation with SFWIFT II software [START_REF] Oueslati | Heterogeneous hydrolytic features for OXA-48-like β-lactamases[END_REF]. Fifty percent inhibitory concentrations (IC 50 ) for imipenem were determined in 100 mM sodium phosphate buffer (pH 7) and 100 µM of benzylpenicillin as a reporter substrate.

Crystallization and data collection. Crystal conditions were screened with 10 mg/mL of concentrated protein using a range of commercially available screens. OXA-232 crystals grew in 0. Structure determination and refinement. X-ray diffraction data sets were processed using XDS [START_REF] Kabsch | XDS[END_REF] and AIMLESS [START_REF] Evans | How good are my data and what is the resolution?[END_REF]; the crystals diffracted to 2.2 Å resolution (Table 4). Initial phases were obtained by molecular replacement with MOLREP (23) using OXA-48 (PDB ID 3HBR) as a search probe. Refinement was performed by successive and alternate rounds of refinement with BUSTER (24) and model improvement using Coot [START_REF] Emsley | Features and development of Coot[END_REF]. The final model was evaluated using MolProbity [START_REF] Williams | MolProbity: More and better reference data for improved all-atom structure validation[END_REF]. Data collection and refinement statistics are provided in Table 4.

Molecular modelling. Molecular modelling was performed to evaluate the effects of the mutations on OXA-232 β-lactamase. Substitution of amino acid was performed by mutating the OXA-232 structure (PDB code 5HFO) in silico using the Dunbrack rotamer library (swapaa command) as a part of the UCSF Chimera software [START_REF] Dunbrack | Rotamer libraries in the 21st century[END_REF][START_REF] Pettersen | UCSF Chimera--a visualization system for exploratory research and analysis[END_REF]. The Dunbrack rotamer library predicts the conformation of the amino acid side-chain based on the global energy minimum of the protein. The identification of interatomic clashes based on VDW radii [START_REF] Tsai | The packing density in proteins: standard radii and volumes[END_REF] was performed with UCSF Chimera software [START_REF] Dunbrack | Rotamer libraries in the 21st century[END_REF][START_REF] Pettersen | UCSF Chimera--a visualization system for exploratory research and analysis[END_REF]. Three-dimensional structures of the βlactam ligands were generated using Corina 3.60 (Molecular Networks GmbH, Erlangen, Germany). Molecular docking calculations were performed using Gold (Cambridge Crystallographic Data Centre, Cambridge, UK) (29) and the GoldScore scoring function. The binding site, defined as a 20 Å radius sphere, was centered on the OG oxygen atom of Ser70.

All other parameters had default values. The receptor-ligand complex images were produced using UCSF Chimera [START_REF] Dunbrack | Rotamer libraries in the 21st century[END_REF][START_REF] Pettersen | UCSF Chimera--a visualization system for exploratory research and analysis[END_REF]. Alternative docking conformations of imipenem (colored in magenta (A) and yellow (B), closed form of the β-lactam ring) as non-covalent (Michaelis) complexes with the OXA-48 R214E mutant (colored in cyan, generated from structure 6P97). Hydrogen bonds and favorable ionic interactions are shown as blue springs. Protein hydrogens are hidden for clarity. D Docking conformation of temocillin (colored in green, closed form of the β-lactam ring) as non-covalent (Michaelis) complex with OXA-48 (colored in pink, PDB 6P97), superposed with the OXA-48 R214E mutant (colored in purple, generated from structure 6P97, only Gln214 residue shown). Hydrogen bonds and favorable ionic interactions are represented as blue springs, and unfavorable ionic interactions are represented as yellow thick lines. Protein hydrogens are hidden for clarity.

2 M

 2 Lithium sulfate, 0.1M Tris pH 8.5 and 30% w/v PEG 3000. Crystals were dipped briefly in cryoprotectant (20% v/v glycerol in the reservoir solution) and flash-frozen in liquid nitrogen prior to be transferred in a stream of nitrogen at 100K (delivered by a Rigaku Xtream cryosystem). A 2.21-Å resolution data set was collected on a Rigaku MicroMax-007-HF rotating-anode generator with Cu Kα radiation, Varimax HF mirror focusing optics (Rigaku) and MAR345 image-plate detector (MAR Research).

Figure 1 .

 1 Figure 1. Amino acid sequence alignment of OXA-48 variants. Asterisks indicate identical residues in all the three sequences, colon indicate a substitution with another amino acid but with the same proprieties. Amino acid motif that are well conserved among class D lactamases are indicated by grey boxes, and the black-outlined box corresponds to the 5-6 loop. Positions D159 and R/S214 are indicated by a thin and large arrow, respectively. CLUSTAL O (1.2.4) multiple sequence alignment.Figure 2. A Superposition of crystal structures and of OXA-48 (yellow, PDB: 3HBR) and OXA-232 (grey; PDB: 5HFO). The 5-6 loop is delimited by the circle. B Partial active-site closeup of sulfate 302A-bound OXA-232. Atoms are colored by type (C, white;N, blue; O, red; S, yellow). Hydrogen bonding and electrostatic interactions are depicted as blue dashes. CAlternative docking conformations of imipenem (colored in magenta (A) and yellow (B), closed form of the β-lactam ring) as non-covalent (Michaelis) complexes with the OXA-48 R214E mutant (colored in cyan, generated from structure 6P97). Hydrogen bonds and favorable ionic interactions are shown as blue springs. Protein hydrogens are hidden for clarity. D Docking conformation of temocillin (colored in green, closed form of the β-lactam ring) as non-covalent (Michaelis) complex with OXA-48 (colored in pink, PDB 6P97), superposed with the OXA-48 R214E mutant (colored in purple, generated from structure 6P97, only Gln214 residue shown). Hydrogen bonds and favorable ionic interactions are represented as blue springs, and unfavorable ionic interactions are represented as yellow thick lines. Protein hydrogens are hidden for clarity.

Figure 2 .

 2 Figure 1. Amino acid sequence alignment of OXA-48 variants. Asterisks indicate identical residues in all the three sequences, colon indicate a substitution with another amino acid but with the same proprieties. Amino acid motif that are well conserved among class D lactamases are indicated by grey boxes, and the black-outlined box corresponds to the 5-6 loop. Positions D159 and R/S214 are indicated by a thin and large arrow, respectively. CLUSTAL O (1.2.4) multiple sequence alignment.Figure 2. A Superposition of crystal structures and of OXA-48 (yellow, PDB: 3HBR) and OXA-232 (grey; PDB: 5HFO). The 5-6 loop is delimited by the circle. B Partial active-site closeup of sulfate 302A-bound OXA-232. Atoms are colored by type (C, white;N, blue; O, red; S, yellow). Hydrogen bonding and electrostatic interactions are depicted as blue dashes. CAlternative docking conformations of imipenem (colored in magenta (A) and yellow (B), closed form of the β-lactam ring) as non-covalent (Michaelis) complexes with the OXA-48 R214E mutant (colored in cyan, generated from structure 6P97). Hydrogen bonds and favorable ionic interactions are shown as blue springs. Protein hydrogens are hidden for clarity. D Docking conformation of temocillin (colored in green, closed form of the β-lactam ring) as non-covalent (Michaelis) complex with OXA-48 (colored in pink, PDB 6P97), superposed with the OXA-48 R214E mutant (colored in purple, generated from structure 6P97, only Gln214 residue shown). Hydrogen bonds and favorable ionic interactions are represented as blue springs, and unfavorable ionic interactions are represented as yellow thick lines. Protein hydrogens are hidden for clarity.

Table 1 :

 1 Minimum inhibitory concentrations (MICs) of -lactams for E. coli Top 10 (pTOPO-OXA-48 and variants), E. coli Top 10 (pTOPO-

	OXA-232 and variants) and E. coli Top 10.
							MIC (mg/L)
	β-lactams	E. coli Top 10 pTOPO-OXA-48 wt R214S R214G R214L R214K R214D R214E	E. coli Top 10 pTOPO-OXA-232 wt S214R S214G S214L S214K S214D S214E	E. coli Top 10
	Temocillin >256 32	32	64	128	8	8	32	256	32	64	128	8	8	8
	Imipenem 0.75 0.25	0.25	0.25	0.38	0.25	0.25	0.38 0.75	0.25	0.25	0.5	0.25	0.25	0.25

Table 2 :

 2 Minimum inhibitory concentrations (MICs) of -lactams for E. coli Top 10 (pTOPO-OXA-48), E. coli Top 10 (pTOPO-OXA-

	181), E. coli Top 10 (pTOPO-OXA-232 and its variants), and E. coli Top10.
					MIC (mg/L)
		E. coli Top 10 (pTOPO-)					E. coli Top 10
	β-lactams	OXA-232 S214G	OXA-232 S214E	OXA-232 S214K	OXA-232 OXA-181	OXA-48
	Benzypenicillin	>256	>256	>256	>256	>256	>256	64
	Temocillin	32	8	128	32	>256	>256	8
	Cephalothin	16	16	16	16	16	16	4
	Imipenem	0.25	0.25	0.5	0.38	0.75	0.75	0.25
	Meropenem	0.047	0.032	0.094	0.047	0.125	0.094	0.016
	Ertapenem	0.012	0.012	0.094	0.012	0.19	0.094	0.004

Table 3 :

 3 Kinetic parameters of OXA-232 mutants NH, no detectable hydrolysis was observed with 1 μM of purified enzyme and 500 μM of substrate; ND, not determined. a Standard deviation values are indicated after the kinetic parameter values.

			Km (M) a			k cat (s -1 ) a			k cat /K m (mM -1 .s -1 ) a
	Substrate	OXA-232	OXA-232	OXA-232	OXA-232	OXA-232	OXA-232	OXA-232	OXA-232	OXA-232
		S214E	S214G	S214K	S214E	S214G	S214K	S214E	S214G	S214K
	Benzypenicillin	140 ± 5	95± 6.5	95 ± 5	235 ± 20.6	130 ± 11.6	115 ± 9	1678 ± 91.6	1369 ± 166.5	1189 ± 72.2
	Temocillin	NH	63 ± 6.7	28 ± 3.6	NH	0.01 ± 0.007 0.03 ± 0.007	NH	0.15 ± 0.07	1 ± 0.29
	Cephalothin	100±7	133 ± 12.1 160 ± 13.8	8 ± 0.6	4 ± 0.6	5 ± 0.6	81 ± 10.6	30 ± 3.8	31 ±1.7
	Imipenem	>2000	5 ± 2,6	6 ± 0.7	ND	0.1 ± 0.07	0.5 ± 0.014	ND	20 ± 6.9	83 ± 8.9

Acknowledgments

All the authors are members of the Laboratory of Excellence in Research on Medication and Innovative Therapeutics (LERMIT). This work was partially funded by the University Paris-Saclay, and by grants from the French National Research Agency (ANR-10-LABX-33, and ANR-17-ASTR-0018).

Transparency declarations

We have no competing interests to declare.

PDB accession codes.

The OXA-232 structure atomic coordinates were deposited in the Protein Data Bank with the accession code 5HFO.