Round 2, Reviewer 6 comments

[labarba]

The authors present a python PB solver based on a fast multipole method library and a Galerkin boundary element method package. The method integrates a python interface with optimized computational kernels, and can be run interactively via Jupyter notebooks for ease of prototyping. Besides describing the mathematical frameworks, they further present results demonstrating the capability of the software, confirming code correctness, and assessing performance with between 8,000 and 2 million elements. To illustrate the interactive computing platform, they studied the conditioning of two variants of the boundary integral formulation with just a few lines of code. Mesh-refinement studies confirm convergence as 1/N, for N boundary elements, and a comparison with results from the APBS code using various proteins shows agreement. Performance results include timings, breakdowns, and computational complexity. To highlight the high performance of the method, they illustrate the compution for the solvation free energy of a Zika virus structure, represented by 1.6 million atoms and 10 million boundary elements.

A. On the intended audience, I think the intentions of the authors to provide a new platform for “high productivity and high performance” are reasonable for a journal interested in publishing articles reporting new resources.

B. On the intended functions, the authors are proposing a python-notebook-based workflow for “high productivity and high performance”. The innovation is not in the underlining algorithms that are already published elsewhere, but in the development of a workflow, which is a reasonable effort worth reporting.

C. On the agreement with existing software packages. There are some numerical disagreements between this implementation and those in the literature. It is possible that these are due to the different treatments of the salt effect. The boundary-element methods model the salt term all the way to the molecular surface in the solvent region, while finite-difference methods usually model the salt term a single ion layer (roughly 2 Angstrom) away from the molecular surface. If the salt concentration is not very low and the solute net charges are high enough, the differences in their solvation energies can be large. As the detail on the salt term is not given for the molecular calculations, it is unclear whether this could account for the observed differences. Another possibility lies in the different numerical surfaces used in different methods. Unless the exact numerical surfaces are used, it is unlikely to have very high numerical consistency even if the same atomic parameters are
used. Overall, PB solvation energies are very sensitive to where the molecular surface is located.

C. Other comments.

1. The major issue of this reviewer is that the manuscript does not appear to be one reporting/showcasing computing resources. The presentation of the manuscript is not consistent with the claimed goals in providing a python workflow. Instead, it is more like a regular manuscript reporting a new PB solver even after multiple rounds of revisions. A big chunk of the manuscript is on validation of the algorithm/implementation. These materials dilute the main theme of “high productivity and high performance” and are better be moved to appendices. It is recommended that the authors refer to those resource articles in Bioinformatics and Nucleic Acid Research on how to showcase their software functionalities.

The bottom line is that the authors need to talk more about the workflow functionalities. You do not just say your software has certain functionalities without showing. Specifically,

• The authors need to show the ease of using the python Jupyter notebooks by using figures/screenshots how prototyping is done with actual lines of code in the text for playing with different methods as the authors mentioned. Of course, more samples can be posted online.

• The authors need to show how to use their Jupyter notebooks to compute energies for multiple structures from MD as claimed as they agree a single solvation energy does not make much sense for a virus particle. Note that they need to conduct statistics as well. If possible, visualization of the data would be great too.

2. To highlight the scalability of their PB solver, a single virus particle is not enough. Instead time scaling for a range of smaller to larger proteins/nucleic acids, molecular complexes, maybe all the way to virus will be very helpful. If they can compare with APBS with its default recommended setting, it will be great.

3. Section “Mesh refinement study using a spherical molecule” is better labeled as comparison with analytical solutions. Otherwise, it is more appropriate to be combined with the next section on “Mesh refinement study using 5PTI”.

4. Section “Comparison with trusted community software” is partially overlapped with the Mesh Refinement sections as both sections addressed convergence. The authors can just show the extrapolated energy values from both methods to save space and stay focused on the main point. If they want to compare with MIBPB, all three values can be put into the same table, saving the space for more important topics.

5. In section “Comparison with trusted community software” it is difficult to assess the quality of the numerical calculations as the authors do not provide actual grid spacing used as shown in their BEM method. Please talk about the finite-difference resolution in term of grid spacing (Angstrom) or how many grids per Angstrom. Number of grid points are not very useful.

6. It is also surprising to see that the APBS convergence scaling is far from the second order, which should be the case for most widely used finite-difference PB solvers. Not sure this is because the grid spacing is still too large or some other reasons. They may want to refer to the original APBS papers to see why this is the case.

7. There should be a section on how authors prepared all molecules. Only the section on 5PTI has such discussion, but it lacks the mentioning of dielectric constants (solute/solvent regions) and ionic strength/kappa for the salt term.

8. There are no citations for pdb2pqr and the charm force field. There are no citations for Bempp and Exafmm packages.

[tingyu66]

2. To highlight the scalability of their PB solver, a single virus particle is not enough. Instead time scaling for a range of smaller to larger proteins/nucleic acids, molecular complexes, maybe all the way to virus will be very helpful. If they can compare with APBS with its default recommended setting, it will be great.

In this revision, we added the performance comparison with APBS using 1RCX, a medium-sized protein. We want to point out that making fair performance comparison is not trivial. We cannot simply run both software with default settings and report the timings since the results might be at different accuracies. Moreover, we cannot force them to achieve the same accuracy by simply matching the computed solvation energies. It is because the two codes carry different modeling errors and will converge to different solutions. Therefore, the only fair way to compare the performance of each code is 1) to perform grid convergence study to obtain the extrapolated solution of each code when the discretization error diminishes to zero, 2) to compute the error based on this extrapolated solution, respectively 3) identify mesh spacings that match the accuracy of the two solvers. Only at this matched accuracy is it appropriate to compare runtimes. Following this procedure, we measured the performance of Bempp using 4 different mesh densities and the performance of APBS using 3 grid spacings respectively, to compute the solvation energy of 1RCX.
We would like to be able to do such a comparison for a spectrum of proteins of different sizes. However, using APBS with structures larger than 1RCX on our cluster node (with 192GB memory) resulted in a runtime error (segmentation fault 11). For problems with smaller structures, most PB solvers can finish within minutes. We think showing these cases (e.g., 60 vs 90 seconds) would not add much to our paper. Instead, the ability of Bempp to handle larger structures with decent accuracy is what we want to demonstrate here.

[labarba]

Section “Mesh refinement study using a spherical molecule”

We changed the heading to say "analytical solution."

[labarba]

In point 1, the reviewer also comments:

The authors need to show how to use their Jupyter notebooks to compute energies for multiple structures from MD as claimed as they agree a single solvation energy does not make much sense for a virus particle. Note that they need to conduct statistics as well. If possible, visualization of the data would be great too.

The reviewer is suggesting that we could demonstrate a workflow where we compute salvation energy for a virus, using multiple structures (conformations) of the viral particle. Running MD simulations to obtain new conformations of a virus particle requires leadership computing resources, and moreover we are not MD experts, so we could only feasibly do what the reviewer suggests if publicly available structures existed for the Zika virus. We could not find any other structures in public databases. Since we do not have an ensemble of structures, we cannot do statistics. Note that running each case for a virus-scale structure takes in the order of 3 hours or more (compute time) on our university cluster. The reviewer's idea is interesting, but it would take enormous computing resources to achieve, even if we had the structures.

This reviewer's comment seems to be in response to our reply to Reviewer #1 in the first round. We don't address this in the paper, that is, we don't make promises to the reader that we're not fulfilling. Of course it is possible to run many structures in a Jupyter notebook in some research scenario: the point we're making is that you can call our solver from a notebook.

[labarba]

Point 4:

Section “Comparison with trusted community software”[…] stay focused on the main point ...

Note that the Results sub-section on "Comparison with trusted community software" was added aft the first round of reviews, in response to Reviewer #3's comments, who asked for comparison with other methods in the literature. We agree with Reviewer #6 that it is a side track, so we did the following: moved the results with APBS and our software on many molecules to an Appendix (this could be Supplementary Material), and removed the results with MIBPB. In the latter case, the fact is that iterative grid refinement does not show monotonically changing salvation energy (in other words, convergence towards a value). These results do not add anything to our main point.

More details:
The reason for removing the MIBPB results in this revision is that we did not observe convergence of MIBPB using a few proteins. Taking 1A63 for example, we generated the pqr file using pdb2pqr and the charmm force field. We then ran MIBPB on its online server using three grid spacing: 1.0, 0.75, 0.5 Angstrom. Ionic strength, relative permittivities are the same as what we reported in the manuscript, and the grid extension is at its default value of 4. The three computed solvation energies are: -586.50, -585.52, -587.43 kcal/mol. The trend is not monotonic and thus does not show convergence. We also downloaded the compiled version to test and ended up with the same results. We have tried our best to compare our code with MIBPB, however, at this point, we feel that it is no longer our jobs to justify the MIBPB results.

[labarba]

5. Please talk about the finite-difference resolution in term of grid spacing (Angstrom) or how many grids per Angstrom. Number of grid points are not very useful.

We moved these results to an Appendix (or Supplementary Information) now. Nevertheless, the reviewer's criticism stands. We should explain that in conducting iterative grid refinement, it is particularly hard to achieve a constant refinement ratio, given that APBS uses a multigrid method. In the results with multiple molecules, we tackled this by choosing a constant number of points in each dimension, which resulted in different grids spacings dx, dy, dz in the rectangular computational domain chosen by the solver. This is just one example of the hidden complexities of performing grid-convergence analysis.

As requested, we add the grid spacing to the table (now in Supplementary Information).

[labarba]

It is also surprising to see that the APBS convergence scaling is far from the second order, which should be the case for most widely used finite-difference PB solvers. Not sure this is because the grid spacing is still too large or some other reasons. They may want to refer to the original APBS papers to see why this is the case.

The order of convergence in APBS is likely affected by their approximations at the molecular boundary. The reviewer says that second order should be expected, but we could not find in the APBS papers that they claim or observe second order. We do see in the paper by Geng and Krasny, J. Comput. Phys. 2013 that they also observed a lower than second order. On their Table 4, you can see that the cases with APBS (last three rows) show that when you halve the grid spacing, the error is also decreasing by 1/2. This means an observed first-order of convergence.

[tingyu66]

7. There should be a section on how authors prepared all molecules. Only the section on 5PTI has such discussion, but it lacks the mentioning of dielectric constants (solute/solvent regions) and ionic strength/kappa for the salt term.

We mentioned the dielectric constants and kappa at the beginning of the Results section, and they apply to all test cases. We also mentioned that we parameterized 5PTI using charmm and 6CO8 using amber force field, respectively. The rest of cases only use spherical molecules where we don't need to apply a force field to generate the solvent-excluded surface.

[tingyu66]

8. There are no citations for pdb2pqr and the charm force field. There are no citations for Bempp and Exafmm packages.

We added citations to charmm and amber force field in commit c820722. We did cite Bempp and Exafmm packages already in the method section, and cite pdb2pqr in the results section.