Computer implemented method for calculating a bound-like structure of a protein

Abstract

An in-silico method to estimate one or more bound-like structure is presented, taking into account quasi-rigid domains and electrically charged residues of the binding site(s).

Description

[0001] P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI COMPUTER IMPLEMENTED METHOD FOR CALCULATING A BOUND-LIKE STRUCTURE OF A PROTEIN DESCRIPTION TECHNICAL FIELD The present invention refers to a computer implemented method for calculating a (plurality of) boundary structure(s) of a protein. STATE OF THE ART Knowledge of bound-like structure of a protein is achieved with accuracy by in vitro experiments and measures. SCOPES AND BRIEF DESCRIPTION OF THE INVENTION The scope of the present invention is to provide an in-silico prediction of the bound-like conformation of a protein. The scope of the present invention is achieved by a computer implemented method according to claim 1. In particular, checking the number of quasi-rigid domains of the protein, is an important step to increase the precision of the calculation and, hence, the prediction of bound-like structures. At the same time, the inventors spotted that including electrically charged and uncharged residues of the binding site, provides additional precision to the prediction. Such predictions, providing a reasonable number of putative bound-like conformations of a given protein, greatly improves the accuracy and ultimately the speed of docking calculations so as to significantly reduce the lead time for the design process of a drug. P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI BRIEF DESCRIPTION OF THE DRAWINGS The invention is described by non limiting embodiments shown in the drawings, which respectively refer to: : - Fig. 1 a simplified flowchart of the computer implemented method according to the invetion; - Fig. 2 information about exemplary residues of a benchamrk protein used to test the present invention - Fig.3 shows an example of principal inertial planes applied to a bonding site; - Fig. 4 shows a sketch to understand (pseudo) contacts among residues; - Fig. 5 show data about a numerical analysis according to the present invention where symbols are used to understand the graphs in b / w; - Fig. 6 is a visualization of a protein in apo configuration (4AKE) and holo configuration (1AKE), having three quasi-rigid domains, elaborated from output data of the present invention. Inhibitor molecule AP5 is also shown (sticks) and the arrows indicate the hinge-like movements between adjacent quasi-rigid domains. DETAILED DESCRIPTION OF THE INVENTION A general flowchart embodying the present invention is sketched in Figure 1. First, we identify the putative binding sites on a target protein. Figure 2 shows, as a mere example, a list of protein residues defining these sites. In real cases of interest, where the binding site is not known and only the apo structure of the protein (that is not bound to any drug-like molecule), the putative binding sites could be identified using publicly available site detection packages / webservers. For example, for the most targets including those used to validate the protocol, there is a good agreement between the experimental binding sites and those identified by the COACH-D webserver (REF 1). Such site finding algorithm P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI (other ones are known in the art) receives as input a model of the native structure of a protein and provides in output, according to a known approach, a list of putative binding amino acids. We then calculate the “inertia planes” of this selection of binding amino acids outputted by the site finding algorithm. Inertia planes are planes orthogonal to the corresponding inertia axes and passing through the geometrical center of the binding amino acids (Figure 3 – planes in phantom). According to a preferred embodiment, inertia axes are calculated via existing libraries in tcl scripting language. According to an alternative embodiment of the site-finding algorithm, the X- ray apo structure was identified from the Protein Data Bank (PDB). We selected, as a non limiting example, the structure with PDB ID 4AKE82. The structure was resolved at 2.2 Å resolution and did not feature any missing residue. The binding region was identified from the apo structure via the freely available site-finder software COACH-D (https: / / doi.org / 10.1093 / bioinformatics / btt447), without exploiting any experimental information. The software outputs a set of 10 possible binding sites, together with a score (C) reflecting the confidence of the prediction. The score ranges from 0 to 1, where a higher score indicates a more reliable prediction. The software identified three binding regions displaying a C- score respectively of 0.99 (BS1 COACH), 0.89 (BS2COACH), and 0.67 (BS3COACH), while the remaining predictions featured C-scores around 0.01. Notably, the second site was virtually equivalent to the union of the others. In general, we use as a binding site the one associated with the highest C score. Importantly, our choice was not biased towards any ligand-specific binding region; in fact, our method aims to be predictive without exploiting any prior knowledge about the target complex. To test the accuracy of the above protocol with respect to several experimental binding sites, we selected four different complex structures, each bearing a P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI different ligand and displaying full sequence identity to the apo one. The PDB IDs of these complexes are 1AKE125, 2ECK104, 1ANK80, 6F7U - https: / / doi.org / 10.1073 / pnas.1721508115). In 1AKE, the protein was resolved in complex with the inhibitor P1,P5-bis(adenosine-5’-)pentaphosphate (hereafter AP5), mimicking the presence of two physiological substrates and binding across the interfaces between the quasi-rigid domains that will be described in greater detail in a few paragraphs. In 2ECK, ADK bore two physiological substrates, ADP and AMP, while in 1ANK the protein was complexed with AMP, a physiological substrate, and with ANP, a non-hydrolysable ATP analog. Another non-hydrolysable ATP analog, GCP, was bound to the LC interface of the protein in the 6F7U structure. Four different experimental binding regions, labeled BSAP5, BSGCP, BSADP, and BSAMP, were defined by taking all the residues within 3.5 Å from the corresponding ligands AP5, GCP, ADP, and AMP in the structures identified by PDB IDs 1AKE, 6F7U, 2ECK, and 1ANK respectively. BSAP5 span both the NC and LC interfaces, while BSGCP and BSADP identify two slightly different regions across the LC interface, and BSAMP is located within the NC interface. Note that while enhanced conformational sampling will be performed on BSCOACH, the sampling performance will be assessed with respect to the experimental binding regions, which are truly relevant to ligand binding. Benchmarking our method in this way is crucial since not all the residues included in the experimental binding sites are retrieved in BSCOACH. Furthermore, collective variables of the residues (CVs) are spotted via the inertia planes, namely: three (pseudo)contacts across inertia plane (CIP) variables, each defined as the number of contacts between residues, i.e. binding amino acids, of the binding site on opposite sides of the corresponding inertia plane (Figure 4), and the gyration radius of the residues (RoGBS). The number of collective variables considered by the present method depends on the complexity of the input protein. A parameter significant to structure complexity is the P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI number of quasi-rigid dynamical domains of the protein structure (that is, groups of amino acids behaving as approximately rigid units in the course of protein equilibrium fluctuations, as established by a Normal Mode Analysis on the protein structure – REF: doi: 10.1093 / bioinformatics / btp512). In a simpler case where the binding site of the protein belongs to a single quasi-rigid domain, the above-mentioned 4 collective variables provide good results. For example, in such a case, then we perform relatively short bias-exchange, well-tempered metadynamics simulations (MD) of the apo protein using the set of four collective variables (CVs). In particular, standard all-atom MD simulations were carried out using the pmemd module of the AMBER16 molecular modeling software developed by Case, D. A.; Betz, R. M.; Cerutti, D. S.; Cheatham, T. E., III; Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W.; Homeyer, N.; Izadi, S.; et al; University of California: San Francisco, in 2016. Topology files were created for each system using the LEaP module of AmberTools17 starting from the experimental structures available in the Protein Data Bank, using standard forcefields. Missing parameters for the latter were generated using the antechamber module of AmberTools17. In particular, atomic restrained electrostatic potential charges were derived after a structural optimization performed with Gaussian 09 software by Gaussian Inc.. Each structure was solvated with the explicit TIP3P water model, and its net charge was neutralized with the required number of randomly placed K+ or Cl- ions. The total number of atoms was ∼86 000 for BGT / BGT-UDP, ∼54 000 for RIC / RIC-NEO, and ∼62 000 for ABP / ABP-ALL. Periodic boundary conditions were employed, and long-range electrostatics was evaluated through the particle-mesh Ewald algorithm using a real-space cutoff of 12 Å and a grid spacing of 1 Å per grid point in each dimension. The van der Waals interactions were treated by a Lennard-Jones potential using a smooth cutoff (switching radius 10 Å, cutoff radius 12 Å). The initial distance between P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI the protein and the edge of the box was set to be at least 16 Å in each direction. Multistep energy minimization with a combination of the steepest-descent and conjugate-gradient methods was carried out to relax internal constraints of the systems by gradually releasing positional restraints. Following this, the systems were heated from 0 to 310 K in 10 ns of constant-pressure heating (NPT) using the Langevin thermostat (collision frequency of 1 ps-1) and the Berendsen barostat. After equilibration, four production runs of 2.5 µs each (for a total of 10 µs for each system) were performed for the apo systems, while a single 1 µs-long simulation was performed for each complex. A time step of 2 fs was used for preproduction runs, while equilibrium MD simulations were carried out with a time step of 4 fs in the NPT ensemble (using a MC barostat) after hydrogen mass repartitioning. Coordinates from production trajectories were saved every 100 and 10 ps for MDapo and MDholo, respectively. Bias-exchange well-tempered metadynamics simulations were performed on the three apo proteins using known software e.g. the GROMACS 2016.5 package and the PLUMED 2.3.5 plugin. The last conformation saved from the equilibration step from MDapowas used as the starting structure for each simulation. AMBER parameters were ported to GROMACS using the acpype parser. To enhance the sampling of different binding site shapes, we used the following four CVs defined by including all heavy atoms of the residues lining the binding site itself: the radius of gyration of the binding site (RoGBS) calculated using the gyration built-in function of PLUMED and the numbers of (pseudo)contacts across the “inertia planes” (CIP1,2,3) of the binding site, defined as the planes orthogonal to the three principal inertia axes and passing through the center of mass of the binding site. Binding site residues were defined as those within 3 Å (BGT and RIC) or 4 Å (ABP) of the ligand in the experimental structure of the complex. The cutoff was increased for ABP-ALL because of the low number of residues (seven) found P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI when a 3 Å cutoff was used. Very similar definitions were found using the COACH-D Web server with the apo structures. The CVs were calculated by an in-house tcl script based on the VMD orient function. Namely, residues lining the binding site were split into two lists A and B according to the positions of the geometrical centers of their backbones on each of the two sides of the inertia plane, and the overall number of pseudocontacts Ncbetween the two groups was calculated through the coordination keyword of PLUMED, which implements a switching function such as the following: replica was simulated for 100 ns (as our aim is primarily to enhance sampling of different shapes of the binding site and not to obtain converged free energy profiles), so that each window accumulated 400 ns of simulation time. Coordinates were saved every 10 ps. The height w was set to 0.6 kcal / mol for all systems, while the widths si of the Gaussian hills were set according to established prescriptions107 to 0.15, 0.05, and 0.08 Å (RoGBS), 5.4, 4.8, and 1.6 (CIP1), 5.1, 3.2, and 4.9 (CIP2), and 5.3, 3.1, and 6.0 (CIP3) for BGT, RIC, and ABP, respectively. Hills were added every 2 ps, while the bias exchange frequency was set to 20 ps. The bias factor for well-tempered metadynamics was set to 10. The “windows” approach briefly described in the Results and Discussion was implemented using RoGBS as the control parameter. Namely, we applied restraints (force constants set to 50 and 10 kcal mol-1Å-2for the upper and lower walls, respectively, as we seek for compression rather than enlargement of the binding site) at values of RoGBS that were 7.5% higher and lower than the value measured in the apo X-ray structure (RoGapoX-ray). Then, P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI from the trajectory corresponding to this first window, we selected a random conformation of the protein whose RoGBS was 5% lower than RoGapoX-ray and performed another simulation with walls centered at ±7.5% RoGapoX-rayfrom this new center, repeating this procedure so as to simulate a total of four windows (see Figure 5). It should be noted that the walls were set to allow partial overlap between adjacent windows, which indeed occurred in all cases (Figure 5). We repeat this procedure to generate up to four windows including the first one. This leads to an overall reduction of RoGBS of 15% relative to the center of the first window RoGapoX-ray. Despite the arbitrariness of our choice, the performance of EDES is not very sensitive to the exact choice of three or four windows (and thus to the exact extent of the collapse induced at the binding site, amounting to 10% or 15% of the initial value, respectively). Moreover, although the imposed ∆ RoGBS correspond only to the true change seen for BGT, EDES performed comparably well for all of the systems investigated here. In the following we refer to the cumulative set of three- and four-window simulations as EDES3w and EDES4w, respectively. However, there may be structures where the residues of the binding site belong to at least a first and second quasi-rigid domain respectively. Grouping residues of the of the binding site on the first or the second quasi-rigid domain is performed by identifying which residue is on the first quasi-rigid domain and which one is on the second quasi-rigid domain. Afterwards, sub-groups of proximal residues are identified, e.g. residues of adjacent quasi-rigid domains that are within a threshold distance, e.g. 8 Å, from one another. Then a further collective variable of contacts is generated, i.e. ‘contacts between quasi-rigid domain cQR CV. For example, in the case of ADK, it is possible to implement a known analysis e.g. SPECTRUM https: / / doi.org / 10.1016 / j.str.2015.05.022 developed by Ponzoni, Polles, Carnevale, Micheletti to identify three quasi-rigid domains, namely the P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI CORE, LID, and NMP domains known from previous literature (Figure 6). To account for hinge-like motions between the two LID-CORE (LC) and NMP- CORE (NC) interfaces, we implemented the “contacts between quasi-rigid domains” (cQR) CVs representing the number of (pseudo)contacts across these interfaces. Clearly, with the aim of not relying on any experimental knowledge of the binding region, we only exploited the knowledge of BSCOACHto identify residues lining the NC and LC (sub)pockets. Namely, the following procedure was adopted: i) we selected all BSCOACH residues belonging to the NMP(LID) and being within a distance threshold e.g. 8 Å from any residue of the CORE; ii) a second specular selection was made by taking all BSCOACH residues belonging to the CORE that are within 8 Å from any residue of the NMP(LID). The union of selections i) and ii) defined the cumulative list of residues used to setup the cQRLC(NC) CV. iii) the residues associated with the NC and LC (sub)pockets were split into two lists via domain assignment, for which the number of (pseudo)contacts was calculated via the coordination keyword of PLUMED. Therefore, in the case of ADK, in view of the above, a total of six CV were considered so far. The inventors further identified the electric charge of residues as a condition to generate a further collective variable. In particular, within each subpocket, the charged amino acids were separated from the non-charged (others) ones, so as to obtain two CVs (cQRNC(LC)c and cQRNC(LC)o) specifically enhancing the conformational sampling of charged amino acids in targets such as ADK that contain many of them within the BS. Pointing to the general applicability of the method, we automatized the workflow so that the latter cQR subdivision occurs only if: i) the binding site presents more than the 25% of charged residues and ii) each residue group defining a cQR variable is composed of at least 2 non-adjacent residues. In this case, 11 out of the 32 residues (~34 %) composing BSCOACHare charged and this P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI subdivision was applied. The cutoff of 8 Å was chosen so that none of the BSCOACH residues got associated with more than one cQR CV and to ensure that each list contained a minimum of 4 residues belonging to each quasi-rigid domain (which limits the onset of large structural distortions within secondary structure elements). Thus, when more than two quasi-rigid domains within the same binding site and / or a sufficient number of charged residues is present in the BS, bias- exchange well-tempered metadynamics simulations were performed on the apo protein, biasing three CIPs, one gyration radius and one or more pseudo contact variable for each subgroup of residues within the threshold distance and respectively belonging to two adjacent quasi-rigid domains. We used the GROMACS 2020.4 packageand the PLUMED 2.6.2 plugin. AMBER parameters were ported to GROMACS using the acpype parser. According to the original implementation of the method, we defined four CVs considering the residues lining BSCOACH: i) the radius of gyration (RoGBS) calculated using the gyration built-in function of PLUMED; ii) the number of (pseudo)contacts across three orthogonal “inertia planes” (CIPs), calculated through a switching function implemented in the coordination keyword of PLUMED. In a first step, an unbiased MD shall be performed. In particular, standard all- atom MD simulations of the apo protein (hereafter MDstd) embedded in a 0.15 KCl water solution (~46.000 atoms in total) and under periodic boundary conditions were carried out using the pmemd module of the AMBER20 package by University of California, San Francisco. The initial distance between the protein and the edge of the box was set to be at least 16 Å in each direction. The topology file was created using the LEaP module of AmberTools starting from the apo structure (PDB ID: 4AKE). The AMBER-14SB force field was used for the protein, the TIP3P model was used for water, and the parameters for the ions were obtained from literature. Long-range electrostatics was evaluated through P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI the particle-mesh Ewald algorithm using a real-space cutoff of 12 Å and a grid spacing of 1 Å in each dimension. The van der Waals interactions were treated by a Lennard-Jones potential using a smooth cutoff (switching radius 10 Å, cutoff radius 12 Å). Multistep energy minimization with a combination of the steepest- descent and conjugate-gradient methods was carried out to relax the internal constraints of the systems by gradually releasing positional restraints. Following this, the system was heated from 0 to 310 K in 10 ns of constantpressure heating (NPT) MD simulation using the Langevin thermostat (collision frequency of 1 ps–1) and the Berendsen barostat. After equilibration, four production runs of 2.5 µs each were performed, for a total of 10 µs. A time step of 2 fs was used for pre- production runs, while equilibrium MD simulations were carried out with a time step of 4 fs in the NPT ensemble (using an MC barostat) after hydrogen mass repartitioning. Coordinates from production trajectory were saved every 100 ps. Furthermore, we employed variable (moving) restraints on this CV allowing smoother variations within a single simulation. Starting from the last conformation sampled along the pre-production step of the unbiased MD, each replica was simulated without restraints for the first 10 ns. Next, an upper restraint centered at the value of RoGX-rayapowas imposed with a force constant increasing linearly from 10 to 50 kcal mol-1Å-2in 0 ns. This preliminary phase is needed to push the system towards a structure featuring a RoGBSvalue close to RoGX-rayapo, disfavoring conformations with an overly enlarged BS. Then the center of the RoGBSrestraint is decreased linearly (every ns) to a value corresponding to 85% of RoGX-rayapoin 400 ns. As remarked in literature focusing on conformations featuring collapsed binding sites is based on the evidence that the binding of ligands to enzymes is most often associated with such structural changes. Finally, the RoGBS restraint is kept at its lowest value for further 100 ns. Thus, the cumulative simulation time of each replica (hereafter referred to as EDEStraj) amounts to 550 ns. Note that the relatively soft upper restraints used in P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI our protocol still allow the sampling of structures with a larger RoGBSwith respect to the value at which the restraint is set. The height w of the Gaussian hills was set to 0.6 kcal / mol, while their widths were set as reported in Table 1 on the basis of the fluctuations recorded during a short (~200 ps) unbiased MD run. The bias factor for well-tempered metadynamics was set to 10. Hills were added every 2.5 ps, while the bias- exchange frequency was set to 50 ps. Further details on CV definitions are reported in Table 1. Coordinates of the system were saved every 10 ps. We used the software fpocket to assess the druggability of the binding site within the ensembles of conformations generated by EDES. For each conformation, we evaluated the druggability score D, which ranges between 0 and 1 with higher values identifying more druggable geometries. It is customary to associate scores >0.5 to putative binding sites. Table 1 shows that EDES generated a much larger set of druggable structures than MDapofor BGT and ABP (actually, no druggable conformation was generated from this trajectory for the latter system), while the performance of the two sets was similar for RIC, as expected. In particular, (i) the EDES-derived ensembles have a higher percentage of structures associated with D > 0.5 than those derived from MDapo, and (ii) the percentage of structures with D > 0.9 is not much lower than that obtained from the MDholoset for BGT and RIC.

[0002] P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI Table I The cluster analysis, preferably an unsupervised clustering, was performed in the CVs space (for both the gEDES and the unbiased runs) as follows. The distribution of RoGBS values sampled during the MD simulation was binned into 30 equally wide slices, and the built-in hclust module was used to perform a hierarchical agglomerative clustering within each slice, setting the number of generated clusters in that slice to xi = (Ni / Ntot)·Nc, where Ni, Ntot, and Nc are the number of structures within the ith slice, the total number of structures, and the total number of clusters, respectively. In our case, Nc was set to 100. We then imposed as an additional requirement to have at least two clusters within each of the 30 slices. This was implemented by iteratively increasing Nc by 10 units until the number of clusters within each of the RoGBS slices was equal or higher than two. The resulting clusters were used as starting points to perform a second cluster analysis with the K-means method (maximum number of iterations set to 10000) and generating the same number of clusters. This multi-step strategy of clustering on the CV space outperforms a more standard RMSD-based approach in generating a maximally diverse ensemble of protein conformations. This P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI approach resulted respectively in 130 and 160 clusters for the gEDES and unbiased runs. Docking was performed with the software HADDOCK and AutoDock according to protocols known in literature. Ligand conformations were extracted from the relative complex structures and prepared according to the standard procedure of each software. The binding region selected for the calculations was BSCOACH, independently from the ligand and thus from its real binding site. For HADDOCK, a single docking run was performed per case, starting from the various ensembles of Nc conformations, with increased sampling (10000 / 400 / 400 models for it0, it1 and wat steps, respectively referring to rigid- body docking, semiflexible and final refinement in explicit solvent). Namely, during it0 the protein BS residues were defined as “active”, effectively drawing the rigid ligand into the BS without restraining its orientation. For the subsequent stages only the ligand was active, improving its exploration of the binding site while maintaining at least one contact with its interacting residues. The weight of the intermolecular van der Waals energy used in it0 was increased to 1.0 (from the default value of 0.01), and RMSD-based clustering was selected with a cutoff of 1 Å. Docking was guided by ambiguous distance restraints defined for the BS residues for and the ligand. For AutoDock, each ligand conformation was rigidly docked on each of the desired protein structures using the Lamarckian Genetic Algorithm (LGA). The grid density and number of energy evaluations were both increased from default values (respectively by decreasing the spacing parameter from 0.375 to 0.25 Å and by increasing the ga_num_evals parameter by a factor of 10) to avoid repeating each calculation several times to obtain converged results. An adaptive grid was used, enclosing all of the residues belonging to the BS in each different protein conformation. Finally, an additional step consisting in the relaxation of the docking poses by means of a multi-step structural relaxation performed with P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI AMBER20 software. Systems were optimized in vacuum through three consecutive cycles of restrained structural relaxation (1000 cycles of steepest descent followed by up to 24000 cycles of conjugate gradients) followed by an unrestrained optimization (2000 cycles of steepest descent followed by up to 8000 cycles of conjugate gradients). During restrained relaxation harmonic forces of 0.3, 0.2, and 0.1 kcal⋅mol-1⋅Å-1 respectively for the first, second, and third cycles were applied on all non-hydrogenous atoms of the system. Long-range electrostatics was evaluated directly using a cutoff of 99 Å, as for the Lennard- Jones potential. The AMBER-14SB force field was used for the protein, while the parameters of the ligands were derived from the GAFF force field1 using the antechamber module of AmberTools. In particular, bond-charge corrections (bcc) charges were assigned to ligand atoms following structural relaxation under the “Austin Model 1 (AM1)” approximation. After this step, poses were scored according to AutoDock’s energy function. Next, the top poses (in total Nc, one for each docking run performed on a different receptor structure) were clustered using the cpptraj module of AmberTools with a hierarchical agglomerative algorithm. Namely, after structural alignment of the BS for the different complex conformations, ligand poses were clustered using a distance RMSD (dRMSD) cutoff dc =0.075·Nnh, where Nnh is the number of non-hydrogenous atoms of the ligand. This choice was made to tune the cutoff to the molecular size of each compound. Finally, clusters were ordered according to the top score (lowest binding free energy) within each cluster. RESULTS AND DISCUSSION The workflow of gEDES protocol is sketched in Figure 1. The main new ingredients can be summarized as follows: i) the introduction of a new class of CVs (cQR) to specifically enhance the relative motions between quasi-rigid protein domains; ii) since the BS of ADK is rich in charged residues, which could be easily trapped into specific conformations by forming strong interactions, each P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI cQR variable was further split into two new ones, one encompassing only charged residues and one all the others; iii) the use of soft moving restraints to allow for a smooth decrease of RoGBSin a single run. In this section, we report the performance of our method in generating holo-like conformations of ADK, both in terms of reproducing the overall protein structure and the fine geometry of the binding sites of four ligands. Notably, sampling results are obtained without using any a priori experimental knowledge of the bound conformations. Moreover, the definition of the binding site used to drive the sampling (BSCOACH) was obtained using a freely available software without putting any bias toward specific chemotypes. The resulting BSCOACH comprised 32 residues (of which 12 are charged). After analyzing the sampling performance, we discuss how accounting for the plasticity of the protein does lead to improved docking outcomes. The present approach was tested on a flexible protein i.e. having three quasi-rigid domains and many electrically charged residues in two adjacent binding sites, is able to estimate many bound-like structures in a reasonable computational time i.e. less than 10 days using a modern GPU RTX4080 type.

Claims

P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI CLAIMS 1. Computer implemented method for predicting bound-like and "druggable" conformation of the binding-site of a protein, comprising the steps of: - receiving data representing a native structure of the protein - receiving or finding via a site-finding algorithm a plurality of residues of one or more binding sites of the proteing - calculating quasi-rigid domains of the native structure - in case of a single quasi-rigid domain, defining a first, second and third collective variable for contacts of residues of the binding site with respect to a first, second and third principal inertia plane of the binding site, and a collective variable about a gyration radius of the residues - performing a metadynamic analysis to the first, second, third and fourth collective variable to provide an output parameter representative of at least a bound-like structure of the protein.

2. Method according to claim 1, wherein in case of a first and second quasi-rigid domains within the binding site are provided, then: - identifying a first group of residues of the binding site belonging to the first quasi-rigid domain and a second group of residues of the binding site belonging to the second quasi-rigid domain; - performing a fifth collective variable including contacts between residues of the first and second group located at a distance below a pre- defined threshold distance; - applying a metadynamic analysis to the first, second, third, fourth and fifth collective variable to provide an output parameter representative of at least a bound-like structure of the protein.P7629PC00 UNIVERSITA’ DEGLI STUDI DI CAGLIARI 3. Method according to any of the preceding claims, further comprising the steps of: - identifying a first further collective variable including electrically charged residues of the binding site and a second further variable including non-charged residues of the binding site; wherein the step of performing includes said first and second further collective variable.

4. Method according to any of the preceding claims, wherein the step of performing generates a plurality of bound-like structures and further comprising the steps of: - clustering via an algorithm said plurality of bound-like structures to group the plurality of bound-like structures in a plurality of clusters - docking via an algorithm one or more ligands on a target bound-like structure extracted from one of the clusters.

5. Method according to claim 4, wherein the ligand is representative of a pharmaceutical active compound.

Images