Interactions between the plasmons of noble metal nanoparticles and non-absorbing biomolecules forms the basis of the plasmonic sensors, which have received much attention. Studying these interactions can help to exploit the full potentials of plasmonic sensors in quantification and analysis of biomolecules. Here, a quasi-static continuum model is adopted for this purpose. We present a boundary-element method for computing the optical response of plasmonic particles to the molecular binding events by solving the Poisson equation. The model represents biomolecules with their molecular surfaces, thus accurately accounting for the influence of exact binding conformations as well as structural differences between different proteins on the response of plasmonic nanoparticles. The linear systems arising in the method are solved iteratively with Krylov generalized minimum residual algorithm, and the acceleration is achieved by applying precorrected-Fast Fourier Transformation technique. We apply the developed method to investigate interactions of biotinylated gold nanoparticles (nanosphere and nanorod) with four different types of biotin-binding proteins. The interactions are studied at both ensemble and single-molecule level. Computational results demonstrate the ability of presented model for analyzing realistic nanoparticle-biomolecule configurations. The method can provide comprehensive study for wide variety of applications, including protein structures, monitoring structural and conformational transitions, and quantification of protein concentrations. In addition, it is suitable for design and optimization of the nano-plasmonic sensors.