Bayesian machine learning analysis of singlemolecule fluorescence colocalization images
Abstract
Multiwavelength singlemolecule fluorescence colocalization (CoSMoS) methods allow elucidation of complex biochemical reaction mechanisms. However, analysis of CoSMoS data is intrinsically challenging because of low image signaltonoise ratios, nonspecific surface binding of the fluorescent molecules, and analysis methods that require subjective inputs to achieve accurate results. Here, we use Bayesian probabilistic programming to implement Tapqir, an unsupervised machine learning method that incorporates a holistic, physicsbased causal model of CoSMoS data. This method accounts for uncertainties in image analysis due to photon and camera noise, optical nonuniformities, nonspecific binding, and spot detection. Rather than merely producing a binary ‘spot/no spot’ classification of unspecified reliability, Tapqir objectively assigns spot classification probabilities that allow accurate downstream analysis of molecular dynamics, thermodynamics, and kinetics. We both quantitatively validate Tapqir performance against simulated CoSMoS image data with known properties and also demonstrate that it implements fully objective, automated analysis of experimentderived data sets with a wide range of signal, noise, and nonspecific binding characteristics.
Editor's evaluation
Using a Bayesian machine learning approach, the authors of this paper have developed a tool for the analysis of singlemolecule fluorescence colocalization microscopy images. The authors develop the algorithm, generate an associated software program, and then benchmark the algorithm and software using both simulated and experimental data. The results provide an important, validated tool for use by the singlemolecule fluorescence microscopy community.
https://doi.org/10.7554/eLife.73860.sa0Introduction
A central concern of modern biology is understanding at the molecular level the chemical and physical mechanisms by which protein and nucleic acid macromolecules perform essential cellular functions. The operation of many such macromolecules requires that they work not as isolated molecules in solution but as components of dynamic molecular complexes that selfassemble and change structure and composition as they function. For more than two decades, scientists have successfully explored the molecular mechanisms of many such complex and dynamic systems using multiwavelength single molecule fluorescence methods such as smFRET (singlemolecule fluorescence resonance energy transfer) (Roy et al., 2008) and multiwavelength singlemolecule colocalization methods (CoSMoS, colocalization single molecule spectroscopy) (Larson et al., 2014; van Oijen, 2011; Friedman and Gelles, 2012).
CoSMoS is a technique to measure the kinetics of dynamic interactions between individual molecules. The CoSMoS method has been used for elucidating the mechanisms of complex biochemical processes in vitro. Examples include cell cycle regulation (Lu et al., 2015b), ubiquitination and proteasomemediated protein degradation (Lu et al., 2015a), DNA replication (Geertsema et al., 2014; Ticau et al., 2015), transcription (Zhang et al., 2012; Friedman and Gelles, 2012; Friedman et al., 2013), microRNA regulation (Salomon et al., 2015), premRNA splicing (Shcherbakova et al., 2013; Krishnan et al., 2013; Warnasooriya and Rueda, 2014), ribosome assembly (Kim et al., 2014), translation (Wang et al., 2015; Tsai et al., 2014; O’Leary et al., 2013), signal recognition particlenascent protein interaction (Noriega et al., 2014), and cytoskeletal regulation (Smith et al., 2013; Breitsprecher et al., 2012).
Figure 1A illustrates an example CoSMoS experiment to measure the interaction kinetics of RNA polymerase II molecules with DNA. In the experiment (Rosen et al., 2020), we first measured the locations of individual DNA molecules (the ‘targets’) tethered to the surface of an observation chamber at low density. Next, a cell extract solution containing fluorescent RNA polymerase II molecules (the ‘binders’) was added to the solution over the surface and the chamber surface was imaged by total internal reflection fluorescence (TIRF) microscopy. When the binder molecules are freely diffusing in solution, they are not visible in TIRF. In contrast, when bound to a target, a single binder molecule is detected as a discrete fluorescent spot colocalized with the target position (Friedman et al., 2006; Friedman and Gelles, 2015).
Effective data analysis is a major challenge in the use of the CoSMoS technique. The basic goal is to acquire information at each time point about whether a binder molecule fluorescence spot is observed at the image position of a target molecule (e.g., whether a colocalized greendyelabeled RNA polymerase II is observed at the surface location of a bluedyelabeled DNA spot in Figure 1B). Although CoSMoS images are conceptually simple – they consist only of diffractionlimited fluorescent spots collected in several wavelength channels – efficient analysis of the images is inherently challenging. The number of photons emitted by a single fluorophore is limited by fluorophore photobleaching. Consequently, it is desirable to work at the lowest feasible excitation power in order to maximize the duration of experimental recordings and to efficiently capture relevant reaction events. Achieving higher time resolution divides the number of emitted photons between a larger number of images, so that photon shot noise ordinarily dominates the data statistics. Furthermore, the required concentrations of binder molecules can sometimes create significant background noise (Peng et al., 2018; van Oijen, 2011), even with zeromode waveguide instruments (Chen et al., 2014). These technical difficulties frequently result in CoSMoS images that have low signaltonoise ratios (SNR), making discrimination of colocalized fluorescence spots from noise a significant challenge. In addition, there are usually nonspecific interactions of the binder molecule with the chamber surface, and these artefacts can give rise to both false positive and false negative spot detection (Friedman and Gelles, 2015). Together, these defects in analyzing spot colocalization interfere with the interpretation of CoSMoS data to measure reaction thermodynamics and kinetics and to infer molecular mechanisms.
Most CoSMoS spot detection methods are based on integrating the binder fluorescence intensity by summing the pixel values in small regions of the image centered on the location of individual target molecules, and then using crossings of an intensity threshold to score binder molecule arrival and departure, e.g., (Friedman and Gelles, 2012; Shcherbakova et al., 2013). However, integration discards data about the spatial distribution of intensity that can (and should) be used to distinguish authentic ontarget spots from artefacts caused by noise or offtarget binding. More recently, improved methods (Friedman and Gelles, 2015; Smith et al., 2019) were developed that directly analyze TIRF images, using the spatial distribution of binder fluorescence intensity around the target molecule location. All these methods, whether image or integrated intensitybased, make a binary decision about the presence or absence of a binder spot at the target location. Treating all such binary decisions as equal neglects differences in the confidence of each spot detection decision caused by variations in noise, signal intensity, and nonspecific binding. Failure to account for spot confidence decreases the reliability of downstream thermodynamic and kinetic analysis.
In this paper, we describe a qualitatively different Bayesian machine learning method for analysis of CoSMoS data implemented in a computer program, Tapqir (Kazakh: clever, inventive; pronunciation: tapkeer). Tapqir analyzes twodimensional image data, not integrated intensities. Unlike prior methods, our approach is based on an explicit, global causal model for CoSMoS image formation and uses variational Bayesian inference (KinzThompson et al., 2021; Gelman et al., 2013) to determine the values of model parameters and their associated uncertainties. This model, which we call ‘cosmos’, implements timeindependent analysis of singlechannel (i.e., onebinder) data sets. The cosmos model is physicsinformed and includes realistic shot noise in fluorescent spots and background, camera noise, the size and shape of spots, and the presence of both targetspecific and nonspecific binder molecules in the images. Most importantly, instead of yielding a binary spot/nospot determination, the algorithm calculates the probability of a targetspecific spot being present at each time point and target location. The calculated probability can then be used in subsequent analyses of the molecular thermodynamics and kinetics. Unlike alternative approaches, Tapqir and cosmos do not require subjective threshold settings so they can be used effectively and accurately by nonexpert analysts. The program is implemented in the Pythonbased probabilistic programming language Pyro (Bingham et al., 2019), which enables efficient use of graphics processing unit (GPU)based hardware for rapid parallel processing of data and facilitates future modifications to the model.
Results
Data analysis pipeline
The initial steps in CoSMoS data analysis involve preprocessing the data set (Figure 1B) to map the spatial relationship between target and binder images, correct for microscope drift (if any) and list the locations of target molecules. Software packages that perform these preprocessing steps are widely available (e.g., Friedman and Gelles, 2015; Smith et al., 2019).
The input into Tapqir consists of the time sequence of images (Figure 1B, right). For colocalization analysis, it is sufficient to consider the image area local to the target molecule. This analyzed area of interest (AOI) needs to be several times the diameter of a diffractionlimited spot to include both the spot and the surrounding background (Figure 1C).
In addition to AOIs centered at target molecules, it is useful to also select negative control AOIs from randomly selected sites at which no target molecule is present (Figure 1B and D). In Tapqir, such offtarget control data is analyzed jointly with ontarget data and serves to estimate the background level of targetnonspecific binding.
Once provided with the preprocessing data and image sequence, Tapqir computes for each frame of each AOI the probability, $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$, that a targetspecific fluorescence spot is present. The $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ values that are output can then be used to extract information about the kinetics and thermodynamics of the targetbinder interaction.
Bayesian image classification analysis
Tapqir calculates $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ values using an objective image classification method built on a rigorous Bayesian statistical approach to the CoSMoS image analysis problem. The Bayesian approach has three components. First, we define a probabilistic model of the CoSMoS images. The probabilistic model, cosmos, is a mathematical formalism that describes the AOI images in terms of a set of parameter values. The model is probabilistic in that each parameter is specified to have a probability distribution that defines the likelihood that it can take on particular values. Model parameters describe physically realistic image features such as the characteristic fluorescence spot width. Second, we specify prior distributions for the parameters of the model. These priors embed preexisting knowledge about the CoSMoS experiment, such as the fact that targetspecific spots will be close to the target molecule locations. Third, we infer the values of the model parameters, including $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$, using Bayes’ rule (Bishop, 2006; KinzThompson et al., 2021). The cosmos model is ‘timeindependent’, meaning that we ignore the time dimension of the recording – the order of the images does not affect the results.
Probabilistic image model and parameters
A single AOI image from a CoSMoS data set is a matrix of noisy pixel intensity values. In each image, multiple binder molecule fluorescence spots can be present. Figure 2A shows an example image where two spots are present; one spot is located near the target molecule at the center of the image and another is offtarget.
The probabilistic model mathematically generates images $D$ as follows. We construct a noisefree AOI image ${\mu}^{I}$ as a constant average background intensity $b$ summed with fluorescence spots modeled as 2D Gaussians ${\mu}^{S}$, which accurately approximate the microscope point spread function (Zhang et al., 2007; Figure 2B). Each 2D Gaussian is described by parameters integrated intensity $h$, width $w$, and position ($x$, $y$). We define $K$ as the maximum number of spots that can be present in a single AOI image. For the data we typically encounter, $K$ = 2 is sufficient. Since the spots may be present or not in a particular image, we define the $K$ = 2 binary indicators $m}_{\mathsf{s}\mathsf{p}\mathsf{o}\mathsf{t}(1)$ and $m}_{\mathsf{s}\mathsf{p}\mathsf{o}\mathsf{t}(2)$. Each indicator can take a value of either 0 denoting spot absence or 1 denoting spot presence.
The resulting mixture model has four possible combinations for $m}_{\mathsf{s}\mathsf{p}\mathsf{o}\mathsf{t}(1)$ and $m}_{\mathsf{s}\mathsf{p}\mathsf{o}\mathsf{t}(2)$: (1) a nospot image that contains only background (Figure 2B, top left), (2) a singlespot image that contains the first binder molecule spot superimposed on background (Figure 2B, bottom left), (3) a singlespot image that contains the second binder molecule spot superimposed on background (Figure 2B, top right), and (4) a twospot image that contains both binder molecule spots superimposed on background (Figure 2B, bottom right).
Among the spots that are present in an AOI image, by assumption at most only one can be targetspecific. We use a state parameter $z$ to indicate targetspecific spot absence ($z$ = 0) or presence ($z$ = 1) in an AOI image. We also introduce an index parameter $\theta$ that identifies which of the spots is the targetspecific spot when it is present ($z$ = 1) (e.g., Figure 2C, middle and right have $\theta$ = 1 and $\theta$ = 2, respectively) and equals zero when it is absent ($z$ = 0) (e.g., Figure 2C, left). Since the offtarget control AOIs by definition contain only nonspecific binding, $z$ = 0 and $\theta$ = 0 for all offtarget AOIs.
Finally, to construct realistic noisy AOI images $D$ from the noisefree images ${\mu}^{I}$, the model adds intensitydependent noise to each pixel. For cameras that use chargecoupled device (CCD) or electronmultiplier CCD (EMCCD) sensors, each measured pixel intensity in a singlemolecule fluorescence image has a noise contribution from photon counting (shot noise) and can also contain additional noise arising from electronic amplification (van Vliet et al., 1998). The result is a characteristic linear relationship between the noise variance and mean intensity with slope defining the gain $g$. This relationship is used to compute the random pixel noise values (see Materials and methods).
The resulting probabilistic image model can be interpreted as a generative process that produces the observed image data $D$. A graphical representation of the probabilistic relationships in the model is shown in Figure 2D. A complete description of the model is given in Materials and methods and Figure 2—figure supplement 1.
Parameter prior distributions
Specifying prior probability distributions for model parameters is essential for Bayesian analysis and allows us to incorporate preexisting knowledge about the experimental design. For most model parameters, there is no strong prior information so we use uninformative prior distributions (see Materials and methods). However, we have strong expectations for the positions of specific and nonspecific binder molecules that can be expressed as prior distributions and used effectively to discriminate between the two. Nonspecific binding can occur anywhere on the surface with equal probability and thus has a uniform prior distribution across the AOI image. Targetspecific binding, on the other hand, is colocalized with the target molecule and thus has a prior distribution peaked at the AOI center (Figure 2—figure supplement 2). The width of this peak, proximity parameter ${\sigma}^{xy}$, depends on multiple features of the experiment such as the spot localization accuracy and the mapping accuracy between target and binder imaging channels. Prior distributions for parameters $\theta$ and $m$ are defined in terms of the average number of targetspecific and target nonspecific spots per AOI image, $\pi$ and $\lambda$, respectively. To facilitate convenient use of the algorithm, it is not necessary to prespecify values of ${\sigma}^{xy}$, $\pi$, and $\lambda$. Instead, values of these parameters appropriate to a given data set are calculated automatically using a hierarchical Bayesian analysis (see Materials and methods; for hierarchical modeling see Chapter 5 of Gelman et al., 2013).
Bayesian inference and implementation
Tapqir calculates posterior distributions of model parameters conditioned on the observed data by using Bayes’ theorem. In particular, Tapqir approximates posterior distributions using a variational inference approach implemented in Pyro (Bingham et al., 2019). Complete details of the implementation are given in Materials and methods.
Tapqir analysis
In initial tests, we used Tapqir to analyze simulated CoSMoS image data with a comparatively high SNR of 3.76 as well as data from the experiment shown in Figure 1B–D, which has a lower SNR of 1.61. The simulated data were generated using the same cosmos model (Figure 2D) that was used for analysis. Tapqir correctly detects fluorescent spots in both simulated and experimental images (compare ‘AOI images’ and ‘Spotdetection’ rows in Figure 3). The program precisely calculates the position ($x$, $y$), intensity ($h$), and width ($w$) for each spot and also determines the background intensity ($b$) for each image without requiring a separate analysis. These parameters confirm the desired behavior of the model and could be used in further calculations. However, the most important output of the analysis is assessment of the presence of targetspecific binding. For each AOI image, we calculate $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})\equiv p(z=1)$ (Figure 3, green), the probability that any targetspecific spot is present. Spots determined as likely targetspecific ($p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ > 0.5) are represented as filled circles in the spot detection row of Figure 3. For a particular spot to have high $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$, it must have a high spot probability and be colocalized with the target molecule at the center of the AOI (Figure 3—figure supplement 1).
Tapqir robustly fits experimental data sets with different characteristics
Next, we evaluated how well the model fits data sets encompassing a range of characteristics found in typical CoSMoS experiments. We analyzed four experimental data sets with varying SNR, frequency of targetspecific spots, and frequencies of nonspecific spots (Table 1). We then sampled AOI images from the posterior distributions of parameters (a method known as posterior predictive checking Gelman et al., 2013). These posterior predictive simulations accurately reproduce the experimental AOI appearances, recapitulating the noise characteristics and the numbers, intensities, shapes, and locations of spots (Figure 3—figure supplement 2, images). The distributions of pixel intensities across the AOI are also closely reproduced (Figure 3—figure supplement 2, histograms) confirming that the noise model is accurate. Taken together, these results confirm that the model is rich enough to accurately capture the full range of image characteristics from CoSMoS data sets taken over different experimental conditions. Importantly, all the results on different experimental data sets were obtained using the same model (Figure 2D) and the same priors (Materials and methods). No tuning of the algorithm or prior measurement of datasetspecific properties was needed to achieve good fits for all data sets.
Tapqir accuracy on simulated data with known global parameter values
Next, we evaluated Tapqir’s ability to reliably infer the values of global model parameters. To accomplish this, we generated simulated data sets using a wide range of randomized parameter values and then fit the simulated data to the model (Supplementary file 2). Fit results show that global model parameters (i.e., average specific spot probability $\pi$, nonspecific binding density $\lambda$, proximity ${\sigma}^{xy}$, and gain $g$; see Figure 2D) are close to the simulated values (Figure 3—figure supplement 3 and Supplementary file 2). This suggests that CoSMoS data contains enough information to reliably infer global model parameters and that the model is not obviously overparameterized.
Tapqir classification accuracy
Having tested the basic function of the algorithm, we next turned to the key question of how accurately Tapqir can detect targetspecific spots in data sets of increasing difficulty.
We first examined the accuracy of targetspecific spot detection in simulated data sets with decreasing SNR (Supplementary file 3). By eye, spots can be readily discerned at SNR >1 but cannot be clearly seen at SNR <1 (Figure 4A). Tapqir gives similar or better performance: if an image contains a targetspecific spot, Tapqir correctly assigns it a targetspecific spot probability $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ that is on average close to one as long as SNR is adequate (i.e., SNR >1) (Figure 4B). In contrast, mean $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ sharply decreases at SNR <1, consistent with the subjective impression that no spot is recognized under those conditions. In particular, images that contain a targetspecific spot are almost always assigned a high $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ for high SNR data and almost always assigned low $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ for low SNR data (Figure 4C, green). At marginal SNR ≈ 1, these images are assigned a broad distribution of $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ values, accurately reflecting the uncertainty in classifying such data. Just as importantly, images with no targetspecific spot are almost always assigned $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ < 0.5, correctly reflecting the absence of the spot (Figure 4C, gray).
Ideally, we want to correctly identify targetspecific binding when it occurs but also to avoid incorrectly identifying targetspecific binding when it does not occur. To quantify Tapqir’s classification accuracy, we next examined binary image classification statistics. Binary classification predictions were obtained by thresholding $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ at 0.5. We then calculated two complementary statistics: recall and precision (Fawcett, 2006; Figure 4D; see Materials and methods). Recall is defined as the fraction of true targetspecific spots that are correctly predicted. Recall is high at high SNR and decreases at lower SNR. Recall is a binary analog of the mean $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ for the subset of images containing targetspecific spots; as expected the two quantities have similar dependencies on SNR (compare Figure 4B and D, black). Precision is the fraction of predicted targetspecific spots that are correctly predicted. Precision is near one at all SNR values tested (Figure 4D, red); this shows that the algorithm rarely misclassifies an image as containing a targetspecific spot when none is present.
In order to quantify the effects of both correctly and incorrectly classified images in a single statistic, we used the binary classification predictions to calculate the Matthews Correlation Coefficient (MCC) (Matthews, 1975; see Materials and methods). The MCC is equivalent to the Pearson correlation coefficient between the predicted and true classifications, giving 1 for a perfect match, 0 for a random match, and –1 for complete disagreement. The MCC results (Figure 4D, blue) suggest that the overall performance of Tapqir is excellent at SNR ≥ 1: the program rarely misses targetspecific spots that are in reality present and rarely falsely reports a targetspecific spot when none is present.
The analyses of Figure 4B–D examined Tapqir performance on data in which the rate of targetnonspecific binding is moderate ($\lambda$ = 0.15 nonspecific spots per AOI image on average). We next examined the effects of increasing the nonspecific rate. In particular, we used simulated data (Supplementary file 1) with high SNR = 3.76 to test the classification accuracy of Tapqir at different nonspecific binding densities up to $\lambda$ = 1, a value considerably higher than typical of usable experimental data (the experimental data sets in Table 1 have $\lambda$ ranging from 0.04 to 0.30). In analysis of these data sets, a few images with targetspecific spots are misclassified as not having a specific spot ($p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ near zero) or as being ambiguous ($p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ near 0.5) (Figure 4F, green bars), and a few images with targetnonspecific spots are misclassified as having specific spot ($p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ near or above 0.5) (Figure 4F, gray bars), but these misclassifications only occurred at the unrealistically high $\lambda$ value. Even in the simulation with this highest $\lambda$ value, Tapqir accurately identified targetspecific spots (Figure 4E and F) and returned excellent binary classification statistics (Figure 4G).
A weakness of some existing imagebased CoSMoS spot discrimination methods is that targetnonspecific binding adjacent to a targetspecific spot can interfere with correctly identifying the latter as targetspecific. The very high recall values obtained at $\lambda$ = 1 (Figure 4G) confirm that there are few such misidentifications by Tapqir even at high nonspecific binding densities. This good performance is likely facilitated by the feature of the Tapqir model that explicitly includes the possibility that both a specifically and a nonspecifically bound spot may occur simultaneously in the same AOI. Consistent with this interpretation, we see effective detection of the specific and nonspecific spots even in example AOIs in which the two spots are so closely spaced that they are not completely resolved (Figure 4H). In contrast, tests of existing CoSMoS image classification methods show that images with targetnonspecific spots are prone to misclassification. As discussed previously (Friedman and Gelles, 2015), methods based on thresholding of integrated AOI intensities are prone to incorrectly classify targetnonspecific spots as targetspecific. Conversely, an existing ‘spotpicker’ method based on empirical binary classification of 2D AOI images (Friedman and Gelles, 2015) is much more likely than Tapqir to fail to detect target specific spots when there is a nearby nonspecific spot (Figure 4—figure supplement 1). This contributes to the superior overall performance we see for Tapqir vs. spotpicker on the $\lambda$ = 1 data set (recall 0.993 vs 0.919; precision 0.943 vs 0.873; MCC 0.961 vs 0.874).
To further evaluate whether Tapqir is prone to misidentifying targetnonspecific spots as specific, we simulated data sets with no targetspecific binding at both low and high nonspecific binding densities (Supplementary file 4). Analysis of such data (Figure 4I) shows that no targetspecific binding (i.e., $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ > 0.6) was detected even under the highest nonspecific binding density, demonstrating that Tapqir is robust to falsepositive targetspecific spot detection even under these extreme conditions.
Since targetnonspecific spots are built into the cosmos model, there is no need to choose excessively small AOIs in an attempt to exclude nonspecific spots from analysis. We found that reducing AOI size (from 14 × 14 to 6 x 6 pixels) did not appreciably affect analysis accuracy on simulated data, when the width ($w$) of the spots was equal to 1.4 pixels (Table 2). In analysis of experimental data, smaller AOI sizes caused occasional changes in calculated $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ values reflecting apparent missed detection of a few spots (Figure 3—figure supplement 4). Out of caution, we therefore used 14 × 14 pixel AOIs routinely, even though the larger AOIs somewhat reduced computation speed (Table 2 and Figure 3—figure supplement 4).
Kinetic and thermodynamic analysis of molecular interactions
The most widespread application of CoSMoS experiments is to measure rate and equilibrium constants for the binding interaction of the target and binder molecules being studied. We next tested whether these constants can be accurately determined using Tapqircalculated posterior predictions.
We first simulated CoSMoS data sets (Supplementary file 5) that reproduced the behavior of a onestep association/dissociation reaction mechanism (Figure 5A and B, blue). Simulated data were analyzed with Tapqir yielding $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ values for each frame (e.g., Figure 5B, green). We wanted to estimate rate constants using the full information contained in the $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ probabilities, so we did not threshold $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ for this analysis. Instead, from each singleAOI $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ time record we constructed a family of binary time records (Figure 5B, black) by Monte Carlo sampling according to the $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ time series. Each family member has welldefined targetspecific binderpresent and binderabsent intervals $\mathrm{\Delta}{t}_{\mathsf{o}\mathsf{n}}$ and $\mathrm{\Delta}{t}_{\mathsf{o}\mathsf{f}\mathsf{f}}$, respectively. Each of these time records was then analyzed with a twostate hidden Markov model (HMM) (see Materials and methods), producing a distribution of inferred rate constants from which we calculated mean values and their uncertainties (Figure 5C and D). Comparison of the simulated and inferred values shows that both $k}_{\mathsf{o}\mathsf{n}$ and $k}_{\mathsf{o}\mathsf{f}\mathsf{f}$ rate constants are accurate within 30% at nonspecific binding densities typical of experimental data ($\lambda$ ≤ 0.5). At higher nonspecific binding densities, rare interruptions caused by falsepositive and falsenegative spot detection shorten $\mathrm{\Delta}{t}_{\mathsf{o}\mathsf{n}}$ and $\mathrm{\Delta}{t}_{\mathsf{o}\mathsf{f}\mathsf{f}}$ distributions, leading to moderate systematic overestimation of the association and dissociation rate constants.
From the same simulated data, we calculated the equilibrium constant $K}_{\mathsf{e}\mathsf{q}$ and its uncertainty. This calculation does not require a timedependent model and can be obtained directly from the posterior distribution of the average specificbinding probability $\pi$. The estimated equilibrium constants are highly accurate even at excessively high values of $\lambda$ (Figure 5E). The high accuracy results from the fact that equilibrium constant measurements are in general much less affected than kinetic measurements by occasional false positives and false negatives in spot detection.
The forgoing analysis shows that Tapqir can accurately recover kinetic and thermodynamic constants from simulated CoSMoS data. However, experimental CoSMoS data sets can be more diverse. In addition to having different SNR and nonspecific binding frequency values, they also may have nonidealities in spot shape (caused by optical aberrations) and in noise (caused by molecular diffusion in and out of the TIRF evanescent field). In order to see if Tapqir analysis is robust to these and other properties of real experimental data, we analyzed several CoSMoS data sets taken from different experimental projects. Analysis of each data set took a few hours of computation time on a GPUequipped desktop computer or cloud computing service (Table 1). We first visualized the results as probabilistic rastergrams (Figure 6A, Figure 6—figure supplement 1A, Figure 6—figure supplement 2A, and Figure 6—figure supplement 3A), in which each horizontal line represents the time record from a single AOI. Unlike the binary spot/nospot rastergrams in previous studies (e.g., Friedman et al., 2013; Rosen et al., 2020) we plotted the Tapqircalculated spot probability $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ using a color scale. This representation allows a more nuanced understanding of the data. For example, Figure 6A reveals that while the longduration spot detection events typically are assigned a high probability (yellow), some of the shortest duration events have an intermediate $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ (green) indicating that the assignment of these as targetspecific is uncertain.
To demonstrate the utility of Tapqir for kinetic analysis of real experimental data, we measured binder association rate constants in previously published experimental data sets (Table 1). We employed our previous strategy (Friedman and Gelles, 2012; Friedman and Gelles, 2015) of analyzing the duration of the binderabsent intervals that preceded the first binding event. Such timetofirst binding analysis improves the accuracy of association rate constant estimates relative to those obtained by analyzing all $\mathrm{\Delta}{t}_{\mathsf{o}\mathsf{f}\mathsf{f}}$ values by minimizing the effects of target molecules occupied by photobleached binders, dye blinking and false negative dropouts that occur within a continuous binder dwell interval. To perform a timetofirstbinding analysis using Tapqir, we used the posterior sampling method (as in Figure 5B, black records) to determine the initial $\mathrm{\Delta}{t}_{\mathsf{o}\mathsf{f}\mathsf{f}}$ in each AOI record. These data were fit to a kinetic model (Friedman and Gelles, 2012; Friedman and Gelles, 2015) in which only a fraction of target molecules $A}_{\mathsf{f}$ were binding competent and which includes both exponential targetspecific association with rate constant $k}_{\mathsf{a}$, as well as exponential nonspecific association with rate constant $k}_{\mathsf{n}\mathsf{s}$ (Figure 6B, Figure 6—figure supplement 1B, Figure 6—figure supplement 2B, and Figure 6—figure supplement 3B). The Tapqirderived fits showed excellent agreement with the kinetic model.
To further assess the utility of the Tapqir method, we used experimental data sets and compared the Tapqir association kinetics results with those from the previously published empirical binary ‘spotpicker’ method (Friedman and Gelles, 2015; Figure 6C, Figure 6—figure supplement 1C, Figure 6—figure supplement 2C, and Figure 6—figure supplement 3C). The values of the association rate constant $k}_{\mathsf{a}$ obtained using these two methods are in good agreement with each other (Figure 6D, Figure 6—figure supplement 1D, Figure 6—figure supplement 2D, and Figure 6—figure supplement 3D). We emphasize that while Tapqir is fully objective, achieving these results with the spotpicker method required optimization by subjective adjustment of spot detection thresholds. We noted some differences between the two methods in the nonspecific association rate constants $k}_{\mathsf{n}\mathsf{s}$. Differences are expected because these parameters are defined differently in the different nonspecific binding models used in Tapqir and spotpicker (see Materials and methods).
Discussion
A broad range of physical processes contribute to the formation of CoSMoS images. These include camera and photon noise, targetspecific and nonspecific binding, and time and positiondependent variability in fluorophore imaging and image background. Unlike prior CoSMoS analysis methods, Tapqir considers these aspects of imaging in a single, holistic model. This cosmos model explicitly includes the uncertainties due to photon noise, camera gain, and spatial variability in intensity offset. The model also includes the possibility of multiple binder molecule fluorescence spots being present in the vicinity of the target, including both targetspecific binding and targetnonspecific interactions of binder molecules with the coverslip surface. This explicit modeling of targetnonspecific spots makes it possible to include offtarget control data as a part of the experimental data set. Similarly, all AOIs and frames in the data set are simultaneously fit to the global model in a way that allows for realistic frametoframe and AOItoAOI variability in image formation caused by variations in laser intensity, fluctuations in background, and other nonidealities. The global analysis based on a single, unified model enables the final results (e.g., kinetic and thermodynamic parameters) to be estimated in a way that is cognizant of the known sources of uncertainty in the data.
Previous approaches to CoSMoS data analysis, including our spotpicker method (Friedman and Gelles, 2015), did not employ a holistic modeling approach and instead relied on a multistep process that includes a separate binary classification step. These prior methods require subjective setting of classification thresholds. Because they are not fully objective, such methods cannot reliably account for uncertainties in spot classification, which compromises error estimates in the analysis pipeline downstream of spot classification. One recent approach (Smith et al., 2019; Smith et al., 2015), which like spotpicker and Tapqir analyzes 2D images instead of integrated intensities, used a Bayesian kinetic analysis but a frequentist hypothesis test (a generalized likelihood ratio test) for spot detection. The frequentist method lacks a key advantage of Tapqir’s modelbased Bayesian approach that here enables prediction of targetspecific spot presence probabilities $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ for each image, rather than a binary ‘spot/no spot’ classification. In general, previous approaches in essence assume that spot classifications are correct, and thus the uncertainties in the derived molecular properties (e.g., equilibrium constants) are systematically underestimated because the errors in spot classification, which can be large, are not accounted for. By performing a probabilistic spot classification, Tapqir enables reliable inference of molecular properties, such as thermodynamic and kinetic parameters, and allows statistically welljustified estimation of parameter uncertainties. This more inclusive error estimation likely accounts for the generally larger kinetic parameter error bars obtained from Tapqir compared to those from the existing spotpicker analysis method (Figure 6, Figure 6—figure supplement 1, Figure 6—figure supplement 2, and Figure 6—figure supplement 3). Even though existing analysis methods take advantage of subjective tuning by a human analyst, our comparisons show that Tapqir performs at least comparably to (Figure 6, Figure 6—figure supplement 1, Figure 6—figure supplement 2, and Figure 6—figure supplement 3) and under some conditions much better than (Figure 4—figure supplement 1) the existing spotpicker method.
The Tapqir cosmos model includes parameters of mechanistic interest, such as the average probability of targetspecific binding, as well as ‘nuisance’ parameters that are not of primary interest but nevertheless essential for image modeling. In previous imagebased methods for CoSMoS analysis (e.g., Friedman and Gelles, 2015; Smith et al., 2019), nuisance parameters were either measured in separate experiments (e.g., gain was determined from calibration data), set heuristically (e.g., a subjective choice of userset thresholds for spot intensity and proximity in colocalization detection), or determined at a separate analysis step (e.g., rate of nonspecific binding). In contrast, Tapqir directly learns parameters from the full set of experimental data, thus eliminating the need for additional experiments, subjective adjustment of tuning parameters, and postprocessing steps.
Bayesian analysis has been used previously to analyze data from singlemolecule microscopy experiments (e.g., KinzThompson et al., 2021 and references cited therein). A key feature of Bayesian analysis is that the extent of prior knowledge of all model parameters is explicitly incorporated. Where appropriate, cosmos uses relatively uninformative priors that only weakly specify information about the value of the corresponding parameters. In these cases, cosmos mostly infers parameter values from the data. In contrast, some priors are more informative. For example, binder molecule spots near the target molecule are more likely to be targetspecific rather than targetnonspecific, so we use this known feature of the experiment by encoding the likely position of targetspecific binding as a databased prior. This tactic effectively enables probabilistic classification of spots as either targetspecific or targetnonspecific, which would be difficult using other inference methodologies, while still accommodating data sets with different accuracies of mapping between binder and target channels.
Tapqir is implemented in Pyro, a Pythonbased probabilistic programming language (PPL) (Bingham et al., 2019). Probabilistic programming is a relatively new paradigm in which probabilistic models are expressed in a highlevel language that allows easy formulation, modification, and automated inference (van de Meent et al., 2018). In this work we focused on developing an image model for colocalization detection in a relatively simple bindertarget singlemolecule experiment. However, Tapqir can be used with more complex models. For example, the cosmos model could be naturally extended to multistate and multicolor analysis. Furthermore, with the development of more efficient sequential hidden Markov modeling algorithms (Särkkä and GarcíaFernández, 2019; Obermeyer et al., 2019b) Tapqir can potentially be extended to directly incorporate kinetic processes, allowing direct inference of kinetic mechanisms and rate constants.
Tapqir is free, opensource software. Tapqir is available at https://github.com/gellesbrandeis/tapqir. The results presented here were obtained using release 1.0 of the program (https://github.com/gellesbrandeis/tapqir/releases/tag/v1.0). The Tapqir documentation, which contains tutorials on program use, is at https://tapqir.readthedocs.io/en/stable/. Source data including Figures, Figure supplements, Supplementary files, manuscript text, and the scripts and data used to generate them are available at https://github.com/ordabayevy/tapqiroverleaf.
Materials and methods
Notation
Request a detailed protocolIn the Materials and methods section, we adopt a mathematical notation for multidimensional arrays from the field of machine learning (Chiang et al., 2021). The notation uses named axes and incorporates implicit broadcasting of arrays when their shapes are different.
Extracting image data
Request a detailed protocolRaw input data into Tapqir consists of (1) binder channel images ($D}^{\mathsf{r}\mathsf{a}\mathsf{w}$), each $W\times H$ pixels in size, for each time point (Figure 1B, right) and (2) lists of locations, corrected for microscope drift if necessary (Friedman and Gelles, 2015), of target molecules and of offtarget control locations (Friedman and Gelles, 2015) within the raw images. For simplicity, we use the same notation ($x}^{\mathsf{t}\mathsf{a}\mathsf{r}\mathsf{g}\mathsf{e}\mathsf{t},\mathsf{r}\mathsf{a}\mathsf{w}$, $y}^{\mathsf{t}\mathsf{a}\mathsf{r}\mathsf{g}\mathsf{e}\mathsf{t},\mathsf{r}\mathsf{a}\mathsf{w}$) both for target molecule locations and offtarget control locations. Tapqir extracts a $P\times P$ AOI around each target and offtarget location and returns (1) the extracted data set $D$ consisting of a set of $P\times P$ grayscale images, collected at $N$ ontarget AOI sites and $N}_{\mathsf{c}$ offtarget AOI sites for a range of $F$ frames (Figure 1C and D; Figure 7), and (2) new target (and offtarget) locations ($x}^{\mathsf{t}\mathsf{a}\mathsf{r}\mathsf{g}\mathsf{e}\mathsf{t}$, $y}^{\mathsf{t}\mathsf{a}\mathsf{r}\mathsf{g}\mathsf{e}\mathsf{t}$) adjusted relative to extracted images $D$ where $x}^{\mathsf{t}\mathsf{a}\mathsf{r}\mathsf{g}\mathsf{e}\mathsf{t}$ and $y}^{\mathsf{t}\mathsf{a}\mathsf{r}\mathsf{g}\mathsf{e}\mathsf{t}$ both lie within the $(P/21,P/2)$ central range of the image. For the data presented in this article, we used $P$ = 14. Cartesian pixel indices (i, $j$) are integers but also represent the center point of a pixel on the image plane. While experimental intensity measurements are integers, we treat them as continuous values in our analysis.
The cosmos model
Our intent is to model CoSMoS image data by accounting for the significant physical aspects of image formation, such as photon noise and binding of targetspecific and targetnonspecific molecules to the microscope slide surface. A graphical representation of the Tapqir model for CoSMoS data similar to that in Figure 2D but including probability distributions and other additional detail is shown in Figure 2—figure supplement 1. The corresponding generative model represented as pseudocode is shown in Figure 8. All variables with short descriptions and their domains are listed in Table 3. Below, we describe the model in detail starting with the observed data and the likelihood function and then proceed with model parameters and their prior distributions.
Image likelihood
Request a detailed protocolWe model the image data $D$ as the sum of a photonindependent offset $\delta$ introduced by the camera and the noisy photondependent pixel intensity values $I$:
In our model, each pixel in the photondependent image $I$ has a variance which is equal to the mean intensity ${\mu}^{I}$ of that pixel multiplied by the camera gain $g$, which is the number of camera intensity units per photon. This formulation is appropriate for cameras that use chargecoupled device (CCD) or electronmultiplier CCD (EMCCD) sensors. (The experimental CoSMoS datasets we analyzed (Table 1) were collected with EMCCD cameras.) It accounts for both photon shot noise and additional noise introduced by EMCCD camera amplification (van Vliet et al., 1998) and is expressed using a continuous Gamma distribution:
The Gamma distribution was chosen because we found it to effectively model the image noise, which includes both Poissonian (shot noise) and nonPoissonian contributions. The Gamma distribution used here is parameterized by its mean and standard deviation. The functional forms of the Gamma distribution and all other distributions we use in this work are given in Table 4.
A competing camera technology based on scientific complementary metaloxide semiconductor (sCMOS) sensors produces images that have also successfully been modeled as having a combination of Poissonian and nonPoissonian (Gaussian, in this case) noise sources. However, sCMOS images have noise characteristics that are considerably more complicated than CCD/EMCCD images, because every pixel has its own characteristic intensity offset, Gaussian noise variance, and amplification gain. Additional validation will be required to determine whether the existing cosmos model requires modification or inclusion of additional prior information (e.g., pixelbypixel calibration data as in Huang et al., 2013) to optimize its performance with sCMOS CoSMoS data.
Image model
Request a detailed protocolThe idealized noisefree image ${\mu}^{I}$ is represented as the sum of a background intensity $b$ and the intensities from fluorescence spots modeled as 2D Gaussians ${\mu}^{S}$:
For simplicity, we allow at most $K$ number of spots in each frame of each AOI. (In this article, we always use $K$ equal to 2.) The presence of a given spot in the image is encoded in the binary spot existence parameter $m$, where $m$ = 1 when the corresponding spot is present and $m$ = 0 when it is absent.
The intensities for a 2D Gaussian spot at each pixel coordinate (i, $j$) is given by:
with spot parameters total integrated intensity $h$, width $w$, and center ($x$, $y$) relative to the target (or offtarget control) location ($x}^{\mathsf{t}\mathsf{a}\mathsf{r}\mathsf{g}\mathsf{e}\mathsf{t}$, $y}^{\mathsf{t}\mathsf{a}\mathsf{r}\mathsf{g}\mathsf{e}\mathsf{t}$).
Our primary interest is whether a targetspecific spot is absent or present in a given AOI. We encode this information using a binary state parameter $z$ with 0 and 1 denoting targetspecific spot absence and presence, respectively. To indicate which of the $K$ spots is targetspecific, we use the index parameter $\theta$ which ranges from 0 to $K$. When a targetspecific spot is present ($z$ = 1), $\theta \in \{1,\dots ,K\}$ specifies the index of the targetspecific spot, while $\theta$ = 0 indicates that no targetspecific spot is present ($z$ = 0). For example, $\{{m}_{\mathsf{s}\mathsf{p}\mathsf{o}\mathsf{t}(1)}=1,{m}_{\mathsf{s}\mathsf{p}\mathsf{o}\mathsf{t}(2)}=1,z=1,\theta =2\}$ means that both spots are present and spot 2 is targetspecific. A combination like $\{{m}_{\mathsf{s}\mathsf{p}\mathsf{o}\mathsf{t}(1)}=0,{m}_{\mathsf{s}\mathsf{p}\mathsf{o}\mathsf{t}(2)}=1,z=1,\theta =1\}$ is impossible (i.e., has zero probability) since spot 1 cannot be absent and targetspecific at the same time. For offtarget control data, in which no spots are targetspecific by definition, $z$ and $\theta$ are always set to zero.
Prior distributions
Request a detailed protocolThe prior distributions for the model parameters are summarized in Figure 2—figure supplement 1 and detailed below. Unless otherwise indicated we assume largely uninformative priors (such as the HalfNormal distribution with large mean).
Background intensity $b$ follows a Gamma distribution:
where the mean $\mu}^{b}\in {\mathbb{R}}_{>0}^{\mathsf{A}\mathsf{O}\mathsf{I}[N]$ and standard deviation $\sigma}^{b}\in {\mathbb{R}}_{>0}^{\mathsf{A}\mathsf{O}\mathsf{I}[N]$ of the background intensity describe the irregularity in the background intensity in time and across the field of view of the microscope. Priors for ${\mu}^{b}$ and ${\sigma}^{b}$ are uninformative:
The targetspecific presence parameter $z$ has a Bernoulli prior parameterized by the average targetspecific binding probability $\pi$ for ontarget AOIs and zero probability for control offtarget AOIs:
The prior distribution for the index of the targetspecific spot $\theta$ is conditional on $z$. When no specifically bound spot is present (i.e., $z$ = 0), $\theta$ always equals 0. Since spot indices are arbitrarily assigned, when the targetspecific spot is present (i.e., $z$ = 1) $\theta$ can take any value between 1 and $K$ with equal probability. We represent the prior for $\theta$ as a Categorical distribution of the following form:
The average targetspecific binding probability $\pi$ has an uninformative Jeffreys prior (Gelman et al., 2013) given by a Beta distribution:
The prior distribution for the spot presence indicator $m$ is conditional on $\theta$. When $\theta$ corresponds to spot index $k$, i.e., $\theta =k$, then $m}_{\mathsf{s}\mathsf{p}\mathsf{o}\mathsf{t}(k)$ = 1. When $\theta$ does not correspond to a spot index $k$, that is, $\theta \ne k$, then either spot $k$ is targetnonspecific or a spot corresponding to $k$ does not exist. Consequently, for $\theta \ne k$ we assign $m}_{\mathsf{s}\mathsf{p}\mathsf{o}\mathsf{t}(k)$ to either 0 or 1 with a probability dependent on the nonspecific binding density $\lambda \in {\mathbb{R}}_{>0}$:
The mean nonspecific binding density $\lambda$ is expected to be much less than two nonspecifically bound spots per frame per AOI; therefore, we use an Exponential prior of the form
The prior distribution for the integrated spot intensity $h$ is chosen to fall off at a value much greater than typical spot intensity values
In CoSMoS experiments, the microscope/camera hardware is typically designed to set the width $w$ of fluorescence spots to a typical value in the range of 1–2 pixels (Ober et al., 2015). We use a Uniform prior confined to the range between 0.75 and 2.25 pixels:
Priors for spot position ($x$, $y$) depend on whether the spot represents targetspecific or nonspecific binding. Nonspecific binding to the microscope slide surface can occur anywhere within the image and therefore has a uniform distribution (Figure 2—figure supplement 2, red). Spot centers may fall slightly outside the AOI image yet still affect pixel intensities within the AOI. Therefore the range for ($x$, $y$) is extended one pixel wider than the size of the image, which allows a spot center to fall slightly beyond the AOI boundary.
In contrast to nonspecifically bound molecules, specifically bound molecules are colocalized with the target molecule with a precision that can be better than one pixel and that depends on various factors including the microscope pointspread function and magnification, accuracy of registration between binder and target image channels, and accuracy of drift correction. For targetspecific binding, we use an AffineBeta prior with zero mean position relative to the target molecule location ($x}^{\mathsf{t}\mathsf{a}\mathsf{r}\mathsf{g}\mathsf{e}\mathsf{t}$, $y}^{\mathsf{t}\mathsf{a}\mathsf{r}\mathsf{g}\mathsf{e}\mathsf{t}$), and a ‘proximity’ parameter ${\sigma}^{xy}$ which is the standard deviation of the AffineBeta distribution (Figure 2—figure supplement 2, green). We chose the AffineBeta distribution because it models a continuous parameter defined on a bounded interval.
We give the proximity parameter ${\sigma}^{xy}$ a diffuse prior, an Exponential with a characteristic width of one pixel:
Tests on data simulated with increasing proximity parameter values ${\sigma}^{xy}$ (true) (i.e., with decreasing precision of spatial mapping between the binder and target image channels) confirm that the cosmos model accurately learns ${\sigma}^{xy}$ (fit) from the data (Figure 3—figure supplement 3D; Table 5). This was the case even if we substituted a lessinformative ${\sigma}^{xy}$ prior (Uniform vs. Exponential; Table 5).
The CoSMoS technique is premised on colocalization of the binder spots with the known location of the target molecule. Consequently, for any analysis method, classification accuracy declines when the images in the target and binder channels are less accurately mapped. For the Tapqir cosmos model, low mapping precision has little effect on classification accuracy at typical nonspecific binding densities ($\lambda$ = 0.15; see MCC values in Table 5).
Gain $g$ depends on the settings of the amplifier and electron multiplier (if present) in the camera. It has a positive value and is typically in the range between 5 and 50. We use a HalfNormal prior with a broad distribution encompassing this range:
The prior distribution for the offset signal $\delta$ is empirically measured from the output of camera sensor regions that are masked from incoming photons. Collected data from these pixels are transformed into a density histogram with intensity step size of 1. The resulting histogram typically has a long right hand tail of low density. For computational efficiency, we shorten this tail by binning together pixel intensity values from the upper 0.5% percentile. Since $D=\delta +I$ (Equation 1) and photondependent intensity $I$ is positive, all $D$ values have to be larger than the smallest offset intensity value. If that is not the case we add a single value $\mathrm{min}(D)1$ to the offset empirical distribution which has a negligible effect on the distribution. Bin values $\delta}_{\mathsf{s}\mathsf{a}\mathsf{m}\mathsf{p}\mathsf{l}\mathsf{e}\mathsf{s}$ and their weights $\delta}_{\mathsf{w}\mathsf{e}\mathsf{i}\mathsf{g}\mathsf{h}\mathsf{t}\mathsf{s}$ are used to construct an Empirical prior:
All simulated and experimental data sets in this work were analyzed using the prior distributions and hyperparameter values given above, which are compatible with a broad range of experimental conditions (Table 1). Many of the priors are uninformative and we anticipate that these will work well with images taken on variety of microscope hardware. However, it is possible that highly atypical microscope designs (e.g., those with effective magnifications that are suboptimal for CoSMoS) might require adjustment of some fixed hyperparameters and distributions (those in Eqs. 6a, 6b, 11, 12, 13, 15, and 16). For example, if the microscope point spread function is more than 2 pixels wide, it may be necessary to increase the range of the $w$ prior in Eq. 13. The Tapqir documentation (https://tapqir.readthedocs.io/en/stable/) gives instructions for changing the hyperparameters.
Joint distribution
Request a detailed protocolThe joint distribution of the data and all parameters is the fundamental distribution necessary to perform a Bayesian analysis. Let $\varphi $ be the set of all model parameters. The joint distribution can be expressed in a factorized form:
The Tapqir generative model is a stochastic function that describes a properly normalized joint distribution for the data and all parameters (Figure 8). In Pyro this is called ‘the model’.
Inference
Request a detailed protocolFor a Bayesian analysis, we want to obtain the posterior distribution for parameters $\varphi $ given the observed data $D$. There are three discrete parameters $z$, $\theta$, and $\delta$ that can be marginalized out exactly so that they do not appear expilictly in either the joint posterior distribution or the likelihood function. Computationally efficient marginalization is implemented using Pyro’s enumeration strategy (Obermeyer et al., 2019a) and KeOps’ kernel operations on the GPU without memory overflows (Charlier et al., 2021). Let ${\varphi}^{\prime}=\varphi \{z,\theta ,\delta \}$ be the rest of the parameters. We obtain posterior distributions of ${\varphi}^{\prime}$ using Bayes’ rule:
Note that the integral in the denominator of this expression is necessary to calculate the posterior distribution, but it is usually analytically intractable. However, variational inference provides a robust method to approximate the posterior distribution $p({\varphi}^{\mathrm{\prime}}\mid D)$ with a parameterized variational distribution $q({\varphi}^{\prime})$ (Bishop, 2006).
$q({\varphi}^{\prime})$ has the following factorization:
The variational distribution $q({\varphi}^{\prime})$ is provided as pseudocode for a generative stochastic function (Figure 9). In Pyro this is called ‘the guide’. Variational inference is sensitive to initial values of variational parameters. In Figure 9, step 1 we provide the initial values of variational parameters used in our analyses.
Calculation of spot probabilities
Request a detailed protocolVariational inference directly optimizes $q(m)\equiv {m}_{\mathsf{p}\mathsf{r}\mathsf{o}\mathsf{b}}$ (see Eq. 21 and Figure 9), which approximates $p(m\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}D)$. To obtain the marginal posterior probabilities $p(z,\theta \phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}D)$, we use a Monte Carlo sampling method:
In our calculations, we used $S$ = 25 as the number of Monte Carlo samples. Marginal probabilities $p(z\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}D)$ and $p(\theta \phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}D)$ are calculated as:
The probability, $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$, that a targetspecific fluorescence spot is present in a given image by definition is:
For simplicity in the main text and figures we suppress the conditional dependency on $D$ in $p(\theta \phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}D)$ and $p(m\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}D)$ and instead write them as $p(\theta )$ and $p(m)$, respectively.
Tapqir implementation
Request a detailed protocolThe model and variational inference method outlined above are implemented as a probabilistic program in the Pythonbased probabilistic programming language (PPL) Pyro (Foerster et al., 2018; Bingham et al., 2019; Obermeyer et al., 2019a). We use a variational approximation because exact inference is not analytically tractable for a model as complex as cosmos. As currently implemented in Pyro, variational inference is significantly faster than Monte Carlo inference methods. In Tapqir, the objective that is being optimized is the evidence lower bound (ELBO) estimator that provides unbiased gradient estimates upon differentiation. At each iteration of inference procedure we choose a random subset of AOIs and frames (minibatch), compute a differentiable ELBO estimate based on this minibatch and update the variational parameters via automatic differentiation. We use PyTorch’s Adam optimizer (Kingma and Ba, 2014) with the learning rate of $5\times {10}^{3}$ and keep other parameters at their default values.
Credible intervals and confidence intervals
Request a detailed protocolCredible intervals were calculated from posterior distribution samples as the highest density region (HDR), the narrowest interval with probability mass 95% using the pyro.ops.stats.hpdi Pyro function. Confidence intervals were calculated from bootstrap samples as the 95% HDR.
Data simulation
Request a detailed protocolSimulated data were produced using the generative model (Figure 8). Each simulation has a subset of parameters ($\pi ,\phantom{\rule{thinmathspace}{0ex}}\lambda$, $g$, ${\sigma}^{xy}$, $b$, $h$, $w,\phantom{\rule{thinmathspace}{0ex}}\delta$) set to desired values while the remaining parameters ($z,\theta$, $m$, $x$, $y$) and resulting noisy images ($D$) are sampled from distributions. The fixed parameter values and data set sizes for all simulations are provided inSupplementary file 1; Supplementary file 2; Supplementary file 3; Supplementary file 4; Supplementary file 5; Supplementary file 6.
For kinetic simulations (Figure 5, Supplementary file 5), $z$ was modeled using a discrete Markov process with the initial probability and the transition probability matrices:
where $k}_{\mathsf{o}\mathsf{n}$ and $k}_{\mathsf{o}\mathsf{f}\mathsf{f}$ are transition probabilities that numerically approximate the pseudofirstorder binding and firstorder dissociation rate constants in units of $\mathsf{s}}^{1$, respectively, assuming 1 s/frame. We assumed that the Markov process is at equilibrium and initialized the chain with the equilibrium probabilities.
Posterior predictive sampling
Request a detailed protocolFor posterior predictive checking, sampled images ($\stackrel{~}{D}$) were produced using Tapqir’s generative model (Figure 8) where model parameters were sampled from the posterior distribution $p(\varphi D)$, which was approximated by the variational distribution $q(\varphi )$:
Signaltonoise ratio
Request a detailed protocolWe define SNR as:
where ${\sigma}_{\mathsf{b}\mathsf{a}\mathsf{c}\mathsf{k}\mathsf{g}\mathsf{r}\mathsf{o}\mathsf{u}\mathsf{n}\mathsf{d}}^{2}=b\cdot g$ the variance of the background intensity, $\sigma}_{\mathsf{o}\mathsf{f}\mathsf{f}\mathsf{s}\mathsf{e}\mathsf{t}}^{2$ the variance of the offset intensity, and the mean is taken over all targetspecific spots. For experimental data, $\mathsf{s}\mathsf{i}\mathsf{g}\mathsf{n}\mathsf{a}\mathsf{l}$ is calculated as
where $\mathsf{w}\mathsf{e}\mathsf{i}\mathsf{g}\mathsf{h}\mathsf{t}$ is
For simulated data theoretical $\mathsf{s}\mathsf{i}\mathsf{g}\mathsf{n}\mathsf{a}\mathsf{l}$ is directly calculated as:
Classification accuracy statistics
Request a detailed protocolAs a metric of classification accuracy we use three commonly used statistics – recall, precision, and Matthews Correlation Coefficient (Matthews, 1975)
where TP is true positives, TN is true negatives, FP is false positives, and FN is false negatives.
Kinetic and thermodynamic analysis
Request a detailed protocolTo estimate simple binding/dissociation kinetic parameters (Figure 5C and D), we sample binary time records $z$ from the inferred $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ time records for all AOIs. For a twostate hidden Markov model, the maximumlikelihood estimates of $k}_{\mathsf{o}\mathsf{n}$ and $k}_{\mathsf{o}\mathsf{f}\mathsf{f}$ are given by:
Repeating this procedure 2,000 times gave the distributions of $k}_{\mathsf{o}\mathsf{n}$ and $k}_{\mathsf{o}\mathsf{f}\mathsf{f}$ from which we compute mean and 95% credible interval.
Similarly, to estimate mean and 95% CI of $K}_{\mathsf{e}\mathsf{q}$ (Figure 5E) we sampled π from $q(\pi )$ and for each sampled value of $\pi$ calculated $K}_{\mathsf{e}\mathsf{q}$ as:
To calculate timetofirst binding kinetics from the Tapqirderived $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ (Figure 6B, Figure 6—figure supplement 1B, Figure 6—figure supplement 2B, and Figure 6—figure supplement 3B), 2,000 binary time records $z$ were sampled from the $p(\mathsf{s}\mathsf{p}\mathsf{e}\mathsf{c}\mathsf{i}\mathsf{f}\mathsf{i}\mathsf{c})$ time record for each AOI. For each sampled time record initial absent intervals were measured and analyzed using Equation 7 in Friedman and Gelles, 2015, yielding distributions of $k}_{\mathsf{a}$, $k}_{\mathsf{n}\mathsf{s}$, and $A}_{\mathsf{f}$. Mean value and 95% credible intervals were calculated from these distributions. Initial absent intervals from ‘spotpicker’ analysis (Figure 6C, Figure 6—figure supplement 1C, Figure 6—figure supplement 2C, and Figure 6—figure supplement 3C) were analyzed as described in Friedman and Gelles, 2015, except that ontarget and offtarget data were here analyzed jointly instead of being analyzed sequentially (Friedman and Gelles, 2015). Note that the $k}_{\mathsf{n}\mathsf{s}$ values determined using the two methods are not directly comparable for several reasons, including that the nonspecific binding frequencies are effectively measured over different areas. For Tapqir, the target area is approximately $\pi {\left({\sigma}^{xy}\right)}^{2}$ (which is between 0.3 and 0.8 pixels^{2} in the different experimental data sets) and for spotpicker the area is subjectively chosen as $\pi \cdot {1.5}^{2}=7$ pixels^{2}.
Data availability
All data generated or analyzed for this study will be available at https://github.com/ordabayevy/tapqiroverleaf. That repository also includes all figures, figure supplements, and the scripts and data used to generate them. It also contains the supplemental data files and preprint manuscript text.

GithubID ordabayevy/tapqiroverleaf. Simulated and experimental data used for "Bayesian machine learning analysis of singlemolecule fluorescence colocalization images".
References

Pyro: Deep universal probabilistic programmingJournal of Machine Learning Research: JMLR 20:1–6.

Rocket launcher mechanism of collaborative actin assembly defined by singlemolecule imagingScience (New York, N.Y.) 336:1164–1168.https://doi.org/10.1126/science.1218062

Kernel operations on the GPU, with autodiff, without memory overflowsJournal of Machine Learning Research: JMLR 22:1–6.

An introduction to ROC analysisPattern Recognition Letters 27:861–874.https://doi.org/10.1016/j.patrec.2005.10.010

Multiwavelength singlemolecule fluorescence analysis of transcription mechanismsMethods (San Diego, Calif.) 86:27–36.https://doi.org/10.1016/j.ymeth.2015.05.026

BookChapman and Hall/CRC Texts in Statistical SciencePhiladelphia: CRC Press.https://doi.org/10.1201/b16018

Bayesian Inference: The Comprehensive Approach to Analyzing SingleMolecule ExperimentsAnnual Review of Biophysics 50:191–208.https://doi.org/10.1146/annurevbiophys082120103921

Biased Brownian ratcheting leads to premRNA remodeling and capture prior to firststep splicingNature Structural & Molecular Biology 20:1450–1457.https://doi.org/10.1038/nsmb.2704

Substrate degradation by the proteasome: a singlemolecule kinetic analysisScience (New York, N.Y.) 348:1250834.https://doi.org/10.1126/science.1250834

Specificity of the anaphasepromoting complex: a singlemolecule studyScience (New York, N.Y.) 348:1248737.https://doi.org/10.1126/science.1248737

Comparison of the predicted and observed secondary structure of T4 phage lysozymeBiochimica et Biophysica Acta 405:442–451.https://doi.org/10.1016/00052795(75)901099

Quantitative Aspects of Single Molecule MicroscopyIEEE Signal Processing Magazine 32:58–69.https://doi.org/10.1109/MSP.2014.2353664

Dynamic recognition of the mRNA cap by Saccharomyces cerevisiae eIF4EStructure (London, England 21:2197–2207.https://doi.org/10.1016/j.str.2013.09.016

Breaking the Concentration Barrier for SingleMolecule Fluorescence MeasurementsChemistry (Weinheim an Der Bergstrasse, Germany) 24:1002–1009.https://doi.org/10.1002/chem.201704065

A practical guide to singlemolecule FRETNature Methods 5:507–516.https://doi.org/10.1038/nmeth.1208

Probabilitybased particle detection that enables thresholdfree and robust in vivo singlemolecule trackingMolecular Biology of the Cell 26:4057–4062.https://doi.org/10.1091/mbc.E15060448

The dynamics of SecMinduced translational stallingCell Reports 7:1521–1533.https://doi.org/10.1016/j.celrep.2014.04.033

Singlemolecule approaches to characterizing kinetics of biomolecular interactionsCurrent Opinion in Biotechnology 22:75–80.https://doi.org/10.1016/j.copbio.2010.10.002

BookDigital fluorescence imaging using cooled CCD array cameras invisibleIn: Celis JE, editors. Cell Biology. New York: Academic Press. pp. 109–120.

Singlemolecule fluorescencebased studies on the dynamics, assembly and catalytic mechanism of the spliceosomeBiochemical Society Transactions 42:1211–1218.https://doi.org/10.1042/BST20140105

Structural basis of transcription initiationScience (New York, N.Y.) 338:1076–1080.https://doi.org/10.1126/science.1227786
Article and author information
Author details
Funding
National Institute of General Medical Sciences (R01GM121384)
 Jeff Gelles
 Douglas L Theobald
National Institute of General Medical Sciences (R01GM081648)
 Jeff Gelles
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This work was supported by grants R01GM121384 and R01GM081648 from the National Institute of General Medical Sciences, NIH. We thank Alex Okonechnikov for his work on an earlier version of this project. We thank Jane Kondev, Timothy M Lohman, and Timothy O Street for helpful comments on the manuscript.
Version history
 Received: September 14, 2021
 Preprint posted: October 1, 2021 (view preprint)
 Accepted: March 19, 2022
 Accepted Manuscript published: March 23, 2022 (version 1)
 Version of Record published: June 9, 2022 (version 2)
Copyright
© 2022, Ordabayev 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.
Metrics

 1,545
 views

 286
 downloads

 2
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Biochemistry and Chemical Biology
 Structural Biology and Molecular Biophysics
Vitamin B6 deficiency has been linked to cognitive impairment in human brain disorders for decades. Still, the molecular mechanisms linking vitamin B6 to these pathologies remain poorly understood, and whether vitamin B6 supplementation improves cognition is unclear as well. Pyridoxal 5’phosphate phosphatase (PDXP), an enzyme that controls levels of pyridoxal 5’phosphate (PLP), the coenzymatically active form of vitamin B6, may represent an alternative therapeutic entry point into vitamin B6associated pathologies. However, pharmacological PDXP inhibitors to test this concept are lacking. We now identify a PDXP and agedependent decline of PLP levels in the murine hippocampus that provides a rationale for the development of PDXP inhibitors. Using a combination of smallmolecule screening, protein crystallography, and biolayer interferometry, we discover, visualize, and analyze 7,8dihydroxyflavone (7,8DHF) as a direct and potent PDXP inhibitor. 7,8DHF binds and reversibly inhibits PDXP with low micromolar affinity and submicromolar potency. In mouse hippocampal neurons, 7,8DHF increases PLP in a PDXPdependent manner. These findings validate PDXP as a druggable target. Of note, 7,8DHF is a wellstudied molecule in brain disorder models, although its mechanism of action is actively debated. Our discovery of 7,8DHF as a PDXP inhibitor offers novel mechanistic insights into the controversy surrounding 7,8DHFmediated effects in the brain.

 Biochemistry and Chemical Biology
 Stem Cells and Regenerative Medicine
Molecules that facilitate targeted protein degradation (TPD) offer great promise as novel therapeutics. The human hepatic lectin asialoglycoprotein receptor (ASGR) is selectively expressed on hepatocytes. We have previously engineered an antiASGR1 antibodymutant RSPO2 (RSPO2RA) fusion protein (called SWEETS) to drive tissuespecific degradation of ZNRF3/RNF43 E3 ubiquitin ligases, which achieved hepatocytespecific enhanced Wnt signaling, proliferation, and restored liver function in mouse models, and an antibody–RSPO2RA fusion molecule is currently in human clinical trials. In the current study, we identified two new ASGR1 and ASGR1/2specific antibodies, 8M24 and 8G8. Highresolution crystal structures of ASGR1:8M24 and ASGR2:8G8 complexes revealed that these antibodies bind to distinct epitopes on opposing sides of ASGR, away from the substratebinding site. Both antibodies enhanced Wnt activity when assembled as SWEETS molecules with RSPO2RA through specific effects sequestering E3 ligases. In addition, 8M24RSPO2RA and 8G8RSPO2RA efficiently downregulate ASGR1 through TPD mechanisms. These results demonstrate the possibility of combining different therapeutic effects and degradation mechanisms in a single molecule.