A new method was proposed by Mongan et al for constant pH molecular dynamics simulation and was implemented in AMBER 8 package. Protonation states are modeled with different charge sets, and titrating residues are sampled from a Boltzmann distribution of protonation states generated by Monte Carlo sampling based on Generalized Born (GB) derived energies. However, when this approach was applied to a bio-toxin, Botulinum Neurotoxin Type A (BoNT/A) at pH 4.4, 4.7, 5.0, 6.8 and 7.2, the pKa predictions yielded by the method were inconsistent with the experimental values. The systems being simulated were divergent. Furthermore, the system behaviors in a very weak acidic solution (pH6.8) and in a very weak basic solution (pH7.2) were significantly different from the neutral case (pH7.0). Hence, we speculate this method may require further study for modeling large biomolecule. Key words: Constant pH molecular dynamics; Botulinum Neurotoxin; Generalized Born Method