Inferring Latent States and Refining Force Estimates via Hierarchical Dirichlet Process Modeling in Single Particle Tracking Experiments

Understanding the basis for intracellular motion is critical as the field moves toward a deeper understanding of the relation between Brownian forces, molecular crowding, and anisotropic (or isotropic) energetic forcing. Effective forces and other parameters used to summarize molecular motion change over time in live cells due to latent state changes, e.g., changes induced by dynamic micro-environments, photobleaching, and other heterogeneity inherent in biological processes. This study discusses limitations in currently popular analysis methods (e.g., mean square displacement-based analyses) and how new techniques can be used to systematically analyze Single Particle Tracking (SPT) data experiencing abrupt state changes in time or space. The approach is to track GFP tagged chromatids in metaphase in live yeast cells and quantitatively probe the effective forces resulting from dynamic interactions that reflect the sum of a number of physical phenomena. State changes can be induced by various sources including: microtubule dynamics exerting force through the centromere, thermal polymer fluctuations, and DNA-based molecular machines including polymerases and protein exchange complexes such as chaperones and chromatin remodeling complexes. Simulations aiming to show the relevance of the approach to more general SPT data analyses are also studied. Refined force estimates are obtained by adopting and modifying a nonparametric Bayesian modeling technique, the Hierarchical Dirichlet Process Switching Linear Dynamical System (HDP-SLDS), for SPT applications. The HDP-SLDS method shows promise in systematically identifying dynamical regime changes induced by unobserved state changes when the number of underlying states is unknown in advance (a common problem in SPT applications). We expand on the relevance of the HDP-SLDS approach, review the relevant background of Hierarchical Dirichlet Processes, show how to map discrete time HDP-SLDS models to classic SPT models, and discuss limitations of the approach. In addition, we demonstrate new computational techniques for tuning hyperparameters and for checking the statistical consistency of model assumptions directly against individual experimental trajectories; the techniques circumvent the need for “ground-truth” and/or subjective information.

In this article, we demonstrate how the Hierarchical Dirichlet Process Switching Linear Dynamical System (HDP-SLDS) framework developed by Fox and co-workers [32] can be used to deduce the direction and magnitude of different forces that contribute to molecular motion in living cells [23]. The utility of combining the HDP-SLDS with SPT was motivated by experiments aiming to quantify the time varying forces driving chromosome dynamics. The approach presented shows promise in both (I) accelerating the scientific discovery process (i.e., statistically significant changes in dynamics can be reliably detected) and (II) automating preprocessing tasks required when analyzing and segmenting large SPT data sets.
The technique introduced is applicable to various scenarios where SPT trajectories are sampled frequently in time and particles can be accurately tracked over multiple frames, e.g. [15,16,23,26,28]. Extracting accurate and reliable force estimates from noisy position vs. time data in the aforementioned setting requires one to account for numerous complications inherent to experimental SPT data in living cells. For example, nonlinear and/or time changing systematic forces need to be differentiated from thermal fluctuations (i.e., random forces), both of which contribute to motion at the length and time scales measurable in living systems [23,24,31]. Furthermore, additional measurement noise (consisting of localization error amongst other factors [18,23,[33][34][35]) induced by the optical measurement apparatus must be systematically accounted for since this noise source varies substantially between and within single trajectories; inaccurate effective measurement noise estimates can appreciably influence estimates of kinetic parameters as well as statistical decisions about the underlying physical system [18,23,28,36]. Finally, extracting forces from position vs. time data requires one to explicitly or implicitly make numerous assumptions about the underlying effective dynamics. We believe these assumptions should be systematically tested directly against experimental data before one trusts kinetic quantities inferred from experimental data [23,37,38]. However, in live cell SPT studies, reference "ground-truth" is rarely available. Hence, techniques checking statistical assumptions directly against data are attractive (e.g., through goodness-of-fit hypothesis testing [23]).
The HDP-SLDS approach combines Hidden Markov Modeling (HMM), Kalman filtering, [39,40], and more recent ideas from Dirichlet Process modeling [41]. We demonstrate how the HDP-SLDS method can be used to reliably identify the time at which a state change occurs as well as the number of states implied by a specific time series. In the HDP-SLDS approach [32], the number of underlying states are inferred from the data via nonparametric Bayesian techniques [41,42]. Inferring the number of states jointly with the parameters determining the dynamics (i.e., in a single fully Bayesian computation) is useful because the number of underlying effective states is rarely known a priori in live cell SPT applications due to inherent heterogeneity between and within trajectories [21][22][23]. Since the HDP-SLDS framework directly infers the number of states from observed data, the user does not need to provide an accurate upper bound on the number of states or worry about a posteriori model selection issues [21,43,44]. This allows the HDP-SLDS approach to readily identify a wide range of distinct dynamical regimes which may occur within a single trajectory. Once the experimental trajectory is segmented into distinct kinetic states, one can use classic maximum likelihood estimation (MLE) techniques to infer kinetic parameters [23,45]. Carrying out MLE estimation after segmentation mitigates the bias induced by prior assumptions. It is emphasized throughout that the procedure described in this work can systematically account for measurement noise and complex spatial-temporal variations in thermal fluctuations & force/velocity fields.
In the particular application that we studied, the heuristics suggested for generating HDP-SLDS prior parameter input required modification before accurate results could be obtained in SPT applications. Specifically, prior specification heuristics employed by the original HDP-SLDS Bayesian analysis (using guidelines outlined in Ref. [32]) required substantial modification (prior specification substantially influences both state inference and molecular motion parameter estimation). In particular, the "hyperparameters" governing the measurement and thermal noise (parameterizing components of the "base-measure" [32,46]) required more careful calibration/tuning before high accuracy results could be obtained. In S1 Text, we discuss how ideas in Refs. [18,23] can be used for this type of hyperparameter tuning.
We note that the experimental applications are meant to serve as a proof-of-concept. The results presented are not intended to provide an exhaustive study of the forces involved in chromosome movement in metaphase. The problem is amenable to the HDP-SLDS approach since the GFP tagged experiments experienced abrupt state changes in the experimental trajectories before other downstream SPT computations could be reliably carried out (e.g., force and diffusion coefficient estimation). Not all "abrupt" were visually obvious to human observers, hence the approach is promising in both data-mining and preprocessing contexts. It should also be noted that two recent articles have discussed the utility of nonparametric Bayesian ideas in single-molecule data analysis [47,48]; Ref. [47] is similar in spirit to this manuscript except the focus is exclusively on simulation data mimicking features of SPT data and Ref. [48] focuses on using nonparametric Bayesian modeling to detect the number of states in HMM models. Note that in Ref. [48] each observation is assumed to be drawn independently from a distribution indexed by a latent state and the approach employed does not assume explicit temporal dependence between successive SPT time series observations whereas this work accounts for such temporal dependence (allowing quantitative approximations of spatially or temporally dependent effective forces and/or velocities).
This article is organized as follows: the mathematical models and experimental methods used are presented in Sec. 1; Results obtained when analyzing chromatid dynamics in yeast are presented in Sec. 2. Conclusions are presented in Sec. 3. We have also included Supporting Material (S1 and S5 Texts) where additional technical details are provided and control simulation studies are analyzed (results shown in S1-S7 Figs.). A companion publication [47] compares the HDP-SLDS approach to an open-source SPT Hidden Markov Model assuming a finite number of states through simulation studies.

Background and Models Considered
We begin by reviewing the main technical ideas underlying the HDP-SLDS introduced in Ref. [32]. Note that the SLDS models assumed allow researchers to decouple noise induced by thermal fluctuations, measurement noise, and spatially anisotropic velocity/forces. Forces are computed by leveraging the overdamped Langevin model [23] (this relationship is made explicit in Sec. 1.1.1 and S3 Text).
In the early SPT works, the spatial and temporal resolution afforded by the measurement device led researchers to focus mainly on Mean-Square-Displacement (MSD) type analyses to analyze singlemolecule data [49][50][51][52]. MSD approaches have many undesirable features, namely they tend to introduce unnecessary temporal averaging (i.e., they ignore the natural time ordering of the trajectory measurements) and they have a difficult time accounting for spatially varying forces (a common occurrence in live cells [23]). Advances in spatial and temporal resolution have inspired many researchers to develop new techniques aiming at more reliably extracting single-molecule level information out of measurements [19, 21, 22, 27, 36-38, 44, 53-55]. The previously cited works are most similar in spirit to the work presented, but all of the works encounter technical difficulties when the number of "states" is not known in advance. Additional technical complications arise when position estimates are obscured by non-negligible "measurement noise" [18,37] or force/velocity fields exhibiting temporal dependence [23].
The method of Fox et al. [32] overcomes the difficulties mentioned above by assuming that a discrete time series model of the form: can be used to describe the dynamics of each unique state. The position of the molecule or particle at time t i is denoted by the vector r i and the measured value of the position at this same time is denoted by ψ i (subscripts are used to index time); the position is not directly measurable due to "localization noise" and other artifacts induced by the experimental apparatus. The term "effective measurement noise" is meant to include the net measurement noise induced by finite photon counts, background fluorescence, motion blur, etc. [18,23,35,37]. Effective measurement noise is modeled as a mean zero normal random variable with covariance R; the expression ∼ N (0, R) conveys that the random vector, , is distributed according to the normal distribution N (0, R). Techniques for checking the validity of this modeling assumption via hypothesis testing technique are discussed in Ref. [23]. The term µ represents a (constant) "velocity vector" experienced by the particle multiplied by the observation time ∆t; the matrix F accounts for systematic spatial variations in forces or velocity; random thermal fluctuations are modeled by η. The "drift terms" [56,57], i.e. µ and F r i , can be used to quantify active (linear) force and velocity fields; spatial anisotropy in the force and velocity fields are accounted for by F r i . The physical assumptions we employ for inferring forces from position vs. time data are discussed in Sec. 1.1.1. The net parameter vector characterizing the (discrete) dynamics of a single state is given by θ = ( µ, F , R, Σ). Note that the HDP-SLDS presented in Ref. [32] assumed that all observations are uniformly spaced by ∆t time units.
Although the molecular position, r, is not directly observable, the discrete state-space model above allows one to readily plug directly into the established Kalman filtering equations [39,40] to infer the dynamics of r using the observations ψ. The Kalman filtering framework allows one to systematically treat measurement noise as well as spatially dependent particle forces and/or velocities (this spatial variation can induce classic "confinement" or "corralling" effects or be indicative of molecular binding in SPT data [23,36]). The importance of properly accounting for measurement noise is demonstrated in S2 Text and is also presented elsewhere [18,23,36,38].
An illustrative trajectory of the (observable) measurements ψ and (unobservable) r is shown in Fig. 1; this figure will also be used to also illustrate the "latent state" modeling discussed in the next paragraph. Note that the presence of F slightly complicates physical interpretation of the parameters [23], but the discrete model above can be readily mapped to a continuous time stochastic differential equation (SDE) studied in Ref. [23] (where diffusion coefficients, effective friction, and instantaneous force terms associated with an overdamped Langevin equation can be readily extracted). The equations mapping discrete time parameters used in the HDP-SLDS framework to continuous time SDE parameters are presented in S3 Text for the reader's convenience. The ability to map between discrete and continuous time models is important for both specifying physically inspired priors and for interpreting parameter estimates since classical quantities like diffusion coefficients and forces are typically defined via continuous time SDE models [17,18,23].
In the discrete time model above, θ contains the parameters (vectors and matrices) required to specify the stochastic dynamics defining the evolution equations of the system in a single "state". However, simple linear evolution equations are not expected to be valid for the entire duration of the trajectory in live cell measurements [23]. In this article, we assume that θ can change abruptly over time or space and hence the system's "state" can change over time. We will use the notation z i to denote the state at time i and will use θ z to denote the parameter vector characterizing the dynamics of r when it is in state z. The standard goal of an HMM inference procedure is to infer the state sequence {z 1 , z 2 , . . . z T } from an observation sequence { ψ 1 , ψ 2 , . . . ψ T }.
The typical HMM framework assumes temporal transitions between the states are governed by Markovian transition probabilities. If there are K fixed states, the transition vector associated with state i in the traditional HMM framework is prescribed by the vector π (i) = π K . Note that each component provides the probability of the state transition from state i to state j and the sum of components of this K dimensional vector is one (hence defines a proper discrete probability distribution); the collection { π (1) , π (2) , . . . π (K) } defines the classical HMM "transition matrix". In Bayesian inference of standard HMMs, the Dirichlet distribution is sometimes used as a prior for π (i) [21,55]. Additional details on the Dirichlet distribution and various infinite dimensional extensions, namely the Dirichlet Process and the Hierarchical Dirichlet Process, are provided in S4 Text.
The HDP-SLDS framework of Fox et al. comprehensively addresses issues not accounted for in Refs. [21-23, 27, 58]. The following modified quote from Ref. [32] captures the essence of the HDP-SLDS method: "The HDP-SLDS is an [infinite discrete state space] extension of hidden Markov models (HMMs) in which each HMM state, or mode, is associated with a linear dynamical process." As in Ref. [23], the statistical influence of thermal fluctuations and effective measurement noise on the dynamics is accounted for using the Kalman filter framework (this allows linear spatial variations in force and/or velocity); however, the HDP-SLDS also provides a mechanism for segmenting trajectories into chunks where dynamics are distinct. In addition, the HDP-SLDS framework advocated in Ref. [32] introduces a "sticky parameter" encouraging temporal state persistence; this feature can be advantageous in many SPT applications [47]. The ability to avoid model selection nuances [21,43,44] and infer the number of states in a data-driven fashion within a single Bayesian inference while also rigorously accounting for spatially or temporally varying stochastic dynamics are major advantages of the HDP-SLDS method.
Before discussing the major technical weaknesses of the HDP-SLDS method, we discuss some strengths/weaknesses of published SPT modeling approaches. Wavelet methods have the ability to identify sharp and abrupt changes in time without making too many assumptions about the underlying stochastic process [22], but they do not readily account for subtle spatial variations in the dynamics or noise (e.g., temporally varying velocity fields or diffusion coefficients) and can have difficulty in separating diffusive noise from measurement noise. The HDP-SLDS overcomes the aforementioned challenges by assuming a specific parametric model (however, if anomalous diffusive noise is deemed important to quantify system dynamics [22,26], this can be a problem for the HDP-SLDS approach). Ref. [58] is one of the pioneering efforts attempting to account for spatial and temporal variations in SPT signals, but the approach neglected to explicitly account for the effects of measurement noise (the approach also focused on comparative hypothesis tests and model selection as opposed to goodness-of-fit tests directly checking the consistency of a model's distributional assumptions to experimental data). The method in Ref. [27] attempts to account for spatial variations (allowing for correlated 2D forces) and the influence of measurement noise, but appeals to ad hoc statistical approximations which can adversely affect state and parameter inference [18]; these approximations (which can be avoided using Kalman filtering ideas [23]) can substantially complicate downstream analysis where one would like to check the assumptions against experimental data [36] (also Ref. [27] did not employ any formal hypothesis testing procedures). The "Windowed local MLE" method in Refs. [37,45] can account for spatial variations in force and measurement noise [23,28] and provides a procedures for goodness-of-fit testing, but does not prescribe a systematic and automated method for segmenting time series data into distinct dynamical segments.
Open problems facing the HDP-SLDS are associated with prior specification and model validation [46]. The priors used are selected primarily for computational convenience. For example, the so-called matrix normal inverse-Wishart (MNIW) and other priors proposed allow for exact Markov chain Monte Carlo (MCMC) sampling [32], but such priors have no real connection to the physical mechanisms governing molecular motion. The priors and their associated hyperparameters used by the HDP-SLDS method are discussed further in S5 Text. In the results, we demonstrate how a variant of the method in Ref. [23] can be used to correct for some artifacts introduced by "bad priors" (i.e., the parametric prior distributions assumed do not accurately reflect the true underlying process). Quantitatively prescribing a "good prior" is difficult because single-molecules trajectories contain a high degree of heterogeneity induced by the local micro-environment, conformational fluctuations of the biomolecules, etc. Hence finding an accurate prior representative of a single trajectory is non-trivial (finding a prior governing a population of trajectories in the spirit of Ref. [55] is even harder in SLDS modeling). Fortunately, the nonparametric Bayesian HDP-SLDS from Ref. [32] combined with frequentist ideas [23] provides researchers in SPT analysis a set of tools which can be used to more reliably and systematically extract information from complex in vivo measurements as we demonstrate in the Results and Supporting Text.

Inferring Instantaneous Forces from Noisy Position Measurements
In what follows, we present a technical discussion explicitly pointing out how we infer forces from a sequence of position measurements. Particle position is modeled using the overdamped limit of the Langevin equation (i.e., a first order SDE). The general form of the corresponding first order SDE is: where γ is a local effective "friction matrix" associated with the Langevin model and f (·) represents the force vector [23,27]; within the SLDS the force is a function of both the position r and the latent state.
The notation above was selected to facilitate comparison to Refs. [27] and [59]. The above model assumes particle inertia can be ignored (i.e., the dynamics are in the "overdamped" regime or the Smoluchowski limit of the full Langevin equation [59]) and local diffusivity is constant. In Text S3, the mapping between the continuous and discrete HDP-SLDS time series models is provided (note that in S3 Text, On the time and length scale of SPT measurements, it is generally accepted that particle inertia is unimportant. However, goodness-of-fit tests can be used to check if unobserved and unmodeled momentum dynamics affect position measurements [28,60]. If inertia is believed or empirically determined to be important, more complex (second order SDE) models for fitting forces from position data can be entertained [61]. The subtler assumption implicit in our estimate of the forces from the above model is that the Einstein relationship can be used to relate local effective friction to the local effective diffusion coefficient via the relationship: The above relationship is a classic microscopic relationship relating molecular diffusivity to local friction [59,62,63]. On longer time scales, the microscopic parameters become effective parameters [23,59] and the validity of the above relationship can become questionable due to phenomena like molecular crowding [64] or other factors (e.g., transient binding) inducing more complex dynamics. The temporal segmentation approach proposed aims to detect when a sudden event alters the dynamics in a statistically significant fashion, hence the odds of "mixing" different dynamical regimes is reduced (e.g., if the obstacle density experienced in one part of the cell is vastly different than another, then the effective dynamics in each regime will be different [64]). The HDP-SLDS's ability to identify regimes, within a single trajectory, where locally linear SDE models are consistent with the data makes the validity of Eq. 4 more plausible.
It should be noted that "effective force" data has been backed out from noisy position measurements in multiple in vitro single-molecule manipulation experiments [37,38,57]. In the aforementioned works, force measurements obtained via an external probe (e.g., AFM) could be compared to internal forces inferred from noisy position vs. time data (more precisely molecular extension vs. time) obtained by appealing to relationships similar to those in Eqns. 3 and 4. Despite the system's being far out of equilibrium in Refs. [37,38,57], appealing to the relationship in Eq. 4 yielded effective "internal forces" consistent with those measured by external probes. However, testing the validity of the above expression directly against in vivo SPT data is non-trivial since directly measuring fores in the crowded and heterogeneous environment of live cells is difficult. The formal statistical hypothesis testing procedure employed can only judge the consistency of the assumed drift and diffusion terms of the SDE in Eq. 3 against observational data; said differently, only the product γ −1 f ( r t ) affects the test (the validity of Eq. 4 does not influence the hypothesis test in any way). As we will discuss later, our force estimates obtained (without assuming or inputting a local effective cellular viscosity) are consistent with other in vivo force measurements [65]. Even if a researcher does not believe Eq. 4 should hold, our fitted models can still produce local effective velocity fields [59] which can aid in our quantitative understanding of intracellular dynamics.
Before proceeding, we explicitly mention a few more technical notes about our nominal force estimates. Our models make no assumptions about the existence of a gradient potential field [27] or the stationarity (temporal or spatial) of statistics characterizing the random process [59]. Transient structures and phenomena, which are often resolvable in live cell SPT measurements [23], can invalidate the validity of a force obtained by taking the gradient of an energy potential field and can also complicate nonparametric estimators [59] (nonparametric estimators often make implicit stationarity assumptions). Our models sample data on millisecond timescales, but allow for the local effective forces and diffusion to change over time and/or space (i.e., statistics are allowed to be non-stationary); we believe this feature will be critical in accurately characterizing heterogeneous high resolution optical microscopy data.

Experimental Methods
We examined the in vivo dynamics of chromatin during mitosis to determine the behavior of a region of the chromosome visualized through the binding of lac repressor fused to GFP (LacI-GFP) to lac operator (lacO). The lacO is a repeated array of operator DNA sequences (256 repeats, 10 kilobase prs.) integrated 6.8 kb from the centromere on chromosome XV. The spindle pole bodies (sites of microtubulare nucleation) are visualized through a fusion protein between a spindle pole component (Spc29) and RFP (red fluorescent protein). Cells were grown to logarithmic phase at 24 • C in rich media. Images were acquired on a Nikon Eclipse Ti wide-field inverted microscope with a 100x Apo TIRF 1.49 NA objective (Nikon, Melville, New York, USA) and Andor Clara CCD camera (Andor, South Windsor, Connecticut, USA) using continuous laser illumination. Images were streamed at the net effective camera acquisition rate of 22 frames/sec. The CCD camera's exposure time was set nominally to be the inverse of the net frame rate. Images were acquired at room temperature with Nikon NIS Elements imaging software (Nikon, Melville, New York, USA). The program "Speckle Tracker" and other methods outlined in Ref. [66] were used to estimate the centroid of GFP and RFP spots for position measurements.

Results
In Fig. 2 we display a histogram of the estimated instantaneous force magnitudes as well as the force magnitude trajectory computed from a time series of measured X/Y chromatid positions (the X/Y data and state estimates are shown in Fig. 3). The force was estimated by first determining the number of states implied by the observed trajectory and using the method described in S1 Text to segment the trajectory. For the experimental trajectories shown in this paper, the initial prior mean of the measurement noise covariance matrix, R, was assumed to be the identity matrix multiplied by a scalar, σ 2 , where σ = 40nm. This value was inspired by the fact that the effective measurement noise standard deviation was found to be in the 10 − 60nm range for these SPT trajectories (with a mode at 40nm). We specified K = 10 for the so-called "weak limit approximation" of the HDP-SLDS [67]; this term is discussed further in S5 Text and the relative insensitivity of the HDP-SLDS method to the K parameter is demonstrated in Fig S5. Additional parameters required to run the HDP-SLDS segmentation are reported in S5 Text. After the HDP-SLDS segmentation was obtained, the estimated MLE parameter vector of each unique state was used to obtain molecular position estimates via the Kalman filter [32,39]. Finally, the position estimates along with the MLE parameters were plugged into the equations shown in S3 Text to evaluate the instantaneous effective force (note: this provides a collection of 2D force vectors).
The force vector magnitudes reported in Fig. 2 observed are representative of other chromatid data sets studied using this analysis. However, some datasets exhibited more interesting "force regime changes" (as shown and discussed in the final experimental SPT trajectory studied). Two things should be emphasized in the analysis of this trajectory: (I) previously published works have reported that chromatids in metaphase-like conditions experience forces in the (relatively low) 0.1-0.2 pN range [65,68]. Our findings are consistent with these previous results, except the model considered here does not require an estimate of the effective viscosity (a quantity difficult to estimate in live cells) to infer force vectors. Local effective forces are estimated using time series analysis techniques applied to the so-called overdamped Langevin equation outlined in Sec. 1.1.1. (II) The HDP-SLDS method aided in accurately determining a transition between two states (physical interpretation of the identified states discussed below).
The top left panel of Fig. 3 displays the white light image of a trajectory obtained from the yeast chromatid SPT experiments. Note that we use the phrase "white light image" throughout to indicate a single diffraction limited image of the yeast cell obtained with white light illumination; the trajectory obtained using a laser frequency tuned to enhance GFP excitation was then overlaid upon this single image (all images are recorded at different times, but the yeast cell is not expected to move substantially during the experiment, so the white light illumination image gives one an idea of the spatial environment explored by the molecule). The top right panel displays the trajectory in X/Y space with time information color coded. The red vertical line in the bottom left panel of Fig. 3 shows the time at which a state change point was detected by the HDP-SLDS method. The bottom right panel of Fig. 3 shows that state estimates of the vbSPT [21] and HDP-SLDS methods both identified two states, though the latter captures the state persistence more accurately.
After the temporal segmentation was formally estimated using the HDP-SLDS, we computed the posterior mode of θ implied by the assumed HDP-SLDS model for each of the states in the two identified temporal segments. These two parameter vectors (one vector drawn for each of the two states) were used along with the experimental data to compute the Q test statistic [69] and the corresponding p−values. The null hypothesis of the one-sided goodness-of-fit test using the Q statistic assumes that the time series was generated by a Markovian SDE (evaluated at an "optimal" parameter estimate [69]) and the alternative hypothesis assumes that the time series was produced by any other process [23,60,69]. "Optimal" parameter estimates were obtained by two means: (i) drawing the parameter vector obtained after 10 4 iterations of the MCMC sampler in the HDP-SLDS inference; note that this approach included prior bias and (ii) via MLE estimation applied to HDP-SLDS segmented trajectories; with this approach, bias is mitigated since priors do not influence the parameter estimates.
In the first segment, using the observed data and the HDP-SLDS parameter vector estimate mentioned above, resulted in a p−value of 0.32 (little evidence for rejection), but in the second segment a p−value < 0.01 was computed (large evidence for rejection). Visual inspection of Data S1 (see rightmost spot) shows the effective measurement noise has abruptly changed, however the prior assumed by the HDP-SLDS segmentation only allowed the components of R to change modestly. The default settings of the MATLAB code associated with Ref. [32] assumes little dispersion about the nominally known mean R (all segments yielded an effective measurement noise having ≈ 40.0nm on all diagonal components via the HDP-SLDS method). Using a so-called "non-informative" or "uniform" prior can potentially remedy the situation, but such a prior can introduce both technical and computational problems in nonparametric Bayesian methods [42]; see S5 Text for further discussions on this issue.
Alternatively, if one uses the segmentation afforded by the HDP-SLDS, but instead computes the MLE using a variant of the technique of Ref. [23] (discussed in S1 Text) as opposed to using the posterior sample provided by the HDP-SLDS (which contains biases induced by priors), the p−values were found to be 0.31 and 0.45 for the left and right segments, respectively. The primary difference in the estimation results was in the noise components estimated, R and D, the effective measurement noise and diffusion coefficient matrices, respectively (see S3 Text). For example, the MLEsR = diag([19.5 2 , 39. s] were obtained in the right segment (i.e., the segment occurring later in time) using the model in Ref. [23]. Hats are used to denote MLE estimates and diag(·) denotes the square diagonal matrix formed by the arguments. Note that the diffusion coefficient estimate, D, is near machine single-precision zero in the right segment.
The GFP image stack associated with this video readily shows photobleaching occurs (see rightmost point spread function spot in S1 Data). The MLE predicted effectively zero diffusion for the right segment (i.e., the state appearing temporally after the first segment), quantitatively indicating that significant photobleaching had occurred and that the second state was mainly background photon noise not corresponding to an individual molecule. The quantitative evidence for this statement is the MLE parameter estimates (effective measurement noise near diffraction limit in the "right segment") and the corresponding p−values, reported in two paragraphs above, resulting from testing the fitted model against the data. In this experiment, the HDP-SLDS approach was able to automatically detect this measurement noise transition despite a poorly tuned prior; the state change was identified without requiring subject matter expertise or manual image inspection. Fig. 4 displays another experimental trajectory with a different type of state switching which was readily identified by the HDP-SLDS procedure. The posterior mode computed using Gibbs sampling of the HDP-SLDS model contained assignments with two and three states; in total 10 4 posterior MCMC samples were sampled. A majority of the MCMC draws implied two-state switching. However, one does not necessarily need to be concerned about the number of states if one is only interested in determine change points. The change points predicted by the HDP-SLDS are denoted by vertical lines (change points were determined to occur if a state change occurred at time i in over 25% of the posterior samples). Accurately identifying change points can help in quantifying the transition times between distinct physical states [21]. However the state switching in this particular system was likely caused by vibrations in the microscope's piezo stage. The HDP-SLDS can be used to automatically identify physically relevant state changes or state changes induced by experimental artifacts. As we show in Fig. S5, isolated outliers (e.g., caused by intense fluorescent background "flashes" in some SPT experiments) are also readily identified by the HDP-SLDS framework.
Next we turn to an example exhibiting subtle state changes (identified using the procedure outlined in S1 Text) in a more biologically relevant study. Fig. 5 shows two GFP tagged chromatid trajectories during the metaphase stage of mitosis. The white light image of these trajectories is shown in Fig. 6 where the spatial location of the RFP tagged mitotic spindle pole bodies (SPB) as well as the GFP tagged chromatid trajectories are shown. HDP-SLDS state segmentation was carried out for each of the trajectories; the state segmentation is denoted by vertical lines in the bottom panel of plots in Fig. 5. The procedure used to infer forces in Fig. 2 was carried out again for the pair of chromatids undergoing mitosis. The histogram of the force magnitudes is not atypical for trajectories observed, however the evolution of the force vectors was suggestive of a dynamic regime shift.
In order to probe this regime shift further, we computed the eigen-values and eigen-vectors of the MLE of B (a matrix associated with the continuous-time overdamped Langevin model; details provided in S3 Text). The eigen-analysis shown in Fig. 6 provides quantitative information about variations in the effective force vector as a function of position (i.e., the eigen-analysis provides both magnitude and directional information about a linear vector field). The upper-right panel shows that the largest magnitude eigen-value consistently decreased over time for each state identified by the HDP-SLDS segmentation. The rate of decrease was similar for the pair of sister chromatids. This phenomenon was originally identified using the crude time window segmentation of Ref. [23], but the HDP-SLDS provides a framework for more systematically dividing the trajectory into different dynamical segments.
In terms of physical interpretations, several molecular components can induce "force cross-talk" between sister chromatids. A potential explanation of the eigen-value phenomenon described above is that microtubules of the mitotic spindle connected to the common centromere of the sister chromatids are simultaneously inducing a common force change to both chromatids. Alternatively (or perhaps in addition), a common tension signal induced by chromatin loops and a network of pericentric proteins (cohesin and condensin [24,70]) might be influencing the forces experienced by the pair of sister chromatids; recall that the GFP dyes used are in close proximity to the centromere so molecular changes in the pericentric region can potentially influence the fluctuations observed in the tagged sister chromatids [24,70].
Regardless of the underlying biological cause of the phenomenon observed, the downstream eigenanalysis substantially benefited from the temporal segmentation provided by the HDP-SLDS framework. Using the segmentation afforded by the HDP-SLDS method resulted in the eigen-vectors shown in the bottom panel of Fig. 6; in this figure, the eigen-vectors associated with the largest magnitude eigen-value is displayed as a solid line and the eigen-vector associated with the smaller eigen-value of the estimated 2D B matrix is shown by a dashed line (the origin for the eigen-vectors correspond to the empirical average of the identified states). The (real-valued) eigen-values and eigen-vectors provide spatial information about restoring force directions experienced by the tagged chromatid. Observe how the pair of chromatids both exhibit a state change occurring around 16s. Note that the time of the change point occurring around 16s coincides for the sister chromatids, but each trajectory was processed independently by separate HDP-SLDS analyses. At 16s, the dominant eigen-vectors of the chromatids change from pointing towards the spindles to a direction nearly orthogonal to the plane connecting the two mitotic spindle poles. The eigen-values for both chromatids are decreasing at roughly the same rate (i.e., their effective "stiffness" is increasing). Also, recall that the pair of chromatids are subject to forces from multiple sources that can change in magnitude and direction; for example, tugging induced by the microtubules connected to the centromere shared by the chromatids, forces associated with nucleosome remodeling (e.g., RNA polymerase "tugging" on the strand as it makes RNA), as well as other spring like forces induced by the shared centromere and percentric components [24,70].
The various explanations posited above are admittedly speculative. What these results unambiguously show is that sister chromatids experience a force regime shift whereby the direction and magnitude of the forces (inferred by position vs. time data) co-occur in both direction and magnitude. This information gives new insight into the effective forces experienced by particles in time changing, complex live cellular environment. Researchers can get a more reliable "force map" that contains rich dynamical information beyond that contained in a mean square displacement or correlation-based analysis. Note that this "force map" allows a spatially dependent vector field to change over time (an important feature for detecting dynamical changes in single-molecule measurements). A dominant force direction pointing orthogonal to the spindle axis is highly suggestive of a new dominant force (not induced by the microtubule pulling) being felt by both sister chromatids. Deciding the correct physical explanation of the phenomenon observed requires additional experiments, but the results obtained on the metaphase trajectories are presented in order to show the power of SPT analysis coupled with HDP-SLDS segmentation and downstream likelihood-based time series analysis [23]. Currently we are attempting to use this type of analysis to more systematically and quantitatively probe forces associated with the various physical phenomena described here using experiments which more carefully control the types of forces that can be experienced by the GFP tagged chromatids.
Before proceeding, we should note that the GFP array used to tag the chromatids induces a gradual continuous trend in the effective measurement noise (the finding can readily be observed in S1-S3 Data). A gradual time trend induced in large part by different fluorophores bleaching out has been observed and quantitatively studied in previous works (see Supp. Mat. of Ref. [23]). However, the HDP-SLDS formulation considered assumes a constant measurement noise for the entire duration of the trajectory. S6 -S7 Figs. explore various simulations mimicking features in the experimental data; specifically we analyze data where the effective measurement changes linearly from 15nm to 35nm over the course of 400 measurements (this induces model misspecification since the HDP-SLDS model assumes fixed covariance). The linear change in measurement noise is consistent with the MLE range of values observed in S3 Data. Figs. S6 -S7 indicate that if one uses a prior near the temporal midpoint of the time changing measurement noise that the HDP-SLDS segmentation can achieve accuracy similar to the ideal situation where the HDP priors match the DGP exactly. For whatever system is being studied, it is highly recommended to employ simulations mimicking features implied by the data in order to be more confident that the segmentation output does not contain artifacts induced by features not explicitly accounted for in the underlying model. Techniques accounting for gradual time trends in certain parameters extending the HDP-SLDS model would be an interesting future research topic.
In terms of more general applicability of this type of analysis, systematically segmenting SPT trajectories into distinct dynamical regimes and then using classical statistical physics models is one alternative to using so-called anomalous diffusion models commonly used in SPT analysis [20,26,71]. This is advantageous since fitted overdamped Langevin model output can readily be interpreted in terms of force, velocity, and molecular friction. Furthermore, an overdamped Langevin SDE model can be tested directly against experimental data and the effects of measurement noise can be systematically accounted for by likelihood methods [23,28,37,38]. Anomalous diffusion models are more difficult to physically interpret and rigorous time series analysis (accounting for measurement noise, sampling noise, spatially dependent forces, etc.) is challenging for various technical reasons [23].

Discussion
We have demonstrated the utility of the HDP-SLDS method of Fox et al. [32] in automatically segmenting SPT data into different dynamical regimes. The approach shows great promise in systematically processing a variety of SPT data sets. When applied to experimental data, we demonstrated how new quantitative information can be extracted about transient forces experienced during mitosis in live yeast cells; the method explicitly accounts for the statistical effects of measurement noise on top of "thermal" or "process" noise (both sources are non-negligible in many SPT applications). Simulations results reported in S2 Text were used to exhibit strengths and weaknesses of the HDP-SLDS approach within an SPT context. Work presented in Ref. [47] compared the HDP-SLDS approach to classic (finite state space) HMM models.
In terms of strengths, the HDP-SLDS approach can account for an unknown number of states, position dependent anisotropic active forces, and can account for the statistical effects of measurement noise by exploiting Kalman filtering modeling in conjunction with HDP ideas [41]; these features are not accounted for in any HMM approach currently used in SPT analysis and hence the HDP-SLDS represents a "state-of-the-art" technology within this domain. The approach provides a more systematic means for the "time window size" selection problem discussed in Ref. [23]. The HDP-SLDS approach was also confirmed to be fairly robust to many modeling assumptions and hyperparameters associated with HDP modeling [32]. In situations where a collection of transient dynamical responses are experienced within a single trajectory and the "states" causing the different kinetics are experimentally resolvable, the data processing procedure discussed here provides an attractive alternative to anomalous diffusion modeling [20,23,26,36,71] since the procedure can mitigate artifacts induced by aggregating distinct kinetic states [20,23]. Furthermore, the HDP-SLDS approach can be modified to produce output that can be readily physically interpreted in terms of classic SPT models (see Sec.1.1.1).
A substantial weakness of the HDP-SLDS approach (also discussed in Ref. [46]) is associated with its dependence on reliable priors for the parameters determining governing the Kalman filter parameters (i.e., hyperparameters specifying the "base measure" [46]). The data-driven scheme advocated in Ref. [32] for selecting the prior mean of R exhibits undesirable properties in many SPT applications (results presented in S2 Text). A variant of the MLE-based approach of Ref. [23,37] was discussed (S1 Text) and was demonstrated (S2 Text) to be capable of providing pilot estimates needed to tune the base measure hyperparameters priors in SPT applications. The formal hypothesis testing methods utilized in Ref. [23,37] also allow one to test the validity of a candidate segmentation and assumed dynamical model against experimental data in situations where reference (or "ground-truth") data is unavailable.
In summary, combining recent nonparametric Bayesian modeling ideas, like the HDP-SLDS, with frequentist ideas show great potential in quantitatively analyzing complex SPT data sets. In the future, we aim to provide user-friendly software making the tools discussed herein readily accessible to the diverse set of researchers involved in SPT. The development of new, statistically rigorous data processing algorithms capable of reliably extracting information about molecular motion from live cell data is becoming ever more important as in vivo microscopy techniques continue to improve in spatial and temporal resolution [10,[13][14][15][16]. particle is denoted by thin red line (i.e., the trajectory of r) and the simulated discrete measurements, ψ, are represented by the grey lines with symbols. The parameter determining the linear dynamical system is denoted by θz. At random times, the latent state z changes and hence the underlying parameter specifying the dynamics, θz, also changes. Note that θz models both thermal and measurement noise (parameters defined in Eqn. 1). The latter is crucial for identifying subtle changes in the underlying stochastic process. For example, the dramatic change in thermal fluctuations shown around observations 300-400 would be difficult to detect without jointly modeling thermal and measurement noise. To achieve this joint modeling, one must also account for temporal correlations induced by measurement noise [18,36,37].