Phenylalanine hydroxylase (PAH) is an enzyme that catalyzes the conversion of phenylalanine to tyrosine, and functions in humans to control free phenylalanine (Phe), an essential amino acid that is neurotoxic at elevated levels. Mutations in PAH can result in phenylketonuria, which is the most common inborn error of amino acid metabolism. The transition from resting-state PAH to activated PAH requires formation of a new intersubunit interface that can be stabilized by Phe binding in an allosteric manner; formation of the new protein:protein interface is coupled to exposure of the active site, thus activating the enzyme. A recent crystal structure of this interface bound to Phe enables new insight into phenylalanine binding mechanism from molecular dynamics (MD) simulations. We performed massively parallel explicit-solvent simulations on the [email protected] distributed computing platform to elucidate pathways and rates of binding of Phe to the interface, which is comprised of a dimer of ACT domains. Time-lagged independent component analysis (tICA) of binding trajectories suggest a conformational selection mechanism. Markov state models (MSMs) constructed from the trajectory data reveal a key loop motion which acts as a “gatekeeper” of allosteric ligand binding. Binding rates estimated by different methods (MSMs, Transition Path Theory and Bayesian inference) agree well with each other. These results warrant further MD studies of the conformational dynamics involved in regulation of PAH activity, the effect of disease-associated mutations, and suggest future directions in simulation-based drug discovery.