Merveille Eguida

Didier Rognan
email: rognan@unistra.fr

A computer vision approach to align and compare protein cavities: Application to fragment-based drug design

Identifying local similarities in ligand binding sites from distant proteins is still a major hurdle to rational drug design approaches. We herewith present a novel method, borrowed from computer vision, particularly adapted to mine fragment subpockets and compare them to entire ligand-binding sites.

Pockets are represented by pharmacophore-annotated point clouds mimicking ideal ligands or fragments. Point cloud registration is used to find the transformation enabling an optimal overlap of points sharing similar topological neighborhoods and pharmacophoric features. Importantly, the local environment of randomly sampled cavity points is used to generate a preliminary alignment that is next refined by an iterative closest point algorithm. A scoring function has been tuned to quantify the degree of shape and pharmacophoric overlap, and assess pocket pairwise similarity. The method (ProCare) was first calibrated on a large set of known similar and dissimilar druggable cavities, shown to be insensitive to moderate variations in atomic coordinates, and then applied to the specific problem of comparing fragment subpockets with entire cavities. Starting from a unique set of cavity points, point cloud registration outperformed a state-of-the-art computational method (shape-based similarity search) in detecting local similarities between fragment subsites and entire cavities. A collection of 33,953 subpockets annotated with their bound fragments was screened for local similarity to cavities from three recently described protein X-ray structures. ProCare was able to detect local similarities between remote pockets and transfer the corresponding fragments to the query cavity space. Interestingly, fragments selected from totally unrelated proteins, nicely overlap substructures of the masked original ligand co-crystallized with the target query, thereby proposing an automated first step to an automated fragment-based design approach targeting ligand-orphan cavities.

INTRODUCTION

Three-dimensional (3D) structures of protein-ligand complexes are the corner stones of structurebased rational approaches to ligand design. [START_REF] Goodsell | RCSB Protein Data Bank: Enabling biomedical research and drug discovery[END_REF] Among the many computational methods [START_REF] Rognan | Chemogenomic approaches to rational drug design[END_REF] to infer putative relationships between ligand and target spaces, detection and pairwise comparison of protein-ligand binding sites have gained considerable popularity in the last decade. [START_REF] Perot | Druggable pockets and binding site centric chemical space: a paradigm shift in drug discovery[END_REF][START_REF] Volkamer | Analyzing the topology of active sites: on the prediction of pockets and subpockets[END_REF][START_REF] Ehrt | Impact of Binding Site Comparisons on Medicinal Chemistry and Rational Molecular Design[END_REF] Potential cavities can be first detected at the surface of macromolecules using a myriad of computational tools, [START_REF] Ehrt | Impact of Binding Site Comparisons on Medicinal Chemistry and Rational Molecular Design[END_REF] classically grouped in three categories: geometry-based (e.g. CavBase, [START_REF] Schmitt | A new method to detect related function among proteins independent of sequence and fold homology[END_REF] VolSite, [START_REF] Desaphy | Comparison and druggability prediction of protein-ligand binding sites from pharmacophore-annotated cavity shapes[END_REF] Fpocket 8 ), energybased (e.g. GRID, [START_REF] Goodford | A computational procedure for determining energetically favorable binding sites on biologically important macromolecules[END_REF] Q-SiteFinder [START_REF] Laurie | Q-SiteFinder: an energy-based method for the prediction of proteinligand binding sites[END_REF] ) and evolutionary-based (e.g.SURFNET-ConSurf [START_REF] Glaser | A method for localizing ligand binding pockets in protein structures[END_REF] ), although some methods may combine different approaches (e.g. Ligsitecsc [START_REF] Huang | LIGSITEcsc: predicting ligand binding sites using the Connolly surface and degree of conservation[END_REF] , SiteMap [START_REF] Halgren | Identifying and characterizing binding sites and assessing druggability[END_REF] ). Whereas geometry-based approaches rely on the prior calculation of the target's molecular surface to identify accessible pockets, energy-based methods compute interaction energies on a 3D lattice between the target protein and several probe atoms. Last, evolutionary-based tools require a multiple sequence or structural alignment of targets from the same family to pinpoint evolutionary conserved motifs that can be linked to the recognition of peculiar ligand structures. Interestingly, structural druggability or ligandability, [START_REF] Vukovic | Quantitative metrics for drug-target ligandability[END_REF] the propensity to accommodate high-affinity drug-like ligands, can be computed on the fly using machine-learning models [START_REF] Guilloux | Fpocket: an open source platform for ligand pocket detection[END_REF][START_REF] Desaphy | Comparison and druggability prediction of protein-ligand binding sites from pharmacophore-annotated cavity shapes[END_REF] trained on sets of known druggable and undruggable sites. Once pockets have been detected, they can be systematically compared at a high-throughput to detect global similarities even in absence of fold conservation. [START_REF] Ehrt | Impact of Binding Site Comparisons on Medicinal Chemistry and Rational Molecular Design[END_REF] Many descriptors (fingerprints, distance counts, pharmacophoric triplets, grid points, point clouds, graphs, and shapes) of protein-ligand binding pockets can be used by geometric hashing [START_REF] Shulman-Peleg | Recognition of functional sites in protein structures[END_REF] or clique detection [START_REF] Schmitt | A new method to detect related function among proteins independent of sequence and fold homology[END_REF] algorithms to find the most prominent shared features guiding the structural alignment of protein cavities.

Following the basic principle that similar cavities recognize similar ligands, protein-ligand binding site comparison methods have been successfully used in many drug discovery scenarios: (i) assigning a function from a target's 3D structure, [START_REF] Kinoshita | Identification of protein biochemical functions by similarity search using the molecular surface database eF-site[END_REF][START_REF] Tseng | Classification of protein functional surfaces using structural characteristics[END_REF][START_REF] Konc | Structure-based function prediction of uncharacterized protein using binding sites comparison[END_REF] (ii) finding hits for a novel target, [START_REF] Willmann | Impairment of prostate cancer cell growth by a selective and reversible lysine-specific demethylase 1 inhibitor[END_REF] (iii) prioritizing compound library design, [START_REF] Al-Gharabli | An efficient method for the synthesis of peptide aldehyde libraries employed in the discovery of reversible SARS coronavirus main protease (SARS-CoV Mpro) inhibitors[END_REF] (iv) repurposing ancient drugs for new targets, 21-23 (v) explaining the poypharmacological profile of known drugs, [START_REF] Xie | Drug Discovery Using Chemical Systems Biology: Identification of the Protein-Ligand Binding Network To Explain the Side Effects of CETP Inhibitors[END_REF] (vi) predicting unexpected off-targets [START_REF] Defranchi | Binding of protein kinase inhibitors to synapsin I inferred from pair-wise binding site similarity measurements[END_REF][START_REF] Cleves | Chemical and protein structural basis for biological crosstalk between PPARalpha and COX enzymes[END_REF][START_REF] Kakisaka | A Novel Antiviral Target Structure Involved in the RNA Binding, Dimerization, and Nuclear Export Functions of the Influenza A Virus Nucleoprotein[END_REF][START_REF] Schirris | Mitochondrial ADP/ATP exchange inhibition: a novel off-target mechanism underlying ibipinabant-induced myotoxicity[END_REF] and extending potential binding sites to new areas of target space. [START_REF] Babu | A dual function of the CRISPR-Cas system in bacterial antivirus immunity and DNA repair[END_REF][START_REF] Da Silva | Exhaustive Repertoire of Druggable Cavities at Protein-Protein Interfaces of Known Three-Dimensional Structure[END_REF] A practical guide to navigate across all available methods and benchmarking data sets has been recently described. [START_REF] Ehrt | A benchmark driven guide to binding site comparison: An exhaustive evaluation using tailor-made data sets (ProSPECCTs)[END_REF] Most of above-described methods consider pocket similarity from a global and not a local point of view. The consequences are two-fold. First, the proposed initial 3D alignment of both pockets will prioritize global properties (e.g. molecular shape, principle axes and moments of inertia) over peculiar microenvironments. A wrong preliminary misalignment will not be corrected after refinement and therefore leads to erroneous similarity estimates. Second, inferring ligand information from pocket similarity searches (e.g. merging ligand coordinates from one reference pocket to a target cavity) will address the entire ligand structure as a whole, without any obvious clues about which ligand substructure ideally fits which subpocket. Therefore, most existing computational methods are well suited to repurpose existing ligands for new pockets, [START_REF] Weber | Unexpected nanomolar inhibition of carbonic anhydrase by COX-2-selective celecoxib: new pharmacological opportunities due to related binding site recognition[END_REF][START_REF] Kinnings | Drug Discovery Using Chemical Systems Biology: Repositioning the Safe Medicine Comtan to Treat Multi-Drug and Extensively Drug Resistant Tuberculosis[END_REF][START_REF] Yang | Computational discovery and experimental verification of tyrosine kinase inhibitor pazopanib for the reversal of memory and cognitive deficits in rat model neurodegeneration[END_REF] but not to prioritize ligand fragments for specific protein subsites, a very important process in fragment-based drug discovery. [START_REF] Erlanson | Twenty years on: the impact of fragments on drug discovery[END_REF] Fewer examples of subpocket comparisons are available to date. [START_REF] Schmitt | A new method to detect related function among proteins independent of sequence and fold homology[END_REF][START_REF] Ramensky | A novel approach to local similarity of protein binding sites substantially improves computational drug design results[END_REF][START_REF] Wallach | Prediction of sub-cavity binding preferences using an adaptive physicochemical structure representation[END_REF][START_REF] Durrant | CrystalDock: A Novel Approach to Fragment-Based Drug Design[END_REF][START_REF] Wood | Pharmacophore Fingerprint-Based Approach to Binding Site Subpocket Similarity and Its Application to Bioisostere Replacement[END_REF][START_REF] Jalencas | Chemoisosterism in the proteome[END_REF][START_REF] Kalliokoski | Subpocket Analysis Method for Fragment-Based Drug Discovery[END_REF][START_REF] Tang | Knowledge-based fragment binding prediction[END_REF][START_REF] Bartolowits | Considerations of Protein Subpockets in Fragment-Based Drug Design[END_REF] Existing approaches follow a common flowchart made of four steps: (i) fragmentation of protein-bound PDB ligands into smaller pieces; (ii) registration of protein-ligand non covalent interactions; (iii) definition of protein microenvironments interacting with above-reported ligand chemical moieties; (iv) mathematical representation of the microenvironment into a graph, pharmacophore or fingerprint; (v) pairwise similarity calculation between a reference and a query microenvironment.

Reported methods differ in the level of ligand fragmentation (few connected atoms, [START_REF] Ramensky | A novel approach to local similarity of protein binding sites substantially improves computational drug design results[END_REF] chemical group, [START_REF] Wallach | Prediction of sub-cavity binding preferences using an adaptive physicochemical structure representation[END_REF] fragment [START_REF] Durrant | CrystalDock: A Novel Approach to Fragment-Based Drug Design[END_REF][START_REF] Wood | Pharmacophore Fingerprint-Based Approach to Binding Site Subpocket Similarity and Its Application to Bioisostere Replacement[END_REF][START_REF] Jalencas | Chemoisosterism in the proteome[END_REF][START_REF] Kalliokoski | Subpocket Analysis Method for Fragment-Based Drug Discovery[END_REF][START_REF] Tang | Knowledge-based fragment binding prediction[END_REF] ), the atomic definition of protein microenvironments (atom [START_REF] Ramensky | A novel approach to local similarity of protein binding sites substantially improves computational drug design results[END_REF] or residue [START_REF] Durrant | CrystalDock: A Novel Approach to Fragment-Based Drug Design[END_REF] based, surface feature pseudoatoms [START_REF] Weber | Unexpected nanomolar inhibition of carbonic anhydrase by COX-2-selective celecoxib: new pharmacological opportunities due to related binding site recognition[END_REF][START_REF] Jalencas | Chemoisosterism in the proteome[END_REF][START_REF] Tang | Knowledge-based fragment binding prediction[END_REF] ), the computational representation of the subpocket (graph, 33, 36- 38 fingerprint [START_REF] Wallach | Prediction of sub-cavity binding preferences using an adaptive physicochemical structure representation[END_REF][START_REF] Tang | Knowledge-based fragment binding prediction[END_REF] ), the alignment method (clique detection, [41][START_REF] Shindyalov | Protein structure alignment by incremental combinatorial extension (CE) of the optimal path[END_REF][START_REF] Berman | The Protein Data Bank[END_REF][START_REF] Desaphy | sc-PDB-Frag: A Database of Protein-Ligand Interaction Patterns for Bioisosteric Replacements[END_REF] rigid-body transformation, [START_REF] Wallach | Prediction of sub-cavity binding preferences using an adaptive physicochemical structure representation[END_REF] rmsd alignment [START_REF] Durrant | CrystalDock: A Novel Approach to Fragment-Based Drug Design[END_REF] ) and the scoring function (simple Tanimoto or cosine metric, [START_REF] Wood | Pharmacophore Fingerprint-Based Approach to Binding Site Subpocket Similarity and Its Application to Bioisostere Replacement[END_REF][START_REF] Jalencas | Chemoisosterism in the proteome[END_REF][START_REF] Kalliokoski | Subpocket Analysis Method for Fragment-Based Drug Discovery[END_REF][START_REF] Tang | Knowledge-based fragment binding prediction[END_REF] shape and/or pharmacophore overlap, [START_REF] Ramensky | A novel approach to local similarity of protein binding sites substantially improves computational drug design results[END_REF][START_REF] Wallach | Prediction of sub-cavity binding preferences using an adaptive physicochemical structure representation[END_REF][START_REF] Kalliokoski | Subpocket Analysis Method for Fragment-Based Drug Discovery[END_REF] rmsd of key atoms [START_REF] Durrant | CrystalDock: A Novel Approach to Fragment-Based Drug Design[END_REF] ) to estimate pairwise pocket similarity. To the best of our knowledge, only retrospective validation of subpocket comparisons have been proposed, one of the most impressive being the a posteriori molecular explanation to the unexpected cross-reactivity of cyclooxygenase-2 inhibitors with human carbonic anhydrase. [START_REF] Weber | Unexpected nanomolar inhibition of carbonic anhydrase by COX-2-selective celecoxib: new pharmacological opportunities due to related binding site recognition[END_REF] Moreover, most approaches are focusing on fragment-bound sub-cavities and cannot easily predict local similarities between a novel entire cavity and a collection of microenvironments. Last, the lack of availability of most methods (KRIPO [START_REF] Wood | Pharmacophore Fingerprint-Based Approach to Binding Site Subpocket Similarity and Its Application to Bioisostere Replacement[END_REF] being a noticeable exception) hampers the usage of above-described tools.

There is therefore still a need for novel computational methods, notably those relying on novel cavity representations and alternative alignment methods, applicable at a high throughout to compare entire cavities to fragment-annotated protein microenvironment collections. Following the above guidelines, we herewith present a novel pocket comparison method (ProCare: Potein Cavity registration), particularly adapted to detect local similarity between entire cavities and fragment subpockets, that significantly differs from existing computational tools. ProCare utilizes the concept of point cloud registration, widely used in computer vision to compare and align 2D/3D images. We first describe the implementation of the method to align and compare entire cavities. After parameter optimization and fine-tuning a scoring function to evaluate pocket similarity, we then apply the new method to the comparison of fragment subpockets to full cavities, thereby enabling to fill new binding pockets with complementary fragments.

RESULTS AND DISCUSSION

In computer vision, pattern recognition, and robotics, point cloud registration [START_REF] Cheng | Registration of Laser Scanning Point Clouds: A Review[END_REF][START_REF] Chui | A new point matching algorithm for non-rigid registration[END_REF][START_REF] Jian | Robust Point Set Registration Using Gaussian Mixture Models[END_REF] is the process of finding the best spatial transformation (e.g., scaling, rotation and translation) that aligns two point clouds (Figure 1). Figure 1. Schematic representation of the point cloud registration process. The red cloud is rotated and translated along its three main axes until the optimal alignment to the green cloud is found.

Interestingly, point cloud registration has rarely been used to overlay molecular surfaces of proteins 48- 49 and ligands. [START_REF] Douguet | SENSAAS (SENsitive Surface As A Shape): utilizing open-source algorithms for 3D point cloud alignment of molecules[END_REF] With respect to a previously described approach locating point clouds on binding site atoms, [START_REF] Hoffmann | A new protein binding pocket similarity measure based on comparison of clouds of atoms in 3D: application to ligand prediction[END_REF] we here take advantage of our previous work describing a protein pocket by a point cloud located in ligand space. [START_REF] Desaphy | Comparison and druggability prediction of protein-ligand binding sites from pharmacophore-annotated cavity shapes[END_REF] The cloud is described as an ensemble of 3D points regularly filling the pocket, each point having a specific pharmacophoric property ("color") complementary to that of the nearby protein residue. [START_REF] Desaphy | Comparison and druggability prediction of protein-ligand binding sites from pharmacophore-annotated cavity shapes[END_REF] The cloud is therefore bigger (200-300 points), regular and complementary in shape and pharmacophoric properties to flanking protein residues. The registration is classically achieved in two steps. First, a rough alignment of both molecular objects is proposed using a fast point feature histogram (FPFH) of each cloud (See Computational methods) that sums up the 3D space distribution of points in a fixed-sized fingerprint. The preliminary alignment is then refined by utilizing an iterative closest point (ICP) algorithm, aiming at minimizing the root-mean square distance between pairs of points sharing the same protein microenvironment.

We will first demonstrate the proof-of-concept of applying this computational method to the problem of protein cavity alignments, next fine tune a set of parameters enabling an optimal performance on a large dataset of known cavities, and then propose a physicochemically relevant score to quantify the alignment and pocket similarity. Last, we will apply the optimized method to the specific problem of finding local similarities between fragment subpockets and entire cavities.

Proof-of-concept in applying point cloud registration to align and compare protein-ligand binding sites. A preliminary attempt on PDB copies of the same protein-ligand complexe showed that an initial registration applying the Random Sample Consensus (RANSAC) method [START_REF]How to use Random Sample Consensus model[END_REF][START_REF] Fischler | Random Sample Consensus -a Paradigm for Model-Fitting with Applications to Image-Analysis and Automated Cartography[END_REF] on standard FPFH descriptors followed by ICP refinement (nr= 2 Å, fr = 10 Å; see parameter meaning Table 3 of the Computational methods section) yielded a correct alignment of cavity points and corresponding protein-ligand coordinates to which the resulting transformation has been applied (Figure 2A). In other words, a 2 Å radius was enough to compute meaningful normal vectors to cavity points, knowing that the distance between two adjacent points is 1.5 Å. Moreover, a 10 Å radius seems large enough to capture local shape information relative to each point. Interestingly, the above-described expected success could be extended to matching subpockets of bioisosteric fragments originating from very different proteins (beta2 adrenergic receptor, cyclin-dependent kinase 2; Figure 2B), thus supporting our hypothesis that the point cloud registration might be appropriate for detecting local similarities for unrelated protein pockets. Parameter optimization. Cavity point clouds are obviously different from 2D/3D images classically used in computer vision. Although the above reported preliminary attempts were promising, we soon noticed that many parameters strongly influence the quality of the alignment. We therefore systematically enumerated 157,465 parameter combinations (Table 3, Computational methods) in order to consider the effect of 15 key parameters as well as their interdependencies. These values were chosen to be lower, equal and higher than numbers previously tested and based on the known distance between two adjacent cavity points (1.5 Å). To test all these conditions, a very simple data set of five similar pairs completed by five dissimilar cavity pairs (EASY1 set, Table 1, Computational methods) was designed, just to filter out those parameter combinations that failed in either producing any kind of alignment (fitness = 0), or could not perfectly discriminate similar from dissimilar pairs (ROC AUCs < 1, ph4-strict Tanimoto and Tversky similarity metrics). These two simple filters enabled to decrease the number of potential combinations from 157,465 to 20,181 (Figure 3). Figure 3. Selection procedure to determine the best alignment parameters. 157,465 different conditions (a set of parameters) were initially enumerated and non-relevant conditions filtered-out with the EASY1 set. The 314 remaining alignment conditions were evaluated with the BO1 set and the best one selected by its discrimination performance (high ROC AUC) and manual inspection.

For the remaining possibilities, the output transformation matrices were applied to the protein coordinates of the similar pairs to ensure whether the 3D structures were correctly aligned (rmsd on backbone heavy atoms < 2 Å) or not. A total number of 314 combinations (0.2 % of the total number) still fulfilled the above-described requirements. In order to benchmark the 314 remaining alignment conditions, we designed a larger and much more diverse data set (BO1 set, see Computational methods) of similar pairs and dissimilar pairs of cavities starting from the sc-PDB archive of 17,000 druggable-protein-ligand complexes. [START_REF] Desaphy | sc-PDB: a 3D-database of ligandable binding sites-10 years on[END_REF] Because the notion of similarity and dissimilarity of protein pockets is context-dependent, we defined two similar cavities as deriving from pairs of different proteins (different Uniprot accession numbers) that are similar in terms of sequence (50-100% identity), structure (rmsd on backbone atoms ≤ 5 Å) and functions (Uniprot keywords annotation). No constraint was applied on the bound-ligand chemical similarity, so that different cases are represented (0 ≤ chemical similarity < 1; see Computational methods for similarity calculation).Conversely, pairs of dissimilar cavities were formed from the same target space, but need be different in terms of function and bound ligands (0 ≤ chemical similarity ≤ 0.4) in order to rule out potential wrong class annotations.

The final sets of similar and dissimilar cavities have comparable distribution of size (i.e. number of points) difference between members of each pair, with the aim of eliminating possible biases in results due to alignment of differently-sized objects (Figure 4). The BO1 data set consists of 766 pairs of nonredundant VolSite cavities (383 similar pairs, 383 dissimilar pairs) covering 507 different proteins (460 in the set of similar, 178 in the set of dissimilar), 62 different sets of Uniprot functional annotations for similar pairs and 38 for dissimilar pairs. The 314 pre-selected conditions were used to align cavity pairs from the BO1 set. The area under the ROC curve (ROC AUC) of a binary classification (similar, dissimilar) was calculated to rank each condition using three possible scoring functions (ph4-strict, ph4-rules and ph4-ext; see Computational methods). We selected the best alignment condition (see parameters in Table S4) that yielded the highest ROC AUC ( ) based on the ph4-ext scoring. Indeed, ph4-ext scoring allows compensating 0.87 0.85 0.89 for slight variations in point coordinates which may cause two property-corresponding points in reference and target cavities to be in adjacent voxels, and/or allows compensating for the transformation residual errors. We further observed that for almost all the conditions, the ph4-rules scoring yielded the worst classification as it increases the probability of assigning a high score to an irrelevant match by almost 3 times, compared to ph4-strict (p ≈ 0.34 and p ≈ 0.13 respectively). We must point out that the fitness score provided by the original point cloud registration is not accounted for in our final scoring scheme since it is hard to appreciate and correctly quantify the relevance of correspondences. In any case, the property-based score is dependent on the shape matching.

From here on, the following abbreviations will be used to label combination of feature-based registration and refinement methods: FPFH-icp: FPFH-based rough registration followed by ICP refinement FPFH-colored-icp: FPFH-based rough registration followed by colored ICP refinement c-FPFH-icp: c-FPFH-based rough registration followed by ICP refinement Although the current approach was successful in aligning and ranking cavity pairs from a large and diverse data set, we observed that some pairs of similar cavities still remained misaligned (see example in Figure S1). Constraining the alignment to consider both shape and color might solve the problem. However, given our system, the colored-ICP algorithm [START_REF] Park | Colored Point Cloud Registration Revisited[END_REF] which aims to optimize both geometric (shape) and photometric (colors) terms is not suitable for two reasons: (i) ICP requires a starting point close to the optimal solution, meaning that ICP would not rescue initial FPFH feature-based misalignments;

(ii) the meaning and assignment of color in a pharmacophoric context (Table 2, Computational methods) do not correspond to that utilized in image processing (RGB primary colors). Using the optimal set of parameters on the BO1 set, but refining the rough RANSAC alignment with the FPFHcolored-icp method confirmed our initial hypothesis, as the corresponding AUC (ROC AUC = ) 0.83 0.81 0.86 was inferior to that reported above. We have therefore implemented a new descriptor to improve the correspondences estimation during the feature-based alignment.

Alignment with feature-based histograms encoding shape and pharmacophoric properties. In light of the interesting results we previously obtained with the FPFH-icp routine and regarding the misalignment issues that arose, we have modified the FPFH descriptor implemented in the Open3D library. Similarly to the way that shape information is binned to form a normalized 33-bin histogram, we encoded the distribution of eight pharmacophoric features (Table 2; The improvement of the discrimination with c-FPFH descriptors is due to the correction of alignment errors previously reported, which are consequently reflected on scores. Comparison of the ranking of pairs based on the similarity scores outputted by the 4 different methods (c-FPFH-icp, FPFH-icp, FPFHcolored-icp, and Shaper) showed the highest Spearman correlation between c-FPFH-icp and Shaper (ρ = 0.85). Differences in the ranking between methods is partially explained by misalignment of some similar pairs, and by the different fuzziness level of the utilized scoring functions. In quite a few cases, alignments of similar cavities were well approximated when evaluating the consequent alignment of the corresponding proteins, while the scores were inferior to the median score obtained for similar pairs. For those misaligned pairs, we did not find any correlation between alignment scores and chemical similarity of the cavity-bound ligands (Tversky on Morgan fingerprint and MCS uniformly ranged from 0 to 1). Another reason for misalignments is the difference in shape (globular vs. planar) observed between the two cavities, rendering neighborhood similarities of randomly sampled points difficult to catch. Of course, we cannot exclude the possibility to have wrongly annotated BO1 pairs, particularly those predicted dissimilar. However, observing a similarity between binding sites of functionally unrelated proteins is a very rare event [START_REF] Kalliokoski | Subpocket Analysis Method for Fragment-Based Drug Discovery[END_REF] so that, even if present in the data set, such errors are negligible.

Statistical evaluation of ProCare score distributions. The ability of the method combining c-FPFH descriptors for aligning and ph4-ext for scoring, was first assessed by its ability to discriminate similar and dissimilar cavities of the BO1 set, using incremental variations of the ph4-ext score (from here on ProCare score). The optimal discriminative power (recall = precision = F-measure = 0.85) is obtained at a threshold value of 0.39 for the investigated data set (Figure 6A). To check whether this threshold value is data set-dependent, we next generated a background distribution of 2.5 million alignments (510 non-redundant BO1 cavities vs. 4,223 sc-PDB cavities). 100 statistically representative samples of 100,000 values each, could be fitted to a generalized extreme value (GEV) distribution (Figure 6B) according to the Kolmogorov-Smirnov test (D = 0.046, P-value = 0.0292,  =0.02) with a probability density function of the type:

k ≠ 0 (8) f(x) = 1 σ exp( -(1 + 𝑘𝑧) -1/𝑘 ) (1 + 𝑘𝑧) -1 -1/𝑘 k = 0 f(x) = 1 σ exp( -z -exp( -z))
with k = -0.15024,  = 0.08338,  = 0.24475, 𝑧 =

𝑥 -𝜇 𝜎

The significance level p of the detected similarity represents the probability of obtaining the same or higher similarity score Z > z by chance is:

k ≠ 0 (9) 𝑝(𝑍 > 𝑧) = 1 -exp ( -(1 + 𝑘𝑧) -1 - 1 𝑘 k=0 𝑝(𝑍 > 𝑧) = 1 -exp ( -exp ( -𝑧)
From the background distribution, a statistically significant threshold for the ProCare score was set at a value of 0.47, which corresponds to a p-value of 0.05. At this threshold, the classification of the previous BO1 set yields to a lower recall (0.72) but a much better precision (0.95). From here on, ProCare will be used with the above-reported best set of parameters, combining c-FPFH descriptors for aligning and ph4-ext for scoring pocket alignments. Sensitivity of the ProCare score with variation in pocket atomic coordinates. We wanted the method needs to be insensitive to variations in atomic coordinates of the pocket, frequently observed upon ligand binding and experimental details of the structural determination method (e.g. X-ray diffraction, single-particle cryo-electron microscopy, homology modeling). We therefore designed two data sets (MD-PLA2, Apo-Holo, See Computational methods) to assess ProCare robustness to align and score identical cavities exhibiting small to large variations in atomic coordinates. In the first set (MD-PLA2), the phospholipase A2-atropine complex was subjected to a 10 ns molecular dynamics (MD) simulation in explicit water, and 1,000 MD snapshots of the atropine-bound cavity were retained for pairwise similarity calculations. The second set (Apo-Holo) is composed of 10 pairs of pockets in a ligand-bound (holo) and ligand-free (apo) form, showing from small (rmsd < 1.0 Å) to large (rmsd > 4.0 Å) variations in the atomic coordinates of cavity-lining heavy atoms.

For both sets, ProCare still detected cavity similarity up to variations in atomic coordinates located in a grey zone around 2.5-3.0 Å rms deviations for heavy atoms (Figure 6C,D). It is a challenge for binding site comparison methods to be able to catch similarity from static structures with a level of fuzziness accounting for experimentally observed variations. The usually admitted 2.0 Å rmsd in posing ligands fits with the robustness of the point cloud registration presented here. Given the right set of parameters for geometrical information, the binning setup of the histogram descriptors allows some flexibility as different angle values can be represented by a same bin. Although the RANSAC algorithm implemented checks whether the relative distances between pairs of points of a correspondence set are consistent, the regular organization of points in VolSite cavities can tolerate some inconsistencies. This can be problematic, but the positive aspect is that an alignment solution is always proposed, to be one of the best possible matches provided the key-points are correctly sampled. At the same time, the local description of each point up to a 6 Å-neighborhood allows a good precision in the proposed alignment. We therefore expect this approach to be pertinent for local alignment of sub-cavities in proteins.

Detecting similarity between fragment subpockets and entire protein cavities. As demonstrated in the previous section, point cloud registration can be successfully applied to align and compare entire protein cavities. Is it still applicable to smaller objects (fragment-binding sites), a notoriously difficult problem in cavity comparisons? [START_REF] Ehrt | Impact of Binding Site Comparisons on Medicinal Chemistry and Rational Molecular Design[END_REF] To answer this question, we systematically aligned cavity pairs from the Frag-Lig set [START_REF] Drwal | Structural Insights on Fragment Binding Mode Conservation[END_REF] (see Computational methods) in which the same protein is bound to either a druglike ligand or a substructural fragment of the later ligand (see Computational Methods). A correct subpocket to full cavity alignment can therefore be easily deduced after applying the ProCare transformation matrix to the corresponding protein-fragment complex and computing two properties:

(i) the rmsd of the fragment-bound protein to the full ligand-bound target, (ii) the similarity of interactions observed between fragment-subpocket and ligand-full cavity complexes.

Examination of pocket sizes, expressed as the number of points in the corresponding clouds, confirmed that the fragment-bound subpockets are much smaller than the entire cavities to which the corresponding full ligands bind to (Figure S2). In 91% of the cases, a structural alignment of both protein structures, performed by the combinatorial extension (CE) method, [START_REF] Shindyalov | Protein structure alignment by incremental combinatorial extension (CE) of the optimal path[END_REF] yields to a rmsd on Calpha atoms below 2 Å, illustrating that no major conformational changes occurs upon ligand binding, when compared to the original fragment-bound protein structure (Figure 7A). In this context, ProCare clearly outperforms Shaper in proposing reliable alignments (rmsd of protein backbone atoms ≤ 2 Å) in 42% of cases vs. 34% for the Gaussian-based Shaper method (Figure 7A). For those structurally wellaligned pockets, the ProCare score was higher than the previously defined threshold (score 0.47, pvalue = 0.05) in 98% of the cases, suggesting that scores obtained by aligning full cavities can be translated to the comparison of pockets of very different sizes.

We next looked whether the better alignments proposed by ProCare, corresponds to a better positioning of the fragments after rotation/translation to the full cavity. Since fragments were not always real substructures of the full drug-like ligand counterpart (but sometimes just bioisosteric substructural parts), we could not compute rms deviations on fragment atomic coordinates. We therefore estimated the similarity of interactions between the fragment subpocket and either the ProCare-aligned fragment of the native drug-like ligand, using a Tanimoto coefficient calculated on molecular interaction fingerprints (IFP). 58 Considering a conserved binding mode for IFP similarities higher than 0.6, [START_REF] Marcou | Optimizing fragment and scaffold docking by use of molecular interaction fingerprints[END_REF] the CE structural alignment indicates that the fragment binding mode is conserved in the full ligand in 53% of cases (Figure 7B).

Provided with this baseline, ProCare succeeded in correctly positioning the fragment in the full pocket in 35% of cases whereas Shaper was only successful in 28% of cases (Figure 7B), thereby confirming that the better cavity alignments provided by ProCare also translates into better poses of the corresponding fragment. In many examples, Shaper misalignments were indeed rescued by the herein described point cloud registration (Figures 7C,D).

Virtual screening of fragment subpockets to assist fragment-based drug design: a first proof-ofconcept. We next extended the concept of fragment positioning inferred from binding sites alignments, to pairs of unrelated proteins. In this fragment-based drug design exercise, we took highresolution X-ray structures of protein-ligand complexes recently disclosed for the first time in the Protein Data Bank, and check whether screening a collection of fragment subpockets for similarity to the novel query cavities (Table 4), could help reconstitute, even partly, the masked query-bound ligands. A collection of 33,953 fragment subpockets was obtained by fragmenting all sc-PDB-bound ligands (sc-PDB fragment set, Computational methods) using a previously reported protocol, [START_REF] Desaphy | sc-PDB-Frag: A Database of Protein-Ligand Interaction Patterns for Bioisosteric Replacements[END_REF] while keeping protein-bound 3D coordinates. The fragment subpocket collection was then screened for ProCare similarity to the three novel cavities whose structure had recently been disclosed and therefore not present in the sc-PDB archive. After point cloud registration, the corresponding fragments were merged into the coordinate frame of the query cavity using the optimal transformation matrix, and filtered according to two criteria: (i) compliance to the fragment rule-of-three [START_REF] Congreve | A 'rule of three' for fragment-based lead discovery?[END_REF] (hence, our fragmentation protocol may find no possible fragmentation of the sc-PDB ligand), (ii) ProCare score > 0.47. Remaining fragments hits were then ranked by a composite score (FragScore, eq.10) taking into account both pocket similarity and interaction fingerprint similarity when comparing selected fragments with the masked ligand co-crystallized with the target query.

𝐹𝑟𝑎𝑔𝑆𝑐𝑜𝑟𝑒 = 𝑃𝑟𝑜𝑐𝑎𝑟𝑒 𝑠𝑐𝑜𝑟𝑒 + 𝐼𝐹𝑃 𝑠𝑖𝑚 + 1 2 𝐼𝐹𝑃_𝑝𝑜𝑙𝑎𝑟 𝑠𝑖𝑚
where IFP sim is the similarity of full interaction fingerprints and IFP_polar sim is similarity of polar interaction fingerprints (10) The first query cavity is small-sized (335 Å 3 ) and was retrieved from the recently published muscarinic M5 receptor structure bound the tiotropium inverse agonist. [START_REF] Vuckovic | Crystal structure of the M5 muscarinic acetylcholine receptor[END_REF] It is intended to be an easy challenge since the same ligand bound to three related muscarinic receptor subtypes (M1, M3 and M4) from five sc-PDB entries. Therefore, this first query was meant as a quality control of the ProCare alignment protocol and subsequent scoring function. Hence, three tiotropium-based fragments are ranked among the top 33 th fragments (Table S5) and nicely posed with respect to the true M5-bound tiotopium pose (Figure 8A). Interestingly, highly ranked fragments derived from ligands bound to M5-bound tiotropium and proposes suitable starting points for fragment growing and/or linking. Of course, visual inspection of the merged fragments into the query cavity space remains necessary to optimize fragment hits (e.g. JH2 fragment lacks the necessary ammonium group for π-cation interaction to Tyr481) for the intended cavity. The second query cavity (681 Å [START_REF] Perot | Druggable pockets and binding site centric chemical space: a paradigm shift in drug discovery[END_REF] ) is present at the interface of an asymmetrical tumor necrosis factor-alpha (TNF-alpha) trimer. This unique inhibitorbound TNF conformation has very recently been reported [START_REF] O'connell | Small molecules that inhibit TNF signalling by stabilising an asymmetric form of the trimer[END_REF] and has no comparable structure in the sc-PDB archive. Nevertheless, several sc-PDB fragments (e.g. 4KZ0_1UJ, 3R04_UNQ; see list of top 100 scorers in Table S6) selected from unrelated proteins, appear among the top ProCare scorers, and are true bioisosteres of the benzimidazole moiety of the TNF-alpha inhibitor (Figure 8C,D). The ProCare poses of the selected fragments nicely overlaps that of the true ligand, and recapitulates aromatic interactions exhibited by the bicyclic benzimidazole ring and a hydrogen bond to Tyr151 side chain of the TNF-alpha cavity. Likewise, the disubstituted aromatic substituent of the true TNF-alpha inhibitor is also mimicked by one of the top scoring aromatic fragment (3N6U_NSU, Figure 8E).

The last query used for this preliminary proof-of-concept comes from the structure of an antagonistbound cysteinyl leukotriene type 2 receptor (CysLTR2, PDB ID 6RZ8). [START_REF] Gusach | Structural basis of ligand selectivity and disease mutations in cysteinyl leukotriene receptors[END_REF] S7). Importantly, subpocket similarity and fragment posing were found for very different reasons ranging from salt bridge mimicry to the conservation of hydrogen bonds and hydrophobic/aromatic interactions.

We acknowledge that the empirical FragScore, used in the present exercise, can only be used in case the query cavity is already filled with a ligand. It enables to retrieve either apolar/aromatic fragments exhibiting a high interaction fingerprint similarity score (IFP), or polar/charged fragments with a high polar interaction fingerprint similarity value (IFP_polar). Cavity pairwise similarity, expressed by the ProCare score remains however the main driver for fragment selection, and can be used to query cavities in the apo-state. The accompanying p-value gives a statistical support to the predictions and can be used as a surrogate to the ProCare similarity value.

CONCLUSIONS

We herewith present a novel method (ProCare), inspired from computer vision, to align and compare protein cavities. Cavities are represented as 3D point clouds annotated by pharmacophoric properties mimicking that of an ideal ligand, and aligned by the point cloud registration. Importantly, Procare takes advantage of a novel point feature histogram to encode cavity microenvironments, thereby favoring the overlay of supockets sharing similar geometrical and physicochemical properties. The new method is able to align either entire pockets, subpockets, and compares subsites to full cavities. It exhibits a comparable performance to a state-of-the-art method (shape-based similarity search) when tested across a variety of benchmarking data sets, and is relatively insensitive to moderate variations of atomic coordinates. A key feature of ProCare is its unique ability to detect local similarities and thereby compare cavities of quite different sizes (e.g. fragment-bound subpockets vs. full ligand-bound cavities). We herewith provide the proof-of-concept of its application in a fragment-based drug design scenario in which cavities from recently described X-ray structures have been compared to a collection of fragment-bound subpockets. Local similarities undetectable with standard cavity comparison tools are found by ProCare, and enable after cavity overlay, to directly locate the corresponding fragments in the query cavity. Interestingly, proposed fragments are derived from remote targets that are totally different from the query, and proved to be identical or bioisosteric to susbtructures of the unmasked query cavity-bound ligand. Of course, designing a full ligand still requires to either grow and/or link ProCare-aligned fragments with any of existing computational fragment linking tool. [START_REF] Lauri | Caveat -a Program to Facilitate the Design of Organic-Molecules[END_REF][START_REF] Maass | Recore: a fast and versatile method for scaffold hopping based on small molecule crystal structure conformations[END_REF][START_REF] Wang | Ligand based lead generation -considering chemical accessibility in rescaffolding approaches via BROOD[END_REF] Nevertheless, the novel method enables to elaborate a fragment-based drug design strategy from the simple knowledge of a cavity 3D structure, by simple detection of local similarities to a large collection of fragment-bound subpockets.

In its current implementation, ProCare can still be optimized with respect to speed and completeness.

A pairwise similarity search can be achieved in a couple of seconds, but the cpu cost could be significantly reduced by optimizing the nearest neighbor search and excluding irrelevant points in the preliminary RANSAC alignment procedure. Moreover, usage of the RANSAC algorithm does not guarantee to find the best possible solution to the registration. Deterministic algorithms able to find the absolute minimum have recently been proposed [START_REF] Yang | Go-ICP: A Globally Optimal Solution to 3D ICP Point-Set Registration[END_REF] and should be tested further on. Last, the method could also be applied to align ligands to cavity points, and propose a computer vision approach to the protein-ligand docking problem. ProCare is freely available upon request to authors.

COMPUTATIONAL METHODS

Data Sets.

EASY1 set. This data set consists of five pairs of known similar cavities and five pairs of known dissimilar cavities (Table 1). Pairing of dissimilar cavities was obtained by mixing PDB entries of different protein names and functions (2RH1-2OUZ; 2RH1-2B7Z, 2RH1-1DM2, 1DM2-2OUZ, 1DM2-2B7Z). Protein-ligand

x-ray structures were extracted from the sc-PDB database (http://bioinfo-pharma.ustrasbg.fr/scPDB) [START_REF] Desaphy | sc-PDB: a 3D-database of ligandable binding sites-10 years on[END_REF] . Cavities were computed from ligand-free sc-PDB protein input (mol2 file format)

with using default parameters of the VolSite 7 algorithm within the IChem v. 5.2.9 toolkit. [START_REF] Da Silva | IChem: A Versatile Toolkit for Detecting, Comparing, and Predicting Protein-Ligand Interactions[END_REF] Cavity points, located on a 1.5-Å three-dimensional (3D) lattice and annotated by pharmacophoric properties, [START_REF] Desaphy | Comparison and druggability prediction of protein-ligand binding sites from pharmacophore-annotated cavity shapes[END_REF] were placed within 6 Å of heavy atoms of the corresponding hidden ligand, and visually checked with Pymol v.2.1.0. 68 Keywords") and singletons were discarded. For each cluster, the proteins sequences in fasta format were retrieved from the UniprotKB API and gathered to form a multi-fasta alignment file of the cluster.

In case several isoforms were available for one protein, only the first one (default) has been considered. Then, multiple sequence alignments were performed with Clustal Omega [START_REF] Sievers | Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega[END_REF] ("biopolymer align_structure" method, default parameters). Pairs of proteins for which the root-mean square deviation (rmsd) of main chain coordinates is higher than 5 Å were discarded. For 30 pairs, a manual structural alignment was performed with Maestro v.11.9.011 [START_REF] Shrödinger | The UniProt Consortium, UniProt: a worldwide hub of protein knowledge[END_REF] to rescue SYBYL misalignments.

For each of the remaining 643 pairs, corresponding cavities were computed from the position of their bound ligands, as described above for the EASY1 set. The transformation matrix used to align the proteins was applied to their corresponding cavities using the realign module of the IChem toolkit.

Pairs of cavity points were next analyzed for co-localization, by measuring all possible pairwise distances. A pair was kept if three conditions were verified: (i) at least 45% of all pairwise distances are lower than 10 Å, (ii) any cavity point in one pair member has more than 50 unique neighbors (d < 1.5 Å) in the cognate pair member; (iii) bound ligands according to Morgan fingerprints (radius = 2) [START_REF]RDKit: Open-Source Cheminformatics Software[END_REF] were not identical (Tanimoto coefficient Tc≠ 1). Finally, a set of 383 pairs (Table S1) was annotated as "similar".

An equally-sized set of dissimilar pairs (Table S2) was defined from the above described clustering of UniprotKB keywords, as protein pairs sharing no single functional keyword and different ligands HET codes with a chemical similarity, expressed by a Tanimoto coefficient on Morgan fingerprints (radius =2), below 0.4. Finally, an equivalent number of 383 dissimilar pairs was retrieved randomly from that list, with the constraint that the distribution of differences in cavity volumes between dissimilar pairs matches that of similar pairs. MD-PLA2 set. X-ray coordinates of phospholipase A2 in complex with atropine was retrieved from the RCSB PDB (PDB ID: 2ARM). The protein-ligand complex was protonated with Protoss 74 to optimize the protonation/ionization/tautomeric state of each molecule, including protein, ligand and waters.

Atropine parameters for the GAFF force field were generated with the Antechamber routine in AMBER18, [START_REF] Case | AMBER 2018[END_REF] using AM1-BBC partial charges. The ff14SB force field was used to model protein coordinates. The complex was solvated with a 10-Å thick TIP3P water box and neutralized with 5 Clions. After energy refinement for 10,000 steps (1000 steps of steepest descent followed by 9,000 steps of conjugate-gradient minimization steps), the system was equilibrated by a 10 ps constant volume molecular dynamics (MD) simulation. The temperature was raised from 100 to 300K over 10 ps, using a weak coupling to a heat bath (temperature coupling of 2 ps) and periodic boundary conditions. The SHAKE algorithm was utilized to constrain bond lengths involving hydrogen atoms, thereby allowing a timestep of 2 fs. The equilibrium phase was continued by a 100 ps constant pressure MD with isotropic pressure scaling and a Berendsen barostat. Equilibrium was followed by a 10-ns constant pressure MD.

Throughout the different MD steps, electrostatics was treated with the Particle Mesh Ewald (PME) method. A snapshot of the MD trajectory was saved every ps, fitted to the starting conformation (heavy chain atoms) and outputted in pdb file format. The 1,000 PDB files were then converted to a mol2 file format using SYBYL-X and used as VolSite input for the atropine-bound cavity detection as described above. Atomic coordinates of the initial input protein-ligand complex (frame 1) and corresponding cavity were randomly translated by 10 Å along the three axes x-y-z and rotated by 180° along the x-axis, in order to put reference and target complexes in different coordinate frames.

Apo-Holo set. A previously reported set of 10 pairs of protein X-ray structures [START_REF] Desaphy | Comparison and druggability prediction of protein-ligand binding sites from pharmacophore-annotated cavity shapes[END_REF] [START_REF] Shrödinger | The UniProt Consortium, UniProt: a worldwide hub of protein knowledge[END_REF] Cofactors, additives and water molecules were removed. The H-bonds were optimized at pH 7 using PROPKA rules. Complexes were split into proteins and ligands and all structures were saved as mol2 files. Full cavities were computed with IChem VolSite v.5.2.9 as described before for proteinligand complexes and using only the protein as query for apostructures. For apostructures, the cavity overlapping the bound-ligand of the holostructure was further kept.

Frag-Lig set. This data set is a subset of the previously reported PDBmob data set, [START_REF] Drwal | Structural Insights on Fragment Binding Mode Conservation[END_REF] and consists of 578 pairs of cavities from the same protein (same Uniprot AC), bound to a drug-like ligand and a substructural fragment of the latter ligand. The data set provides already aligned proteinligand/fragment complexes for each target set. For each unique protein of the PDBmob data set, all possible pairs of protein-fragment and protein-druglike ligand were formed. The Tversky similarity of the paired fragments and ligands were calculated using RDKit Morgan fingerprints (radius = 2) and maximum common substructures (RDKit FindMCS default parameters). A first selection conserved pairs with both similarity metrics superior to 0.6. The corresponding cavities were computed with IChem VolSite v.5.2.9 using default parameters. For fragment-bound structures, only the close vicinity (4 Å) of the fragment was considered for cavity detection (VolSite CAVITY_4 output). For ligand-bound structures, the entire cavity (VolSite CAVITY_ALLoutput) was retrieved. This preliminary list was then filtered to remove drug-like-bound cavities of smaller volume that that of the fragment counterpart.

Then, fragment/ligand occupancy in their cognate cavities was inspected to ascertain that any heavy atom has a cavity point within a 2 Å distance. Last cavity overlap (fragment-bound vs. ligand-bound) was computed by estimating the number of fragment-bound cavity points with a close neighbor (≤ 2Å) in druglike-bound cavity points. Only pairs with 100% overlap were finally retained to yield 578 pairs (Table S3). For each pair, atomic coordinates of the fragment-protein complex were randomly translated by 10 Å along the three axes x-y-z and rotated by 180° along the x-axis, in order to put reference and target complexes in different coordinate frames.

sc-PDB fragment set.

For each of the 16,034 entries of the sc-PDB data set, [START_REF] Desaphy | sc-PDB: a 3D-database of ligandable binding sites-10 years on[END_REF] the corresponding 3D structure of the ligand was fragmented using a previously-described protocol [START_REF] Desaphy | sc-PDB-Frag: A Database of Protein-Ligand Interaction Patterns for Bioisosteric Replacements[END_REF] in three steps. First, a ring perception algorithm is used to detect aromatic and aliphatic rings of the ligand. Second, acyclic atoms are then parsed to assign either a linker or substituent label, as whether to the corresponding bonds are connecting two rings or not. Linker atoms are left unchanged. In case of substituent atoms, single bonds involving the closest apolar carbon (in terms of bond distance) to any ring are later cleaved at the condition that the cleaved bond is at least three bonds away from the cyclic root atom. Third, fragments are kept at the condition that they make at least 4 interactions (including ≥1 polar or aromatic) with the target. The fragment set contains 33,953 fragments out of which 7,294 are unique.

For each of the 33,953 protein-bound fragments, the 4 Å-surrounding cavity was computed in IChem VolSite as described above.

Point Cloud registration. The herein described method relies on Open3D v.0.5.0, 76 a library for point cloud processing. The library is available in C++ programming language but provides a python interface with pybind11, and allows parallel computing via the OpenMP environment. For the sake of efficiency, Open3D was compiled and installed from source in conda environment following the provided guidelines. Protein cavitiesfiles computed with VolSite (mol2 format) were converted into PCD (Point Cloud Data) file format version 0.7. The Header was kept as default unless the "WIDTH" and "POINT" sections that were updated with the cavity size (number of cavity points). The "DATA ascii" section contained the x, y, z coordinates of the mol2 file and a fourth column assigning a color to each of the eight VolSite pharmacophoric properties. [START_REF] Desaphy | Comparison and druggability prediction of protein-ligand binding sites from pharmacophore-annotated cavity shapes[END_REF] Normal vectors and fast point feature histograms (FPFH) [START_REF] Rusu | Fast Point Feature Histograms (FPFH) for 3D Registration[END_REF] were computed for the source cloud and the target cloud. A first rough alignment was performed based on FPFH descriptors with the Random Sample Consensus (RANSAC) method [START_REF]How to use Random Sample Consensus model[END_REF][START_REF] Fischler | Random Sample Consensus -a Paradigm for Model-Fitting with Applications to Image-Analysis and Automated Cartography[END_REF] The neighborhood is perceived without considering the points colors. Considering a point P q (green) whose FPFH is to be computed, its neighbor points = {1, 2, 3} within a radius r are determined (green circle). For each Ƥ 𝑘 neighbor in , their respective neighbors are also determined within the radius r; B) Between a point Ƥ 𝑘 and each of its neighbor, an ensemble of θ, α and φ angular values are computed to reflect the local environment of each point; C) Each of the θ, α and φ computed values for the point Pq and its normal n q are respectively binned into 11-bin histograms with regular intervals deduced from minimal and maximal distances. The resulting 33-bin histogram forms the simplified point feature histogram (SPFH) of the point Pq. Similarly, the SPFH is computed for each point in ; D) The FPFH of the point P q is the Ƥ 𝑘 sum of its SPFH and the distance-weighted average of its neighbors' SPFHs; E) Simplified schematic representation of a cloud of points with perception of point colors. Considering a point P q (green) whose c-FPFH is to be computed, its neighbor points = {1, 2, 3} within a radius r are determined Ƥ 𝑘 (green circle). For each neighbor in , their respective neighbors are also determined within radius r; Ƥ 𝑘 F) The 33-bin histogram SPFH is computed for the point P q , in addition to eight bins coding for the eight pharmacophoric features respectively, encompassing the percentage of each pharmacophoric feature in . The final 41-bin histogram forms the c-SPFH of the point P q . Similarly, the c-SPFH is computed Ƥ 𝑘 for each point in ; G) The c-FPFH of the point P q is the sum of its c-SPFH and the distance-weighted Ƥ 𝑘 average of its neighbors' c-SPFHs.

ProCare scoring. The quality of the alignment was estimated by two scores (fitness, RMSE) in Open3D.

The fitness score (eq. 1) measures the overlap of source and target clouds as the ratio of the number of inlier correspondences (points in the source cloud that are fitted to the target cloud, based on a nearest neighbor search on coordinates after transformation) to the total number of points in the source cloud.

𝑓𝑖𝑡𝑛𝑒𝑠𝑠 =

𝑁𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑖𝑛𝑙𝑖𝑒𝑟 𝑐𝑜𝑟𝑟𝑒𝑠𝑝𝑜𝑛𝑑𝑒𝑛𝑐𝑒𝑠 𝑇𝑜𝑡𝑎𝑙 𝑛𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑝𝑜𝑖𝑛𝑡𝑠 𝑖𝑛 𝑠𝑜𝑢𝑟𝑐𝑒 𝑐𝑙𝑜𝑢𝑑 (1) RMSE (eq. 2) is the root-mean square error between corresponding pairs of points in source and target clouds.

𝑅𝑀𝑆𝐸 = 1 𝑁 𝑁 ∑ 𝑖 = 1 (𝑃𝑠 𝑖 -𝑃𝑡 𝑖 ) 2 , 𝑃𝑠 ∈ Ƥ 𝑠𝑜𝑢𝑟𝑐𝑒 , 𝑃𝑡 ∈ Ƥ 𝑡𝑎𝑟𝑔𝑒𝑡 (2)
We then implemented 3 additional scoring functions (ph4-strict, ph4-rules, ph4-ext) to evaluate the alignment of pharmacophoric properties. The ph4-strict scoring method, relies on the ball-tree algorithm implemented in scikit-learn, [START_REF]scikit-learn: Machine learning in Python[END_REF]

Where c is the number of fitted points of identical pharmacophoric properties, a and b are number of points of the smallest and the largest cavity, respectively, c i is the number of points of property i aligned, f i is the average frequency of points with property i in all sc-PDB cavities.

The ph4-rules scoring function is defined as the ph4-strict, with c equal to the number of fitted points of similar pharmacophoric properties (Table 2). The ph4-ext scoring function is defined as the ph4strict, with c as the number of points in the smallest cloud which has a point of the same property of any of its neighbors in the target cloud. As for the ph4-strict scoring method, the Tanimoto, Tversky and frequency-weighted metrics are calculated. All possible combinations were tested on the EASY1 data set and their performance evaluated in three steps. First, parameter sets having rough and refined alignment fitness values higher than 0 were retrieved and their corresponding alignments were rescored with the above-described ph4-strict scoring scheme. Second, the area under the receiver operating characteristic (ROC) curve was assessed using either the Tanimoto or the Tversky metric to rank alignment similarity values. Corresponding parameter sets were conserved only if the ROC AUC was equal to 1. Finally, the target protein structures were aligned with UCSF Chimera v.1.12 80 using the cavity transformation matrix previously generated by ProCare for three EASY1 pairs (HIV protease: residues 1-99, 100-198; beta-2 adrenergic receptor: residues 1-202, 363-44; cyclin-dependent kinase 2: 2c6t-residues 1-35, 45-148; 1dm2 residues 1-35, 36-139, 140-272). Only parameter sets leading to a mean rmsd (backbone heavy atoms) below 2 Å were kept for further analysis on the BO1 data set.

Alignment of cavities with Shaper. Previously reported cavity pairs were also aligned and compared with default parameters of the previously described Shaper method. [START_REF] Desaphy | Comparison and druggability prediction of protein-ligand binding sites from pharmacophore-annotated cavity shapes[END_REF] Starting the from the same set of point clouds, Shaper relies on the OpenEye ShapeTK toolkit 41 and a smooth Gaussian function to maximize the overlap of both cavity shapes and colors (pharmacophoric properties). The alignment between cavities A and B was scored as the higher of two Tversky metrics (eq. 6-7).

𝑆 𝐴,𝐵 = 𝑂 𝐴,𝐵 0.95 𝐼 𝐴 + 0.05 𝐼 𝐵 + 𝑂 𝐴,𝐵

𝑆 𝐴,𝐵 = 𝑂 𝐴,𝐵 0.05 𝐼 𝐴 + 0.95 𝐼 𝐵 + 𝑂 𝐴,𝐵

where O A,B is the overlap between colors of cavities A and B, and I non-overlapped colors of each entity A and B.

Running times. Cavity alignments were run on a 64-bit Intel Core i5-4590 @ 3. List of BO1 dissimilar pairs; Frag-Lig set of 578 pairs of protein-fragment and related protein-ligand and complexes; Optimal Open3D parameters to align cavities from the BO1 set; Fragment hits for the muscarinic M5 receptor; Fragment hits for the TNF-alpha trimer; Fragment hits for the cysteinyl leukotriene receptor 2 (PDB ID 6RZ8).

Figure 2 .

 2 Figure 2. Proof-of-concept of applying point cloud registration to protein cavities. A) ProCare overlay of cavities from two structures of the beta2 adrenergic receptor bound to carazolol (PDB IDs: 2RH1, 5D6L), starting from different coordinate frames. The original point cloud registration of cavity points (left panel) is applied to the corresponding protein coordinates (middle panel, rmsd on backbone heavy atoms = 0.44 Å) and the bound ligands (right panel, rmsd on heavy atoms: 0.25 Å).The protein is displayed by white and cyan ribbons, respectively. The bound ligand is represented by sticks with a cpk color-coding (nitrogen, blue; oxygen, red; carbon: white and cyan, respectively); B) ProCare overlay of cavities from unrelated proteins (beta2 adrenergic receptor, PDB ID 2RH1; cyclin-dependent kinase 2, PDB ID 1DM2). The original point cloud registration of cavity points (left panel) is applied to the corresponding protein coordinates and the bound ligands (right panel, rmsd on pyrrole/imidazole heavy atoms: 0.88 Å). The protein is displayed by white and cyan ribbons, respectively. The bound ligand is represented by sticks with a cpk color-coding (nitrogen, blue; oxygen, red; carbon: white and cyan, respectively).

Figure 4 .

 4 Figure 4. Properties of the BO1 data set of 766 protein-ligand cavity pairs (338 similar, 338 dissimilar). A) Protein-bound ligand similarity (Tanimoto coefficient from Morgan fingerprints) vs. protein sequence identity (PIM of ClustalΩ alignment with default parameters) for similar (green) and dissimilar pairs (redbrick); B) Distribution of the difference in the size of cavity point clouds for similar pairs; C) Distribution of the difference in the size of cavity point clouds for dissimilar pairs.

  Computational methods) in the neighborhood of a point into an eight-bin histogram, each bin corresponding to one of the eight pharmacophoric features. The final 41-bin histogram, termed c-FPFH (Figure 9, Computational methods,) was next utilized to improve RANSAC preliminary alignments of BO1 cavity pairs. Obtained results were compared to that obtained using the standard FPFH descriptor and to the alignments obtained our previously-reported Shaper 7 tool that uses a smooth Gaussian function to optimize the shape overlap of cavity points. Using the ph4-ext scoring function to score alignment of BO1 cavity pairs, the novel c-FPFH appears clearly superior to the standard one (c-PFPH, ROC AUC= ; performance of the novel descriptor was almost similar to that obtained with the state-of-the art Shaper alignment (ROC AUC = ) on the same data set. 0.92 0.90 0.93

Figure 5 .

 5 Figure 5. Evaluation of ProCare scoring in comparing cavities from the BO1 set. A) Receiver operating characteristics (ROC) plot in ranking BO1 cavity pairs with the ph4-ext scoring function, using ProCare (standard FPFH descriptor, new c-FPFH descriptor) and Shaper; B) Distribution of ph4-ext scores after ProCare overlay with FPFH-icp refinement; C) Distribution of ph4-ext scores after ProCare overlay with c-FPFH-icp refinement; D) Distribution of scores after Shaper overlay.

Figure 6 .

 6 Figure 6. Statistical evaluation and sensitivity of ProCare to variations in atomic coordinates. A)Variation of statistical parameters (recall, precision, F-measure) of a binary classification model (similar/dissimilar) of BO1 cavity pairs for increasing ProCare similarity score thresholds; B) Fitting randomly sampled ProCare scores to a generalized extreme value (GEV) distribution. Repeated random samples (n = 100) showed to be representative of the whole population of scores (Scipy combined pvalue for the 100 Kolgomorov-Smirnov p-values with Fisher's method: 0.90). GEV parameters were estimated with EasyFit; 56 C) Sensitivity of the ProCare score to variations in atomic coordinates of the pocket (rmsd on heavy atoms to the fist snapshot), induced by molecular dynamics simulation of phospholipase A2 in complex with atropine (PDB ID: 2ARM). A score of 0.47 (dotted line) corresponds a statistically significant threshold (p-value = 0.05) to discriminate similar from dissimilar cavities; D) Sensitivity of the ProCare score to ligand-induced variations in atomic coordinates of test pockets (rmsd on heavy atoms of pocket-lining residues) of the Apo-Holo set (cell division protein kinase 2, CDK2, PDB IDs: 1DM2, 2JGZ; HIV-1 protease, HIVP1, PDB IDs: 1QBS, 1G6L; estrogen-related receptor gamma, ERRγ, PDB IDs: 2ZKC, 2ZBS; aldose reductase, AR, PDB IDs: 1ADS, 2NVD; hexokinase, Hex, PDB IDs: 2E2O, 2E2N; alginate-binding protein, ALGI, PDB IDs:1Y3N, 1Y3Q; Osmo-protection protein, OSMO, PDB IDs: 1SW2, 1SW5; D-allose binding protein, ALLO, PDB IDs:1RPJ, 1GUD; guanylate kinase, GUA, PDB IDs: 1EX7, 1EX6; 5-enolpyruvylshikimate-3-phosphate synthase, ESP, PDB IDs:1RF4, 1RF5). A score of 0.47 (dotted line) corresponds a statistically significant threshold (p-value = 0.05) to discriminate similar from dissimilar cavities.

Figure 7 .

 7 Figure 7. Evaluation of ProCare alignment of fragment supockets to full cavities. A) Proportion of pairs of proteins poses yielding rmsd on main chain atoms falling into the following intervals (Å) [0;2[, [2;4[, [4;6[, [6;10[, [10;∞[ after applying the transformation matrix derived from ProCare and Shaper alignments. The values were compared to the original structural alignments of the proteins obtained by the CE algorithm; 42 B) Proportion of pairs of fragment poses yielding IFP similarity with their paired ligands which falls into the following intervals [0;0.2], ]0.2;0.4], ]0.4;0.6], ]0.6;0.8], ]0.8;1.0]; C) Example of Shaper misalignment of cavities from cytochrome P121 bound to fragment 1G9 (PDB ID 4IQ7) and ligand YTT (PDB ID 3G5H; rmsd of proteins backbone heavy atoms: 22 Å; rmsd of ligands matching substructure: 5.4 Å); D) ProCare correct alignment of the same cavity pair (rmsd of proteins backbone heavy atoms: 0.45 Å; rmsd of ligands matching substructure: 0.59 Å).

Figure 8 .

 8 Figure 8. ProCare positioning of sc-PB fragments in novel cavities. Atoms are colored using a cpk colorcoding (nitrogen: blue; oxygen: red; sulfur; yellow; carbon of fragment: cyan/rosy salmon, green; carbon of true ligand, white). A) Placing a fragment derived from a muscarinic M1 receptor-bound ligand (PDB ID: 5CXV; HET: 0HK, rank = 1, , FragScore = 1.61, ProCare similarity = 0.82 [p-value = 2.04e-

  Again, this structure has no similar homologue in the sc-PDB archive, such that the ProCare search for potential subpocket matching has no obvious bias. The CysLTR2 pocket is wider (813 Å 3 ) than the two previous ones, and is fully occupied by a high molecular weight ligand (ONO-2080365, HET: KNZ) filling three separate subsites, thereby challenging ProCare for finding local similarity to each of the three subpockets and finding appropriate fragments. The benzoxazine dicarboxylic acid-binding subpocket in CysLTR2 is found similar to that of two adenosine-3',5'-cyclic-monophosphate (cAMP) pockets from unrelated proteins (catabolite gene activator protein, PDB ID: 1RUO; Potassium/sodium hyperpolarizationactivated cyclic nucleotide-gated channel 2, PDB ID: 1Q43) with the cyclic phosphate group mimicking each of the two carboxylic acids of the CysLTR2 antagonist (Figure8F) and interacting with a basic residue (Arg82 for 1RUO, Arg591 for 1Q43) that drives the subpocket similarity to the CysLTR2 cavity (FigureS3). Local similarity to the central phenoxy-binding subsite is also found in a subpocket from a receptor tyrosine phosphatase (PDB ID: 3QCH, Figure8F) with a nice overlap of the corresponding dichlorophenyl fragment to the fluorophenyl substructure of the CysLTR2 ligand. Another fragment mimicking both the benzoxazine and the central fluophenyl CysLTR2 antagonist is selected by ProCare from remote pocket similarity to that of a glucokinase pocket (PDB ID 3F9M, Figure8G). Last the hydrophobic CysLTR2 subsite accommodating the terminal difluorophenyl ring of the bound inhibitor is found similar to that of a MAP kinase 14 subpocket (PDB ID: 3DT1) with a nice overlap of the cognate phenyl fragment to the terminal aromatic ring of the CysLTR2 ligand (Figure8G). Altogether, ProCare managed to find subpocket similarity between each of the three CysLTR2 subsites with totally unrelated subpockets and proposes reliable fragments for a structure-based fragment linking strategy (see the list of 100 top fragments in Table

  via the EMBL-EBI web services REST API[START_REF] Madeira | The EMBL-EBI search and sequence analysis tools APIs in 2019[END_REF] using default parameters, and outputted in ClustalW format. The Percent Identity Matrix (PIM) files were processed to retrieve pairs of proteins having different Uniprot AC and a sequence identity between 50 and 100%. For enzymes (Function-Keywords containing one of the 6 enzyme classes), the Enzyme Classification (E.C.) number was fetched from UniprotKB and additional filtering was performed to discard pairs having different E.C. numbers and pairs in which at least one partner is not annotated with E.C. number (e.g. TrEMBL entries). At this stage, PDB atomic coordinates of ligand-bound protein chains were extracted and structurally aligned with Sybyl-X 2.1.172

  in an iterative way (registration_ransac_based_on_feature_matching function). The rough alignment was subsequently refined with an Iterative Closest Point algorithm[START_REF] Besl | A Method for Registration of 3-D Shapes[END_REF] (registration_icp function) starting from the transformation matrix of the rough alignment. Alternative to registration_icp is registration_colored_icp, which is a function considering the color of points to compute transformation matrices. We further implemented a new descriptor, the colored-FPFH (c-FPFH). c-FPFH consists of 41 bins: the 33 FPFH bins, with eight additional normalized bins accounting for the distribution of the eight colors (pharmacophoric properties) in the neighborhood of the point (Figure9).

Figure 9 .

 9 Figure 9. Fast point feature histogram (FPFH) and colored fast point feature histogram (c-FPFH) computation. A) Simplified schematic representation of a cloud of points. The neighborhood is perceived without considering the points colors. Considering a point P q (green) whose FPFH is to be computed, its neighbor points = {1, 2, 3} within a radius r are determined (green circle). For each Ƥ 𝑘 neighbor in , their respective neighbors are also determined within the radius r; B) Between a point Ƥ 𝑘 and each of its neighbor, an ensemble of θ, α and φ angular values are computed to reflect the local environment of each point; C) Each of the θ, α and φ computed values for the point Pq and its normal

,

  and searches for the nearest neighbor point in the largest cavity for each point of the smallest cavity, within a maximum distance d ( d = 1.5 Å by default). Three similarity indices (Tanimoto, Tversky, Wei) are computed for each alignment (eq. 3-5). (α=0.95, β=0.05) 𝑇𝑣𝑒𝑟𝑠𝑘𝑦(𝛼,𝛽) = 𝑐 𝛼(𝑎 -𝑐) + 𝛽(𝑏 -𝑐) + 𝑐 {𝐶𝐴, 𝐶𝑍, 𝑂,𝑁,𝑂𝐷1,𝑁𝑍,𝑂𝐺,𝐷𝑈}

Table 4 .

 4 Binding site comparison of three protein-ligand complexes recently released in the PDB.

	Target	PDB ID Ligand a Resolution, Å	Release date Cavity size b
	M5 muscarinic receptor	6OL9	0HK	2.5	2019-12-11	99
	TNF-alpha trimer	6OOY	A7M	2.5	2019-12-25	208
	Cysteinyl leukotriene receptor 2 6RZ8	KNZ	2.7	2019-12-11	241

a Ligand chemical component HET code b number of cavity points. The volume of cavity (in Å

[START_REF] Perot | Druggable pockets and binding site centric chemical space: a paradigm shift in drug discovery[END_REF]

) is the number of points x 3.375 (third power of the grid resolution in Å)

Table 1 :

 1 Entries of the EASY1 data set Starting from all 16,034 sc-PDB protein-ligand complexes, unique proteins were retrieved and clustered according to UniProt 69 keywords. Proteins without keywords (cluster "No

	PDB ID a	HET code b	Protein name	Family
	2RH1	CAU	Beta-2 adrenergic receptor	G protein-coupled receptor
	5D6L	CAU	Beta-2 adrenergic receptor	G protein-coupled receptor
	2OUZ	CD3	Estrogen receptor α	Nuclear hormone receptor
	3ERT	OHT	Estrogen receptor α	Nuclear hormone receptor
	2C6T	DT5	Cyclin-dependent kinase 2	Protein kinase
	1DM2	HMD	Cyclin-dependent kinase 2	Protein kinase
	1C6X	3IN	HIV-1 protease	Protease
	2B7Z	MKA	HIV-1 protease	Protease
	1FTL	DNQ	Glutamate receptor 2	Ion channel
	1LB9	DNQ	Glutamate receptor 2	Ion channel
	a Protein Data Bank (PDB) 43 identifier

b Chemical component identifier (https://www.rcsb.org/ligand/XXX, XXX= HET code) BO1 diverse set.

Table 2 :

 2 Pharmacophoric matching rules used by the ph4-rules scoring function.ProCare parameters optimization.A set of values were rationally defined for 15 Open3D parameters (Table3). A combination of these values led to 157,464 different alignment conditions.

	Property	Definition	Compatible pharmacophoric properties
	CA	Hydrophobic	CA, CZ
	CZ	Aromatic	CZ, CA
	N	H-bond donor	N, NZ, OG
	NZ	Positive	NZ, N, OG
	O	H-bond acceptor	O, OD1, OG
	OD1	Negative	OD1, O, OG
	OG	H-bond acceptor & donor	OG, N, O, OD1, NZ
	DU	Dummy atom	DU

Table 3 :

 3 Open3D parameters values for ProCare alignment (default values are underlined)

	Parameter

  30 GHz processor with 4 threads,16 Go RAM. Average running times are 2.17 and 0.16 s for ProCare and Shaper, respectively. Statistical analysis. Data analysis was performed with in-house python scripts. The 90 % confidence calculated with the NumPy 81 package to be the 95 th and the 5 th percentiles. Sampling fitting to the generalized extreme value (GEV) distribution and statistical tests were performed with EasyFit[START_REF]EasyFit version 5.6, MathWave Technologies[END_REF] and The supporting information is available free of charge on the ACS Publications website at DOI:xxx Example of misalignment for a pair of similar cavities from the BO1 set; Distribution of pocket size for fragments and full ligands; ProCare overlay of cavities from unrelated targets; List of BO1 similar pairs;

	intervals for	𝐴𝑈𝐶 𝑖 𝑢𝑝𝑝𝑒𝑟 𝑖 𝑙𝑜𝑤𝑒𝑟	were obtained with 5,000 bootstrap samples, where	𝑖 𝑢𝑝𝑝𝑒𝑟	and	𝑖 𝑙𝑜𝑤𝑒𝑟	were
	Scipy. 82
	ASSOCIATED CONTENT
	Supporting Information.

ACS Paragon Plus EnvironmentJournal of Medicinal Chemistry

This work was funded by a grant of the Ministère de l'Enseignement supérieur, de la Recherche et de l'Innovation to M.E. The Calculation Center of the IN2P3 (CNRS, Villeurbanne, France) is acknowledged for allocation of computing time and excellent support. We sincerely thank Prof. E. Kellenberger and Guillaume Bret (UMR7200, University of Strasbourg) as well as Prof. F. Hetroy-Wheeler (ICube, University of Strasbourg, France) for helpful discussions. ABBREVIATIONS USED 3D, three-dimensional; PDB, Protein Data Bank; rmsd, root-mean square deviation; RANSAC, Random Sample Consensus; ICP, iterative closest point algorithm; FPFH, Fast Point Feature Histogram.

AUTHOR INFORMATION

Notes

The authors declare no competing financial interest

ACKNOWLEDGMENTS Table of Contents Graphics

A computer vision approach to align and compare protein cavities: Application to fragment based drug design.

Merveille Eguida and Didier Rognan


Catalina Torres

Arianna Caporali

Gilles Pison

The Human Multiple Births Database (HMBD) An international, open access database on twin and other multiple births

Keywords: multiple births, twin births, twinning rate, demographic data, database

The twinning rate has increased dramatically over the last four decades in many countries. This increase constitutes a public health challenge because multiple babies are more fragile than singleton ones and are exposed to higher infant mortality risks. Thus, following the trends and the spatial variations of twin and other multiple deliveries is valuable and relevant to public health. The Human Multiple Births Database (HMBD) facilitates access to statistics on the annual number of deliveries by multiplicity (singletons, twins, triplets, quadruplets and more) and the annual twinning and multiple rates by country. As of May 2021, the database comprises 14 countries with reliable vital statistics. The data series go back in time as far as possible, depending on the availability of statistics; they start in the 19 th century for some countries, until the most recent year available. In the short term, the database will comprise approximately 30 countries. In the long term, we aim to extend its geographical coverage, as the quality of vital statistics improves in other regions of the world. The HMBD is designed for anyone interested in the topic of twin and other multiple births. All data, metadata and other material are freely available at https://www.twinbirths.org

Introduction

The twinning rate (defined as the annual number of twin deliveries divided by the total number of deliveries registered in a country, per 1,000) has increased dramatically over the last four decades in nearly all countries for which we have vital statistics information [START_REF] Monden | Twin Peaks: more twinning in humans than ever before[END_REF][START_REF] Pison | Twinning Rates in Developed Countries: Trends and Explanations[END_REF][START_REF] Hoekstra | Dizygotic twinning[END_REF][START_REF] Pison | Frequency of twin births in developed countries[END_REF]. For example, it has increased in the United States from 9.5 twin deliveries per 1,000 deliveries in 1975 to 16.7 in 2015 [START_REF] Martin | Is Twin Childbearing on the Decline? Twin Births in the United States, 2014-2018[END_REF]. Over the same period, it has practically doubled in many developed countries, e.g. from 9.9 to 16.1 in England and Wales, from 9.2 to 18.4 in Germany, from 9.3 to 17.5 in France, from 9.6 to 16.7 in Denmark, and from 5.9 to 9.9 in Japan (data based on national vital statistics). This significant increase constitutes an important public health challenge because multiple babies are more fragile than singleton ones. Compared to singletons, they tend to have lower birth weight, more complications at birth, and are more often born premature, which are all associated with many longterm health implications [START_REF] Johnson | Early-life origins of adult disease: national longitudinal population-based study of the United States[END_REF][START_REF] Delobel-Ayoub | Behavioral problems and cognitive performance at 5 years of age after very preterm birth: the EPIPAGE Study[END_REF][START_REF] Larroque | Survival of very preterm infants: Epipage, a population based cohort study[END_REF]. Stillbirth and infant mortality rates are also much higher among twins and other multiple children [START_REF] Pison | Twins in Sub-Saharan Africa: Frequency, Social Status and Mortality, in Mortality and society in Sub-Saharan Africa[END_REF][START_REF] Guo | Child mortality among twins in less developed countries[END_REF][START_REF] Monden | Mortality among twins and singletons in sub-Saharan Africa between 1995 and 2014: a pooled analysis of data from 90 Demographic and Health Surveys in 30 countries[END_REF][START_REF] Fresson | Analyse de la mortinatalité par le Programme de médicalisation du système d'information (PMSI): premiers résultats nationaux[END_REF]. Multiple-child pregnancies can have negative impacts on the parents as well. For instance, these pregnancies are associated with higher risk of gestational diabetes, pre-eclampsia, post-partum depression, and increased divorce risk [START_REF] Jena | Association between the birth of twins and parental divorce[END_REF][START_REF] Choi | Multiple births are a risk factor for postpartum maternal depressive symptoms[END_REF][START_REF] Rauh-Hain | Risk for developing gestational diabetes in women with twin pregnancies[END_REF][START_REF] Bdolah | Twin pregnancy and the risk of preeclampsia: bigger placenta or relative ischemia?[END_REF]. A better understanding of the variations and trends in twinning (and other multiple) rates is therefore urgently needed.

The Human Multiple Births Database (HMBD) facilitates access to aggregate, demographic data on the births of twins and other multiple children. The database provides annual data-series at the national level on the number of deliveries by multiplicity (singletons, twins, triplets, quadruplets and more); the annual twinning and multiple rates are also provided.

The HMBD is the continuation of a database on multiple births created by Gilles Pison, with the collaboration of Agata D'Addato [START_REF] Pison | Frequency of twin births in developed countries[END_REF], as well as Christiaan Monden and Jeroen Smits [START_REF] Monden | Twin Peaks: more twinning in humans than ever before[END_REF][START_REF] Pison | Twinning Rates in Developed Countries: Trends and Explanations[END_REF]. Formerly available upon request, this dataset has been updated, extended back in time, harmonized, and made freely available online as the HMBD, available at: https://www.twinbirths.org

This website gives open access to the full dataset, country-specific input data and metadata files, country-specific R-codes, interactive and user-friendly web applications that facilitate a rapid visualisation of the database (for example, the HMBD Data explorer), a glossary of relevant terms regarding demographic data on multiple births, explanations of key points, and a list of references where further, more disaggregated data on multiple births can be found.

The HMBD is not only for researchers but also for anyone interested in multiple deliveries (e.g. NGOs and associations of parents of multiple children). The HMBD is funded by the French National Research Agency (Agence nationale de la recherche, contract ANR -18-CE36-0007-01). The members of the operations team are affiliated with the French Museum of Natural History and the French Institute for Demographic Studies (INED).

Data collected: coverage and content

The input data for the HMBD come from the vital statistics of each country included. The data cover total national populations (unless otherwise indicated in the metadata). The main criteria for including a country in the HMBD are that its civil registration is complete (or nearly complete) and that vital statistics are regularly published, with births specified by multiplicity. As of May 2021, the online database comprises 14 national populations (Figure 1). For each country, the period covered extends from the earliest until the most recent year, using available annual vital statistics on births. In the short term, the database will comprise approximately 30 countries. In the long term, we aim to extend its geographical coverage, given the improvement of vital statistics in other regions of the world.

The HMBD provides aggregate, annual counts and rates of deliveries by multiplicity, i.e. single and multiple deliveries, specifying, whenever possible, the number of children among the latter (twin deliveries, triplet deliveries, etc.). Twin deliveries are by far the most frequent among multiple deliveries; deliveries of triplets or more children represent only a small fraction of the latter (for example, in France in 2010-2015, 97.9% of all multiple deliveries were those involving twins; the proportion of multiple deliveries of triplets and quadruplets and more was respectively 2.0% and 0.1%). Consequently, the HMBD focuses on twin deliveries, although it does provide information on other multiple deliveries whenever available. For each country, the database offers the most complete possible annual series of the variables shown in Table 1. Data on multiple deliveries by certain characteristics (such as the sex-composition of the newborns, the age of the mother, etc.) may also be available for some countries over specific periods. As these details are rare, they are not included in the annual series provided in the HMBD. Nevertheless, a list of sources where such information can be found is provided on the HMDB website. However, some factors may slightly compromise data comparability, especially the treatment of stillbirths in the statistics. For example, for a given country, stillbirths may be included in data for recent years, but not in the historical data. In cases of twin births involving only one live birth, the delivery may be counted as a single delivery or as a twin delivery, depending on the definitions and the reporting practices in each country. If both children are stillborn, the delivery may or may not be counted. Moreover, the definition of what counts as a stillbirth may vary over time within the same country, or between countries. All these issues are fully documented in the country-specific metadata documents, which contain relevant notes and definitions regarding the treatment of the stillbirths in the statistics, a crucial issue because stillbirths are more common in pregnancies involving multiple children. Additionally, in the database, the variable Stillbirths indicates whether these are included in the data for a given country and year. Other factors that may affect data coverage (e.g. territorial changes and other country-specific historical processes) are also explained in the metadata.

Main research results and current projects using the HMBD

Using the HMBD complemented with data from other sources, Monden et al. [START_REF] Monden | Twin Peaks: more twinning in humans than ever before[END_REF] found that since the 1980s, the global twinning rate has increased by a third, from 9.1 to 12.0 twin deliveries per 1,000 deliveries. In 2010-2015, the number of twin deliveries was higher than ever, both at the world level and for all global regions (except South America). While the global deliveries have increased by only 8%, the number of twin deliveries has increased by 42%.

Using a precursor of the HMBD, Pison et al. (2015) [START_REF] Pison | Twinning Rates in Developed Countries: Trends and Explanations[END_REF] showed an unprecedented and rapid increase in the twinning rate from the mid-1970s in various developed countries. Figure 2 shows this pattern, using the updated data from the HMBD. Two factors are mainly responsible for this increase: the delay in childbearing and the rise in medically assisted reproduction (MAR). MAR, which has increased substantially since the 1970s [START_REF] Wyns | ART in Europe, 2016: results generated from European registries by ESHRE[END_REF], explains between 22% and 87% of the total change over the period 1970-2005, depending on the country [START_REF] Pison | Twinning Rates in Developed Countries: Trends and Explanations[END_REF].

The increase in twinning rates due to MAR has raised concerns among governments and medical authorities, which has led to a change in MAR regulations and practices in most developed countries. For instance, the mean number of embryos transferred has been reduced [START_REF] Wyns | ART in Europe, 2016: results generated from European registries by ESHRE[END_REF]. These changes are probably responsible for the finding that in 2015, in about a quarter of the countries studied, the twinning rates reached a plateau in the early 2000s and decreased afterwards [START_REF] Hoekstra | Dizygotic twinning[END_REF]. This recent development can be observed using the updated data of the HMBD (Figure 2).

Ongoing specific data analysis projects using the HMBD include documenting recent levels and trends in twinning rates in different countries and attempting to disentangle the effects of various factors in the twinning rate by examining two specific questions: (a) Age/parity/fertility control. Does parity really influence twinning, or is it just an artefact? If so, what mechanism is responsible for the increase in twinning with higher birth orders? and (b) Delayed childbearing / medically assisted reproduction. What share of the recent increase in twinning rates in developed countries is due to each factor? Is this twin boom ending?

Figure 1 .

 1 Figure 1. Coverage of the HMBD according to the availability of twinning rates, by country

Figure 2 .

 2 Figure 2. Twinning rates in selected countries

Table 1 .

 1 Variables included in the Human Multiple Births Database (HMBD)

	Variable name	Description
	Country	Country name
	Source	Source name (abbreviation or code)
	Year	Year of reference of the data
	Stillbirths	Indicates whether stillborn children are counted or not
	Singletons	Number of single deliveries
	Twin_deliveries	Number of twin deliveries
	Triplet_deliveries*	Number of triplet deliveries
	Quadruplet_plus_deliveries*	Number of deliveries involving four or more children
	Multiple_deliveries	Total number of multiple deliveries
	Multiple_children	Number of children born from multiple deliveries (twins, triplets, etc.)
	Total_deliveries	Total deliveries (numbers of single and multiple deliveries combined)
	Twinning_rate	Number of twin deliveries / total deliveries, per 1000.
	Multiple_rate	Number of multiple deliveries / total deliveries, per 1000.

twinbirths.org/en/data-metadata/ input data files and R-codes can be downloaded at https://www.twinbirths.org/en/data-metadata/