A quantitative study of the interplay between the structure, energy, and dynamics of heme models is presented. Our calculations are based on density functional theory (DFT) combined with molecular dynamics (MD), within the Car-Parrinello scheme. The systems analyzed are the five-coordinated FeP-AB (FeP = iron-porphyrin, AB = CO, NO, O-2) and Fe(Heme)-O-2, the six-coordinated FeP(Im)-AB (Im = imidazole), and the picket fence-oxygen complex, Fe(TpivP)(2me-Im)-O-2. Starting from the optimized structures of these systems, we analyze the trends in their ligand binding properties. Our calculations show the peculiar trans repulsive effect of the NO ligand when binding to iron-porphyrin (it elongates and weakens the trans axial Ligand bond), in contrast with the synergistic effect of CO and O-2. Energy variations associated to the rotation of the O-2 molecule around the Fe-O bond in Fe(TpivP)(2me-Im)-O-2 are analyzed, in comparison with that of a simpler FeP(Lm)-O-2 model. A small energy barrier (<2 kcal/mol) is found in both cases. Room temperature molecular dynamics simulations reveal how the O-2 Ligand overcomes this barrier at room temperature: It undergoes large-amplitude oscillations within one porphyrin quadrant, jumping to another one in the picosecond timescale. In contrast, the CO dynamics is characterized by small, albeit extremely complex, displacements around its equilibrium position. Small fluctuations of the FeCO tilt (<delta>) and bend (theta) angles (delta less than or equal to 8 degrees, theta less than or equal to 13 degrees) are observed. (C) 2000 John Wiley & Sons, Inc.