Abstract Constant pH molecular dynamics (MD) simulations sample protonation states on the fly according to the conformational environment and user specified pH condition; however, the current accuracy is limited due to the use of implicit-solvent models or a hybrid solvent scheme. Here we report the first GPU-accelerated implementation, parameterization, and validation of the all-atom continuous constant pH MD (CpHMD) method with particle-mesh Ewald (PME) electrostatics in the Amber22 pmemd. cuda engine. The titration parameters for Asp, Glu, His, Cys, and Lys were derived for the CHARMM c22 and Amber ff14sb and ff19sb force fields. We then evaluated the PME-CpHMD method using the asynchronous pH replica-exchange titration simulations with the c22 force field for six benchmark proteins, including BBL, hen egg white lysozyme (HEWL), staphylococcal nuclease (SNase), thioredoxin, ribonuclease A (RNaseA), and human muscle creatine kinase (HMCK). The root-mean-square deviation from the experimental p K a ’s of Asp, Glu, His, and Cys is 0.76 pH units, and the Pearson’s correlation coefficient for the p K a shifts with respect to model values is 0.80. We demonstrated that a finite-size correction or much enlarged simulation box size can remove a systematic error of the calculated p K a ’s and improve agreement with experiment. Importantly, the simulations captured the relevant biology in several challenging cases, e.g., the titration order of the catalytic dyad Glu35/Asp52 in HEWL and the coupled residues Asp19/Asp21 in SNase, the large p K a upshift of the deeply buried catalytic Asp26 in thioredoxin, and the large p K a downshift of the deeply buried catalytic Cys283 in HMCK. We anticipate that PME-CpHMD offers proper pH control to improve the accuracies of MD simulations and enables mechanistic studies of proton-coupled dynamical processes that are ubiquitous in biology but remain poorly understood due to the lack of experimental tools and limitation of current MD simulations.

Broad-spectrum antiviral drugs are urgently needed to stop the COVID-19 pandemic and prevent future ones. The novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is related to SARS-CoV and Middle East respiratory syndrome coronavirus (MERS-CoV), which have caused the previous outbreaks. The papain-like protease (PLpro) is an attractive drug target due to its essential roles in the viral life cycle. As a cysteine protease, PLpro is rich in cysteines and histidines and their protonation/deprotonation modulates catalysis and conformational plasticity. Here we report the pKa calculations and assessment of the proton-coupled conformational dynamics of SARS-CoV-2 in comparison to SARS-CoV and MERS-CoV PLpros using a newly developed GPU-accelerated implicit-solvent continuous constant pH molecular dynamics method with an asynchronous replica-exchange scheme. The calculated pKa's support the catalytic roles of the Cys-His-Asp triad. We also found that several residues can switch protonation states at physiological pH, among which is C270/271 located on the flexible blocking loop 2 (BL2) of SARS-CoV-2/CoV PLpro. Simulations revealed that the BL2 conformational dynamics is coupled to the titration of C271/270, in agreement with the crystal structures of SARS-CoV-2 PLpro. Simulations also revealed that BL2 in MERS-CoV PLpro is very flexible, sampling both open and closed states despite the lack of an analogous cysteine. Our work provides a starting point for more detailed mechanistic studies to assist structure-based design of broad-spectrum inhibitors against CoV PLpros.

Abstract Many membrane channels, transporters, and receptors utilize a pH gradient or proton coupling to drive functionally relevant conformational transitions. Conventional molecular dynamics simulations employ fixed protonation states, thus neglecting the coupling between protonation and conformational equilibria. Here we describe the membrane-enabled hybrid-solvent continuous constant pH molecular dynamics method for capturing atomic details of proton-coupled conformational dynamics of transmembrane proteins. Example protocols from our recent application studies of proton channels and ion/substrate transporters are discussed.

During the continuing evolution of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the Omicron variant of concern emerged in the second half of 2021 and has been dominant since November of that year. Along with its sublineages, it has maintained a prominent role ever since. The Nsp5 main protease (Mpro) of the Omicron virus is characterized by a single dominant mutation, P132H. Here we determined the X-ray crystal structures of the P132H mutant (or O-Mpro) as a free enzyme and in complex with the Mpro inhibitor, the alpha-ketoamide 13b-K, and we conducted enzymological, biophysical as well as theoretical studies to characterize the O-Mpro. We found that O-Mpro has a similar overall structure and binding with 13b-K; however, it displays lower enzymatic activity and lower thermal stability compared to the WT-Mpro (with "WT" referring to the prototype strain). Intriguingly, the imidazole ring of His132 and the carboxylate plane of Glu240 are in a stacked configuration in the X-ray structures determined here. Empirical folding free energy calculations suggest that the O-Mpro dimer is destabilized relative to the WT-Mpro due to less favorable van der Waals interactions and backbone conformations in the individual protomers. All-atom continuous constant-pH molecular dynamics (MD) simulations reveal that His132 and Glu240 display coupled titration. At pH 7, His132 is predominantly neutral and in a stacked configuration with respect to Glu240 which is charged. In order to examine whether the Omicron mutation eases the emergence of further Mpro mutations, we also analyzed the P132H+T169S double mutant, which is characteristic of the BA.1.1.2 lineage. However, we found little evidence of a correlation between the two mutation sites.

ABSTRACT The molecular mechanisms governing the human voltage-gated proton channel hH v 1 remain elusive. Here we used membrane-enabled hybrid-solvent continuous constant pH molecular dynamics (CpHMD) simulations with pH replica exchange to further evaluate the recently obtained structural models of hH v 1 in hyperpolarized (closed channel) and depolarized (open channel) states (Geragotelis, Tobias et al., Proc. Natl. Acad. Sci. USA 2020) and explore potential pH-sensing residues. The CpHMD titration at a set of symmetric pH conditions revealed three residues that can gain or lose protons upon channel depolarization. Among them residue H168 at the intracellular end of the S3 helix switches from the deprotonated to the protonated state and its protonation is correlated with the increased tilting of the S3 helix during the transition from the hyper-to the depolarized state. Thus, the simulation data suggest H168 as an interior pH sensor, in support of a recent finding based on electrophysiological experiments of H v 1 mutants (Cherny, DeCoursey et al., J. Gen. Physiol. 2018). Our work represents an important step towards deciphering the pH-dependent gating mechanism of hH v 1. SIGNIFICANCE The human voltage-gated proton channel hH v 1 is comprised of a proton-selective voltage sensing domain and responsible for cellular pH homeostasis. Despite intense experimental and theoretical investigations, its pH-dependent gating mechanism is not understood. Our simulation data offer strong evidence supporting the role of H168 as a pH i sensor (Cherny, DeCoursey et al., J. Gen. Physiol. 2018). Deciphering the interior pH sensor moves us a step closer to elucidating the structure-function relationship of hH v 1.

Abstract Nirmatrelvir is an orally available inhibitor of SARS-CoV-2 main protease (Mpro) and the main ingredient of PAXLOVID, a drug approved by FDA for high-risk COVID-19 patients. Recently, a rare natural mutation, H172Y, was found to significantly reduce nirmatrelvir’s inhibitory activity. As the COVID-19 cases skyrocket in China and the selective pressure of antiviral therapy builds up in the US, there is an urgent need to characterize and understand how the H172Y mutation confers drug resistance. Here we investigated the H172Y Mpro’s conformational dynamics, folding stability, catalytic efficiency, and inhibitory activity using all-atom constant pH and fixed-charge molecular dynamics simulations, alchemical and empirical free energy calculations, artificial neural networks, and biochemical experiments. Our data suggests that the mutation significantly weakens the S1 pocket interactions with the N-terminus and perturbs the conformation of the oxyanion loop, leading to a decrease in the thermal stability and catalytic efficiency. Importantly, the perturbed S1 pocket dynamics weakens the nirma-trelvir binding in the P1 position, which explains the decreased inhibitory activity of nirmatrelvir. Our work demonstrates the predictive power of the combined simulation and artificial intel-ligence approaches, and together with biochemical experiments they can be used to actively surveil continually emerging mutations of SARS-CoV-2 Mpro and assist the discovery of new antiviral drugs. The presented workflow can be applicable to characterize mutation effects on any protein drug targets.

In 2019, drug overdose has claimed over 70,000 lives in the United States. More than half of the deaths are related to synthetic opioids represented by fentanyl which is a potent agonist of mu-opioid receptor (mOR). In recent years, the crystal structures of mOR in complex with morphine derivatives have been determined; however, structural basis of mOR activation by fentanyl-like synthetic opioids remains lacking. Exploiting the X-ray structure of mOR bound to a morphinan ligand and several state-of-the-art simulation techniques, including weighted ensemble and continuous constant pH molecular dynamics, we elucidated the detailed binding mechanism of fentanyl with mOR. Surprisingly, in addition to forming a salt-bridge with Asp147

All-atom constant pH molecular dynamics simulations offer a powerful tool for understanding pH-mediated and proton-coupled biological processes. As the protonation equilibria of protein sidechains are shifted by electrostatic interactions and desolvation energies, pKa values calculated from the constant pH simulations may be sensitive to the underlying protein force field and water model. Here we investigated the force field dependence of the all-atom particle mesh Ewald (PME) continuous constant pH (PME-CpHMD) simulations of a mini-protein BBL. The replica-exchange titration simulations based on the Amber ff19SB and ff14SB force fields with the respective water models showed significantly overestimated pKa downshifts for a buried histidine (His166) and for two glutamic acids (Glu141 and Glu161) that are involved in salt-bridge interactions. These errors (due to undersolvation of neutral histidines and overstabilization of salt bridges) are consistent with the previously reported \pka's based on the CHARMM c22/CMAP force field, albeit in larger magnitudes. The pKa calculations also demonstrated that ff19SB with OPC water is significantly more accurate than ff14SB with TIP3P water, and the salt-bridge related pKa downshifts can be partially alleviated by the atom-pair specific Lennard-Jones corrections (NBFIX). Together, these data suggest that the accuracies of the protonation equilibria of proteins from constant pH simulations can significantly benefit from improvements of force fields.

Machine learning (ML) identification of covalently ligandable sites may significantly accelerate targeted covalent inhibitor discoveries and expand the druggable proteome space. Here we report the development of the tree-based models and convolutional neural networks trained on a newly curated database (LigCys3D) of over 1,000 liganded cysteines in nearly 800 proteins represented by over 10,000 X-ray structures as reported in the protein data bank (PDB). The unseen tests yielded 94% AUC (area under the receiver operating characteristic curve), demonstrating the highly predictive power of the models. Interestingly, application to the proteins evaluated by the activity-based protein profiling (ABPP) experiments in cell lines gave a lower AUC of 72%. Analysis revealed significant discrepancies in the structural environment of the ligandable cysteines captured by X-ray crystallography and those determined by ABPP. This surprising finding warrants further investigations and may have implications for future drug discoveries. We discuss ways to improve the models and project future directions. Our work represents a first step towards the ML-led integration of big genome data, structure models, and chemoproteomic experiments to annotate the human proteome space for the next-generation drug discoveries.

Abstract Renin is a pepsin-like aspartyl protease and an important drug target for the treatment of hypertension; despite three decades’ research, its pH-dependent structure-function relationship remains poorly understood. Here we employed the continuous constant pH molecular dynamics (CpHMD) simulations to decipher the acid/base roles of renin’s catalytic dyad and the conformational dynamics of the flap, which is a common structural feature among aspartyl proteases. The calculated p K a ’s suggest that the catalytic Asp38 and Asp226 serve as the general base and acid, respectively, in agreement with experiment and supporting the hypothesis that renin’s neutral optimum pH is due to the substrate-induced p K a shifts of the aspartic dyad. The CpHMD data confirmed our previous hypothesis that hydrogen bond formation is the major determinant of the dyad p K a order. Additionally, our simulations showed that renin’s flap remains open regardless of pH, although a Tyr-inhibited state is occasionally formed above pH 5. These findings are discussed in comparison to the related aspartyl proteases, including β -secretases 1 and 2, capthepsin D, and plasmepsin II. Our work represents a first step towards a systematic understanding of the pH-dependent structure-dynamics-function relationships of pepsin-like aspartyl proteases that play important roles in biology and human disease states.