Search

A major challenge in drug design is to obtain compounds that bind selectively to their target receptors and do not cause side-effects by binding to other similar receptors. Here, we investigate strategies for applying COMBINE (COMparative BINding Energy) analysis, in conjunction with PIPSA (Protein Interaction Property Similarity Analysis) and ligand docking methods, to address this problem. We evaluate these approaches by application to diverse sets of inhibitors of three structurally related serine proteases of medical relevance: thrombin, trypsin, and urokinase-type plasminogen activator (uPA). We generated target-specific scoring functions (COMBINE models) for the three targets using training sets of ligands with known inhibition constants and structures of their receptor-ligand complexes. These COMBINE models were compared with the PIPSA results and experimental data on receptor selectivity. These scoring functions highlight the ligand-receptor interactions that are particularly important for binding specificity for the different targets. To predict target selectivity in virtual screening, compounds were docked into the three protein binding sites using the program GOLD and the docking solutions were re-ranked with the target-specific scoring functions and computed electrostatic binding free energies. Limits in the accuracy of some of the docking solutions and difficulties in scoring them adversely affected the predictive ability of the target specific scoring functions. Nevertheless, the target-specific scoring functions enabled the selectivity of ligands to thrombin versus trypsin and uPA to be predicted.

The factors that determine the extent to which diffusion and thermal activation processes govern electron transfer (ET) between proteins are debated. The process of ET between plastocyanin (PC) and cytochrome f (CytF) from the cyanobacterium Phormidium laminosum was initially thought to be diffusion-controlled but later was found to be under activation control (Schlarb-Ridley, B. G.; et al. Biochemistry 2005, 44, 6232). Here we describe Brownian dynamics simulations of the diffusional association of PC and CytF, from which ET rates were computed using a detailed model of ET events that was applied to all of the generated protein configurations. The proteins were modeled as rigid bodies represented in atomic detail. In addition to electrostatic forces, which were modeled as in our previous simulations of protein-protein association, the proteins interacted by a nonpolar desolvation (hydrophobic) force whose derivation is described here. The simulations yielded close to realistic residence times of transient protein-protein encounter complexes of up to tens of microseconds. The activation barrier for individual ET events derived from the simulations was positive. Whereas the electrostatic interactions between P. laminosum PC and CytF are weak, simulations for a second cyanobacterial PC-CytF pair, that from Nostoc sp. PCC 7119, revealed ET rates influenced by stronger electrostatic interactions. In both cases, the simulations imply significant contributions to ET from both diffusion and thermal activation processes.

Beta-Lactamases are responsible for bacterial resistance to beta-lactams and are thus of major clinical importance. However, the identity of the general base involved in their mechanism of action is still unclear. Two candidate residues, Glu166 and Lys73, have been proposed to fulfill this role. Previous studies support the proposal that Glu166 acts during the deacylation, but there is no consensus on the possible role of this residue in the acylation step. Recent experimental data and theoretical considerations indicate that Lys73 is protonated in the free beta-lactamases, showing that this residue is unlikely to act as a proton abstractor. On the other hand, it has been proposed that the pKa of Lys73 would be dramatically reduced upon substrate binding and would thus be able to act as a base. To check this hypothesis, we performed continuum electrostatic calculations for five wild-type and three beta-lactamase mutants to estimate the pKa of Lys73 in the presence of substrates, both in the Henri-Michaelis complex and in the tetrahedral intermediate. In all cases, the pKa of Lys73 was computed to be above 10, showing that it is unlikely to act as a proton abstractor, even when a beta-lactam substrate is bound in the enzyme active site. The pKa of Lys234 is also raised in the tetrahedral intermediate, thus confirming a probable role of this residue in the stabilization of the tetrahedral intermediate. The influence of the beta-lactam carboxylate on the pKa values of the active-site lysines is also discussed.

Enzyme kinetic parameters can differ between different species and isoenzymes for the same catalysed reaction. Computational approaches to calculate enzymatic kinetic parameters from the three-dimensional structures of proteins will be reviewed briefly here. Enzyme kinetic parameters may be derived by modelling and simulating the rate-determining process. An alternative, approximate, but more computationally efficient approach is the comparison of molecular interaction fields for experimentally characterized enzymes and those for which parameters should be determined. A correlation between differences in interaction fields and experimentally determined kinetic parameters can be used to determine parameters for orthologous enzymes from other species. The estimation of enzymatic kinetic parameters is an important step in setting up mathematical models of biochemical pathways in systems biology.

Recent studies provide evidence that wild-type Cu/Zn-superoxide dismutase (SOD1(hWT)) might be an important factor in mutant SOD1-mediated amyotrophic lateral sclerosis (ALS). In order to investigate its functional role in the pathogenesis of ALS, we designed fusion proteins of two SOD1 monomers linked by a polypeptide. We demonstrated that wild-type-like mutants, but not SOD1(G85R) homodimers, as well as mutant heterodimers including SOD1(G85R)-SOD1(hWT) display dismutase activity. Mutant homodimers showed an increased aggregation compared with the corresponding heterodimers in cell cultures and transgenic Caenorhabditis elegans, although SOD1(G85R) heterodimers are more toxic in functional assays. Our data show that (i) toxicity of mutant SOD1 is not correlated to its aggregation potential; (ii) dismutase-inactive mutants form dismutase-active heterodimers with SOD1(hWT); (iii) SOD1(hWT) can be converted to contribute to disease by forming active heterodimers. Therefore, we conclude that toxicity of mutant SOD1 is at least partially mediated through heterodimer formation with SOD1(hWT) in vivo and does not correlate with the aggregation potential of individual mutants.

In protein unfolding simulations, elevated temperature, significantly exceeding the melting temperature Tm, provides an important means to accelerate unfolding to a computationally accessible time range. This procedure is based on the assumption that protein thermal unfolding has Arrhenius behavior and therefore that increasing temperature does not alter the protein unfolding pathways. However, in nature, proteins can show non-Arrhenius behavior and, in practice, overly fast unfolding in high-temperature simulations can result in difficulties in identifying unfolding intermediates and distinguishing their relative stabilities. In this paper, we describe simulations of two WW domains, small protein domains that have a three-stranded β-sheet structure. Simulations were carried out at several temperatures ranging from 300 K to 500 K, starting from folded structures. The results demonstrate the temperature dependence of the unfolding pathways, showing that to obtain unfolding pathways corresponding to those observed in experiments, the elevation of the simulation temperature has to be controlled. Based on trajectory analysis, we proposed a qualitative criterion for judging when an elevated temperature is acceptable or not, namely, that the temperature must be such that the native folded state is sampled substantially before protein unfolding begins. While depending on force field parameters and protein fold complexity, this criterion can be quantified to obtain the upper bound of an "acceptable elevated temperature", which was observed to be dependent on the thermostabilities of the two WW domain proteins.

How can we make the connection between the three-dimensional structures of individual proteins and understanding how complex biological systems involving many proteins work? The modelling and simulation of protein structures can help to answer this question for systems ranging from multimacromolecular complexes to organelles and cells. On one hand, multiscale modelling and simulation techniques are advancing to permit the spatial and temporal properties of large systems to be simulated using atomic-detail structures. On the other hand, the estimation of kinetic parameters for the mathematical modelling of biochemical pathways using protein structure information provides a basis for iterative manipulation of biochemical pathways guided by protein structure. Recent advances include the structural modelling of protein complexes on the genomic level, novel coarse-graining strategies to increase the size of the system and the time span that can be simulated, and comparative molecular field analyses to estimate enzyme kinetic parameters.

The active site of cytochromes P450 is situated deep inside the protein next to the heme cofactor. Consequently, enzyme specificity and kinetics can be influenced by how substrates pass through the protein to access the active site and how products egress from the active site. We previously analysed the channels between the active site and the protein surface in P450 crystal structures available in October 2003 [R.C. Wade, P.J. Winn, I. Schlichting, Sudarko, A survey of active site access channels in cytochromes P450, J. Inorg. Biochem. 98 (2004) 1175-1182]. Since then, 52 new P450 structures have been made available, including entries for ten isozymes for which structures were not previously available. We present an updated survey covering all P450 crystal structures available in March 2006. This survey shows channels not observed earlier in crystal structures, some of which were identified in previous molecular dynamics simulations. The crystal structures demonstrate how some of the channels can merge when the protein structure opens up resulting in a wide cleft to the active site, caused largely by movements of the F-G helix-loop-helix and the B-C loop. Significant differences were observed between the channels in the crystal structures of the mammalian and bacterial enzymes. The multiplicity of channels suggests possibilities for substrate channelling to and from the P450s.

The small guanosine triphosphate (GTP)-binding proteins of the Ras family are involved in many cellular pathways leading to cell growth, differentiation, and apoptosis. Understanding the interaction ofRaswith other proteins is of importance not only for studying signalling mechanisms but also, because of their medical relevance as targets, for anticancer therapy. To study their selectivity and specificity, which are essential to their signal transfer function, we performed COMparative BINding Energy (COMBINE) analysis for 122 differentwild- type andmutant complexes between the Ras proteins, Ras and Rap, and their effectors, Raf and RalGDS. The COMBINE models highlighted the amino acid residues responsible for subtle differences in binding of the same effector to the two different Ras proteins, aswell asmore significant differences in the binding of the two different effectors (RalGDS and Raf) to Ras. The study revealed that E37, D38, and D57 in Ras are nonspecific hot spots at its effector interface, important for stabilization of both the RalGDS– Ras and Raf– Ras complexes. The electrostatic interaction between a GTP analogue and the effector, either Raf or RalGDS, also stabilizes these complexes. The Raf– Ras complexes are specifically stabilized by S39, Y40, and D54, and RalGDS– Ras complexes by E31 and D33. Binding of a small molecule in the vicinity of one of these groups of amino acid residues could increase discrimination between the Raf– Ras and RalGDS– Ras complexes. Despite the different size of the RalGDS– Ras and Raf– Ras complexes, we succeeded in building COMBINE models for one type of complex that were also predictive for the other type of protein complex. Further, using system-specific models trained with only five complexes selected according to the results of principal component analysis, we were able to predict binding affinities for the other mutants of the particular Ras-effector complex. As the COMBINE analysis method is able to explicitly reveal the amino acid residues that have most influence on binding affinity, it is a valuable aid for protein design.

International audience; Conservation of function across families of orthologous enzymes is generally accompanied by conservation of their active site electrostatic potentials. To study the electrostatic conservation in the highly conserved essential enzyme, thymidylate synthase (TS), we conducted a systematic species-based comparison of the electrostatic potential in the vicinity of its active site. Whereas the electrostatics of the active site of TS are generally well conserved, the TSs from minimal organisms do not conform to the overall trend. Since the genomes of minimal organisms have a high thymidine content compared to other organisms, the observation of non-conserved electrostatics was surprising. Analysis of the symbiotic relationship between minimal organisms and their hosts, and the genetic completeness of the thymidine synthesis pathway suggested that TS from the minimal organism Wigglesworthia glossinidia (W.g.b.) must be active. Four residues in the vicinity of the active site of Escherichia coli TS were mutated individually and simultaneously to mimic the electrostatics of W.g.b TS. The measured activities of the E. coli TS mutants imply that conservation of electrostatics in the region of the active site is important for the activity of TS, and suggest that the W.g.b. TS has the minimal activity necessary to support replication of its reduced genome. The electrostatic potential of a protein plays a crucial role in steering ligands to their binding sites, and orienting them correctly for binding 1. In enzymes, the active site electrostatic potential is important for stabilizing the transition state and thereby catalyzing the reaction 2. Therefore, conservation of protein function across a protein family is often accompanied by conservation of the electrostatic potential in the active site region, even though the rest of the protein may lack a conserved electrostatic potential 3,4. Consequently, comparison of protein electrostatic potentials has been employed as a tool to predict protein function and to derive similarities in protein function across protein families 5–7. Optimizing the electrostatic complementarity between a ligand and the binding site of a protein is also an important aspect in drug design 8,9 and may provide a route to gain target selectivity 10 .

Fluorescent labels are often attached to proteins to monitor binding and adsorption processes. Docking simulations for native hen egg white lysozyme (HEWL) and HEWL labeled with fluorescein isothiocyanate show that these adsorb differently on charged surfaces. Attachment of even a small label can significantly change the interaction properties of a protein. Thus, the results of experiments with fluorescently labeled proteins should be interpreted by modeling the structures and computing the interaction properties of both labeled and unlabeled species.

Background Intellectual disability (ID) is a neurodevelopmental disorder affecting 1%–3% of the population worldwide. It is characterised by high phenotypic and genetic heterogeneity and in most cases the underlying cause of the disorder is unknown. In our study we investigated a large consanguineous family from Baluchistan, Pakistan, comprising seven affected individuals with a severe form of autosomal recessive ID (ARID) and epilepsy, to elucidate a putative genetic cause. Methods and results Whole exome sequencing (WES) of a trio, including a child with ID and epilepsy and its healthy parents that were part of this large family, revealed a homozygous missense variant p.R53Q in the lectin mannose-binding 2-like ( LMAN2L ) gene. This homozygous variant was co-segregating in the family with the phenotype of severe ID and infantile epilepsy; unaffected family members were heterozygous variant carriers. The variant was predicted to be pathogenic by five different in silico programmes and further three-dimensional structure modelling of the protein suggests that variant p.R53Q may impair protein–protein interaction. LMAN2L (OMIM: 609552) encodes for the lectin, mannose-binding 2-like protein which is a cargo receptor in the endoplasmic reticulum important for glycoprotein transport. Genome-wide association studies have identified an association of LMAN2L to different neuropsychiatric disorders. Conclusion This is the first report linking LMAN2L to a phenotype of severe ARID and seizures, indicating that the deleterious homozygous p.R53Q variant very likely causes the disorder.

Macromolecular interactions play a crucial role in biological systems. Simulation of diffusional association (SDA) is a software for carrying out Brownian dynamics simulations that can be used to study the interactions between two or more biological macromolecules. webSDA allows users to run Brownian dynamics simulations with SDA to study bimolecular association and encounter complex formation, to compute association rate constants, and to investigate macromolecular crowding using atomically detailed macromolecular structures. webSDA facilitates and automates the use of the SDA software, and offers user-friendly visualization of results. webSDA currently has three modules: 'SDA docking' to generate structures of the diffusional encounter complexes of two macromolecules, 'SDA association' to calculate bimolecular diffusional association rate constants, and 'SDA multiple molecules' to simulate the diffusive motion of hundreds of macromolecules. webSDA is freely available to all users and there is no login requirement. webSDA is available at http://mcm.h-its.org/webSDA/.

The conformational properties of unbound multi-Cys2 His2 (mC2H2) zinc finger proteins, in which zinc finger domains are connected by flexible linkers, are studied by a multiscale approach. Three methods on different length scales are utilized. First, atomic detail molecular dynamics simulations of one zinc finger and its adjacent flexible linker confirmed that the zinc finger is more rigid than the flexible linker. Second, the end-to-end distance distributions of mC2H2 zinc finger proteins are computed using an efficient atomistic pivoting algorithm, which only takes excluded volume interactions into consideration. The end-to-end distance distribution gradually changes its profile, from left-tailed to right-tailed, as the number of zinc fingers increases. This is explained by using a worm-like chain model. For proteins of a few zinc fingers, an effective bending constraint favors an extended conformation. Only for proteins containing more than nine zinc fingers, is a somewhat compacted conformation preferred. Third, a mesoscale model is modified to study both the local and the global conformational properties of multi-C2H2 zinc finger proteins. Simulations of the CCCTC-binding factor (CTCF), an important mC2H2 zinc finger protein for genome spatial organization, are presented.

Protein aggregates are the hallmark of stressed and ageing cells, and characterize several pathophysiological states1, 2. Healthy metazoan cells effectively eliminate intracellular protein aggregates3, 4, indicating that efficient disaggregation and/or degradation mechanisms exist. However, metazoans lack the key heat-shock protein disaggregase HSP100 of non-metazoan HSP70-dependent protein disaggregation systems5, 6, and the human HSP70 system alone, even with the crucial HSP110 nucleotide exchange factor, has poor disaggregation activity in vitro4, 7. This unresolved conundrum is central to protein quality control biology. Here we show that synergic cooperation between complexed J-protein co-chaperones of classes A and B unleashes highly efficient protein disaggregation activity in human and nematode HSP70 systems. Metazoan mixed-class J-protein complexes are transient, involve complementary charged regions conserved in the J-domains and carboxy-terminal domains of each J-protein class, and are flexible with respect to subunit composition. Complex formation allows J-proteins to initiate transient higher order chaperone structures involving HSP70 and interacting nucleotide exchange factors. A network of cooperative class A and B J-protein interactions therefore provides the metazoan HSP70 machinery with powerful, flexible, and finely regulatable disaggregase activity and a further level of regulation crucial for cellular protein quality control. published

The experimentally determined structures of human dipeptidyl peptidase III (DPP III) for the wild-type protein and for the complex of its E451A mutant with the peptide substrate, tynorphin, differ significantly in their overall shape. The two domains of the enzyme are separated by a wide cleft in the structure of the ligand-free enzyme, while in the ligand-bound mutant they are very close to each other, and the protein structure is extremely compact. Here, we applied a range of molecular dynamics simulation techniques to investigate the DPP III conformational landscape and the influence of ligand binding on the protein structure and dynamics. We used conventional, accelerated and steered methods to simulate DPP III and its complexes with tynorphin and with the preferred, synthetic, substrate Arg-Arg-2-naphthylamide. We found that DPP III can adopt a number of different forms in solution. The compact forms are more stable, but the open and partially closed states, spanning a wide range of conformations, can more effectively recognize the substrate which preferentially binds to the five-stranded β-core of the lower DPP III domain. The simulations indicated the existence of a dynamic equilibrium between open and semi-closed states and revealed two ways that the protein can close, leading to two distinct compact structures. The way in which the protein closes depends on the presence of the ligand.

We present an analysis of trajectories from Brownian dynamics simulations of diffusional protein-protein encounter for the well-studied system of barnase and barstar. This analysis reveals details about the optimal association pathways, the regions of the encounter complex, possible differences of the pathways for dissociation and association, the coupling of translational and rotation motion, and the effect of mutations on the trajectories. We found that a small free-energy barrier divides the energetically most favorable region into a region of the encounter complex above the barnase binding interface and a region around a second energy minimum near the RNA binding loop. When entering the region of the encounter complex from the region near the RNA binding loop, barstar has to change its orientation to increase the electrostatic attraction between the proteins. By concentrating the analysis on the successful binding trajectories, we found that the region of the second minimum is not essential for the binding of barstar to barnase. Nevertheless, this region may be helpful to steer barstar into the region of the encounter complex. When applying the same analysis to several barnase mutants, we found that single mutations may drastically change the free-energy landscape and may significantly alter the population of the two minima. Therefore, certain protein-protein pairs may require careful adaptation of the positions of encounter and transition states when interpreting mutation effects on kinetic rates of association and/or dissociation.

The oxidation of melatonin by the mammalian myeloperoxidase (MPO) provides protection against the damaging effects of reactive oxygen species. Indole derivatives, such as melatonin and serotonin, are also substrates of the plant horseradish peroxidase (HRP), but this enzyme exhibits remarkable differences from MPO in the specificity and reaction rates for these compounds. A structural understanding of the determinants of the reactivity of these enzymes to indole derivatives would greatly aid their exploitation for biosynthetic and drug design applications. Consequently, after validation of the docking procedure, we performed computational docking of melatonin and serotonin to structural models of the ferric and compound I and II (co I and co II, respectively) states of HRP and MPO. The substrates dock at the heme edge on the distal side, but with different orientations in the two proteins. The distal cavity is larger in MPO than in HRP; however, in MPO, the substrates make closer contacts with the heme involving ring stacking, whereas in HRP, no ring stacking is observed. The observed differences in substrate binding may contribute to the higher reaction rates and lower substrate specificity of MPO relative to those of HRP. The docking results, along with the previously measured heme-protein reduction potentials, suggest that the differentially lowered reaction rates of co II of HRP and MPO with respect to those of co I could stem from as yet undetermined conformational or electrostatic differences between the co I and co II states of MPO, which are absent in HRP.

We report on our work with vibrational absorption, vibrational circular dichroism, Raman scattering, Raman optical activity, and surface-enhanced Raman spectroscopy to study protein and DNA structure, hydration, and the binding of ligands, drugs, pesticides, or herbicides via a combined theoretical and experimental approach. The systems we have studied systematically are the amino acids (L-alanine, L-tryptophan, and L-histidine), peptides (N-4271 acetyl L-alanine N′-methyl amide, N-acetyl L-tryptophan N′-methyl amide, N-acetyl L-histidine N′-methyl amide, L-alanyl L-alanine, tri-L-serine, N-acetyl L-alanine L-proline L-tyrosine N′-methyl amide, Leu-enkephalin, cyclo-(gly-L-pro)3, N-acetyl (L-alanine)nN′-methyl amide), 3-methyl indole, and a variety of small molecules (dichlobenil and 2,6-dochlorobenzamide) of relevance to the protein systems under study. We have used molecular mechanics, the SCC-DFTB, SCC-DFTB+disp, RHF, MP2, and DFT methodologies for the modeling studies with the goal of interpreting the experimentally measured vibrational spectra for these molecules to the greatest extent possible and to use this combined approach to understand the structure, function, and electronic properties of these molecules in their various environments. The application of these spectroscopies to biophysical and environmental assays is expanding, and therefore a thorough understanding of the phenomenon from a rigorous theoretical basis is required. In addition, we give some exciting and new preliminary results which allow us to extend our methods to even larger and more complex systems. The work presented here is the current state of the art to this ever and fast changing field of theoretical spectroscopic interpretation and use of VA, VCD, Raman, ROA, EA, and ECD spectroscopies. © 2005 Wiley Periodicals, Inc. Int J Quantum Chem, 2006

Understanding substrate binding and product release in cytochrome P450 (CYP) enzymes is important for explaining their key role in drug metabolism, toxicity, xenobiotic degradation and biosynthesis. Here, molecular simulations of substrate and product exit from the buried active site of a mammalian P450, the microsomal CYP2C5, identified a dominant exit channel, termed pathway (pw) 2c. Previous simulations with soluble bacterial P450s showed a different dominant egress channel, pw2a. Combining these, we propose two mechanisms in CYP2C5: (i) a one-way route by which lipophilic substrates access the enzyme from the membrane by pw2a and hydroxylated products egress along pw2c; and (ii) a two-way route for access and egress, along pw2c, for soluble compounds. The proposed differences in substrate access and product egress routes between membrane-bound mammalian P450s and soluble bacterial P450s highlight the adaptability of the P450 fold to the requirements of differing cellular locations and substrate specificity profiles.

Accounting for protein flexibility in protein-protein docking algorithms is challenging, and most algorithms therefore treat proteins as rigid bodies or permit side-chain motion only. While the consequences are obvious when there are large conformational changes upon binding, the situation is less clear for the modest conformational changes that occur upon formation of most protein-protein complexes. We have therefore studied the impact of local protein flexibility on protein-protein association by means of rigid body and torsion angle dynamics simulation. The binding of barnase and barstar was chosen as a model system for this study, because the complexation of these 2 proteins is well-characterized experimentally, and the conformational changes accompanying binding are modest. On the side-chain level, we show that the orientation of particular residues at the interface (so-called hotspot residues) have a crucial influence on the way contacts are established during docking from short protein separations of approximately 5 A. However, side-chain torsion angle dynamics simulations did not result in satisfactory docking of the proteins when using the unbound protein structures. This can be explained by our observations that, on the backbone level, even small (2 A) local loop deformations affect the dynamics of contact formation upon docking. Complementary shape-based docking calculations confirm this result, which indicates that both side-chain and backbone levels of flexibility influence short-range protein-protein association and should be treated simultaneously for atomic-detail computational docking of proteins.

The E2 enzymes are key enzymes in the ubiquitin and ubiquitin-like protein ligation pathways. To understand the functionality of the different E2 enzymes, we analyzed 190 protein sequences and 211 structures and electrostatic potentials. Key findings include: The ScUbc1 orthologs are defined by a C-terminal UBA domain. An N-terminal sequence motif that is highly conserved in all E2s except for Cdc34 orthologs is important for the stabilization of the L7 loop and is likely to be involved in E1 binding. ScUbc11p has a different electrostatic potential from E2-Cp and other proteins with which it has high sequence similarity but different functionality. All the E2s known to ubiquitinate histones have a negative potential. The members of the NCUBE family have a positive electrostatic potential, although its form is different from that of the SUMO conjugating E2s. The specificities of only the ScUbc4/Ubc5 and ScUbc1p orthologs are reflected in their L4 and L7 loops.

Structure-based drug design of species-specific inhibitors generally exploits structural differences in proteins from different organisms. Here, we demonstrate how achieving specificity can be aided by targeting differences in the dynamics of proteins. Thymidylate synthase (TS) is a good target for anticancer agents and a potential target for antibacterial agents. Most inhibitors are folate-analogs that bind at the folate binding site and are not species specific. In contrast, α156 is not a folate-analog and is specific for bacterial TS; it has been shown crystallographically to bind in a nonconserved binding site. Docking calculations and crystal structure-based estimation of the essential dynamics of TSs from five different species show that differences in the dynamics of TSs make the active site more accessible to α156 in the prokaryotic than in the eukaryotic TSs and thereby enhance the specificity of α156.

A major problem in mathematical modeling of the dynamics of complex biological systems is the frequent lack of knowledge of kinetic parameters. Here, we apply Brownian dynamics simulations, based on protein three-dimensional structures, to estimate a previously undetermined kinetic parameter, which is then used in biochemical network simulations. The peroxidase-oxidase reaction involves many elementary steps and displays oscillatory dynamics important for immune response. Brownian dynamics simulations were performed for three different peroxidases to estimate the rate constant for one of the elementary steps crucial for oscillations in the peroxidase-oxidase reaction, the association of superoxide with peroxidase. Computed second-order rate constants agree well with available experimental data and permit prediction of rate constants at physiological conditions. The simulations show that electrostatic interactions depress the rate of superoxide association with myeloperoxidase, bringing it into the range necessary for oscillatory behavior in activated neutrophils. Such negative electrostatic steering of enzyme-substrate association presents a novel control mechanism and lies in sharp contrast to the electrostatically-steered fast association of superoxide and Cu/Zn superoxide dismutase, which is also simulated here. The results demonstrate the potential of an integrated and concerted application of structure-based simulations and biochemical network simulations in cellular systems biology.

The activity of a biological compound is dependent both on specific binding to a target receptor and its ADME (Absorption, Distribution, Metabolism, Excretion) properties. A challenge to predict biological activity is to consider both contributions simultaneously in deriving quantitative models. We present a novel approach to derive QSAR models combining similarity analysis of molecular interaction fields (MIFs) with prediction of logP and/or logD. This new classification method is applied to a set of about 100 compounds related to the auxin plant hormone. The classification based on similarity of their interaction fields is more successful for the indole than the phenoxy compounds. The classification of the phenoxy compounds is however improved by taking into account the influence of the logP and/or the logD values on biological activity. With the new combined method, the majority (8 out of 10) of the previously misclassified derivatives of phenoxy acetic acid are classified in accord with their bioassays. The recently determined crystal structure of the auxin-binding protein 1 (ABP1) enabled validation of our approach. The results of docking a few auxin related compounds with different biological activity to ABP1 correlate well with the classification based on similarity of MIFs only. Biological activity is, however, better predicted by a combined similarity of MIFs + logP/logD approach.

We evaluate the applicability of automated molecular docking techniques and quantum mechanical calculations to the construction of a set of structures of enzyme-substrate complexes for use in Comparative binding energy (COMBINE) analysis to obtain 3D structure-activity relationships. The data set studied consists of the complexes of eighteen substrates docked within the active site of haloalkane dehalogenase (DhlA) from Xanthobacter autotrophicus GJ10. The results of the COMBINE analysis are compared with previously reported data obtained for the same dataset from modelled complexes that were based on an experimentally determined structure of the DhlA-dichloroethane complex. The quality of fit and the internal predictive power of the two COMBINE models are comparable, but better external predictions are obtained with the new approach. Both models show a similar composition of the principal components. Small differences in the relative contributions that are assigned to important residues for explaining binding affinity differences can be directly linked to structural differences in the modelled enzyme-substrate complexes: (i) rotation of all substrates in the active site about their longitudinal axis, (ii) repositioning of the ring of epihalohydrines and the halogen substituents of 1,2-dihalopropanes, and (iii) altered conformation of the long-chain molecules (halobutanes and halohexanes). For external validation, both a novel substrate not included in the training series and two different mutant proteins were used. The results obtained can be useful in the future to guide the rational engineering of substrate specificity in DhlA and other related enzymes.

The suitability of three implicit solvent models for flexible protein-protein docking by procedures using molecular dynamics simulation is investigated. The three models are (i) the generalized Born (GB) model implemented in the program AMBER6.0; (ii) a distance-dependent dielectric (DDD) model; and (iii) a surface area-dependent model that we have parameterized and call the NPSA model. This is a distance-dependent dielectric model modified by neutralizing the ionizable side-chains and adding a surface area-dependent solvation term. These solvent models were first tested in molecular dynamics simulations at 300 K of the native structures of barnase, barstar, segment B1 of protein G, and three WW domains. These protein structures display a range of secondary structure contents and stabilities. Then, to investigate the performance of the implicit solvent models in protein docking, molecular dynamics simulations of barnase/barstar complexation, as well as PIN1 WW domain/peptide complexation, were conducted, starting from separated unbound structures. The simulations show that the NPSA model has significant advantages over the DDD and GB models in maintaining the native structures of the proteins and providing more accurate docked complexes.

The three-dimensional structures of proteins are being solved apace, yet this information is often underused in quantitative structure-activity relationship (QSAR) studies. Here, we describe and compare methods for exploiting protein structures to derive 3D-QSARs. These methods can facilitate molecular design and lead optimization and should increasingly become a standard component of the drug designer's repertoire.

Many signaling events require the binding of cytoplasmic proteins to cell membranes by recognition of specific charged lipids, such as phosphoinositol-phosphates. As a model for a protein-membrane binding site, we consider one charged phosphoinositol phosphate (PtdIns(3)P) embedded in a phosphatidylcholine bilayer. As the protein-membrane binding is driven by electrostatic interactions, continuum solvent models require an accurate representation of the electrostatic potential of the phosphoinositol phosphate-containing membrane. We computed and analyzed the electrostatic potentials of snapshots taken at regular intervals from molecular dynamics simulations of the bilayer. We observe considerable variation in the electrostatic potential of the bilayer both along a single simulation and between simulations performed with the GAFF or CHARMM c36 force fields. However, we find that the choice of GAFF or CHARMM c36 parameters has little effect on the electrostatic potential of a given configuration of the bilayer with a PtdIns(3)P embedded in it. From our results, we propose a remedian averaging method for calculating the electrostatic potential of a membrane system that is suitable for simulations of protein-membrane binding with a continuum solvent model.

Despite high similarity in sequence and catalytic properties, the l-lactate dehydrogenases (LDHs) in lactic acid bacteria (LAB) display differences in their regulation that may arise from their adaptation to different habitats. We combined experimental and computational approaches to investigate the effects of fructose 1,6-bisphosphate (FBP), phosphate (Pi), and ionic strength (NaCl concentration) on six LDHs from four LABs studied at pH 6 and pH 7. We found that 1) the extent of activation by FBP (Kact) differs. Lactobacillus plantarum LDH is not regulated by FBP, but the other LDHs are activated with increasing sensitivity in the following order: Enterococcus faecalis LDH2

This paper presents a quantum description of L -alanine in aqueous solution predicting the Raman and Raman optical activity (ROA) spectra. We have investigated theoretically the chiral sensitive spectroscopic method of ROA as a probe for secondary structural features. We have utilized the Becke 3LYP/6-31G* description to determine the geometries and Hessians [Theoretical and Computational Genome Research, Plenum Press, New York, 1997, p. 225; J. Phys. Chem. B 102 (1998) 5899; Chem. Phys. 225 (2000) 165], and have calculated the Becke 3LYP/6-31G* level electric dipole–electric dipole polarizability derivatives. The electric dipole–magnetic dipole polarizability derivatives and electric dipole–electric quadrupole polarizability derivatives further required to obtain the ROA spectra have also been calculated. The Raman scattering and ROA spectral simulations for the various structures are compared with the experimentally measured spectra and previously reported spectral simulations. The combination of Raman and ROA spectroscopy is shown to be as sensitive to secondary structural changes as its sister combination of vibrational absorption and vibrational circular dichroism and thus offers a complementary source of information when investigating the secondary structural and chiral features of biomolecules.

We describe how similarity indices may be applied to semiquantitatively compare the molecular interaction fields, such as electrostatic and hydrophobic potentials, of proteins. We then review ways in which this approach has been used to investigate protein interaction properties in diverse applications ranging from in-depth comparison of pairs of proteins to large-scale analysis of hundreds of modeled protein structures. © 2001 John Wiley & Sons, Inc. Int J Quant Chem 83: 122–127, 2001

Cytochrome P450s form a ubiquitous protein family with functions including the synthesis and degradation of many physiologically important compounds and the degradation of xenobiotics. Cytochrome P450cam from Pseudomonas putida has provided a paradigm for the structural understanding of cytochrome P450s. However, the mechanism by which camphor, the natural substrate of cytochrome P450cam, accesses the buried active site is a long-standing puzzle. While there is recent crystallographic and simulation evidence for opening of a substrate-access channel in cytochrome P450BM-3, for cytochrome P450cam, no such conformational changes have been observed either in different crystal structures or by standard molecular dynamics simulations. Here, a novel simulation method, random expulsion molecular dynamics, is presented, in which substrate-exit channels from the buried active site are found by imposing an artificial randomly oriented force on the substrate, in addition to the standard molecular dynamics force field. The random expulsion molecular dynamics method was tested in simulations of the substrate-bound structure of cytochrome P450BM-3, and then applied to complexes of cytochrome P450cam with different substrates and with product. Three pathways were identified, one of which corresponds to a channel proposed earlier on the basis of crystallographic and site-directed mutagenesis data. Exit via the water-filled channel, which was previously suggested to be a product exit channel, was not observed. The pathways obtained by the random expulsion molecular dynamics method match well with thermal motion pathways obtained by an analysis of crystallographic B-factors. In contrast to large backbone motions (up to 4 A) observed in cytochrome P450BM-3 for the exit of palmitoleic acid, passage of camphor through cytochrome P450cam only requires small backbone motions (less than 2.4 A) in conjunction with side-chain rotations. Concomitantly, in almost all the exit trajectories, salt-links that have been proposed to act as ionic tethers between secondary structure elements of the protein, are perturbed.

Three possible channels by which substrates and products can exit from the buried active site of cytochrome P450cam have been identified by means of random expulsion molecular dynamics simulations. In the investigation described here, we computed estimates of the relative probabilities of ligand passage through the three channels using steered molecular dynamics and adiabatic mapping. For comparison, the same techniques are also applied to investigate substrate egress from cytochrome P450-BM3. The channel in cytochrome P450cam, for which there is the most supporting evidence from experiments (which we name pathway 2a), is computed to be the most probable ligand exit channel. It has the smallest computed unbinding work and force. For this channel, the ligand exits between the F/G loop and the B' helix. Two mechanistically distinct, but energetically similar routes through this channel were observed, showing that multiple pathways along one channel are possible. The probability of ligand exit via the next most probable channel (pathway 3), which is located between the I helix and the F and G helices, is estimated to be less than 1/10 of the probability of exit along pathway 2a. Low-frequency modes of the protein extracted from an essential dynamics analysis of a 1 ns duration molecular dynamics simulation of cytochrome P450cam with camphor bound, support the opening of pathway 2a on a longer timescale. On longer timescales, it is therefore expected that this pathway becomes more dominant than estimated from the present computations.

The Poisson-Boltzmann method was used to compute the pKa values of titratable residues in a set of class C β-lactamases. In these calculations, the pKa of the phenolic group of residue Tyr150 is the only one to stand out with an abnormally low value of 8.3, more than one pKa unit lower than the measured reference value for tyrosine in solution. Other important residues of the catalytic pocket, such as the conserved Lys67, Lys315, His314, and Glu272 (hydrogen-bonded to the ammonium group of Lys315), display normal protonation states at neutral pH. pKa values were also computed in catalytically impaired β-lactamase mutants. Comparisons between the relative kcat values and the Tyr150 pKa value in these mutants revealed a striking correlation. In active enzymes, this pKa value is always lower than the solution reference value while it is close to normal in inactive enzymes. These results thus support the hypothesis that the phenolate form of Tyr150 is responsible for the activation of the nucleophilic serine. The possible roles of Lys67 and Lys315 during catalysis are also discussed. Proteins 2000;40:23–28. © 2000 Wiley-Liss, Inc.

Specific binding of transcription factors to DNA is crucial for gene regulation. We derived models for the binding specificity of transcription factors of the nuclear receptor family to DNA using two QSAR methods: a Free-Wilson-like method and COMparative BINding Energy (COMBINE) analysis. The analysis is based on experimental data for the interaction of 20 mutant glucocorticoid receptor DNA-binding domains with 16 different response elements in a total of 320 complexes (Zilliacus, J.; Wright, A. P.; Carlstedt-Duke, J.; Nilsson, L.; Gustafsson, J. A. Proteins 1995, 21, 57-67). The predictive abilities of the models obtained by the two methods are similar. The COMBINE analysis indicates that the most important properties for determining binding specificity for this dataset are the changes upon binding of the solvation free energies of the bases that are mutated in the dataset and the electrostatic interactions of the mutated nucleotides with certain charged amino acids. Further important descriptors are the changes of solvation free energy and surface area of the side chain of the mutated residue. It is clear, however, that there are additional features important for the specificity of binding that are not included in the models, such as differences in interfacial hydration of the complexes.

Cisproline(i−1)-aromatic(i) interactions have been detected in several short peptides in aqueous solution by analysis of anomalous chemical shifts measured by 1H-NMR spectroscopy. This formation of local structure is of importance for protein folding and binding properties. To obtain an atomic-detail characterisation of the cisproline(i−1)-aromatic(i) interaction in terms of structure, energetics and dynamics, we studied the minimal peptide unit, blocked Ala-cisPro-Tyr, using computational and experimental techniques. Structural database analyses and a systematic search revealed two groups of conformations displaying a cisproline(i−1)-aromatic(i) interaction. These conformations were taken as seeds for molecular dynamics simulations in explicit solvent at 278 K. During a total of 33.6 ns of simulation, all the `folded' conformations and some `unfolded' states were sampled. 1H- and 13C-chemical shifts and 3J-coupling constants were measured for the Ala-Pro-Tyr peptide. Excellent agreement was found between all the measured and computed NMR properties, showing the good quality of the force field. We find that under the experimental and simulation conditions, the Ala-cisPro-Tyr peptide is folded 90% of the time and displays two types of folded conformation which we denote `a' and `b'. The type a conformations are twice as populated as the type b conformations. The former have the tyrosine ring interacting with the alanine α proton and are enthalpically stabilised. The latter have the aromatic ring interacting with the proline side chain and are entropically stabilised. The combined and complementary use of computational and experimental techniques permitted derivation of a detailed scenario of the `folding' of this peptide.

Protein electrostatics plays a key role in ligand binding and protein–protein interactions. Therefore, similarities or dissimilarities in electrostatic potentials can be used as indicators of similarities or dissimilarities in protein function. We here describe a method to compare the electrostatic properties within protein families objectively and quantitatively. Three-dimensional structures are built from database sequences by comparative modeling. Molecular potentials are then computed for these with a continuum solvation model by finite difference solution of the Poisson-Boltzmann equation or analytically as a multipole expansion that permits rapid comparison of very large datasets. This approach is applied to 104 members of the Pleckstrin homology (PH) domain family. The deviation of the potentials of the homology models from those of the corresponding experimental structures is comparable to the variation of the potential in an ensemble of structures from nuclear magnetic resonance data or between snapshots from a molecular dynamics simulation. For this dataset, the results for analysis of the full electrostatic potential and the analysis using only monopole and dipole terms are very similar. The electrostatic properties of the PH domains are generally conserved despite the extreme sequence divergence in this family. Notable exceptions from this conservation are seen for PH domains linked to a Dbl homology (DH) domain and in proteins with internal PH domain repeats. Proteins 1999;37:379–387. ©1999 Wiley-Liss, Inc.

Electrostatic interactions play a key role in many aspects of protein engineering. Consequently, much effort has been put into the design of software for calculating electrostatic fields around macromolecules. We show that optimization of hydrogen bonding networks can improve both the results of pK(a) calculations and the results of electrostatic calculations performed by commonly used programs such as DelPhi. Further optimization can often be achieved by flipping the side chains of asparagine, histidine and glutamine around their chi2, chi2 and chi3 torsion angles, respectively, when this improves the local hydrogen bonding network. These optimizations are applied to some well characterized proteins: BPTI, hen egg white lysozyme and superoxide dismutase. A search for flipped residues in the PDB revealed that significant improvements in electrostatic calculations in or near the active site of enzymes can be expected for about one quarter of all enzymes in the PDB.

When two proteins diffuse together to form a bound complex, an intermediate is formed at the end-point of diffusional association which is called the encounter complex. Its characteristics are important in determining association rates, yet its structure cannot be directly observed experimentally. Here, we address the problem of how to construct the ensemble of three-dimensional structures which constitute the protein–protein diffusional encounter complex using available experimental data describing the dependence of protein association rates on mutation and on solvent ionic strength and viscosity. The magnitude of the association rates is fitted well using a variety of definitions of encounter complexes in which the two proteins are located at up to about 17 A root-mean-squared distance from their relative arrangement in the bound complex. Analysis of the ionic strength dependence of bimolecular association rates shows that this is determined to a greater extent by the (protein charge) – (salt ion) separation distance than by the protein–protein charge separation distance. Consequently, ionic strength dependence of association rates provides little information about the geometry of the encounter complex. On the other hand, experimental data on electrostatic rate enhancement, mutation and viscosity dependence suggest a model of the encounter complex in which the two proteins form a subset of the contacts present in the bound complex and are significantly desolvated. Copyright © 1999 John Wiley & Sons, Ltd.

Heparin−protein interactions play an important role in many steps of the immune system. Here, we evaluated the search capabilities of three widely used programsGRID, DOCK, and AutoDockfor heparin binding sites. Because of the weak surface complementarity and the high charge density of the sulfated sugar chain, the docking of heparin to its protein partners presents a challenging task for computational docking. Our protocols were tested on antithrombin and acidic and basic fibroblast growth factor, the only three proteins for which structures of their complexes with heparin are available. With all three programs, the heparin binding site for these test cases was determined correctly. We then used these protocols to predict the heparin binding site on Interleukin-8, a chemokine with a central role in the human immune response. The results indicate that His18, Lys20, Arg60, Lys64, and Arg68 in interleukin-8 bind to heparin.

The structures and relative energies for L-alanyl-L-alanine (LALA) in the presence of explicit water molecules have been determined by using the density functional theory (DFT) Becke3–Lee–Yang–Parr functional and the 6-31G* basis set (B3LYP/6-31G*). In aqueous solution the dominant state of LALA is the zwitterionic form, while its neutral form is dominant in vacuo. Attempts to locate or determine a gas-phase zwitterionic species failed. That is, on the B3LYP/6-31G* potential energy surface, there is no barrier to proton transfer from the positively charged ammonium group to the negatively charged carboxylate group, or from the ammonium group to the adjacent carbonyl oxygen and from the amide nitrogen to the carboxylate group. To stabilize the zwitterion, we modelled the system by adding explicit water molecules and by placing the zwitterion within a sphere surrounded by a medium with a dielectric constant of 78.5, that is, within the Onsager continuum model, where the recommended cavity radius is obtained from a solute volume calculation. The zwitterionic species is only stable in the presence of water at the B3LYP/6-31G* level. This makes it imperative to include water molecules to model the zwitterionic species of LALA, peptides and amino acids at the B3LYP/6-31G* level. Finally, the zwitterionic structure stabilized by explicit water molecules has also been modelled within the Onsager theory. Here the Onsager model represents the effects due to the bulk water and the explicit water molecules stand for the effect due to direct H-bonding between the zwitterion and the solvent, that is, the first solvation shell. We used molecular dynamics simulations utilizing the CHARMm force field to produce structural input for the subsequent quantum-mechanical simulations. The structures determined using various methods to model the LALA zwitterionic form in aqueous solution were compared. We were able to find additional stable structures for LALA by adding water molecules and optimizing it which could not be obtained by using the Onsager theory. This shows that one must be careful when using continuum models to study peptides and proteins or other methods which do not take into account the explicit interactions between the solute and the first solvent shell.

Hydration of protein cavities influences protein stability, dynamics, and function. Protein active sites usually contain water molecules that, upon ligand binding, are either displaced into bulk solvent or retained to mediate protein-ligand interactions. The contribution of water molecules to ligand binding must be accounted for to compute accurate values of binding affinities. This requires estimation of the extent of hydration of the binding site. However, it is often difficult to identify the water molecules involved in the binding process when ligands bind on the surface of a protein. Cytochrome P450cam is, therefore, an ideal model system because its substrate binds in a buried active site, displacing partially disordered solvent, and the protein is well characterized experimentally. We calculated the free energy differences for having five to eight water molecules in the active site cavity of the unliganded enzyme from molecular dynamics simulations by thermodynamic integration employing a three-stage perturbation scheme. The computed free energy differences between the hydration states are small (within 12 kJ mol-1) but distinct. Consistent with the crystallographic determination and studies employing hydrostatic pressure, we calculated that, although ten water molecules could in principle occupy the volume of the active site, occupation by five to six water molecules is thermodynamically most favorable.

NMR studies of the YTGP tetrapeptide in aqueous solution show a large structural shift on the glycine peptide proton, indicating that the peptide adopts structures in which the aromatic ring lies close to this proton in an aromatic−(i+2) amine interaction (Kemmink and Creighton, J. Mol. Biol. 1993, 234, 861). Here, molecular dynamics simulations are used to obtain details at the atomic level of the structural properties of this system. To maximize the volume of conformational space searched, 13 separate trajectories of 300 ps were calculated at 300 K, with starting structures chosen from a conformational search, protein crystal structures, and modeling. Analysis shows that together the trajectories sample all the important regions of conformational space. Structures from the trajectories were clustered into conformational states, defined as regions in coordinate space around maxima in the probability distribution. Nine conformational states with significant populations were identified in which the glycine...

Triose phosphate isomerase (TIM) is a diffusion-controlled enzyme whose rate is limited by the diffusional encounter of the negatively charged substrate glyceraldehyde 3-phosphate (GAP) with the homodimeric enzyme's active sites. Translational and orientational steering of GAP toward the active sites by the electrostatic field of chicken muscle TIM has been observed in previous Brownian dynamics (BD) simulations. Here we report simulations of the association of GAP with TIMs from four species with net charges at pH 7 varying from -12e to +12e. Computed second-order rate constants are in good agreement with experimental data. The BD simulations and computation of average Boltzmann factors of substrate–protein interaction energies show that the protein electrostatic potential enhances the rates for all the enzymes. There is much less variation in the computed rates than might be expected on the basis of the net charges. Comparison of the electrostatic potentials by means of similarity indices shows that this is due to conservation of the local electrostatic potentials around the active sites which are the primary determinants of electrostatic steering of the substrate. Proteins 31:406–416, 1998. © 1998 Wiley-Liss, Inc.

The ability to reliably compute accurate protein−ligand binding affinities is crucial to understanding protein−ligand recognition and to structure-based drug design. A ligand's binding affinity is specified by its absolute binding free energy, ΔGbind, the free energy difference between the bound and unbound states. To compute accurate free energy differences by free energy perturbation (FEP), “alchemical” rather than physical processes are usually simulated by molecular dynamics simulations so as to minimize the perturbation to the system. Here, we report a novel “alchemistic” application of the FEP methodology involving a large perturbation. By mutating a ligand with 11 non-hydrogen atoms into six water molecules in the binding site of a protein, we computed a ΔGbind within 3 kJ/mol of the experimental value. This is the first successful example of the computation of ΔGbind for a protein:ligand pair with full treatment of the solvent degrees of freedom.

Although auxins were the first type of plant hormone to be identified, little is known about the molecular mechanism of this important class of plant hormones. We present a classification of a set of about 50 compounds with measured auxin activities, according to their interaction properties. Four classes of compounds were defined: strongly active, weakly active with weak antiauxin behaviour, inactive and inhibitory. All compounds were modeled in two low-energy conformations, 'P' and 'T', so as to obtain the best match to the 'planar' and 'tilted' conformations, respectively, of indole 3-acetic acid. Each set of conformers was superimposed separately using several different alignment schemes. Molecular interaction energy fields were computed for each molecule with five different chemical probes and then compared by computing similarity indices. Similarity analysis showed that the classes are on average distinguishable, with better differentiation achieved for the T conformers than the P conformers. This indicates that the T conformation might be the active one. Further, a screening was developed which could distinguish compounds with auxin activity from inactive compounds and most antiauxins using the T conformers. The classifications rationalize ambiguities in activity data found in the literature and should be of value in predicting the activities of new plant growth substances and herbicides.

The threat of a catastrophic outbreak of influenza is ever present. Vaccines are only partially effective and the two compounds, amantidine and rimantidine, used clinically against influenza A cause side-effects and rapid viral resistance. Recent advances bring hope that specific and potent drugs against influenza may soon be available in the clinic. These compounds were designed to inhibit influenza neuraminidase (NA), one of the viral coat glycoproteins, using the crystal structure of NA which was first published in 1983. In this review, the application of structure-based drug design approaches to the design of anti-influenza agents targeted at NA and haemagglutinin (HA), the other viral surface glycoprotein, is discussed.

Cytochrome P450 enzymes are hemoprotein monooxygenases that catalyse the oxidation of a variety of compounds. The mechanism by which camphor, the natural substrate of Cytochrome P450cam (P450cam), accesses the active site is a long-standing puzzle, although putative access channels have been proposed. A thermal motion pathway (TMP) analysis was performed on the crystal structure of P450cam with camphor bound. Hereby, three distinct thermal motion pathway families (TMPFs) were found. Possible substrate access channels obtained by this analysis based on B-factors are compared with exit channels explored by molecular dynamics simulations (MDS) by imposing an artificial expulsion force on the substrate in addition to the standard MD force field. Two out of three TMPFs are supported by results obtained with the random expulsion MDS method. However, the pathway found by the TMP method to have the highest average B-factor could not be observed by MDS. The pathway proposed from crystallographic data, which is a small opening above the active site located near residues 185, 87 and 395 corresponds to the TMPF with the second highest average B-factor.