Note from referee 1 This paper presents many results on Monte Carlo methods for solving systems of PDEs used to model biological electrostatics problems. These involve various implicit solvent models, the Poisson, linearized Poisson-Boltzmann (PB) and the fully nonlinear PB equations. The paper reviews many of the latest results, but has two major contributions to the literature. The first is the use of more sophisticated computational geometry algorithms in the linear case. Here, the fundamental computation is geometric, where a query point becomes the center of the largest sphere that can be constructed outside the domain of interest. The authors use Edelsbrunner's power diagram as the key to this computation, as power diagrams are implemented in the Computational Geometry Algorithms Library (CGAL) library. The second contribution is the use of branching processes to solve the nonlinear PB equation. Here, the paper presents only preliminary results for geometries where Kirkwood published exact solutions: one and two spheres. This is a very nicely written paper. The organization is very nice, the English is impeccable, and the results convincing. While the work is a collaboration between mathematicians and computer scientists, the paper itself is more heavily weighted towards the mathematics. Thus, the main criticisms of the paper really deal with computer science issues that are not adequately discussed. In addition, the paper assembles many results together into a nice story, as well as providing new results. As a paper for applied mathematicians, it is quite nice, but for computer scientists and biochemists, these results should really be extended and published in a manner that is more suitable to these other professionals. M. Mascagni (Florida State University) Note from referee 2 This articles proposes a new method to solve a non-linear Poisson-Bolzmann PDE, here in the context of computing the eletrostatic potential close to a biomolecule. This is done through a probabilistic representation in which the positions of particles moving randomly in the medium (and sometimes killed) are averaged. While in the free space, the particles are moved using the random walk on spheres methods for computational efficiency, there are several issues this articles deals with, and this article completes and expand an previous article by some of the authors (ESAIM M2AN, vol 45, 2010). The first one is to describe and implement the proper behavior of the particle at the interface between the biomolecule and the external medium. Several replacement schemes are proposed and studied. The second one is the non-linearity, for which a new probabilistic representation is given through a branching mechanism, together with the corresponding algorithm. In addition to new theoretical considerations, simulation are performed both in the linearized and the non-linear case, with a discussion on their efficiency of the method and the computational geometry problems one faces using the replacement scheme. The non-linear part could be improved through variance reduction techniques in futures works.