The invention relates to the field of localized surface plasmon resonance (LSPR) spectrometry. More specifically, the invention relates to a system and method for determining the peak wavelength in a spectrum sensed by an LSPR spectrometer and novel method for reducing the noise of the peak wavelength signal.
A localized surface plasmon resonance (LSPR) spectrometer is a chemical analysis spectrometer in which ligand protein molecules are immobilized onto nanoparticles such as gold nanoparticles. The molecule to be analyzed, known as the analyte, binds to the ligand, causing a shift in LSPR resonant frequency of the nanoparticle. This resonant frequency is probed using absorbance/reflectance spectrometry, and is seen as a peak in the frequency/wavelength of the absorbance/reflectance. The peak wavelength signal as a function of time can be analyzed to evaluate the binding kinetics and other chemical parameters associated with the analyte and the ligand.
The invention provides systems and methods for finding the peak wavelength of the spectrum sensed by an LSPR spectrometer is described herein. The method comprises reading an image representing the reflected/absorbed spectrum, using a mathematical model of the LSPR spectrometer system to estimate a parametric curve representing the absorbance/reflectance spectrum, and adjusting or optimizing the parameters of the parametric curve so as to increase the likelihood of the parametric curve representing the sensed spectrum. Also described herein is a novel method to achieve LSPR peak wavelength signal noise reduction using an adaptive regularization algorithm.
1.1. Acronyms
“Absorbance spectrum” is the ratio of the spectrum of light absorbed by the LSPR gold nanoparticles to the uniform input light spectrum incident on them.
“Core-shell nanoparticles” means nanoparticles that consist of a core particle encapsulated by a shell.
“Localized surface plasmon resonance” means the collective oscillation of electrons at the interface of metallic structures.
“Nanoparticles” means particles with one or more dimensions less than 100 nm.
“Particles” means particles with one or more dimensions greater than 100 nm.
“Reflectance spectrum” means the ratio of the spectrum of light reflected from the LSPR gold nanoparticles to the uniform input light spectrum incident on them.
The accompanying drawings illustrate various embodiments and are a part of the specification. The illustrated embodiments are merely examples and do not limit the scope of the disclosure. Throughout the drawings, identical or similar reference numbers designate identical or similar elements.
The invention provides systems and methods for finding the peak wavelength of the spectrum sensed by an LSPR spectrometer is described herein. In one aspect, the method includes:
In an embodiment, the image processor 110 may be a dedicated hardware designed to perform the image processing task. In another embodiment, the image processor 110 may be a computer running a program which performs the computations for estimating the peak reflectance/absorbance wavelength of the LSPR sensor 104.
In an embodiment, the reflectance/absorbance spectrum of the LSPR sensor 104 may be modeled as a function of the model 204:
R=a+bƒ(c,d),
where a, b, c and d are parameters of the model and f is a known function involving the peak wavelength. The peak of R comes from the peak of ƒ(c,d) which encompasses the peak wavelength 202. The maximum value of R will depend on parameters a, b, c and d. In an embodiment, the parameters a and b represent an arbitrary base and scale respectively.
In a preferred embodiment, the function ƒ(c,d) is determined using a computer to calculate the Mie scattering theory applied to AuNPs. Mie scattering theory is a solution to the Maxwell's equations that describes the scattering of an electromagnetic wave by a homogeneous spherical particle having a refractive index (Ri) different from that of the medium surrounding it. Mie theory can also be applied to non-homogeneous core-shell spheres. In the LSPR case, the AuNP forms the core, while the ligand and bound analyte together form the shell. As the shell thickness or refractive index changes, Mie theory predicts a change in the scattering properties of the nanoparticles. It also predicts a shift in the wavelength peak. The parameters c and d represent shell thickness and shell refractive index or some functions of these quantities.
In another embodiment, the function ƒ is modeled as a log-normal function. A log-normal function is the probability distribution function of a random variable where the logarithm of that function is normally distributed. So, if Y is normally distributed, then X, such that Y=ln(X), is log-normally distributed. Function ln( ) represents the natural log. The mathematical expression for log normal function is as follows:
For log-normal function, it's maxima can be found at:
ωmax=e(μ−σ
Considering parameters c=μ and d=σ, the function ƒ can be re-written as
In an embodiment, the parameters of ƒ are defined as: c=eμ and d=eσ. The function ƒ can then be written as
Source spectrum 206 represents the spectrum provided by the illumination source 102 on the LSPR sensor 104. The reflectance/absorbance spectrum R (model 204) is multiplied with the source spectrum 206 to give a reflected/absorbed spectrum of light. Thus, reflected/absorbed spectrum of light is given by:
r=S(a+bƒ(c,d)),
where S represents the illumination source spectrum.
The effects of dispersive optics 106 are modeled by the diffraction model 208. In an embodiment, the diffraction model 208 is a convolution. In an embodiment, the convolution kernel is modeled as the area of a circle convolved with a 1D Gaussian function. The circle represents the exit surface of the fiber optics and the Gaussian represents the dispersion due to the optical components. Thus, the effect of the dispersive optics on the reflected/absorbed spectrum of light can be represented as a 2D convolution operator G:
m=G*S(a+bƒ(c,d)), (1)
where G represents the convolution kernel of the dispersive optics. In an embodiment, G is modeled as the convolution of a discrete 2D function of the area of a circle, which is 1 inside a circle and 0 outside the circle, with a Gaussian distribution.
The parameter m represents the modeled spectral image as opposed to the actual spectral image recorded by the imaging sensor 108. Due to physical limitations as well as the limitations of the imaging sensor 108, the model spectral image m, in equation (1) above, is perturbed by an imaging noise 210. The imaging noise 210 comprises one or more of photon noise, quantization noise, dark noise, additive white noise, etc.
In an embodiment, the imaging noise 210 comprises photon noise only. The number of photons falling on each pixel of the imaging sensor 108 is different in each exposure interval. This variation, also called the photon noise, typically follows a Poisson distribution. For a Poisson distributed photon noise and an expected pixel value of mi, the probability of observing pixel value xi in data is given by:
where k is the number of digital levels per photon—also known as analog gain of the imaging sensor 108.
The dispersed and noise perturbed reflected/absorbed light data is stored in the imaging sensor 108 as spectral image 212. The data stored in spectral image 212 are the pixel intensity values obtained from the image sensor 108 by the image processor 110. Even if values do not change in subsequent m (the modeled spectral image), the actual image pixel values will be different in each subsequent image owing to random photon noise. In the sections below we describe how the maximum likelihood formulation attempts to find the modeled spectral image, m, which increases, e.g., maximizes the likelihood of seeing the observed pixel intensity values of the actual image.
Step 402 is followed by a step 404 of constructing a maximum likelihood function for obtaining the read spectral image, given the initial parameter values found in step 402. In an embodiment, the photon noise is the most dominant sensor noise and the other forms of noise are ignored. In this case, m=G*S(a+bƒ(c,d)) is the expected value of the 2D spectral image. The actual pixel values would be perturbed by a Poisson distributed photon noise.
As the noise follows a scaled Poisson distribution, the probability of observing the particular pixel value xi, given expected value mi is given as follows:
where k=analog gain and i iterates over pixels.
Thus, the total probability of observing the sensor image, given expected value mi at pixel i is given by the product of probability for each pixel as shown in equation (2) below.
As the logarithm function is monotonic, taking the logarithm of a function preserves its maxima, i.e, argmaxxƒ(x)=argmaxxln(ƒ(x)).
Taking the natural logarithm on both sides of the equation (2) results in the following equation:
The last term in the summation does not change with any of the parameters a, b, c or d. Since the constant term will not affect the maxima, it can be omitted. Thus, the likelihood function for maximization would be equation (3) shown below.
In an embodiment, the likelihood formulation derived above is directly used as the objective function for parameter estimation. Optionally, the likelihood function could be formulated such that the estimated parameter values and hence peak wavelength values are more regular. Here regular implies that the values are not completely devoid of any structure. Depending on the amount of noise, estimated parameters and hence estimated peak wavelengths could vary quite a bit from one spectral image to another, with the variation showing up as noise in the wavelength peak signal. A structure could be defined which does not allow the peak wavelengths to vary significantly from one frame to the next. This reduces the noise in the temporal peak wavelength signal. Such a structure could be added to the likelihood formulation itself.
In an embodiment, the option is available to add a temporal regularization term to the likelihood term, in step 406. The Bayesian estimation technique allows the specification of a distribution on the expected pixel values m. According to the Bayesian formulation:
Taking the natural logarithm of both sides of the above equation produces the following equation:
ln(P(m|x))=ln(p(x|m))+ln(P(m))−ln(P(x)).
The last term does not depend on the parameters and can be ignored for optimization e.g., maximization. Thus, it can be written as:
argmaxm[ln(p(m|x))]=argmaxm[ln(p(x|m))+ln(p(m))] (4).
The first term on the RHS is equal to the LHS of equation (3). The second term specifies a structure on m and in turn a structure on peak wavelength signal. In an embodiment, the Wiener-Khinchin-Einstein theorem is used to formulate this structure.
The Wiener-Khinchin-Einstein theorem states: Given a wide-sense stationary random process M, there exists a causal filter g and a stationary random process W, such that M=W*g and W is white noise.
In an embodiment, the peak wavelength signal is assumed to change very gradually. It is expected that the high frequency content in the peak wavelength signal is very low. So, the suitable filter to be used would be a low pass filter. A pure autoregressive filter with a pole at 0 frequency has a suitable frequency response and also simplifies the math. Consider an autoregressive filter g with the following recurrence equation:
m[n]=w[n]−(c1m[n−1]+c2m[n−2] . . . ckm[n−k]) Thus, w[n]=m[n]+c1m[n−1]+c2m[n−2] . . . ckm[n−k]w[n]=Σicim[n−i]
In an embodiment, the white noise process, W, is an IID (independent and identically distributed) Gaussian white noise process, where in the gaussian distributed white noise process equation, −w2/n represents (x−μ)2 with μ=0. The probability distribution function would be expressed as:
Thus,
Thus,
ln(pM(m))=C1+C2Σ1(cim[n−i])2
Constant C1 does not affect the maxima and can be ignored. Substituting the above relationships in equations (3) and (4),
Here, index i iterates over spectral image pixels, while j iterates over the temporal spectral image sequence.
Thus, the objective function to be increased, e.g., maximized, is:
C2 represents the relative scale between the regularization term and the likelihood term.
In an embodiment, the value of C2 is fixed. In an embodiment, C2 is adapted based on the expected trend of the peak wavelength signal.
The objective function derived after step 404 or optionally, after step 406 having the maximum likelihood is obtained in step 408. Algorithms known in the art, such as gradient descent, conjugate gradient, newton method, etc. can be used for finding the maxima. The peak wavelength corresponding to the parameter values at maxima represents the peak wavelength for the spectral image.
The greater the value of C2, the smaller is the variation in the estimated peak wavelength signal. Thus, high C2 is desirable when no binding events are happening at the LSPR sensor and the variation in peak wavelength signal is purely stochastic. On the other hand, a low C2 is desirable when the estimated peak wavelength signal needs to track fast changes in the peak wavelength due to binding events at the LSPR sensor.
In another embodiment, the approximate value of peak wavelength is found by fitting an analytical function to the spectral image and then finding the maxima of that analytical function. In an embodiment, the analytical function is the Gaussian function. In another embodiment, the analytical function is a log-normal function.
In an embodiment, approximate peak wavelength is found by successively relaxing parameters a, b and d from equation (3). Successive relaxation involves modifying a parameter value for minimizing the sum of the squared error between the estimate and spectral image. In the implementation, the best fit for parameters a and b is obtained using the algorithm as described for
In an embodiment the approximate peak wavelength values are stored for future use.
After finding the approximate peak wavelength values for a lookahead number of spectral images, the approximate slope of the peak wavelength signal is found in step 504. The approximate slope is the slope of the best fit (least squares) line passing through the lookahead number of approximate peak wavelength values.
A step 506 computes the adaptive regularization scale by evaluating a particular function of the slope from step 504. In an embodiment the function is a binary function giving two distinct values of regularization scale across a range of values of slope. In an embodiment, the function is a discrete function giving n distinct values of regularization scale across a range of values of slope. In an embodiment, the function is a linear function from slope to regularization scale. In an embodiment, the function from slope to regularization scale is bounded by an experimentally calculated quadratic function, using a computer, defined as follows:
C2=Cmax2(1−25*max(slope2, 0.04)), where Cmax2 is the maximum value that C2 can take. This formula ensures that the maximum value of C2 is equal to Cmax2. In addition, this formula ensures that the C2 saturates at a slope equal to 0.2 and is zero for all higher slopes. The slope of 0.2 was calculated experimentally using a computer and determined to be a good saturation level. Finally, C2 depending quadratically on slope performed better than linear dependence.
By chain rule,
and finally, the gradient:
Method 600 further comprises a step 604 of taking second order partial derivatives of the objective function with respect to the parameters, to find the hessian matrix of the objective function. The chain rule can be applied to equations 6a-f to obtain second order derivatives. Finally, the hessian matrix can be written as follows:
Using the gradient and hessian 4×4 matrix, a step 606 is performed to compute the change in parameters using the Newton method as described below.
[δa,δb,δc,δd]=step*hessian−1(Plog)∇Plog
where step is a real number constant chosen so as to prevent the parameter values from diverging. In an embodiment, the constant step is equal to 1. A step 608 updates the parameter values as follows:
anew=aold+δa,bnew=bold+δb,Cnew=/Cold+δc and dnew=dold+δd.
The procedure for updating the parameter values is iterated until the stopping condition is reached. In an embodiment, the stopping condition is defined as follows:
where cstop is a predefined constant. In an embodiment, cstop=1×10−4.
Once the optimal parameters for the mathematical model of the LSPR spectrometer system have been obtained using method 600, a separate step is performed where the peak wavelength signal with reduced noise is analyzed by downstream software to gauge the reaction kinetics parameters for the ligand-analyte pair under test.
Various modifications and variations of the disclosed methods, compositions and uses of the invention will be apparent to the skilled person without departing from the scope and spirit of the invention. Although the invention has been disclosed in connection with specific preferred embodiments, it should be understood that the invention as claimed should not be unduly limited to such specific embodiments. Indeed, various modifications of the disclosed modes for carrying out the invention, which are obvious to the skilled person in view of this specification are intended to be within the scope of the following claims.
The present invention may be implemented using hardware, software, or a combination thereof and may be implemented in one or more computer systems or other processing systems. In one aspect, the invention is directed toward one or more computer systems capable of carrying out the functionality described herein.
Unless specifically stated otherwise, terms such as “receiving,” “routing,” “updating,” “providing,” or the like, refer to actions and processes performed or implemented by computing devices that manipulates and transforms data represented as physical (electronic) quantities within the computing device's registers and memories into other data similarly represented as physical quantities within the computing device memories or registers or other such information storage, transmission, or display devices. Also, the terms “first,” “second,” “third,” etc., as used herein are meant as labels to distinguish among different elements and may not necessarily have an ordinal meaning according to their numerical designation.
This application is a 371 application of International Application No. PCT/CA2022/050291, filed Mar. 1, 2022, which claims benefit of priority to U.S. Pat. App. Ser. No. 63/158,214, filed on Mar. 8, 2021, which are incorporated herein by reference.
Filing Document | Filing Date | Country | Kind |
---|---|---|---|
PCT/CA2022/050291 | 3/1/2022 | WO |
Publishing Document | Publishing Date | Country | Kind |
---|---|---|---|
WO2022/187931 | 9/15/2022 | WO | A |
Number | Name | Date | Kind |
---|---|---|---|
7217574 | Pien et al. | May 2007 | B2 |
10254216 | Sieben | Apr 2019 | B2 |
20090004757 | Yguerabide | Jan 2009 | A1 |
20130003070 | Sezaki | Jan 2013 | A1 |
20150247846 | Gerion et al. | Sep 2015 | A1 |
20170089832 | Uemura | Mar 2017 | A1 |
20170219488 | Fukuda | Aug 2017 | A1 |
Entry |
---|
Abumazwed, Enhancing data extraction from localized surface plasmon resonance biosensors. McGill University (Canada). PhD Thesis 182pgs. (2017). |
PCT/CA2022/050291 International Search Report and Written Opinion dated May 10, 2022. |
Number | Date | Country | |
---|---|---|---|
20240175749 A1 | May 2024 | US |
Number | Date | Country | |
---|---|---|---|
63158214 | Mar 2021 | US |