On Inactivation of the Coronavirus Main Protease.

A deeper understanding of the inactive conformations of the coronavirus main protease (MPro) could inform the design of allosteric drugs. Based on extensive molecular dynamics simulations, we built a Markov State Model to investigate structural changes that can inactivate the SARS-CoV-2 MPro. In a subset of structures, one subunit of the homodimer assumes an inactive conformation that resembles an inactive crystal structure. However, contradicting the widely held half-of-sites activity hypothesis, the most populated enzyme structures have two active subunits. We then used transition path theory (TPT) and the Jensen-Shannon Divergence (JSD) to pinpoint residues involved in the inactivation process. A π stack between Phe140 and His163 is a key feature that can distinguish active and inactive conformations of MPro. Each subunit has unique inactive conformations stabilized by π stacking interactions involving residues Phe140, Tyr118, His163, and His172, a hydrogen bonding network centered around His163 and His172, and a modified network of interactions in the dimer interface. The importance of these residues in maintaining an active structure explains the sensitivity of enzymatic activity to site-directed mutagenesis.

Learned mappings for targeted free energy perturbation between peptide conformations.

Targeted free energy perturbation uses an invertible mapping to promote configuration space overlap and the convergence of free energy estimates. However, developing suitable mappings can be challenging. Wirnsberger et al. [J. Chem. Phys. 153, 144112 (2020)] demonstrated the use of machine learning to train deep neural networks that map between Boltzmann distributions for different thermodynamic states. Here, we adapt their approach to the free energy differences of a flexible bonded molecule, deca-alanine, with harmonic biases and different spring centers. When the neural network is trained until "early stopping"-when the loss value of the test set increases-we calculate accurate free energy differences between thermodynamic states with spring centers separated by 1 Å and sometimes 2 Å. For more distant thermodynamic states, the mapping does not produce structures representative of the target state, and the method does not reproduce reference calculations.

Bayesian Regression Quantifies Uncertainty of Binding Parameters from Isothermal Titration Calorimetry More Accurately Than Error Propagation.

We compare several different methods to quantify the uncertainty of binding parameters estimated from isothermal titration calorimetry data: the asymptotic standard error from maximum likelihood estimation, error propagation based on a first-order Taylor series expansion, and the Bayesian credible interval. When the methods are applied to simulated experiments and to measurements of Mg(II) binding to EDTA, the asymptotic standard error underestimates the uncertainty in the free energy and enthalpy of binding. Error propagation overestimates the uncertainty for both quantities, except in the simulations, where it underestimates the uncertainty of enthalpy for confidence intervals less than 70%. In both datasets, Bayesian credible intervals are much closer to observed confidence intervals.

Inclusion of Control Data in Fits to Concentration-Response Curves Improves Estimates of Half-Maximal Concentrations.

Concentration-response curves, in which the effect of varying the concentration on the response of an assay is measured, are widely used to evaluate biological effects of chemical compounds. While National Center for Advancing Translational Sciences guidelines specify that readouts should be normalized by the controls, recommended statistical analyses do not explicitly fit to the control data. Here, we introduce a nonlinear regression procedure based on maximum likelihood estimation that determines parameters for the classical Hill equation by fitting the model to both the curve and the control data. Simulations show that the proposed procedure provides more precise parameters compared with previously prescribed practices. Analysis of enzymatic inhibition data from the COVID Moonshot demonstrates that the proposed procedure yields a lower asymptotic standard error for estimated parameters. Benefits are most evident in the analysis of the incomplete curves. We also find that Lenth's outlier detection method appears to determine parameters more precisely.

H-Bonds in Crambin: Coherence in an α-Helix.

We applied coherence analysis-used by engineers to identify linear interactions in stochastic systems-to molecular dynamics simulations of crambin, a thionin storage protein found in Abyssinian cabbage. A key advantage of coherence over other analyses is that it is robust, independent of the properties, or even the existence of probability distributions often relied on in statistical mechanics. For frequencies between 0.391 and 5.08 GHz (corresponding reciprocally to times of 2.56 and 0.197 ns), the displacements of oxygen and nitrogen atoms across α-helix H-bonds are strongly correlated, with a coherence greater than 0.9; the secondary structure causes the H-bonds to effectively act as a spring. Similar coherence behavior is observed for covalent bonds and other noncovalent interactions including H-bonds in ß-sheets and salt bridges. In contrast, arbitrary pairs of atoms that are physically distant have uncorrelated motions and negligible coherence. These results suggest that coherence may be used to objectively identify atomic interactions without subjective thresholds such as H-bond lengths angles and angles. Strong coherence is also observed between the average position of adjacent leaves (groups of atoms) in an α-helix, suggesting that the harmonic analysis of classical molecular dynamics can successfully describe the propagation of allosteric interactions through the structure.

Bayesian regression and model selection for isothermal titration calorimetry with enantiomeric mixtures.

Bayesian regression is performed to infer parameters of thermodynamic binding models from isothermal titration calorimetry measurements in which the titrant is an enantiomeric mixture. For some measurements the posterior density is multimodal, indicating that additional data with a different protocol are required to uniquely determine the parameters. Models of increasing complexity-two-component binding, racemic mixture, and enantiomeric mixture-are compared using model selection criteria. To precisely estimate one of these criteria, the Bayes factor, a variation of bridge sampling is developed.

Identification of the riboflavin cofactor-binding site in the Vibrio cholerae ion-pumping NQR complex: A novel structural motif in redox enzymes.

The ion-pumping NQR complex is an essential respiratory enzyme in the physiology of many pathogenic bacteria. This enzyme transfers electrons from NADH to ubiquinone through several cofactors, including riboflavin (vitamin B2). NQR is the only enzyme reported that is able to use riboflavin as a cofactor. Moreover, the riboflavin molecule is found as a stable neutral semiquinone radical. The otherwise highly reactive unpaired electron is stabilized via an unknown mechanism. Crystallographic data suggested that riboflavin might be found in a superficially located site in the interface of NQR subunits B and E. However, this location is highly problematic, as the site does not have the expected physiochemical properties. In this work, we have located the riboflavin-binding site in an amphipathic pocket in subunit B, previously proposed to be the entry site of sodium. Here, we show that this site contains absolutely conserved residues, including N200, N203, and D346. Mutations of these residues decrease enzymatic activity and specifically block the ability of NQR to bind riboflavin. Docking analysis and molecular dynamics simulations indicate that these residues participate directly in riboflavin binding, establishing hydrogen bonds that stabilize the cofactor in the site. We conclude that riboflavin is likely bound in the proposed pocket, which is consistent with enzymatic characterizations, thermodynamic studies, and distance between cofactors.

Electrostatics and water occlusion regulate covalently-bound flavin mononucleotide cofactors of Vibrio cholerae respiratory complex NQR.

Proteins like NADH:ubiquinone oxidoreductase (NQR), an essential enzyme and ion pump in the physiology of several pathogenic bacteria, tightly regulate the redox properties of their cofactors. Although flavin mononucleotide (FMN) is fully reduced in aqueous solution, FMN in subunits B and C of NQR exclusively undergo one-electron transitions during its catalytic cycle. Here, we perform ab initio calculations and molecular dynamics simulations to elucidate the mechanisms that regulate the redox state of FMN in NQR. QM/MM calculations show that binding site electrostatics disfavor anionic forms of FMNH2 , but permit a neutral form of the fully reduced flavin. The potential energy surface is unaffected by covalent bonding between FMN and threonine. Molecular dynamics simulations show that the FMN binding sites are inaccessible by water, suggesting that further reductions of the cofactors are limited or prohibited by the availability of water and other proton donors. These findings provide a deeper understanding of the mechanisms used by NQR to regulate electron transfer through the cofactors and perform its physiologic role. They also provide the first, to our knowledge, evidence of the simple concept that proteins regulate flavin redox states via water occlusion.

A Structural Model for Bax∆2-Mediated Activation of Caspase 8-Dependent Apoptosis.

Bax∆2 is a pro-apoptotic anti-tumor protein in the Bax family. While most of the Bax family causes cell death by targeting mitochondria, Bax∆2 forms cytosolic aggregates and activates caspase 8-dependent cell death. We previously showed that the Bax∆2 helix α9 is critical for caspase 8 recruitment. However, the interaction between these two proteins at the structural level is unknown. In this in silico study, we performed molecular dynamics (MD) simulations and protein-protein docking on Bax∆2 variants. The results suggest that the Bax∆2 variants have different stable states. Mutating the Baxα mitochondria-targeting signal [L26P/L27P] appears to introduce a kink into helix α1. Protein-protein docking suggests that helices α9 of both wild-type Bax∆2 and Bax∆2 caspase 8 binding-deficient mutant [L164P] can fit in the same caspase 8 binding site, but the mutant is unable to fit as well as wild-type Bax∆2. Together, these data point to a structural basis for explaining Bax∆2 function in caspase 8-dependent cell death.

On the polarization of ligands by proteins.

Although ligand-binding sites in many proteins contain a high number density of charged side chains that can polarize small organic molecules and influence binding, the magnitude of this effect has not been studied in many systems. Here, we use a quantum mechanics/molecular mechanics (QM/MM) approach, in which the ligand is the QM region, to compute the ligand polarization energy of 286 protein-ligand complexes from the PDBBind Core Set (release 2016). Calculations were performed both with and without implicit solvent based on the domain decomposition Conductor-like Screening Model. We observe that the ligand polarization energy is linearly correlated with the magnitude of the electric field acting on the ligand, the magnitude of the induced dipole moment, and the classical polarization energy. The influence of protein and cation charges on the ligand polarization diminishes with the distance and is below 2 kcal mol-1 at 9 Å and 1 kcal mol-1 at 12 Å. Compared to these embedding field charges, implicit solvent has a relatively minor effect on ligand polarization. Considering both polarization and solvation appears essential to computing negative binding energies in some crystallographic complexes. Solvation, but not polarization, is essential for achieving moderate correlation with experimental binding free energies.

Robosample: A rigid-body molecular simulation program based on robot mechanics.

BACKGROUND: Compared with all-atom molecular dynamics (MD), constrained MD methods allow for larger time steps, potentially reducing computational cost. For this reason, there has been continued interest in improving constrained MD algorithms to increase configuration space sampling in molecular simulations. METHODS: Here, we introduce Robosample, a software package that implements high-performance constrained dynamics algorithms, originally developed for robotics, and applies them to simulations of biomolecular systems. As in the gMolmodel package developed by Spiridon and Minh in 2017, Robosample uses Constrained Dynamics Hamiltonian Monte Carlo (CDHMC) as a Gibbs sampling move - a type of Monte Carlo move where a subset of coordinates is allowed to change. In addition to the previously described Cartesian and torsional dynamics moves, Robosample implements spherical and cylindrical joints that can be distributed along the molecule by the user. RESULTS: In alanine dipeptide simulations, the free energy surface is recovered by mixing fully flexible with torsional, cylindrical, or spherical dynamics moves. Ramachandran dynamics, where only the two key torsions are mobile, accelerate the slowest transition by an order of magnitude. We also show that simulations of a complex glycan cover significantly larger regions of the configuration space when mixed with constrained dynamics. MAJOR CONCLUSIONS: Robosample is a tool of choice for efficient conformational sampling of large biomolecules. GENERAL SIGNIFICANCE: Robosample is intended as a reliable and user-friendly simulation package for fast biomolecular sampling that does not require extensive expertise in mechanical engineering or in the statistical mechanics of reduced coordinates.

Implicit ligand theory for relative binding free energies: II. An estimator based on control variates.

Implicit ligand theory describes the relationship between the noncovalent binding free energy and the binding free energy between a ligand and multiple rigid receptor conformations. We have previously shown that if the receptor conformations are sampled from or reweighed to a holo ensemble, the binding free energy relative to the ligand that defines the ensemble can be calculated. Here, we apply a variance reduction technique known as control variates to derive a new statistical estimator for the relative binding free energy. In applications to a data set of 6 reference ligands and 18 test ligands, statistically significant differences between the estimators are not observed for most systems. However, in cases where such differences are observed, the new estimator is more accurate, precise, and converges more quickly. Performance improvements are most consistent where there is a clear correlation, with a correlation coefficient greater than 0.3, between the control variate and the statistic being averaged.

On Restraints in End-Point Protein-Ligand Binding Free Energy Calculations.

The impact of harmonic restraints on protein heavy atoms and ligand atoms on end-point free energy calculations is systematically characterized for 54 protein-ligand complexes. We observe that stronger restraints reduce the equilibration time and statistical inefficiency, suppress conformational sampling, influence correlation with experiment, and monotonically decrease the estimated loss of entropy upon binding, leading to stronger estimated binding free energies in most systems. A statistical estimator that reweights for the biasing potential and includes data prior to the estimated equilibration time has the highest correlation with experiment. A spring constant of 20 cal mol-1 Å-2 maintains a near-native energy landscape and suppresses artifactual energy minima while minimally limiting thermal fluctuations about the crystal structure. © 2019 Wiley Periodicals, Inc.

Alchemical Grid Dock (AlGDock): Binding Free Energy Calculations between Flexible Ligands and Rigid Receptors.

Alchemical Grid Dock (AlGDock) is open-source software designed to compute the binding potential of mean force-the binding free energy between a flexible ligand and a rigid receptor-for a small organic ligand and a biological macromolecule. Multiple BPMFs can be used to rigorously compute binding affinities between flexible partners. AlGDock uses replica exchange between thermodynamic states at different temperatures and receptor-ligand interaction strengths. Receptor-ligand interaction energies are represented by interpolating precomputed grids. Thermodynamic states are adaptively initialized and adjusted on-the-fly to maintain adequate replica exchange rates. In demonstrative calculations, when the bound ligand is treated as fully solvated, AlGDock estimates BPMFs with a precision within 4 kT in 65% and within 8 kT for 91% of systems. It correctly identifies the native binding pose in 83% of simulations. Performance is sometimes limited by subtle differences in the important configuration space of sampled and targeted thermodynamic states. © 2019 Wiley Periodicals, Inc.

Role of Subunit D in Ubiquinone-Binding Site of Vibrio cholerae NQR: Pocket Flexibility and Inhibitor Resistance.

The ion-pumping NADH: ubiquinone dehydrogenase (NQR) is a vital component of the respiratory chain of numerous species of marine and pathogenic bacteria, including Vibrio cholerae. This respiratory enzyme couples the transfer of electrons from NADH to ubiquinone (UQ) to the pumping of ions across the plasma membrane, producing a gradient that sustains multiple homeostatic processes. The binding site of UQ within the enzyme is an important functional and structural motif that could be used to design drugs against pathogenic bacteria. Our group recently located the UQ site in the interface between subunits B and D and identified the residues within subunit B that are important for UQ binding. In this study, we carried out alanine scanning mutagenesis of amino acid residues located in subunit D of V. cholerae NQR to understand their role in UQ binding and enzymatic catalysis. Moreover, molecular docking calculations were performed to characterize the structure of the site at the atomic level. The results show that mutations in these positions, in particular, in residues P185, L190, and F193, decrease the turnover rate and increase the Km for UQ. These mutants also showed an increase in the resistance against the inhibitor HQNO. The data indicate that residues in subunit D fulfill important structural roles, restricting and orienting UQ in a catalytically favorable position. In addition, mutations of these residues open the site and allow the simultaneous binding of substrate and inhibitors, producing partial inhibition, which appears to be a strategy used by Pseudomonas aeruginosa to avoid autopoisoning.

Nonequilibrium path-ensemble averages for symmetric protocols.

According to the nonequilibrium work relations, path-ensembles generated by irreversible processes in which a system is driven out of equilibrium according to a predetermined protocol may be used to compute equilibrium free energy differences and expectation values. Estimation has previously been improved by considering data collected from the reverse process, which starts in equilibrium in the final thermodynamic state of the forward process and is driven according to the time-reversed protocol. Here, we develop a theoretically rigorous statistical estimator for nonequilibrium path-ensemble averages specialized for symmetric protocols, in which forward and reverse processes are identical. The estimator is tested with a number of model systems: a symmetric 1D potential, an asymmetric 1D potential, the unfolding of deca-alanine, separating a host-guest system, and translocating a potassium ion through a gramicidin A ion channel. When reconstructing free energies using data from symmetric protocols, the new estimator outperforms existing rigorous unidirectional and bidirectional estimators, converging more quickly and resulting in a smaller error. However, in most cases, using the bidirectional estimator with data from a forward and reverse pair of asymmetric protocols outperforms the corresponding symmetric protocol and estimator with the same amount of simulation time. Hence, the new estimator is only recommended when the bidirectional estimator is not feasible or is expected to perform poorly. The symmetric estimator shows similar performance to a unidirectional protocol of half the length and twice the number of trajectories.

Conserved residue His-257 of Vibrio cholerae flavin transferase ApbE plays a critical role in substrate binding and catalysis.

The flavin transferase ApbE plays essential roles in bacterial physiology, covalently incorporating FMN cofactors into numerous respiratory enzymes that use the integrated cofactors as electron carriers. In this work we performed a detailed kinetic and structural characterization of Vibrio cholerae WT ApbE and mutants of the conserved residue His-257, to understand its role in substrate binding and in the catalytic mechanism of this family. Bi-substrate kinetic experiments revealed that ApbE follows a random Bi Bi sequential kinetic mechanism, in which a ternary complex is formed, indicating that both substrates must be bound to the enzyme for the reaction to proceed. Steady-state kinetic analyses show that the turnover rates of His-257 mutants are significantly smaller than those of WT ApbE, and have increased Km values for both substrates, indicating that the His-257 residue plays important roles in catalysis and in enzyme-substrate complex formation. Analyses of the pH dependence of ApbE activity indicate that the pKa of the catalytic residue (pKES1) increases by 2 pH units in the His-257 mutants, suggesting that this residue plays a role in substrate deprotonation. The crystal structures of WT ApbE and an H257G mutant were determined at 1.61 and 1.92 Å resolutions, revealing that His-257 is located in the catalytic site and that the substitution does not produce major conformational changes. We propose a reaction mechanism in which His-257 acts as a general base that deprotonates the acceptor residue, which subsequently performs a nucleophilic attack on FAD for flavin transfer.

Alchemical Grid Dock (AlGDock) calculations in the D3R Grand Challenge 3 : Binding free energies between flexible ligands and rigid receptors.

We participated in Subchallenges 1 and 2 of the Drug Design Data Resource (D3R) Grand Challenge 3. To prepare our submissions, we performed molecular docking with UCSF DOCK 6 and binding potential of mean force (BPMF) calculations-free energy calculations between flexible ligands and rigid receptors-using our open-source software package Alchemical Grid Dock (AlGDock). For each system, submissions were based on the minimum BPMF calculated for a selected set of crystal structures. In Subchallenge 1, our workflow performed poorly. Possible reasons for the poor performance include the neglect of cooperative ligands and limited sampling of ligand binding poses. In Subchallenge 2, our workflow led to some of most highly correlated submissions (Pearson R = 0.5) for vascular endothelial growth factor receptor 2. However, our results were poorly correlated for Janus Kinase 2 and Mitogen-activated protein kinase 14. Affinity prediction could potentially be improved by systematic selection of more diverse receptor configurations.

Design, Synthesis, and Biological Evaluation of Polyaminocarboxylate Ligand-Based Theranostic Conjugates for Antibody-Targeted Cancer Therapy and Near-Infrared Optical Imaging.

We report the design, synthesis, and evaluation of polyaminocarboxylate ligand-based antibody conjugates for potential application in targeted cancer therapy and near-infrared (NIR) fluorescence imaging. We synthesized a new polyaminocarboxylate chelate (CAB-NE3TA) as a potential anticancer agent. CAB-NE3TA displayed potent inhibitory activities against various cancer cell lines. We then designed a multifunctional theranostic platform (CAB-NE3TA-PAN-IR800) constructed on an epidermal growth factor receptor (EGFR)-targeted antibody (panitumumab, PAN) labeled with a NIR fluorescent dye. We also built the first atomistic model of the EGFR-PAN complex and loaded it with the cytotoxic CAB-NE3TA and the NIR dye. The therapeutic (CAB-NE3TA-PAN) and theranostic (CAB-NE3TA-PAN-IR800) conjugates were evaluated using an EGFR-positive A431 (human skin cancer) cell xenograft mouse model. Biodistribution studies using NIR fluorescence imaging demonstrated that the CAB-NE3TA-PAN labeled with the IR800 dye selectively targeted the A431 tumors in mice and resulted in prolonged retention in the tumor tissue and displayed excellent clearance in blood and normal organs. The therapeutic conjugate was capable of significantly inhibiting tumor growth, leading to nearly complete disappearance of tumors in the mice. The results of our pilot in vivo studies support further evaluation of the novel ligand-based therapeutic and theranostic conjugates for targeted iron chelation cancer therapy and imaging applications.