Resolution limit of image analysis algorithms

Bioinformatics

Preparing HMEC-1 cells for fluorescence imaging

HMEC-1 cells were fixed using 1.7% (w/v) paraformaldehyde (Electron Microscopy Sciences) at room temperature and permeabilized by incubation with 0.02% (w/v) saponin in phosphate-buffered saline (PBS) for 10 min at room temperature. Cells were then preblocked with 3% bovine serum albumin (BSA) in PBS, incubated with anti-Clathrin primary antibody (mouse monoclonal X22, Catalog # ab2731, Abcam, diluted 1000-fold in 1% BSA/PBS) for 25 min at room temperature, and treated with goat serum diluted 50-fold. Bound primary antibody was detected by treatment with Alexa 555-labeled anti-mouse IgG (Catalog # A-21424, Thermo Fisher Scientific, diluted 750-fold in 1% BSA/PBS) for 25 min at room temperature. Cells were washed twice with PBS between each incubation and finally immersed in 1.5 ml of 1% BSA/PBS prior to imaging.

Fluorescence microscopy imaging

Fixed HMEC-1 cells were imaged with a Zeiss (Axiovert 200M) inverted epifluorescence microscope fitted with a 63× (1.4 NA) Plan Apo objective (Carl Zeiss) using a CCD camera (Orca ER, Hamamatsu). The sample was illuminated using a broadband LED illumination (X-Cite 110LED, Excelitas Technologies) filtered through a standard Cy3 filterset (Cy3-4040C-ZHE M327122 Brightline, Semrock). Signal from the sample was also filtered through this filterset before being acquired by the camera.

Preparing BS-C-1 cells for superresolution imaging

BS-C-1 cells (CCL-26, American Type Culture Collection) were plated on custom, glass-bottom dishes (Catalog # PG35G-10C-NON, MatTek Corporation) fitted with high-performance Zeiss coverglasses (Catalog # 474030-9000-00, Thickness 1.5) at a density of 20,000 cells and incubated for 16–24 h. The coverglasses were cleaned by sonicating them sequentially in 50% HPLC grade ethanol, 1 mM HCL in 50% HPLC grade ethanol, and 1 M KOH in 50% HPLC grade ethanol for 20 min each with extensive washing using MilliQ water after every sonication. The cells were fixed using 3.4% paraformaldehyde (Electron Microscopy Sciences) in PBS for 10 min at 37 °C. The fixed cells were then permeabilized using 0.1% Triton X-100 in PBS for 10 min at room temperature. Cells were then preblocked with 5% BSA (Catalog # BP1600, Fisher Scientific), incubated with mouse α-tubulin antibody (Catalog # A11126, Thermo Fisher Scientific, diluted 200-fold in 1% BSA/PBS) for 30 min at room temperature, and treated with goat serum (Catalog # G6767, Sigma, diluted 50-fold in 1% BSA/PBS). Bound primary antibody was counter-stained with Alexa 647-labeled goat anti-mouse IgG (Catalog # A-21236, Thermo Fisher Scientific, diluted 750-fold in 1% BSA/PBS) for 30 min at room temperature.

Optical setup for superresolution imaging

Images of BS-C-1 cells were acquired using a standard inverted epifluorescence microscope (Observer A1, Zeiss) fitted with a 63× (1.4 NA) Plan-Apochromat oil-immersion objective (Carl Zeiss). The sample was illuminated using 635 nm and 450 nm diode lasers (OptoEngine) for the excitation of Alexa 647 and for photoswitching respectively. The illumination light was directed towards the microscope using custom laser excitation optics, reflected into the sample using a dichroic filter (Di01-R405/488/561/635-25×36, Semrock), and focused onto the back focal plane of the objective lens. Signal from the sample was filtered using a bandpass filter (FF01-676/29-25, Semrock) and acquired by an electron multiplying charge coupled device (EMCCD) camera (iXon DU897-BV, Andor) set to conventional readout mode. Axial drift during the course of the acquisition was corrected in real time using a custom focus stabilization system comprising an 850 nm diode laser (PI, USA), a quadrant position detector (Thorlabs), and an XYZ piezo positioner (PI, USA). All devices including lasers, shutters, and cameras were controlled and synchronized using custom-written software in the C programming language.

Superresolution image acquisition

Buffer solutions in which the BS-C-1 cell samples were prepared were replaced with an imaging buffer immediately prior to the imaging of each sample. This imaging buffer consisted of 50 mM cysteamine, 0.5 mg ml−1 glucose oxidase, 40 μg ml−1 catalase, and 10% glucose (w/v) in PBS, pH 7.4. (All buffer ingredients were purchased from Thermo Fisher Scientific.) The imaging buffer was freshly prepared for each acquisition. A coverslip was placed over the sample and sealed using vacuum grease. For each acquisition, the sample was briefly illuminated (approximately 5–10 min) using the maximum-power setting of the 635 nm laser to switch off most of the Alexa 647 fluorophores and achieve sparse single-molecule activation conditions. Acquisition was initiated once steady-state, single-molecule blinking was observed with periodic illumination using the 405 nm laser to reactivate fluorophores and maintain sufficient activation density. The activation laser power was kept low throughout the acquisition. For each dataset, typically 50,000 images of stochastically activated fluorophores were acquired with 50 ms exposure time for each image and without EM gain resulting in an acquisition rate of 20 frames per second.

Generating completely spatially random locations

For a completely spatially random distribution of locations in an S μm × S μm region, the (x, y) coordinate for each location was obtained by drawing realizations of independent random variables X and Y, each uniformly distributed with probability density functions pX(x) = pY(y) = 1/S, 0 ≤ x, y ≤ S. For the simulated image in Fig. 2, S = 50 resulting in a 50 μm × 50 μm region.

Generating locations for deterministic structures

Deterministic structures consist of D molecules located at evenly spaced points on the circumference of a circle of radius r. The location of the dth point ((x_0^d,y_0^d)) is given by (x_0^d = r,{mathrm{cos}}left( {2pi d/D} right)) and (y_0^d = r,{mathrm{sin}},left( {2pi d/D} right)), where d = 1,…,D.

Generating locations with a preferred spacing (clustering)

A set of N random locations, denoted by Δ, such that Nc pairs of those locations are spaced at distances between rmin and rmax is generated by combining three subsets of locations so that Δ = ΔcΔc’Δu. The subset Δ({varDelta }^{mathrm{c}}: = { d_1^{mathrm c}, ldots ,d_{N_{mathrm{c}}}^{mathrm{c}}}) consists of completely spatially random locations (d_i^{mathrm{c}}: = (x_i^{mathrm{c}},y_i^{mathrm{c}})). The subset Δ({varDelta }^{{mathrm{c}}’ }: = { d_1^{{mathrm{c}}’ }, ldots ,d_{N_{mathrm{c}}}^{{mathrm{c}}’ }}) consists of locations corresponding to Δc, where each location (d_i^{{mathrm{c}}’ }: = (x_i^{{mathrm{c}}’ },y_i^{{mathrm{c}}’ })), is calculated as

$$begin{array}{*{20}{c}} {x_i^{{mathrm{c}}’ } = x_i^{mathrm{c}} + r_i,{mathrm{cos}},theta _i,} \ {y_i^{{mathrm{c}}’ } = y_i^{mathrm{c}} + r_i,{mathrm{sin}},theta _i.} end{array}$$

The distance ri between (d_i^{mathrm{c}}) and (d_i^{{mathrm{c}}’ }) is uniformly distributed between rmin and rmax, and θi is uniformly distributed between 0 and 2π. The subset Δu is an additional completely spatially random distribution of Nu = N − 2Nc locations. For the simulated images analyzed to obtain the results in Fig. 3, the following values were used: N = 2500, Nc = 250, rmin = 2990 nm, and rmax = 3010 nm.

Generating locations avoiding specific spacings (inhibition)

A set of locations, Δ {d1,…,dN}, in which no two locations are spaced between rmin and rmax of each other is generated as follows. For i = 1,…,N, the ith location di is drawn from a completely spatially random distribution. The ith point is not added as a location if rmin ≤ dij ≤ rmax for some 1 ≤ j ≤ i, where dij denotes the distance between the ith and jth locations. For the simulated images analyzed to obtain the results in Fig. 3, the following values were used: N = 2500, rmin = 2990 nm and rmax = 3010 nm.

Simulating an image of clathrin-coated pits

When simulating an image of clathrin-coated pits, the detector is modeled as a set of pixels {C1,…,CK}. The photon count detected in the kth pixel is modeled as ({cal T}_k: = S_k + B_k), where Sk and Bk are both Poisson random variables21. The total photons detected at the kth pixel from all clathrin-coated pits within the region represented by the image is denoted by Sk. The background photon count Bk has a mean of B = 100 photons per pixel.

The mean of Sk is given by

$$mu _k: = mathop {sum}limits_{d = 1}^D left( {N_d{int}_{C_k} f_d(r){mathrm d}r} right),$$

where D is the total number of clathrin-coated pits in the image, Nd denotes the total number of photons detected from the dth clathrin-coated pit, Ck denotes the area of the kth pixel, and fd denotes the photon distribution profile for the dth clathrin-coated pit. For the simulated image of clathrin-coated pits in Fig. 1, D = 419 and Nd had values uniformly distributed between 500 to 2000 photons.

For each clathrin-coated pit, fd is modeled as a Gaussian profile given by

$$f_d(r): = frac{1}{{2pi M^2sigma ^2}} cdot exp{left( { – frac{{left( {x – Mx_0^d} right)^2}}{{2(Msigma )^2}} – frac{{left( {y – My_0^d} right)^2}}{{2(Msigma )^2}}} right)},$$

where M denotes magnification, σ denotes the width of the Gaussian profile, and ((x_0^d,y_0^d)) denotes the center of the dth clathrin-coated pit. The coordinate ((x_0^d,y_0^d)) is drawn from a completely spatially randomly distributed set of D locations generated as described above. For the simulated image in Fig. 1, M = 63 and σ = 120 nm.

Simulating images of single molecules

When simulating an image of single molecules, the detector is again modeled as a set of pixels {C1,…,Ck}. The photon count detected at the kth pixel is modeled as a Poisson random variable with mean given by,

$$mu _k: = mathop {sum}limits_{d,=,1}^D N{int}_{C_k} f_d(r){mathrm d}r + B,$$

where N denotes the total number of photons detected from the molecule, Ck denotes the area of the kth pixel, fd denotes the photon distribution profile for the dth molecule, and B denotes the uniform number of background photons detected in each pixel. For each molecule, fd is modeled as an Airy profile given by

$$f_d(r): = frac{{left[ {J_1left( {frac{kappa }{M}||r – r_0^d||} right)} right]^2}}{{pi ||r – r_0^d||^2}},$$

where J1 denotes the first-order Bessel function of the first kind, (||r – r_0^d|| = left((x – Mx_0^d)^2 + (y – My_0^d)^2right)^{1/2}), ((x_0^d,y_0^d)) denotes the location of the dth molecule, and M denotes the magnification of the optical system. κ is calculated as κ = 2πNa/λ, where Na denotes the numerical aperture, and λ denotes the wavelength of the photons emitted by the molecule. The following values were used for simulating all images of single molecules: Na = 1.3, λ = 525 nm, M = 100, and N = 1000 photons. Images that were analyzed for all figures were simulated using a background of B = 0 photons per pixel.

Simulating PALM images of tubulin

When simulating images of fluorescently labeled tubulin molecules that are stochastically photoactivated and detected, each image is taken as an image of single molecules and is simulated as described above. The positions of molecules within each image are obtained by taking a spatial point pattern consisting of N points and uniformly distributing them among nF = q−1 images, where q is the probability that a molecule is visible in any particular frame. For all simulated tubulin datasets, the spatial point pattern provided by Sage et al.10 was used.

List of image analysis approaches

The following is a list of the image analysis approaches that were used:

• Algorithm 1 (Wavelet): Detects molecules or clathrin-coated pits using wavelet-filtering22 and estimates their locations by fitting Airy profiles to the detected molecules or Gaussian profiles to the detected pits. Further details are provided below.

• Algorithm 2 Detects and localizes molecules using the default settings of the software package taken from ref. 23.

• Algorithm 3 Detects and localizes molecules using the default settings of the software package taken from ref. 24.

• Algorithm 4 (Global-Thresholding): Detects molecules or clathrin-coated pits by identifying pixels above a threshold value and estimates their locations by fitting Airy profiles to the detected molecules or Gaussian profiles to the detected pits. Further details are provided below.

• Algorithm 5 Detects and localizes molecules using the default settings of the software package taken from ref. 25.

• Algorithm 6 Detects and localizes molecules using the default settings of the multiemitter version of the software package taken from ref. 23.

For Algorithms 2, 3, 5, and 6, the software packages provided by the authors of the algorithms were used. We used the default settings of those software packages and did not attempt to optimize them for use here. The purpose of this paper is not to investigate the different algorithms. We used several algorithms solely to illustrate that different algorithms can have different algorithmic resolutions limits. It is quite likely that an expert could obtain very different results from those shown here. We would therefore like to emphasize that we are not interested in an evaluation of these algorithms per se. To arrive at such an evaluation in a rigorous manner would require a very different approach. For example, the involvement of the authors of the software packages, as is routinely done in software evaluation challenges, is the preferred approach. It is under these provisos that we somewhat reluctantly mention the references to the software packages that were used.

Identifying objects by wavelet-filtering

The image was filtered using the product of two consecutive wavelet transforms as described in ref. 22. Each isolated set of one or more edge-connected pixels obtained from the filtering was identified as a region of the image containing an object, i.e., an individual molecule or clathrin-coated pit. For the subsequent localization of that object, a 5 × 5 pixel region centered on the average pixel coordinate of the corresponding set of identified pixels was used.

Identifying objects by global-thresholding

The image was thresholded using 25% of the maximum pixel intensity in the dataset as the threshold value. Each isolated set of one or more edge-connected pixels obtained from the thresholding was identified as a region of the image containing an object, i.e., an individual molecule or clathrin-coated pit. For the subsequent localization of that object, a 5 × 5 pixel region centered on the average pixel coordinate of the corresponding set of identified pixels was used.

Localizing objects detected by wavelet/threshold approaches

Each clathrin-coated pit or single molecule was localized by fitting a Gaussian or Airy profile, respectively, to the corresponding 5 × 5 pixel image identified using either wavelet-filtering or global-thresholding as described above. An initial location estimate and an initial value for the σ parameter denoting the width of a Gaussian profile was calculated for each clathrin-coated pit or single molecule by applying the approach described in ref. 26 to the corresponding image. An initial value for the κ parameter of the Airy profile was calculated as 1.323/σ. The background associated with each clathrin-coated pit or single molecule was taken as the median of the intensities in the edge pixels of the corresponding 5 × 5 pixel image. An initial estimate of the photon count detected from each molecule or clathrin-coated pit was taken as the sum of the pixel intensities in the corresponding image after subtracting the background. Airy or Gaussian profiles with initial values for the various parameters calculated as described above were fitted to each 5 × 5 pixel image using a least-squares estimator to obtain the final location estimates. The location parameters (x0, y0), width parameter (σ when fitting Gaussian profiles and κ when fitting Airy profiles), and the total photon count were estimated for each clathrin-coated pit or single molecule.

Estimating Ripley’s K-function from one image

When estimating L(r) − r for a set of localizations obtained by analyzing an image of either clathrin-coated pits or single molecules, (L(r) = sqrt {K(r)/pi }) for r > 0. K(r) denotes the Ripley’s K-function, defined as

$$K(r): = lambda ^{ – 1}E{ {mathrm{number}},{mathrm{of}},{mathrm{events}},{mathrm{within}},{mathrm{a}},{mathrm{distance}},r,{mathrm{of}},{mathrm{an}},{mathrm{arbitrary}},{mathrm{event}}} .$$

The estimator for K(r) is given by

$$hat K(r) = frac{{S^2}}{{D(D – 1)}}mathop {sum}limits_{i = 1}^D mathop {sum}limits_{j = 1}^D w_{ij}I(0 < d_{ij} < r),$$

where S2 denotes the area in the object space corresponding to the image being analyzed, wij denotes the Ripley’s isotropic edge correction weights27, D denotes the total number of localizations, and dij denotes the distance between the ith and jth localizations. The indicator function is defined as

$$I(0 < d_{ij} < r): = left{ {begin{array}{*{20}{c}} {1} & {{mathrm{if}},0 < ,di_j < 1} \ {0} & {{mathrm{otherwise}}} end{array}}. right.$$

Estimating Ripley’s K-function using multiple images

When estimating L(r) − r for a particular image analysis approach from a total of D localizations distributed among M images,

$$hat L(r) – r = mathop {sum}limits_{m = 1}^M left( {frac{{D^m}}{D}} right)hat K_m(r) – r,$$

where Dm denotes the number of localizations obtained from the mth image and (hat K_m(r)) denotes estimates of the Ripley’s K-function calculated using the localization obtained from that image.

Estimating pair-correlations for an image analysis approach

Estimates of the pair-correlation results for an image analysis approach, denoted as a, were calculated by a weighted averaging of pair-correlation estimates from multiple simulated images as follows. A total of M images containing D single molecules were simulated. A set of localizations of size (D_a^m) were obtained by applying analysis approach a to the mth image, for m = 1,…,M. Pair-correlations estimates (hat g_a^m(r)) were calculated independently for each set of (D_a^m) localizations using a MATLAB implementation of the approach in ref. 28. The weighted-average pair-correlation estimates for each analysis approach a is then calculated as

$$hat g_a(r) = mathop {sum}limits_{m = 1}^M left( {frac{{D_a^m}}{{D_a}}} right)hat g_a^m(r),$$

where (D_a = D_a^1 + D_a^2 + cdots + D_a^M). The pair-correlation results shown in Fig. 2 were calculated using M = 2000 images containing D = 250,000 molecules.

Determining the algorithmic resolution limit

The radius of correlation for analysis approach a is defined as

$$rho : = inf_{r> 0}left{ {left. {r:g_a(rprime) = 1,{mathrm{for}},{mathrm{all}},rprime , in ,left[ {r,infty} right.} right)} right},$$

when a analyzes completely spatially random data. To estimate ρ from (hat g_a(r)), the scheme presented in Supplementary Note 9 was implemented. For a set ({cal R} = { r_1, ldots ,r_m}) of finely, regularly spaced proposal values of ρ,

$$hat rho = operatorname*{argmin}_ {r_i in {cal R}} T(r_i)$$

where

$$T(r_i) = (m – i + 1)^{ – 1}(m – i)^{ – 1}mathop {sum}limits_{j = i}^m (hat g_a(r_j) – bar g_i)^2$$

and (bar g_i = (m – i + 1)^{ – 1}mathop {sum}nolimits_{j = i}^m hat g_a(r_j)). Algorithmic resolution limit α is then estimated as (hat rho /2).

For the algorithmic resolution limit α for a specific image analysis approach determined as described above, L2α(r) − r is calculated as

$$L_{2alpha }(r) – r = sqrt {frac{{K_{2alpha }(r) + 4pi alpha ^2}}{pi }} – r,$$

where K2α(r) = K(r) − K(2α). See Supplementary Note 9 for details regarding the determination of 2α for each image analysis approach from the corresponding pair-correlation results.

Analyzing experimental superresolution data

Let ({cal D}_i) denote the ith dataset consisting of Mi images ({ {cal I}_1^i, ldots ,{cal I}_{M_i}^i}) for i = 1,…,N. Let ({cal D}_1) denote the set of acquired experimental images. The remaining datasets ({cal D}_2, ldots ,{cal D}_N) were generated by summing images from ({cal D}_1) as follows:

$${cal I}_m^i = mathop {sum}limits_{j = m_0^prime}^{m^prime } {cal I}_j^1,$$

where (m_0^prime = m(i – 1)) + 1, m′ = mi for m = 1,…,Mi. Here, Mi denotes the number of images in dataset ({cal D}_i) and is calculated as Mi = M1/i. Each of these datasets were analyzed independently using the image analysis approaches described above. The experimental dataset analyzed for Fig. 4 consisted of M1 = 50,000 images.

Reconstructing superresolution images

Superresolution images were reconstructed from a set of localizations by simulating Gaussian profiles (as described above in the simulation of clathrin-coated pits) centered at each localization. The superresolution image was generated using a pixel size of 1/63 μm in the object space which corresponds to 1/16th the object space pixel size of the acquired images. Each Gaussian profile was simulated using σ = 17 nm for the width parameter and 100 detected photons. The value for σ corresponds to the limit of the localization accuracy calculated using the average number of photons detected from each localized fluorophore and the average background noise associated with each localization.

Software

Region of interest identification using wavelet-filtering or global-thresholding followed by fitting with either Airy or Gaussian profiles was performed using custom programs developed with the MIATool software framework29 in Java. The ThunderSTORM23 and SimpleFit24 software packages were used with default settings for the various options within the software. The QuickPALM25 software was used with the Full Width at Half Maximum  = 2 setting to match the width of the single molecule or clathrin-coated pit being localized. Calculations for L(r) − r, pair-correlations, and L2α(r) − r were performed using custom-developed scripts in the MATLAB programming environment (The MathWorks, Inc., Natick, MA). All figures were similarly prepared using MATLAB.