An explicitly polarizable force field based exclusively on quantum data is applied to calculations of relative binding affinities of ligands to proteins. Five ligands, differing by replacement of an atom or functional group, in complexes with three serine proteases-trypsin, thrombin, and urokinase-type plasminogen activator-with available experimental binding data are used as test systems. A special protocol of thermodynamic integration was developed and used to provide sufficiently low levels of systematic error along with high numerical efficiency and statistical stability. The calculated results are in excellent quantitative (rmsd ؍ 1.0 kcal/mol) and qualitative (R 2 ؍ 0.90) agreement with experimental data. The potential of the methodology to explain the observed differences in the ligand affinities is also demonstrated. molecular dynamics simulation ͉ serine protease ͉ drug design T he ability to accurately calculate the binding affinity, or equivalently binding free energy, of a ligand for a protein would be highly useful in the field of drug design for lead selection and optimization. Although screening and docking methods (1, 2) have been successful in filtering large chemical databases, they cannot provide definitive calculations of binding energy because of the simplified scoring functions used and the restricted number of states tested. Hence, more accurate calculations have generally been based on molecular mechanics models (3, 4). In these methods the required properties of statistical ensembles are determined by molecular dynamical or Monte Carlo simulations of systems by using a physically grounded description of interactions between particles. The theoretical thermodynamic foundation of such methods is clear, simple, and well established (5).Despite these advantages, as well as initial optimism and long development, only a limited number of successful simulations have been published during the past decade (e.g., refs. 6-14; for review of early results, see refs. 3 and 4). In part, this is explained by many methodological difficulties that must be overcome, especially in accurate and efficient description of long-range interactions and adequate sampling of the conformational space. Major efforts devoted to solving these problems have resulted both in partial success (15, 16) and the recognition of some principal limitations. Notably, several theoretically based techniques have been developed, such as umbrella sampling, the concept of potential of mean force, and artificial restraining potentials, which restrict or decompose the conformational space and simplify adequate sampling (e.g., 11, 17-20).The other principal problem has been the quality of the model potentials or force fields (FFs) used to describe atom-atom interactions. The most questionable point in this respect is the role of nonadditive effects, particularly electronic polarizability. Widely used FFs such as MMFF (21), AMBER (22), OPLS (23), CHARMM (24), and GROMOS (25) are not explicitly polarizable but rather include polarizabil...