We considered how to infer dynamic subpopulations from the short, fragmented trajectories produced by sptPALM in a manner robust to the effects of localization error and defocalization (Figure 1).

A simple and popular approach to this problem is to make a separate estimate for the parameters of each trajectory, then compile a histogram of the results. In the case of Brownian motion, we refer to this method as the ‘MSD histogram’ approach since the MSD is the maximum likelihood estimator for the diffusion coefficient of a Brownian motion with no localization error.

Real estimates of a particle’s position, however, are invariably associated with localization error. In sptPALM, this problem is more significant due to motion blur, which increases the magnitude of the error (Figure 1—figure supplement 1). To incorporate these effects, we refer to the combination of regular Brownian motion with normally distributed, mean-zero localization error as ‘RBME’ (‘Materials and methods’). Each RBME is characterized by two parameters: the diffusion coefficient and the localization error variance. (For brevity, we refer to the latter simply as ‘localization error.’) Importantly, the increments of RBME are only Markovian when the localization error is zero (Martin et al., 2002; Figure 1—figure supplement 1).

Because individual trajectories produced by sptPALM are usually too short to estimate localization error, and because it does not take into account other effects like defocalization, the MSD histogram approach is prone to large systematic biases (Michalet and Berglund, 2012, Hansen et al., 2018). While techniques exist to mitigate some biases of MSD fitting (Kepten et al., 2015), most are difficult to apply at the single trajectory level due to the small number of points per trajectory.

A distinct approach is represented by Bayesian finite state mixture models (Marin et al., 2005, McLachlan et al., 2019; Figure 2A, Figure 2—figure supplement 1A). Such models are comprised of a collection of states labeled $k=1,\mathrm{\dots },K$. Each state is associated with an occupation ${\tau }_k$ (describing the probability to observe trajectories from that state) and a vector of state parameters ${\mathit{\theta }}_k$ (describing the kind of trajectories produced by that state). Importantly, ${\mathit{\theta }}_k$ can also incorporate measurement parameters like the localization error. The probability to observe a particular trajectory $x$ is then ${\displaystyle \sum _{k=1}^K{\tau }_k{p}_X(x|{\mathit{\theta }}_k)}$, where $p_X(x|{\mathit{\theta }}_k)$ is a distribution over trajectories produced by state $k$ and depends on the type of motion being considered. The goal is to infer ${\tau }_k$ and ${\mathit{\theta }}_k$ for each state given some observed set of trajectories $\mathit{F}$. A challenge with such methods is choosing the number of states $K$ as well as the high computational cost when $p_X(x|\mathit{\theta })$ is nonconjugate to the prior over $\mathit{\theta }$.

Figure 2 with 1 supplement see all

Download asset Open asset

Potential solutions can be found in the Bayesian nonparametric class of methods. These approaches begin with a single model comprising a very large or infinite collection of states. A Bayesian inference algorithm is then used to prune away superfluous complexity, leaving a sparse subset of states sufficient to explain the observed trajectories. The foundational example is the DPMM (Ferguson, 1973), which has the distinct advantage of being able to approximate essentially any mixture of states, discrete or continuous (Neal, 1992, Teh, 2010; Figure 2B). Its disadvantage is the high computational cost associated with inference, which becomes especially severe when considering types of motion with multiple parameters (such as RBME) (Neal, 2000, Andrieu et al., 2003).

We considered two responses to this challenge. First, we constructed a DPMM that uses a cheap approximation to RBME by treating the RBME as a Markov process (Matsuda et al., 2018; Figure 3C). This assumption is strictly true only when the localization error is zero and is the same assumption used to estimate diffusion coefficient via the MSD (Michalet and Berglund, 2012). Because localization error is never actually zero, we were curious to see when and how this method breaks down.

Figure 3 with 2 supplements see all

Download asset Open asset

The second approach we explored is a model we refer to as a ‘state array’ (SA). This model is a special case of the finite state mixture, obtained by selecting a large number of states $K$ and fixing the state parameters to the vertices of an ‘array’ that spans some target parameter space (Figure 2C, Figure 2—figure supplement 1). For example, the array for RBME might span a range of biologically plausible diffusion coefficients and localization error variances. An array for an anomalous diffusion model may also incorporate one or more anomaly parameters. The occupation of each ‘state’ in this array is inferred through a variational Bayesian algorithm, driving the occupation of most states to zero to leave a minimal set sufficient to explain the observations (‘Materials and methods’). Importantly, SAs jointly infer a ‘global’ distribution over the state parameters along with ‘individual’ distributions for each trajectory. The nature of the variational inference algorithm means that the ‘global’ distribution is always a weighted mean of these ‘individual’ distributions. We focus our attention on the global distribution in this article, with some consideration of the individual distributions for each trajectory at the end.

Because the parameters for each state in an SA are fixed, the most expensive computations can be cached and reused throughout inference. As a result, SAs can handle more complex models than DPMMs. In this article, we use a 2D SA for RBME spanning a range of diffusion coefficients and localization error variances. After inference, we marginalize out the localization error part to yield 1D functions of the diffusion coefficient (Figure 3B). This procedure naturally incorporates uncertainty about localization error variance, rendering SAs more robust to variations in localization error than DPMMs (Figure 3—figure supplement 1).

DPMMs and SAs work best with thousands to tens of thousands of trajectories. This often requires aggregating trajectories across multiple cells, which can mask cell-to-cell variability. To assess cell-to-cell variability, we also found it useful to have a ‘cheap and dirty’ estimate of state occupation that works with a smaller number (100 s) of trajectories. This is derived from the SA calculation and is simply the sum of the normalized RBME likelihood function across all of the trajectories observed in a cell. We refer to this as the ‘naive occupation estimate.’ Functionally it behaves like a less precise version of the SA method (‘Materials and methods’).

Finally, to account for defocalization we developed a method applicable to the posterior distributions of both DPMMs and SAs (Figure 3—figure supplement 2, ‘Materials and methods’).

As the target for inference, we considered a mixture of RBMEs enclosed in a spherical membrane with a thin focal volume bisecting the sphere, with dimensions similar to a mammalian cell nucleus. Emitters photoactivate and photobleach throughout the sphere and are only observed when their positions coincide with the focal volume. Because no gaps are allowed during tracking, the result is a highly fragmented set of trajectories with mean length 3–5 frames. We chose simulation settings to approximate real sptPALM experiments, with bleaching rates ≥10 Hz, diffusion coefficients in the range 0–100 μm² s^-1, and localization error variances between 0² and 0.06²μm².

We compared the ability of DPMMs, SAs, and MSD histograms to recover the underlying distribution of diffusion coefficients from this data. We divided these simulations into four classes with increasing difficulty. In class 1, localization error for all states was provided as a known constant to the algorithms (Figure 4A, Figure 4—figure supplement 1A). In class 2, localization error was held constant for all states but was unknown to the algorithms (Figure 4B, Figure 4—figure supplement 1B). In class 3, localization error was allowed to vary between diffusive states and was also unknown to the algorithms (Figure 4C, Figure 4—figure supplement 1C). Finally, for class 4 we simulated full sptPALM-like movies that incorporate heterogeneous localization error, motion blur, camera noise, tracking errors, and defocus (Figure 4—figure supplement 5, Figure 4—figure supplement 6). In these simulations, the localization error is unique for each emitter and depends on the emitter’s axial position, the stochastic number of photons it emits during each integration, and its pattern of motion blur (Video 4, Video 5).

Figure 4 with 12 supplements see all

Download asset Open asset

Video 4

Download asset

Video 5

Download asset

DPMMs and SAs both recovered the dynamic profile for simulations in class 1 with a resolution that exceeded the MSD histogram approach. With large samples of trajectories, DPMMs and SAs inferred even nondiscrete distributions of states (Figure 4A, Figure 4—figure supplement 1A).

When knowledge of the localization error was removed (classes 2 and 3), the SA approach outperformed both the MSD and DPMM approaches. The DPMM’s performance was especially poor when the contributions of diffusion and error to jump variance were similar ($D\mathrm{\Delta }t\approx {\sigma }_{\text{loc}}^2$), likely due to its simplistic treatment of localization error. Meanwhile, the dynamic profile estimated by SAs was unperturbed by variations in the localization error (Figure 4B and C, Figure 4—figure supplement 1B and C). Comparing the results from simulations in class 3 numerically, we found that the root mean squared deviation of the estimated CDF from the true CDF was ≤ 5% for SAs, while it was 5–20% for both the MSD histogram and DPMM approaches (Figure 4—figure supplement 2).

The dynamic profiles produced by the MSD, DPMM, and SA approaches can be integrated to yield occupation estimates over particular diffusion coefficient ranges. We compared the accuracy and precision of these estimates with discrete two-, three-, or four-state models (Figure 4D, Figure 4—figure supplement 3, Figure 4—figure supplement 4). As the number of trajectories increased, occupations estimated by DPMMs and SAs converted to within 3% of the true values. In contrast, the MSD approach was associated with large systematic errors, an effect previously reported (Mazza et al., 2012, Hansen et al., 2018).

On full optical and dynamic simulations in class 4, SAs also outperformed the DPMM approach (Figure 4—figure supplement 5, Figure 4—figure supplement 6). Again, the difference was particularly pronounced for small diffusion coefficients, for which the DPMM state occupation estimates were severely inaccurate. Both methods had difficulty recovering the fastest diffusion coefficient tested (Figure 4—figure supplement 5B), possibly due to the restrictive conditions on the maximum jump distance used during tracking.

A central limitation of DPMMs and SAs is that they do not account for transitions between diffusive states. To determine the effect of state transitions on the output of these algorithms, we simulated mixtures of two diffusive states with increasing transition rates (Figure 4G, Figure 4—figure supplement 7). While slow transition rates had a negligible effect on the estimated state profile, transition rates approaching the frame interval appeared as single state with intermediate diffusion coefficient (Figure 4—figure supplement 7C), consistent with a result from reaction-diffusion systems (Crank, 1975). The shift from the two-state to single-state regime occurred in a narrow window of mean state dwell times between 0.05 and 0.5 frame intervals.

In this article, we restricted DPMM/SA inference to a range of diffusion coefficients from 10^-2 to 10² μm² s^-1. We also explored what happens when the true diffusion coefficient lies outside this range. DPMMs and SAs still recovered the correct state occupations by using the closest diffusion coefficient in their respective supports (Figure 4—figure supplement 8).

In the presence of multiple diffusing states with similar diffusion coefficients, both DPMMs and SAs tended to identify a single population with occupation equal to the sum of the occupations for each true state (Figure 4F, Figure 4—figure supplement 9).

We compared the performance of SAs and vbSPT (Persson et al., 2013) using simulated SPT movies with different dynamic models (Figure 4—figure supplement 10). Both methods had comparable accuracy on simple two-state models (Figure 4—figure supplement 10B). On more complex models (Figure 4—figure supplement 10C), both methods encountered distinct difficulties, with vbSPT tending to overestimate and SAs tending to underestimate the number of states. For clusters of states with similar parameters (Figure 4—figure supplement 10C, bottom), SAs tend to produce a ‘smear’ of state occupations over a range of diffusion coefficients, while vbSPT tended to produce a different cluster of states in the same region of parameter space. vbSPT was noticeably less accurate at recovering slow-moving states with small diffusion coefficients (<0.1 μm² s^-1). We concluded that both approaches are useful and may provide complementary information.

While our investigation focused primarily on Brownian motion, SAs can be applied to any motion model parameterized by a likelihood function. To explore applications of SAs outside of Brownian motion, we applied it to fractional Brownian motion (FBM), a generalization of Brownian motion capable of producing anomalous diffusion (Mandelbrot and Van Ness, 1968). Whereas Brownian motion’s sole parameter is the diffusion coefficient, FBM parameterizes both the magnitude (via a scaling coefficient) and the temporal correlations (via the Hurst parameter) of a particle’s increments. As with Brownian motion, we simulated sptPALM movies with fraction Brownian particles with variable diffusion coefficient and Hurst parameter (Video 6). To construct a state array for FBM, we used a 3D array over scaling coefficient, Hurst parameter, and localization error variance (Figure 4—figure supplement 11C). As with the RBME array, we marginalized out localization error after inference. While the SA accurately recovered the diffusion coefficient and Hurst parameter for multistate FBM models (Figure 4—figure supplement 11D), we noted a systematic error in the estimation of low (subdiffusive) Hurst parameters due to motion blur (Figure 4—figure supplement 12).

Video 6

Download asset

After observing that SAs outperformed DPMMs on simulations, we proceeded to evaluate SAs on real data. We acquired an sptPALM dataset in U2OS osteosarcoma nuclei with endogenously tagged retinoic acid receptor-α-HaloTag (RARA-HT) (Pontén and Saksela, 1967, Los et al., 2008; (Figure 5—figure supplement 1). RARA-HT is a type II nuclear receptor that heterodimerizes via its ligand-binding domain (LBD) with the retinoid X receptor (RXR) to form a complex competent to bind chromatin and regulate target genes Giguere et al., 1987, Petkovich et al., 1987, Brand et al., 1988, Yu et al., 1991, Bugge et al., 1992, Marks et al., 1992, Leid et al., 1992; reviewed in Evans and Mangelsdorf, 2014). In addition, association of coregulator complexes with the RAR/RXR heterodimer has been shown to influence the dimer’s dynamics in FCS studies (Brazda et al., 2011, Brazda et al., 2014). As such, RARA-HT is expected to inhabit a variety of dynamic states in sptPALM.

For comparison, we also performed identical sptPALM experiments with histone H2B-HaloTag (H2B-HT), a protein with a high-occupation immobile state (Hansen et al., 2017, McSwiggen et al., 2019), as well as HaloTag and HaloTag-NLS (HT and HT-NLS), which are fast-diffusing proteins with low immobile fractions.

The four proteins presented distinct dynamic profiles (Figure 5A). For both HT and HT-NLS, the SA identified a single highly mobile state. In agreement with previous reports (Xiang et al., 2020), we observed that addition of the NLS reduces HaloTag’s diffusion coefficient by two- to threefold. In contrast, both RARA-HT and H2B-HT had substantial immobile fractions, accounting for roughly 40 and 70% of their total populations, respectively (Figure 5C). SAs identified stark differences in the mobile subpopulations for RARA-HT and H2B-HT. Whereas H2B-HT presented a fast population at 8–10 μm² s^-1, RARA-HT inhabited a broad spectrum of diffusing states ranging from 0.3 to 10.0 μm² s^-1. Biological replicates gave similar results (Figure 5—figure supplement 2A).

Figure 5 with 2 supplements see all

Download asset Open asset

Figure 5—source data 1

Raw and labeled RARA-HaloTag Western blots used in Figure 5.

https://cdn.elifesciences.org/articles/70169/elife-70169-fig5-data1-v2.zip

Download elife-70169-fig5-data1-v2.zip

To determine the origins of the dynamic states observed for RARA-HT, we performed domain deletions (Figure 5B). Removal of either the DNA-binding domain (DBD) or LBD resulted in loss of the immobile population. Because both the DBD and LBD are required for chromatin binding by the RAR/RXR heterodimer, this suggests that the immobile fraction represents chromatin-bound molecules. To confirm this, we introduced a point mutation (C88G) in the zinc fingers for the RARA-HT DBD that abolishes DNA-binding in vitro (Zhu et al., 1999). This led to loss of the immobile fraction (Figure 5B). Deletion of the unstructured N-terminal domain (NTD) or C-terminal domain (CTD) had a milder effect, suggesting that these domains are not the primary determinants of the dynamic behavior of RARA-HT.

To understand the origins of heterogeneity in the diffusive profile, we performed three variants of bootstrap aggregation (Figure 5—figure supplement 2B). The primary origins of variability for both DPMMs and SAs were cell-to-cell rather than clone-to-clone variability or intrinsic variability due to finite sample sizes.

In the process of inferring the global distribution over state parameters for an sptPALM dataset, SAs jointly infer individual distributions for each trajectory. Up to this point, we have analyzed the global distribution. However, it is also possible to aggregate the individual distribution for each trajectory as a function of space or time, yielding, for instance, separate dynamic profiles for every spatial location in an experiment. This approach offers a potential route to understand spatiotemporal variation in the dynamics of a protein target.

We explored this aspect of SAs with a U2OS nucleophosmin-HaloTag (NPM1-HT) sptPALM dataset. NPM1-HT exhibits partial nucleolar localization (Figure 6—figure supplement 1B) and distinct dynamic behavior inside and outside nucleoli (Mitrea et al., 2018). The SA identified a broad range of diffusion coefficients for NPM1-HT, with three modes including an effectively immobile population (Figure 6A). Selecting four ranges of diffusion coefficients for analysis (Figure 6A), we visualized the posterior distribution as a function of space, calculating local fractional occupations for each range (Figure 6B, Figure 6—figure supplement 1C). This analysis revealed that some populations (including a slow-moving mobile population at 0.23 μm² s^-1) are enriched in nucleoli, while others (for instance, a fast-moving population at 4 μm² s^-1) are depleted and still others show no preference (Figure 6C). Notably, these preferences are apparent even in the naive occupations for trajectories in each compartment (Figure 6—figure supplement 1D).

Figure 6 with 1 supplement see all

Download asset Open asset

The NPM1-HT tracking experiments were performed with an acquisition sequence comprising several phases with distinct levels of photoactivation. As a result, the localization density varied temporally in each movie. To understand the effect of localization density on the diffusion coefficient likelihoods, we aggregated the naive state occupations over 100-frame temporal blocks (Figure 6D). These experiments demonstrated that high localization densities led to a deflation in the occupation of slower-moving states, probably due to tracking errors. As a result, only phases with low localization density were used for posterior estimation. This demonstrates how the temporal perspective on the posterior may be useful as a guide for subsequent analysis, including quality control in SPT experiments.

Unless otherwise noted, all PCRs were performed with New England Biosciences Phusion High-Fidelity DNA polymerase (M0530S), and Gibson assemblies (Gibson et al., 2009) were performed with New England Biosciences Gibson Assembly Master Mix (E2611S) following the manufacturer’s instructions. Cloning and expression of plasmids was performed in E. coli DH5α using the Inoue protocol (Im et al., 2011). Plasmids used for nucleofections were purified with Zymo midiprep kit (Zymo D4200) and concentrations were quantified by absorption at 260 nm. Cloning primers were synthesized by Integrated DNA Technologies as 25 nmol DNA oligos with standard desalting, and sequences were verified by Sanger sequencing at the UC Berkeley DNA Sequencing Facility. A complete list of the primers used in this article is provided in Supplementary file 1, and a complete list of the plasmids used in this article is provided in Supplementary file 2.

We produced the vector PB PGKp-PuroR L30p MCS-GDGAGLIN-HaloTag-3xFLAG by amplifying the human L30 promoter with prAH675 and prAH676 and assembling into AsiSI- (NEB R0630) and XbaI- (NEB R0145) digested PB PGKp-PuroR EF1a MCS-GDGAGLIN-HaloTag-3xFLAG. For the expression plasmid PB PGKp-PuroR EF1a 3x-FLAG-HaloTag-GDGAGLIN, we cloned three tandem copies of the SV40 nuclear localization sequence into XbaI- and BamHI-HF (NEB R3136)-digested PB PGKp-PuroR EF1a 3xFLAG-HaloTag-MCS using Gibson assembly.

For constructs expressing RARA-HaloTag domain deletions and point mutations, we first cloned the RARA coding sequence out of U2OS cDNA by extracting RNA from cycling U2OS cells with a QIAGEN RNeasy kit (QIAGEN 74104), preparing cDNA with the iScript Reverse Transcription Supermix (Bio-Rad 1708840), amplifying the CDS with prAH495 and prAH496, then assembling into an XbaI- and NotI-HF- (NEB R3189) digested PB PGKp-PuroR EF1a MCS-GDGAGLIN-HaloTag-3xFLAG using Gibson assembly. Next, to produce the mutants, we amplified parts of the RARA coding sequence in PCR fragments while introducing point mutations or domain deletions at the intersections of the fragments. PCR fragments were assembled into XbaI- and BamHI-HF-digested PB PGKp-PuroR L30p-MCS-GDGAGLIN-HaloTag-3xFLAG using Gibson assembly. The primers used for each construct were as follows: for PB PGKp-PuroR EF1a RARA[ΔNTD]-HaloTag-GDGAGLIN-3xFLAG, PCR fragment 1 was produced with prAH1111 and prAH1112; for PB PGKp-PuroR EF1a RARA[ΔCTD]-HaloTag-GDGAGLIN-3xFLAG, PCR fragment 1 was produced with prAH1113 and prAH1114; for PB PGKp-PuroR EF1a RARA[ΔNTD,ΔCTD]-HaloTag-GDGAGLIN-3xFLAG, PCR fragment 1 was produced with prAH1111 and prAH1114; for PB PGKp-PuroR EF1a RARA[C88G]-HaloTag-GDGAGLIN-3xFLAG, PCR fragment 1 was produced with prAH1113 and prAH1069 and PCR fragment 2 was produced with prAH1112 and prAH1070; for PB PGKp-PuroR EF1a RARA[ΔDBD]-HaloTag-GDGAGLIN-3xFLAG, PCR fragment 1 was produced with prAH596 and prAH704 and PCR fragment 2 was produced with prAH597 and prAH705; for PB PGKp-PuroR EF1a RARA[ΔLBD]-HaloTag-GDGAGLIN-3xFLAG, PCR fragment 1 was produced with prAH596 and prAH706 and PCR fragment 2 was produced with prAH597 and prAH707.

To generate the plasmid-based homology repair donor for gene editing at the human RARA exon 9 locus, we assembled the following fragments by Gibson assembly. For fragment 1, we digested the pUC57 vector with EcoRI and HindIII. For fragment 2, we amplified the left homology arm out of U2OS genomic DNA with prAH599 and prAH600. For fragment 3, we amplified the GDGAGLIN-HaloTag-3xFLAG insert out of the plasmid PB PGKp-PuroR L30p MCS-GDGAGLIN-HaloTag-3xFLAG with prAH601 and prAH602. For fragment 4, we amplified the right homology arm out of U2OS genomic DNA with prAH603 and prAH604.

To generate guide RNA/Cas9 expression plasmids for gene editing at the human RARA exon 9 locus, we cloned the two guide RNA sequences under a U6 promoter in a vector that coexpresses the sgRNA, mVenus, and S. pyogenes Cas9, which has been previously described (Hansen et al., 2017).

In luciferase assays, we used the retinoic acid-responsive firefly luciferase expression vector pGL3-RARE-luciferase (Addgene plasmid #13458; http://n2t.net/addgene:13458; RRID:Addgene_13458), a gift from T. Michael Underhill (Hoffman et al., 2006). Renilla luciferase was expressed from pRL CMV Renilla (Promega E2261).

Human U2OS cells (female, 15 years old, osteosarcoma) obtained from the UC Berkeley Cell Culture Facility were cultured under 5% CO₂ at 37°C in DMEM containing 4.5 g/L glucose supplemented with 10% fetal bovine serum and 10 U/mL penicillin-streptomycin. Cells were subpassaged at a ratio of 1:6 every 3–4 days. The stable cell line expressing H2B-HaloTag-SNAPf was described previously (Hansen et al., 2017, McSwiggen et al., 2019). We induced exogenous expression of HaloTag, HaloTag-NLS, and point mutants and domain deletions of RARA-HaloTag by nucleofection of PiggyBac vectors containing the proteins under EF1a promoters. Expression of wildtype RARA-HaloTag and NPM1-HaloTag was induced by endogenous gene editing, as described in the ‘CRISPR/Cas9-mediated gene editing’ section. The U2OS cell line used here was validated by whole-genome sequencing as described in Hansen et al., 2017, and mycoplasma testing was performed by DAPI staining.

For sptPALM experiments, cells were grown on 25 mm circular No. 1.5 H coverglasses (Marienfeld, Germany, High-Precision 0117650) that were first sonicated in ethanol for 10 min, plasma-cleaned, then stored in isopropanol until use. U2OS cells were grown directly on the coverglasses in regular culture medium. The medium was changed after dye labeling and immediately before imaging into phenol red-free medium.

For all imaging experiments involving exogenous expression, we used the Lonza Amaxa II Nucleofector System with Cell Line Nucleofector Kit V reagent (Lonza VCA-1003). Briefly, U2OS cells were grown in 10 cm plates (Thermo Fisher) for 2 days prior to nucleofection, trypsinized, spun down at 1200 rpm for 5 min, combined with vector and Kit V reagent according to the manufacturer’s instructions, and nucleofected with program X-001 on an Lonza Amaxa II Nucleofector. After nucleofection, cells were immediately resuspended in regular culture medium at 37°C and plated onto coverslips. In all imaging experiments involving nucleofection, imaging was performed within 24 hr of plating.

Endogenous tagging of RARA in U2OS cells was performed with a protocol roughly following Hansen et al., 2017 with some modifications. A complete list of the plasmids used in gene editing is provided in Supplementary file 2, and a list of the guide sequences is provided in Supplementary file 3.

For U2OS cells, we nucleofected cells with plasmid expressing 3xFLAG-SV40NLS-pSpCas9 from a CBh promoter (Ran et al., 2013), mVenus from a PGK promoter, and guide RNA from a U6 promoter (pU6_sgRNA_CBh_Cas9_PGK_Venus_anti-RARA-C_terminus_1 and pU6_sgRNA_CBh_Cas9_PGK_Venus_anti-RARA-C_terminus_2), along with a second plasmid encoding the homology repair donor (pUC57_homRep_RARA-HaloTag). The homology repair donor was built in a pUC57 backbone modified to contain HaloTag-3xFLAG with ~500 base pairs of homologous genomic sequence on either side. Synonymous mutations were introduced at the cut site to prevent retargeting by Cas9. Each of the two guide RNA plasmids were nucleofected into separate populations of cells to be pooled for subsequent analysis. Then, 24 hr after the initial nucleofection, we screened for mVenus-expressing cells using FACS and pooled these mVenus-positive cells in 10 cm plates. Then, 5 days after plating, we labeled cells with HTL-TMR (Promega G8251) and screened for TMR-positive, mVenus-negative cells. Cells were diluted to single clones and plated in 96-well plates for a 2–3-week outgrowth step, during which the medium was replaced every 3 days. The 96-well plates were then screened for wells containing single colonies of U2OS cells, which were split by manual passage into two replicate wells in separate 96-well plates. One of these replicates was used to subpassage, while the other was used to harvest genomic DNA for PCR and sequencing-based screening for the correct homology repair product. In PCR-based genotyping, we used three primer sets: (A) primers external to the homology repair arms, expected to amplify both the wildtype allele and the edited allele (‘PCR1’), (B) a primer internal to HaloTag and another external to it on the 5′ side, expected to amplify only the edited allele (‘PCR2’), and (C) a primer internal to HaloTag and another external to it on the 3′ side, expected to amplify only the edited allele (‘PCR3’). The primer sets for each target were the following. For RARA-GDGAGLIN-HaloTag-3xFLAG, we used prAH586 and prAH761 for PCR1, prAH761 and prAH762 for PCR2, and prAH763 and prAH764 for PCR3. For NPM1-GDGAGLIN-HaloTag-3xFLAG, we used prAH1092 and prAH1093 for PCR1, prAH1093 and prAH377 for PCR2, and prAH1092 and prAH373 for PCR3. U2OS cDNA from selected clones was isolated with DirectPCR Lysis Reagent (Viagen 101T), treated with 0.5 mg/mL proteinase K for 15 min, incubated at 95°C for 1 hr, then subjected to PCRs 1–3 using Phusion polymerase in the presence of 5% DMSO. Amplicons from candidate clones were gel-purified (QIAGEN 28704) and Sanger sequenced; only clones with the correct target sequence were kept for continued screening. A subset of these clones were chosen for characterization by Western blot, imaging, and luciferase assays.

For NPM1-GDGAGLIN-HaloTag-3xFLAG knock-in cell lines, we used a different strategy relying on nucleofected Streptococcus pyogenes Cas9 sgRNPs and linear dsDNA homology repair donors. The target insert (GDGAGLIN-HaloTag-3xFLAG from the vector PB PGKp-PuroR L30p MCS-GDGAGLIN-HaloTag-3xFLAG) was first amplified with ultramers encoding 120 bp homology arms (prAH867 and prAH868; IDT) using KAPA2G Robust HotStart polymerase (Kapa Biosystems KR0379) for 12 cycles. A small volume of this reaction was then used to seed a PCR reaction using primers prAH869 and prAH870 in Q5 High-Fidelity 2X Master Mix (QIAGEN M0492). Products were purified by RNAClean XP magnetic beads (Beckman-Coulter A63987) and further cleaned by ethanol precipitation, followed by resuspension in a small volume of RNase-free water. For guides, we performed a three-primer PCR using prAH2000 and prAH2001 along with a unique oligo encoding the spacer (either prAH979 or prAH980) to produce a linear dsDNA product encoding the sgRNA preceded by a T7 promoter. We then used T7 RNA polymerase (NEB E2040S) to transcribe sgRNA from this template and purified the sgRNA with RNAClean XP magnetic beads according to the manufacturer’s instructions. To assemble the sgRNP, we incubated 80 pmol sgRNA with 40 pmol purified SpyCas9-NLS (UC Berkeley Macrolab) for 15 min at 37°C in 20 mM HEPES pH 7.5, 150 mM KCl, 10 mM MgCl₂, and 5% glycerol. sgRNPs were subsequent kept on ice and combined with donor immediately before nucleofection. For each nucleofection, we used 40 pmol sgRNP and 5 pmol dsDNA donor template suspended in <10 μL with Lonza Amaxa Nucleofector II protocol X-001 in Lonza Kit V reagent. Roughly 1 million cells were used for nucleofection. Sorting for labeled cells, subcloning, and genotyping proceeded as previously described for RARA-GDGAGLIN-HaloTag-3xFLAG.

Antibodies were as follows (the ratio indicates the dilution factors used for Western blot): human TBP, Abcam Ab51841, 1:2500 (mouse); FLAG, Sigma-Aldrich F3165, 1:2000 (mouse).

For Western blots, cells were scraped from plates in ice-cold PBS then pelleted. Pellets were resuspended in lysis buffer (0.15 M NaCl, 1% NP-40, 50 mM Tris–HCl [pH 8.0], and a cocktail of protease inhibitors [Sigma-Aldrich 11697498001 dissolved in PBS with supplemented PMSF, aprotinin, and benzamidine]), agitated for 30 min at 4°C, then centrifuged for 20 min at 12,000 rpm, 4°C. The supernatant was then mixed with 2× Laemmli (to final 1×), boiled for 5 min, then run on 12.5% SDS-PAGE. After transfer to nitrocellulose, the membrane was blocked with 10% condensed milk in TBST (500 mM NaCl, 10 mM Tris–HCl [pH 7.4], 0.1% Tween-20) for 1 hr at room temperature. Antibodies were suspended in 5% condensed milk in TBST at the dilutions indicated above and incubated rocking at 4°C overnight. After primary hybridization, the membrane was washed three times for 10 min with TBST at room temperature, hybridized with an anti-mouse HRP secondary antibody in 5% condensed milk in TBST for 60 min at room temperature, washed three more times with TBST for 10 min, then visualized with Western Lightning Plus-ECL reagent (PerkinElmer NEL103001) according to the manufacturer’s instructions and imaged on a Bio-Rad ChemiDoc imaging system. Different exposure times were used for each antibody. The raw Western blots images for RARA-HaloTag and NPM1-HaloTag are provided as Figure 5—source data 1 and Figure 6—figure supplement 1—source data 1 , respectively.

All luciferase assays used pGL3-RARE-luciferase, a reporter containing firefly luciferase driven by an SV40 promoter with three retinoic acid response elements (RAREs). pGL3-RARE-luciferase was a gift from T. Michael Underhill (Addgene plasmid 13458; http://n2t.net/addgene:13458; RRID:Addgene_13458; Hoffman et al., 2006). Luciferase assays were performed on cells cultivated in 6-well plates. Cells were transfected with 100 ng pGL3-RARE-luciferase and 10 ng pRL Renilla (Promega E2261) using Mirus TransIT-2020 Transfection Reagent (Mirus MIR 5404). Transfection was performed 1 day before assaying luciferase expression with the Dual-Luciferase Reporter Assay System (Promega E1910) according to http://n2t.net/addgene:13458 the manufacturer’s instructions. Readout was performed on a GloMax luminometer (Promega).

For sptPALM experiments, cells were labeled with one of two methods, depending on the dye. For non-photoactivatable fluorescent dyes including TMR-HTL (tetramethylrhodamine-HaloTag ligand; Promega G8251), we stained cells with 100 nM dye in regular culture medium for 10 min, then performed three 10 min incubations in dye-free culture medium separated by PBS washes. All PBS and culture medium was incubated at 37°C between medium changes and washes.

For experiments with photoactivatable dyes, which have lower cell permeability and slower wash-in/wash-out kinetics, we labeled cells with 100 nM dye in regular culture medium for 10–20 min, followed by four 30 min incubations in dye-free culture medium at 37°C. Between each incubation, we washed twice with PBS at 37°C. After the final incubation, cells were changed into phenol red-free medium for imaging.

sptPALM experiments were performed with a custom-built Nikon TI microscope equipped with a ×100/NA 1.49 oil-immersion TIRF objective (Nikon apochromat CFI Apo TIRF 100X Oil), an EMCCD camera (Andor iXon Ultra 897), a perfect focus system to account for axial drift, an incubation chamber maintaining a humidified 37°C atmosphere with 5% CO₂, and a laser launch with 405 nm (140 mW, OBIS, Coherent), 488 nm, 561 nm, and 633 nm (all 1 W, Genesis Coherent) laser lines. Laser intensities were controlled by an acousto-optic Tunable Filter (AA Opto-Electronic, AOTFnC-VIS-TN) and triggered with the camera TTL exposure output signal. Lasers were directed to the microscope by an optical fiber, reflected using a multi-band dichroic (405 nm/488 nm/561 nm/633 nm quad-band, Semrock) and focused in the back focal plane of the objective. The angle of incident laser was adjusted for highly inclined laminated optical sheet (HiLo) conditions (Tokunaga et al., 2008). Emission light was filtered using single band-pass filters (Semrock 593/40 nm for PAJFX549 and Semrock 676/37 nm for PAJF646). Hardware was controlled with the Nikon NIS-Elements software.

For stroboscopic illumination, the excitation laser (561 nm or 633 nm) was pulsed for 1–2 ms (most commonly 1 ms) at maximum (1 W) power at the beginning of the frame interval, while the photoactivation laser (405 nm) was pulsed during the ~447 μs camera transition time, so that the background contribution from the photoactivation laser is not integrated. For all sptPALM, we used an EMCCD vertical shift speed of 0.9 μs and conversion gain setting 2. On our setup, the pixel size after magnification is 160 nm and the photon-to-grayscale gain is 109. A total of 15,000–30,000 frames with this sequence were collected per nucleus, during which the 405 nm intensity was manually tuned to maintain low density of fluorescent particles per frame.

To produce trajectories from raw sptPALM movies, we used a custom sptPALM tracking pipeline publicly available on GitHub (https://github.com/alecheckert/quot, swh:1:rev:1adf7a0574c62f38140f1dec2d14555bfc03b9a7, Heckert, 2022b). All localization and tracking for this article was performed with the following settings:

• Detection: Generalized log likelihood ratio test (Sergé et al., 2008) with a 2D Gaussian kernel with $\sigma =190$ nm (detection method ”llr” with k = 1.2, a 15 pixel window size [w = 15], and a log ratio threshold of 16.0 [t = 16.0]).

• Subpixel localization: Levenberg–Marquardt fitting of a 2D integrated Gaussian point spread function model (localization method “ls_int_gaussian”) with fixed $\sigma =190$ nm, window size 9 pixels, maximum 20 iterations per PSF, with a damping term of 0.3 for parameter updates. The 2D integrated Gaussian PSF model is described in Smith et al., 2010 and the Levenberg–Marquardt routine in Laurence and Chromy, 2010. We used the radial symmetry method (Parthasarathy, 2012) to make the initial guess used to start the Levenberg–Marquardt algorithm.

• Tracking: We used the tracking algorithm “conservative” in quot with a 1.2 µm search radius. This simple algorithm searches for particle–particle reconnections that are ‘unambiguous’ in the sense that no other reconnections are possible within the specified search radius. These reconnections are then used to synthesize trajectories, while ‘ambiguous’ connections are discarded.

After localization and tracking, all trajectories in the first 1000 frames of each movie were discarded. Localization density tends to be high in these frames, so they can contribute tracking errors that compromise accuracy. The mean localization density for most movies in the remaining set of frames was less than one emitter per frame.

For experiments involving HaloTag or HaloTag-NLS, which have high mobility, we used a broader search radius at 2.5 μm. All other settings were kept the same.

All trajectories from real sptPALM experiments used in this article are publicly accessible as a Dryad dataset (https://doi.org/10.6078/D13H6N).

Experiments using spinning disk confocal imaging were performed at the UC Berkeley High-Throughput Screening Facility on a Perkin Elmer Opera Phenix equipped with a controller for 37°C and 5% CO₂, using a built-in 40x water immersion objective.

All simulations in this article belong to one of two categories:

• Trajectory simulations: Individual trajectories are simulated in a HiLo-like geometry, then localization error and defocalization are injected as detailed below. The output of each simulation is a set of trajectories. This type of simulation does not incorporate tracking errors. These include the simulations in Figure 3, Figure 4, Figure 4—figure supplement 1, Figure 4—figure supplement 2, Figure 4—figure supplement 3, Figure 4—figure supplement 4, Figure 4—figure supplement 7, Figure 4—figure supplement 8, and Figure 4—figure supplement 9.

• Optical–dynamical simulations: Starting from a dynamical model and an approximation to sptPALM imaging system, we simulate full SPT movies. The output is a stack of images similar to that acquired on a real sptPALM system. Analysis follows the same steps as for real SPT data: we recover trajectories using a localization and tracking algorithm, which are subjected to the relevant downstream analyses. These include the simulations in Figure 1—figure supplement 1, Figure 4—figure supplement 5, Figure 4—figure supplement 6, Figure 4—figure supplement 10, Figure 4—figure supplement 11, and Figure 4—figure supplement 12.

Both types of simulation are important. Trajectory simulations allow us to separate the accuracy of the tracking algorithm from the accuracy of the SA/DPMM algorithm in a tightly controlled setting, while optical–dynamical simulations are ‘end-to-end’ tests that also incorporate realistic features such as motion blur, camera noise, and tracking errors.

This section describes the SA and DPMM used in this article. We begin with a classic Bayesian finite state mixture model, then introduce modifications that lead to SAs and DPMMs. The finite state mixture has been reviewed in detail elsewhere (Marin et al., 2005, McLachlan et al., 2019), so here we keep details to a minimum.

A finite state mixture is a collection of ‘states’ $k=1,\mathrm{\dots },K$, each of which is associated with an occupation ${\tau }_k$ and a vector of state parameters ${\mathit{\theta }}_k$. (Where convenient, we let $\mathbf{\theta }=({\mathit{\theta }}_1,\mathrm{\dots },{\mathit{\theta }}_K)$ be the collection of parameters for all states.) Each state generates trajectories $X$ according to some distribution $p_X(x|{\mathit{\theta }}_k)$. The overall generative process for each trajectory is:

1. Randomly select a state $k$ with probability ${\tau }_k$.

2. Randomly generate a trajectory $X$ from that state according to $p_X(x|{\mathit{\theta }}_k)$.

The probability to generate a particular trajectory $X$ is then

To represent the origin state for each trajectory, we use a 1-of-$K$ encoding ${\mathit{Z}}_i\in {\{0,1\}}^K$ so that $Z_{ik}=1$ if trajectory $i$ originates from state $k$ and $Z_{ik}=0$ otherwise. For a dataset with $N$ trajectories, we let $\mathit{Z}\in {\{0,1\}}^{N\times K}$ be the matrix such that the ith row is ${\mathit{Z}}_i$.

Finally, we specify priors over $\mathit{\tau }$ and ${\mathit{\theta }}_k$. The full Bayesian finite state mixture can then be written as

where $H$ is the prior over the parameters ${\displaystyle {\mathit{\theta }}_k}$, usually chosen to be conjugate to $p_X(x|{\mathit{\theta }}_k)$.

This corresponds to the first graphical model in Figure 2—figure supplement 1. The objective is to infer the posterior distribution $p(\mathit{Z},\mathit{\tau },\mathbf{\theta }|\mathit{F})$, where $\mathit{F}$ represents some observed set of trajectories.

Three common challenges with the finite state mixture 4 are:

1. Choosing $K$, the number of states. Because $K$ is a hyperparameter in 4, some kind of meta-algorithm is required to infer it, and this process can be fraught (Marin et al., 2005).

2. Choosing $H$, the prior over ${\mathit{\theta }}_k$. Ideally the prior is chosen to be conjugate to $p_X(x|{\mathit{\theta }}_k)$, but this is only possible for the simplest forms of $p_X(x|{\mathit{\theta }}_k)$.

3. Computing $p_X(x|{\mathit{\theta }}_k)$ is often expensive, especially if it needs to be evaluated repeatedly during inference.

SAs are a special case of finite mixture models designed in response to these issues. Rather than equating $K$ with the true number of states, SAs instead choose a large, fixed value of $K$ and constant values for each ${\mathit{\theta }}_k$. A Bayesian routine is then used to drive the occupation of most states to zero, leaving minimal models sufficient to explain the observations. (The ability of Bayesian inference to identify sparse explanatory models in the presence of more complex alternatives is the same property that drives automatic relevance determination [ARD] in machine learning with Bayesian models.)

Because the state parameters are constant, the only parameters left to infer are $\mathit{Z}$ and $\mathit{\tau }$. Together, this simplified model is

This corresponds to the third graphical model shown in Figure 2A. Notice that since each $X_i$ and ${\mathit{\theta }}_k$ are constant, $p_X(X_i|{\mathit{\theta }}_k)$ is also constant and only needs to be evaluated once during inference.

To infer the posterior distribution $p(\mathit{Z},\mathit{\tau }|\mathit{F})$, we take a variational approach, constructing an approximation $q(\mathit{Z},\mathit{\tau })\approx p(\mathit{Z},\mathit{\tau }|\mathit{F})$ such that

where $L[q]$ is the variational lower bound:

Motivation for the variational lower bound is discussed in detail elsewhere (Bishop, 2006). Here, we only remark that maximization of ${\displaystyle L[q]}$ minimizes the Kullback–Leibler divergence between the approximation and the true posterior. The factorability criterion in 6 enables an expectation-maximization routine (Dempster et al., 1977) by iteratively evaluating.

The constants are chosen so that each factor, $q(\mathit{Z})$ or $q(\mathit{\tau })$, is normalized. Combining Equations 8 for model 5 yields the solution

where $L_i$ is the number of jumps in trajectory $i$ and $\psi (n)$ is the digamma function. For brevity here, the derivation of Equation 9 is placed in its own section below.

$q(\mathit{Z},\mathit{\tau })$ is parameterized by $\mathit{n}$ and $\mathit{r}$. These can be inferred with a simple EM algorithm:

1. Initialize a matrix $\mathit{A}\in R^{N\times K}$ by setting $A_{ik}=p(X_i|{\mathit{\theta }}_k)$.

2. Initialize a matrix ${\mathit{r}}^{(0)}\in R^{N\times K}$ such that $r_{ik}^{(0)}=\frac{A_{ik}}{\sum _{j=1}^K{A}_{ij}}$.

3. For each iteration $t=1,2,\mathrm{\dots }$:

   1. For each state $k$, evaluate $n_k^{(t)}=\frac{\alpha }{K}+\sum _{i=1}^N{L}_i{r}_{ik}^{(t-1)}$.

   2. Evaluate the matrix ${\mathit{r}}^{(t)}$ such that

• $r_{ik}^{(t)}=\frac{A_{ik}{e}^{\psi (n_k^{(t)})}}{\sum _{j=1}^K{A}_{ij}{e}^{\psi (n_j^{(t)})}}$

1. After convergence of $\mathit{n}$ and $\mathit{r}$, we can summarize the posterior distribution by taking its mean:

   • (10) ${\displaystyle \begin{array}{ll}E\left[Z_{ik}\right]& =r_{ik}\\ E\left[{\tau }_k\right]& =\frac{n_k}{\sum _{j=1}^K{n}_j}\end{array}}$

2. Finally, we perform two postprocessing steps on the posterior mean:

   1. If localization error is a parameter, we marginalize it out by projecting through that axis of the array.

   2. We adjust the posterior mean to account for defocalization biases, as described in ‘Defocalization’.

Throughout this article, we always report occupations for the SA model as the mean of $q(\mathit{\tau })$ according to Equation 10, twith localization error marginalized out and the appropriate defocalization correction applied.

Inference of the SA posterior works optimally with thousands to tens of thousands of trajectories. We also found it useful to have a cheap, dirty estimate for state occupations that can be evaluated on a small number of trajectories to visualize nuclei-to-nuclei variability (for instance, in Figure 5A).

For these purposes, we define the ‘naive occupation estimate’ ${\mathit{\tau }}_{\text{naive}}$ such that

Notice that this is just the posterior occupations based on the initial value for $\mathit{r}$ in the algorithm for SA inference. We use the same postprocessing steps for ${\mathit{\tau }}_{\text{naive}}$ as for SAs, including marginalizing out localization error and correcting for defocalization.

In the above section, we have left $p_X(x|{\mathit{\theta }}_k)$ unspecified as it depends on the type of motion being considered. This section states the form of $p_X(x|{\mathit{\theta }}_k)$ for RBME, the type of motion considered in this article.

Suppose that trajectory $i$ is constructed by measuring the position of a Brownian particle over sequential frame intervals of duration $\mathrm{\Delta }t$, and that each measured position has some error associated with it. We assume that this error is normally distributed with mean zero and variance ${\sigma }_{\text{loc}}^2$.

We refer to the change in the particle’s position over each frame interval as a ‘jump.’ If there are $L_i$ total jumps, let $\mathit{x},\mathit{y}\in R^{L_i}$ be the displacements of these jumps along the $x$ and $y$ axes, respectively. Then, the probability density over $\mathit{x}$ and $\mathit{y}$ is

where ${\displaystyle \mathbf{\Gamma }\in R^{L_i\times L_i}}$ is the covariance matrix defined by

where $D$ is the diffusion coefficient and ${\sigma }_{\text{loc}}^2$ is the localization error (Michalet and Berglund, 2012). Due to the contribution of the localization error to the off-diagonal terms of the covariance matrix, the jumps of an RBME are not a Markov process except when ${\sigma }_{\text{loc}}^2=0$.

The SA for RBME uses a 2D grid of diffusion coefficients and localization errors. In this grid, the diffusion coefficients $D$ are log-spaced between 10^-2 and 10² μm² s^-1, while the localization errors ${\sigma }_{\text{loc}}$ are linearly spaced between 0 and 0.06 μm.

In Figure 4—figure supplement 11 and Figure 4—figure supplement 12, we consider a generalization of RBME that we refer to as fractional Brownian motion with localization error (FBME). This is a simple modification of Mandelbrot and Van Ness’s FBM (Mandelbrot and Van Ness, 1968) that incorporates localization error.

We define 1D FBME as a mean-zero Gaussian process $X_t$ with the covariance function

where $S$ is the scaling coefficient, $H$ is the Hurst parameter (${\displaystyle 0<H<1}$), ${\sigma }_{\text{loc}}^2$ is the variance of the localization error, and $I_{t=s}$ is the indicator function (1 if $t=s$ and 0 otherwise). Because we always measure the position at regular frame intervals of duration $\mathrm{\Delta }t$, we let $t=i\mathrm{\Delta }t$ and $s=j\mathrm{\Delta }t$ so that this can be written as

The corresponding increment process ${\tilde{X}}_i=X_{i\mathrm{\Delta }t}-X_{(i-1)\mathrm{\Delta }t}$ is a mean-zero Gaussian process with the covariance function

2D and 3D FBMEs are constructed with independent 1D FBMEs along each spatial axis.

In Equation 13, the scaling coefficient has units of ${\displaystyle \mu {\mathrm{m}}^2{\mathrm{s}}^{-2\mathrm{H}}}$. As a result, its magnitude is highly dependent on $H$. Because we often want to parameterize the magnitude of the particle’s jumps separately from the covariance between jumps, in this article we use a ‘modified’ scaling parameter $\overline{S}$ defined by

As a result, the jump variance is $\text{Var}({\tilde{X}}_i)=2\overline{S}\mathrm{\Delta }t$, regardless of the Hurst parameter. While $\overline{S}$ is much easier to work with for one dataset, since it is dependent on $\mathrm{\Delta }t$ it must not be compared across datasets with different frame intervals and should first be converted to $S$ with Equation 14.

Here, we derive the SA posterior (Equation 9) by substituting model 5 into Equation 8 and imposing some additional physical constraints.

First, let $A_{ik}=p_X(X_i|{\mathit{\theta }}_k)$. Then factor $\log p(\mathit{F},\mathit{Z},\mathit{\tau })$ as

where the constant accounts for normalization factors. Plugging this into the second equation in Equation 8, we have

We have collected terms that do not depend on $\mathit{\tau }$ into the constant. In this article, we choose to weight the contribution of each trajectory to $\log q(\mathit{\tau })$ by the number of jumps in the trajectory. This is equivalent to treating jumps (rather than trajectories) as individual observations and is more robust to issues arising from the shallow observation depth of most sptPALM setups. It results in the modified equation

where $L_i$ is the number of jumps in trajectory $i$. We recognize this as a log Dirichlet distribution, so that

Next, we substitute Equation 15 into the first equation in Equation 8, giving

Since $q(\mathit{\tau })$ is the Dirichlet distribution given by Equation 16,

where $\psi (x)$ is the digamma function. Normalizing over the states for each trajectory $i$, we have

Together, Equations 16 and 17 constitute the result in Equation 9.

As mentioned above, a fundamental challenge with the finite state mixture (Equation 4) is determining the number of states. SAs deal with this issue by selecting a large, finite value for $K$ and relying on an inference routine that selects sparse subsets of states from a $K$-dimensional initial model.

DPMMs are more extreme, taking the limit $K\to \mathrm{\infty }$ (Ferguson, 1973). In this limit, the discrete vector of state occupations is replaced by a continuous distribution over the entire space of state parameters. The generative process for each trajectory is,

1. Randomly draw some state parameters ${\mathit{\theta }}_i\sim H$, where $H$ is a continuous distribution over the space of state parameters.

2. Randomly generate a trajectory $X$ from that state according to $p_X(x|{\mathit{\theta }}_i)$.

This process is formalized by replacing the Dirichlet distribution in Equation 4 with the Dirichlet process $\text{DP}(\alpha ,H)$, its infinite-dimensional analog. Here, α has the same function as in the finite mixture (defining the relative strength of the prior) and $H$ is the ‘base distribution’ over state parameters. The full DPMM is then

This corresponds to the second graphical model in Figure 2A. Each draw $G$ is a discrete probability distribution over part of the parameter space (Blackwell, 1973). This formalism is discussed in detail in Teh, 2010 or Neal, 1992. Here, we only remark that recovering the posterior $p(\mathit{\theta }|\mathit{F})$ requires marginalizing over $G$, yielding a continuous distribution over the parameter space.

To estimate the posterior distribution $p(\mathit{\theta }|\mathit{F})$, we take the Gibbs sampling approach introduced by Neal (Algorithm 8 in Neal, 2000). This involves sampling each ${\mathit{\theta }}_i$ while hold the other ${\mathit{\theta }}_{j\ne i}$ constant, yielding samples from the posterior distribution (Geman and Geman, 1984). To counter autocorrelation in the samples, Neal also endowed the sampler with additional Metropolis–Hastings nudges to the candidate state parameters between rounds of Gibbs sampling. For these nudges, we use a Gaussian proposal distribution.

In the case of RBME, the state parameters are $\mathit{\theta }=(D,{\sigma }_{\text{loc}}^2)$. Even with Neal’s sampler, a large number of samples are required to estimate the posterior over this 2D space, potentially requiring hours of computational time per dataset.

To make the problem more tractable, we replace this 2D space with a 1D approximation by neglecting the off-diagonal terms in the covariance matrix for RBME (Equation 12). With this approximation, Equation 12 can be rewritten as the log gamma density as

where $\varphi =\log \left[4(D\mathrm{\Delta }t+{\sigma }_{\text{loc}}^2)\right]$, $S_i$ is the sum of squared 2D jumps in trajectory $i$, and $L_i$ is the number of jumps. Notice that we cannot distinguish the contributions of $D$ and ${\sigma }_{\text{loc}}^2$ to $\varphi$ without measuring ${\sigma }_{\text{loc}}^2$ by some other method, such as averaging the negative sequential jump covariance across all trajectories in the dataset. This is the price we pay for a tractable DPMM and is the major disadvantage of this model (see, for instance, Figure 3A).

The complete Gibbs sampling routine for our DPMM is the following, which is essentially a modified version of Algorithm 8 from Neal, 2000:

1. Draw a random sample ${\mathit{\varphi }}^{(0)}=({\varphi }_1,\mathrm{\dots },{\varphi }_{m_0})$ from a uniform distribution on the interval $[{\varphi }_{\text{min}},{\varphi }_{\text{max}}]$, where the interval is selected to span the parameter space of interest. Each element of the vector ${\mathit{\varphi }}^{(0)}$ represents a candidate ‘state.’ At each iteration, we will add or remove states from this vector as the sampler explores the posterior.

2. Assign each trajectory $i$ to a state $k\in \{1,\mathrm{\dots },m_0\}$ with log probability proportional to $\log p_X(X_i|{\varphi }_k)$. Let this assignment be $Z_i^{(0)}$.

3. For each iteration $t=1,2,\mathrm{\dots }$

   1. For each trajectory $i=1,2,\mathrm{\dots }$, either set $Z_i^{(t)}$ to a state in the current set ${\mathit{\varphi }}^{(t-1)}$ with probability $(N-1)/(\alpha +N-1)$, or create a new state with probability $\alpha /(\alpha +N-1)$.

      1. If setting to an existing state, choose state $k$ with log probability proportional to $\log n_k+\log p_X(X_i|{\varphi }_k)$, where n_k is the number of jumps already assigned to state $k$.

      2. If creating a new state, pick m₀ values of $\varphi$ from the interval $[{\varphi }_{\text{min}},{\varphi }_{\text{max}}]$. Among these, accept a particular value ${\varphi }^{\prime }$ with log probability proportional to $\log p_X(X_i|{\varphi }^{\prime })$. Add a new state with this parameter to the set of current states ${\mathit{\varphi }}^{(t-1)}$.

   2. For each state $k$, if there are no trajectories currently assigned to it, remove it from consideration. Otherwise add it to ${\mathit{\varphi }}^{(t)}$, the next set of states, and update its parameter according to a Metropolis–Hastings step as follows:

      1. Propose a new ${\displaystyle {\varphi }^{\mathrm{\prime }}\sim \mathcal{N}\left({\varphi }_k^{(t-1)},{\nu }^2\right)}$ (resampling if ${\varphi }^{\prime }$ is outside the range $[{\varphi }_{\text{min}},{\varphi }_{\text{max}}]$).

      2. Evaluate the likelihood ratio

         • $r=\frac{\prod _{i=1}^N{p}_X(X_i|{\varphi }^{\mathrm{\prime }}{)}^{I_{Z_i=k}}}{\prod _{i=1}^N{p}_X(X_i|{\varphi }_k^{(t-1)}{)}^{I_{Z_i=k}}}\frac{\mathrm{\Phi }\left(\frac{{\varphi }_{\text{max}}-{\varphi }_k^{(t-1)}}{\nu }\right)-\mathrm{\Phi }\left(\frac{{\varphi }_{\text{min}}-{\varphi }_k^{(t-1)}}{\nu }\right)}{\mathrm{\Phi }\left(\frac{{\varphi }_{\text{max}}-{\varphi }^{\mathrm{\prime }}}{\nu }\right)-\mathrm{\Phi }\left(\frac{{\varphi }_{\text{min}}-{\varphi }^{\mathrm{\prime }}}{\nu }\right)}$

      3. Draw $u\sim \text{Uniform}(0,1)$. If ${\displaystyle r<u}$, set ${\varphi }_k^{(t)}={\varphi }^{\prime }$. Otherwise set ${\varphi }_k^{(t)}={\varphi }_k^{(t-1)}$.

4. The posterior mean can be estimated by making a histogram of the samples ${\mathit{\varphi }}^{(t)}$ weighted by their occupations ${\mathit{n}}^{(t)}$, where $n_k^{(t)}$ is the number of jumps assigned to the state with parameter ${\varphi }_k^{(t)}$ at iteration $t$.

5. Finally, we account for defocalization as discussed in ‘Defocalization.’

In this algorithm, $\mathrm{\Phi }(x)$ is the unit Gaussian CDF and its contribution to $r$ is required to make an unbiased proposal distribution for the Metropolis–Hastings updates given that $\varphi$ is confined to the range ${\displaystyle [{\varphi }_{\text{min}},{\varphi }_{\text{max}}]}$. ${\displaystyle I_{Z_i=k}}$ is the indicator function and is 1 if ${\displaystyle Z_i=k}$ and 0 otherwise.

While the gamma approximation 19 is what makes DPMMs computationally scalable, it also means that in order to disambiguate the contributions of diffusion and localization error to $\varphi$ we need to measure localization error by a different method. This is particularly relevant when accounting for defocalization, which relies on knowledge of $D$ independent of ${\sigma }_{\text{loc}}^2$. In this article, we always use the mean negative covariance between sequential jumps to estimate localization error prior to launching the Gibbs sampler above. However, this means that the DPMM is only as good as our estimate of ${\sigma }_{\text{loc}}^2$ – and as demonstrated in Figure 3 and Figure 3—figure supplement 1, our estimate of ${\sigma }_{\text{loc}}^2$ can be quite noisy with small numbers of trajectories and starts to fail completely when localization error varies a lot between states. SAs, although they require discretizing the parameter space, handle the problem of localization error in a more graceful manner than DPMMs.

We use ‘defocalization’ to refer to the axial movement of fluorescent emitters out of the microscope’s focus during an sptPALM acquisition. Because fluorescent emitters move quickly, defocalization is rapid and often limits trajectory length to a few frames. Due to defocalization, the probability to observe a jump from a fast-moving particle is less than that of a slow-moving particle because the jumps of a fast-moving particle are more likely to land outside the microscope’s focus.

Defocalization was considered as an experimental avenue to measure diffusion by Kues and Kubitscheck, 2002. In the jump histogram modeling frameworks of Mazza et al., 2012 and Hansen et al., 2018, who investigated the effect in detail, it appears as a correction term. The latter two sets of authors evaluated the defocalization probability by treating the microscope’s focal volume as a slab with absorbing boundaries and using the solution to the diffusion equation within these boundaries. Because the boundaries for the focal volume are not actually absorbing, both sets of authors then applied a correction term derived from Monte Carlo simulations of regular Brownian motion to ‘correct’ their correction.

Here, we provide a simpler alternative that is not based on Monte Carlo simulations, enables nonuniform probabilities of detection in the axial detection, and extends to a broader class of diffusion processes than regular Brownian motion. Although the framework can be extended to tracking with gaps, here we consider the case without gaps in tracking (all jumps are strictly between sequential frames).

Let $f(z,t=0)$ be the initial profile of particles in the axial direction of the microscope, and let $g(z,\mathrm{\Delta }t)$ be the Green’s function for the diffusion process at this frame interval. For regular Brownian motion, $g(z,\mathrm{\Delta }t)=e^{-z^2/4D\mathrm{\Delta }t}/\sqrt{4\pi D\mathrm{\Delta }t}$. Then the axial probability density for the particle after one frame interval can be obtained by convolving its initial profile with the Green’s function:

To account for defocalization, we multiply this density with an appropriate transmission function. For example, if our focal volume is a slab with depth $\mathrm{\Delta }z$, infinite XY extent, and perfect recall at any point inside the slab (i.e., all particles inside the slab are detected and no particles outside are detected), then our transmission function $T$ is

(This is the transmission function considered by Mazza et al., 2012 and Hansen et al., 2018.) The resulting axial profile is

To calculate the axial profile after $n$ frame intervals, we repeat this process iteratively:

where ${\text{Diffuse}}^{(n)}$ denotes $n$ sequential applications of the function

This scheme is illustrated in Figure 3—figure supplement 2A. The fraction of particles remaining in focus after $n$ frame intervals can be found by integrating this density:

In the SA and DPMM algorithms, we use this method to account for defocalization in the following way. Suppose that ${\tau }_k$ is the estimated occupation and $D_k$ is the estimated diffusion coefficient for state $k$. Then, we define the corrected state occupations ${\mathit{\tau }}^{\prime }$ such that

where ${\eta }_k$ is the probability for a Brownian motion to remain in focus after one frame interval and $\mathrm{\Delta }z$ is the focal depth. While defocalization can be incorporated explicitly into the models for SAs or DPMMs, in practice we find it makes little difference if it used as a final postprocessing step after inferring the posterior mean occupations.

To determine the focal depth $\mathrm{\Delta }z$, we used the method described in Hansen et al., 2017.

We have implemented SAs as a simple tool (https://github.com/alecheckert/saspt, Heckert, 2022c), available on the Python Package Index (PyPI) as https://pypi.org/project/saspt/. Documentation is also available at https://saspt.readthedocs.io/en/latest/.

DPMMs for SPT analysis have a publicly accessible implementation at https://github.com/alecheckert/dpsp, (copy archived at swh:1:rev:2f5196e4cae5943a5822be7c4493df50cd564a0c, Heckert, 2022a). As a result of the investigation in this article, we recommend that researchers looking to try these methods start with SAs due to their superior performance.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

We thank Luke Lavis for generously synthesizing the Janelia Fluor dyes used in these experiments, Thomas Graham for his creative suggestions about the defocalization problem, Ana Robles for the monumental task of keeping our SPT microscopes in working order, and Anders Hansen and Maxime Woringer for insights and advice at the outset of the project. Claudia Cattoglio provided indispensable advice on Western blots. Portions of this work were performed on shared instrumentation at the UC Berkeley Cancer Research Laboratory Molecular Imaging Center (MIC), supported by The Gordon and Betty Moore Foundation. We thank MIC gurus Holly Aaron and Feather Ives for their assistance with shared microscopes. Sanger sequencing was performed at the UC Berkeley DNA Sequencing Facility. This work was supported by NIH grant 1U54CA231641 (XD) and the Howard Hughes Medical Institute (CC34430, RT). AH was supported by the NIH Stem Cell Biological Engineering predoctoral fellowship T32 GM098218. We would like to thank David Schaffer for his support and feedback throughout the project as a program director for the T32 fellowship.

© 2022, Heckert et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.