Preparing HMEC1 cells for fluorescence imaging
HMEC1 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 phosphatebuffered saline (PBS) for 10 min at room temperature. Cells were then preblocked with 3% bovine serum albumin (BSA) in PBS, incubated with antiClathrin primary antibody (mouse monoclonal X22, Catalog # ab2731, Abcam, diluted 1000fold in 1% BSA/PBS) for 25 min at room temperature, and treated with goat serum diluted 50fold. Bound primary antibody was detected by treatment with Alexa 555labeled antimouse IgG (Catalog # A21424, Thermo Fisher Scientific, diluted 750fold 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 HMEC1 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 (XCite 110LED, Excelitas Technologies) filtered through a standard Cy3 filterset (Cy34040CZHE M327122 Brightline, Semrock). Signal from the sample was also filtered through this filterset before being acquired by the camera.
Preparing BSC1 cells for superresolution imaging
BSC1 cells (CCL26, American Type Culture Collection) were plated on custom, glassbottom dishes (Catalog # PG35G10CNON, MatTek Corporation) fitted with highperformance Zeiss coverglasses (Catalog # 474030900000, 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 X100 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 200fold in 1% BSA/PBS) for 30 min at room temperature, and treated with goat serum (Catalog # G6767, Sigma, diluted 50fold in 1% BSA/PBS). Bound primary antibody was counterstained with Alexa 647labeled goat antimouse IgG (Catalog # A21236, Thermo Fisher Scientific, diluted 750fold in 1% BSA/PBS) for 30 min at room temperature.
Optical setup for superresolution imaging
Images of BSC1 cells were acquired using a standard inverted epifluorescence microscope (Observer A1, Zeiss) fitted with a 63× (1.4 NA) PlanApochromat oilimmersion 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 (Di01R405/488/561/63525×36, Semrock), and focused onto the back focal plane of the objective lens. Signal from the sample was filtered using a bandpass filter (FF01676/2925, Semrock) and acquired by an electron multiplying charge coupled device (EMCCD) camera (iXon DU897BV, 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 customwritten software in the C programming language.
Superresolution image acquisition
Buffer solutions in which the BSC1 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 maximumpower setting of the 635 nm laser to switch off most of the Alexa 647 fluorophores and achieve sparse singlemolecule activation conditions. Acquisition was initiated once steadystate, singlemolecule 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 p_{X}(x) = p_{Y}(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 N_{c} pairs of those locations are spaced at distances between r_{min} and r_{max} 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 r_{i} between (d_i^{mathrm{c}}) and (d_i^{{mathrm{c}}’ }) is uniformly distributed between r_{min} and r_{max}, and θ_{i} is uniformly distributed between 0 and 2π. The subset Δ^{u} is an additional completely spatially random distribution of N_{u} = N − 2N_{c} locations. For the simulated images analyzed to obtain the results in Fig. 3, the following values were used: N = 2500, N_{c} = 250, r_{min} = 2990 nm, and r_{max} = 3010 nm.
Generating locations avoiding specific spacings (inhibition)
A set of locations, Δ ≔ {d_{1},…,d_{N}}, in which no two locations are spaced between r_{min} and r_{max} of each other is generated as follows. For i = 1,…,N, the ith location d_{i} is drawn from a completely spatially random distribution. The ith point is not added as a location if r_{min} ≤ d_{ij} ≤ r_{max} for some 1 ≤ j ≤ i, where d_{ij} 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, r_{min} = 2990 nm and r_{max} = 3010 nm.
Simulating an image of clathrincoated pits
When simulating an image of clathrincoated pits, the detector is modeled as a set of pixels {C_{1},…,C_{K}}. The photon count detected in the kth pixel is modeled as ({cal T}_k: = S_k + B_k), where S_{k} and B_{k} are both Poisson random variables^{21}. The total photons detected at the kth pixel from all clathrincoated pits within the region represented by the image is denoted by S_{k}. The background photon count B_{k} has a mean of B = 100 photons per pixel.
The mean of S_{k} 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 clathrincoated pits in the image, N_{d} denotes the total number of photons detected from the dth clathrincoated pit, C_{k} denotes the area of the kth pixel, and f_{d} denotes the photon distribution profile for the dth clathrincoated pit. For the simulated image of clathrincoated pits in Fig. 1, D = 419 and N_{d} had values uniformly distributed between 500 to 2000 photons.
For each clathrincoated pit, f_{d} 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 clathrincoated 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 {C_{1},…,C_{k}}. 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, C_{k} denotes the area of the kth pixel, f_{d} 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, f_{d} 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 J_{1} denotes the firstorder 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πN_{a}/λ, where N_{a} 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: N_{a} = 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 n_{F} = ⌈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 clathrincoated pits using waveletfiltering^{22} 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 (GlobalThresholding): Detects molecules or clathrincoated 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 waveletfiltering
The image was filtered using the product of two consecutive wavelet transforms as described in ref. ^{22}. Each isolated set of one or more edgeconnected pixels obtained from the filtering was identified as a region of the image containing an object, i.e., an individual molecule or clathrincoated 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 globalthresholding
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 edgeconnected pixels obtained from the thresholding was identified as a region of the image containing an object, i.e., an individual molecule or clathrincoated 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 clathrincoated pit or single molecule was localized by fitting a Gaussian or Airy profile, respectively, to the corresponding 5 × 5 pixel image identified using either waveletfiltering or globalthresholding 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 clathrincoated 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 clathrincoated 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 clathrincoated 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 leastsquares estimator to obtain the final location estimates. The location parameters (x_{0}, y_{0}), width parameter (σ when fitting Gaussian profiles and κ when fitting Airy profiles), and the total photon count were estimated for each clathrincoated pit or single molecule.
Estimating Ripley’s Kfunction from one image
When estimating L(r) − r for a set of localizations obtained by analyzing an image of either clathrincoated pits or single molecules, (L(r) = sqrt {K(r)/pi }) for r > 0. K(r) denotes the Ripley’s Kfunction, 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 S^{2} denotes the area in the object space corresponding to the image being analyzed, w_{ij} denotes the Ripley’s isotropic edge correction weights^{27}, D denotes the total number of localizations, and d_{ij} 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 Kfunction 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 D^{m} denotes the number of localizations obtained from the mth image and (hat K_m(r)) denotes estimates of the Ripley’s Kfunction calculated using the localization obtained from that image.
Estimating paircorrelations for an image analysis approach
Estimates of the paircorrelation results for an image analysis approach, denoted as a, were calculated by a weighted averaging of paircorrelation 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. Paircorrelations 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 weightedaverage paircorrelation 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 paircorrelation 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).
Calculating resolution adjusted Ripley’s Kfunction
For the algorithmic resolution limit α for a specific image analysis approach determined as described above, L_{2α}(r) − r is calculated as
$$L_{2alpha }(r) – r = sqrt {frac{{K_{2alpha }(r) + 4pi alpha ^2}}{pi }} – r,$$
where K_{2α}(r) = K(r) − K(2α). See Supplementary Note 9 for details regarding the determination of 2α for each image analysis approach from the corresponding paircorrelation results.
Analyzing experimental superresolution data
Let ({cal D}_i) denote the ith dataset consisting of M_{i} 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′ = m ⋅ i for m = 1,…,M_{i}. Here, M_{i} denotes the number of images in dataset ({cal D}_i) and is calculated as M_{i} = ⌊M_{1}/i⌋. Each of these datasets were analyzed independently using the image analysis approaches described above. The experimental dataset analyzed for Fig. 4 consisted of M_{1} = 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 clathrincoated 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 waveletfiltering or globalthresholding followed by fitting with either Airy or Gaussian profiles was performed using custom programs developed with the MIATool software framework^{29} in Java. The ThunderSTORM^{23} and SimpleFit^{24} software packages were used with default settings for the various options within the software. The QuickPALM^{25} software was used with the Full Width at Half Maximum = 2 setting to match the width of the single molecule or clathrincoated pit being localized. Calculations for L(r) − r, paircorrelations, and L_{2α}(r) − r were performed using customdeveloped scripts in the MATLAB programming environment (The MathWorks, Inc., Natick, MA). All figures were similarly prepared using MATLAB.