In drug discovery process, improvement of ADME/Tox properties of lead compounds including metabolic stability is critically important. Cytochrome P450 (CYP) is one of the major metabolizing enzymes and the prediction of sites of metabolism (SUM) on the given lead compounds is key information to modify the compounds to be more stable against metabolism. There are two factors essentially important in SUM prediction. First is accessibility of each substrate atom to the oxygenated Fe atom of heme in a CYP protein, and the other is the oxidative reactivity of each substrate atom. To predict accessibility of substrate atoms to the heme iron, conventional protein-rigid docking simulations have been applied. However, the docking simulations without consideration of protein flexibility often lead to incorrect answers in the case of very flexible proteins such as CYP3A4. In this study, we demonstrated an approach utilizing molecular dynamics (MD) simulation for SUM prediction in which multiple MD runs were executed using different initial structures. We applied this strategy to CYP3A4 and carbamazepine (CBZ) complex. Through 10 no MD simulations started from five different CYP3A4-CBZ complex models, our approach correctly predicted SUM observed in experiments. The experimentally known epoxidized sites of CBZ by CYP3A4 were successfully predicted as the most accessible sites to the heme iron that was judged from a numerical analysis of calculated Delta G(binding) and the frequency of appearance. In contrast, the predictions using protein-rigid docking methods hardly provided the correct SUM due to protein flexibility or inaccuracy of the scoring functions. Our strategy using MD simulation with multiple initial structures will be one of the reliable methods for SUM prediction. (C) 2011 Elsevier Ltd. All rights reserved.