# Computational multifocal microscopy

Kuan He  1,* Zihao Wang  1 Xiang Huang  2 Xiaolei Wang  3 Seunghwan Yoo  1
Pablo Ruiz  1 Itay Gdor  3 Alan Selewa  3 Nicola J. Ferrier  2 Norbert Scherer
3,4 Mark Hereld  2 Aggelos K. Katsaggelos  1 and Oliver Cossairt1  1Department
of Electrical Engineering and Computer Science, Northwestern University,
Evanston, IL 60208, USA
2Mathematics and Computer Science, Argonne National Laboratory, 9700 South
Cass Avenue, Lemont, IL 60439, USA
3Department of Chemistry, The University of Chicago, 5801 South Ellis Avenue,
Chicago, IL 60637, USA
4James Franck Institute, The University of Chicago, 5801 South Ellis Avenue,
Chicago, IL 60637, USA<EMAIL_ADDRESS>
###### Abstract

Despite recent advances, high performance single-shot 3D microscopy remains an
elusive task. By introducing designed diffractive optical elements (DOEs), one
is capable of converting a microscope into a 3D "kaleidoscope", in which case
the snapshot image consists of an array of tiles and each tile focuses on
different depths. However, the acquired multifocal microscopic (MFM) image
suffers from multiple sources of degradation, which prevents MFM from further
applications. We propose a unifying computational framework which simplifies
the imaging system and achieves 3D reconstruction via computation. Our optical
configuration omits chromatic correction grating and redesigns the multifocal
grating to enlarge the tracking area. Our proposed setup features only one
single grating in addition to a regular microscope. The aberration correction,
along with Poisson and background denoising, are incorporated in our
deconvolution-based fully-automated algorithm, which requires no empirical
parameter-tuning. In experiments, we achieve the spatial resolutions of
$0.35$um (lateral) and $0.5$um (axial), which are comparable to the resolution
that can be achieved with confocal deconvolution microscopy. We demonstrate a
3D video of moving bacteria recorded at $25$ frames per second using our
proposed computational multifocal microscopy technique.

††journal: boe††articletype: Research Article

##  1 Introduction

Typically, an optical microscope focuses on one focal plane at a time. Thus,
recovering 3D information with sufficient resolution usually requires
sequentially z-scanning the focal planes. This process is time consuming and
requires precise mechanical control. The long acquisition time of focal
scanning microscopy fundamentally prevents the applications in _in vivo_
imaging, especially when the 3D movement of biomedical objects is desired. To
realize _single-shot_ 3D microscopy, a standard microscope has to be modified
in such a way that 3D information can be encoded onto a 2D plane.

One of the few modifications to achieve single-shot 3D microscopy, for
example, is by placing a diffractive optical element (DOE) in the Fourier
plane of a standard microscope (Fig. 1a). The DOE is designed as a distorted
phase gratings [1, 2, 3], also referred to as a multifocal gratings (MFG). The
MFG projects different depth layers in a 3D volume onto different sub-regions
of the imaging plane simultaneously. In particular, if the imaging plane is
divided into $l\times l$ tiles and each tile is focused on one focal plane by
design, then a 3D volume of $l^{2}$ depths can be retrieved by z-stacking all
the tiles. An $l=3$ case is shown in Fig. 1(c). This imaging technique is also
referred as multifocal microscopy (MFM) [4].

Figure 1: Single-DOE multifocal microscopy (MFM) setup (a) and computational
3D reconstruction pipeline (d-f). In (a), a conventional microscope is
augmented by a $4f$ system. An MFG (b) is inserted at the Fourier plane of the
$4f$ system to produce an array of $l\times l$ differently focused tile images
in a single exposure (c; $l=3$). Note that our CMFM system discards CA
corrective optics, significantly reducing the system complexity and cost
compared to conventional MFM. We correct for CA computationally rather than
optically. (d-f) The pipeline of proposed computational framework: by
capturing a z-stack 3D PSFs (d), the algorithm can simultaneously recover the
background noise $b^{*}$, the optimal regularizer parameter $\lambda^{*}$ and
a high resolution 3D image (f) from a single captured 2D MFM image (e).

MFM stands out among other single-shot 3D microscopy, e.g. light fields
microscopy [5, 6] and lensless 3D coded microscopy [7, 8, 9], for its
capability of achieving diffraction-limited 3D resolution comparable to
conventional focal scanning microscopy. However, as an imaging system, MFM
faces several imaging issues, some of which are even more severe that prevents
MFM from broader applications:

  * •

Out-of-focus blur. MFM captures an array of multifocal planes. Like standard
wide-field microscopy, each focal plane image is convolved by out-of-focus
light. Estimating and removing out-of-focus blur requires implementation of a
3D deconvolution algorithm. _Our proposed algorithm builds on the regularized
Rechardson-Lucy (RL) deconvolution, and features automatic estimation of the
background noise and the optimal regularizer parameter._

  * •

Chromatic aberration (CA). Conventionally, CA is optically corrected by adding
a chromatic correction gratings (CCG) and a multifacet prism after the MFG
[4], resulting in excessive hardware cost. In our work, we propose a
simplified version by abandoning the CCG and the prism. Instead, we place a
$10$ nm bandpass filter in front of the detector to mitigate CA, as shown in
Fig. 1(a). We demonstrate that the resolution loss due to CA can be
_computationally_ compensated. This comes from a fact that CA from MFG is
directional due to tile geometry. In each tile, the PSF is stretched mainly in
one direction, but preserves imaging quality in the orthogonal direction.
_Therefore, by jointly incorporating all the PSFs in each tile, our proposed
3D deconvolution algorithm enables restoring the resolution loss due to CA._

  * •

Poisson noise due to limited photon budget. There are three factors contribute
to the limited photon budget of the MFM. First, as mentioned before, a
narrower bandpass filter of $10$nm is used to mitigate CA, and thus, the total
number of emission photons arriving at the detector is smaller than when using
CA corrective optics. Second, the theoretical maximum efficiency of the MFG is
smaller than one. For example, the theoretical limits are $68\%$ for a binary
phase-only ($0$ or $\pi$) MFG with $3\times 3$ diffraction orders, and $78\%$
for a MFG with $5\times 5$ diffraction orders, which means $32\%$ or $22\%$ of
the total emission light will be lost and cannot be collected by the detector.
The efficiency can be improved by using a multi-phase MFG [10] at the expense
of the complicate fabrication process. Third, each tile only shares about
$1/l^{2}$ of the total photons. The limited photon budget requires a noise
modeling of Poisson process. _We incorporate the Poisson noise modeling for
MFM, for the first time._

  * •

Field-of-view (FOV). As described above, MFM trades lateral FOV for depth
information. Conventional MFM designs use large tile spacings for large
objects. We address here that this design is not optimal for tracking dynamic
objects, which are usually small in size but require large FOV. We analyze the
FOV performance theoretically, and demonstrate that _by modifying the MFG
design to create a smaller tile spacing, small objects (e.g. bacteria) can be
tracked over a larger lateral area than the tile width._

###  1.1 Related works

The 2D sensing limitation imposed by optical sensors makes 3D imaging an
interesting and active research topic both in macroscopic and microscopic
regime. In microscopy, several approaches have been proposed for a variety of
imaging conditions. One scheme for achieving 3D imaging is to use coherent
illumination, e.g. digital holographic microscopy [11, 12, 13], variably-
controlled light-emitting diodes (LED) arrays [14, 15], etc. Yet the coherent
modeling of light propagation does not apply to incoherent/fluorescent
objects. For the incoherent case, one seminal idea is to introduce self-
interference by projecting a set of Fresnel patterns [16, 17, 18]. However,
most of the existing methods, especially for incoherent cases, require
multiple exposures and massive processing time. This limitation prevents
broader applications such as _in vivo_ imaging as the motion of the objects
during capture process is hard to circumvent and thus, deteriorates the image
quality [19, 20]. Recent works have proposed to reduce the measuring
requirement by computationally exploiting spatial-temporal redundancy of the
scenes [21, 22]. Here, we review two types of computational microscopy that
enable single-shot 3D fluorescence imaging.

_Light field microscopy (LFM)._ By introducing a microlens array on the
primary image plane of a microscopy, LFM encodes 4D (spatial-angular)
information into a 2D image and computationally recover the 3D objects [5, 6].
Because LFM trades spatial resolution for single-shot capture, its spatial
resolution is lower than that of the conventional microscope. Recently, a RL
deconvolution method based on the wave optics theory [23] is derived and
demonstrated to improve the resolution for LFM. So far, the best resolutions
experimentally achieved by LFM are $\sim 1.4$um and $2.6$um [5] in the lateral
and axial dimensions respectively. In addition, due to the variation on
sampling density, the lateral resolution decreases to $\sim 3.75$um at the
focal plane [5].

_Lensless 3D microscopy._ Traditionally, lens-based cameras/microscopes map a
point in the scene to a pixel on the detector. Lensless imaging architectures,
instead, replace the lens with other encoding elements such that a point is
mapped to many points on the detector and thus, require computation to
recover. Several seminal literatures have proposed to place a single encoding
element, such as a coded mask [7] or diffuser [9], directly in front of a
detector. Therefore, the imaging system is compact and cost-effective, and has
large FOV. A computational algorithm is then designed to make use of the
physical effects, e.g. the Point Spread Functions (PSF), to recover the 3D
information of the scene. However, the spatial resolution of lensless 3D
microscopes is restricted to the pixel pitch of the detector. In particular, a
lens-free fluorescent imaging platform was reported to achieve the spatial
resolution of $\sim 10$um based on a compressive sensing algorithm for sparse
objects [8]. The recovered 3D volume consisted of two or three depth layers
with intervals of $50$um or $100$um.

###  1.2 Computational multifocal microscopy

Our computational imaging framework, which we call computational multifocal
microscopy (CMFM), balances and optimizes the processing capabilities of
optics and computation. We demonstrate a CMFM system [Fig.1(a)] that utilizes
simplified optics, as well as algorithms that correct for CA and low-photon
counts. The pipeline of our computational framework is shown in Fig.1(d-f).

We perform two experiments to demonstrate the effectiveness of our CMFM
system. First, we capture a static image of several frozen periplasms in 3D to
demonstrate the spatial resolutions of $0.35$um and $0.5$um in the lateral and
axial dimensions [see Fig. 7(d)], which are comparable to those achievable
with confocal deconvolution microscopy [see Fig. 7(b)]. Note that we compare
$0.5$s captures with our CMFM instrument to a $20$s confocal scan taken with a
dual spinning disk confocal microscope (Model: CSU-W1) made by Yokogawa
Electric Corporation. Our CMFM results show similar 3D image quality, but
achieve a 40x reduction in acquisition time. _However, we kindly remind
readers that the principle of the axial resolution improvement is different
for the two techniques; Further discussion will appear in section 4.1_.
Second, we record a video of _in vivo_ bacterium at $25$ frames per second
(see Visualization 1) and perform high resolution 3D reconstruction with
tracking results (see Visualization 2) using our CMFM technique.

##  2 Methods

###  2.1 Image formation model

The image formation of our CMFM fluorescence system can be modeled as the
axial integral of the 2D convolution of the PSF with its corresponding depth
layer. In the noiseless case, such model can be written as:

$g(x,y)=\int_{-\infty}^{\infty}o_{z}(x,y)*h_{z}(x,y)dz,$ (1)

where $g$ is the observed image, $o$ is the 3D object we want to recover, and
$h$ is a z-stack 3D PSF. $o_{z}$ and $h_{z}$ are 2D slices of $o$ and $h$ at
axial depth $z$, respectively.

Here, we discretize $o$ into $N_{x}\times N_{y}\times N_{z}$ voxels in three
dimensions. Each voxel has size $\Delta_{x}\times\Delta_{y}\times\Delta_{z}$.
The 2D image $g$ on the detector is sampled with $M_{x}\times M_{y}$ pixels.
Each pixel has size $\Delta^{\prime}_{x}\times\Delta^{\prime}_{y}$. For a
band-limited system, the Nyquist sampling rate has to be satisfied in order to
avoid aliasing. In MFM, the lateral cut-off frequency is $f_{c}=2NA/\lambda$,
and therefore $\Delta^{\prime}_{x}\leq 1/(2f_{c})$ and
$\Delta^{\prime}_{y}\leq 1/(2f_{c})$. For simplicity, we set
$\Delta_{x}=\Delta^{\prime}_{x}$ and $\Delta_{y}=\Delta^{\prime}_{y}$ in
practice. In order to match $o$ and $g$ dimensions, $h$ is therefore of
$M_{x}\times M_{y}\times N_{z}$ voxels, with each voxel size of
$\Delta^{\prime}_{x}\times\Delta^{\prime}_{y}\times\Delta_{z}$ . Then the
discrete form of Eq. (1) can be written as a matrix-vector multiplication
form:

$\mathbf{g}=[\mathbf{H}_{1},\mathbf{\cdots},\mathbf{H}_{N_{z}}][\mathbf{o}_{1},\mathbf{\cdots},\mathbf{o}_{N_{z}}]^{T}=\mathbf{H}\mathbf{o}$
(2)

where $\mathbf{g}$ is an $M\times 1$ column vector, in which $M=M_{x}\times
M_{y}$, $\mathbf{o}$ is an $N\times 1$ column vector, in which $N=N_{x}\times
N_{y}\times N_{z}$, and $\mathbf{H}$ is the sensing matrix of a size $M\times
N$. Note that each $\mathbf{H}_{z}$ is a Toeplitz matrix representing 2D
convolution, and can be constructed from $h_{z}$.

Note that, in our case, two types of noise are considered. First, our imaging
process has Poisson noise due to the limited photon budget resulted from (1)
the narrower bandpass filter, (2) the lower diffraction efficiency of the MFG
and (3) splitting of light. Second, the background noise (room stray light, or
the sample itself) is also considered. We assume a uniform background photon
noise across all the pixels in the acquired MFM image. Therefore, the image
formation model under the Poisson noise and additive background noise model
can be expressed as

$\mathbf{g}=\mathcal{P}\left\\{\mathbf{H}\mathbf{o}+\mathbf{b}\right\\},$ (3)

where $\mathcal{P}$ represents Poisson statistics originated from signal
photons, and $\mathbf{b}$ models the uniform background noise. Note that
$\mathbf{b}$ is an $M\times 1$ column vector and each entry of $\mathbf{b}$ is
a same constant, denoted by $b$.

###  2.2 Joint regularized RL deconvolution

Equation (3) leads to the following likelihood function according to Poisson
distribution

$p(\mathbf{g}|\mathbf{o},\mathbf{b})=\prod_{i=1}^{M}\frac{{\left(\mathbf{H}\mathbf{o}+\mathbf{b}\right)_{i}^{\mathbf{g}_{i}}}e^{-\left(\mathbf{H}\mathbf{o}+\mathbf{b}\right)_{i}}}{\mathbf{g}_{i}!},$
(4)

where $i$ stands for the pixel coordinate in $\mathbf{g}$. In RL
deconvolution, the optimal solution $\mathbf{o}^{*}$ is found from the
observation $\mathbf{g}$ by maximizing Eq. (4), or equivalently minimizing its
negative logarithm, subject to all the voxels of the restored image have non-
negative values, that is,

$\mathbf{o}^{*}=\arg\;\max_{\mathbf{o}\geq\mathbf{0}}p(\mathbf{g}|\mathbf{o},\mathbf{b})=\arg\;\min_{\mathbf{o}\geq\mathbf{0}}f(\mathbf{o},\mathbf{b}),$
(5)

where

$f(\mathbf{o},\mathbf{b})=\sum_{i=1}^{M}\left[\left(\mathbf{H}\mathbf{o}+\mathbf{b}\right)_{i}-\mathbf{g}_{i}\log{\left(\mathbf{H}\mathbf{o}+\mathbf{b}\right)_{i}}\right]$
(6)

is the negative Poisson log-likelihood of
$p(\mathbf{g}|\mathbf{o},\mathbf{b})$ in Eq. (4), in which a constant term
$\log{\left(\mathbf{g}_{i}!\right)}$ is omitted.

In regularized RL algorithm, a regularizer term is added to the objective
function in Eq. (5). Here, we consider total variation (TV) regularization
because it can avoid the noise amplification problem in RL deconvolution by
allowing to recover a smooth and stable solution with sharp edges. The
objective function in the TV regularized RL algorithm is therefore written as

$\Phi(\mathbf{o},\mathbf{b},\lambda)=f(\mathbf{o},\mathbf{b})+\lambda\text{TV}\left(\mathbf{o}\right),$
(7)

where $\lambda$ is the regularization parameter and
$\text{TV}(\mathbf{o})=\sum_{j=1}^{N}{\lvert\mathbf{\nabla o}_{j}\rvert}$, in
which $j$ denotes the voxel coordinate in the $\mathbf{o}$. Our joint
algorithm estimates the 3D image, the background noise and the regularization
parameter simultaneously by minimizing $\Phi$ in Eq. (7) subject to all the
voxels of the restored 3D image have non-negative values, that is,

$\\{\mathbf{o}^{*},\mathbf{b}^{*},\lambda^{*}\\}=\arg\;\min_{\mathbf{o}\geq\mathbf{0},\mathbf{b},\lambda}\Phi(\mathbf{o},\mathbf{b},\lambda).$
(8)

The minimization of $\Phi(\mathbf{o},\mathbf{b},\lambda)$ in Eq. (8) with
respect to all three unknown variables can be performed by an iterative
alternating gradient descent method. More specifically, at each iteration, the
three variables are updated sequentially, and when one variable is updated,
the other two variables are fixed as constants.

####  2.2.1 Estimation of 3D image

To update $\mathbf{o}$ for $(k+1)$-th iteration, we substitute
$\mathbf{o}^{k},\mathbf{b}^{k}$ and $\lambda^{k}$ estimated from $k$-th
iteration in Eq. (7), and take partial derivatives with respect to each voxel
$\mathbf{o}^{k}_{j}$ in $\mathbf{o}^{k}$

$\frac{\partial\Phi(\mathbf{o}^{k};\mathbf{b}^{k},\lambda^{k})}{\partial\mathbf{o}^{k}_{j}}=\sum_{i=1}^{M}h_{i,j}-\left[\mathbf{H}^{T}\frac{\mathbf{g}}{\mathbf{H}\mathbf{o}^{k}+\mathbf{b}^{k}}\right]_{j}-\lambda^{k}div\left(\frac{\mathbf{\nabla
o}^{k}_{j}}{\lvert\mathbf{\nabla o}^{k}_{j}\rvert}\right),$ (9)

where the divisions are element wise, $h_{i,j}$ is the entry in the $i$-th row
and the $j$-th column of $\mathbf{H}$, and $\sum_{i=1}^{M}h_{i,j}=1$ if
$\mathbf{H}$ is column normalized (which is equivalent to normalizing PSFs
$h_{z}$), and $div$ stands for the divergence operator. Given the partial
derivative in Eq. (9), We update $\mathbf{o}$ using the gradient-based
algorithm (or equivalently by using EM algorithm), defined by [24, 25]

$\displaystyle\mathbf{o}^{k+1}_{j}=\left[\mathbf{H}^{T}\frac{\mathbf{g}}{\mathbf{H}\mathbf{o}^{k}+\mathbf{b}^{k}}\right]_{j}\frac{\mathbf{o}^{k}_{j}}{1-\lambda^{k}div\left(\frac{\mathbf{\nabla
o}^{k}_{j}}{\lvert\mathbf{\nabla o}^{k}_{j}\rvert}\right)},$ (10)
$\displaystyle\qquad\qquad\qquad\qquad\mathbf{o}^{k+1}_{j}\geq 0,$

where the second line in Eq. (10) is to enforce the non-negativity on $j$-th
reconstructed voxel. Note that for $1\leq j\leq N$, all the voxels in
$\mathbf{o}$ can be updated simultaneously by storing them in a vector.

####  2.2.2 Estimation of the uniform background

When we have $\mathbf{o}^{k+1}$, we substitute
$\mathbf{o}^{k+1},\mathbf{b}^{k}$ and $\lambda^{k}$ in Eq. (7), and take
partial derivatives with respect to $b$

$\frac{\partial\Phi(\mathbf{b}^{k};\mathbf{o}^{k+1},\lambda^{k})}{\partial
b^{k}}=\sum_{i=1}^{M}\left\\{1-\left[\frac{\mathbf{g}}{\mathbf{H}\mathbf{o}^{k+1}+\mathbf{b}^{k}}\right]_{i}\right\\}.$
(11)

Note that each pixel in the captured MFM image is assumed to have the same
background value $b$. Then we update $b$ using the same gradient-based
algorithm as is in 3D image estimation

$b^{k+1}=\left\\{\frac{1}{M}\sum_{i=1}^{M}\left[\frac{\mathbf{g}}{\mathbf{H}\mathbf{o}^{k+1}+\mathbf{b}^{k}}\right]_{i}\right\\}b^{k}.$
(12)

####  2.2.3 Estimation of optimal regularization parameter

Algorithm 1 Joint regularized RL algorithm for MFM

1:inputs: $\mathbf{g}$ and $\mathbf{H}$ $\triangleright$ 2D MFM data and 3D
PSFs

2:initialize: $\mathbf{o}^{0},\mathbf{b}^{0},$ and $\lambda^{0}$

3:set $k=0$

4:while not convergence do

5: update $\mathbf{o}^{k+1}$ using $\mathbf{o}^{k},\mathbf{b}^{k}$ and
$\lambda^{k}$ in Eq. (10)

6: update $\mathbf{b}^{k+1}$ using $\mathbf{o}^{k+1},$ and $\mathbf{b}^{k}$ in
Eq. (12)

7: update $\lambda^{k+1}$ using $\mathbf{o}^{k+1},$ and $\mathbf{b}^{k+1}$ in
Eq. (14)

8: $k=k+1$

9:return: $\mathbf{o}^{*},\mathbf{b}^{*},$ and $\lambda^{*}$ $\triangleright$
joint estimations of 3 unknowns

The approach of estimating optimal $\lambda$ is based on the fact that when
the optimal solution $\mathbf{o}^{*}$ is found, the partial derivatives with
respect to each restored voxel $\mathbf{o}^{k}_{j}$ in Eq. (9) should be equal
or close to zero. So the $\lambda$ is obtained by minimizing the sum of the
square of Eq. (9) over all the reconstructed voxels [26]

$\lambda^{k+1}=\arg\;\min_{\lambda}\sum_{j=1}^{N}\|\frac{\partial\Phi(\mathbf{o}^{k+1};\mathbf{b}^{k+1},\lambda)}{\partial\mathbf{o}^{k+1}_{j}}\|^{2}.$
(13)

Because the objective function in Eq. (13) is a quadratic form of $\lambda$,
the close form solution of $\lambda$ can be found by equating the derivatives
of Eq. (13) with respect to $\lambda$ to zero

$\lambda^{k+1}=\frac{\sum_{j=1}^{N}\left\\{1-\left[\mathbf{H}^{T}\frac{\mathbf{g}}{\mathbf{H}\mathbf{o}^{k+1}+\mathbf{b}^{k+1}}\right]_{j}\right\\}div\left(\frac{\mathbf{\nabla
o}^{k+1}_{j}}{\lvert\mathbf{\nabla
o}^{k+1}_{j}\rvert}\right)}{\sum_{j=1}^{N}{\left[div\left(\frac{\mathbf{\nabla
o}^{k+1}_{j}}{\lvert\mathbf{\nabla o}^{k+1}_{j}\rvert}\right)\right]}^{2}}.$
(14)

The joint regularized RL deconvolution algorithm is summarized in Algorithm 1.
We performed a series of simulations to evaluate the effectiveness and
convergence of the joint regularized RL algorithm for our CMFM system. The
simulation results and convergence plots are shown in Appendix A.

##  3 System analysis

A schematic of our CMFM imaging setup is shown in Fig. 1(a). A regular
microscope is augmented by a $4f$ relay system (L1 = 200mm and L2 = 400mm).
The total magnification of the whole system is $120\times$ with the use of
$60\times$ objective lens. An electron multiplying charge coupled device
(EMCCD) has $1024\times 1024$ pixels with pixel pitch $13\text{um}\times
13\text{um}$. Thus, each pixel corresponds to a size of $108\text{nm}\times
108\text{nm}$ in the object domain, which satisfies the Nyquist sampling
criterion. A multifocal gratings (MFG) is placed at the Fourier plane of the
$4f$ relay system. As an example, schematics of the MFG used in our first
experiment is shown in Fig. 1(b). This MFG is optimized and designed to
produce $3\times 3$ differently focused tile images on a single 2D detector as
shown in Fig. 1(c), with an equal focal step $\Delta z$ of $0.25\text{um}$
between adjacent tiles and almost even light energy distribution of
$[7.56,7.48,7.21,7.47,7.62,7.47,7.21,7.48,7.56]$ percent for each tile. A
total diffraction efficiency of the MFG is therefore $67.06\%$, which is close
to the theoretical maximum efficiency. To avoid the overlapping between tile
images, a field stop slightly smaller than the tile width is placed at the
intermediate image plane to reduce the lateral FOV [see Fig. 1(a)].

Next, we provide analysis on the Point Spread Function (PSF), the chromatic
aberration (CA), and the FOV of our CMFM system.

###  3.1 Point Spread Function (PSF) and resolution

We image a z-stack of a small fluorescent bead (about $170$nm in diameter),
which approximates an ideal light-emitting point, to measure PSF. Figure 2
(first row) shows CMFM lateral PSF imaged under five different axial positions
(columns). Each lateral PSF consists of $3\times 3$ tiles, with one tile in
focus (outlined by a green box; also zoomed in second row) and eight tiles out
of focus. From left to right, as the focus depth increases, we can observe a
focal shift on the axial direction of the x-z plots of PSFs (third row). The
Optical Transfer Functions (OTFs) are plotted as reference (bottom two rows).

We measure the full-width half-magnitude (FWHM) values to quantify the
resolution. Note that the PSFs are affected by chromatic aberration and the
lateral resolution is anisotropic. Therefore, we perform a principle component
analysis (PCA) of xy PSFs. The FWHM values of those five tiles’ PSFs (second
row) are $\\{0.32,0.87,1.17,0.87,1.17\\}\text{um}$ in the CA direction, and
$[0.32,0.37,0.38,0.37,0.36]\text{um}$ in the perpendicular direction. The
axial FWHM values are $\\{1.03,1.39,1.50,1.39,1.55\\}\text{um}$. The variation
of the in-focus PSF sizes results from directional chromatic aberration. Based
on the PCA results, each horizontal and vertical diffraction order tile has a
resolution loss of $0.55\text{um}$ due to CA. Each diagonal diffraction order
tile has a resolution loss of $0.85\text{um}$ due to geometry. We show, in the
next sub-section, that this observation is consistent with the theoretical
calculation.

Figure 2: A z-stack 3D PSFs of CMFM, measured from a $170$nm fluorescent bead.
Top row: CMFM lateral PSFs imaged under five different axial positions
(columns). each PSF consists of one focused image (outlined by a green box)
and eight out-of-focus version images of the bead. xy and xz PSFs (second and
third rows) and corresponding OTFs (bottom two rows) of five differently
focused tiles (columns). The focal shift property of CMFM can be observed from
xz PSFs (third row), verifying that CMFM is capable of capturing a focal stack
instantaneously. Although the central tile’s PSF CA-free (first column), the
off-axis tiles’ PSFs suffer from directional CA (second to last columns) due
to geometry. The lateral spatial frequencies that are lost by CA are shown in
xy OTFs (fourth row).

###  3.2 Chromatic aberration (CA)

Due to the diffractive nature of the MFG, the CMFM system suffers from obvious
CA effect, which causes the loss of the intensity and resolution.
Mathematically, the lateral CA per tile can be expressed as [27]:

$\delta_{x}=\frac{mf\Delta\lambda}{d_{x}};\delta_{y}=\frac{nf\Delta\lambda}{d_{y}},$
(15)

where $m$ and $n$ are horizontal and vertical diffraction orders of the tile,
$f$ is the focal length of the second lens of the $4f$ system, $d_{x}$ and
$d_{y}$ are periods of the MFG in the $x$ and $y$ directions, and
$\Delta\lambda$ is the bandpass spectrum of the emission collected by the
detector. The axial CA for each tile, i.e., different wavelengths are focused
at different distances from the lens, can also be expressed as [27]:

$\delta_{z}=\frac{\Delta\lambda}{\lambda_{c}}f_{m,n},$ (16)

where $\lambda_{c}$ the wavelength used to design the MFG, and
$f_{m,n}=(m+ln)\Delta z$ is the focusing distance of each individual tile.

From Eq. (15) and Eq. (16), we can see that for the central tile where
$m=n=0$, $\delta_{x}=\delta_{y}=\delta_{z}=0$, which means the central tile’s
PSF and image are free of CA. However, for off-axis tile, both lateral and
axial CA exist. Each horizontal and vertical diffraction oder tile has a
dispersion of $mf\Delta\lambda/d_{x}$, which is equal to $0.59\text{um}$ when
$m=1$ for our CMFM system of $f=400\text{mm}$, $\Delta\lambda=10$nm,
$d_{x}=d_{y}=56$um and $\hat{M}=120$. Each diagonal diffraction order tile has
a dispersion of $\sqrt{m^{2}+n^{2}}f\Delta\lambda/d_{x}$. In addition, the
dispersion direction varies for different diffraction order tiles due to
geometry. If left uncorrected, the CA would cause a loss of the the intensity
and resolution in the chromatic dispersion direction. In principle, CA can be
optically corrected by using CA corrective optics. However, the corrective
optics increases system cost and complexity[4]. Here, we show that CA can be
effectively compensated by _computation_. This is possible because the CA is
directional due to geometry of the tile distribution in our CMFM system. For
example, in Fig. 2, the eight non-centric tiles exhibit directional stretching
towards the image center. At depth $z=0$ (first column), the PSF is CA-free.
On the other focal planes, from $\Delta z$ to $4\Delta z$ (second column to
last column), even though the green boxes indicate in-focus tiles, the PSFs
stretch in different directions at different diffraction orders. However, the
stretching does not occur at the perpendicular direction of CA. This implies
that although a certain tile is affected by CA on one direction, the
perpendicular direction is less affected and can be used to compensate another
tile which suffers CA at this direction. Therefore, by jointly utilizing all
the tiles in the model, 3D deconvolution can preserve 2D spatial frequencies
that are lost by CA.

Figure 3: CA blur vs defocus blur. An in-focus tile with CA blur (red) and a
defocus blur (blue) are highlighted in (a), whose PSFs and OTFs are shown in
(b). (c) plots a comparison of linecuts indicated by blue and magenta lines in
(b). For reference, a linecut in CA-free central tile’s OTF (shown in first
column and fourth row of Fig. 2) is also plotted (red). A reconstruction
comparison is shown in the right panel. (d) Object image. (e) Observation
image (for visualization purpose, each tile image is cropped). (f)
Reconstruction using only in-focus PSF. (g) Reconstruction using all the PSFs.

_CA blur vs defocus blur._ Note here that the blur resulted from CA is
different than the defocus blur. Fig. 3 presents a comparison of the two types
of blurs. In Fig. 3(a), the red box highlights the in-focus tile. This in-
focus tile is not located at the image center and has CA blur in horizontal
dimension. Meanwhile, the blue box highlights the tile that suffers both
defocus blur and CA blur but in vertical dimension. However, in horizontal
dimension, this out-of-focus tile only has defocus blur. An Optical Transfer
Function (OTF) comparison is shown in Fig. 3(b). The defocus blur size is
smaller than the CA blur size in horizontal direction, resulting in a wider
OTF lobe, which implies that higher resolution reconstruction can be achieved
by using all the information provided from the PSFs. On the right panel, we
present numerical simulation results. The simulated resolution target [Fig.
3(d)] contains bars pointing four directions. By using only the in-focus PSF,
as shown in Fig. 3(f), only horizontal bars can be resolved, as the high
frequency has only been preserved in vertical direction [Fig. 3(b)]. On the
other hand, by making full use of the PSFs, signals on different directions
can be well-recovered [Fig. 3(g)].

Figure 4: (a) Conventionally, a MFG is designed to project the diffraction
pattern occupying the entire imaging space. (b) We propose to use a smaller
tile spacing, so as to achieve a small lateral FOV that can be tracked over a
large area for MFM tracking applications.

###  3.3 Field-of-View (FOV)

In conventional MFM designs [Fig. 4(a)], the lateral FOV is the same as the
lateral size of each tile image, which can be expressed in the sample domain
as:

$FOV_{x}=\frac{L_{x}}{l\hat{M}};FOV_{y}=\frac{L_{y}}{l\hat{M}},$ (17)

where $L_{x}$ and $L_{y}$ are detector size in the $x$ and $y$ directions, $l$
is the number of tiles in each dimension, and $\hat{M}$ is the magnification
of the MFM. The axial FOV is limited by the focused imaging range of the
designed MFG as

$FOV_{z}=(l^{2}-1)\Delta z.$ (18)

From Eq. (17) and Eq. (18), we can see that by increasing the number of tiles
$l^{2}$, the axial FOV will increase but at the expense of the reduced lateral
FOV for a fixed detector size. Thus, MFM trades its lateral FOV for the depth
resolution. The maximum recovered 3D volume or tracking space in MFM is
$FOV_{x}\times FOV_{y}\times FOV_{z}$.

####  3.3.1 Enlarging lateral tracking space

In conventional MFM designs, large tile spacings are used to ensure that large
objects can be imaged [see Fig. 4(a)]. However, this comes at the cost of a
large reduction in lateral tracking area. This, however, is not optimal for
the MFM 3D tracking applications where the objects of interest (bacteria or
molecule, etc.) are quite small but move freely in a large 3D space. In single
molecule tracking, larger tracking space is more desirable than the lateral
FOV.

By modifying our MFG design to produce a smaller tile spacing, we can achieve
a small lateral FOV that can be tracked over a large area [see Fig. 4(b)]. The
comparison of MFM image diagrams obtained by two different design methods are
shown in Fig. 4. When optimizing the MFG for tracking applications, the
dimensions of the trackable area

$T_{x}=\frac{L_{x}-(l-1)s_{x}}{\hat{M}};T_{y}=\frac{L_{y}-(l-1)s_{y}}{\hat{M}},$
(19)

in $x$ and $y$ dimensions, where $s_{x}$ and $s_{y}$ are tile spacing in x and
y dimensions. Note that Eq. (19) is a general form of Eq. (17). From Eq. (19),
we can clearly see that the lateral tracking area $T_{x}$ and $T_{y}$ can be
enlarged by designing a smaller tile spacing $s_{x}$ and $s_{y}$. Therefore,
our method can achieve a larger lateral tracking area for MFM tracking
applications.

####  3.3.2 Numerical validation

We simulated a 3D tracking space of $1024\times 1024\times 71$ voxels with
each voxel size of $108\text{nm}\times 108\text{nm}\times 200\text{nm}$. The
scene consists of an ellipsoid [Fig. 5(a)] with a diameter of $4.32$um in $x$
and $y$ directions and $8$um in $z$ direction, which is similar to the size of
a bacterium in our experiment [see Fig. 8]. The $xz$ slice and 1D axial
profile of the ellipsoid are shown in Fig. 5(a) (right) and Fig. 5(d) (red
line), respectively. The center of the ellipsoid is placed at $35.8$um
horizontally from the center of the detector.

Figure 5: Simulations showing the capability of the proposed MFM of achieving
larger lateral tracking space than conventional MFM does. (a) The synthetic 3D
ground truth of an ellipsoid (left) and its $xz$ slice (right). The center of
the ellipsoid is $35.8$um away from the center of the detector in $x$
direction. MFM measurements (e-f) and corresponding reconstructions (b-c) by
different design methods. (d) 1D axial profile comparison between ground truth
(red) and reconstructions (black and blue). It is clear that the ellipsoid is
reconstructed poorly from conventional MFM method while our design provides a
good reconstruction. Signal loss as a function of lateral position of the
tracked object for conventional (g) and our designed MFM (h). Similar to
vignetting effect, the signal falls off when approaching the edges. Our
proposed design alleviates peripheral signal loss and achieves an enlarged
lateral tracking area.

For conventional MFM, we used an experimentally captured z-stack PSF for the
simulation. The PSF consists of $5\times 5$ tiles with a tile spacing of $205$
pixels in both $x$ and $y$ dimensions. The 2D MFM measurement [Fig. 5(e)] was
generated based on the forward model of Eq. (3). The maximum signal and
background photon counts are set to be $100$ and $5$, respectively. The
Poisson noise was then added to the measurement by using Matlab’s Poisson
random number generator. For the proposed MFM, the PSF was generated by
cropping the central $101\times 101$ pixels region of each tile of the
experimentally acquired MFM PSF and then recombining the cropped small tiles.
As a result, the tile spacing is $101$ pixels instead of $205$ pixels. The
measurement [Fig. 5(f)] of this new MFM design was also generated based on Eq.
(3) and then corrupted by background and Poisson noise.

In conventional MFM, since the horizontal position of the synthetic ellipsoid
exceeds the tile spacing, the left two-column tiles produced by MFM which
contains axial depth information from $-12\Delta z$ to $-3\Delta z$ are
shifted out of the detector’s FOV and therefore can not recorded by the
detector [Fig. 5(e)]. However, in our new MFM design, the $5\times 5$ tiles
are recorded in their entirety by the detector [Fig. 5(f)]. For
reconstruction, we recovered a 3D space on a grid of $1024\times 1024\times
71$ with each grid size being $108\text{nm}\times 108\text{nm}\times
200\text{nm}$ for both methods. The reconstructed ellipsoid (cropped from the
larger 3D recovery space) and its $xz$ slice by different design methods are
shown in Figs. 5(b) and (c), respectively. The 1D axial profiles of the
reconstructed ellipsoids are also compared in Fig. 5(d) (black and blue
lines). We can clearly see that since the left two-column tiles of the
ellipsoid are missing in the conventional MFM measurement, the ellipsoid is
reconstructed poorly while our design provides a good reconstruction. We also
compare the signal loss as a function of lateral position of the tracked
object in Figs. 5(g) and (h). Similar to vignetting effect, the signal falls
off when approaching the edges. Our proposed design alleviates peripheral
signal loss and preserves full-sized PSF array over an enlarged lateral space.

##  4 Experimental results

We present two experimental results with different biological samples to test
the spatial and temporal performance of our system. As shown in Fig. 6, we
divide the whole image into different sets of tiles by modifying the MFG
pattern. However, the focal step $\Delta z$ is designed to be $0.25$um for
both cases. In Fig. 6(a), the image is divided into $3\times 3$ tiles so that
9 focal planes are in focus, while in Fig. 6(b), $5\times 5$ focal planes are
obtained for the purpose of tracking.

Figure 6: Two experimental snapshot MFM raw images. (a) Experiment 1: snapshot
captured 2D MFM image of multiple static periplasms by using an MFG with $9$
focal planes under exposure time of $0.5$s. (b) Experiment 2: a frame from an
MFM video of a moving bacterium captured at $25$ fps by using an MFG with $25$
focal planes. The raw MFM video is shown in Visualization 1.

###  4.1 Static scene: CMFM _vs_ confocal microscopy

Figure 7: Proposed computational 3D reconstruction of CMFM image in comparison
with confocal deconvolution results. (a) Confocal raw data and (b) its
deconvolution results. (c) CMFM raw data and (d) its computational
reconstruction results. In (a), confocal scan is taken with a dual spinning
disk confocal microscope (Model: CSU-W1) with the total acquisition time of
$20$s, while in (c), the CMFM raw data is captured in a single exposure of
$0.5$s. The lateral resolution of the proposed computational reconstruction is
about $0.35$um (second row of d) and the axial resolution is about $0.5$um
(fourth row of d), which are comparable with those achievable with confocal
deconvolution microscopy (second and fourth rows of b).

We prepare several periplasms (static) as our imaging sample for 3D spatial
resolution characterization. The captured raw image is shown in Fig. 6(a). The
exposure time is $0.5$s. The 3D PSFs are measured by recording a z-stack of a
small fluorescent bead with z focal step of $50$nm. The 3D image are
reconstructed on a $300\times 300\times 51$ grid with each grid size being
$108\text{nm}\times 108\text{nm}\times 50\text{nm}$. Figure 7 shows our
computational 3D reconstruction [Fig. 7(d)] in comparison with the captured
focal stack [Fig. 7(c)] assembled by z-stacking $9$ tiles from the 2D raw MFM
image. For ground truth, we captured the 3D image of periplasms using a
scanning confocal microscope [Fig. 7(a)] and performed deconvolution of the
confocal image [Fig. 7(b)] by using a commercial Huygens software. For each 3D
image, a XY slice at $z=0$ (first row) and a YZ slice at $x=-6.5$um (third
row) are displayed. The comparisons of vertical line profiles indicated by the
red solid line in XY slice, and comparison of axial profiles indicated by the
red dotted line in YZ slice are shown in the second and fourth rows,
respectively. From Fig. 7(c), we can see that without computational
reconstruction, the 3D image assembled by z-stacking nine tile images from 2D
MFM measurement is degraded by noise, out-of-focus blur and low spatial
resolution in all three dimensions. The outer membranes are blurred and the
empty space between them cannot be recognized. In addition, due to the out-of-
focus blur, the FWHM of the axial profile is about $2$um. However, after
deconvolution, the reconstructed 3D image [Fig. 7(d)] is much cleaner, less
out-of-focus blur, and more importantly, has high resolution in three
dimensions. The empty space between outer membranes can be clearly recognized.
The outer membranes which are vertically blurred to a single peak can be well
separated by a dip of $84\%$ between two peaks with the lateral FWHM of
$0.35$um [second row of Fig. 7(d)], which is consistent with FWHM of the
measured PSF. Due to removal of out-of-focus blur, the axial FWHM is about
$0.5$um [fourth row of Fig. 7(d)], with about two times improvement over that
of the raw MFM image. Note that we compare $0.5$s captures with our CMFM
instrument to a $20$s confocal scan taken with a dual spinning disk confocal
microscope (Model: CSU-W1) made by Yokogawa Electric Corporation. Our CMFM
results show similar 3D image quality, but achieve a 40x reduction in
acquisition time. The automatically recovered background noise and optimal
regularizer parameter at each iteration of the joint RL-TV deconvolution
process is shown in Appendix B Fig. 11.

_Cautionary remarks._ Note that although both our CMFM and confocal
deconvolution microscopy achieve $0.5$um axial resolution in the experiment,
the principle of the axial resolution improvement is different for two
techniques. Confocal microscopy increases axial resolution by means of using a
(or multiple) spatial pinhole(s) to block out-of-focus light in image
detection process. Therefore, the axial extent of confocal PSF is narrower
than that in the widefield microscope, and thus a high contrast and resolution
image can be obtained. However, for CMFM, the axial resolution improvement is
based on sparsity-based super-resolution microscopy techniques [28, 29, 30] by
utilizing the signal sparsity in the arbitrary known domain (i.e., gradient
domain in our case). Recently, a new method called SPARCOM [31]: sparsity-
based super-resolution correlation microscopy by utilizing sparsity in the
correlation domain, is also reported to achieve spatial resolution comparable
to PALM and STORM.

###  4.2 3D tracking of moving bacterium

In the second experiment, we demonstrate the 3D video reconstruction of a
moving bacterium (see Visualization 2). The MFG is modified to focus on
$5\times 5$ planes, as shown in Fig. 6(b). A video is captured at $25$ frame
per second (fps) (see Visualization 1). According to our first experimental
reconstruction where the axial FWHM can achieve $0.5$um, we reconstruct the 3D
image of the bacterium on a grid of $150\times 150\times 41$ with each grid
size being $108\text{nm}\times 108\text{nm}\times 200\text{nm}$, which
corresponds to a FOV of $16.2\text{um}\times 16.2\text{um}\times 8\text{um}$
in 3D space for each video frame. Figures 8(a-e) show five out of fifty frames
3D reconstructions. Note that for the visualization purpose, we cropped the
reconstructed 3D volume and just showed $4\text{um}\times 6\text{um}\times
4\text{um}$ region around the bacterium. By computing the center of mass for
each frame reconstruction, a 3D trajectory of the moving bacterium can be
tracked [Fig. 8(f)]. The colorbar in Fig. 8(f) indicates the frame index over
time. The automatically recovered background value and the optimal regularizer
parameter for each video frame is shown in Appendix B Fig. 12.

Figure 8: Experimental 3D reconstructions of a movable bacterium. A raw MFM
video (shown in Visualization 1) was captured at $25$fps as the bacterial
moves in 3D space. The computational 3D reconstruction was performed for each
video frame. Five out of sixty frames reconstruction is shown in (a-e). (f) 3D
trajectory of the bacterium by computing and tracking its center of mass for
each frame reconstruction. The colorbar indicate the frame index over time.
The complete 3D video reconstruction from the first frame to the last frame is
shown in Visualization 2.

##  5 Conclusion

We have presented computational multi-focus microscopy (CMFM), a framework
that balances the imaging system design and computational processing to
achieve single-shot 3D microscopy. Our optical design significantly reduces
the system complexity and experimental alignment by discarding CA correction
optics. We correct for CA computationally rather than optically. This comes
from the fact that CA occurs in different directions at different diffraction
order tiles. Therefore, by jointly utilizing all the tiles in the model, 3D
deconvolution can preserve 2D spatial frequencies that are lost by CA. By
incorporating TV regularization, our algorithm can not only compensate for CA,
but also perform high quality 3D reconstructions from noisy data. We build on
a joint regularized RL deconvolution algorithm and incorporate two types of
noise. Notice that our proposed algorithm is free of any parameter tuning by
automatically estimating the background noise and the optimal reularization
parameter using alternating gradient descent.

We experimentally demonstrate that the out-of-focus blur along $z$ can be
significantly suppressed and the axial resolution as high as $0.5$um is
achievable. The lateral resolution of $0.35$um, which is consistent with the
diffraction-limited resolution, is also experimentally demonstrated. A high
resolution 3D video of a movable bacterium at 25fps is also computationally
reconstructed to verify the proposed CMFM framework.

Finally, we propose a new design method of MFG to enlarge the lateral tracking
area of MFM tracking applications without sacrificing its axial FOV and
single-shot capture speed. The benefits of the simple system design and high
resolution image recovery offered by the proposed CMFM will broaden its
applications in 3D single-shot imaging.

## Appendix A: algorithm evaluation

To verify the effectiveness of the joint regularized RL algorithm for our CMFM
and also quantify the reconstruction quality, we performed a series of
simulations by using an experimentally captured CMFM PSFs. The synthetic 3D
image [Fig. 9(top left)] was obtained from the confocal microscope image of
the 3D bacterial that was imaged under the CMFM. The CMFM measurement image
was generated based on image formation model of Eq. (3) and then degraded with
background and Poisson noise. The maximum signal and background photon counts
are set to be $50$ and $5$, respectively. The Poisson noise was then added to
the measurement by using Matlab’s Poisson random number generator.

Figure 9: Simulations that demonstrate the capability of the joint RL-TV
algorithm to simultaneously recover the 3D image, background noise and the
optimal regularizer parameter for CMFM. Top left: a ground truth image. Top
right: the standard RL deconvolution without TV regularizer ($\lambda=0$).
Bottom left: RL-TV deconvolution by using an incorrect background values
($b=10$). Bottom right: joint RL-TV deconvolution can simultaneously recover a
3D image, background noise and the optimal regularizer parameter. The PSNRs
for three methods are $34.2\text{dB},27.5\text{dB},40.2\text{dB}$, and
I-divergences are $204.3,517$, and $96.7$, respectively.

To quantify the reconstruction quality, we computed peak signal-to-noise ratio
(PSNR) and I-divergence [25] between ground truth image $u$ and the
reconstructed image $v$:

$PSNR=10\log_{10}\left(\frac{\text{MAX}^{2}_{u}}{\text{MSE}}\right),\\\
I_{u,v}=\sum_{i=1}^{N}\left[u_{i}\text{ln}\frac{u_{i}}{v_{i}}-\left(u_{i}-v_{i}\right)\right],$
(20)

where $\text{MAX}_{u}$ is the maximum pixel value of the image $u$, MSE is the
mean square error (MSE) between two images $u$ and $v$, and $i$ stands for the
voxel index.

Figure 9(top right) shows the reconstruction of the standard RL deconvolution
without TV constraint by setting $\lambda=0$. Because TV prior is not used,
the reconstruction can not preserve the smooth edge of the original image and
contains the artifacts. To investigate the effect of the background noise on
the CMFM reconstruction quality, we performed RL-TV deconvolution but with an
incorrect background value of $10$ (the truth value is 5). The reconstructed
image is shown in Fig. 9(bottom left). Since the used background value is two
times bigger than the real one, the recovered 3D image suffers from
substantial signal loss. The reconstructed 3D image by joint RL-TV algorithm
is shown in Fig. 9(bottom right). The PSNRs between the original image and
reconstructed images by three methods are
$34.2\text{dB},27.5\text{dB},40.2\text{dB}$, and I-divergences are
$204.3,517$, and $96.7$, respectively. Clearly, the joint RL-TV algorithm,
which simultaneously recovers 3D image, the background noise value and optimal
regularizer parameter, gave the best reconstruction quality of CMFM among
three methods. Figures 10(a-d) show the reconstructed background value $b$,
the optimal regularizer parameter $\lambda$, PSNR and I-divergence at each
iteration of the deconvolution process. Note that both the background value
and regularizer parameter were initialized to be $100$, but they gradually
converged to $4.99$ and $2.2\times 10^{-3}$, respectively. From Fig. 10(c-d),
we can also see that the image quality is dramatically improved during the
iterative reconstruction process, proving the effectiveness of joint RL-TV
algorithm for the 3D reconstruction of our snapshot CMFM system.

Figure 10: The convergence analysis of the joint RL-TV algorithm for the
simulated CMFM. (a) The reconstructed background noise value and (b) optimal
regularizer parameter during the iteration of joint RL-TV deconvolution
process. (c) PSNR and (d) I-divergence between the ground truth image and
reconstructed image at each iteration of the deconvolution process.

## Appendix B: experimental results

We also tested the algorithm on two types of real data. The 3D reconstruction
of static periplasms is shown in Fig. 7. The automatically recovered
background noise and optimal regularizer parameter at each iteration of the
joint RL-TV deconvolution process is shown in Fig. 11.

Figure 11: The recovered background values (left) and optimal regularizer
parameter $\lambda$ (right) at each iteration of the joint RL-TV deconvolution
process for the first experiment.

For the second experiment, we captured an MFM video of a moving bacterium at
25 frames per second (fps) using our CMFM techniques. The complete 3D video
reconstruction from the first frame to the last frame is shown in
Visualization 2. The automatically recovered background value and the optimal
regularizer parameter for each video frame is shown in Fig. 12.

Figure 12:  The recovered background values (left) and optimal regularizer
parameter $\lambda$ (right) for each video frame of a movable bacterium.

## Funding

Biological Systems Science Division, Office of Biological and Environmental
Research, Office of Science, U.S. Dept. of Energy, under Contract DE-
AC02-06CH11357; NSF CAREER grant IIS-1453192; ONR award N00014-15-1-2735.

## Disclosures

The authors declare that there are no conflicts of interest related to this
article.

## References

  * [1] P. M. Blanchard and A. H. Greenaway, “Simultaneous multiplane imaging with a distorted diffraction grating,” Appl. Opt. 38, 6692–6699 (1999).
  * [2] S. Abrahamsson, R. Ilic, J. Wisniewski, B. Mehl, L. Yu, L. Chen, M. Davanco, L. Oudjedi, J.-B. Fiche, B. Hajj, X. Jin, J. Pulupa, C. Cho, M. Mir, M. E. Beheiry, X. Darzacq, M. Nollmann, M. Dahan, C. Wu, T. Lionnet, J. A. Liddle, and C. I. Bargmann, “Multifocus microscopy with precise color multi-phase diffractive optics applied in functional neuronal imaging,” Biomed. Opt. Express 7, 855–869 (2016).
  * [3] L. Oudjedi, J.-B. Fiche, S. Abrahamsson, L. Mazenq, A. Lecestre, P.-F. Calmon, A. Cerf, and M. Nöllmann, “Astigmatic multifocus microscopy enables deep 3d super-resolved imaging,” Biomed. Opt. Express 7, 2163–2173 (2016).
  * [4] S. Abrahamsson, B. H. J. Chen, S. Stallinga, A. Y. Katsov, J. Wisniewski, G. Mizuguchi, P. Soule, F. Mueller, C. D. Darzacq, X. Darzacq, C. Wu, C. I. Bargmann, D. A. Agard, M. Dahan, , and M. G. L. Gustafsson, “Fast multicolor 3d imaging using aberration-corrected multifocus microscopy,” Nature methods 10, 60–63 (2013).
  * [5] R. Prevedel, Y.-G. Yoon, M. Hoffmann, N. Pak, G. Wetzstein, S. Kato, T. Schrödel, R. Raskar, M. Zimmer, E. S. Boyden, , and A. Vaziri, “Simultaneous whole-animal 3d imaging of neuronal activity using light-field microscopy,” Nature methods 11, 727–730 (2014).
  * [6] N. C. Pégard, H.-Y. Liu, N. Antipa, M. Gerlock, H. Adesnik, and L. Waller, “Compressive light-field microscopy for 3d neural activity recording,” Optica 3, 517–524 (2016).
  * [7] J. K. Adams, V. Boominathan, B. W. Avants, D. G. Vercosa, F. Ye, R. G. Baraniuk, J. T. Robinson, and A. Veeraraghavan, “Single-frame 3d fluorescence microscopy with ultraminiature lensless flatscope,” Science Advances 3 (2017).
  * [8] A. F. Coskun, I. Sencan, T.-W. Su, and A. Ozcan, “Lensless wide-field fluorescent imaging on a chip using compressive decoding of sparse objects,” Opt. Express 18, 10510–10523 (2010).
  * [9] N. Antipa, G. Kuo, R. Heckel, B. Mildenhall, E. Bostan, R. Ng, and L. Waller, “Diffusercam: lensless single-exposure 3d imaging,” Optica 5, 1–9 (2018).
  * [10] J. N. Mait, “Understanding diffractive optic design in the scalar domain,” J. Opt. Soc. Am. A 12, 2145–2158 (1995).
  * [11] G. Pedrini and S. Schedin, “Short coherence digital holography for 3d microscopy,” Optik-International Journal for Light and Electron Optics 112, 427–432 (2001).
  * [12] B. W. Schilling, T.-C. Poon, G. Indebetouw, B. Storrie, K. Shinoda, Y. Suzuki, and M. H. Wu, “Three-dimensional holographic fluorescence microscopy,” Optics Letters 22, 1506–1508 (1997).
  * [13] Z. Wang, D. Ryu, K. He, O. Cossairt, and A. K. Katsaggelos, “4d tracking of biological samples using lens-free on-chip in-line holography,” in _Digital Holography and Three-Dimensional Imaging,_ (Optical Society of America, 2017), pp. Tu2A–4.
  * [14] L. Tian, J. Wang, and L. Waller, “3d differential phase-contrast microscopy with computational illumination using an led array,” Optics letters 39, 1326–1329 (2014).
  * [15] L. Tian and L. Waller, “3d intensity and phase imaging from light field measurements in an led array microscope,” optica 2, 104–111 (2015).
  * [16] J. Rosen and G. Brooker, “Digital spatially incoherent fresnel holography,” Optics letters 32, 912–914 (2007).
  * [17] J. Rosen and G. Brooker, “Non-scanning motionless fluorescence three-dimensional holographic microscopy,” Nature Photonics 2, 190 (2008).
  * [18] O. Cossairt, K. He, R. Shang, N. Matsuda, M. Sharma, X. Huang, A. Katsaggelos, L. Spinoulas, and S. Yoo, “Compressive reconstruction for 3d incoherent holographic microscopy,” in _Image Processing (ICIP), 2016 IEEE International Conference on,_ (IEEE, 2016), pp. 958–962.
  * [19] D. Ryu, Z. Wang, K. He, G. Zheng, R. Horstmeyer, and O. Cossairt, “Subsampled phase retrieval for temporal resolution enhancement in lensless on-chip holographic video,” Biomedical optics express 8, 1981–1995 (2017).
  * [20] Z. Wang, Q. Dai, D. Ryu, K. He, R. Horstmeyer, A. Katsaggelos, and O. S. Cossairt, “Dictionary-based phase retrieval for space-time super resolution using lens-free on-chip holographic video,” in _Computational Optical Sensing and Imaging,_ (Optical Society of America, 2017), pp. CTu2B–3.
  * [21] H. Y. Liu, J. Zhong, and L. Waller, “Multiplexed phase-space imaging for 3d fluorescence microscopy,” Optics Express 25, 14986–14995 (2017).
  * [22] Z. Wang, L. Spinoulas, K. He, L. Tian, O. Cossairt, A. K. Katsaggelos, and H. Chen, “Compressive holographic video,” Optics express 25, 250–262 (2017).
  * [23] M. Broxton, L. Grosenick, S. Yang, N. Cohen, A. Andalman, K. Deisseroth, and M. Levoy, “Wave optics theory and 3-d deconvolution for the light field microscope,” Opt. Express 21, 25418–25439 (2013).
  * [24] P. J. GREEN, “On the use of em algorithm for penalized likelihood estimation,” J. Royal Statist. Soc. B 52, 443–452 (1990).
  * [25] N. Dey, L. Blanc-Feraud, C. Zimmer, P. Roux, Z. Kam, J. C. Olivo-Marin, and J. Zerubia, “Richardson-lucy algorithm with total variation regularization for 3d confocal microscope deconvolution,” Microsc. Res. Tech. 4, 260–266 (2006).
  * [26] M. Laasmaa, M. Vendelin, and P. Peterson, “Application of regularized richardson-lucy algorithm for deconvolution of confocal microscopy images,” J. Microsc 2, 124–140 (2011).
  * [27] K. He, X. Huang, X. Wang, S. Yoo, P. Ruiz, I. Gdor, N. J. Ferrier, N. Scherer, M. Hereld, A. K. Katsaggelos, and O. Cossairt, “Design and simulation of a snapshot multi-focal interferometric microscope,” Opt. Express 26 (2018).
  * [28] S. Hugelier, J. J. De Rooi, R. Bernex, S. Duwé, O. Devos, M. Sliwa, P. Dedecker, P. H. Eilers, and C. Ruckebusch, “Sparse deconvolution of high-density super-resolution images,” Scientific reports 6, 21413 (2016).
  * [29] S. Hugelier, P. Eilers, O. Devos, and C. Ruckebusch, “Improved superresolution microscopy imaging by sparse deconvolution with an interframe penalty,” Journal of Chemometrics 31, e2847 (2017).
  * [30] A. Szameit, Y. Shechtman, E. Osherovich, E. Bullkich, P. Sidorenko, H. Dana, S. Steiner, E. B. Kley, S. Gazit, T. Cohen-Hyams _et al._ , “Sparsity-based single-shot subwavelength coherent diffractive imaging,” Nature materials 11, 455 (2012).
  * [31] O. Solomon, M. Mutzafi, M. Segev, and Y. C. Eldar, “Sparsity-based super-resolution microscopy from correlation information,” Opt. Express 26, 18238–18269 (2018).


# Analyzing the contribution of different passively collected data to predict
Stress and Depression

Irene Bonafonte1,7, Cristina Bustos1,4, Abraham Larrazolo2,4, Gilberto Lorenzo
Martínez Luna2,
Adolfo Guzmán Arenas2, Xavier Baró1, Isaac Tourgeman3, Mercedes Balcells4,6
and Agata Lapedriza1,5 Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>1 Universitat Oberta de Catalunya, Barcelona, Spain  2
Instituto Politecnico Nacional, CDMX, Mexico  3 Albizu University, Miami, FL,
USA  4 Massachusetts Institute of Technology, Cambridge, MA, USA
5 Northeastern University, Boston, MA, USA  6 IQS School of Engineering,
Universitat Ramon Llull, Barcelona, Spain
7 Helmholtz Munich, Neuherberg, Germany

###### Abstract

The possibility of recognizing diverse aspects of human behavior and
environmental context from passively captured data motivates its use for
mental health assessment. In this paper, we analyze the contribution of
different passively collected sensor data types (WiFi, GPS, Social
interaction, Phone Log, Physical Activity, Audio, and Academic features) to
predict daily self-report stress and PHQ-9 depression score. First, we compute
$125$ mid-level features from the original raw data. These $125$ features
include groups of features from the different sensor data types. Then, we
evaluate the contribution of each feature type by comparing the performance of
Neural Network models trained with all features against Neural Network models
trained with specific feature groups. Our results show that WiFi features
(which encode mobility patterns) and Phone Log features (which encode
information correlated with sleep patterns), provide significative information
for stress and depression prediction.

###### Index Terms:

depression prediction, stress prediction, Digital Phenotyping, feature
importance

##  I Introduction

Wearable devices and smartphones have enabled the collection of large amounts
of data that reflect human behavior patterns. This has given rise to a field
called Digital Phenotyping, which attempts to quantify behavior using
passively captured mobile data [1]. This field has gained interest in mental
health diagnosis due to its ability to capture behaviors and habits like
sociability, physical activity, mobility, sleep, among others [2]. Mental
health conditions are highly prevalent. However, the diagnosis and management
of these conditions face limitations like barriers in seeking help or reliance
on self-reports during short clinical visits [3]. Hence, Digital Phenotyping
can be a potential tool to support decisions in a clinical setting.

Extensive research in Digital Phenotyping aims to predict mental states or
mental disorders, but there are still a lot of open questions. One of them is
understanding well what type of data is actually relevant for mental health
assessment. From the Machine Learning perspective, most of the solutions do
not provide sufficient information about which factors or what types of
features contribute to the model’s inference of certain mental state or
disorder. As a consequence, it is difficult for mental health specialists to
rely on model’s predictions to evaluate or intervene in patient diagnosis [2].

In this work, we study the contribution of different data types, coming from
different sensors and sources of information, like GPS, WiFi, or phone logs
for predicting mental health states using deep learning models. Concretely, we
focus on two particular tasks: prediction of the depression score and
prediction of the self-report stress level. We use the StudentLife dataset [4]
for this study. StudentLife is one of the few public and available datasets
for Digital Phenotyping.

Inspired by the feature analysis done in [5], we study the importance of
features by analyzing the performance of a neural network trained with all the
features extracted from passive sensing data or trained just using a group of
features. We can gain insight into which features are more relevant for
predicting depression or stress. Until now, we have not found a previous study
that evaluates the importance of features for prediction of the depression
score and prediction of the stress level using the StudentLife dataset. In our
analysis, we find the WiFi features, which represent the mobility patterns of
a student, are the most discriminant for the neural models for both prediction
tasks. Additionally, phone log features show higher discrimination information
for stress level classification. For reproducibility purposes, our feature
extraction and models can be found in our public code
repository111https://github.com/cristinabustos16/digitalPhenotyping_featuresAnalysis.

###  I-A Related Work

Our study is conducted on the StudentLife dataset [4], which includes data
from 48 students from Dartmouth College passively monitored during a 10-week
academic term. The dataset contains passive sensing data from different
sensors and an Ecological Momentary Assessment (EMA) component that provides
self-reported data such as clinical survey responses, daily stress, sleep
hours, mood, and exercise, among others. The StudentLife dataset has been
widely used to study depression diagnosis [6, 7] sedentary behavior prediction
and activity recommendation [8], food purchase prediction [9], recommendations
for improving academic performance [10], mood monitoring [11] and stress level
prediction [12, 13, 14]

Stress level prediction on the StudentLife dataset has been approached in
different studies. It was first explored in [15], authors categorized the
stress level into below-median, median or above-median based on the user’s
median, and use a Fully Connected Network with mobility and temporal features.
Later, Shaw et al. [12] utilized 24-hour histograms of passive data and Long-
Short-Term-Memory (LSTM) autoencoder to transform temporal data into a vector
representation, then trained a personalized multi task neural network on a
reduced sample of 23 students. Furthermore, stress prediction has been
addressed as a binary classification problem, as seen in [13], where authors
used LSTMs on the last 2-12h of sensor data to predict between stress / not
stress.

Predicting depression using the StudentLife dataset has been addressed by
dividing the depression score into two classes (depressed/not depressed). For
example, in [16, 17] authors computed features from the GPS sensor data (e.g.
location variance, speed mean, total distance, cluster of locations, circadian
movement, entropy and heuristic that indicates the time spent at home). Saeb
et al. [17] found the mobility features grouped in blocks of 2-weeks are
significantly correlated with the depression score. In [16], authors extract
mobility features with an autoencoder and then use an SVM classifier on top of
these features. Colbaugh et al. [7] used GPS, WiFi, screen lock, light sensor,
and microphone data, while Kim et al. [6] used the screen lock and unlock data
for the last 2-weeks.

##  II From Passive Data to Features

Inspired by the hierarchical framework introduced by Mohr et al. [18], which
transforms raw data into informative mid-level features (e.g. location,
exercise, bedtime, phone usage patterns, and social activity), we compute a
collection of $125$ informative features, representing the student’s daily
information divided into three periods as defined by Wang et al. [4]: day (9
am - 6 pm), night (12 am - 9 am), and evening (6 pm - 12 am). We divided the
features into groups according to the sensor or/and behavior patterns
representation. We excluded data contained in the EMA like self-reported sleep
hours, exercise patterns, because we analyze the model performance mostly
using passive data. The rest of this section describes the features we
computed:

WiFi (36 features): WiFi data tracks the student’s connection time to
different locations on the campus (over 60 different locations), reflecting
their mobility patterns. We compute two summary variables: the number of
different locations visited and the variance in time spent at each location.
Additionally, we include 3 variables representing the time spent in the
student’s three most visited locations, typically reflecting their homestay
and classroom building. Furthermore, we include 7 variables that indicate the
time spent in the 7 most frequently visited locations for all students.

GPS (30 features): From the GPS data, we compute and aggregate the maximum
distance, the total distance, distance variance, mean speed, speed variance,
area of the convex hull that contains all the visited coordinates, and the
total time spent indoors and outdoors. Notice that this information is
complementary to just identifying frequented locations, which are already
encoded with the WiFi features.

Social (9 features): We computed number of SMS, calls, app usages, Bluetooth
contacts, and the total time in seconds the student spend in a phone
conversation.

Phone Log (14 features): These features provide information about the duration
in seconds of the smartphone in certain state: phone charging time, phone lock
time, and light sensor (time in a dark environment).

Activity (12 features): These features encode the time spent doing physical
activity: running, walking, and stationary.

Audio (9 features): These features represent the time the student spent in
silent, conversational or noisy environment.

Academic (13 features): We include the academic performance information (GPA),
the interaction with the campus website like the number of views,
contributions, questions, notes and answers. We also add variables such as
days to the nearest deadline, daily class hours, and the day of the week.

##  III Neural Network Architectures

TABLE I: Stress Level Classification. Results (F1 and accuracy) for Binary L-H
, Binary LM-H, and Multiclass with FCN -top- and LSTM -bottom- models. Best
and second best results per each model (row) and marked in bold face. FCN
---
| metric | most freq. | all features | WiFi | GPS | Social | phone log | Activity | Audio | Academic
L-H | F1 | 48,3 |  61,1$\pm$ 3,5 | 53,7 $\pm$ 3,9 | 46,8 $\pm$ 3,6 | 45,6 $\pm$ 3,2 |  55,6 $\pm$ 3,9 | 50,6 $\pm$ 2,9 | 50,4 $\pm$ 3,4 | 49,0 $\pm$ 2,7
Acc. | 52,6 |  60,8 $\pm$ 3,2 | 53,4 $\pm$ 3,9 | 46,8 $\pm$ 3,4 | 45,5 $\pm$ 3,5 |  55,4 $\pm$ 3,8 | 50,4 $\pm$ 3,2 | 50,1 $\pm$ 3,6 | 48,8 $\pm$ 2,6
LM-H | F1 | 38,8 |  61,4 $\pm$ 3,8 | 55,8 $\pm$ 2,9 | 52,5 $\pm$ 3,2 | 50,6 $\pm$ 3,7 |  56,8 $\pm$ 3,0 | 53,7 $\pm$ 3,6 | 55,8 $\pm$ 3,6 | 52,4 $\pm$ 2,9
Acc. | 40,1 |  64,7 $\pm$ 3,0 | 59,8 $\pm$ 2,4 | 58,6 $\pm$ 1,8 | 56,1 $\pm$ 2,6 |  60,5 $\pm$ 2,4 | 58,8 $\pm$ 2,7 | 59,9 $\pm$ 2,4 | 58,6 $\pm$ 2,1
Multiclass | F1 | 24,2 |  50,9 $\pm$ 3,2 |  43,8 $\pm$ 3,9 | 43,7 $\pm$ 2,4 | 40,0 $\pm$ 2,8 | 43,5 $\pm$ 2,5 | 43,5 $\pm$ 2,6 | 43,6 $\pm$ 2,3 | 42,9 $\pm$ 3,5
Acc. | 32,3 |  51,9 $\pm$ 3,2 |  44,9 $\pm$ 2,4 | 44,8 $\pm$ 2,4 | 41,3 $\pm$ 2,9 | 44,8 $\pm$ 2,5 | 44,5 $\pm$ 3,0 | 44,5 $\pm$ 2,3 | 44,1 $\pm$ 3,4
LSTM
| metric | most freq. | all features | WiFi | GPS | Social | Phone log | Activity | Audio | Academic
L-H | F1 | 48,3 |  61,4 $\pm$ 5,0 |  57,5 $\pm$ 2,9 | 49,7 $\pm$ 2,9 | 48,8 $\pm$ 2,9 | 54,5 $\pm$ 4,9 | 51,8 $\pm$ 2,7 | 49,0 $\pm$ 3,1 | 48,8 $\pm$ 2,8
Acc. | 52,6 |  61,5 $\pm$ 2,8 |  57,8 $\pm$ 2,7 | 49,7 $\pm$ 3,4 | 49,1 $\pm$ 3,0 | 55,9 $\pm$ 4,3 | 52,7 $\pm$ 3,0 | 48,7 $\pm$ 3,2 | 48,8 $\pm$ 2,6
LM-H | F1 | 38,8 |  58,5 $\pm$ 4,7 | 58,4 $\pm$ 4,2 | 53,7 $\pm$ 3,0 | 52,5 $\pm$ 2,6 |  59,0 $\pm$ 3,1 | 56,6 $\pm$ 4,2 | 55,3 $\pm$ 3,2 | 52,6 $\pm$ 3,5
Acc. | 40,1 |  62,1 $\pm$ 2,1 |  62,3 $\pm$ 2,2 | 58,9 $\pm$ 2,2 | 57,4 $\pm$ 3,0 | 61,6 $\pm$ 1,6 | 61,1 $\pm$ 2,1 | 60,9 $\pm$ 1,7 | 58,6 $\pm$ 2,1
Multiclass | F1 | 24,2 |  52,4 $\pm$ 5,0 | 48,0 $\pm$ 2,8 | 43,6 $\pm$ 2,7 | 42,6 $\pm$ 2,4 |  48,5 $\pm$ 3,1 | 44,7 $\pm$ 3,0 | 44,1 $\pm$ 3,2 | 43,6 $\pm$ 3,9
Acc. | 32,3 |  53,9 $\pm$ 2,8 | 49,9 $\pm$ 2,4 | 46,5 $\pm$ 2,7 | 46,4 $\pm$ 3,0 |  50,4 $\pm$ 1,6 | 47,5 $\pm$ 2,6 | 46,8 $\pm$ 3,4 | 44,1 $\pm$ 3,4
TABLE II: Root Mean Square Error (RMSE) for the regression of PHQ-9 depression score using the FCN model. Best and second best results and marked in bold face. | mean | all features | WiFi | GPS | Social | Phone log | Activity | Audio | Academic
---|---|---|---|---|---|---|---|---|---
PHQ-9 regression | 4.4 |  4,1 $\pm$ 4,4 |  3,5 $\pm$ 3,4 | 4,9 $\pm$ 3,9 | 4,7 $\pm$ 4,3 | 5 $\pm$ 4,1 | 4,9 $\pm$ 4,7 | 5,1 $\pm$ 4,2 | 4,9 $\pm$ 4,1

###  III-A Modeling daily self-reported stress level

We predict daily stress levels based on current (and past) passive data. We
rescale the original stress data that comes in a 1-5 scale into 3 stress
levels: below-median (low stress), median (medium stress) and above-median
(high stress). We compute the median separately for each individual, as done
in [15]. We test 3 different problem configurations: binary low versus high
classification (L-H); binary low and medium versus high classification (LM-H);
and multiclass classification (low/medium/high). We train two types of neural
networks: a Fully Connected Network (FCN) that takes passive data from a
single day as input; and a Recurrent Neural Network (LSTM) that takes data
from 5 past days and predicts the stress of the last day. The FCN is composed
by 3 layers (57, 35, 3), each one followed by ReLU activation and a dropout
layer (0.35, 0.15, 0.15). The recurrent network is composed by a LSTM layer
with 50 nodes and a dropout of 0.2, followed by a fully connected layer of 15
neurons with ReLU activation. We order chronologically the data, used the 80%
first days of the study for training and the last 20% for testing. For this
task, we also add an extra one-hot encoded variable indicating the user.

###  III-B Modeling depression score (PHQ-9)

For predicting depression, we used the PHQ-9 (Patient Health Questionnaire), a
self-report measure of depression symptom severity, with scores ranging from 0
to 27 (each item scored from 0 to 3). In total, 38 students applied for the
PHQ-9 at the end of the study. For this purpose, we created a FCN with Root
Mean Square Error (RMSE) loss to predict the depression score. The FCN is
composed by 3 layers (128 neurons) each one followed by ReLU activation
function and a dropout layer of 0.3. We perform leave-one-out evaluation,
training using 37 subjects and testing with one subject. The input data
consists of aggregated features representing the mean and standard deviation
of each feature over a period of two weeks leading up to the end of the study.

##  IV Results and Discussion

The accuracy and F1 score for stress level classification are reported in
Table I. In turn, Table II reports the Root Mean Square Error (RMSE) for the
depression score regression. We train each model with all $125$ features, as
well as with each group of features separately. In all the cases we perform
$50$ rounds of training and testing with different random initializations of
the model. The tables report the mean and standard deviation results obtained
per each performance measure. As a baseline, we include the accuracy and F1
score for a classifier that always predicts the most frequent stress level
value for each student and RMSE for a regressor that always predicts the mean
of all the PHQ-9 scores.

Performance results for stress level classification. The results in Table I
show that, in general, the best performance is obtained when using all the
features. Surprisingly, for the binary predictions the FCN model outperforms
the LSTM, suggesting that using passive sensor data from the same day is more
informative than incorporating data from the past 5 days. For the multi-class
classification, the LSTM performs better than FCN, indicating that more
context from past days is needed for a finer-grained classification problem.

Performance results for PHQ-9 regression. Table II shows that the model
trained with mobility features embedded in WiFi data outperforms the model
using all features. WiFi features, in particular, also exhibit a lower
standard deviation, which supports the finding that these features have
significant relevance for predicting depression. This can be attributed to the
correlation between the time a person spends at home.

Feature importance. We observe that WiFi and Phone log features exhibit higher
performance when used as individual group features, when compared to the other
feature groups. WiFi data captures the mobility patterns of students within
the experiment, while phone log features indirectly represent sleep patterns
through indicators such as phone charging and time spent in the dark
environment. These features also have lower standard deviations, indicating
their consistency and relevance. On the contrary, social features show the
worst performance overall, suggesting that metrics like the number of calls or
SMS may not be highly influential for stress level prediction in this dataset.
However, we notice that social features represent the smallest feature group,
and their limited impact could be due to their lack of information for stress
prediction. The results suggest that mobility patterns and sleep habits have
relevant information for stress prediction.

Comparison of stress level classification with previous works. Previous
studies on stress level classification with the StudentLife dataset report
results that are comparable with the results obtained by our approach. For
example, for the binary stress level classification, In [13], an accuracy of
62.8% is reported, while our model reaches from 60% to almost 65%. Kadri et al
[14] reports 93% but they are using active data that can have a direct
correlation with stress, so it is not comparable with our approach. In
particular, the state-of-the-art results in multi class stress prediction with
the StudentLife dataset is obtained with a complex multitask architecture
[12], difficult to generalize to new unseen participants and using a reduced
subset of students. Furthermore, we have observed a strong sensitivity to
model initialization. For example, notice that the standard deviations range
from $1,6$ to $3,9$ in the case of accuracy for the binary predictions.
Unfortunately, none of these previous works has published their code, trained
models, training-test data partitions or details on hyperparameters, which
does not allow us to reproduce their results and perform a fair comparison. We
encourage the practice of releasing code, so the community can fairly
reproduce the reported results.

Comparison of depression prediction with previous works. Previous studies
primarily focused on detecting depression using the Studentlife dataset, and
they mainly relied on classification techniques [7]. However, our research
takes a complementary approach by focusing on regression analysis. Our
findings align with the results reported in [17], which emphasize the
significance of mobility features for detecting depression. Again, we also
notice that none of these previous works released their code, trained models,
or data partition details, making reproducibility and comparison not possible.

Weaknesses. In our study, the number of features differs depending on the
group, and this might have an impact on the amount of information each group
of features might encode. The encoded features are inspired by previous works,
but other additional features could be computed from the raw data, which might
also have an impact on the results.

##  V Conclusions

In this study we use the StudentLife dataset to compare the contribution of
various passive data features for stress level prediction and depression score
regression. We found that different feature sets exhibited varying performance
levels when tested with Neural Network models. WiFi features were particularly
informative for stress and depression prediction, highlighting the importance
of mobility patterns. Phone log features has relevance in stress level
prediction, indicating the relation of sleep patterns and mental health
states.

## Ethical Impact Statement

Our work focuses on evaluating the importance of different sensor data
features in predicting depression and stress in students. Our analysis aims to
contribute the understanding of mental health states and enhancing
interventions. However, we acknowledge the limitations of our study and
emphasize the need for further research and validation before implementing it
in real-world diagnosis. It is crucial to recognize the complexity of mental
health disorders, which may not be fully reflected in the data collected from
wearable devices or smartphones. In our experiments, we use a publicly
available and anonymized dataset where participants gave their consent for the
collection of personal data. However, this dataset is very small and does not
represent global student’s conditions and general human behavior patterns.
Therefore, larger, more balanced, and diverse datasets are necessary for
future research. Furthermore, considerations such as privacy, data protection,
and participant well-being should be carefully addressed when creating these
datasets and training machine learning models. Lastly, machine learning models
should not be seen as a replacement of the diagnosis provided by clinicians,
but rather serve as a supportive tool.

## Acknowledgments

This work is partially supported by the Spanish Ministry of Science,
Innovation and Universities RTI2018-095232-B-C22 to A.L and by Fundacion
Bancaria La Caixa (ref. 2019-198898) to M.B. We also thank NVIDIA for their
generous hardware donations.

## References

  * [1] I. Perez-Pozuelo, D. Spathis, E. A. Clifton, and C. Mascolo, “Wearables, smartphones, and artificial intelligence for digital phenotyping and health,” in _Digital Health_. Elsevier, 2021.
  * [2] I. Moura, A. Teles, D. Viana, J. Marques, L. Coutinho, and F. Silva, “Digital phenotyping of mental health using multimodal sensing of multiple situations of interest: A systematic literature review,” _Journal of Biomedical Informatics_ , 2022.
  * [3] S. Z. Williams, G. S. Chung, and P. A. Muennig, “Undiagnosed depression: A community diagnosis,” _SSM-population health_ , 2017.
  * [4] R. Wang, F. Chen, Z. Chen, T. Li, G. Harari, S. Tignor, X. Zhou, D. Ben-Zeev, and A. T. Campbell, “Studentlife: assessing mental health, academic performance and behavioral trends of college students using smartphones,” in _Proceedings of the 2014 ACM international joint conference on pervasive and ubiquitous computing_ , 2014.
  * [5] T. Umematsu, A. Sano, S. Taylor, and R. W. Picard, “Improving students’ daily life stress forecasting using lstm neural networks,” in _2019 IEEE EMBS international conference on biomedical & health informatics (BHI)_. IEEE, 2019.
  * [6] J. Kim, J. Hong, and Y. Choi, “Automatic depression prediction using screen lock/unlock data on the smartphone,” in _2021 18th International Conference on Ubiquitous Robots (UR)_. IEEE, 2021.
  * [7] R. Colbaugh, K. Glass, and V. Global, “Detecting and monitoring brain disorders using smartphones and machine learning,” _medRxiv_ , 2020.
  * [8] J. A. R. Rojas, J. Rosas, Y. Shen, H. Jin, and A. K. Dey, “Activity recommendation: Optimizing life in the long term,” in _IEEE International Conference on Pervasive Computing and Communications_ , 2020.
  * [9] F. Chen, R. Wang, X. Zhou, and A. T. Campbell, “My smartphone knows i am hungry,” in _Proceedings of the 2014 workshop on physical analytics_ , 2014\.
  * [10] H. Djeghri, H. Seridi, M. Messaadia, and M. T. Khadir, “A recommendation system for improving student life based on neural network and matrix vectorization,” in _2021 1st International Conference On Cyber Management And Engineering (CyMaEn)_. IEEE, 2021.
  * [11] X. Ma, X. Yang, J. Gao, and C. Xu, “Health status prediction with local-global heterogeneous behavior graph,” _ACM Transactions on Multimedia Computing, Communications, and Applications (TOMM)_ , 2021.
  * [12] A. Shaw, N. Simsiri, I. Deznaby, M. Fiterau, and T. Rahaman, “Personalized student stress prediction with deep multitask network,” _arXiv preprint arXiv:1906.11356_ , 2019.
  * [13] Y. Acikmese and S. E. Alptekin, “Prediction of stress levels with lstm and passive mobile sensors,” _Procedia Computer Science_ , 2019.
  * [14] N. Kadri, S. H. Turki, A. Ellouze, and M. Ksantini, “An hybrid deep learning approach for prediction and binary classification of student’s stress,” in _Intelligent Systems and Pattern Recognition: Second International Conference_. Springer, 2022.
  * [15] G. Mikelsons, M. Smith, A. Mehrotra, and M. Musolesi, “Towards deep learning models for psychological state prediction using smartphone data: Challenges and opportunities,” _arXiv:1711.06350_ , 2017.
  * [16] W. Gerych, E. Agu, and E. Rundensteiner, “Classifying depression in imbalanced datasets using an autoencoder-based anomaly detection approach,” in _2019 IEEE 13th International Conference on Semantic Computing (ICSC)_. IEEE, 2019.
  * [17] S. Saeb, E. G. Lattie, S. M. Schueller, K. P. Kording, and D. C. Mohr, “The relationship between mobile phone location sensor data and depressive symptom severity,” _PeerJ_ , 2016.
  * [18] D. C. Mohr, M. Zhang, and S. M. Schueller, “Personal sensing: understanding mental health using ubiquitous sensors and machine learning,” _Annual review of clinical psychology_ , 2017.