Reviewer 3 comments

[labarba]

There are a number of concerns regarding this work.

• [x] 1) There is little innovation in this work. The formulation is well-known. There have been numerous papers on this approach.

• [x] 2) It is not clear what is the convergence of the methods. Note that high-order boundary element based Poisson-Boltzmann (PB) solvers have been developed by many authors, including Krasny and Geng.

• [x] 3) The convergence of the software was not carefully validated with problems of known solutions. This validation is required to publish in a normal journal, not to mention a Nature one.

• [x] 4) There is no systematical comparison with other methods in the literature for either simple geometry or complex geometry.

• [x] 5) How does the package do for geometric singularities, which are common in protein surfaces?

• [x] 6) I have also a concern of the perspective of the field of this paper. Many important progresses were not mentioned, for example: Improvements to the APBS biomolecular solvation software suite, E Jurrus, D Engel, K Star, K Monson, J Brandi, LE Felberg, DH Brookes, Protein Science 27 (1), 112-128.

• [x] 7) Many online PBE solvers available. A user just needs to give a PBD ID to get result. However, the speed of solver is important. I am not sure if Python codes are the best of speed.

[labarba]

item 1

There is little innovation in this work. The formulation is well-known. There have been numerous papers on this approach.

We do not claim novelty in regards to the model formulation or numerical method. We claim to be providing the field with a new workflow that combines high researcher productivity (via the ability to run the solvers from Jupyter), with high-performance execution (via the binding with a performant fast multipole library). The paper presents evidence for the correctness (solution verification via grid-convergence analysis), and results of performance with the fast algorithm (FMM) that provides O(N) computational complexity. The software infrastructure is state-of-the-art in terms of performance, while providing high researcher productivity via the ability to run computational experiments interactively from the Jupyter environment. We assert that the combination of high performance and high researcher productivity is a paradigm for computational science that is emerging and will take hold. Moreover, our work is presented with utmost care for reproducible research, which also is a model to follow.

[labarba]

item 2

It is not clear what is the convergence of the methods. Note that high-order boundary element based Poisson-Boltzmann (PB) solvers have been developed by many authors, including Krasny and Geng.

We observe second order convergence of the energy (evaluated at the interior point charges) with respect to the grid spacing h (which appears as slope -1 with respect to the number of elements, N) in the computational grid convergence study. This is the same convergence rate seen by Krasny and Geng in "A Treecode-Accelerated Boundary Integral Poisson-Boltzmann Solver for Electrostatics of Solvated Biomolecules" (2013) J. Comput. Phys. using piecewise constant elements, and slightly better than the N^(1/2) convergence reported by Lu and McCammon in "Improved Boundary Element Methods for Poisson-Boltzmann Electrostatic Potential and Force Calculations" (2007) J. Chem. Theory Comput. using a node-patch approach that yields linear boundary elements.

[labarba]

item 3

The convergence of the software was not carefully validated with problems of known solutions. This validation is required to publish in a normal journal, not to mention a Nature one.

We provide detailed convergence studies. The first uses a spherical molecule, which has an analytical solution (Kirkwood's well-known solution). The second uses a small biomolecule (5PTI), where we study convergence using Richardson extrapolation. However, we can perform additional convergence analyses using other biomolecules, for extra evidence, as requested.

We included 8 more molecules in our convergence study in this revision, see item 4 below.

---

From the editor, after clarification with the reviewer:

Concern 3) Grid refinement validation should be carried out for a large set of realistic biomolecules in terms of various evaluation metrics.

[labarba]

item 4

There is no systematical comparison with other methods in the literature for either simple geometry or complex geometry.

Since our paper is not introducing a new method that should be systematically compared with others to be shown credible, we interpret this comment as referring to the software comparison. We are using an established mathematical model (linearized Poisson-Boltzman representation with implicit solvent) and established numerical methods (Galerkin boundary element method, and fast multipole methods). The paper does includes formal solution verification via grid-convergence analysis with analytical solutions and real-world structures using Richardson extrapolation. However, we had not included comparison with other trusted software, which has now been added.

In this revision, we supplemented the grid-convergence study with 8 more proteins, and compared our results with both APBS and MIBPB. We performed a grid-convergence study for each molecule on APBS as well, to obtain an approximate exact solution using Richardson extrapolation. We compared our extrapolated solutions to those from APBS, and these results can be found in the result section. In addition, we added the comparison with MIBPB only using finer meshes in the Appendix.

New results: grid-convergence analysis with 9 different molecules, with our software and two others (APBS and MIBPB), obtaining the estimated exact value with Richardson extrapolation, reported on table 2 and the appendix.

---

From the editor:

Concern 4) Examples of other methods:

Li, A. and Gao, K., 2016. Accurate estimation of electrostatic binding energy with Poisson-Boltzmann equation solver DelPhi program. Journal of Theoretical and Computational Chemistry, 15(08), p.1650071. Nguyen, D.D., Wang, B. and Wei, G.W., 2017. Accurate, robust, and reliable calculations of Poisson–Boltzmann binding energies. Journal of computational chemistry, 38(13), pp.941-948.

[labarba]

item 5

How does the package do for geometric singularities, which are common in protein surfaces?

We suspect the reviewer is addressing the effect of using a different representation of the molecular surface: solvent-excluded, or solvent-accessible surface. The latter can produce singularities (cusps), leading to slight changes in the mathematical formulation. Bempp uses a Galerkin formulation, which is well suited for any globally Lipschitz, piecewise smooth surface; we always use solvent-excluded surfaces, which satisfy this condition.

[labarba]

item 6

I have also a concern of the perspective of the field of this paper. Many important progresses were not mentioned, for example: Improvements to the APBS biomolecular solvation software suite, E Jurrus, D Engel, K Star, K Monson, J Brandi, LE Felberg, DH Brookes, Protein Science 27 (1), 112-128.

We have added this citation.

• fixed in 115bc2c57dd09a4b2791988a3af9bf32c4bfc4e4

[labarba]

item 7

Many online PBE solvers available. A user just needs to give a PBD ID to get result. However, the speed of solver is important. I am not sure if Python codes are the best of speed.

The Python environment introduces an overhead that is noticeable with the small-scale problems, but at the larger scale, this cost is minor, while the benefit and convenience of Python provides the researcher with ample opportunity to be more productive.

The intense computations are being performed at high performance with Exafmm (in C/C++) and low-level OpenCL optimized code in the case of Bempp. Take 1F6W, a large protein with over 8k atoms, for example, our code took 10mins to solve on a workstation (with two 20-core Intel Xeon Gold 6148 CPUs) while MIBPB took 15mins for a medium-sized mesh. We reported that we computed the solvation energy of a Zika virus with 10mil elements within 1.5 hrs on a single node, with 4 digits of accuracy in FMM. As a rough comparison, DASHMM-accelerated AFMPB solved a 20mil-element Dengue virus using 1.9 hrs on a single node, with 3 digits of accuracy in FMM. (This is calculated base on the time-to-solution being 30s on 512 nodes with a 45% parallel efficiency). We decided against adding these comparisons in the paper because performance comparisons often bring about concerns about the "fairness" of the comparisons. These reference timings are all obtained with different hardware, and thus easily criticized as "unfair."

Ref: Bo Zhang, Jackson DeBuhr, Drake Niedzielski, Silvio Mayolo, Benzhuo Lu & Thomas Sterling. (2020). DASHMM Accelerated Adaptive Fast Multipole Poisson-Boltzmann Solver on Distributed Memory Architecture. Communications in Computational Physics. 25 (4). 1235-1258. doi:10.4208/cicp.OA-2018-0098