In this talk I will give an overview of the recent development of boundary integral methods for solving Poisson-Boltzmann equations (PBE) for the electrostatic potential in a solvation system of biomolecules. Electrostatic forces are crucial in determining the structure and dynamics of biomolecules and their interactions with solvents. As a mean-field approximation, the PBE has proved to be a very useful model of such electrostatics. However, numerically solving the PBE in a very efficient and accurate way is challenging. I will first describe the method and present some of my computational results. I will then compare different numerical methods for solving the PBE with complicated geometry. In particular, I will present the newly developed fast multipole method for the boundary integral discretization of the PBE. I will finally mentions some open questions. The results presented in talk are mainly from the work done in McCammon's group at UCSD.