Advancing image quantification methods and tools for analysis of nanoparticle electrokinetics

Advancing image quantification methods and tools for analysis of nanoparticle electrokinetics D. J. Bakewell,1,a J. Bailey,2,3 and D. Holmes2,a 1Department of Electrical Engineering and Electronics, University of Liverpool, Brownlow Hill, Liverpool, L69 3GJ, UK 2London Centre for Nanotechnology, University College London, 17-19 Gordon Street, London, WC1H 0AH, UK 3CoMPLEX: Centre for Mathematics and Physics in the Life Sciences and Experimental Biology, University College London, Gower Street, London, WC1E 6BT, UK


I. INTRODUCTION
Nanoparticle technologies have been developing rapidly over the last decade with many practical applications in nano-medicine (e.g.diagnostics and drug delivery).As a result of this growth, there is a particular interest in the development of methods and tools for the rapid and accurate characterisation of the electrical properties of biological particles suspended in aqueous solution.The advent of microfluidic and lab-on-a-chip devices has enabled the capability of manipulating small volumes of liquid sample; and the inclusion of integrate sensors within these devices allows for rapid and controlled measurements to be performed on sub-microlitre volumes.In this paper we describe an AC electrokinetic method, with associated imaging and data processing, capable of measuring the dielectric properties of nanoparticles.
Dielectrophoresis (DEP) is an electrokinetic technique that enables manipulation of microscopic and nanosized objects suspended in solution.][3] DEP is the translational movement of an electrically polarisable body by the action of a non-uniform electric field.Typically, this is implemented by applying a low voltage radio frequency (RF) electrical signal to microfabricated electrodes immersed in a liquid (typically of low conductivity).Many investigations of DEP, particularly for miniaturised applications, involve video recording the motion of biological particles and colloids using light microscopy.Indeed fluorescence microscopy has become a standard tool in many electrokinetic laboratories for real-time imaging of fluorescently labelled nanosized particles.Subsequent to recording the experiment, video images of microscopic particle dynamics are typically processed and quantified on a frame-by-frame basis, often by bespoke software developed by investigators and their groups.An important class of imaging is quantifying the motion of particles onto a planar object, such as, an electrode array whereby the entire array lies in focus for the duration of an experiment.Planar arrays often involve symmetric and periodic design features and widely used examples are interdigitated and castellated geometries. 1,3 is paper describes the significant advancement of image processing methods for DEP with examples of experiments using 200 nm diameter fluorescent nanoparticles manipulated using planar castellated arrays.A schematic of a typical DEP experimental setup for manipulating nanoparticles is shown in Fig. 1(a) and 1(b).The methods presented in this paper concern both (1) the actual process of image quantification of DEP collections, and (2) providing a robust foundation for investigating new imaging processing techniques that inform about the nanoparticle electrokinetics.The two imaging processing and quantification techniques described in this work represent novel approaches to quantifying data from DEP experiments and considerably advance the literature.As illustrated in Fig. 1(c), they entail (1) a geometric method that is based on two-dimensional (2D) spatial filtering predicated on the electrode design, and (2) a statistical filtering method that selects image pixels according to their levels of time-dependent variation in fluorescence.5][6] This process involves the repeated application of pDEP for durations, ranging from milliseconds to seconds.The applied All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported license.See: http://creativecommons.org/licenses/by/3.0/Downloaded to IP: 144.82.107.91 On: Wed, 27 May 2015 13:35 :27  signal can be varied in terms of signal amplitude (i.e. by modulating the carrier envelope) and frequency.
Applications of pDEP using planar electrode arrays include measurements of collection rates for characterising cells and their constituents, e.g.][6] Consequently, DEP video recorded fluorescence microscope dielectrophoretic spectroscopy is becoming an important tool for laboratory investigations of AC electrokinetic nanoparticle behaviour and could also advance novel investigations, e.g.][23][24] Although bespoke image processing methods have been reported for some DEP applications, e.g. cells suspended in circular apertures 25 and single feature structures, [26][27][28] the literature is sparse in terms DEP imaging involving spatially periodic structures and with capability for high resolution quantification.Consequently, there is considerable need for development of unsupervised, rapid image processing and quantification methods for DEP collections on periodic features and a theoretical foundation to support the techniques.This paper attempts to remedy this deficiency by advancing and demonstrating new, first principle, mathematical and statistical models for image processing.The paper briefly reviews DEP and important statistical concepts and parameters relevant to the work.A geometric method, suitable for castellated electrode planar array design, is described that achieves quantification of the time dependent DEP nanoparticle collections by spatial averaging each frame of a video.Two statistic-based methods, that do not require any information about the electrode array geometry, are introduced with examples that consider the application of the quantile and standard deviation functions for quantifying video images.Quantitative comparisons are made between the two methods of DEP collection experiments using 200 nm latex particles.

II. IMAGE PROCESSING METHODS
Quantifying DEP collections requires video filtering to achieve a reasonable signal-to-noise ratio (SNR).The two types of filtering, spatial geometric and statistic-based quantile and standard deviation, are detailed after a brief review of basic concepts and parameters.

A. DEP nanoparticle collection and release
The system of on-off switching of pDEP has been described in previous experimental and theoretical work. 4,5,14,29,30 A crtoon showing the repeated movement of nanoparticles attracted by pDEP towards, and away from, a planar array is shown in Supplementary Information, Fig. S.1(a) and S.1(b). 36The system essentially alternates between two states, 'off' and 'on'.Before switching on the DEP force at the start of the collection phase at t = t c , Fig The process is repeated; repetition can be regular (cyclic) or irregular (acyclic).

B. Space-time ordinates and probabilistic framework
A probabilistic description of nanoparticle movement is underpinned by a number of random processes, including the fact that the stochastic effect of Brownian thermal motion is significant for nanosized particles. 31The distribution of nanoparticles over space and time is described by the nanoparticle concentration, c(x, t) where x = (x, y, z).The concentration equals the product of the probability density function (pdf), p(x, t), and the total number of nanoparticles that remains constant throughout the experiment, i.e. nanoparticles do not enter or escape.Thus, a simple correspondence exists between the experimental and intuitive nanoparticle concentration, and the theoretically useful pdf.The imaging process involves optical fluorescence, emitted from nanoparticles located in 3D space at a moment in time, being detected on a CCD, 2D pixel array.
In terms of a probabilistic framework, the random variable (RV) in this context, is a process whereby a pixel attribute, a monochromatic intensity, is assigned a greyscale integer value, i g , (subscript 'g' denotes greyscale).The RV is a function, I g , that is dependent on outcomes, ω, of simultaneous light detection by each pixel element in the array. 32The pixel values range from black to white, i g ∈ {0, 1, . . ., i gx }, with 8-bit greyscale maximum value, i gx = 255.The pixel positions on the array are indexed according to a video-frame image matrix with matrix row, r, and column, c, so that the set of detection outcomes is abbreviated, ω = {ω 11 , ω 12 , . . .ω rc } where {r, c : 1 ≤ r ≤ r mx , 1 ≤ c ≤ c mx } with r mx and c mx , being the respective maxima, e.g.r mx = 512, c mx = 640.The sample space is the entire video, V , that can be partitioned into a set of frames, . ., F kx } with each frame associated with a discrete time, t = t 1 , t 2 , .., t k x .Thus, for a given frame sample space, F , the measure of probability, P, that the RV is a discrete greyscale value can be written in terms of a set of pixel detection outcomes, P{Ig(ω) = i g }.In the following sections, for convenience, the outcome of the RV is understood to be implicit, I g ≡ I g ( r, c ω ) and the set notation is simplified; 33,34 thus, P(Ig = i g ) reads as the probability of the RV, I g , being a particular discrete greyscale value, i g .
An example of a frame image segment is shown in Fig. 2(a).The image intensity is related to the optical fluorescence intensity, I (x, t), so that the RV at discrete time, t = t k , can be represented, where y is the vertical depth-of-focus, the coordinate basis is the same as adopted previously 4,5,11,12,14 and the pixilation is denoted by the subscript 'discrete'.The 2D CCD camera pixel array is parallel to the DEP electrode array so there are simple scalar relations between the pixel indices and continuous space coordinates, (r, c) and (x, −z), respectively.The detected light enables quantification of the nanoparticle concentration, e.g. for sufficiently dilute suspensions, I (x, t) ∝ c(x, t).There are also a number of phenomena that introduce randomizing processes, e.g.scattered light (optical noise), thermal noise from the camera, etc.For example, a pixel of an image of an electrode that ideally should be black with value, i g = 0, due to noise, can have values representing a range of grey shades, i.e. {i g : 0 ≤ i g ≤ i gx }.
Normalising a frequency histogram, f (i g ), of pixel values from a sample space, , with n pixels yields the relative frequency, f r (i g ), (subscript 'r' denotes relative).The grouping, according to greyscale value, is equivalent to the discrete probability distribution, where the relevant sample space is implicitly understood, P(Ig = i g ) = P(Ig = i g | ) e.g.= V , for the video.The expected value, or mean, of the RV, E{I g }, is given by the first moment, In Eq. (3), Eq. ( 2) is applied to the greyscale grouping on the left hand side (lhs) to regroup pixels according to a suitable independent variable (e.g.physical space and time parameters, r, c, t k ) to yield the right hand side (rhs); the rhs is useful for the geometric method.

C. Spatial geometric filtering
The pDEP nanoparticle collections occur near the castellated electrode edges so to improve the SNR pixels near and across the edges are selectively filtered.In previous work using an interdigitated electrode design, the coordinate location of edges was entered manually using a screen cursor. 11,12,14 Stial filtering requires a 2D grid that points to the locations of electrode edges and is determined using a template image.The 2D matrix template frame image is transformed into two, 1D vectors and feature recognition methods are applied to extract grid coordinates.We found this method for locating electrode edges to be robust against noise.

Template image
The template frame, I g (r, c), is found by time averaging a sequence of frames with beads undergoing Brownian motion (i.e.DEP switched off).Using Eq. ( 1) and Eq.(3), where k m need be only a few frames to be effective, i.e. 3 ≤ k m < k x .The averaging effect tends to smear any motion of the beads and reinforces the pattern of the stationary array, thus, improving the SNR for edge detection.Edge detection, alone, often incurs two types of errors, i.e. detecting an edge when none exists, and not detecting an edge that actually exists.Instead, error-free edge detection, which is a local property, is enabled by preserving information about the global periodic and symmetric properties of the castellated design.Global design properties, in each of the longitudinal and transverse directions, are preserved by taking a 1D projection of the 2D template matrix in each of the two orthogonal directions.The simplest projection that preserves the length of the vector, is the mean of pixel values of rows and columns.The 2D template matrix is, therefore, transformed into two, 1D vectors.The periodic and symmetric properties of the annotated part-frame castellated design, shown in Fig. 2(b), means that 2D electrode edge positions can be represented by a number of parameters, as labelled.In the longitudinal direction (denoted by 'l' subscript), o l is the longitudinal pixel offset and w l is the average number of pixels to specify the width between electrode edges.
In the transverse direction (denoted by the first 't' subscript), o t is the offset, w tc is the width between the bottom of the castellated bay and tip, g t is the inter-electrode gap that includes tips of the edges, and w te is the width between the bay edges across the electrode surface (the second subscripts 'c' and 'e' denote as italicised).The parameters sum to the average number of pixels, in the transverse direction, per cell pair, w tp = 2(2w tc + g t + w te ), as shown in the figure.The offset pixels are a set of positive real integers, {o l , o t } ∈ J + , and the set of widths, specifying the average number of pixels, are positive real numbers, {w l , w tc , g t , w te , w tp } ∈ R + .

Transverse mean
Using Eq. ( 3), the mean of transverse (row) pixel values of the 2D template matrix given by Eq. ( 4), I g (r, c), is the longitudinal 1D vector quantity, where the over-bar denotes a 1D spatial mean and the greyscale value is implicit with subscript 'g' being replaced by 'l' that denotes 'longitudinal' direction.Referring to the annotated template partframe, Fig. 2(b), the associated plot of Īl (c) is shown in Fig. 3(a).Optical effects of the fluorescence microscope can locally bias Īl (c) in certain locations, making some areas bright and others dull.Therefore, to reduce unwanted biasing, our interest lies in the short-scale greyscale variation, v l , over tens of pixels typifying the pattern expected from a castellated design, where b l (c ) is the baseline, or large-scale, spline-fit and it is set below min{ν l (c)} so that the variation is always positive, v l (c) > 0 ∀ c.The design of the castellated array means that the edges coincide with the periodic minima, min{v t (c)}.A robust method of finding the edge positions, is to re-parameterise the minimisation problem into 2D: o l and w l .Values for the offset and electrode width parameters are estimated by minimising a periodic test function, where [ ] I denotes the integer function and k ml is selected, such that, [w l k ml + o l ] I ≤ c mx .Optimum estimates for the longitudinal offset, ôl and width, ŵl ('ˆ' denotes estimate) occur at the minimum of the test function, ṽl , where the '+' subscript denotes the column edge ordinates shown by the grid locations superposed over the array in Fig. 2(a) and 2(b).In this particular example the bias is evident but not detrimental, so that it is sufficient, b l (c ) = 0.

Longitudinal mean
Using Eq. ( 3), the mean of longitudinal (column) pixel values of the 2D template matrix given by Eq. ( 4), I g (r, c), is given by the transverse 1D vector quantity, where 't' denotes 'transverse' direction.To accommodate optical biasing effects in the transverse direction, the variation of pixel values is, where b t (r ) is the large-scale baseline spline fit.The design of the castellated array means that four parameters are required to specify the edges.Referring to the annotated template part-frame, Fig. 2(b), the associated plot of Īt (r ) is shown in Fig. 3(b).The pixel parameters in Fig. 3(b) (with the first subscript 't' denoting transverse direction) are o t , and widths w tc , g t , and w te .Since the parameters sum to w tp , only three out of the four parameters need estimation.The procedure for finding the edges to the 'bottle' shape peak features in Fig. 3(b) involves fitting four parameters and requires a different approach compared with transverse mean.It entails fitting four line segments to short, straight sections, of the 'neck' and 'shoulder' profiles, in each of the 'bottle' shape peak features as labelled in the Fig. 3(b), centred on initial estimates.The details are given in Appendix A. The transverse row ordinates, denoted by '+' subscript, are calculated from the estimates (denoted by 'ˆ'), where u( j − a) = 1, j > a 0, j ≤ a is the unit step function, k = {1, 2, . . ., 2n tp } is the 'bottle peak' feature index, j = {1, 2, 3, 4 } is the fitted line segment number for each of the features, n tp is the number of transverse cell pairs (typically n tp = 2), and [ ] I is the integer function.It is understood the segment number cycles through all four elements, j = {1, . . ., 4 } for each increment of the feature index, k.The placement of the row edge ordinates is illustrated by the '+' grid superimposed over the array image in Fig. 2(a) and 2(b).In Fig. 3(b) the far right marker is located at r In terms of a statistical framework, the rectangular '+' grid region enables partitioning of the 2D frame sample into physically separate regions, F = { S , X }: the pixels within the rectangular grid constitutes the sample space for geometric processing, S , whereas pixels out-with the grid constitute the excluded space, X , that are not used in the geometric approach, as depicted in Fig. 2(a).Using the same parameter values as the examples in Fig. 2 and 3, S is 22 unit cells by 2 pairs, or 22w l × 2w tp .

Characteristic intensity
Applying the template grid to frames in the video enables harnessing the symmetric and periodic properties of the castellated array to generate a 2D greyscale surface of an array cell that characterises the light intensity of the frame region of interest; typically representative of the entire frame.Advancing the much simpler procedure for interdigitated electrodes, 12 castellated array cells exhibit small variations in the number of pixels both in their longitudinal and transverse widths and quantization errors result.To circumvent this problem, the frames are interpolated and re-sampled with slightly more pixels than the maximum number of pixels required to ensure all the frame matrixes have the same size.The 2D interpolation and re-sampling process is expressed, reconstructed intensity (13)   where I c denotes pixel intensities for a cell pair and subscript 's' denotes 're-sample'.A typical cell pair is shown in Fig. 2(b).Cell pairs have the same spatial orientation are added and divided by the total number, and the process is known as periodic averaging. 12Periodic averaging is expressed in Eq. ( 14) where the re-sample subscript is understood and dropped.The cells exhibit symmetry about either the longitudinal or transverse axis, so that two unique spatial orientations result from periodic averaging.Since n tp = 2, in this example, four cells result from periodic averaging and are shown in Fig. S.3. 36They are symmetric so they can be further averaged by rotating about either the longitudinal or transverse axis (or rotating 180 • ), so that they superpose, and thus, be averaged.The result is called symmetric averaging and yields a single characteristic cell shown in Fig. 4(a).The resultant light intensity of a unit cell, shown by the contours in Fig. 4(a), characterises the optical intensity for the entire electrode array -except the unimportant electrode interior regions that have been excluded for clarity.The methods and algorithms used for 2D are an advancement of the 1D relations described in Ref. 12, hence a summary is given.The example contour plot is shown in Fig. 4(a) has been superposed with coordinates of the electrode castellation.Selecting a spatial rectangular Region of Interest (RoI), the greyscale intensity of each pixel is averaged over selected cell rows and columns resulting in a greyscale number for the frame.The process is summarised with time made explicit, so that a sequence of re-sampled frames, is progressively averaged, starting with the result of Eq. (13) Thus, the 2D fluorescence region of interest on each video frame becomes converted into a number.Consequently, a sequence of 2D frames yields the time dependent intensity, with double bars overhead, denoting 2D spatial mean, I gm (t), with time made continuous t k → t in the last step.In Eq. ( 14) the subscripts 'gm' denotes geometric-based method and it is understand that it is evaluated for the gap region.
An example of time dependent intensity generated by the geometric method is shown in Fig. 4(b).The characteristic cell shown in Fig. 4(a) is in fact a time average of three frames, as encircled on the plot before the first DEP collection occurs.Time averaging makes the contours smooth and ideal for template parameter estimation, as described earlier.An example of contours that are spatially integrated to yield the second DEP collection peak time-point is shown by the right RoI in Fig. 4(b).The contours clearly indicate the accumulation of fluorescent nanoparticles near the electrode tips that is expected of pDEP.

D. Spatial statistical filtering
Statistical filtering of each image frame in a video sequence does not require explicit information about the DEP electrode geometry.Clearly, avoiding feature extraction about the 2D electrode array pattern is an advantage compared with the geometric method.

Frame statistics and sample space partitioning
Referring to the probabilistic framework introduced in section II.B, in addition to the frequency and discrete probability, the action of pDEP collection and release is also conveniently shown using the discrete relative cumulative sum, where the number of pixels in the frame is n p = r mx c mx , (e.g. using values given in section II.B, n p = 327680) and other parameters have been previously defined.
Within the geometric sample space, the 2D planar electrode and inter-electrode gap is further partitioned into sample subsets, written as, S = { G , E } and E = { C , I } (subscripts 'G' and 'E' denote 'gap' and 'electrode'; and 'C' and 'I' denote 'castellated' and 'interior' areas of the electrodes, respectively.Referring to Fig. 2 and Fig. 3, for example, the pixel sample 2D areas for the unit cell pair detailed in Fig. 2(b), are: G is w l × 2g t , C is w l × 4w tc , and I is w l × 2w te .The discrete distributions are related to each other using conditional notation, 32 where P{Ig = i g | G } reads as the probability conditioned on a given event occurring, G , and ¯ G is the complementary area of F that includes { C , I , X }.Letting P( G ) = γ where γ is numerically equivalent to the proportion of frame pixels imaging the gap, n g , then axiomatically, P( ¯ G ) = 1 − γ .Using the same data as for the examples described in Fig. 2 to 5, γ = n g /n p = 35840/327680 ∼ = 0.11.The discrete distribution is sufficiently resolved and smooth that we consider the interpolating between discrete data values, i k , so that they become continuous, i v , and the same applies for time, t k → t.This means that the relative frequency is easily converted into a probability density function (pdf), f r (i g , t k ) → p(i v , t) and the cumulative relative frequency becomes the cumulative density function (cdf), F r (i k , t k ) → C(i v , t).Two examples of pdfs and cdfs are shown in Fig. 5 corresponding to nanoparticle fluorescence before (t = 0.3 s) and after (t = 47.2 s) DEP collection shown in Fig. 4(b).Using (16) it can be demonstrated that for both time points, the left mode of the pdfs arises from the large number of pixels that constitute the dark electrode areas, whereas the right mode indicates fluorescence over the glass, as labelled.Thus, it is the right mode that changes during DEP collections, whereas the left exhibits spiky variation attributed solely to noise.The cdfs corresponding to the pdfs are illustrated in Fig. 5 inset and it is clear that DEP collections result in changes in the upper tails of the cdfs as they approach unity, i.e. as C(i v , t) ∼ = 1, as shown.

Video statistical filtering
Our interest lies in variations of the intensity distributions over the duration of the video (i.e. the experiment) without any information about the electrode geometry being needed.Of particular interest is finding a good phenomenological approximation for the gap mean.An example of frame pdf and cdf surfaces, p(i v , t) and C(i v , t), using the frame sequence that generated Fig. 4 36 The surfaces show the action of pulsed pDEP.This means that by selecting suitable value for the greyscale, i v , as shown, the cdf can yield a pulsed pDEP signal, C(t).However, controlling the greyscale intensity is not desirable and there are two more effective ways of finding a pulsed pDEP signal that approximates the gap mean.
a. Approximating the gap mean by the quantile function.The first method is motivated by the differences between the cdfs in Fig. 5 being pronounced for cdf values near unity, i.e. selecting intensity values, with ordered greyscale values upwards from zero, for almost all of the pixel population.This is straightforwardly achieved by selecting a chosen quantile value, q v , of the cdf, q v = C(i v ), e.g.q v = 99%.The inverse cdf of a distribution is the quantile function, Q, 35 Q(x) = C −1 (x) and it is supported in standard software, e.g.Mathematica TM 8 (Wolfram Research, Campaign, IL) and Matlab TM 7.14 (Mathworks Inc., Natick, MA, USA).By definition, the Q function FIG. 6. Frame quantile function (inverse of cumulative sum) versus time (frame no.) and listed quantile constants, Q F (q ν , t), showing pulsed DEP collection and release for the same video sequence as Fig. 4(b).Quantile plots for constant values from 0.98 to 0.995 are the most satisfactory, showing patterns consistent with on-off DEP collections and release, whereas lower values are too noisy.In terms of frame statistics, this means that about 99% of frame pixels, ordered from lowest to highest, account for DEP behavior.A percentage too low resolves poorly, and very close to 100% they saturate.
for a specific cdf, operating for a particular quantile value, maps to the greyscale value, and we can also write Q [C(i v , t)] = i v (t).A typical frame, for example, for an inverse cdf at a particular moment in time and using a quantile value q v = 0.95, yields Q(0.95) =i v = 100.Our interest in the use of the quantile function arises from the property shown in Appendix B, that for two normal distributions with means μ and μ , under certain conditions relating to their standard deviations (stds), an estimate of their difference, μ, can be approximated by applying the quantile operation to the cdfs, The conditions include, e.g. that the means and stds are approximately linearly correlated, or indeed the stds are time invariant.The estimation is a good approximation for normal pdfs that are truncated the lower tail -as is relevant to our DEP collection experiments.Since the incremental change in mean occurs on a frame-by-frame basis (with a constant frame rate), then Eq. ( 18) can be expressed with respect to time as, where the subscript 'F' denotes that the quantile function is applied to each, entire, frame -thus avoiding any need for inputting geometric information, e.g.feature extraction.Pixel greyscale intensity values for the frame sequence that generated Fig.  36 are illustrated in Fig. 6 for {q v : 0.900 ≤ q v ≤ 0.995}.The intensity profiles show that for very high quantile q v values selected, Q F (q v , t) is similar to I gm (t) in Fig. 4(b).In the following sections initial collection rate ratio parameters are experimentally evaluated and circumvent the need to concern ourselves with scale and absolute level; linear proportionality is therefore sufficient for the application of Eq. ( 19) in the next section.Consequently, the statistic-based quantile method can be directly compared with the geometric method for quantifying DEP collections.
b. Approximating the gap mean by the frame standard deviation.The pdf of the frame pixel values in Fig. 5 is bimodal with tail of the smaller right mode being of interest.The features of the distribution have a number of consequences: (i) the frame mean is noisy due to the large proportion of frame pixels representing the electrodes (left mode), (ii) the mean and std of the gap are positively, approximately linearly correlated, and (iii) the std of the left mode is small so that including its component in the total frame std appears to add very little noise.Consequently the pulsed DEP activity tail of the right mode can be indicated using the std and not the mean of the frame pixels.
The spatial variance of the frame, V ar{I g }, can be evaluated by finding the second moment, using a similar method described for the first moment in Eq. ( 3).Applying Eq. ( 20) and Eq. ( 3), for a large number of pixels, the asymptotically unbiased estimator 34 of the variance is Hence, the time varying spatial standard deviation for the entire frame, (hence not requiring any geometric information to be inputted), denoted by subscript 'F', is given by, As stated previously, the application of Eq. ( 22) in the next section avoids the need to concern ourselves with scale and level, so that σ F (t) can be compared with I gm (t).

III. EXPERIMENTAL COMPARISON OF GEOMETRIC AND STATISTICAL FILTERING
A dielectrophoretic DEP initial collection rate ratio method, using dual-cycle for characterising the dielectric properties of 200 nm diameter fluorescent latex nanospheres by estimating nanoparticle conductivity, was demonstrated by our group recently. 6In this paper we use the same method in combination with automated image analysis.

A. Laboratory setup
An arbitrary function generator (Tektronix AFG 3022B, OR, USA) was used to provide 2 V peak-to-peak, variable duty cycle, square wave enveloped, sinusoidal signals over a carrier frequency range of 10 kHz-10 MHz to the microelectrodes.Pulse duration, amplitude and applied frequencies were controlled by custom software written in LabVIEW TM 2011 (National Instruments Corp., Austin, TX, USA).
Castellated geometry interdigitated microelectrodes with critical feature sizes of 5 microns were fabricated using standard photolithography and lift-off techniques.A 100 nm thick layer of platinum was lithographically patterned on 500 μm thick Pyrex wafers.Individual devices were cut from the wafer.The devices were mounted on Veroboard and wire bonded to the copper strips on the Veroboard allowing robust electrical connection to the signal generator.A 3 mm diameter sample reservoir was fabricated in poly-dimethylsiloxane (PDMS).The PDMS and glass chip were exposed to an O 2 plasma to facilitate bonding of the PDMS to the glass/platinum substrate.The completed device was then sealed with a cover-slip to prevent sample evaporation.
Ultra-pure water having a resistance of 18.2 M cm (Purelab ultra, Elga process water, Buckinghamshire, UK) was used to prepare KCl (Sigma-Aldrich R Inc., ST Louis, MO) buffer solutions with conductivity of 2 × 10 three times in KCl buffer (2 × 10 −4 S/m) and suspended in the same buffer at a concentration of 5 × 10 10 spheres/ml (diluted 1:100 from 2% w/v stock solution).The concentration and monodispersity of the nanospheres was verified using nanoparticle tracking analysis (NanoSight LM10, Wiltshire, UK).
The motion of the nanospheres was observed using an inverted microscope (Axiovert 200 M, Zeiss, Germany) with epi-fluorescent illumination (HBO100, Zeiss), imaged with x40 objective and recorded with a digital camera (Thorlabs USB 2.0, Newton, NJ) at 10 frames/s.Videos were analysed using software written in Mathematica TM 8 (Wolfram Research, Campaign, IL) and Matlab TM 7.14 (Mathworks Inc., Natick, MA).Nanosphere collections at the high-field gradient regions of the castellated electrodes (i.e.areas where particles collect under the influence of pDEP) were quantified by measuring the fluorescence intensity at the electrode tips.

B. Experimental results and analysis
In each set of collection experiments, approximately 50μl aliquots of nanosphere suspension was pipetted into the PDMS reservoir on top of the device and sealed with a cover slip.The continuously pulsed (or AM) dual-cycle DEP frequency collections used an 'on' control pulse set at f = 0.7 MHz and the probe at f = 3.0 MHz, with equal dual-cycle periods, T = 17 s comprising 7 s on, 10 s off.The ground-to-peak voltage, during the DEP 'on' phase, was 1.0 V for all experiments.Each experiment recorded approximately five dual-cycles in 170 s with frame rate 10 frames/s.Collections for the first cycle were transient whereas the remaining cycles exhibited cyclic steady state, concurring with theoretical predictions. 4,6 hus, the choice of parameters values achieved (1) cyclic steady state and (2) an adequate fluorescence amplitude or nanoparticle fluctuation suitable for quantification.The geometric-based and statistic-based quantile and standard deviation fluorescence intensities, I gm (t), Q F (q v , t) with q v = 0.99 and σ F (t), respectively, are plotted using greyscale units of intensity in Fig. 7.The plots demonstrate the high level of correlation between the methods; the data from the quantile method has been offset vertically by about 11 units to allow for clearer visualization.
The plots comprise five dual-cycles, each comprising a control and probe cycle, as labelled in Fig. 7, where each cycle itself consists of a DEP collection (DEP switched on) followed by release (DEP switched off).The use of the control cycle mitigates against variations in the initial nanoparticle concentration, so that ratios of initial collection rates can be used to estimate relative measures of the frequency dependent DEP force, and consequently infer the nanoparticle conductivity DEP collection experiments. 6The mean ratio of the ith dual-cycle probe and control initial collection rates is given by, where t c i is the initial time point for the control and probe collection rates, e.g. at the start of the control cycle of the second dual-cycle in Fig. 7, t c i ∼ = 40 s.The first dual-cycle in Fig. 7 is transient whereas the remaining four indicate dual-cyclic steady state.Therefore, the first dual-cycle is not analysed and only the collection rates ratios of the following n = 4 dual-cycles are analysed.Collection rates for the probe and control cycles were each estimated by linear fitting to half the collection phase points (35 five data values).Hence, the example plotted in Fig. 7, yielded mean ratios for the quantile and std statistic and geometric methods, ρ = 0.45 ρ = 0.48 and ρ = 0.52, respectively.Therefore the three different approaches concurred to within 15% of each other.The geometric-based gap mean and statistic-based frame standard deviation intensities, I gm (t) and σ F (t) plotted in Fig. 7 practically overlay each other and data values used in the initial collection rate are plotted for the four dual (control and probe) cycles at cyclic steady state in Fig. 8.The relationship between frame standard deviation and gap mean data points, shown by the fitted scatterplot, shows two linear relations, with the control and probe points separately plotted, that are practically co-linear.Scatterplots for all data points for the gap mean and frame standard deviation also exhibited a very good co-linear relationship (data not shown).Thus, under the caveat that scale FIG. 7. Plot of the inter-electrode gap mean using geometric method, I gm (t), and frame standard deviation and frame quantile (upper trace; using q v = 0.99 value) using the statistical method, Q F (q ν , t)with q ν = 0.99 and σ F (t), respectively.All curves have been re-scaled vertically for illustration; the quantile has been also shifted vertically (by about 11 units) with respect to the other two traces to allow for better visualization of the curves.The geometric gap mean and statistical frame standard deviation are shown to almost overlay each other and are highly correlated for both the DEP collection and release phases.The statistical frame quantile is strongly correlated with the other two methods for the collection phase.Video frame rate was 10 frames/s.and level are not important quantification parameters, the behaviour of the statistic-based frame standard deviation can be used to approximate the gap mean.

IV. DISCUSSION
The geometric method achieves spatial filtering of the array fluorescence by recognising electrode edge locations by reducing the dimensions from 2D to two 1D.The techniques, implemented by simply taking pixel means in each of the orthogonal directions, combine local edge features with global periodic and symmetry properties inherent in the design of the castellated planar array.The fusion of these two properties significantly improves the feature recognition edge detection SNR so that the detection of electrode edge features is demonstrated to be error free.Specialised methods, described in section II and Appendix A, are then applied to find parameter values in each of the orthogonal directions, e.g.2D array region offsets, castellation electrode widths, inter-electrode gap, etc.The geometric filter parameters are inputted by algorithms using almost fully-automated software implemented in Matlab and from a template that is a time average of a few frames, or more.Once the parameter values for the edge positions are inputted, the software automatically evaluates the entire video and outputs the time dependent intensity, I gm (t).
Importantly, the frames used for the template are part of a video sequence briefly recorded just before the first DEP collection and release pulse of an experiment.That is, the template is created without the need to adjust the experiment or microscope working conditions.Some other trials have also successfully used the later stage of release phase for template making.This is an important benefit because adjustments can compromise parameter values due to small changes in array focus, mechanical vibrations, etc.To explain further; applying broad spectrum back-lighting for the fluorescence microscope, for example, gave more ideal longitudinal and transverse profiles with less noise, than exhibited in Fig. 3(a) and 3(b) (data not shown).However, backlighting while helpful, is not at all necessary for successful parameter estimation.Fast Fourier Transforms (FFTs) were also applied to frame images, and to the reduced dimension, 1D, vectors.However, the multiple frequencies resulting from the FFTs proved problematic; instead the specialised methods described in section II exhibited better robustness and viability.The global periodic and symmetry properties of the array were exploited to improve SNR, although, it could be used for aperiodic and/or asymmetric if appropriate electrode spatial design rules were available.
The geometric method for evaluating the filtered spatial 2D mean, on a frame-by-frame basis, essentially performs three roles.First and primary, as image processing and quantification for DEP induced nanoparticle manipulation.It is important to note that it is not necessary for the DEP to be positive, or be pulsed, since the template, which is fundamental for automated processing, uses the fluorescence from a few frames without the DEP switched on.The random, Brownian, motion of the beads, smeared by time averaging, provides enough optical scattered light for the electrode edges with periodic features, to be distinguished well above noise.The second role the geometric spatial method performs is, more subtly, a basis for comparing other methods -including the quantile and standard deviation statistic-based approaches that are applied to each, entire, frame of the video and, hence, avoid the need to input geometric information about the array.The third role lies in using the geometric method, in synergic conjunction with statistical methods, for exploring nanoparticle electrokinetics and other fluid dynamics.
One of the costs of the geometric method is the amount of time and effort required to setup geometric information, either by developing feature recognition as described, and/or manually by hand.As a consequence, we investigated statistic-based methods.The statistic-based quantile method, ideally, assumes normally distributed light intensity to estimate the spatial mean for each frame -as shown in Appendix B. However, a normal distribution truncated in the lower tail for the 'gap' (i.e.not the electrode region), as shown in Fig. 5, has also satisfactorily produced the time dependent intensity.Ideally, the estimation of the spatial mean could be achieved by any quantile value from above zero to unity.Instead, quantile values of limited range 0.98 to 0.995 produced a satisfactory time dependent intensity -and consistently across a range of DEP experiments.Lower quantile values tended to be confounded with noise.
The relationship between the geometric and statistic-based approaches is illustrated by Fig. 4 and 5.In Fig. 4(b), the fluorescence intensity is represented by the contours of the two gap rectangular RoIs.They show nanoparticle distributions before and after DEP collections coinciding with increases in spatial mean and spread (or standard deviation) as shown by the associated pdfs and cdfs in Fig. 5.The approximate linear correlation between gap spatial mean and frame measures of spread (quantile and standard deviation) is attributed to specific regional nanoparticle pDEP collection at the electrode tips.It is important to realise that there can be linear correlation between the mean and standard deviation, particularly for example, when these statistics include evaluating a bimodal distribution.Simple models show that a linear correlation between mean and standard deviation, for example, occurs when there is a relatively large pixel population that does not change with pDEP, i.e. left mode of pdfs in Fig. 5 representing dark electrodes, and small (γ ∼ = 11%) proportion of pixels that brighten considerably when pDEP is switched on, i.e. right mode of pdfs in Fig. 5 representing gaps.The results give support for the viability of the statistical approaches, thereby providing an exciting opportunity for image quantification that does not require any input (e.g.feature recognition) about electrode array edge locations.
Dual-cycle pulsed pDEP collection, a technique recently developed for determining the dielectric properties of nanoparticles, [4][5][6] was used as a method for an experimental example with AM control and probe frequencies, 0.7 MHz and 3.0 MHz, respectively; and we demonstrated agreement between the geometric method and the statistic-based quantile and standard deviation methods to within 15%.Good agreement between the methods was also found for using other probe frequencies (data not shown), nevertheless, more analyses are needed to establish regimes of operation, bounds for accuracy and consistency, etc.
The statistic-based approach offers a promising rapid and supervisor-free quantification method, making it suitable for nanoparticle analysis under a range of investigations.Fast, supervisor-free, 'on-the-fly' quantification is especially important for initial analysis and optimizing laboratory experiments and projects.There is scope for extending the methods to include negative DEP and for further investigating pDEP collection and release by varying the inter-electrode gap RoI geometry in conjunction with the use of statistical analyses.There is also scope for minimizing the excluded sample space so more of the array area is evaluated, exploring the effects of different electrode designs and experimental setups, and advancing the geometric method towards full automation -as required.
In this paper and earlier work, [4][5][6] the samples of nanoparticles are monodisperse and the dualcycle method has been shown to accurately predict the dielectric properties.Polydispersity in the size of the nanoparticles within the sample, on the other hand, will tend to introduce uncertainty in the accuracy of the measurement.As such we envisage our current method be used to: 1) assay monodisperse samples, or 2) assay a polydisperse sample, as part of an integrated Lab-ona-Chip system, with integrated upstream processing that involves size-based fractionation prior to measurement.
Although beyond the scope of current work, the effect of polydispersity on measurement accuracy requires further comment.We refer to the dual-cycle method, described earlier, for determining nanoparticle dielectric properties where the frequency of the probe cycle DEP collection phase is varied.The cycles could be influenced by the nanoparticle size in several ways.During the collection phase, the nanoparticle DEP velocity is proportional to the square of the nano-particle diameter so that large nanoparticles move more quickly than small nanoparticles.In addition, the nanoparticle conductivity tends to be inversely proportional to the nanoparticle diameter, so that changes in the frequency dependent DEP collections will occur.The effect has been observed experimentally, for example, 450 nm and 110 nm diameter beads were shown in 13 to exhibit crossover frequencies of about 2 and 6 MHz, respectively.During the release phase, the movement of nanoparticles from the electrode tips is driven by diffusion and tends to be inversely proportional to the bead diameter.This means that larger nanoparticles are released more slowly than smaller beads, and has been quantified in earlier work. 4,5 mbining these effects means that the nature of the dual-cycle process excludes the effect of very small and very large nanoparticles.Estimation of nanoparticle conductivity for monodisperse nanoparticles typically fits one frequency dependent Clausius-Mossotti (CM) curve [4][5][6]13 to the ratio of the measured probe and control initial collection rate data values, given by equation (23). Wihin the nanoparticle size range that moves in accordance with dual-cycling, the effect of polydispersity would "smear" the frequency dependent response of the data ratios and consequently make the estimate of conductivity less certain.
Errors in the measurement process due to variation in the fluorescence intensity of the sample are mitigated against due to the differential nature of the dual-cycle measurement process and the resulting normalisation of the fluorescence signal.Any effects due to photo-bleaching and variation of the fluorescence from the sample can easily be monitored, quantified and compensated through the use of the dual-cycle technique, as the method allows monitoring of variation in absolute fluorescent intensity over time as well as performing a differential measurement over shorter timescales.Such an analysis is therefore reasonably robust against the effect of photo-bleaching.In general, the reported work considerably advances the arsenal of DEP quantification and analysis techniques and tools available to researchers that can be implemented on standard academic and scientific computing platforms.

V. CONCLUDING REMARKS
Two important advancements in image quantification tools and methods for DEP-driven nanoparticle movement using planar castellated electrode arrays have been reported in detail.The first of the two quantification methods for estimation 2D spatial mean details a geometric approach that entails spatial filtering using almost fully-automated edge feature detection algorithms derived by first principles.The second approach uses statistic-based quantile and standard deviation methods, applied to each frame, that avoid the need for geometric feature recognition and are therefore automatic and supervisor-free.Evaluating initial collection rate ratios using the dual-cycle method, a technique previously developed for determining the dielectric properties of nanoparticles, the two methods were shown to concur to within 15% of each other.The geometric and statistic-based methods, hence, advance pDEP quantification techniques.They can also complement each other, acting not only to enhance current DEP image quantification methods, but also expand the scope of electrokinetic nanoparticle investigations.The system described used fluorescence microscope dielectrophoretic spectroscopy to image continuously cyclically pulsed DEP.However, there is ample scope for further implementation of our method to a range of other systems (e.g.other device geometries and non-DEP systems that rely on monitoring variation in optical intensities over time) and these are being explored at present.returns the x-ordinate, x q .Applying Eq. (B2) to (B3) yields ( The inverse error function, −1 (x), can be readily applied to the rhs of Eq. (B4) for finding, x q , −1 [2c q − 1] = k q (B5) where k q is a constant because c q is constant.By definition, the left hand side (lhs) of Eq. (B4) leads to −1 [ ( Thus, Eq. (B5) and (B6) lead to the relation for finding x q or horizontal abscissa, The cdf coordinates (c q , x q ) are shown by the 'plus', '+' printed on the blue and red curves in plotted in Fig. S.I.B.1(b). 36uppose that the mean and std are approximately linearly correlated σ = α + βμ where α and β are respective estimates for the abscissa and gradient, and k q remain constant (since c q is fixed).Substituting the linear relation Allowing the mean assume a new incremental value, μ = μ so that Eq. (B8) yields x q = μ (1 + βk q √ 2) + k q α√ 2, then the difference cancels the abscissa component and leads to the useful relation, Note that the relation μ ∝ x q applies to all values of the quantiles -as is demonstrated in . 36Applying Eq. (B3) to (B9) and Eq. ( 19) follows.That is, if we are not concerned with scale and the contribution from the std (e.g. in evaluating DEP collection rate ratios), the difference of the mean values of two normal distributions is proportional to the difference of the quantile functions applied to the two cdfs.The result is easily generalised to more points, and hence, can progress with time.Similar arguments apply where σ remains time invariant, and so forth.
It is useful to note that Eq. (B10) relies on the linear relation Eq. (B8) that, in turn, is consequence of expected value, E{X } = μ and variance, V ar{X } = σ 2 .In the case of a log-normal pdf, a linear relation can be obtained by transforming the data on the log-scale.

FIG. 1 .
FIG. 1. Scheme of image processing described for pDEP-driven nanoparticle collection and release onto planar castellated electrodes (a) experimental setup: microelectrode array, signal source (including amplifier and monitoring oscilloscope), fluorescence microscope and camera (b) close up of pDEP nanoparticle collection experiment using castellated planar microelectrodes with opposite polarity as shown (c) schematic showing two image processing methods (geometric and statistic-based) that lead to quantification of pDEP dynamics.
. S.1(a)(i) and S.1(b), the nanoparticles are initially uniformly distributed.Applying an RF voltage to the castellated electrodes, the action of the pDEP force causes downward nanoparticle movement, particularly near the electrode array where the DEP force is strong, Fig. S.1(a)(ii).The initial rate of collection, with respect to time, is the nanoparticle collection rate, and it progressively becomes limited, Fig. S.1(a)(iii) and S.1(b).Switching off the AC potential at t = t r initiates the release phase since there is no longer any pDEP force to trap the nanoparticles, and they diffuse into the bulk medium, Fig. S.1(a)(iv) and S.1(b).

FIG. 2 .
FIG. 2. (a)Image of part of a 5 μm planar castellated electrode array frame, and (b) expanded view showing unit cell pair as a white partly transparent oblong with dimensions as shown and x and −z coordinate system labeled (compatible with Matlab ij axis).The unit cell pair can be split into two cells since the top grid, 4 × 3 '+', of the unit cell pair, superposed on the castellation, is symmetric.This means it can be rotated about either the longitudinal or transverse axis, and forms the same pattern as the bottom grid of 4 × 3 '+' and castellation; likewise for the electrode interior.Longitudinal (direction) bay widths are slightly less than electrode tip 'crest' widths, but there is little advantage computationally to distinguish them so the cell is evenly split longitudinally to two 1/2 half-cells each with width w l /2, as illustrated.

I 1 )
(x, y, −z, t k )dy discrete (All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported license.See: http://creativecommons.org/licenses/by/3.0/Downloaded to IP: 144.82.107.91 On: Wed, 27 May 2015 13:35:27 An example of a 2D surface variation function is shown in Fig. S.2 36 where {w l , o l : 15 ≤ w l ≤ 35, 0 ≤ o l ≤ 20} pixels.Hence, edge locations in the longitudinal direction are generated from the pair of estimates, FIG. 3. (a) Mean of the 2D image in transverse x direction leads to fluorescence intensity in longitudinal −z direction, I l (c), with edges (dips) located, as shown, by triangles for each of 22 cells.Consecutive triangles, in the same direction, represent one longitudinal cell width.Estimates of w l and o l , are 26.3 and 14 pixels, respectively.(b) Mean of the 2D image in longitudinal −z direction leads to intensity in transverse x direction (__), I t (r ), with numerical labels on far left 'bottle' peak shaped feature.Subtracting the 'basin' fluorescence from the original profile (. . ..), levels out the valleys (representing electrode interiors), shoulders (labeled '1' and '4', castellation bays), necks (labeled '2' and '3', electrode tips) and peaks (gaps).Markers for gap/electrode tip interfaces (circle and square) and electrode bay edges (triangle and diamond) are found by software algorithms that automatically best-fit line segments along neck and shoulders as labeled.Lengths between triangles represent one cell, thus, consecutive right pointing triangles indicate two transverse cell pairs.

FIG. 4 .
FIG. 4. (a) Characteristic cell contour plot and (b) associated two DEP dual (control and probe) cycles, each cycle comprises a collection and release phase, as shown.The characteristic in (a) excludes the interior electrode region that is not of interest and uses 4 × 22 = 88 cells.Spatial 2D integration of the rectangular region of interest (RoI denoted by hashed outline) and enlarged left rectangle in (b), yields a numerical value I ∼ = 78 as shown.The template frame is an average of three frames (represented by hashed circle) before the first DEP collection.Contours extending between electrode tips, shown by the right RoI example in (b), represent DEP collections, e.g.associated with the second peak, located by small full circle shown (t = 47.2 s or frame 472).Repeating the operation, a sequence of 2D frames is imaged processed to form a sequence of 2D RoIs, each of which is then spatially integrated, hence forms a series of points, or time profile, I gm (t).

FIG. 5 .
FIG. 5. Comparison of frame probability density functions (pdfs), and inset, frame cumulative density functions (cdfs), associated with two corresponding time points shown in the DEP collection and release dual-cycles in Fig. 4(b).Each of the frame pdfs is comprised of two distributions as shown: the left mode is attributed to pixels associated with mainly in the interior electrode region, and the right mode is due to pixels located over the inter-electrode gap.Comparing two frame pdfs and cdfs for the points associated with template at t = 0.3 s and second dual-peak at t = 47.2 s shows that the effect of DEP nanoparticle collection is manifested in mainly in the tails of the right pdf, and hence the upper tail of the cdf (shown in inset).

FIG. 8 .
FIG.8.Initial collection rate ratio scatterplot showing almost linear correlation between frame standard deviation, σ F (t), and gap mean, I gm (t), for both control (solid) and probe (dash dot).The number of data points for each phase was 35.