Analysis of Brownian motion trajectories of non-spherical nanoparticles using deep learning Open Access

To evaluate the properties of heterogeneous nanoparticles with mixed particle types and sizes, it is necessary to measure and analyze multiple physical properties for each identified particle. For example, as new therapeutic and diagnostic technologies using extracellular vesicles and artificial nanoparticles attract attention, evaluation techniques for nanoparticles and nanocolloids that support safety and reliability are becoming increasingly important, but there are many challenges and limitations to measurement in the 10 to 100 nm range.^1,2

Nano-tracking analysis (NTA),^3,4 which captures the Brownian motion (BM) through scattered light imaging and obtains the particle size from each trajectory, is widely used as a technique for measuring micro- to nano-sized single-particles in liquids due to its simplicity and the time and cost advantages of the measurement. However, the properties that can be extracted with NTA are generally limited, such as the inability to directly see the shape of the particles. It is well known that nanoparticle shape strongly influences optical properties, such as plasmon resonance^5–7 and Raman enhancement,⁸ but it is also an important factor indicating the degree of aggregation,^9,10 chemical properties such as solubility,¹¹ drug delivery performance,¹² and material toxicity.^13–15 We, therefore, tried to extract shape information from NTA data without significantly changing the experimental method.

It is also known that BM is affected by the shape of nanoparticles by previous studies on micrometer-sized particles¹⁶ and stochastic simulations¹⁷ in the long history of BM research from basic physics and analytical methods to the applied technology.^16–28 Therefore, in principle, it is expected that particle shape information can be detected by analyzing the trajectory of BM (time-series coordinate data). However, it is practically difficult to detect shape anisotropy of nanoparticles using NTA measurement data because nanoparticles have relaxation times that are many orders of magnitude shorter than micrometer-sized particles (theoretical rotational relaxation time is nine orders of magnitude smaller from the following definition: τ_r = πηd³/2k_BT), and the actual measured data contain experimental errors.

The analysis method also has points to be discussed. Although Stokes’ law is applied to the analysis of nanoparticles by assuming a rigid sphere, actual nanoparticles are often not spherical. In the case of non-spherical particles, a correction factor called the dynamic shape factor¹⁸ is introduced to consider translational BM, and the volume equivalent diameter was used as a solution to apply Stokes’ law. Both theoretical and experimental studies have been conducted to obtain the dynamic shape factor.^19–22 For rotational BM, the similar analysis can be performed as for spherical particles by including a factor of shape in the rotational friction coefficient to represent the rotational diffusion coefficient.^27–29 This means that if the optical measurement results for non-spherical particles are analyzed using the “spherical” Stokes–Einstein equation, the analytical values will include errors due to the shape factor not applied. Since the NTA method does not allow for direct observation of particle shape, it is impossible to determine the shape of each particle and then apply an appropriate analytical formula.

A promising solution is to apply deep learning (DL) to the analysis.^30,31 It is superior in extracting hidden correlations in large datasets, even when the data acquired are averaged, and is used to analyze the diffusion characteristics of nanoparticles, such as the Hurst exponent of fractional BM and the diffusion coefficient of BM,³² classification of diffusion modes,³³ and classification of particle trajectories in heterogeneous media from a number of diffusive states.³⁴ Deep learning analysis may also be able to extract shape information from statistically similar BM trajectory patterns, as in the example of classifying particles by shape based on their scattered light patterns.³⁵ When extracting shape information from the trajectory of BM, there is no explicit formula that directly links the two. In previous studies of BM of non-spherical particles using DL, shape is usually given as a known factor, and no previous studies have identified the particle shape from the BM trajectory measured by NTA. Therefore, in this study, we have developed a model that can estimate shape classification from time-series trajectory data of BM using DL.

Development of the shape estimation model consisted of three stages: NTA measurements for raw data acquisition, creation of datasets and models for DL, and DL.

Brownian motion trajectories of spherical gold nanoparticles (AuNP, Citrate capping, 80 nm diameter, Sigma-Aldrich) and gold nanorods (AuNR, Citrate capping, 40 nm diameter × 180 nm length, Nanopartz) were measured by NTA. A microcapillary chip-based particle analysis system³⁶ based on the NTA method developed for characterizing nanoparticles in liquid, which can simultaneously measure particle size and zeta potential, is shown in Fig. 1(a). The nanoparticles in the microfluidic channels made of cyclo-olefin polymers [Fig. 1(b)] are irradiated with a sheet laser (405 nm, 100 mW) from the side of the channel, and the scattered light is focused by an objective lens below the channel and formed by a high-sensitivity camera to obtain a high-contrast dark field image every 10 ms (100 fps) for 100 frames [Fig. 1(c)]. The center of gravity of the scattered light from each nanoparticle in consecutive 100-frame images was tracked [Fig. 2(a)] to obtain a trajectory [Fig. 2(b)].

The difference between adjacent frames of the time-series 2D coordinate data was converted into a displacement vector. Statistics for each BM trajectory in the displacement vector is shown in Table I. The time series trajectory data were given a correct answer label (AuNP: 0, AuNR: 1) indicating whether the trajectory related to spherical or rod-shaped particles and was used as the DL dataset [Fig. 2(c)]. In evaluating all DL models, the dataset was randomly divided into training data (80%), validation data (10%), and test data (10%). For shape prediction, the dataset was randomly divided into training data (90%) and validation data (10%). As the test data, another dataset was created and loaded.

Here, we designed and built DL models using Keras based on three orthodox techniques: multilayer perceptron (MLP), long short-term memory (LSTM),³⁷ and one-dimensional convolutional neural network (1D CNN).³⁸ MLP is a basic Fully Connected Neural Network (FCNN)³⁹ and was used to compare the performance of the models. LSTM is a derivative model of Recurrent Neural Networks (RNNs)^40,41 and can learn features of long-term time series data by introducing memory cells (LSTM block) into RNN and widely used for analysis and prediction of continuous data.⁴² Three gates (input, forgetting, and output) work together to learn how information is put in and out of the cell, control the balance between long-term memory and short-term memory, and retain information for the long term.⁴³ 1D CNN performs convolution processing in the time axis direction, extracts local features in the convolution layer, and summarizes them in the pooling layer.⁴⁴ Compared to FCNN, it is possible to recognize spatial features with high accuracy with fewer parameters. Furthermore, in order to consider both the time series and the correlation with the surroundings, we designed and constructed a 1D CNN+Bi-LSTM composite model that can extract local features by one-dimensional convolution and accumulate trends in a bidirectional LSTM structure.

The structure of our 1D CNN+Bi-LSTM model and a typical flowchart when performing shape prediction is shown in Figs. 3(a) and 3(b), respectively. Input data are the array with 100 × 1 elements. A filter with 7 × 1 elements is employed in first two convolution layers and that with 5 × 1 elements for the third convolution layer. The number of channels of all filters is set to be 16. Bidirectional LSTM blocks with five layers are installed instead of pooling layers. As the activation functions, we used ReLU for 1D conversion, hyperbolic tangent for the LSTM block, and softmax for the output layer. Binary cross-entropy loss is employed for the loss function, and adaptive moment estimation (ADAM)⁴⁵ is employed for the optimization of the loss function. The dropout⁴⁶ with the rate of 0.5 is employed in the affine layer to avoid the overfitting. The architecture and flowchart of MLP, LSTM, and 1D CNN models we used in this study are shown in  Appendix A. To implement the deep learning models, we used Python 3.6.9 as the programming language and implemented the shape estimation model using the deep learning libraries TensorFlow and Keras 2.2.0.

where k_B is the Boltzmann constant (J/K), T is the temperature (K), η is the viscosity (Pa s), L is the rod length (m), and d is the rod diameter (m).

The median particle size determined by NTA measurements was 113.3 nm for AuNP and 123.5 nm for AuNR in this experiment [Figs. 4(a) and 4(b)]. The hydrodynamic diameter d_h of particles less than 100 nm in diameter has been measured to be about 5%–10% larger,⁴⁹ but a larger difference was found in our case. Considering above factors, trajectories of 5000 AuNP and 5000 AuNR particles of 95–125 nm and 110–140 nm in diameter, respectively, were used in datasets for model development.

Using time series data with varying trajectory lengths (20, 40, 60, 80, and 100 frames), we verified the convergence of learning by changing the hyperparameters for each number of frames for the four models (MLP, LSTM, 1D CNN, and 1D CNN+Bi-LSTM). Typical training/validation loss and accuracy curves at 100 frames for each model are shown in  Appendix B. The conditions showing good evaluation indices stably were used as the optimal hyperparameters for the following calculations. The shape classification accuracy of the four models was evaluated at their respective optimal hyperparameters [Fig. 5(a)]. The accuracy increased with trajectory length for all models, indicating that data length is an important factor in revealing features depending on the shape. The area under the curves (AUCs) of receiver operating characteristic (ROC) curves for four models [Fig. 5(b)] and the corresponding ROC curves [Fig. 5(c)] showed similar trends as accuracy. Here, the true positive rate represents the proportion of AuNR (rod shaped particles) correctly determined as AuNR (non-spherical), and false positive rate represents the proportion of AuNP (sphere particles) falsely determined as AuNR (non-spherical). The high accuracy above 80% at 100 frames for the LSTM and 1D CNN model [Fig. 5(a)] indicates that both the extraction of local features by convolution and accumulation of temporal dynamics are effective for extracting shape characteristics. At the same time, the high accuracy indicates that shape classification of nanoparticles in liquids has reached a realistic level for single-particle analysis using NTA with DL.

As seen above, the 1D CNN+Bi-LSTM model achieved accuracy, AUC, and ROC comparable to those of the 1D CNN and LSTM models, suggesting that it is possible to combine the characteristics of both models in one model. It means that local features can be extracted while capturing features of temporal dynamics. In actual BM, multiple factors, such as the particle size, shape, and surface condition, are involved, and it is considered that a DL model specialized for a single feature value cannot fully analyze the actual phenomenon. Therefore, rather than selecting a model only based on the high accuracy rate, we adopted the 1D CNN+Bi-LSTM model, which can consider multiple physical effects with the same accuracy as the conventional architectures, because it is more suitable for shape prediction.

Mixed suspensions of AuNP and AuNR in different mixing ratios (7:3, 5:5, 3:7 by number) were prepared for NTA. Figures 6(a)–6(c) show the particle size distribution of mixtures. From the trajectory data of each mixture, a test dataset (total 500 time-series displacement vectors) was made and labeled spheres. The training/validation dataset was created from measurements of 9000 AuNP and 9000 AuNR particles without removing data containing errors that appear as size distribution, unlike model evaluation described above. It is for including equivalent error in the training/validation data and test data, as it is difficult to get a correct error reading (i.e., determine the optimal size range) for the mixed samples that serve as the test samples. To check the effect of experimental errors, such as agglomeration, we also generated and verified a pseudosample from unmixed AuNP and AuNR data according to the experimental mixing ratio [Fig. 6(d)]. Using the data of the mixed sample and the pseudosample at each mixture ratio as the test sample, shape prediction was performed with a 1D CNN+Bi-LSTM model to obtain the content of particles predicted to be non-spherical [Fig. 6(e)]. The difference between the value of predicted content (average of eight runs) for the mixed sample and generated pseudosample was from 0.2% to 5.6%. The two regression lines are in good agreement, and the difference between the values of predicted content and that calculated from the regression line is within 2.6%. This indicates that the regression line is valid as a calibration curve and allows for the estimation of the mixing ratio from the predicted content. For example, in the present case, a mixed solution predicted to have 50% non-spherical content would be considered to contain ∼43% non-spherical particles. We have shown that the shape of nanoparticles, for which information cannot be obtained only by ordinary NTA, can be distinguished when DL is used for analysis.

With the DL method, particles with different shapes, even if they have the same hydrated diameter, can be distinguished from a mixture based on their trajectories. In comparison, conventional NTA methods cannot distinguish between particle shapes, and nanoparticles of different shapes in a mixture cannot be distinguished unless there is a large difference in hydrated particle size or scattered light intensity. This method has the potential to open a new approach to the fundamental study of BM of nanoparticles in liquids.

We attempted to determine the shapes of two types of particles in this paper, but considering the types of shapes of commercially available nanoparticles, we think that this method can be used in practical applications, such as the detection of foreign substances in homogeneous systems. It would be interesting to extend the measurement object to particles of various shapes and materials, and it is a future research topic to examine the applicability of the DL+NTA method.

Here, we demonstrated the effectiveness of a new nanoparticle characterization method that combines DL and NTA measurements. For applying DL to general NTA data, we constructed a 1D CNN+Bi-LSTM model that predicts the shape of single nanoparticles in liquid with an accuracy of more than 80%. Using this model, the mixing ratio of the AuNP and AuNR mixed sample can be determined with a very small error: the difference between the values of predicted content and that calculated from the regression line is within 2.6%. We demonstrated that the shape of nanoparticles could be distinguished in a homogeneously dispersed system with little difference in particle size. It was shown that the shape of nanoparticles, which cannot be obtained by conventional NTA alone, can be identified by applying DL to BM trajectory analysis.

The datasets created for this paper are available within the supplementary material.

The authors are grateful for the technical assistance given by Ms. H. Kishita. This research was financially supported by the Japan Science and Technology Agency (JST) through the open innovation platform for industry-academia co-creation (COI-NEXT) program (Grant No. JPMJPF2202).

The authors have no conflicts to disclose.

Hiroaki Fukuda and Hiromi Kuramochi contributed equally to this work.

Hiroaki Fukuda: Formal analysis (supporting); Investigation (equal); Methodology (equal); Software (lead); Writing – original draft (supporting); Writing – review & editing (equal). Hiromi Kuramochi: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Yasushi Shibuta: Methodology (equal); Supervision (supporting); Writing – review & editing (equal). Takanori Ichiki: Conceptualization (equal); Funding acquisition (lead); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.

The architecture of MLP, LSTM, and 1D CNN models is as follows: Typical flowcharts are shown in Fig. 7.

MLP: Input data of MLP are the array with 200 × 1 elements. The 2D output was obtained from a quadruple stack of 1000 layers, 500 layers, and two 100 layers. We used ReLU for activation in all layers except output layer (softmax). The dropout with the rate of 0.5 is employed in all stages. Binary cross-entropy loss is employed for the loss function.

LSTM: Input data of LSTM are the array with 200 × 1 elements. The 2D output was obtained with a structure in which ten layers of LSTM are stacked in three stages. A dropout rate of 0.05 is adopted for all layers. We used hyperbolic tangent as an activation function, sigmoid for recurrent activation, and softmax for the output layer. Binary cross-entropy loss is employed for the loss function.

1D CNN: Input data of 1D CNN are the array with 100 × 1 elements. A filter with 7 × 1 elements is employed in first two convolution layers and that with 5 × 1 elements for the third convolution layer. The number of channels of all filters is set to be 16. Max pooling with a pool size of five is employed. As the activation functions, we used ReLU for 1D conversion and softmax for the output layer. Binary cross-entropy loss is employed for the loss function, and ADAM is employed for the optimization of loss function. The dropout with the rate of 0.5 is employed in the affine layer to avoid the overfitting.

As shown in Fig. 8, for three models (LSTM, 1D CNN, and 1D CNN+Bi-LSTM), the gap between the training and validation curves of both accuracy and loss is small; the training was performed properly and converged within the implemented epochs. MLP has a large gap, and accuracy and loss are unstable even if epochs are increased.