Dynamics and Allosteric Information Pathways of Unphosphorylated c-Cbl.

Human c-Cbl is a RING-type ligase and plays a central role in the protein degradation cascade. To elucidate its conformational changes related to substrate binding, we performed molecular dynamics simulations of different variants/states of c-Cbl for a cumulative time of 68 µs. Our simulations demonstrate that before the substrate binds, the RING domain samples a broad set of conformational states at a biologically relevant salt concentration, including the closed, partially open, and fully open states, whereas substrate binding leads to a restricted conformational sampling. Phe378 and the C-terminal region play an essential role in stabilizing the partially open state. To visualize the allosteric signal transmission pathways from the substrate-binding site to the 40 Å apart RING domain and identify the critical residues for allostery, we have created a subgraph from the optimal and suboptimal paths. Redundant paths are seen in the SH2 domain where the substrate binds, while the major bottlenecks are found at the junction between the SH2 domain and the linker helix region as well as that between the SH2 domain and the 4H bundle. These bottlenecks separate the paths into two overall routes. The nodes/residues at the bottlenecks on the subgraph are considered allosteric hot spots. This subgraph approach provides a general tool for network visualization and determination of critical residues for allostery. The structurally and allosterically critical residues identified in our work are testable and would provide valuable insights into the emerging strategies for drug discovery, such as targeted protein degradation.

Approximating dynamic proximity with a hybrid geometry energy-based kernel for diffusion maps.

The diffusion map is a dimensionality reduction method. The reduction coordinates are associated with the leading eigenfunctions of the backward Fokker-Planck operator, providing a dynamic meaning for these coordinates. One of the key factors that affect the accuracy of diffusion map embedding is the dynamic measure implemented in the Gaussian kernel. A common practice in diffusion map study of molecular systems is to approximate dynamic proximity with RMSD (root-mean-square deviation). In this paper, we present a hybrid geometry-energy based kernel. Since high energy-barriers may exist between geometrically similar conformations, taking both RMSD and energy difference into account in the kernel can better describe conformational transitions between neighboring conformations and lead to accurate embedding. We applied our diffusion map method to the ß-hairpin of the B1 domain of streptococcal protein G and to Trp-cage. Our results in ß-hairpin show that the diffusion map embedding achieves better results with the hybrid kernel than that with the RMSD-based kernel in terms of free energy landscape characterization and a new correlation measure between the cluster center Euclidean distances in the reduced-dimension space and the reciprocals of the total net flow between these clusters. In addition, our diffusion map analysis of the ultralong molecular dynamics trajectory of Trp-cage has provided a unified view of its folding mechanism. These promising results demonstrate the effectiveness of our diffusion map approach in the analysis of the dynamics and thermodynamics of molecular systems. The hybrid geometry-energy criterion could be also useful as a general dynamic measure for other purposes.

Inherent structure versus geometric metric for state space discretization.

Inherent structure (IS) and geometry-based clustering methods are commonly used for analyzing molecular dynamics trajectories. ISs are obtained by minimizing the sampled conformations into local minima on potential/effective energy surface. The conformations that are minimized into the same energy basin belong to one cluster. We investigate the influence of the applications of these two methods of trajectory decomposition on our understanding of the thermodynamics and kinetics of alanine tetrapeptide. We find that at the microcluster level, the IS approach and root-mean-square deviation (RMSD)-based clustering method give totally different results. Depending on the local features of energy landscape, the conformations with close RMSDs can be minimized into different minima, while the conformations with large RMSDs could be minimized into the same basin. However, the relaxation timescales calculated based on the transition matrices built from the microclusters are similar. The discrepancy at the microcluster level leads to different macroclusters. Although the dynamic models established through both clustering methods are validated approximately Markovian, the IS approach seems to give a meaningful state space discretization at the macrocluster level in terms of conformational features and kinetics.

Long-Range Reactivity Modulations in Geranyl Chloride Derivatives.

Derivatives of geraniol are versatile synthetic intermediates that are useful for synthesizing a variety of terpenoid natural products; however, the results presented herein show that subtle differences in the structures of functionalized geranyl chlorides can significantly impact their abilities to function as effective electrophiles in synthetic reactions. A series of focused kinetics experiments identify specific structure-activity relationships that illustrate the importance not only of steric bulk, but also of electronic effects from distant regions of the molecules that contribute to their overall levels of reactivity. Computational modeling suggests that destabilization of the reactant by filled-filled orbital mixing events in some, but not all, conformations may be a critical contributor to these important electronic effects.

Network representation of conformational transitions between hidden intermediates of Rd-apocytochrome b562.

The folding kinetics of Rd-apocytochrome b562 is two-state, but native-state hydrogen exchange experiments show that there are discrete partially unfolded (PUF) structures in equilibrium with the native state. These PUF structures are called hidden intermediates because they are not detected in kinetic experiments and they exist after the rate-limiting step. Structures of the mimics of hidden intermediates of Rd-apocytochrome b562 are resolved by NMR. Based upon their relative stability and structural features, the folding mechanism was proposed to follow a specific pathway (unfolded → rate-limiting transition state → PUF1 → PUF2 → native). Investigating the roles of equilibrium PUF structures in folding kinetics and their interrelationship not only deepens our understanding of the details of folding mechanism but also provides guides in protein design and prevention of misfolding. We performed molecular dynamics simulations starting from a hidden intermediate and the native state of Rd-apocytochrome b562 in explicit solvent, for a total of 37.18 µs mainly with Anton. We validated our simulations by detailed comparison with experimental data and other computations. We have verified that we sampled the post rate-limiting transition state region only. Markov state model was used to analyze the simulation results. We replace the specific pathway model with a network model. Transition-path theory was employed to calculate the net effective flux from the most unfolded state towards the most folded state in the network. The proposed sequential folding pathway via PUF1 then more stable, more native-like PUF2 is one of the routes in our network, but it is not dominant. The dominant path visits PUF2 without going through PUF1. There is also a route from PUF1 directly to the most folded state in the network without visiting PUF2. Our results indicate that the PUF states are not necessarily sequential in the folding. The major routes predicted in our network are testable by future experiments such as single molecule experiment.

Euclidean sections of protein conformation space and their implications in dimensionality reduction.

Dimensionality reduction is widely used in searching for the intrinsic reaction coordinates for protein conformational changes. We find the dimensionality-reduction methods using the pairwise root-mean-square deviation (RMSD) as the local distance metric face a challenge. We use Isomap as an example to illustrate the problem. We believe that there is an implied assumption for the dimensionality-reduction approaches that aim to preserve the geometric relations between the objects: both the original space and the reduced space have the same kind of geometry, such as Euclidean geometry vs. Euclidean geometry or spherical geometry vs. spherical geometry. When the protein free energy landscape is mapped onto a 2D plane or 3D space, the reduced space is Euclidean, thus the original space should also be Euclidean. For a protein with N atoms, its conformation space is a subset of the 3N-dimensional Euclidean space R(3N). We formally define the protein conformation space as the quotient space of R(3N) by the equivalence relation of rigid motions. Whether the quotient space is Euclidean or not depends on how it is parameterized. When the pairwise RMSD is employed as the local distance metric, implicit representations are used for the protein conformation space, leading to no direct correspondence to a Euclidean set. We have demonstrated that an explicit Euclidean-based representation of protein conformation space and the local distance metric associated to it improve the quality of dimensionality reduction in the tetra-peptide and ß-hairpin systems.

Graph representation of protein free energy landscape.

The thermodynamics and kinetics of protein folding and protein conformational changes are governed by the underlying free energy landscape. However, the multidimensional nature of the free energy landscape makes it difficult to describe. We propose to use a weighted-graph approach to depict the free energy landscape with the nodes on the graph representing the conformational states and the edge weights reflecting the free energy barriers between the states. Our graph is constructed from a molecular dynamics trajectory and does not involve projecting the multi-dimensional free energy landscape onto a low-dimensional space defined by a few order parameters. The calculation of free energy barriers was based on transition-path theory using the MSMBuilder2 package. We compare our graph with the widely used transition disconnectivity graph (TRDG) which is constructed from the same trajectory and show that our approach gives more accurate description of the free energy landscape than the TRDG approach even though the latter can be organized into a simple tree representation. The weighted-graph is a general approach and can be used on any complex system.

Observation of two families of folding pathways of BBL.

BBL is an independent folding domain of a large multienzyme complex, 2-oxoglutarate dehydrogenase. The folding mechanism of BBL is under debate between the views of noncooperative downhill-type and classical two-state. Extensive replica exchange molecular dynamics simulations of BBL in explicit solvent have shown some non-two-state behaviors despite no definitive evidence of downhill folding. In this work, we postprocess the replica exchange data using our roadmap-based MaxFlux reaction path algorithm to reveal atomically detailed folding pathways. A connected graph is used to organize and visualize the folding pathways initiated from random coils. High structural and transition heterogeneity is seen in the early stage of folding. Two main parallel folding pathways emerge in the later stage; one path shows that tertiary contact and helix formation develop at different stages of folding, whereas the other path exhibits concurrence of secondary and tertiary structure formation to some extent. Because the native state of BBL is sensitive to experimental conditions, we speculate that the relative predominance of the two pathways may vary with the protein construct and solvent conditions, possibly leading to the seeming discrepancy of experimental results. Our roadmap-based reaction path algorithm is a general tool to extract path information from replica exchange.

Interplay between human islet amyloid polypeptide aggregates and micro-heterogeneous membranes.

Human islet amyloid polypeptides (hIAPP) aggregate into amyloid deposits in the pancreatic islets of Langerhans, contributing to the loss of ß-cells of patients with type 2 diabetes. Despite extensive studies of membrane disruption associated with hIAPP aggregates, the molecular details regarding the complex interplay between hIAPP aggregates and raft-containing membranes are still very limited. Using all-atom molecular dynamics simulations, we investigate the impact of hIAPP aggregate insertion on lipid segregation. We have found that the domain separation of 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) and 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) is enhanced upon hIAPP membrane permeabilization in the absence of cholesterol, while within our simulation timescale, we cannot provide definitive evidence regarding the impact of hIAPP insertion on domain segregation in the ternary mixture (DOPC/DPPC/cholesterol). When the lipid domains are perturbed, their restoration occurs rapidly and spontaneously in the presence of hIAPP aggregates. hIAPP insertion affects membrane thickness in its immediate surroundings. On average, hIAPP causes the fluidity of lipids to increase and even cholesterol shows enhanced diffusivity. The acyl chain packing of the lipids near hIAPP is disrupted as compared to that further away from it. Cholesterol not only modulates membrane mobility and ordering but also hIAPP aggregates' structure and relative orientation to the membrane. Our investigations on the interaction between hIAPP aggregates and raft-containing membranes could lead to a better understanding of the mechanisms of amyloid cytotoxicity.

Quantum mechanical studies on model alpha-pleated sheets.

Pauling and Corey proposed a pleated-sheet configuration, now called alpha-sheet, as one of the protein secondary structures in addition to alpha-helix and beta-sheet. Recently, it has been suggested that alpha-sheet is a common feature of amyloidogenic intermediates. We have investigated the stability of antiparallel beta-sheet and two conformations of alpha-sheet in solution phase using the density functional theoretical method. The peptides are modeled as two-strand acetyl-(Ala)(2)-N-methylamine. Using stages of geometry optimization and single point energy calculation at B3LYP/cc-pVTZ//B3LYP/6-31G* level and including zero-point energies, thermal, and entropic contribution, we have found that beta-sheet is the most stable conformation, while the alpha-sheet proposed by Pauling and Corey has 13.6 kcal/mol higher free energy than the beta-sheet. The alpha-sheet that resembles the structure observed in molecular dynamics simulations of amyloidogenic proteins at low pH becomes distorted after stages of geometry optimization in solution. Whether the alpha-sheets with longer chains would be increasingly favorable in water relative to the increase in internal energy of the chain needs further investigation. Different from the quantum mechanics results, AMBER parm94 force field gives small difference in solution phase energy between alpha-sheet and beta-sheet. The predicted amide I IR spectra of alpha-sheet shows the main band at higher frequency than beta-sheet.

Formation and growth of oligomers: a Monte Carlo study of an amyloid tau fragment.

Small oligomers formed early in the process of amyloid fibril formation may be the major toxic species in Alzheimer's disease. We investigate the early stages of amyloid aggregation for the tau fragment AcPHF6 (Ac-VQIVYK-NH2) using an implicit solvent all-atom model and extensive Monte Carlo simulations of 12, 24, and 36 chains. A variety of small metastable aggregates form and dissolve until an aggregate of a critical size and conformation arises. However, the stable oligomers, which are beta-sheet-rich and feature many hydrophobic contacts, are not always growth-ready. The simulations indicate instead that these supercritical oligomers spend a lengthy period in equilibrium in which considerable reorganization takes place accompanied by exchange of chains with the solution. Growth competence of the stable oligomers correlates with the alignment of the strands in the beta-sheets. The larger aggregates seen in our simulations are all composed of two twisted beta-sheets, packed against each other with hydrophobic side chains at the sheet-sheet interface. These beta-sandwiches show similarities with the proposed steric zipper structure for PHF6 fibrils but have a mixed parallel/antiparallel beta-strand organization as opposed to the parallel organization found in experiments on fibrils. Interestingly, we find that the fraction of parallel beta-sheet structure increases with aggregate size. We speculate that the reorganization of the beta-sheets into parallel ones is an important rate-limiting step in the formation of PHF6 fibrils.

Predicting the folding pathway of engrailed homeodomain with a probabilistic roadmap enhanced reaction-path algorithm.

To predict a protein-folding pathway, we present an alternative to the time-consuming dynamic simulation of atomistic models. We replace the actual dynamic simulation with variational optimization of a reaction path connecting known initial and final protein conformations in such a way as to maximize an estimate of the reactive flux or minimize the mean first passage time at a given temperature, referred to as MaxFlux. We solve the MaxFlux global optimization problem with an efficient graph-theoretic approach, the probabilistic roadmap method (PRM). We employed CHARMM19 and the EEF1 implicit solvation model to describe the protein solution. The effectiveness of our MaxFlux-PRM is demonstrated in our promising simulation results on the folding pathway of the engrailed homeodomain. Our MaxFlux-PRM approach provides the direct evidence to support that the previously reported intermediate state is a genuine on-pathway intermediate, and the demand of CPU power is moderate.

A novel mutant cardiac troponin C disrupts molecular motions critical for calcium binding affinity and cardiomyocyte contractility.

Troponin C (TnC) belongs to the superfamily of EF-hand (helix-loop-helix) Ca(2+)-binding proteins and is an essential component of the regulatory thin filament complex. In a patient diagnosed with idiopathic dilated cardiomyopathy, we identified two novel missense mutations localized in the regulatory Ca(2+)-binding Site II of TnC, TnC((E59D,D75Y)). Expression of recombinant TnC((E59D,D75Y)) in isolated rat cardiomyocytes induced a marked decrease in contractility despite normal intracellular calcium homeostasis in intact cardiomyocytes and resulted in impaired myofilament calcium responsiveness in Triton-permeabilized cardiomyocytes. Expression of the individual mutants in cardiomyocytes showed that TnC(D75Y) was able to recapitulate the TnC((E59D,D75Y)) phenotype, whereas TnC(E59D) was functionally benign. Force-pCa relationships in TnC((E59D,D75Y)) reconstituted rabbit psoas fibers and fluorescence spectroscopy of TnC((E59D,D75Y)) labeled with 2-[(4'-iodoacetamide)-aniline]naphthalene-6-sulfonic acid showed a decrease in myofilament Ca(2+) sensitivity and Ca(2+) binding affinity, respectively. Furthermore, computational analysis of TnC showed the Ca(2+)-binding pocket as an active region of concerted motions, which are decreased markedly by mutation D75Y. We conclude that D75Y interferes with proper concerted motions within the regulatory Ca(2+)-binding pocket of TnC that hinders the relay of the thin filament calcium signal, thereby providing a primary stimulus for impaired cardiomyocyte contractility. This in turn may trigger pathways leading to aberrant ventricular remodeling and ultimately a dilated cardiomyopathy phenotype.

Structural and pathway complexity of beta-strand reorganization within aggregates of human transthyretin(105-115) peptide.

Interstrand conformational rearrangements of human transthyretin peptide (TTR(105-115)) within dimeric aggregates were simulated by means of molecular dynamics (MD) with implicit solvation model for a total length of 48 micros. The conformations sampled in the MD simulations were clustered to identify free energy minima without any projections of free energy surface. A connected graph was constructed with nodes (=clusters) and edges corresponding to free energy minima and transitions between nodes, respectively. This connected graph which reflects the complexity of the free energy surface was used to extract the transition disconnectivity graph, which reflects the overall free energy barriers between pairs of free energy minima but does not contain information on transition paths. The routes of transitions between important free energy minima were obtained by further processing the original graph and the MD data. We have found that both parallel and antiparallel aggregates are populated. The parallel aggregates with different alignment patterns are separated by nonnegligible free energy barriers. Multiroutes exist in the interstrand conformational reorganization. Most visited routes do not dominant the kinetics, while less visited routes contribute a little each but they are numerous and their total contributions are actually dominant. There are various kinds of reptation motions, including those through a beta-bulge, side-chain aided reptation, and flipping or rotation of a hairpin formed by one strand.

Evaluation of configurational entropy methods from peptide folding-unfolding simulation.

A 4-micros molecular dynamics simulation of the second beta-hairpin of the B1 domain of streptococcal protein G is used to characterize the free energy surface and to evaluate different configurational entropy estimators. From the equilibrium folding-unfolding trajectory, 200 000 conformers are clustered according to their root-mean-square deviation (RMSD). The height of the free energy barrier between pairs of clusters is found to be significantly correlated with their pairwise RMSD. Relative free energies and relative configurational entropies of the clusters are determined by explicit evaluation of the partition functions of the different clusters. These entropies are used to evaluate different entropy estimators for the largest 20 clusters as well as a subensemble comprising exclusively extended conformers. It is found that the quasi-harmonic entropy estimator operating in dihedral angle space performs better than the one using Cartesian coordinates. A recent generalization of the quasi-harmonic approach that computes Shannon entropies of probability distributions obtained by projecting the conformers along the eigenvectors of the covariance matrix performs similarly well. For the best entropy estimators, a linear correlation coefficient between 0.92 and 0.97 is found. Unexpectedly, when correlations between dihedral angles are neglected, the agreement with the reference entropies improved.

Observations on AMBER Force Field Performance under the Conditions of Low pH and High Salt Concentrations.

Molecular dynamics simulations have become an important tool for the study of structures, dynamics, and functions of biomolecules. Time scales and force field accuracy are key factors for the reliability of these simulations. With the advancement of computational platforms and simulation technologies, all-atom simulations of proteins in explicitly represented aqueous solutions can reach as long as milliseconds. However, there are indications of minor force field flaws in literature. In this work we present our observations on force field accuracy under uncommon conditions. We performed molecular dynamics simulations of BBL refolding in aqueous solutions of acidic pH and high salt concentrations (up to 6 M) with the AMBER99SB-ILDN force field for a microsecond time scale. The reliability of the simulations relies on the accuracy of the physical models of protein, water, and ions. Our simulations show the same trend as experiments: higher salt concentration facilities refolding. However, we have observed the presence of ß-sheet in the native helical region as well as α-helix and ß-sheet in the native loop region. Some of the nonnative secondary structures are even more stable than native helices. Aside from the secondary structure issue under the uncommon conditions, we have found that the rigidity of glycine dihedral angles in the loop region leads to low root-mean-square deviations but large topological differences from the native structure. Whether this is due to a force field deficiency or not needs further investigations. Recently developed ion parameters exhibit evident liquid features even in the 6 M LiCl solution. However, cation-anion interactions in TIP3P water still seem too strong, leading to high fractions of contact ion pairs. We do not find any specific ion-binding motif, thus we conclude that the effect of salt is a nonspecific electrostatic screening in our simulations. Our observations on the AMBER force field performance under acidic conditions and high salt concentrations show that simulations under extreme conditions can provide informative tests of physical models.

Peptide plane can flip in two opposite directions: implication in amyloid formation of transthyretin.

Transthyretin (TTR) is one of the known 20 or so human proteins that form fibrils in vivo, which is a hallmark of amyloid diseases. Recently, molecular dynamics simulations using ENCAD force field have revealed that under low pH conditions, the peptide planes of several amyloidogenic proteins can flip in one direction to form an alpha-pleated structure which may be a common conformational transition in the fibril formation. We performed molecular dynamics simulations with AMBER force fields on a recently engineered double mutant TTR, which was shown experimentally to form amyloid fibrils even under close to physiological conditions. Our simulations have demonstrated that peptide-plane flipping can occur even under neutral pH and room temperature for this amyloidogenic TTR variant. Unlike previously reported peptide-plane flipping of TTR using ENCAD force field, we have found two-way flipping using AMBER force field. We propose a new mechanism of amyloid formation based on the two-way flipping, which gives a better explanation of various experimental and computational results. In principle, the residual dipolar and hydrogen-bond scalar coupling techniques can be applied to the wild-type TTR and the variant to study the peptide-plane flipping of amyloidogenic proteins.

Potential influence of Asp in the Ca2+ coordination position 5 of parvalbumin on the calcium-binding affinity: a computational study.

Parvalbumins (PV) are calcium-binding proteins, all sharing the common helix-loop-helix (EF-hand) motif. This motif contains a central twelve-residue Ca(2+)-binding loop with the flanking helices positioned roughly perpendicular to each other. The precise role of these coordination residues has been the subject of intense studies. In this work, we focus on the coordination position 5 in the CD Ca(2+)-binding site of silver hake parvalbumin isoform B (SHPV-B). The most common residue at site 5 of calcium-binding loop in canonical EF-hands is Asp [B.J. Marsden, G.S. Shaw, B.D. Sykes, Biochem. Cell Biol. 68 (1990) 587-601], but in the CD site of PV, this position is almost always serine (Ser). The substitution of Ser with Asp will add the 5th carboxylate residue in the CD coordination sphere. However, as predicted by the acid pair hypothesis, the Ca(2+)-binding affinity would be maximized in an EF-hand motif that has four carboxylate ligands paired along the +/-x, and +/-z-axes [R.E. Reid, R.S. Hodges, J. Theor. Biol. 84 (1980) 401-444]. Molecular dynamics simulations and free energy calculations were employed to investigate the influence of Ser to Asp mutation at position 5 on calcium-binding affinity. We found that the Asp variant exhibited remarkable stability during the entire molecular dynamics simulation, with not only the retention of the Ca(2+)-binding site, but also increased compactness in the coordination sphere. The S55D fragment also accommodated Ca(2+) well. We conclude that the reason why Asp which is the most common residue at site 5 of calcium-binding loop in canonical EF-hands has never been identified at this position experimentally for PVs might be related to its physiological functions.

Why is Leu55-->Pro55 transthyretin variant the most amyloidogenic: insights from molecular dynamics simulations of transthyretin monomers.

Transthyretin (TTR) is one of the known human amyloidogenic proteins. Its native state is a homotetramer with each monomer having a beta-sandwich structure. Strong experimental evidence suggests that TTR dissociates into monomeric intermediates and that the monomers subsequently self-assemble to form amyloid deposits and insoluble fibrils. However, details on the early steps along the pathway of TTR amyloid formation are unclear, although various experimental approaches with resolutions at the molecular or residue level have provided some clues. It is highly likely that the stability and flexibility of monomeric TTR play crucial roles in the early steps of amyloid formation; thereby, it is essential to characterize initial conformational changes of TTR monomers. In this article we probe the possibility that the differences in the monomeric forms of wild-type (WT) TTR and its variants are responsible for differential amyloidogenesis. We begin with the simulations of WT, Val30-->Met (V30M), and Leu55-->Pro (L55P) TTR monomers. Nanosecond time scale molecular dynamics simulations at 300 K were performed using AMBER. The results indicate that the L55P-TTR monomer undergoes substantial structural changes relative to fluctuations observed in the WT and V30M TTR monomers. The observation supports earlier speculation that the L55P mutation may lead to disruption of the beta-sheet structure through the disorder of the "edge strands" that might facilitate amyloidogenesis.

Molecular dynamics and free energy analyses of cathepsin D-inhibitor interactions: insight into structure-based ligand design.

In this study, we compare the calculated and experimental binding free energies for a combinatorial library of inhibitors of cathepsin D (CatD), an aspartyl protease. Using a molecular dynamics (MD)-based, continuum solvent method (MM-PBSA), we are able to reproduce the experimental binding affinity for a set of seven inhibitors with an average error of ca. 1 kcal/mol and a correlation coefficient of 0.98. By comparing the dynamical conformations of the inhibitors complexed with CatD with the initial conformations generated by CombiBuild (University of California, San Francisco, CA, 1995), we have found that the docking conformation observed in an X-ray structure of one peptide inhibitor bound to CatD (Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 6796-6800) is in good agreement with our MD simulation. However, the DOCK scoring function, based on intermolecular van der Waals and electrostatics, using a distance-dependent dielectric constant (J. Comput. Chem. 1992, 13, 505-524), poorly reproduces the trend of experimental binding affinity for these inhibitors. Finally, the use of PROFEC (J. Comput.-Aided Mol. Des. 1998, 12, 215-227) analysis allowed us to identify two possible substitutions to improve the binding of one of the better inhibitors to CatD. This study offers hope that current methods of estimating the free energy of binding are accurate enough to be used in a multistep virtual screening protocol.