Extracting nanoscale membrane morphology from single-molecule localizations

Membrane surface reconstruction at the nanometer scale is required for understanding mechanisms of subcellular shape change. This historically has been the domain of electron microscopy, but extraction of surfaces from specific labels is a difficult task in this imaging modality. Existing methods for extracting surfaces from fluorescence microscopy have poor resolution or require high-quality super-resolution data that are manually cleaned and curated. Here, we present NanoWrap, a new method for extracting surfaces from generalized single-molecule localization microscopy data. This makes it possible to study the shape of specifically labeled membranous structures inside cells. We validate NanoWrap using simulations and demonstrate its reconstruction capabilities on single-molecule localization microscopy data of the endoplasmic reticulum and mitochondria. NanoWrap is implemented in the open-source Python Microscopy Environment.


S.1 Point cloud simulation
In order to demonstrate the e ectiveness of the algorithm, it was necessary to generate biologically-motivated test structures.To do this, we established a constructive solid geometry library where each primitive is represented as a signed distance function (SDF).An SDF describes an object by taking a point location as input and returning the signed distance between the point and the surface of the object.An example SDF is that of a sphere of radius ' 2 R, 3 ( ?) = || ?|| ', where ? 2 R 3 .3 ( ?) = 0 indicates that point ? is on the surface of the sphere.Points located within the sphere have a negative distance to the sphere surface and points outside of the sphere have a positive distance to the sphere surface.
Points were generated from a test structure using Monte-Carlo sampling of an octree of the structure's SDF.The octree subdivides space that straddles where the SDF goes to zero until a fixed sampling rate 3G is reached.The center of each octree box with an absolute SDF value of less than 3G is kept with a probability of ? and rejected otherwise.This results in a point set of size (.
Next, non-specific localizations are simulated.For a given noise fraction =, = =( 1 = additional points are generated with uniform randomness over the bounding box of the point set.This ensures /) = = for ) = + ( total points. To simulate localization precision, each point 8 in the resulting point set is jittered by sampling values from a Gaussian with mean `8 equivalent to the location of point 8 2 0 . . ) 1 and standard deviation 2 8,4 for 4 2 G, H, I. Up to # samples are generated per point.The total point set is then clipped to size ) via random uniform selection of points.
The uncertainty 2 8,4 is calculated by simulating an exponential distribution with a mean value of d, corresponding to the average photon count of a simulated localization.This distribution is clipped to only contain values greater than a background photon count 1, simulating the noise floor for data collection.For the remaining counts, the first ( are selected and where A 4 is the resolution of the simulated system along axis 4.

S.2 NanoWrap parameters
To generate the parameter space for evaluation of SPR meshes derived from simulated data in Figure 2, we followed recommendations in (18).For SPR we estimated normals from 10, 30, 50 and 100 nearest neighbors, set the smoothing parameter U = 0, 1, 2 and 4, and set the samples-per-node to 1.5, 5, 10, 20 and 30.We used an octree depth of 8, 8 Gauss-Seidel relaxations, and a scale factor of 1.1.

S.3 Parameter sensitivity
The search of the parameter space used to quantify RMS surface error (Figure 2) also permits an analysis of the sensitivity of that error to parameter choices.As the starting surface is an isosurface on density, the optimal choice of the threshold parameter for this is expected to correlate to labeling density.In practice, however, as long as the starting estimate is not too far o , the subsequent optimisation steps will ensure good surface quality.Because the number of localizations pulling on a mesh face also correlates with density, the curvature weight required to balance the point attraction force is also expected to increase with increasing density.Figure S6 shows this is the case, although the minimum in RMS error is fairly broad indicating that the algorithm is not overly sensitive to the exact curvature weight.Increasing the background at a constant density shifts the optimal curvature weight to slightly higher values, but is not as significant as a change in density.The underlying curvature of the structure itself has only a very slight e ect on the optimal curvature weight (Figure S7) the localization errors are already considered in the algorithm we do not expect a strong change in the optimal curvature weight with changing localization precision.
Unlike the RMS surface error, which can tolerate quite large variations in the starting surface, estimating the correct genus (number of holes) depends strongly on the threshold used to generate the starting surface.This is consistent with other shrink-wrapping approaches, which rely on manual thresholds to set hole size (21,22).We attempt to address this automatically in NanoWrap, but to avoid obvious misrepresentation of underlying structure, the topology modification routines used during NanoWrap's shrink-wrapping iterations are quite conservative about when they will change the topology -i.e. they will only change topology if it is very obvious the current topology is wrong.They are also more reliable at removing spurious connections and holes than adding missing connections or holes.Because we bias the starting surface to be slightly outside the structure of interest, we are more likely to miss holes in areas of low density than to generate spurious holes.To assess the ability of our method to accurately capture membrane topology, we performed simulation of a synthetic structure with a range of hole sizes (Supplementary Figure S8).To avoid gaming the analysis we used an automated algorithm to set the threshold for the starting surface (see below).We also chose the curvature weight so as to minimise RMS error, not hole count.Under these conditions, the ability to capture a hole depends on a combination of localisation density, precision, and background level.At high densities and low levels of background, we can typically capture holes that are larger than around 3 times the localization precision, although we occasionally miss a hole we really should have captured.Manual tweaking of the parameters could have achieved improved the results.Correctly estimating surface genus is an area where we believe there is scope for improvement in our algorithm, especially in deciding when to modify the topology as we iterate.

Automatic threshold selection for starting surfaces
We expect that the starting surface should be close to the localizations used to generate it, and the median distance from each localization to the starting surface can be used as quality metric for the surface.By trying di erent thresholds within an optimisation routine (we use bisection) a threshold value leading to a good starting surface can be chosen.It practice it is helpful to aim for a distance (e.g. 10 to 40 nm) which puts the localizations just inside the surface and avoids the generation of spurious holes.This algorithm can be activated by selecting the 'Auto threshold' box in the Dual Marching Cubes settings, and will almost always generate good starting surfaces.Manual tweaking can, however, be useful if there are topological features (e.g.holes) which are right at the limit of being resolvable with the given localisation precision and densities.

Rational choice of curvature weighting
Similar to choosing a threshold for the starting surface, the distance between localizations and the reconstructed membrane can be a useful metric to determine whether an appropriate curvature weight has been selected.This time, however, we look at the distribution of distances rather than their median.Localizations are expected to scatter about the true membrane position due to their localisation error, with the distance between a localisation and the true position of that molecule expected to follow a scaled Chi(3) distribution.The choice of curvature weight is appropriate if the distribution of distances between localizations and the membrane is similar to that predicted by the localisation error.If the spread is larger, the curvature weight is too high, but if the spread is less than that predicted by localization error the algorithm is over-fitting to the points and the curvature weight should be increased.For performance reasons, we have not automated this, and setting the curvature weight remains a manual process.A function to plot the distance distributions is available under Mesh->Show shrinkwrap residuals within PYMEVis.
These distributions are also useful for showing if the number of iterations is suitable -if the residual distribution is heavily weighted to negative distances (inside the surface) this indicates that either more iterations are likely needed to allow the surface to converge.We manually experimented with the input hyperparameters of each algorithm to achieve the best surface in each case.A grid search may yield better surfaces for each of these methods, but we do not expect significant improvement, due to the nature of the input data.Both SPR and Point2Mesh require point normals as input.Point normals were estimated using Meshlab's (48) estimator.Note that these methods may perform well in the presence of true normals for a point cloud, but normal estimation usually relies on neighboring points, including noise points that confound the estimate.Point2Mesh requires the user to provide an initial starting surface.We provided it the same starting surface we used for our approach, which is a dual marching cubes (30) approximation on the density of the point cloud.The localization point cloud is shown in green and the surface fitting these points is shown in magenta.B 1 has area-minimizing tendencies and the surface pulls through the point cloud.B 2 's preservation of area means that areas of the mesh far from point influence can establish stable "blebs" protruding from the true structure.The combination of forces subject to a point influence weighting U creates an accurate reconstruction of a surface.See the methods section for more details.

CGAL Alpha Wrap 3 Rhinoceros 8 Shrink Wrap Point2Mesh
Figure S3: A NanoWrap surface fit over 39 iterations, remeshed every 0, 5, 10, 15, 20 iterations.The structure does not change significantly depending on remesh frequency.Figure S7: To test whether the optimal curvature weight depends on the underlying structure curvature we simulated a scaled version of the double-torus structure that had about twice the radius of curvature, whilst holding localisation precision, density, and noise levels constant in the middle of our range.There is only a very slight change in the optimal curvature weight, although the flatter structure shows a faster degradation in performance when the curvature weight is decreased below the optimum.This can be attributed to an additional regularising e ect of the finite mesh size at small radii of curvature.

SPR
Figure S8: Ability of NanoWrap to correctly estimate the number of holes (genus) for a test structure using an automatically thresholded starting surface.Our ability to correctly estimate hole count is broadly comparable to SPR, but degrades a bit more gracefully at poor localization density.The genus estimation is very sensitive to the quality of the starting surface, and in the few cases where SPR correctly estimates the genus and we do not, the correct genus is easily obtained by manually tweaking the parameters for starting surface generation.This implies that there is still scope to improve the algorithm we use for intial surface estimation.

Figure S1 :
Figure S1: A comparison of our method to other common methods used to generate surfaces from a point cloud.Top left is an example point cloud from a subregion of an endoplasmic reticulum.A surface is generated from this point cloud with NanoWrap, Screened Poisson Reconstruction (18), CGAL's Alpha Wrap (22), Rhinoceros 8's shrink wrap algorithm (47) and Point2Mesh (23).NanoWrap allows the surface to adapt per-point based on localization positions and achieves a better representation of the underlying structure than other methods, which strictly adhere to the point cloud or generate a surface at a constant o set from each point.We manually experimented with the input hyperparameters of each algorithm to achieve the best surface in each case.A grid search may yield better surfaces for each of these methods, but we do not expect significant improvement, due to the nature of the input data.Both SPR and Point2Mesh require point normals as input.Point normals were estimated using Meshlab's (48) estimator.Note that these methods may perform well in the presence of true normals for a point cloud, but normal estimation usually relies on neighboring points, including noise points that confound the estimate.Point2Mesh requires the user to provide an initial starting surface.We provided it the same starting surface we used for our approach, which is a dual marching cubes (30) approximation on the density of the point cloud.

Figure S2 :
Figure S2: This figure demonstrates the independent influence of curvature forces B 1 and B 2 on the final NanoWrap surface.The localization point cloud is shown in green and the surface fitting these points is shown in magenta.B 1 has area-minimizing tendencies and the surface pulls through the point cloud.B 2 's preservation of area means that areas of the mesh far from point influence can establish stable "blebs" protruding from the true structure.The combination of forces subject to a point influence weighting U creates an accurate reconstruction of a surface.See the methods section for more details.

Figure S4 :
Figure S4:The radius of curvature of a reconstructed surface of a tapered cylinder as a function of density.The tapered cylinder's diameter scales between 30 and 200 nm over a 2000 nm distance as the square of the distance along the cylinder.The left column shows a cross section of the fit point cloud (blue) and the reconstructed surface (orange).The right column shows the estimated radius of curvature for each vertex of the fit surface along with an error bar indicating the full range of curvature estimates extracted from the vertices for each expected radius of curvature.As the density increases, the minimum radius of curvature that can be estimated decreases.The point cloud localization precision is fixed at 11 nm.

Figure S5 :
Figure S5:The radius of curvature of a reconstructed surface of a tapered cylinder as a function of localization precision.The tapered cylinder's diameter scales between 30 and 200 nm over a 2000 nm distance as the square of the distance along the cylinder.The left column shows a cross section of the fit point cloud (blue) and the reconstructed surface (orange).The right column shows the estimated radius of curvature for each vertex of the fit surface along with an error bar indicating the full range of curvature estimates extracted from the vertices for each expected radius of curvature.As the precision increases, the minimum radius of curvature that can be estimated scales to match.The point cloud density is fixed at 218 points µm 2 .

Figure S6 :
Figure S6: NanoWrap's sensitivity of fit performance to the chosen value of curvature weight at di erent values of structure density and background level.The optimal curvature weight increases with both increasing point density and background, but in each case there is a broad minimum.

Figure S9 :
Figure S9: Left, Tubular ER from the lower left portion of Figure 3, colored by mean curvature.Lookup table:-0.01 to 0.01 nm 1 .Scale bar is 200 nm.Right, Histogram of ⇡ = 2/: max,; where : max,; is the maximal principal curvature of vertex ;.The mean and standard error of the mean for the diameter distribution are reported.

Figure S10 :
Figure S10: Distance between two surfaces.Surface generated from Sec61V (ER membrane) localizations (red-white-blue) and surface generated from TOMM20 (outer mitochondrial membrane) localizations (magenta).The Sec61V surface is colored by its distance from the TOMM20 surface (0 to 100 nm; red to blue).A, Full field of view, displaying both surfaces.B, Region of interest, showing both surfaces.C, The Sec61V surface alone.Scale bars are 1 µm.

Figure S11 :
Figure S11: Data from user-guide-data.hdfdisplayed as 10 nm point sprites with 0.2 alpha.The points are colored by their location in the z-direction according to the lookup table on the right.

Figure S12 :
Figure S12: Ideal parameters to use to create an isosurface based on user-guide-data.hdf.

Figure S13 :
FigureS13: The isosurface generated using the parameters in S12.The isosurface is displayed using the wireframe method to make it easy to assess how well it fits the point cloud.

Figure S14 :
Figure S14: Ideal parameters to create a surface based on user-guide-data.hdf.

Figure S15 :
Figure S15: The surface generated using the parameters in S14.