The protonation state of the catalytic aspartates of HIV-1 protease (HIVPR) is atypical and as a result is the subject of much debate. Modeling of the correct protonation state of the aspartates is vital in computational drug design. Using pH replica-exchange molecular dynamics, we simulated the apo and bound forms of HIV-1 protease with 12 different protease inhibitors to investigate the pK(a) of not only the catalytic dyad but also the other titrating residues in HIVPR. The results obtained from these simulations are compared to experiment where possible. This study provides evidence that the catalytic aspartates are primarily in a monoprotonated state for both the apo and bound forms of HIVPR in the pH range where generally most experiments and computational simulations occur.