Vibrational relaxation of CO bound to myoglobin (MbCO) following photoexcitation is investigated using nonequilibrium molecular dynamics (MD) simulations. It is found that harmonic potential energy functions for bond vibrations are not suited to simultaneously and accurately describe vibrational de-excitation and the vibrational spectroscopy of the bound ligand. Only when anharmonic (e.g. Morse) potentials are introduced for both the C-O and the adjacent Fe-C-CO bonds to allow anharmonic coupling, rapid (tens of ps) relaxation of the vibrationally excited CO is possible. To capture both relaxation and vibrational spectroscopy, the parameters of the potential energy functions are fitted by an interactive, nonlinear least-squares procedure using averages over multiple MD trajectories. The sensitivity of cooling rate to the difference in vibrational frequency between coupled modes is demonstrated. Potential cooling mechanisms are suggested, based on the sensitivity of the CO relaxation rate to changes in the force field parameters of local degrees of freedom. Accounting for quantum correction leads to relaxation rates around 20 ps, in good agreement with experiment. Finally, the importance of electronic effects is explored by fitting a 2D potential energy surface to ab initio data to describe the strengthening and weakening of the CO bond as a function of Fe-C-CO bond length, and vice versa.