G protein-coupled receptors (GPCRs) continue to hold leading positions as drug targets. However, many GPCR drug candidates fail in clinical trials because of limited in vivo efficacy. While binding affinity, i.e. the strength of association of a drug to its receptor, has traditionally been viewed as an appropriate surrogate for in vivo efficacy, retrospective analyses of marketed drugs suggest that kinetic parameters, such as the rates at which the drug associates with or dissociates from the target, may play a role that is as important as, or even more important than, binding affinity in determining in vivo efficacy. Thus, predictions of both kinetic and thermodynamic parameters of ligand binding to GPCRs are highly desirable because they may inform the rational discovery of improved therapeutics.The main challenge in studying molecular recognition by a GPCR using molecular dynamics (MD) is that ligand binding and unbinding are rare events on microscopic time-scales, and as such, they are difficult to observe using unbiased simulations. In this work, we present a general strategy that employs biased MD simulations to build Markov State Models (MSMs) of the binding of small molecules to a prototypic family A GPCR. Using Perron cluster analysis and transition path theory we were able to identify the kinetic basins of the binding process, and to characterize both the metastable and the transition states between bound and unbound conformations. By investigating the role of hydrophobic interactions, dewetting, and conformational changes in the binding pocket, we were able to characterize the microscopic determinants that influence association and dissociation rates of the ligand. This information has a direct applicability in rational drug design approaches.