A Hidden Markov Model for Detecting Confinement in Single-Particle Tracking Trajectories.

State-of-the-art single-particle tracking (SPT) techniques can generate long trajectories with high temporal and spatial resolution. This offers the possibility of mechanistically interpreting particle movements and behavior in membranes. To this end, a number of statistical techniques have been developed that partition SPT trajectories into states with distinct diffusion signatures, allowing a statistical analysis of diffusion state dynamics and switching behavior. Here, we develop a confinement model, within a hidden Markov framework, that switches between phases of free diffusion and confinement in a harmonic potential well. By using a Markov chain Monte Carlo algorithm to fit this model, automated partitioning of individual SPT trajectories into these two phases is achieved, which allows us to analyze confinement events. We demonstrate the utility of this algorithm on a previously published interferometric scattering microscopy data set, in which gold-nanoparticle-tagged ganglioside GM1 lipids were tracked in model membranes. We performed a comprehensive analysis of confinement events, demonstrating that there is heterogeneity in the lifetime, shape, and size of events, with confinement size and shape being highly conserved within trajectories. Our observations suggest that heterogeneity in confinement events is caused by both individual nanoparticle characteristics and the binding-site environment. The individual nanoparticle heterogeneity ultimately limits the ability of interferometric scattering microscopy to resolve molecule dynamics to the order of the tag size; homogeneous tags could potentially allow the resolution to be taken below this limit by deconvolution methods. In a wider context, the presented harmonic potential well confinement model has the potential to detect and characterize a wide variety of biological phenomena, such as hop diffusion, receptor clustering, and lipid rafts.

Introduction

Single-particle tracking (SPT) experiments directly observe the motion of single molecules and hence offer a powerful method to analyze the membrane environment. For instance, detection and characterization of heterogenous diffusion behaviors yields information on membrane structure (1, 2). However, SPT methods require the molecule of interest to be tagged with a trackable label that is imaged over a number of time steps. A number of experimental design limitations constrain the amount of information that can be extracted from such data, including spatial accuracy, temporal resolution, and the tracking period. New technologies are capable of extending the trajectory length while retaining high sampling rates and high spatial resolution. For example, interferometric scattering microscopy (iSCAT) can generate very long (50,000 step) trajectories, with high spatial (<2 nm) and temporal (up to 500 kHz) resolution (3, 4, 5, 6). However, a fundamental problem that impacts interpretation is the effect of the tag itself (7). This is particularly relevant for iSCAT because the gold nanoparticle (AuNP) tags are 20–40 nm in diameter, whereas spatial resolution is estimated to be ∼2 nm for a 20-nm AuNP (4, 5); relative movements between the AuNP and the bound GM1 will thus convolve with the movement of the GM1. For example, an iSCAT study on model membranes demonstrated both Gaussian-like and ring-like confinement events, which was ascribed to transient multivalent binding of the tag (4). Thus, to extend this technique to in vivo experiments, there is a need to deconvolve the tag signature from the environment signal. Failure to achieve this separation means that interpretation of the high-resolution dynamics measured by these techniques may be limited to the order of the tag’s size.

Analysis of SPT data is not straightforward primarily because of the stochastic nature of diffusion. This has led to the development of a range of statistical methods that detect deviations from Brownian motion, such as mean-square displacement (MSD) (8, 9, 10, 11, 12, 13) and confinement (14, 15, 16, 17, 18, 19) analyses. A new breed of methods model switching of the movement dynamics between various dynamic states (20, 21, 22, 23, 24), often within a hidden Markov chain framework (25, 26, 27, 28, 29, 30, 31). For high-resolution data, the latter techniques can utilize the high level of information present in the trajectory to extract detailed motion characteristics and potentially infer underlying biophysical mechanisms.

However, the majority of existing hidden Markov approaches only incorporate changes in the particle diffusivity and (or) drift. Such methods can only approximate confinement through a change in the particle’s effective diffusion coefficient. To our knowledge, the only exception is the work of Bernstein et al. (31) who explicitly model confinement using a hidden Markov model (HMM) within a maximal likelihood framework. In this article, we develop an HMM harmonic potential well (HPW) confinement analysis method using a Bayesian approach. Specifically, the particle moves between two states hidden to the observer: free diffusion with (to be determined) diffusion coefficient D and confinement in an HPW (center location and well strength to be determined). We developed a Markov chain Monte Carlo (MCMC) algorithm to infer model parameters and hidden states from a single trajectory. We tested the algorithm on simulated data, then applied it to previously published experimental iSCAT trajectories of GM1 lipids diffusing in model membranes (4). Specifically, a (20 or 40 nm) AuNP was coated in cholera toxin B subunits (CTxBs) by streptavidin binding, with each CTxB then binding 5 GM1 molecules in the lipid membrane to form an AuNP/CTxB/GM1 complex. In trajectories of 20 nm AuNP/CTxB/GM1 diffusing in model membranes on a glass substrate, we detected clear periods of trapping in wells of a mean radius of 18 nm with a mean trapping time of 0.024 s. However, we also observed inherent heterogeneities in both AuNP/CTxB/GM1 particles and trapping sites, which ultimately affect trajectory characteristics.

This article is organized as follows. In Methods, we introduce the HPW confinement model and an associated inference (MCMC) algorithm. The full derivation of the MCMC algorithm is described in Note S1 of the Supporting Materials and Methods. In Results, we demonstrate accurate inference of model parameters and hidden states on simulated trajectories, then apply the algorithm to iSCAT trajectories of AuNP/CTxB/GM1 diffusing in model membranes.

Methods

We implemented the following methods in MATLAB (The MathWorks, Natick, MA). The source code, documentation, trajectory data, and working examples are freely available (https://doi.org/10.5281/zenodo.1405647).

HPW model

We developed a model for a particle that switches between a freely diffusing state and a confinement state localized around a slowly diffusing center. The state is encoded by a hidden variable z, with if the particle is freely diffusing at time and if confined, where denotes the time point (i.e., frame). The state depends only on with transition probabilities (constant frame rate):where and are the per-frame probabilities of switching into and out of confinement, respectively. The probability of being in state given state is therefore as follows:where Bernoulli denotes the Bernoulli probability distribution with variable x and parameter p. In the free state, the particle diffuses freely with diffusion coefficient D. In the confined state, the particle is assumed to have a directed component to its diffusive motion, proportional to the distance from the well center , i.e., the force is proportional to where is the particle position at time . (Note that and are two-dimensional (2D) vectors). During confinement, the center diffuses much slower than the particle itself (diffusion coefficient ). When the particle is free, C diffuses with diffusion coefficient , where is sufficiently high that the center can relocate between different confinement sites. The center is thus still present even when it is not affecting the particle. The stochastic differential equations for this model are as follows:where are independent Weiner processes. During confinement, has Ornstein-Uhlenbeck (OU) dynamics with center . We assume that switching can only occur at the sampling points. We also assume that is slowly varying and therefore ignore its time dependence over the time step . The frame-to-frame dynamics are hence,

See Note S1 of the Supporting Materials and Methods for full details. If the step size is sufficiently small relative to the confinement strength , an Euler-Maruyama approximation is justified, but if the particle explores the well over , this OU solution is required. We refer to this discrete-time stochastic model as the HPW confinement model.

The model has two hidden states to be inferred at all trajectory time points : the state (confined or free) and the position of the HPW center when confined. There are also five parameters to be inferred: two diffusion coefficients (D and ), the strength of the HPW (κ), and two transition probabilities ( and ). is treated separately because it only weakly affects the trajectory and does not affect the likelihood or parameter estimates provided it is sufficiently high. Fig. 1 A shows a simulated HPW model trajectory. In the simulation, we include a drift term for the center so that C tracks X when not confined. This ensures that C is close to X when the particle switches from free diffusion to confinement, and therefore confinement zones remain within a reasonable field of view. This tracking of X by C is not included in the inference algorithm because diffusion alone is sufficient to allow the Markov chain to find high-probability paths.

Figure 1

Simulated harmonic potential well (HPW) model trajectory. (A) A simulated trajectory colored by state. Model parameters are , , , , , , time step s, and frames. The simulation was performed using Eq. 5 and a modified version of Eq. 6 as discussed in the main text. Trajectory colored blue, denoting free diffusion or yellow, denoting the confined state. Color bar length, 0.1 μm. (B) A schematic of AuNP/CTxB/GM1 complex in 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) lipid bilayer, based on figure 5 in (4).

MCMC sampler

There are a number of MCMC samplers for linear switching models in the literature; the main distinction is whether variables are integrated out using an inverse Wishart prior (32), or a Markov chain incorporating all variables is used. The latter approach allows greater control of prior information, including use of uninformative priors, whereas the Wishart distribution, motivated by computational convenience, imposes a dependence between variable correlations and scale, which is a concern for inference (33). We developed an MCMC algorithm (Supporting Materials and Methods, Note S1) for the full system of variables to fit the HPW model (Eqs. 5 and 6) to 2D trajectory data, . We chose uninformative priors for all parameters except for the transition probabilities, where we use an informative prior to restrict rapid switching between states (details in Note S1 of the Supporting Materials and Methods). For an SPT trajectory, the algorithm samples the posterior distribution, , giving K samples of the parameters and hidden states . Here, for each sample k, and are the set of hidden states and center locations (2D vectors) throughout the trajectory.

We determined convergence of the MCMC sampler by calculating the Gelman potential scale reduction factor (PSRF) (34), considering a run converged provided the PSRF was below a threshold in all variables, set to 1.2 on experimental trajectories. The MCMC run length was increased up to a maximum of 4 × 105 steps on trajectories that failed the convergence criteria on shorter runs.

GM1 molecules diffusing in model membranes

We applied the MCMC algorithm to previously published iSCAT SPT data (4), where CTxB-coated AuNPs were introduced to a 1,2-dioleoyl-sn-glycero-3-phosphocholine lipid bilayer containing 0.03% GM1 lipids (Fig. 1 B). A confinement event corresponds to an interaction between an AuNP/CTxB/GM1 complex on the upper leaflet with a lower leaflet GM1 that is immobilized on a glass surface. This was previously referred to as “interleaflet coupling and molecular pinning” (4). Both Gaussian and non-Gaussian confinement events were observed; we investigated these events in greater detail using our HPW model.

The data set includes 71 trajectories of 20-nm AuNP/CTxB/GM1 diffusing in a model membrane on a glass substrate and 18 trajectories of 40-nm AuNP/CTxB/GM1 in a model membrane on a mica substrate. There is a dynamic error in the localization accuracy at the 50-kHz sampling rate resulting in apparent superdiffusive behavior, which we removed by subsampling down to 5 kHz (Figs. S1 and S2; Note S2 of the Supporting Materials and Methods). We also removed trajectory artifacts because of multiple AuNPs in the focal area (Supporting Materials and Methods, Note S2). The MCMC algorithm did not converge on five 20-nm AuNP/CTxB/GM1 trajectories (PSRF convergence criteria of 1.2) leaving a set of 66 trajectories for further analysis. MCMC runs on all 18 40-nm AuNP/CTxB/GM1 on mica trajectories converged.

Thresholding hidden states for lifetime analysis

For each trajectory, at each time point i, we computed the probability of confinement from the MCMC posterior distribution samples. This probability distribution is concentrated near 0 and 1 on experimental trajectories (only 2.7% of the state probabilities were between 0.2 and 0.8; Fig. S3), indicating high confidence in confinement state estimates. To annotate the trajectory by state, we define the binary signal, or , for free diffusion and confinement respectively, using a threshold of 0.5 on . We then identify confinement events as a series of ones in the (posterior) binary state vector and free diffusions as a series of zeros, allowing event lifetimes (and per-event spatial statistics) to be computed. When considering event lifetimes, we exclude those containing either the first or last time point of the trajectory, because the full event is not witnessed, hence the state lifetime is unknown.

Confinement event profiling

To analyze confinement events in 20-nm AuNP/CTxB/GM1 trajectories, we utilized spatial statistics (including the mean confinement radius and radial skewness, defined in Table 1) based on the Euclidean distance between the particle and the confinement center. We calculate these statistics for all events of at least 0.01 s (50 frames). Unlike the event lifetime analysis, we allow events that contain either the first or last (or both) time points. Furthermore, we compute statistics including events revisiting a previous confinement zone where applicable (details in Table 1). These restrictions left a set of 271 confinement events when excluding repeat events and 427 when including them. The number of events within a trajectory ranges from 1 (there were six examples in which the particle remained trapped for the entire trajectory) to 11 (without repeats) or 25 (with repeats).

Table 1

Calculation of Confinement Event Statistics

Statistic	Calculation
Confinement radius	R_lm={R_i}_i∈T_lm,R_i=‖X_i−C¯_lm‖
Mean confinement radius	R¯_lm=1M_lm∑i∈T_lmR_i
Radial skewness	S_lm=1M_lm∑_i∈T_lm(R_i−R¯_lm)^3[1M_lm∑_i∈T_lm(R_i−R¯_lm)^2]^3
Radial mean-median distance	|R¯_lm−Rˆ_lm|
Radial SD	var[R_lm]

We denote the time points of the mth trapping event in the lth trajectory . Events have associated particle positions and harmonic well center positions . The mean posterior harmonic well center is given by , where is the number of time points in . To remove events that revisit a previous trapping zone, we did not include events if was within 30 nm of a previous confinement center (, ) within trajectory l. denotes the Euclidean distance, and denotes the median.

Results

MCMC on simulated data

The HPW model sampler was extensively tested on simulated data. Figs. 2 and 3 show an MCMC run on the simulated trajectory of Fig. 1 A. The parameter posteriors are consistent with the true (i.e., simulation) values (Fig. 2). When confined, the inferred center closely tracks the simulated center (Fig. 3, A and B), and every confinement event is accurately inferred (Fig. 3 C). The inferred model parameters are independent of the algorithm parameter in both the noiseless case (Fig. S5) and in the presence of static localization error at the same level as the 20-nm AuNP/CTxB/GM1 trajectories (Fig. S6); are typically underestimated because of their informative priors. There also appears to be a small underestimation in the confinement strength κ in the presence of localization error, which we further investigated by increasing the localization error’s SD (Fig. S7). This revealed a clear trend, with an increase in confinement strength underestimation bias with localization error. However, the effect is relatively small at the noise level (SD 2.7 nm) of the experimental data presented in this article. Performance was robust to changes in trajectory length and number of events (Figs. S8 and S9). Estimation of D, , κ, , and switching rates were robust to the time series subsampling rate (Figs. S10 and S11); in particular, most events were still detected even with a 10-fold subsampling (Fig. S11). Parameter estimation is also robust to changes in confinement strength κ (Fig. S12). However, we found that confinement strength estimation accuracy decreases dramatically when the confinement center diffusivity approaches 10% of the free diffusivity D (Fig. S13), reflecting the degeneracy of the problem when .

Figure 2

Posterior parameter distributions of the HPW model for a simulated trajectory. (A) The posterior distribution for D (blue, solid line) and (red, dashed line), with simulation values indicated (circles) (B) The posterior for κ and simulation value (circle). (C) The posterior for (blue, solid line) and (red, dashed line), with simulation values (cross, circle respectively). The trajectory is as in Fig. 1. MCMC priors are as in Note S1 of the Supporting Materials and Methods. Corresponding MCMC runs are shown in Fig. S4. Data are based on the pooling of five independent chains of 2000 steps with a 1000-step burn-in. MCMC priors are as in Note S1 of the Supporting Materials and Methods. To see this figure in color, go online.

Figure 3

Hidden state inference for the HPW model for a simulated trajectory. (A and B) The mean inferred position of the harmonic potential center in x and y directions (black) and simulated (true) center (red). The colored line at the top represents the following particle states: free diffusion (blue) and confinement (yellow). (C) The inferred confinement probability (black line) and simulated (true) confinement state (yellow area). (D) The trajectory is colored by mean inferred confinement state, from (blue, free) to (yellow, confined). Color bar length, 0.1 μm. MCMC is as in Fig. 2.

MCMC on 20-nm AuNP/CTxB/GM1 on glass trajectories

An example of the model fit is shown in Fig. 4, with the segmented trajectory shown in Fig. 4 E. Fig. 5 shows the associated parameter posterior estimates. Particle state (confined or free diffusion) is well determined, with state probabilities near zero or one (Fig. 4 C). The parameter estimates for D and κ are tight (low relative SD), whilst the diffusion coefficient of the center is very low, (mean ± SD) compared to , indicating near complete immobilization of the well. The inferred position of the well center is also practically stationary in both coordinates during periods of confinement consistent with immobilization (Fig. 4, A and B). As an independent measure of changes in mobility, we estimated the effective local diffusion coefficient (Fig. 4 D), which demonstrates a clear shift at around 0.5 s (i.e., the first inferred switch point). By color coding the trajectory according to the probability of being confined per frame (Fig. 4 E), we can extract periods of confinement with non-Gaussian occupation profiles (Fig. 4, F and H). In this trajectory, we observed that one confinement zone is visited twice (Fig. S15) and that the repeat confinements had remarkably similar occupation profiles (Fig. 4, G and H). The probability per frame of switching is reasonably well inferred (Fig. 5 C) despite the small number of events. The probability of escape from a confinement zone is smaller than the probability of trapping, reflecting the short periods of time that the AuNP/CTxB/GM1 complex undergoes free diffusion.

Figure 4

Hidden state inference for the HPW model applied to a 20-nm AuNP/CTxB/GM1 trajectory. (A and B) The mean inferred position of the HPW center C (x, y components) and upper colored bar representing , (color scale goes from (blue, free) to (yellow, confined)). (C) The inferred mean confinement state and (D) moving average of local maximal likelihood diffusion coefficient estimate (window size: 100 subsampled frames). (E) The trajectory colored by mean inferred confinement state (color bar length, 0.1 μm). (F–H) Density-colored 2D spatial histograms of confinement events. The two events in (G) and (H) are spatially colocated. Data is based on the pooling of 10 independent chains. MCMC priors and convergence criteria are as in Note S1 of the Supporting Materials and Methods.

Figure 5

Posterior parameters of the HPW model applied to a 20-nm AuNP/CTxB/GM1 trajectory. (A) Posterior distributions for D (blue, solid line) and (red, dashed line). (B) The posterior for κ. (C) The posteriors for (blue, solid line) and (red, dashed line). Distributions consist of samples pooled from 12 independent runs. Corresponding MCMC chains are shown in Fig. S14. Trajectory and MCMC runs are as in Fig. 4. To see this figure in color, go online.

Applying our MCMC algorithm to the 66 trajectories, we obtain parameter estimates across the population (Fig. S16). The mean value of D over all trajectories was 1.15 ± 0.106 μm2s−1 (mean ± standard error (SE); population SD 0.86 μm2s−1); MSD analysis (using the @msdanalyzer package (35)) gave a smaller estimate, 0.0525 ± 0.017 μm2s−1 (mean ± SE). This difference reflects the fact that MSD does not account for confinement, which is the dominant state, whereas our diffusion coefficient estimate does. Our Bayesian analysis provides estimates of parameter confidence per trajectory, which are in fact substantially smaller than the spread between trajectories (Fig. S17); specifically, the ratio of the population variances of D and κ are 257 and 59 times larger than the average trajectory posterior variances, respectively. This indicates the presence of system variability, giving rise to trajectory heterogeneity. To understand its cause, we investigate whether heterogeneity is manifest in the confinement events of individual trajectories, specifically the size, shape, and lifetime of these events.

Lifetime and shape analysis of confinement events

The mean confinement state lifetime (as defined in Methods) is 0.024 s, but there is a large variation in event lifetimes across trajectories (Fig. 6, A and B) and significant heterogeneity across trajectories (p = 0.02, Kruskal-Wallis test, 1779 events across 60 trajectories). The lifetimes of free-diffusion events (mean 0.002 s) did not show significant heterogeneity across trajectories, (p = 0.86, Kruskal-Wallis test on 60 trajectories, 1770 events). Further, we examined if the population of lifetimes across trajectories conform to an exponential waiting time model, i.e., whether switching between states obey first-order kinetics. A quantile-quantile plot demonstrates that there is a distinct deviation from an exponential distribution fit (mean event time μ = 0.024 s); specifically, there are a far higher proportion of longer trapping events, indicative of heterogeneity. A mixture of two exponentials is a better fit (Fig. 6 D) suggesting that the confinement events derive from a heterogeneous population with at least two components with short and long average lifetimes. The minor population of long lifetime events are dispersed over trajectories (Fig. 6 B); in particular, trajectories are not split into two groups with long and short mean confinement times. In contrast, the free diffusion state lifetimes closely follow a mono-exponential distribution (Fig. 6, E and F).

Figure 6

Confinement event lifetimes are not exponentially distributed. (A) A histogram of all confinement lifetimes (n = 1959 events). (B) A scatterplot of confinement lifetimes against trajectories ordered by mean confinement lifetimes. (C–F) Quantile-quantile plots of state lifetimes against exponential fits. (C) Confinement events against the exponential distribution ( s, ) and (D) confinement events against samples () from a mixture of two exponentials ( s, s; weights 0.80 and 0.20, respectively; ) are shown. (E) Free diffusion lifetimes (n = 2011 events) against the exponential distribution ( s, ) and (F) free diffusion lifetimes against samples from a mixture of two exponentials ( s, s; weights 0.99 and 0.01, respectively; ) are shown. The red line is an extrapolated linear fit to the first and third quantiles. Plots include all confinement events except those that contained the trajectories’ first or last time point. To see this figure in color, go online.

We next analyzed confinement event shape using the spatial statistics defined in Table 1. The mean confinement radius over all trajectories is 18 nm, comparable to the size of the AuNP, although estimator inflation is likely to be present (36). The mean radial skewness is 0.88; for comparison, a 2D Gaussian distribution gives a radial displacement (from the mean) that is Rayleigh distributed with skew . Mean confinement radius and radial skewness show a wide distribution of values across confinement events (Fig. 7, A and B), with significant (one-way analysis of variance: mean confinement radius ; radial skewness ; 271 events grouped by 66 trajectories) heterogeneity across trajectories (Fig. 7, C and D). Confinement event spatial histograms for all 66 20-nm AuNP/CTxB/GM1 trajectories, ordered by the average within-trajectory mean confinement radius (Fig. S18) and average radial skewness (Fig. S19) demonstrate the wide variety of confinement shapes.

Figure 7

Shape statistics for confinement events in 20-nm AuNP/CTxB/GM1 trajectories. (A and B) Histograms over confinement events. (C and D) Spatial statistics for all confinement events, ordered by the average within trajectory statistic. Plots include all confinement events of at least 0.01 s, with events revisiting a previous trapping zone removed (giving 271 events). (E and F) Scatterplots of state lifetime against the given spatial statistic for the same 271 confinement events. To see this figure in color, go online.

The observed heterogeneity in confinement time, size, and (previously reported (4)) shape raises two key questions:

Does the shape of confinement events determine their lifetime?

Does heterogeneity predominately arise from a mechanism operating at individual confinement sites (local environment dependent) or at the level of trajectories (AuNP/CTxB/GM1 nanoparticle dependent)?

To probe the relationship between confinement event shape and lifetime (the first question), we examined their correlation. We found no correlation between confinement state lifetime and mean confinement radius (Fig. 7 E; Table S1) but a weak (and significant) negative correlation between lifetime and radial skewness (Fig. 7 F; Table S1). This suggests that the mixed-exponential nature of the binding lifetime is only weakly related to the shape of the binding event, i.e., these arise from different physical mechanisms.

Regarding the second question, the heterogeneity analysis above (see also Fig. 7, C and D) indicates that confinement events are statistically more similar within trajectories than across trajectories. Additionally, the ratio of the mean variance within trajectories to the variance across all events is 0.6 for both mean confinement radius and radial skewness (Table S2). To determine which confinement statistic is most strongly conserved within trajectories, we clustered events by each confinement statistic and quantified the similarity of events within single trajectories (Fig. 8). Confinement size is the most conserved, followed by confinement lifetime when excluding events revisiting a previous confinement zone (Fig. 8 A). Incorporating revisiting events dramatically improves the conservation of confinement size and shape statistics relative to lifetime (Fig. 8 B); this suggests that, although shape is conserved, confinement time is variable between events at the same location. This shape conservation at the same site is evident from the confinement event spatial histograms for long events (Fig. 9). Of note, the mean-median distance statistic failed to show significant heterogeneity across trajectories (p = 0.06, one-way analysis of variance) reflecting its lack of conservation when excluding events revisiting a previous confinement zone (Fig. 8 A) but was as conserved as lifetime when all events were included (Fig. 8 B). Thus, in answer to the second question, heterogeneity arises at both the trajectory- and confinement-site level, with different effects on confinement lifetime and shape. This suggest that nanoparticle confinement events are described by two degrees of freedom.

Figure 8

Clustering of confinement event statistics across individual confinement events were clustered (k-means++ algorithm (42) with squared Euclidean distance metric) based on event statistics. For each trajectory, l, the Shannon diversity index, , was calculated ( is the proportion of the events in trajectory l that appeared in cluster j). The sum of the Shannon diversity index over all trajectories is then a measure of the dissimilarity of events within trajectories (the lower the Shannon index the higher the similarity). The event statistics (shown in the legend) are defined in Table 1. For each choice of , 50 separate clusterings were performed (because the k-means++ algorithm stochastically assigns initial values for cluster centroids), and the sum of the Shannon diversity index was averaged over these clusterings. (A) The clustering of events was obtained as described in “Confinement Event Profiling” in the main text, except with events containing the first or last time points excluded (214 events total). (B) The clustering of events is the same as (A), except with events that revisited a previous confinement zone included. To see this figure in color, go online.

Figure 9

Spatial conservation of confinement events revisiting the same site. Each column shows particle position histograms for two spatially colocated confinement events and one event at a different location in the same trajectory. The spatially colocated confinement events are distinct, i.e., the particle moved away from the trapping zone between the displayed events. Each plot has a side length of 0.1 μm.

Analysis of 40-nm AuNP/CTxB/GM1 trajectories on mica

As a control, we analyzed 18 trajectories of 40-nm AuNP/CTxB/GM1 diffusing in supported lipid bilayers (SLBs) on a mica substrate. The previous analysis demonstrated that no confinement for this treatment was present (4). We applied our HPW model MCMC algorithm to this data and detected no confinement events (Fig. S20); the posterior confinement probability was <0.01 for all time, in all trajectories. The mean D was 1.2048 ± 0.09 μm2s−1 (mean ± SE), comparable to 20-nm AuNP/CTxB/GM1 on glass trajectories (1.15 ± 0.20 μm2s−1). The mean MSD-derived (with @msdanalyzer (35)) D was 0.87 ± 0.12 μm2s−1. These values are in closer agreement than for the 20-nm AuNP/CTxB/GM1 on glass data set, which is expected because of the lack of confinement. However, we again observed trajectory heterogeneity in the diffusion coefficients with a ratio of population variance/mean trajectory variance of 123, indicative of individual AuNP-dependent diffusion coefficients.

Discussion

We developed a Bayesian algorithm to infer an HPW confinement HMM and used it to partition SPT trajectories into periods of free diffusion and confinement. When applied to experimental AuNP/CTxB/GM1 trajectories, we detected clear periods of confinement and free diffusion (Fig. 4). It was previously proposed that confinement event shape heterogeneity (Gaussian versus non-Gaussian confinement) in this data set was due to transient multivalent binding of the tag (4). Our analysis of confinement events attained using the HPW model revealed the following heterogeneity trends:

Confinement size and shape are conserved within trajectories (Fig. 7, C and D), and repeat events at the same site show similarities (Figs. 5, G and H and 9).

Confinement event lifetimes are heterogeneous across trajectories and comprise a mixture of at least two exponentials with short (4 ms) and long (100 ms) mean lifetimes (Fig. 6, C and D).

Spatial heterogeneity and lifetime heterogeneity are effectively uncorrelated, suggesting they arise from different mechanisms.

Based on these observations, we propose a refinement to the transient multivalent tag-binding hypothesis. Namely, the characteristics of individual confinement events are determined by the following two factors: the size and geometry of the GM1 platform on the lower leaflet (determining residence times) and the number and distribution of CTxB complexes bound to the surface of the AuNP (determining size and shape of confinement event) (Fig. 1 B).

These dependencies are consistent with CTxBs remaining attached to the surface via GM1s throughout the entire trajectory (Fig. 10); this is supported by the high affinity of the CTxB/GM1 bond with a dissociation rate in SLBs of s−1, giving a mean binding lifetime of s (37). We propose that differences in the geometry of bound CTxB on the surface of the nanoparticle causes trajectory-conserved variation in the observed confinement as follows: tightly packed (or single) CTxBs have more freedom to “wobble” (Fig. 10 B), and broadly spaced, multiple (bound) CTxBs have less freedom (Fig. 10 C), giving a large, respectively small confinement radius for binding events. Additionally, non-Gaussian confinement events occur when there is a second (or potentially multiple) CTxB/GM1 attachment that is not immobilized, which restricts movement to a rotation or nonuniform “wobbling” around the immobilized binding site (Fig. 10 D). These hypotheses are consistent with the fact that there are around 25 CTxBs per 20-nm AuNP (4); it is expected that there will be variability in both their number and spatial distribution. We observe a strong correlation of the diffusion coefficient with the mean confinement radius (; Fig. 11 A), but not with the confinement event lifetime (; Fig. 11 B). This is consistent with the hypothesis that with more attachments, the AuNP experiences higher drag, whereas the mean confinement radius decreases because of stronger geometric constraints. Our analysis thus suggests that variation in the number and spatial configuration of bound CTxB contributes to the nature of the AuNP interaction with the upper leaflet of the bilayer, thereby giving each individual AuNP/CTxB/GM1 complex a confinement signature (Figs. 8, S18, and S19) and diffusion coefficient, with the latter also being evident in the confinement-free 40-nm AuNP data.

Figure 10

Schematic of AuNP/CTxB/GM1 structures leading to Gaussian and non-Gaussian confinement profiles. (A) Free diffusion, (B) wide Gaussian-like confinement, (C) narrow Gaussian-like confinement, and (D) non-Gaussian confinement. Insets in (B)–(D) are example histograms of particle positions pooled over confinement events within selected trajectories (e.g., Figs. S18 and S19). Insets have a side length of 0.1 μm. The schematic is based on a figure in (4). To see this figure in color, go online.

Figure 11

Correlation between D and confinement statistics. (A) The mean confinement radius and (B) confinement lifetime. Plots include all confinement events of at least 0.01 s. To see this figure in color, go online.

The characteristics of the lower-leaflet GM1 platform contribute a confinement site dependence, with conservation of shape and size upon revisiting the same site (Figs. 8 and 9). On the other hand, confinement lifetime at the same site is variable. The GM1 in the lower leaflet is immobilized by hydroxyl pinning sites on the glass surface, with these sites having an estimated size of <10 nm (4). However, aggregation of GM1 with domain sizes of 15–60 nm in SLBs has been observed in atomic force microscopy experiments (38). Large sites consist of more aggregated GM1 in the lower leaflet. Our mean confinement radius is 18 nm, which would comprise both AuNP/GM1/CTxB nanoparticle degrees of freedom around the binding site and displacements of the GM1 platforms between the leaflets. This suggests that either pinning sites are small, or no relative movement is possible. Larger pinning sites may trap multiple CTxB molecules on the AuNP, leading to more Gaussian behavior as rotational degrees of freedom are lost and possibly longer (on average) trapping times. This could be the cause of the negative correlation between non-Gaussian confinement shape and event lifetime (the lower the radial skewness statistic, the longer the typical confinement time (Fig. 7 F)). However, the double exponential mixture distribution of confinement lifetimes cannot be explained by these mechanisms. Mean lifetime does not partition by trajectories (with long event lifetimes being distributed throughout the trajectories (Fig. 6 B)) suggesting a random process is responsible. The simplest explanation is that binding of the AuNP/CTxB/GM1 at the pinning sites is heterogeneous, e.g. there could be a multistep binding sequence with the second (long lifetime) step proceeding in only a fraction of the binding events. We note that the six trajectories that remain confined throughout are a third population, because even on the long-lifetime distribution, observing binding events of 1 s or longer is negligible (probability ).

Outlook and future work

Analysis of SPT trajectories with HMMs has advantages over other methods for detecting confinement in single trajectories. In particular, they do not rely on tuning algorithm parameters through a comparison with Brownian motion. Additional parameters (such as the confinement strength κ, center C, and switching times as inferred here) can also be extracted, which allows for interpretation and comparison of confinement event characteristics across and within trajectories. However, appropriate HMMs are necessary for successful analysis. Specifically, models must approximate well the behavior of different dynamic states in the data. For instance, confinement is often associated with a decrease in the effective diffusion coefficient, suggesting that models that switch diffusivities (25, 26, 27, 28, 29, 30) should also be able to detect confinement in these iSCAT particle trajectories. However, we found that a two-state diffusion coefficient switching HMM (30) could not segment these trajectories (data not shown). This implies that the effective diffusion coefficient of the AuNP/CTxB/GM1 complex does not change sufficiently during confinement events. The effective diffusion coefficient under confined Brownian motion is dependent on the temporal and spatial resolution of the data, suggesting that the high temporal sampling rate and positional accuracy of iSCAT data do not reduce the effective diffusion coefficient during confinement. In fact, the sampling timescale leads to frame-to-frame displacements that are on the order of the confinement radius, which required us to use OU dynamics in our MCMC algorithm. Failing to account for the effect of the frame-to-frame displacement/confinement size ratio on the integration scheme accuracy leads to a bias in the parameter estimates (data not shown).

This suggests that specifically modeling confinement, rather than modelling as a change in the diffusion coefficient, is necessary for data such as this. We believe the only other published confinement HMM for SPT analysis is by Bernstein et al. (31). Our models are similar, with both switching between free diffusion and confinement in a HPW. However, there are key differences. Firstly, we incorporate diffusion of the harmonic well to relax the constraint of a circular potential. The importance of this will depend on the spatial-temporal resolution of the data and whether the confinement zone is static or has time-dependent shape variation or drift. For instance, transient confinement in lipid rafts, which both diffuse and have an irregular shape, has been hypothesized. Secondly, the inference frameworks are different; we use a Bayesian approach to determine full posterior distributions for parameters and hidden states per trajectory, thus quantifying the level of uncertainty in these estimates. Bernstein et al. (31) use a maximal likelihood approach inferring a point estimate, using a particle filter within an expectation-maximization algorithm. Using a Bayesian analysis was essential in our study to show that significant intertrajectory heterogeneity was present. Such information was key to identifying the AuNP signatures and separating environment effects from movement characteristics specific to individual AuNPs. In absence of this information, we would only be able to pool the trajectory estimates to define the population statistics. Furthermore, Bernstein et al. (31) infer multiple, static confinement zones separately, as opposed to our single, moving confinement-zone solution. The merits of these distinct approaches to the problem of multiple confinement zones in different experimental applications is a rich area for future work. Finally, we comment that we have implemented a “hop diffusion” model for multiple adjacent confinement zones using reversible-jump MCMC (39), which may be more appropriate for some systems.

As with all analysis, model accuracy needs to be balanced against computational complexity. In practice, this balance is also affected by the data; data at higher spatial and temporal resolution allows more detailed models to be inferred because subtle model differences can then be distinguished. We demonstrated our algorithm on high-resolution iSCAT data, with static localization error 2.7 nm, mean particle diffusivity 1.15 ± 0.106 μm2s−1, and exposure time equal to s (4). By subsampling the data at rate 10 (Figs. S1 and S2; Note S2 of the Supporting Materials and Methods), we mitigate any effects due to dynamic localization error. Under these conditions, confinement models and switching-diffusivity models can be distinguished. For other data types with lower localization accuracy and shorter trajectories, such as fluorescent probes, explicit consideration of measurement noise may be necessary. This is suggested by our simulations, in which we observed a bias in confinement strength κ as the measurement error increased beyond that of iSCAT data (Fig. S7). Our MCMC algorithm could be extended using an HMM with an unknown (hidden) particle position to model measurement noise as in (30); static and dynamic errors could also be separated (21, 40).

The observation that individual AuNP/CTxB/GM1 complexes have a specific spatial signature means that distinguishing the effects of the tag from other factors, such as the cell membrane environment, is difficult. It follows that homogenous tags should improve characterization of the membrane environment. Because the variability in the tag signature presumably arises from the random placing of CTxB molecules on the AuNP surface, using particles with a structured surface is predicted to reduce or potentially eliminate this problem. Virions are ideal, given their highly geometric 3D structure. Interferometric label-free tracking of virions has been demonstrated at 3 s temporal resolution (41); thus, achieving the high spatial and temporal resolution of recent iSCAT (such as in the data set explored here) with viral particle tags is a distinct possibility. The resolution of the tag’s signature and the length of events will then determine the resolution of that trajectory and thus the length scale to which SPT can discriminate different types of particle movement. Whether this can be taken below the size of the tag, reminiscent of super resolution, remains to be ascertained.

Conclusions

We use an HMM-based analysis to partition SPT trajectories into periods of free diffusion and confinement. Our algorithm infers the switching times between these two states, the diffusion coefficient D, and the characteristics of the following confinement events: the HPW strength κ, the position of the HPW center C, and the center diffusion coefficient . We demonstrate the utility of the method on simulated and experimental data; on simulated data, confinement zones were accurately detected, and HPW centers accurately tracked while experimental trajectories were partitioned with high confidence. The model could potentially detect various biological phenomena such as lipid microdomains (or “rafts”), receptor clustering, and hop diffusion.

Author Contributions

Analyzed the data, P.J.S. and N.J.B.; Wrote the article, P.J.S. and N.J.B.; Developed the software, P.J.S.