Understanding the binding mode of agonists to adrenergic receptors is crucial to enabling improved rational design of new therapeutic agents. However, so far the high conformational flexibility of G protein-coupled receptors has been an obstacle to obtaining structural information on agonist binding at atomic resolution. In this study, we report microsecond classical molecular dynamics simulations of b1 and b2 adrenergic receptors bound to the full agonist isoprenaline and in their unliganded form. These simulations show a novel agonist binding mode that differs from the one found for antagonists in the crystal structures and from the docking poses reported by in silico docking studies performed on rigid receptors. Internal water molecules contribute to the stabilization of novel interactions between ligand and receptor, both at the interface of helices V and VI with the catechol group of isoprenaline as well as at the interface of helices III and VII with the ethanolamine moiety of the ligand. Despite the fact that the characteristic N-C-C-OH motif is identical in the co-crystallized ligands and in the full agonist isoprenaline, the interaction network between this group and the anchor site formed by Asp(3.32) and Asn(7.39) is substantially different between agonists and inverse agonists/antagonists due to two water molecules that enter the cavity and contribute to the stabilization of a novel network of interactions. These new binding poses, together with observed conformational changes in the extracellular loops, suggest possible determinants of receptor specificity.