PHYSICAL REVIEW E 91, 042308 (2015)

Principal-component analysis of particle motion H. Y. C h e n ,12 Raphael L iegeois,3 John R. de Bruyn,2’* and A ndrea S oddu2i*-t 1Department o f Nuclear Science and Technology, Fudan University, Shanghai 200433, China 1Department o f Physics and Astronomy, University o f Western Ontario, London, Ontario, Canada N6A 3K7 3Montefiore Institute, Universite de Liege, 4000 Liege, Belgium (Received 22 July 2014; published 15 April 2015) We demonstrate the application of principal-component analysis (PCA) to the analysis of particle motion data in the form of a time series of images. PCA has the ability to resolve and isolate spatiotemporal patterns in the data. Using simulated data, we show that this translates into the ability to separate individual frequency components of the particle motion. We also show that PCA can be used to extract the fluid viscosity from images of particles undergoing Brownian motion. PCA thus provides an efficient alternative to more traditional particle-tracking methods for the analysis of microrheological data. DOI: 10.1103/PhysRevE.91.042308

PACS number(s): 83.85.Ns, 83. lO.Pp, 05.45.Tp, 02.70.Rr

I. INTRODUCTION Principal-com ponent analysis (PCA) [1,2], also know n as K arhunen-L oeve transform ation, is a statistical method for extracting patterns from com plex data sets and reducing the dim ension o f com plex signals. PCA involves transform ing a set o f data into a linear superposition o f orthogonal com ponents, arranged such that the first principal com ponent has the largest possible variance; i.e., it accounts for the largest contribution to the variation in the data. Successive principal com ponents sim ilarly have the m axim um possible variance subject to the orthogonality condition. By construction, the first few principal com ponents typically contain m ost o f the inform ation em bedded in the data, m aking PCA the basis o f very efficient data com pression and spatiotem poral signal­ processing m ethods [3], PCA has im portant applications in fields such as biom edical science [4,5], neuroscience [6,7], environm ental science [8-10], im age com pression [11,12], and image analysis [13], In this paper, we dem onstrate the application o f PCA to the analysis o f particle m otion, prim arily in the context o f m icrorheological m easurem ents. M icrorheology involves the analysis o f therm ally driven or externally forced m otion o f m icron-scale tracer particles suspended in a com plex fluid to determ ine the local viscoelastic properties o f the fluid or to study its m icrostructure [14-16], In an active m icrorheological experim ent, an external force is applied to tracer particles using, for exam ple, optical tw eezers, and the resulting response o f the particles is analyzed to give the inform ation about the viscoelastic properties o f the fluid [14,16-21]. In a passive m ultiple-particle-tracking experim ent, the positions o f many tracer particles are recorded over tim e and the resulting time series of im ages is analyzed to determ ine the m ean squared displacem ent o f the particles as a function o f lag tim e r [14-16,22,23], The local frequency-dependent viscous and elastic m oduli can be determ ined from the mean squared displacem ent using a generalized Stokes-Einstein rela­ tion [15,24], Both types o f experim ent produce a set o f

spatiotem poral data in the form o f a time series o f im ages w hich is w ell suited for analysis using PCA. PCA differs from m ost com m only used particle-tracking analysis m ethods in that the entire spatiotem poral data set is analyzed as a whole, as opposed to identifying individual particles in each image and reconstructing their trajecto­ ries [25,26]. Using sim ulated data, we dem onstrate in Sec. IIA that PCA can resolve periodic particle m otion with a range o f frequencies, and in Secs. IIB and IIC we show that PCA provides an alternative to conventional particle-tracking algorithm s for determ ining the viscosity o f N ew tonian fluids from m icrorheological data and for studying non-N ew tonian behavior, respectively. O ur results dem onstrate that PCA is a prom ising tool for particle m otion analysis. D etailed introductions to PCA can be found in several excellent references [1,3,27,28], The decom position o f a data set into its principal com ponents is illustrated schem atically in Fig. 1. We consider spatiotem poral data in the form o f a tim e series o f im ages. Each image consists o f m pixels, and im ages are recorded at discrete tim es t = 1 , 2 , . . . ,7". W hile in this paper we focus on the analysis o f tw o-dim ensional particle motion, there is no restriction on the spatial dim ension o f the data, and it is straightforw ard to apply our analysis to three- or higher-dim ensional images. The entire data set can be represented as an m x T m atrix X, each colum n o f w hich is a vector x, containing the im age data recorded at a particular time. Let W be a m atrix w hose rows are the n principal-com ponent vectors w,-, norm alized so that the m odulus squared |w ,|2 = 1. A matrix Y can be calculated via the transform ation

Y =

= wx =

W| w2

and y,-, the / th row o f Y, is the tim e-dependent am plitude of the principal com ponent w,-, as illustrated in Fig. 1. The principal com ponents w,- are the eigenvectors o f C = X X 1, w hich is proportional to the covariance m atrix o f X. The m atrix D, given by

These two authors contributed equally to this work. [email protected] 1539-3755/2015/91 (4)/042308(6)

"yi" y2

D = Y Y r = W X X 7W r = W C W r , 042308-1

(2)

©2015 American Physical Society

PHYSICAL REVIEW E 91, 042308 (2015)

CHEN, LIEGEOIS, DE BRUYN, AND SODDU

y.

X

0

0



7 0 0 I 0

0



0









0





0

W ,„

(

W 31

)

W 2„

)

W 3*

)

T

FIG. 1. (Color online) A schematic illustration of the decompo­ sition of an m x T data matrix X into its principal components w;. The vectors y, are the time courses of the corresponding principal components.

600 x 800 pixels, with a simulated frame rate of 10 s-1. Five groups of 10 particles were placed randomly in the first image. The particles moved in simple harmonic motion with random orientation and an amplitude of 50 pixels, with each group oscillating at a different frequency. The data were decomposed using PCA [29] into principal components w; and their respective time courses y,-. The results are summarized in Fig. 2. W|, the principal component with the largest variance, is shown in the left-hand panel of Fig. 2(a). Its time course yi is shown in the inset in the right panel of Fig. 2(a), and the

is diagonal, and the diagonal elements of D are its eigenvalues. The ith diagonal element of D, which we label as \ is the variance associated with w,-. The principal components are sorted according to magnitude of V (,\ with = max['D(,)} [1], For low-noise signals, the amount of information embed­ ded in each principal component is related to the corresponding variance: components with a larger variance carry a larger portion of the information from the original signal. Most of the information in the signal is therefore preserved by keeping only the first few principal components. X is conventionally centered on zero by subtracting its mean before applying PCA, but this process inevitably results in the loss of the physical information embedded in the average signal. To avoid this, we did not center X in the present work. In this case, the first principal component is just the average signal. The PCA calculations in this paper were performed using software [29] written in Matlab [30], For concreteness in the following discussion, we assume that each image in the data set has a fixed number of particles n in the field of view, and that each particle occupies k pixels. Assuming that the particle concentration is small enough that particles do not overlap (as would normally be the case in a microrheological experiment [23]), p = kn pixels will be occupied. Occupied pixels have a value /, and vacant pixels are 0, as illustrated in Fig. 1. The pixel values change with time as the tracer particles move. As noted above, each twodimensional image is reshaped to form one column of the data matrix X, so that each row of X corresponds to a given spatial location (pixel) and each column to a particular time. -

II. RESULTS A. Periodic motion In a typical active microrheology experiment, tracer parti­ cles would be driven, and would all respond, at a single known frequency. In other applications, however, nonlinear coupling between individual oscillators can lead them to respond at different frequencies. One example comes from neuroscience, in which coupled oscillations in cortical circuits have been studied using the Kuramoto model [31]. We demonstrate here that such systems can be studied by taking advantage of the frequency-resolving ability of PCA, which we illustrate by analyzing images of particles undergoing periodic motion. We created a simulated data set comprising 300 images of

0.03 - 0.02 - 0.01

0

0.01 0.02 0.03

Frequency (Hz)

FIG. 2. (Color online) The results of PCA of simulated data for particles undergoing simple harmonic motion. The left-hand panels of (a)-(e) show the first, second, fourth, fifth, and ninth principal components, respectively. The inset on the corresponding right-hand panel shows the time course of the principal component, and the power spectrum is plotted in the main graph of the right-hand panel, with all y-axis units arbitrary, (a) The first principal component is the time average of the data. Its time course is roughly constant, and there is no strong peak (other than that at zero frequency) in the power spectrum. The other principal components shown in (b)-(e) selectively pick out oscillations at frequencies 0.5 Hz, 0.167 Hz, 0.2 Hz, and 0.1 Hz respectively. Oscillations at 0.056 Hz were not isolated by any principal component.

042308-2

PRINCIPAL-COMPONENT ANALYSIS OF PARTICLE MOTION

Fourier power spectrum of jq in the main plot on the right of Fig. 2(a). This principal component represents the average state of the system and, as discussed above, could be eliminated by subtracting the time average of X from X. The probability distribution function of the position z of a particle undergoing harmonic oscillation about the origin with an amplitude A is (32] f ( z ) = 7T~i(A2 —z2)~'/2, and the image of W| reflects this distribution for each particle in the data set. Since it includes information from all particles, its time course yi is roughly constant; its power spectrum is peaked at 0 frequency and shows no other significant features. The other leading principal components tend to pick out oscillations of a certain frequency. As illustrated in Fig. 2(b)—2(e), Wi. w_t, W5, and W9 are dominated by motion at frequencies 0.5, 0.167, 0.2, and 0.1 Hz, respectively, these being the oscillation frequencies of four of the five groups of particles in the data. The most prominent features in the images of these principal components are the end points of the trajectories of particles oscillating at the selected frequency, and the power spectra of the corresponding time courses y2, y4, ys. and y9 have single strong peaks at the selected frequency. For our simulated data, PCA is able to cleanly separate motion at frequencies that differ by only 20%, and its frequency selectivity would presumably be higher for a longer data set. The fifth set of particles in the data oscillated at 0.056 Hz. This low frequency was not detected by our analysis, presumably because the length of the data set, which corresponded to 1.68 oscillation periods at this frequency, was too short. Other low-index principal components not shown in Fig. 2 contain either harmonics or combinations of the fundamental particle frequencies. As seen in the results shown in Fig. 2, PCA tends to rank the principal components from high to low frequency, indicating that higher-frequency oscillations contribute more variance to the data. Similar results were obtained from data in which small stochastic perturbations were added to the oscillatory motion, indicating that the analysis is robust with respect to small amounts of noise.

PHYSICAL REVIEW E 91. 042308 (2015)

probability given by the distribution m ,n j)

jij)2'

(5)

2a 2

p W|,- =

{Hi - t t j f 2a2

exp i= 1

(6)

where a is a normalization factor to be determined. From Eq. (2), the variance associated with the first principal component is then £>(1) = w,XXTw[ = I 2

( j 2 WR t — 1 \i=occ

where the second sum is over the p occupied pixels in the image. For simplicity we now consider the motion of a single particle occupying a single pixel, in which case p = 1 and we can drop the j subscript. The variance in the first principal component is then ( H - v ) 2'

a exp

£>(') = I 2 J 2

(H

I 2T

/{““p[-

i zt

— (TV

r L

(8 )

2a2 H?

-

|

r

!

f{H,n)dH

2(7 2

( ? - m )2 i 2

3

d'H

a2

e x p

(9) ( 10)

l 2T or GO



The normalization condition is

PCA can be used to determine the viscosity of a Newtonian fluid from a series of images of tracer particles diffusing in the fluid. These tracers undergo Brownian motion, and their mean squared displacement (A r2) grows linearly with time t : (Ar 2(t)) = 2clDt, where d is the dimensionality and D is the diffusion coefficient. For Brownian motion, the distribution of particle displacements Ar over a time interval period At is Gaussian with a standard deviation a = s/ldDAt.

ts/2n

-

where a is a function of A t as in Eq. (3). To maximize the variance V (U, the elements of the first principal component w, should have values determined by Eq. (5), summed over the p occupied pixels. If the position £ corresponds to the i th row of the data matrix X, then the f th component of W| will be

V3

B. Viscosity determination

(Hi

exp

|W | | 2 —

(Hi

f

-

I a exp

)2

(12)

2a 2

1=1 ‘

{H -ri2 2a2

dH

oTa^frt = 1,

(13) (14)

from which a % [1/( ct7TI/2)]i/2. This finally gives

(3)

I 2T

T)(l)

For a single particle, the viscosity 77 is related to D by the Stokes-Einstein relation,

(15)

r^/37r More generally, when p >

t] — kBT / 2 n a d D,

(4)

with k B the Boltzmann constant, T the absolute temperature, and a the radius of the particle [33]. Although hydrodynamic interactions between particles causep to increase with particle volume fraction [34], we neglect any such effects in this work. From Eq. (3), a particular particle j whose original position is jij will appear at position £,■ in a subsequent frame with a

1 4. p - l \ % ( _2_ v£ 3

2

I 2T

(16)

J a J7t

The total variance of all m principal components is

042308-3

m

v

101 =

m

=

t r ( c

!=l

}

=

E

T E

4

1= 1 t = 1

=

/ 2 ^

,

(17)

PHYSICAL REVIEW E 91, 042308 (2015)

CHEN, LIEGEOIS, DE BRUYN, AND SODDU

where x,; is the (i,t) element of X. X>tot is constant for a given data set, regardless of the details of the motion of the particles. Using this, the fraction of the total variance accounted for by the first principal component is £>(') S(,) =

± + P^l V3

V™

2

(18)

p o ^ /l T

In a typical multiple-particle-tracking experiment [23], p implying that

Principal-component analysis of particle motion.

We demonstrate the application of principal-component analysis (PCA) to the analysis of particle motion data in the form of a time series of images. P...
4MB Sizes 0 Downloads 10 Views