Despite intense interest in designing positive allosteric modulators (PAMs) as selective drugs of the adenosine A1 receptor (A1AR), structural binding modes of the receptor PAMs remain unknown. Using the first X-ray structure of the A1AR, we have performed all-atom simulations using a robust Gaussian accelerated molecular dynamics (GaMD) technique to determine binding modes of the A1AR 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 A1AR. 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 A1AR alone. In summary, the GaMD simulations elucidated structural binding modes of the PAMs and provided important insights into allostery in the A1AR, which will greatly facilitate the receptor structure-based drug design.