Merveille Eguida

Christel Schmitt-Valencia

Marcel Hibert

Pascal Villa

Target-focused library design by pocket-applied computer vision and fragment deep generative linking

We here describe a computational approach (POEM: Pocket Oriented Elaboration of Molecules) to drive the generation of target-focused libraries while taking advantage of all publicly available structural information on protein-ligand complexes. A collection of 31 384 PDB-derived images with key shapes and pharmacophoric properties, describing fragment-bound microenvironments, is first aligned to the query target cavity by a computer vision method. The fragments of the most similar PDB subpockets are then directly positioned in the query cavity using the corresponding image transformation matrices. Lastly, suitable connectable atoms of oriented fragment pairs are linked by a deep generative model to yield fully connected molecules. POEM was applied to generate a library of 1.5 million potential cyclin-dependent kinase 8 inhibitors. By synthesizing and testing as few as 43 compounds, a few nanomolar inhibitors were quickly obtained with limited resources in just two iterative cycles.

Abs trac t

INTRODUCTION

Fragment-based drug design (FBDD) [START_REF] Erlanson | Twenty Years On: The Impact of Fragments on Drug Discovery[END_REF] has gain considerable popularity in the last 20 years for identifying new lead compounds and guiding the optimization towards drug candidates, even up to the market with four recently approved drugs. [START_REF] Li | Application of Fragment-Based Drug Discovery to Versatile Targets[END_REF] Common FBDD programs starts by screening libraries of low molecular weight compound (fragments) [START_REF] Troelsen | Library Design Strategies to Accelerate Fragment-Based Drug Discovery[END_REF] by multiple biophysical methods such as nuclear magnetic resonance spectroscopy (NMR), surface plasmon resonance (SPR), isothermal titration calorimetry (ITC) or mass spectroscopy (MS) to cite just a few. [START_REF] Coyle | Applied Biophysical Methods in Fragment-Based Drug Discovery[END_REF] Key advantages of FBDD with respect to biochemical high-throughput screening (HTS) are the sampling of a much larger chemical space as well as higher hit rates, even for difficult targets for which other approaches failed. Despite low affinities, fragment hits can be progressed to leads by linking, merging or growing approaches. [START_REF] Erlanson | Fragment-Based Drug Discovery[END_REF] Although not necessary, it is usually advisable to start from high quality X-ray diffraction data to position fragment hits in their cognate target. [START_REF] Erlanson | Fragment-Based Drug Discovery: Advancing Fragments in the Absence of Crystal Structures[END_REF] Even if FDBB is now widely used for hit identification, not all targets and fragments are suitable to X-ray diffraction. One the one hand, some targets still proved to be hard to isolate, purify in large scale and produce high-quality crystals for X-ray diffraction.

On the other hand, some fragments cannot be detected by the latter technique because of poor physicochemical properties or too low affinities. In such cases, computational approaches are the only alternatives to predict the most viable positions of fragment hits identified experimentally [START_REF] Bian | Computational Fragment-Based Drug Design: Current Trends, Strategies, and Applications[END_REF] or to identify new hits by in silico screening. [START_REF] De Graaf | Crystal Structure-Based Virtual Screening for Fragment-Like Ligands of the Human Histamine H(1) Receptor[END_REF] Three computational approaches can be used to predict the relative orientation of a fragment in a target cavity: molecular docking, functional group mapping and deconstruction-reconstruction. Molecular docking [START_REF] Brooijmans | Molecular Recognition and Docking Algorithms[END_REF] is by far the most popular structure-based approach and aims at identifying both the bound conformation and the orientation of the ligand in a target cavity from their respective stereochemical and topological complementarities. Although it has mostly been applied to drug-like compounds, docking can be used to pose fragments with an accuracy comparable to that of lead-like compounds. [START_REF] Sandor | Virtual Fragment Docking by Glide: A Validation Study on 190 Protein-Fragment Complexes[END_REF][START_REF] Verdonk | Docking Performance of Fragments and Druglike Compounds[END_REF] Docking is the computational method that is the closest to experimental fragment screening, and can be directly applied to any fragment library. In addition to potential hit identification, the fragment position in the target cavity is also predicted. Unfortunately, scoring weak-binding fragments remains a challenge and requires an efficient post-processing, e.g. knowledge-based protein-ligand interaction rescoring. [START_REF] Marcou | Optimizing Fragment and Scaffold Docking by Use of Molecular Interaction Fingerprints[END_REF][START_REF] Jacquemard | Binding Mode Information Improves Fragment Docking[END_REF][START_REF] Chachulski | Leads-Frag: A Benchmark Data Set for Assessment of Fragment Docking Performance[END_REF] Functional group mapping [START_REF] Guvench | Computational Functional Group Mapping for Drug Discovery[END_REF] uses probe atoms or groups to map a protein cavity at their preferential location. Probes can be positioned according to protein-ligand interaction energies at regular points of a three-dimensional (3D) lattice [START_REF] Kozakov | The Ftmap Family of Web Servers for Determining and Characterizing Ligand-Binding Hot Spots of Proteins[END_REF][START_REF] Radoux | Identifying Interactions That Determine Fragment Binding at Protein Hotspots[END_REF] or by molecular dynamics (MD) sampling. [START_REF] Guvench | Computational Fragment-Based Binding Site Identification by Ligand Competitive Saturation[END_REF] Interestingly, exhaustive all-atom MD better captures protein flexibility and solvation issues, and may also unmask transient cavities hidden to conventional docking protocol. Key drawback is the computational burden limiting a wide applicability for virtual screening. Moreover, reconstructing a fully connected ligand from several discontinuous propensity maps is not straightforward.

Last, deconstruction-reconstruction approaches [START_REF] Chen | Evolutions in Fragment-Based Drug Design: The Deconstruction-Reconstruction Approach[END_REF] aim at computationally splitting protein-bound ligand X-ray structures into fragments according to well-known retrosynthetic organic chemistry rules. [START_REF] Lewell | Recap--Retrosynthetic Combinatorial Analysis Procedure: A Powerful New Technique for Identifying Privileged Molecular Fragments with Useful Applications in Combinatorial Chemistry[END_REF][START_REF] Degen | On the Art of Compiling and Using 'Drug-Like' Chemical Fragment Spaces[END_REF] Resulting fragments can then be recombined into new chemical entities while taking into account the protein environment. The method still suffers from the tricky recombination step (linking, merging, scaffold hopping) [START_REF] Maass | Recore: A Fast and Versatile Method for Scaffold Hopping Based on Small Molecule Crystal Structure Conformations[END_REF] that may disturb the original fragment binding modes or generate conformational strains. Interestingly, deep generative models [START_REF] Imrie | Deep Generative Models for 3d Linker Design[END_REF][START_REF] Yang | Syntalinker: Automatic Fragment Linking with Deep Conditional Transformer Neural Networks[END_REF][START_REF] Imrie | Deep Generative Design with 3d Pharmacophoric Constraints[END_REF] for linking disconnected fragments have shown some promises as they learn from millions of existing bioactive ligands. Deconstructionreconstruction is mainly target-specific and applicable to targets for which numerous co-crystallized ligands are already available, although docking poses may be used in principle. None of the above-reported method really takes profit of the increasing amount of structural data on protein-ligand complexes and their druggable pockets. [START_REF] Bhagavat | An Augmented Pocketome: Detection and Analysis of Small-Molecule Binding Pockets in Proteins of Known 3d Structure[END_REF] Since low molecular weight fragments have been shown to bind to preferential protein microenvironments regardless of their evolutionary relationship, [START_REF] Gao | A Comprehensive Survey of Small-Molecule Binding Pockets in Proteins[END_REF] a FBDD approach considering the whole universe of druggable ligands and pockets is desired. Capitalizing on our recent numerical image processing tool to describe and align protein cavities, [START_REF] Eguida | A Computer Vision Approach to Align and Compare Protein Cavities: Application to Fragment-Based Drug Design[END_REF] we here propose to pose fragments according to the local similarity of their respective subpockets to the target cavity. Applying the transformation matrix leading to the optimal subpocketcavity alignment, the corresponding fragments are directly positioned into the target cavity and connected, under topological constraints, by a deep generative linker to yield fully connected molecules. Applying the method to the catalytic site of human cyclin dependent kinase 8 (CDK8), a focused library of 1.5 million chemical entities could be quickly generated. Interestingly, most newly generated compounds exhibited unprecedented structures. In vitro biological evaluation of 43 carefully selected compounds identified several nanomolar inhibitors within just two design iterations and limited experimental efforts.

RESULTS AND DISCUSSION

Setting the scene. We herein present a novel method to design target cavity-focused libraries based on predicted similarities between the target cavity and a library of PDB fragment-bound subpockets (Figure 1). The underlying idea is to locate the most complementary fragments in the target cavity based on the estimated similarity of their corresponding subpockets, and then to link the prepositioned fragments into drug-like compounds using a deep generative linker. Accordingly, this approach can be implemented even in the absence of known ligands for the target protein. To assess its applicability and limits in a real-life drug design project, the method is here applied to CDK8, a target of pharmaceutical interest [START_REF] Dale | A Selective Chemical Probe for Exploring the Role of Cdk8 and Cdk19 in Human Disease[END_REF] and known X-ray structure. [START_REF] Mallinger | Discovery of Potent, Selective, and Orally Bioavailable Small-Molecule Modulators of the Mediator Complex-Associated Kinases Cdk8 and Cdk19[END_REF] In the following sections, we will describe, step by step, each part of the workflow until the experimental validation of newly generated inhibitors. Alignment of fragments to the target cavity. Subpockets, defined as the immediate protein environment around bound fragments of druggable protein-ligand complexes (sc-PDB dataset), [START_REF] Desaphy | Sc-Pdb: A 3d-Database of Ligandable Binding Sites--10 Years On[END_REF] were compared and aligned to the ATP pocket of CDK8 with the aim to use the hidden bound fragments for library design. The rationale of this implementation is that according to the similarity principle, fragments originating from similar subpockets are likely to reproduce favorable interactions with the target pocket. The term 'fragment' here refers to the molecular moieties obtained after interactionaware 3D fragmentation of ligands bound to proteins so that each fragment exhibits at least one polar interaction and at least four interactions with its target. [START_REF] Desaphy | Sc-Pdb-Frag: A Database of Protein-Ligand Interaction Patterns for Bioisosteric Replacements[END_REF] The query CDK8 pocket and the sc-PDB subpockets are represented as a cloud of 1.5 Å-spaced points annotated by eight pharmacophoric properties (hydrophobic, aromatic, H-bond acceptor, H-bond donor, H-bond acceptor and donor, positive ionizable, negative ionizable, null). [START_REF] Desaphy | Comparison and Druggability Prediction of Protein-Ligand Binding Sites from Pharmacophore-Annotated Cavity Shapes[END_REF] The term 'pocket' describes the full druggable cavity available at the surface of the protein while a subpocket is defined from its bound fragment. Since we aimed at targeting the ATP binding site in its type-I 'DMG in' conformation, the druggable pockets were first detected from 19 available CDK8 structures (Table S1). The largest pocket (830.3 Å 3 ) selected as representative was retrieved from the 5HBH 30 PDB entry (Figure 2). This pocket incorporates regions around the hinge, the gatekeeper F97, extends to a solvent exposed area near the αD helix whereas on an opposite side. It covers the DMG motif and reaches the αC-helix (Figure 2A). It thus spans several already described kinase subpockets: the adenine pocket, the front pockets FP-I and FP-II, the back pockets BP-I-A and BP-I-B in the gate area. [START_REF] Van Linden | Klifs: A Knowledge-Based Structural Database to Navigate Kinase-Ligand Interaction Space[END_REF] The 31 384 sc-PDB subpockets were compared and aligned to the CDK8 cavity with the in-house ProCare method (Figure S1). [START_REF] Eguida | A Computer Vision Approach to Align and Compare Protein Cavities: Application to Fragment-Based Drug Design[END_REF] Briefly, ProCare finds the best possible local alignment of cavity-defining points using a point cloud registration algorithm [START_REF] Rusu | 3d Is Here: Point Cloud Library (Pcl)[END_REF][START_REF] Zhou | Open3d: A Modern Library for 3d Data Processing[END_REF] and scores the alignment according to the overlap of pharmacophoric properties of the aligned points.

According to a preliminary study on the set of CDK8 structures, the original ProCare alignment c-FPFH fingerprint was modified to account only for the spatial distribution of pharmacophoric features (Figure S2-S3), a modification leading to a better alignment of CDK8 subpockets and fragments to the corresponding full cavities. Noteworthy, the novel c-FH fingerprint does not change the distribution of ProCare scores with respect to alignments generated by the original c-FPFH fingerprint. Once transformation matrices of the alignment of sc-PDB subpockets to the target cavity were obtained, the same rotation/translation matrices were applied to the corresponding sc-PDB fragments to position them in the CDK8 cavity. Posed fragments were then filtered according to five criteria (Figure 2B). Fragments originating from subpockets exhibiting a similarity score to the CDK8 pocket above a threshold value of 0.39 (previously shown to optimally discriminate known similar from known dissimilar binding sites) [START_REF] Eguida | A Computer Vision Approach to Align and Compare Protein Cavities: Application to Fragment-Based Drug Design[END_REF] were first selected, leading to a set of 12 661 fragments. Remaining fragments were further pruned according to three criteria: (i) belonging to a cofactor (therefore avoiding purine-base fragments), (ii) insufficient buriedness in the target cavity, (iii) no compliance to the fragment rule-of-three. [START_REF] Sterling | Zinc 15--Ligand Discovery for Everyone[END_REF] Remaining fragments were then annotated by one of the six CDK8 areas in which they were positioned: hinge (H), gate (GA1, GA2), solvent-accessible (SE1, SE2), αC helix (AC) (Table 1, Figure 3). 4 152 fragments could be unambiguously assigned to one CDK8 area: H (1.4%), GA1

(2.7%), GA2 (22.5%), SE1 (61.9%), SE2 (2.8%) and AC (8.7%) (Figure 3A). We next analyzed the origin of the sc-PDB ligands these fragments were derived from. As to be expected, 70% of fragments assigned to the hinge area (H) come from protein kinase inhibitors, the remaining 30% originating from a ligand co-crystallized with a protein that belong to a non-kinase family (Figure 3B). However, it should be noted that fragments from known CDK8 inhibitors were not selected as occupying the hinge region. Two simple reasons explain this absence: (i) the seven CDK8 ligands in the sc-PDB dataset are type II inhibitors binding to a DMG-out conformation and occupy the back pocket, (ii) the only CDK8 ligand (3RGF) that binds to the hinge could not be fragmented by our protocol and therefore did not pass our filters. The other areas (GA1, GA2, SE1, SE2, AC) were assigned fragments from both kinase (~25%) and non-kinase ligands (~75%). While the initial sc-PDB subpocket database contains 16% of entries from protein kinases, the enrichment observed for hinge-selected fragments (4.4) is logically due to the specific stereoelectronic features of the hinge area, notably the hydrogen bonding capacity of Asp98 and Ala100 backbone heteroatoms imposing complementary features on the ligand side. To limit the size of the library, all fragments were not considered for full enumeration of complete molecules. Whereas all fragments bound to H (n=57), GA1 (n=111) and SE2 (n=117) subpockets were selected, only 100 GA2-bound fragments were randomly chosen. Duplicates, in other words 2D identical fragments were kept as they do not originate from the same 3D subpocket, therefore resulted in different alignments that may differently impact molecules design.

Comprehensive statistics of the pairwise fragment similarity (Figure S4) and the observed distribution of their physicochemical properties (Figure S5) clearly evidence their chemical diversity. 385 fragments were selected at this stage for the next linking stage.

Round-1 library generation. The DeLinker deep generative model [START_REF] Imrie | Deep Generative Models for 3d Linker Design[END_REF] was used to link the aboveselected fragments. Briefly, DeLinker uses a graph-based deep generative model, trained on the ZINC [START_REF] Sterling | Zinc 15--Ligand Discovery for Everyone[END_REF] or PDBbind [START_REF] Wang | The Pdbbind Database: Collection of Binding Affinities for Protein-Ligand Complexes with Known Three-Dimensional Structures[END_REF] databases, to expand bond by bond the two fragments to be connected until final SMILES strings are generated by a variational autoencoder while keeping 3D constraints through a set of distances and angles between connectable atoms. [START_REF] Imrie | Deep Generative Models for 3d Linker Design[END_REF] In the current work, all possible connectable atoms of hinge-annotated fragments (H) were used as seeds to find potential connectable atoms in fragments filling three remaining subpockets (GA1, GA2, SE2) (Figure S6).

An atom is considered connectable if it is a heavy atom covalently bonded to a hydrogen, that bond being used as exit vector for the linking. Pairs of atoms belonging to different fragments are then associated by restricting the angle between the exit vectors and distances between the corresponding heavy atoms (see Methods) in order to avoid pointless connections and lower the number of combinations (Figure S7). Starting from 385 fragments, 1 517 488 SMILES strings were generated by linking fragment pairs with DeLinker. 15% of the proposed solutions were discarded since they correspond to uncomplete molecules where the SMILES consisted of a linker moiety attached to only one of the two fragments (Figure 4). (2) Druglikeness is defined by customized OpenEye Filter rules available in Table S2.

(3) Synthetic accessibility score. [START_REF] Ertl | Estimation of Synthetic Accessibility Score of Drug-Like Molecules Based on Molecular Complexity and Fragment Contributions[END_REF] (4) Filter to remove unwanted aliphatic linkers.

The remaining molecules were filtered for drug-likeness (Table S2) resulting in 566 989 unique SMILES.

Although the redundant SMILES per pair of connectable atoms were removed during the linking process, duplicated molecules still arose when connecting the same 3D fragments via equivalent exit atoms (symmetry cases) or connecting the same duplicated fragments originating from different subpockets. After keeping only molecules that are likely to be synthesized (SAscore [START_REF] Ertl | Estimation of Synthetic Accessibility Score of Drug-Like Molecules Based on Molecular Complexity and Fragment Contributions[END_REF] ≤ 3), only those having a linker compliant with defined rules (Figure S8) were finally kept. The remaining 141 125 molecules composed the first-round R1 library (Figure 4). A majority of the generated molecules arose from combining the hinge and the solvent-exposed SE2 fragments which account for more than 50% of the sets (Figure 5). Indeed, the average number of generated SMILES strings per pair of H-SE fragments is higher than for the two other areas, a consequence of having more pairs of connectable atoms and more generated linkers per connectable atoms for the H-SE subpockets. While it was expected that kinase-derived fragments would contribute to most of the generated molecules, only 14% of SMILES strings were generated by linking two kinase-bound fragments. Interestingly, around 26% of the molecules were made of two fragments originating from a non-kinase protein. Interestingly, the observed proportions do not vary between the full set, the drug-like subset and the R1 set (Figure 5). Most of the generated molecules (> 90 %) were already compliant with the Lipinski's rule of five [START_REF] Gaulton | Chembl: A Large-Scale Bioactivity Database for Drug Discovery[END_REF] (Figure S9). Albeit two fragments were assembled, many generated molecules still remained in the fragment space with around 10 % of SMILES strings being compliant with the fragment rule-of-three [START_REF] Congreve | A 'Rule of Three' for Fragment-Based Lead Discovery?[END_REF] (Figure S9). Filtering the designed molecules to R1 library members did not bias our selection towards molecules with particular properties as the distribution of the molecular properties, although reported individually, remained comparable among the sets (full, drug-like and R1; Figure S9). To give insights on the chemical space covered by R1 library members, we further assessed its overlap with either a broad purpose bioactive chemical space [START_REF] Gaulton | Chembl: A Large-Scale Bioactivity Database for Drug Discovery[END_REF] (1.7 million ChEMBL compounds) or a recently described kinasefocused ligand space (6.7 million KinFragLib library members). [START_REF] Sydow | Kinfraglib: Exploring the Kinase Inhibitor Space Using Subpocket-Focused Fragmentation and Recombination[END_REF] 259 unique R1 library molecules were exactly found in ChEMBL among which only a few have been assayed against protein kinases, while only five R1 library compounds were identical to KinFragLib molecules. Considering similarity, only 0.85% and 13% of R1 library members were found similar to KingFragLib and ChEMBL molecules, respectively, according to a Tanimoto coefficient, computed from Morgan2 fingerprints 43 higher than 0.60. The herein proposed computational workflow is therefore able to generate really new chemical entities, the chemical diversity of the generated molecules stemming from the diversity of the seed fragments pool, the connectivity and the possible linkers.

As a first validation of the structure-based workflow, we verified whether the drug-like subset contains molecules highly similar to 302 submicromolar human CDK8 inhibitors retrieved from the ChEMBL database. Using the similarity search protocol described in the methods section, we found 44 molecules that matched with 35 unique known CDK8 inhibitors (representing three series of congeneric molecules). While these molecules were built with fragments from all possible areas, most of them were assembled from hinge-fragments originally co-crystallized with protein kinases, linked to fragments originally co-crystallized with non-kinase proteins.

The round-1 library contains novel and potent CDK8 inhibitors. To identify chemically novel hits, we filtered first-round R1 library members by dissimilarity (Tanimoto coefficient < 0.5, RDKit7 fingerprints) to all CDK8 compounds available in ChEMBL [START_REF] Gaulton | Chembl: A Large-Scale Bioactivity Database for Drug Discovery[END_REF] and to all seed sc-PDB fragments. With respect to the previously used Morgan 2 fingerprint that is best suited to ligand-based virtual screening, the RDKit7 descriptor [START_REF] Riniker | Open-Source Platform to Benchmark Fingerprints for Ligand-Based Virtual Screening[END_REF][START_REF] O'boyle | Comparing Structural Fingerprints Using a Literature-Based Similarity Benchmark[END_REF] was here chosen for its ability to account for substructure similarities between library members and known inhibitors. [START_REF] Riniker | Open-Source Platform to Benchmark Fingerprints for Ligand-Based Virtual Screening[END_REF][START_REF] O'boyle | Comparing Structural Fingerprints Using a Literature-Based Similarity Benchmark[END_REF] Hits were then searched for availability among 8.2 million commercially available drug-like compounds (Table S3) to select 37 compounds that are identical or very similar (Tanimoto coefficient > 0.90, RDKit7 fingerprints) to their queries (Table S4). These compounds were purchased and tested for CDK8 inhibited the CDK8 kinase by more than 50% at the single concentration of 10 μM (Figure 6). Notably two related compounds (12 and 37), exhibiting more than 80% inhibition were assembled from the same pair of 3D fragments by just inverting the ester linkage (Figure 6). They differ from the original R1 library members by just a carbon atom (methoxy for ethoxy substitution, Table S4). Round-2 library design by fragment hit growing. The most potent hit (12) from round-1 library , generated by linking a H-area pyridine fragment to a GA2-area p.methoxyphenyl fragment, is still a fragment-like compound (MW = 229 g.mol -1 ) that can be optimized by growing towards the nearby and yet unexploited SE2 and GA1 subpockets. We thus explored the possible connections between the hinge-binding fragment of 12 and all remaining SE2 or GA1-anchored fragments, to generate a secondround library R2 of 5 700 compounds. R2 library members were filtered by physicochemical properties (number of rotatable bonds ≤ 6, no chiral centers) and synthetic accessibility (SAscore ≤ 3) to yield a final set of 151 candidates (Table S5). Six representative compounds (Table 2) were chosen for their ease of synthesis (i.e. availability of building blocks, costs of goods, number of synthetic steps) and predicted buriedness upon preliminary docking to CDK8. Three linkers (urea, piperidine, pyrazole)

were chosen for their capacity to connect the H-anchoring pyridine ring to a SE2-anchored phenyl fragment. Two positions of the pyridine ring (ortho and meta position to the benzoyl ester) were predicted compatible, therefore leading to six possible analogs (Table 2).

Table 2. Round-2 library of optimized hits and their CDK8 inhibitory potency. The six compounds were synthesized (Scheme S1), checked for purity (Figures S10-S15) and tested for in vitro CDK8 inhibition using the same HTRF assay as described above, to build concentration-response curves (Figure 7). The ProCare pose of the three fragment-bound subpockets used to yield compound 49 (Figure S16)

Compound

illustrate the alignment of key residues within each subpocket to the target cavity, and the capacity of ProCare to detect local similarities extending to 5-7 conserved amino acid pairs. As part as the fragment selection workflow, only fragments bound to subpockets whose pharmacophoric points share enough polar matches with that of the target query were considered for linking. In the present cases, several polar matches (h-bond donor/acceptor, positive/negative ionizable, aromatic) could be detected between each of the three subpockets and the CDK8 cavity (Figure S16). The topological relationships (distance between connectable atoms, angles between exit vectors) between all selected fragments, notably that leading to the final hits 12 49 are given in Figure S17.

TheIts putative binding mode of compound 49, deduced form molecular docking, suggests that the pyridine nitrogen atom h-bonds to the hinge backbone atoms (E98, A100) while the ethoxyphenyl and the newly introduced pyrazole moieties exhibit π-π interactions to H106 (SE2 subpocket) and the gatekeeper F97 (GA1 subpocket). Last, the terminal phenyl ring is oriented towards K52 (GA2 subpocket) for a putative π-cation interaction (Figure 8). While the parent hit 12 showed two possible docking poses (ethoxyphenyl towards GA2 or SE2), growing by a pyrazole prioritized the SE2 orientation, still with exhibited interactions compatible with the rationale of the initial fragment alignments. At this point, we should recall that neither early safety (e.g. kinase selectivity) nor pharmacokinetic properties (e.g. metabolic stability) have been considered in either generating or post-processing the target-focused library members. Although technically feasible, target selectivity assessment requires applying the same workflow to different cavities and prioritizing compounds generated only for the target of interest. This approach is feasible for a comparing a few targets but is rapidly impracticable at a larger scale (e.g. whole kinome). It has not been applied in the current study aimed at demonstrating the proof-of-concept of the structure-based workflow.

CONCLUSION

We herewith propose a novel fragment-based library design method to generate target-focused compound libraries. The originality of the approach is that seed fragments are chosen from a large repertoire of protein-bound fragmented ligand X-ray structures, and positioned in the target according to the local similarity of their protein subpocket to the target cavity. This ligand-agnostic posing protocol does not require scoring protein-ligand interactions and is fuzzy enough to transfer ligand information across unrelated target spaces. Once fragments have been posed, they are linked by a deep generative model to enumerate full molecules which are later post-processed to account for druglikeness and synthetic accessibility. The linking step still deserves improvement, notably to enumerate candidate molecules directly in the original target 3D coordinate frame. Hence, the variational autoencoder used here generates SMILES strings and just accounts for the target binding site topology in the form of topological relationships between fragment atoms to be connected. A true 3D deep generative model [START_REF] Li | Structure-Based De Novo Drug Design Using 3d Deep Generative Models[END_REF] considering complementarity to the binding site shape and the ligand conformational freedom would be highly desirable to link subpocket-selected seed fragments. It would avoid a tedious post-processing of unrealistic solutions and the necessary docking of candidates to verify whether the starting binding hypothesis of the seed fragments is conserved.

When applied to the test case of the CDK8 kinase, the method was able to quickly suggest potential inhibitors. Within two iterations and 43 compounds, a single digit nanomolar inhibitor could be identified thereby demonstrating a first proof-of-concept of the underlying methodology. Interestingly, the POEM method is technically applicable to any target of known 3D structure and does not require prior ligand knowledge.

MATERIAL AND METHODSEXPERIMENTAL SECTION

CDK8 cavity detection. All publicly available X-ray structures of human CDK8 (UniProt accession number P49336; Table S1) were downloaded from the Protein Data Bank [START_REF]The Protein Data Bank[END_REF][START_REF] Berman | The Protein Data Bank[END_REF] . Type I structures (DMGin, α-C helix-out) were put in the same coordinates frame by subsequent structural alignment to the 4F7S reference with Maestro v.2019-3 (Schrödinger, New York, NY 10036, U.S.A.) and refinement to ensure that the hinge residue Ala100 heavy atoms were fitted. Aligned structures (proteins, co-factors, ligands) were then protonated with Protoss v.4.0, [START_REF] Bietz | Protoss: A Holistic Approach to Predict Tautomers and Protonation States in Protein-Ligand Complexes[END_REF] while optimizing the intra and inter-molecular hydrogen bond network. After discarding crystallization additives, each PDB entry was split to afford a protein (no water molecules) and a ligand in separate mol2 files using SYBYL-X 2.1.1 (Certara USA, Inc., Princeton, NJ 08540, U.S.A.). For each protein file, entire cavities ("CAVITY_ALL" output) were next computed with the VolSite 33 module of the IChem v.5.2.9 package, [START_REF] Da Silva | Ichem: A Versatile Toolkit for Detecting, Comparing, and Predicting Protein-Ligand Interactions[END_REF] using default parameters and saved as point clouds annotated by pharmacophoric features. Only cavities corresponding to the catalytic site were retained for the next steps. Upon visual inspection, the corresponding three clouds for PDB entry 5HBH were merged into a single cavity in mol2 file, yielding the reference pocket for CDK8.

sc-PDB subpocket-fragment database. 16 034 drug-like ligands in their protein-bound X-ray structure were retrieved from the sc-PDB database [START_REF] Desaphy | Sc-Pdb: A 3d-Database of Ligandable Binding Sites--10 Years On[END_REF] of druggable protein-ligand complexes and fragmented in three dimensional (3D) space within their protein binding site using the IChem fragmentation tool. [START_REF] Desaphy | Sc-Pdb-Frag: A Database of Protein-Ligand Interaction Patterns for Bioisosteric Replacements[END_REF] Only fragments exhibiting at least 4 non-covalent interactions [START_REF] Marcou | Optimizing Fragment and Scaffold Docking by Use of Molecular Interaction Fingerprints[END_REF] (out of which one is polar, hydrogenbond or electrostatic interaction) with the protein target were retained. The fragments exit bonds (dummy atoms 'Z') were converted into hydrogen atoms. The immediate protein environment of each of the c-FPFH descriptor. [START_REF] Eguida | A Computer Vision Approach to Align and Compare Protein Cavities: Application to Fragment-Based Drug Design[END_REF] For each subpocket-cavity pair, the optimal alignment matrix was used to position the corresponding sc-PDB fragment into the CDK8 cavity. The comparison protocol was validated by successful cross-comparison of CDK8 subpockets from type I PDB entries (Figure S3).

In the second stage, aligned sc-PDB fragments were filtered according to their subpocket similarity to the CDK8 cavity (ProCare score ≥ 0.39), their compliance to the fragment rule-of-three, [START_REF] Congreve | A 'Rule of Three' for Fragment-Based Lead Discovery?[END_REF] and their embedding into the CDK8 cavity such that at least half of the fragment atoms are less than 1.5 Å away to the closest CDK8 cavity point. Fragments originating from the sc-PDB list of cofactors were excluded.

Resulting fragments were further annotated with the CDK8 cavity area to which they have been aligned

based on their distance (closest heavy atom should be within 6 Å) to subpocket-specific preliminary defined atom centers (hinge H area, Asp98 O atom and Ala100 N and O atoms; gate area 1 GA1, Phe97

CA atom; gate area 2 GA2, Lys52 NZ atom; solvent-exposed area 1 SE1, Arg356 CZ atom; solventexposed area 2 SE2 subpocket, His106 CE1 atom; αC area AC, Ser62 CA atom). For selecting hingebinding fragments, hydrogen bonds to Asp98 O or Ala100 N or O was mandatory. Since a few fragments were assigned to multiple subpockets, the following prioritization scheme was applied: H annotation takes precedence over all the other annotations, therefore a fragment interacting with the hinge centers is only annotated as such. SE1 and SE2 were defined compatible so that fragments annotated as from both areas were automatically assigned only SE2. Similarly, fragments annotated as from both AC and GA2 areas were automatically assigned only GA2. In any other case of combination (e.g. fragments annotated as from GA2 and SE1), the annotations were considered ambiguous and the fragments were discarded.

In the third stage, H fragments were defined connectable to either GA1, GA2 or SE2 fragments (in the current work, although other connections are possible). Selected fragments were converted into sdf format with OpenEye v.2.5.1.4. toolkit. 51 For each pair of fragments with hydrogen atoms connected, pairs of connectable atoms were searched based on their respective orientation as follows. A right circular cone (half-angle=π/4) is projected along the bond axis between any heavy atom A i and a bound hydrogen atom H i . A connectable atom pair A 1 A 2 is selected if heavy atoms A 1 and A 2 are located in the projection cone of their counterpart (Figure S7).

In the fourth stage, the recently-described DeLinker [START_REF] Imrie | Deep Generative Models for 3d Linker Design[END_REF] deep learning method was employed to generate linkers between above-described connectable atom pairs using the default model distributed with the package and a batch size of 1. Input data were prepared as ZINC atom types features to be ready for DeLinker using the 'prepare_data' module and by setting the 'test' parameter of the 'preprocess' function to 'True' as molecules are to be found. The linker length was set to a minimum of on and a maximum of six heavy atoms. Other parameters were kept by default. Generated molecules were saved as SMILES strings and further processed to remove redundancy for each connectable atom pair.

In the final stage, unsuccessful linking attempts where only a single fragment is attached to the linker were removed using the function 'get_linker' in the 'frag_utils' utility. The remaining SMILES were filtered to keep only drug-like compounds according to in-house rules (Table S2). Next, the synthetic accessibility scores were computed with the the SAscore 40 method distributed with RDKit 52 to remove molecules with SAscore higher that three. Finally, molecules made of long flexible linkers were discarded according to our in-house filtering workflow (Figure S8), resulting in the first-round library (R1).

Comparison with ChEMBL and KinFraglib ligands. Standardized ChEMBL (1.7 million compounds) and KinFragLib (6.7 million) data were retrieved from the KinFragLib website. [START_REF] Kinfraglib | [END_REF] Pairwise 2D fingerprint similarity to R1 molecules were assessed with RDKit 52 Morgan (radius = 2) topological fingerprint (default parameters, maximum path = 7).

Comparison to known CDK8 inhibitors.

A search in the ChEMBL database [START_REF] Chembl | [END_REF][START_REF] Gaulton | Chembl: A Large-Scale Bioactivity Database for Drug Discovery[END_REF] for human CDK8 target assays resulted in three target report cards (CHEMBL3038474, CHEMBL5719 and CHEMBL3885556) from which bioassay data were joined and processed to keep compounds with a half maximal inhibitory concentration IC 50 inferior or equal to 1 μM. Duplicates were then removed according to and the SMILES were standardized with OpenEye Filter v.3.0.1.2 (OpenEye Scientific Software, Santa Fe, NM 87508, U.S.A.). The final list of 302 inhibitors was searched in the generated drug-like subset described above for substructure 2D similarity using both RDKit Morgan (radius = 2) and topological (maximum path = 7) fingerprints and a combination of Tanimoto (Tc) and Tversky (Tv) metrics. Pairs were reported when morgan2 Tc ≥ 0.6 or morgan2 Tv ≥ 0.8 or RDKit7 Tc ≥ 0.75 or RDKit7 Tv ≥ 0.9.

Search for new potential CDK8 inhibitors. R1 library members were considered as potentially new at the condition that their similarity to any of 946 unique human CDK8-tested compounds (both active and inactive) reported in ChEMBL (target card reports CHEMBL3038474, CHEMBL5719 and CHEMBL3885556) and any of the 31 384 sc-PDB fragment is inferior to 0.50 (Tanimoto coefficient from RDKit topological fingerprints). Last, the subsequent list was searched for substructure similarity (RDKit topological fingerprint Tanimoto ≥ 0.90) to an in-house library of 8 280 193 commercially available drug-like compounds (Supporting Table S3).

Molecular docking. Virtual hits were drawn as 2D sketches with ChemAxon MarvinSketch v.16.10.17, tested compounds used in the biological assays were determined by HPLC/MS using the area percentage method on the UV trace recorded at a wavelength of 254 nm. All compounds were found to have >95% purity.

(
1-(3-hydroxypyridin-4-yl)3-phenylurea (38). To a stirred solution of phenylisocyanate (0.4 g, 3.4 mmol) in DMF (5 ml) was added a solution of 4-aminopyridin-3-ol hydrochloride (0.5 g, 3.4 mmol) in DMF (5 ml) followed by the addition of triethylamine (1.4 ml, 10.2 mmol) at room temperature (r.t.). The resulting mixture was stirred at room temperature overnight. The reaction mixture was concentrated under reduced pressure and the crude residue was purified by HPLC to afford 50 mg (6%) of the 1-(3hydroxypyridin-4-yl)-3-phenylurea 38 as a white solid which was used for the next step without further purification. 1-(5-hydroxypyridin-3-yl)-3-phenylurea (40). Compound 40 was prepared as described above for compound 38, starting from 5-aminopyridin-3-ol hydrobromide (0.65 g, 3.4 mmol). The reaction mixture was concentrated under reduced pressure and the crude residue was purified by HPLC (method B) to afford 60 mg (8%) of 1-(3-hydroxypyridin-5-yl)-3-phenylurea 40 as a white solid which was used for the next step without further purification. (41). Compound 41 was prepared as described above for compound 39, starting from 1-(5-hydroxypyridin-3-yl)-3-phenylurea 40 (60 mg, 0.264 mmol). The residue was purified by HPLC (method B) to afford compound 41 (36 mg, 45%) as a white solid. [START_REF] Erlanson | Twenty Years On: The Impact of Fragments on Drug Discovery[END_REF] 4-(1-phenyl-3,6-dihydro-2H-pyridin-4-yl)pyridin-3-ol (42). To a stirred solution of 4-iodopyridin-3-ol (0.63 g, 2.86 mmol, 1.1 eq.) and 1-phenyl-4-(4,4,5,5-tetramethyl-1,3,2-dioxaborolan-2-yl)-3,6-dihydro-2H-pyridine (0.74 g, 2.6 mmol, 1 eq.) in a mixture of 1,4-dioxane and water (20 ml, v/v=4:1), K 2 CO 3 (1.8 g, 13 mmol, 5 eq.) was added and purged with argon for 30 min followed by the addition of Pd(dppf)Cl 2 (0.1 g, 0.05 eq.) and stirred at 90°C overnight. After completion, the reaction mixture was cooled to room temperature, diluted with ethyl acetate and water. The organic layer was washed with water and brine, dried over anhydrous sodium sulfate and evaporated under reduced pressure. The crude product was purified by column chromatography on silica gel (hexane/EtOAc) to afford 42 (251 mg, 38%). (43). Compound 42 (251 mg, 1 mmol) was dissolved in MeOH (20 ml), followed by addition of Pd (10 wt % on activated carbon, 50 mg), and then the resulting suspension was stirred at room temperature under 1 atm. hydrogen pressure overnight. The resulting reaction was filtered, concentrated under reduced pressure, and dried under vacuum, to afford 43 (201 mg, 79%) which was used for the next step without further purification.

4-(3-phenylureido)pyridin-3-yl 4-ethoxybenzoate (39

5-(3-phenylureido)pyridin-3-yl 4-ethoxybenzoate

4-(1-phenyl-4-piperidyl)pyridin-3-ol

[4-(1-phenyl-4-piperidyl)-3-pyridyl] 4-ethoxybenzoate (44). A solution of compound 43 (201 mg, 1 eq.),

4-ethoxybenzoic acid (131 mg, 1 eq.), Et 3 N (0.27 ml, 2.5 eq.) and HATU (360 mg, 1.2 eq.) in dry DMSO

5-(1-phenyl-3,6-dihydro-2H-pyridin-4-yl)pyridin-3-ol (45).

To a stirred solution of 5-iodopyridin-3-ol (0.63 g, 2.86 mmol, 1.1 eq.) and 1-phenyl-4-(4,4,5,5-tetramethyl-1,3,2-dioxaborolan-2-yl)-3,6-dihydro-2H-pyridine (0.74 g, 2.6 mmol, 1 eq.) in a mixture of 1,4-dioxane and water (20 ml, v/v=4:1), K 2 CO 3 (1.8 g, 13 mmol, 5 eq.) was added and purged with argon for 30 min followed by the addition of Pd(dppf)Cl 2 (0.1 g, 0.05 eq.) and stirred at 90 °C overnight. After completion, the reaction mixture was cooled to room temperature, diluted with ethyl acetate and water. The organic layer was washed with water and brine, dried over anhydrous sodium sulfate and evaporated under reduced pressure. The crude product was purified by column chromatography on silica gel (hexane/EtOAc) to afford compound 45 (326 mg, 49%). (46). Compound 45 (251 mg, 1 mmol) was dissolved in MeOH (20 ml), followed by addition of Pd (10 wt% on activated carbon, 50 mg), and then the resulting suspension was stirred at room temperature under 1 atm. hydrogen pressure overnight. The resulting reaction was filtered, concentrated under reduced pressure, and dried under vacuum, to afford compound 46 (220 mg, 86%) which was used for the next step without further purification.

4-(1-phenyl-4-piperidyl)pyridin-3-ol

[5-(1-phenyl-4-piperidyl)-3-pyridyl] 4-ethoxybenzoate (47). A solution of compound 46 (200 mg, 1 eq.),

4-ethoxybenzoic acid (131 mg, 1 eq.), Et 3 N (0.27 mL, 2.5 eq.) and HATU (360 mg, 1.2 eq.) in dry DMSO ( 4-bromopyridin-3-yl 4-ethoxybenzoate (48). A solution of 4-bromopyridin-3-ol (300 mg, 1.7 mmol, 1 eq.), 4-ethoxybenzoic acid (310 mg, 1.87 mmol, 1.1 eq.), DIPEA (0.89 ml, 5.1 mmol, 3 eq.) and HATU (760 mg, 2 mmol, 1.2 eq.) in DMF (10 ml) was stirred at 25°C for 16 h. The reaction mixture was poured into 50 ml of water and extracted with ethyl acetate (3x15 ml). The combined organic layers were washed with saturated ammonium chloride solution (50 ml) and brine (50 ml), dried over anhydrous sodium sulfate, and concentrated under reduced pressure to afford compound 48 as a brown solid (320 mg, purity 85%), which was used in the next step without further purification.

4-(1-phenyl-1H

-pyrazol-4-yl)pyridin-3-yl 4-ethoxybenzoate (49). A mixture of compound 48 (200 mg, 0.62 mmol, 1 eq.), 1-(phenylpyrazol-4-yl)boronic acid (130 mg, 0.68 mmol, 1.1 eq.), cesium carbonate (400 mg, 1.24 mmol, 2 eq.) and Pd(dppf)Cl 2 (25 mg, 0.03 mmol, 0.05 eq.) in dioxane/water (5 ml, 10:1 v/v) was degassed and stirred at 105°C for 16 h under inert atmosphere. After cooling, the reaction acceptor fluorophore on the kinase tracer. Binding of an inhibitor to CDK8/CyclinC competes for binding with the tracer, resulting in a loss of FRET. Binding assay was performed into 384-well small volume plates (CORNING 3824) using kinase buffer provided by supplier (HEPES 50mM pH7.5, MgCl2 10mM, EGTA 1mM, Brij-35 0.01%) in a final volume of 15 µL. Briefly, 5µL of 3X compound (increasing concentrations from 3.10 -11 to 3.10 -5 M) prepared in kinase buffer are added to 5µL of 3X kinase/Ab solution (15nM kinase, 6nM biotin anti-His-tag antibody, 6nM Eu-streptavidin) and 5µL of 30nM kinase tracer236 (Kd 8 nM). The plate was incubated 1h at room temperature before reading with a TRFcompatible multi-well plate reader (Envision, PerkinElmer) using a classic TRF reading protocol (excitation at 337 nm; donor emission measured at 620 nm; acceptor emission measured at 665 nm).

The TR-FRET signal was collected both at 665 and 620 nm, and TR-FRET ratios were calculated (acceptor signal value divided by donor signal value). IC 50 and K i values of the tested compounds were determined from competitive binding curves using GraphPad Prism software (version 6.07) as follows:

𝑆 𝑆 𝑆 𝑆 1

REFERENCES Table of content graphic

Target-focused library design by pocket-applied computer vision and fragment deep generative linking.

Merveille Eguida, Christel Schmitt-Valencia, Marcel Hibert, Pascal Villa, and Didier Rognan.

Figure 1 .

 1 Figure 1. Overall workflow of the POEM computational method including in vitro experimental validation.

Figure 2 .

 2 Figure 2. Seed fragments selection to fill the CDK8 query cavity. A) Description of the reference CDK8 pocket (PDB ID: 5HBH). Cavity points (grey dots, 246 points) delineate a ligand-accessible envelope (solid surface, 830.3 Å 3 ) and areas (hinge, H; gate area 1, GA1; gate area 2, GA2; solvent-exposed area 1, SE1; solvent-exposed area 2, SE2; αc area, AC) according to the distance to key CDK8 atoms (spheres). B) Fragments selection workflow. (1) A list of cofactors (PDB HET code) is provided in the sc-PDB database. (2) Fragments buriedness is approximated as the percentage of heavy atoms within 1.5 Å of one CDK8 cavity point. (3) fragment rule-of-three: 37 molecular weight ≤ 300 g.mol-1, logP ≤ 3, H-bond donor count ≤ 3 and H-bond acceptor count ≤ 3. (4) ambiguous annotation denotes assignment of two or more incompatible areas (Methods section) out of the six possible areas. (5) All annotated fragments from H, GA1, SE2 areas and a random sampling of 100 fragments from GA2 were selected.

Figure 4 .

 4 Figure 4. Focused library design via linking selected fragments. Fragments aligned in the H area were paired with fragments from GA1, GA2 and SE2 areas. SMILES were generated by linking fragment pairs with DeLinker 23 and filtered to compose the first-round library R1. (1) Successful linking signifies that both fragments have been attached to the linker whereas cases where only one of the fragments was linked were considered unsuccessful. (2) Druglikeness is defined by customized OpenEye Filter rules available in Table S2. (3) Synthetic accessibility score. 40 (4) Filter to remove unwanted aliphatic linkers.

Figure 5 .

 5 Figure 5. Protein origin of fragments pairs in newly generated molecules. From left to right, the full set after cleaning unsuccessful generation out (GEN), the drug-like subset (DL) and round-1 library (R1). The distribution is given for the combinations annotated by the targeted CDK8 area (H, hinge; GA1, gate area 1; GA2, gate area 2, SE2, solvent-exposed area 2) and color-coded according to the protein origin (co-crystallized target) of the two connected fragments (K, protein kinase; O, other; K-K, both fragments were derived from a protein kinase structure; K-O, H-fragment derived from a protein kinase and the other fragment from a non-kinase protein structure; O-K, H-fragment derived from a non-protein kinase and the other fragment from a kinase protein structure; O-O, both fragments were derived from a non-kinase protein structure).

  inhibition in a homogeneous time-resolved fluorescence (HTRF) assay aimed at measuring the FRET signal between a fluorescent-labelled ATP competitive inhibitor and the fluorescent-tagged CDK8 soluble kinase (see Methods). Six out of the 37 tested molecules (compounds 9, 11, 12, 29, 32, 37)

Figure 6 .

 6 Figure 6. CDK8 inhibition (LanthaScreen Eu kinase competitive binding assay) by 37 commercially available compounds identical or very similar to R1 library members. Results are expressed as mean ± SEM of two independent experiments using a 10 μM concentration of competitor (STA, staurosporine control).

  moiety (blue) is attached via different linkers (red) to round-1 compound 12. b Inhibition of CDK8 measured in a LanthaScreen Eu kinase competitive binding assay. Results are expressed as mean ± SEM of three independent experiments. c confidence interval at a 95% confidence level

Figure 7 .

 7 Figure 7. Inhibition of human CDK8 by six selected round-2 library compounds. Concentration-response curves are derived from three independent experiments with duplicates per experiment.

Figure 8 .

 8 Figure 8. PLANTS docking pose of compound 49 (green sticks) to the catalytic site of CDK8 (PDB ID 5HBH, solid surface). H-bond to the hinge (E98, A100) and π-π interactions to F97, H106 are displayed by yellow broken bonds.

  selected fragment was considered to compute VolSite point clouds, keeping only those with at least 3 points, each being closer than 4.0 Å from any fragment heavy atom ("CAVITY_4" output), thereby defining a subpocket point cloud in mol2 file format for 31 384 fragments. CDK8-focused library design. In the first stage, 31 384 sc-PDB subpocket point clouds (Figure S1) were aligned to the reference 5HBH CDK8 cavity point clouds with ProCare 28 v.0.1.1 using default parameters and the c-FH color-based descriptor (Figure S2) corresponding to the eight terminal bins

( 2

 2 ml) was stirred at room temperature for 12h. The completion of the reaction was monitored by LCMS. The mixture was purified by HPLC (Method A) to give compound 44 (120 mg, 38% yield) as a white solid. 1H NMR (400 MHz, DMSO-d6). δ 8.46 (d, J = 5.4 Hz, 2H), 8.15 -8.09 (m, 2H), 7.50 (d, J = 5.1 Hz, 1H), 7.22 -7.10 (m, 4H), 6.93 (d, J = 8.2 Hz, 2H), 6.75 (t, J = 7.3 Hz, 1H), 4.16 (q, J = 6.9 Hz, 2H), 3.78 (d, J = 12.3 Hz, 2H), 2.87 -2.79 (m, 1H), 2.63 (t, J = 10.0 Hz, 2H), 1.82 (t, J = 5.1 Hz, 4H), 1.37 (t, J = 7.0 Hz, 3H). LC-MS (ESI) m/z 403.2 [(M+H)+, calcd. C 25 H 27 N 2 0 3 , 403.2].

10 S

 10 ABBREVIATIONS

Table 1 .

 1 Annotation of the CDK8 target cavity by key pharmacophoric atoms.

	Area	Label Key CDK8 atoms		KLIFS subpockets a
	Hinge area	H	Asp98.O, Ala100.N, Ala100.O AP
	Gate area 1	GA1 Phe97.CA	(gatekeeper	AP, BP-I-A, BP-I-B
			residue)
	Gate area 2	GA2 Lys52.NZ		AP, FP-I, FP-II
	Solvent-accessible area 1	SE1	Arg366.CZ		-
	Solvent-accessible area 2	SE2	His106.CE1		-
	αC helix area	AC	Ser62.CA		-

a Full

or partial overlap with KLIFS 34 subpockets: AP: adenine pocket, BP: back pocket, FP: front pocket Figure 3. CDK8

  subpocket occupancy of sc-PDB fragments. A) Assignment of CDK8 pocket areas to 4 152 sc-PDB fragments. B) Proportion of sc-PDB fragments per area.

  ChemAxon Ltd., 1031 Budapest, Hungary) saved in sdf file format, ionized at physiological pH with The search space was set at 13 Å from the binding site center with a search speed of 1 (highest accuracy). 10 poses were generated per ligand, scored by the ChemPLP scoring function and clustered using a root-mean square deviations (RMSD) of 2 Å on ligand heavy atoms. The flipped/rotated side chains were reconstructed in the protein structure for each corresponding PLANTS pose when applicable. Molecular descriptors (molecular weight (g.mol -1 ), the count of heavy atoms (non-hydrogen atoms), logP, polar surface area (Å), count of H-bond acceptor, count of H-bond donor, count of rotatable bonds, count of ring systems, count of heteroatoms, bonds) were computed with RDKit. Data were processed with Python v.3.7.

	Method B) 2-10 min 0-50% acetonitrile, 30 ml/min ((loading pump 4 ml acetonitrile); column: SunFire
	C18; 100 mm x 19 mm; 5 µm)
	Analytical RP-HPLC-MS was performed using Agilent Technologies 1260 Infinity LC/MSD system with
	DAD\ELSD Alltech 3300 and Agilent LC\MSD G6120B mass-spectrometer using the following
	OpenEye Filter v.2.5.1.4 and finally converted in 3D structures (mol2 file) with Corina v.3.40 (Molecular Networks GmbH, 90411 Nürnberg, Germany), generating all possible stereoisomers and ring conformers simultaneously. The prepared molecules were docked into the above-described CDK8 Chemistry. All reactions were carried out under usual atmosphere unless otherwise stated. Chemicals and solvents were purchased from Enamine (LV-1035 Riga, Latvia) and were used without further purification. Yields refer to isolated compounds, estimated to be >95% pure as determined by 1H NMR or HPLC. 1H NMR spectra were recorded at 298 K on Bruker Avance III Spectrometer operating at 400 MHz. All chemical shift values δ and coupling constants J are quoted in ppm and in Hz, respectively, multiplicity (s = singulet, d = doublet, t = triplet, q = quartet, quin = quintet, sex = sextet m = multiplet, br = broad). Preparative HPLC was performed using two methods: Method A) 2-10 min 30-70% acetonitrile, 30 ml/min ((loading pump 4 ml acetonitrile); column: YMC-ACTUS TRIART (C18; 100 mm x 20 mm; 5 µm); acquisition parameters: column, Agilent Poroshell 120 SB-C18 4.6x30mm 2.7 μm with UHPLC Guard Infinity Lab Poroshell 120 SB-C18 4.6x 5mm 2.7 μm; Temperature 60 C; Mobile phase А -acetonitrile : water (99:1%), 0.1% formic acid, В -water (0.1% formic acid); Flow rate 3 ml/min; Gradient : 0.01 min -99% B, 1.5 min -0% B, 1.73 min -0% B, 1.74 min -99% B; Injection volume 0.5μl; Ionization mode cavity using PLANTS 55 v.1.2 Molecular data analysis. Plots were generated with Matplotlib v3.0.2 56 in Python v.3.7. Electrospray ionization (ESI); Scan range m/z 83-600; DAD 215 nm, 254nm, 280 nm. Purities of all

Data visualization. Molecules were drawn in 2D with RDKit and MarvinSketch v.16.10.17, (ChemAxon Ltd., 1031 Budapest, Hungary). Three-dimensional structures were analyzed with Maestro v.2019-3 (Schrödinger, New York, NY 10036, U.S.A.) and Pymol v.2.1 (Schrödinger, New York, NY 10036, U.S.A.).

  ml) was stirred at room temperature for 12h. The completion of the reaction was monitored by LCMS. The mixture was purified by HPLC (Method B) to give compound 47 (140 mg, 44% yield) as a

	white solid. 1H NMR (400 MHz, DMSO-d6) δ 8.48 (d, J = 1.8 Hz, 1H), 8.41 (d, J = 2.4 Hz, 1H), 8.12 -8.05
	(m, 2H), 7.71 (t, J = 2.2 Hz, 1H), 7.21 (dd, J = 8.6, 7.1 Hz, 2H), 7.15 -7.09 (m, 2H), 6.98 (d, J = 7.8 Hz,
	2H), 6.76 (t, J = 7.3 Hz, 1H), 4.16 (q, J = 7.0 Hz, 2H), 3.82 (d, J = 12.1 Hz, 2H), 2.88 -2.71 (m, 2H), 2.54
	(d, J = 1.0 Hz, 1H), 1.92 (d, J = 11.8 Hz, 2H), 1.81 (qd, J = 12.4, 3.9 Hz, 2H), 1.37 (t, J = 7.0 Hz, 3H). LC-MS
	(ESI) m/z 403.2 [(M+H)+, calcd. C 25 H 27 N 2 0 3 , 403.2].

DATA AVAILABILITY

Data. All scripts and input data necessary to run the workflow (CDK8 structures, ProCare subpockets alignment and scoring, fragment selection, DeLinker generative linking, library processing) have been made available at https://github.com/kimeguida/POEM. Software. ProCare (version 0.1.1) is available at https://github.com/kimeguida/ProCare. IChem (version 5.2.9) was downloaded from http://bioinfo-pharma.u-strasbg.fr/labwebsite/download.html.

mixture was poured into 30 ml of water and extracted with ethyl acetate (4x10 ml). The combined organic layers were washed with brine (20 ml), dried over anhydrous sodium sulfate, and concentrated under reduced pressure. The crude material was purified by HPLC (Method A) to afford compound 49 as a white solid (235 mg, 36% yield after 2 steps). [START_REF] Erlanson | Twenty Years On: The Impact of Fragments on Drug Discovery[END_REF] H NMR (400 MHz, DMSO-d6). (50). Compound 50 was prepared as described above for compound 48, starting from 5-bromopyridin-3-ol (300 mg, 1.7 mmol, 1 eq.) to afford a yellow solid (260 mg, purity 90%), which was used in the next step without further purification.

5-bromopyridin-3-yl 4-ethoxybenzoate

5-(1-phenyl-1H-pyrazol-4-yl)pyridin-3-yl 4-ethoxybenzoate (51).

Compound 51 was prepared as described above for compound 49, starting from 5-bromopyridin-3-yl 4-ethoxybenzoate 50 (200 mg, 0.62 mmol, 1eq.). The crude material was purified by HPLC (method B) to afford compound 51 as a white solid (50 mg, 8% yield after 2 steps). [START_REF] Erlanson | Twenty Years On: The Impact of Fragments on Drug Discovery[END_REF] H NMR (400 MHz, DMSO-d6). δ 9.21 (s, 1H), 8.95 (d, J = 1.9 Hz, 1H), 8.44 (d, J = 2.6 Hz, 1H), 8.39 (s, 1H), 8.17 In vitro CDK8 inhibition. Inhibitory activity of compounds was tested by using the LanthaScreen® Eu kinase binding assay optimized for CDK8/CyclinC (Invitrogen). This assay is based on the binding and displacement of an Alexa Fluor® 647-labeled ATP-competitive kinase inhibitor scaffold (kinase tracer) to the kinase. Binding of the tracer to the kinase is detected using a europium-labeled anti-tag antibody, which binds to the tagged CDK8/CyclinC. Simultaneous binding of both the tracer and antibody to the kinase results in a close proximity suitable for a high degree of FRET (fluorescence resonance energy transfer) from the europium (Eu) donor fluorophore to the Alexa Fluor® 647

ASSOCIATED CONTENT SUPPORTING INFORMATION

Supplementary Methods section and additional figures and tables including the comparison and alignment of sc-PDB subpocket and fragments to CDK8 ATP binding site, the colored feature histogram (c-FH descriptor) used to align sc-PDB subpockets to the target cavity, the validation of the subpocket comparison protocol, the pairwise similarity of selected fragments, the properties of selected fragments, the definition of connectable fragments, the topological requirements to connect fragment atoms by a linker, the filters for DeLinker-generated linkers, the properties of generated molecules, the LC-MS analysis of compounds 39, 41, 44, 47, 49 and 51, the ProCare alignment of selected fragment-bound subpockets to the CDK8 cavity, the topological relationships between the 385 fragments selected for the linking stage the synthesis of round-2 library compounds, the list of CDK8 X-ray structures, the filtering rules to select drug-like compounds, the in-house catalog of commercially available drug-like compounds, the list of 37 commercially available compounds structurally similar or identical to round-1 library members, the list of 151 round-2 library members (PDF).

Molecular formula strings-SMILES codes (CSV)

Docking pose of compound 49 (PDB file format) This material is available free of charge on the ACS Publications website at http://pubs.cas.org


Sandra Profice

Didier Lasseux

Yves Jannot

Naime Jebara

Gérald Hamon

Permeability, Porosity and Klinkenberg Coefficient Determination on Crushed Porous Media 1

Permeability estimation of poor permeable formations like tight or gas-shale reservoirs using a pulse-decay experiment performed on crushed samples has been shown in earlier works to be an interesting alternative for it is faster and less expensive than traditional transient tests performed on carefully prepared core plugs, although it is restricted to measurement in the absence of overburden pressure. Due to reservoir depletion during production, sample characterization over a wide range of pore-fluid pressure is essential. If the Darcy-Klinkenberg model is thought to be a satisfactory gas-flow model for these tight formations, the full characterization can be achieved by determining both the intrinsic permeability, k l , and Klinkenberg coefficient, b.

In this work, the conditions under which reliable estimates of k l , b and porosity, ϕ can be expected from this type of measurement are carefully analyzed. Considering a bed of monodisperse-packed spheres and a complete physical model to carry out direct simulations and inversion of the pressure decay, important conclusions are drawn opening wide perspectives for significant operational improvement of the method. In particular, it is shown that:

i) The particle size of the crushed sample must be well selected for a reliable pressure-decay signal record.

ii) The simultaneous determination of both k l and b by inversion of the pressure-decay signal is very difficult because the sensitivities of the pressure decay to both coefficients are correlated.

iii) The porosity of the particles can be accurately estimated when the experimental setup has been properly calibrated (volumes of the chambers and of the porous sample). The precision on the estimation of this parameter is however strongly dependent on a bias on the crushed sample volume.

iv) When identification of k l and b is possible, a very significant error may occur in the determination of the intrinsic permeability due to a bias on the porous sample volume. Errors on the estimated values of ϕ and k l due to a bias on the chamber volume are not very significant. Moreover, b remains insensitive to bias on both the chamber and porous sample volumes.

INTRODUCTION

Gas shales have become a topic of major interest over the past 10 years due to the large natural gas resource they represent. However, reliable characterization of gas shales is a challenging task since these unconventional reservoirs have very tight pores and therefore specific flow properties, in particular very significant Klinkenberg effects, at least after a sufficiently long period of reservoir depletion. These special features call upon a careful identification of tractable and reliable methods to determine single-phase flow properties that are of major importance to evaluate accumulation and production capabilities of the reservoir. Depletion of the reservoir may take place over a wide range of pore pressure. Assuming the Darcy-Klinkenberg model is valid for accounting for single-phase gas flow in such a porous medium, a characterization of porosity, ϕ, intrinsic permeability, k l , and Klinkenberg coefficient, b, should provide relevant insights in the reservoir properties.

Since the permeability of shales ranges from tens of microdarcies to nanodarcies, the time needed to reach stationary flow while performing a steady-state experiment on such samples is relatively long (it typically varies as the inverse of k l ) even on short plugs. Whereas associated flow rates are small (and difficult to measure), it must be noted that the separate determination of k l and b requires several different experiments when performed under steadystate conditions. Although this technique is widely used and despite an alternative quasi-steady method proposed elsewhere [START_REF] Jannot | A New Method to Measure Gas Permeabilty of Weakly Permeable Porous Media[END_REF], methods based on fully unsteady techniques are appealing. In the late 1960s, [START_REF] Brace | Permeability of Granite Under High Pressure[END_REF] elaborated the pulse-decay technique, which consists of applying a pressure pulse at one end of a core plug while recording the response at the opposite end. Since then, numerous studies have sought to improve the investigative setup and associated data interpretation methods. More recently, [START_REF] Luffel | Matrix Permeability Measurement of Gas Productive Shales[END_REF] proposed the use of pulsedecay methods on crushed samples. It involves recording the unsteady step of a pycnometry-type experiment using a device schematically represented in Fig. 1.

This option is particularly appealing as it shortens the test duration by increasing the medium exchange area with the penetrating gas while decreasing the penetration depth. In addition, the cost of the test is reduced because a small volume of cuttings recovered during the drilling is sufficient for the analysis. However, experiments under realistic overburden pressures are impossible.

Several variations of this technique were envisaged. For instance, [START_REF] Egermann | A Fast and Direct Method of Permeability Measurement on Drill Cuttings[END_REF] proposed a technique based on the flow of a viscous oil within saturated cuttings leading to a compression of the residual-gas content. [START_REF] Luffel | Matrix Permeability Measurement of Gas Productive Shales[END_REF] and [START_REF] Lenormand | Permeability Measurement on Small Rock Samples[END_REF], among others, introduced similar methods relying on the gas desorption of entire plugs subjected to a suction or vacuum. Recently, [START_REF] Wang | Permeability Characterization on Tight Gas Samples Using Pore Pressure Oscillation Method[END_REF] worked on a pore-pressure oscillation setup based on the generation of a sinusoidal excitation at the upstream side of a sample that is preconditioned at a given pore pressure.

All these transient tests are highly sensitive to the porosity of the medium that governs the transient gas flow. If this parameter is considered as a given datum in the interpretative model, it must be known with extreme accuracy, otherwise, estimates of the remaining parameters will be severely biased. As a consequence, a method allowing the determination of all the three parameters (k l , b and ϕ) from a single experiment is highly desirable.

The purpose of the present work is to analyze the capabilities of a pulse-decay experiment performed on a crushed sample to achieve this goal. To do so, a physical description on a model configuration is implemented and validated by an approximated analytical solution. Both the direct and inverse procedures are derived, which allows a) investigation of the available information in the pressure-decay record that is strongly related to the radius of the particles of the crushed sample, b) analysis of the possibility of extracting the three parameters k l , b and ϕ from the pressure-decay signal using a sensitivity analysis and c) study of the impact of the chamber and sample volume bias on the estimated parameters when inversion is possible.

PROCEDURE Modeling

The pulse-decay experiment carried out on the crushed sample must be first described by the flow model.

Direct Model

The physical model built to describe the pulse-decay experiment on a crushed sample relies on several hypotheses. First of all, the gas is supposed to be ideal, a commonly used approximation that is appropriate for gases like N 2 or He at operating pore pressures. Furthermore, the flow is assumed to be isothermal and slow enough to prevent development of inertial (Forchheimer) effects. At the pore level, a firstorder slip flow is considered, yielding a Darcy-Klinkenberg flow model at the macroscale. Lastly, the crushed medium is assimilated to a monodisperse-sphere pack, the porous material being rigid, homogeneous isotropic and dry while pressure is assumed to be uniform in the spaces between spherical particles at any time. Under these circumstances, the flow within the pack is one-dimensional along the radial direction of each sphere of radius R 0 . By combining the mass-conservation equation with the ideal gas law and the momentum conservation equation (Darcy-Klinkenberg), it is described by the following equation on the pressure P at the position r in each particle at time t

P t 1 1 1 r 2 2 r r k b P P P ϕ . ( 1
)
The associated initial and boundary conditions are given by:

P 0 (t = 0) = P 0i , (2)
P (r R 0 , t = 0) = P 1i , ( 3
) P r r = 0 = 0 , ( 4
) P t r = R 0 P r r = R 0 (V 0 + DV ) = k 1 P 0 P 0 b S T 1 + . (5)
Equation 3 indicates that the system is initially at equilibrium with the surrounding atmosphere at P li prior to the emission of the pressure pulse at t = 0 (Eq. 2). The condition (Eq. 4) simply results from the symmetry of the problem while Eq. 5 expresses the balance between the mass flowrate coming in all the particles of total surface area S T and the mass flow-rate out of the volume of gas available to the particles, V 0 +DV. The interest is focused on P 0 (t) = P(r = R 0 ,t) which is the pressure in the reservoir V 0 +DV as measured in the experiment.

This initial boundary value problem under this complete form has no analytical solution. It was hence solved by means of a numerical procedure, specifically an explicit finite difference scheme of first-order in time and second-order in space, inspired from previous works on the pulse-decay experiment on core plugs [START_REF] Jannot | A Detailed Analysis of Permeability and Klinkenberg Coefficient Estimation from Unsteady-State Pulse-Decay or Draw-Down Experiments[END_REF]. Using ψ = (P+b) 2 and the notation ψ i n for the nodal value of ψ at r = (i-1)∆r and t = (n-1)∆t where ∆r and ∆t are the grid spacing and time step respectively, the scheme is given by

t r 2 3 n+1 1 n 1 n 1 n 2 n 1 , ( 6
) t r n+1 i n i n i n i-1 1 1 r i r 1 1 r i r 2 r n i+1 n i [ ] ; 2 i m -1 (7) , n+1 m 1 1 4 n m r 2 t n m-3 n m 1 4 n m-1 n m-2 (8) , 0 i (P 1 + b) 2 1 i m -1 , ( 9
)
0 m (P 0i + b) 2 . ( 10
)
This scheme, in which we have denoted , and

k l V 0 R 0 2 r S T k l 1 2 3 n m
is subjected to a stability criterion on the time step, classical for a diffusion-like equation, such that

min i i r t 2 2
.

Approximated Analytical Solution

The numerical scheme was validated on a simplified configuration for which an analytical solution can be determined. Such a solution is achieved after linearization of Eqs. 1 to 5. Noting that for perfect gas, 1/P represents the fluid compressibility, denoted β f in what follows, and assuming that gas experiences only slight pressure fluctuation during the test, β f remains quasi constant. Under these circumstances, the apparent permeability, k = k l (1+bβ f ), is treated as a constant as well and the problem can be solved analytically.

The sequence of derivation of the solution splits in two parts. The set of equations is first Laplace transformed and solved in the Laplace domain. The solution is inverted making use of the Cauchy theorem, as detailed by [START_REF] Hsieh | A Transient Laboratory Method for Determining the Hydraulic Properties of Tight Rocks-1[END_REF] in a second step. Following this procedure, the analytical solution on P 0 (t) is obtained as

3 KR 0 2 + R 0 3 2 K P 0 (t) R 0 R 0 R 0 R 0 2 t K 3 K 2 m m m m m m cos sin sin Kexp 2 K 1/ 2 P 1i 2 (12) , with: R 0 2 (P 0i 2 _ P 1i 2 ) , ( 13
)
(V 0 + DV ) S T k f , ( 14
) and k K f . ( 15
)
The ξ m are the roots of:

tan m m m R 0 1 2 K . ( 16
)
In Fig. 2, we have represented the evolution of P 0 obtained from the numerical method along with the evolution predicted by the analytical solution given by Eq. 12. Parameters used for this comparison are k l = 10 -18 m 2 , b = 0 Pa, ϕ = 10%, R 0 = 10 -2 m, V 0 = 20 x 10 -6 m 3 , V 1 = 75.4 x 10 -6 m 3 , DV = 26.5 x 10 -6 m 3 , P 0i = 10 x 10 5 Pa, P li = 10 5 Pa, μ = 1.8 x 10 -5 Pa.s, t f = 60 s, NP = 1200, m = 500 (see nomenclature for the significance of the different parameters). The crushed sample volume V 1 was fixed to that of the standard core plug and the space trapped between the spheres, i.e. the dead volume (DV ), was calculated assuming that the spherical particles are packed according to a perfect facecentered cubic lattice whose compactness reaches 74%. Note that for real crushed samples for which particle shapes may strongly deviate from sphericity, larger compactness, i.e. smaller DV, may be achieved.

In order to estimate the deviation introduced by the hypothesis of a constant compressibility, the numerical procedure was first run without any particular assumption using the direct model Eqs. (1) to (5) and the scheme described above, yielding P 0 (t) represented as (P 0 ) 1 in Fig. 2a. In a second step, the numerical procedure was run using the same parameters with a constant β f equal to 1/P m , P m being the time average over t f of (P 0 ) 1 . The corresponding result is represented as (P 0 ) 2 in Fig. 2a. The analytical solution in this last situation is represented as (P 0 ) 3 in Fig. 2a.

The numerical and analytical results for a constant β f are in perfect agreement since the residue (Fig. 2b) remains globally close to 1 mbar for the chosen number of space nodes. The matching can even be improved when m is increased. This comparison clearly validates our numerical procedure. The comparison between the analytical solution and the numerical one simulated with a variable β f shows a much more significant discrepancy that increases with time as suggested in [START_REF] Wu | Gas Flow in Porous Media with Klinkenberg Effects[END_REF]. Hence, even if the pressure decay remains small during the test, the assumption of constant gas compressibility may introduce a significant error and should be used with great care for further parameter identification since uncontrolled bias may be introduced. The numerical procedure is now used to analyze the impact of the different experimental parameters on the pulse-decay test.

RESULTS AND DISCUSSION

In this section, the numerical procedure is used to investigate the impact of the particle size of the crushed sample on the pressure decay and to perform a sensitivity analysis of P 0 (t) to the parameters k l , b and ϕ. An inverse procedure based on the direct numerical model is further employed to analyze the bias on V 0 and V 1 on the estimated values of k l , b and ϕ when the estimation is possible.

Duration of the Experiment

Several numerical tests performed with a particle size R 0 of the crushed sample ranging from 1 mm to 1 cm were performed keeping all other parameters the same, namely, k l = 10 -21 m 2 , b = 68.6 x 10 5 Pa, ϕ = 1%, V 0 = 20 x 10 -6 m 3 , V 1 = 75.4 x 10 -6 m 3 , DV = 26.5 x 10 -6 m 3 , P 0i = 50 x 10 5 Pa, P 1i = 10 5 Pa, μ = 1.8 x 10 -5 Pa.s, t f = 200 s, NP = 2000 and m = 200. Results of the simulations are reported in Fig. 3.

From these results, it must be emphasized that the period over which the pressure decays and contains useful information about the physical parameters k l , b and ϕ is strongly related to the grain size. In fact, this characteristic dimension determines the test duration as it is reduced by a factor of ~10 when the particle radius is decreased from 1 cm to 1 mm. Thus, the smaller the spheres, the faster the decay, an acceleration that can be easily understood as the result of an increase of the exchange surface S T in conjunction with a reduction of the characteristic dimension to be penetrated by the pressure signal while reducing R 0 , V 1 remaining the same. However, just a few seconds of recording (as for R 0 = 1mm, see Fig. 3) are insufficient to guarantee an accurate post-estimation of the medium properties, taking into account actual pressure fluctuations at short times induced by valve opening and induced temperature effects resulting from the rapid gas expansion from V 0 to V 0 +DV which are not considered in this model. Consequently, for a given rock and predefined experimental conditions, the mean grain size must be carefully selected to provide a usable signal for the estimation of k l , b and ϕ by signal inversion. This dimension must be considered even more crucial when considering real particle shapes. Indeed, the hypothesis of spherical particles implies that the surfaceto-volume ratio is at minimum, hence the time to reach final pressure equilibrium is maximum. With particles of arbitrary shape having the same average dimension the available pressure-decay signal is expected to be even shorter.

Sensitivity Analysis

We now analyze the sensitivity of P 0 (t) to the parameters that are to be estimated, namely k l , b and ϕ. Several runs were performed with ϕ = 10%, R 0 = 10 -2 m, V 0 = 20 x 10 -6 m 3 , V 1 = 75.4 x 10 -6 m 3 , DV = 26.5 x 10 -6 m 3 , P 0i = 50 x 10 5 Pa, P 1i = 10 5 Pa, μ = 1.8 x 10 -5 Pa.s, t f = 12600 s, NP = 12600 and m = 200. The remaining parameters k l and b were those reported in Table 1. Numerical results obtained on P 0 (t) for the different (k l , b) pairs are represented in Fig. 4.

Figure 4, in which all pressure signals are superimposed, highlights a key point related to the pulse-decay experiment on crushed core plugs, namely, the difficulty to discriminate the independent role of k l and b on a given pressure decay. Indeed, the direct simulations show that the same signal can be obtained from tests performed with different combinations of the two parameters k l and b that can be varied over an order of magnitude (and even more). The physical explanation of this behavior is as follows. Because the pressure decay remains small, the Klinkenberg term b/P experiences small variations so that k l (1+b⁄P) remains quasi-constant during the test, making impossible the identification of the separate role of k l and b on P 0 (t). In other words, this indicates that k l and b appear 1.

as correlated parameters for the signal P 0 (t), a feature that is more pronounced when k l (and ϕ) are small, as for instance in the range of interest for tight or gas shale reservoirs. A clear illustration can be made from the sensitivity of P 0 (t) to each of the three parameters, k l , b and ϕ.

The reduced sensitivity (S r ) i of the measured pressure P 0 to the parameter θ i is defined by:

P 0 i i (S r ) i . ( 17
)
It basically expresses how the evolution of P 0 is modified upon a slight modification of the parameter θ i , all other parameters remaining the same. It should be noted that the identification of θ i from P 0 (t) is possible only if (S r ) i is larger than the sensitivity of the pressure sensor used to measure P 0 (t). An optimal estimation of θ i also requires the experiment be designed to maximize (S r ) i and minimize the θ i to θ j correlation.

Two different parameters θ i and θ j are correlated when their respective sensitivities (S r ) i and (S r ) j are proportional. In Fig. 5a, we have represented the sensitivity curves of P 0 (t) to k l , b and ϕ using k l = 10 -21 m 2 , b = 68.6 x 10 5 Pa, ϕ = 1%, R 0 = 10 -2 m, V 0 = 20 x 10 -6 m 3 , V 1 = 75.4 x 10 -6 m 3 , DV = 26.5 x 10 -6 m 3 , P 1i = 10 5 Pa, μ = 1.8 x 10 -5 Pa.s, t f = 1400 s, NP = 14000 and m = 200.

Clearly, the sensitivity of P 0 (t) to ϕ greatly exceeds those to k l and b, the latter being the lowest one. The large sensitivity to ϕ is indeed expected since the pulse-decay test in that case is a pycnometry-type experiment, a result that contrasts with the classical pulse decay test performed on an entire plug. All the three sensitivities remain large compared to the precision of a standard pressure sensor.

Figure 5b points out the correlation of k l and b for moderate pressure pulses P 0 i (~50 bars and less) through the quite perfect straight lines of the graph representing the sensitivity to b versus that to k l . Using larger pressure pulses, near 100 bars and more, slightly weakens the correlation in accordance with the fact that the apparent local permeability k l (1+b/P) experiences larger variations over the pressuredecay period. However, with such an experimental design, the correlation remains so significant that the simultaneous identification of the three parameters k l , b and ϕ is seriously compromised as further analyzed in the next section.

Inverse Procedure

The pressure-decay signal inversion is operated using a Levenberg-Marquardt numerical algorithm as in [START_REF] Jannot | A Simultaneous Determination of Permeability and Klinkenberg Coefficient from an Unsteady-State Pulse-Decay Experiment[END_REF]. This method which is commonly used to treat least square curve fitting problems returns the vector of parameters Θ ensuring the best correspondence between a set of measured data (t i , p i ) (here they are generated with the direct model) and the estimated values P(t i ,Θ) derived from the solution to the direct problem. Optimal estimated parameters are found by minimizing the quadratic residues:

S i = 1 2 p i _ P(t i , ) . ( 18
)
The inverse procedure was run for two configurations of permeability and porosity. In both cases, the pressure-decay signal was generated using the direct simulation detailed above while superimposing a Gaussian noise which amplitude was fixed according to the pressure-decay magnitude (see Tables 2 and3) to guarantee the convergence of the algorithm. Inversion was started giving two different initial sets (k li n , b i n , ϕ i n ) j . The (b i n ) j were chosen according to the correlation proposed by [START_REF] Jones | A Rapid and Accurate Unsteady-State Klinkenberg Permeameter[END_REF]:

(b in ) j = 0.189 (k 1in ) j -0.36 , ( 19
)
whereas the (ϕ in ) j were taken equal to the nominal value ϕ.

Note that since the final equilibrium pressure P f satisfies:

P 0 i V 0 + P 1i (DV + V 1 ) V 0 + DV + V 1 ) P f . ( 20
)
ϕ can be effectively pre-estimated with accuracy if the various pressures and volumes needed for the calculation are measured with care and the experiment is long enough to reach final equilibrium. Direct simulations and inversions were carried out with the following parameters: R 0 = 10 -2 m, V 0 = 40 x 10 -6 m 3 , V 1 = 75.4 x 10 -6 m 3 , DV = 40 x 10 -6 m 3 , P 1i = 10 5 Pa, μ = 1.8 x 10 -5 Pa.s and m = 200. The results of the inversion are gathered in tables 2 and 3.

As expected, the strong correlation between k l and b explained above justifies the fact that the estimated values are dependent upon the starting data. Despite the relatively small amplitude of the superimposed noise, errors on k l and b are quite large. This analysis clearly shows that: a) convergence can only be achieved if initial values of k l and b are close enough to the actual values, and b) when convergence is achieved, uncertainty of the estimated values of the parameters can be very large.

Although larger pressure pulses could be used to induce larger pressure decays, slightly reducing the correlation between k l and b, this may cause experimental artifacts due to pressure fluctuations and thermal effects associated to an abrupt pressure drop at the valve opening as well as mechanical effects on the porous particles. Since the small variation of the apparent permeability leading to correlation between k l and b is due to a large volume V 0 +DV compared to the pore volume ϕV 1 , an alternative is to considerably reduce V 0 and/or DV. In the remaining part of this work, V 0 is considered to be zero (more precisely, in the above numerical model, DV is set to zero while V 0 is taken as the actual value of DV) so that V s is directly connected to the gas supply. The actual experimental situation would be that of Fig. 1, for which the valve between V 0 and V s is opened in order to supply gas and generate the pulse and then closed so that the available gas volume is DV only during the entire test. This configuration is now considered and to reduce V 0 (i.e. DV) even more, a sample made of a single spherical particle having the volume of the standard core plug is chosen. Under these circumstances, reliable estimates of k l , b (and ϕ) are obtained since inversion identical results whatever the starting values of the parameters, provided that they lead to convergence. In this situation, the objective is now to investigate the error on the parameters k l , b and ϕ estimated by inversion due to bias in the volume of the chamber V 0 and in the sample volume V 1 that may always occur due to measurement errors in these two parameters.

To do so, pressure-decay signals were generated with the nominal values of V 0 and V 1 while respective bias of ±2 % in V 0 and ±6 % in V 1 were introduced to carry out the inversion. It should be noted that a bias in V 1 is unavoidably reflected in V 0 (V 0 = DV = V s -V 1 ). Tests were performed with k l = 10 -17 m 2 , b = 2.49 x 10 5 Pa, ϕ = 10%, R 0 = 2.6 x 10 -2 m, V 0 = 9.8 x 10 -6 m 3 , V 1 = 75.4 x 10 -6 m 3 , DV = 0 m 3 , P 0i = 50 x 10 5 Pa, P 1i = 10 5 Pa, t f = 30 s, NP = 300, μ = 1.8 x 10 -5 Pa.s and m = 200.

Results of the inverse procedure are reported in Table 4. They clearly show that the estimation of the Klinkenberg coefficient is insensitive to bias in both V 0 and V 1 . In addition, a bias in V 0 seems to be directly reflected in k l and ϕ, which means that the expected error due to this bias is not very significant provided the volume chamber is carefully characterized. However, the results show that the bias in V 1 has a very strong impact on both k l and ϕ, requiring that the volume of the crushed sample must be very carefully measured if one wishes to extract reliable information on both k l and ϕ from the pulse-decay signal.

CONCLUSION

A direct model and the related inverse procedure have been proposed to characterize in depth the physical properties of a crushed porous medium using a pulse-decay test.

From this analysis, it first appears that the available information in the pressure-decay signal is intimately related to the particle size of the crushed sample. If small particle radii are used, the pulse relaxation time, which is even reduced for increasing permeabilities, may be so small that no interpretation of the pressure-decay signal is possible. This suggests that a pre-estimation (for instance on an entire core plug) of the order of magnitude of both k l and ϕ would provide an advantage in defining the appropriate particle mean radius for a further pulse-decay experiment on the crushed sample. A sensitivity analysis performed on the pressure decay indicates that, for a standard experimental design, the identification of both k l and b is seriously compromised because these two parameters are strongly correlated. The physical origin of this lies in the restricted variation of the apparent gas permeability during the test. Without any pre-estimation (for instance on an entire core plug) of the order of magnitude of both k l and b, their estimation on crushed rock might be off by an order of magnitude. Conversely, and as expected, if volumes of both the chambers V 0 (and DV) and porous sample (V 1 ) are well characterized, a precise determination of the porosity can be performed as in a classical pycnometry-type experiment. Finally, a quick analysis of the impact of the bias in both V 0 and V 1 shows that when the inversion can be performed the resulting estimated values of b are unaffected by these bias. The effect of bias in V 0 on both k l and ϕ may not be very significant since it seems to be directly reflected as an error in the estimated values of these coefficients. In contrast, a bias in V 1 induces very significant errors in both k l and ϕ, suggesting that the volume of the crushed porous sample must be precisely measured for reliable identification of these two parameters from the pulse-decay signal.

Δk l /k l (%) ΔΦ/Φ (%) V 0 +ΔV 0 V 0 -ΔV 0 V 1 +ΔV 1 V 1 -

ACKNOWLEDGEMENTS

TOTAL is gratefully acknowledged for its support. She obtained an engineering degree from the National School of Electricity and Mechanics (Nancy), in Fluid Mechanics in 2010. Her research, carried out in collaboration with TOTAL, deals with the characterization of the physical properties of low permeable porous media, such as tight gas and gas shale reservoirs.

NOMENCLATURE

Didier Lasseux is a senior scientist at the Centre National de la Recherche Scientifique (Institute of Mechanics and Engineering in Bordeaux). He received an engineering degree from the Ecole Nationale Supérieure des Arts et Métiers and a PhD from the University of Bordeaux. He has been working on transport in porous media, more specifically on single and multiphase flow, for about 20 years. His recent research deals with the modeling of non Darcy flow and characterization of porous media properties. He has published more than 90 papers in journals and conference proceedings including SPE and SCA and has been granted 3 patents. He has co-founded the International Society for Porous Media (INTERPORE).

Yves Jannot is a senior scientist at the Centre National de la Recherche Scientifique (Laboratory of Theoretical and Applied Energetic and Mechanics). He received an engineering degree from the Ecole des Mines de Nancy and a PhD from the University of Nancy. He has been working on characterization of thermal properties of porous media for about 15 years. His recent research deals with the adaptation of methods for the characterization of thermal properties to the determination of the parameters related to fluid flow in porous media. He has published more than 60 papers in journals and conference proceedings and has been granted one patent.

Naime Jebara is a reservoir engineer at "Total SA" main technical and scientific research center in Pau, France. He is in charge of experimental pertrophysical characterization of conventional and unconventional special core analysis. He graduated from the Institut Français du Pétrole "IFP-School" with a MSc in petroleum engineering in 2009 and earned a degree in mechanical engineering from Institut National des Sciences Appliquées de Lyon (France) in 2008.

Gerald Hamon is currently an expert for Petrophysics with TOTAL. He earned his PhD in Fluid Mechanics in 1980 from the Institut National Polytechnique de Grenoble (France). He joined the oil industry in 1980 and worked in a variety of senior technical positions, including core analysis, field development in the North Sea and special expertise in the evaluation and numerical simulation of naturally fractured reservoirs. He authored or co-authored more than 60 papers for symposiums/journals of the Society of Core Analysts, the Society of Petroleum Engineers and the Society of Professional Well Log Analysts. In 2007, he received the Darcy Award from the Society of Core Analysts for lifetime technical achievement.

Fig. 1 -

 1 Fig. 1-Schematic experimental setup for a pulse-decay investigation of a crushed sample.

Fig. 2 -

 2 Fig. 2-(a)Numerical and analytical pressure-decay curves, and (b) residues drawn from the comparison of (P 0 ) 1 and (P 0 ) 2 to (P 0 ) 3 . (P 0 ) 1 , numerical pressure decay keeping a variable β f ; (P 0 ) 2 , numerical pressure decay for a constant β f ; (P 0 ) 3 , analytical pressure decay for a constant β f ).

Fig. 3 -

 3 Fig. 3-Pressure decay curves for different R 0 .

Table 1 -Fig. 4 -

 14 Fig. 4-Pressure decay curves for the different (k l ,b) pairs reported in Table1.

(

  P 0i = 5 x 10 5 Pa, t f = 100 s, NP = 1000) (pressure-decay amplitude: 18.7 mbars; -1.5 mbars <noise amplitude < 1.5 mbars) of k l , b and ϕ in Case 2 (P 0i = 5 x 10 5 Pa, t f = 100 s, NP = 1800). (pressure-decay amplitude: 42 mbars; -3 mbars <noise amplitude < 3 mbars) in ) 2 = 157.2 x 10 5 (Φ in ) 1 = 0.05 (Φ in ) 2 = 0.05

Fig. 5 -

 5 Fig. 5-(a) Sensitivity curves to k l , b and ϕ for P 0i = 50 x 10 5 Pa. (b) k l -b correlation curves (in absolute value) for different pressure pulses P 0i . (b)

Table 2 -

 2 Estimates of k l , b and ϕ in Case 1

Table 4 -

 4 Errors in the Estimates of k l , b and ϕ, Due to a Bias of 2 % in V 0 and 6 % in V l .

	Estimated Values	Errors
	k l (m 2 )	b (Pa)	Δb/b (%)

  Sandra Profice is a PhD student at the University of Bordeaux (Institute of Mechanics and Engineering).

		ABOUT THE AUTHORS	Δr	node spacing (m)
			Δt	time step (s)
			ΔX/X relative bias on X (-)
			ϕ	porosity (-)
			Θ	vector of parameters
			θ i	i th parameter
	b	Klinkenberg coefficient (Pa)
	DV	dead volume (m 3 )
	k	apparent gas permeability k = k l (1+b ⁄P) (m 2 )
	k l	intrinsic permeability (m 2 )
	m	number of space nodes (-)
	NP	number of measurement points (-)
	P/p i	gas pressure/i th experimental datum (Pa)
	P 0 i	initial pressure pulse (Pa)
	P 0 (t) pressure at r = R 0 and t (Pa)
	P 1i	initial steady-state pressure (Pa)
	P m	mean value of P 0 over t f (Pa)
	r	radial coordinate (m)
	R 0	particles radius (m)
	S T	particles total exchange area (m 2 )
	t	time (s)
	t i	date of measurement (s)
	t f	duration of the experiment (s)
	V s	sample cell volume (m 3 )
	V 0	reference cell volume (m 3 )
	V 1	crushed sample volume (m 3 )
	β f	ideal gas compressibility (Pa -1 )