Bayesian object classification of gold nanoparticles

Nanoparticles (NPs) are tiny particles of matter with diameters typically ranging from a few nanometers to a few hundred nanometers which possess distinctive properties. These particles, larger than typical molecules but too small to be considered bulk solids, can exhibit hybrid physical and chemical properties which are absent in the corresponding bulk material. The particles in their nano regime exhibit special properties which are not found in the bulk properties, for example, catalysis [Kundu et al. (2003)], electronic properties [Jana, Sau and Pal (1999)] and size and shape dependent optical properties [Jana and Pal (1999)], which have potential ramifications in medicinal applications and optical devices [Link and El-Sayed (1999), Kamat (1993)]. The current challenge is to develop capabilities to understand and synthesize materials at the nano stage, instead of the bulk stage.

Among the various NPs studied, colloidal gold (Au) NPs were found to have tremendous importance due to their unique optical, electronic and molecular-recognition properties [Hirsch et al. (2003) and Gaponik et al. (2000)]. For example, selective optical filters, bio-sensors, are among the many applications that use optical properties of gold NPs related to surface plasmon resonances which depend strongly on the particle shape and size [Yu et al. (1997)]. Moreover, there is an enormous interest in exploiting gold NPs in various biomedical applications since their scale is similar to that of biological molecules (e.g., proteins, DNA) and structures (e.g., viruses and bacteria) [Chitrani, Ghazani and Chan (2006)].

In recent years it has become possible to investigate the dependency of chemical and physical properties on size and shape of NPs, due to Transmission Electron Microscopy (TEM) images. Sau, Pal and Pal (2000) and Kundu, Lau and Liang (2009), respectively, showed size and shape dependence of synthesis and catalysis reaction where they observed different rates. They also observed that circular gold NPs are better catalysts compared to triangular NPs for a specific reaction. The development of new pathways for the systematic manipulation of size and shape over different dimensions is thus critical for obtaining optimal properties of these materials. In this paper we develop novel, model-based image analysis tools that classify and characterize the images of the NPs which provide their morphological characteristics to enable a better understanding of the underlying physical and chemical properties. Once we are able to accurately characterize the shapes of NPs by using this method, we can develop different techniques to control these shapes to extract useful material properties.

Substantial work in estimating the closed contours of objects in an image has been done by Blake and Yuille (1992), Qiang and Mardia (1995), Pievatolo and Green (1998), Jung, Ko and Nam (2008), Kothari, Chaudhry and Wang (2009), among others. Imaging processing tools, especially for cell segmentation, also exist; for instance, ImageJ [ImageJ (2004)] is a tool recommended by the National Institute of Health (NIH). However, the features of the data we are dealing with are quite different from those considered in the literature reviewed, as there are various degrees of overlapping of the NPs differing in shapes and sizes, as well as a significant number of NPs lying along the image boundaries.

High-level statistical image analysis techniques model an image as a collection of discrete objects and are used for object recognition [Baddeley and van Lieshout (1993)]. In images with object overlapping, Bayesian approaches have been preferred over maximum likelihood estimators (MLE). The unrestricted MLE approaches tend to contain clusters of identical objects allowing one object to sit on the top of the other, whereas the Bayesian approaches mitigate this problem by penalizing the overlapping as part of the prior specification [Ripley and Kelly (1977), Baddeley and van Lieshout (1993)], offering flexibility over controlling the overlapping or the touching.

In Mardia et al. (1997), a Bayesian approach using a prior which forbids objects to overlap completely is proposed to capture predetermined shapes (mushrooms, circular in shape). Inference is carried out by finding the Maximum A Posteriori (MAP) estimates and the prior parameters are chosen by simulation experience, in effect, fixing the parameters that define the penalty terms. Rue and Hurn (1999) also used a similar framework to handle the unknown number of objects but introduce polygonal templates to model the objects. However, their application is restricted to cell detection problems, where the objects do not overlap but barely touch each other and the method works more like a segmentation technique than as a classification technique. Moreover, the success of this approach depends on prior parameters, which are assumed known throughout the simulation. Al-Awadhi, Jenninson and Hurn (2004) used the same model except that they considered elliptical templates instead of polygonal templates and applied their method to similar cell images. All the above methods take advantage of the Marked Point Process (MPP), in particular, the Area Interaction Process Prior (AIPP), or any other prior that penalizes the overlapping or touching, which we explain later in the paper.

Since the structure of the data we are analyzing is different from the literature, we adapt object representation strategies discussed above to the problem at hand. When we refer to a shape, we refer to a family of geometrical objects which share certain features, for example, an isosceles and a right triangle both belong to the triangle family. There are five types of possible shapes of the NPs in our problem. The scientific reason is that the final shape of the particle is dominated by the potential energy and the growth kinetics. There is a balance between surface energy and bulk energy once a nucleus is formed. The arrangement of atoms in a crystal determines those energies such that only one of these specified shapes can be formed. We use similar scientific reasons to construct shape templates. These templates are determined by the parameters which vary from shape to shape.

Since there is a difference in the degree of overlapping from image to image, we assume that the parameters of the AIPP are unknown and ought to be inferred. This leads to a hierarchical model setting where the prior distribution has an intractable normalizing constant. As a result, the posterior is doubly intractable and we use the Markov chain Monte Carlo (MCMC) framework to carry out the inference. Simulating from distributions with doubly intractable normalizing constants has received much attention in the recent literature, but most of these methods consider the normalizing constant in the likelihood and not in the hierarchical prior; Møller et al. (2006), Murray, Ghahramani and MacKay (2006), and Liang (2010), among others. In this paper, we borrow the idea of Liang and Jin (2011), which is a modified version of the reweighting mixtures given in Chen and Shao (1998) and Geyer and Møller (1994), which can deal with doubly intractable normalizing constants in the hierarchical prior as well. The MCMC algorithm used can be described as a two-step MCMC algorithm. We first sample the parameters from the pseudo posterior distribution which is a part of the posterior that does not contain the AIPP normalizing constant—and then an additional Monte Carlo Metropolis–Hastings (MCMH) step that accounts for this normalizing constant.

Sampling from the pseudo posterior distribution is also quite challenging. Inferring the unknown number of objects with undetermined shapes is a complex task. We propose Reversible Jumps MCMC (RJ-MCMC) type of moves to handle both the tasks [Green (1995)]. Birth, death, split and merge moves have been designed based on the work of Ripley (1977), Rue and Hurn (1999). We also propose RJ-MCMC moves to swap (switch) the shape of an object. Using the above mentioned computational scheme, we obtain the posterior distributions for all the parameters which characterize the NPs: number, shape, size, center, rotation, mean intensity, etc. Owing to the model specification and the computational engine for inferring the model parameters, our approach extracts the morphological information of NPs, detects NPs laying on the boundaries, quantities uncertainty in shape classification, and successfully deals with the object overlapping, when most of the existing shape analysis methods fail.

The rest of the paper is organized as follows: Section 2 describes the TEM images, Section 3 deals with the object specification procedure, Section 4 describes the model specification, Section 5 describes the MCMC algorithm, Section 6 describes a simulation study and Section 7 applies the method to the real data. Conclusions are presented in Section 8.

In this paper we analyze a mixture of gold NPs in a water solution. In order to analyze the morphological characteristics, NPs are sampled from this solution onto a very thin layer of carbon film. After the water evaporates, the two-dimensional morphology of NPs is measured using an Electron microscopy such as TEM. In our case, a JEOL $2010$ high resolution TEM operating at $200$ kV accelerating voltage was used, which has $0.27$ nm of point resolution. The TEM shoots a beam of electrons onto the materials embedded with NPs and captures the electron wave interference by using a detector on the other side of the material specimen, resulting in an image. The electrons cannot penetrate through the NPs, resulting in a darker area in that part of the image. The output from this application will be an eight bit gray scale image where darker parts indicate the presence of a nanoparticle. The gray scale intensity is varying as an integer between $1$ and $256$. Refer to Figure 1 for examples of TEM images.

Due to the absorption of electrons by the gold atoms, the regions occupied by the NPs look darker in the image. The darkness pattern may vary according to specific arrangements of the atoms inside any single nanoparticle. Additionally, one can see many tiny dark dots in the background, which are uniformly distributed throughout the image region. These dark dots are generated because the carbon atoms of the carbon film also absorb electrons. One may also notice a white thin aura wrapping around the whole or partial boundary of a particle. This is the result of having surfactants on the rim of the particles. The surfactants are added to keep the particles from aggregating in the process of making colloidal gold. Analyzing the shapes of the NPs in a TEM image is primarily based on modeling them as objects, whose shapes are parametrized. Treating a nanoparticle as an object is the critical component of our modeling framework, which we discuss in the next section.

An object is specified in a series of steps that allow us to model a wide variety of shapes. They are as follows: (a) template, (b) shift, scale and rotate operators, (c) object multiplicity. We discuss each of them in detail below.

The MCMC algorithm used in this paper can be described as a two-stage Metropolis–Hastings algorithm. We first sample the parameters from the pseudo posterior distribution followed by a Monte Carlo Metropolis–Hastings step to account for $A^{*}$ [Liang and Jin (2011), Liang, Liu and Carroll (2010)].

The MCMC algorithm will have the following form:

• •

  Given the current state $\bolds{\Theta}^{k},\bolds{\gamma}^{k}$ draw $\bolds{\Theta}^{\prime},\bolds{\gamma}^{\prime}$ from $p^{*}$ using any standard MCMC sampler.

• •

  Given all the parameters, simulate auxiliary variables $z_{1},\ldots,z_{M}$ from the likelihood $z\sim f(z;\bolds{\Theta}^{\prime})$ using an exact sampler.

• •

  Estimate $R=\frac{A(\bolds{\eta}^{\prime},m^{\prime},\bolds{\gamma}^{\prime})}{A(\bolds{\eta}^{k},m^{k},\bolds{\gamma}^{k})}$ as

  $$\hat{R}=\frac{1}{M}\sum_{i=1}^{M}\frac{f(z;\bolds{\Theta}^{\prime})}{f(z;\bolds{\Theta}^{k})}\frac{\pi(\bolds{\Theta}^{\prime}|\bolds{\gamma}^{\prime})}{\pi(\bolds{\Theta}^{k}|\bolds{\gamma}^{k})}\frac{\pi(\bolds{\gamma}^{\prime})}{\pi(\bolds{\gamma}^{k})},$$

  which is also known as the importance sampling (IS) estimator of $R$.

• •

  Compute (estimate) the MH rejection ratio $\alpha$ as

  $$\hat{\alpha}=\frac{1}{\hat{R}}\frac{p^{*}(\bolds{\Theta}^{\prime},\bolds{\gamma}^{\prime})}{p^{*}(\bolds{\Theta}_{k},\bolds{\gamma}_{k})}\frac{Q(\bolds{\Theta}^{\prime},\bolds{\gamma}^{\prime}\rightarrow\bolds{\Theta}^{k},\bolds{\gamma}^{k})}{Q(\bolds{\Theta}^{k},\bolds{\gamma}^{k}\rightarrow\bolds{\Theta}^{\prime},\bolds{\gamma}^{\prime})}=\frac{1}{\hat{R}}.$$

  (9)

  The last equation is true since $Q=p^{*}$. So, the above approximates the normalizing constant of the posterior.

• •

  Accept $\bolds{\Theta}^{\prime},\bolds{\gamma}^{\prime}$ with probability $\min(1;\hat{\alpha})$.

Simulating auxiliary variables $z_{i}$ from the likelihood is straightforward and simply requires us to sample from normal distribution with parameters defined at the proposed state of the sampler. The challenge lies in drawing from the pseudo posterior.

A generalized Metropolis-within-Gibbs sampling with a reversible jump step is used to simulate from the pseudo posterior distribution with known number of objects. Additionally, a reversible jump MCMC (RJ-MCMC) with spatial birth-death as well as merge-split move is invoked to sample the number of objects and their corresponding parameters.

We draw from the joint pseudo posterior $p^{*}(\bolds{\mu},\bolds{\sigma}^{2},\bolds{\eta},\bolds{\gamma},m|y)$ by alternately drawing from the conditional pseudo posteriors of $\bolds{\mu},\bolds{\sigma}^{2}\bolds{\eta}|m,y,\bolds{\gamma}$, $\bolds{\gamma}|\bolds{\mu},\bolds{\sigma}^{2}\bolds{\eta},m,y$ and $m|\bolds{\eta},\bolds{\mu},\bolds{\sigma}^{2},\bolds{\gamma},y$ as follows:

• •

  Draw $\bolds{\eta}^{k+1},\bolds{\mu}^{k+1},\bolds{\sigma}^{k+1}$ from $p^{*}(\bolds{\eta},\bolds{\mu},\bolds{\sigma}|m^{k},\bolds{\gamma}^{k},y)$ using a Metropolis-within-Gibbs sampler.

• •

  Draw $m^{k+1}$ from the pseudo posterior $p^{*}(m|\bolds{\mu}^{k+1},\bolds{\sigma}^{k+1},\bolds{\eta}^{k+1},\bolds{\gamma}^{k},y)$ using a RJ-MCMC.

• •

  Draw $\gamma_{1}^{(k+1)},\gamma_{2}^{(k+1)}$ from the distribution $p^{*}(\bolds{\gamma}|y,\bolds{\Theta})$ using an M–H step.

We explain these steps in detail, in the following paragraphs.

In this section we use a simulation study to evaluate the performance of our proposed MCMC method. We consider two examples wherein a $200\times 200$ image with ten objects each are generated from the prior distributions described in Section 4.2 with area interaction parameter $\gamma_{2}=40$ and $\gamma_{2}=10$, respectively. The pixels inside each object have constant mean, which is different from object to object. The covariance matrix is chosen from a CAR model with parameters very close to the extreme dependence. Objects in both the example images have different morphological properties and belong to the five different shape families described in Section 2. The image used in the first example is shown in Figure 2(a) and the second is shown in Figure 2(b). For the

example images, the MCMC samples drawn from the posterior distribution of $\gamma_{2}$ are given in Figure 3. From these simulations, we can see that the Markov chain mixes well and the posterior mean is close to the true values we used to simulate the data. Values close to $40$ are drawn in example 1 [Figure 1(a)], while values close to $10$ are drawn in example 2 [Figure 2(b)]. A general observation in the

simulations is that the variance of the posterior distribution of $\gamma_{2}$ depends on the value of $\gamma_{2}$. For large values of $\gamma_{2}$ we observe relatively larger posterior variance than for small values. Another significant observation is that there is a dependence on the accuracy and the variance of the posterior distribution of $\gamma_{2}$ on the number and size of objects. To investigate this phenomenon, we fixed the value of $\gamma_{2}$ but simulated images with a different number of objects and sizes. As we increase the number and the size of objects, the posterior distribution of $\gamma$ will be closer to the true value. Below, we discuss two features of our method using the two examples.

Using the MCMC samples, we can obtain the distribution of the particle size, which is characterized by the area of the nanoparticle and the distribution of the particle shape. The aspect ratio, defined as the length of the perimeter of a boundary divided by the area of the same boundary, can be derived from the combination of size, shape and the pure parameters. The statistics of size, shape and aspect ratio are widely adopted in nano science and engineering to characterize the morphology of NPs, and are believed to strongly affect the physical or chemical properties of the NPs [El-Sayed (2001), Nyiro-Kosa, Nagy and Posfai (2009)]. For example, the aspect ratio is considered as an important parameter relevant to certain macro-level material properties because physical and chemical reactions are believed to frequently occur on the surface of molecules so that as the aspect ratio of a nanoparticle gets larger, those reactions are more active.

We apply our method to three different TEM images. The parameters that maximize the posterior distribution (MAP) obtained from the (MCMC) are presented in detail. Our classification results of particular type are verified by our collaborators with domain expertise; this manual verification appears the only valid way for the time being. More than $95\%$ of the NPs in those images are classified correctly. This also includes the particles in the boundary as well as having overlapping regions. For completely observed objects, there is almost $100\%$ correct classification.

We start our application with the image in Figure 1(a). Morphological image processing operations, such as watershed transformation and erosion, can be used to get an approximate count of the number of NPs in the model [Gonzalez and Woods (2007)]. They also can be used in initializing the MCMC chains and in constructing proposal distributions required by the MCMC sampler. The morphological image processing we used in this dissertation has the following steps: (1) image filtering and segmentation, (2) determining the number of objects, (3) estimating location, size and rotation parameters. We first transform the image from grey to a binary image and then apply watershed transformation to partition the image into subimages. In each binary subimage we apply erosion and dilation operations to find initial values for the parameters inside of each subimage. Because this morphological processing is not the subject of the present work, it is not presented in more detail. After the initial values are obtained from the preprocessing step, all five templates are randomly assigned for starting template specifications. From the MCMC sampler described in Section 5 we obtain a random sample of the posterior distribution for all the parameters which characterize the NPs, namely, the shape $\mathbf{T}$, the size $\mathbf{s}$, the rotation $\bolds{\theta}$, the random pure parameter $\mathbf{g}^{r}$, the mean intensity $\bolds{\mu}$ and the variance $\bolds{\sigma}^{2}$. We use this posterior sample for inferring the model parameters and extracting the morphological information of NPs with uncertainty in shape size and classification. To better present our results, we chose to work with the Maximum a posteriori (MAP) estimations of these parameters.

In Figure 7 we show the TEM image and MAP estimates of the parameters for $20\mbox{,}000$ MCMC sample. In Figure 8 we present the parameters of $\mathbf{s}$, $\mathbf{g}^{r}$ and $\bolds{\mu}$ that correspond to the MAP estimate for all the number of objects, $m$, corresponding to that value. Summary statistics of the shape parameters are given in Table 1. From the table and the histogram it is clear that the mean intensity is different from nanoparticle to nanoparticle, justifying our assumption of different means in (4.1). We also obtain the posterior probability of the classification for each of the objects. This probability depends on the complexity of the shape of the object. For example, object $2$ has been classified as an ellipse with probability $0.98$, whereas object $20$ has been classified as an ellipse with probability $0.68$ (circle with probability $0.32$). In Table 1 (and in all the following tables of this chapter) we presented the classification with the highest posterior probability of some of the nanopartiles. In this example we successfully deal with the object overlapping and objects laying on the boundaries.

Our second application deals with a more complex image shown in Figure 9. In this image at least $6$ overlapping areas and at least $6$ nanoparticles laying in the boundary are observed. More specifically, nanoparticles $1$, $2$, $3$, $14$, $15$, $16$, $18$, and $19$ lay in the boundary of the image while pairs 2–4, 3–4, 9–10, 10–11, 17–18, and 10–12 overlap. In this example, the overlapping is more complex and existing methods fail to represent the real situation. A number of nanoparticles are overlapped together forming a groups such as nanoparticles 9–10–11–12. MAP estimate values for all the parameters are obtained after $20\mbox{,}000$ MCMC iterations. Complex shapes have been classified accurately; see Figure 9. For example, nanoparticle $18$ has an incomplete image and it has been classified as a circle with posterior probability $0.77$. The MAP estimates of the parameters drawn from MCMC, namely, shape $\mathbf{T}$, size $\mathbf{s}$, rotation $\bolds{\theta}$, random pure parameter $\mathbf{g}^{r}$ and mean intensity $\bolds{\mu}$ are presented for the first six objects in Table 2. In this application, $11$ out of the $17$ objects are ellipses (E) and $6$ are circles (C) and one is a triangle (TR). We also present the histogram of the MAP estimates of parameters $\mathbf{s}$, $\mathbf{g}^{r}$ and $\bolds{\mu}$ in Figure 10. Summary statistics of various shape parameters are given in Table 2. We see from the table that our proposed algorithm captures triangles, circles etc. quite accurately.

Our next application deals with an image with $76$ nanoparticles with $4$ shapes; see Figure 1(b). In this image, few objects have overlapping areas and at least $10$ nanoparticles are laying in the boundary. Some objects do not have very clear shape like objects $29$ and $31$.

Different shapes are captured with different templates with the proposed method. In addition to the circles and ellipses which were successfully captured in the previous images, the triangles and squares are also captured accurately. Nanoparticles denoted by $29$ and $31$ are classified correctly, even if they have vague shapes; see Figure 11. In this example, out of $76$ nanoparticles, $47$ are classified as a circle, $23$ as an ellipse, $4$ as a triangle and $2$ as a square. Distribution of the various parameters of the identified objects are shown in Figure 12. In Table 3 we present all the triangular shapes in order to compare the pure parameter $h_{1}$. As we can see from the table, triangular shape nanoparticles $4$ and $12$ are closer to the equilateral triangle, with value close to $h_{1}=2.33$, while triangular shape nanoparticles $51$ and $57$ have wider sides, since their $h_{1}<2.3$.

In this image we can see more than $85\%$ percent of the nanoparticles are in the same shapes like circular or slightly tilted like ovals. Normally when we do shape controlled synthesis, we called it nano spheres or circular nanoparticles. Approximately five to ten percent of the other shapes or slight changes we usually neglect because in solution synthesis routes it is very difficult to synthesis $100\%$ of the same size and same shapes. However, if we consider critically the reason of shape evolution or statistical analysis of different shapes, then this small difference might be considered. We classify this particular example as spherical gold nanoparticles having almost the same size and shapes.

As a part of the verification process, we compare the accuracy of our method with that of the current practice used in nanoscience. In brief, the current practice is largely a manual process with support of image processing tools such as ImageJ Particle Analyzer (http://rsbweb.nih.gov/ij) and AxioVision (http://www.zeiss.com/), which have been popularly used for biomedical image processing. The results are shown in Figures 13 and 14.

The manual counting process, subject to the application of the above imaging tools, is necessitated by the low accuracy of the autonomous procedures. For three TEM images with overlaps among particles, our procedure recognized $95\%$ of the total articles compared to the 20–50% recognition rate of the ImageJ. Considering frequent occurrence of overlaps in the TEM images of nanoparticles, the existing software cannot be used as more than a supporting tool. We have also applied our method to other images with the same success, encouraging its applicability.

We adopted a Bayesian approach to image classification and segmentation simultaneously and applied it in TEM images of gold nanoparticles. The merit of our development is to provide a tool for nanotechnology practitioners to recognize the majority of the nanoparticles in a TEM image so that the morphology analysis can be performed subsequently. This can evaluate how well the synthesis process of nanoparticles is controlled, and may even be used to explain or design certain material properties. Several factors like kinetic and thermodynamic parameters, flux of growing material, structure of the support, presence of defects and impurities can affect the morphology of NPs. In the future, we are planning to perform a factorial type experiment to identify the significant factors for morphological study. These significant factors can be properly controlled to develop NPs of required shapes.

From the experimental point of view, several improvements of existing techniques will be helpful to characterize the shape of the NPs. One is TEM tomography that allows to image an object in three dimensions, by automatically taking a series of pictures of the same particle at different tilt angles [Midgley and Weyland (2003)]. Another improvement of TEM is environmental HRTEM that is able to image nanoparticles, with atomic lattice resolution, at various temperatures and pressures [Hansen et al. (2002)].

From the modeling point of view, we used marked point process to represent the NPs in the image, where points represent the location of NPs and marks represent their geometrical features. More specifically, we treated the NPs in the image as objects, wherein the geometrical properties of the object were largely determined by templates and the interaction between the objects was modeled using the area interaction process prior. By varying the template parameters and applying operators such as scaling, shifting and rotation to the template, we modeled different shapes very realistically. In our current applications, we chose circle, triangle, square and ellipse as our templates. Other templates can be also constructed in the same framework. To solve the intractability of the posterior distribution, we proposed a complex Markov chain Monte Carlo (MCMC) algorithm which involves Reversible Jump, Metropolis–Hastings, Gibbs sampling and a Monte Carlo Metropolis–Hastings (MCMH) for the intractable normalizing constants in the prior. The first steps deal with simulating from a pseudo posterior distribution without involving the random normalizing constant. A generalized Metropolis-within-Gibbs sampling with a reversible jump step is used to simulate from a pseudo posterior distribution given the number of objects. Additionally, a reversible jump MCMC with the use of birth-death and merge-split moves is invoked on moving from a state with a different number of objects. Finally, we simulate from the intractable normalizing constant posterior using Monte Carlo Metropolis–Hastings where the acceptance ratio of the sample taken from the pseudo posterior is estimated by simulating from an auxiliary variable. We reported the posterior summary statistics of the shapes and the number of objects in the image. We successfully applied this algorithm to real TEM images, outperforming convention tools aided by manual screening. Our proposed methodology can help practitioners to associate morphological characteristics to physical and chemical properties of the NPs, and in synthesizing materials that have potential applications in optics and medical electronics, to name a few.