Structural Fluctuation, Relaxation, and Folding of Protein: An Approach Based on the Combined Generalized Langevin and RISM/3D-RISM Theories.

In 2012, Kim and Hirata derived two generalized Langevin equations (GLEs) for a biomolecule in water, one for the structural fluctuation of the biomolecule and the other for the density fluctuation of water, by projecting all the mechanical variables in phase space onto the two dynamic variables: the structural fluctuation defined by the displacement of atoms from their equilibrium positions, and the solvent density fluctuation. The equation has an expression similar to the classical Langevin equation (CLE) for a harmonic oscillator, possessing terms corresponding to the restoring force proportional to the structural fluctuation, as well as the frictional and random forces. However, there is a distinct difference between the two expressions that touches on the essential physics of the structural fluctuation, that is, the force constant, or Hessian, in the restoring force. In the CLE, this is given by the second derivative of the potential energy among atoms in a protein. So, the quadratic nature or the harmonicity is only valid at the minimum of the potential surface. On the contrary, the linearity of the restoring force in the GLE originates from the projection of the water's degrees of freedom onto the protein's degrees of freedom. Taking this into consideration, Kim and Hirata proposed an ansatz for the Hessian matrix. The ansatz is used to equate the Hessian matrix with the second derivative of the free-energy surface or the potential of the mean force of a protein in water, defined by the sum of the potential energy among atoms in a protein and the solvation free energy. Since the free energy can be calculated from the molecular mechanics and the RISM/3D-RISM theory, one can perform an analysis similar to the normal mode analysis (NMA) just by diagonalizing the Hessian matrix of the free energy. This method is referred to as the Generalized Langevin Mode Analysis (GLMA). This theory may be realized to explore a variety of biophysical processes, including protein folding, spectroscopy, and chemical reactions. The present article is devoted to reviewing the development of this theory, and to providing perspective in exploring life phenomena.

A theory of chemical reactions in biomolecules in solution: Generalized Langevin mode analysis (GLMA).

The generalized Langevin mode analysis (GLMA) is applied to chemical reactions in biomolecules in solution. The theory sees a chemical reaction in solution as a barrier-crossing process, similar to the Marcus theory. The barrier is defined as the crossing point of two free-energy surfaces that are attributed to the reactant and product of the reaction. It is assumed that both free-energy surfaces are quadratic or harmonic. The assumption is based on the Kim-Hirata theory of structural fluctuation of protein, which proves that the fluctuation around an equilibrium structure is quadratic with respect to the structure or atomic coordinates. The quadratic surface is a composite of many harmonic functions with different modes or frequencies. The height of the activation barrier will be dependent on the mode or frequency-the less the frequency, the lower the barrier. Hence, it is essential to decouple the fluctuational modes into a hierarchical order. GLMA is impeccable for this purpose. It is essential for a theoretical study of chemical reactions to choose a reaction coordinate along which the reaction proceeds. We suppose that the mode whose center of coordinate and/or the frequency changes most before and after the reaction is the one relevant to the chemical reaction and choose the coordinate as the reaction coordinate. The rate of reaction along the reaction coordinate is krate=ν⁡exp-ΔF(†)/kBT, which is similar to the Marcus expression for the electron transfer reaction. In the equation, ΔF(†) is the activation barrier defined by ΔF(†)≡F(r)Q†-F(r)(Qeq (r)), where F(r)(Qeq (r)) and F(r)Q† denote the free energies at equilibrium Qeq (r) and the crossing point Q†, respectively, both on the free energy surface of the reactant.

Realization of the structural fluctuation of biomolecules in solution: Generalized Langevin mode analysis.

A new theoretical method, referred to as Generalized Langevin Mode Analysis (GLMA), is proposed to analyze the mode of structural fluctuations of a biomolecule in solution. The method combines the two theories in the statistical mechanics, or the Generalized Langevin theory and the RISM/3D-RISM theory, to calculate the second derivative, or the Hessian matrix, of the free energy surface of a biomolecule in aqueous solution, which consists of the intramolecular interaction among atoms in the biomolecule and the solvation free energy. The method is applied to calculate the wave-number spectrum of an alanine dipeptide in water for which the optical heterodyne-detected Raman-induced spectroscopy (RIKES) spectrum is available to compare with. The theoretical analysis reproduced the main features of the experimental spectrum with respect to the peak positions of the four bands around ~90 cm-1 , ~240 cm-1 , ~370 cm-1 , and 400 cm-1 , observed in the experimental spectrum, in spite that the physics involved in the two spectrum was not exactly the same: the experimental spectrum includes the contributions from the dipeptide and the water molecules interacting with the solute, while the theoretical one is just concerned with the solute molecule, influenced by solvation. Two major discrepancies between the theoretical and experimental spectra, one in the band intensity around ~100 cm-1 , and the other in the peak positions around ~370 cm-1 , are discussed in terms of the fluctuation mode of water molecules interacting with the dipeptide, which is not taken explicitly into account in the theoretical analysis.

Differential mode of cholesterol inclusion with 2-hydroxypropyl-cyclodextrins increases safety margin in treatment of Niemann-Pick disease type C.

BACKGROUND AND PURPOSE: Niemann-Pick disease type C (NPC) is a lysosomal storage disorder with disrupted intracellular cholesterol trafficking. A cyclic heptasaccharide, 2-hydroxypropyl-ß-cyclodextrin (HP-ß-CD), is a cholesterol solubilizer that is being developed to treat NPC, but its ototoxicity and pulmonary toxicity remain important issues. We have characterized 2-hydroxypropyl-γ-cyclodextrin (HP-γ-CD), a cyclic octasaccharide with a larger cavity than HP-ß-CD, as a candidate drug to treat NPC. However, the molecular target of HP-γ-CD with respect to NPC and its potential for clinical application are still unclear. EXPERIMENTAL APPROACH: We investigated the mode of interaction between HP-γ-CD and cholesterol by phase-solubility analysis, proton NMR spectroscopy and molecular dynamics simulations. We then evaluated the therapeutic effects of HP-γ-CD compared with HP-ß-CD using cellular and murine NPC models. Mouse auditory and pulmonary function tests were also conducted. KEY RESULTS: HP-γ-CD solely formed a 1:1 inclusion complex with cholesterol with an affinity similar to that of HP-ß-CD. In vitro, HP-γ-CD and HP-ß-CD amelioration of NPC-related manifestations was almost equivalent at lower concentrations. However, at higher concentrations, the cholesterol inclusion mode of HP-ß-CD shifted to the highly soluble 2:1 complex whereas that of HP-γ-CD maintained solely the 1:1 complex. The constant lower cholesterol solubilizing ability of HP-γ-CD conferred it with significantly reduced toxicity compared with HP-ß-CD, but equal efficacy in treating a mouse model of NPC. CONCLUSIONS AND IMPLICATIONS: HP-γ-CD can serve as a fine-tuned cholesterol solubilizer for the treatment of NPC with a wider safety margin than HP-ß-CD in terms of ototoxicity and pulmonary toxicity.

Computational Screening of a Functional Cyclodextrin Derivative for Suppressing a Side Effect of Doxorubicin.

The binding affinity of the beta-cyclodextrin (ß-CyD) derivatives with Doxorubicin (Dox) is evaluated by means of the 3D-RISM/KH theory combined with the molecular dynamics simulation in order to screen the compounds for suppressing a side-effect of the cancer drug. A protocol revised for the external and conformational entropies of the host-guest system is employed to calculate the binding free energy. It is found that the direct interactions of CyD with Dox and the desolvation free-energies of the both compounds largely cancel out to leave moderate contributions to the affinity, which are comparable to those from the entropies. The results shed light on the entropy terms for determining the binding affinity, although the external-entropy terms are essentially constant over all the compounds examined and do not affect the screening. The theoretical result is compared with the experimental data of the association constant for a CyD derivative which was predicted to be the best compound by the preliminary calculation without the entropy terms.

Molecular Recognition and Self-Organization in Life Phenomena Studied by a Statistical Mechanics of Molecular Liquids, the RISM/3D-RISM Theory.

There are two molecular processes that are essential for living bodies to maintain their life: the molecular recognition, and the self-organization or self-assembly. Binding of a substrate by an enzyme is an example of the molecular recognition, while the protein folding is a good example of the self-organization process. The two processes are further governed by the other two physicochemical processes: solvation and the structural fluctuation. In the present article, the studies concerning the two molecular processes carried out by Hirata and his coworkers, based on the statistical mechanics of molecular liquids or the RISM/3D-RISM theory, are reviewed.

New Protocol for Predicting the Ligand-Binding Site and Mode Based on the 3D-RISM/KH Theory.

An efficient algorithm to find the binding position and mode of small ligands bound at an active site of protein is proposed based on the spatial distribution function (SDF) obtained from the three-dimensional reference interaction site model (3D-RISM) theory with the Kovalenko-Hirata (KH) closure relation. The ligand examined includes hydrophobic, acidic, and basic molecules and zwitterions. Eighteen different types of proteins, which serve as targets for those ligands, are selected to examine the robustness of the algorithm. An imaginary atom, referred to as an "anchor site", is defined at the center of geometry of a ligand molecule that serves as a center for searching the binding position and mode of the ligand molecule in the translational and rotational spaces. The probable binding sites (PBSs) are identified based on the SDFs of the ligand molecules around the protein, and the PBS is ranked according to the peak height of SDF. The deviations from the mean height of the peak values of SDFs for 50 PBSs are analyzed based on the z-score, which is a measure of prominence of the site. The PBS found at the closest distance from the anchor site of the crystal structure is referred to as the "nearest site". The orientation of the ligand molecule at each PBS is explored by changing the Euler angles, and the most probable binding mode is determined based on the superposition approximation. The binding position of ligand molecules is successfully predicted as one of the distinct peaks in SDF of the anchor site, with a few exceptions. The binding mode of the ligand molecule predicted based on the superposition approximation is consistent with the X-ray crystal structure in nine systems, a half of the systems investigated. The significance of the results is discussed in detail. An application of the new protocol to fragment-based drug discovery is suggested.

Role of Mg2+ Ions in DNA Hydrolysis by EcoRV, Studied by the 3D-Reference Interaction Site Model and Molecular Dynamics.

The role of Mg2+ ions during precursor formation in DNA hydrolysis by the homodimeric restriction enzyme EcoRV was elucidated based on the 3D-reference interaction site model (RISM) theory and the molecular dynamics (MD) simulation. From an analysis of the spatial distribution of Mg2+ in an active site using 3D-RISM, we identified a new position for Mg2+ in the X-ray EcoRV-DNA complex structure ( 1rvb ). We refer to the position as site IV†. Site IV† is almost at the same position as that of a Ca2+ ion in the superimposed X-ray crystallographic active-site structure of the PvuII-DNA complex ( 1f0o ). 3D-RISM was also used to locate the position of water molecules, including the water nucleophile at the active site. MD simulations were carried out with the initial structure having two Mg2+ ions at site IV† and at site I*, experimentally identified by Horton et al., to find a stable complex structure in which the DNA fragment was rearranged to orient the scissile bond direction toward the water nucleophile. The equilibrium active-site structure of the EcoRV-DNA complex obtained from the MD simulation was similar to the superimposed X-ray crystallographic structure of the BamHI-DNA complex ( 2bam ). In the active-site structure, two metal ions have almost the same position (≤1.0 Å) as that of 2bam , and the scissile phosphate is twisted to orient the scissile bond toward the water nucleophile, as is the case in 2bam . We propose the equilibrium active-site structure obtained in this study as a precursor for the hydrolysis reaction of EcoRV.

Predicting the Binding Mode of 2-Hydroxypropyl-ß-cyclodextrin to Cholesterol by Means of the MD Simulation and the 3D-RISM-KH Theory.

It has been found that a cyclodextrin derivative, 2-hydroxypropyl-ß-cyclodextrin (HPßCD), has reasonable therapeutic effect on Niemann-Pick disease type C, which is caused by abnormal accumulation of unesterified cholesterol and glycolipids in the lysosomes and shortage of esterified cholesterol in other cellular compartments. We study the binding affinity and mode of HPßCD with cholesterol to elucidate the possible mechanism of HPßCD for removing cholesterol from the lysosomes. The dominant binding mode of HPßCD with cholesterol is found based on the molecular dynamics simulation and a statistical mechanics theory of liquids, or the three-dimensional reference interaction site model theory with Kovalenko-Hirata closure relation. We examine the two types of complexes between HPßCD and cholesterol, namely, one-to-one (1:1) and two-to-one (2:1). It is predicted that the 1:1 complex makes two or three types of stable binding mode in solution, in which the ßCD ring tends to be located at the edge of the steroid skeleton. For the 2:1 complex, there are four different types of the complex conceivable, depending on the orientation between the two HPßCDs: head-to-head (HH), head-to-tail (HT), tail-to-head (TH), and tail-to-tail (TT). The HT and HH cyclodextrin dimers show higher affinity to cholesterol compared to the other dimers and to all the binding modes of 1:1 complexes. The physical reason why the HT and HH dimers have higher affinity compared to the other complexes is discussed based on the consistency with the 1:1 complex. On the one hand, in case of the HT and HH dimers, the position of each CD in the dimer along the cholesterol chain comes right on or close to one of the positions where a single CD makes a stable complex. On the other hand, one of the CD molecules is located on unstable region along the cholesterol chain, for the case of TH and TT dimers.

Perspective: Structural fluctuation of protein and Anfinsen's thermodynamic hypothesis.

The thermodynamics hypothesis, casually referred to as "Anfinsen's dogma," is described theoretically in terms of a concept of the structural fluctuation of protein or the first moment (average structure) and the second moment (variance and covariance) of the structural distribution. The new theoretical concept views the unfolding and refolding processes of protein as a shift of the structural distribution induced by a thermodynamic perturbation, with the variance-covariance matrix varying. Based on the theoretical concept, a method to characterize the mechanism of folding (or unfolding) is proposed. The transition state, if any, between two stable states is interpreted as a gap in the distribution, which is created due to an extensive reorganization of hydrogen bonds among back-bone atoms of protein and with water molecules in the course of conformational change. Further perspective to applying the theory to the computer-aided drug design, and to the material science, is briefly discussed.

A Systematic Analysis of the Binding Affinity between the Pim-1 Kinase and Its Inhibitors Based on the MM/3D-RISM/KH Method.

A systematic study of the binding affinities of 16 lead compounds targeting the Pim-1 kinase based on the 3D-RISM/KH theory and MD simulations is reported. The results show a correlation coefficient R = 0.69 between the theoretical and experimental values of the binding free energy. This demonstrates that the method is applicable to the problem of compound screening and lead optimization, for which relative values of the free energy among the compounds have significance. We elucidate the contribution of the ligand fragments to the binding free energy. Our results indicate that the interactions between the residues and the triazolo[4,3-b]pyridazine scaffold as well as the phenyl ring of the ligand molecule make significant contributions to stabilization of the complex. Using the 3D-RISM/KH theory, we further analyze the probability distribution of a ligand fragment around the protein-ligand complex in which the substituent around the phenyl ring is removed from the ligand. The results demonstrate that the 3D-RISM/KH theory is capable of predicting the position of substitution on a ligand that has a higher affinity to a target protein.

A molecular theory of the structural dynamics of protein induced by a perturbation.

An equation to describe the structural dynamics of protein molecule induced by a perturbation such as a photo-excitation is derived based on the linear response theory, which reads 𝐑α(t)=𝐑α(t=∞)-1kBT∑γ⟨Δ𝐑α(t)Δ𝐑γ⟩eq(0)⋅𝐟γ(0). In the equation, α and γ distinguish atoms in protein, 𝐟γ(0) denotes a perturbation at time t = 0, 𝐑α(t) the average position (or structure) of protein atom α at time t after the perturbation being applied, and 𝐑a(t=∞) the position at t=∞. ⟨Δ𝐑α(t)Δ𝐑γ⟩eq(0) is a response function in which Δ𝐑α(t) is the fluctuation of atom α at time t in the equilibrium system. The perturbation is defined in terms of the free energy difference between perturbed and unperturbed equilibrium-states, which includes interactions between solute and solvent as well as those among solvent molecules in a renormalized manner. The response function signifies the time evolution of the variance-covariance matrix of the structural fluctuation for the unperturbed system. A theory to evaluate the response function ⟨Δ𝐑α(t)Δ𝐑γ⟩eq(0) is also proposed based on the Kim-Hirata theory for the structural fluctuation of protein [B. Kim and F. Hirata, J. Chem. Phys. 138, 054108 (2013)]. The problem reduces to a simple eigenvalue problem for a matrix which includes the friction and the second derivative of the free energy surface of protein with respect to its atomic coordinates.

Predicting the binding free energy of the inclusion process of 2-hydroxypropyl-ß-cyclodextrin and small molecules by means of the MM/3D-RISM method.

A protocol to calculate the binding free energy of a host-guest system is proposed based on the MM/3D-RISM method, taking cyclodextrin derivatives and their ligands as model systems. The protocol involves the procedure to identify the most probable binding mode (MPBM) of receptors and ligands by means of the umbrella sampling method. The binding free energies calculated by the MM/3D-RISM method for the complexes of the seven ligands with the MPBM of the cyclodextrin, and with the fluctuated structures around it, are in agreement with the corresponding experimental data in a semi-quantitative manner. It suggests that the protocol proposed here is promising for predicting the binding affinity of a small ligand to a relatively rigid receptor such as cyclodextrin.

A 3D-RISM/RISM study of the oseltamivir binding efficiency with the wild-type and resistance-associated mutant forms of the viral influenza B neuraminidase.

The binding affinity of oseltamivir to the influenza B neuraminidase and to its variants with three single substitutions, E119G, R152K, and D198N, is investigated by the MM/3D-RISM method. The binding affinity or the binding free energy of ligand to receptor was found to be determined by a subtle balance of two major contributions that largely cancel out each other: the ligand-receptor interactions and the dehydration free energy. The theoretical results of the binding affinity of the drug to the mutants reproduced the observed trend in the resistivity, measured by IC50 ; the high-level resistance of E119G and R152K, and the low-level resistance of D198N. For E119G and R152K, reduction of the direct drug-target interaction, especially at the mutated residue, is the main source of high-level oseltamivir resistance. This phenomenon, however, is not found in the D198N strain, which is located in the framework of the active-site.

Water Turns the "Non-biological" Fluctuation of Protein into "Biological" One.

Structural fluctuation of protein is not just an mechanical "oscillation," but an event induced by an interplay of mechanical and thermodynamic processes in which water plays crucial role. The chapter is devoted to provide a theoretical description concerning the concept of structural fluctuation of protein, based on methods of the statistical mechanics.

Structural fluctuation of proteins induced by thermodynamic perturbation.

A theory to describe structural fluctuations of protein induced by thermodynamic perturbations, pressure, temperature, and denaturant, is proposed. The theory is formulated based on the three methods in the statistical mechanics: the generalized Langevin theory, the linear response theory, and the three dimensional interaction site model (3D-RISM) theory. The theory clarifies how the change in thermodynamic conditions, or a macroscopic perturbation, induces the conformational fluctuation, which is a microscopic property. The theoretical results are applied, on the conceptual basis, to explain the experimental finding by Akasaka et al., concerning the NMR experiment which states that the conformational change induced by pressure corresponds to structural fluctuations occurring in the ambient condition. A method to evaluate the structural fluctuation induced by pressure is also suggested by means of the 3D-RISM and the site-site Kirkwood-Buff theories.

Cavity as a source of conformational fluctuation and high-energy state: high-pressure NMR study of a cavity-enlarged mutant of T4 lysozyme.

Although the structure, function, conformational dynamics, and controlled thermodynamics of proteins are manifested by their corresponding amino acid sequences, the natural rules for molecular design and their corresponding interplay remain obscure. In this study, we focused on the role of internal cavities of proteins in conformational dynamics. We investigated the pressure-induced responses from the cavity-enlarged L99A mutant of T4 lysozyme, using high-pressure NMR spectroscopy. The signal intensities of the methyl groups in the (1)H/(13)C heteronuclear single quantum correlation spectra, particularly those around the enlarged cavity, decreased with the increasing pressure, and disappeared at 200 MPa, without the appearance of new resonances, thus indicating the presence of heterogeneous conformations around the cavity within the ground state ensemble. Above 200 MPa, the signal intensities of >20 methyl groups gradually decreased with the increasing pressure, without the appearance of new resonances. Interestingly, these residues closely matched those sensing a large conformational change between the ground- and high-energy states, at atmospheric pressure. (13)C and (1)H NMR line-shape simulations showed that the pressure-induced loss in the peak intensity could be explained by the increase in the high-energy state population. In this high-energy state, the aromatic side chain of F114 gets flipped into the enlarged cavity. The accommodation of the phenylalanine ring into the efficiently packed cavity may decrease the partial molar volume of the high-energy state, relative to the ground state. We suggest that the enlarged cavity is involved in the conformational transition to high-energy states and in the volume fluctuation of the ground state.

Determination of the formylglycinamide ribonucleotide amidotransferase ammonia pathway by combining 3D-RISM theory with experiment.

Molecular tunnels in enzyme systems possess variable architecture and are therefore difficult to predict. In this work, we design and apply an algorithm to resolve the pathway followed by ammonia using the bifunctional enzyme formylglycinamide ribonucleotide amidotransferase (FGAR-AT) as a model system. Though its crystal structure has been determined, an ammonia pathway connecting the glutaminase domain to the 30 Å distal FGAR/ATP binding site remains elusive. Crystallography suggested two purported paths: an N-terminal-adjacent path (path 1) and an auxiliary ADP-adjacent path (path 2). The algorithm presented here, RismPath, which enables fast and accurate determination of solvent distribution inside a protein channel, predicted path 2 as the preferred mode of ammonia transfer. Supporting experimental studies validate the identity of the path, and results lead to the conclusion that the residues in the middle of the channel do not partake in catalytic coupling and serve only as channel walls facilitating ammonia transfer.

Massively parallel implementation of 3D-RISM calculation with volumetric 3D-FFT.

A new three-dimensional reference interaction site model (3D-RISM) program for massively parallel machines combined with the volumetric 3D fast Fourier transform (3D-FFT) was developed, and tested on the RIKEN K supercomputer. The ordinary parallel 3D-RISM program has a limitation on the number of parallelizations because of the limitations of the slab-type 3D-FFT. The volumetric 3D-FFT relieves this limitation drastically. We tested the 3D-RISM calculation on the large and fine calculation cell (2048(3) grid points) on 16,384 nodes, each having eight CPU cores. The new 3D-RISM program achieved excellent scalability to the parallelization, running on the RIKEN K supercomputer. As a benchmark application, we employed the program, combined with molecular dynamics simulation, to analyze the oligomerization process of chymotrypsin Inhibitor 2 mutant. The results demonstrate that the massive parallel 3D-RISM program is effective to analyze the hydration properties of the large biomolecular systems.

Theoretical characterization of the "ridge" in the supercritical region in the fluid phase diagram of water.

The density fluctuation of water in the supercritical region was investigated theoretically using the reference interaction site model theory combined with the Kovalenko-Hirata closure relation, the so-called RISM-KH theory. The density fluctuation was evaluated by the numerical differentiation of density with respect to pressure at constant temperature. The density fluctuations plotted against density show finite maxima along a line slightly off from the critical isochore, in accordance with experimental results. The microscopic structures of water on both regions that were separated by the line were investigated by analyzing the site-site radial distribution functions. The analysis clearly indicates that the structure is determined by the two effects featuring liquid states: the packing or volume exclusion effect and the screening of the Coulomb interaction or the hydrogen bond, both becoming more important at higher densities. An interplay of the two effects creates maxima of the density fluctuation in the supercritical region of water.