Detection of Transition Times from Singleparticletracking Trajectories
Abstract
In heterogeneous environments, the diffusivity is not constant but changes with time. It is important to detect changes in the diffusivity from singleparticletracking trajectories in experiments. Here, we devise a novel method for detecting the transition times of the diffusivity from trajectory data. A key idea of this method is the introduction of a characteristic time scale of the diffusive states, which is obtained by a fluctuation analysis of the timeaveraged mean square displacements. We test our method in silico by using the Langevin equation with a fluctuating diffusivity. We show that our method can successfully detect the transition times of diffusive states and obtain the diffusion coefficient as a function of time. This method will provide a quantitative description of the fluctuating diffusivity in heterogeneous environments and can be applied to time series with transitions of states.
The mean square displacement (MSD) is one of the most popular observables for quantifying the diffusivity. In Brownian motion, the MSD increases linearly with time, and the diffusivity can be quantified by the slope of the MSD, i.e., the diffusion coefficient. The diffusion coefficient is determined by the surrounding environment including the viscosity of the medium and the properties of the diffusing particle, e.g., the shape of the Brownian particle. When there are no fluctuations in the properties of the surrounding environment and the diffusing particle, no intrinsic differences arise between the diffusivities for shorttime and longtime measurements except for fluctuations of the diffusivity due to the finite measurement times.
In heterogeneous environments such as amorphous materials and living cells, diffusion often becomes anomalous; that is, the MSD does not increase linearly with time Scher and Montroll (1975); Golding and Cox (2006); Manzo et al. (2015). The local diffusivities in these environments are highly heterogeneous. These heterogeneities are sometimes static or fluctuating. For example, the charge transport in amorphous materials Scher and Montroll (1975) as well as the diffusion of proteins on DNA Granéli et al. (2006); Wang et al. (2006) can be modeled by a quenched trap model, where a random walker jumps in static random energy landscape Bouchaud and Georges (1990). In other words, the characteristic time scale of a change in the energy landscape is much longer than that of random walkers. On the other hand, in supercooled liquids, mobile and immoblie particles are distributed in space, and the diffusive properties (mobile and immobile properties) will change with time, i.e., dynamic heterogeneity Yamamoto and Onuki (1998a, b); Richert (2002). Moreover, transmembrane proteins Sergé et al. (2008); Weron et al. (2017) and membranebound proteins on biological membranes Yamamoto et al. (2017) exhibit a temporally heterogeneous diffusivity. In these systems, the diffusivity for shorttime measurements is intrinsically different from that for longtime measurements.
One of the most important issues in heterogeneous environments is to uncover the local diffusivity from singleparticle trajectories. However, there is a crucial difficulty in extracting the local diffusivity in both spatially and temporally heterogeneous environments. In particular, one cannot know the boundaries of regions with the same diffusivities and transition times when the diffusive states change in spatially and temporally heterogeneous environments, respectively. In previous studies, maximum likelihood estimators were proposed to determine the dynamic changes of the diffusivity Montiel et al. (2006); Koo and Mochrie (2016), where the key idea is to detect the transition times when the diffusivity changes drastically. However, an empirical parameter is necessary to implement the method.
The detection of the transition times is also important in statetransition processes, e.g., channel gating Wang et al. (2016), the conformational transition of proteins Chung et al. (2012), the rotation of F1ATPase Noji et al. (1997), the fluorescence of quantum dots Stefani et al. (2009), and nanopore sensing of single molecules Howorka and Siwy (2009). A method for detecting transition times using only trajectories without prior knowledge and empirical parameters is desired in timeseries analysis.
Here, we devise an estimation method for characterizing the shorttime diffusivity from trajectory data without knowing the transition times of the diffusive states. In our method, there are no parameters that are determined empirically. Thus, our method can be applied when many singleparticletracking trajectories are obtained. We show that our method can successfully detect the transition times of the diffusivity and estimate the local diffusivity in the (overdamped) Langevin equation with a fluctuating diffusion coefficient.
We assume that there are many trajectories for the same system and that the system can be described by the Langevin equation with a fluctuating diffusivity (LEFD):
(1) 
where is the position of a particle at time , is dimensional white Gaussian noise with and , and is the diffusion coefficient at time , which is a stochastic process independent of . Although we do not assume any condition on , i.e., the diffusion coefficient may be nonMarkov and depend on the position , we assume that the variance of is sufficiently large. In particular, it is much greater than the variance of the diffusion coefficients obtained by the timeaveraged MSD defined by Eq. (2) when the measurement time is the same as the characteristic time scale of the diffusive state.
In our setting, we do not know

the number of diffusive states and

the time scales of the diffusive states.
This is because we do not know the transition times when a diffusive state changes in singleparticletracking trajectories. This is one of the most difficult issues when estimating the fluctuating diffusivity. To overcome this difficulty, we apply a fluctuation analysis of the timeaveraged MSDs to obtain a characteristic time scale of the diffusive states. The timeaveraged MSD is defined as
(2) 
To characterize the fluctuations in the timeaveraged MSDs, we use the relative standard deviation (RSD) of the timeaveraged MSDs, defined as
(3) 
as a function of the measurement time ( is fixed). This type of quantity is widely used to investigate the ergodic property He et al. (2008); Deng and Barkai (2009) as well as the characteristic time of the system Akimoto et al. (2011); Miyaguchi and Akimoto (2011); Uneyama et al. (2012, 2015); Miyaguchi et al. (2016). In fact, the RSD analysis provides the longest relaxation time in the reptation model, which is a model of entangled polymers Uneyama et al. (2012, 2015).
When is a stationary stochastic process, i.e., the characteristic time of the stochastic process is finite, the general formula for the RSD is derived as Uneyama et al. (2015)
(4) 
where is the normalized correlation function of the diffusion coefficient, i.e., . Thus, if the relaxation time of the system is (roughly speaking, the correlation function decays as ), the asymptotic form of the RSD becomes
(5) 
Therefore, is obtained by the crossover time from the plateau to the decay in the RSD. In particular, when the correlation function decays exponentially, the crossover time in the RSD is given by . From many singleparticletracking trajectories, one can calculate the timeaveraged MSDs. Taking the ensemble average of the timeaveraged MSDs gives us the RSD. In this way, one can obtain the characteristic time scale of from singleparticletracking trajectories.
Here, we devise a novel method to detect the changes in states from a singleparticletracking trajectory. First, we define the timeaveraged diffusion coefficient (TDC) at time by
(6) 
There are two parameters, and , in the TDC. We set as the minimal time step of the trajectory; thus, it is not necessary to tune this parameter. On the other hand, we have to tune the parameter by introducing a tuning parameter as . Since is of the same order as the system’s characteristic time, can be smaller than one. In what follows, we use .
Second, using the effective diffusion coefficient , which is obtained by the ensemble average of the timeaveraged MSD, i.e., , we define the crossing points as the points at which the TDC crosses , i.e., and or and , satisfying , where is the time step of the trajectories (see Fig. LABEL:LEFD_two_statesA). Note that the crossing points are not exact points representing changes in the diffusive states because different diffusive states coexist in a time window of . Therefore, we define the transition time as . The term is not exact when the threshold is not at the middle of two successive diffusive states. If only one transition occurs in the time interval , which is a physically reasonable assumption, the transition times represent the points of changes in the diffusive states. Note that some transition times of the diffusive states will be still missing.
To correct the transition times obtained above, we test whether successive diffusive states are significantly different. Since we know the transition times of the diffusive states, we can estimate the diffusion coefficient in the time interval : the diffusion coefficient of the th diffusive state is given by
(7) 
Since we consider a situation that is sufficiently large (), fluctuations of can be approximated as a Gaussian distribution by the central limit theorem. According to a statistical test, the th and th states can be considered as the same state if there exists such that both the and states satisfy
(8) 
where is the variance of the TDC with the time window and the diffusion coefficient , which is given by , and is determined by the level of statistical significance, e.g., when the value is 0.05. Therefore, the transition times can be corrected if the two successive diffusion states are the same. We repeat this procedure: Eq. (7) will be calculated again after correcting the transition times , and the above test will be repeated to correct the transition times.
Furthermore, one can improve the transition times by changing the thresholds around the transition times. The detailed procedure and flowchart of our method are given in the Supplemental Material SM ().
Here, we test our method with the trajectories of three different LEFD models, where the number of diffusive states is two, three, and uncountable. The crossover times in the RSD are finite for all models.
In the Langevin equation with the twostate diffusivity, Fig. 1B shows the diffusion coefficient obtained by our method. Almost all diffusive states can be classified into two states according to the condition (8) with . Moreover, the deviations in the transition times from the actual transition times are within 0.25. Thus, we successfully extract the underlying diffusion process from a single trajectory after obtaining the characteristic time scale of the diffusive states.
We introduce different relaxation times in the two sojourntime distributions (we use the exponential distribution for both sojourntime distributions) and examine the effects of the tuning parameter. Figure 2 shows the TDCs for different tuning parameters , , and 0.01, corresponding to , 1.6, and 0.16, respectively. As clearly seen, when the tuning parameter is small, the fluctuations in the TDC become large. Therefore, inaccurate transition times may be detected when is too small. On the other hand, the actual transition times may not be detected when is too large. In fact, the transition times around cannot be detected in the case of the green dotted line (). As a result, the tuning parameter can be set to or between 0.1 and 0.01.
Next, we analyze the LEFD with the threestate diffusivity. The sojourntime distributions are exponential distributions, and their relaxation times are the same in each state (). For the threestate LEFD, one can obtain several diffusive states from a single trajectory by our method after calculating the crossover time in the RSD using many trajectories. Figure 3 shows the diffusion coefficient obtained by our method, where we revised the threshold using the procedures described in the Supplemental Material SM (). As shown in Fig. 3A, the transition times are correctly detected. Moreover, almost all diffusive states belong to the three diffusive states using , and the distribution of the estimated diffusion coefficients has three peaks corresponding to the exact diffusion coefficient (see Fig. 3B).
Finally, we apply our method to a diffusion process with an uncountable number of diffusive states. In particular, we use the annealed transit time model (ATTM) Massignan et al. (2014); Akimoto and Yamamoto (2016). The ATTM was proposed to describe heterogeneous diffusion in living cells Massignan et al. (2014); Manzo et al. (2015). The diffusion process is described by the LEFD where is coupled to the sojourn time. When the sojourn time is , the diffusion coefficient is given by . Here, we assume that the sojourntime distribution follows an exponential distribution . One can obtain by the RSD analysis Akimoto and Yamamoto (2016).
Figure 4A shows the diffusion coefficient obtained by our method. Because the variance of is not large, the transition times are not correctly detected compared with the other two models. However, the transition times for the highly diffusive states can be detected correctly. Moreover, Fig. 4B shows the relation between the obtained diffusion coefficient and the sojourn times, which exhibits a powerlaw relation . Therefore, our method can also be applied to systems with an uncountable number of diffusive states.
Diffusivity changes with time in temporally/spatially heterogeneous environments such as cells and supercooled liquids. It is difficult to estimate such a fluctuating diffusivity from singleparticle trajectories because one does not have information about the transition times when the diffusivity changes. In this paper, we have proposed a new method for detecting the transition times from single trajectories. Our method is based on a fluctuation analysis of the timeaveraged MSD to extract information on the characteristic time scale of the system. We have applied this method to three different diffusion processes, i.e., the LEFD with two states, the LEFD with three states, and the ATTM, which has an uncountable number of diffusive states. Our method successfully extracts the transition times of the diffusivities and estimates the fluctuating diffusion coefficients in the three models. Since our method can be conducted with singleparticle trajectories, the application will be useful and of importance in experiments. Furthermore, a slight modification of this method will be also applied to the time series of statetransition processes.
E.Y. was supported by an MEXT (Ministry of Education, Culture, Sports, Science and Technology) GrantinAid for the “Building of Consortia for the Development of Human Resources in Science and Technology.”
References
 Scher and Montroll (1975) H. Scher and E. W. Montroll, Phys. Rev. B 12, 2455 (1975).
 Golding and Cox (2006) I. Golding and E. C. Cox, Phys. Rev. Lett. 96, 098102 (2006).
 Manzo et al. (2015) C. Manzo, J. A. TorrenoPina, P. Massignan, G. J. Lapeyre, M. Lewenstein, and M. F. Garcia Parajo, Phys. Rev. X 5, 011021 (2015).
 Granéli et al. (2006) A. Granéli, C. C. Yeykal, R. B. Robertson, and E. C. Greene, Proc. Natl. Acad. Sci. USA 103, 1221 (2006).
 Wang et al. (2006) Y. M. Wang, R. H. Austin, and E. C. Cox, Phys. Rev. Lett. 97, 048302 (2006).
 Bouchaud and Georges (1990) J. Bouchaud and A. Georges, Phys. Rep. 195, 127 (1990).
 Yamamoto and Onuki (1998a) R. Yamamoto and A. Onuki, Phys. Rev. Lett. 81, 4915 (1998a).
 Yamamoto and Onuki (1998b) R. Yamamoto and A. Onuki, Phys. Rev. E 58, 3515 (1998b).
 Richert (2002) R. Richert, J. Phys. Cond. Matt. 14, R703 (2002).
 Sergé et al. (2008) A. Sergé, N. Bertaux, H. Rigneault, and D. Marguet, Nat. Methods 5, 687 (2008).
 Weron et al. (2017) A. Weron, K. Burnecki, E. J. Akin, L. Solé, M. Balcerek, M. M. Tamkun, and D. Krap, Sci. Rep. 7, 5404 (2017).
 Yamamoto et al. (2017) E. Yamamoto, T. Akimoto, A. C. Kalli, K. Yasuoka, and M. S. P. Sansom, Sci. Adv. 3, e1601871 (2017).
 Montiel et al. (2006) D. Montiel, H. Cang, and H. Yang, J. Phys. Chem. B 110, 19763 (2006).
 Koo and Mochrie (2016) P. K. Koo and S. G. J. Mochrie, Phys. Rev. E 94, 052412 (2016).
 Wang et al. (2016) S. Wang, R. Vafabakhsh, W. F. Borschel, T. Ha, and C. G. Nichols, Nat. Struct. Mol. Biol. 23, 31 (2016).
 Chung et al. (2012) H. S. Chung, K. McHale, J. M. Louis, and W. A. Eaton, Science 335, 981 (2012).
 Noji et al. (1997) H. Noji, R. Yasuda, M. Yoshida, and K. Kinosita Jr, Nature 386, 299 (1997).
 Stefani et al. (2009) F. D. Stefani, J. P. Hoogenboom, and E. Barkai, Phys. Today 62, 34 (2009).
 Howorka and Siwy (2009) S. Howorka and Z. Siwy, Chem. Soc. Rev. 38, 2360 (2009).
 He et al. (2008) Y. He, S. Burov, R. Metzler, and E. Barkai, Phys. Rev. Lett. 101, 058101 (2008).
 Deng and Barkai (2009) W. Deng and E. Barkai, Phys. Rev. E 79, 011112 (2009).
 Akimoto et al. (2011) T. Akimoto, E. Yamamoto, K. Yasuoka, Y. Hirano, and M. Yasui, Phys. Rev. Lett. 107, 178103 (2011).
 Miyaguchi and Akimoto (2011) T. Miyaguchi and T. Akimoto, Phys. Rev. E 83, 062101 (2011).
 Uneyama et al. (2012) T. Uneyama, T. Akimoto, and T. Miyaguchi, J. Chem. Phys. 137, 114903 (2012).
 Uneyama et al. (2015) T. Uneyama, T. Miyaguchi, and T. Akimoto, Phys. Rev. E 92, 032140 (2015).
 Miyaguchi et al. (2016) T. Miyaguchi, T. Akimoto, and E. Yamamoto, Phys. Rev. E 94, 012109 (2016).
 (27) See Supplementary Material for the detailed procedures and flowchart of our method.
 Massignan et al. (2014) P. Massignan, C. Manzo, J. A. TorrenoPina, M. F. GarcíaParajo, M. Lewenstein, and G. J. Lapeyre, Phys. Rev. Lett. 112, 150603 (2014).
 Akimoto and Yamamoto (2016) T. Akimoto and E. Yamamoto, J. Stat. Mech. 2016, 123201 (2016).