Computational nanoplasmonics in the quasistatic limit for biosensing applications

TL;DR

This paper extends the open-source PyGBe software to compute the localized surface plasmon resonance (LSPR) response of biosensors with metallic nanoparticles, demonstrating accurate, scalable simulations of sensor-analyte interactions in the quasistatic limit.

Contribution

It introduces a method for modeling biosensing using nanoplasmonics with an open-source tool, including validation and sensitivity analysis of nanoparticle-protein systems.

Findings

PyGBe accurately computes extinction cross-section with grid convergence.

A 0.5 nm red-shift in resonance frequency due to BSA proteins was observed.

The software handles large boundary element problems efficiently.

Abstract

This work uses the long-wavelength limit to compute LSPR response of biosensors, expanding the open-source PyGBe code to compute the extinction cross-section of metallic nanoparticles in the presence of any target for sensing. The target molecule is represented by a surface mesh, based on its crystal structure. PyGBe is research software for continuum electrostatics, written in Python with computationally expensive parts accelerated on GPU hardware, via PyCUDA. It is also accelerated algorithmically via a treecode that offers O(N log N) computational complexity. These features allow PyGBe to handle problems with half a million boundary elements or more. Using a model problem consisting of an isolated silver nanosphere in an electric field, our results show grid convergence as 1/N, and accurate computation of the extinction cross-section as a function of wavelength (compared with an analytical solution). For a model of a sensor-analyte system, consisting of a spherical silver nanoparticle and a set of bovine serum albumin (BSA) proteins, our results again obtain grid convergence as 1/N (with respect to the Richardson extrapolated value). Computing the LSPR response as a function of wavelength in the presence of BSA proteins captures a red-shift of 0.5 nm in the resonance frequency due to the presence of the analytes at 1-nm distance. The final result is a sensitivity study of the biosensor model, obtaining the shift in resonance frequency for various distances between the proteins and the nanoparticle. All results in this paper are fully reproducible, and we have deposited in archival data repositories all the materials needed to run the computations again and re-create the figures. PyGBe is open source under a permissive license and openly developed. Documentation is available at http://barbagroup.github.io/pygbe/docs/.

Figures

Tables

Table 1. Table 1: Grid-convergence study: Gauss quadrature points; K 𝐾 K and K f ​ i ​ n ​ e subscript 𝐾 𝑓 𝑖 𝑛 𝑒 K_{fine} are per element; N k subscript 𝑁 𝑘 N_{k} is per element edge (semi-analytical integration).

distant elements:	$K=4$
near-singular integrals:	$K_{fine}=37$
singular elements:	$N_k=9$

Table 2. Table 2: Grid-convergence study: treecode and solver parameters.

treecode order of expansion:	$P=15$
MAC	$\theta =0.5$
GMRES tolerance	${10}^{-5}$

Table 3. Table 3: Percentage error in the grid-convergence cases with an isolated silver nanosphere.

N	% error
$512$	$29.86$
$2048$	$7.33$
$8192$	$1.9$
$32768$	$0.52$

Table 4. Table 4: Verification: Gauss quadrature points; K 𝐾 K and K f ​ i ​ n ​ e subscript 𝐾 𝑓 𝑖 𝑛 𝑒 K_{fine} are per element; N k subscript 𝑁 𝑘 N_{k} is per element edge (semi-analytical integration).

distant elements:	$K=4$
near-singular integrals:	$K_{fine}=19$
singular elements:	$N_k=9$

Table 5. Table 5: Verification: treecode and solver parameters.

treecode order of expansion:	$P=6$
MAC	$\theta =0.5$
GMRES tolerance	${10}^{-3}$

Table 6. Table 6: Estimated percentage error of the BSA-sensor system (Fig.  3 ), with respect to the extrapolated value (using Richardson extrapolation).

N	% error
$512$	$29.39$
$2048$	$7.13$
$8192$	$1.82$
$32768$	$0.46$

Equations

$$\nabla \cdot {\mathbf{E}}_{1s}$$

$$\nabla \cdot {\mathbf{E}}_{2s}$$

$$({ϵ}_1{\mathbf{E}}_{1s}-{ϵ}_2{\mathbf{E}}_{2s})\cdot \mathbf{n}$$

$${\nabla }^2{\varphi }_{1s}$$

$${ϵ}_1\frac{\partial {\varphi }_{1s}}{\partial \mathbf{n}}-{ϵ}_2\frac{\partial {\varphi }_{2s}}{\partial \mathbf{n}}$$

$${\mathbf{E}}_{2s}=\frac{1}{4\pi {ϵ}_2}{k}^2\frac{e^{ikr}}{r}(\hat{\mathbf{r}}\times \mathbf{p})\times \hat{\mathbf{r}}.$$

$${\mathbf{E}}_{2s}(\mathbf{r}{)}_{r\to \infty }=\frac{e^{ikr}}{r}\mathbf{F}(\mathbf{k},{\mathbf{k}}_0),$$

$$C_{\text{ext}}=\frac{4\pi }{k^{\prime }}\mathrm{Im}\left[\frac{{\hat{\mathbf{e}}}_i}{|{\mathbf{E}}_i|}\mathbf{F}(\mathbf{k}={\mathbf{k}}_0,{\mathbf{k}}_0)\right].$$

$$k=k^{\prime }+i{k}^{\prime \prime }=\frac{2\pi }{\lambda }n,$$

$$\frac{{\varphi }_{1s,\Gamma }}{2}+K_L^{\Gamma }({\varphi }_{1s,\Gamma })-V_L^{\Gamma }\left(\frac{\partial }{\partial \mathbf{n}}{\varphi }_{1s,\Gamma }\right)=0$$

$$\frac{{\varphi }_{2s,\Gamma }}{2}-K_L^{\Gamma }({\varphi }_{2s,\Gamma })+V_L^{\Gamma }\left(\frac{\partial }{\partial \mathbf{n}}{\varphi }_{2s,\Gamma }\right)=0$$

$$V_L^{\Gamma }(\psi ({\mathbf{r}}_{\Gamma }))={\oint }_{\Gamma }\psi ({\mathbf{r}}_{\Gamma }^{\prime })G_L({\mathbf{r}}_{\Gamma },{\mathbf{r}}_{\Gamma }^{\prime })\text{d}{\Gamma }^{\prime },$$

$$K_L^{\Gamma }(\psi ({\mathbf{r}}_{\Gamma }))={\oint }_{\Gamma }\psi ({\mathbf{r}}_{\Gamma }^{\prime })\frac{\partial }{\partial \mathbf{n}}{G}_L({\mathbf{r}}_{\Gamma },{\mathbf{r}}_{\Gamma }^{\prime })\text{d}{\Gamma }^{\prime }.$$

$$G_L(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{1}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}$$

$$\frac{{\varphi }_{1s,\Gamma }}{2}+K_L^{\Gamma }({\varphi }_{1s,\Gamma })-V_L^{\Gamma }\left(\frac{\partial }{\partial \mathbf{n}}{\varphi }_{1s,\Gamma }\right)$$

$$\frac{{\varphi }_{1s,\Gamma }}{2}-K_L^{\Gamma }({\varphi }_{1s,\Gamma })+\frac{{ϵ}_1}{{ϵ}_2}{V}_L^{\Gamma }\left(\frac{\partial }{\partial \mathbf{n}}{\varphi }_{1s,\Gamma }\right)$$

$$\frac{{ϵ}_2-{ϵ}_1}{{ϵ}_2}{V}_L^{\Gamma }\left(\frac{\partial }{\partial \mathbf{n}}{\varphi }_{i,\Gamma }\right)\text{on Γ.}$$

$${\nabla }^2{\varphi }_{1s}$$

$${\nabla }^2{\varphi }_{3s}$$

$${ϵ}_1\frac{\partial {\varphi }_{1s}}{\partial \mathbf{n}}-{ϵ}_2\frac{\partial {\varphi }_{2s}}{\partial \mathbf{n}}$$

$${ϵ}_3\frac{\partial {\varphi }_{3s}}{\partial \mathbf{n}}-{ϵ}_2\frac{\partial {\varphi }_{2s}}{\partial \mathbf{n}}$$

$$\frac{{\varphi }_{1s,{\Gamma }_1}}{2}+K_{L,{\Gamma }_1}^{{\Gamma }_1}({\varphi }_{1s,{\Gamma }_1})-V_{L,{\Gamma }_1}^{{\Gamma }_1}\left(\frac{\partial }{\partial \mathbf{n}}{\varphi }_{1s,{\Gamma }_1}\right)$$

$$\frac{{\varphi }_{2s,{\Gamma }_1}}{2}-K_{L,{\Gamma }_1}^{{\Gamma }_1}({\varphi }_{2s,{\Gamma }_1})+V_{L,{\Gamma }_1}^{{\Gamma }_1}\left(\frac{\partial }{\partial \mathbf{n}}{\varphi }_{2s,{\Gamma }_1}\right)-K_{L,{\Gamma }_2}^{{\Gamma }_1}({\varphi }_{2s,{\Gamma }_2})+V_{L,{\Gamma }_2}^{{\Gamma }_1}\left(\frac{\partial }{\partial \mathbf{n}}{\varphi }_{2s,{\Gamma }_2}\right)$$

$$\frac{{\varphi }_{2s,{\Gamma }_2}}{2}-K_{L,{\Gamma }_1}^{{\Gamma }_2}({\varphi }_{2s,{\Gamma }_1})+V_{L,{\Gamma }_1}^{{\Gamma }_2}\left(\frac{\partial }{\partial \mathbf{n}}{\varphi }_{2s,{\Gamma }_1}\right)-K_{L,{\Gamma }_2}^{{\Gamma }_2}({\varphi }_{2s,{\Gamma }_2})+V_{L,{\Gamma }_2}^{{\Gamma }_2}\left(\frac{\partial }{\partial \mathbf{n}}{\varphi }_{2s,{\Gamma }_2}\right)$$

$$\frac{{\varphi }_{3s,{\Gamma }_2}}{2}+K_{L,{\Gamma }_2}^{{\Gamma }_2}({\varphi }_{3s,{\Gamma }_2})-V_{L,{\Gamma }_2}^{{\Gamma }_2}\left(\frac{\partial }{\partial \mathbf{n}}{\varphi }_{3s,{\Gamma }_2}\right)$$

$$\frac{{\varphi }_{1s,{\Gamma }_1}}{2}$$

$$-K_{L,{\Gamma }_2}^{{\Gamma }_1}({\varphi }_{3s,{\Gamma }_2})+V_{L,{\Gamma }_2}^{{\Gamma }_1}\left(\frac{{ϵ}_3}{{ϵ}_2}\frac{\partial }{\partial \mathbf{n}}{\varphi }_{3s,{\Gamma }_2}\right)-V_{L,{\Gamma }_2}^{{\Gamma }_1}\left(\frac{{ϵ}_2-{ϵ}_3}{{ϵ}_2}\frac{\partial }{\partial \mathbf{n}}{\varphi }_{i,{\Gamma }_2}\right)=0$$

$$\frac{{\varphi }_{3s,{\Gamma }_1}}{2}$$

$$-K_{L,{\Gamma }_2}^{{\Gamma }_2}({\varphi }_{3s,{\Gamma }_2})+V_{L,{\Gamma }_2}^{{\Gamma }_2}\left(\frac{{ϵ}_3}{{ϵ}_2}\frac{\partial }{\partial \mathbf{n}}{\varphi }_{3s,{\Gamma }_2}\right)-V_{L,{\Gamma }_2}^{{\Gamma }_2}\left(\frac{{ϵ}_2-{ϵ}_3}{{ϵ}_2}\frac{\partial }{\partial \mathbf{n}}{\varphi }_{i,{\Gamma }_2}\right)=0$$

$$\frac{{\varphi }_{3s,{\Gamma }_2}}{2}$$

$$V_{L,\text{disc}}^{{\mathbf{r}}_{\Gamma }}\left(\frac{\partial }{\partial \mathbf{n}}\varphi ({\mathbf{r}}_{\Gamma })\right)$$

$$K_{L,\text{disc}}^{{\mathbf{r}}_{\Gamma }}(\varphi ({\mathbf{r}}_{\Gamma }))$$

$$\left[\begin{array}{cc}\frac{1}{2}+K_L^{\Gamma }& -V_L^{\Gamma }\text{vspace}0.2cm\\ \frac{1}{2}-K_L^{\Gamma }& \frac{{ϵ}_1}{{ϵ}_2}{V}_L^{\Gamma }\text{vspace}0.2cm\end{array}\right]\left[\begin{array}{c}{\varphi }_{1s,\Gamma }\text{vspace}0.2cm\\ \frac{\partial }{\partial \mathbf{n}}{\varphi }_{1s,\Gamma }\text{vspace}0.2cm\end{array}\right]=\left[\begin{array}{c}0\\ V_L^{\Gamma }\left(\frac{{ϵ}_2-{ϵ}_1}{{ϵ}_2}\frac{\partial {\varphi }_i}{\partial \mathbf{n}}\right)\text{vspace}0.2cm\end{array}\right]$$

$$\left[\begin{array}{cccc}\frac{1}{2}+K_{L,{\Gamma }_1}^{{\Gamma }_1}& -V_{L,{\Gamma }_1}^{{\Gamma }_1}& 0& 0\text{vspace}0.2cm\\ \frac{1}{2}-K_{L,{\Gamma }_1}^{{\Gamma }_1}& \frac{{ϵ}_1}{{ϵ}_2}{V}_{L,{\Gamma }_1}^{{\Gamma }_1}& -K_{L,{\Gamma }_2}^{{\Gamma }_1}& \frac{{ϵ}_3}{{ϵ}_2}{V}_{L,{\Gamma }_2}^{{\Gamma }_1}\text{vspace}0.2cm\\ -K_{L,{\Gamma }_1}^{{\Gamma }_2}& \frac{{ϵ}_1}{{ϵ}_2}{V}_{L,{\Gamma }_1}^{{\Gamma }_2}& \frac{1}{2}-K_{L,{\Gamma }_2}^{{\Gamma }_2}& \frac{{ϵ}_3}{{ϵ}_2}{V}_{L,{\Gamma }_2}^{{\Gamma }_2}\text{vspace}0.2cm\\ 0& 0& \frac{1}{2}+K_{L,{\Gamma }_2}^{{\Gamma }_2}& -V_{L,{\Gamma }_2}^{{\Gamma }_2}\text{vspace}0.2cm\end{array}\right]\cdot \left[\begin{array}{c}{\varphi }_{1,{\Gamma }_1}\text{vspace}0.2cm\\ \frac{\partial }{\partial \mathbf{n}}{\varphi }_{1,{\Gamma }_1}\text{vspace}0.2cm\\ {\varphi }_{3,{\Gamma }_2}\text{vspace}0.2cm\\ \frac{\partial }{\partial \mathbf{n}}{\varphi }_{3,{\Gamma }_2}\text{vspace}0.2cm\end{array}\right]$$

$$=\left[\begin{array}{c}0\text{vspace}0.2cm\\ V_{L,{\Gamma }_1}^{{\Gamma }_1}\left(\frac{{ϵ}_2-{ϵ}_1}{{ϵ}_2}\frac{\partial }{\partial \mathbf{n}}{\varphi }_{i,{\Gamma }_1}\right)+V_{L,{\Gamma }_2}^{{\Gamma }_1}\left(\frac{{ϵ}_2-{ϵ}_3}{{ϵ}_2}\frac{\partial }{\partial \mathbf{n}}{\varphi }_{i,{\Gamma }_2}\right)\text{vspace}0.2cm\\ V_{L,{\Gamma }_1}^{{\Gamma }_2}\left(\frac{{ϵ}_2-{ϵ}_1}{{ϵ}_2}\frac{\partial }{\partial \mathbf{n}}{\varphi }_{i,{\Gamma }_1}\right)+V_{L,{\Gamma }_2}^{{\Gamma }_2}\left(\frac{{ϵ}_2-{ϵ}_3}{{ϵ}_2}\frac{\partial }{\partial \mathbf{n}}{\varphi }_{i,{\Gamma }_2}\right)\text{vspace}0.2cm\\ \frac{1}{4\pi {ϵ}_3}{\sum }_{k=0}^{N_q}\frac{q_k}{|{\mathbf{r}}_{{\Gamma }_2}-{\mathbf{r}}_k|}\text{vspace}0.2cm\end{array}\right]$$

$$V_{L,ij}^{\Gamma }$$

$$K_{L,ij}^{\Gamma }$$

$$\mathbf{p}={\int }_{\Omega }\mathbf{r}\rho \text{d}\Omega ,$$

$$\mathbf{p}=-{ϵ}_2{\int }_{\Omega }\mathbf{r}{\nabla }^2{\varphi }_{2s}\text{d}\Omega .$$

Full text

I Introduction

Localized surface plasmon resonance (LSPR) is an optical effect where an electromagnetic wave excites the free electrons on the surface of a metallic nanoparticle. The vibrations of the electron cloud are known as plasmons, and in LSPR they resonate with the incoming field (see Figure 1). When this happens, most of the incoming energy is either absorbed by the nanoparticle, or scattered in different directions, both effects creating a shadow behind the scatterer (a.k.a., extinction). In the case of nanoparticles smaller than 20 nm, absorption dominates and scattering contributions are negligible Petryayeva and Krull (2011); Olson et al. (2015). In LSPR, the wavelength of the incoming wave is often much larger than the size of the nanoparticle, which allows for valid approximations that simplify the mathematical model.

The phenomenon of LSPR can be used for biosensing, as the resonance frequency is highly dependent on the dielectric environment around the scatterer. The resonance frequency shifts whenever an analyte binds to the nanoparticle, resulting in a very sensitive means of detecting its presence Haes and Van Duyne (2002); Haes et al. (2004).

Numerical models for LSPR generally rely on the solution of Maxwell’s equations in some form, using finite difference time-domain (FDTD), boundary element, or finite element methods Solis et al. (2014). These methods have been used to study the optical properties of dielectric or metallic nanoparticles Hohenester (2018); Hohenester and Trugler (2012); Jung et al. (2010); Videen and Sun (2003); Mayergoyz et al. (2005); Mayergoyz and Zhang (2007), interactions between nanoparticles and electron beams García de Abajo and Aizpurua (1997); de Abajo and Howie (2002), and surface plasmon resonance sensors. In the latter application, researchers have used simple mathematical models for the interaction between a metallic nanoparticle and biomolecules, like representing the medium and the dissolved analytes with an effective permittivity Jung et al. (1998); Willets and Van Dyune (2007); Phan et al. (2013), or representing the target molecules as spheres Davis et al. (2010); Antosiewicz et al. (2011).

Progress in biosensor research is still predominantly made through experimental investigations, which can often be costly and time consuming. Computational approaches could assist the design process and play a role in optimizing biosensors, giving access to details that are not available in experimental settings. For example, empirical studies showed that the sensitivity of the sensor is highly dependent on the distance between the nanoparticle and the analyte Haes et al. (2004). These studies were complemented with models using a discrete dipole approximation (DDA), which includes the effect of the analyte through the effective permittivity. Other experimental studies complemented by modeling fully ignore the presence of the target molecules. For example, Beuwer et al. Beuwer et al. (2018) and Henkel et al. Henkel et al. (2018) used a boundary element method (BEM) in studies of the sensitivity of plasmonic sensors relying on (at least) two metallic nanoparticles (one on the sensor and one attached to the analyte). Explicitly including the target molecules in the model may be needed in some cases, however. For instance, despite experimental evidence showing that LSPR sensors are sensitive enough to detect conformational changes of the analytes Hall et al. (2011), these simplified models are not able to capture such details.

Even though LSPR is an optical effect, electrostatic theory provides a good approximation in the long-wavelength limit. This work uses the boundary integral electrostatics solver PyGBe Cooper et al. (2016) to compute the extinction cross-section of metallic nanoparticles, and to study how LSPR response changes in the presence of a biomolecule. We treat Maxwell’s equations quasi-statically Mayergoyz and Zhang (2007) and explicitly represent the target biomolecules by a surface mesh built from the crystal structures.

PyGBe is a Python implementation of continuum electrostatic theory, used for computing solvation energy of biomolecular systems. It has also been used to study protein orientation near charged nanosurfaces Cooper et al. (2015). The code was recently extended to allow for complex dielectric constants Clementi et al. (2017), aiming towards the LSPR biosensing application. The boundary element solver in PyGBe is accelerated algorithmically via a treecode—an $\mathcal{O}(N\log N)$ fast-summation method—and on hardware by taking advantage of graphic processing units (GPUs). With these features, PyGBe is able to easily handle problems with in the order of half a million boundary elements, or more, allowing for the explicit representation of the biomolecular surface. Other research software that could be used in this setting includes BEM++ Śmigaj et al. (2015) and a Matlab toolbox called MNPBEM Hohenester and Trugler (2012), which have the capability to solve the full Maxwell’s equations and the electrostatic approximation in the long-wavelength limit. We believe in both cases the size of problems they can solve, in terms of number of boundary elements, may not be enough to resolve the details of target biomolecules from their crystal structure.

The software is shared under the BSD 3-clause license and is openly developed via its repository on Github (https://github.com/barbagroup/pygbe). This study also follows careful reproducibility practices, and all materials necessary to reproduce the results are publicly available in reproducibility packages. We use the Figshare and Zenodo services to deposit the computational meshes, input and configuration files, and file bundles corresponding to the main figures in the paper. See the figure captions for references to the open data artifacts.

II Methods

The original implementation of PyGBe used continuum electrostatic theory to compute the solvation energy of biomolecular systems. In that setting, biomolecules are modeled as dielectric cavities inside an infinite continuum solvent, leading to a Poisson equation inside the molecules and Laplace or Poisson-Boltzmann in the solvent medium (with appropriate boundary conditions). This set of partial differential equations can be expressed with the corresponding boundary integral equation along the molecular interface, which PyGBe solves using a boundary element method Cooper et al. (2014, 2015).

The present work extends PyGBe to the LSPR biosensing application. In the long-wavelength limit, Maxwell’s equations can be approximated by a Laplace equation, which permits using the methods implemented in PyGBe, with modifications to allow for complex-valued permittivities, and to include the effect of an external electric field. This section describes the mathematical formulation for computing electromagnetic scattering in the long-wavelength setting, and develops the associated boundary integral equations and their discretized form.

II.1 Scattering of small particles

Electromagnetic scattering is usually modeled with Maxwell’s equations. When the wavelength of the incoming wave is much larger than the scatterer, these can be reduced to a quasi-static first-order approximation Mayergoyz and Zhang (2007):

[TABLE]

In Equation (II.1), $\mathbf{E}_{1s}$ and $\mathbf{E}_{2s}$ are the electric fields of the scattered wave in the nanoparticle and host regions, respectively (see Figure 2), $\mathbf{E}_{i}$ is the field of the incoming wave, and $\epsilon_{1}$ and $\epsilon_{2}$ are the permittivities. This approximation decouples the electric and magnetic fields, neglects the magnetic field, and describes the electric field as a curl-free vector field. Hence, we can reformulate Equation (II.1) with a scalar potential ($-\nabla\phi_{js}=\mathbf{E}_{js}$), as follows:

[TABLE]

Equation (II.1) is an electrostatic equation with an imposed electric field $\mathbf{E}_{i}=-\nabla\phi_{i}$, where $\Gamma$ is the boundary between regions $\Omega_{1}$ and $\Omega_{2}$.

II.2 Far-field scattering

In LSPR, the scattered electromagnetic wave is measured by a detector located far away from the scatterer (nanoparticle), and plasmon resonance is identified when the energy detected is minimum. In the far-field limit, the scattered field in the outside region ($\Omega_{2}$) is given by:

[TABLE]

where $k=2\pi/\lambda$ is the wave number and $\lambda$ the wavelength, $\mathbf{\hat{r}}$ is a unit vector in the direction of the observation point, and $\mathbf{p}$ is the dipole moment. We can obtain the scattered field using the scattering amplitude Jackson (1998):

[TABLE]

where $\mathbf{F}$ is the scattering amplitude, $\mathbf{k}$ is the scattered wave vector in the direction of propagation, and $\mathbf{k}_{0}$ the wave vector of the incident field.

II.3 Extinction cross-section and optical theorem

The extinction cross-section ($C_{\text{ext}}$) is a measure of the energy that does not reach the detector, either because of scattering in other directions, or absorption. This quantity is defined as the ratio between the lost energy and the intensity of the incoming wave, and has units of area. The extinction cross-section peaks at resonance of plasmons.

The extinction cross-section is related to the forward-scattering amplitude via the optical theorem. The traditional expression for this relationship applies for non-absorbing media Mayergoyz and Zhang (2007); Jackson (1998); Mishchenko Mishchenko (2007) corrected it for absorbing media, giving an expression that can be re-written using Jackson’s notation Jackson (1998) as follows:

[TABLE]

Here, $k^{\prime}$ is the real part of the complex wave number,

[TABLE]

and $n$ is the refraction index of the host medium.

Combining Equations (3) and (4), we can compute the scattering amplitude to then obtain the extinction cross-section with Equation (5).

II.4 The boundary element method

II.4.1 Electrostatic potential of a nanoparticle under an electric field

Integral formulation

Using Green’s second identity, the system of partial differential equations in Equation (II.1) can be rewritten as a system of boundary integral equations Brebbia and Dominguez (1992). Evaluating on the surface $\Gamma$, this becomes

[TABLE]

where $V$ and $K$ are the single- and double-layer operators, respectively:

[TABLE]

[TABLE]

Here, $G_{L}$ is the free-space Green’s function of the Laplace equation:

[TABLE]

Applying the interface conditions of Equation (II.1), leads to:

[TABLE]

II.4.2 Analyte-sensor electrostatic potential under an electric field

The sketch in Figure 3 shows a metallic nanoparticle ($\Omega_{1}$) interacting with an analyte ($\Omega_{3}$), under an external electric field. Mathematically, this situation can be modeled as

[TABLE]

where $q_{k}$ are the point charges of the atoms inside the protein, located at $\mathbf{r}_{k}$.

Integral formulation

Similar to Equation (II.4.1), we can write the system of partial differential equations in (II.4.2) as

[TABLE]

where $V$ and $K$ are the single- and double-layer operators in equations (8) and (9). In this case, we distinguish between the surface where the integrals run (subindex), and the surface that contains the evaluation point (superindex).

Applying the interface conditions of equation (II.4.2), leads to:

[TABLE]

Discretization and linear system

We discretize the surface into flat triangles, and assume that $\phi$ and $\partial\phi/\partial\mathbf{n}$ are constant within each element. We can then write the layer operators in their discretized form as follows:

[TABLE]

where $N_{p}$ is the number of discretization elements on $\Gamma$, and $\phi(\mathbf{r}_{\Gamma_{j}})$ and $\frac{\partial}{\partial\mathbf{n}}\phi(\mathbf{r}_{\Gamma_{j}})$ are the values of $\phi$ and $\frac{\partial\phi}{\partial\mathbf{n}}$ on panel $\Gamma_{j}$. Using centroid collocation, we can write equation (II.4.1) in matrix form as:

[TABLE]

Equation (II.4.2) can be represented as:

[TABLE]

where the elements of the matrix are

[TABLE]

with $\mathbf{r}_{\Gamma_{i}}$ being at the center of panel $\Gamma_{i}$.

Integral evaluation

We evaluate the integrals in Equation (II.4.2) with Gauss quadrature rules. The $1/r$ singularity of the Green’s function poses a problem to obtaining good accuracy when the integral is singular or near-singular. Therefore, we define three different regions, as follows.

Singular integrals:

If the collocation point is in the integration element, the singularity is difficult to resolve with standard Gauss integration schemes. In this case, we use a semi-analytical technique Hess and Smith (1967); Zhu et al. (2001) that places $N_{k}$ quadrature nodes on the edges of the triangle.

Near-singular integrals:

If the collocation point is close to the integration element, the integrand has a high gradient, and high-order quadrature rules are required. We use the representative length of the integrated triangle ($L=\sqrt{2\cdot\text{Area}}$) to define a threshold of the nearby region, for example, when the integration panel is $2L$ or less away from the collocation point. For near-singular integrals, we use $K_{fine}=19,25\text{ or }37$ points per triangle.

Far-away integrals:

When the distance between the collocation point and the integration element is beyond the threshold, they are considered to be far-away. At this point, the integrand is smooth enough that we obtain good accuracy with low-order integration, for example, with $K=1,3\text{ or }4$ Gauss quadrature points per boundary element.

II.4.3 Boundary integral expression of the dipole moment

As shown in Equation (3), the scattered electric field in the far-away limit depends on the dipole moment. The dipole moment is defined as

[TABLE]

and rewriting this equation using Gauss’ law, we obtain

[TABLE]

For component $i$, this becomes:

[TABLE]

Using the identity

[TABLE]

with $f=x_{i}$ and $\mathbf{v}=\nabla\phi_{2s}$, we can rewrite Equation (21) as

[TABLE]

and applying the divergence theorem

[TABLE]

Using the identity (22) again in Equation (23), this time taking $f=\phi_{2s}$ and $\mathbf{v}=\nabla x_{i}$, we get:

[TABLE]

Throughout this derivation, the normals are pointing into $\Omega_{1}$. However, in our implementation all normals are pointing outwards, and we need to include an extra negative sign, yielding:

[TABLE]

Using BEM, we obtain the electrostatic potential and its normal derivative, on the surface of the nanoparticle, which we use in Equation (25) to get the dipole moment, and in Equation (3) to obtain the scattered electric field. We can then use Equation (4) and Equation (5) to get the extinction cross section.

II.5 Acceleration strategies

One disadvantage of the Boundary Element Method (BEM) is that it generates dense matrices after discretization. Solving the resulting linear system using Gaussian elimination would require $\mathcal{O}(N^{3})$ computations and $\mathcal{O}(N^{2})$ storage, whereas for a Krylov-subspace iterative solver, like the Generalized Minimal Residual Method (GMRES), computations drop to $\mathcal{O}(N^{2})$ because they are dominated by dense matrix-vector products. This makes BEM inefficient with more than a few thousand boundary elements, which are the mesh sizes required for real applications.

In our formulation with Gaussian quadrature and collocation, the matrix-vector product becomes an $N$-body problem, with Gauss nodes acting as centers of mass (sources), and the collocation points acting as evaluation points for the potential (targets). To overcome the unfavorable scaling, we accelerate the matrix-vector product using a treecode algorithm Barnes and Hut (1986); Duan and Krasny (2001), which is a fast-summation algorithm capable of reducing $\mathcal{O}(N^{2})$ computational patterns like

[TABLE]

to a computational complexity of $\mathcal{O}(N\log N)$. In Equation (26) $q_{j}$ is the weight, $\psi$ the kernel, $\mathbf{y}_{j}$ the locations of sources and $\mathbf{x}_{i}$ the locations of targets.

The treecode groups sources geometrically in boxes of an octree, built ensuring that no box in the lowest level has more than $N_{\text{crit}}$ sources. If a group of sources is far away from a target, their influence is aggregated at an expansion center, and the target interacts with the box, rather than with each source independently. If the group of targets is close, the treecode queries the child boxes. If the box has no children and still is not far enough, the interaction is performed directly via (26). The threshold to decide if a box is far enough is called the multipole- acceptance criterion (MAC), defined as:

[TABLE]

where $r_{b}$ is the box size and $r$ the distance between the box center and the target. Common values of $\theta$ are $1/2$ and $2/3$. To approximate the contribution of the sources, we use Taylor expansions of order $P$. The treecode allows us to control the accuracy of the approximation by modifying $\theta$ and $P$. Further details of the treecode implementation in PyGBe can be found in Cooper and Barba (2013); Cooper et al. (2014).

II.6 Code modifications and added features

As mentioned at the beginning of this section, the present work extends the PyGBe code to allow its application to nano-plasmonics. The code required the following modifications and added features:

• •

Re-writing the GMRES solver to accept complex numbers.

• •

Splitting treecode calculations into real and imaginary parts.

• •

Re-formatting configuration files to include electric field intensity and wavelength.

• •

Adding the new function read_electric_field, to read the electric field intensity and its wavelength from configuration files.

• •

Adding the new function dipole_moment to compute numerically the dipole moment by Equation (25).

• •

Adding a new function to compute the extinction cross section (extinction_cross_section).

• •

Organizing LSPR computations on a different main script (called lspr.py).

For information about how to use the code, run examples and tests, see the PyGBe documentation at http://barbagroup.github.io/pygbe/docs/

II.7 Protein mesh preparation

In Figure 3, $\Omega_{3}$ is a region that represents the analyte molecule, which contains a point charge distribution of the partial charges, and is interfaced with the solvent by $\Gamma_{2}$, the solvent excluded surface (SES). The SES is generated by rolling a spherical probe of the size of a water molecule ($1.4$Å  radius) around the analyte, and tracking the points where the probe and molecule make contact. The open-source software Nanoshaper Decherchi and Rocchia (2013) uses the molecular structure to produce a triangulation of the SES, which can be read by our software. In particular, Nanoshaper takes as inputs the atomic coordinates, obtained from the Protein Data Bank, and radii, which were extracted from a pqr file generated with pdb2pqr Dolinsky et al. (2004). We obtained the charge and van der Waals parameters of the analyte from pdb2pqr using the built-in amber force field. In support of the reproducibility of our results, we deposited the final meshes in the Zenodo data repository. See section III.3 for details.

III Results

We present results for two kinds of problems. The first is a model problem for which an analytical solution is available, allowing for a grid-refinement study and code verification using that solution. It consists of a spherical nanoparticle in a constant electric field, where the extinction cross-section can be derived in closed form. The second set of results use a model for a biosensor detecting a target molecule, via frequency shifts in the plasmon resonance of a metallic nanoparticle. In this case, since an analytical solution is not available, we can use Richardson extrapolation to estimate the errors in a grid-refinement study. We also computed the variation of the extinction cross-section with respect to wavelength for the isolated nanoparticle, and in the presence of bovine serum albumin (BSA) proteins, varying the location of the analytes. The final result is a sensitivity study of the biosensor model, looking at how the peak in frequency response varies with distance of the protein to the nanoparticle.

All results were obtained on a lab workstation, built from parts. Hardware specifications are as follows:

• •

CPU: Intel Core i7-5930K Haswell-E 6-Core 3.5GHz LGA 2011-v3

• •

RAM: G.SKILL Ripjaws 4 series 32GB (4 x 8GB)

• •

GPU: Nvidia Tesla K40c (with 12 GB memory)

III.1 Grid convergence and verification with an isolated silver nanoparticle

In the long-wavelength limit, the electrostatic approximation applies and the electromagnetic scattering of a small spherical particle can be modeled by a sphere in a constant electric field. Figure 4 illustrates this scenario.

This model problem has an analytical solution, which allows us to compare with the numerical calculations of the extinction cross-section obtained with PyGBe, for code verification and grid-convergence analysis. Mishchenko Mishchenko (2007) derived the following analytical result, valid for lossy mediums:

[TABLE]

Here, $a$ is the radius of the sphere, $k$ is the complex wave number ($k=k^{\prime}+ik^{\prime\prime}$), $\epsilon_{p}$ is the dielectric constant of the particle, and $\epsilon_{m}$ is the dielectric constant of the host medium. If the medium is not lossy, then $k^{\prime\prime}=0$ and $k=k^{\prime}$.

We completed a grid-convergence study of PyGBe for the extinction cross-section of a spherical silver nanoparticle of radius 8 nm immersed in water, under a $z$-polarized electric field with a wavelength of 380 nm and intensity of $-0.0037e/({\AA}^{2}\,\epsilon_{0})$. In these conditions, water has a dielectric constant of $1.7972\,+\,8.5048^{-09}i$ Johnson and Christy (1972) and silver of $-3.3877\,+\,0.1922i$ Hale and Querry (1972). Table 1 lists the Gauss quadrature points used for each type of boundary element. The threshold parameter defining the near-singular region was 0.5 (refer to the PyGBe documentation, under “Parameter file format”). Table 2 shows the treecode and solver parameters for this grid-convergence study.

The results are shown in Figure 5, where the mesh sizes are 512, 2048, 8192, and 32768 elements. The analytical solution with equation (28) is $C_{ext}=1854.48$ nm2, and the computed errors are as shown in Table 3. The dashed line in Figure 5 shows a $1/N$ slope, and the observed order of convergence is $0.98$, evidence that the meshes are correctly resolving the numerical solutions with PyGBe.

As another verification test of PyGBe in the LSPR setting, we computed the extinction cross-section of an isolated sphere for a range of wavelengths. The results are shown in Figure 6, comparing with the analytical solution. The values of the dielectric constant for each wavelength were obtained by interpolation of experimental data Johnson and Christy (1972); Hale and Querry (1972). For reproducibility of these results, we provide a Jupyter notebook with the code used for this interpolation step. See section III.3 for details. We used a mesh with $N=32,768$, and relaxed some parameters compared with the grid-convergence results shown previously, still yielding errors below $1\%$ at all frequencies. This results in a $12\times$ decrease in the runtime for each case. The parameters used are shown in Tables 4 and 5. Figure 6 shows good agreement between the computed and analytical results, evidence that PyGBe can accurately represent the mathematical model.

III.2 LSPR response to bovine serum albumin (BSA)

Localized Surface Plasmon Resonance (LSPR) biosensors detect a target molecule by monitoring frequency shifts in the plasmon resonance of metallic nanoparticles, in presence of an analyte Willets and Van Dyune (2007). In this section, we model the response of LSPR biosensors using the expanded capacity of PyGBe. We consider a spherical silver nanoparticle, and compute the extinction cross-section placing bovine serum albumin (BSA) proteins (PDB code: 4FS5, a BSA dimer) in different locations. We placed two BSA dimers opposite to each other in three configurations ($\pm z$, $\pm y$, and $\pm x$), as shown by figures 8 and 12. As a point of comparison, experiments by Teichroeb and co-workers Teichroeb et al. (2008) find a coverage of $2\times 10^{12}\quad\text{molecules}/cm^{2}$ or $3.3\times 10^{12}\quad\text{molecules}/cm^{2}$, with a gold sphere 15-nm in diameter. In that work, the molecular size reported is 5.5 nm$\times$5.5 nm$\times$9 nm, resulting in a number of attached molecules between 4 and 6. The BSA molecule used in our work corresponds to a dimer, i.e., approximately double the size of that in Teichroeb et al.’s experiment. With two BSA dimers in the proximity of the sensor, the volume fractions in the near-by region are comparable.

III.2.1 Grid-convergence study

We performed a grid-convergence analysis of the system sketched in Figure 3. Since we compute the extinction cross-section of the spherical nanoparticle only, we set a fixed mesh density for the protein and refined the mesh of the sphere (meshes of 512, 2048, 8192 and 32768 elements). We found that the protein meshed with two triangles per ${\AA}^{2}$ was fine enough for the convergence analysis, resulting in $N_{prot}=98116$ elements.

We used the same physical conditions as in the grid convergence with an isolated silver nanoparticle, and the same numerical parameters, presented in Tables 1 and 2. For the protein dielectric constant, we used $2.7514+0.2860i$, obtained from the functional relationship provided by Phan, et al. Phan et al. (2013). The distance between the sensor and the analyte was $d=1$ nm, and the BSA protein was oriented such that its dipole moment was aligned with the $y$-axis. To obtain the error estimates shown in Figure 7 and Table 6, we used the Richardson extrapolated value of extinction cross-section as a reference, $C_{ext}=1778.73$ nm2.

The observed order of convergence is $0.99$, and Figure 7 shows that the error decays with the number of boundary elements ($1/N$), which is consistent with our verification results in Section III.1. This provides evidence that the numerical solutions computed with PyGBe are correctly resolved by the meshes. The percentage errors for the different meshes are presented in Table 6.

III.2.2 Resonance frequency shift

We computed the LSPR response as a function of the wavelength in the presence of the BSA protein. To optimize run-times without compromising accuracy, we used a relaxed set of parameters, where the protein mesh density was one element per ${\AA}^{2}$ ($N_{prot}=45140$) and the sphere mesh had $N_{sensor}=32768$ elements. These calculations used the same parameters as shown in Tables 4 and 5. This parameter choice resulted in a percentage error below 1%, with respect to the Richardson-extrapolated value. The run time for each one of these cases was approximately $7.5$ min using one NVIDIA Tesla K40c GPU. When two proteins are present, the run time per case is approximately $15$ min.

Figure 8 shows a visualization of the meshes for these calculations, with two BSA proteins placed at a distance $d=1$ nm away from a spherical silver nanoparticle, along the $z$-axis. The surface-mesh data, plotting scripts and figure are available openly on Figshare, in support of the paper’s reproducibility Clementi et al. (2018c). The position of the BSA molecule in the $+z$ axis was the same as in the convergence analysis in Section III.2.1, whereas the BSA in the $-z$ position is a 180∘ solid rotation about the $y$-axis of the BSA in $+z$. We performed calculations for wavelengths between $382$ nm and $387$ nm, every $0.25$ nm, which are around the peak seen in Figure 6.

Figure 9 shows the variation of the extinction cross-section with respect to wavelength for the isolated nanoparticle ($d=\infty$) and with BSA proteins placed $d=1$ nm away. The result shows a red-shift ($0.5$ nm) in the resonance frequency due to the presence of the BSA analytes.

To study the effect of location of the analytes, we re-computed the result placing the BSA proteins along the $x$- and $y$-axis, at $\pm 1$ nm, as shown in Figure 12. These configurations were obtained via a 90-degree solid rotation of the $z$-configuration (Figure 8) along the $x$- and $y$-axis, respectively. Figure 10 shows the results, in each case.

III.2.3 Sensitivity calculations

The sensitivity of an LSPR biosensor corresponds to the relationship between the size of the resonance frequency shift and the number of analytes bound to the sensor (through a ligand). Experiments show that the distance between the nanoparticle and the analyte affects the sensitivity of the sensor, to the point that targets placed $15$ nm away from the surface are very hard to detect Haes et al. (2004). This is a critical issue, considering that common ligands (for example, antibodies) can be larger than $15$ nm. Figure 11 shows how the peak varies with the distance at which the analytes ($+z$ and $-z$) are placed. In particular, we see a shift of $0.25$ nm when $d=2$ nm to $0.75$ nm when the analytes are placed at $d=0.5$ nm. The parameters used in this case remain the same as the ones used in Figures 9 and 10 .

III.3 Reproducibility and data management

To facilitate the reproducibility and replication of our results, we consistently release our research code and data with every publication. PyGBe is openly developed and shared under the BSD3-clause license via its repository at https://github.com/barbagroup/pygbe.

We also release all of the data and scripts needed to run the calculations reported in this work, as well as the post-processing scripts to reproduce the figures in this paper. All the input files necessary to reproduce the computations are available in one Zenodo data set Clementi et al. (2018d). Each problem corresponds to a folder, wherein the user can find geometry files (surface meshes), configuration files, parameter files, and when it applies, the protein charges (.pqr). The scripts and auxiliary files needed to run PyGBe to re-generate every result in the paper are collected in another Zenodo deposit Clementi et al. (2018e). After execution, the resulting data needed to re-create the figures in the paper will be saved in the running folder and the input files (the first Zenodo set) can at that point be deleted. (For more details, the reader can consult a README file in the Zenodo archive.) Reproducibility packages to reproduce the figures in the paper are deposited on Figshare, including the figures, plotting scripts and Jupyter notebooks that organize and re-create the results Clementi et al. (2018a, b, c, f).

IV Discussion

Extending PyGBe to the LSPR biosensing application required considerable code modifications and added functionality. The results presented in the previous section offer evidence to build confidence on the suitability of the mathematical model and the correctness of the code. The grid-convergence study with a nanosphere under a constant electric field shows a $1/N$ rate of convergence, consistent with convergence results in previous work using PyGBe Cooper et al. (2014). Further verification of PyGBe’s new ability to compute extinction cross-section of a scatterer in the long-wavelength limit is provided in Figure 6. The computed extinction cross-section of a silver nanoparticle in a range of frequencies is within 1% of the analytical value, with the numerical parameters chosen. This level of accuracy is likely sufficient, given that experimental uncertainty in the values of the dielectric constant for silver is in the order of 1%, also Johnson and Christy (1972).

Figure 9 shows a red shift of the plasmon resonance frequency peak in presence of the BSA proteins. Experimental observations of Tang, et al. Tang et al. (2010) with silver nanoparticles of approximately 17 nm in diameter and BSA proteins in solution revealed a red shift upon adding the proteins. Similar to the effect we see with our model, they observed as well a decrement of the peak amplitude. Moreover, recent experiments Pu et al. (2018) report a resonance frequency for a silver nanoparticle in the presence of BSA proteins of between 380 and 400 nm, which is consistent with our results. Other experiments Raphael et al. (2013) also report a red shift in the resonance frequency in the presence of (different) proteins. Our boundary element method approach using electrostatic approximation is thus able to capture the characteristic resonance-frequency shift of LSPR biosensors.

With the electric field aligned in the $z$-direction, placing the proteins at a distance in the $x$ or $y$ directions from the nanoparticle shows a negligible shift in the resonance peak: the shifts in Figure 10 are smaller than the resolution between wavelengths ($<0.25$ nm). This finding is consistent with the free electrons oscillating in the $z$ direction under a $z$-polarized electric field, and not in the $x$ and $y$ directions (see Figure 1). The analytes have a marked effect when placed in the $z$ direction, where they can interfere with the free oscillating electrons.

Figure 11 shows how the shift in resonance frequency varies with the distance between the sensor and the analyte. As expected, the shift decays as the BSA moves away from the sensor, to the point that if the BSA proteins are placed $d=2$ nm away, the shift is only $0.25$ nm. This result shows the potential of PyGBe and the electrostatic approach to study biosensor sensitivity with distance. Note that possible quantum effects (e.g., tunneling) at $d=0.5$nm are ignored with our classical model. Even if this distance could be close to or in the quantum regime, evidence that classical theory is valid at this distance in similar systems has been reported Savage et al. (2012); Esteban et al. (2012).

Even though there is evidence that techniques such as Plasmon Enhanced Raman Scattering are capable of detecting all the way to single molecules Zhang et al. (2013), as far as we know, there is no evidence of purely LSPR approaches that can sense such low concentration of analytes. These computational studies can shine light on potential improvements that would enhance sensitivity of LSPR biosensors, for example, by using smaller ligands.

We are not aware of other LSPR simulations where the molecular details of the analyte are considered, however, similar calculations could be performed with other software. For example, BEM++ Śmigaj et al. (2015) also models the system as a set of boundary integral equations, discretized in flat triangular panels. This software uses the Galerkin approach and algorithmic acceleration via hierarchical matrices, which is slower and less memory efficient than the treecode and limits the accessible problem sizes. The Matlab toolbox MNPBEM Hohenester and Trugler (2012) is another alternative software designed to simulate scattering of metallic nanoparticles. Its BEM implementation is similar to PyGBe as it uses a centroid collocation scheme on flat triangular panels, but differs in the algorithmic acceleration technique, which is also based on hierarchical matrices rather than a treecode. This results in higher memory usage compared to our code, making it harder to simulate large analytes in detail. Commercial finite-element or finite-difference solvers could also be used in this application, for example, COMSOL. These volumetric approaches, however, struggle to correctly impose the zero boundary condition at infinity, which is exactly met for a BEM formulation.

V Conclusion

In this work, we combined the implicit-solvent model of electrostatics interactions in PyGBe with a long-wavelength representation of LSPR response in nanoparticles. We extended PyGBe to work with complex-valued quantities, and added functionality to include an imposed electric field and compute relevant quantities (dipole moment, extinction cross-section). Previous work with PyGBe showed its suitability for computing biomolecular electrostatics considering solvent-filled cavities and Stern layers Cooper et al. (2014), and for protein-surface electrostatic interactions Cooper and Barba (2016). This latest extension can offer a valuable computational approach to study nanoplasnomics and aid in the design of LSPR biosensors. Thanks to algorithmic acceleration with a treecode, and hardware acceleration with GPUs, PyGBe is able to compute problems with half a million elements, or more, which is required to represent the molecular surface accurately.

Acknowledgements.

CDC acknowledges the financial support from CONICYT through projects FONDECYT Iniciación 11160768 and Basal FB0821.