Mechanosensitive Piezo1 and Piezo2 channels transduce various forms of mechanical forces into cellular signals that play vital roles in many important biological processes in vertebrate organisms. Besides mechanical forces, Piezo1 is selectively activated by micromolar concentrations of the small molecule Yoda1 through an unknown mechanism. Here, using a combination of all-atom molecular dynamics simulations, calcium imaging and electrophysiology, we identify an allosteric Yoda1 binding pocket located in the putative mechanosensory domain, approximately 40 Å away from the central pore. Our simulations further indicate that the presence of the agonist correlates with increased tension-induced motions of the Yoda1-bound subunit. Our results suggest a model wherein Yoda1 acts as a molecular wedge, facilitating force-induced conformational changes, effectively lowering the channel’s mechanical threshold for activation. The identification of an allosteric agonist binding site in Piezo1 channels will pave the way for the rational design of future Piezo modulators with clinical value.

Piezo proteins are transmembrane ion channels which transduce many forms of mechanical stimuli into electrochemical signals. Their pore, formed by the assembly of three identical subunits, opens by an unknown mechanism. Here, to probe this mechanism, we investigate the interaction of Piezo1 with the small molecule agonist Yoda1. By engineering chimeras between mouse Piezo1 and its Yoda1-insensitive paralog Piezo2, we first identify a minimal protein region required for Yoda1 sensitivity. We next study the effect of Yoda1 on heterotrimeric Piezo1 channels harboring wild type subunits and Yoda1-insensitive mutant subunits. Using calcium imaging and patch-clamp electrophysiology, we show that hybrid channels harboring as few as one Yoda1-sensitive subunit exhibit Yoda1 sensitivity undistinguishable from homotrimeric wild type channels. Our results show that the Piezo1 pore remains fully open if only one subunit remains activated. This study sheds light on the gating and pharmacological mechanisms of a member of the Piezo channel family.

Mechanosensitive Piezo1 channels are essential mechanotransduction proteins in eukaryotes. Their curved transmembrane domains, called arms, create a convex membrane deformation, or footprint, which is predicted to flatten in response to increased membrane tension. Here, using a hyperbolic tangent model, we show that, due to the intrinsic bending rigidity of the membrane, the overlap of neighboring Piezo1 footprints produces a flattening of the Piezo1 footprints and arms. Multiple all-atom molecular dynamics simulations of Piezo1 further reveal that this tension-independent flattening is accompanied by gating motions that open an activation gate in the pore. This open state recapitulates experimentally obtained ionic selectivity, unitary conductance, and mutant phenotypes. Tracking ion permeation along the open pore reveals the presence of intracellular and extracellular fenestrations acting as cation-selective sites. Simulations also reveal multiple potential binding sites for phosphatidylinositol 4,5-bisphosphate. We propose that the overlap of Piezo channel footprints may act as a cooperative mechanism to regulate channel activity.

Reversible covalent inhibitors have many clinical advantages over noncovalent or covalent drugs. However, apart from selecting a warhead, substantial efforts in design and synthesis are needed to optimize noncovalent interactions to improve target-selective binding. Computational prediction of binding affinity for reversible covalent inhibitors presents a unique challenge since the binding process consists of multiple steps, which are not necessarily independent of each other. In this study, we lay out the relation between relative binding free energy and the overall reversible covalent binding affinity using a two-state binding model. To prove the concept, we employed free energy perturbation (FEP) coupled with λ-exchange molecular dynamics method to calculate the binding free energy of a series of α-ketoamide analogs relative to a common warhead scaffold, in both noncovalent and covalent bond states, and for two highly homologous proteases, calpain-1 and calpain-2. We conclude that covalent binding affinity alone, in general, can be used to predict reversible covalent binding selectivity. However, exceptions may exist. Therefore, we also discuss the conditions under which the noncovalent binding step is no longer negligible and propose a novel approach that combines the relative FEP calculations with a single QM/MM calculation of warhead to predict the binding affinity and binding kinetics for a large number of reversible covalent inhibitors. Our FEP calculations also revealed that covalent and noncovalent states of an inhibitor do not necessarily exhibit the same selectivity. Thus, investigating both binding states, as well as the kinetics will provide extremely useful information for optimizing reversible covalent inhibitors.

Continuum solvent models have been widely used in biomolecular modeling applications. Recently much attention has been given to inclusion of implicit membrane into existing continuum Poisson-Boltzmann solvent models to extend their applications to membrane systems. Inclusion of an implicit membrane complicates numerical solutions of the underlining Poisson-Boltzmann equation due to the dielectric inhomogeneity on the boundary surfaces of a computation grid. This can be alleviated by the use of the periodic boundary condition, a common practice in electrostatic computations in particle simulations. The conjugate gradient and successive over-relaxation methods are relatively straightforward to be adapted to periodic calculations, but their convergence rates are quite low, limiting their applications to free energy simulations that require a large number of conformations to be processed. To accelerate convergence, the Incomplete Cholesky preconditioning and the geometric multi-grid methods have been extended to incorporate periodicity for biomolecular applications. Impressive convergence behaviors were found as in the previous applications of these numerical methods to tested biomolecules and MMPBSA calculations.

Membrane proteins, due to their roles as cell receptors and signaling mediators, make prime candidates for drug targets. The computational analysis of protein-ligand binding affinities has been widely employed as a tool in rational drug design efforts. Although efficient implicit solvent-based methods for modeling globular protein-ligand binding have been around for many years, the extension of such methods to membrane protein-ligand binding is still in its infancy. In this study, we extended the widely used Amber/MMPBSA method to model membrane protein-ligand systems, and we used it to analyze protein-ligand binding for the human purinergic platelet receptor (P2Y12R), a prominent drug target in the inhibition of platelet aggregation for the prevention of myocardial infarction and stroke. The binding affinities, computed by the Amber/MMPBSA method using standard parameters, correlate well with experiment. A detailed investigation of these parameters was conducted to assess their impact on the accuracy of the method. These analyses show the importance of properly treating the non-polar solvation interactions and the electrostatic polarization in the binding of nucleotide agonists and non-nucleotide antagonists to P2Y12R. Based on the crystal structures and the experimental conditions in the binding assay, we further hypothesized that the nucleotide agonists lose their bound magnesium ion upon binding to P2Y12R, and our computational study supports this hypothesis. Ultimately, this work illustrates the value of computational analysis in the interpretation of experimental binding reactions.

To understand how protein function changes upon an allosteric perturbation, such as ligand binding and mutation, significant progress in characterizing allosteric network from molecular dynamics (MD) simulations has been made. However, determining which amino acid(s) play an essential role in the propagation of signals may prove challenging, even when the location of the source and sink is known for a protein or protein complex. This challenge is mainly due to the large fluctuations in protein dynamics that cause instability of the network topology within a single trajectory or between multiple replicas. To solve this problem, we introduce the current-flow betweenness scheme, originated from electrical network theory, to protein dynamical network analysis. To demonstrate the benefit of this new method, we chose a prototypic allosteric enzyme (IGPS or HisH-HisF dimer) as our benchmark system. Using multiple replicas of simulations and multiple network topology comparison metrics (edge ranking, path length, and node frequency), we show that the current-flow betweenness provides a significant improvement in the convergence of the allosteric networks. The improved stability of the network topology allows us to generate a delta-network between the apo and holo forms of the protein. We illustrated that the delta-network is a more rigorous way to capture the subtle changes in the networks that would otherwise be neglected by comparing node usage frequencies alone. We have also investigated the use of a linear smoothing function to improve the stability of the contact map. The methodology presented here is general and may be applied to other topology and weighting schemes. We thus conclude that, for determining protein signaling pathways between the pair(s) of source and sink, multiple MD simulation replicas are necessary, and the current-flow betweenness scheme introduced here provides a more robust approach than the geodesic scheme based on correlation edge weighting.

Electrostatic interactions play crucial roles in biophysical processes such as protein folding and molecular recognition. Poisson-Boltzmann equation (PBE)-based models have emerged as widely used in modeling these important processes. Though great efforts have been put into developing efficient PBE numerical models, challenges still remain due to the high dimensionality of typical biomolecular systems. In this study, we implemented and analyzed commonly used linear PBE solvers for the ever-improving graphics processing units (GPU) for biomolecular simulations, including both standard and preconditioned conjugate gradient (CG) solvers with several alternative preconditioners. Our implementation utilizes standard Nvidia® CUDA libraries cuSPARSE, cuBLAS, and CUSP. Extensive tests show that good numerical accuracy can be achieved given that the single precision is often used for numerical applications on GPU platforms. The optimal GPU performance was observed with the Jacobi-preconditioned CG solver, with a significant speedup over standard CG solver on CPU in our diversified test cases. Our analysis further shows that different matrix storage formats also considerably affect the efficiency of different linear PBE solvers on GPU, with the diagonal format best suited for our standard finite-difference linear systems. Further efficiency may be possible with matrix-free operations and integrated grid stencil setup specifically tailored for the banded matrices in PBE-specific linear systems.