Neurocomputing 451 (2021) 279–289

Contents lists available at ScienceDirect

Neurocomputing
journal homepage: www.elsevier.com/locate/neucom

Fast and robust multiplane single-molecule localization microscopyusing a deep neural network
Toshimitsu Aritake a,⇑, Hideitsu Hino b, Shigeyuki Namiki c, Daisuke Asanuma c, Kenzo Hirose c,Noboru Murata aabcSchool of Advanced Science and Engineering, Waseda University, 3-4-1 Okubo, Shinjuku-ku, Tokyo 169-8555, JapanThe Institute of Statistical Mathematics 10-3 Midori-cho, Tachikawa, Tokyo 190-8565, Japan
Graduate School of Medicine, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8643, Japana r t i c l e

i n f oArticle history:
Received 6 March 2020
Revised 21 August 2020
Accepted 14 April 2021
Available online 20 April 2021
Communicated by Zidong Wang
Keywords:
3D single-molecule localization microscopy
Multi-focal plane microscopy
Lateral drift
Convolutional neural networka b s t r a c t
Single-molecule localization microscopy is a widely used technique in biological research for measuring
the nanostructures of samples smaller than the diffraction limit. This study uses multifocal plane microscopy and addresses the three-dimensional (3D) single-molecule localization problem, where lateral andaxial locations of molecules are estimated. However, when multifocal plane microscopy is used, the estimation accuracy of 3D localization is easily deteriorated by the small lateral drifts of camera positions. A
3D molecule localization problem was presented along with the lateral drift estimation as a compressedsensing problem. A deep neural network (DNN) was applied to solve this problem accurately and efficiently. The results show that the proposed method is robust to lateral drift and achieves an accuracyof 20 nm laterally and 50 nm axially without an explicit drift correction.
Ó 2021 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND
license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction
Fluorescence microscopy is a widely used technique in biological research to analyze in vivo structures of samples. However,owing to the diffraction limit of light, the resolution of conventional fluorescence microscopy is limited to approximately
200 nm laterally and 500 nm axially. To overcome this diffraction
limit, several super-resolution microscopy methods [1], including
single-molecule localization microscopy (SMLM) [2,3], have been
proposed. The fundamental problem in super-resolution microscopy is estimating the true molecular distribution from an
observed image. In SMLM, only a few molecules are activated at
one time using photoactivatable molecules. Therefore, the positions of the activated molecules can be accurately estimated using
a localization algorithm, such as Gaussian fitting. Then, a highresolution image can be obtained by integrating the localization
results from many frames into one image.
In many biological studies, three-dimensional (3D) imaging
techniques are important for observing the 3D structures of samples and various 3D fluorescence microscopy techniques have been
proposed [4]. In the past decade, SMLM methods have also been
extended to achieve 3D super-resolution. The most commonly⇑ Corresponding author.

used method for 3D SMLM is point spread function (PSF) engineering. In particular, several types of filters [5–7] have been proposedto achieve 3D localization. In these methods, the 3D locations ofmolecules are estimated from the differences in the shape of the
PSFs. Hence, additional instruments, such as a cylindrical lenses
or phase masks, are required for the optical system, and, the reconfiguration of the optical system for other applications is difficult.
Meanwhile, multifocal plane microscopy (MUM) [8,9], which
we used in this work, is a simple extension of two-dimensional(2D) SMLM to 3D using multiple cameras. Fig. 1(a) shows the
MUM using four cameras. The 3D locations of the molecules are
estimated by the images obtained from multiple focal planes. Additional instruments are not necessary for MUM; thus, the opticallayout is simple, and the reconfiguration of the optical system forother applications is relatively easy. However, another problem
arises in 3D SMLM using MUM: any lateral drift in the camera positions affects the localization quality. When we use MUM, the camera positions may have sub-pixel-sized lateral drift as shown in
Fig. 1(b). Because images from the focal planes are obtained
through the same objective lens, the focal planes do not share identical coordinate systems owing to the lateral drifts. The locations of
the molecules are estimated from the observations from each focal
plane, which is affected by the lateral drift; hence, the drift in thecamera positions make the estimation less accurate.

E-mail address: toshimitsu.aritake@ruri.waseda.jp (T. Aritake).
https://doi.org/10.1016/j.neucom.2021.04.050
0925-2312/Ó 2021 The Authors. Published by Elsevier B.V.
This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

T. Aritake, H. Hino, S. Namiki et al.

Neurocomputing 451 (2021) 279–289

Fig. 1. (a) Ideal multi-focal plane microscopy. Images are simultaneously obtained at each focal plane using an objective lens, and the lateral coordinates of the focal planes
are identical. (b) Lateral drift in MUM. Owing to the lateral drift, the lateral coordinates of the focal planes are not identical.

pose a CNN which stacks several deconvolution layers which areused in fast super-resolution convolutional neural network(FSRCNN) [23] for the single image super-resolution of natural
images. The architecture of this network is based on the iterativealgorithm to solve the compressed sensing problem as unrolling
methods. Namely, when we used iterative optimization to solve
the compressed sensing problem of molecule localization, the resolution of an input image was enhanced at each iteration. Therefore, the proposed network stacks several deconvolution layers
so that the resolution of an input image is enhanced at each layer.
As opposed to the super-resolution of natural images, which
restores the edges or the textures of an image, the molecule localization only requires to estimate the locations of molecules in an
observed image. Therefore, we reformulated the molecule localization problem as a binary classification problem where the network
estimates whether a molecule exists in a given pixel of the target
resolution. In addition, this network is trained to be robust to lateral drifts in camera positions using the dataset which considers
the effect of lateral drifts so that the network can localize molecules accurately without explicit drifts correction. Then, we use
the trained network to estimate molecule locations in an observed
image.
Our contributions are summarized as follows. First, we formulate the 3D molecule localization problem using MUM as a compressed sensing problem to clarify the inverse problem that we
attempt to address in this study. To cope with lateral drift in
MUM, the problem aims to estimate the locations of molecules
and the amount of lateral drifts of cameras. Then, we point out thatthe problem becomes intractable owing to the high computational
cost when the size of the input images becomes larger or the target
resolution of the localization becomes higher. Second, we use CNN
to solve the inverse problem approximately, which significantly
reduces the computational cost of solving the problem. In addition,
this network is trained to be robust to the lateral drift in camera
positions. Therefore, the network correctly estimates the molecule
distribution from the given input image without explicitly estimating the extent of lateral drift. To the best of the authors’ knowledge,
this work is the first attempt to tackle the problem of lateral drifts
in MUM for SMLM.
The rest of this paper is organized as follows. The observation
model and the formulation of the single-molecule localization
problem using MUM are explained in Section 2. Then, the problemof the molecule localization using MUM and the details of the

Therefore, for an accurate localization of the molecules, it isessential to compensate for the effect of lateral drift. The MUM
proposed in [9] used a cross-correlation approach [10,11] to correct drift between cameras. However, the cross-correlation
approach corrects the drift of specimens from a set of localization
results at certain time intervals. Then, the molecules must be localized from the observed image obtained from each focal plane independently. Therefore, MUM information is not completely used to
localize molecules. Moreover, additional computational cost isrequired to correct the drift when localizing the molecules.
In this study, we formulated the simultaneous 3D molecule
localization and lateral drift correction as a compressed sensing
problem [12]. Our formulation uses observed images that are
obtained from different focal planes simultaneously and takes
advantage of the benefits of MUM for molecule localization anddrift correction. By addressing this problem, the locations of molecules are accurately estimated. However, as the size of the input
images becomes larger or the target resolution of the localization
becomes higher, the problem becomes intractable owing to its highcomputational cost.
Recently, deep neural networks (DNNs) have been successfully
applied to a variety of inverse problems in imaging [13]. Conventionally, inverse problems in imaging were formulated as compressed sensing problems and solved by iterative optimization,
which has a high computational cost. Instead of solving the optimization problem iteratively, DNNs use deeply layered networks
and approximate the inverse of the observation process by learning
the mapping from observations to their corresponding inputs.
Intuitively, many iterations of optimization are replaced by a smaller number of trained layers as unrolling methods [14,15].
Recently, DNNs have also been applied to inverse problems of
SMLM, such as molecule localization [16–20], colorization [21],and background estimation [22]. Specifically, convolutional neuralnetworks (CNNs) have achieved a remarkable speedup in 2D and
3D molecule localization. Although training a neural network
requires a large amount of training data, and the training process
can take several hours to several days, the trained network canestimate the molecule locations accurately and efficiently. Moreover, an infinite amount of training data can be generated using
an approximated PSF; hence, DNNs are suitable for SMLM.
Here, a CNN was also used to estimate the locations of molecules, along with the lateral drift of cameras instead of solving acompressed sensing problem by an iterative algorithm. We pro280

Neurocomputing 451 (2021) 279–289

T. Aritake, H. Hino, S. Namiki et al.

proposed method are presented in Section 3. In the experiments,
we used quad-plane microscopy as our MUM method. The detailed
experimental setup and experimental results that validate the proposed algorithm are discussed in Section 4. Finally, our concluding
remarks and the discussion are provided in Section 5.

and corresponding molecule distribution c ¼ ðc1 ; c2 ; . . . ; cm Þ 2 Rm
to approximate the observation yl . Here, most of the elements ofc are zeros, and weight cj is nonzero if a molecule is contained in
the j-th voxel. We denote the coordinate of the center of the j-thvoxel by xhj ðj ¼ 1; 2; . . . mÞ.

2. Formulation

Then, the ði; jÞ-element of matrix H is the fluorescence from amolecule at xhj , observed at xli , and can be written as hðxli ; xhj Þ, where
h is a PSF.

2

2.1. Observation model

We assumed that an observed image contains both shot noise
and additive observation noise. The shot noise follows a Poisson
distribution, whereas the observation noise follows a Gaussian distribution for each observation independently. Hence, the observation can be modeled asyl ¼ Hc þ :

sity at the i-th observation coordinate xli ¼ ðxli1 ; xli2 ; xli3 Þ 2 R3 , andn ¼ kW l Hl . Note that xli is a coordinate on one of the k focal planes;
hence xi3 2 fd1 ; d2 ; . . . ; dk g, where dz ðz ¼ 1; 2; . . . ; kÞ is the depth of
the z-th focal plane. The observed image yl is a convolution of thetrue molecule density and the PSF which depends on the opticalsystem of a microscope. Therefore, yl can be approximated by a linear equation as

ð2Þ

where  is composed of shot noise and observation noise.
Now, the problem of molecule localization is to estimate theweights cj for all j ¼ 1; 2; . . . m from a low-resolution image set yl .
Here, the observation matrix H is an overcomplete matrix
ðn < mÞ; thus, coefficient vector c cannot be recovered by minimizing the noise in Eq. (2). We assum that c is a sparse vector, wheremost of the elements are zeros; therefore, c can be recovered bysolving the following problem:

ð1Þminimize kyl  Hck22 þ kkck1 ;where H 2 R
is an observation matrix. The j-th column of H isfluorescence from the j-th molecule, and K is the number of molecules in an image set. The vector c  ¼ ðc1 ; c2 ; . . . ; cK Þ 2 RK is a molecule distribution, where cj ðj ¼ 1; . . . ; KÞ represents the weight of
the intensity of the j-th molecule.
In theory, the true coordinates of molecules can be used to generate H  ; however, they are not given or known in advance.
Instead, we assumed that the target resolution is given. Hence,
we divided a target 3D space into m ¼ W h  Hh  Dh voxels at
the target resolution, as shown in Fig. 2 (b). We considered allmolecules located at the center of voxels and generated observation matrix H 2 Rnm , instead of H  , where m  K. We used H
nK

3

7
6
6 hðxl2 ; xh1 Þ    hðxl2 ; xhm Þ 776H¼6
7 2 Rnm
..
..
..
7
6
.
.
.
5
4
l
h
l
h
hðxn ; x1 Þ    hðxn ; xm ÞIn this study, we used MUM to observe a sample in the target
3D space. The MUM allows the simultaneous imaging at k focal
planes in the target 3D space using an objective lens. A conceptual
diagram of MUM when k ¼ 4 is shown in Fig. 2(a). The intensitiesof fluorescence are observed as a W l  Hl image at each focal planein the 3D space. Owing to the diffraction limit, low-resolution
images are obtained from the focal planes. Hence, we call a set ofk observed images as a low-resolution image set.
Let yl ¼ ðyl1 ; yl2 ; . . . ; yln Þ 2 Rn be a low-resolution image set
obtained using MUM, where yli is the observed fluorescence inten-yl  H  c  ;hðxl1 ; xh1 Þ    hðxl1 ; xhm Þcð3Þwhich is known as Lasso (see [24] and the references therein). This
type of inverse problem is known as compressed sensing [12].
2.2. 3D PSF
In the above formulation, the true PSF h is also not known gen^ to approximate
erally; therefore, we used a parametric function, h,^ to solve Eq. (3). Theh, and an approximated observation matrix, H,
PSF of MUM is modeled based on the following function:

Fig. 2. (a) Illustration of multi-focal plane microscopy observing a sample in 3D space. (b) 3D space divided into voxels at the target resolution. Molecules located at thecenters of the voxels are used to approximate an observation. The molecules of the frontmost voxels are shown by the blue circles.
281

T. Aritake, H. Hino, S. Namiki et al.

Neurocomputing 451 (2021) 279–289^ l ; xh Þ ¼ hðx^ l ; xl ; xl ; xh ; xh ; xh Þ
hðx
i
j
i1 i2
i3
j1 j2
j3
 l h 2 l h 2ð4Þ
ðxi1 xj1 Þ þðxi2 xj2 Þlhþ b;¼ aðxi3  xj3 Þ exp 2wðxl xh Þ
i3

j3

where a is the peak fluorescence intensity, w is the squared width of
the fluorescence, and b is the background fluorescence intensity.
This PSF is similar to the PSF used in [25] for biplane microscopy.
The squared width wðx3 Þ of the PSF varies depending on the distance between a molecule and the focal planes; it is described by
the following defocus curve:


x 2
x 3
x 4 3
3
3;
wðx3 Þ ¼ w20 1 þ
þA
þB
d
d
d
where w0 is the width of the PSF when a molecule is on the focal
plane and d is the focus depth of the microscope. Peak a of the
PSF depends on the width wðx3 Þ and is modeled as:

aðx3 Þ ¼

a02pwðx3 Þ2:

Fig. 3. Width of the observed fluorescent beads and values of the defocus curve. The
observed widths of the fluorescent beads at each depth are indicated by circles. The
red line shows the value of the defocus curve that approximates the width using theparameters in the Table 1.

ð6Þ

The parameters a0 ; b; d; A, and B are determined as listed in Table 1
using a set of images of fluorescent beads obtained from different
depths. Namely, these parameters are fitted to a small number ofreal observations, so that Eq. (4) is close to the real observations.
Fig. 3 shows the width of the observed fluorescent beads and thevalues of the defocus curve (5).

Further, based on Eq. (7), if Dz is the same for all z ¼ 1; 2; . . . ; k,we cannot distinguish whether the amount of lateral drifts is Dz ortrue location of the molecule is xhj  Dz , even when we use theimages of all planes. Instead, we estimated the relative lateral drift
D0z ¼ Dz  D1 ðz ¼ 2; . . . ; kÞ from a reference plane z ¼ 1. In thisstudy, the amount of relative lateral drift, as well as the molecule
positions, were estimated at a high-resolution voxel level.

2.3. Lateral drift of the focal planes

3. Method

When we used MUM to localize the molecules, we must consider the lateral drift of the focal planes. Otherwise, the estimation
accuracy of the localization becomes deteriorated because theappropriate observation matrix H varies depending on the drift.
Let Dz1 and Dz2 2 Z ðz ¼ 1; 2; . . . ; kÞ be the amount of lateral drift
along the horizontal and vertical axes of the z-th focal plane,
respectively. In this study, we only considered high-resolution
voxel-level lateral drifts. The drift vector is written asDz  ðDz1 ; Dz2 ; 0Þ. The PSF was defined in Eq. (4), and the following^ l ; xh Þ:
equation holds for hðxi3.1. Compressed sensing with lateral drift estimation
When solving the molecule localization problem stated in Eq.
(3), the lateral drift of the focal planes should be considered to estimate the molecule location accurately. The approximated observa^ varies depending on the lateral drift; hence, thetion matrix H^ can be modified so that theapproximated observation matrix H
molecule location is correctly estimated.
^ based on the
By ordering the rows of the observation matrix Hj^ l þ Dz ; xh Þ ¼ hðx^ l ; xh  Dz Þ;
hðx
i
j
i
j

ð7Þ

^
axial position xli3 of the observation coordinate xli ði ¼ 1; 2; . . . nÞ, H
becomes a block matrix

where the observation coordinates xli is on the z-th plane. This equation implies that the observation from the molecule at xhj with a lateral drift of Dz is identical to the observation from the molecule at
xhj  Dz without any lateral drift. As the observation from each plane
is affected by each drift in camera position, the locations of a molecule estimated from the images obtained from each focal plane are
not identical. However, the problem Eq. (3) considers all focal
planes to estimate the location of a molecule. Hence, unless wequantify the amount of lateral drift for each plane, the true molecule position cannot be correctly estimated even when a single
molecule exists in the 3D space.

2^ 3
H1
6^ 7
6 H2 7^ ¼6
7H
6 .. 7;4 . 5^kH^ z 2 Rn0 m represents the observation matrix
where each submatrix Hof the z-th focal plane and n0 ¼ nk. As the lateral drift affects eachplane independently, the drift of each block can be considered individually. Submatrices shifted by D0z ðz ¼ 2; . . . ; kÞ are written as
follows:3^ l  D0 ; xh Þ    hðx^ l  D0 ; xh Þ
hðx
1
1
1
m
z
z
7
6
^ l  D0 ; xh Þ    hðx^ l  D0 ; xh Þ 7
6 hðx
z
z
2
1
2
m 7
^ z ðD0 Þ  6
H
7:
6
z
..
..
7
6
...
5
4
.
.
^ l 0  D0 ; xh Þ
^ l 0  D0 ; xh Þ    hðx
hðx
z
z
n
1
n
m
2

Table 1
PSF parameters.
Parameter

Value

a0

5:00  107
0

b
w0

ð8Þ

1:33  102d

3:02  102

A

7:37  104Here, we consider the maximum amount of shift is given as a
hyperparameter.
282

Neurocomputing 451 (2021) 279–289

T. Aritake, H. Hino, S. Namiki et al.

layer at the last layer, the first convolution layer extracts featuresfrom an input image, followed by the three deconvolution layers.
As we assume that the scaling factor is 8, based on the doubling
the lateral resolution of an input image at each layer, we used threedeconvolution layers. For these layers, we also used a rectified linear unit, also known as ReLU, as our activation function, followedby batch normalization layers [26], to enhance the training speed
and the estimation accuracy for these layers. For the last layer ofthe network, a convolution layer followed by a sigmoid function
is used to obtain a set of images for the target axial resolution.
Unlike a natural image super-resolution task, the localization problem in Eq. (3) aims to estimate the existence of a molecule in eachhigh-resolution voxel. Hence, the network outputs the probabilityof the existence of a molecule in each voxel using a sigmoid function as an activation function. In other words, at the last layer ofthe network, the binary classification problem, where the network
estimates whether a molecule exists or not, is solved for eachvoxel.
The architecture of the proposed network is based on iterative
optimization to solve the compressed sensing problem Eq. (9),and the proposed network can be seen as one of the unrolling
methods [14,15]. When we used iterative optimization to solve
the compressed sensing problem, as in Eq. (9), each iteration
enhances the resolution of the estimated output. Namely, a deconvolution problem is solved and an input image is upsampled anddeblurred at each iteration. Therefore, to mimic the iterative
deconvolution process, we stacked deconvolution layers so thatthe resolution of a low-resolution input image is doubled at eachdeconvolution layer while features for the localization are
extracted. Intuitively, many iterations of optimization algorithm
are replaced by a smaller number of trained layers as unrolling
methods.
The input of the network is the observed low-resolution image
that is possibly affected by the drift in camera positions, and theoutput is the estimated molecule distribution. The network wastrained to minimize the error between the true molecule distribution and the distribution which is estimated from the observed
image. In addition, the network was trained to be robust to the lateral drifts in camera positions so that the trained network can correctly estimate the molecule distribution without explicitly

Now, the problem is to estimate both the lateral drifts and themolecule locations from the observation yl , which can be formulated as:

minimize
k

T 

X
2
0 k
^
kylt  HðfDz gz¼2 Þc t kF þ kkc t k1 ;

ð9Þ

fc t gTt¼1 ;fD0z gz¼2 t¼1

where

2

^1
H

3

7
6^
6 H 2 ðD02 Þ 7
7
6
0 k
^
HðfD
g
Þ
¼
7:
6
z z¼2
..
7
6
.
5
4
0
^
H k ðDk Þ

ð10Þ

Here, we considered T images simultaneously to ensure that theamount of lateral drifts is correctly estimated from the images. As
we assumed that only a small number of molecules exist in the target 3D space, estimating all the lateral drifts from a single image is
not always possible.
However, the optimization problem above has a high computational cost because we need to consider all the possible pairs of lateral drifts to obtain an optimal solution. Although this problem cankbe solved by alternating the optimization of fcgTt¼1 and fD0z gz¼2 toobtain a sub-optimal solution, the optimization of fc t gTt¼1

has anincreasingly high computational cost as the input image size
becomes larger or the target resolution of the molecule localization
becomes higher. Moreover, the optimization with respect to fc t gTt¼1
is solved using an iterative algorithm. The slow computational
speed prevents us from applying 3D SMLM using MUM for interactive observation. Therefore, a faster method for solving the optimization problem is required to obtain a super-resolution image
within a reasonable computational time.
3.2. Convolutional neural network (CNN)
In this study, we used a CNN to reduce the computational costto solve the optimization problem Eq. (9). We assumed that theresolution of the input and output images is given, and the scaling
factor is 8 along each axis and the number of focal planes is k ¼ 4
henceforth. As the other SMLM methods using DNNs [16–20], theCNN was used to reduce the computational cost of the optimization with respect to c t .
We propose a CNN which is composed of several convolution
layers and deconvolution layers, as shown in Fig. 4. The mainbuilding blocks of the network are the deconvolution layers, which
are used in FSRCNN [23] for single image super-resolution of natural images. Unlike FSRCNN, which uses only one deconvolutionkestimating the amount of lateral drifts fD0z gz¼2 . To train the robust
network, artificial data which are affected by the random lateraldrifts are used and parameters of the network are optimized to
localize molecules from images that are affected by the lateraldrifts.
Let p ¼ ðp1 ; p2 ; . . . pm Þ be the ground-truth molecule existence
probability for each voxel, where

Fig. 4. Network architecture of the proposed model.
283

T. Aritake, H. Hino, S. Namiki et al.


pj ¼

Neurocomputing 451 (2021) 279–289

0 cj ¼ 0;
1 cj > 0;ð11Þ

and let q ¼ ðq1 ; q2 ; . . . ; qm Þ be the probability of the existence of amolecule as estimated by the CNN. As the network solves the binaryclassification problem at each voxel, we used the sum of the binary
cross-entropy (BCE):‘ðp; qÞ ¼mX
pj logðqj Þ  ð1  pj Þ logð1  qj Þ ;j¼1

as the loss function to train the network.
To train the network, we need a training dataset. However, thetrue molecule density of the real samples is not known. Therefore,
we used artificial low-resolution images generated from artificial
distributions to train the network.
We generated random molecule distributions that contain K
molecules in a 3D space. In this study, K ¼ 3 was used; the sizeof the target 3D space is 3072 nm  3072 nm  1600 nm, andthe coordinates of the molecules xj ðj ¼ 1; 2; . . . ; KÞ were independently drawn from a uniform distribution on this space. The weightck ðj ¼ 1; 2; . . . ; KÞ of each molecule was also independently drawn
from a continuous uniform distribution on ½0:3; 1:0 , and the target
value p was generated as in Eq. (11).
We assumed that the size of the low-resolution image obtained
from each focal plane is 16  16 and the size of each pixel ofimages is 192 nm  192 nm. The size of the high-resolution voxel
is 24 nm  24 nm  50 nm, and a high-resolution image set is composed of 32 images of size 128  128. The relative lateral drift of
the focal planes D0z1 ; D0z2 ðz ¼ 2; . . . ; 4Þ are randomly chosen as 24s
nm independently, where s is drawn from a discrete uniform distribution on ½2; 2
Z. Then, low-resolution images are generatedPK ^hðx; xh Þ on the low-resolution
by calculating the values of
k¼1

Fig. 5. Optical layout of the quad-plane microscope. The intermediate image isrelayed onto each camera via a pair of lenses (L1, f = 125.0 mm; L2, f = 100.0 mm).
TL, tube lens; M1, 1:1 beam-splitter mirror. The inset shows the focusing planes offour cameras.

a bandpass mirror (672–712 nm, FF01-692/40, Semrock) was used
to separate the excitation and emission light. The fluorescence
image formed by the internal tube lens of the inverted microscope
was relayed by an achromatic lens (f = 125.0 mm, Thorlabs), split
twice by 1:1 beam-splitter mirrors (BSW29R, Thorlabs), and refocused onto the four cameras via achromatic lenses (f= 100.0 mm, Thorlabs). The axial positions of the achromatic lenses
in front of the cameras were adjusted so that the four planes at
400 nm intervals in the Z-axis direction of the specimen corresponding to the conjugate planes of the sensor surface of therespective camera. The relative distance among planes was estimated based on a shift in Z-axis position dependence of PSFs,
which were determined by imaging fluorescent beads (FluoSphere
Carboxylate-Modified Microspheres, 0.2 lm, Invitrogen) whilevarying Z-axis positions of the objective using a piezo positioner(P725.1, PI). The difference in the field of view of the cameras
was corrected by coordinate registration using affine transformation, the parameters of which were determined by images of multiple fluorescent beads captured using different cameras.
Methanol-fixed COS7 cells were used for STORM imaging of
tubulin molecules expressing inside the cells as described previously [27]. The primary and secondary antibodies were an antitubulin antibody (YL1/2, Abcam) and an Alexa Fluor 647-labeledanti-rat IgG antibody, respectively. The specimen was mounted
in a STORM buffer (10 mM NaCl, 60% sucrose, 10% glucose, 0.1%
b-mercaptoethanol, 0.5 mg/mL glucose oxidase, 0.04 mg/mL catalase, and 50 mM HEPES, pH 8.0) and then subjected to imaging.
Images were acquired at 22 Hz with 20 ms exposure.

kgrids x ¼ xli ði ¼ 1; 2; . . . ; nÞ.
To generate a training dataset, the coordinates of the moleculesfxj gKj¼1 and the relative lateral drifts of the focal planes
44

fD0z1 gz¼2 ; fD0z2 gz¼2 were randomly drawn, as above, for each framet ¼ 1; 2; . . . ; T independently. By training the network with thisdataset, it is expected that the trained network becomes robust
to the lateral drifts in camera positions within ½2; 2 highresolution voxels.
4. Experiments
We showed the experimental results of the localization by our
proposed method using both artificial images and the real microscopy images. In the following experiments, images observed by
the above microscope were processed on an Nvidia Tesla V100
32 GB GPU to localize molecules in the images.
4.1. Optical layout of the microscope
In this section, the experimental settings of the quad-plane
microscope used in the experiments were presented. A multifocus microscope equipped with four EM-CCD cameras (iXon
897, Andor) was constructed based on a commercial inverted
microscope (ECLIPSE Ti, Nikon) (Fig. 5). A 640 nm laser beam
(HL6366DG, Thorlabs) that passed through a cleanup filter(LD01-640/8, Semrock) was focused on the back focal plane of a100 oil immersion objective (Plan Apo VC 100X/1.40, Nikon) toilluminate an Alexa Fluor 647-stained specimen at an excitation
intensity of approximately 5 kW/cm2. The fluorescence emitted
from the specimen was collected by the same objective. A filter
cube consisting of an excitation filter (608–648 nm, FF02-628/40,
Semrock), a dichroic mirror (669 nm, FF660-Di02, Semrock), and4.2. Experiments with artificial images
First, we used two different datasets to assess the effects of lateral drifts for molecule localization. One dataset was affected by
the random lateral drifts D0k ¼ ðD0k1 ; D0k2 ; 0Þ ðk ¼ 2; 3; 4Þ for each
image independently, and the other dataset did not consider thelateral drifts when generating the artificial images. Then, we usedeach dataset to train the network. Each dataset consists of 90,000
90000

low-resolution images fylt gt¼1

and the corresponding molecule

generated from the random moleexistence probabilities fpt g90000
t¼1
284

Neurocomputing 451 (2021) 279–289

T. Aritake, H. Hino, S. Namiki et al.

cule distribution. The parameters of the experiments are presented
in Table 2. An example of training data is shown in Fig. 6. We used
Adam [28] as an optimizer, with the parameters b1 ¼ 0:9; b2 ¼ 0:99
and the initial learning rate was set to 1:0  103 , and the batch
size was 100. The dataset was randomly shuffled at the end of eachone of 30 epochs.
To assess the effect of lateral drifts, we evaluated the estimation
accuracy of the two trained networks. We generated 500 artificial
images that contain only one molecule in the 3D space, and theother parameters were the same as the parameters in Table 2. In
this experiment, we assumed that a molecule exists in a voxel
where the network outputs the highest molecule existence probability. Fig. 7 shows a 2D histogram of the estimation error along x
and y axes by two networks. The error is discretized at the resolution of the high-resolution voxels, and the element at (0, 0) shows
the number of correctly localized molecules. Further, the Jaccard
indices of the classification of the voxels that contain a molecule
are shown in the caption. As we can see from the figure, the network that is trained to be robust to the lateral drifts localizes mostof the molecules correctly. Meanwhile, the error of the networkthat is not robust to the lateral drifts distributed to a broader area,
and the number of the correctly localized molecules is less than in
Fig. 7(a). This result suggests that if we do not consider the effect of
the lateral drifts, the accuracy deteriorates even when only one
molecule exists in the target space.
We also evaluated the computational speed and localization
accuracy of the compressed sensing method and robustly trained
network using 500 artificial images that contain only one moleculein the target 3D space. As the previous experiment, the generated
images were affected by the lateral drifts. Also, we assumed thatthe number of molecules in the target space is known for the compressed sensing method. Table 3 shows the average lateral andaxial error between the estimated voxel and the true voxel. The lateral error of the CNN is smaller than the lateral error of the compressed sensing. On the other hand, the axial error of thecompressed sensing method is smaller than the axial error of
CNN, however, this is due to the assumption that the number ofmolecules is known in advance. Most importantly, when estimating the molecule distribution, CNN process shows much higher
fps than the compressed sensing method.
Next, we validated the accuracy of the robustly trained network
against lateral drifts. As the previous experiment, we generated
artificial images that contain only one molecule in the 3D spaceand estimated that a molecule exists in the voxel where the network outputs the highest molecule-existence probability. In
Fig. 8, the mean localization accuracy along the horizontal (X), vertical (Y), and axial (Z) directions with 95% confidence intervals areshown for each true molecule depth. The figure indicates that theerror along each axis is within a high-resolution voxel on averagealong each axis at all depths.
Fig. 9 shows the estimation results with multiple molecules.
The molecules are sampled from the helix curve (red line), andtheir high-resolution coordinates are shown by red circles. In thisexperiment, we estimated that the molecules exist in the voxels

Fig. 6. (a) An example of artificial molecule distribution and a part of the target 3D
space is shown for visualization purpose. True locations of molecules are shown by
red circle and the voxels which contain the molecules are shown by green cuboid.
(b) The observed image generated from the artificial distribution. Note that theimages are normalized for visual emphasis.

whose molecule-existence probability exceeds a certain threshold.
The thresholding value is common for all locations and needs to be
specified by a user. Here, we selected 0.1 for the thresholding
value. As demonstrated in Fig. 9(a), three molecules distributedin the 3D space are also accurately localized. Although, as Fig. 9(b) indicates, the closely located molecules are difficult to localize,
the estimated locations are close to the ground-truth locations. Byplotting all the detected molecules from all frames, the helix curve
structure behind the molecules can be observed in Fig. 9(c).
The processing speed of the network is presented in Fig. 10. As
shown in the figure, the computational speed decreases as the sizeof the image increases and is inversely proportional to the numberof pixels of an image. Because 22 images are obtained by our
microscopy at every second, further improvement of the processing speed is required to process the large images in real-time. Still,
the computational time is significantly reduced in comparison to
the compressed sensing method. By solving the problem Eq. (9)by an alternating minimization, the processing speed is only

Table 2
Parameters for experiments with artificial data.
Parameter

Value

Target space size
Resolution of a low-resolution voxel
Resolution of a high-resolution voxel
Size of observed images
Number of molecules
Variance of Gaussian noise
Rate parameter of Poisson noise

3072  3072  1,600 [nm]192  192  400 [nm]24  24 50 [nm]16  16 4
3
r2 ¼ 9k¼4

1:45  103 fps even for a small 16  16  4 input.
285

T. Aritake, H. Hino, S. Namiki et al.

Neurocomputing 451 (2021) 279–289

Fig. 7. Histogram of the localization error of a single molecule. (a) the network that is trained to be robust against the lateral drift and (b) the network that is trained without
considering the lateral drifts.

4.3. Experiments with real images

Table 3
Computational speed and accuracy of the CNN and compressed sensing method.

Lateral error [nm]Axial error [nm]Frame per second

CNN

Compressed sensing

16.7
8.90
125

57.6
0.400
0.129

In this section, we show the experimental results with real datathat observed microtubules by our microscopy. In this experiment,
there were no ground-truth results; hence, we used the trainedneural network from the previous subsection to localized
molecules.
The resolution of the low- and high-resolution images are thesame as in the previous subsection. The size of the input images
for each frame is 256  256  4, and the target image size is2048  2048  32. An example of a frame is shown in Fig. 11.
The dataset contains 39,000 frames of images, and each frame
was processed independently to localize the molecules. In thisexperiment, we estimated that a molecule exists in a voxel if themolecule-existence probability exceeds 0:05.
Fig. 12 shows the estimated high-resolution image generated bymerging localization results of all the frames. Each pixel of theimage is a binary value, which indicates that the voxel contains amolecule in more than one frame. From the figures, we can observea tubular structure of the microtubules, which varies depending on
the depth.

5. Discussion
This study presents the 3D molecule localization problem using
quad-plane microscopy. The problem with using MUM is lateraldrift in the camera positions, which makes the localization lessFig. 8. Average estimation error along the horizontal (X), vertical (Y), and axial (Z)axes with their 95% confidence intervals.

Fig. 9. Localization results using artificial data. The red lines shows the true molecular structure from which molecules were sampled from. The red circles show true highresolution molecule coordinates and the blue triangles show a estimated molecule locations. Figures (a) and (b) shows the estimation results of selected frames, and figure (c)shows the image reconstructed from 300 frames.
286

Neurocomputing 451 (2021) 279–289

T. Aritake, H. Hino, S. Namiki et al.

Fig. 10. Computational speed (fps) of estimations by the trained network forimages of the size 16 px  16 px, 32 px  32px, . . . ,256 px  256 px. The estimation
speed of our network is inversely proportional to the number of pixels of an image.

Fig. 12. Estimated high-resolution image of the microtubules data. The depthdependent tubular structure of the sample is visualized by the colors.

the molecule location without estimating the amount of lateraldrifts. We expect that this technique can be used to broaden theapplicability of MUM to 3D imaging because an explicit drift correction is not required.
However, some limitations are worth noting. Although our proposed method significantly reduces the computational cost of solving the localization problem, large images are still difficult toprocess with real-time processing speed. Therefore, future work
should include further improvement in computational speed.
Using a faster method to extract possible molecule locations andlocalizing the molecule by the proposed method may further
improve the computational efficiency.
CRediT authorship contribution statement
Toshimitsu Aritake: Conceptualization, Methodology, Software, Validation, Formal analysis, Writing - original draft, Visualization. Hideitsu Hino: Formal analysis, Conceptualization.
Shigeyuki Namiki: Investigation. Daisuke Asanuma: Investigation. Kenzo Hirose: Supervision. Noboru Murata: Project
administration.

Fig. 11. A frame which observed and microtubules. The observed images is a
256  256  4 image and the difference in the field of view of the cameras wascorrected by affine transformation.

Declaration of Competing Interest
accurate. We formulated the localization problem as a compressedsensing problem that consists of the molecule localization and an
estimation of the amount of drifts. However, the computational
cost to solve this problem is high, and the optimal solution cannot
be obtained within a reasonable computational time. In this study,
we used CNN to solve this problem accurately and efficiently. The
network is trained to be robust against the sub-pixel-sized lateraldrift in the camera locations. Therefore, the trained network canaccurately estimate the locations of molecules in an observed
image without explicitly estimating the amount of drift in thecameras.
Experiments with both artificial data and real data were presented. The results suggest that the network achieves 3D localization of the molecules with an average lateral resolution of 25 nm
and an axial resolution of 50 nm. The trained network is also
robust to the lateral drift in the camera positions and estimates

The authors declare the following financial interests/personal
relationships which may be considered as potential competinginterests: Dr. Hideitsu Hino reports grants from Japan Society for
the Promotion of Science KAKENHI Grant Numbers 17H01793 andgrants from Japan Science and Technology Agency CREST Grant
Number JPMJCR1761 during the conduct of the study.
Dr. Noboru Murata reports grants from Japan Society for the Promotion of Science KAKENHI Grant Numbers 17H01793 , 18H03291
and grants from Japan Science and Technology Agency CREST Grant.
Acknowledgments
The authors would like to thank M. Tanaka for technical assistance. This work was partially supported by JSPS KAKENHI (Grant
Nos. 17H01793 and 18H03291) and JST CREST (Grant Nos.
JPMJCR1761 and JPMJCR14D7).
287

T. Aritake, H. Hino, S. Namiki et al.

Neurocomputing 451 (2021) 279–289[25] L. Gu, Y. Sheng, Y. Chen, H. Chang, Y. Zhang, P. Lv, W. Ji, T. Xu, High-density 3D
single molecular analysis based on compressed sensing, Biophysical Journal
106 (11) (2014) 2443–2449, https://doi.org/10.1016/j.bpj.2014.04.021.
[26] S. Ioffe, C. Szegedy, Batch normalization: Accelerating deep network training
by reducing internal covariate shift, in: Proceedings of the 32nd International
Conference on Machine Learning - vol. 37, ICML’15, JMLR.org, 2015, pp. 448–456.
[27] D.W. Cleveland, K.F. Sullivan, Molecular biology and genetics of tubulin,
Annual Review of Biochemistry 54 (1) (1985) 331–366, https://doi.org/
10.1146/annurev.bi.54.070185.001555, pMID: 3896122.
[28] D.P. Kingma, J. Ba, Adam: A method for stochastic optimization, in: 3rdInternational Conference on Learning Representations, ICLR 2015, San Diego,
CA, USA, May 7–9, 2015, Conference Track Proceedings, 2015.

References[1] L. Schermelleh, A. Ferrand, T. Huser, C. Eggeling, M. Sauer, O. Biehlmaier, G.P.
Drummen, Super-resolution microscopy demystified (jan 2019). doi:10.1038/
s41556-018-0251-8.
[2] E. Betzig, G.H. Patterson, R. Sougrat, O.W. Lindwasser, S. Olenych, J.S.
Bonifacino, M.W. Davidson, J. Lippincott-Schwartz, H.F. Hess, Imaging
intracellular fluorescent proteins at nanometer resolution, Science 313(5793) (2006) 1642–1645, https://doi.org/10.1126/science.1127344.
[3] S.T. Hess, T.P. Girirajan, M.D. Mason, Ultra-high resolution imaging byfluorescence photoactivation localization microscopy, Biophysical Journal 91(11) (2006) 4258–4272, https://doi.org/10.1529/biophysj.106.091116.
[4] W. Liu, K.C. Toussaint, C. Okoro, D. Zhu, Y. Chen, C. Kuang, X. Liu, Breaking theaxial diffraction limit: a guide to axial super-resolution fluorescencemicroscopy, Laser and Photonics Reviews 12 (8) (2018) 1–29, https://doi.org/
10.1002/lpor.201700333.
[5] B. Huang, W. Wang, M. Bates, X. Zhuang, Three-dimensional super-resolution
imaging by stochastic optical reconstruction microscopy, Science 319 (5864)(2008) 810–813, https://doi.org/10.1126/science.1153529.
[6] S.R.P. Pavani, M.A. Thompson, J.S. Biteen, S.J. Lord, N. Liu, R.J. Twieg, R. Piestun,
W.E. 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 106 (9) (2009) 2995–2999,
https://doi.org/10.1073/pnas.0900245106.
[7] Y. Shechtman, S.J. Sahl, A.S. Backer, W.E. Moerner, Optimal point spreadfunction design for 3d imaging, Physical Review Letters 113 (2014) 133902,
https://doi.org/10.1103/PhysRevLett.113.133902.
[8] S. Ram, P. Prabhat, J. Chao, E.S. Ward, R.J. Ober, High accuracy 3D quantum dot
tracking with multifocal plane microscopy for the study of fast intracellular
dynamics in live cells, Biophysical Journal 95 (12) (2008) 6025–6043, https://doi.org/10.1529/biophysj.108.140392.
[9] H.P. Babcock, Multiplane and spectrally-resolved single molecule localization
microscopy with industrial grade CMOS cameras, Scientific Reports 8 (1)(2018) 1726, https://doi.org/10.1038/s41598-018-19981-z.
[10] M.J. Mlodzianoski, J.M. Schreiner, S.P. Callahan, K. Smolková, A. Dlasková, J.
Šantorová, P. Ježek, J. Bewersdorf, Sample drift correction in 3d fluorescencephotoactivation localization microscopy, Optics Express 19 (16) (2011)
15009–15019, https://doi.org/10.1364/OE.19.015009.
[11] Y. Wang, J. Schnitzbauer, Z. Hu, X. Li, Y. Cheng, Z.-L. Huang, B. Huang,
Localization events-based sample drift correction for localization microscopy
with redundant cross-correlation algorithm, Optics Express 22 (13) (2014)
15982–15991, https://doi.org/10.1364/OE.22.015982.
[12] D.L. Donoho, Compressed sensing, IEEE Transactions on Information Theory 52(4) (2006) 1289–1306, https://doi.org/10.1109/TIT.2006.871582.
[13] A. Lucas, M. Iliadis, R. Molina, A.K. Katsaggelos, Using deep neural networks forinverse problems in imaging: Beyond analytical methods, IEEE Signal
Processing Magazine 35 (1) (2018) 20–36, https://doi.org/10.1109/
MSP.2017.2760358.
[14] K. Gregor, Y. LeCun, Learning fast approximations of sparse coding, in:Proceedings of the 27th International Conference on Machine Learning,
ICML’10, Omnipress, Madison, WI, USA, 2010, p. 399–406.
[15] J.R. Hershey, J.L. Roux, F. Weninger, Deep unfolding: Model-based inspirationof novel deep architectures, ArXiv abs/1409.2574.
[16] P. Zelger, K. Kaser, B. Rossboth, L. Velas, G.J. Schütz, A. Jesacher, Threedimensional localization microscopy using deep learning, Optics Express 26(25) (2018) 33166–33179, https://doi.org/10.1364/OE.26.033166.
[17] W. Ouyang, A. Aristov, M. Lelek, X. Hao, C. Zimmer, Deep learning massively
accelerates super-resolution localization microscopy, Nature Biotechnology 36(5) (2018) 460–468, https://doi.org/10.1038/nbt.4106.
[18] E. Nehme, L.E. Weiss, T. Michaeli, Y. Shechtman, Deep-storm: super-resolution
single-molecule microscopy by deep learning, Optica 5 (4) (2018) 458–464,
https://doi.org/10.1364/OPTICA.5.000458.
[19] N. Boyd, E. Jonas, H. Babcock, B. Recht, Deeploco: Fast 3d localization
microscopy using neural networks, bioRxivdoi:10.1101/267096.
[20] P. Zhang, S. Liu, A. Chaurasia, D. Ma, M.J. Mlodzianoski, E. Culurciello, F. Huang,
Analyzing complex single-molecule emission patterns with deep learning,
Nature Methods 15 (11) (2018) 913–916, https://doi.org/10.1038/s41592018-0153-5.
[21] E. Hershko, L.E. Weiss, T. Michaeli, Y. Shechtman, Multicolor localization
microscopy and point-spread-function engineering by deep learning, Optics
Express 27 (5) (2019) 6158–6183, https://doi.org/10.1364/OE.27.006158.
[22] L. Möckl, A.R. Roy, P.N. Petrov, W.E. Moerner, Accurate and rapid backgroundestimation in single-molecule localization microscopy using the deep neuralnetwork bgnet, Proceedings of the National Academy of Sciences 117 (1)(2020) 60–67, https://doi.org/10.1073/pnas.1916219117.
[23] C. Dong, C.C. Loy, X. Tang, Accelerating the super-resolution convolutional
neural network, in: B. Leibe, J. Matas, N. Sebe, M. Welling (Eds.), Computer
Vision – ECCV 2016, Springer International Publishing, Cham, 2016, pp. 391–
407.
[24] T. Hastie, R. Tibshirani, M. Wainwright, Statistical Learning with Sparsity: The
Lasso and Generalizations, Chapman & Hall/CRC, 2015.

Toshimitsu Aritake received the B.Eng. and M.Eng.
from Waseda University in 2012 and 2014. He joined
Hitachi’s Central Research Laboratory and worked as a
research staff from April 2014 to September 2017. Since
October 2017, he has been working toward the Doctor’s
degree at Waseda University. His research interests
include machine learning such as sparse modeling,
matrix factorization and deep learning.

Hideitsu Hino received his bachelor’s degree in engineering in 2003, and a master’s degree in Applied
Mathematics and Physics in 2005 from Kyoto University. He was with Hitachi’s Systems Development Laboratory from April 2005 to August 2007. He earned D.
Eng. in 2010 from Waseda University. From April 2013
to 2018, he was an Assistant and Associate Professor at
University of Tsukuba, and he is currently a Professor at
The Institute of Statistical Mathematics. His researchinterests include analysis of learning algorithms and theapplication of statistical methods to various problems.

Shigeyuki Namiki received the B.Sc. and M.Sc. from the
University of Tsukuba and received a Ph.D. from the
University of Tokyo in 2002. He was an assistant professor in 2007 in the Department of Cell Physiology at
Nagoya University. He was an assistant professor in2009 and is currently a lecturer in the Department of
Neurobiology at the University of Tokyo. His researchinterest is the relationship between the dynamics ofsynaptic molecules and synaptic plasticity.

Daisuke Asanuma is an Assistant Professor in the
Department of Pharmacology at the University of Tokyo
since 2018. He is also a Japan Science and Technology
Agency (JST)-PRESTO researcher since 2017. From 2011 to
2018, he served as an Assistant Professor in the Department of Neurobiology at the University of Tokyo. He
received his Ph.D. in 2011 from the University of Tokyo. He
was a Research Fellow of the Japan Society for the Promotion of Science (JSPS) from 2008 to 2011. He has been
involved in development of chemical tools for biologicaland medical applications such as fluorescence imaging of
bioactive molecules and in vivo tumor imaging.

288

Neurocomputing 451 (2021) 279–289

T. Aritake, H. Hino, S. Namiki et al.
Kenzo Hirose is a Professor in the Department ofPharmacology at the University of Tokyo since 2016.
From 2008 to 2018 he served as a professor in the
Department of Neurobiology. From 2005 to 2008 heserved as a professor in the Department of Cell Physiology at Nagoya University. He earned his M.D. in 1992
and received his Ph.D. in 1996 in the University ofTokyo. He is interested in the cellular signaling and has
been developing new technologies to visualize signaling
molecules in living cells.

Noboru Murata received the B. Eng., M. Eng. and Dr.
Eng. in Mathematical Engineering and Information
Physics from the University of Tokyo in 1987, 1989 and1992, respectively. After working at the University ofTokyo, GMD FIRST in Germany, and RIKEN in Japan, in
April of 2000, he joined Waseda University in Japan
where he is presently a professor. His research interest
includes the theoretical aspects of learning machines
such as neural networks, focusing on the dynamics andstatistical properties of learning.

289