(PRE) Reviewer 1 comments

The manuscript presents a software package for boundary element method (BEM) calculations that has been modified to calculate optical extinction spectra for models of biological molecules in close proximity to plasmonic metal nanoparticles. These systems require advances in modelling to address questions in the plasmonic response changes used in widespread optical sensing applications. The software uses treecode to reduce computational costs. The changes to the code from its previous version include edits to account for the required complex number input and to enable the calculations as a function of wavelength and to calculate optical exctinction. The manuscript presents several examples illustrating convergence of the code for problems with analytical solutions and then presents an application to a system containing a small silver particle with two protein molecules in close proximity.

The manuscript is generally well written and presents the changes to the code and theoretical background with sufficient clarity. The challenge of representing complex geometries in BEM calculations with sufficient mesh density is important for improving modelling capabilities for realistic representations of plasmonic particles and for understanding subtleties of their optical effects.

[x] 0) The manuscript may therefore merit publication in Physical Review E, although it can be improved significantly to assist readers in understanding the broader context of BEM approaches. The explanation of how the mesh is used to model biological proteins also requires revision to clarify the physical constraints or interpretation of the resulting calculations.

I have identified several key questions and comments below:

[x] 1) The manuscript would be significantly improved if the BEM formulation could be compared with alternative implementations in more detail. This would assist with comparisons between implementations which is ultimately linked to aspirations for checks on reproducibility and validation of the calculations beyond analytical models. The question of validation between codes is beyond the scope of the manuscript, but pointing readers to key similarities or key differences in BEM implementations would provide a more complete introduction and background framework on the use of BEM in solving quasi-static as well as full Maxwell’s equations.
[x] 2) The relaxation of some parameters for calculations shown in Figure 6 is not inherently a problem, but I would have thought that a figure showing the wavelength dependence for the silver sphere with decreasing meshes and for directly comparable sets of treecode parameters used in Figure 5 would provide another important consistent comparison.
[x] 3) For BSA and other proteins, is it physical to represent them as an unconstrained mesh? What I mean is, there are presumably regions that are (chemically) motivated to be charged in particular ways which would seem to possibly contribute dipoles at relevant length scales? It is not clear to me in this manuscript how the dielectric function value and its local or total dipole moment is represented on the mesh used for the protein molecules. This may be discussed in other work, but should be explained within this manuscript as it is central to the calculations reported here.
[x] 4) Is a classical treatment appropriate for distances of 1 nm (or less than 1 nm)? How do results compare with 5 nm away? And how would the calculations change for a gold particle rather than a silver one?
[x] 5) Why have the authors elected to show calculations for two BSA proteins not just one? Are the positions of the two particles opposite each other selected because this relative positioning enhances the reported effects?
[x] 6) What are the relevant distances and concentrations for experimental comparisons? Can this code handle those? (THE REPLY NEEDS SOME MORE WORK)
[x] 7) Why do the authors report only the extinction cross section Cext? Are there easy ways to calculate other optical properties from the solved system (e.g. Cabs, Csca, near fields, scattering matrix)? Perhaps a comment on the feasibility or computational obstacles to these alternative outputs could be offered.

A few minor points:

[x] 8) In the introduction of extinction, it is unclear if ‘a.k.a.’ is referring to both scattering and absorption or just to the second part of sentence (scattering).
[x] 9) Is Gamma defined at its first appearance on p. 2?

Comment 0

The manuscript may therefore merit publication in Physical Review E, although it can be improved significantly to assist readers in understanding the broader context of BEM approaches. The explanation of how the mesh is used to model biological proteins also requires revision to clarify the physical constraints or interpretation of the resulting calculations.

Ideas

"Help reader in understanding BEM in a broader context": we should explain why BEM is an appropriate tool for this problem.
"The explanation of how the mesh is used...": we should explain what the solvent excluded surface is, and the pipeline from PDB to the mesh.
Add a subsection under Methodology, in the end, explaining the pipeline to obtain the mesh and explaining the solvent excluded surface. Tentative name: Computational workflow
Reference Fig. 3 to explain the point charges
Move explanation about Nanoshaper in results (right before section III.B.1.) to the new section.

Reply

The way biomolecules are modeled with the boundary element method is not detailed the original manuscript, which led to this comment by the reviewer. This is not the main effort of the paper, and has been done previously with our code in reference [26] (Cooper, Bardhan, Barba (2014)), nonetheless, we should clarify for completion. We have added a subsection under Methods called 'Protein mesh preparation', where we explain how we start from a molecular structure from the Protein Data Bank, and represent it as a surface mesh, clarifying the physical meaning of the molecular surface, and how the biomolecule is parameterized in terms of van der Waals parameters and charges.

Modifications

We added a subsection called 'Protein mesh preparation' in commit fb1785c and ab38e80

Comment 1

1) The manuscript would be significantly improved if the BEM formulation could be compared with alternative implementations in more detail. This would assist with comparisons between implementations which is ultimately linked to aspirations for checks on reproducibility and validation of the calculations beyond analytical models. The question of validation between codes is beyond the scope of the manuscript, but pointing readers to key similarities or key differences in BEM implementations would provide a more complete introduction and background framework on the use of BEM in solving quasi-static as well as full Maxwell’s equations.

Ideas

Expand discussion on comparison with BEM++ and MNPBEM, and place it earlier in the text.
Move the discussion on BEM++ and MNPBEM to the introduction. Place it before the sentence 'The software is ahred under BSD...', then we could add a paragraph break.

Reply

There are two research software codes that could be used in this setting: BEM++ and MNPBEM. They can solve either full Maxwell's equations or the quasi-static electrostatic approximation. However, they haven't been used in the context of biomolecules, and are not able to model with all the way to hundred-thousand boundary elements, which is required to represent molecular surfaces accurately.

Part of this discussion was present in the conclusion of the original manuscript, however, we've moved it to the introduction, and expanded it. This way, the reader will have a clearer idea of differences with similar available software.

Modifications

We have moved the discussion on similar software (BEM++ and MNPBEM) to the introduction and expanded it in commit d44c84e

Comment 2

2) The relaxation of some parameters for calculations shown in Figure 6 is not inherently a problem, but I would have thought that a figure showing the wavelength dependence for the silver sphere with decreasing meshes and for directly comparable sets of treecode parameters used in Figure 5 would provide another important consistent comparison.

Ideas

Report timings of relaxed and non-relaxed parameters to show the importance of relaxing the parameters
Perhaps we're not being clear on the difference of a mesh convergence study and its difference with the Treecode parameters. We should explain more clearly that the errors introduced by Treecode/integration parameters need to be below the meshing error.
We could mention the timings for calculations in Figs. 5 and 6

Reply

[LB] Figure 6 shows a visually undetectable error compared with the analytical solution. Repeating these calculations with a finer mesh or tighter parameters for the treecode will decrease the errors, but the figure will look the same. Currently, errors are everywhere smaller than 1%. These finer mesh and parameters, however, increase the time to solution by 12X.

Modifications

added runtime advantage of relaxed parameters in 91e73b5

Comment 3

3) For BSA and other proteins, is it physical to represent them as an unconstrained mesh? What I mean is, there are presumably regions that are (chemically) motivated to be charged in particular ways which would seem to possibly contribute dipoles at relevant length scales? It is not clear to me in this manuscript how the dielectric function value and its local or total dipole moment is represented on the mesh used for the protein molecules. This may be discussed in other work, but should be explained within this manuscript as it is central to the calculations reported here.

Ideas

This is related with Comment 0. It might help to add an explanation of the model with an isolated protein only, or explain eq 13, and how its RHS projects the charges to the surface mesh.
Mention after Eq.(17) that the last term accounts for the point charges in the protein.
We never introduce the variable q_k

Reply

[LB] The molecular charges that appear in the last term of Eq. (17) are considered explicitly, and their effect is not transferred to the mesh in the model. The effect of these charges on the boundary (mesh) is captured by this mathematical term. The dielectric function is a uniform complex value inside the protein, and appears as epsilon_3 in Eq. (17). Its value depends on the wavelength.

[CC] Perhaps this was not clear as we didn't introduce the physical meaning of the q_k variable in the original manuscript. This is fixed in the revised version.

Modifications

To explain how the BSA molecules is modeled, we added a subsection called 'Protein mesh preparation' in commit fb1785c and ab38e80. We also added an explanation on the meaning of q_k in commit f648cd1. The molecule contains point charges at the local partial charges from applying pdb2pqr with an Amber force field. This model is broadly used to represent biomolecules.

Comment 4

4) Is a classical treatment appropriate for distances of 1 nm (or less than 1 nm)? How do results compare with 5 nm away? And how would the calculations change for a gold particle rather than a silver one?

Ideas

Related with Comment 1 of reviewer 2.
At 5nm the system is barely capable of detecting (shift is very small).

Is a classical treatment appropriate for distances of 1 nm (or less than 1 nm)? Naty found this paper that is an MD simulation (coarse grained) with NP and an electric field https://pubs.rsc.org/en/content/articlepdf/2016/nr/c6nr02051h

"The electrostatic interactions were calculated using the particle mesh Ewald (PME) method with a real space cut off length of 1.2 nm and a fast Fourier-transform grid spacing of 0.24 nm"
Fig 5 shows distances between nanoparticle and membrane and they have values smaller than 1nm, even close to 0nm.

Chris found other papers where there is no electric field but the nanoparticles or nanotubes are smaller than 1nm:

https://pubs.acs.org/doi/abs/10.1021/acs.jpcc.8b02582 1.4 nm diameter but the particles get closer than that distance. Distances between particles is smaller.

Reply

Classical treatment is appropriate for this small distances, a point of comparison is this example of coarse-grained MD (classic limit) where metallic nanoparticles are at similar distances from a membrane under an electric field. In this paper (https://pubs.rsc.org/en/content/articlepdf/2016/nr/c6nr02051h)
- "The electrostatic interactions were calculated using the particle mesh Ewald (PME) method with a real space cut off length of 1.2 nm and a fast Fourier-transform grid spacing of 0.24 nm"
- Fig 5 shows distances between nanoparticle and membrane and they have values smaller than 1nm, even close to 0nm.
At 5 nm we won't see the effect of the proteins on the NP, the shift will be negligible. At 2 nm we already see a smaller shift, reported in section Results-Sensitivity calculations (Fig 12).
In order to answer fully this question we would need to perform multiple calculations using gold, which is out of the scope of this paper. However, we can say that the shift magnitude would be smaller due gold Figure of Merit (how a nanoparticle's sensing capabilities is characterized). The reason why we chose silver instead of gold it is because the FOM of silver is better Ref FarooqAraujo2018.

Modifications

no modifications needed.

Comment 5

5) Why have the authors elected to show calculations for two BSA proteins not just one? Are the positions of the two particles opposite each other selected because this relative positioning enhances the reported effects?

Ideas

One BSA molecule generates a very small effect.
The reason why we placed them in opposite positions was that this was the configuration where the shift was the greatest.
I think this is not a problem considering this is a proof-of-concept work
Reply
Simulations were performed using two BSA dimers, hence, effectively, there are four BSA molecules. There is experimental evidence that supports this choice (Teichroeb et al 2008, reference [41] in the paper, see Comment 6). Moreover, with less analytes the red shift would be smaller, and we would need to examine a finer range of wavelengths to see a shift.

Modifications

no modifications needed

Comment 6

6) What are the relevant distances and concentrations for experimental comparisons? Can this code handle those?

Ideas

We should add some justification of the parameters used: sizes, number of molecules, etc.

Look for size of sdAb in raphael paper for reference. There are sdAB used in Biosensor development that have a size of 4x2.5x3.5 nm (Saerens etal 2008) this can be used of an estimate in the the order of magnitud of the distances we are using. https://www.mdpi.com/1424-8220/8/8/4669/htm

try to make an estimate based on some volume to get count of BSA

Update Based on the concentrations reported on certain papers in micro moles/L I've tried to compute amount of BSA molecules in a volume (estimated layer of BSA) but this calculations tells us that the number of molecules is less than 1. After reading more papers, where they report amount of molecules per surface area we found that the number of molecules attached to some nanoparticles is between 1-6 depending on the case. This numbers are base on a monomer BSA (in our case we have a dimer.) Citations and more explanation in the possible reply section.

Dasgupta et al 2016 Case of silver NP of 30nm radius they report (based on a theoretical calculation) a maximum of 2.7 BSA per particles, while experimentally they see 1 or 2. In this case the size of BSA is a monomer, so it's about half the size of our protein.
Teichroeb et al 2008 This paper reports on a concentration per surface area of 3.2x10^12 molecules /cm^2 (they say is on end on absorption, not sure what that means. This numbers is 2x10^12 molecules /cm^2 for side-on absorption) which gives 0.032 molecules/ nm^2 (0.02 for side on) for molecules of diameter approx 15 nm. In our case diameter is 16 nm. Doing the math we an area of 201.06 nm^2 what gives approximately 6 (4 for value of side-on) molecules. This is assuming the dimenssions of BSA are (5.5 × 5.5 × 9 nm) In our case the the molecules are a bit bigger since it is a dimer.

Reply

Our simulations were performed under realistic conditions. Similar to our model with BSAs placed very close to the nanoparticle, there is experimental work where BSA is directly adsorbed on the metallic surface (without a layer of ligand in between), and the distance between them is effectively zero (Tang et al 2010 and Nghiem et al 2012). This has also been modeled computationally by other researchers (Phan et al 2013 and Dan and Hu 2014).

In our simulations, we looked at the resonance shift of two BSA dimers near a 16nm-diameter nanoparticle, which agrees with Teichroeb et al 2008. This work reports nanoparticles of 15nm diameter, with a concentration per surface area of 3.2x10^12 molecules /cm^2 (0.032 molecules/ nm^2). This gives approximately 4 to 6 molecules ( BSA-monomers) per nanoparticle.

Modifications

Add citation to [Teichroeb et al 2008]
Corrected citation of [Teichroeb et al 2008], since this paper only talks about concentration no distances. Commit number 0eb2637
Add more details on paper regarding Teichroeb results. Commit number aaa9e79, fixed wording and typos in eef190d3 and 14a5e06. Removed sentence regarding Nanoshaper in commit 789d931

Comment 7

7) Why do the authors report only the extinction cross section Cext? Are there easy ways to calculate other optical properties from the solved system (e.g. Cabs, Csca, near fields, scattering matrix)? Perhaps a comment on the feasibility or computational obstacles to these alternative outputs could be offered.

Ideas

It's common to report the Cext quantity.
Regarding of "easy ways to calculate other optical properties":
- in our case due to the sizes of the NP, I understand that absorption dominates, so extinction and absorption are really close.
- The code spits out potential and normal field on the surface. That could be called 'near field', but you could also compute it anywhere with BEM
- The forward scattering amplitude is one of the terms of the scattering matrix (the dominant one). Using the optical theorem you could calculate it (eq 5). We could mention this on the paper.
We can mention on what are other outputs you get from the code.

Reply

The extinction cross section is the sum of the absorption and scattering cross sections. The reason why we report the extinction cross section only is that with our model, it is not possible to separate these two components, however, for such nanoparticles smaller than 20 nm absorption dominates over the scattering (as detailed in Petryayeva, E., & Krull, U. J. (2011).. and Olson et al 2015), making Cext approximately equal to Cabs. We will mention this fact in the manuscript.

By default, the code returns a file with the values of the potential and its derivative on the surface of the scatterer. This could be used to compute the field anywhere in the domain, including the region close to the nanoparticle, however, it would include some modifications on the code.

Also, in the computation of the extinction cross-section, we use the forward scattering amplitude, as detailed in Equation 5. This is one entry of the scattering matrix, and can be back-calculated from Equation 5. With the current status of the code, it is not possible to compute other entries of the scattering matrix.

Modifications

@labarba This modification might need a check on the wording used.
Mention that absorption dominates when NP are smaller than 20 nm on 252807e

Comment 8

8) In the introduction of extinction, it is unclear if ‘a.k.a.’ is referring to both scattering and absorption or just to the second part of sentence (scattering).

Ideas

The phrase in the paper he/she is referring to is "When this happens, most of the incoming energy is either absorbed by the nanoparticle, or scattered in different directions, creating a shadow behind the scatterer (a.k.a., extinction)"

Reply

The "a.k.a., extinction" refers to the "shadow," produced by the combination of absorption and scattering.

Rewrite as: "When this happens, most of the incoming energy is either absorbed by the nanoparticle, or scattered in different directions, both creating a shadow behind the scatterer (a.k.a., extinction)"

Modifications

Add modification on commit ae87a75

Comment 9

9) Is Gamma defined at its first appearance on p. 2?

Ideas

It seems that \Gamma appears first in eq. 2 in page 2 but it's not defined and it's not in the text or picture what it is. We should say that it's the boundary.

Reply

Yes, the first appearance is on page two. We added an explanation after it is first introduced saying that it is a boundary between regions \Omega_1 and \Omega_2 .

Modifications

Fixed on commit 75c74a0

barbagroup / pygbe_lspr_paper