Chasing allosteric inhibition of the SARS-CoV-2 PLpro via Molecular Dynamics Simulations with Flooding Fragments (MDFFr)

Chasing allosteric inhibition of the SARS-CoV-2 PLpro via Molecular Dynamics Simulations with Flooding Fragments (MDFFr) | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Chasing allosteric inhibition of the SARS-CoV-2 PLpro via Molecular Dynamics Simulations with Flooding Fragments (MDFFr) Jason Pattis, Khaled Elokely, Eleonora Gianti This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-7499914/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 19 Dec, 2025 Read the published version in Journal of Computer-Aided Molecular Design → Version 1 posted 9 You are reading this latest preprint version Abstract The SARS-CoV-2 papain-like protease (PLpro) represents a crucial therapeutic target due to its dual role in viral polyprotein processing and suppression of host immune responses through de-ubiquitination and de-ISGylation activities. To identify novel allosteric druggable sites on PLpro, we developed a molecular dynamics approach with flooding fragments (MDFFr), advancing previous Molecular Dynamics (MD) flooding simulation methods. Using MDFFr, we evaluated interactions of known phenolic inhibitors with SARS-CoV-2 PLpro and discovered multiple allosteric druggable sites. Our simulations not only validated experimentally known binding sites but also revealed previously unidentified hotspots, including protein-protein interaction sites for ubiquitin and ISG-15 (Interferon-Stimulated Gene 15). The MDFFr approach demonstrates robust capability for physics-based druggability assessment of biological targets using only protein 3D structure, while providing detailed insights into fragment-protein interactions at both druggable sites and protein-protein interfaces. These findings unveil new opportunities for allosteric inhibition of PLpro, potentially advancing therapeutic strategies against SARS-CoV-2 and other coronavirus-related diseases. Molecular Dynamics Simulations Flooding Simulations Binding Site Identification Cosolvent Mapping Protein-protein Interaction Sites Sites of Allosteric Modulation Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 INTRODUCTION The SARS-CoV-2 virus, a member of the beta Coronavirus (CoV) family,[ 1 – 3 ] emerged in 2019 and precipitated the COVID-19 pandemic, one of the most devastating global health crises in recent history. With over 780M documented cases and > 7.1M deaths worldwide (as of August 2025), its impact has been profound and far-reaching. Despite the World Health Organization's declaration ending the international public health emergency in May 2023, SARS-CoV-2 continues to evolve, generating variants of interest (VOIs) such as JN.1 and variants under monitoring (VUMs) including XFG, NB.1.8.1, and others that pose ongoing risks for new epidemic waves, as seen in the summer 2025 surges reported across the United States, Europe, and other regions. While current countermeasures include variant-targeted vaccines and FDA-approved therapeutics (both antiviral drugs and monoclonal antibodies),[ 4 , 5 ] the persistent threat of SARS-CoV-2 to public health systems necessitates the development of novel therapeutic strategies. Particularly crucial is the need for broad-spectrum interventions that can effectively target multiple circulating variants and mitigate the risk of future outbreaks. Among potential therapeutic targets against SARS-CoV-2, the papain-like protease (PLpro) emerges as a particularly promising candidate (Fig. 1 ), given its critical role in modulating host antiviral immune responses against emerging viral variants. Encoded by the nsp3 gene,[ 6 ] PLpro functions alongside the main protease (3CLpro/Mpro, encoded by nsp5) as one of two essential and highly conserved coronavirus cysteine-proteases. PLpro's primary function involves the processing of viral polyproteins (pp1a and pp1ab), which are synthesized from viral RNA by host ribosomes.[ 7 ] This crucial proteolytic processing generates 16 mature nonstructural proteins (nsps) that subsequently form the viral replication transcription complex (RTC). Given these vital functions, inhibition of PLpro's catalytic activity represents a compelling strategy for antiviral drug development.[ 8 ] Beyond its proteolytic function, PLpro plays a crucial role in suppressing host innate immunity through its deubiquitinase (DUB) activity. This immunomodulatory function is achieved by cleaving ubiquitin (Ub) and ubiquitin-like (UbL) proteins, particularly interferon stimulating gene product 15 (ISG15), from host proteins.[ 9 , 10 ] Specifically, PLpro recognizes and cleaves the C-terminal RLRGG motif of Ub/UbLs, thereby removing these modifications from cellular proteins. This DUB activity significantly impacts host immune responses by dysregulating both inflammatory and interferon-mediated pathways, a mechanism that has been directly linked to increased SARS-CoV-2 mortality rates.[ 11 , 12 ] SARS-CoV-2 PLpro (SCoV2-PLpro) exhibits remarkable conservation among coronaviruses, sharing 83% sequence identity with its SARS counterpart (SCoV1-PLpro) and maintaining 100% conservation of catalytic-site residues crucial for Ub interactions.[ 13 ] A key distinguishing feature of SCoV2-PLpro is its preferential activity toward ISG15, particularly in cleaving ISG15 from interferon responsive factor 3 (IRF3), whereas SCoV1-PLpro primarily targets ubiquitin chains. This specificity directly impacts type I interferon responses, making SCoV2-PLpro an attractive therapeutic target with potential dual benefits: inhibition of viral infection and enhancement of antiviral immunity.[ 12 ] Structural studies have revealed that CoV PLpro comprises four distinct domains: a ubiquitin-like domain (UBL) and three domains (thumb, fingers, and palm) that together form the ubiquitin-specific proteases (USP) domain (Fig. 1 ). The enzyme's active site specifically recognizes the LXGG/XX substrate motif at the nsp1/2, nsp2/3, and nsp3/4 junctions,[ 9 , 14 , 15 ] utilizing a catalytic triad of Cys111, His272, and Asp286.[ 16 , 17 ] This site catalyzes two crucial reactions: (1) cleavage of viral polyprotein after the P1 glycine, releasing nsp1, nsp2, and nsp3 for viral replication, and (2) hydrolysis of isopeptide bonds following the LXGG motif at the C-terminus of cellular Ub and ISG15,[ 9 , 10 ] thereby modulating innate immune responses during infection and pathogenesis.[ 18 – 20 ] The protein features two distinct ubiquitin-binding sites: the S1 site (Ub-S1), positioned between the fingers and palm domains, which accommodates ubiquitin and the N-terminal UbL domain of ISG15 (Fig. 2 C), and the S2 site (Ub-S2), located at the α2 helix-thumb domain interface, which binds the C-terminal UbL domain of ISG15 (Fig. 2 B). The development of effective PLpro inhibitors has seen significant progress, though the number of validated compounds remains limited.[ 21 ] The first notable inhibitor was GRL0617, a potent non-covalent naphthalene-based compound that targets the active site.[ 16 ] X-ray crystallography revealed its unique mechanism of action, where inhibition occurs through induced loop closure of the active site.[ 22 ] A key challenge in PLpro inhibitor development has been the limited druggability of the 'site I' and 'site II' sub-pockets, primarily due to their featureless Gly-Gly recognition pattern (LXGG substrate).[ 21 ] Recent advances have addressed this challenge through multiple strategies: (1) leveraging binding cooperativity across multiple sites, leading to 2-phenylthiophenes with nanomolar potency,[ 21 ] (2) development of selective covalent inhibitors,[ 23 ] and (3) design of novel compounds targeting the ubiquitin-binding sites. Recent breakthroughs include the discovery of orally bioavailable inhibitors with efficacy against multiple SARS-CoV-2 variants and related coronaviruses, demonstrating PLpro's continued promise as a therapeutic target. As of August 2025, PLpro inhibitors for SARS-CoV-2 are still in preclinical or early development stages and lack clinical approval. Approved outpatient therapies mainly target alternative viral proteins, such as Mpro with nirmatrelvir-ritonavir, or employ nucleoside analogs. These include nirmatrelvir-ritonavir (Paxlovid), an oral antiviral for adults and children aged ≥ 12 years (≥ 40 kg), initiated within 5 days of symptom onset; remdesivir (Veklury), an intravenous antiviral for adults and children aged ≥ 28 days (≥ 3 kg), given over 3 days and started within 7 days of symptoms; molnupiravir (Lagevrio), an oral antiviral for adults, begun within 5 days of symptoms (with potential expiration of distributed supplies); and convalescent plasma for immunocompromised patients under emergency use authorization. Symptomatic supportive care, such as antipyretics and analgesics, is additionally advised. An alternative approach to address the limited druggability of sub-pockets involves targeting allosteric sites through multiple mechanisms, including disruption of protein-protein interactions and interference with Ub/UbL substrate binding sites.[ 13 , 24 ] This strategy has yielded promising results, as demonstrated by the discovery of two fragment-like natural compounds through X-ray crystallography-based high-throughput screening: 4-(2-hydroxyethyl)phenol (YRL) and 4-hydroxybenzaldehyde (HBA). These compounds bind at the interface between the thumb and UBL domains (Fig. 2 A), where they are thought to disrupt the S2 ubiquitin binding site (Ub-S2).[ 25 ] Additionally, a third compound, 3,4-dihydroxybenzoate (HE9), was co-crystallized in a distinct location exposed to solvent at the back of the α2-helix (Fig. 2 A). Both HE9 and YRL demonstrate significant inhibitory activity against PLpro, with EC50 values of 0.13 µM and 10 µM respectively, establishing them as valuable starting points for further drug development.[ 25 ] Fragment-based drug discovery has emerged as a powerful approach for identifying novel therapeutic compounds through high-throughput screening of fragment libraries against biological targets. This method typically employs biophysical techniques, including differential scanning fluorimetry and surface plasmon resonance, followed by detailed structural characterization using NMR spectroscopy or X-ray crystallography. The identified fragments are then optimized through growing, linking, or merging strategies. However, experimental characterization of protein-fragment interactions presents significant challenges and limitations. To address these constraints, numerous computational approaches,[ 26 – 29 ] including our previous work,[ 30 – 32 ] have been developed to enhance various aspects of fragment-screening campaigns. While current computational methods, primarily based on molecular docking and virtual screening, excel at identifying active scaffolds and chemical leads, they have important limitations. These approaches often provide insufficient insight into two critical aspects: (1) the dynamic effects of ligand binding on functional or allosteric protein pathways, and (2) the identification of druggable pockets on protein surfaces. These limitations can only be effectively addressed through molecular simulations,[ 33 , 34 ] which provide crucial dynamic information about protein-ligand interactions. In this study, we developed and implemented Molecular Dynamics simulations with Flooding Fragments (MDFFr), a computational protocol applied to SCoV2-PLpro over a multi-microsecond timescale (2.7 µs total). This approach serves two primary purposes: (1) computational characterization of fragment ligand binding at experimentally validated locations, and (2) discovery of novel allosteric binding sites. Our method builds upon existing protocols that incorporate protein flexibility in hot-spot mapping, such as those using 'standard probes' (isopropanol, phenol, benzene, etc.)[ 35 – 38 ] and the MixMD technique, which analyzes cosolvent density distributions across MD trajectories.[ 39 – 43 ] We also drew inspiration from "MD flooding" methodology, previously employed to study volatile anesthetic binding pathways in membrane proteins.[ 44 – 47 ] MDFFr represents an advancement of these approaches, offering a generalized MD flooding protocol applicable to both globular proteins and enzymes. Using this protocol, we conducted multiple unbiased simulations of SCoV2-PLpro using high concentrations of experimentally validated co-crystal fragment ligands ("flooding fragments") to accelerate binding events and enhance ligand sampling efficiency.[ 44 – 47 ] Our analysis revealed several key findings: Multiple sites suitable for allosteric ligand binding exist on SCoV2-PLpro, including both experimentally confirmed sites[ 25 ] and previously unidentified druggable spots. The flooding fragments showed preferential binding to locations that align with experimental screening results.[ 25 ] Additional "ligand hot spots" were identified that correspond to ubiquitin and ISG-15 interaction sites.[ 9 , 10 ] These results not only unveil new opportunities for allosteric inhibition of SCoV2-PLpro—a crucial strategy in combating COVID-19—but also establish MDFFr as a powerful tool for computational fragment-based screening and binding site identification. MATERIALS AND METHODS Structural Models . Three high resolution crystal structures of SARS-CoV-2 PLpro 7OFS (1.9 Å), 7OFT (1.95 Å), and 7OFU (1.72 Å),[ 25 ] co-crystallized with phenolic fragment-like compound inhibitors YRL (4-(2-hydroxyethyl)phenol), HBA (p-hydroxybenzaldehyde), and HE9 (methyl 3,4-bis(oxidanyl)benzoate), respectively, were prepared for unbiased, all-atom molecular dynamics (MD) simulations. Co-crystallized water molecules and fragment-ligands were removed, before the protein structures were protonated using the pdb2pqr program.[ 48 ] Cystines 189, 192, 224, and 226 were left deprotonated because they coordinate a zinc molecule. Each protein structure was then placed in a cubic box with a 10 Å buffer. The GROMACS[ 49 ] insert-molecule tool was used to place 10 millimolar (mM), 50 mM, or 100 mM concentration of non-overlapping, individual fragment ligands in the box. (The “insert-molecule” tool places a molecule randomly in the box but only accepts the placement if there is no overlap with other molecules.) Each box was solvated with TIP3P[ 50 ] explicit solvent before adding 150 mM of NaCl. A total of 9 systems were prepared, leading to a cubic periodic simulation cell of \(\:113\times\:\:113\times\:\:113\:\) Å for a total of 144,343 atoms for 7OFS, 10 mM; \(\:113\times\:\:113\times\:\:113\) Å with a total of 144,244 atoms for 7OFS, 50 mM; \(\:113\times\:\:113\times\:\:113\:\) Å with a total of 144,101 atoms for 7OFS, 100 mM; 111 x 111 x 111 Å with a total of 137,557 atoms for 7OFT, 10 mM; 111 x 111 x 111 Å with a total of 137,476 atoms for 7OFT, 50 mM; 111 x 111 x 111 Å with a total of 137,269 atoms for 7OFT, 100 mM; 111 x 111 x 111 Å with a total of 137,306 atoms for 7OFU, 10 mM; 111 x 111 x 111 Å with a total of 137,182 atoms for 7OFU, 50 mM; and 111 x 111 x 111 Å with a total of 137,022 atoms for 7OFU, 100 mM. Molecular Dynamics Simulations. Systems were minimized with 5000 steeps of steepest decent energy minimization. Position restraints of 1000 KJ/mol were placed on each protein and fragment-ligand heavy atoms before conducting 0.5 ns of restrained NPT equilibration runs at 1 bar with the Parrinello-Rahman barostat[ 51 ] and 300 K with the Nose-Hoover thermostat. One additional run of equilibration for 0.5 ns under restrained NVT conditions was performed. Simulations were carried out using GROMACS[ 49 ] version 2021.2. The Charmm36m[ 52 ] force field was used to simulate all systems with the CHARMM generalizable force field (CGenFF)[ 53 – 55 ] used to parameterize the fragment-ligands. Production runs were carried out in the NVT ensemble at 300 K, reaching 0.3 µs for each system, for a total of 2.7 µs of simulation (i.e., 0.3 µs × 9 systems). Coordinates of the systems were collected every 10 ps for a total of 30,000 frames for each run. Analysis of the Results RMSD and RMSF . The Python library MDTraj,[ 56 ] version 1.9.3, was used to compute the root-mean-square deviation (RMSD) of protein and ligand atoms relative to their respective reference experimental structure over the simulation time.[ 16 ] Similarly, root-mean-square fluctuation (RMSF) of protein and ligand atoms were calculated using GROMACS.[ 49 ] All plots were generated using the Python library Matplotlib,[ 57 ] version 3.5.0. PCA . The essential motions of the simulated systems were captured by the global principal component analysis (PCA), performed using PyEMMA,[ 58 ] version 2.5.7, on the Cartesian coordinates of protein backbone atoms. PCA was fit on all nine trajectories to capture the highest variance motions of any system. In PCA plots, clusters of samples were color coded by the density of frames. A local PCA of the YRL binding site was performed (using PyEMMA)[ 58 ] by first selecting all alpha- and beta-carbons located within 9 Å of the YRL-ligand as in the experimental structure. Pairwise distances between these alpha- and beta-carbons were measured in all 3 systems with YRL (i.e., the protein extracted from PDB-ID 7OFS and YRL ligands at 10-, 50- and 100-mM concentrations). Then, PCA was performed on these distances, extracting the largest variance motions. Principal components are the eigenvectors of the data’s covariance matrix. Domain Angle Analysis. To monitor the extent of the movement of the UBL domain with respect to the remaining ones, domain angles were defined using MDTraj[ 56 ] by calculating the angle (in radians) between the alpha carbons for three-residue sets. Two sets of residues were used (Figures S4-S5), each with one residue located in the UBL domain and the remaining two, including the vertex, in highly flexible regions of the thumb domain, i.e. – the \(\:\alpha\:2\) - and \(\:\alpha\:3\) -helix. Clustering and Occupancy . To find fragment binding spots, first, the trajectory was aligned to the protein backbone using the starting structure as the reference. The center of mass (COM) of each fragment was calculated using MDTraj.[ 56 ] Unsupervised clustering was performed on the fragment COM using the Density-Based Spatial Clustering of Applications with Noise algorithm, DBSCAN,[ 59 , 60 ] from the machine learning scikit-learn python package,[ 61 ] version 0.22.1. To belong to the same cluster, a 1.5-Å maximum distance between the COM of two fragments was set as a threshold, for a minimum of 60 ns over 300-ns individual simulations. The center of each cluster is rendered in Visual Molecular Dynamics (VMD) molecular visualization program[ 62 ] as a red van der Waals bead. Lastly, the equation below (Eq. 1) was then used to convert the total number of frames in a cluster into occupancy, or total amount of time spent in that spot. Global occupancy was calculated over all trajectories (2.7 \(\:\mu\:\) s). \(\:Occupancy\:\%=\:\frac{\#\:Frames\:in\:a\:cluster}{Total\:Trajectory\:Frames}\:\) (Eq. 1) Cavity Detection . Pocket volume on MD trajectories was measured using the cavity detection and characterization tool MDpocket.[ 63 ] Images were produced using VMD,[ 62 ] version 1.9.3 or with The PyMOL Molecular Graphic System, version 2.5.2. (Schrödinger, LLC). Water thermodynamics . Water mapping calculations were carried out using SZMAP Version 1.6.3.0 from OpenEye Scientific software.[ 64 – 66 ] The simulation frames containing a fragment-ligand in a given spot were collected, and the computation was performed on the simulation frame that was closest, in RMSD space, to the average structure of a particular spot. SPRUCE 1.4.0.0 (Spruce 1.6.0.3. OpenEye, Cadence Molecular Sciences, Santa Fe, NM) was used to generate biologically relevant units, add missing sidechains, adjust protonation states, enumerates alternate locations and tautomerization states for protein, ligands and co-factors. The output file has had charges and radii assigned. SZMAP grids were calculated for the binding pocket regions around the fragment. Then the stabilization calculations were performed for the complex, apo and ligand grids to describe where water molecules are influencing the binding process. The neutral difference in the free energy between the water and uncharged hypothetical water probes was computed to map the polarity throughout the binding site. Contact analysis . The number of contact pairs was defined as the number of protein residue pairs with at least one heavy atom within 4 Å from each other. Contacts were calculated for all frames in which at least one fragment-ligand was found at any given binding spot, then averaged and rounded up to the nearest whole residue using PyEMMA.[ 58 ] Exposons . Cooperative changes in the solvent exposure of protein residues were calculated using the exposons method,[ 67 ] which identifies functional conformational changes in protein structure ensembles. The method relies on identifying clusters of protein residues that are subjected to cooperative changes in their solvent exposure – i.e., the exposons – and deciphering the hierarchy of interactions that occur between them. Hence, in exposons calculations, the solvent accessible surface area (SASA) per residue-sidechain is computed. Then, correlations between pairs of residues are quantified by calculating a pairwise mutual information matrix. Next, the mutual information matrix is clustered into groups using affinity propagation clustering[ 68 ] with a damping value of 0.9. These groups of residues are often cryptic pockets , pockets which only open in rare events. Exposons was calculated using the enspara python package.[ 67 ] The exposons calculation was performed on the combined set of all nine trajectories, analyzing the protein only. Exposons were ranked by the number of residues involved in the exposon. RESULTS Conformational dynamics of SARS-CoV-2 PLpro via MD “flooding” simulations . We performed nine independent 300-ns "flooding" MD simulations of SCoV2-PLpro bound to different fragment-ligands, accumulating a total simulation time of 2.7 µs. Three distinct fragments were investigated:[ 25 ] YRL (4-(2-hydroxyethyl)phenol) from PDB ID: 7OFS, HBA (p-hydroxybenzaldehyde) from PDB ID: 7OFT, and HE9 (methyl 3,4-bis(oxidanyl)benzoate) from PDB ID: 7OFU. The chemical structures of these fragments are illustrated in Fig. 2 D. Each fragment was simulated at three different concentrations: 10 mM, 50 mM, and 100 mM (Table 1 ). For clarity, we designated these as different systems, where systems 1–3 correspond to YRL at increasing concentrations, systems 4–6 to HBA, and systems 7–9 to HE9. The protein dynamics remained largely stable across most systems, with root-mean-square deviation (RMSD) values typically ranging between 2–3 Å (Figure S1). However, systems 5 and 6 (HBA) and 8 (HE9) exhibited higher RMSD values, though this increased mobility appeared independent of both concentration and fragment type. Root-mean-square fluctuation (RMSF) analysis revealed a similar pattern: while systems 1–3 (YRL) maintained consistent residue deviations across increasing ligand concentrations, systems 5–6 (HBA) and 8 (HE9) displayed notably higher RMSFs in the UBL domain (residues 1–62) (Figure S2). Additionally, systems 5 and 8 showed elevated RMSFs in the thumb domain (residues 120–170), though this behavior appeared independent of fragment type or concentration. To further explore the protein dynamics and capture essential motions, we conducted principal component analysis (PCA) on the Cartesian coordinates of the protein backbone atoms (Fig. 3 ). The analysis revealed that the first three principal components accounted for 59% of the total variance, with PC1 alone explaining 34% (Fig. 3 B). The strongest correlations with PC1 corresponded to motions of residues 23, 46, 26, 25, 24, and 45, all within the UBL domain. Examination of frames with minimum and maximum PC1 values showed nearly perfect alignment of residues 63–315, excluding the UBL domain, indicating that PC1 primarily represents the forward movement of the UBL domain (Fig. 3 A, left panels). The second principal component (PC2) accounted for 18% of the variance (Fig. 3 B), with strongest correlations among residues 18, 17, 19, 9, and 31, also located in the UBL domain. Analysis of minimum and maximum PC2 structures revealed that this component corresponds to a downward movement and forward twisting of both the UBL domain and the adjacent helix-3 in the thumb domain (Fig. 3 A, right panels). While systems 1, 2, 3, 4, 7, and 9 maintained stability in a single energy basin within the PC1/PC2 space, systems 5, 6, and 8 explored multiple energy basins, moving in negative PC1 and PC2 directions, indicating substantial UBL domain movement. The third principal component (PC3) contributed 7% of the variance (Figure S3), with strongest correlations among residues 226, 225, 227, 224, and 191. These residues, located in the fingers domain, form part of the loops containing four cysteine residues responsible for tetrahedral coordination of a zinc ion (Zn 2+ ) in a zinc finger motif. This motif has been established as essential for maintaining structural integrity and protease activity in CoV PLpro enzymes.[ 9 , 69 ] The motion captured by PC3, which was evenly sampled across all systems (Figure S3), corresponds to an inward curling of the fingers domain toward the active site. Table 1 Summary of MDFFr simulations and spots identified (total time 2.7 \(\:\varvec{\mu\:}\varvec{s})\) . System PDB-ID Fragment Concentration (mM) Simulation time (ns) Spot # Sys-1 7OFS YRL 10 300 S11 Sys-2 7OFS YRL 50 300 - Sys-3 7OFS YRL 100 300 S1, S4, S5, S7, S11, S12 Sys-4 7OFT HBA 10 300 - Sys-5 7OFT HBA 50 300 S4, S9, S10 Sys-6 7OFT HBA 100 300 S11 Sys-7 7OFU HE9 10 300 S2, S11 Sys-8 7OFU HE9 50 300 S6, S8, S9 Sys-9 7OFU HE9 100 300 S3, S11 To validate the UBL-related motions identified by PC1 and PC2, we analyzed specific inter-domain angles throughout the simulation trajectories. We monitored two key angles: (1) between the alpha carbons of residues 134→84 (thumb domain)→8 (UBL domain) to capture the forward motion of the UBL domain (Figure S4), and (2) between residues 88→81 (thumb domain)→5 (UBL domain) to track the downward UBL movement (Figure S5). The forward motion analysis revealed distinct behavioral patterns among the systems. Systems 1–4, 7, and 9 maintained stable domain angles throughout the simulations (Figure S4B). In contrast, system 5 exhibited consistently large increases in domain-angle amplitude across the entire trajectory, while systems 6 and 8 showed significant angle increases toward the end of their simulation periods. The downward motion analysis showed similar system-dependent behavior. The angle tracking the UBL domain's downward movement remained stable in systems 1–4, 7, and 9, while systems 5, 6, and 8 demonstrated pronounced downward UBL motion in the latter portions of their simulations (Figure S5B). These angular analyses provide quantitative support for the domain motions identified through PCA, confirming the significant conformational flexibility of the UBL domain in specific systems. Identification and Characterization of Binding Spots for Fragment Ligands. We identified and characterized potential fragment binding sites on SARS-CoV2-PLpro through systematic analysis of our MDFFr simulations. The binding sites were determined by clustering the center of mass positions of flooded fragments throughout the simulation trajectories, revealing a total of 12 distinct binding spots (Fig. 4 ). Each identified spot underwent comprehensive evaluation using multiple parameters: Simulation occupancy: The duration of fragment binding during individual simulations. Global occupancy: The percentage of total simulation time (2,700 ns) during which a fragment remained bound. Contact analysis: The number of protein residues within 4 Å heavy atom distance of the bound fragment. Water thermodynamics: The local solvent behavior and energetics at each binding site This multi-parameter assessment approach enabled us to characterize not only the location of potential binding sites but also their physicochemical properties and binding stability, providing insights into their potential druggability. Occupancy and Contact Residues. Spot 1 (S1) , located on the back of the fingers domain between β8 and β12 (Fig. 4 ), exhibits the lowest global occupancy at 2.41% (Table 2 ) and was exclusively observed in system 3 (Table 1 , Table S1). Fragments binding at S1 established contacts with an average of 8 residues (Table S2 and S3). Table 2 Comparison of MDFFr spots by global occupancy (calculated by Eq. 1) Spot # Experimental Global Occupancy (%) S1 no 2.41 S2 no 4.17 S3 no 8.91 S4 no 7.16 S5 no 6.17 S6 no 5.36 S7 no 3.77 S8 no 4.64 S9 no 9.50 S10 no 4.51 S11 Yes (YRL and HBA) 49.60 S12 no 7.47 Spot 2 (S2) is positioned on the back of the palm domain (Fig. 4 ), making contact with the loop between β13 and β14. With a low global occupancy of 4.17% and exclusive presence in system 7, fragments at S2 maintained contacts with an average of 8 residues (Table S2 and S3). Notably, this site represents a previously unidentified allosteric binding location. Spot 3 (S3) is situated at the base of the fingers where it meets the palm domain, contacting residues on β10 and α7. With the third-highest global occupancy at 8.91% (Table 2 ), S3 was uniquely observed in system 9 and established an average of 6 contacts (Table S2 and S3). Significantly, four of its contact residues (ARG166, GLU167, GLU203, and MET208) overlap with the experimentally observed UBS1 site (PDB ID: 6XAA,[ 70 ] UBS1 bound to UB; and PDB ID: 6YVA,[ 12 ] UBS2 bound to ISG15; Table S4), representing the most extensive overlap with these functionally crucial protein-protein interaction sites. Spot 4 (S4) is positioned at the fingers-thumb interface, establishing contacts with β8 and β12 at their respective termini (Fig. 4 ). With a global occupancy of 7.16%, S4 appears in systems 3 and 5, maintaining an average of 4 contacts (Table S2 and S3). Spot 5 (S5) , located on the thumb domain between α3 and α6, shows a global occupancy of 6.17% and appears exclusively in system 8. S5 maintains an average of 5 contacts, with two residues (ASN156, LYS157) overlapping with the experimentally observed UBS1 site (Table S1V). Spot 6 (S6) is positioned behind the thumb domain between α5, α6, and the UBL domain. With a global occupancy of 5.36% and exclusive presence in system 8, S6 establishes an average of 9 contacts (Table S2 and S3). Spot 7 (S7) , situated on the front of α3 in the thumb domain, displays the second-lowest global occupancy at 3.77% and appears only in system 3. Along with S4, it maintains the lowest average number of contacts at 4 (Tables S2). Spot 8 (S8) is located on the thumb domain's posterior surface, just below α2 (Fig. 4 ). With a low global occupancy of 4.64% and exclusive presence in system 8 (Table S1), S8 maintains 6 contacts, including GLU70, which overlaps with the experimentally observed UBS2 site and establishes contacts with ISG15 (6YVA) (Table S2 and S3). Spot 9 (S9) occupies a deep cavity between the thumb and UBL domains, positioned directly below S11, the nearest site to the experimental fragment binding pocket (Figs. 1 and 4 ). With the second-highest global occupancy at 9.50%, S9 appears in systems 5 and 8 and maintains the highest number of contacts at 12 (Table S2 and S3). Spot 10 (S10) is positioned beneath α2, toward the posterior interface between thumb and UBL domains. With a low global occupancy of 4.51%, S10 appears only in system 5 and maintains an average of 5 contacts (Tables S2). Spot 11 (S11) resides in a deep cavity between the thumb and UBL domains. It demonstrates the highest global occupancy by far at 49.60% and appears in systems 1, 3, 6, 7, 8, and 9. S11 lies closest to the experimental YRL and HBA binding site (Fig. 2 A) and maintains the second-highest number of contacts, with 8 residues in each system where it appears (Tables S2). Additional experimental structures with co-crystallized ligands are known.[ 71 ] Spot 12 (S12) is located between α3 helix in the thumb domain and the UBL domain. With the fourth-highest global occupancy at 7.47%, S12 appears only in system 3 and maintains an average of 6 contacts. S12 lies nearest to S11, with the spots nearly overlapping (Fig. 4 and Table S2 and S3). Both spots include contacts with residues found in the experimentally known Ub-S2 site (6YVA and SARS 5E6J, Table S3), with S11 contacting ARG65, PHE69, and PRO77, the latter also being part of S12's contacts (Table S2 and S3). Water Thermodynamics. To characterize the physicochemical properties of the binding spots identified through MDFFr, we performed SZMAP calculations on representative frames closest to the average structure of each spot (Fig. 5 ). SZMAP analyzes water thermodynamics to map binding site polarity, where negative regions (yellow, blue, and green) indicate favorable interactions with polar groups, while positive regions (red) favor non-polar interactions. The analysis of spots S1-S4 revealed distinct polarity patterns. In S1, while the phenol hydroxyl occupies an unfavorable hydrophobic pocket, the hydroxyethyl hydroxyl extends into a hydrophilic region. S2 features a hydrophobic pocket accommodating the fragment's bulk, with its carboxyl group engaging in polar interactions. S3 comprises a large hydrophobic pocket where the presence of oxygen atoms creates a polarity mismatch, though potential expansion into adjacent polar regions appears possible through extended polar groups. In S4, the binding mode combines a hydrophobic pocket housing the aromatic ring with a hydrophilic pocket accommodating the phenol hydroxyl. Spots S5-S6 show complementary physicochemical properties. S5 positions YRL in a hydrophobic pocket flanked by hydrophilic regions near both hydroxyl groups. In S6, UBL domain displacement creates a large hydrophobic pocket where HE9 binds between α5, α6, and the UBL loop, with its alcohol groups extending into solvent-exposed regions and its carbonyl group engaging with polar moieties. The remaining spots (S7-S12) display varied thermodynamic profiles. S7 presents a narrow hydrophobic channel along α3, while S8 forms a weakly hydrophobic, solvent-exposed pocket near α2. S9 features an extensive hydrophobic cavity extending from S11, offering considerable space for fragment expansion. S10, characterized by significant UBL displacement, shares similarities with S6 in its interaction pattern with HBA engaging the loop between β1 and β2, but differs in interacting with α2 rather than α5's bottom region. The highly occupied S11 forms a large hydrophobic pocket between UBL, α2, and α3, with strategic hydrophilic contact points that could guide fragment optimization. S12 presents a compact hydrophobic pocket with a hydrophilic patch accommodating the phenol hydroxyl, and its proximity to S11 suggests potential for fragment linking strategies. Comparison of computational versus experimental binding sites & cryptic pockets Experimental Binding . Our MDFFr simulations revealed significant correlations with experimental binding data. Spot 11 (S11), while positioned slightly above the crystallographic pose, emerged as the closest fragment-binding hot-spot to the experimentally determined positions of YRL and HBA. Similarly, Spot 9 (S9) closely approximates these crystallographic poses, though located just below the experimental binding position. The significance of these sites is reinforced by their occupancy statistics: S11 shows dominant binding with 49.60% global occupancy (Table 2 ), while S9 follows with 9.50%. Their combined occupancy of 59.20% substantially exceeds that of all other identified binding spots, strongly validating their biological relevance. The binding dynamics at S11 showed system-dependent behavior. In systems 1, 3, and 6, fragments exhibited rapid binding to S11 and maintained stable occupancy throughout the simulation (Figure S8). However, several systems showed no fragment binding at S11, primarily due to conformational changes leading to pocket collapse (Figure S9). Specific examples include sys-2, where movements of PRO59, LEU80, and TYR72 caused the pocket's upper region to collapse; sys-4, where similar collapse prevented fragment access; and sys-5, where ubiquitin domain rotation resulted in pocket closure (Figure S5). Regarding HE9 binding, our simulations revealed interesting disparities with experimental data. While none of our identified spots directly corresponded to the crystallographic HE9 pose, which uniquely occupies a solvent-exposed position at the ISG15 interface, our simulations showed HE9 preferentially binding to S11, near the crystallographic poses of YRL and HBA (Figure S10). This finding is particularly noteworthy given that all three compounds (YRL, HBA, and HE9) experimentally bind at the ISG15/Ub-S2 protein-protein interaction site, though HE9 binds on the opposite face of the α2 helix compared to the concave recession where YRL and HBA bind.[ 16 ] Importantly, the experimental HE9 binding location coincides with the outer, solvent-exposed region of a continuous binding surface formed by the merger of S11, S9, and S10. The absence of spots corresponding to the experimentally known GRL binding site likely stems from simulation time limitations. With individual simulations of 300 ns (totaling 2.7 µs), the BL2 loop between β14 and β15 may not have had sufficient time to adopt an open conformation, thus preventing access to the GRL binding pocket. Cryptic Pockets . To investigate potential functional changes in PLpro conformational dynamics revealed by MDFFr, we employed "exposons" calculations.[ 67 ] This analytical method identifies groups of residues that exhibit strong correlations in their solvent accessible surface area changes. These exposons, which include the previously mentioned GRL binding site, typically manifest only during rare conformational events. Our analysis revealed three significant exposons with distinct characteristics. The largest exposon (shown in red in Fig. 6 ) corresponds to the BL2 loop opening that exposes the active site. This finding is particularly significant as it demonstrates that our MDFFr trajectories, despite individual simulation times of 300 ns, can collectively (2.7 µs) capture cooperative residue movements associated with functional events. The second largest exposon (shown in blue in Fig. 6 ) maps to the Ub-S2 binding site at the base of the ISG15 interaction site, overlapping with our identified spots S5, S7, S9, S11, and S12. This exposon correlates with the UBL domain's movement away from the thumb domain, as confirmed by our domain angle analysis (Figs. S4 and S5). Notably, this region encompasses both the experimentally validated binding sites for YRL and HBA, and includes S11 and S9, our highest occupancy binding spots (Table 2 ). The third largest exposon (shown in green in Fig. 6 ) coincides with our S2 spot (Table 2 ), positioned on the thumb domain's posterior surface, beneath the UBL domain (Fig. 4 ). The functional significance of this concerted amino acid movement warrants further investigation, as it may reveal previously unknown allosteric regulatory mechanisms. DISCUSSION AND CONCLUSIONS We have developed a computational framework, called Molecular Dynamics Simulations with Flooding Fragments (MDFFr), to accurately identify druggable hot-spots for fragment-ligand binding based solely on the 3D structure of the protein, without prior knowledge of binding site(s). This physics-based approach, which evolves from previous studies conducted on membrane proteins,[ 44 – 47 ] has been significantly modified and generalized for use with globular proteins, as demonstrated with the viral PLpro enzyme of SARS-CoV2 (our case study). From a biological perspective, MDFFr identified 12 putative ligand binding sites (S1 to S12; Fig. 4 and Table 2 ) on PLpro. Two sites exhibited notably high overall occupancy: S11 (~ 50%) and S9 (~ 10%), with their individual and combined occupancies significantly exceeding all other spots (Table 2 ). Notably, S11 corresponds to the experimentally validated binding site for two fragment-probes (YRL and HBA) used in our MD studies (Figs. 2 A and 4 ; Table 2 ), providing proof-of-concept for MDFFr's application in blind-challenge MD studies (with random fragment placement and no prior binding knowledge). Analysis of the overall site mapping on PLpro (Fig. 4 ) revealed that S11, S9, and S10 merge to form a continuous binding pocket (combined global occupancy ~ 64%) anchored within the thumb domain at the interface of helices α3, α5, and α6 (deeply buried), and the bordering α2 helix (solvent-exposed). While this pocket partially overlaps with the shallow Ub-S2 interaction site, the YRL/HBA binding site, and the thiosemicarbazones derivatives binding site (PDB ID: 7QCI)[ 71 ] (S10-S11), it extends inward to form a previously unobserved buried groove (S9). A particularly significant finding is the S2 fragment-binding spot (third highest occupancy, ~ 4%; Table 2 ), positioned near but not overlapping with the experimentally known Ub-S1 protein-protein interaction site (Figs. 2 C and 4 ). This newly identified site coincides with the third largest cluster of residues showing cooperative changes in solvent accessibility (exposon), suggesting potential roles as a conformational switch, allosteric site, or cryptic pocket.[ 67 ] Interestingly, S2 is proximal to (but not overlapping with) an allosteric site recently mapped on PLpro.[ 72 ] This site was identified through computational structural analyses of PLpro experimental structures and validated through enzyme titration experiments of computationally predicted protein mutants.[ 72 ] While complete validation of S2 as a new allosteric regulatory site requires additional computational and experimental work, these findings support S2's significance and demonstrate MDFFr's utility in physics-based mapping of novel allosteric/functional sites. Methodologically, MDFFr adapts the "flooding MD simulations"[ 44 – 47 ] concept of using known or putative binder molecules to reveal binding sites or access pathways (e.g., cocrystal-fragments to PLpro or volatile anesthetics to ion channels). Unlike the original approach,[ 44 – 47 ] MDFFr offers broader applicability to all proteins, including non-membrane bound ones. Our PLpro application demonstrates MDFFr's potential as a general MD-based fragment mapping approach, suitable for "computational-fragment screening" to complement experimental high-throughput screening and support hit-to-lead optimizations. Current MD-based binding-site mapping approaches[ 26 , 35 – 39 , 42 , 43 , 43 , 73 ] typically employ large-scale simulations with classical or accelerated MD implementations, using cosolvent "probes" to map protein surface features. While effective, these methods face challenges in cosolvent selection and computational costs. In contrast, MDFFr offers straightforward implementation with any forcefield and MD engine, utilizing a "building block" approach with short simulation runs (300 ns per system) and flexible fragment concentrations. Important considerations for MDFFr include: (i) lower cosolvent fragment concentrations compared to membrane-specific "flooding MD simulations",[ 44 – 47 ] and (ii) focus on binding hot spot location prediction rather than binding affinity correlation. Future work will expand MDFFr applications to diverse protein sets with varying fragment types and concentrations, while evaluating its effectiveness in computational fragment screening and protein-protein interaction site identification, including systems with additional biomacromolecules like DNA/RNA. Declarations ACKNOWLEDGMENTS This research was supported in part by The RF Research Foundation and the City University of New York (CUNY) Queens College startup funds (E.G.); the University of Wyoming (UWYO) and School of Pharmacy (SOP) startup funds (K.E.). Computational resources were provided by Temple University's high-performance computing facility, supported in part by the National Science Foundation through major research instrumentation grant number 1625061. The authors gratefully acknowledge OpenEye Scientific Software for providing a free academic license. AUTHORS CONTRIBUTIONS Conceptualization: Jason Pattis, Khaled Elokely, Eleonora Gianti. Methodology: Jason Pattis, Khaled Elokely, Eleonora Gianti. Formal analysis and investigation: Jason Pattis, Khaled Elokely, Eleonora Gianti. Writing: Jason Pattis, Khaled Elokely, Eleonora Gianti. COMPLIANCE WITH ETHICAL STANDARDS Competing Interests : The authors declare no conflict of interest. DATA AVAILABILITY All data produced by this study are disclosed in the publication and its associated supporting information (SI) file. The code for these analyses and input files are available at: https://github.com/jgpattis/MDFFr References Zhou P, Yang X-L, Wang X-G, et al (2020) A pneumonia outbreak associated with a new coronavirus of probable bat origin. Nature 579:270–273. https://doi.org/10.1038/s41586-020-2012-7 Wu F, Zhao S, Yu B, et al (2020) A new coronavirus associated with human respiratory disease in China. Nature 579:265–269. https://doi.org/10.1038/s41586-020-2008-3 Lu R, Zhao X, Li J, et al (2020) Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding. The Lancet 395:565–574. https://doi.org/10.1016/S0140-6736(20)30251-8 Vegivinti CTR, Evanson KW, Lyons H, et al (2022) Efficacy of antiviral therapies for COVID-19: a systematic review of randomized controlled trials. BMC Infect. Dis. 22:1–45 Hwang Y-C, Lu R-M, Su S-C, et al (2022) Monoclonal antibodies for COVID-19 therapy and SARS-CoV-2 detection. J. Biomed. Sci. 29:1–50 Lei J, Kusov Y, Hilgenfeld R (2018) Nsp3 of coronaviruses: Structures and functions of a large multi-domain protein. Antiviral Res 149:58–74. https://doi.org/10.1016/j.antiviral.2017.11.001 Hilgenfeld R (2014) From SARS to MERS: crystallographic studies on coronaviral proteases enable antiviral drug design. FEBS J 281:4085–4096. https://doi.org/10.1111/febs.12936 Narayanan A, Narwal M, Majowicz SA, et al (2022) Identification of SARS-CoV-2 inhibitors targeting Mpro and PLpro using in-cell-protease assay. Commun. Biol. 5:1–17 Barretto N, Jukneliene D, Ratia K, et al (2005) The Papain-Like Protease of Severe Acute Respiratory Syndrome Coronavirus Has Deubiquitinating Activity Lindner HA, Lytvyn V, Qi H, et al (2007) Selectivity in ISG15 and ubiquitin recognition by the SARS coronavirus papain-like protease. Arch Biochem Biophys 466:8–14 Del Valle DM, Kim-Schulze S, Huang H-H, et al (2020) An inflammatory cytokine signature predicts COVID-19 severity and survival. Nat. Med. 26:1636–1643 Shin D, Mukherjee R, Grewe D, et al (2020) Papain-like protease regulates SARS-CoV-2 viral spread and innate immunity. Nature 587:657–662 Rut W, Lv Z, Zmudzinski M, et al (2020) Activity profiling and crystal structures of inhibitor-bound SARS-CoV-2 papain-like protease: A framework for anti-COVID-19 drug design. Sci Adv 6 Harcourt BH, Jukneliene D, Kanjanahaluethai A, et al (2004) Identification of Severe Acute Respiratory Syndrome Coronavirus Replicase Products and Characterization of Papain-Like Protease Activity. J Virol 78:13600–13612. https://doi.org/10.1128/jvi.78.24.13600-13612.2004 Báez-Santos YM, St. John SE, Mesecar AD (2015) The SARS-coronavirus papain-like protease: Structure, function and inhibition by designed antiviral compounds. Antiviral Res 115:21–38. https://doi.org/10.1016/j.antiviral.2014.12.015 Osipiuk J, Azizi SA, Dvorkin S, et al (2021) Structure of papain-like protease from SARS-CoV-2 and its complexes with non-covalent inhibitors. Nat Commun 12:743 Gao X, Qin B, Chen P, et al (2021) Crystal structure of SARS-CoV-2 papain-like protease. Acta Pharm Sin B 11:237–245. https://doi.org/10.1016/j.apsb.2020.08.014 Devaraj SG, Wang N, Chen Z, et al (2007) Regulation of IRF-3-dependent Innate Immunity by the Papain-like Protease Domain of the Severe Acute Respiratory Syndrome Coronavirus *. J Biol Chem 282:32208–32221. https://doi.org/10.1074/jbc.M704870200 Clementz MA, Chen Z, Banach BS, et al (2010) Deubiquitinating and Interferon Antagonism Activities of Coronavirus Papain-Like Proteases. J Virol 84:4619–4629. https://doi.org/10.1128/jvi.02406-09 Frieman M, Ratia K, Johnston RE, et al (2009) Severe Acute Respiratory Syndrome Coronavirus Papain-Like Protease Ubiquitin-Like Domain and Catalytic Domain Regulate Antagonism of IRF3 and NF-κB Signaling. J Virol 83:6689–6705. https://doi.org/10.1128/jvi.02220-08 Shen Z, Ratia K, Cooper L, et al (2021) Design of SARS-CoV-2 PLpro Inhibitors for COVID-19 Antiviral Therapy Leveraging Binding Cooperativity. J Med Chem Ratia K, Pegan S, Takayama J, et al (2008) A noncovalent class of papain-like protease/deubiquitinase inhibitors blocks SARS virus replication. Proc Natl Acad Sci 105:16119–16124. https://doi.org/10.1073/pnas.0805240105 Sanders BC, Pokhrel S, Labbe AD, et al (2023) Potent and selective covalent inhibition of the papain-like protease from SARS-CoV-2. Nat Commun 14:1733. https://doi.org/10.1038/s41467-023-37254-w Patchett S, Lv Z, Rut W, et al (2021) A molecular sensor determines the ubiquitin substrate specificity of SARS-CoV-2 papain-like protease. Cell Rep 36:109754 Srinivasan V, Brognaro H, Prabhu PR, et al (2022) Antiviral activity of natural phenolic compounds in complex at an allosteric site of SARS-CoV-2 papain-like protease. Commun. Biol. 5:805 Brenke R, Kozakov D, Chuang G-Y, et al (2009) Fragment-based identification of druggable ‘hot spots’ of proteins using Fourier domain correlation techniques. Bioinformatics 25:621–627. https://doi.org/10.1093/bioinformatics/btp036 Erlanson DA, McDowell RS, O’Brien T (2004) Fragment-Based Drug Discovery. J Med Chem 47:3463–3482. https://doi.org/10.1021/jm040031v Erlanson DA, Jahnke W (2015) Fragment-based Drug Discovery: Lessons and Outlook. John Wiley & Sons Bottegoni G, Favia AD, Recanatini M, Cavalli A (2012) The role of fragment-based and computational methods in polypharmacology. Drug Discov Today 17:23–34. https://doi.org/10.1016/j.drudis.2011.08.002 Zauhar RJ, Moyna G, Tian L, et al (2003) Shape Signatures: A New Approach to Computer-Aided Ligand- and Receptor-Based Drug Design. J Med Chem 46:5674–5690. https://doi.org/10.1021/jm030242k Gianti E, Sartori L (2008) Identification and Selection of “Privileged Fragments” Suitable for Primary Screening. J Chem Inf Model 48:2129–2139. https://doi.org/10.1021/ci800219h Gianti E, Messick TE, Lieberman PM, Zauhar RJ (2016) Computational analysis of EBNA1 “druggability” suggests novel insights for Epstein-Barr virus inhibitor design. J Comput Aided Mol Des 30:285–303. https://doi.org/10.1007/s10822-016-9899-y Bissaro M, Sturlese M, Moro S (2020) The rise of molecular simulations in fragment-based drug design (FBDD): an overview. Drug Discov Today 25:1693–1701. https://doi.org/10.1016/j.drudis.2020.06.023 de Souza Neto LR, Moreira-Filho JT, Neves BJ, et al (2020) In silico Strategies to Support Fragment-to-Lead Optimization in Drug Discovery. Front Chem 8:. https://doi.org/10.3389/fchem.2020.00093 Yanagisawa K, Moriwaki Y, Terada T, Shimizu K (2021) EXPRORER: Rational Cosolvent Set Construction Method for Cosolvent Molecular Dynamics Using Large-Scale Computation. J Chem Inf Model 61:2744–2753. https://doi.org/10.1021/acs.jcim.1c00134 Bakan A, Nevins N, Lakdawala AS, Bahar I (2012) Druggability Assessment of Allosteric Proteins by Dynamics Simulations in the Presence of Probe Molecules. J Chem Theory Comput 8:2435–2447. https://doi.org/10.1021/ct300117j Zariquiey FS, Souza JV de, Bronowska AK (2019) Cosolvent Analysis Toolkit (CAT): a robust hotspot identification platform for cosolvent simulations of proteins to expand the druggable proteome. Sci Rep 9:19118. https://doi.org/10.1038/s41598-019-55394-2 Raman EP, Yu W, Lakkaraju SK, MacKerell ADJr (2013) Inclusion of Multiple Fragment Types in the Site Identification by Ligand Competitive Saturation (SILCS) Approach. J Chem Inf Model 53:3384–3398. https://doi.org/10.1021/ci4005628 Ghanakota P, Carlson HA (2016) Moving Beyond Active-Site Detection: MixMD Applied to Allosteric Systems. J Phys Chem B 120:8685–8695. https://doi.org/10.1021/acs.jpcb.6b03515 Ghanakota P, Carlson HA (2016) Driving Structure-Based Drug Discovery through Cosolvent Molecular Dynamics. J Med Chem 59:10383–10399. https://doi.org/10.1021/acs.jmedchem.6b00399 Ghanakota P, DasGupta D, Carlson HA (2019) Free Energies and Entropies of Binding Sites Identified by MixMD Cosolvent Simulations. J Chem Inf Model 59:2035–2045. https://doi.org/10.1021/acs.jcim.8b00925 Smith RD, Carlson HA (2021) Identification of Cryptic Binding Sites Using MixMD with Standard and Accelerated Molecular Dynamics. J Chem Inf Model 61:1287–1299. https://doi.org/10.1021/acs.jcim.0c01002 Graham SE, Leja N, Carlson HA (2018) MixMD Probeview: Robust Binding Site Prediction from Cosolvent Simulations. J Chem Inf Model 58:1426. https://doi.org/10.1021/acs.jcim.8b00265 Vemparala S, Domene C, Klein ML (2008) Interaction of Anesthetics with Open and Closed Conformations of a Potassium Channel Studied via Molecular Dynamics and Normal Mode Analysis. Biophys J 94:4260–4269. https://doi.org/10.1529/biophysj.107.119958 Raju SG, Barber AF, LeBard DN, et al (2013) Exploring Volatile General Anesthetic Binding to a Closed Membrane-Bound Bacterial Voltage-Gated Sodium Channel via Computation. PLoS Comput Biol 9:e1003090. https://doi.org/10.1371/journal.pcbi.1003090 Yang E, Granata D, Eckenhoff RG, et al (2018) Propofol inhibits prokaryotic voltage-gated Na+ channels by promoting activation-coupled inactivation. J Gen Physiol 150:1299–1316. https://doi.org/10.1085/jgp.201711924 Barber AF, Carnevale V, Klein ML, et al (2014) Modulation of a voltage-gated Na+ channel by sevoflurane involves multiple sites and distinct mechanisms. Proc Natl Acad Sci 111:6726–6731. https://doi.org/10.1073/pnas.1405768111 Dolinsky TJ, Czodrowski P, Li H, et al (2007) PDB2PQR: expanding and upgrading automated preparation of biomolecular structures for molecular simulations. Nucleic Acids Res 35:W522–W525. https://doi.org/10.1093/nar/gkm276 Abraham MJ, Murtola T, Schulz R, et al (2015) GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 1–2:19–25. https://doi.org/10.1016/j.softx.2015.06.001 Jorgensen WL, Chandrasekhar J, Madura JD, et al (1983) Comparison of simple potential functions for simulating liquid water. J Chem Phys 79:926. https://doi.org/10.1063/1.445869 Parrinello M, Rahman A (1981) Polymorphic transitions in single crystals: A new molecular dynamics method. J Appl Phys 52:7182–7190. https://doi.org/10.1063/1.328693 Huang J, Rauscher S, Nawrocki G, et al (2017) CHARMM36m: an improved force field for folded and intrinsically disordered proteins. Nat Methods 14:71–73. https://doi.org/10.1038/nmeth.4067 Vanommeslaeghe K, Hatcher E, Acharya C, et al (2009) CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J Comput Chem NA-NA. https://doi.org/10.1002/jcc.21367 Vanommeslaeghe K, MacKerell ADJr (2012) Automation of the CHARMM General Force Field (CGenFF) I: Bond Perception and Atom Typing. J Chem Inf Model 52:3144–3154. https://doi.org/10.1021/ci300363c Vanommeslaeghe K, Raman EP, MacKerell ADJr (2012) Automation of the CHARMM General Force Field (CGenFF) II: Assignment of Bonded Parameters and Partial Atomic Charges. J Chem Inf Model 52:3155–3168. https://doi.org/10.1021/ci3003649 McGibbon RT, Beauchamp KA, Harrigan MP, et al (2015) MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories. Biophys J 109:1528–32 Hunter JD (2007) Matplotlib: A 2D Graphics Environment. Comput. Sci. Eng. 9:90–95 Scherer MK, Trendelkamp-Schroer B, Paul F, et al (2015) PyEMMA 2: A Software Package for Estimation, Validation, and Analysis of Markov Models. J Chem Theory Comput 11:5525–42 Ester M, Kriegel H-P, Xu X A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise Schubert E, Sander J, Ester M, et al (2017) DBSCAN Revisited, Revisited: Why and How You Should (Still) Use DBSCAN. ACM Trans Database Syst 42:19:1-19:21. https://doi.org/10.1145/3068335 Pedregosa F, Varoquaux G, Gramfort A, et al (2011) Scikit-learn: Machine Learning in Python. J Mach Learn Res 12:2825–2830 Humphrey W, Dalke A, Schulten K (1996) VMD: visual molecular dynamics. J Mol Graph 14:33–38, 27–28 Schmidtke P, Bidon-Chanal A, Luque FJ, Barril X (2011) MDpocket: open-source cavity detection and characterization on molecular dynamics trajectories. Bioinformatics 27:3276–3285. https://doi.org/10.1093/bioinformatics/btr550 Tanger JCI, Pitzer KS (1989) Calculation of the thermodynamic properties of aqueous electrolytes to 1000.degree.C and 5000 bar from a semicontinuum model for ion hydration. J Phys Chem 93:4941–4951. https://doi.org/10.1021/j100349a053 Rashin AA, Bukatin MA (1991) Continuum based calculations of hydration entropies and the hydrophobic effect. J Phys Chem 95:2942–2944. https://doi.org/10.1021/j100161a002 Grant JA, Pickup BT, Nicholls A (2001) A smooth permittivity function for Poisson–Boltzmann solvation methods. J Comput Chem 22:608–640. https://doi.org/10.1002/jcc.1032 Porter JR, Moeder KE, Sibbald CA, et al (2019) Cooperative Changes in Solvent Exposure Identify Cryptic Pockets, Switches, and Allosteric Coupling. Biophys. J. 116:818–830 Frey BJ, Dueck D (2007) Clustering by Passing Messages Between Data Points. Science 315:972–976. https://doi.org/10.1126/science.1136800 Herold J, Siddell SG, Gorbalenya AE (1999) A Human RNA Viral Cysteine Proteinase That Depends upon a Unique Zn2+-binding Finger Connecting the Two Domains of a Papain-like Fold *. J Biol Chem 274:14918–14925. https://doi.org/10.1074/jbc.274.21.14918 Klemm T, Ebert G, Calleja DJ, et al (2020) Mechanism and inhibition of the papain‐like protease, PLpro, of SARS‐CoV‐2. EMBO J 39:e106275. https://doi.org/10.15252/embj.2020106275 Ewert W, Günther S, Miglioli F, et al (2022) Hydrazones and Thiosemicarbazones Targeting Protein-Protein-Interactions of SARS-CoV-2 Papain-like Protease. Front Chem 10:. https://doi.org/10.3389/fchem.2022.832431 Ferreira JC, Villanueva AJ, Adem KA, et al (2024) Identification of novel allosteric sites of SARS-CoV-2 papain-like protease (PLpro) for the development of COVID-19 antivirals. J Biol Chem 300:. https://doi.org/10.1016/j.jbc.2024.107821 Vithani N, Zhang S, Thompson JP, et al (2024) Exploration of Cryptic Pockets Using Enhanced Sampling Along Normal Modes: A Case Study of KRAS G12D. J Chem Inf Model. https://doi.org/10.1021/acs.jcim.4c01435 Supplemental Figures Supplemental Figures S1-S5 are not available with this version. Additional Declarations No competing interests reported. Cite Share Download PDF Status: Published Journal Publication published 19 Dec, 2025 Read the published version in Journal of Computer-Aided Molecular Design → Version 1 posted Editorial decision: Revision requested 05 Oct, 2025 Reviews received at journal 05 Oct, 2025 Reviews received at journal 22 Sep, 2025 Reviewers agreed at journal 18 Sep, 2025 Reviewers agreed at journal 14 Sep, 2025 Reviewers invited by journal 12 Sep, 2025 Editor assigned by journal 12 Sep, 2025 Submission checks completed at journal 02 Sep, 2025 First submitted to journal 31 Aug, 2025 You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-7499914","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":509156419,"identity":"9673b75f-743b-4e31-baf9-93c9767fc75d","order_by":0,"name":"Jason Pattis","email":"","orcid":"","institution":"Temple University","correspondingAuthor":false,"prefix":"","firstName":"Jason","middleName":"","lastName":"Pattis","suffix":""},{"id":509156420,"identity":"0462ac83-ce53-4dba-b7f8-6aeadcf9c723","order_by":1,"name":"Khaled Elokely","email":"","orcid":"","institution":"University of Wyoming","correspondingAuthor":false,"prefix":"","firstName":"Khaled","middleName":"","lastName":"Elokely","suffix":""},{"id":509156421,"identity":"0cfc9c66-f249-438a-a8a0-9dc612945423","order_by":2,"name":"Eleonora Gianti","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABGklEQVRIiWNgGAWjYBAC+3YYi/1gA5CUkEOSTMCqxYAZxuJJBGpJkDAmQYsERBFIIyEtzAcfF1Tck2eQYG78XPjDIr3v+AHGz7w7DjPws+cYYPULM1uy8YwzxYYNEozN0jMSJHJnnklgluY9c5hBsucNVi0GzDxm0rxtCYwNEg8bpHmAWjbcYGBj5m07zGBwA7stBsz830Ba7BskEpt/A7WkG8C02OPUwsMG0pII1NIGsiUBrsVAApcWNmNjnjMJyW08iW3WPGkShjPPJDZLzm1L55E486wAq/fbmx8+5qlIsO1nP/74No9NnTzf8cMHP7xts5bjb0/egDWUYYANzjrA2AAkm3nwKkcFB8BkHQk6RsEoGAWjYJgDAOnOVZyuqYCEAAAAAElFTkSuQmCC","orcid":"","institution":"Queens College, CUNY","correspondingAuthor":true,"prefix":"","firstName":"Eleonora","middleName":"","lastName":"Gianti","suffix":""}],"badges":[],"createdAt":"2025-08-31 10:23:10","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-7499914/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-7499914/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1007/s10822-025-00730-0","type":"published","date":"2025-12-19T15:56:56+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":91375296,"identity":"d5e0dae5-7051-49b1-8887-d64537089be4","added_by":"auto","created_at":"2025-09-15 19:55:40","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":936405,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eDomain organization and structural features of SARS-CoV-2 PLpro. \u003c/strong\u003eThe protein consists of four distinct domains: the UBL domain (residues 1-62, yellow), the thumb domain (residues 63-182, green), the fingers domain (residues 183-241, cyan), and the palm domain (residues 242-315, orange). The thumb, fingers, and palm domains together form the ubiquitin-specific proteases (USP) domain. The catalytic triad (C111, H272, and D286) is displayed in licorice representation with carbon atoms in grey and oxygen atoms in red. The protein backbone is shown in new cartoon representation. Inset: Close-up view of the catalytic active site showing the spatial arrangement of the catalytic residues. Abbreviations: UBL, ubiquitin-like domain; USP, ubiquitin-specific proteases domain.\u003c/p\u003e","description":"","filename":"floatimage1.png","url":"https://assets-eu.researchsquare.com/files/rs-7499914/v1/8623cc07425b8573f4635001.png"},{"id":91375307,"identity":"e281b485-f1bb-40bf-92b5-f9a13995058d","added_by":"auto","created_at":"2025-09-15 19:55:41","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":459496,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eStructural characterization of PLpro binding interactions. \u003c/strong\u003e(A) Crystal structures of PLpro in complex with three phenolic inhibitors: YRL (purple), HBA (pink), and HE9 (grey), shown in licorice representation. (B) Complex structure of PLpro with ISG15 (grey), highlighting the protein-protein interaction interface. (C) PLpro-ubiquitin complex showing two distinct ubiquitin binding sites: Ub-S1 (magenta) and Ub-S2 (silver). (D) Chemical structures of the three experimental ligands (YRL, HBA, HE9) used as flooding fragments in MDFFr simulations. In panels A-C, PLpro is displayed in new cartoon representation with domain coloring consistent with Figure 1 (UBL in yellow, thumb in green, fingers in cyan, and palm in orange).\u003c/p\u003e","description":"","filename":"floatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-7499914/v1/7e16e5ac7448a297cb7289dd.png"},{"id":91375294,"identity":"ca09db88-b44f-4c62-83da-3e0d5369b672","added_by":"auto","created_at":"2025-09-15 19:55:40","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":994167,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003ePrincipal Component Analysis (PCA) of PLpro global motions. \u003c/strong\u003e(A) Structural representation of the extreme conformations (minimum and maximum frames) for the first two principal components. Left panels show PC1, which captures the forward movement of the UBL domain toward the catalytic site. Right panels show PC2, which describes the combined downward motion of the UBL domain and helix-3. (B) Conformational sampling density plots for all nine systems projected onto the PC1/PC2 subspace. The analysis was performed using PyEmma on the Cartesian coordinates of protein backbone atoms across all trajectories. The observed motions show no clear dependence on fragment concentration or type. Additional motion analysis is provided in Supplemental Figure 3.\u003c/p\u003e","description":"","filename":"floatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-7499914/v1/8164218bc16de026d0406d5e.png"},{"id":91375311,"identity":"d66f009c-7da4-45f8-a49e-4346c408d13e","added_by":"auto","created_at":"2025-09-15 19:55:41","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":986832,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eIdentification of fragment binding sites on PLpro through MDFFr analysis. \u003c/strong\u003eThe twelve identified binding spots (S1-S12) are represented as red spheres, with each sphere marking the center of mass of clustered fragment positions throughout the simulations. Main panels show front and back views of PLpro, with insets providing rotated perspectives to highlight binding spots that are partially obscured in the main views. The protein structure is displayed using the same domain coloring scheme as in Figure 1 (UBL in yellow, thumb in green, fingers in cyan, and palm in orange). The binding spots were identified through clustering analysis of fragment positions in Cartesian coordinate space across all simulation trajectories.\u003c/p\u003e","description":"","filename":"floatimage4.png","url":"https://assets-eu.researchsquare.com/files/rs-7499914/v1/d0a0dd2b39384b450712aa35.png"},{"id":91375309,"identity":"e702a703-fc03-4e5d-88b8-cc7cc8c8b899","added_by":"auto","created_at":"2025-09-15 19:55:41","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":1646857,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eWater thermodynamics analysis of PLpro binding spots. \u003c/strong\u003eSZMAP calculations for all twelve binding spots (S1-S12) showing the energetic preferences for water molecules in each APO binding site. The thermodynamic mapping is represented by colored spheres: red spheres indicate regions of positive (unfavorable) free energy where water displacement would be energetically favorable, while blue spheres indicate regions of negative (favorable) free energy where water molecules are energetically stabilized. This analysis reveals the physicochemical nature of each binding site and potential opportunities for fragment optimization through strategic placement of polar and nonpolar groups.\u003c/p\u003e","description":"","filename":"floatimage5.png","url":"https://assets-eu.researchsquare.com/files/rs-7499914/v1/e3c6f10913f4771acac3869e.png"},{"id":91375299,"identity":"61ce410f-6b3d-4834-8575-95753c2443e8","added_by":"auto","created_at":"2025-09-15 19:55:40","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":1247593,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eAnalysis of cooperative solvent exposure changes in PLpro through exposon identification. \u003c/strong\u003eThe three largest exposons identified during MDFFr simulations are visualized: first exposon (red, BL2 loop region), second exposon (blue, Ub-S2 binding site), and third exposon (green, S2 spot region). (A) Front view and (B) back view of PLpro. The protein backbone is displayed in new cartoon representation, while residue side chains within each exposon are shown in licorice representation. These exposons represent groups of residues that undergo coordinated changes in their solvent accessibility, potentially indicating functionally relevant conformational changes or cryptic binding sites.\u003c/p\u003e","description":"","filename":"floatimage6.png","url":"https://assets-eu.researchsquare.com/files/rs-7499914/v1/b5c2617f33f33e9cb6d3bf8b.png"},{"id":98813823,"identity":"9da7e3ae-498c-4fcc-b85a-a60ffe5880af","added_by":"auto","created_at":"2025-12-22 16:02:04","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":7119045,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7499914/v1/c9ffe48f-6133-49f5-9547-7c0aaca25da8.pdf"}],"financialInterests":"No competing interests reported.","formattedTitle":"Chasing allosteric inhibition of the SARS-CoV-2 PLpro via Molecular Dynamics Simulations with Flooding Fragments (MDFFr)","fulltext":[{"header":"INTRODUCTION","content":"\u003cp\u003eThe SARS-CoV-2 virus, a member of the beta Coronavirus (CoV) family,[\u003cspan additionalcitationids=\"CR2\" citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e] emerged in 2019 and precipitated the COVID-19 pandemic, one of the most devastating global health crises in recent history. With over 780M documented cases and \u0026gt;\u0026thinsp;7.1M deaths worldwide (as of August 2025), its impact has been profound and far-reaching. Despite the World Health Organization's declaration ending the international public health emergency in May 2023, SARS-CoV-2 continues to evolve, generating variants of interest (VOIs) such as JN.1 and variants under monitoring (VUMs) including XFG, NB.1.8.1, and others that pose ongoing risks for new epidemic waves, as seen in the summer 2025 surges reported across the United States, Europe, and other regions. While current countermeasures include variant-targeted vaccines and FDA-approved therapeutics (both antiviral drugs and monoclonal antibodies),[\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e, \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e] the persistent threat of SARS-CoV-2 to public health systems necessitates the development of novel therapeutic strategies. Particularly crucial is the need for broad-spectrum interventions that can effectively target multiple circulating variants and mitigate the risk of future outbreaks. Among potential therapeutic targets against SARS-CoV-2, the papain-like protease (PLpro) emerges as a particularly promising candidate (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e1\u003c/span\u003e), given its critical role in modulating host antiviral immune responses against emerging viral variants. Encoded by the nsp3 gene,[\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e] PLpro functions alongside the main protease (3CLpro/Mpro, encoded by nsp5) as one of two essential and highly conserved coronavirus cysteine-proteases. PLpro's primary function involves the processing of viral polyproteins (pp1a and pp1ab), which are synthesized from viral RNA by host ribosomes.[\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e] This crucial proteolytic processing generates 16 mature nonstructural proteins (nsps) that subsequently form the viral replication transcription complex (RTC). Given these vital functions, inhibition of PLpro's catalytic activity represents a compelling strategy for antiviral drug development.[\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e]\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eBeyond its proteolytic function, PLpro plays a crucial role in suppressing host innate immunity through its deubiquitinase (DUB) activity. This immunomodulatory function is achieved by cleaving ubiquitin (Ub) and ubiquitin-like (UbL) proteins, particularly interferon stimulating gene product 15 (ISG15), from host proteins.[\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e, \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e] Specifically, PLpro recognizes and cleaves the C-terminal RLRGG motif of Ub/UbLs, thereby removing these modifications from cellular proteins. This DUB activity significantly impacts host immune responses by dysregulating both inflammatory and interferon-mediated pathways, a mechanism that has been directly linked to increased SARS-CoV-2 mortality rates.[\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e, \u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e]\u003c/p\u003e\u003cp\u003eSARS-CoV-2 PLpro (SCoV2-PLpro) exhibits remarkable conservation among coronaviruses, sharing 83% sequence identity with its SARS counterpart (SCoV1-PLpro) and maintaining 100% conservation of catalytic-site residues crucial for Ub interactions.[\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e] A key distinguishing feature of SCoV2-PLpro is its preferential activity toward ISG15, particularly in cleaving ISG15 from interferon responsive factor 3 (IRF3), whereas SCoV1-PLpro primarily targets ubiquitin chains. This specificity directly impacts type I interferon responses, making SCoV2-PLpro an attractive therapeutic target with potential dual benefits: inhibition of viral infection and enhancement of antiviral immunity.[\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e] Structural studies have revealed that CoV PLpro comprises four distinct domains: a ubiquitin-like domain (UBL) and three domains (thumb, fingers, and palm) that together form the ubiquitin-specific proteases (USP) domain (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e1\u003c/span\u003e). The enzyme's active site specifically recognizes the LXGG/XX substrate motif at the nsp1/2, nsp2/3, and nsp3/4 junctions,[\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e, \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e, \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e] utilizing a catalytic triad of Cys111, His272, and Asp286.[\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e, \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e] This site catalyzes two crucial reactions: (1) cleavage of viral polyprotein after the P1 glycine, releasing nsp1, nsp2, and nsp3 for viral replication, and (2) hydrolysis of isopeptide bonds following the LXGG motif at the C-terminus of cellular Ub and ISG15,[\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e, \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e] thereby modulating innate immune responses during infection and pathogenesis.[\u003cspan additionalcitationids=\"CR19\" citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e] The protein features two distinct ubiquitin-binding sites: the S1 site (Ub-S1), positioned between the fingers and palm domains, which accommodates ubiquitin and the N-terminal UbL domain of ISG15 (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e2\u003c/span\u003eC), and the S2 site (Ub-S2), located at the α2 helix-thumb domain interface, which binds the C-terminal UbL domain of ISG15 (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e2\u003c/span\u003eB).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eThe development of effective PLpro inhibitors has seen significant progress, though the number of validated compounds remains limited.[\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e] The first notable inhibitor was GRL0617, a potent non-covalent naphthalene-based compound that targets the active site.[\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e] X-ray crystallography revealed its unique mechanism of action, where inhibition occurs through induced loop closure of the active site.[\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e] A key challenge in PLpro inhibitor development has been the limited druggability of the 'site I' and 'site II' sub-pockets, primarily due to their featureless Gly-Gly recognition pattern (LXGG substrate).[\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e] Recent advances have addressed this challenge through multiple strategies: (1) leveraging binding cooperativity across multiple sites, leading to 2-phenylthiophenes with nanomolar potency,[\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e] (2) development of selective covalent inhibitors,[\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e] and (3) design of novel compounds targeting the ubiquitin-binding sites. Recent breakthroughs include the discovery of orally bioavailable inhibitors with efficacy against multiple SARS-CoV-2 variants and related coronaviruses, demonstrating PLpro's continued promise as a therapeutic target.\u003c/p\u003e\u003cp\u003eAs of August 2025, PLpro inhibitors for SARS-CoV-2 are still in preclinical or early development stages and lack clinical approval. Approved outpatient therapies mainly target alternative viral proteins, such as Mpro with nirmatrelvir-ritonavir, or employ nucleoside analogs. These include nirmatrelvir-ritonavir (Paxlovid), an oral antiviral for adults and children aged\u0026thinsp;\u0026ge;\u0026thinsp;12 years (\u0026ge;\u0026thinsp;40 kg), initiated within 5 days of symptom onset; remdesivir (Veklury), an intravenous antiviral for adults and children aged\u0026thinsp;\u0026ge;\u0026thinsp;28 days (\u0026ge;\u0026thinsp;3 kg), given over 3 days and started within 7 days of symptoms; molnupiravir (Lagevrio), an oral antiviral for adults, begun within 5 days of symptoms (with potential expiration of distributed supplies); and convalescent plasma for immunocompromised patients under emergency use authorization. Symptomatic supportive care, such as antipyretics and analgesics, is additionally advised.\u003c/p\u003e\u003cp\u003eAn alternative approach to address the limited druggability of sub-pockets involves targeting allosteric sites through multiple mechanisms, including disruption of protein-protein interactions and interference with Ub/UbL substrate binding sites.[\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e, \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e] This strategy has yielded promising results, as demonstrated by the discovery of two fragment-like natural compounds through X-ray crystallography-based high-throughput screening: 4-(2-hydroxyethyl)phenol (YRL) and 4-hydroxybenzaldehyde (HBA). These compounds bind at the interface between the thumb and UBL domains (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e2\u003c/span\u003eA), where they are thought to disrupt the S2 ubiquitin binding site (Ub-S2).[\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e] Additionally, a third compound, 3,4-dihydroxybenzoate (HE9), was co-crystallized in a distinct location exposed to solvent at the back of the α2-helix (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e2\u003c/span\u003eA). Both HE9 and YRL demonstrate significant inhibitory activity against PLpro, with EC50 values of 0.13 \u0026micro;M and 10 \u0026micro;M respectively, establishing them as valuable starting points for further drug development.[\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e]\u003c/p\u003e\u003cp\u003eFragment-based drug discovery has emerged as a powerful approach for identifying novel therapeutic compounds through high-throughput screening of fragment libraries against biological targets. This method typically employs biophysical techniques, including differential scanning fluorimetry and surface plasmon resonance, followed by detailed structural characterization using NMR spectroscopy or X-ray crystallography. The identified fragments are then optimized through growing, linking, or merging strategies. However, experimental characterization of protein-fragment interactions presents significant challenges and limitations. To address these constraints, numerous computational approaches,[\u003cspan additionalcitationids=\"CR27 CR28\" citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e] including our previous work,[\u003cspan additionalcitationids=\"CR31\" citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e] have been developed to enhance various aspects of fragment-screening campaigns. While current computational methods, primarily based on molecular docking and virtual screening, excel at identifying active scaffolds and chemical leads, they have important limitations. These approaches often provide insufficient insight into two critical aspects: (1) the dynamic effects of ligand binding on functional or allosteric protein pathways, and (2) the identification of druggable pockets on protein surfaces. These limitations can only be effectively addressed through molecular simulations,[\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e, \u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e] which provide crucial dynamic information about protein-ligand interactions.\u003c/p\u003e\u003cp\u003eIn this study, we developed and implemented Molecular Dynamics simulations with Flooding Fragments (MDFFr), a computational protocol applied to SCoV2-PLpro over a multi-microsecond timescale (2.7 \u0026micro;s total). This approach serves two primary purposes: (1) computational characterization of fragment ligand binding at experimentally validated locations, and (2) discovery of novel allosteric binding sites. Our method builds upon existing protocols that incorporate protein flexibility in hot-spot mapping, such as those using 'standard probes' (isopropanol, phenol, benzene, etc.)[\u003cspan additionalcitationids=\"CR36 CR37\" citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e] and the MixMD technique, which analyzes cosolvent density distributions across MD trajectories.[\u003cspan additionalcitationids=\"CR40 CR41 CR42\" citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e] We also drew inspiration from \"MD flooding\" methodology, previously employed to study volatile anesthetic binding pathways in membrane proteins.[\u003cspan additionalcitationids=\"CR45 CR46\" citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e] MDFFr represents an advancement of these approaches, offering a generalized MD flooding protocol applicable to both globular proteins and enzymes. Using this protocol, we conducted multiple unbiased simulations of SCoV2-PLpro using high concentrations of experimentally validated co-crystal fragment ligands (\"flooding fragments\") to accelerate binding events and enhance ligand sampling efficiency.[\u003cspan additionalcitationids=\"CR45 CR46\" citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e] Our analysis revealed several key findings:\u003c/p\u003e\u003cp\u003e\u003col\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eMultiple sites suitable for allosteric ligand binding exist on SCoV2-PLpro, including both experimentally confirmed sites[\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e] and previously unidentified druggable spots.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eThe flooding fragments showed preferential binding to locations that align with experimental screening results.[\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e]\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eAdditional \"ligand hot spots\" were identified that correspond to ubiquitin and ISG-15 interaction sites.[\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e, \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e]\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003c/ol\u003e\u003c/p\u003e\u003cp\u003eThese results not only unveil new opportunities for allosteric inhibition of SCoV2-PLpro\u0026mdash;a crucial strategy in combating COVID-19\u0026mdash;but also establish MDFFr as a powerful tool for computational fragment-based screening and binding site identification.\u003c/p\u003e"},{"header":"MATERIALS AND METHODS","content":"\u003cp\u003e\u003cb\u003eStructural Models\u003c/b\u003e. Three high resolution crystal structures of SARS-CoV-2 PLpro 7OFS (1.9 \u0026Aring;), 7OFT (1.95 \u0026Aring;), and 7OFU (1.72 \u0026Aring;),[\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e] co-crystallized with phenolic fragment-like compound inhibitors YRL (4-(2-hydroxyethyl)phenol), HBA (p-hydroxybenzaldehyde), and HE9 (methyl 3,4-bis(oxidanyl)benzoate), respectively, were prepared for unbiased, all-atom molecular dynamics (MD) simulations. Co-crystallized water molecules and fragment-ligands were removed, before the protein structures were protonated using the pdb2pqr program.[\u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e] Cystines 189, 192, 224, and 226 were left deprotonated because they coordinate a zinc molecule.\u003c/p\u003e\u003cp\u003eEach protein structure was then placed in a cubic box with a 10 \u0026Aring; buffer. The GROMACS[\u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e] insert-molecule tool was used to place 10 millimolar (mM), 50 mM, or 100 mM concentration of non-overlapping, individual fragment ligands in the box. (The \u0026ldquo;insert-molecule\u0026rdquo; tool places a molecule randomly in the box but only accepts the placement if there is no overlap with other molecules.) Each box was solvated with TIP3P[\u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e50\u003c/span\u003e] explicit solvent before adding 150 mM of NaCl. A total of 9 systems were prepared, leading to a cubic periodic simulation cell of \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:113\\times\\:\\:113\\times\\:\\:113\\:\\)\u003c/span\u003e\u003c/span\u003e\u0026Aring; for a total of 144,343 atoms for 7OFS, 10 mM; \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:113\\times\\:\\:113\\times\\:\\:113\\)\u003c/span\u003e\u003c/span\u003e \u0026Aring; with a total of 144,244 atoms for 7OFS, 50 mM; \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:113\\times\\:\\:113\\times\\:\\:113\\:\\)\u003c/span\u003e\u003c/span\u003e\u0026Aring; with a total of 144,101 atoms for 7OFS, 100 mM; 111 x 111 x 111 \u0026Aring; with a total of 137,557 atoms for 7OFT, 10 mM; 111 x 111 x 111 \u0026Aring; with a total of 137,476 atoms for 7OFT, 50 mM; 111 x 111 x 111 \u0026Aring; with a total of 137,269 atoms for 7OFT, 100 mM; 111 x 111 x 111 \u0026Aring; with a total of 137,306 atoms for 7OFU, 10 mM; 111 x 111 x 111 \u0026Aring; with a total of 137,182 atoms for 7OFU, 50 mM; and 111 x 111 x 111 \u0026Aring; with a total of 137,022 atoms for 7OFU, 100 mM.\u003c/p\u003e\u003cp\u003e\u003cb\u003eMolecular Dynamics Simulations.\u003c/b\u003e Systems were minimized with 5000 steeps of steepest decent energy minimization. Position restraints of 1000 KJ/mol were placed on each protein and fragment-ligand heavy atoms before conducting 0.5 ns of restrained NPT equilibration runs at 1 bar with the Parrinello-Rahman barostat[\u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e51\u003c/span\u003e] and 300 K with the Nose-Hoover thermostat. One additional run of equilibration for 0.5 ns under restrained NVT conditions was performed. Simulations were carried out using GROMACS[\u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e] version 2021.2. The Charmm36m[\u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e52\u003c/span\u003e] force field was used to simulate all systems with the CHARMM generalizable force field (CGenFF)[\u003cspan additionalcitationids=\"CR54\" citationid=\"CR53\" class=\"CitationRef\"\u003e53\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e55\u003c/span\u003e] used to parameterize the fragment-ligands. Production runs were carried out in the NVT ensemble at 300 K, reaching 0.3 \u0026micro;s for each system, for a total of 2.7 \u0026micro;s of simulation (i.e., 0.3 \u0026micro;s \u0026times; 9 systems). Coordinates of the systems were collected every 10 ps for a total of 30,000 frames for each run.\u003c/p\u003e\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e\u003ch2\u003eAnalysis of the Results\u003c/h2\u003e\u003cp\u003e\u003cb\u003eRMSD and RMSF\u003c/b\u003e. The Python library MDTraj,[\u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e56\u003c/span\u003e] version 1.9.3, was used to compute the root-mean-square deviation (RMSD) of protein and ligand atoms relative to their respective reference experimental structure over the simulation time.[\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e] Similarly, root-mean-square fluctuation (RMSF) of protein and ligand atoms were calculated using GROMACS.[\u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e] All plots were generated using the Python library Matplotlib,[\u003cspan citationid=\"CR57\" class=\"CitationRef\"\u003e57\u003c/span\u003e] version 3.5.0.\u003c/p\u003e\u003cp\u003e\u003cb\u003ePCA\u003c/b\u003e. The essential motions of the simulated systems were captured by the global principal component analysis (PCA), performed using PyEMMA,[\u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e58\u003c/span\u003e] version 2.5.7, on the Cartesian coordinates of protein backbone atoms. PCA was fit on all nine trajectories to capture the highest variance motions of any system. In PCA plots, clusters of samples were color coded by the density of frames. A local PCA of the YRL binding site was performed (using PyEMMA)[\u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e58\u003c/span\u003e] by first selecting all alpha- and beta-carbons located within 9 \u0026Aring; of the YRL-ligand as in the experimental structure. Pairwise distances between these alpha- and beta-carbons were measured in all 3 systems with YRL (i.e., the protein extracted from PDB-ID 7OFS and YRL ligands at 10-, 50- and 100-mM concentrations). Then, PCA was performed on these distances, extracting the largest variance motions. Principal components are the eigenvectors of the data\u0026rsquo;s covariance matrix.\u003c/p\u003e\u003cp\u003e\u003cb\u003eDomain Angle Analysis.\u003c/b\u003e To monitor the extent of the movement of the UBL domain with respect to the remaining ones, domain angles were defined using MDTraj[\u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e56\u003c/span\u003e] by calculating the angle (in radians) between the alpha carbons for three-residue sets. Two sets of residues were used (Figures S4-S5), each with one residue located in the UBL domain and the remaining two, including the vertex, in highly flexible regions of the thumb domain, i.e. \u0026ndash; the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:2\\)\u003c/span\u003e\u003c/span\u003e- and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:3\\)\u003c/span\u003e\u003c/span\u003e-helix.\u003c/p\u003e\u003cp\u003e\u003cb\u003eClustering and Occupancy\u003c/b\u003e. To find fragment binding spots, first, the trajectory was aligned to the protein backbone using the starting structure as the reference. The center of mass (COM) of each fragment was calculated using MDTraj.[\u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e56\u003c/span\u003e] Unsupervised clustering was performed on the fragment COM using the Density-Based Spatial Clustering of Applications with Noise algorithm, DBSCAN,[\u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e59\u003c/span\u003e, \u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e60\u003c/span\u003e] from the machine learning scikit-learn python package,[\u003cspan citationid=\"CR61\" class=\"CitationRef\"\u003e61\u003c/span\u003e] version 0.22.1. To belong to the same cluster, a 1.5-\u0026Aring; maximum distance between the COM of two fragments was set as a threshold, for a minimum of 60 ns over 300-ns individual simulations. The center of each cluster is rendered in Visual Molecular Dynamics (VMD) molecular visualization program[\u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e62\u003c/span\u003e] as a red van der Waals bead. Lastly, the equation below (Eq.\u0026nbsp;1) was then used to convert the total number of frames in a cluster into occupancy, or total amount of time spent in that spot. Global occupancy was calculated over all trajectories (2.7 \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\mu\\:\\)\u003c/span\u003e\u003c/span\u003es).\u003c/p\u003e\u003cp\u003e\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:Occupancy\\:\\%=\\:\\frac{\\#\\:Frames\\:in\\:a\\:cluster}{Total\\:Trajectory\\:Frames}\\:\\)\u003c/span\u003e\u003c/span\u003e (Eq.\u0026nbsp;1)\u003c/p\u003e\u003cp\u003e\u003cb\u003eCavity Detection\u003c/b\u003e. Pocket volume on MD trajectories was measured using the cavity detection and characterization tool MDpocket.[\u003cspan citationid=\"CR63\" class=\"CitationRef\"\u003e63\u003c/span\u003e] Images were produced using VMD,[\u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e62\u003c/span\u003e] version 1.9.3 or with The PyMOL Molecular Graphic System, version 2.5.2. (Schr\u0026ouml;dinger, LLC).\u003c/p\u003e\u003cp\u003e\u003cb\u003eWater thermodynamics\u003c/b\u003e. Water mapping calculations were carried out using SZMAP Version 1.6.3.0 from OpenEye Scientific software.[\u003cspan additionalcitationids=\"CR65\" citationid=\"CR64\" class=\"CitationRef\"\u003e64\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR66\" class=\"CitationRef\"\u003e66\u003c/span\u003e] The simulation frames containing a fragment-ligand in a given spot were collected, and the computation was performed on the simulation frame that was closest, in RMSD space, to the average structure of a particular spot. SPRUCE 1.4.0.0 (Spruce 1.6.0.3. OpenEye, Cadence Molecular Sciences, Santa Fe, NM) was used to generate biologically relevant units, add missing sidechains, adjust protonation states, enumerates alternate locations and tautomerization states for protein, ligands and co-factors. The output file has had charges and radii assigned. SZMAP grids were calculated for the binding pocket regions around the fragment. Then the stabilization calculations were performed for the complex, apo and ligand grids to describe where water molecules are influencing the binding process. The neutral difference in the free energy between the water and uncharged hypothetical water probes was computed to map the polarity throughout the binding site.\u003c/p\u003e\u003cp\u003e\u003cb\u003eContact analysis\u003c/b\u003e. The number of contact pairs was defined as the number of protein residue pairs with at least one heavy atom within 4 \u0026Aring; from each other. Contacts were calculated for all frames in which at least one fragment-ligand was found at any given binding spot, then averaged and rounded up to the nearest whole residue using PyEMMA.[\u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e58\u003c/span\u003e]\u003c/p\u003e\u003cp\u003e\u003cb\u003eExposons\u003c/b\u003e. Cooperative changes in the solvent exposure of protein residues were calculated using the exposons method,[\u003cspan citationid=\"CR67\" class=\"CitationRef\"\u003e67\u003c/span\u003e] which identifies functional conformational changes in protein structure ensembles. The method relies on identifying clusters of protein residues that are subjected to cooperative changes in their solvent exposure \u0026ndash; i.e., the exposons \u0026ndash; and deciphering the hierarchy of interactions that occur between them. Hence, in exposons calculations, the solvent accessible surface area (SASA) per residue-sidechain is computed. Then, correlations between pairs of residues are quantified by calculating a pairwise mutual information matrix. Next, the mutual information matrix is clustered into groups using affinity propagation clustering[\u003cspan citationid=\"CR68\" class=\"CitationRef\"\u003e68\u003c/span\u003e] with a damping value of 0.9. These groups of residues are often \u003cem\u003ecryptic pockets\u003c/em\u003e, pockets which only open in rare events. Exposons was calculated using the enspara python package.[\u003cspan citationid=\"CR67\" class=\"CitationRef\"\u003e67\u003c/span\u003e] The exposons calculation was performed on the combined set of all nine trajectories, analyzing the protein only. Exposons were ranked by the number of residues involved in the exposon.\u003c/p\u003e\u003c/div\u003e"},{"header":"RESULTS","content":"\u003cp\u003e\u003cb\u003eConformational dynamics of SARS-CoV-2 PLpro via MD \u0026ldquo;flooding\u0026rdquo; simulations\u003c/b\u003e. We performed nine independent 300-ns \"flooding\" MD simulations of SCoV2-PLpro bound to different fragment-ligands, accumulating a total simulation time of 2.7 \u0026micro;s. Three distinct fragments were investigated:[\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e] YRL (4-(2-hydroxyethyl)phenol) from PDB ID: 7OFS, HBA (p-hydroxybenzaldehyde) from PDB ID: 7OFT, and HE9 (methyl 3,4-bis(oxidanyl)benzoate) from PDB ID: 7OFU. The chemical structures of these fragments are illustrated in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e2\u003c/span\u003eD. Each fragment was simulated at three different concentrations: 10 mM, 50 mM, and 100 mM (Table\u0026nbsp;\u003cspan refid=\"Tab1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). For clarity, we designated these as different systems, where systems 1\u0026ndash;3 correspond to YRL at increasing concentrations, systems 4\u0026ndash;6 to HBA, and systems 7\u0026ndash;9 to HE9. The protein dynamics remained largely stable across most systems, with root-mean-square deviation (RMSD) values typically ranging between 2\u0026ndash;3 \u0026Aring; (Figure S1). However, systems 5 and 6 (HBA) and 8 (HE9) exhibited higher RMSD values, though this increased mobility appeared independent of both concentration and fragment type. Root-mean-square fluctuation (RMSF) analysis revealed a similar pattern: while systems 1\u0026ndash;3 (YRL) maintained consistent residue deviations across increasing ligand concentrations, systems 5\u0026ndash;6 (HBA) and 8 (HE9) displayed notably higher RMSFs in the UBL domain (residues 1\u0026ndash;62) (Figure S2). Additionally, systems 5 and 8 showed elevated RMSFs in the thumb domain (residues 120\u0026ndash;170), though this behavior appeared independent of fragment type or concentration. To further explore the protein dynamics and capture essential motions, we conducted principal component analysis (PCA) on the Cartesian coordinates of the protein backbone atoms (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e3\u003c/span\u003e). The analysis revealed that the first three principal components accounted for 59% of the total variance, with PC1 alone explaining 34% (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e3\u003c/span\u003eB). The strongest correlations with PC1 corresponded to motions of residues 23, 46, 26, 25, 24, and 45, all within the UBL domain. Examination of frames with minimum and maximum PC1 values showed nearly perfect alignment of residues 63\u0026ndash;315, excluding the UBL domain, indicating that PC1 primarily represents the forward movement of the UBL domain (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e3\u003c/span\u003eA, left panels). The second principal component (PC2) accounted for 18% of the variance (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e3\u003c/span\u003eB), with strongest correlations among residues 18, 17, 19, 9, and 31, also located in the UBL domain. Analysis of minimum and maximum PC2 structures revealed that this component corresponds to a downward movement and forward twisting of both the UBL domain and the adjacent helix-3 in the thumb domain (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e3\u003c/span\u003eA, right panels). While systems 1, 2, 3, 4, 7, and 9 maintained stability in a single energy basin within the PC1/PC2 space, systems 5, 6, and 8 explored multiple energy basins, moving in negative PC1 and PC2 directions, indicating substantial UBL domain movement. The third principal component (PC3) contributed 7% of the variance (Figure S3), with strongest correlations among residues 226, 225, 227, 224, and 191. These residues, located in the fingers domain, form part of the loops containing four cysteine residues responsible for tetrahedral coordination of a zinc ion (Zn\u003csup\u003e2+\u003c/sup\u003e) in a zinc finger motif. This motif has been established as essential for maintaining structural integrity and protease activity in CoV PLpro enzymes.[\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e, \u003cspan citationid=\"CR69\" class=\"CitationRef\"\u003e69\u003c/span\u003e] The motion captured by PC3, which was evenly sampled across all systems (Figure S3), corresponds to an inward curling of the fingers domain toward the active site.\u003c/p\u003e\u003cp\u003e\u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab1\" border=\"1\"\u003e\u003ccaption language=\"En\"\u003e\u003cdiv class=\"CaptionNumber\"\u003eTable 1\u003c/div\u003e\u003cdiv class=\"CaptionContent\"\u003e\u003cp\u003eSummary of MDFFr simulations and spots identified (total time 2.7\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\varvec{\\mu\\:}\\varvec{s})\\)\u003c/span\u003e\u003c/span\u003e.\u003c/p\u003e\u003c/div\u003e\u003c/caption\u003e\u003ccolgroup cols=\"6\"\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e\u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e\u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c5\" colnum=\"5\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c6\" colnum=\"6\"\u003e\u003c/div\u003e\u003cthead\u003e\u003ctr\u003e\u003cth align=\"left\" colname=\"c1\"\u003e\u003cp\u003eSystem\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c2\"\u003e\u003cp\u003ePDB-ID\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c3\"\u003e\u003cp\u003eFragment\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c4\"\u003e\u003cp\u003eConcentration (mM)\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c5\"\u003e\u003cp\u003eSimulation time (ns)\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c6\"\u003e\u003cp\u003eSpot #\u003c/p\u003e\u003c/th\u003e\u003c/tr\u003e\u003c/thead\u003e\u003ctbody\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eSys-1\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e7OFS\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eYRL\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e\u003cp\u003e10\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e\u003cp\u003e300\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003eS11\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eSys-2\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e7OFS\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eYRL\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e\u003cp\u003e50\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e\u003cp\u003e300\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e-\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eSys-3\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e7OFS\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eYRL\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e\u003cp\u003e100\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e\u003cp\u003e300\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003eS1, S4, S5, S7, S11, S12\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eSys-4\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e7OFT\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eHBA\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e\u003cp\u003e10\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e\u003cp\u003e300\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e-\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eSys-5\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e7OFT\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eHBA\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e\u003cp\u003e50\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e\u003cp\u003e300\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003eS4, S9, S10\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eSys-6\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e7OFT\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eHBA\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e\u003cp\u003e100\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e\u003cp\u003e300\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003eS11\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eSys-7\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e7OFU\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eHE9\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e\u003cp\u003e10\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e\u003cp\u003e300\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003eS2, S11\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eSys-8\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e7OFU\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eHE9\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e\u003cp\u003e50\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e\u003cp\u003e300\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003eS6, S8, S9\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eSys-9\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e7OFU\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eHE9\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e\u003cp\u003e100\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e\u003cp\u003e300\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003eS3, S11\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003c/tbody\u003e\u003c/colgroup\u003e\u003c/table\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eTo validate the UBL-related motions identified by PC1 and PC2, we analyzed specific inter-domain angles throughout the simulation trajectories. We monitored two key angles: (1) between the alpha carbons of residues 134\u0026rarr;84 (thumb domain)\u0026rarr;8 (UBL domain) to capture the forward motion of the UBL domain (Figure S4), and (2) between residues 88\u0026rarr;81 (thumb domain)\u0026rarr;5 (UBL domain) to track the downward UBL movement (Figure S5). The forward motion analysis revealed distinct behavioral patterns among the systems. Systems 1\u0026ndash;4, 7, and 9 maintained stable domain angles throughout the simulations (Figure S4B). In contrast, system 5 exhibited consistently large increases in domain-angle amplitude across the entire trajectory, while systems 6 and 8 showed significant angle increases toward the end of their simulation periods. The downward motion analysis showed similar system-dependent behavior. The angle tracking the UBL domain's downward movement remained stable in systems 1\u0026ndash;4, 7, and 9, while systems 5, 6, and 8 demonstrated pronounced downward UBL motion in the latter portions of their simulations (Figure S5B). These angular analyses provide quantitative support for the domain motions identified through PCA, confirming the significant conformational flexibility of the UBL domain in specific systems.\u003c/p\u003e\u003cp\u003e\u003cb\u003eIdentification and Characterization of Binding Spots for Fragment Ligands.\u003c/b\u003e We identified and characterized potential fragment binding sites on SARS-CoV2-PLpro through systematic analysis of our MDFFr simulations. The binding sites were determined by clustering the center of mass positions of flooded fragments throughout the simulation trajectories, revealing a total of 12 distinct binding spots (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003e). Each identified spot underwent comprehensive evaluation using multiple parameters:\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003col\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eSimulation occupancy: The duration of fragment binding during individual simulations.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eGlobal occupancy: The percentage of total simulation time (2,700 ns) during which a fragment remained bound.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eContact analysis: The number of protein residues within 4 \u0026Aring; heavy atom distance of the bound fragment.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eWater thermodynamics: The local solvent behavior and energetics at each binding site\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003c/ol\u003e\u003c/p\u003e\u003cp\u003eThis multi-parameter assessment approach enabled us to characterize not only the location of potential binding sites but also their physicochemical properties and binding stability, providing insights into their potential druggability.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003cb\u003eOccupancy and Contact Residues.\u003c/b\u003e\u003c/p\u003e\u003cp\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eSpot 1 (S1)\u003c/span\u003e, located on the back of the fingers domain between β8 and β12 (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003e), exhibits the lowest global occupancy at 2.41% (Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e) and was exclusively observed in system 3 (Table\u0026nbsp;\u003cspan refid=\"Tab1\" class=\"InternalRef\"\u003e1\u003c/span\u003e, Table S1). Fragments binding at S1 established contacts with an average of 8 residues (Table S2 and S3).\u003c/p\u003e\u003cp\u003e\u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab2\" border=\"1\"\u003e\u003ccaption language=\"En\"\u003e\u003cdiv class=\"CaptionNumber\"\u003eTable 2\u003c/div\u003e\u003cdiv class=\"CaptionContent\"\u003e\u003cp\u003eComparison of MDFFr spots by global occupancy (calculated by Eq.\u0026nbsp;1)\u003c/p\u003e\u003c/div\u003e\u003c/caption\u003e\u003ccolgroup cols=\"3\"\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e\u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e\u003cthead\u003e\u003ctr\u003e\u003cth align=\"left\" colname=\"c1\"\u003e\u003cp\u003eSpot #\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c2\"\u003e\u003cp\u003eExperimental\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c3\"\u003e\u003cp\u003eGlobal Occupancy (%)\u003c/p\u003e\u003c/th\u003e\u003c/tr\u003e\u003c/thead\u003e\u003ctbody\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eS1\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eno\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e\u003cp\u003e2.41\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eS2\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eno\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e\u003cp\u003e4.17\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eS3\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eno\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e\u003cp\u003e8.91\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eS4\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eno\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e\u003cp\u003e7.16\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eS5\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eno\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e\u003cp\u003e6.17\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eS6\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eno\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e\u003cp\u003e5.36\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eS7\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eno\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e\u003cp\u003e3.77\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eS8\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eno\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e\u003cp\u003e4.64\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eS9\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eno\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e\u003cp\u003e9.50\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eS10\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eno\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e\u003cp\u003e4.51\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eS11\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eYes (YRL and HBA)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e\u003cp\u003e49.60\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eS12\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eno\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e\u003cp\u003e7.47\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003c/tbody\u003e\u003c/colgroup\u003e\u003c/table\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eSpot 2 (S2)\u003c/span\u003e is positioned on the back of the palm domain (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003e), making contact with the loop between β13 and β14. With a low global occupancy of 4.17% and exclusive presence in system 7, fragments at S2 maintained contacts with an average of 8 residues (Table S2 and S3). Notably, this site represents a previously unidentified allosteric binding location.\u003c/p\u003e\u003cp\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eSpot 3 (S3)\u003c/span\u003e is situated at the base of the fingers where it meets the palm domain, contacting residues on β10 and α7. With the third-highest global occupancy at 8.91% (Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e), S3 was uniquely observed in system 9 and established an average of 6 contacts (Table S2 and S3). Significantly, four of its contact residues (ARG166, GLU167, GLU203, and MET208) overlap with the experimentally observed UBS1 site (PDB ID: 6XAA,[\u003cspan citationid=\"CR70\" class=\"CitationRef\"\u003e70\u003c/span\u003e] UBS1 bound to UB; and PDB ID: 6YVA,[\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e] UBS2 bound to ISG15; Table S4), representing the most extensive overlap with these functionally crucial protein-protein interaction sites.\u003c/p\u003e\u003cp\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eSpot 4 (S4)\u003c/span\u003e is positioned at the fingers-thumb interface, establishing contacts with β8 and β12 at their respective termini (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003e). With a global occupancy of 7.16%, S4 appears in systems 3 and 5, maintaining an average of 4 contacts (Table S2 and S3).\u003c/p\u003e\u003cp\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eSpot 5 (S5)\u003c/span\u003e, located on the thumb domain between α3 and α6, shows a global occupancy of 6.17% and appears exclusively in system 8. S5 maintains an average of 5 contacts, with two residues (ASN156, LYS157) overlapping with the experimentally observed UBS1 site (Table S1V).\u003c/p\u003e\u003cp\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eSpot 6 (S6)\u003c/span\u003e is positioned behind the thumb domain between α5, α6, and the UBL domain. With a global occupancy of 5.36% and exclusive presence in system 8, S6 establishes an average of 9 contacts (Table S2 and S3).\u003c/p\u003e\u003cp\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eSpot 7 (S7)\u003c/span\u003e, situated on the front of α3 in the thumb domain, displays the second-lowest global occupancy at 3.77% and appears only in system 3. Along with S4, it maintains the lowest average number of contacts at 4 (Tables S2).\u003c/p\u003e\u003cp\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eSpot 8 (S8)\u003c/span\u003e is located on the thumb domain's posterior surface, just below α2 (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003e). With a low global occupancy of 4.64% and exclusive presence in system 8 (Table S1), S8 maintains 6 contacts, including GLU70, which overlaps with the experimentally observed UBS2 site and establishes contacts with ISG15 (6YVA) (Table S2 and S3).\u003c/p\u003e\u003cp\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eSpot 9 (S9)\u003c/span\u003e occupies a deep cavity between the thumb and UBL domains, positioned directly below S11, the nearest site to the experimental fragment binding pocket (Figs.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e1\u003c/span\u003e and \u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003e). With the second-highest global occupancy at 9.50%, S9 appears in systems 5 and 8 and maintains the highest number of contacts at 12 (Table S2 and S3).\u003c/p\u003e\u003cp\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eSpot 10 (S10)\u003c/span\u003e is positioned beneath α2, toward the posterior interface between thumb and UBL domains. With a low global occupancy of 4.51%, S10 appears only in system 5 and maintains an average of 5 contacts (Tables S2).\u003c/p\u003e\u003cp\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eSpot 11 (S11)\u003c/span\u003e resides in a deep cavity between the thumb and UBL domains. It demonstrates the highest global occupancy by far at 49.60% and appears in systems 1, 3, 6, 7, 8, and 9. S11 lies closest to the experimental YRL and HBA binding site (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e2\u003c/span\u003eA) and maintains the second-highest number of contacts, with 8 residues in each system where it appears (Tables S2). Additional experimental structures with co-crystallized ligands are known.[\u003cspan citationid=\"CR71\" class=\"CitationRef\"\u003e71\u003c/span\u003e]\u003c/p\u003e\u003cp\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eSpot 12 (S12)\u003c/span\u003e is located between α3 helix in the thumb domain and the UBL domain. With the fourth-highest global occupancy at 7.47%, S12 appears only in system 3 and maintains an average of 6 contacts. S12 lies nearest to S11, with the spots nearly overlapping (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003e and Table S2 and S3). Both spots include contacts with residues found in the experimentally known Ub-S2 site (6YVA and SARS 5E6J, Table S3), with S11 contacting ARG65, PHE69, and PRO77, the latter also being part of S12's contacts (Table S2 and S3).\u003c/p\u003e\u003cp\u003e\u003cb\u003eWater Thermodynamics.\u003c/b\u003e To characterize the physicochemical properties of the binding spots identified through MDFFr, we performed SZMAP calculations on representative frames closest to the average structure of each spot (Fig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e5\u003c/span\u003e). SZMAP analyzes water thermodynamics to map binding site polarity, where negative regions (yellow, blue, and green) indicate favorable interactions with polar groups, while positive regions (red) favor non-polar interactions.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eThe analysis of spots S1-S4 revealed distinct polarity patterns. In S1, while the phenol hydroxyl occupies an unfavorable hydrophobic pocket, the hydroxyethyl hydroxyl extends into a hydrophilic region. S2 features a hydrophobic pocket accommodating the fragment's bulk, with its carboxyl group engaging in polar interactions. S3 comprises a large hydrophobic pocket where the presence of oxygen atoms creates a polarity mismatch, though potential expansion into adjacent polar regions appears possible through extended polar groups. In S4, the binding mode combines a hydrophobic pocket housing the aromatic ring with a hydrophilic pocket accommodating the phenol hydroxyl.\u003c/p\u003e\u003cp\u003eSpots S5-S6 show complementary physicochemical properties. S5 positions YRL in a hydrophobic pocket flanked by hydrophilic regions near both hydroxyl groups. In S6, UBL domain displacement creates a large hydrophobic pocket where HE9 binds between α5, α6, and the UBL loop, with its alcohol groups extending into solvent-exposed regions and its carbonyl group engaging with polar moieties.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eThe remaining spots (S7-S12) display varied thermodynamic profiles. S7 presents a narrow hydrophobic channel along α3, while S8 forms a weakly hydrophobic, solvent-exposed pocket near α2. S9 features an extensive hydrophobic cavity extending from S11, offering considerable space for fragment expansion. S10, characterized by significant UBL displacement, shares similarities with S6 in its interaction pattern with HBA engaging the loop between β1 and β2, but differs in interacting with α2 rather than α5's bottom region. The highly occupied S11 forms a large hydrophobic pocket between UBL, α2, and α3, with strategic hydrophilic contact points that could guide fragment optimization. S12 presents a compact hydrophobic pocket with a hydrophilic patch accommodating the phenol hydroxyl, and its proximity to S11 suggests potential for fragment linking strategies.\u003c/p\u003e\n\u003ch3\u003eComparison of computational versus experimental binding sites \u0026 cryptic pockets\u003c/h3\u003e\n\u003cp\u003e\u003cb\u003eExperimental Binding\u003c/b\u003e. Our MDFFr simulations revealed significant correlations with experimental binding data. Spot 11 (S11), while positioned slightly above the crystallographic pose, emerged as the closest fragment-binding hot-spot to the experimentally determined positions of YRL and HBA. Similarly, Spot 9 (S9) closely approximates these crystallographic poses, though located just below the experimental binding position. The significance of these sites is reinforced by their occupancy statistics: S11 shows dominant binding with 49.60% global occupancy (Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e), while S9 follows with 9.50%. Their combined occupancy of 59.20% substantially exceeds that of all other identified binding spots, strongly validating their biological relevance.\u003c/p\u003e\u003cp\u003eThe binding dynamics at S11 showed system-dependent behavior. In systems 1, 3, and 6, fragments exhibited rapid binding to S11 and maintained stable occupancy throughout the simulation (Figure S8). However, several systems showed no fragment binding at S11, primarily due to conformational changes leading to pocket collapse (Figure S9). Specific examples include sys-2, where movements of PRO59, LEU80, and TYR72 caused the pocket's upper region to collapse; sys-4, where similar collapse prevented fragment access; and sys-5, where ubiquitin domain rotation resulted in pocket closure (Figure S5).\u003c/p\u003e\u003cp\u003eRegarding HE9 binding, our simulations revealed interesting disparities with experimental data. While none of our identified spots directly corresponded to the crystallographic HE9 pose, which uniquely occupies a solvent-exposed position at the ISG15 interface, our simulations showed HE9 preferentially binding to S11, near the crystallographic poses of YRL and HBA (Figure S10). This finding is particularly noteworthy given that all three compounds (YRL, HBA, and HE9) experimentally bind at the ISG15/Ub-S2 protein-protein interaction site, though HE9 binds on the opposite face of the α2 helix compared to the concave recession where YRL and HBA bind.[\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e] Importantly, the experimental HE9 binding location coincides with the outer, solvent-exposed region of a continuous binding surface formed by the merger of S11, S9, and S10.\u003c/p\u003e\u003cp\u003eThe absence of spots corresponding to the experimentally known GRL binding site likely stems from simulation time limitations. With individual simulations of 300 ns (totaling 2.7 \u0026micro;s), the BL2 loop between β14 and β15 may not have had sufficient time to adopt an open conformation, thus preventing access to the GRL binding pocket.\u003c/p\u003e\u003cp\u003e\u003cb\u003eCryptic Pockets\u003c/b\u003e. To investigate potential functional changes in PLpro conformational dynamics revealed by MDFFr, we employed \"exposons\" calculations.[\u003cspan citationid=\"CR67\" class=\"CitationRef\"\u003e67\u003c/span\u003e] This analytical method identifies groups of residues that exhibit strong correlations in their solvent accessible surface area changes. These exposons, which include the previously mentioned GRL binding site, typically manifest only during rare conformational events.\u003c/p\u003e\u003cp\u003eOur analysis revealed three significant exposons with distinct characteristics. The largest exposon (shown in red in Fig.\u0026nbsp;\u003cspan refid=\"Fig12\" class=\"InternalRef\"\u003e6\u003c/span\u003e) corresponds to the BL2 loop opening that exposes the active site. This finding is particularly significant as it demonstrates that our MDFFr trajectories, despite individual simulation times of 300 ns, can collectively (2.7 \u0026micro;s) capture cooperative residue movements associated with functional events. The second largest exposon (shown in blue in Fig.\u0026nbsp;\u003cspan refid=\"Fig12\" class=\"InternalRef\"\u003e6\u003c/span\u003e) maps to the Ub-S2 binding site at the base of the ISG15 interaction site, overlapping with our identified spots S5, S7, S9, S11, and S12. This exposon correlates with the UBL domain's movement away from the thumb domain, as confirmed by our domain angle analysis (Figs. S4 and S5). Notably, this region encompasses both the experimentally validated binding sites for YRL and HBA, and includes S11 and S9, our highest occupancy binding spots (Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eThe third largest exposon (shown in green in Fig.\u0026nbsp;\u003cspan refid=\"Fig12\" class=\"InternalRef\"\u003e6\u003c/span\u003e) coincides with our S2 spot (Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e), positioned on the thumb domain's posterior surface, beneath the UBL domain (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003e). The functional significance of this concerted amino acid movement warrants further investigation, as it may reveal previously unknown allosteric regulatory mechanisms.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e"},{"header":"DISCUSSION AND CONCLUSIONS","content":"\u003cp\u003eWe have developed a computational framework, called Molecular Dynamics Simulations with Flooding Fragments (MDFFr), to accurately identify druggable hot-spots for fragment-ligand binding based solely on the 3D structure of the protein, without prior knowledge of binding site(s). This physics-based approach, which evolves from previous studies conducted on membrane proteins,[\u003cspan additionalcitationids=\"CR45 CR46\" citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e] has been significantly modified and generalized for use with globular proteins, as demonstrated with the viral PLpro enzyme of SARS-CoV2 (our case study).\u003c/p\u003e\u003cp\u003eFrom a biological perspective, MDFFr identified 12 putative ligand binding sites (S1 to S12; Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003e and Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e) on PLpro. Two sites exhibited notably high overall occupancy: S11 (~\u0026thinsp;50%) and S9 (~\u0026thinsp;10%), with their individual and combined occupancies significantly exceeding all other spots (Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e). Notably, S11 corresponds to the experimentally validated binding site for two fragment-probes (YRL and HBA) used in our MD studies (Figs.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e2\u003c/span\u003eA and \u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003e; Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e), providing proof-of-concept for MDFFr's application in blind-challenge MD studies (with random fragment placement and no prior binding knowledge). Analysis of the overall site mapping on PLpro (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003e) revealed that S11, S9, and S10 merge to form a continuous binding pocket (combined global occupancy\u0026thinsp;~\u0026thinsp;64%) anchored within the thumb domain at the interface of helices α3, α5, and α6 (deeply buried), and the bordering α2 helix (solvent-exposed). While this pocket partially overlaps with the shallow Ub-S2 interaction site, the YRL/HBA binding site, and the thiosemicarbazones derivatives binding site (PDB ID: 7QCI)[\u003cspan citationid=\"CR71\" class=\"CitationRef\"\u003e71\u003c/span\u003e] (S10-S11), it extends inward to form a previously unobserved buried groove (S9).\u003c/p\u003e\u003cp\u003eA particularly significant finding is the S2 fragment-binding spot (third highest occupancy, ~\u0026thinsp;4%; Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e), positioned near but not overlapping with the experimentally known Ub-S1 protein-protein interaction site (Figs.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e2\u003c/span\u003eC and \u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003e). This newly identified site coincides with the third largest cluster of residues showing cooperative changes in solvent accessibility (exposon), suggesting potential roles as a conformational switch, allosteric site, or cryptic pocket.[\u003cspan citationid=\"CR67\" class=\"CitationRef\"\u003e67\u003c/span\u003e] Interestingly, S2 is proximal to (but not overlapping with) an allosteric site recently mapped on PLpro.[\u003cspan citationid=\"CR72\" class=\"CitationRef\"\u003e72\u003c/span\u003e] This site was identified through computational structural analyses of PLpro experimental structures and validated through enzyme titration experiments of computationally predicted protein mutants.[\u003cspan citationid=\"CR72\" class=\"CitationRef\"\u003e72\u003c/span\u003e] While complete validation of S2 as a new allosteric regulatory site requires additional computational and experimental work, these findings support S2's significance and demonstrate MDFFr's utility in physics-based mapping of novel allosteric/functional sites.\u003c/p\u003e\u003cp\u003eMethodologically, MDFFr adapts the \"flooding MD simulations\"[\u003cspan additionalcitationids=\"CR45 CR46\" citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e] concept of using known or putative binder molecules to reveal binding sites or access pathways (e.g., cocrystal-fragments to PLpro or volatile anesthetics to ion channels). Unlike the original approach,[\u003cspan additionalcitationids=\"CR45 CR46\" citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e] MDFFr offers broader applicability to all proteins, including non-membrane bound ones. Our PLpro application demonstrates MDFFr's potential as a general MD-based fragment mapping approach, suitable for \"computational-fragment screening\" to complement experimental high-throughput screening and support hit-to-lead optimizations.\u003c/p\u003e\u003cp\u003eCurrent MD-based binding-site mapping approaches[\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e, \u003cspan additionalcitationids=\"CR36 CR37 CR38\" citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e, \u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e, \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e, \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e, \u003cspan citationid=\"CR73\" class=\"CitationRef\"\u003e73\u003c/span\u003e] typically employ large-scale simulations with classical or accelerated MD implementations, using cosolvent \"probes\" to map protein surface features. While effective, these methods face challenges in cosolvent selection and computational costs. In contrast, MDFFr offers straightforward implementation with any forcefield and MD engine, utilizing a \"building block\" approach with short simulation runs (300 ns per system) and flexible fragment concentrations. Important considerations for MDFFr include: (i) lower cosolvent fragment concentrations compared to membrane-specific \"flooding MD simulations\",[\u003cspan additionalcitationids=\"CR45 CR46\" citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e] and (ii) focus on binding hot spot location prediction rather than binding affinity correlation.\u003c/p\u003e\u003cp\u003eFuture work will expand MDFFr applications to diverse protein sets with varying fragment types and concentrations, while evaluating its effectiveness in computational fragment screening and protein-protein interaction site identification, including systems with additional biomacromolecules like DNA/RNA.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eACKNOWLEDGMENTS\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis research was supported in part by The RF Research Foundation and the City University of New York (CUNY) Queens College startup funds (E.G.); the University of Wyoming (UWYO) and School of Pharmacy (SOP) startup funds (K.E.). Computational resources were provided by Temple University\u0026apos;s high-performance computing facility, supported in part by the National Science Foundation through major research instrumentation grant number 1625061. The authors gratefully acknowledge OpenEye Scientific Software for providing a free academic license.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAUTHORS CONTRIBUTIONS\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eConceptualization: Jason Pattis, Khaled Elokely, Eleonora Gianti. Methodology: Jason Pattis, Khaled Elokely, Eleonora Gianti. Formal analysis and investigation: Jason Pattis, Khaled Elokely, Eleonora Gianti. Writing: Jason Pattis, Khaled Elokely, Eleonora Gianti.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCOMPLIANCE WITH ETHICAL STANDARDS\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting Interests\u003c/strong\u003e: The authors declare no conflict of interest.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eDATA AVAILABILITY\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAll data produced by this study are disclosed in the publication and its associated supporting information (SI) file. The code for these analyses and input files are available at: https://github.com/jgpattis/MDFFr\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eZhou P, Yang X-L, Wang X-G, et al (2020) A pneumonia outbreak associated with a new coronavirus of probable bat origin. Nature 579:270\u0026ndash;273. https://doi.org/10.1038/s41586-020-2012-7\u003c/li\u003e\n\u003cli\u003eWu F, Zhao S, Yu B, et al (2020) A new coronavirus associated with human respiratory disease in China. Nature 579:265\u0026ndash;269. https://doi.org/10.1038/s41586-020-2008-3\u003c/li\u003e\n\u003cli\u003eLu R, Zhao X, Li J, et al (2020) Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding. The Lancet 395:565\u0026ndash;574. https://doi.org/10.1016/S0140-6736(20)30251-8\u003c/li\u003e\n\u003cli\u003eVegivinti CTR, Evanson KW, Lyons H, et al (2022) Efficacy of antiviral therapies for COVID-19: a systematic review of randomized controlled trials. BMC Infect. Dis. 22:1\u0026ndash;45\u003c/li\u003e\n\u003cli\u003eHwang Y-C, Lu R-M, Su S-C, et al (2022) Monoclonal antibodies for COVID-19 therapy and SARS-CoV-2 detection. J. Biomed. Sci. 29:1\u0026ndash;50\u003c/li\u003e\n\u003cli\u003eLei J, Kusov Y, Hilgenfeld R (2018) Nsp3 of coronaviruses: Structures and functions of a large multi-domain protein. Antiviral Res 149:58\u0026ndash;74. https://doi.org/10.1016/j.antiviral.2017.11.001\u003c/li\u003e\n\u003cli\u003eHilgenfeld R (2014) From SARS to MERS: crystallographic studies on coronaviral proteases enable antiviral drug design. FEBS J 281:4085\u0026ndash;4096. https://doi.org/10.1111/febs.12936\u003c/li\u003e\n\u003cli\u003eNarayanan A, Narwal M, Majowicz SA, et al (2022) Identification of SARS-CoV-2 inhibitors targeting Mpro and PLpro using in-cell-protease assay. Commun. Biol. 5:1\u0026ndash;17\u003c/li\u003e\n\u003cli\u003eBarretto N, Jukneliene D, Ratia K, et al (2005) The Papain-Like Protease of Severe Acute Respiratory Syndrome Coronavirus Has Deubiquitinating Activity\u003c/li\u003e\n\u003cli\u003eLindner HA, Lytvyn V, Qi H, et al (2007) Selectivity in ISG15 and ubiquitin recognition by the SARS coronavirus papain-like protease. Arch Biochem Biophys 466:8\u0026ndash;14\u003c/li\u003e\n\u003cli\u003eDel Valle DM, Kim-Schulze S, Huang H-H, et al (2020) An inflammatory cytokine signature predicts COVID-19 severity and survival. Nat. Med. 26:1636\u0026ndash;1643\u003c/li\u003e\n\u003cli\u003eShin D, Mukherjee R, Grewe D, et al (2020) Papain-like protease regulates SARS-CoV-2 viral spread and innate immunity. Nature 587:657\u0026ndash;662\u003c/li\u003e\n\u003cli\u003eRut W, Lv Z, Zmudzinski M, et al (2020) Activity profiling and crystal structures of inhibitor-bound SARS-CoV-2 papain-like protease: A framework for anti-COVID-19 drug design. Sci Adv 6\u003c/li\u003e\n\u003cli\u003eHarcourt BH, Jukneliene D, Kanjanahaluethai A, et al (2004) Identification of Severe Acute Respiratory Syndrome Coronavirus Replicase Products and Characterization of Papain-Like Protease Activity. J Virol 78:13600\u0026ndash;13612. https://doi.org/10.1128/jvi.78.24.13600-13612.2004\u003c/li\u003e\n\u003cli\u003eB\u0026aacute;ez-Santos YM, St. John SE, Mesecar AD (2015) The SARS-coronavirus papain-like protease: Structure, function and inhibition by designed antiviral compounds. Antiviral Res 115:21\u0026ndash;38. https://doi.org/10.1016/j.antiviral.2014.12.015\u003c/li\u003e\n\u003cli\u003eOsipiuk J, Azizi SA, Dvorkin S, et al (2021) Structure of papain-like protease from SARS-CoV-2 and its complexes with non-covalent inhibitors. Nat Commun 12:743\u003c/li\u003e\n\u003cli\u003eGao X, Qin B, Chen P, et al (2021) Crystal structure of SARS-CoV-2 papain-like protease. Acta Pharm Sin B 11:237\u0026ndash;245. https://doi.org/10.1016/j.apsb.2020.08.014\u003c/li\u003e\n\u003cli\u003eDevaraj SG, Wang N, Chen Z, et al (2007) Regulation of IRF-3-dependent Innate Immunity by the Papain-like Protease Domain of the Severe Acute Respiratory Syndrome Coronavirus *. J Biol Chem 282:32208\u0026ndash;32221. https://doi.org/10.1074/jbc.M704870200\u003c/li\u003e\n\u003cli\u003eClementz MA, Chen Z, Banach BS, et al (2010) Deubiquitinating and Interferon Antagonism Activities of Coronavirus Papain-Like Proteases. J Virol 84:4619\u0026ndash;4629. https://doi.org/10.1128/jvi.02406-09\u003c/li\u003e\n\u003cli\u003eFrieman M, Ratia K, Johnston RE, et al (2009) Severe Acute Respiratory Syndrome Coronavirus Papain-Like Protease Ubiquitin-Like Domain and Catalytic Domain Regulate Antagonism of IRF3 and NF-\u0026kappa;B Signaling. J Virol 83:6689\u0026ndash;6705. https://doi.org/10.1128/jvi.02220-08\u003c/li\u003e\n\u003cli\u003eShen Z, Ratia K, Cooper L, et al (2021) Design of SARS-CoV-2 PLpro Inhibitors for COVID-19 Antiviral Therapy Leveraging Binding Cooperativity. J Med Chem\u003c/li\u003e\n\u003cli\u003eRatia K, Pegan S, Takayama J, et al (2008) A noncovalent class of papain-like protease/deubiquitinase inhibitors blocks SARS virus replication. Proc Natl Acad Sci 105:16119\u0026ndash;16124. https://doi.org/10.1073/pnas.0805240105\u003c/li\u003e\n\u003cli\u003eSanders BC, Pokhrel S, Labbe AD, et al (2023) Potent and selective covalent inhibition of the papain-like protease from SARS-CoV-2. Nat Commun 14:1733. https://doi.org/10.1038/s41467-023-37254-w\u003c/li\u003e\n\u003cli\u003ePatchett S, Lv Z, Rut W, et al (2021) A molecular sensor determines the ubiquitin substrate specificity of SARS-CoV-2 papain-like protease. Cell Rep 36:109754\u003c/li\u003e\n\u003cli\u003eSrinivasan V, Brognaro H, Prabhu PR, et al (2022) Antiviral activity of natural phenolic compounds in complex at an allosteric site of SARS-CoV-2 papain-like protease. Commun. Biol. 5:805\u003c/li\u003e\n\u003cli\u003eBrenke R, Kozakov D, Chuang G-Y, et al (2009) Fragment-based identification of druggable \u0026lsquo;hot spots\u0026rsquo; of proteins using Fourier domain correlation techniques. Bioinformatics 25:621\u0026ndash;627. https://doi.org/10.1093/bioinformatics/btp036\u003c/li\u003e\n\u003cli\u003eErlanson DA, McDowell RS, O\u0026rsquo;Brien T (2004) Fragment-Based Drug Discovery. J Med Chem 47:3463\u0026ndash;3482. https://doi.org/10.1021/jm040031v\u003c/li\u003e\n\u003cli\u003eErlanson DA, Jahnke W (2015) Fragment-based Drug Discovery: Lessons and Outlook. John Wiley \u0026amp; Sons\u003c/li\u003e\n\u003cli\u003eBottegoni G, Favia AD, Recanatini M, Cavalli A (2012) The role of fragment-based and computational methods in polypharmacology. Drug Discov Today 17:23\u0026ndash;34. https://doi.org/10.1016/j.drudis.2011.08.002\u003c/li\u003e\n\u003cli\u003eZauhar RJ, Moyna G, Tian L, et al (2003) Shape Signatures: A New Approach to Computer-Aided Ligand- and Receptor-Based Drug Design. J Med Chem 46:5674\u0026ndash;5690. https://doi.org/10.1021/jm030242k\u003c/li\u003e\n\u003cli\u003eGianti E, Sartori L (2008) Identification and Selection of \u0026ldquo;Privileged Fragments\u0026rdquo; Suitable for Primary Screening. J Chem Inf Model 48:2129\u0026ndash;2139. https://doi.org/10.1021/ci800219h\u003c/li\u003e\n\u003cli\u003eGianti E, Messick TE, Lieberman PM, Zauhar RJ (2016) Computational analysis of EBNA1 \u0026ldquo;druggability\u0026rdquo; suggests novel insights for Epstein-Barr virus inhibitor design. J Comput Aided Mol Des 30:285\u0026ndash;303. https://doi.org/10.1007/s10822-016-9899-y\u003c/li\u003e\n\u003cli\u003eBissaro M, Sturlese M, Moro S (2020) The rise of molecular simulations in fragment-based drug design (FBDD): an overview. Drug Discov Today 25:1693\u0026ndash;1701. https://doi.org/10.1016/j.drudis.2020.06.023\u003c/li\u003e\n\u003cli\u003ede Souza Neto LR, Moreira-Filho JT, Neves BJ, et al (2020) In silico Strategies to Support Fragment-to-Lead Optimization in Drug Discovery. Front Chem 8:. https://doi.org/10.3389/fchem.2020.00093\u003c/li\u003e\n\u003cli\u003eYanagisawa K, Moriwaki Y, Terada T, Shimizu K (2021) EXPRORER: Rational Cosolvent Set Construction Method for Cosolvent Molecular Dynamics Using Large-Scale Computation. J Chem Inf Model 61:2744\u0026ndash;2753. https://doi.org/10.1021/acs.jcim.1c00134\u003c/li\u003e\n\u003cli\u003eBakan A, Nevins N, Lakdawala AS, Bahar I (2012) Druggability Assessment of Allosteric Proteins by Dynamics Simulations in the Presence of Probe Molecules. J Chem Theory Comput 8:2435\u0026ndash;2447. https://doi.org/10.1021/ct300117j\u003c/li\u003e\n\u003cli\u003eZariquiey FS, Souza JV de, Bronowska AK (2019) Cosolvent Analysis Toolkit (CAT): a robust hotspot identification platform for cosolvent simulations of proteins to expand the druggable proteome. Sci Rep 9:19118. https://doi.org/10.1038/s41598-019-55394-2\u003c/li\u003e\n\u003cli\u003eRaman EP, Yu W, Lakkaraju SK, MacKerell ADJr (2013) Inclusion of Multiple Fragment Types in the Site Identification by Ligand Competitive Saturation (SILCS) Approach. J Chem Inf Model 53:3384\u0026ndash;3398. https://doi.org/10.1021/ci4005628\u003c/li\u003e\n\u003cli\u003eGhanakota P, Carlson HA (2016) Moving Beyond Active-Site Detection: MixMD Applied to Allosteric Systems. J Phys Chem B 120:8685\u0026ndash;8695. https://doi.org/10.1021/acs.jpcb.6b03515\u003c/li\u003e\n\u003cli\u003eGhanakota P, Carlson HA (2016) Driving Structure-Based Drug Discovery through Cosolvent Molecular Dynamics. J Med Chem 59:10383\u0026ndash;10399. https://doi.org/10.1021/acs.jmedchem.6b00399\u003c/li\u003e\n\u003cli\u003eGhanakota P, DasGupta D, Carlson HA (2019) Free Energies and Entropies of Binding Sites Identified by MixMD Cosolvent Simulations. J Chem Inf Model 59:2035\u0026ndash;2045. https://doi.org/10.1021/acs.jcim.8b00925\u003c/li\u003e\n\u003cli\u003eSmith RD, Carlson HA (2021) Identification of Cryptic Binding Sites Using MixMD with Standard and Accelerated Molecular Dynamics. J Chem Inf Model 61:1287\u0026ndash;1299. https://doi.org/10.1021/acs.jcim.0c01002\u003c/li\u003e\n\u003cli\u003eGraham SE, Leja N, Carlson HA (2018) MixMD Probeview: Robust Binding Site Prediction from Cosolvent Simulations. J Chem Inf Model 58:1426. https://doi.org/10.1021/acs.jcim.8b00265\u003c/li\u003e\n\u003cli\u003eVemparala S, Domene C, Klein ML (2008) Interaction of Anesthetics with Open and Closed Conformations of a Potassium Channel Studied via Molecular Dynamics and Normal Mode Analysis. Biophys J 94:4260\u0026ndash;4269. https://doi.org/10.1529/biophysj.107.119958\u003c/li\u003e\n\u003cli\u003eRaju SG, Barber AF, LeBard DN, et al (2013) Exploring Volatile General Anesthetic Binding to a Closed Membrane-Bound Bacterial Voltage-Gated Sodium Channel via Computation. PLoS Comput Biol 9:e1003090. https://doi.org/10.1371/journal.pcbi.1003090\u003c/li\u003e\n\u003cli\u003eYang E, Granata D, Eckenhoff RG, et al (2018) Propofol inhibits prokaryotic voltage-gated Na+ channels by promoting activation-coupled inactivation. J Gen Physiol 150:1299\u0026ndash;1316. https://doi.org/10.1085/jgp.201711924\u003c/li\u003e\n\u003cli\u003eBarber AF, Carnevale V, Klein ML, et al (2014) Modulation of a voltage-gated Na+ channel by sevoflurane involves multiple sites and distinct mechanisms. Proc Natl Acad Sci 111:6726\u0026ndash;6731. https://doi.org/10.1073/pnas.1405768111\u003c/li\u003e\n\u003cli\u003eDolinsky TJ, Czodrowski P, Li H, et al (2007) PDB2PQR: expanding and upgrading automated preparation of biomolecular structures for molecular simulations. Nucleic Acids Res 35:W522\u0026ndash;W525. https://doi.org/10.1093/nar/gkm276\u003c/li\u003e\n\u003cli\u003eAbraham MJ, Murtola T, Schulz R, et al (2015) GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 1\u0026ndash;2:19\u0026ndash;25. https://doi.org/10.1016/j.softx.2015.06.001\u003c/li\u003e\n\u003cli\u003eJorgensen WL, Chandrasekhar J, Madura JD, et al (1983) Comparison of simple potential functions for simulating liquid water. J Chem Phys 79:926. https://doi.org/10.1063/1.445869\u003c/li\u003e\n\u003cli\u003eParrinello M, Rahman A (1981) Polymorphic transitions in single crystals: A new molecular dynamics method. J Appl Phys 52:7182\u0026ndash;7190. https://doi.org/10.1063/1.328693\u003c/li\u003e\n\u003cli\u003eHuang J, Rauscher S, Nawrocki G, et al (2017) CHARMM36m: an improved force field for folded and intrinsically disordered proteins. Nat Methods 14:71\u0026ndash;73. https://doi.org/10.1038/nmeth.4067\u003c/li\u003e\n\u003cli\u003eVanommeslaeghe K, Hatcher E, Acharya C, et al (2009) CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J Comput Chem NA-NA. https://doi.org/10.1002/jcc.21367\u003c/li\u003e\n\u003cli\u003eVanommeslaeghe K, MacKerell ADJr (2012) Automation of the CHARMM General Force Field (CGenFF) I: Bond Perception and Atom Typing. J Chem Inf Model 52:3144\u0026ndash;3154. https://doi.org/10.1021/ci300363c\u003c/li\u003e\n\u003cli\u003eVanommeslaeghe K, Raman EP, MacKerell ADJr (2012) Automation of the CHARMM General Force Field (CGenFF) II: Assignment of Bonded Parameters and Partial Atomic Charges. J Chem Inf Model 52:3155\u0026ndash;3168. https://doi.org/10.1021/ci3003649\u003c/li\u003e\n\u003cli\u003eMcGibbon RT, Beauchamp KA, Harrigan MP, et al (2015) MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories. Biophys J 109:1528\u0026ndash;32\u003c/li\u003e\n\u003cli\u003eHunter JD (2007) Matplotlib: A 2D Graphics Environment. Comput. Sci. Eng. 9:90\u0026ndash;95\u003c/li\u003e\n\u003cli\u003eScherer MK, Trendelkamp-Schroer B, Paul F, et al (2015) PyEMMA 2: A Software Package for Estimation, Validation, and Analysis of Markov Models. J Chem Theory Comput 11:5525\u0026ndash;42\u003c/li\u003e\n\u003cli\u003eEster M, Kriegel H-P, Xu X A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise\u003c/li\u003e\n\u003cli\u003eSchubert E, Sander J, Ester M, et al (2017) DBSCAN Revisited, Revisited: Why and How You Should (Still) Use DBSCAN. ACM Trans Database Syst 42:19:1-19:21. https://doi.org/10.1145/3068335\u003c/li\u003e\n\u003cli\u003ePedregosa F, Varoquaux G, Gramfort A, et al (2011) Scikit-learn: Machine Learning in Python. J Mach Learn Res 12:2825\u0026ndash;2830\u003c/li\u003e\n\u003cli\u003eHumphrey W, Dalke A, Schulten K (1996) VMD: visual molecular dynamics. J Mol Graph 14:33\u0026ndash;38, 27\u0026ndash;28\u003c/li\u003e\n\u003cli\u003eSchmidtke P, Bidon-Chanal A, Luque FJ, Barril X (2011) MDpocket: open-source cavity detection and characterization on molecular dynamics trajectories. Bioinformatics 27:3276\u0026ndash;3285. https://doi.org/10.1093/bioinformatics/btr550\u003c/li\u003e\n\u003cli\u003eTanger JCI, Pitzer KS (1989) Calculation of the thermodynamic properties of aqueous electrolytes to 1000.degree.C and 5000 bar from a semicontinuum model for ion hydration. J Phys Chem 93:4941\u0026ndash;4951. https://doi.org/10.1021/j100349a053\u003c/li\u003e\n\u003cli\u003eRashin AA, Bukatin MA (1991) Continuum based calculations of hydration entropies and the hydrophobic effect. J Phys Chem 95:2942\u0026ndash;2944. https://doi.org/10.1021/j100161a002\u003c/li\u003e\n\u003cli\u003eGrant JA, Pickup BT, Nicholls A (2001) A smooth permittivity function for Poisson\u0026ndash;Boltzmann solvation methods. J Comput Chem 22:608\u0026ndash;640. https://doi.org/10.1002/jcc.1032\u003c/li\u003e\n\u003cli\u003ePorter JR, Moeder KE, Sibbald CA, et al (2019) Cooperative Changes in Solvent Exposure Identify Cryptic Pockets, Switches, and Allosteric Coupling. Biophys. J. 116:818\u0026ndash;830\u003c/li\u003e\n\u003cli\u003eFrey BJ, Dueck D (2007) Clustering by Passing Messages Between Data Points. Science 315:972\u0026ndash;976. https://doi.org/10.1126/science.1136800\u003c/li\u003e\n\u003cli\u003eHerold J, Siddell SG, Gorbalenya AE (1999) A Human RNA Viral Cysteine Proteinase That Depends upon a Unique Zn2+-binding Finger Connecting the Two Domains of a Papain-like Fold *. J Biol Chem 274:14918\u0026ndash;14925. https://doi.org/10.1074/jbc.274.21.14918\u003c/li\u003e\n\u003cli\u003eKlemm T, Ebert G, Calleja DJ, et al (2020) Mechanism and inhibition of the papain‐like protease, PLpro, of SARS‐CoV‐2. EMBO J 39:e106275. https://doi.org/10.15252/embj.2020106275\u003c/li\u003e\n\u003cli\u003eEwert W, G\u0026uuml;nther S, Miglioli F, et al (2022) Hydrazones and Thiosemicarbazones Targeting Protein-Protein-Interactions of SARS-CoV-2 Papain-like Protease. Front Chem 10:. https://doi.org/10.3389/fchem.2022.832431\u003c/li\u003e\n\u003cli\u003eFerreira JC, Villanueva AJ, Adem KA, et al (2024) Identification of novel allosteric sites of SARS-CoV-2 papain-like protease (PLpro) for the development of COVID-19 antivirals. J Biol Chem 300:. https://doi.org/10.1016/j.jbc.2024.107821\u003c/li\u003e\n\u003cli\u003eVithani N, Zhang S, Thompson JP, et al (2024) Exploration of Cryptic Pockets Using Enhanced Sampling Along Normal Modes: A Case Study of KRAS G12D. J Chem Inf Model. https://doi.org/10.1021/acs.jcim.4c01435\u003c/li\u003e\n\u003c/ol\u003e"},{"header":"Supplemental Figures","content":"\u003cp\u003eSupplemental Figures S1-S5 are not available with this version.\u003c/p\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"journal-of-computer-aided-molecular-design","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"jcam","sideBox":"Learn more about [Journal of Computer-Aided Molecular Design](http://link.springer.com/journal/10822)","snPcode":"10822","submissionUrl":"https://submission.nature.com/new-submission/10822/3","title":"Journal of Computer-Aided Molecular Design","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"Springer Hybrid","inReviewEnabled":true,"inReviewRevisionsEnabled":false},"keywords":"Molecular Dynamics Simulations, Flooding Simulations, Binding Site Identification, Cosolvent Mapping, Protein-protein Interaction Sites, Sites of Allosteric Modulation","lastPublishedDoi":"10.21203/rs.3.rs-7499914/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-7499914/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eThe SARS-CoV-2 papain-like protease (PLpro) represents a crucial therapeutic target due to its dual role in viral polyprotein processing and suppression of host immune responses through de-ubiquitination and de-ISGylation activities. To identify novel allosteric druggable sites on PLpro, we developed a molecular dynamics approach with flooding fragments (MDFFr), advancing previous Molecular Dynamics (MD) flooding simulation methods. Using MDFFr, we evaluated interactions of known phenolic inhibitors with SARS-CoV-2 PLpro and discovered multiple allosteric druggable sites. Our simulations not only validated experimentally known binding sites but also revealed previously unidentified hotspots, including protein-protein interaction sites for ubiquitin and ISG-15 (Interferon-Stimulated Gene 15). The MDFFr approach demonstrates robust capability for physics-based druggability assessment of biological targets using only protein 3D structure, while providing detailed insights into fragment-protein interactions at both druggable sites and protein-protein interfaces. These findings unveil new opportunities for allosteric inhibition of PLpro, potentially advancing therapeutic strategies against SARS-CoV-2 and other coronavirus-related diseases.\u003c/p\u003e","manuscriptTitle":"Chasing allosteric inhibition of the SARS-CoV-2 PLpro via Molecular Dynamics Simulations with Flooding Fragments (MDFFr)","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-09-15 19:55:35","doi":"10.21203/rs.3.rs-7499914/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"decision","content":"Revision requested","date":"2025-10-06T01:22:31+00:00","index":"","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-10-05T22:35:13+00:00","index":"hide","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-09-23T00:44:43+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"162512526350332730846730284941833619245","date":"2025-09-18T13:24:10+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"281108622354171261379005865597830994278","date":"2025-09-14T07:36:43+00:00","index":"hide","fulltext":""},{"type":"reviewersInvited","content":"","date":"2025-09-12T23:51:30+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2025-09-12T23:50:20+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2025-09-02T13:56:46+00:00","index":"","fulltext":""},{"type":"submitted","content":"Journal of Computer-Aided Molecular Design","date":"2025-08-31T10:09:33+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"journal-of-computer-aided-molecular-design","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"jcam","sideBox":"Learn more about [Journal of Computer-Aided Molecular Design](http://link.springer.com/journal/10822)","snPcode":"10822","submissionUrl":"https://submission.nature.com/new-submission/10822/3","title":"Journal of Computer-Aided Molecular Design","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"Springer Hybrid","inReviewEnabled":true,"inReviewRevisionsEnabled":false}}],"origin":"","ownerIdentity":"b2f486b4-fdfa-45bc-9007-5f0be9d2af11","owner":[],"postedDate":"September 15th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[],"tags":[],"updatedAt":"2025-12-22T15:59:02+00:00","versionOfRecord":{"articleIdentity":"rs-7499914","link":"https://doi.org/10.1007/s10822-025-00730-0","journal":{"identity":"journal-of-computer-aided-molecular-design","isVorOnly":false,"title":"Journal of Computer-Aided Molecular Design"},"publishedOn":"2025-12-19 15:56:56","publishedOnDateReadable":"December 19th, 2025"},"versionCreatedAt":"2025-09-15 19:55:35","video":"","vorDoi":"10.1007/s10822-025-00730-0","vorDoiUrl":"https://doi.org/10.1007/s10822-025-00730-0","workflowStages":[]},"version":"v1","identity":"rs-7499914","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-7499914","identity":"rs-7499914","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}