In the present study, we have performed combined molecular dynamics and potential of mean force (PMF) simulations to determine the enzyme-substrate (ES) binding pathway and the corresponding free energy profiles for wildtype butyrylcholinesterase (BChE) binding with (-)/(+)-cocaine and for the A328W/Y332G mutant binding with (-)-cocaine. According to the PMF simulations, for each ES binding system, the substrate first binds with the enzyme at a peripheral anionic site around the entrance of the active-site gorge to form the first ES complex (ES1-like) during the binding process. Further evolution from the ES1-like complex to the nonprereactive ES complex is nearly barrierless, with a free energy barrier lower than 1.0 kcal/mol. So, the nonprereactive ES binding process should be very fast. The rate-determining step of the entire ES binding process is the subsequent evolution from the nonprereactive ES complex to the prereactive ES complex. Further accounting for the entire ES binding process, the PMF-based simulations qualitatively reproduced the relative order of the experimentally derived binding free energies (Delta G(bind)), although the simulations systematically overestimated the magnitude of the binding affinity and systematically underestimated the differences between the Delta G(bind) values. The obtained structural and energetic insights into the entire ES binding process provide a valuable base for future rational design of high-activity mutants of BChE as candidates for an enzyme therapy for cocaine overdose and abuse.