New Methodology for Computer-Aided Modelling of Biomolecular Structure and Dynamics 1. Non-Cyclic Structures

Abstract

A general methodology is proposed for the conformational modelling of biomolecular systems. The approach allows one: (i) to describe the system under investigation by an arbitrary set of internal variables, i.e., torsion angles, bond angles, and bond lengths; it offers a possibility to pass from the free structure to a completely fixed one with the number of variables from 3N to zero, respectively, where N is the number of atoms; (ii) to consider both, a single molecule and a complex of many molecules, (e.g., proteins, water, ligands, etc.) in terms of one universal model; (iii) to study the dynamics of the system using explicit analytical Lagrangian equations of motion, thus opening up possibilities for investigations of slow concerted motions such as domain oscillations in proteins etc.; (iv) to calculate the partial derivatives of various functions of conformation, e.g., the conformatinal energy or external constraints imposed, using a standard efficient procedure regardless of the variables and the structure of the system. The approach is meant to be used in various investigations concerning the conformations and dynamics of biomacromolecules.     

Normal mode analysis as a method to derive protein dynamics information from the Protein Data Bank

Normal mode analysis (NMA) can facilitate quick and systematic investigation of protein dynamics using data from the Protein Data Bank (PDB). We developed an elastic network model-based NMA program using dihedral angles as independent variables. Compared to the NMA programs that use Cartesian coordinates as independent variables, key attributes of the proposed program are as follows: (1) chain connectivity related to the folding pattern of a polypeptide chain is naturally embedded in the model; (2) the full-atom system is acceptable, and owing to a considerably smaller number of independent variables, the PDB data can be used without further manipulation; (3) the number of variables can be easily reduced by some of the rotatable dihedral angles; (4) the PDB data for any molecule besides proteins can be considered without coarse-graining; and (5) individual motions of constituent subunits and ligand molecules can be easily decomposed into external and internal motions to examine their mutual and intrinsic motions. Its performance is illustrated with an example of a DNA-binding allosteric protein, a catabolite activator protein. In particular, the focus is on the conformational change upon cAMP and DNA binding, and on the communication between their binding sites remotely located from each other. In this illustration, NMA creates a vivid picture of the protein dynamics at various levels of the structures, i.e., atoms, residues, secondary structures, domains, subunits, and the complete system, including DNA and cAMP. Comparative studies of the specific protein in different states, e.g., apo- and holo-conformations, and free and complexed configurations, provide useful information for studying structurally and functionally important aspects of the protein.     

Mechanical force can fine-tune redox potentials of disulfide bonds

Mechanical force applied along a disulfide bond alters its rate of reduction. We here aimed at quantifying the direct effect of force onto the chemical reactivity of a sulfur-sulfur bond in contrast to indirect, e.g., steric or mechanistic, influences. To this end, we evaluated the dependency of a disulfide bond's redox potential on a pulling force applied along the system. Our QM/MM simulations of cystine as a model system take conformational dynamics and explicit solvation into account and show that redox potentials increase over the whole range of forces probed here (30-3320 pN), and thus even in the absence of a significant disulfide bond elongation (<500 pN). Instead, at low forces, dihedrals and angles, as the softer degrees of freedom are stretched, contribute to the destabilization of the oxidized state. We find physiological forces to be likely to tune the disulfide's redox potentials to an extent similar to the tuning within proteins by point mutations.     

Improving consensus structure by eliminating averaging artifacts

Background  

Common structural biology methods (i.e., NMR and molecular dynamics) often produce ensembles of molecular structures. Consequently, averaging of 3D coordinates of molecular structures (proteins and RNA) is a frequent approach to obtain a consensus structure that is representative of the ensemble. However, when the structures are averaged, artifacts can result in unrealistic local geometries, including unphysical bond lengths and angles.     

Interdependence of conformational variables in double-helical DNA.

DNA exhibits conformational polymorphism, with the details depending on the sequence and its environment. To understand the mechanisms of conformational polymorphism and these transitions, we examine the interrelationships among the various conformational variables of DNA. In particular, we examine the stress-strain relation among conformational variables, describing base-pair morphology and their effects on the backbone conformation. For the calculation of base pairs, we use the method previously developed to calculate averages over conformational variables of DNA. Here we apply this method to calculate the Boltzmann averages of conformational variables for fixed values of one particular conformational variable, which reflects the strain in the structure responding to a particular driving stress. This averaging over all but one driving variable smooths the usual rough energy surface to permit observation of the effects of one conformational variable at a time. The stress-strain analyses of conformational variables of base pair slide, twist, and roll, which exhibit characteristic changes during the conformational transition of DNA, have shown that the conformational changes of base pairs are strongly correlated with one another. Furthermore, the stress-strain relations are not symmetrical with respect to these variables, i.e., the response of one coordinate to another is different from the reverse direction. We also examine the effect of conformational changes in base-pair variables on the sugar-backbone conformation by using the minimization method we developed. The conformational changes of base pairs affect the sugar pucker and other dihedral angles of the backbone of DNA, but each variable affects the sugar-backbone differently. In particular, twist is found to have the most influence in affecting the sugar pucker and backbone conformation. These calculated conformational changes in base pairs and backbone segments are consistent with experimental observations and serve to validate the calculation method.     

Conformational dynamics of helix S6 from Shaker potassium channel: simulation studies

Prolines in transmembrane (TM) alpha-helices are believed to play an important structural and/or functional role in membrane proteins. At a structural level a proline residue distorts alpha-helical structure due to the loss of at least one stabilizing backbone hydrogen bond, and introduces flexibility in the helix that may result in substantial kink and swivel motions about the effective "hinge." At a functional level, for example in Kv channels, it is believed that proline-induced molecular hinges may have a direct role in gating, i.e., the conformational change linked to opening/closing the channel to movement of ions. In this article we study the conformational dynamics of the S6 TM helix from of the Kv channel Shaker, which possesses the motif PVP--a motif that is conserved in Kv channels. We perform multiple molecular dynamics simulations of single S6 helices in a membrane-mimetic environment in order to effectively map the kink-swivel conformational space of the protein, exploiting the ability of multiple simulations to achieve greater sampling. We show that the presence of proline locally perturbs the helix, disrupting local dihedral angles and producing local twist and unwinding in the region of the hinge--an effect that is relaxed with distance from the PVP motif. We furthermore show that motions about the hinge are highly anisotropic, reflecting a preferred region of kink-swivel conformation space that may have implications for the gating process.     

Structural and functional investigations of cholinesterases by means of affinity electrophoresis

1. After a brief survey of the basic affinity electrophoresis concepts, the usual ways for preparing affinity electrophoresis ligands are examined. 2. Then results obtained on cholinesterases are reviewed. This section includes (a) structural and functional investigations on anionic sites, i.e., study of ligand-induced conformational change, organophosphate-induced "aging," genetic variants, and active-site topology; and (b) characterization of cholinesterase conjugates (hybrid proteins) and glycoinositol phospholipid-anchored cholinesterases. 3. The future prospects of affinity electrophoresis, e.g., investigations on the esteratic site and exploration of the carbohydrate moiety, are emphasized in the concluding section.     

Algorithmic generation of freely jointed hard sphere chains and properties of their inertial tensors

A statistical algorithm, capable of generating a large number of freely jointed hard sphere chains, is presented. This is the first of a series of algorithms being developed to model unfolded proteins by different modes of hard sphere chains. The aim of these studies is to systematically investigate the effects of different factors, such as atomic radii, bond angles, torsion angles, chain length, etc., on the conformation of unfolded proteins and other random polymers. As continuous models, various types of hard sphere chains enable one to isolate the aforementioned factors one at a time for investigation and thus are advantageous over discrete lattice models. In particular, the freely jointed hard sphere chain model allows one to evaluate the excluded volume effect. As a first step in this endeavor, the average determinant D(N, r) and the average trace T(N, r) of the inertial tensor A of the random chains were calculated at various sphere radii r and chain lengths N. It is found that both the average determinant D(N, r) and the average trace T(N, r) scale linearly with chain length N after logarithmic transformation. However, the critical exponent of D(N, r) increases with r faster than that of T(N, r) as a result of the non-commutativity between the det operator and the average operator < >. The significance of the algorithm and the results obtained on understanding random polypeptide chains are discussed.     

Pex, analytical tools for PDB files. II. H-Pex: noncanonical H-bonds in alpha-helices

We use the H-Pex (Thomas et al., this issue) to analyze the main chain interactions in 131 proteins. In antiparallel beta-sheets, the geometry of the N...O bond is: median N...O distances, 2.9 SA, C==O...N angles at 154 degrees and the C alpha--C==O...H angles are dispersed around 3 degrees. In some instances, the other side of the C==O axis is occupied by a HC alpha. As recently supported by Vargas et al. (J Am Chem Soc 2000;122:4750-4755) C alpha H...O and NH...O could cooperate to sheet stability. In alpha-helices, the main chain C==O interact with the NH of their n + 4 neighbor on one side, and with a C beta H or C gamma H on the other side. The median O...N distance (3.0 A) and C==N angle (147 degrees) suggest a canonical H-bond, but the C alpha--C==O...H dihedral angle invalidates this option, since the hydrogen attacks the oxygen at 122 degrees, i.e., between the sp(2) and pi orbitals. This supports that the H-bond is noncanonical. In many instances, the C gamma H or the C beta H of the n + 4 residue stands opposite to the NH with respect to the oxygen. Therefore, we propose that, in alpha-helices, the C gamma H or C beta H and the NH of the n + 4 residue hold the oxygen like an electrostatic pincher. Proteins 2001;43:37-44.     

Molecular dynamics simulation of Lewis blood groups and related oligosaccharides.

Molecular dynamics simulations without explicit inclusion of solvent molecules have been performed to study the motions of Lewisa and Lewisb blood group oligosaccharides, and two blood group A tetrasaccharides having type I and type II core chains. The blood group H trisaccharide has also been studied and compared with the blood group A type II core chain. The potential energy surface developed by Rasmussen and co-workers was used with the molecular mechanics code CHARMM. The lowest energy minima of the component disaccharide fragments were obtained from conformational energy mapping. The lowest energy minima of these disaccharide fragments were used to build the tri- and tetrasaccharides that were further minimized before the actual heating/equilibration and dynamics simulations. The trajectories of the disaccharide fragments, e.g., Fuc alpha- (1----4)GlcNAc, Gal beta-(1----4)GlcNAc, etc., show transitions among various minima. However, the oligosaccharides were found to be dynamically stable and no transitions to other minimum energy conformations were observed in the time series of the glycosidic dihedral angles even during trajectories as long as 300 ps. The stable conformations of the glycosidic linkages in the oligosaccharides are not necessarily the same as the minimum energy conformation of the corresponding isolated disaccharides. The average fluctuations of the glycosidic angles in the oligosaccharides were well within the range of +/- 15 degrees. The results of these trajectory calculations were consistent with the relatively rigid single-conformation models derived for these oligosaccharides from 1H-nmr data.     

An analysis of the class of gene regulatory functions implied by a biochemical model

Understanding the integrated behavior of genetic regulatory networks, in which genes regulate one another's activities via RNA and protein products, is emerging as a dominant problem in systems biology. One widely studied class of models of such networks includes genes whose expression values assume Boolean values (i.e., on or off). Design decisions in the development of Boolean network models of gene regulatory systems include the topology of the network (including the distribution of input- and output-connectivity) and the class of Boolean functions used by each gene (e.g., canalizing functions, post functions, etc.). For example, evidence from simulations suggests that biologically realistic dynamics can be produced by scale-free network topologies with canalizing Boolean functions. This work seeks further insights into the design of Boolean network models through the construction and analysis of a class of models that include more concrete biochemical mechanisms than the usual abstract model, including genes and gene products, dimerization, cis-binding sites, promoters and repressors. In this model, it is assumed that the system consists of N genes, with each gene producing one protein product. Proteins may form complexes such as dimers, trimers, etc. The model also includes cis-binding sites to which proteins may bind to form activators or repressors. Binding affinities are based on structural complementarity between proteins and binding sites, with molecular binding sites modeled by bit-strings. Biochemically plausible gene expression rules are used to derive a Boolean regulatory function for each gene in the system. The result is a network model in which both topological features and Boolean functions arise as emergent properties of the interactions of components at the biochemical level. A highly biased set of Boolean functions is observed in simulations of networks of various sizes, suggesting a new characterization of the subset of Boolean functions that are likely to appear in gene regulatory networks.     

The conformational analysis of adenosine triphosphate by classical potential energy calculations

The conformational analysis of adenosine triphosphate was conducted by using classical potential energy calculations. All rotatable bonds were examined, i.e., no dihedral angles were fixed at predetermined conformations except for the ribofuranose ring, which was held in the C(3′)-endo conformation—the conformation observed for adenosine in the crystal state. The energy terms included in the total energy expression consist of nonbonded pairwise interaction, electrostatic pairwise interaction, free energy of solvation, and torsional bond potentials. Two separate approaches were used in the conformational analyses. The first consisted of a sequential fragment approach were four bonds were rotated simultaneously at 30° increments. Each fragment overlapped the preceding one by at least one bond. All rotors were then simultaneously examined at their minima and at ±15°. The second approach consisted of a coarse grid search where all rotors were examined simultaneously, but only at staggered positions. The low-energy conformations thus obtained were then used as starting conformations for a minimization routine based on the method of conjugate directions. The first approach required about 40 hr of central processing unit (CPU) computer time, while the coarse grid/minimization approach required about 4 hr of CPU time. Both the sequential fragment approach and the minimization approach yielded lowest-energy conformations which are remarkably similar to the solid-state conformation of C(3′)-endo ATP.     

Conformational analysis of lignin models: a chemometric approach

In the present work, conformational analysis of lignin models was accomplished by considering four cross-link types (3–5′, β-5′, α-O-4 and β-O-4) and three monomer units [guaiacyl (G), p-hydroxyphenyl (H) and syringyl (S)]. Analysis involving the 3–5′ and β-5′ dimers was conducted following the standard procedure, i.e., rotating the monomers around the single bond. On the other hand, analysis of α-O-4 and β-O-4 dimers followed a distinct protocol with the aid of an interesting chemometric tool called Box-Behnken (BB) design. This methodology was applied with the aim of screening the most relevant dihedral angles. The results show that the conformational space for large systems with several dihedral angles can be mapped satisfactorily through the BB approach, reducing the number of dimensions to be treated at the quantum mechanical level. Furthermore, the quantum mechanics-chemometry-quantum mechanics (QM/BB/QM) method proposed here allows us to determine calculated torsional angles for lignin models in good agreement with crystallographic data for some model compounds.     

Dynamically Driven Protein Allostery Exhibits Disparate Responses for Fast and Slow Motions

There is considerable interest in the dynamic aspect of allosteric action, and in a growing list of proteins allostery has been characterized as being mediated predominantly by a change in dynamics, not a transition in conformation. For considering conformational dynamics, a protein molecule can be simplified into a number of relatively rigid microdomains connected by joints, corresponding to, e.g., communities and edges from a community network analysis. Binding of an allosteric activator strengthens intermicrodomain coupling, thereby quenching fast (e.g., picosecond to nanosecond) local motions but initiating slow (e.g., microsecond to millisecond), cross-microdomain correlated motions that are potentially of functional importance. This scenario explains allosteric effects observed in many unrelated proteins.     

Combined use of molecular dynamics simulations and NMR to explore peptide bond isomerization and multiple intramolecular hydrogen-bonding possibilities in a cyclic pentapeptide, cyclo(Gly-Pro-D-Phe-Gly-Val).

The conformational behavior of a model cyclic pentapeptide--cyclo(Gly-L-Pro-D-Phe-Gly-L-Val)--has been explored through the combined use of in vacuo molecular dynamics simulations and a range of nmr experiments (preceding paper). The molecular dynamics analysis suggests that, despite the conformational constraints imposed by formation of the pentapeptide cycle, this pentapeptide undergoes conformational transitions between various hydrogen-bonded conformations, characterized by low energy barriers. An inverse gamma turn with Pro in position i + 1 and a gamma turn with D-Phe in position i + 1 are two alternatives occurring frequently. Like other DLDDL cyclic pentapeptides, cyclo(Gly-Pro-D-Phe-Gly-Val) is also stabilized by an inverse gamma-turn structure with the beta-branched Val residue in position i + 1, and this hydrogen bond is retained in the different conformational families. The gamma-turn around D-Phe3 and the inverse gamma turn around Val5 are consistent with the nmr observations. 3JNH-CH alpha coupling constants of the all-trans forms were calculated from one of the molecular dynamics trajectories and are comparable to nmr experimental data, suggesting that the conformational states visited during the simulation are representative of the conformational distribution in solution. In addition to the equilibrium among various hydrogen-bonded all-trans conformers, the observation in nmr spectra of two sets of resonances for all peptide protons indicated a slow conformational interconversion of the Gly-Pro peptide bond between trans and cis isomers. The activation energy between these two conformers was determined experimentally by magnetization transfer and was calculated by high temperature constrained molecular dynamics simulation. Both methods yield a free energy of activation of ca. 20 kcal/mol. Furthermore, the free energy of activation is dependent on the direction of rotation of the Gly-Pro peptide bond.     

Nonorthogonal tight-binding model with H–C–N–O parameterisation

A parametric nonorthogonal tight-binding model (NTBM1) with the set of parameters for H–C–N–O systems is presented. This model compares well with widely used semi-empirical AM1 and PM3/PM7 models but contains less fitting parameters per atom. All NTBM1 parameters are derived based on a criterion of the best agreement between the calculated and experimental values of bond lengths, valence angles and binding energies for various H–C–N–O molecules. Results for more than 200 chemical compounds are reported. Parameters are currently available for hydrogen, carbon, nitrogen, oxygen atoms and corresponding interatomic interactions. The model has a good transferability and can be used for both relaxation of large molecular systems (e.g., high-molecular compounds or covalent cluster complexes) and long-timescale molecular dynamics simulation (e.g., modelling of thermal decomposition processes). The program package based on this model is available for download at no cost from http://ntbm.info.     

Dynamic properties of the first enzymatic reaction steps of porcine pancreatic elastase. How rigid is the active site of the native enzyme? Molecular dynamics simulation

Two molecular dynamics simulations (100 and 50 ps) of native porcine pancreatic elastase i.e., without bound substrate and with the active site hydrated by a dome of water (630 molecules) have been performed. Dynamical properties of the catalytic tetrad have been examined. While relative conformations of the Asp 102, His 57, and Ser 214 are rather stable in time, the side chain of Ser 195 undergoes several conformational changes. No preferences are observed for the formation of a hydrogen bond between the O gamma-H group (Ser 195) and nitrogen N, (His 57). A cluster of ordered water molecules effectively competes with the H-O gamma group (Ser 195) and thereby prevents the formation of this H bond, which is generally agreed to be crucial for catalysis.     

Learning generative models of molecular dynamics

We introduce three algorithms for learning generative models of molecular structures from molecular dynamics simulations. The first algorithm learns a Bayesian-optimal undirected probabilistic model over user-specified covariates (e.g., fluctuations, distances, angles, etc). L1 regularization is used to ensure sparse models and thus reduce the risk of over-fitting the data. The topology of the resulting model reveals important couplings between different parts of the protein, thus aiding in the analysis of molecular motions. The generative nature of the model makes it well-suited to making predictions about the global effects of local structural changes (e.g., the binding of an allosteric regulator). Additionally, the model can be used to sample new conformations. The second algorithm learns a time-varying graphical model where the topology and parameters change smoothly along the trajectory, revealing the conformational sub-states. The last algorithm learns a Markov Chain over undirected graphical models which can be used to study and simulate kinetics. We demonstrate our algorithms on multiple molecular dynamics trajectories.     

Serine and threonine residues bend alpha-helices in the chi(1) = g(-) conformation

The relationship between the Ser, Thr, and Cys side-chain conformation (chi(1) = g(-), t, g(+)) and the main-chain conformation (phi and psi angles) has been studied in a selection of protein structures that contain alpha-helices. The statistical results show that the g(-) conformation of both Ser and Thr residues decreases their phi angles and increases their psi angles relative to Ala, used as a control. The additional hydrogen bond formed between the O(gamma) atom of Ser and Thr and the i-3 or i-4 peptide carbonyl oxygen induces or stabilizes a bending angle in the helix 3-4 degrees larger than for Ala. This is of particular significance for membrane proteins. Incorporation of this small bending angle in the transmembrane alpha-helix at one side of the cell membrane results in a significant displacement of the residues located at the other side of the membrane. We hypothesize that local alterations of the rotamer configurations of these Ser and Thr residues may result in significant conformational changes across transmembrane helices, and thus participate in the molecular mechanisms underlying transmembrane signaling. This finding has provided the structural basis to understand the experimentally observed influence of Ser residues on the conformational equilibrium between inactive and active states of the receptor, in the neurotransmitter subfamily of G protein-coupled receptors.     

ATP-induced conformational transition of denatured proteins.

Although the conformational change occurring in proteins upon ATP binding is important in many biological reactions, the mechanism by which ATP binding induces the conformational change is unknown. We found that ATP induces acid-unfolded (pH 2) ferricytochrome c or apomyoglobin to adopt a compact structure with a significant amount of alpha-helix and increased hydrophobicity. A very similar conformational transition was observed at neutral pH for an amphiphilic model polypeptide. The effectiveness of various adenine nucleotides in inducing the conformational transition was found to be proportional to their phosphate group contents, i.e., adenosine tetraphosphate greater than ATP greater than ADP greater than AMP. These results should be important when considering the mechanism of the ATP-induced conformational change in proteins during various biological reactions.