
==== Front
Commun Biol
Commun Biol
Communications Biology
2399-3642
Nature Publishing Group UK London

6772
10.1038/s42003-024-06772-8
Article
MultiMatch: geometry-informed colocalization in multi-color super-resolution microscopy
http://orcid.org/0000-0003-3766-0945
Naas Julia 12
Nies Giacomo 34
http://orcid.org/0000-0002-2434-9878
Li Housen 34
Stoldt Stefan 456
Schmitzer Bernhard 7
http://orcid.org/0000-0002-8028-3121
Jakobs Stefan 4568
http://orcid.org/0000-0002-9181-9331
Munk Axel munk@math.uni-goettingen.de

34
1 grid.22937.3d 0000 0000 9259 8492 Center for Integrative Bioinformatics Vienna (CIBIV), Max Perutz Labs, University of Vienna and Medical University of Vienna, Vienna, Austria
2 grid.22937.3d 0000 0000 9259 8492 Vienna Biocenter PhD Program, a Doctoral School of the University of Vienna and Medical University of Vienna, Vienna, Austria
3 https://ror.org/01y9bpm73 grid.7450.6 0000 0001 2364 4210 Institute for Mathematical Stochastics, University of Göttingen, Göttingen, Germany
4 https://ror.org/01y9bpm73 grid.7450.6 0000 0001 2364 4210 Cluster of Excellence ‘Multiscale Bioimaging: from Molecular Machines to Networks of Excitable Cells’ (MBExC), University of Göttingen, Göttingen, Germany
5 https://ror.org/03av75f26 Department of NanoBiophotonics, Max Planck Institute for Multidisciplinary Sciences, Göttingen, Germany
6 https://ror.org/021ft0n22 grid.411984.1 0000 0001 0482 5331 Clinic of Neurology, University Medical Center Göttingen, Göttingen, Germany
7 https://ror.org/01y9bpm73 grid.7450.6 0000 0001 2364 4210 Institute for Computer Science, University of Göttingen, Göttingen, Germany
8 https://ror.org/01s1h3j07 grid.510864.e Fraunhofer Institute for Translational Medicine and Pharmacology ITMP, Translational Neuroinflammation and Automated Microscopy TNM, Göttingen, Germany
13 9 2024
13 9 2024
2024
7 113928 2 2024
22 8 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
With recent advances in multi-color super-resolution light microscopy, it is possible to simultaneously visualize multiple subunits within biological structures at nanometer resolution. To optimally evaluate and interpret spatial proximity of stainings on such an image, colocalization analysis tools have to be able to integrate prior knowledge on the local geometry of the recorded biological complex. We present MultiMatch to analyze the abundance and location of chain-like particle arrangements in multi-color microscopy based on multi-marginal optimal unbalanced transport methodology. Our object-based colocalization model statistically addresses the effect of incomplete labeling efficiencies enabling inference on existent, but not fully observable particle chains. We showcase that MultiMatch is able to consistently recover existing chain structures in three-color STED images of DNA origami nanorulers and outperforms geometry-uninformed triplet colocalization methods in this task. MultiMatch generalizes to an arbitrary number of color channels and is provided as a user-friendly Python package comprising colocalization visualizations.

MultiMatch performs geometry-informed colocalization analysis for multi-color microscopy based on multi-marginal optimal unbalanced transport methodology.

Subject terms

Software
Image processing
Statistical methods
Super-resolution microscopy
Fluorescence imaging
501100002428 Austrian Science Fund (Fonds zur Förderung der Wissenschaftlichen Forschung) F78 (to Arndt von Haeseler) 501100001659 Deutsche Forschungsgemeinschaft (German Research Foundation) Germany's Excellence Strategy EXC 2067/1-390729940 501100001659 Deutsche Forschungsgemeinschaft (German Research Foundation) CRC 1456 792 'Mathematics of Experiment' (Project Number B04, C06) European Research Council (ERC), Advanced Grant, Call: ERC-2018-AdG, Grant: 835102, 2019-2024 (to Stefan Jakobs)issue-copyright-statement© Springer Nature Limited 2024
==== Body
pmcIntroduction

Colocalization analysis aims to unravel the interconnection and interaction network between two or more groups of particles based on their spatial proximity in a microscopy image. By visualizing biological structures, like DNA, RNA and proteins, that are only a few nanometers in size, colocalization analysis makes it possible to study a wide range of biological processes, such as DNA replication and the transcription of genes1, nuclear import of splicing factors2 or the dynamics of cargo sorting zones in the trans-Golgi networks of plants3, to name only a few.

In the following, we will denote any objects of interest that are depicted within a microscopy image, e.g., proteins as well as loci on DNA or RNA strands, as particles. In fluorescence light microscopy, such particles are stained, i.e., in case they do not already intrinsically fluoresce, they are labeled with fluorophores, which in turn are excited by an external light source. The emitted fluorescence radiation then can be imaged via several microscopy technologies.

Diffraction unlimited super-resolution fluorescence microscopy technologies, also called nanoscopy, are classified into two broad concepts4:

In coordinate-stochastic microscopy, fluorophores within the sample are stochastically excited resulting in a temporally resolved blinking dynamic5–7, which allows to spatially separate fluorophores. Their coordinates are estimated by means of the detected radiation peak, yielding a list of coordinates of detected fluorophores as output data. If only one fluorophore is detected for one particle, the output translates into a list of particle coordinates. Else, fluorophore coordinates can be aggregated in order to localize the particle of interest in the imaged biological sample.

In scanning-based microscopy methods such as Stimulated Emission Depletion (STED)8–10, the fluorescence distribution is stored as an intensity matrix, in which every entry encodes the detected radiation within a respective pixel of the microscopy image. To obtain coordinate estimates of particle positions, object detection algorithms have to be applied to the intensity matrix.

In order to study possible particle interactions or connections, stainings with different fluorescent markers are recorded in different color channels. Particles colocalize, if they are spatially closer than or equal to a colocalization distance, which heavily depends on the underlying biological setting and might be unknown prior to colocalization analysis11.

Colocalization methods are divided in two categories based on the input data format they require:

Pixel-based colocalization methods take an intensity matrix as input and compare the pixel intensities across color channels, e.g., by utilizing overlap, correlation or intensity transport analysis. Such approaches are thus only applicable for scanning-based images and examples for well-established methods are Mander’s Colocalization Coefficient12,13, Pearson’s Correlation Coefficient14, BlobProb15, SACA16, and OTC curves17.

Object-based colocalization methods, which our method MultiMatch classifies as, require the coordinates of particles and evaluate their distances, where pairwise particle distances can be defined in several ways18. Examples for other object-based tools are ConditionalColoc18 and Ripley’s K based methods19,20 as SODA21.

While nanoscopy for dual-color stainings is well studied for a long time, multi-color imaging including three or more stainings has received increased attention more recently since it allows simultaneous measurements of multiple particle types. There is a steadily increasing number of published multi-color STED microscopy datasets22–28, of other super-resolution microscopy methods29,30 and the development of appropriate labeling methods allowing for an ever-increasing number of channels is ongoing24,30–33.

However, most pixel- and object-based colocalization tools are designed for and therefore limited to the analysis of two-color stainings. Applying them to multi-color images is not an obvious task: Particle arrangements with more than two different particle types can occur in different configurations and depending on the biological context, some may be of interest and others may simply not exist in the imaged sample. A geometry-uninformed, pairwise analysis of all possible channel combinations34, as well as the few established methods that are explicitly presented as multi-color pixel-based15,35–37 and object-based18,21,38 colocalization tools are prone to overestimate colocalization, as soon as the biological complex of interest has a fixed geometry and stoichiometry, as we can show in a simulation study. To exploit the full potential of multi-color microscopy imaging in such a situation, it is therefore beneficial to actively incorporate prior knowledge of the local geometry into the colocalization analysis.

To this end, we introduce MultiMatch, a widely applicable colocalization methodology based on optimal transport theory, which is especially tailored to detect chain-like, one-to-one particle arrangements. Integrating this type of colocalization geometry optimizes the multi-color colocalization analysis of quadruples, triplets, pairs, and singlets, as they appear when marking different loci of a chain-like molecule with multi-color stainings. MultiMatch is able to statistically address the effect of incomplete labeling efficiencies on the detection results and includes statistical guarantees on the estimated number of structures of interest. It is provided as computationally efficient Python package allowing for a user-friendly visualization of colocalization results via colocalization curves and the exploratory napari viewer39.

Results

Chain-like particle assembly detection with MultiMatch

One exemplary biological framework, in which the localization of chain-like particle arrangements is especially insightful, is the highly condensed mammalian mitochondrial genome: It is transcribed from both strands of the mitochondrial DNA as long polycistronic transcripts that have to undergo multiple processing steps, including endonucleolytic cleavage, in order to get to the different functional RNA species. Transcription of the heavy strand leads to polycistronic primary transcripts containing the premature mRNAs of 12 of the 13 oxidative phosphorylation (OXPHOS) subunits encoded in the mitochondrial genome. Labeling more than two of the mRNAs within such a primary construct, in combination with our colocalization approach, can significantly contribute to our understanding of the post-transcriptional processing steps and their dynamics, that lead to the generation of matured mRNA molecules40,41.

We consider a particle arrangement as chain-like when exactly one particle of each type is stringed together in an ordered fashion and pairwise distances of chain-neighbors are smaller than or equal to a maximal colocalization threshold t. In MultiMatch we implemented the distance between reference points, i.e., the center of detected particles, as t by default (Fig. 1A). This approach is especially suited for particles of small size or in case the center of the particle is a suitable representation for its location on the microscopy image. However, we allow the user to also input arbitrary particle-to-particle distance matrices18 as they are output by alternative particle detection and segmentation algorithms.Fig. 1 MultiMatch workflow to detect chain-like particle arrangements in multi-color microscopy images.

A After microscopy imaging (0) and object detection (1), the distances between channel-specific lists of reference points or a user-defined distance matrix are input to the optimal matching procedure. Restricted on particle pairs with distance smaller or equal than colocalization distance t (2), MultiMatch either outputs the maximal number of triplets and subsequently pairs (Mode I) or simultaneously searches for triplets and pairs (Mode II) (3). MultiMatch provides the localization and number of detected chains for a known or abundance curves for a range of colocalization distances t (4). For known incomplete labeling efficiencies true abundances can be estimated with confidence statements (5) (“Methods” section). B If more than two different particle types are involved, multiple geometric colocalization patterns can emerge. In case the chain is a substructure of the colocalization geometry of interest, its detection will help to localize and quantify colocalization events. C Structures of interest in three-color colocalization analysis for chain-like, one-to-one particle interactions and fixed particle type order. All pairwise distances between neighboring particles in a chain are smaller or equal than colocalization distance t. D Exemplary MultiMatch output for an experimental STED image of DNA origami nanoruler structures (as sketched in C) in the interactive napari viewer39.

Even if the biological complex of interest itself is not chain-like, chain detection still can give substantial insights on the abundance and location of colocalization events inside a microscopy image as soon as the chain is a substructure of the colocalization geometry (Fig. 1B). The converse, on the other hand, does not hold true in general.

To fix the chain order of particles, we will refer to color channels, in which the respective particle type was imaged, as channel A, B, C, D etc. For simplicity, we will explain the main methodology for a three-color setting in what follows, but MultiMatch is applicable to an arbitrary number of color channels, which we showcase in the evaluation of simulated four-color STED images. We stress, that our software is already designed to process any number of channels (“Methods” section).

All configurations resulting from a three-color staining of an chain-like molecule are sketched in Fig. 1C, where we assume the following unknown abundances n=nABC,nAB,nBC,nA,nB,nC of chain-like assemblies, wherenABC is the number of true ABC triplets,

nAB, nBC is the numbers of true AB and BC pairs,

nA, nB, nC is the numbers of true A, B, and C singlets.

Optimal transport (OT) theory42 has a wide range of applications throughout statistics43, data science, and machine learning44. Generally, OT aims to allocate (transport plan) one mass distribution into another by minimizing the transportation cost arising from moving one mass unit from one location to another. Applied to fluorescence intensity distributions on a pixel grid and using the euclidean distances between pixels as transportation cost, OT introduces an intuitive distance between two microscopy images and could already successfully be utilized in the context of pixel-based, dual-color colocalization methods17,45.

For object-based analysis, reference points of detected particles can also be interpreted as support points of mass one of a (discrete) two-dimensional distribution. For only two color channels with the same number of particles the standard OT problem simply assigns each particle from the first channel to one particle from the second channel while minimizing the total sum of Euclidean matching costs. This can be also generalized for other particle-to-particle distance matrices, in case the Euclidean distance between particle reference points is not suitable to represent particle proximity. We can obtain an optimal matching between more than two particle types by multi-marginal OT46,47 and at the same time account for the not necessarily equal numbers of support points per channel by utilizing an unbalanced OT formulation48 (“Methods” section). A combination of both OT generalizations, i.e., multi-marginal optimal unbalanced transport problems, have been recently discussed in the literature49–52.

In this manner, the basic concept of MultiMatch can be interpreted as a linear assignment problem as described, e.g., in the field of object tracking53–56. In contrast to methods of this research field, we explicitly formulate the matching problem as a function of the colocalization threshold, allowing to plot the chain abundances dependent on a range of t (“Methods” section). We utilize the equivalence of the optimal transport methodology to a network flow problem to overcome the otherwise prohibitively high computational complexity of its corresponding linear program formulation57 (“Methods” section, Supplementary Note 1, and Supplementary Fig. 1).

MultiMatch provides two different modes to solve the particle matching problem (Fig. 1A(3)):

Mode I: By restricting a k-marginal optimal unbalanced transport problem to particle pairs with a distance smaller than t and introducing a chain-cost that only considers distances between neighboring particle types (“Methods” section), the resulting OT plan encodes the maximal number of, for k = 3, triplets within the nanoscopy image. If requested, the matching process is subsequently repeated on the remaining particles to detect yet unresolved AB and BC pairs, respectively.

Mode II: This mode only detects AB, BC, etc. pairs by solving respective two-marginal unbalanced OT problems. Subsequently, the two-marginal OT matchings are coupled to chain structures: For k = 3, all pairs occupying the same intermediate particle are redefined as respective ABC triplet.

Depending on the underlying biological experiment, the user can select the appropriate mode for colocalization analysis: Mode I prioritizes the detection of a predefined chain structure of choice. For example, if a user aims to analyze triplets, Mode I will detect a triplet as soon as three particles A, B, and C are close enough to each other – even if another particle A or C is nearby that would allow to match two pairs instead of one triplet (as depicted in Fig. 1A(3)). If k > 3 and the user wants to detect multiple chain structures, one needs to set a prioritization order for Mode I. For example, for k = 4 and after ABCD quadruplet detection, one can search either for ABC or BCD triplets next. Depending on the order, the final matching results may change as soon as some particles cannot be uniquely assigned to one particle arrangement.

Mode II, on the other hand, does not need a predefined prioritization order of structures for subsequent matching steps, hence it does not overemphasize structures that are matched in the earlier steps. It is useful in case we do not have any prior knowledge on which structures might appear in the microscopy image and we do not want to prioritize any chain structures.

In the evaluation of experimental and simulated three-color STED microscopy images we show that for sparse particle distributions and mixed singlets, pairs, and triplet ratios the differences in detected abundances between the two modes is neglectable (Supplementary Note 2, Supplementary Fig. 2). However, in case of dense particle distributions (Supplementary Note 3 and Supplementary Figs. 3–5), or in case we know in advance that only one chain structure exists in the biological context, the multi-marginal approach of Mode I, which is also the default setting in the MultiMatch tool, outperforms the pairwise matching approach of Mode II.

MultiMatch outputs detected abundances w = (wABC, wAB, wBC, wA, wB, wC) for a known colocalization distance t and depicts configuration positions on the respective microscopy image allowing further investigation on the spatial distribution of recorded biological complexes. If t is unknown (optionally channel-wise scaled) abundance curves w(t) are output for a user-defined range of t values. MultiMatch is compatible with the interactive Graphical User Interface of napari (Fig. 1D) enabling the visual evaluation of structure locations for different t values in form of a colocalization threshold slider.

The differentiation between triplets, pairs, and singlets within a microscopy image is additionally hindered by incomplete labeling efficiencies and point detection artifacts. This is a notorious problem in fluorescence microscopy, e.g., described in Hummert et al.58, and missing detections can add an unpredictable bias toward systematic underestimation of triplet numbers and overestimation of singlet abundances, if not corrected. Currently, the problem of incomplete labeling efficiency is barely addressed in the field of colocalization analysis. Therefore, we propose a statistical framework to correct for incomplete labeling efficiencies and introduce an unbiased estimator n^(t) of true chain-structure abundances and confidence statements on the estimated quantities (“Methods” section, Supplementary Note 4, Supplementary Fig. 6, Supplementary Table 1).

An overview on the full workflow of MultiMatch from microscopy image to abundance curves is depicted in Fig. 1A.

Simulation study

To systematically evaluate the performance of MultiMatch against compatible colocalization methods, we simulated 100 microscopy images for each of three scenarios with different combinations of singlets, pairs, and triplet abundances. For this simulation study, we decreased the noise level to a minimum to allow a fair comparison despite different point detection tools implemented in the respective colocalization tools. Also, we amplified simulating linear triplet structures over randomly folded triplets (“Methods” section). For every simulated image,Scenario 1: 50 singlets of each type A, B, and C were simulated.

Scenario 2: 50 A, B, and C singlets and 50 AB and BC pairs were simulated, respectively.

Scenario 3: 100 triplets and 50 AB and BC pairs and 50 A, B, and C singlets were simulated, respectively.

Exemplary, simulated images and the results of the simulation study for a fixed colocalization threshold of t = 5 pixels are shown in Fig. 2A, B. Analysis results for all considered methods across a range of colocalization thresholds are presented in Supplementary Note 3 and Supplementary Fig. 3.Fig. 2 Simulation study for three-color microscopy images with three combinations of chain structures.

In each scenario 100 independent STED images and different abundances of triplets, pairs, and singlets were simulated with 100% labeling efficiency. A Method specific boxplots of the errors in detected relative (scaled by the total number of points in channel B) structure abundances are displayed. The error is computed by subtracting true relative abundance from detected relative abundances. In Scenario 1 only A, B, and C singlets, in Scenario 2 all possible singlets as well as AB and BC pairs and in Scenario 3 ABC triplets, AB, BC pairs and A, B, and C singlets were simulated. B Simulated STED images for Scenarios 1, 2, and 3 with respective image details. For visualization purposes, contrast stretching and increasement of image brightness was applied.

As a representative of pixel-based methods, we include BlobProb15, which counts the number of colocalized intensity blobs, i.e., groups of neighboring pixels with high intensity. In each channel, blobs are detected via image segmentation and for each blob the local intensity maximum is defined as reference particle coordinate. A blob pair colocalizes if the first reference point lies within the second blob and vice versa. Triplet colocalization is detected if all involved reference points are included in all three blobs. SODA21 is an object-based method, which uses the Ripley’s K function19 and computes the coupling probability of point pairs based on marked-point process theory. In the most recently published method ConditionalColoc18 particles are defined as colocalized as soon as their distance is below a maximal colocalization radius. Then, utilizing Bayes’ Theorem, (conditional) probabilities are computed and assigned for triplet and pair colocalization. We experienced that ConditionalColoc, although aiming to output probabilities, in some cases yields values greater than one and hence the errors in relative abundance detection are not bounded by one as well. For a better comparison, we restricted the respective results to values between -0.5 and 1 in Fig. 2A and show ConditionalColoc outliers in Supplementary Note 5, Supplementary Fig. 7.

In none of the above methods triplet colocalization is restricted to one-to-one interactions. This has barely any negative effect on the detection of singlets in Scenario 1, where no additional pairs and triplets occur. Apart from few outliers of overestimation in pairs and triplet abundances in ConditionalColoc and SODA, all considered colocalization measures show consistently low errors with small variability. The maximal median error in relative abundances of 0.03 in Scenario 1 is obtained by ConditionalColoc in the detection of AB as well as BC pairs.

In Scenarios 2 and 3 on the other hand, we observe a consistent overestimation of relative pairs and triplet abundances in object-based methods SODA and ConditionalColoc, since one particle can be included in several structures at the same time. Additionally, in Scenario 2 SODA exhibits a larger variation in pairs abundances, resulting in median errors 0.14 in both AB and BC pairs with interquartile ranges of 0.16, respectively. In Scenario 3 the variation in abundance detection decreased and median errors are 0.1 for ABC triplets and 0.04 for AB as well as BC pairs. ConditionalColoc performances worst in Scenario 3 yielding a median error of 0.48 for ABC triplets.

The pixel-based method BlobProb mostly obtains zero relative abundances of triplets and pairs across all three scenarios and hence severely underestimates the triplet and pair configurations within the simulated images. This is due to the high resolution in the simulation setup, which was chosen to mimic conventional STED imaging. If particles are small and their respective intensity blobs do not overlap, BlobProb does not detect any colocalization.

MultiMatch on the other hand searches for optimal matches on a global scale while considering the local geometry of chain-like particle assemblies. It consistently recovers the ground truth abundances for each simulation scenario. The maximal median error across all scenarios and chain structures for both Modes of MultiMatch is 0.03 with a maximal interquartile range in errors of 0.04.

Apart from above considered, already established colocalization methods, we also implemented a Nearest Neighbor Matching as comparable object-based method. We can show that greedily matching particle pairs based on local optima leads to underestimation of ABC triplets in dense particle distributions (Supplementary Note 3 and Supplementary Fig. 4).

Incomplete labeling efficiencies and point detection errors

In experimental STED microscopy, typically it is impossible to record all existing particles of interest. This can, for example, be due to the fluorescent marker not being successfully attached to the probe or a flawed point detection. All such scenarios resulting in a failure of particle detection for simplicity will be summarized under incomplete labeling efficiency hereafter.

If only singlets were to be counted in multi-color images with the same labeling efficiency across channels, the relative abundance could still be estimated consistently. However, as soon as configurations of two or more particle types are to be recovered, incomplete labeling efficiencies can lead to under- and overestimation of structures. Figure 3A shows that a triplet can be erroneously detected as pair or singlet or not at all, which can introduce a severe bias. However, if the labeling efficiencies are known, the detection success of a particle can be modeled with a Bernoulli distribution, which allows the definition of an unbiased estimator n^ for the vector of true chain structures abundances n. This approach allows for constructing multi-dimensional joint confidence ellipsoids covering n with a given significance level, e.g., α = 0.1 (Fig. 3B, C). The multi-dimensional confidence ellipsoid then can be respectively projected onto one dimension to obtain structure-specific confidence intervals or bands for a range of t values, while fixing the estimated abundances of all other considered structures (“Methods” section).Fig. 3 Chain-like particle structure detection is influenced by incomplete labeling efficiencies and structure rotation.

A Because of channel-specific incomplete labeling efficiencies, triplets and pairs can erroneously counted to other structure abundances. B For entrywise large enough n, estimator n^ is approximately multi-dimensional normally distributed: Estimated abundances of 10,000 independent simulations with labeling efficiencies sA = sB = sC = 0.95 and true abundances nABC = 500, nAB = nBC = nA = nB = nC = 50 (“Methods” section). The respective 3-dimensional, normal 90% quantile ellipsoid is plotted. C Estimated abundance curves for one of the experimental multi-color STED images in Setting 3 with additional confidence bands for significance level α = 0.1. D Restricted image resolution and 3-dimensional rotation of particle arrangements lead to variability in the observed colocalization thresholds: simulation study of 100 independent images only containing one triplet with pairwise distances set to 70 nm = 2.8 pixels per image (100% complete labeling efficiency, ”Methods” section).

Note that microscopy images are also influenced by other sources of noise that complicates the detection of chain-like particles as we show in Fig. 3D: In this small study we simulated 100 STED images containing only one triplet (nABC = 1) and observe that the discrete nature of the pixel grid can effect on the accuracy of the measured distance between particles and hence the stabilization behavior of colocalization curves. For square pixels with side length l the worst case for pairwise particle comparisons is 2l.

Evaluation of experimental STED images

Chain-like particle structures occur within several biological complexes. To showcase the performance of our method on experimentally retrieved data we used one-, two-, and three-color nanorulers. Nanorulers are DNA-origamis with a predefined distance between spots at which 20 fluorophores are attached and hence, as their name suggests, can be used as rulers inside a microscopy image1,59–61. For this experimental setup, we chose nanorulers with pairwise distances between neighboring spots of 70 nanometers (nm). For each chain structure (as depicted in Fig. 1C), respective nanoruler origamis are available in separate solutions, which allows us to control whether in an experiment we record singlets, pairs or triplets only or a combination of those structures. We performed three experiments:Setting 1: The experiment consists of all three single marker nanorulers (22 images in total). We expect to detect no pairs or triplets, i.e., wABC = wAB = wBC = 0.

Setting 2: The experiment consists of all three singlets, two pairs and triplet marker nanoruler solutions (22 images in total). We expect to detect all possible configurations, i.e., A, B, and C singlets, AB and BC pairs as well as ABC triplets.

Setting 3: The experiment consists of only triplet marker nanorulers (12 images in total). We expect to detect ABC triplets only, i.e., wAB = wAB = wA = wB = wC = 0.

For each experimental setting we recorded STED images of size 400 × 400 pixels with a pixel size of 25 × 25 nm. In channel A, stainings with Star Red 640 nm are recorded, in channel B, stainings with Alexa 488 and in channel C, stainings with Alexa 594. Note, however, that the exact numbers of nanorulers within a recorded STED image is unknown. Due to misfolding and clumping of nanorulers and different nanoruler immobilization rates for each STED image one cannot compute a fixed unit of nanorulers per microscopy image and experiment.

The results of the colocalization analysis for all three settings (with default MultiMatch Mode I) are shown in Fig. 4A via relative abundance curves with standard deviation bands quantifying variation across images within the same setting. Here, we used MultiMatch Mode I and included the analysis with Mode II showing comparable results, but slightly underestimating the number of triplets in Setting 3, in Supplementary Note 2, Supplementary Fig. 2. Exemplary images for each setting are shown in Fig. 4B.Fig. 4 MultiMatch Mode I relative abundance curves w(t) for experimental STED images.

For each setting the solid curves are mean relative abundances with standard deviation bands across a range of colocalization threshold t from 0 to 10 pixels (25 nm = 1pixel). The abundances are scaled by the total number of points detected in channel B. Additionally, incomplete labeling efficiency (90% in each channel) corrected abundances are plotted as dotted curves. The true colocalization distance of 70 nm within nanoruler structures is depicted as vertical line. A Setting 1: mean abundance curves for only singlets consistently show the expected 0% relative triplet and pair abundances (22 independent experimental STED images). Setting 2: triplets, pairs, and singlet nanoruler are detected with stable abundances for ~t ≥ 4 pixels (22 independent experimental STED images). Setting 3: mean abundance curves for analyzing the triplet nanoruler solution only. The incorporation of incomplete labeling efficiency clearly corrects the relative triplet abundance towards the in this setup expected 100% (12 independent experimental STED images). B Representative STED images for Settings 1, 2, and 3 with image details. For visualization purposes, contrast stretching and increasement of image brightness was applied.

For Setting 1 we can appreciate that, as expected, across a range of t values only a few pairs and triplets are detected (Fig. 4A). The rise of relative abundance curves is unavoidable for large t, since the probability increases that randomly scattered particles are matched. In Setting 2, despite experimental variation, we clearly recover all supplied nanoruler structures. Even more, colocalization curves are still stabilizing for a colocalization threshold t greater than approximately 4 pixels (=100 nm): For t > 100 nm ABC triplets are approximately detected with relative abundance of 0.32, AB pairs with 0.16 and BC pairs with 0.42 relative abundance, yielding a relative amount of 0.1 unmatched B singlets. The relative abundance curves of all structures reach a plateau at approximately t ≥ 4 pixels (= 100 nm), i.e., the slope of all curves within the same setting decreases rapidly. In Setting 3, as expected, the relative abundances of AB and BC pairs converge to zero while triplets are the dominantly detected structure for t ≥ 4 pixels.

Notably, in Settings 2 and 3 stable abundance curves are reached at around 100 nm, which is 30 nm more than the experimentally fixed, maximal distance between neighboring fluorophore spots in the nanoruler structures. This effect can be explained by the still limited resolution in the microscopy image and can be reproduced via simulation (Fig. 3D).

Limited resolution alone does not explain why 20%–30% of detected B particles (for t ≥ 5 pixels) are not matched to a triplet in Setting 3: The attachment of a single fluorophore to a nanoruler spot is expected to have a success probability of 85% to 90% and hence at least one fluorophore should be attached to each spot in almost 100% of all cases. Still, due to the above-described experimental variation in nanoruler imaging and additional errors in point detection, especially due to nanoruler clumping, the overall success rate of fluorophore spot detection is incomplete. Hence, we erroneously detect pairs instead of triplets or singlets due to noise. As in Setting 1 those artifacts will be matched into triplets for large enough t.

For simplicity, we model a 90% labeling efficiency across all three-color channels in the experimental STED setup. The estimated abundance curves n^(t) (dotted lines in Fig. 4A), in Setting 3 visibly correct the measurements towards the expected relative abundances. Additional confidence bands around n^ allow to infer on the robustness of the abundance estimation as presented in (Fig. 3C) for one of the experimental STED images of Setting 3.

Evaluation of simulated four-color STED images

MultiMatch is applicable to an arbitrary number of color channels, which we showcase in a second simulation study with an adapted simulation setup for quadruples, triplets, pairs, and singlets in simulated four-color STED microscopy images. In contrast to the first simulation study simulating triplets, tuples, and singlets, we additionally challenged our MultiMatch tool with an increased noise level and by allowing arbitrarily curled chain structures (“Methods” section). In Fig. 5A–D we show the colocalization analysis results of two simulation scenarios:Scenario I: We simulated 50 ABCD quadruples, 30 ABC triplets, 20 AB pairs and 30 C and D singlets, respectively, to mimic a chain-like molecule being split at loci C and D.

Scenario II: We simulated 100 ABCD quadruples and no triplets, pairs nor singlets

Fig. 5 MultiMatch Mode II absolute abundance curves w(t) and estimation results n^(t) for simulated four-color STED images.

For each simulated scenario 100 independent images were simulated with complete labeling efficiency (sA = sB = sC = sD = 1) and with incomplete labeling efficiency (sA = sB = sC = sD = 0.95), respectively. Solid curves are mean absolute detected abundances with standard deviation bands across a range of colocalization thresholds t from 0 to 10 pixels (25 nm = 1pixel). Corrected abundances are plotted as dotted curves. A Scenario I: a mixture of ABCD quadruplets, ABC triplets, AB pairs and C, D singlets were simulated. All curves stabilize at approximately t = 4 pixel close to the true simulated number of structures. For images with incomplete labeling efficiency uncorrected detected abundances plus standard deviation bands are plotted as solid curves showing consistent underestimation of quadruples. Corrected abundances recover the true number of simulated structures. B For one exemplary STED image of Scenario I simulated with incomplete labeling efficiency, corrected abundance curves and corresponding confidence bands are shown. C, D show the same analysis as shown in A and B but for Scenario II: only ABCD quadruplets were simulated. E Representative STED images for Scenarios I and Scenario II with image details. For visualization purposes, contrast stretching and increasement of image brightness was applied.

Exemplary images of both simulation scenarios are shown in Fig. 5E and three additional simulations setups are shown in Supplementary Note 3, Supplementary Fig. 5. For each scenario, we simulated 100 images with full labeling efficiencies (sA = sB = sC = sD = 1) and 100 images with incomplete labeling efficiencies (sA = sB = sC = sD = 0.95) by randomly deleting 5% of points simulated in the prior, full labeling efficiency simulation in each channel.

For this analysis we applied MultiMatch Mode II, i.e., allowing the detection of both ABC as well as BCD triplets and AB, BC and CD pairs without any prioritization order of chain structures. Again, also in the case of four-color images, we can appreciate that MultiMatch consistently recovers true abundances of quadruplets in case of full labeling efficiencies. Absolute abundance curves, as also described in the analysis of our experimental dataset in Fig. 4, stabilize for approximately t = 4 pixels. For images simulated with incomplete labeling efficiencies, the colocalization curves show underestimation of quadruplets as expected. With our statistical framework we again can visibly correct the colocalization curves towards the true, simulated structures abundances and additionally gain confidence bands confirming the stability of our estimator.

For denser distributions, as shown in Supplementary Note 3, Supplementary Fig. 5C-F, we can observe that 1. MultiMatch II misses quadruples for the sake of closer particle pairs, and 2. similar to the experimental nanoruler analysis depends on the performance of the point detection and hence the noise level in the microscopy image. If consistent noise challenges the point detection, abundance curves still stabilize, but the plateau shows a smaller number of matched quadruples than simulated in absolute numbers. Hence, we advise user of MultiMatch to check the noise level of the microscopy image and the point detection result with the interactive napari viewer (Fig. 1D and Supplementary Fig 5G) and if necessary evaluate channel-wise scaled, relative instead of absolute abundances.

Discussion

In this article we introduce multi-marginal optimal unbalanced transport methodology for geometry-informed, multi-color colocalization analysis. We are able to show, that for the analysis of more than two color channels, it is crucial to take into account the colocalization geometry of the biological complex.

By either choosing chain costs in a multi-marginal OT problem (Mode I) or coupling consecutive two-marginal OT matchings (Mode II), MultiMatch successfully detects k-chain particle assemblies such as quadruples, triplets, pairs, and singlets, as they appear when staining multiple loci on chain-like molecules like DNA or RNA strands. Both modes have their advantages, which depend on the number of particles imaged and prior knowledge on the biological context: Mode I is best for detecting one chain structure of choice and is more robust in dense particle distributions. When the particle distribution is sparser and multiple chain structures in the imaged biological setting are of interest, Mode II is suited to detect them without any predefined prioritization order.

Since often the true colocalization distance is unknown, MultiMatch results can be output as structure-wise relative or absolute abundance curves across a range of colocalization thresholds t. In our simulation studies as well as our experimental settings we could show, that output curves stabilize close to ground truth abundances.

However, as for all object-based colocalization methods, the performance MultiMatch scales with the noise level of the microscopy image, the performance of the object detection and the resolution of the microscopy. Abundance curve plateaus can be less clear in case the microscopy image contains detected singlets of different particles types. In this case, the larger t, the more far away singlets are matched. In such cases it might be unclear, whether singlets truly exist in the biological sample or whether they are an artifact of the experiment and image processing. For such cases, we advise to observe the quality of the microscopy image with the MultiMatch compatible, interactive napari viewer.

Our network flow implementation significantly decreases computational costs compared to standard approaches solving comparable OT problems and comparable colocalization tools (“Methods” section). The simulation studies show that as soon as we have prior knowledge on the chain colocalization geometry, MultiMatch, in contrast to other triplet colocalization methods, is robust against overestimation of triplets with chain geometry since it only considers one-to-one interactions. MultiMatch is also tested on experimental STED images of different nanoruler combinations and can correct structure abundances for predefined incomplete labeling efficiencies and point detection errors, where confidence bands allow further inference on the estimated abundances.

All experimental studies have been performed for k = 3 color channels. However, in many scientific fields the detection of k-chains for larger k is of interest. The mathematical and statistical frameworks allow straight-forward generalization (“Methods” section) and we exemplarily show successful detection results for simulated four-color STED images. With current technical standards, the experimental setup of multi-color nanoscopy imaging is still challenging, costly and time consuming, but in view of further technological improvements our algorithm is already applicable for the evaluation of this type of experimental setups, and especially promising in view of recent developments in super-resolution microscopy with a resolution of a few nanometers and below62,63.

In the same way channel specific colocalization thresholds as tAB, tBC and tCD can be considered within the OT problem. Although we only present the evaluation of 2D STED images with constant labeling efficiencies across channels, our software package can directly be applied to multi-color 3D microscopy images with channel-specific labeling efficiencies.

Limitations: If the microscopy image shows especially dense point clouds, MultiMatch necessarily will have difficulties in differentiating between random and biological reasonable proximity. Note, however, that this is not a specific weakness of MultiMatch, but any other method will face this identifiability problem, which is caused by missing linkage information. It can only be overcome with additional prior information of the underlying biological sample. However, MultiMatch Mode I is especially robust against dense particle distribution in comparison to pairwise matching approaches as implemented in MultiMatch Mode II or greedy Nearest Neighbor Matchings. An adaption to tree like particle arrangements and the inclusion of additional constraints, e.g., incorporating regions of interest are future research objectives.

Methods

Optimal chain-matching

In the following we will denote the sets of two-dimensional particle coordinates in the image domain for each of the k color channels as1 X(1):=xl(1)l=1n1,…,X(k):=xl(k)l=1nk⊆R2,

where number of particles nj∈N≥0 for j ∈ {1, …, k}. For simplicity and related to the considered data in this article, we will only consider the cases k = 2, 3 in the following. Generalization to larger k is straight-forward. In a chain-like particle arrangement of the form x(1),…,x(k) with x(j) ∈ X(j), all neighbors x(j), x(j+1) have to be closer than the colocalization threshold t and we will denote according tuples as dtk-chains:

Definition 1

(dtk-chain). Fix k ≥ 2. For sets X(1), …, X(k), a distance d:R2×R2→R≥0 and a predefined maximal threshold t ≥ 0 a tuple of k points2 x(1),…,x(k)∈R2×kwithx(j)∈X(j)forj∈{1,…,k}

is a dtk-chain, if pairwise point distances along the fixed tuple point order are smaller or equal than t, i.e.,3 dx(j),x(j+1)≤tforj∈{1,…,k−1}.

In the context of our colocalization problem, d is the Euclidean distance (this can easily be generalized), a dt3-chain is a triplet and a dt2-chain a pair. For given t, we now aim to detect as many dtk-chains as possible:

Definition 2

(Optimal dtk-matching). A collection of pairwise disjoint dtk-chains is called dtk-matching. It is called optimal if its number of chains is maximal among all matchings.

Such an optimal dtk-matching can be found by utilizing a multi-marginal and unbalanced formulation of OT. For example, if k = 3, for each channel i = 1, 2, 3, we interpret coordinates of detected particles as support points with mass 1 of a respective discrete distribution. Due to this discrete structure, the resulting optimization problem will be finite-dimensional. Since in our measurements the number of detected particles per channel might differ, we require an unbalanced formulation to compare distributions with different total masses. A wide variety of penalty terms for mass discrepancies has been studied in the literature, see for instance64. Our problem formulation is closely related to an ℓ1-penalty for unmatched particles, see also52. We first consider the problem of finding optimal dt2-matchings between two point clouds, i.e. k = 2. This can be solved via the following optimization problem:

Definition 3

(Optimal dt2-matchings via unbalanced optimal transport). Let λ∈R≥0, set the cost function4 c:R2×R2→R≥0∪{∞},(x1,x2)↦d(x1,x2)−λifd(x1,x2)≤t,+∞otherwise,

and c∈Rn1×n2 the pairwise cost between all points in X(1) and X(2), defined by ci1,i2=c(xi1(1),xi2(2)). The optimal unbalanced transport problem of interest can now be stated as the following linear program5 argminπ∈Rn1×n2×n3∑i1=1n1∑i2=1n2ci1i2πi1i2s.t.∑i2=1n2πi1i2≤1foralli1=1,…,n1∑i1=1n1πi1i2≤1foralli2=1,…,n2πi1i2≥0forall(i1,i2)∈{1,…,n1}×{1,…,n2}.

Entries of an optimal π indicate which particles have been matched. The constraints enforces that each particle can at most be part of one matching, but it may also be discarded. By the definition of the cost vector c, the solution of Equation (5) does not match points x(1) and x(2) as soon as they are farther apart than t, but for each matching below distance t there is an incentive by the parameter λ. For λ sufficiently large in comparison to t one can show that the solution yields an optimal dt2-matching. Among all optimal matchings the above problem prefers one with the lowest sum of pairwise particle distances among matched particles.

We now generalize this to k = 3 via a multi-marginal transport problem.

Definition 4

(Optimal dt3-matchings via unbalanced multi-marginal optimal transport). Let λ∈R≥0, set the cost function6 c:R2×R2×R2→R≥0∪{∞},(x1,x2,x3)↦d(x1,x2)+d(x2,x3)−λifd(x1,x2)≤t∧d(x2,x3)≤t,+∞otherwise,

and let c∈Rn1×n2×n3, be the cost tensor between all triplets in (X(1), X(2), X(3)), defined by ci1i2i3=c(xi1(1),xi2(2),xi3(3)). Then the unbalanced multi-marginal OT problem can be stated as the following linear program:7 argminπ∈Rn1×n2×n3∑i1=1n1∑i2=1n2∑i3=1n3ci1i2i3πi1i2i3s.t.∑i2=1n2∑i3=1n3πi1i2i3≤1foralli1∈[n1]∑i1=1n1∑i3=1n3πi1i2i3≤1foralli2∈[n2]∑i1=1n1∑i3=1n3πi1i2i3≤1foralli3∈[n1]πi1i2i3≥0forall(i1,i2,i3)∈[n1]×[n2]×[n3],

where we used the notation [n] = {1, …, n}. As mentioned above, note that per the marginal constraints, particles may be matched at most once and can also be discarded. Likewise, by definition of the cost vector c only allows matchings between points that are valid dt3-chains. Analogously there is a matching incentive via the parameter λ and for sufficiently high values (relative to t) one can show that the above problem provides an optimal dt3-matching. Among all these matchings, one with minimal sum of pairwise distances is selected by the problem.

Generalization of Definition 4 to arbitrary k is now obvious, leading to a multi-marginal problem with k marginals. In general, multi-marginal problems quickly become numerically impractical due to the large number of variables. The cost function c in (6) has a chain structure, i.e. it can be written as a sum of functions only depending on (x1, x2) and (x2, x3). This chain structure allows the reformulation of the problem as a much more compact network flow problem (see Section below), and it implies the existence of optimal binary matchings. Problems where the cost exhibits a tree-structure can still be solved efficiently, see ref. 51 and references therein, but they cannot be formulated as network flow problems and do not exhibit binary minimizers in general.

Network flow formulation

In this section, we show that the multi-marginal optimal unbalanced transport problem corresponds to a min cost flow problem if the cost function has a chain structure as in (6). This has two relevant consequences:It guarantees that (7) has integer solutions and thus indeed corresponds to a matching problem, which in general does not hold true for discrete OT problems;

It allows us to solve the multi-marginal optimal unbalanced transport problem efficiently.

Definition 5

Let (V, E) be a directed graph with a source node S ∈ V, a target node T ∈ V, an edge capacity function lE:E→R∪∞ and an edge cost function cE:E→R∪∞. Then we call (V, E, cE, lE) a flow network. Given an amount of flow, m∈{R}_{+} the min cost flow problem consists in finding a function f:E→R that solves the following optimization problem:minf∑(u,v)∈Ef(u,v)cE(u,v)s.t.0≤f(u,v)≤lE(u,v)forall(u,v)∈E(capacityconstraints)∑{u:(u,v)∈E}f(u,v)−∑{w:(v,w)∈E}f(v,w)=0forallv≠S,T(flowconservation)∑{u:(S,u)∈E}f(S,u)−∑{v:(v,S)∈E}f(v,S)=m(flowsource)∑{u:(u,T)∈E}f(u,T)−∑{v:(T,v)∈E}f(T,v)=m(flowsink).

Notably, due to the total unimodularity of the constraint matrix, the min cost flow problem with integer total flow m and integer capacity function lE has an integer solution (Theorem 13.11 in65). In the following, we recast (7) to a min cost flow problem (see sketch in Supplementary Fig. 1):Node set V: Define source node S ∈ V and target node T ∈ V and add two nodes vl(j) and v^l(j) for each detected particle position xl(j) in Equation (1).

Edge set E:Connect nodes referring to the same detected point and set edge costs cE(vl(j),v^l(j))=−λk where k is the number of point clouds as in (1).

Add all possible edges of form (v^(j),v(j+1))∈E for j = 1, …, k − 1. Set edge costscE(v^(j),v(j+1))=∞,ifd(x(j),x(j+1))>td(x(j),x(j+1)),otherwise.

Include source and target nodes via edges of form (S,v(1)),(v^(k),T),(S,T)∈E, and set its costs to 0.

Define edge capacitieslE(vi,vj)=∞,ifvi=Sandvj=T1,otherwise.

Proposition 1

Let f:E→R be an integer solution of the min cost flow problem for the flow network (V, E, cE, lE) defined above with transported mass m=min(n1,n2,n3). Then one of the optimal solutions π* of the multi-marginal optimal unbalanced transport problem (7) is given by,8 πi1i2i3*=f(v^i1(1),vi2(2))f(v^i2(2),vi3(3)),

for i1 ∈ [n1], i2 ∈ [n2] and i3 ∈ [n3] with notation [n] = {1, …, n}.

Proof

First we show that π* as defined in Eq. (8) is in fact a valid transport plan for Eq. (7). For any i3 ∈ [n3] we have that, using the conservation constraint,∑i1=1n1∑i2=1n2πi1i2i3*=∑i2=1n2f(v^i2(2),vi3(3))∑i1=1n2f(v^i1(1),vi2(2))=∑i2=1n2f(v^i2(2),vi3(3))f(v^i2(2),vi2(2))≤∑i2=1n2f(v^i2(2),vi3(3))=f(v^i3(3),vi3(3))≤1

Analogously it is easy to verify that π* satisfies∑i1=1n1∑i3=1n3πi1i2i3*≤1foralli2∈[n2]∑i2=1n2∑i3=1n3πi1i2i3*≤1foralli1∈[n1].

Hence, π* is a feasible solution of (7). Further, since the source node S is directly connected to the target node T with an edge of infinite capacity and finite cost, the total flow cost must be finite. This implies that for any i1 ∈ [n1], i2 ∈ [n2] and i3 ∈ [n3], we have that f(v^i1(1),vi2(2))=0 if d(xi1(1),xi2(2))>t and f(v^i2(2),vi3(3))=0 if d(xi2(2),xi3(3))>t. Hence, using the shorthand notation⟨c,π⟩=∑i1=1n1∑i2=1n2∑i3=1n3ci1i2i3πi1i2i3,

we can rewrite the total cost of the transport problem as⟨c,π*⟩=∑i1=1n1∑i2=1n2∑i3=1n3d(xi1(1),xi2(2))+d(xi2(2),xi3(3))−λ⋅f(v^i1(1),vi2(2))f(v^i2(2),vi3(3)).

By the flow conservation constraints and the fact that f is an integer solution, we can simply reformulate the sum above in terms of the network flow cost function to obtain⟨c,π*⟩=∑(u,v)∈EcE(u,v)f(u,v).

Let us now assume that there exists a feasible solution of (7), π~, such that⟨c,π~⟩<⟨c,π*⟩.

Then we can define the flow f~:E→R by setting:f~(v^i1(1),vi2(2))=∑i3=1n3π~i1i2i3andf~(v^i2(2),vi3(3))=∑i1=1n1π~i1i2i3,

for i1 ∈ [n1], i2 ∈ [n2] and i3 ∈ [n3]. The value of the flow on the remaining nodes of E can then be determined by the conservation constraints. In particular, we have f~(S,T)=min{n1,n2,n3}−∑i1=1n1∑i2=1n2∑i3=1n3πi1i2i3. This flow is a feasible solution of the given min cost flow problem and hence, by the definition of the cost function for the edges we can derive a contradiction:∑(u,v)∈EcE(u,v)f~(u,v)=⟨c,π~⟩<∑(u,v)∈EcE(u,v)f(u,v).□

As a result of Proposition 1, we immediately obtain that the multi-marginal optimal unbalanced transport problem Eq. (7) has an integer solution and hence provides one-to-one point matchings.

Another significant consequence of Proposition 1 is that we can solve the unbalanced optimal transport problem given in Eq. (7) efficiently. While it is often unfeasible to compute directly the solution of the n1 ⋅ n2 ⋯ ⋅ nk-dimensional linear programming problem in Eq. (7), the min cost flow problem can be solved by the Scaling Minimum-Cost Flow Algorithm in ref. 66 in O(∣V∣2∣E∣log(∣V∣)) elementary operations, where ∣V∣ is the number of nodes, ∣E∣ is the number of edges. In our case the number of nodes is of the order O(n1 ⋅ ⋯ ⋅ nk) and the number of edges can be upper bounded by an expression of the order O(n1⋅n22⋅⋯⋅nk−12⋅nk). In practice, it is further possible to omit all edges with infinite cost, since the source S and the sink T are connected through an edge of cost 0 and with infinite capacity. This implies that for small t much fewer edges to the network are added which results in better computational performance.

For an image containing around 1,000 points in each color channel, a solution of the min cost flow problem can be computed for about 10 different values of t in ~1 s on a standard laptop.

Estimating the true chain-like particle abundances

The quality of fluorescence microscopy suffers from non-optimal labeling efficiencies and point detection errors. This will be addressed by a statistical framework to infer on how many of the detected structures in the image actually concur with the ground truth biological structure and how many detections represent only incomplete parts of the underlying particle assembly. For color channels i ∈ {1, …, k} let9 ξj(i)j=1ni⊂R2

be the pairs of coordinates of all particles that lie within the scope of the microscope. Note that these point clouds do not necessarily equal those defined in Eq. (1) describing the coordinates of detected particles, since we might not be able to measure all of the existing particles to do unsuccessful labeling or point detection errors.

Definition 6

(Labeling Efficiency). For each color channel i ∈ {1, …, k} we assume that there is a specific probability si ∈ (0, 1] quantifying whether a particle of this channel is successfully imaged and detected. For simplicity in the following we will always call probabilities silabeling efficiencies.

We further assume that the random event of successful detection is statistically independent for each point. Accordingly, the detection success can be described by independent Bernoulli variables10 Zj(i)j=1ni~Ber(si),

where si ∈ (0, 1] and ξj(i) is detectable, if and only if Zj(i)=1.

If there exists a true dtk-chain of form (ξ(1),…,ξ(k)), then this can only be correctly identified as such, if each of the included particles was detected, i.e., if and only if ∏i=1kZ(i)=1. From independence it follows that11 ∏i=1kZ(i)~Ber∏i=1ksi.

Detecting an ABC triplet correctly is Ber(sAsBsC) distributed. Therefore, all possible substructures that can be detected conditioned on the true underlying ABC triplet, i.e.,ABC triplet, if we see all particles

AB pair, if we do not see C

BC pair, if we do not see A

AC substructure, if we do not see B – which is detected as A and C singlets

A singlet, if we do not see B and C

B singlet, if we do not see A and C

C singlet, if we do not see A and B

∅, if we do not see A,B and C which can not be detected at all,

can accordingly be modeled as Multinomial random variable12 W⋅∣ABC=WABC∣ABCWAB∣ABCWBC∣ABCWAC∣ABCWA∣ABCWB∣ABCWC∣ABCW∅∣ABC.

This can be done in the same manner for all other structures of interest, i.e., true AB and BC pairs and A, B, and C singlets (and their respective substructures) yielding random variables W⋅∣AB, W⋅∣BC, W⋅∣A, W⋅∣B, W⋅∣C. The actual detectable numbers of those structures are13 WABC= ∑WABC∣⋅,WAB= ∑WAB∣⋅,WBC= ∑WBC∣⋅,WA= ∑WA∣⋅+∑WAC∣⋅,WB= ∑WB∣⋅,WC= ∑WC∣⋅+∑WAC∣⋅,

which define a random variable W=(WABC,WAB,WBC,WA,WB,WC)T. This leads to a statistical framework, that allows us to estimate the true underlying structures abundances from the detected number of structures.

Theorem 2

Let known, positive labeling efficiencies sA > 0, sB > 0 and sC > 0 and unknown structure abundances n=(nABC,nAB,nBC,nA,nB,nC)T and define N = ∑i∈{ABC, …, C}ni. Assume the multinomial model as described in Equation (12) and Equation (13).

Part 1: An unbiased estimator n^ of true abundances n is given as14 n^=1sAsBsC00000sC−1sAsBsC1sAsB0000sA−1sAsBsC01sBsC000sB−1sAsBsB−1sAsB01sA00(1−sA)(1−sC)sAsBsCsA−1sAsBsC−1sBsC01sB0sB−1sBsC0sB−1sBsC001sCW.

Part 2: For n → ∞ entrywise, nj/N → fj with ∞ > fj > 0 constant for each j ∈ {ABC, …, C}, and ΘΣ(n^)ΘT invertible,15 PΞ≤χ6,α2≤1−α,

where16 Ξ=(n^−n)T(Θμ)TΘΣ(n^)ΘT−1(Θμ)(n^−n)

and χ6,α2 is the α-quantile of a chi-squared distribution with 6 degrees of freedom and with Θ, μ and Σ(n^) defined as in the following proof. If ΘΣ(n^)ΘT−1 does not exist, we get Equation (15) with χr,α2 plugging its pseudoinverse ΘΣ(n^)ΘT+ in Equation (16), where r=rankΘΣ(n^)ΘT.

Proof

Part 1: conditioned on a true ABC triplet, the number of (mis)specifications resulting from incomplete labeling efficiencies is multinomially distributed:17 W⋅∣ABC=WABC∣ABCWAB∣ABCWBC∣ABCWAC∣ABCWA∣ABCWB∣ABCWC∣ABCW∅∣ABC~Mnom(nABC,pABC)

with probability vector18 pABC=sAsBsCsAsB(1−sC)(1−sA)sBsCsA(1−sB)sCsA(1−sB)(1−sC)(1−sA)sB(1−sC)(1−sA)(1−sB)sC(1−sA)(1−sB)(1−sC),

where ∑j=18pABC[j]=1. Accordingly, the abundances of (mis)detections of a true AB pair are19 WABC∣ABWAB∣ABWBC∣ABWAC∣ABWA∣ABWB∣ABWC∣ABW∅∣AB~Mnom(nAB,pAB)

with20 pAB=0sAsB00sA(1−sB)(1−sA)sB0(1−sA)(1−sB).

This can be done accordingly for all other structures of interest, i.e. BC pairs and A, B, and C singlets yielding21 W⋅∣ABC~Mnom(nABC,pABC)W⋅∣AB~Mnom(nAB,pAB)W⋅∣BC~Mnom(nBC,pBC)W⋅∣A~Mnom(nA,pA)W⋅∣B~Mnom(nB,pB)W⋅∣C~Mnom(nC,pC)

with22 pBC=00sBsC00sB(1−sC)(1−sB)sC(1−sB)(1−sC),pA=0000sA00(1−sA),pB=00000sB0(1−sB),pC=000000sC(1−sC).

Note, that ∅ can not be detected at all and substructure AC is counted as a separate A and C singlet Hence, the total numbers of detected triplets, pairs and singlets are defined as the following sums23 WABC=WABC∣ABCWAB=WAB∣ABC+WAB∣ABWBC=WBC∣ABC+WBC∣BCWA=WA∣ABC+WAC∣ABC+WA∣AB+WA∣AWB=WB∣ABC+WB∣AB+WB∣BC+WB∣BWC=WC∣ABC+WAC∣ABC+WC∣BC+WC∣C.

This can be rewritten as24 W=WABCWABWBCWAWBWC=ΘW⋅∣ABC+W⋅∣AB+W⋅∣BC+W⋅∣A+W⋅∣B+W⋅∣C,

using the transformation matrix25 Θ=100000000100000000100000000110000000010000010010∈R6×8.

With this definition of Θ we delete the last entry in each binomial distributed vector and add an AC substructure appearance to singlet detections A and B. By Eq. (24) we get that26 E[W]=Θμn

with27 μ=pABCpABpBCpApBpC∈R8×6.

Hence, with positive labeling efficiencies sA > 0, sB > 0 and sC > 0, multiplying28 (Θμ)−1=1sAsBsC00000sC−1sAsBsC1sAsB0000sA−1sAsBsC01sBsC000sB−1sAsBsB−1sAsB01sA00(1−sA)(1−sC)sAsBsCsA−1sAsBsC−1sBsC01sB0sB−1sBsC0sB−1sBsC001sC

with W introduces an unbiased estimator n^.

Part 2: we utilize that by the central limit theorem for a multinomially distributed random variable M ~ Mnom(m, p) with probability vector p=(p1,p2,…,pk)T29 1mM−mp→DNk0k,diag(p)−ppTform→∞,

where30 diag(p)=p10⋯0p2⋯⋮⋮pk

and 0k=(0,...,0)T∈Rk (see, e.g., ref. 67). Hence, for n entrywise large enough, we can approximate properly scaled independent, multinomial random vectors31 W⋅∣ABC,W⋅∣AB,W⋅∣BC,W⋅∣A,W⋅∣B,W⋅∣C

with multi-dimensional normal distributions, respectively. In the following assume n → ∞ entrywise and nj/N → fj with ∞ > fj > 0 constant for each j ∈ {ABC, …, C}, where N = ∑i∈{ABC, …, C}ni. Then, it holds that32 ∑i∈{ABC,…,C}niN1niW⋅∣i−nipi=1N∑i∈{ABC,…,C}W⋅∣i−nipi→DN808,∑i∈{ABC,…,C}fidiag(pi)−pipiT.

For now, suppose ∑nidiag(pi)−pipiT is invertible. Then in the limit33 ∑fidiag(pi)−pipiT−1/21N ∑W⋅∣i−nipi=∑Nfidiag(pi)−pipiT−1/2 ∑W⋅∣i−nipi=∑nidiag(pi)−pipiT−1/2 ∑W⋅∣i−nipi

and hence34 ∑nidiag(pi)−pipiT−1/2 ∑W⋅∣i−nipi→DN808,I8×8,

where I8×8 is the 8-dimensional identity matrix. In the following we denote35 Σ(n)=∑nidiag(pi)−pipiT.

Multiplying (Θμ)−1Θ with Eq. (32) consequently yields36 (Θμ)−1ΘΣ(n)ΘT(Θμ)−1T−1/2(Θμ)−1Θ ∑W⋅∣i−(Θμ)−1Θ ∑nipi=(Θμ)−1ΘΣ(n)ΘT(Θμ)−1T−1/2n^−n→DN606,I6×6

with n^=(Θμ)−1Θ∑W⋅∣i and n = (Θμ)−1Θμn = (Θμ)−1Θ∑nipi. By law of large numbers, it holds that37 1Nn^−n=n^N−nN→P06.

and hence for all j ∈ {ABC, …, C}38 n^jN→Pfj.

By Slutsky’s Lemma we can use Eq. (38) to replace n in Σ(n) with n^. For n → ∞ entrywise this yields39 Ξ=(n^−n)T(Θμ)TΘΣ(n^)ΘT−1(Θμ)(n^−n)→Dχ62.

In case ΘΣ(n^)ΘT is not invertible, one can use its pseudoinverse yielding convergence to a chi-square distribution with r degrees of freedom, i.e., χr2 in Equation (39), where r=rankΘΣ(n^)ΘT□.

With Part 2 of Theorem 2 we can construct a confidence ellipsoid around n^ in a straight-forward manner. To show that Ξ in our setting is approximately chi-square distributed for finite sample sizes and to compare simulated and theoretical coverages of n^, we performed a simulation study as described in the following section.

Simulation study setup

In the first simulation study a predefined number of triplets, pairs, and singlets are generated as follows:Step 1: Draw the coordinate for channel B as b~U([0,400⋅r]2), where U is the continuous uniform distribution.

Step 2a: Draw angle α~U[0,2π] and normally distributed distance dA~N(t,0.5). Set a=bcos(α)dA+sin(α)dA.

Step 2b: Draw ϵ~N(0,0.2) and set angle β = α + π + ϵ. Draw dC~N(t,0.5) and set c=bcos(β)dC+sin(β)dC.

Step 3: Round a, b and c to match the pixel grid [0,400]2⊆N≥02.

This design favors to simulate triplets of an approximately linear structure. Pairs are simulated by skipping either Step 2a or 2b. Singlets are drawn as in Step 1.

In the second simulation study quadruples, triplets, pairs, and singlets n are generated similarly, but replacing and addingStep 2b: Draw angle β~U[0,2π] and dC~N(t,0.5) and set d=bcos(β)dC+sin(β)dC.

Step 2c: Draw angle γ~U[0,2π] and dD~N(t,0.5) and set d=ccos(γ)dD+sin(γ)dD.

This simulation setup allows arbitrarily curved chain-structures. The distance threshold is always fixed to t = 70 nm.

To obtain intensity images close to an experimental STED setup from the simulated point sets we followed the simulation setup introduced in Tameling et al.17, to mimic experimental STED images of 400  × 400 pixels with full-width at half-maximum (FWHM) value of 40 nm (approximately the resolution of the STED microscope) and pixel size 25 nm = 1 pixel). In the second simulation study (including quadruples) the Poisson noise level was on average increased by a factor of 10.

Methods included in the simulation study

For the Ripley’s K based Statistical Object Distance Analysis (SODA)21 we used the triplet colocalization protocol SODA 3 Colors in ICY68 (version 2.4.0.0). For the analysis we used default input parameters and set scale threshold per channel to be 100. The plugin BlobProb15 was called in ImageJ/Fiji69 (version 2.3.0/1.53q) and the number of colocalized blobs were considered. We set voxel size to 25 nm in every dimension and the threshold per channel to 100. The ConditionalColoc18 from GitHub (https://github.com/kjaqaman/ConditionalColoc) was executed on MATLAB (version R2023a). Particles were detected using the “point-source detection” algorithm provided via the integrated u-track package (https://github.com/DanuserLab/u-track).

For all implementations but ConditionalColoc the detected chain-structure abundances were output as integers. Therefore, we scaled abundances, i.e., divided them by the total number of particles detected in channel B. ConditionColoc already aims to output probabilities that are scaled by detected particles per channel, hence no further transformation of the output was performed by us. Since for all simulated Scenarios the same number of particles was generated in every channel, we ensured that both scaling procedures are comparable. The maximal colocalization threshold is set to t = 5 pixels = 125 nm throughout all considered methods.

Nanoruler samples

Custom-made DNA nanoruler samples featuring one, two, or three fluorophore spots, each consisting of 20 fluorophores (Alexa Fluor488, Alexa Fluor594, Star Red), with a distance between the spots of 70 nm, were purchased from Gattaquant - DNA Nanotechnologies (Gräfelfing, Germany). The biotinylated nanorulers were immobilized on a BSA-biotin-neutravidin surface according to the manufacturer’s specifications.

Stimulated emission depletion super-resolution light microscopy

Image acquisition was done using a quad scanning STED microscope (Abberior Instruments, Göttingen, Germany) equipped with a UPlanSApo 100x/1,40 Oil objective (Olympus, Tokyo, Japan). Excitation of Alexa Fluor 488, Alexa Fluor 594 and Star Red was achieved by laser beams featuring wave lengths of 485 nm, 561 nm, and 640 nm, respectively. For STED imaging, a laser beam with an emission wavelength of 775 nm was applied. For all experimental STED images, a pixel size of 25 nm was utilized. For visualization purposes, contrast stretching and increasement of image brightness was applied to exemplary STED images within the figures of this manuscript. No image processing was applied prior to the application of the MultiMatch analysis workflow.

Statistics and reproducibility

The statistical framework developed and applied in this manuscript and the settings of simulation studies performed are presented in the Method sections. All sample sizes and significance levels of the confidence bands are listed in the respective figure legends. Experimental and simulated data and analysis scripts to reproduce results and figures are provided on Zenodo (10.5281/zenodo.7221879)70.

Supplementary information

Supplementary Information

Supplementary information

The online version contains supplementary material available at 10.1038/s42003-024-06772-8.

Acknowledgements

We would like to thank Leo Lehmann for his help in software implementation, Jan-Niklas Dohrke for the graphic illustration of the nanorulers in Fig. 1C, and Christiane Elgert and Arndt von Haeseler for constructive criticism on the manuscript. J.N. is supported by the Austrian Science Fund (FWF) project number F78 to Arndt von Haeseler. This work was supported by the European Research Council (ERC AdG no. 835102) (to S.J.). G.N., H.L., S.J., and A.M. are supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy–EXC 2067/1-390729940 and H.L., B.S., S.J., and A.M. by the DFG CRC 1456 “Mathematics of Experiment” (Project Number B04, C06).

Author contributions

S.S. and S.J. designed the experimental work and S.S. acquired the experimental data. J.N. designed the MultiMatch algorithm with support from H.L., B.S., and A.M. G.N. implemented the MultiMatch Python package. J.N. and G.N. performed data analysis and method evaluation and developed the statistical methodology with input from H.L. and A.M.; A.M. initiated and S.J. and A.M. supervised the project. J.N. wrote the manuscript with contributions from all co-authors. All authors read and approved the final manuscript.

Peer review

Peer review information

Communications Biology thanks Thibault Lagache and the other, anonymous, reviewer for their contribution to the peer review of this work. Primary Handling Editor: Ophelia Bu.

Funding

Open Access funding enabled and organized by Projekt DEAL.

Data availability

Datasets generated and analyzed in this manuscript can be accessed via Zenodo (10.5281/zenodo.7221879)70.

Code availability

The Python package MultiMatch is available on GitHub repository https://github.com/gnies/multi_match. All scripts used to create the main and Supplementary Figs. are implemented in R (version 4.1.0) and Python (version 3.8.5) and are available via Zenodo (10.5281/zenodo.7221879)70. In order to locate the positions of the particles in STED images, we perform point detection via the Python package scikit-image71 (version 0.19.1). This is provided as an optional analysis step in our MultiMatch implementation for the evaluation of intensity matrices. Multi-color microscopy images, point detection results and MultiMatch output can be loaded into the interactive napari viewer. MultiMatch is compatible with Python package napari39 (version 0.4.18) and an exemplary use-case is described on our repository https://github.com/gnies/multi_match. We utilize the minimum-cost flow solver provided in the package ortools72 (version 9.4.1874).

Competing interests

The authors declare no competing interests.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Cainero I Measuring nanoscale distances by structured illumination microscopy and image cross-correlation spectroscopy (SIM-ICCS) Sensors 2021 21 2010 10.3390/s21062010 33809144
Cainero, I. et al. Measuring nanoscale distances by structured illumination microscopy and image cross-correlation spectroscopy (SIM-ICCS). Sensors 21, 2010 (2021).33809144 10.3390/s21062010
2. Costa R Morphological study of TNPO3 and SRSF1 interaction during myogenesis by combining confocal, structured illumination and electron microscopy analysis Mol. Cell. Biochem. 2021 476 1797 1811 10.1007/s11010-020-04023-y 33452620
Costa, R. et al. Morphological study of TNPO3 and SRSF1 interaction during myogenesis by combining confocal, structured illumination and electron microscopy analysis. Mol. Cell. Biochem. 476, 1797–1811 (2021).33452620 10.1007/s11010-020-04023-y
3. Shimizu Y Cargo sorting zones in the trans-Golgi network visualized by super-resolution confocal live imaging microscopy in plants Nat. Commun. 2021 12 1901 10.1038/s41467-021-22267-0 33772008
Shimizu, Y. et al. Cargo sorting zones in the trans-Golgi network visualized by super-resolution confocal live imaging microscopy in plants. Nat. Commun. 12, 1901 (2021).33772008 10.1038/s41467-021-22267-0
4. Sahl SJ Hell SW Jakobs S Fluorescence nanoscopy in cell biology Nat. Rev. Mol. Cell Biol. 2017 18 685 701 10.1038/nrm.2017.71 28875992
Sahl, S. J., Hell, S. W. & Jakobs, S. Fluorescence nanoscopy in cell biology. Nat. Rev. Mol. Cell Biol. 18, 685–701 (2017).28875992 10.1038/nrm.2017.71
5. Betzig E Imaging intracellular fluorescent proteins at nanometer resolution Science 2006 313 1642 1645 10.1126/science.1127344 16902090
Betzig, E. et al. Imaging intracellular fluorescent proteins at nanometer resolution. Science 313, 1642–1645 (2006).16902090 10.1126/science.1127344
6. Hess ST Girirajan TPK Mason MD Ultra-high resolution imaging by fluorescence photoactivation localization microscopy Biophys. J. 2006 91 4258 4272 10.1529/biophysj.106.091116 16980368
Hess, S. T., Girirajan, T. P. K. & Mason, M. D. Ultra-high resolution imaging by fluorescence photoactivation localization microscopy. Biophys. J. 91, 4258–4272 (2006).16980368 10.1529/biophysj.106.091116
7. Rust MJ Bates M Zhuang X Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM) Nat. Methods 2006 3 793 796 10.1038/nmeth929 16896339
Rust, M. J., Bates, M. & Zhuang, X. Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM). Nat. Methods 3, 793–796 (2006).16896339 10.1038/nmeth929
8. Hell SW Wichmann J Breaking the diffraction resolution limit by stimulated emission: Stimulated-emission-depletion fluorescence microscopy Opt. Lett. 1994 19 780 782 10.1364/OL.19.000780 19844443
Hell, S. W. & Wichmann, J. Breaking the diffraction resolution limit by stimulated emission: Stimulated-emission-depletion fluorescence microscopy. Opt. Lett. 19, 780–782 (1994).19844443 10.1364/OL.19.000780
9. Hell SW Far-field optical nanoscopy Science 2007 316 1153 1158 10.1126/science.1137395 17525330
Hell, S. W. Far-field optical nanoscopy. Science 316, 1153–1158 (2007).17525330 10.1126/science.1137395
10. Klar TA Jakobs S Dyba M Egner A Hell SW Fluorescence microscopy with diffraction resolution barrier broken by stimulated emission Proc. Natl Acad. Sci. USA 2000 97 8206 8210 10.1073/pnas.97.15.8206 10899992
Klar, T. A., Jakobs, S., Dyba, M., Egner, A. & Hell, S. W. Fluorescence microscopy with diffraction resolution barrier broken by stimulated emission. Proc. Natl Acad. Sci. USA 97, 8206–8210 (2000).10899992 10.1073/pnas.97.15.8206
11. Malkusch S Coordinate-based colocalization analysis of single-molecule localization microscopy data Histochem. Cell Biol. 2012 137 1 10 10.1007/s00418-011-0880-5 22086768
Malkusch, S. et al. Coordinate-based colocalization analysis of single-molecule localization microscopy data. Histochem. Cell Biol. 137, 1–10 (2012).22086768 10.1007/s00418-011-0880-5
12. Manders EMM Verbeek FJ Aten JA Measurement of co-localization of objects in dual-colour confocal images J. Microsc. 1993 169 375 382 10.1111/j.1365-2818.1993.tb03313.x 33930978
Manders, E. M. M., Verbeek, F. J. & Aten, J. A. Measurement of co-localization of objects in dual-colour confocal images. J. Microsc. 169, 375–382 (1993).33930978 10.1111/j.1365-2818.1993.tb03313.x
13. Xu L Resolution, target density and labeling effects in colocalization studies - suppression of false positives by nanoscopy and modified algorithms FEBS J. 2016 283 882 898 10.1111/febs.13652 26756570
Xu, L. et al. Resolution, target density and labeling effects in colocalization studies - suppression of false positives by nanoscopy and modified algorithms. FEBS J. 283, 882–898 (2016).26756570 10.1111/febs.13652
14. Adler J Parmryd I Quantifying colocalization by correlation: the Pearson correlation coefficient is superior to the Mander’s overlap coefficient Cytom. Part A 2010 77A 733 742 10.1002/cyto.a.20896
Adler, J. & Parmryd, I. Quantifying colocalization by correlation: the Pearson correlation coefficient is superior to the Mander’s overlap coefficient. Cytom. Part A 77A, 733–742 (2010).10.1002/cyto.a.20896
15. Fletcher PA Scriven DRL Schulson MN Moore EDW Multi-image colocalization and its statistical significance Biophys. J. 2010 99 1996 2005 10.1016/j.bpj.2010.07.006 20858446
Fletcher, P. A., Scriven, D. R. L., Schulson, M. N. & Moore, E. D. W. Multi-image colocalization and its statistical significance. Biophys. J. 99, 1996–2005 (2010).20858446 10.1016/j.bpj.2010.07.006
16. Wang S Spatially adaptive colocalization analysis in dual-color fluorescence microscopy IEEE Trans. Image Process. 2019 28 4471 4485 10.1109/TIP.2019.2909194
Wang, S. et al. Spatially adaptive colocalization analysis in dual-color fluorescence microscopy. IEEE Trans. Image Process. 28, 4471–4485 (2019).10.1109/TIP.2019.2909194
17. Tameling C Colocalization for super-resolution microscopy via optimal transport Nat. Comput. Sci. 2021 1 199 211 10.1038/s43588-021-00050-x 35874932
Tameling, C. et al. Colocalization for super-resolution microscopy via optimal transport. Nat. Comput. Sci. 1, 199–211 (2021).35874932 10.1038/s43588-021-00050-x
18. Vega-Lugo J da Rocha-Azevedo B Dasgupta A Jaqaman K Analysis of conditional colocalization relationships and hierarchies in three-color microscopy images J. Cell Biol. 2022 221 e202106129 10.1083/jcb.202106129 35552363
Vega-Lugo, J., da Rocha-Azevedo, B., Dasgupta, A. & Jaqaman, K. Analysis of conditional colocalization relationships and hierarchies in three-color microscopy images. J. Cell Biol. 221, e202106129 (2022).35552363 10.1083/jcb.202106129
19. Ripley BD The second-order analysis of stationary point processes J. Appl. Probab. 1976 13 255 266 10.2307/3212829
Ripley, B. D. The second-order analysis of stationary point processes. J. Appl. Probab. 13, 255–266 (1976).10.2307/3212829
20. Mukherjee S Gonzalez-Gomez C Danglot L Lagache T Olivo-Marin J-C Generalizing the statistical analysis of objects’ spatial coupling in bioimaging IEEE Signal Process. Lett. 2020 27 1085 1089 10.1109/LSP.2020.3003821
Mukherjee, S., Gonzalez-Gomez, C., Danglot, L., Lagache, T. & Olivo-Marin, J.-C. Generalizing the statistical analysis of objects’ spatial coupling in bioimaging. IEEE Signal Process. Lett. 27, 1085–1089 (2020).10.1109/LSP.2020.3003821
21. Lagache T Mapping molecular assemblies with fluorescence microscopy and object-based spatial statistics Nat. Commun. 2018 9 698 10.1038/s41467-018-03053-x 29449608
Lagache, T. et al. Mapping molecular assemblies with fluorescence microscopy and object-based spatial statistics. Nat. Commun. 9, 698 (2018).29449608 10.1038/s41467-018-03053-x
22. Winter FR Multicolour nanoscopy of fixed and living cells with a single STED beam and hyperspectral detection Sci. Rep. 2017 7 46492 10.1038/srep46492 28417977
Winter, F. R. et al. Multicolour nanoscopy of fixed and living cells with a single STED beam and hyperspectral detection. Sci. Rep. 7, 46492 (2017).28417977 10.1038/srep46492
23. Spahn C Grimm JB Lavis LD Lampe M Heilemann M Whole-cell, 3D, and multicolor STED imaging with exchangeable fluorophores Nano Lett. 2019 19 500 505 10.1021/acs.nanolett.8b04385 30525682
Spahn, C., Grimm, J. B., Lavis, L. D., Lampe, M. & Heilemann, M. Whole-cell, 3D, and multicolor STED imaging with exchangeable fluorophores. Nano Lett. 19, 500–505 (2019).30525682 10.1021/acs.nanolett.8b04385
24. Butkevich AN Photoactivatable fluorescent dyes with hydrophilic caging groups and their use in multicolor nanoscopy J. Am. Chem. Soc. 2021 143 18388 18393 10.1021/jacs.1c09999 34714070
Butkevich, A. N. et al. Photoactivatable fluorescent dyes with hydrophilic caging groups and their use in multicolor nanoscopy. J. Am. Chem. Soc. 143, 18388–18393 (2021).34714070 10.1021/jacs.1c09999
25. Glogger M Synergizing exchangeable fluorophore labels for multitarget STED microscopy ACS Nano 2022 16 17991 17997 10.1021/acsnano.2c07212 36223885
Glogger, M. et al. Synergizing exchangeable fluorophore labels for multitarget STED microscopy. ACS Nano 16, 17991–17997 (2022).36223885 10.1021/acsnano.2c07212
26. Gonzalez Pisfil M Stimulated emission depletion microscopy with a single depletion laser using five fluorochromes and fluorescence lifetime phasor separation Sci. Rep. 2022 12 14027 10.1038/s41598-022-17825-5 35982114
Gonzalez Pisfil, M. et al. Stimulated emission depletion microscopy with a single depletion laser using five fluorochromes and fluorescence lifetime phasor separation. Sci. Rep. 12, 14027 (2022).35982114 10.1038/s41598-022-17825-5
27. Wang J Fan Y Sanger JM Sanger JW STED analysis reveals the organization of nonmuscle muscle II, muscle myosin II, and F-actin in nascent myofibrils Cytoskeleton 2022 79 122 132 10.1002/cm.21729 36125330
Wang, J., Fan, Y., Sanger, J. M. & Sanger, J. W. STED analysis reveals the organization of nonmuscle muscle II, muscle myosin II, and F-actin in nascent myofibrils. Cytoskeleton 79, 122–132 (2022).36125330 10.1002/cm.21729
28. Saal KA Heat denaturation enables multicolor X10-STED microscopy Sci. Rep. 2023 13 5366 10.1038/s41598-023-32524-5 37005431
Saal, K. A. et al. Heat denaturation enables multicolor X10-STED microscopy. Sci. Rep. 13, 5366 (2023).37005431 10.1038/s41598-023-32524-5
29. Andronov L Genthial R Hentsch D Klaholz BP SplitSMLM, a spectral demixing method for high-precision multi-color localization microscopy applied to nuclear pore complexes Commun. Biol. 2022 5 1 13 10.1038/s42003-022-04040-1 34987157
Andronov, L., Genthial, R., Hentsch, D. & Klaholz, B. P. SplitSMLM, a spectral demixing method for high-precision multi-color localization microscopy applied to nuclear pore complexes. Commun. Biol. 5, 1–13 (2022).34987157 10.1038/s42003-022-04040-1
30. Unterauer EM Spatial proteomics in neurons at single-protein resolution Cell 2024 187 1785–1800.e16 10.1016/j.cell.2024.02.045 38552614
Unterauer, E. M. et al. Spatial proteomics in neurons at single-protein resolution. Cell 187, 1785–1800.e16 (2024).38552614 10.1016/j.cell.2024.02.045
31. Beater S Holzmeister P Lalkens B Tinnefeld P Simple and aberration-free 4color-STED - multiplexing by transient binding Opt. Express 2015 23 8630 8638 10.1364/OE.23.008630 25968701
Beater, S., Holzmeister, P., Lalkens, B. & Tinnefeld, P. Simple and aberration-free 4color-STED - multiplexing by transient binding. Opt. Express 23, 8630–8638 (2015).25968701 10.1364/OE.23.008630
32. Willig KI Rizzoli SO Westphal V Jahn R Hell SW STED microscopy reveals that synaptotagmin remains clustered after synaptic vesicle exocytosis Nature 2006 440 935 939 10.1038/nature04592 16612384
Willig, K. I., Rizzoli, S. O., Westphal, V., Jahn, R. & Hell, S. W. STED microscopy reveals that synaptotagmin remains clustered after synaptic vesicle exocytosis. Nature 440, 935–939 (2006).16612384 10.1038/nature04592
33. Reinhardt SCM Ångström-resolution fluorescence microscopy Nature 2023 617 711 716 10.1038/s41586-023-05925-9 37225882
Reinhardt, S. C. M. et al. Ångström-resolution fluorescence microscopy. Nature 617, 711–716 (2023).37225882 10.1038/s41586-023-05925-9
34. Smallcombe A Multicolor imaging: The important question of co-localization BioTechniques 2001 30 1240 1246 10.2144/01306bt01 11414212
Smallcombe, A. Multicolor imaging: The important question of co-localization. BioTechniques 30, 1240–1246 (2001).11414212 10.2144/01306bt01
35. Sastre, D., Estadella, I., Bosch, M. & Felipe, A. Methods in Molecular Biology (Springer, 2019).
36. Goucher DR Wincovitch SM Garfield SH Carbone KM Malik TH A quantitative determination of multi-protein interactions by the analysis of confocal images using a pixel-by-pixel assessment algorithm Bioinformatics 2005 21 3248 3254 10.1093/bioinformatics/bti531 15947019
Goucher, D. R., Wincovitch, S. M., Garfield, S. H., Carbone, K. M. & Malik, T. H. A quantitative determination of multi-protein interactions by the analysis of confocal images using a pixel-by-pixel assessment algorithm. Bioinformatics 21, 3248–3254 (2005).15947019 10.1093/bioinformatics/bti531
37. Humpert F Yahiatène I Lummer M Sauer M Huser T Quantifying molecular colocalization in live cell fluorescence microscopy J. Biophoton. 2015 8 124 132 10.1002/jbio.201300146
Humpert, F., Yahiatène, I., Lummer, M., Sauer, M. & Huser, T. Quantifying molecular colocalization in live cell fluorescence microscopy. J. Biophoton. 8, 124–132 (2015).10.1002/jbio.201300146
38. Haas KT Peaucelle A Protocol for multicolor three-dimensional dSTORM data analysis using MATLAB-based script package Grafeo STAR Protoc. 2021 2 100808 10.1016/j.xpro.2021.100808 34541556
Haas, K. T. & Peaucelle, A. Protocol for multicolor three-dimensional dSTORM data analysis using MATLAB-based script package Grafeo. STAR Protoc. 2, 100808 (2021).34541556 10.1016/j.xpro.2021.100808
39. napari contributers. Napari: A Multi-dimensional Image Viewer For Python10.5281/zenodo.8115575 (2019).
40. Boettiger AN Super-resolution imaging reveals distinct chromatin folding for different epigenetic states Nature 2016 529 418 422 10.1038/nature16496 26760202
Boettiger, A. N. et al. Super-resolution imaging reveals distinct chromatin folding for different epigenetic states. Nature 529, 418–422 (2016).26760202 10.1038/nature16496
41. Miron E Chromatin arranges in chains of mesoscale domains with nanoscale functional topography independent of cohesin Sci. Adv. 2020 6 eaba8811 10.1126/sciadv.aba8811 32967822
Miron, E. et al. Chromatin arranges in chains of mesoscale domains with nanoscale functional topography independent of cohesin. Sci. Adv. 6, eaba8811 (2020).32967822 10.1126/sciadv.aba8811
42. Villani, C. Optimal Transport (Springer Berlin, 2009).
43. Panaretos VM Zemel Y Statistical aspects of Wasserstein distances Annu. Rev. Stat. Appl. 2019 6 405 431 10.1146/annurev-statistics-030718-104938
Panaretos, V. M. & Zemel, Y. Statistical aspects of Wasserstein distances. Annu. Rev. Stat. Appl. 6, 405–431 (2019).10.1146/annurev-statistics-030718-104938
44. Peyré, G. & Cuturi, M. Computational Optimal Transport: With Applications to Data Science. Foundations and Trends® in Machine Learning 11, 355–607 (2019).
45. Zaritsky A Decoupling global biases and local interactions between cell biological variables eLife 2017 6 e22323 10.7554/eLife.22323 28287393
Zaritsky, A. et al. Decoupling global biases and local interactions between cell biological variables. eLife 6, e22323 (2017).28287393 10.7554/eLife.22323
46. Kim Y-H Pass B A general condition for monge solutions in the multi-marginal optimal transport problem SIAM J. Math. Anal. 2014 46 1538 1550 10.1137/130930443
Kim, Y.-H. & Pass, B. A general condition for monge solutions in the multi-marginal optimal transport problem. SIAM J. Math. Anal. 46, 1538–1550 (2014).10.1137/130930443
47. Pass B Multi-marginal optimal transport: theory and applications ESAIM: Math. Model. Numer. Anal. 2015 49 1771 1790 10.1051/m2an/2015020
Pass, B. Multi-marginal optimal transport: theory and applications. ESAIM: Math. Model. Numer. Anal. 49, 1771–1790 (2015).10.1051/m2an/2015020
48. Chizat L Peyré G Schmitzer B Vialard F-X Unbalanced optimal transport: Dynamic and Kantorovich formulations J. Funct. Anal. 2018 274 3090 3123 10.1016/j.jfa.2018.03.008
Chizat, L., Peyré, G., Schmitzer, B. & Vialard, F.-X. Unbalanced optimal transport: Dynamic and Kantorovich formulations. J. Funct. Anal. 274, 3090–3123 (2018).10.1016/j.jfa.2018.03.008
49. Friesecke G Matthes D Schmitzer B Barycenters for the Hellinger–Kantorovich Distance Over $\mathbb{R}^d$ SIAM Journal on Mathematical Analysis 2021 53 62 110 10.1137/20M1315555
Friesecke, G., Matthes, D. & Schmitzer, B. Barycenters for the Hellinger–Kantorovich Distance Over $\mathbb{R}^d$. SIAM Journal on Mathematical Analysis 53, 62–110 (2021).10.1137/20M1315555
50. Heinemann F Klatt M Munk A Kantorovich-Rubinstein distance and barycenter for finitely supported measures: foundations and algorithms Appl. Math. Optim. 2022 87 4 10.1007/s00245-022-09911-x
Heinemann, F., Klatt, M. & Munk, A. Kantorovich-Rubinstein distance and barycenter for finitely supported measures: foundations and algorithms. Appl. Math. Optim. 87, 4 (2022).10.1007/s00245-022-09911-x
51. Beier, F., von Lindheim, J., Neumayer, S. & Steidl, G. Unbalanced multi-marginal optimal transport. J. Math. Imag. Vis. 65, 394–413 (2023).
52. Le, K., Nguyen, H., Nguyen, K., Pham, T. & Ho, N. On multimarginal partial optimal transport: Equivalent forms and computational complexity. Proceedings of The 25th International Conference on Artificial Intelligence and Statistics 4397–4413 (2022).
53. Schulter, S., Vernaza, P., Choi, W. & Chandraker, M. Deep network flow for multi-object tracking. 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR) 2730–2739 (2017).
54. Chari, V., Lacoste-Julien, S., Laptev, I. & Sivic, J. On pairwise costs for network flow multi-object tracking. 2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR) 5537–5545 (2015).
55. Jaqaman K Robust single-particle tracking in live-cell time-lapse sequences Nat. Methods 2008 5 695 702 10.1038/nmeth.1237 18641657
Jaqaman, K. et al. Robust single-particle tracking in live-cell time-lapse sequences. Nat. Methods 5, 695–702 (2008).18641657 10.1038/nmeth.1237
56. Zhang, L., Li, Y. & Nevatia, R. Global data association for multi-object tracking using network flows. 2008 IEEE Conference on Computer Vision and Pattern Recognition 1–8 (2008).
57. Lin T Ho N Cuturi M Jordan MI On the complexity of approximating multimarginal optimal transport J. Mach. Learn. Res. 2022 23 1 43
Lin, T., Ho, N., Cuturi, M. & Jordan, M. I. On the complexity of approximating multimarginal optimal transport. J. Mach. Learn. Res. 23, 1–43 (2022).
58. Hummert J Tashev SA Herten D-P An update on molecular counting in fluorescence microscopy Int. J. Biochem. Cell Biol. 2021 135 105978 10.1016/j.biocel.2021.105978 33865985
Hummert, J., Tashev, S. A. & Herten, D.-P. An update on molecular counting in fluorescence microscopy. Int. J. Biochem. Cell Biol. 135, 105978 (2021).33865985 10.1016/j.biocel.2021.105978
59. Schmied JJ DNA origami-based standards for quantitative fluorescence microscopy Nat. Protoc. 2014 9 1367 1391 10.1038/nprot.2014.079 24833175
Schmied, J. J. et al. DNA origami-based standards for quantitative fluorescence microscopy. Nat. Protoc. 9, 1367–1391 (2014).24833175 10.1038/nprot.2014.079
60. Schmied JJ Fluorescence and super-resolution standards based on DNA origami Nat. Methods 2012 9 1133 1134 10.1038/nmeth.2254 23223165
Schmied, J. J. et al. Fluorescence and super-resolution standards based on DNA origami. Nat. Methods 9, 1133–1134 (2012).23223165 10.1038/nmeth.2254
61. Rothemund PWK Folding DNA to create nanoscale shapes and patterns Nature 2006 440 297 302 10.1038/nature04586 16541064
Rothemund, P. W. K. Folding DNA to create nanoscale shapes and patterns. Nature 440, 297–302 (2006).16541064 10.1038/nature04586
62. Balzarotti F Nanometer resolution imaging and tracking of fluorescent molecules with minimal photon fluxes Science 2017 355 606 612 10.1126/science.aak9913 28008086
Balzarotti, F. et al. Nanometer resolution imaging and tracking of fluorescent molecules with minimal photon fluxes. Science 355, 606–612 (2017).28008086 10.1126/science.aak9913
63. Gwosch KC MINFLUX nanoscopy delivers 3D multicolor nanometer resolution in cells Nat. Methods 2020 17 217 224 10.1038/s41592-019-0688-0 31932776
Gwosch, K. C. et al. MINFLUX nanoscopy delivers 3D multicolor nanometer resolution in cells. Nat. Methods 17, 217–224 (2020).31932776 10.1038/s41592-019-0688-0
64. Liero M Mielke A Savaré G Optimal transport in competition with reaction: the Hellinger–Kantorovich distance and geodesic curves SIAM J. Math. Anal. 2016 48 2869 2911 10.1137/15M1041420
Liero, M., Mielke, A. & Savaré, G. Optimal transport in competition with reaction: the Hellinger–Kantorovich distance and geodesic curves. SIAM J. Math. Anal. 48, 2869–2911 (2016).10.1137/15M1041420
65. Alexander, S. Combinatorial optimization: Polyhedra and efficiency, 24 edn (Springer, 2003).
66. Goldberg AV An efficient implementation of a scaling minimum-cost flow algorithm J. Algorithms 1997 22 1 29 10.1006/jagm.1995.0805
Goldberg, A. V. An efficient implementation of a scaling minimum-cost flow algorithm. J. Algorithms 22, 1–29 (1997).10.1006/jagm.1995.0805
67. Morris C Central limit theorems for multinomial sums Ann. Stat. 1975 3 165 188 10.1214/aos/1176343006
Morris, C. Central limit theorems for multinomial sums. Ann. Stat. 3, 165–188 (1975).10.1214/aos/1176343006
68. de Chaumont F Icy: an open bioimage informatics platform for extended reproducible research Nat. Methods 2012 9 690 696 10.1038/nmeth.2075 22743774
de Chaumont, F. et al. Icy: an open bioimage informatics platform for extended reproducible research. Nat. Methods 9, 690–696 (2012).22743774 10.1038/nmeth.2075
69. Schindelin J Fiji: an open-source platform for biological-image analysis Nat. Methods 2012 9 676 682 10.1038/nmeth.2019 22743772
Schindelin, J. et al. Fiji: an open-source platform for biological-image analysis. Nat. Methods 9, 676–682 (2012).22743772 10.1038/nmeth.2019
70. Naas, J. et al. Source Data and Scripts - MultiMatch: Geometry-Informed Colocalization in Multi-Color Super-Resolution Microscopy (v0.0.2). Zenodo10.5281/zenodo.7221879 (2024).
71. Walt Svd Scikit-image: Image processing in Python PeerJ 2014 2 e453 10.7717/peerj.453 25024921
Walt, Svd et al. Scikit-image: Image processing in Python. PeerJ 2, e453 (2014).25024921 10.7717/peerj.453
72. Perron, L. & Furnon, V. OR-tools https://developers.google.com/optimization/ (2022).
