Binding free energies between coenzyme (FAD and NADH) and the apoenzyme of NADH-cytochrome b5 reductase (b5R) were estimated by applying the continuum Poisson-Boltzmann (PB) model to structures sampled from molecular dynamics simulations in explicit water molecules. Important residues for the enzymatic catalysis were clarified using a computational alanine scanning method. The binding free energies calculated by applying an alanine scanning method can successfully reproduce the trends of the measured steady-state enzymatic activities k(cat)(NADH)/K-m(NADH). Significant decreases in the binding free energy are expected when one of the four residues Arg91, Lysl 10, Ser127, and Thr181 is mutated into Ala. According to the results of the molecular dynamics simulation, Thr181 is considered to be one of the key residues that helps NADH to approach the isoalloxazine in FAD. Finally, we have constructed very simplified model systems and carried out density functional theory calculations using B3LYP/LANL2DZ//ROHF(or RHF)/LANL2DZ level of theory in order to elucidate a realistic and feasible mechanism of the hydride-ion transfer from NADH to FAD affected by HEME(Fe3+) as an electron acceptor. Our calculated results suggest that the electron and/or hydrideion transfer reaction from NADH to FAD can be accelerated in the presence of HEME(Fe3+).