AbstractDespite intense interest in designing positive allosteric modulators (PAMs) as selective drugs of the adenosine A₁ receptor (A₁AR), structural binding modes of the receptor PAMs remain unknown. Using the first X-ray structure of the A₁AR, we have performed all-atom simulations using a robust Gaussian accelerated molecular dynamics (GaMD) technique to determine binding modes of the A₁AR allosteric drug leads. Two prototypical PAMs, PD81723 and VCP171, were selected. Each PAM was initially placed at least 20 Å away from the receptor. Extensive GaMD simulations using the AMBER and NAMD simulation packages at different acceleration levels captured spontaneous binding of PAMs to the A₁AR. The simulations allowed us to identify low-energy binding modes of the PAMs at an allosteric site formed by the receptor extracellular loop 2 (ECL2), which are highly consistent with mutagenesis experimental data. Furthermore, the PAMs stabilized agonist binding in the receptor. In the absence of PAMs at the ECL2 allosteric site, the agonist sampled a significantly larger conformational space and even dissociated from the A₁AR alone. In summary, the GaMD simulations elucidated structural binding modes of the PAMs and provided important insights into allostery in the A₁AR, which will greatly facilitate the receptor structure-based drug design.