DeepSTORM3D: dense three dimensional localization microscopy and point spread function design by deep learning

DeepSTORM3D: dense thr ee dimensional localization microscopy and point spr ead function design by deep learning Elias Nehme 1,2 , Daniel Freedman 3 , Racheli Gordon 2 , Boris Ferdman 2,4 , Lucien E. W eiss 2 , Onit Alalouf 2 , Reut Orange 2,4 , T omer Michaeli 1 , and Y oav Shechtman 2,4,* 1 Department of Electrical Engineering, T echnion, 32000 Haifa, Israel 2 Department of Biomedical Engineering & Lorry I. Lokey Center for Life Sciences and Engineering, T echnion, 32000 Haifa, Israel 3 Google Research, Haifa, Israel 4 Russel Berrie Nanotechnology Intitute, T echnion, 32000 Haifa, Israel * Corresponding author: yoavsh@bm.technion.ac.il Abstract Localization microscopy is an imaging technique in which the positions of individual nanoscale point emitters (e.g. fluor escent molecules) are determined at high pr ecision from their images. This is the key ingredient in single/multiple-particle-tracking and several super-r esolution microscopy approaches. Localiza- tion in three-dimensions (3D) can be performed by modifying the image that a point-source creates on the camera, namely , the point-spread function (PSF). The PSF is engineered using additional optical elements to vary distinctively with the depth of the point-source. However , localizing multiple adjacent emit- ters in 3D poses a significant algorithmic challenge, due to the lateral overlap of their PSFs. Here, we train a neural network to receive an image containing densely overlapping PSFs of mul- tiple emitters over a large axial range, and output a list of their 3D positions. Furthermore, we then use the network to design the optimal PSF for the multi-emitter case. W e demonstrate our approach numerically as well as experimentally by 3D STORM imaging of mitochondria, and volumetric imaging of dozens of fluorescently-labeled telomer es occupying a mammalian nu- cleus in a single snapshot. 1 Introduction Determining the nanoscale positions of point emitters forms the basis of localization microscopy techniques such as single particle tracking [ 1 , 2 ], (fluor escence) photoactivated localization microscopy (f)P ALM [ 3 , 4 ], stochastic optical reconstruction mi- croscopy (ST ORM) [ 5 ], and related single molecule localization microscopy (SMLM) methods. These techniques have revolu- tionized biological imaging, revealing cellular processes and structures at the nanoscale [ 6 ]. Notably , most samples of interest extend in thr ee dimensions, necessitating thr ee-dimensional (3D) localization microscopy [ 7 ]. In a standard microscope, the pr ecise z position of an emitter is difficult to ascertain because the change of the point-spread function (PSF) near the focus is appr oximately symmetric. Fur - thermore, outside of this focal range ( ≈ ± 350 nm for a high numerical aperture imaging system), the rapid defocusing of the PSF reduces the signal-to-noise ratio causing the localization precision to quickly degrade. One method to extend the useful z-range and explicitly encode the z position is PSF engineering [ 9 – 11 ]. Here, an additional optical element, e.g. a phase mask, - 0 - 0 c b 3 μm 4 μm 5 μm 6 μm - Standard PSF T etrapod PSF a z Laser SLM Sample Objective T ube lens Dichroic mirror Polarizer Mirror Final image plane Intermediate image plane 4f system 2 μm z = 3D high resolution CNN 2D low resolution Fig. 1. Optical setup and approach overview . a The light emitted from a fluorescent micr oscopic particle is collected by the objective and focused thr ough the tube lens into an image at the intermediate image-plane. This plane is extended using a 4f system with a phase mask placed at the Fourier plane in between the two 4f lenses. b The implemented phase mask (using either a Spatial Light Modulator (SLM) or fabricated fused-silica) dictates the shape of the PSF as function of the emitter ’s axial position. c After training, our CNN r eceives a 2D low resol ution image of overlapping PSFs and outputs a 3D high-resolution volume which is translated to a list of 3D localizations. Blue empty spheres denote simulated GT po- sitions along the surface of an ellipsoid. Red spheres denote CNN detections. The T etrapod PSF is depicted here, however the approach is applicable to any PSF , including those opti- mized by the net itself (Fig. 4 ). Scale bars are 3 µ m. 1 a b X[ μ m] Z[ μ m] Y[ μ m] 2 -5 5 4 0 0 6 5 -5 Ground T ruth MP 2 -5 5 4 0 0 6 5 -5 Ground T ruth Deep X[ μ m] Z[ μ m] Y[ μ m] 0 0.1 0.2 0.3 0.4 Density 0 0.5 1 Jaccard Index [a.u.] Deep Matching Pursuit 0 0.1 0.2 0.3 0.4 Density 0 20 40 60 80 Lateral RMSE [nm] 0 0.1 0.2 0.3 0.4 Density 0 20 40 60 80 Axial RMSE [nm] 2 e m i t t e r s μ m       2 e m i t t e r s μ m       2 e m i t t e r s μ m       Fig. 2. Comparison to MP . a The trained CNN is superior to the matching pursuit appr oach in both detectability (Jaccard index) and in accuracy (Lateral\Axial RMSE). Matching of points was computed with a thr eshold distance of 150 nm using the Hungarian algorithm [ 8 ]. b Example of a simulated frame of density 0.124 h emitters µ m 2 i alongside 3D comparisons of the r ecovered positions by MP (middle) and by the CNN (right). Scale bar is 2 µ m. is placed in the emission path of the microscope, modifying the image formed on the detector [ 12 ] (Fig. 1 a); the axial position can then be recover ed via image processing using a theoretical or experimentally-calibrated PSF model [ 11 , 13 – 15 ]. In practically all applications, it is desirable to be able to localize nearby emitters simultaneously . For example, in super- resolution SMLM experiments, the number of emitters localized per frame determines the temporal resolution. In tracking ap- plications, PSF overlap fr om multiple emitters often precludes localization, potentially biasing results in emitter -dense regions. The problem is that localizing overlapping emitters poses a sig- nificant algorithmic challenge even in 2D localization, and much more so in 3D. Specifically , encoding the axial position of an emitter over large axial ranges (>3 µ m) requir es the use of later- ally large PSFs, e.g. the T etrapod [ 11 , 16 ] (Fig. 1 b), increasing the possibility of overlap. Consequently , while a variety of methods have been developed to cope with overlapping emitters for the in-focus, standard-PSF [ 17 – 19 ], the performance in high-density 3D localization situations is far from satisfactory [ 20 ]. Deep learning has proven to be adept at analyzing micro- scopic microscopy data [ 21 – 26 ], especially for single-molecule localization, handling dense fields of emitters over small axial ranges (<1.5 µ m) [ 19 , 27 – 34 ] or sparse emitters spread over larger ranges [ 35 ]. Moreover , an emer ging application is to jointly de- sign the optical system alongside the data processing algorithm, enabling end-to-end optimization of both components [ 33 , 36 – 43 ]. Her e we present DeepSTORM3D, consisting of two funda- mental contributions to high-density 3D localization microscopy over large axial ranges. First, we employ a convolutional neu- ral network (CNN) for analyzing dense fields of overlapping emitters with engineered PSFs, demonstrated with the large- axial-range T etrapod PSF [ 11 , 16 ]. Second, we design an optimal PSF for 3D localization of dense emitters over a large axial range of 4 µ m. By incorporating a physical-simulation layer in the CNN with an adjustable phase modulation, we jointly learn the optimal PSF (encoding) and associated localization algorithm (decoding). This appr oach is highly flexible and easily adapted for any 3D SMLM dataset parameters, i.e. emitter density , SNRs, and z-range. W e quantify the performance of the method by simulation, and demonstrate the applicability to 3D biological samples, i.e. mitochondria and telomer e. 2 Results T o solve the high-density localization problem in 3D, we trained a CNN that r eceives a 2D image of overlapping T etrapod PSFs spanning an axial range of 4 µ m, and outputs a 3D grid with a voxel-size of 27.5 × 27.5 × 33 nm 3 (Fig. 1 c). For architecture details and learning hyper-parameters see Supplementary Infor- mation sections 1.1 and 3. T o compile a list of localizations, we apply simple thresholding, and local maximum finding on the output 3D grid (Supplementary Information section 3.4). W e compare our method to a fit-and-subtract based Match- ing Pursuit (MP) approach [ 44 ] (see Supplementary Information section 4) as we are unaware of any other methods capable of lo- calizing overlapping T etrapod PSFs. T o quantitatively compare 2 xy a c b 4 0 Experimental μ m Overlay Fig. 3. Super -resolution 3D imaging over a 4 µ m z-range . a Super-r esolved image of mitochondria spanning a ≈ 4 µ m z-range r en- dered as a 2D histogram where z is encoded by color . b Representative experimental frame (top), and render ed frame from the 3D recover ed positions by the CNN overlaid on top (bottom). c Diffraction limited (left), super-r esolved (middle), and cross-section of the super-r esolved image at z = 1.5 µ m (right). Scale bars ar e 3 µ m. our method with MP solely in terms of density , we simulated emitters with high signal-to-noise ratio (30K signal photons, 150 background photons per pixel) at 10 dif ferent densities ranging from 1 to 75 emitters per 13 × 13 µ m 2 field-of-view . The results are shown in Fig. 2 . As evident in both the Jaccard index (de- fined as T P T P + F P + F N , wher e TP , FP , FN are true positives, false positives, and false negatives [ 20 ]) and the lateral/axial RMSE (Fig. 2 a) the CNN achieves remarkable performance in localizing high-density T etrapods. In the single-emitter (very low density) case, where the performance of the CNN is bounded by the discretization on the 3D grid, the RMSE of the MP localization is lower (better). This is because for a single-emitter , MP is equiv- alent to a continuous Maximum Likelihood Estimator (MLE) (Supplementary Information section 4), which is asymptotically optimal [ 45 ], whereas the CNN’s precision is bounded by pixila- tion of the grid ( i.e. half voxel of 13.75 nm in xy and 16.5 nm in z). However , quickly beyond the single-emitter case, the CNN drastically outperforms MP . A similar result was obtained when compared to a leading single-emitter fitting method [ 15 ] appli- cable also for the multiple emitter case [ 20 ] (see Supplementary Information section 6). Next, we validated our method for super-r esolution imaging of fluorescently labeled mitochondria in COS7 cells (Fig. 3 ). W e acquired 20K diffraction limited frames of a 50 × 30 µ m 2 FOV and localized them using the CNN in about ≈ 10 hours, resulting in ≈ 360K localizations. The T etrapod PSF was implemented using a fabricated fused-silica phase-mask (see Supplementary Information section 7.1). The estimated resolution was ≈ 40 nm in xy , and ≈ 50 nm in z (see Supplementary Information section 7.2). T o visually evaluate localization performance in a single frame (Fig. 3 b top), we regenerated the corresponding 2D low-resolution image, and overlayed the recovered image with a uniform photons scaling on top of the experimental frame 3 b bottom). As seen in the overlay image, the emitter PSFs (3D positions) are faithfully r ecovered by the CNN. Mor eover , emitters with extremely low number of signal photons wer e ignored. The T etrapod is a special PSF that has been optimized for the single emitter case by Fisher Information maximization [ 11 , 16 ]. However , when considering the multiple-emitter case, an intriguing question arises: What is the optimal PSF for high density 3D localization over a large axial range? T o answer this question we need to r ethink the design metric; extending the Fisher Information criterion [ 11 ] to account for emitter density is not-trivial, and while it is intuitive that a smaller-footprint PSF would be preferable for dense emitters, it is not clear how to 3 Z Y X Simulated emitters in 3D Learned phase mask Simulated 2D image 3D localizations Focal plane Fourier plane T ube lens Image plane Reconstruction volume Decoding Encodin g Learned PSF a b c - CNN 0 0.1 0.2 0.3 0.4 Density 0.5 0.6 0.7 0.8 0.9 Jaccard Index 2 e m i t t e r s μ m       0 0.1 0.2 0.3 0.4 Density 20 30 40 Axial RMSE [nm] 2 e m i t t e r s μ m       0 0.1 0.2 0.3 0.4 20 30 40 Lateral RMSE [nm] T etrapod Learned Learned T etrapod 2 μm z = 3 μm 4 μm 5 μm 6 μm - Y Z X Fig. 4. PSF learning for high density 3D imaging . a Simulated 3D emitter positions ar e fed to the image formation model to sim- ulate their low r esolution CCD image (Encoding). Next, this image is fed to a CNN that tries to recover the simulated emitter po- sitions (Decoding). The difference between the simulated positions and the positions recovered by the CNN is used to jointly op- timize the phase mask at the Fourier plane, and the r ecovery CNN parameters. b Simulation of the learned PSF as function of the emitter axial position (left). 3D isosurface rendering of the learned PSF (right). c Example frame of density 0.197 h emitters µ m 2 i (top) with the same simulated emitter positions, using the T etrapod (left) and the learned PSF (right). Jaccard index (bottom) and lateral \axial RMSE comparison (right) between two CNNs with the same ar chitecture, one trained to r ecover 3D positions fr om 2D im- ages of T etrapod PSF (black), and the second trained to recover 3D positions fr om 2D images of the learned PSF (orange). Scale bars are 3 µ m. mathematically balance this demand with the r equirement for high localization precision per emitter . Our PSF-design logic is based on the following: since we have already established that a CNN yields superior reconstruc- tion for high-density 3D localization, we are interested in a PSF (encoder) that would be optimally localized by a CNN (decoder) . Ther e- fore, in contrast to a sequential paradigm where the PSF and the localization algorithm are optimized separately , we adopt a co-design appr oach (Fig. 4 a). T o jointly optimize the PSF and the localization CNN, we introduce a differ entiable physical simula- tion layer , which is parametrized by a phase mask that dictates the micr oscope’s PSF . This layer encodes 3D point sources to their respective low-resolution 2D image (see Supplementary Information section 2). This image is then fed to the localization CNN which decodes it and r ecovers the underlying 3D source positions. During training, the net is presented with simulated point sources at random locations and, using the differ ence be- tween the CNN recovery and the simulated 3D positions, we optimize both the phase mask and the localization CNN param- eters in an end-to-end fashion. The learned PSF (Fig. 4 b) has a small lateral footprint, which is critical for minimizing overlap at high densities. Moreover , the learned phase mask twists in a spiral trajectory causing the PSF to rapidly r otate throughout the axial range, a trait that was previously shown to be valuable 4 [ μ m] Z [ μ m] Y [ μ m] X a b c d e Experimental X Y Z Y X Z Reconstruction Experimental Reconstruction [ ] μ m phot ons pi xe l s       2.7 3.2 3.7 4.2 Depth 1000 1500 2000 2500 Mean Intensity Measured Fit Recovered Depth [ μ m] Y [ μ m] X [ μ m] Z Coverslip Cell Nucleus T elomeres 3um scale bar 3um scale bar Ground truth T rue positive False positive Fig. 5. Three dimensional imaging of telomeres in a single-snapshot . a Schematic of imaging fixed U2OS cells with fluorescent labeled telomeres inside their nucleus. b Focus slice with the standar d PSF inside a U2OS cell nucleus, obtained via a z-scan. The yellow rectangles mark the same emitter in all three orthogonal planes. c Example fit of the mean intensity in sequential axial slices used to estimate the approximate emitter axial position. d Experimental snapshot with the T etrapod PSF (left), render ed image from the 3D r ecovered positions by the T etrapod CNN (middle), and a 3D comparison of the r ecovered positions and the approximate experimental ground truth (right). e Experimental snapshot with the learned PSF (left), render ed image from the 3D recover ed positions by the learned PSF CNN (middle), and a 3D comparison of the recover ed positions and the appr oximate experimental ground truth (right). Scale bars are 3 µ m. for encoding depth [ 9 ]. T o quantify the improvement introduced by our new PSF , we first compare it to the T etrapod PSF in simulations. Specifically , we train a similar r econstruction net for both the T etrapod and the learned PSF using a matching training set composed of sim- ulated continuous 3D positions along with their corresponding 2D low-resolution images. The learned PSF performs similar to the T etrapod PSF for low emitter densities (Fig. 4 c). However , as the density goes up the learned PSF outperforms the T etrapod PSF in both localization precision and in emitter detectability (Jaccard index) (Fig. 4 c). This result is not surprising, as the learned PSF has a smaller spatial foorprint, and hence it is less likely to overlap than the T etrapod (Fig. 4 c). Next, we demonstrate the superiority of the new PSF ex- perimentally by imaging fluorescently labeled telomeres (TRF1- DsRed) in Fixed U2OS cells. The cell contains tens of telomer es squeezed in the volume of a nucleus with ≈ 20 µ m diameter (Fig. 5 a, b). From a single snapshot focused inside the nucleus, the CNN outputs a list of 3D positions of telomeres spanning an axial range of ≈ 3 µ m. Using the T etrapod PSF snapshot, the T etrapod-trained CNN was able to r ecover 49 out of 62 telomeres with a single false positive, yielding a Jaccard index of 0.77 (Fig. 5 d). In comparison, using the learned PSF snapshot, the corre- sponding CNN was able to recover 57 out of the 62 telomeres with only 2 false positives, yielding a Jaccard index of 0.89 (Fig. 5 e). The recover ed positions were compared to approximated ground-truth 3D positions (Fig. 5 c), obtained by axial scanning and 3D fitting (see Supplementary Information section 9). T o qualitatively compare the recover ed list of localizations to the acquired snapshot, we fed this list to the physical simulation layer and generated the matching 2D low-resolution image (Fig. 5 d,e). As verified by the regenerated images, the 3D positions of the telomeres are faithfully r ecovered by the CNNs. Moreover , the misses in both snapshots were either due to local aberra- 5 tions and/or an extremely low number of signal photons (see Supplementary Information section 10 for more experimental results). 3 Discussion In this work we demonstrated 3D localization of dense emitters over a large axial range both numerically and experimentally . The described network architecture exhibits excellent flexibil- ity in dealing with various experimental challenges, e.g. low signal-to-noise ratios and optical aberrations. This versatility is facilitated in three ways: (1) the net was trained solely on simulated data, thus producing sufficiently large datasets for optimization; (2) the phase mask which governs the PSF was optimized with respect to the implementation in the imaging system, i.e. the pixels of the spatial light modular , rather than over a smaller subspace, e.g. Zernike polynomials [ 11 ]; (3) the CNN localization algorithm was designed in coordination with the development of the PSF , thus the system was optimized for the desired output [ 33 ] rather than a pr oxy . Attaining a sufficiently large training dataset has thus far been a major limitation for most applications of CNNs. W ith this limitation in mind, the application of CNNs to single-molecule localization would seemingly be an ideal one, since each emit- ter ’s behavior should be approximately the same. This uni- formity is broken, however , by spatially-varying background, sample density , and variable emitter size in biological samples (Supplementary Information 3.1), all of which diversify datasets and necessitate relevant training data. By implementing an accu- rate simulator , we have shown that it is possible to build a r obust network entirely in silico, generating arbitrarily large, r ealistic datasets with a known gr ound truth to optimize the nets. This aligns with our previous work in 2D SMLM [ 19 ]). For super-r esolution reconstructions using the T etrapod PSF , the simulator was particularly important due to the highly vari- able SNR of emitters in the sample. Her e, our net was able to selectively localize the emitters even in very dense regions by focusing on those with a high SNR (Fig. 3 ). T o optimize a PSF while simultaneously training the net, the simulator was also essential, as it would be prohibitively time consuming to experi- mentally vary the PSF , while recor ding and analyzing images to train the net. An intriguing aspect of our optimization approach is that the optimized PSF is found by continuously varying the pixels of an initialized mask while evaluating the output of the localization net, thus the final result repr esents a local minimum (Fig. 4 ). By changing the initialization conditions, we have recognized several patterns that indicate how the optimal PSF varies with the experimental conditions: namely , density , axial range, and SNR (see Supplementary Information section 1.2). Some of the recurr ent features are intuitive: for example, in dense fields of emitters with limited SNR, the optimized PSFs have a small footprint over the designed axial range, enabling high density and compacting precious signal photons into as few pixels as possible. What distinguishes the net PSFs over predetermined designs is the utilization of multiple types of depth encoding; namely , simultaneously employing astigmatism, rotation, and side lobe movement (Fig. 4 ), all of which have been conceived of and implemented previously , but never simultaneously! This work, therefor e, triggers many possible questions and resear ch directions r egarding its capabilities and limitations. For example, how globally-optimal is the r esulting PSF? Similarly , how sensitive is the resulting PSF and its performance to dif fer- ent loss functions, CNN ar chitectures, initializations (e.g. with an existing phase mask), and the sampled training set of loca- tions? Curr ently , it is unclear how each of these components af- fects the learning pr ocess although we began to partially answer them in simulations (see Supplementary Information section 1.2). Finally , the co-design approach employed her e paves the way to a wide variety of inter esting applications in microscopy wher e imaging systems have traditionally been designed separately from the pr ocessing algorithm. Funding Information Google; H2020 European Research Council Horizon 2020 (802567); Israel Science Foundation (ISF) (852/17); Ollendorff Foundation; T echnion-Israel Institute of T echnology (Career Ad- vancement Chairship); Zuckerman Foundation. Acknowledgements The authors wish to thank the Garini lab for the U2OS cells. W e also thank Jonas Ries for his help with the application of SMAP- 2018 to T etrapod PSFs. W e gratefully acknowledge the support of NVIDIA Corporation with the donation of the T itan V GPU used for this resear ch. W e thank the staff of the Micro-Nano- Fabrication & Printing Unit (MNF-PU) at the T echnion for their assist with the phase mask fabrication. Finally , we thank Google for the resear ch cloud units provided to accelerate this resear ch. Author contributions EN, DF , TM and YS conceived the approach. EN performed the simulations and analyzed the data with contributions fr om all authors. EN, RG, BF , LEW , and OA took the data. RO fabricated the physical phase mask. EN, DF , LEW , TM and YS wrote the paper with contributions from all authors. Competing Interests The authors declare no competing inter ests. Data availability Data will be made available upon reasonable r equest. Code availability Code will be made publicly available. References [1] Y . Katayama, O. Burkacky , M. Meyer , C. Bräuchle, E. Grat- ton, and D. C. Lamb, “Real-time nanomicroscopy via three-dimensional single-particle tracking,” ChemPhysChem , vol. 10, no. 14, pp. 2458–2464, 2009. [2] C. Manzo and M. F . Garcia-Parajo, “A review of progr ess in single particle tracking: from methods to biophysical insights,” Reports on progr ess in physics , vol. 78, no. 12, p. 124601, 2015. 6 [3] E. Betzig, G. H. Patterson, R. Sougrat, O. W . Lindwasser , S. Olenych, J. S. Bonifacino, M. W . Davidson, J. Lippincott- Schwartz, and H. F . Hess, “Imaging intracellular fluorescent proteins at nanometer r esolution,” Science , 2006. [4] S. T . Hess, T . P . Girirajan, and M. D. Mason, “Ultra-high resolution imaging by fluorescence photoactivation local- ization microscopy ,” Biophysical Journal , vol. 91, no. 11, pp. 4258 – 4272, 2006. [5] M. J. Rust, M. Bates, and X. Zhuang, “Sub-diffraction-limit imaging by stochastic optical r econstruction microscopy (STORM),” Natur e Methods , vol. 3, no. 10, pp. 793–795, 2006. [6] S. J. Sahl and W . Moerner , “Super -resolution fluor escence imaging with single molecules,” Curr ent Opinion in Struc- tural Biology , vol. 23, no. 5, pp. 778 – 787, 2013. [7] A. von Diezmann, Y . Shechtman, and W . Moerner , “Three- dimensional localization of single molecules for super- resolution imaging and single-particle tracking,” Chemical reviews , vol. 117, no. 11, pp. 7244–7275, 2017. [8] H. W . Kuhn, “The hungarian method for the assignment problem,” Naval r esearch logistics quarterly , vol. 2, no. 1-2, pp. 83–97, 1955. [9] S. R. P . Pavani, M. A. Thompson, J. S. Biteen, S. J. Lord, N. Liu, R. J. T wieg, R. Piestun, and W . Moerner , “Three- dimensional, single-molecule fluorescence imaging beyond the diffraction limit by using a double-helix point spread function,” Proceedings of the National Academy of Sciences , vol. 106, no. 9, pp. 2995–2999, 2009. [10] B. Huang, W . W ang, M. Bates, and X. Zhuang, “Three- dimensional super-resolution imaging by stochastic opti- cal reconstr uction microscopy ,” Science , vol. 319, no. 5864, pp. 810–813, 2008. [11] Y . Shechtman, S. J. Sahl, A. S. Backer , and W . Moerner , “Op- timal point spread function design for 3d imaging,” Physical review letters , vol. 113, no. 13, p. 133902, 2014. [12] A. S. Backer and W . Moerner , “Extending single-molecule microscopy using optical fourier processing,” The Journal of Physical Chemistry B , vol. 118, no. 28, pp. 8313–8329, 2014. [13] S. Liu, E. B. Kromann, W . D. Krueger , J. Bewersdorf, and K. A. Lidke, “Three dimensional single molecule localiza- tion using a phase retrieved pupil function,” Optics express , vol. 21, no. 24, pp. 29462–29487, 2013. [14] H. P . Babcock and X. Zhuang, “Analyzing single molecule localization microscopy data using cubic splines,” Scientific reports , vol. 7, no. 1, p. 552, 2017. [15] Y . Li, M. Mund, P . Hoess, J. Deschamps, U. Matti, B. Ni- jmeijer , V . J. Sabinina, J. Ellenberg, I. Schoen, and J. Ries, “Real-time 3d single-molecule localization using experimen- tal point spread functions,” Nature methods , vol. 15, no. 5, p. 367, 2018. [16] Y . Shechtman, L. E. W eiss, A. S. Backer , S. J. Sahl, and W . Moerner , “Pr ecise three-dimensional scan-free multiple- particle tracking over lar ge axial ranges with tetrapod point spread functions,” Nano letters , vol. 15, no. 6, pp. 4194–4199, 2015. [17] J. Min, C. V onesch, H. Kirshner , L. Carlini, N. Olivier , S. Holden, S. Manley , J. C. Y e, and M. Unser , “Falcon: fast and unbiased reconstruction of high-density super- resolution microscopy data,” Scientific r eports , vol. 4, p. 4577, 2014. [18] N. Boyd, G. Schiebinger , and B. Recht, “The alternating descent conditional gradient method for sparse inverse problems,” SIAM Journal on Optimization , vol. 27, no. 2, pp. 616–639, 2017. [19] E. Nehme, L. E. W eiss, T . Michaeli, and Y . Shechtman, “Deep-storm: super -resolution single-molecule microscopy by deep learning,” Optica , vol. 5, no. 4, pp. 458–464, 2018. [20] D. Sage, T .-A. Pham, H. Babcock, T . Lukes, T . Pengo, J. Chao, R. V elmurugan, A. Herbert, A. Agrawal, S. Colabrese, et al. , “Super-r esolution fight club: assessment of 2d and 3d single- molecule localization microscopy software,” Nature methods , vol. 16, no. 5, p. 387, 2019. [21] Y . Rivenson, Y . Zhang, H. Günaydın, D. T eng, and A. Oz- can, “Phase recovery and holographic image reconstr uction using deep learning in neural networks,” Light: Science & Applications , vol. 7, no. 2, p. 17141, 2018. [22] T . Nguyen, Y . Xue, Y . Li, L. T ian, and G. Nehmetallah, “Deep learning approach for fourier ptychography mi- croscopy ,” Optics express , vol. 26, no. 20, pp. 26470–26484, 2018. [23] M. W eigert, U. Schmidt, T . Boothe, A. Müller , A. Dibrov , A. Jain, B. W ilhelm, D. Schmidt, C. Broaddus, S. Culley , et al. , “Content-aware image restoration: pushing the limits of fluorescence micr oscopy ,” Nature methods , vol. 15, no. 12, p. 1090, 2018. [24] Y . Rivenson, T . Liu, Z. W ei, Y . Zhang, K. de Haan, and A. Ozcan, “Phasestain: the digital staining of label-free quantitative phase microscopy images using deep learning,” Light: Science & Applications , vol. 8, no. 1, p. 23, 2019. [25] T . Liu, K. de Haan, Y . Rivenson, Z. W ei, X. Zeng, Y . Zhang, and A. Ozcan, “Deep learning-based super-resolution in coherent imaging systems,” Scientific reports , vol. 9, no. 1, p. 3926, 2019. [26] J. T . Smith, R. Y ao, N. Sinsuebphon, A. Rudkouskaya, J. Mazurkiewicz, M. Barroso, P . Y an, and X. Intes, “Ultra- fast fit-free analysis of complex fluor escence lifetime imag- ing via deep learning,” bioRxiv , p. 523928, 2019. [27] N. Boyd, E. Jonas, H. P . Babcock, and B. Recht, “Deeploco: Fast 3d localization microscopy using neural networks,” BioRxiv , p. 267096, 2018. [28] W . Ouyang, A. Aristov , M. Lelek, X. Hao, and C. Zimmer , “Deep learning massively accelerates super -resolution local- ization microscopy ,” Nature biotechnology , 2018. [29] B. Diederich, P . Then, A. Jügler , R. Förster , and R. Heintz- mann, “cellstorm—cost-effective super -resolution on a cell- phone using dstorm,” PloS one , vol. 14, no. 1, p. e0209827, 2019. [30] J. M. Newby , A. M. Schaefer , P . T . Lee, M. G. For est, and S. K. Lai, “Convolutional neural networks automate detection for tracking of submicron-scale particles in 2d and 3d,” Pro- ceedings of the National Academy of Sciences , vol. 115, no. 36, pp. 9026–9031, 2018. [31] P . Zelger , K. Kaser , B. Rossboth, L. V elas, G. Schütz, and A. Jesacher , “Three-dimensional localization microscopy us- ing deep learning,” Optics express , vol. 26, no. 25, pp. 33166– 33179, 2018. [32] K. Liu, H. Qiao, J. W u, H. W ang, L. Fang, and Q. Dai, “Fast 3d cell tracking with wide-field fluorescence microscopy through deep learning,” arXiv preprint , 2018. [33] E. Hershko, L. E. W eiss, T . Michaeli, and Y . Shecht- man, “Multicolor localization microscopy and point-spread- function engineering by deep learning,” Optics express , vol. 27, no. 5, pp. 6158–6183, 2019. [34] A. Speiser , S. C. T uraga, and J. H. Macke, “T eaching deep neural networks to localize sources in super-r esolution mi- croscopy by combining simulation-based learning and un- 7 supervised learning,” arXiv preprint , 2019. [35] P . Zhang, S. Liu, A. Chaurasia, D. Ma, M. J. Mlodzianoski, E. Culurciello, and F . Huang, “Analyzing complex single- molecule emission patterns with deep learning,” Nature methods , vol. 15, no. 11, p. 913, 2018. [36] A. Chakrabarti, “Learning sensor multiplexing design through back-propagation,” in Advances in Neural Infor- mation Processing Systems , pp. 3081–3089, 2016. [37] R. Horstmeyer , R. Y . Chen, B. Kappes, and B. Judkewitz, “Convolutional neural networks that teach microscopes how to image,” arXiv preprint , 2017. [38] A. T urpin, I. V ishniakou, and J. D. Seelig, “Light scattering control with neural networks in transmission and reflec- tion,” arXiv preprint , 2018. [39] H. Haim, S. Elmalem, R. Giryes, A. M. Bronstein, and E. Marom, “Depth estimation from a single image using deep learned phase coded mask,” IEEE T ransactions on Com- putational Imaging , vol. 4, no. 3, pp. 298–310, 2018. [40] L. He, G. W ang, and Z. Hu, “Learning depth from single images with deep neural network embedding focal length,” IEEE T ransactions on Image Processing , vol. 27, no. 9, pp. 4676– 4689, 2018. [41] V . Sitzmann, S. Diamond, Y . Peng, X. Dun, S. Boyd, W . Hei- drich, F . Heide, and G. W etzstein, “End-to-end optimization of optics and image processing for achr omatic extended depth of field and super -resolution imaging,” ACM T rans- actions on Graphics (TOG) , vol. 37, no. 4, p. 114, 2018. [42] J. Chang and G. W etzstein, “Deep optics for monocular depth estimation and 3d object detection,” arXiv pr eprint arXiv:1904.08601 , 2019. [43] Y . W u, V . Boominathan, H. Chen, A. Sankaranarayanan, and A. V eeraraghavan, “Phasecam3d—learning phase masks for passive single view depth estimation,” [44] Y . Shechtman, L. E. W eiss, A. S. Backer , M. Y . Lee, and W . Moerner , “Multicolour localization micr oscopy by point- spread-function engineering,” Nature photonics , vol. 10, no. 9, p. 590, 2016. [45] P . J. Bickel and K. A. Doksum, Mathematical Statistics: Basic Ideas and Selected T opics, V olumes I-II Package . Chapman and Hall/CRC, 2015. [46] I. Bronshtein, E. Kepten, I. Kanter , S. Berezin, M. Lindner , A. B. Redwood, S. Mai, S. Gonzalo, R. Foisner , Y . Shav- T al, et al. , “Loss of lamin a function increases chromatin dynamics in the nuclear interior ,” Nature communications , vol. 6, p. 8044, 2015. [47] L. Nahidiazar , A. V . Agronskaia, J. Broertjes, B. van den Broek, and K. Jalink, “Optimizing imaging conditions for demanding multi-color super resolution localization mi- croscopy ,” PLoS One , vol. 11, no. 7, p. e0158884, 2016. [48] M. Ovesn ` y, P . K ˇ rížek, J. Borkovec, Z. Švindrych, and G. M. Hagen, “Thunderstorm: a comprehensive imagej plug-in for palm and storm data analysis and super-resolution imaging,” Bioinformatics , vol. 30, no. 16, pp. 2389–2390, 2014. [49] J. Schindelin, I. Arganda-Carreras, E. Frise, V . Kaynig, M. Longair , T . Pietzsch, S. Preibisch, C. Rueden, S. Saalfeld, B. Schmid, et al. , “Fiji: an open-source platform for biological-image analysis,” Nature methods , vol. 9, no. 7, p. 676, 2012. [50] F . Y u and V . Koltun, “Multi-scale context aggregation by dilated convolutions,” arXiv pr eprint arXiv:1511.07122 , 2015. Methods Sample preparation COS7 cells were grown for 24 hr on cleaned 22 × 22 mm, 170 µ m thick coverslips in 6-well plate in Dulbecco’s Modified Eagle Medium (DMEM) W ith 1g/l D-Glucose (Low Glucose), supple- mented with Fetal bovine serum, Penicillin-Streptomycin and glutamine, at 37 o C , and 5% CO 2 . The cells were fixed with 4% paraformaldehyde and 0.2% glutaraldehyde in PBS, pH 6.2, for 45 min, washed and incubated in 0.3M glycine/PBS solution for 10 minutes. The coverslips were transferred into a clean 6-well plate and incubated in a blocking solution for 2 hr (10% goat serum, 3% BSA, 2.2% glycine, and 0.1% T riton-X in PBS, filtered with 0.45 µ m PVDF filter unit, Millex). The cells were then immunostained overnight with anti TOMM20-AF647 (Ab- cam , ab209606) 1:230 diluted in the blocking buffer , and washed X5 with PBS. Cover glasses (22 × 22 mm, 170 µ m thick) were cleaned in an ultrasonic bath with 5% Decon90 at 60 o C for 30 min, then washed with water , incubated in ethanol absolute for 30 min and sterilized with 70% filtered ethanol for 30 min. U2OS cells were grown on cleaned 0.18 µ m coverslips in 12-well plate in Dulbecco’s Modified Eagle Medium (DMEM) W ith 1g/l D-Glucose (Low Glucose), supplemented with Fetal bovine serum, Penicillin-Str eptomycin and glutamine, at 37 o C , and 5% CO 2 . The day after cells were transfected with dsRed- TRF1 plasmid [ 46 ] using Lipofectamine 3000 reagent. 24 hr after transfection cells wer e fixed with 4% paraformaldehyde for 20 minutes, washed 3 times with PBS and attached to a slide together with mounting medium. Optical setup Imaging experiments were performed on the experimental sys- tem shown schematically in Fig. 1 b. The 4f optical processing system was built alongside the side-port of a Nikon Eclipse T i inverted fluorescence micr oscope, with a × 100/1.45 numerical aperture (NA) oil-immersion objective lens (Plan Apo × 100/1.45 NA, Nikon). STORM imaging For super-resolution imaging, a PDMS chamber was attached to a glass coverslip containing fixed COS7 cells. Blinking buffer (100 mM β -mercaptoethylamine hydr ochloride, 20% sodium lac- tate, and 3% OxyFluor (Sigma, SAE0059), modified from [ 47 ], was then added and a glass coverslip was placed on top to prevent evaporation. Low-intensity illumination for r ecording diffraction-limited images was applied using a T opica laser (640 nm), on the Nikon TI imaging setup described pr eviously , and recor ded with an EMCCD (iXon, Andor) in a standar d imag- ing setup. For super-resolution blinking using the T etrapod PSF , high-intensity (1W at the back of the objective lens) 640 nm light was applied using a 638 nm 2000 mW red dot laser module, whose beam shape was cleaned using a 25 µ m pinhole (Thorlabs) in coor dination with low-intensity (<5mW) 405 nm light. Emission light was filtered through a 500 nm Long pass dichroic and a 650 nm long pass (Chroma), projected through a 4f system containing the dielectric T etrapod phase mask (see Sup- plementary Information section 7.1), and imaged on a Prime95b Photometrics camera. 8 Super-resolution image rendering Prior to rendering the super-resolved image (Fig. 3 b), we first corrected for sample drift using the ThunderSTORM [ 48 ] Im- ageJ [ 49 ] plugin. Afterwards, we render ed the 3D localizations as a 2D average shifted histogram, with color encoding the z position. T elomere imaging For telomere imaging, the 4f system consisted of two f=15 cm lenses (Thorlabs), a linear polarizer (Thorlabs) to filter out the light that is polarized in the unmodulated direction of the SLM, a 1920 × 1080 pixel SLM (PLUTO-VIS, Holoeye) and a mirror for beam-steering. A sCMOS camera (Prime95B, Photometrics) was used to recor d the data. The sample was illuminated with 561 nm fiber-coupled laser-light sour ce (iChrome MLE, T optica). The excitation light was reflected up through the microscope objective by a multibandpass dichroic filter (TRF89902-EM - ET - 405/488/561/647nm Laser Quad Band Set, Chroma). Emis- sion light was filtered by the same dichroic and also filtered by another 617 nm band pass filter (FF02-617/73 Semrock). CNN Architecture In a nutshell, our localization CNN architectur e is composed of 3 main modules. First, a multi-scale context aggr egation module process the input 2D low-resolution image and extracts features with a gr owing receptive field using dilated convolutions [ 50 ]. Second, an up-sampling module increase the lateral resolution of the predicted volume by a factor of x4. Finally , the last mod- ule refines the depth and lateral position of the emitters and outputs the predicted vacancy-grid. For more details regar ding the architectur e see Supplementary Information section 1. 9