{'fraction_non_alphanumeric': 0.04334448994865185, 'fraction_numerical': 0.02218156871712783, 'mean_word_length': 4.461398728428701, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Covariance Distributions in Single Particle Tracking

Several recent experiments, including our own in the fission yeast, S. pombe, have characterized the motions of gene loci within living nuclei by measuring the locus position over time, then proceeding to obtain the statistical properties of this motion. To address the question of whether a population of single particle tracks, obtained from many different cells, corresponds to a single mode of diffusion, we derive theoretical equations describing the probability distribution of the displacement covariance, assuming the displacement is a zero-mean multivariate Gaussian random variable. We also determine the corresponding theoretical means, variances, and third central moments. Bolstering the theory is good agreement between its predictions and the results obtained for various simulated and measured data sets, including simulated particle trajectories of simple and anomalous diffusion, and the measured trajectories of an optically-trapped bead in water, and in a viscoelastic solution. We also show that, for sufficiently long tracks, each covariance distribution in these examples is well-described by a skew-normal distribution with mean, variance, and skewness given by theory. For experimental S. pombe gene locus data, however, we find that the first two covariance distributions are wider than predicted, although the third and subsequent covariances are well-described by theory. This suggests that the origin of the theory-experiment discrepancy is associated with localization noise, which influences only the first two covariances. Thus, we hypothesize that the discrepancy is caused by locus-to-locus heterogeneity in the localization noise. Further simulations reveal excess covariance widths can be largely recreated on the basis of heterogeneous noise. We conclude that the motion of gene loci in fission yeast is consistent with a single mode of diffusion.


I. INTRODUCTION
Single particle tracking has long been been applied to elucidate the dynamics of various soft-matter and biological systems [1][2][3][4][5]. Recent advances in fluorescent tagging and imaging now also enable tracking-based studies of single molecules or moieties within living cells. The motion of fluorescently labeled particles is most often analyzed by determining the mean squared displacement (MSD) as a function of time delay between observations. This approach has been applied to a wide variety of macro-and supra-molecular complexes inside cells from diverse organisms [6][7][8][9][10][11][12][13][14][15][16][17]. Examples in mammalian cells include reports of simple diffusion for a transmembrane protein [7] as well as multiple diffusive states within subpopulations (immobilized, subdiffusive and diffusive) for a viral protein [8] or several sub-populations (immobilized, slow and fast diffusion) for a DNA-binding protein [10]. For the DNA itself, subdiffusion has been reported for DNA loci in bacteria [18] and subdiffusion [9] or confined diffusion [6,19] in eukaryotes. * simon.mochrie@yale.edu The question of how to optimally determine the diffusivity from the time-averaged MSD (taMSD) has been investigated thoroughly for the case of simple diffusion [15,16,[20][21][22][23]. Typically, the slope of the taMSD for a limited number of time delays is used to yield an estimate of the diffusivity for each track and distribution of diffusivities across multiple tracks. However, within a track, the displacements measured for different time delays are not statistically independent. This, together with the particle position measurements being imperfect, as they include static localization noise (the uncertainty in particle position due to a limited number of detected photons) [24][25][26] and motion blur (the spatial spread of the signal due to camera integrating particle's positions over exposure time) [25], confounds MSD analysis. Recently, Vestergaard et al. [27] have shown that the optimal way to gauge diffusivity, while accounting for localization noise and motion blur, is to use an estimator based on the covariances of the particle displacements.
Another fundamental question is whether a single or multiple diffusion coefficients are appropriate to describe data collected from heterogeneous biological systems. Several approaches have been developed to sort particle trajectories into different diffusive states [11,15,17,28,29]. However, the simplest way to investigate whether data realize a single diffusivity or not is to compare the width of the measured diffusivity distribution to the width one would expect for a single diffusivity.
Vestergaard et al. [27] derived the expected distribution and variance of diffusivities measured from tracks of finite length in the case of simple diffusion, assuming that the static localization error is Gaussian-distributed. However, this analysis is not applicable to the more general case of anomalous diffusion. In the case of simple diffusion analyzed in Ref. [27], there are only two non-zero covariances, linearly related to the diffusivity. By contrast, there are as many non-zero covariances as there are steps in a track in the case of anomalous diffusion. A biologically important case of anomalous diffusion is the motion of a chromosomal locus. While the mean covariances for anomalous diffusion have been analyzed [30,31], to-date we are unaware of any comparison between the measured and predicted covariance distributions in the context of gene loci motion, nor any consideration of whether the motion of gene loci is homogeneous in time or whether a gene locus may undergo transitions among different modes of diffusion. The goal of this paper is to answer the question: Does the in vivo motion of a gene locus in fission yeast follow a single mode of diffusion, or not?
The paper is organized as follows. In Sec. II A and Sec. II B, we present the probability distribution, the mean, the variance, and the third central moment of the displacement covariance matrix elements for arbitrary modes of diffusion in one and two dimensions (1D and 2D), respectively. In Sec. II C, we specialize to give the diffusivity distribution, its variance, and third central moment in the case of simple 2D diffusion, reproducing some of the results of Ref. [27] via a different route. In Sec. III, we compare theoretical covariance distributions for a zero-mean Gaussian random process, expressed in terms of the mean covariance matrix elements, to a number of simulated and experimentally measured covariances. For simulated simple diffusion (Sec. III A) we find excellent agreement with theory. We also apply our methods to optical tweezers (OT) data (Sec. III B), and again find good agreement with theory. The agreement with theory in this case underscores that the usual procedures for fitting OT measurements should be generally modified to account for the motion blur, as pointed out previously [32]. Finally, we analyze motion of gene loci in living S. pombe yeast cells (Sec. III C), and find an important discrepancy with theory. Our analysis points to the discrepancy stemming from experimental locus-tolocus heterogeneity in localization noise, while the theory assumes the localization error to be drawn from a single normal distribution. After accounting for this heterogeneity with simulations, through implementing heterogeneous localization noise, we conclude that the gene loci in living S. pombe indeed undergo a single mode of diffusion. Finally, in Sec. IV, we summarize and conclude.

A. One-dimensional analysis
We consider data that consist of a collection of single particle tracks, each characterized by the same welldefined diffusive properties. There are at least three important reasons to then consider tracks of finite length. The first is that over time, fluorescently-labeled proteins in vivo may experience changes in their diffusive properties -changes in their diffusive state -such as those caused by binding and unbinding events, for example, to DNA or other complexes in the cell. The second is that, even for proteins that always remain in a single diffusive state, individual fluorescent labels can blink and eventually bleach. Third, in microscopy experiments, proteins can readily diffuse out of the focal volume. All of these processes give rise to tracks that realize a single diffusive state for a finite duration. Therefore, we start off by considering a population of tracks each comprising N + 1 particle coordinates along the x-axis, {x j }, with j = 1 through N , corresponding to N particle displacements along the x-axis, {∆x j = x j+1 −x j }, and we focus on the covariance of these one-dimensional displacements. Each measured track provides an estimate of the covariance matrix elements, which, given that the diffusive properties do not vary within a track, should depend only on time separation n for a stationary process. Each measured track provides an estimate of the covariance matrix elements via: where ∆x j is the displacement of step n and N is the total number of steps in the track. Because of their inherent stochasticity, different tracks yield different values for S n . However, averaging over many tracks yields the underlying mean covariance matrix: where the angular brackets indicate an average over tracks, and we have defined Σ n = S n . In general, the covariance matrix is a symmetric Toeplitz matrix.
EQ. 1 permits us to calculate S n for each individual track, and thus, to determine the experimental distribution S n from a population of experimental or simulated tracks. We can then test our understanding of the underlying diffusive process by comparing these empirical distributions to corresponding theoretical predictions.
To determine theoretical expressions for the distributions of S n , our starting point is the hypothesis that the probability of observing a particular N -step track is given by a multivariate Gaussian distribution, characterized by the N × N covariance matrix Σ (EQ. 2): where ∆x = (∆x 1 , ∆x 2 , ...∆x N ) T is the vector of N successive particle displacements along the x-direction within the track, and |Σ| is the determinant of the co-variance matrix. To proceed, we introduce the (Toeplitz) matrices, [C 0 ] jk = 2 N I and [C n ] jk = 1 (N −n) δ j k±n , which permits us to re-write EQ. 1 as and for n > 0. Using EQ. 3, EQ. 4, and EQ. 5, and the Fourier transform representation of the Dirac delta function, we may express the probability distribution of the random variables, S n , for a given covariance matrix, in terms of a certain matrix determinant: where |Σ −1 +iωC n | and |I+iωΣC n | are the determinants of Σ −1 + iωC n and I + iωΣC n , respectively. EQ. 6 is the probability density of S n , given the covariance martrix Σ. After calculating the determinant as a function of ω, the integral over ω may be performed numerically to obtain p(S n |Σ).
In Appendix A, we show from Eq. 6 that the first three moments of S n are and Combining EQ. 7 and EQ. 8, we find that the variance of S n is Combining EQ. 7, EQ. 8, and EQ. 9, we find that the third central moment of S n is B. Two-dimensional analysis Single-particle tracking often results in twodimensional data, so that two vectors are available for each track, namely ∆x and ∆y. In this case, given that the motion along each dimension is independent of the other, the likelihood of observing a particular trajectory is where Σ x and Σ y are the covariance matrices for x-and y-displacements, respectively. If we redefine S n to corre-spond to two dimensions via then, correspondingly, we can introduce a twodimensional covariance matrix, Σ, via Commonly, the diffusive behavior is isotropic, in which case we furthermore have that In this isotropic case, it follows that the probability density for the two-dimensional covariance is In Appendix B, we calculate S n , S 2 n , and S 3 n in this two-dimensional case. Using these results, it is straightforward to show that the mean, variance, and third central moment of the S n in 2D are and

C. Simple 2D diffusion
Many experimental systems can be expected to realize simple 2D diffusion with experimental errors, corresponding to a tridiagonal covariance matrix, where the only non-zero covariance matrix elements are Σ 0 and Σ 1 , which are related to the diffusion coefficient, D, the static localization noise, σ 2 , the time between camera exposures, ∆t, and the exposure time, ∆t E , via [27] and The static localization noise, σ 2 , is the error in particle localization that results from counting a limited number of photons. The terms involving the exposure time, ∆t E , correspond to motion blur, because the particle position is integrated while the shutter is open. In this case, 1 4 Tr(ΣC n ) 2 and 1 4 Tr(ΣC n ) 3 can be evaluated in closed form with the results that for n > 1, for n ≥ 4. EQs. 22, 23, and 24 reproduce EQ. 9 of Ref. [27], which employs a different method to find the variance of S n . EQs. 20 and 21 can be used to estimate D from each individual track. For a given N -displacement track we have: Because of the physical importance of the diffusion coefficient, we also consider its distribution, mean, variance and third moment for a distribution of D values estimated using EQ. 30.
It follows that the the probability density, the variance, and the third central moment of the diffusion coefficient are respectively.
We show in Appendix E, that EQ. 31 reproduces EQ. C17 of Ref. [27]. EQ. 32 is exact; it reproduces approximated solution EQ. 17 of Ref. [27] to second order in 1 N .

III. COMPARISONS OF THEORY TO SIMULATIONS AND EXPERIMENTS
A. Simulations of simple 2D diffusion First, we compare the theory of Sec. II to simulated particle trajectories undergoing simple 2D diffusion, generated as described in Ref. [29]. We simulated 10 6 steps, and set D = 0.0055 µm 2 s −1 , σ 2 = 7.94 × 10 −5 µm 2 , and ∆t E ∆t = 1. We partitioned these displacements into tracks of varied length, N , ranging from 19 to 419 steps, and then calculated the covariance matrix elements and diffusion coefficient (EQ. 30) for each such track.
Normalized histograms of thus calculated S 0 , S 1 , S 2 , and S 3 are plotted in Fig. 1(a), (b), (c), and (d), respectively, for different track lengths N = 19, 39, and 79. The S 0 -distributions are entirely positive, as required. The mean of the S 1 -distributions are positive, as expected when motion blur dominates static localization error. For n > 1, the means of the S n -distributions appear to be zero, also as expected. The distributions become progressively narrower with increasing track length, because a longer track length represents a more accurate measurement. Also shown in Fig. 1 as the solid red lines are the corresponding theoretical distributions (EQ. 16), calculated using only the experimental mean covariance matrix elements. Clearly, these theoretical predictions closely match the simulated distributions. Logarithmiclinear plots of the same data and model are presented in Fig. 2, demonstrating that the simulated and theoretical distributions continue to agree well, even in the far tails.
The Central Limit Theorem informs us that the distributions must each approach a Gaussian in the limit of large N . However, the simulated and theoretical distributions for N = 19, and even for N = 39, are noticeably skewed. To provide a simple way to empirically gauge the mean, variance, and third central moment of measured covariance distributions by fitting, we approximate them with a skew normal distribution [33]: (1−2δ 2 /π) 3 2 (ρ 2 (1 − 2α 2 π (1 + α))) 3 2 . The black lines in Figs. 1 and 2 are the best fits of a skew normal distribution to distributions sampled from the simulations, varying α, ζ, and ρ as fitting parameters. Near the peak, it is evident that the skew normal fit is able to accurately capture the shape of each distribution. While the bestfit skew normal distribution shows small deviations from both theory and simulated data (Fig. 2) in the far tails, we judge that it provides a good approximation to both.  Distribution of diffusion coefficients, D, for simulated simple 2D diffusion with D = 0.0055 µm 2 s −1 , σ 2 = 7.94 × 10 −5 µm 2 , and ∆t = ∆tE = 0.058s (same data as shown in FIGs. 1-2), for track lengths of (a) 29 steps, (b) 79 steps, and (c) 319 steps, shown using a linear y-axis. The red curve represents the theoretical prediction given by EQ. 31. The black curve corresponds to the best skew normal fit to the simulated distribution.
The distributions of estimated diffusion coefficients are shown in Figs. 3 and 4, together with skew-normal fits. Again close examination reinforces that the theoretical curves very closely match the simulated distributions, and that the skew-normal fits provide an excellent description with only small deviations in the tails. Fig. 5 shows how estimates of the covariance distributions from skew-normal fits depend on the inverse length of the tracks. Estimated mean (Fig. 5a), variance (Fig. 5b) and 3rd moment (Fig. 5d) for S 0 , S 1 , S 2 , and S 3 are compared to the ground-truth mean and to theoretical predictions for the variance and 3rd moment (EQs.  Fig. 5(c), we show specifically the variance of S 0 + 2S 1 . This combination removes static localization noise, leaving a result proportional to the diffusion coefficient. In every case, the best fit parameters match closely the corresponding theoretical values, supporting the utility of the skew normal function as a simple route to describe covariance and diffusivity distributions.

17-19). In
We further evaluate whether a skew normal distribution provides a good approximation to the data by applying the Kolmogorov-Smirnov (KS) test, which uses the difference between the CDFs of two distributions to estimate whether they originate from the same underlying distribution [34]. We performed analysis with track lengths N = 4 − 419. A plot of the KS statis-tic versus inverse track length for fits of S 0 and S 1 covariance distributions on simple 2D diffusion simulations, with D = 0.1389 µm 2 s −1 , σ 2 = 0.0015 µm 2 , and ∆t = ∆t E = 0.01s, and D = 0.0062 µm 2 s −1 , σ 2 = 0.0013 µm 2 , and ∆t = ∆t E = 0.01s, as well as the diffusivity parameters used in Figs. 3-5, is shown in Fig. 6(a). For smaller track lengths (i.e. N =19-49), the KS statistic is small, indicating a better fit, relative to other track lengths. Interestingly, very short track lengths (less than about 19) produce a poor fit using the skew normal distribution. This phenomenon can be explained by a limitation of the skew normal distribution, which is constrained to a maximum/minimum skewness of ±1. For the same three data sets, Fig. 6(b) shows the theoretically calculated skewness values for S 0 and S 1 versus inverse track length. For very small track lengths, the skewness can go beyond the range between ±1. In the case of D = 0.0062 µm 2 s −1 , for example, the skewness of S 1 is more negative than −1 at N = 12, rendering the skew normal distribution a poor descriptor of the data for this and smaller track lengths. To avoid this limitation, we use the exact theory to calculate the covariance distribution for very small track lengths, instead of fitting with the skew normal distribution. We also note that for larger track lengths, the KS statistic increases. While further investigation into this observation is needed, this suggests there may be a range of optimal track lengths to use while performing analysis.

Optically-trapped bead in water
Next, we compare the theory of Sec. II to measurements of confined diffusion of the x-coordinate of a 1 µmdiameter optically-trapped polystyrene bead suspended in water, from Ref. [35]. Specifically, ten independent time series of the bead position were recorded at 20 kHz, each containing 100,000 points, using a National Instruments DAQ PCIe-6343. In an optical trap, a particle experiences thermal Brownian motion, subject to a quadratic confining potential [36]. Theoretical results for the mean-square displacement (MSD) and covariance matrix elements of such a particle are presented in Appendix F, demonstrating that the covariance matrix elements are non-zero for n > 1, in contrast to the case of simple diffusion.
Normalized histograms of the experimental covariance distributions for S 0 , S 1 , S 2 , and S 3 are shown in Fig. 7, along with the corresponding theoretical curves (red lines, EQ. 18) and best fits to the skew normal distribution (black lines). This figure demonstrates good agreement between theory and experiment, and that the skew normal fit describes both well. The same plots are presented on logarithmic-linear axes in Fig. 8, again revealing that the skew normal fits show only small discrepancies in the tails of the distributions. Results from the  skew normal fits are compared to theory in Fig. 9, again showing good agreement. The mean covariance (red) and mean MSD (blue) of the bead determined from 12650 79-step tracks for the first twenty time delays (i.e., n = 1..20) is shown in Fig. 10. The MSD is equivalent, via Fourier transformation, to the power spectral density (PSD), which is usually fit to determine the stiffness of the optical trap. However, the fitting is typically carried out without considering the integration time of the optical detector. To assess accuracy of such an approach, we fit the measured MSD and covariance simultaneously by varying κ, τ , and σ 2 , along with ∆t E (i.e., considering integration time) or by fixing ∆t E = 0 (i.e., neglecting integration time). The results are shown as dashed (∆t E = 0) and solid (varied ∆t E ) lines in Fig. 10. Including an integration time provides a significantly improved description of the data as  judged by χ 2 values (4600 vs 23 for ∆t E = 0 and varied ∆t E , respectively) and by visual inspection (Fig. 10). Importantly, the best value of ∆t E is 0.048 ms, just slightly smaller than the time between successive measurements (0.050 ms). This analysis suggests that it is important to take into consideration a non-zero exposure time and the concomitant motion blur, when calibrating and analyzing optical tweezer data, as previously pointed out in Ref. [32]. Although EQ. F10 reveals that the exponential time dependence of the MSD is unchanged by motion blur, this circumstance is peculiar to this particular (exponential) form of the MSD [37]. In general, the shapes of the MSD and PSD are modified by motion blur (Appendix G).

Optically-trapped bead in a viscoelastic PEO solution
We also compared the theory of Sec. II to other data from Ref. [35], comprising 10 6 measurements of the x-coordinate of a 1 µm-diameter optically-trapped polystyrene bead suspended in a 1.5 wt% viscoelastic solution of long-chain (8 MDa) polyethylene oxide (PEO). Again, ten independent time series were recorded at 20 kHz, each containing 100,000 points, using a National Instruments DAQ PCIe-6343. As shown in Ref. [35], the PSD of a bead within this solution is well-described by a power law versus frequency, suggesting that a power law in time could be an appropriate basis for modeling the corresponding MSD and covariance, i.e. suggesting that the bead undergoes fractional Brownian motion (fBM), characterized by a subdiffusive exponent α < 1. In Appendix G, we derive theoretical expressions for the MSD and covariance for a bead undergoing fractional Brownian motion (fBM) with a subdiffusive exponent (α < 1), in the presence of motion blur and static localization noise. In this case too, the covariance is non-zero for n > 1.
Figs. 11 and 12 show the experimental covariance distributions, represented as histograms, overlaid with the theoretical distributions from EQ. 6 (red lines) and bestfit skew normal curves (black lines). Fig. 13 displays the skew normal fitting results versus track length, together with the correspponding theoretical predictions. For all of these figures, there is good agreement between theory and experiment.
In Fig. 14, we compare the measured covariance and MSD to the corresponding theoretically expected quantities (Appendix G). Plotted in this figure as circles and squares are the mean covariance and mean MSD, respectively, determined from 12650 79-step tracks, for the first twenty time delays. In this case too, we carried out global fits of the theoretically expected forms (Eq. G2, and Eq. G12 through G14), first, allowing for a non-zero ∆t E and varying α, ∆t E , D, and σ 2 , yielding best fits with χ 2 = 110, shown as the solid lines in the figure, and, second, setting ∆t E = 0, varying α, D, and σ 2 as fitting parameters, yielding best fits with χ 2 = 1.7 × 10 4 , shown as the dashed lines in the figure. Evidently, in this example too, including motion blur provides a superior description of the experimental measurements. Here, we also find that the best fit value of ∆t E is close to 0.05 ms.
C. Measurements and simulations of a chromosomal locus in fission yeast

Experimental measurements of a chromosomal locus
Next, we sought to apply the theory of Sec. II to the motion of a fluorescently-labeled genetic locus. Specifically, we examined the motion of a lacO/GFP-lacI tagged DNA sequence [38], adjacent to the mmf1 gene on chromosome II of fission yeast, S. pombe, as described in detail in Ref. [39]. Fluorescence and bright field images were acquired at 30C on a DeltaVision widefield microscope (Applied Precision/GE) equipped with a temperature control chamber, a 1.4 NA, ×100 objective (Olympus), solid-state-based illumination (Lumencor), and an Evolve 512 EMCCD camera (Photometrics). Fluorescence was excited at 488 nm and collected with emission filters passing 500-550 nm.
We analyzed 19 video microscopy movies, each containing 1000 images. Each image was acquired for an integration time, ∆t E = 10 ms and was separated from the next image by ∆t = 58 ms [39]. A total of 157871 time steps were included in the analysis. The position of each labeled locus was then tracked, as described in Ref. [11], and the resultant trajectories were partitioned into tracks of length N = 19, 29, 39, 49, 59, 69, 79, 89, 99, 109, 159, 319, and 419 steps, as above.
As also described in Ref. [39], the locus MSD is welldescribed by a model corresponding to fBm, chartacterized by an exponent α 0.44. In Fig. 15, we compare the measured covariance (blue circles) and MSD (red squares), determined from 5398 29-step tracks to the corresponding theoretically expected quantities (Appendix G), similar to the analysis done on the OT data in Sec. III B 2. We carried out global fits of the theoretically expected forms (Eq. G2 through G14), varying α, D, and σ 2 , yielding best fits with χ 2 = 166, shown as the solid lines in the figure. The resulting best-fit parameter values were α = 0.39 ± 5 × 10 −4 , D = 0.0061 ± 9 × 10 −4 µm 2 s −1 , and σ 2 = 0.0027 ± 0.08 µm 2 .
The measured covariance distributions, together with the corresponding theoretical predictions, based on the mean covariance, and the corresponding best-fits to a skew normal distribution are shown in Fig. 16. In contrast to the previous examples, the theoretical curves (red) for the S 0 distribution predict a significantly narrower distribution, than found experimentally (histograms). A similar trend, albeit less pronounced, can be discerned for the S 1 distributions. By contrast, the theoretical curves seem to accurately represent the experimental S 2 and S 3 distributions. The means and variances, determined from skew normal fits to these distri-butions, are plotted in Fig. 17(a) and 17(b), and provide further insight. As expected, the mean covariance is independent of track length, and the variance in track length increases with decreasing track length. However, in this case, the measured variance of S 0 is much larger than predicted, especially for long track lengths (Fig. 17(b) red circles and line). The measured variance of S 1 is also noticeably larger than predicted, but to a lesser extent. In contrast, the measured variances of S 2 and S 3 seem well-described by theory. It appears that the variances of both S 0 and S 1 display more-or-less constant offsets above their predicted values. To test this idea, we fit the experimental covariances of S 0 and S 1 versus inverse track length to the theoretical form plus a constant. The results of these fits are shown as the dashed lines in Fig. 17(b) and 17(c), and indeed describe the data well.
In general, localization noise contributes 2σ 2 to Σ 0 and −σ 2 to Σ 1 but does not contribute to Σ n for n > 1 (EQs. 20 and 21). This observation suggests that the discrepancy may originate from the localization noise in this experimental system. Given that localization noise has no contribution to the mean of Σ 0 + 2Σ 1 , we compared experimental and predicted variances of Σ 0 + 2Σ 1 to test whether the discrepancy with theory stems from the localization noise. Plotted in Fig. 17(c) versus inverse track length are the experimental variance of Σ 0 + 2Σ 1 , (diamonds) and the corresponding theoretical prediction (solid line). Although a discrepancy between experiment and theory remains, the actual value of the offset is about a factor of four smaller than the offset observed for the variance of S 0 alone, bolstering the idea that static localization noise is the culprit.
In experiments, each gene locus of the same site that is independently measured from cell-to-cell may show a different static localization noise, because of uneven illumination, different focal planes, different optical environments, etc. from other loci, all of which affect its localization noise. Therefore, we should expect locusto-locus heterogeneity in the localization noise, and thus some variation in the covariance matrix Σ underlying the observed locus motion. However, the theory presented in this paper is predicated on EQ. 3, which assumes that all samples are drawn from the same multivariate Gaussian distribution. To investigate the possible role of localization noise heterogeneity, we carried out additional simulations, described in Sec. III C 2.

Simulated fractional Brownian motion data with heterogeneous localization noise
To incorporate localization noise heterogeneity, we carried out simulations of fBm, for which the localization error varies from track to track, determined via two Gaussian random variables as follows. We consider that each track has its own unique localization noise described by a Gaussian (referred to as G1) with zero mean and a standard deviation determined randomly from another Gaussian distribution with the mean equal to the experimentally determined static localization noise σ 2 = 0.0032 µm 2 , and a standard deviation equal to σ σ 2 , which is a variable parameter of the model. Thus, each track has its own unique localization noise. The localization noise heterogeneity may be quantified by the value of σ σ 2 σ 2 , which we term the "localization heterogeneity".
For each of several different inhomogeneities, we simulated 1000 500-step tracks, for D = 0.0055 µm 2 s −0.44 , α = 0.44, σ 2 = 0.0032 µm 2 , and ∆t = ∆t E = 0.058s, comparable to the values determined from the experiment [39]. The resultant trajectories were partitioned as usual into tracks of length 19, 29, 39, 49,..... 319, 419, and the resultant covariance distributions fit by a skew normal function. Fig. 18 shows histograms of simulated S 0 , S 1 , S 2 and S 3 distributions, corresponding to σ σ 2 σ 2 = 0.20. Also shown in Fig. 18 are the corresponding theoretical predictions (red lines) and the skew normal best-fits (black lines), which describe the simulated distributions well in this case too. By contrast, the theoretical curves are significantly narrower, reminiscent of the disparity seen in Fig. 16 for the experimental gene locus data. Also similar to the gene locus data, the simulated and theoretical distributions of S 2 and S 3 agree.
The means, variances and 3rd moment of the covariances calculated from the best-fits for the simulated data are shown in Fig. 19a-d as a function inverse track length. The theory (EQ. 18) and the theory plus offset are shown as solid and dashed lines, respectively. The means reproduce what is expected from the given simulation parameters. Note, that this is true regardless of the values of the localization heterogeneity (data not shown). Strikingly reminiscent of the experimental results shown in Fig. 17, while the variance of S 0 lies above the theoretically predicted values by a relatively large offset, the variance of S 1 also lies above that predicted theoretically, but by a relatively small offset.
The fact that our simple simulations recapitulate important features of the experimental results on gene loci support the hypothesis that the observed theoryexperiment discrepancy is the result of inhomogeneous localization noise. To further test our model, we also determined the simulated and theoretical variance of Σ 0 + 2Σ 1 which are plotted together in Fig. 19(c) as the diamonds and line, respectively. In this case, however, simulation and theory match well, in contrast to the experimenttheory comparison which shows a residual discrepancy (Fig. 17(c)), indicating that our simple model does not perfectly model the experimental situation.
Next, we explore how localization heterogeneity levels affect the discrepancy between theory and the simulations. We vary the heterogeneity level, and, for each simulation, fit S 0 and S 1 by the theory plus offset and plot corresponding offsets vs. σ σ 2 σ 2 (Fig. 20). Fig. 20 plots the offsets in S 0 (red) and S 1 (blue) versus σ σ 2 σ 2 . Unsurprisingly, small locus-to-locus variation in the localization noise results in small offsets; and ramping up the localization noise increases offsets. For increasing inho-mogeneity, i.e. a broader initial Gaussian, G1, the offset increases. Also shown in Fig. 20 are the offsets for the experimental mmf1 gene locus data, where the red horizontal line is the S 0 offset, and the blue horizontal line is the S 1 offset. The lines derived from experiment intersect the simulated curves at comparable values of the heterogeneity, namely σ σ 2 σ 2 = 0.21 for Σ 0 and σ σ 2 σ 2 = 0.16 for Σ 1 . In view of the simplicity of our model, we regard this small discrepancy as satisfactory. Distributions alternative to Gaussians yielded similar results.

IV. CONCLUSION
We have developed and applied a simple method to determine whether a population of finite-length single particle tracks exhibits a single mode of diffusion. First, we derived theoretical equations that describe the distribution of displacement covariance matrix elements for a particle in an arbitrary, but well-defined, diffusive state. For both 1D and 2D, we determined the probability distribution of the covariance matrix elements, S n , given the covariance matrix, and then calculated the first three moments of S n , namely, the mean (trivially equal to Σ n ), variance, and third central moment. To test these theoretical results, we first simulated tracks undergoing simple 2D diffusion, uncovering excellent agreement between theoretical and simulated covariance matrix element distributions. We further showed that least-mean-squares fits of a skew normal function to the simulated distributions were able to accurately capture the shape of the simulated distributions. Best-fit parameters from skew normal fits, carried out for different track lengths also matched well to the theoretical expressions for the variance and third central moment versus track length. In addition to simulated tracks, we also followed this procedure for experimental data, first for an optically-trapped bead in water and then an optically-trapped bead in a viscoelastic high polymer solution. In these cases also, we found that the experimental covariance distributions are well described by the theoretical covariance distributions, using the measured mean covariances as the sole input. In these cases too, skew normal best fits yielded parameters that well described the mean, variance, and third central moment of the experimental covariance distributions versus inverse track length. These collected results give us confidence that the theory, presented in Sec. 2, is correct.
Finally, we analyzed the covariance distributions of the mmf1 gene locus in S. pombe. Here, the analysis revealed that the experimental S 0 and S 1 distributions are significantly wider than predicted on the basis of the mean covariance matrix elements, initially suggesting the presence of more than one diffusive state. However, the fact that the discrepancy is confined only to S 0 and S 1 led us to consider the hypothesis that the discrepancy originates with locus-to-locus heterogeneity in the static localization noise. To test this idea, we simulated fBm tracks with static localization noise of varying degrees of heterogeneity. Analyses of these tracks recreated a similar disparity between the theoretical and measured variance of S 0 and S 1 versus inverse track length. We explored varying the degree of the heterogeneity in the context of simulations in an attempt to more closely recapitulate the discrepancy seen in the gene locus data. While the heterogeneous noise we added largely accounts for the offsets we saw between the measured and theoretical variances of the covariance, it did not fully account for this behavior. Further investigation is needed to completely understand the factors that contribute to this discrepancy. In summary, we found that, for a single mode of diffusion, covariance distributions can be accurately predicted using the presented theory and mean covariances. Comparison of the predicted distributions with the experimentally measured ones enabled us to answer whether particles display a single or multiple diffusive states. Hence, this method is a simple yet powerful tool to determine whether a biological system exhibits a single diffusive state. To calculate the mean, variance, and skewness of S n in the one-dimensional case, we expand |I + iωΣC n | in a small-ω power series to third order in ω, i.e.
It follows that and First, we calculate the mean, S n . Using EQ. 6, we have that Integrating by parts, Eq. A5 becomes where the integral over S n now equals δ(ω). Using Eq. A2, it immediately follows that We calculate S 2 n similarly. Using EQ. 6, we have that Integrating by parts twice then yields where the integral over S n again equals δ(ω). Thus, we find Next, we calculate S 3 n .

(A13)
Appendix B: Moments of Sn (2D case) In the two-dimensional case, we expand |I + i 2 ωΣC n | to the third order, It follows that and Thus, Integrating by parts, Eq. B5 becomes where the integral over S n now equals δ(ω). Using Eq. B2, it follows that We calculate S 2 n similarly: Integrating by parts twice then yields where the integral over S n again equals δ(ω). Thus, we find We also calculate S 3 n : Integrating by parts three times then yields Finally, substituting EQ. B2, B3, B4 into EQ. B12, we find First, we'll focus on the 1-dimensional case. Now we need to expand |I + iωΣC n | to the third order It follows that and The mean value of S 3 n is Integrating by parts three times then yields Substituting EQ. C2, C3, C4 into EQ. C6, The third central moment µ 3 is The skewness (3rd standardized moment) µ 3 is Similarly, for 2-dimensional case, we next expand |I + i 2 ωΣC n | to the third order It follows that and The mean value of S 3 n is Integrating by parts three times then yields Substituting EQ. D3, D4, D5 into EQ. C15, The third central moment µ 3 is The skewness (3rd standardized moment) µ 3 is .

(C18)
Appendix D: Generalizing the calculation of the central moments to k-dimensional space The probability distribution of S n in k-dimensional space is Expand |I + i k ωΣC n | to the third order, It follows that and Using EQ. D6, D7 and D8, we obtain the first central moment the second central moment and the third central moment The skewness (3rd standardized moment) µ 3 is Because all tridiagonal Toeplitz matrices possess the same eigenvectors, the matrices, Σ, C 1 , I + i 2 ωΣC 0 , I + i 2 ωΣC 1 and I+ i 8∆t ω(ΣC 0 +2ΣC 1 ), are all diagonalized by the same orthogonal transformation, given by [U] jk = 2 N sin πjk N +1 . For Σ and C 1 , the eigenvalues are respectively, where k = 1, 2, .. N . It follows that and which is the 2D version of EQ. C17 of Ref. [27].
Appendix F: MSD and covariance matrix elements for a particle in a harmonic potential The usual approach for analyzing optical tweezers data is to calculate the power spectrum of the displacement fluctuations away from the potential minimum, which is equivalent to consideration of the mean-square displacement (MSD) versus time. To determine what new information can be gleaned by considering the covariance matrix elements of optical tweezers data, in this section, we calculate both the MSD and the covariance matrix elements for a particle in a viscous fluid, subject to a harmonic potential.
Our starting point is the theoretical result for the dis-placement correlation function: where κ is the trap stiffness, τ is the characteristic time of the trap, k B is Boltzmann's constant and T is the ab-solute temperature. Thus, for times t 0 and t n , separated by a time n∆t, we have (F2) It follows that the MSD, M n , after n time steps, each of ∆t, is while the covarience matrix elements are and However, data acquisition actually occurs for a time period ∆t E , while the spacing between successive acquisitions is ∆t. Thus, rather than measuring an instantaneous particle position, in general, experiments measure a motion-blurred position, averaged over the duration of data acquisition, i.e., acquisition period n measures where x(n∆t + t) is the particle position at time n∆t + t.
To incorporate motion blur into our calculations of the MSD and the covariance matrix elements, we must replace x 2 0 , etc. in EQs. F3 through F6 by the appropriate time-averaged quantities. Using and and incorporating static localization noise, the motion-blurred MSD becomes: Similarly, the covariance matrix elements become: and for n > 1.
In the limit that ∆t τ and ∆t E τ , using k B T κτ = D, where D is the particle's diffusion coefficient, EQs. F11, F12, and F13 reproduce the corresponding results for free diffusion with motion blur and static localization noise (EQs. 20 and 21), as expected. Generally, powerspectrum-based analyses of optical tweezers data ignore any possible motion blur [40]. Comparison between EQ. F3 and F10 suggests that in the case of exponential time-dependence motion blur changes the interpretation of the fluctuation amplitude without changing the shape of the MSD or PSD. However, this circumstance is special to the case of exponential relaxations and does not hold in general [37]. (See also, for example, Sec. III B 2) Appendix G: MSD and Covariance matrix elements for fractional Brownian motion (1D) In the case of fractional Brownian motion (fBm), our starting point is theoretical mean-square displacement from time t 1 to time t 2 [31], namely where α is the exponent characterizing the fBm, and where we refer to D as the diffusion coefficient, although its dimensions are m 2 s −α . Incorporating motion blur and static localization noise, the MSD becomes which reproduces the result given in Ref. [25]. The calculation for the fBm covariance terms can be performed as follows: Now break EQ. G5 into 3 integrals, and calculate each one. Starting with the left-most integral: Absolute value can be ignored since n∆t + ∆t + t − s > 0. The result is: −2(n∆t + ∆t) α+2 (α + 1)(α + 2) + (n∆t + ∆t − ∆t E ) α+2 (α + 1)(α + 2) + (n∆t + ∆t + ∆t E ) α+2 (α + 1)(α + 2) (G7) Next we have the middle integral: To evaluate this integral, we must take the cases when n = 1 and n > 1 into account. For n = 1, the integrand becomes |t − s| α : For n > 1, n∆t − ∆t + t − s > 0, so we can do the integral normally.
Lastly, the third part of the integral can be done without the absolute value, since n∆t > ∆t E : Thus, the covariance terms are: For n = 2, 3, 4, ...
Experimental MSD (circles) and covariance (squares) versus time delay for an optically-trapped bead in 1.5 wt% PEO solution, compared to corresponding best-fit model MSDs (blue solid and dashed lines) and covariances (red solid and dashed lines). The experimental MSD and covariance were determined from 12650 79-step tracks. Errors (standard errors of the mean) are smaller than the plotted points, and therefore no error bars are plotted. The MSD and the covariance were fit together in two ways. For the fits shown as the solid lines, four fitting parameters were used, namely α, D, σ 2 , and ∆tE, yielding best-fit parameter values of α = 0.20 ± 0.13, D = 3.3 ± 2.4 × 10 −5 µm 2 /s 0.2 , σ 2 = 2.8 ± 2.9 × 10 −7 µm 2 and ∆tE = 0.05 ± 0.03 ms, yielding χ 2 = 110. For the fits shown as dashed lines, ∆tE was set equal to zero, leaving three fitting parameters, α, D, and σ 2 . The best-fit parameter values in this case are α = 0.42, D = 8.7×10 −5 µm 2 /s 0.42 , and σ 2 = 2.8×10 −7 µm 2 , yielding χ 2 = 1.7 × 10 4 . Experimental MSD (circles) and covariance (squares) versus time delay for the gene locus mmf1 in S. pombe, compared to the corresponding best-fit model MSD (blue solid line) and covariance (red solid line). The experimental MSD and covariance were determined from 5398 29step tracks. Errors (standard errors of the mean) are smaller than the plotted points, and therefore no error bars are plotted. The MSD and the covariance were fit together, using three fitting parameters, namely α, D, and σ 2 . The resulting best-fit parameter values were α = 0.39 ± 5 × 10 −4 , D = 0.0061 ± 9 × 10 −4 µm 2 s −1 , and σ 2 = 0.0027 ± 0.08 µm 2 , yielding χ 2 = 166.  The mean of the covariances S0 (red circle), S1 (blue square), S2 (green triangle), and S3 (orange upside-down triangle) vs. inverse track length. Theory (straight lines) is calculated by EQ. 17. (b)Variance of the covariance vs. inverse track size of S0 (red circle), S1 (blue square), S2 (green triangle), and S3 (orange upside-down triangle). The theory is calculated by EQ. 18. Dashed lines correspond to a fit of the experimental data to the theoretical form plus a constant. (c) Variance of the covariance vs. inverse bin size of S0 + 2S1 (diamond). Calculating S0 + 2S1 removes static localization noise. Hence, the data variances are closer in value to the theory, compared the plot of S0. Theory (solid line) is calculated using Tr(Σ(C 0 +2C 1 )) 2

4
. Dashed lines correspond to a fit of the experimental data to the theoretical form plus a constant.    . 19. Dependence of the covariance matrix estimates on the track length for simulated particles undergoing fractional Brownian motion, with noise determined from a Gaussian distribution of widths for each track. The width of the noise distribution G1 is determined using homogeneity = σ σ 2 σ 2 = 0.20. (a) The mean of the covariances S0 (red circle), S1 (blue square), S2 (green triangle), and S3 (orange upside-down triangle) vs. inverse track length. (b)Variance of the covariance vs. inverse track size of S0 (red circle), S1 (blue square), S2 (green triangle), and S3 (orange upside-down triangle). The theory is calculated by EQ. 18. Dotted lines represent a linear fit plus a constant to the data points. (c) Variance of the covariance vs. inverse bin size of S0 + 2S1 (diamond). Theory (solid line) is calculated using Tr(Σ(C 0 +2C 1 )) 2 4 (d) The third central moment of the covariance Σ0 (circle), Σ1 (square), Σ2 (triangle), and Σ3 (upside-down triangle) vs. inverse track length. Theory (lines) is calculated by EQ. 19. Offsets between theoretical and fitted variances of simulated data with inhomogeneous localization noise for S0 (red circle) and S1 (blue square), versus heterogeneity of the localization used to determine the Gaussian noise distribution. Dashed lines represent a guide to the eye for each Sn offset. The horizontal red and blue solid lines are the S0 and S1 offsets for gene locus mmf1, respectively.