{'fraction_non_alphanumeric': 0.032977862246154926, 'fraction_numerical': 0.01967409284482455, 'mean_word_length': 4.455837173579109, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Estimation of diffusive states from single-particle trajectory in heterogeneous medium using machine-learning methods

We propose a novel approach to analyze random walks in heterogeneous medium using a hybrid machine-learning method based on a gamma mixture and a hidden Markov model. A gamma mixture and a hidden Markov model respectively provide the number and the most probable sequence of diffusive states from the time series position data of particles/molecules obtained by single-particle/molecule tracking (SPT/SMT) method. We evaluate the performance of our proposed method for numerically generated trajectories. It is shown that our proposed method can correctly extract the number of diffusive states when each trajectory is long enough to be frame averaged. We also indicate that our method can provide an indicator whether the assumption of a medium consisting of discrete diffusive states is appropriate or not based on the available amount of trajectory data. Then, we demonstrate an application of our method to the analysis of experimentally obtained SPT data.


Introduction
Single-particle/molecule tracking (SPT/SMT) techniques are widely used to investigate intercellular kinetics and biophysical processes at cell membranes [1][2][3][4][5][6][7][8][9] and extract non-bulk information on soft-materials. [10][11][12] As a quantity of interest for SPT/SMT trajectory data, mean-squared displacement (MSD) is widely used. The shape of an MSD indicates the modes of the particle/ molecule motion: normal, anomalous, directed, and confined diffusion. 1,7,13 For example, an MSD for a trajectory of a particle/ molecule in a homogeneous medium is a linear function of time, and the slope of the MSD indicates the diffusion coefficient. Though the transient behavior of MSD contains further information about particle/molecule motion, the MSD analysis is not convenient for practical use due to scattering of MSD plot induced by the limitation of available trajectory data. Thus, MSD is usually adopted to extract the diffusion coefficient of particles/molecule moving in a homogeneous medium. To analyze SPT/SMT trajectories containing adsorption/desorption motion and diffusion anisotropy, several analysis methods have been developed. [14][15][16][17][18][19][20] In the past few years, new approaches based on Bayesian and machine-learning methods have been developed. Metzner et al. 21 developed a Bayesian inference method based on a timediscrete Ornstein-Uhlenbeck process. Their method is specialized for analyzing a Brownian motion with drift and can extract time-dependent statistical parameters. Ott et al. 22 applied a hidden Markov model (HMM) to an analysis of two different modes of diffusion: a fast and a slow diffusive state in heterogeneous medium. The drawback of their method is that the number of diffusive states has to be fixed a priori in their maximum likelihood approach to an HMM analysis. Persson et al. 23 proposed an analysis method, the variational Bayes SPT (vbSPT), based on a variational Bayesian treatment of an HMM. [24][25][26] Their vbSPT can estimate the number of diffusive states and state transition rates from trajectories in a heterogeneous medium. Their vbSPT is a powerful analysis method when a sufficient amount of trajectory data is available. However, it is often pointed out that a proper choice of a prior distribution is difficult and a large amount of data is required for a reliable estimation of parameters in a variational Bayesian treatment. For example, more than 3000 trajectories (one trajectory consists of steps less than 20) were analyzed for the reliable estimation of the number of diffusive states in the vbSPT. 23 Unfortunately, there remains many situations that limited number of samples is available [27][28][29][30] due to photobleaching of the dyes, limited imaging area, and so on. Moreover, since these studies 22,23 provided the discussion in the dimensional form, the users have to adjust parameters to their problem.
In this study, we extend these HMM approaches by introducing a gamma mixture model (GMM). 31 In other word, we propose a hybrid machine learning method based on a GMM and an HMM. Our method is simple and can be easily implemented in a postprocessing program for SPT/SMT. The number of diffusive states is estimated by the GMM with the expectationmaximization (EM) algorithm 32 because the GMM is a more suitable classification algorithm than an HMM. An HMM originally estimate the most probable path of hidden states, where the number of the states are given as a prior knowledge. 32 By combining GMM and HMM, we can estimate most probable sequence of diffusive states from the trajectory data by HMM based on the number of the states estimated by the GMM in our method. We apply our proposed method to a numerically generated Brownian motion to investigate the performance. Since the previous studies such as ref. 20, 22 and 23 and our present method assumed that a medium consists of discrete diffusive states, we also discuss the statistical validity of the assumption from the given number of data points of trajectories. Then, the experimentally obtained SPT data is analyzed as a demonstration.

Numerical model of Brownian motion
We assume that a particle is a single point; that is, the position of a maker and the center of diffusion are indistinguishable. Since the positions of the maker and the diffusion center cannot be distinguished in an ordinary optical microscope due to the diffraction limit, 1,3,4 this assumption is reasonable in most SPT/SMT experiments. We consider the trajectory of a particle in a medium consisting of K media of the diffusion coefficients D m (m = 1, 2, 3,. . .,K), where D m may be expressed as a function of position or time. The diffusion coefficients as functions of position and time correspond to a spatially heterogeneous medium such as a surface having adsorption sites 22,[33][34][35] and to a temporally heterogeneous medium whose property is temporally controlled, 36,37 respectively. We also assume that the particle isotropically diffuses in each medium and the effect of interfaces between the media on the particle motion is negligible. The trajectory is generated by the Wiener process described by the overdamped Langevin equation as follows where r(t) is the position vector of the particle at time t. The variables where the r and t are nondimensionalized quantities, and D R,m is the ratio of D m to D 1 (the subscript R indicates ratio). Eqn (2) is numerically solved by the finite difference method expressed by where t j is the time of j-th time step, and Dt is the discretized time (i.e., t j = jDt). The dimensional time step tDt corresponds to the time interval between each frame of SPT/SMT experiment. In this study, t = 1 to simplify the discussion. To obtain particle trajectories having a few diffusive states, D R,m is varied in the calculation of eqn (3). We treat D R,m as a function of time varying with given state transition probability, because it is difficult to obtain the trajectories under the controlled condition with treating D R,m as a function of position. Since the squared displacement is the only feature and does not contain the position information of the particle, this treatment of D R,m is reasonable in our method. Fig. 1 shows the sample trajectory under the conditions of t = 1.0 Â 10 À2 , D R,1 = 1.0, and D R,2 = 5.0. The total number of time steps is 50. The locations of a particle are classified by color according to the diffusive states.

Gamma mixture model
After experimentally/numerically obtaining trajectory data, a squared displacement d(t j ) at the time of j-th time step is calculated as Considering ensemble average of the squared displacement, one can obtain mean squared displacement (MSD). We use only the squared displacement d(t j ) as the feature in our proposed method; thus, our method can be applied to multi-dimensional trajectories without any modification. Since the probability distribution of r(t j ) À r(t jÀ1 ) follows a Gaussian distribution (see eqn (3)), the probability distribution of d(t j ) follows a gamma distribution. Then, the number of diffusive states and its diffusion coefficients are estimated from the squared displacement d(t j ) using a univariate GMM. 31 We apply the EM algorithm to find maximum likelihood solutions for models having unobserved latent variables. Here, we introduce a GMM following the introduction procedure of a well-known mixture distribution, a Gaussian mixture model. 32 The univariate gamma distribution with a shape parameter a and an inverse scale parameter b is written as where G(x) is the gamma function. In the gamma mixture distribution, the probability distribution p(x) is expressed as a linear superposition of gamma distributions as follows: where p k are the mixing coefficients with 0 r p k r 1 and P K k¼1 p k ¼ 1. The number of gamma distributions having different parameters is represented as K, which also corresponds to the number of diffusive states in our method. Here, we introduce latent variables of a K-dimensional binary (0 or 1) This means that the probability of an observation x comes from z k is represented by p(z k = 1) = p k . Since z is a binary random variable, the distribution p(z k = 1) can be written in the form The conditional distribution of x at given z is also expressed as When we consider a data set X = (x 1 , x 2 ,. . .,x N ) T consisting of N observations, the corresponding latent variables can be written by an N Â K matrix of Z in which with n-th row is given by z T n , where the superscript T means transpose. Using eqn (7) and (8), the likelihood for the data set {X,Z} is expressed as where a = (a 1 , a 2 ,. . ., and z nk is n, k element of Z. By maximizing the likelihood, we can find the best model to fit the data set. The log likelihood is written as The direct maximization of eqn (10) is difficult, because the Z is unobserved variables. Therefore, we consider the maximization problem of the expectation of eqn (10) with respect to the posterior distribution of the latent variables Z. The posterior distribution p(Z|X,a,b,p) is obtained from eqn (7) and (8) with Bayes' theorem as Then the expectation of z nk with respect to the posterior distribution is calculated as Using eqn (12), the expectation of the log likelihood is expressed as As shown in eqn (13), the expectation of the log likelihood consists of the two terms, g(z nk )ln p k and g(z nk ){ln G(x|a k ,b k )}. First, we maximize the expectation of the log likelihood with respect to the mixing coefficients p k . Since there is the constraint condition P K k¼1 p k ¼ 1, we introduce a Lagrange multiplier l and consider the following Lagrange function L(p,l): Then, taking the derivative of L(p,l) with respect to p k equal to 0, we can obtain Multiplying eqn (15) by p k and summing over k with the eqn (15) can be written as To maximize the term g(z nk ){ln G(x|a k ,b k )} in eqn (13), the partial derivatives of it with respect to a k and b k are set to be 0 as follows, Eqn (18) can be easily solved as Using eqn (19), a k is expressed as where c is the digamma function. Since eqn (20) (20), (19) and (16), respectively. The maximum likelihood parameters can be obtained by the iterative calculation of the E and M steps. The number of gamma distributions having different parameters or the number of diffusive states K is extracted based on the Bayesian information criterion (BIC) represented as where m is the number of estimated parameters. We adopt K to minimize BIC in eqn (21) as the number of diffusive states to avoid the over-fitting of GMM to the data. The ratios of diffusion coefficients, D R , m , can be calculated as the mean values of each gamma distribution. When the EM algorithm converge to a local solution or some of the mixing coefficients are degenerate (p k { 1), we seek the maximum likelihood solution by randomly resetting initial conditions of the parameters.

Hidden Markov model
After extracting the number of diffusion states, an HMM is applied to the sequential data of the squared displacements d(t) to extract the most probable path of the diffusive sates. In this study, we use an HMM based on the maximum likelihood approach. Latent variables are also considered in the HMM. An observation d(t j ) is generated through an emission probability (conditional distribution of the observation from a specific state) 32 from the corresponding latent variable as in the case of the GMM. The latent variable transfers from one state to another following the transition probability. The transition and emission probabilities are estimated by the Baum-Welch algorithm. 32 Then, using the estimated transition and emission probabilities, the most likely path is calculated by the Viterbi algorithm. 32 The detail of these algorithm is described in text books such as ref. 32

Results and discussions
We validate our proposed method using numerically generated trajectories (see eqn (3)) in which the number and sequence of diffusive states are given and known in advance. Then, we apply our method to experimentally obtained SPT data as a demonstration.

Classification performance of GMM for trajectories with localization error
We first compare our method with the existing methods. 20 where 2L + 1 is the number of averaged frames and 2L + 1 of the squared displacements are averaged. After calculating and averaging the squared displacement for each trajectory, we estimate the number of diffusive states using all trajectories. It is assumed that the trajectories are independent and identically distributed. In this study, we calculate the correctness of estimation for the number of diffusive states calculated from 50 independent trials. The correctness is 1.0 for L = 1 and L = 2. These L are equivalent to the value adopted in ref. 20. However, the correctness is lower (o0.1) for L = 0 (without a frame average). These results indicate that our method shows better performance when each trajectory is long enough to be frame averaged. This can be achieved under the conditions of relatively small Dt or long diffusion-state lifetime. 11 We investigate the detailed performance of our method based on the analysis of the trajectories consisting of two diffusive states. Hereafter, we consider the trajectories without the localization error of the particles. Since the localization error will increase the apparent diffusion coefficient of the particles and obscure the true value, we compare the estimated diffusion coefficient with the exact one without considering the localization error. The trajectories are generated under the conditions of D R,1 = 1.0, D R,2 = 4.0, and Dt = 4.0 Â 10 À2 . The particles transfer between the two diffusive states with the transition probability of 0.05 or stay in the same states with the probability of 0.95 during the time interval of Dt. We prepare 1000 trajectories with consecutive 30 frames (time steps). These conditions are similar to those discussed in the previous section and the literature 20,22,23 in which the conditions were determined as a model of membrane or protein kinetics. We extract the number of the diffusive states, their most probable sequence, and the transition probability from the trajectories using our method. The squared displacements, d(t), defined by eqn (4) are calculated from the generated trajectory data. Then, the number of diffusive states and each diffusion coefficient are estimated from the probability distribution of the squared displacement using the GMM. We consider a frame average of the squared displacements defined as eqn (22) to reduce the scattering of the squared displacements and improve the estimation by the GMM. The choice of a certain value of L means that the probability of the state transition is sufficiently small in the frames of 2L + 1, and L should be less than 10 in the experimental conditions reported in the literature. 20,22,23 By substituting % d(t j ) for x j , g(z nk ) is calculated by eqn (12) as the E step. Then, as the M step, a new , b new , and p new are calculated from eqn (20), (19) and (16), respectively. In the GMM, the maximum likelihood parameters are obtained by the iterative calculation of the E and M steps. Fig. 2a shows the typical result of the GMM for the trajectories, where consecutive 5 (L = 2) steps of the squared displacements are averaged. The histogram of the squared displacement is also shown in the figure for the comparison. Note that the histogram is not used in the GMM calculation and the bin width does not affect the GMM result. Two peaks corresponding to the two diffusive states have been detected using the GMM. The ratios of diffusion coefficients are estimated as D R,1 = 1.03 and D R,2 = 3.68 in this example. In the HMM calculation, the initial estimation of transfer probability between each state is determined based on the number of frames for averaging (5 in this example). The most likely path for each trajectory is calculated by the Viterbi algorithm. Fig. 2b shows a typical example of an estimation for the most likely path. The correctness of the estimation is defined as C HMM ¼ number of correctly estimated states number of the squared displacements : In this example, C HMM = 0.966, and the reliability of the estimation is reduced near the transition from state 1 to state 2 (time step of 14 in Fig. 2b). The transition probability is estimated as 0.033 from state 1 to 2, where the exact probability is 0.05. This example shows that our proposed method is effective to extract the diffusive states and estimate the transfer between the states.
In this study, we calculate C GMM from 50 independent trials. The tested numbers of averaging frames are from 3, 5, and 7 (L = 1, 2, 3). The correctness, C GMM , drastically increases with increasing number of averaging frames or time steps as shown in Fig. 3a. The correctness, C GMM , is 0 for L = 0, when the time steps of each trajectory is from 15 to 100. This result shows that approaches the assigned value for large numbers of time steps. Though D es R,1 is larger than the assigned value, the deviation of D es R,1 from the assigned value is small (B3%).
In the above discussion, the number of trajectories is fixed at 1000. In the following, we fix the time steps of each trajectory of 30 and vary the number of trajectories. The other conditions remain the same (i.e., D R,1 = 1.0, D R,2 = 4.0, Dt = 4.0 Â 10 À2 , L = 2). Fig. 4a shows the correctness of estimation by the GMM, C GMM . The correctness C GMM monotonically increases and is larger than 0.9 when the number of trajectories is larger than 50. The estimated D R,1 and D R,2 shown in Fig. 4b indicate that reliable estimation is achieved by our method even when the number of trajectories is small. For example, the D R,1 and D R,2 can be estimated with 10% when the number of trajectories larger than 25.  Fig. 5 shows the correctness C GMM , where L = 1, 2, 3. It is considered that the correct estimation will be difficult for small D R,2 . The correctness also drastically increases with increasing D R,2 for L = 2, and the reliable estimation is achieved for L = 3 even when D R,2 is small (C GMM = 0.74 for D R,2 = 2.0). The correctness for L = 1 is low even when D R,2 Z 4. Thus, a larger number of trajectories or time steps of each trajectory is required for reliable estimation when L = 1 (see also Fig. 3).    Fig. 6 shows the correctness, C HMM , and the estimated transfer probability from state 1 to 2. The correctness, C HMM , gradually increases with increasing the time steps of each trajectory. However, C HMM is smaller than C GMM . Though there are a lot of trajectories without state transitions in the given time steps due to small transition probability and number of time steps, the HMM always classifies the data into two states; thus, C HMM is small and the transition probability is over estimated for the small numbers of time steps.

Three diffusive states
We discuss the application of our method to the analysis of the trajectories for particles traveling in a medium having three diffusive states. We generate the trajectories under the following conditions: three diffusive states with D R,1 = 1.0 (state 1), D R,2 = 4.0 (state 2), and D R,3 = 12.0 (state 3) and the discretized time Dt = 3.0 Â 10 À3 . The characteristic time and length are t = 1 s and 0.5 mm, respectively. The particles transfer each state during the time interval of Dt under the following probabilities: the transition probability from state 1 to state 2 is 0.031 and to state 3 is 0.1. The transition probability from state 2 to state 1 is 0.017 and to state 3 is 0.01. The transition probability from state 3 to state 1 is 0.055 and to state 3 is 0.055. These conditions correspond to those discussed in Persson et al.,23 where the parameters were given in dimensional quantities: D 1 = 0.25 mm 2 s À1 , D 2 = 1.0 mm 2 s À1 , D 3 = 3.0 mm 2 s À1 , and Dt = 3.0 Â 10 À3 s.
We estimate the number of diffusive states, where the time steps of each trajectory are fixed at 30. The correctness for the estimation of number of diffusive states is calculated from 50 independent trials as shown in Fig. 7a. The correctness rapidly increases in 100 trajectories and reaches to unity over 750 trajectories. The correctness for this 3 states estimation is smaller than that for 2 states estimation shown in Fig. 4a because the number of data in each state is smaller by construction. Fig. 7b shows the estimated ratio of diffusion coefficients. The estimated D R,1 converges to the input value with 1000 trajectories. On the other hand, the estimated D R,2 and D R,3 are 10% or 20% lager than the input values because the variances of squared displacement generated from D R,2 and D R,3 are larger than that from D R,1 . Our method works well in the 3-states problem.

Distinction of discrete diffusive states
In the above discussions, we have considered the particle trajectories in a medium having discrete diffusive states. The medium is assumed as consisting of several discrete diffusive states because HMM can treat only discrete hidden states. 32 This is one of the limitation of the HMM. It is sometimes convenient to regard a medium consisting of multiple discrete  diffusive states or particle motion is approximately expressed as such, when the underlying process is unclear before the analysis. 28,[40][41][42] Therefore, it is worth clarifying how much amount of data is necessary for treating a medium consisting of discrete diffusive states. Fig. 3a, 4a, 5 and 7a can be understood as showing the required amount of data to determine whether a medium consists of discrete diffusive states or not by the GMM based on BIC.
As a comparison to a medium having several discrete diffusive states, we consider a case where a particle travels in a medium having a continuously varying diffusion coefficient. One thousand trajectories are generated by eqn (3), where D R,m is replaced by the time varying diffusion coefficient D R (t j ). The diffusion coefficient linearly varies from 1 to D max as D R (t j ) = (D max À 1)( j À 1)/(N À 1) + 1, where D max = 4.0 and the time steps N = 30. This condition corresponds to that discussed in Section 4.2.2 except for a continuously varying diffusion coefficient, and the correctness is unity when the number of frame averaging is over 5 (L Z 2) in Fig. 3a. Fig. 8 shows the probability densities for the squared displacements of the trajectories with the continuously varying diffusion coefficient D R (t j ) and the discrete diffusion states (D R,1 = 1.0, and D R,2 = 4.0) for comparison, where the number of frame averaging is 5 (L = 2). The solid and dashed curves are obtained by the GMM in which the number of the mixing distribution is determined by BIC. The GMM returns a single state with the diffusion coefficient of 2.45 corresponding to the time average of D R (t j ) for the trajectory with the continuously varying diffusion coefficient (solid line). For the trajectory traveling in the discrete states, the discrete diffusive states are successfully classified by the GMM (dashed line). As shown in figures from 3 to 7, the GMM can classify the discrete diffusive states when the sufficient number of time steps is available and/or the number of averaging frames is sufficiently large for the trajectories with the discrete diffusive states. By contrast, for the trajectory in the continuously varying diffusion coefficient, the estimated number of diffusive states is determined as one by our method even when the large number of times steps and/or averaging frames are used. Our method provides a statistical validity of the assumption that a medium can be treated as consisting of discrete diffusive states for the given number of data points of trajectory.

Analysis of experimentally obtained SPT data
Now, we apply our proposed method to the analysis of experimentally obtained SPT data. The SPT data was obtained by an optical microscope in the same manner as that of our previous study. 43 We used ZnS-AgInS 2 nanoparticles 44-46 as probe particles, and mixed it into a polydimethylsiloxane (PDMS) layer. An oil immersion objective lens with 100 times magnification with N.A. of 1.4 was used. The time step (frame interval) was Dt = 0.2 s. The SPT images were captured after 20 h adding a curing agent in the PDMS layer at a room temperature of 20 1C. The SPT images were analyzed using IDL-code. 47 For the completely cured PDMS layer, the diffusion coefficient was calculated as O(10 À4 ) mm 2 s À1 . This indicates that the localization error is considered as O(10 2 ) nm. Under this condition, there were some diffusive states in the PDMS layer. 43 Fig. 9a shows the estimation result of the number of diffusive states using the GMM. Fig. 9b shows the single particle trajectory having 185 location data points with the most probable diffusive states. In the GMM analysis, it was found that there were three diffusive states of diffusion coefficients of 4 Â 10 À3 , 8 Â 10 À3 , and 12 Â 10 À3 mm 2 s À1 . These values reasonably agreed with those obtained by MSD analysis. 43 The particle  motions having relatively large displacement can be captured as the fast state (large diffusion coefficient). Though there is the place around x = 0.25, y = 0.12 containing two diffusive states due to the estimation error induced from the small amount of the data, each diffusive state occupies different place in the layer as shown in Fig. 9b. It is considered that the PDMS layer consists of several portions having different viscosity, when we assume discrete diffusive states. We conclude that our proposed method is applicable to the analysis of experimental SPT/SMT data.

Conclusions
In this article, we have proposed a hybrid method of a gamma mixture model (GMM) and a hidden Markov model (HMM) to classify particle trajectory for a particle moving in a heterogeneous medium. We introduce a GMM as an extension of a Gaussian mixture model based on the expectation-maximization (EM) algorithm. The number of diffusive states is estimated by the GMM, and then the HMM is used to extract the most likely path of the diffusive states. The correct estimation of the number of diffusive states can be achieved from small amount of trajectory data by considering frame average of squared displacements. The transition path of diffusive states is estimated by the Viterbi algorithm based on the estimated number of diffusive states. We compare our method with existing methods by calculating trajectory for a particle moving in a medium having two diffusive states. It is shown that our proposed method can extract the number of diffusive states more reliably than existing methods when the number of averaging frames is large. Thus, our method is a powerful tool for the trajectories obtained with relatively large frame rate or having long diffusion-state lifetime. Furthermore, we also indicate that our method can provide an indicator whether the assumption of a medium consisting of discrete diffusive states is appropriate or not based on the amount of the given data. Our hybrid method of the GMM and HMM is promising method for analyzing single-particle/ molecule tracking data when limited number of data is available.

Conflicts of interest
There are no conflicts to declare.