Improved method for estimating adlayer thickness and bulk RI change for gold nanocrescent sensors

This paper presents a novel method employing the maximum likelihood estimation (MLE) technique alongside a nonlinear sensor response model to improve and extract more quantitative sensing results for localized surface plasmon resonance biosensors. The nonlinear response model treats the sensor response as a nonlinear function of the biomolecular adlayer thickness. This method makes use of the multiple resonance characteristic of nanocrescent structures in order to estimate the adlayer thickness and bulk refractive index (RI) change. Nanoimprint lithography is used here to fabricate the nanostructures. The finite element method (FEM) is used to model the nanocrescents and numerically validate the nonlinear-MLE method. Comparing to the established linear model, the proposed nonlinear-MLE method achieves 75% improvement in the limit of detection based on the estimated adlayer thickness and improves the bulk RI resolution by two orders of magnitude.

channel. This platform was shown to reduce the interfering effects, but it does not provide quantitative information about the binding molecules.
U-shaped nanostructures have been previously used in differentiating between specific and nonspecific binding by extracting information from two or more LSPR modes 9 . The study assumed that specific binding mostly occurs at the metal surface, and the nonspecific binding occurs at a distant location (on the substrate), and ignored any non-specific binding that may occur at the metal surface. However, the method requires repetitive simulation in order to determine the sensitivity matrix based on other biological samples, which is a practical limitation of the method. Alternatively, a model that distinguishes bulk RI and adlayer thickness changes is a more practical solution, as it can decouple the effects associated with them. The sensor response-at each resonance-is related to the adlayer thickness and bulk RI changes,and the effects can be determined by solving the two equations (corresponding to the number of resonances). This linear response model has been previously employed for an SPR sensor with a dielectric overlayer to excite two resonances 10 , an SPR sensor based on simultaneous excitation of short and long range SPR modes 11,12 , a dual-resonance SPR sensor with different penetration depths 13 , simultaneous excitation of transverse and longitudinal modes of nanorods 14 , and a sensor based on surface plasmon resonance and plasmon waveguide resonance 15 . This model, however, requires a sensitivity matrix (including the adlayer and bulk RI sensitivities for both resonances) with a low condition number to avoid any numerical errors, which may not be the case for many sensing platforms. We have improved this method by applying the MLE method to a set of results based on three-resonance nanorod structures, estimating the adlayer thickness and the bulk RI change 16 . Although the accuracy and precision have been improved with self-referencing capabilities, the method remains valid only for extremely low adlayer thickness (<l d /10). This article presents a method based on the maximum likelihood estimation (MLE) technique alongside a nonlinear response model to improve the accuracy and precision of the estimated quantities. This indicates a considerable improvement in the RI resolution and the limit of detection and sensing dynamic range (<l d /2). Herein the proposed method is applied to nanocrescent structures as they can be fabricated using cost-effective fabrication methods and they feature a multiple resonance absorption spectrum 17 . However, the method can be applied to any multiple resonance sensor.

Nonlinear Model for Sensor Response
The electromagnetic field decays exponentially from the surface of the nanostructures. If the maximum EM field is located at z = 0, it decays with a factor of exp(−z/l d ) along the z direction, where l d is the electromagnetic (EM) decay length, see Fig. 1 Since the intensity is the square of the electric field strength, it decays with a factor of exp(−2z/l d ) with respect to the nanostructure surface (z). This suggests the following for the effective refractive index Accordingly, the effective refractive index along the z direction is weighted by using the same decay factor exp(−2d/l d ) as shown in Fig. 1. The effective refractive index can be obtained by integrating the refractive indices along z direction as follows 18 where n(z) changes along the z direction as where n a and n B are the analyte and bulk refractive indices, respectively. The resonance shift can be determined by where S B is the bulk RI sensitivity. From equation (2), substituting for n eff , the following is obtained for the resonance shift This can be used to estimate the the adlayer thickness d and bulk RI change Δn B . First, the bulk RI sensitivity S B is measured based on bulk RI change for ethanol solutions of known refractive indices. The nanocrescents support three resonances, but only two resonances are used here. The base solution is injected to obtain the sensing baselines for these resonances, and the shifts in the resonance wavelengths Δλ i are tracked in real time as the biological sample is introduced.
The maximum response can be obtained when a thick adsorbate layer is reached d → ∞, and equation (4) becomes Dividing equation (4) by equation (5), the following is obtained Equation (6) can be rearranged in order to obtain an expression that is linear in adlayer thickness as follows For a nanostructure with multiple resonances, the electromagnetic decay length is dependent on the resonance wavelength, and the bulk RI sensitivity. Considering a two-resonance system and including the errors associated with the measured data, equation (7) can be written as follows , C d i is a coefficient related to the electromagnetic field decay length for the i th resonance calculated as − l 2/ d i , and ε y i represents the error due to noise which is directly related to the variance of the measured quantities R y over the measured shifts in resonance wavelengths. In matrix notation where R 12 represents the covariance accounting for correlated noise. Considering a multiple resonance response, y 1 , y 2 , ..., y i , and a normal distribution for the noise effect, the likelihood of obtaining these measured quantities (y i ) given the true value (d) and the variance R y can be obtained by multiplying the normal distributions  of these estimates. In matrix notation this can be expressed as follows According to the MLE, the maximum likelihood of estimating the true value is obtained when the derivative of the above likelihood with respect to the true value approaches zero 19 . For simplicity, we obtain the log of the above likelihoods as follows Now, the true value can be estimated such that the derivative of the log likelihood with respect to this true value is equal to zero 19 .
This can be solved for the estimate d In the case of the dual-resonance nanocrescents (i = 2), the adlayer thickness can be estimated as follows www.nature.com/scientificreports/  If the noise is uncorrelated (R y12 = 0), the estimate becomes Likewise, the RI change can be estimated with an improved accuracy using the MLE method. From equation (4), the shifts in the resonance wavelengths are given by 2 / di represents the sensitivity coefficient for the i th resonance, and d is the value estimated by equation (13) or (14). ε λ i is the error due to noise, represented by the variance associated with the i th measured wavelength shift. The estimated Δ  n that maximizes the likelihood of the above probability assuming a correlated noise can be obtained using the MLE

Simulated Results and Validation of the Estimation Method
This section compares the nonlinear-MLE method with the established linear response model based on FEM simulated results, and presents a FEM evaluation for the accuracy based on each method with respect to deviated resonance wavelengths (simulating the effect of noise on each resonance wavelength). The simulated results show distinct values for the sensitivity to bulk RI and adlayer thickness changes. Inspecting both modes in terms of the adlayer and bulk RI sensitivities can provide an insight into the sensing efficiency for each mode. It is evident that the adlayer thickness, the adlayer refractive index n a and the bulk refractive index n B contribute to the effective refractive index n eff [Equation (1)]. Therefore, the adlayer sensitivity (S d ) can be given by 20 eff eff where ∂λ/∂d and ∂λ/∂n eff are the adlayer sensitivity S d and bulk RI sensitivity S B , respectively. This equation suggests that the adlayer sensitivity is determined by the bulk RI sensitivity and the rate at which the adlayer thickness contributes to the effective refractive index n eff . Using equation (17) we can define the adlayer sensing efficiency η of a given biosensor as follows This determines how effectively a given biosensor can detect changes in adlayer thickness. Referring to Figs 2 and 3, the first mode (1200 nm) exhibits a lower bulk RI sensitivity, comparing to that associated with the second mode ~1700 nm (325.25 nm/RIU versus 787.35 nm/RIU). These modes yield adlayer sensitivities of 1.47 and 2.2, respectively. Therefore, the calculated adlayer sensing efficiency η is 4.5 × 10 −3 RIU/nm and 2.8 × 10 −3 RIU/nm. In other words, changing the adlayer thickness by 1 nm would alter the effective refractive index by 4.5 × 10 −3 RIU based on the first mode and 2.8 × 10 −3 RIU in the case of the second mode, although the second mode features a higher bulk RI sensitivity than that of the first mode. We can now consider the established linear response model for comparison, the LM method relates each resonance wavelength shift to the changes in RI and thickness of an adsorbed biological material by the following relationship From equation (19), the adlayer thickness and bulk RI change can be calculated as The nonlinear-MLE method is now employed to estimate the bulk RI change and adlayer thickness-that were used in the simulation-by adding uncertainties to the simulated shifts of the resonance wavelengths ( σ ± λ i ) such as λ σ Δ = λ / 10 i i , as shown in Fig. 4(a). The resonance wavelength shifts are used to determine the corresponding values for ln(1−Δλ/Δλ max ), which are then used in equation (14) to determine the adlayer thickness, as shown in Fig. 4(b). The adlayer thickness is then used to determine the sensitivity coefficient, C n = S[1−exp(−2d/l d )], to estimate the bulk RI change using equation (16). The estimated adlayer thickness and bulk RI change are shown in Fig. 4(c). The results obtained based on the LM are also shown in Fig. 4(b,c). The nonlinear-MLE method revealed the following for the estimates: d = 5.95 nm, Δn B ≈ 0, with 0.47 nm and 1.4 × 10 -3 RIU uncertainties, respectively. Under the same conditions, the linear response model revealed 5.5 nm and 1.3 × 10 −3 for the estimated adlayer thickness and bulk RI change with uncertainties of 1.8 nm and 6 × 10 −3 , respectively. This suggests that the nonlinear-MLE can improve the accuracy of the estimated adlayer thickness by one order of magnitude, and the precision by a factor of 4, as shown in Fig. 4. Figure 5 compares both methods in terms of accuracy and precision based on adlayer thickness ranging from 6 nm to 25 nm. The percentage error in the estimated adlayer thickness ranges from 0.83-1.96% based on the nonlinear-MLE method, and 8.3-71.6%, based on the linear response model. This indicates that the nonlinear-MLE method can improve the results by a factor of 36 as compared to those based on the linear response model, when the adlayer thickness approaches ∼l /2 d . The error in the estimated RI change tends to be negligible based on the nonlinear-MLE method and increases drastically based on the linear response model, as shown in Fig. 5(b). The nonlinear-MLE method and the linear response model reveal 5 × 10 −3 and 1.5 × 10 −2 RIU uncertainties when the adlayer thickness approaches 25 nm, as shown Fig. 5(b).

Experimental Results
This section provides the experimental results based on the Bioten-Streptavidin binding events. The samples were prepared according to the established protocol 9,22 . Before being used in the sensing experiments, the fabricated nanocrescent structures were cleaned by DI water and ethanol solution, blown dry with nitrogen,and plasma treated to remove any biological contaminant. The substrate was then functionalized by biotin-hpdp, and streptavidin solution was prepared in tris-buffer according to the established surface chemistry protocol 9,22 . The bulk RI sensitivity at each resonance was first determined based on resonance shift against RI change due to changing the concentration in ethanol solution. The sensitivities to adlayer thickness were then determined based on the measured results by correcting the simulated counterparts based on the measured bulk RI sensitivities. Since the adayer sensitivity is related to the bulk RI sensitivity by S d = S B (n a −n B )/l d . The corrected adlayer sensitivity ′ S d can be determined as follows where S d is the simulated adlayer sensitivity and S B and ′ S B are the simulated and measured buk RI sensitivities, respectively. The bulk RI sensitivity is used along the EM decay length in equations (14) and (16) to estimate the adlayer thickness and bulk RI change, respectively, based on the nonlinear-MLE method. This section also compares these results to those obtained based on the linear model. The linear model uses the bulk and adlayer sensitivity factors in equation (19) to estimate the adlayer thickness and bulk RI change. Figure 6(a) shows the shifts in the resonance wavelengths, based of the measured extinction curves for streptavidin biotin sensing experiments. These shifts were translated into the logarithmic normalized quantity, ln(1 − Δλ/Δλ max ), and used in equations (14) and equation (16) to estimate the adlayer thickness and bulk RI change, as shown in Fig. 6(b,c). Now the linear model is used based on the same measured data, shown in Fig. 6(a) in estimating the adlayer thickness and bulk RI change using the sensitivity factors in equation (19).  Figure 6(c) shows the results based on the linear model. The limit of detection and bulk RI resolution can be determined based on the standard deviation of the estimated adlayer thickness and bulk RI change, respectively. Based on these results, the nonlinear-MLE method improves the accuracy and precision as compared with the linear model. The limit of detection for adlayer thickness is reduced by a factor of 4, and the bulk RI resolution is improved by a factor of 22 based on the bulk RI change. Table 1 compares between the proposed nonlinear-MLE method and the linear model based on the estimated adlayer thickness and measured sensing characteristics.

Discussion
The paper presented a statistical method, combining the the MLE method with a nonlinear model, for extracting the adlayer thickness and bulk RI change with improved accuracy. The method not only provided a quantitative information about the binding events, but also improved the precision of LSPR sensors. The paper adopted a nonlinear model for the sensor response, to adlayer thickness change, which reflects the real scenario of typical sensing experiments. The model follows a similar trend to that based on the association constant in chemical sensors.
The nonlinear model uses the EM field decay length and sensitivity for each resonance to estimate the adlayer and bulk RI change, whilst the linear model uses the sensitivity to bulk RI change and the adlayer sensitivity that needs to be recalculated (corrected) for other target biomolecules. The latter represents a substantial disadvantage of the linear model. In contrast to the linear model, the nonlinear model avoids the redundant calculation of the sensitivity to adlayer thickness for various biomolecular adlayer. Based on the simulated and measured result, and compared to the linear model, the nonlinear-MLE method improved the bulk RI resolution by two orders of magnitude (6 × 10 −5 RIU vs 1.3 × 10 −3 RIU). The method also achieves 75% improvement in the limit of detection based on the adlayer thickness (0.06 nm vs 0.24 nm uncertainty in the estimated adlayer thickness). In this paper, we considered the nanocrescent structures, but the proposed method can be applied to a wide range of different structures and various sensing scenarios.

Methods
The finite element method (FEM) was used to calculate the sensor response to bulk RI and adlayer thickness changes and validate the MLE method based on the nonlinear model. Periodic boundary conditions were enforced along the sides of a hexagonal simulation domain created with commercial COMSOL Multiphysics  simulation package, as shown in Fig. 7(a). The refractive index for gold was obtained from Johnson and Christy 23 . The nanocrescents and the adlayer were discretized using a tetrahedral mesh, and triangular mesh was used with the rest of the domain. The transmission efficiency was obtained from the scattering parameter S 21 , as input and output ports were assigned to the bottom and top surfaces of the simulation domain, shown in Fig. 7(a). The nanostructures were fabricated on a cyclic olefin polymer (COP) substrate by the nanoimprint lithography method and reactive ion etching as explained in a previous work 17 . A scanning electron microscope (SEM) image for the fabricated structures is shown in the inset of Fig. 7(b).
The fluidic channel used in the experiments was fabricated using replica moulding method as explained in 16,24 . Cary 5000 spectrometer was used in the set-up, illustrated in Fig. 7(b), to measure the extinction curves for the nanocrescents while injecting the ethanol and streptavidin solutions.