# End-to-End Simultaneous Learning of Single-particle Orientation and 3D Map
Reconstruction from Cryo-electron Microscopy Data

Youssef S. G. Nashed<EMAIL_ADDRESS>0000-0001-6146-3939 Machine
Learning Initiative, SLAC National Accelerator LaboratoryMenlo
ParkCaliforniaUSA ,  Frederic Poitevin<EMAIL_ADDRESS>0000-0002-3181-8652 Department of LCLS Data Analytics, SLAC National
Accelerator LaboratoryMenlo ParkCaliforniaUSA ,  Harshit Gupta  Machine
Learning Initiative, SLAC National Accelerator LaboratoryMenlo
ParkCaliforniaUSA ,  Geoffrey Woollard  Department of Mathematics, University
of British ColumbiaVancouverCanada ,  Michael Kagan  Fundamental Physics
Directorate, SLAC National Accelerator LaboratoryMenlo ParkCaliforniaUSA ,
Chuck Yoon  Department of LCLS Data Analytics, SLAC National Accelerator
LaboratoryMenlo ParkCaliforniaUSA and  Daniel Ratner  Machine Learning
Initiative, SLAC National Accelerator LaboratoryMenlo ParkCaliforniaUSA

###### Abstract.

Cryogenic electron microscopy (cryo-EM) provides images from different copies
of the same biomolecule in arbitrary orientations. Here, we present an end-to-
end unsupervised approach that learns individual particle orientations from
cryo-EM data while reconstructing the average 3D map of the biomolecule,
starting from a random initialization. The approach relies on an auto-encoder
architecture where the latent space is explicitly interpreted as orientations
used by the decoder to form an image according to the linear projection model.
We evaluate our method on simulated data and show that it is able to
reconstruct 3D particle maps from noisy- and CTF-corrupted 2D projection
images of unknown particle orientations.

##  1\. Introduction

Determining the structure of biomolecules is an important step towards
understanding their functional mechanism or designing new drugs. Structural
determination techniques such as nuclear magnetic resonance spectroscopy,
X-ray diffraction (XRD) crystallography or cryogenic electron microscopy
(cryo-EM) have been very successful over the years (noa, 2021). Of those
techniques, cryo-EM has been used increasingly in recent years to determine
structures of many biomolecules at near-atomic resolution. This revolution has
been possible because of better hardware and software techniques (Kühlbrandt,
2014). Cryo-EM is the fastest growing technique in terms of structures
deposited in the Protein Data Bank, projected to be comparable to XRD within a
few years. While the typical resolution reported for cryo-EM structures is
still worse than those reported for XRD structures, recent developments have
proven possible to distinguish individual atoms in cryo-EM maps (Nakane et
al., 2020). Whereas XRD simultaneously measures millions of copies of a
molecule, cryo-EM is a single-particle imaging (SPI) technique, with each
cryo-EM image corresponding to a single molecular copy. While SPI
reconstruction methods must learn the conformation and orientation of
individual particles, this challenge is also an opportunity: rather than
learning an average structure, SPI methods can resolve biologically relevant
dynamics and avoid blurring from averaging over heterogeneous particles
(Ourmazd, 2019). In this paper, we discuss a new approach to solving the
single-particle orientation problem using a neural-network auto-encoder.

###  1.1. Cryo-EM 3D Reconstruction

In cryo-EM, a large number ($10^{4}-10^{7}$) of identical copies of the same
biomolecule are first frozen in a thin vitreous layer of ice. A beam of
electron then goes through this sample and is acquired by a detector. It is
typically assumed that only the phase of the incoming beam has been changed by
the electrostatic potential of the sample, not its amplitude. Since the sample
is typically prepared so the material that the electrons go through is very
thin, it is also usually assumed that the linear projection approximation can
be used to model the image formation process. The resulting image is called a
micrograph. Individual particles are located and extracted from the micrograph
as square patches of identical size $N\times N$ resulting in a particle stack
which constitutes the starting dataset for tomographic reconstruction of the
particle $N\times N\times N$ volume. Tomographic reconstruction is an ill-
posed problem that faces multiple challenges. The individual pose (orientation
and translational shift) of each particle in the stack is not known and needs
to be estimated from images that have been corrupted in two different ways.
First, the image is convolved with the point-spread function (PSF) of the
microscope which acts as a filter reducing the information content of each
image across its spectrum. Second, the main part of the signal in the image
comes from the surrounding ice which is a major source of noise, in addition
to other sources of noise such as shot noise resulting from dose fractionation
and detector response. Indeed a balance between radiation damage and getting
signal at all needs to be found during data collection. These difficulties
make cryo-EM 3D reconstruction a very challenging inverse problem.

###  1.2. Related work

Many methods have been developed to tackle cryo-EM tomographic reconstruction
(Sorzano et al., 2004; Grigorieff, 2007; Shaikh et al., 2008; Tang et al.,
2007; Hohn et al., 2007; Scheres, 2012b; de la Rosa-Trevín et al., 2016;
Punjani et al., 2017), including Bayesian formulations (Dempster et al., 1977;
Sigworth, 1998; Scheres, 2012b; Punjani et al., 2017), projection-matching
strategies (Penczek et al., 1994; Baker and Cheng, 1996), common-lines
approaches (Vainshtein and Goncharov, 1986) or other techniques (Mallick et
al., 2006; Singer et al., 2010; Wang et al., 2013; Greenberg and Shkolnisky,
2017; Pragier and Shkolnisky, 2019; Zehni et al., 2020). The Bayesian approach
has been popularized by the widely used software RELION (Scheres, 2012a) which
performs maximum-a-posteriori (MAP) optimization through Expectation-
Maximization (EM) to reconstruct a map which maximizes the likelihood of the
acquired data while still meeting some _a priori_ condition about the map. In
EM, during the Expectation-step, a conditional distribution over the poses and
shifts is estimated for each projection using the current estimate of the map.
In the next Maximization-step, these estimated distributions are used to
update the map. These two steps are performed iteratively until meeting a
convergence criterion. MAP optimization does not guarantee global convergence
and therefore its reconstruction is dependent on the quality of the initial
map. First implementations of the method suffered from their lack of
robustness to map initialization, so a competing software cryoSPARC (Punjani
et al., 2017) proposed optimizing a MAP solution _ab initio_ through
stochastic gradient descent (SGD) - however, in doing so, individual poses are
not estimated explicitly, thereby limiting the achievable resolution. The
resulting low-resolution map can then be further refined with the EM
algorithm. Because the EM reconstruction approach requires estimation of
conditional distribution on poses for each measurement, the number of
variables to determine grows directly with the number of measurements.

Neural networks (NNs) with their state-of-the-art performance on various
different inverse problems (Litjens et al., 2017; McCann et al., 2017) have
recently been introduced in the cryo-EM processing pipeline, mainly in pre-
processing steps such as denoising of the micrograph (Bepler et al., 2020) and
particle picking (Wang et al., 2016; Zhu et al., 2017; Tegunov and Cramer,
2019; Wagner et al., 2019; Bepler et al., 2019). More recently a few methods
using NNs have been proposed to solve the cryo-EM reconstruction problem
(Zhong et al., 2021; Rosenbaum et al., 2021; Gupta et al., 2020a, b; Punjani
and Fleet, 2021; Miolane et al., 2021). These methods require no prior
training and are fully unsupervised. In (Zhong et al., 2021) a modified
variational autoencoder (VAE) (Kingma and Welling, 2014) is used to
reconstruct continuous conformations of dynamic biomolecules. In (Gupta et
al., 2020a, b), generative adversarial networks (GANs) (Goodfellow et al.,
2014) are modified to reconstruct structure of biomolecules with single and
continuous conformations, respectively. In (Miolane et al., 2021) a
composition of GANs and VAE is used to find the latent variables that explain
the data, which then can be further used for reconstruction. In contrast to
these methods, (Punjani and Fleet, 2021) argues that amortized inference as
done by VAEs does not yield results that are precise enough and instead
perform direct inference of the conformational coordinate using an auto-
decoder approach. None of these methods use NNs to estimate the poses for each
projections. Instead they either use an external traditional pose estimation
routine (Zhong et al., 2021) or bypass the pose assignment step altogether
(Gupta et al., 2020a, b).

###  1.3. Contributions

In this paper, we propose a method which uses a modified auto-encoder to
simultaneously reconstruct the 3D map and estimate individual poses. Our
method is fully unsupervised, and does not require any prior training nor an
ab-initio solution. We use a combination of an encoder parameterized by an NN
and a cryo-EM physics based decoder parameterized by a learnable 3D structure.
The encoder outputs the estimated pose for a given input measurement, which is
then fed to the decoder that outputs the simulated cryo-EM projection from the
current structure. The encodings are constrained to correspond to the pose by
the physics-based decoder. The weights of the encoder and the 3D structure in
the decoder are simultaneously optimized in order to decrease the error
between the given projection (input of the encoder) and the output of the
decoder.

We leverage this auto-encoder approach to skip the EM-like step in traditional
methods for estimating poses. As opposed to an EM approach, the number of
encoder parameters is fixed and we do not require explicit estimation of pose
distribution for each measurement. Therefore, the number of variables to be
estimated in our method does not grow with the number of measurements.

##  2\. Background and Current Methods

The reconstruction problem of cryo-EM requires estimation of the structure
${\bf{x}}\in{\mathbb{R}}^{N\times N\times N}$ from the measurements
$\\{{\bf{y}}_{1},\ldots,{\bf{y}}_{M}\\}$. The acquisition of each measurement
${\bf{y}}_{m}\in{\mathbb{R}}^{N\times N}$ can be modelled as

(1)
${\bf{y}}_{m}=\underbrace{C_{{\bf{d}}_{m}}*S_{{\bf{t}}_{m}}\\{P_{\boldsymbol{\theta}_{m}}\\{{\bf{x}}\\}\\}}_{{\bf{H}}_{\boldsymbol{\varphi}_{m}}}+{\bf{n}}_{m},$

where ${\bf{n}}_{m}\in{\mathbb{R}}^{N\times N}$ is the additive noise. The
imaging operator ${\bf{H}}_{\boldsymbol{\varphi}_{m}}$ depends on the imaging
parameters
$\boldsymbol{\varphi}_{m}=(\boldsymbol{\theta}_{m},\boldsymbol{t}_{m},\boldsymbol{d}_{m})\in{\mathbb{R}}^{8}$.
The imaging operator consists of the projection operator
$P_{\boldsymbol{\theta}^{m}}$ which outputs the tomographic projection of the
structure rotated by the Euler angles
$\boldsymbol{\theta}_{m}=(\theta_{m,1},\theta_{m,1},\theta_{m,1})$; the shift
operator $S_{{\bf{t}}_{m}}$, which shifts the projection by
${\bf{t}}_{m}=(t_{m,1},t_{m,2})$ consisting of horizontal and vertical
directions, respectively; and the convolution operator $C_{{\bf{d}}_{m}}$
which corrupts the image by the contrast transfer function (CTF) with
parameters ${\bf{d}}_{m}=(d_{m,1},d_{m,2},\alpha_{m})$ that consists of major
defocus, minor defocus, and astigmatism angle, respectively.

As discussed earlier, the imaging parameters, $\boldsymbol{\theta}_{m}$ and
$\boldsymbol{t}_{m}$, are unknown for each measurement. This coupled with the
loss of information due to the CTF, and the high level of noise, makes the
reconstruction of ${\bf{x}}$ a challenging problem.

Maximum Likelihood. A naive approach to solving the reconstruction problem
would consist in searching for the unknown imaging parameters for each
measurement and a global structure which maximizes the likelihood of the
measurements. This quest can be described as

(2) $\displaystyle{\bf{x}}_{\rm
rec},\tilde{\boldsymbol{\varphi}}_{1},\ldots,\tilde{\boldsymbol{\varphi}}_{m}$
$\displaystyle=\arg\max_{{\bf{x}},\tilde{\boldsymbol{\varphi}}_{1},\ldots,\tilde{\boldsymbol{\varphi}}_{M}}\sum_{m=1}^{M}\log
p({\bf{y}}_{m}|{\bf{x}},\tilde{\boldsymbol{\varphi}}_{m}),$ (3)
$\displaystyle=\arg\max_{{\bf{x}},\tilde{\boldsymbol{\varphi}}_{1},\ldots,\tilde{\boldsymbol{\varphi}}_{M}}\sum_{m=1}^{M}\|{\bf{y}}_{m}-{\bf{H}}_{\tilde{\boldsymbol{\varphi}}_{m}}{\bf{x}}\|^{2},$

where $p({\bf{y}}_{m}|{\bf{x}},\tilde{\boldsymbol{\varphi}}_{m})$ denotes the
likelihood of the measurement ${\bf{y}}_{m}$ given the imaging parameter
$\tilde{\boldsymbol{\varphi}}_{m}$ and the structure ${\bf{x}}$. For an
independent white Gaussian noise model, it takes the form (3).

However, since the formulation (3) is highly non-convex and filled with poor
local minima (Ullrich et al., 2019), this naive approach would rarely give a
reasonable solution. A brute-force approach to solving the reconstruction
problem would follow an iterative procedure where instead of estimating all
the parameters at the same time, the structure would be updated using the
current estimates of the poses which would be updated in turn using the
current estimate of the structure. At iteration $k$, this is given by

(4) $\displaystyle\tilde{\boldsymbol{\varphi}}_{m,k}$
$\displaystyle=\arg\max_{\tilde{\boldsymbol{\varphi}}_{m}}\log
p({\bf{y}}_{m}|{\bf{x}}_{k},\tilde{\boldsymbol{\varphi}}_{m})\forall
m\in[1,\ldots,M],$ (5) $\displaystyle{\bf{x}}_{k+1}$
$\displaystyle=\arg\max_{{\bf{x}}}\log
p({\bf{y}}_{m}|{\bf{x}},\tilde{\boldsymbol{\varphi}}_{k,m})$

This approach would not guarantee to find a global optimum either and the
quality of the reconstruction would be highly dependent on the initialization.
Moreover, the high level of noise makes the pose estimation error-prone.

Maximum Marginalized Likelihood. To remedy these issues, current methods use a
marginalized likelihood formulation which instead of estimating a single pose
for each measurement, effectively weighs the contribution of the poses from
the whole search space. This is given by Eq.(6) where
$p({\bf{y}}_{m}|{\bf{x}})$ denotes the probability of the projection given the
structure.

(6) $\displaystyle{\bf{x}}_{\rm rec}$
$\displaystyle=\arg\max_{{\bf{x}}}\sum_{m=1}^{M}\log
p({\bf{y}}_{m}|{\bf{x}}),$
$\displaystyle=\arg\max_{{\bf{x}}}\sum_{m=1}^{M}\log\int
p({\bf{y}}_{m}|{\bf{x}},\boldsymbol{\varphi})p(\boldsymbol{\varphi})\,\mathrm{d}\boldsymbol{\varphi},$

Formulation (6) can be solved in two stages. First, a reasonable approximation
to the global maximum can be found through SGD or some other method. Next,
this ab-initio structure ${\bf{x}}_{0}$ is used to initialize an iterative EM
procedure (also called iterative refinement) which will refine the solution
further, at the cost of doing more operations for pose estimation.

The expectation step of the $k$-th iteration estimates a conditional
distribution on the space of poses for each projection given the current
estimate of the structure ${\bf{x}}_{k}$. This estimate is given by

(7) $\displaystyle
p(\boldsymbol{\varphi}|{\bf{x}}_{k},{\bf{y}}_{m})=\frac{p({\bf{y}}_{m}|{\bf{x}}_{k},\boldsymbol{\varphi})p(\boldsymbol{\varphi})}{\int_{\boldsymbol{\varphi}}p({\bf{y}}_{m}|{\bf{x}}_{k},\boldsymbol{\varphi})p(\boldsymbol{\varphi})\,\mathrm{d}\boldsymbol{\varphi}}.$

The maximization step then uses these poses to update the structures by
solving

(8)
$\displaystyle{\bf{x}}_{k+1}=\arg\max_{{\bf{x}}}\sum_{m=1}^{M}\mathbb{E}_{p(\boldsymbol{\varphi}|{\bf{x}}_{k},{\bf{y}}_{m})}[\log
p({\bf{y}}_{m}|{\bf{x}},\boldsymbol{\varphi})].$

For feasibility, the space of $\boldsymbol{\varphi}$ is discretized to compute
the integrals in (7) and (8). This weighted form of pose estimation is less
sensitive to the initial reference and yields better quality reconstructions.
In most methods, prior knowledge over the structure is used to obtain a
modified MAP formulation. However, the traditional approaches just described
are computationally intensive procedures because estimating poses or
conditional distributions over them against an ever changing reference
structure scales poorly - the number of variables to estimate grows directly
with the size of the data.

##  3\. Proposed Method

To solve the scaling problem, we propose a neural network based representation
of poses where we consider the unknown poses
$\tilde{\boldsymbol{\varphi}}_{m}$ as the output of a tunable function
$E_{\gamma}$ for a given input measurement. In this work, we consider the
shifts ${\bf{t}}_{m}$ and defocus parameters ${\bf{d}}_{m}$ known, but in
principle they could also be added to $E_{\gamma}$ outputs.

For a general error function $R$, we solve (10) which becomes (11) (similar to
(3)) when a Gaussian noise model is assumed

(9) $\displaystyle\tilde{\boldsymbol{\varphi}}_{m}$
$\displaystyle=E_{\gamma}({\bf{y}}_{m})\quad\forall\,m\in[1,\ldots,M].$ (10)
$\displaystyle{\bf{x}}_{\rm rec}$
$\displaystyle=\arg\max_{{\bf{x}},\gamma}\sum_{m=1}^{M}R(y_{m},{\bf{H}}_{E_{\gamma}({\bf{y}}_{m})}{\bf{x}}).$
(11)
$\displaystyle=\arg\max_{{\bf{x}},\gamma}\sum_{m=1}^{M}\|{\bf{y}}_{m}-{\bf{H}}_{E_{\gamma}({\bf{y}}_{m})}{\bf{x}}\|^{2}.$

We minimize (10) using SGD as described in Algorithm 1. At each iteration, for
a batch of measurements $\\{{\bf{y}}_{1},\ldots,{\bf{y}}_{B}\\}$, we get the
empirical estimate of the loss function in (10) by

(12) $\displaystyle
L({\bf{x}},\gamma)=\sum_{b=1}^{B}R(y_{b},{\bf{H}}_{E_{\gamma}({\bf{y}}_{b})}{\bf{x}}).$

The structure and the encoder weights are optimized using the gradients
$\nabla_{{\bf{x}}}L({\bf{x}},\gamma)$ and $\nabla_{\gamma}L({\bf{x}},\gamma)$,
respectively. In summary, simultaneously with structure estimation we tune the
weights of a neural network so that the imaging parameters maximize the
likelihood of the measurements. Our scheme is shown in Figure 1.

Algorithm 1 Reconstruct cryo-EM data

Input: acquired dataset $\\{{\bf{y}}_{1},\ldots,{\bf{y}}_{M}\\}$ ; number of
reconstruction iterations, $n_{\mathrm{rec}}$; size of the batches used for
SGD, $B$; optimizer parameters; Error function $R$;
Initialization: ${\bf{x}}_{\rm rec}$ with a uniform random distribution, and
$E_{\gamma}$ with random weights

1:for $n_{\mathrm{rec}}$ do

2: From acquired dataset, sample a batch$\\{{\bf{y}}^{1}_{\rm
batch},\ldots,{\bf{y}}_{\rm batch}^{B}\\}$

3: Obtain $\tilde{\boldsymbol{\varphi}}_{b}=E_{\gamma}({\bf{y}}_{\rm
batch}^{b})\forall b\in[1,\ldots,B]$

4: Compute $L({\bf{x}},\gamma)=\sum_{b=1}^{B}R({\bf{y}}_{\rm
batch}^{b},{\bf{H}}_{\tilde{\boldsymbol{\varphi}}_{b}}{\bf{x}}_{\rm rec})$

5: Update ${\bf{x}}_{\rm rec}$ using $\nabla_{{\bf{x}}}L({\bf{x}},\gamma)$

6: update $\gamma$ using $\nabla_{\gamma}L_{(}{\bf{x}},\gamma)$

Output: ${\bf{x}}_{\rm
rec},\tilde{\boldsymbol{\varphi}}_{1},\ldots,\tilde{\boldsymbol{\varphi}}_{M}.$

Figure 1. Auto-encoder architecture.The input image ${\bf{y}}$ is encoded into
a pose $\boldsymbol{\varphi}$ by the encoder. The decoder outputs the
tomographic projection of the structure oriented at $\varphi$, and then
corrupted by the CTF. The image formation is implemented using the Fourier-
slice theorem, with the estimated 3D structure ${\bf{x}}$ Fourier transformed,
sliced into a plane normal to $\boldsymbol{\varphi}$, and corrupted by a known
CTF with parameter ${\bf{d}}$. The parameters of the encoder and decoder are
respectively the weights $\boldsymbol{\gamma}$ of the CNN and MLP and the
structure ${\bf{x}}$. They are optimized during training through minimization
of the loss $L$ that measures the dissimilarity between the input image
${\bf{y}}$ and the reconstructed one ${\bf{H}}_{\varphi}{\bf{x}}$.

Encoder. The encoder NN is composed of a convolutional neural network (CNN)
followed by a multilayer percepetron (MLP), also known as a fully-connected
layer. The CNN extracts shift-invariant features from the measurements, which
are then transformed by the MLP into orientation parameters as described
below. We use a standard CNN encoder architecture of a conv2D-¿conv2D-¿maxpool
block repeated three times. All convolution layers use $3\times 3$ filters,
and maxpooling downsamples by a factor of $2$. The number of convolution
filters per block is $32,64,128$, respectively, and the MLP is composed of two
layers, each containing $512$ units/neurons.

Differentiable cryo-EM physics model. The cryo-EM physics operator
${\bf{H}}_{\boldsymbol{\varphi}}$ uses the output of the encoder as the pose
parameters ${\boldsymbol{\varphi}}$. In order to compute the gradients
$\nabla_{{\bf{x}}}L({\bf{x}},\gamma)$ and $\nabla_{\gamma}L({\bf{x}},\gamma)$,
the operator ${\bf{H}}_{\boldsymbol{\varphi}}$ needs to be differentiable with
respect to the structure as well as the imaging parameters. We therefore,
implement a differentiable cryo-EM physics model which enables the learning of
encoder weights using backpropagation. The image formation model is
implemented in reciprocal space, in which the the 3D Fourier transform of
${\bf{x}}$ is computed once per batch. Predicted/Decoded projections are
generated by the 2D inverse Fourier transform of a slice in the 3D Fourier
volume extracted using the corresponding predicted pose. Since the encoder is
optimized using backpropagation, the parameterization of the poses directly
affects the encoder’s performance and thereby, the quality of the
reconstruction. We use the following parameterizations in our experiments:

  * •

Euler Angles: For each measurement, the encoder yields a 3-dimensional output
which is then fed to the imaging model as Euler angles (commonly used in cryo-
EM software). The structure is then rotated using these angles.

  * •

Quaternions: The three-dimensional MLP output is transformed into a
4-dimensional vector using a fixed transformation which has the properties of
a unit quaternion ${\bf{q}}=(q_{1},q_{2},q_{3},q_{4})$. The imaging model then
rotates the structure by $2\cos^{-1}{q_{1}}$, around the axis
$(q_{2},q_{3},q_{4})$.

  * •

$\mathcal{S}^{2}\times\mathcal{S}^{2}$ (s2s2): The MLP outputs two
3-dimensional vectors. These are then orthonormalized to obtain ${\bf{w}}_{1}$
and ${\bf{w}}_{2}$. A third vector is then computed by the cross product of
the first two, ${\bf{w}}_{3}={\bf{w}}_{1}\times{\bf{w}}_{2}$. These three
vectors define a local coordinate system that relates to a $3\times 3$
rotation matrix. This matrix is then used to rotate the structure. Formally,
this case is an algebraic parameterization of the Lie group of 3D rotations
SO(3), using two orthonormal vectors ${\bf{w}}_{1}$ and ${\bf{w}}_{2}$.

NN pose estimation offers many advantages. First, instead of searching for
individual independent poses $\\{\boldsymbol{\varphi}_{m}\\}$, we estimate a
global function that maps the measurements to their respective poses. Since
this mapping is global, it indirectly integrates the information from all the
projections for each pose estimation. Second, since the information of poses
has been condensed in the fixed size $\gamma$, the number of variables to
estimate does not grow with the size of the data. Finally, the universal
approximation property of the neural networks lets us learn complicated
mappings between projections and their poses.

Simulation of datasets. All the experiments presented here used datasets
generated from the atomic model 4AKE111https://www.rcsb.org/structure/4AKE of
_E.coli_ adenylate kinase (Müller et al., 1996), a small 47 kDa protein. TEM
simulator (Rullgård et al., 2011) was used to generate a $128\times 128\times
128$ electrostatic potential map with pixel size 0.8 Å. Images associated with
pose $\varphi$ were obtained by rotating the centered map with $\varphi$ and
projecting it along the $z$-direction in real space after resampling on the
original grid. To account for CTF corruption, the resulting image was Fourier
transformed, padded to double length, multiplied by a pre-computed 2D CTF
image with given defocus $d$. Gaussian white noise is added _a posteriori_.
Datasets are generated by sampling SO(3) uniformally 10,000 times, amounting
to 9,000 images for training, and 1,000 for testing. When relevant, the CTF
range [0.4 $\upmu$m, 1.2 $\upmu$m] is sampled uniformly. The noise is sampled
from the normal distribution and scaled to match the desired signal-to-noise
ratio (SNR). Table 1 lists the datasets generated.

Dataset | SNR (dB) | CTF | size
---|---|---|---
A | - | - | $9,000\times 128\times 128$
B, C, D, E, F | [10, 5, 0, -5, -10] | - | $9,000\times 128\times 128$
G | 0 | [0.4 $\upmu$m, 1.2 $\upmu$m] | $9,000\times 128\times 128$
Table 1. Datasets used in the numerical experiments

Training. For all the described datasets, training was carried with the
following parameters: we use the Adam (Adaptive Moment Estimation) (Kingma and
Ba, 2017) optimizer with a learning rate of
$5\text{\times}{10}^{-10}\text{\,}\mathrm{,}$ minibatch size $B=32$, number of
minibatch update steps $n_{\mathrm{rec}}=10,000$ which is equivalent to
approximately $35.5$ epochs or passes over the whole dataset. The framework is
implemented using Tensorflow and runs on an NVIDIA Tesla V100 GPU in around 5
hours for a full training run.

##  4\. Results

We tested our method on simulated datasets. We detail below how the methods
performed in the absence of noise, in the presence of increasing noise, and in
a more realistic setting with noise and CTF corruption.

SO(3) parameterization. In the absence of noise (dataset A), training
converged in a few thousand steps, as can be seen from the loss curve on
Fig.2-top. Following a similar convergence pattern, the MAE between estimated
and ground truth poses $\varphi$ decreased drastically in the first few
hundred steps (see Fig.2-bottom). Three different approach to SO(3)
parameterization have been tried, of which the Euler representation converges
the slowest and s2s2 the fastest and ultimately most accurate. On close
inspection (data not shown), the main reason for this discrepancy is
attributed to parameterization singularities arising from both the Euler and
quaternion formulations that could affect the numerical stability of gradients
computed around those singularities. For a more in depth proof of those
formulations the reader is referred to these publications (Falorsi et al.,
2018; Kobayashi and Nomizu, [n.d.]). In all the subsequent experiments we
adopt s2s2 as the chosen orientation parameterization.

Figure 2. SO(3) parameterization and convergence. (top) L2 loss convergence
for various representations of SO(3): Euler angles (euler), quaternions (quat)
and s2s2. Each corresponds to 10 independent training runs (average and
standard deviation). For comparison, the loss curves obtained through
automatic differentiation (AD) or with known poses (tomo) are also shown.
(bottom) Evolution of the mean absolute error (MAE) over the poses.

Comparison to Automatic differentiation. To provide points of comparisons with
previous work, we considered the following two scenarios. First, we monitored
the performance of our approach when poses were known (see $tomo$ curve in
Fig.2). In this case, as expected, the network converges almost instantly to
what we consider a lower bound for the loss. We also considered the case where
the poses and map were solved through stochastic gradient descent, similar to
(Punjani et al., 2017). In this case, the loss hardly decreases over a few
thousand steps of optimization, mainly due to attempting to solve for
orientation and structure simultaneously, which is difficult without careful
initialization of the solved variables or alternating the updates between
them. Additionally, as mentioned earlier, the number of orientation variables
to be solved by AD scales with the dataset size.

Reconstruction quality. To measure the reconstruction quality of the estimated
map ${\bf{x}}$, we compute the Fourier shell correlation (FSC) every 200 steps
between the current and the ground truth map (see Fig.3-left). In the first
few hundred steps, the FSC curve intersects the 0.143 cutoff at resolutions
worse than 3.2 Å and the Nyquist-Shannon limit of 1.6 Å is reached after 600
steps. The FSC keeps improving moderately until it reaches a stable curve.

Figure 3. Fourier-Shell Correlation to ground truth. The dotted red line
indicates the FSC threshold used to estimate the map resolution. (left) FSC
curves at regular intervals, from black to white, during optimization in the
absence of noise (dataset A). (center) FSC curves for final models at various
noise levels (datasets B-F), in the absence of CTF corruption. (right) FSC
curves for final models at zero SNR with (dataset D) and without CTF (dataset
G) corruption of the dataset.

Effect of noise on reconstruction. We tested the ability of our method to
handle realistic noise typically encountered in cryoEM data. Figure 3-center
summarizes our findings. As expected, adding noise is detrimental to the
quality of the reconstruction. Yet we see that for the relatively small
dataset size that we tried, the Nyquist-Shannon limit is still hit at SNR
levels of 0 dB. At more realistic value of -10 dB we see that the effective
resolution of the reconstruction is still better than 2.5 Å. Visual inspection
of the reconstructed maps (see Fig.4) is consistent with this observation.

Effect of CTF. Finally, we tested the ability of our method to reconstruct
images that were both degraded by added noise (0 dB) and by CTF corruption
(dataset G). The FSC curve of the resulting reconstruction is very similar to
the one obtained without CTF corruption, both hitting the Nyquist-Shannon
limit (see Fig.3-right). However, the shape differs slightly, with higher
correlation in the highest resolution shells and lower correlation in the
medium-resolution range.Visual inspection of the resulting map is consistent
with this observation (see Fig.4).

Figure 4. Reconstructions. The left column shows samples from their respective
datasets (A-top, F-center, G-bottom). The center column shows the
reconstructed image after training. The right column shows the reconstructed
map at the same contour level.

##  5\. Conclusion

We have presented a new method for solving the pose estimation problem in
cryoEM using an autoencoder architecture employing a differentiable forward
model decoder. We have shown, on simulated data, that the method was able to
reconstruct high-resolution 3D map from simulated data in the presence of
noise or CTF corruption and without any prior structural knowledge. We have
recently become aware of contemporaneous work using a VAE for estimation of
pose and conformation(Rosenbaum et al., 2021). While this approach also uses
NNs to estimate poses, it shares with (Zhong et al., 2021; Punjani and Fleet,
2021) the limitation that it relies on a consensus structure initialization,
the quality of which might impact the ability of the methods to reliably learn
poses and conformations, in particular for highly flexible systems. We
hypothetize that the methods presented here would provide a solution to this
issue. Because the computational cost is independent of the number of images,
we also envision this method to be the basis of a new pipeline for the large
datasets expected to be needed to resolve continuous conformations.

###### Acknowledgements.

We thank Wah Chiu, Khanh Dao Duc, Mark Hunter, TJ Lane, Julien Martel, Nina
Miolane, and Gordon Wetzstein for numerous discussions that helped shape this
project. This work was supported by the U.S. Department of Energy, under DOE
Contract No. DE-AC02-76SF00515. We acknowledge the use of the computational
resources at the SLAC Shared Scientific Data Facility (SDF).

## References

  * (1)
  * noa (2021) 2021\.  A celebration of structural biology.  _Nature Methods_ 18, 5 (May 2021), 427–427.  https://doi.org/10.1038/s41592-021-01163-x
  * Baker and Cheng (1996) Timothy S. Baker and R. Holland Cheng. 1996.  A Model-Based Approach for Determining Orientations of Biological Macromolecules Imaged by Cryoelectron Microscopy.  _Journal of Structural Biology_ 116, 1 (Jan. 1996), 120–130.  https://doi.org/10.1006/jsbi.1996.0020
  * Bepler et al. (2020) Tristan Bepler, Kotaro Kelley, Alex J. Noble, and Bonnie Berger. 2020.  Topaz-Denoise: general deep denoising models for cryoEM and cryoET.  _Nat Commun_ 11, 1 (Oct. 2020), 5208\.  https://doi.org/10.1038/s41467-020-18952-1 Number: 1 Publisher: Nature Publishing Group.
  * Bepler et al. (2019) Tristan Bepler, Andrew Morin, Micah Rapp, Julia Brasch, Lawrence Shapiro, Alex J. Noble, and Bonnie Berger. 2019.  Positive-unlabeled convolutional neural networks for particle picking in cryo-electron micrographs.  _Nat Methods_ 16, 11 (Nov. 2019), 1153–1160.  https://doi.org/10.1038/s41592-019-0575-8 Number: 11 Publisher: Nature Publishing Group.
  * de la Rosa-Trevín et al. (2016) J. M. de la Rosa-Trevín, A. Quintana, L. del Cano, A. Zaldívar, I. Foche, J. Gutiérrez, J. Gómez-Blanco, J. Burguet-Castell, J. Cuenca-Alba, V. Abrishami, J. Vargas, J. Otón, G. Sharov, J. L. Vilas, J. Navas, P. Conesa, M. Kazemi, R. Marabini, C. O. S. Sorzano, and J. M. Carazo. 2016.  Scipion: A software framework toward integration, reproducibility and validation in 3D electron microscopy.  _Journal of Structural Biology_ 195, 1 (July 2016), 93–99.  https://doi.org/10.1016/j.jsb.2016.04.010
  * Dempster et al. (1977) A. P. Dempster, N. M. Laird, and D. B. Rubin. 1977\.  Maximum Likelihood from Incomplete Data Via the EM Algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ 39, 1 (Sept. 1977), 1–22.  https://doi.org/10.1111/j.2517-6161.1977.tb01600.x
  * Falorsi et al. (2018) Luca Falorsi, Pim de Haan, Tim R. Davidson, Nicola De Cao, Maurice Weiler, Patrick Forré, and Taco S. Cohen. 2018.  Explorations in Homeomorphic Variational Auto-Encoding.  (July 2018).  https://arxiv.org/abs/1807.04689v1
  * Goodfellow et al. (2014) Ian J. Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. 2014\.  Generative Adversarial Networks.  _arXiv:1406.2661 [cs, stat]_ (June 2014).  http://arxiv.org/abs/1406.2661 arXiv: 1406.2661.
  * Greenberg and Shkolnisky (2017) Ido Greenberg and Yoel Shkolnisky. 2017.  Common lines modeling for reference free Ab-initio reconstruction in cryo-EM.  _Journal of Structural Biology_ 200, 2 (Nov. 2017), 106–117.  https://doi.org/10.1016/j.jsb.2017.09.007
  * Grigorieff (2007) Nikolaus Grigorieff. 2007\.  FREALIGN: High-resolution refinement of single particle structures.  _Journal of Structural Biology_ 157, 1 (Jan. 2007), 117–125.  https://doi.org/10.1016/j.jsb.2006.05.004
  * Gupta et al. (2020a) Harshit Gupta, Michael T. McCann, Laurène Donati, and Michael Unser. 2020a.  CryoGAN: A New Reconstruction Paradigm for Single-Particle Cryo-EM via Deep Adversarial Learning.  _bioRxiv_ (July 2020), 2020.03.20.001016.  https://doi.org/10.1101/2020.03.20.001016 Publisher: Cold Spring Harbor Laboratory Section: New Results.
  * Gupta et al. (2020b) Harshit Gupta, Thong H. Phan, Jaejun Yoo, and Michael Unser. 2020b.  Multi-CryoGAN: Reconstruction of Continuous Conformations in Cryo-EM Using Generative Adversarial Networks. In _Computer Vision – ECCV 2020 Workshops_ _(Lecture Notes in Computer Science)_ , Adrien Bartoli and Andrea Fusiello (Eds.). Springer International Publishing, Cham, 429–444.  https://doi.org/10.1007/978-3-030-66415-2_28
  * Hohn et al. (2007) Michael Hohn, Grant Tang, Grant Goodyear, P. R. Baldwin, Zhong Huang, Pawel A. Penczek, Chao Yang, Robert M. Glaeser, Paul D. Adams, and Steven J. Ludtke. 2007\.  SPARX, a new environment for Cryo-EM image processing.  _Journal of Structural Biology_ 157, 1 (Jan. 2007), 47–55.  https://doi.org/10.1016/j.jsb.2006.07.003
  * Kingma and Ba (2017) Diederik P. Kingma and Jimmy Ba. 2017.  Adam: A Method for Stochastic Optimization.  _arXiv:1412.6980 [cs]_ (Jan. 2017).  http://arxiv.org/abs/1412.6980 arXiv: 1412.6980.
  * Kingma and Welling (2014) Diederik P. Kingma and Max Welling. 2014.  Auto-Encoding Variational Bayes.  _arXiv:1312.6114 [cs, stat]_ (May 2014).  http://arxiv.org/abs/1312.6114 arXiv: 1312.6114.
  * Kobayashi and Nomizu ([n.d.]) Shoshichi Kobayashi and Katsumi Nomizu. [n.d.].  _Foundations of Differential Geometry, 2 Volume Set | Wiley_.  https://www.wiley.com/en-us/Foundations+of+Differential+Geometry%2C+2+Volume+Set-p-9780470555583
  * Kühlbrandt (2014) Werner Kühlbrandt. 2014\.  The Resolution Revolution.  _Science_ 343, 6178 (March 2014), 1443–1444.  https://doi.org/10.1126/science.1251652 Publisher: American Association for the Advancement of Science Section: Perspective.
  * Litjens et al. (2017) Geert Litjens, Thijs Kooi, Babak Ehteshami Bejnordi, Arnaud Arindra Adiyoso Setio, Francesco Ciompi, Mohsen Ghafoorian, Jeroen A. W. M. van der Laak, Bram van Ginneken, and Clara I. Sánchez. 2017.  A survey on deep learning in medical image analysis.  _Medical Image Analysis_ 42 (Dec. 2017), 60–88.  https://doi.org/10.1016/j.media.2017.07.005
  * Mallick et al. (2006) S.P. Mallick, S. Agarwal, D.J. Kriegman, S.J. Belongie, B. Carragher, and C.S. Potter. 2006\.  Structure and View Estimation for Tomographic Reconstruction: A Bayesian Approach. In _2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’06)_ , Vol. 2. 2253–2260.  https://doi.org/10.1109/CVPR.2006.295 ISSN: 1063-6919.
  * McCann et al. (2017) Michael T. McCann, Kyong Hwan Jin, and Michael Unser. 2017\.  Convolutional Neural Networks for Inverse Problems in Imaging: A Review.  _IEEE Signal Processing Magazine_ 34, 6 (Nov. 2017), 85–95.  https://doi.org/10.1109/MSP.2017.2739299 Conference Name: IEEE Signal Processing Magazine.
  * Miolane et al. (2021) Nina Miolane, Frédéric Poitevin, Yee-Ting Li, and Susan Holmes. 2021.  Estimation of Orientation and Camera Parameters from Cryo-Electron Microscopy Images with Variational Autoencoders and Generative Adversarial Networks.  _arXiv:1911.08121 [cs, eess, q-bio, stat]_ (May 2021).  http://arxiv.org/abs/1911.08121 arXiv: 1911.08121.
  * Müller et al. (1996) C. W. Müller, G. J. Schlauderer, J. Reinstein, and G. E. Schulz. 1996.  Adenylate kinase motions during catalysis: an energetic counterweight balancing substrate binding.  _Structure_ 4, 2 (Feb. 1996), 147–156.  https://doi.org/10.1016/s0969-2126(96)00018-4
  * Nakane et al. (2020) Takanori Nakane, Abhay Kotecha, Andrija Sente, Greg McMullan, Simonas Masiulis, Patricia M. G. E. Brown, Ioana T. Grigoras, Lina Malinauskaite, Tomas Malinauskas, Jonas Miehling, Tomasz Uchański, Lingbo Yu, Dimple Karia, Evgeniya V. Pechnikova, Erwin de Jong, Jeroen Keizer, Maarten Bischoff, Jamie McCormack, Peter Tiemeijer, Steven W. Hardwick, Dimitri Y. Chirgadze, Garib Murshudov, A. Radu Aricescu, and Sjors H. W. Scheres. 2020.  Single-particle cryo-EM at atomic resolution.  _Nature_ 587, 7832 (Nov. 2020), 152–156.  https://doi.org/10.1038/s41586-020-2829-0 Number: 7832 Publisher: Nature Publishing Group.
  * Ourmazd (2019) Abbas Ourmazd. 2019\.  Cryo-EM, XFELs and the structure conundrum in structural biology.  _Nat Methods_ 16, 10 (Oct. 2019), 941–944.  https://doi.org/10.1038/s41592-019-0587-4 Number: 10 Publisher: Nature Publishing Group.
  * Penczek et al. (1994) Pawel A. Penczek, Robert A. Grassucci, and Joachim Frank. 1994.  The ribosome at improved resolution: New techniques for merging and orientation refinement in 3D cryo-electron microscopy of biological particles.  _Ultramicroscopy_ 53, 3 (March 1994), 251–270.  https://doi.org/10.1016/0304-3991(94)90038-8
  * Pragier and Shkolnisky (2019) Gabi Pragier and Yoel Shkolnisky. 2019.  A common lines approach for ab-initio modeling of cyclically-symmetric molecules.  _Inverse Problems_ 35, 12 (Dec. 2019), 124005\.  https://doi.org/10.1088/1361-6420/ab2fb2 arXiv: 1901.10888.
  * Punjani and Fleet (2021) Ali Punjani and David J. Fleet. 2021.  3D Flexible Refinement: Structure and Motion of Flexible Proteins from Cryo-EM.  _bioRxiv_ (April 2021), 2021.04.22.440893.  https://doi.org/10.1101/2021.04.22.440893 Publisher: Cold Spring Harbor Laboratory Section: New Results.
  * Punjani et al. (2017) Ali Punjani, John L. Rubinstein, David J. Fleet, and Marcus A. Brubaker. 2017.  cryoSPARC: algorithms for rapid unsupervised cryo-EM structure determination.  _Nat Methods_ 14, 3 (March 2017), 290–296.  https://doi.org/10.1038/nmeth.4169 Number: 3 Publisher: Nature Publishing Group.
  * Rosenbaum et al. (2021) Dan Rosenbaum, Marta Garnelo, Michal Zielinski, Charlie Beattie, Ellen Clancy, Andrea Huber, Pushmeet Kohli, Andrew W. Senior, John Jumper, Carl Doersch, S. M. Ali Eslami, Olaf Ronneberger, and Jonas Adler. 2021.  Inferring a Continuous Distribution of Atom Coordinates from Cryo-EM Images using VAEs.  _arXiv:2106.14108 [cs, eess]_ (June 2021).  http://arxiv.org/abs/2106.14108 arXiv: 2106.14108.
  * Rullgård et al. (2011) H. Rullgård, L.-G. Öfverstedt, S. Masich, B. Daneholt, and O. Öktem. 2011.  Simulation of transmission electron microscope images of biological specimens.  _Journal of Microscopy_ 243, 3 (2011), 234–256.  https://doi.org/10.1111/j.1365-2818.2011.03497.x _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1365-2818.2011.03497.x.
  * Scheres (2012a) Sjors H.W. Scheres. 2012a.  A Bayesian View on Cryo-EM Structure Determination.  _J Mol Biol_ 415, 2-4 (Jan. 2012), 406–418.  https://doi.org/10.1016/j.jmb.2011.11.010
  * Scheres (2012b) Sjors H. W. Scheres. 2012b.  RELION: Implementation of a Bayesian approach to cryo-EM structure determination.  _Journal of Structural Biology_ 180, 3 (Dec. 2012), 519–530.  https://doi.org/10.1016/j.jsb.2012.09.006
  * Shaikh et al. (2008) Tanvir R. Shaikh, Haixiao Gao, William T. Baxter, Francisco J. Asturias, Nicolas Boisset, Ardean Leith, and Joachim Frank. 2008.  SPIDER image processing for single-particle reconstruction of biological macromolecules from electron micrographs.  _Nat Protoc_ 3, 12 (Dec. 2008), 1941–1974.  https://doi.org/10.1038/nprot.2008.156 Number: 12 Publisher: Nature Publishing Group.
  * Sigworth (1998) F. J. Sigworth. 1998\.  A Maximum-Likelihood Approach to Single-Particle Image Refinement.  _Journal of Structural Biology_ 122, 3 (Jan. 1998), 328–339.  https://doi.org/10.1006/jsbi.1998.4014
  * Singer et al. (2010) Amit Singer, Ronald R. Coifman, Fred J. Sigworth, David W. Chester, and Yoel Shkolnisky. 2010\.  Detecting consistent common lines in cryo-EM by voting.  _Journal of Structural Biology_ 169, 3 (March 2010), 312–322.  https://doi.org/10.1016/j.jsb.2009.11.003
  * Sorzano et al. (2004) C. O. S. Sorzano, R. Marabini, J. Velázquez-Muriel, J. R. Bilbao-Castro, S. H. W. Scheres, J. M. Carazo, and A. Pascual-Montano. 2004\.  XMIPP: a new generation of an open-source image processing package for electron microscopy.  _Journal of Structural Biology_ 148, 2 (Nov. 2004), 194–204.  https://doi.org/10.1016/j.jsb.2004.06.006
  * Tang et al. (2007) Guang Tang, Liwei Peng, Philip R. Baldwin, Deepinder S. Mann, Wen Jiang, Ian Rees, and Steven J. Ludtke. 2007.  EMAN2: An extensible image processing suite for electron microscopy.  _Journal of Structural Biology_ 157, 1 (Jan. 2007), 38–46.  https://doi.org/10.1016/j.jsb.2006.05.009
  * Tegunov and Cramer (2019) Dimitry Tegunov and Patrick Cramer. 2019.  Real-time cryo-electron microscopy data preprocessing with Warp.  _Nat Methods_ 16, 11 (Nov. 2019), 1146–1152.  https://doi.org/10.1038/s41592-019-0580-y Number: 11 Publisher: Nature Publishing Group.
  * Ullrich et al. (2019) Karen Ullrich, Rianne van den Berg, Marcus Brubaker, David Fleet, and Max Welling. 2019\.  Differentiable probabilistic models of scientific imaging with the Fourier slice theorem.  _arXiv:1906.07582 [cs, eess, stat]_ (June 2019).  http://arxiv.org/abs/1906.07582 arXiv: 1906.07582.
  * Vainshtein and Goncharov (1986) B. K. Vainshtein and A. B. Goncharov. 1986.  Determination of the spatial orientation of arbitrarily arranged identical particles of unknown structure from their projections.  _Soviet Physics Doklady_ 31 (April 1986), 278\.  http://adsabs.harvard.edu/abs/1986SPhD...31..278V
  * Wagner et al. (2019) Thorsten Wagner, Felipe Merino, Markus Stabrin, Toshio Moriya, Claudia Antoni, Amir Apelbaum, Philine Hagel, Oleg Sitsel, Tobias Raisch, Daniel Prumbaum, Dennis Quentin, Daniel Roderer, Sebastian Tacke, Birte Siebolds, Evelyn Schubert, Tanvir R. Shaikh, Pascal Lill, Christos Gatsogiannis, and Stefan Raunser. 2019.  SPHIRE-crYOLO is a fast and accurate fully automated particle picker for cryo-EM.  _Commun Biol_ 2, 1 (June 2019), 1–13.  https://doi.org/10.1038/s42003-019-0437-z Number: 1 Publisher: Nature Publishing Group.
  * Wang et al. (2016) Feng Wang, Huichao Gong, Gaochao Liu, Meijing Li, Chuangye Yan, Tian Xia, Xueming Li, and Jianyang Zeng. 2016\.  DeepPicker: A deep learning approach for fully automated particle picking in cryo-EM.  _Journal of Structural Biology_ 195, 3 (Sept. 2016), 325–336.  https://doi.org/10.1016/j.jsb.2016.07.006
  * Wang et al. (2013) Lanhui Wang, Amit Singer, and Zaiwen Wen. 2013.  Orientation Determination of Cryo-EM Images Using Least Unsquared Deviations.  _SIAM J. Imaging Sci._ 6, 4 (Jan. 2013), 2450–2483.  https://doi.org/10.1137/130916436 Publisher: Society for Industrial and Applied Mathematics.
  * Zehni et al. (2020) Mona Zehni, Laurène Donati, Emmanuel Soubies, Zhizhen J. Zhao, and Michael Unser. 2020\.  Joint Angular Refinement and Reconstruction for Single-Particle Cryo-EM.  _IEEE Trans. on Image Process._ 29 (2020), 6151–6163.  https://doi.org/10.1109/TIP.2020.2984313 arXiv: 2003.10062.
  * Zhong et al. (2021) Ellen D. Zhong, Tristan Bepler, Bonnie Berger, and Joseph H. Davis. 2021.  CryoDRGN: reconstruction of heterogeneous cryo-EM structures using neural networks.  _Nat Methods_ 18, 2 (Feb. 2021), 176–185.  https://doi.org/10.1038/s41592-020-01049-4 Number: 2 Publisher: Nature Publishing Group.
  * Zhu et al. (2017) Yanan Zhu, Qi Ouyang, and Youdong Mao. 2017.  A deep convolutional neural network approach to single-particle recognition in cryo-electron microscopy.  _BMC Bioinformatics_ 18, 1 (July 2017), 348\.  https://doi.org/10.1186/s12859-017-1757-y


# A Mass Model for the Lensing Cluster SDSS J1004+4112: Constraints From the
Third Time Delay

R. Forés-Toribio Departamento de Astronomía y Astrofísica, Universidad de
Valencia, E-46100 Burjassot, Valencia, Spain  Observatorio Astronómico,
Universidad de Valencia, E-46980 Paterna, Valencia, Spain  J. A. Muñoz
Departamento de Astronomía y Astrofísica, Universidad de Valencia, E-46100
Burjassot, Valencia, Spain  Observatorio Astronómico, Universidad de Valencia,
E-46980 Paterna, Valencia, Spain  C. S. Kochanek Department of Astronomy, The
Ohio State University, 140 West 18th Avenue, Columbus, OH 43210, USA  Center
for Cosmology and Astroparticle Physics, The Ohio State University, 191 West
Woodruff Avenue, Columbus, OH 43210, USA  E. Mediavilla Instituto de
Astrofísica de Canarias, Vía Láctea S/N, La Laguna, E-38200, Tenerife, Spain
Departamento de Astrofísica, Universidad de la Laguna, La Laguna, E-38200,
Tenerife, Spain Raquel Forés-Toribio<EMAIL_ADDRESS>
###### Abstract

We have built a new model for the lens system SDSS J1004+4112 including the
recently measured time delay of the fourth quasar image. This time delay has a
strong influence on the inner mass distribution of the lensing cluster
($\rho\propto r^{-\alpha}$) allowing us to determine
$\alpha=1.18^{+0.02(+0.11)}_{-0.03(-0.18)}$ at the 68% (95%) confidence level
in agreement with hydrodynamical simulations of massive galaxy clusters. We
find an offset between the brightest cluster galaxy (BCG) and the dark matter
halo of $3.8^{+0.6(+1.4)}_{-0.7(-1.3)}$ kpc at 68% (95%) confidence which is
compatible with other galaxy cluster measurements. As an observational
challenge, the estimated time delay between the leading image C and the faint
(I=24.7) fifth image E is roughly 8 yr.

Gravitational lensing (670): Strong gravitational lensing (1643) — Galaxy
clusters (584) — Dark matter distribution (356) — Quasars (1319): individual:
SDSS J1004+4112

##  1 Introduction

The lens system SDSS J1004+4112 was the first example of a quasar lensed by a
cluster (Inada et al., 2003). The source and lens redshifts are $z_{s}=1.73$
and $z_{l}=0.68$, respectively. This leads to an image separation of 15″ and
large magnified images of the host galaxy. It is the second largest quasar
lens after SDSS J1029+2623 (Inada et al., 2006). Since its discovery, many
additional observational constraints have been obtained: more galaxy cluster
members (Oguri et al., 2004), the central quasar image (Inada et al., 2005,
2008), other background lensed sources (Sharon et al., 2005; Liesenborgs et
al., 2009; Oguri, 2010), the velocity dispersion of the brightest galaxy
cluster (Inada et al., 2008) and the time delays between quasar images A, B,
and C (Fohlmeister et al., 2007, 2008).

As new observational constraints were measured, increasingly sophisticated
mass models of the lens cluster were made using both parametric (Inada et al.,
2003; Oguri et al., 2004; Fohlmeister et al., 2007; Oguri, 2010) and
nonparametric (Williams & Saha, 2004; Saha et al., 2006; Liesenborgs et al.,
2009; Mohammed et al., 2015) approaches. Recently, the time delay of the
fourth quasar image (Muñoz et al., 2022) has been measured at 2457 days
between images D and C, and it is over one year longer than predicted (Oguri,
2010). Indeed, the ability of the models to predict the unmeasured time delays
has been remarkably poor.

Apart from the models of the mass distribution, this system has been widely
studied not only in the optical but also from X-rays to radio (Lamer et al.,
2006; Ota et al., 2006; Ross et al., 2009; Jackson, 2011; Chen et al., 2012;
McKean et al., 2021). However, there are still issues that remain unsolved
like the origin of the enhanced asymmetric wings in the broad emission lines
of quasar image A (Richards et al., 2004; Gómez-Álvarez et al., 2006; Lamer et
al., 2006; Motta et al., 2012; Fian et al., 2018; Popović et al., 2020; Fian
et al., 2021).

Here we revisit the mass model making use of the newly measured time delay.
The new delay should help to constrain the inner mass profile of the cluster
halo and to look for deviations from a Navarro-Frenk-White (NFW) profile
(Kawano & Oguri, 2006; Oguri, 2010). The paper is organized as follows. In
Section 2 we present the available observational constraints. In Section 3 we
describe the modeling method and the components of the mass model. Section 4
presents the results and they are discussed in Section 5. Throughout this
paper we assume a flat $\Lambda$CDM cosmology with $\Omega_{M}=0.26$,
$\Omega_{\Lambda}=0.74$ and $H_{0}$=72 km s${}^{-}1$ Mpc-1.

##  2 Observational Data

In addition to the lensed quasar, SDSS J1004+4112 has seven lensed background
galaxies at three different redshifts. Their redshifts, image positions, and
magnitudes (with respect to quasar image A) are from Oguri (2010) and listed
in Table 1. For the quasar time delays we use the results of Muñoz et al.
(2022) with the uncertainties set to five times the formal errors (see Table
1) since time delay uncertainties are almost always underestimated (see, e.g.,
Tie & Kochanek, 2018). We also adopted the positions, ellipticities, position
angles, and luminosity ratios (with respect to the central galaxy) of 14
cluster galaxy members (see Table 2) from Oguri (2010). For the brightest
cluster galaxy (BCG) we adopt a position of (7$\farcs$114,4$\farcs$409), an
ellipticity of $e=0.30\pm 0.05$, a major axis PA of $\theta_{e}=152\pm 5$°
(Oguri, 2010) and a central velocity dispersion of $352\pm 13$ km s${}^{-}1$
(Inada et al., 2008).

Table 1: Redshift, Positions, Magnitudes, and Time Delays for the Observed
Images of the Lensed Sources

Name | $z_{s}$ |  $x$ (″) |  $y$ (″) | $\Delta m$ |  $\Delta t$ (days)
---|---|---|---|---|---
QSO.A |  | 0.000 | 0.000 |  $\equiv$0 | 825.99$\pm$ | 2.10 |
QSO.B |  |  $-$1.317 | 3.532 | 0.35 | 0.30 | 781.92$\pm$ | 2.20 |
QSO.C | 1.734 | 11.039 |  $-$4.492 | 0.87 | 0.30 |  $\equiv$0 |
QSO.D |  | 8.399 | 9.707 | 1.50 | 0.30 | 2456.99$\pm$ | 5.55 |
QSO.E |  | 7.197 | 4.603 | 6.3 | 0.8 |  $\cdots$ |
A1.1 |  | 3.93 |  $-$2.78 |  $\cdots$ |  $\cdots$ |
A1.2 |  | 1.33 | 19.37 |  $\cdots$ |  $\cdots$ |
A1.3 | 3.33 | 19.23 | 14.67 |  $\cdots$ |  $\cdots$ |
A1.4 |  | 18.83 | 15.87 |  $\cdots$ |  $\cdots$ |
A1.5 |  | 6.83 | 3.22 |  $\cdots$ |  $\cdots$ |
A2.1 |  | 4.13 |  $-$2.68 |  $\cdots$ |  $\cdots$ |
A2.2 |  | 1.93 | 19.87 |  $\cdots$ |  $\cdots$ |
A2.3 | 3.33 | 19.43 | 14.02 |  $\cdots$ |  $\cdots$ |
A2.4 |  | 18.33 | 15.72 |  $\cdots$ |  $\cdots$ |
A2.5 |  | 6.83 | 3.12 |  $\cdots$ |  $\cdots$ |
A3.1 |  | 4.33 |  $-$1.98 |  $\cdots$ |  $\cdots$ |
A3.2 |  | 2.73 | 20.37 |  $\cdots$ |  $\cdots$ |
A3.3 | 3.33 | 19.95 | 13.04 |  $\cdots$ |  $\cdots$ |
A3.4 |  | 18.03 | 15.87 |  $\cdots$ |  $\cdots$ |
A3.5 |  | 6.83 | 3.02 |  $\cdots$ |  $\cdots$ |
B1.1 |  | 8.88 |  $-$2.16 |  $\cdots$ |  $\cdots$ |
B1.2 | 2.74 |  $-$5.45 | 15.84 |  $\cdots$ |  $\cdots$ |
B1.3 |  | 8.33 | 2.57 |  $\cdots$ |  $\cdots$ |
B2.1 |  | 8.45 |  $-$2.26 |  $\cdots$ |  $\cdots$ |
B2.2 | 2.74 |  $-$5.07 | 16.04 |  $\cdots$ |  $\cdots$ |
B2.3 |  | 8.33 | 2.57 |  $\cdots$ |  $\cdots$ |
C1.1 |  | 10.25 |  $-$3.06 |  $\cdots$ |  $\cdots$ |
C1.2 | 3.28 |  $-$7.55 | 15.39 |  $\cdots$ |  $\cdots$ |
C1.3 |  | 8.49 | 2.72 |  $\cdots$ |  $\cdots$ |
C2.1 |  | 9.95 |  $-$3.36 |  $\cdots$ |  $\cdots$ |
C2.2 | 3.28 |  $-$7.30 | 15.44 |  $\cdots$ |  $\cdots$ |
C2.3 |  | 8.49 | 2.72 |  $\cdots$ |  $\cdots$ |

Note. — The position errors adopted for the QSO are 0$\farcs$04 whereas for
the galaxies are 0$\farcs$4.

Table 2: Positions, Ellipticities, Position Angles, and Luminosity Ratios of 14 Cluster Galaxy Members $x$ (″) |  $y$ (″) | e | PA (°) | $L/L_{BCG}$
---|---|---|---|---
30.78 | 4.50 | 0.2723 |  $-$2.90 | 0.5050 |  |
12.14 | 3.67 | 0.2426 |  $-$131.80 | 0.2970 |  |
2.76 | 14.13 | 0.1077 |  $-$161.50 | 0.2250 |  |
25.29 |  $-$9.06 | 0.1028 |  $-$70.00 | 0.0950 |  |
$-$9.22 |  $-$2.53 | 0.0885 |  $-$42.80 | 0.1510 |  |
14.54 | 24.23 | 0.0461 |  $-$82.00 | 0.1830 |  |
24.61 | 4.72 | 0.0433 |  $-$60.40 | 0.1720 |  |
9.36 | 2.41 | 0.2249 |  $-$100.30 | 0.4060 |  |
2.767 |  $-$0.171 | 0.0517 |  $-$21.8 | 0.238 |  |
14.789 |  $-$5.454 | 0.0453 |  $-$127.7 | 0.41 |  |
$-$1.359 | 0.482 | 0.0224 | 51.3 | 0.21 |  |
12.00 | 13.82 | 0.0711 |  $-$92.7 | 0.368 |  |
7.84 | 9.10 | 0.0118 |  $-$17.0 | 0.10 |  |
$-$7.21 |  $-$8.84 | 0.0661 |  $-$139.6 | 0.450 |  |
|  |  |  |  |  |

##  3 System Modeling

We modeled the mass distribution of the lens using both lensmodel (Keeton,
2001) and glafic (Oguri, 2010). Both packages employ parametrized mass
profiles and optimize the model using the downhill simplex method (see, e.g.,
Press et al., 1992). The final results are computed with glafic in order to
ease the comparison with Oguri (2010). The different mass components and any
constraints on their properties are described in the following subsections.

###  3.1 Dark Matter Halo

We modeled the dark matter (DM) mass distribution as a generalized Navarro-
Frenk-White (gNFW) profile (Navarro et al., 1997; Jing & Suto, 2000) with a 3D
density profile

$\rho=\frac{\rho_{s}}{(r/r_{s})^{\alpha}(1+r/r_{s})^{3-\alpha}}.$ (1)

With the addition of its 2D ellipticity, the projected profile is completely
described by seven parameters: $M_{vir}$ (virial mass), $x$ and $y$
(position), e (ellipticity), $\theta_{e}$ (position angle), $c_{-2}$
(concentration parameter), and $\alpha$ (inner slope). The relationships
between ($\rho_{r},r_{s}$) and ($M_{vir},c_{-2}$) can be found in glafic’s
manual111Available at
https://www.slac.stanford.edu/~oguri/glafic/index_v1.html. We do not impose
any restrictions on the parameters, however the inner slope $\alpha$ is set to
unity until the whole model is optimized.

###  3.2 BCG

The BCG is parametrized by a pseudo-Jaffe (pJaffe) profile

$\kappa=2\pi\left(\frac{\sigma}{c}\right)^{2}\frac{D_{ls}}{D_{os}}\left[\frac{1}{R}-\frac{1}{\sqrt{r_{trun}^{2}+R^{2}}}\right].$
(2)

This model is defined by six parameters when ellipticity is added: $\sigma$
(velocity dispersion), $x$ and $y$ (position), e (ellipticity), $\theta_{e}$
(position angle) and $r_{trun}$ (truncation radius, where beyond this scale,
the convergence falls as $R^{-3}$). We fix the position to the observed value
and include Gaussian priors on e and $\theta_{e}$. The observed velocity
dispersion does not directly correspond to the model dispersion and the
Appendix A details how we use the measured dispersion as a constraint to
obtain a model velocity dispersion prior of 325$\pm$20 km s${}^{-}1$. The
truncation radius is loosely constrained using the same prior as Oguri (2010)
($r_{trun}=8\arcsec\pm 4\arcsec$) because of the observed correlation between
the velocity dispersion and truncation radius (Natarajan et al., 2009).

###  3.3 Other Cluster Galaxy Members

The rest of the cluster galaxies are modeled by scaled pseudo-Jaffe ellipsoids
(mass profile ‘gals’ in the glafic software). Their $\sigma$ and $r_{trun}$
are scaled relative to the luminosity of the BCG as

$\sigma=\sigma_{*}\left(\frac{L}{L_{*}}\right)^{1/4}\;\textnormal{ and
}\;r_{trun}=r_{*}\left(\frac{L}{L_{*}}\right)^{1/2}$ (3)

where $\sigma_{*}$ and $r_{*}$ are model parameters. The ratio $L/L_{*}$,
ellipticity, and position angle of each galaxy member are constrained by the
observed values given in Table 2. Since the reference galaxy is the BCG, the
scale parameters $\sigma_{*}$ and $r_{*}$ are initially set equal to the
velocity dispersion and truncation radius of the BCG, but, during the
optimization, both are allowed to vary.

###  3.4 External Perturbations

We also include multipole perturbations to mimic deviations of the DM halo
from a perfect ellipsoid. Like Oguri (2010) we used multipoles of order m=2,
3, 4, and 5 (mpoles model) with the potential

$\phi=-\frac{\epsilon}{m}r^{2}\cos m(\theta-\theta_{\epsilon}-\pi/2).$ (4)

For m=2, the perturbation is an external shear (named ‘pert’ in glafic) so we
denote the m=2 parameters as $\gamma$ and $\theta_{\gamma}$ instead. We center
these potentials on the position of the DM halo.

###  3.5 Background Sources

At higher redshifts than the quasar there are three groups of galaxies (A, B,
and C) at different redshifts. Group A is composed of three galaxies each
lensed into five images. Groups B and C each have two galaxies and each is
split into three images. The observed positions of these images along with
their redshifts and magnitudes are also given in Table 1. All the background
sources are modeled as point sources with fixed redshifts so the fit
parameters are their unlensed positions.

In total we have 77 observational constraints and 37 model parameters leaving
us with $\nu$=40 degrees of freedom.

##  4 Results

The optimization was performed as follows. First we optimized all the
parameters fitting the images on the source plane. Once we found the minimum
$\chi^{2}$, we ran the optimization fitting the images on the image plane.
Next we computed the parameter uncertainties fitting the images on the source
plane in order to maintain a sensible computing time. Finally, we ran the
optimization on the image plane using as initial parameters the model that
produced the minimum $\chi^{2}$ in the uncertainties estimation step.

The central values and the 1$\sigma$ and 2$\sigma$ (corresponding to
$\Delta\chi^{2}=1$ and $\Delta\chi^{2}=4$, respectively) uncertainties are
given in Table 3. We define the uncertainties using the model most distant
from the $\chi^{2}$ minimum that is within the $\Delta\chi^{2}$ limit. The
$\chi^{2}$ contours and profiles for the mass $M_{vir}$, concentration
$c_{-2}$ and inner slope $\alpha$ of the gNFW model are displayed in Figure 1.
The parameters $M_{vir}$ and $c_{-2}$ are strongly correlated because they are
defined on scales much larger than the Einstein ring, an issue we discuss in
more detail below. The parameters do not have smooth $\chi^{2}$ profiles
because the minimizer has difficulties in finding the true minimum.

Table 3: Central Values and Errors at 1-$\sigma$ and 2-$\sigma$ for the Mass Model Parameters Model | Parameters |  | 1$\sigma$ | 2$\sigma$ | Best Fit
---|---|---|---|---|---
gNFW |  $M_{vir}$ ($10^{15}h^{-1}M_{\sun}$) | 1.0 | ${}^{+0.3}_{-0.2}$ | ${}^{+0.7}_{-0.4}$ | 1.0
$x$ (″) | 7.04 | ${}^{+0.07}_{-0.07}$ | ${}^{+0.13}_{-0.14}$ | 7.03
$y$ (″) | 4.95 | ${}^{+0.09}_{-0.10}$ | ${}^{+0.2}_{-0.19}$ | 4.88
e | 0.36 | ${}^{+0.02}_{-0.03}$ | ${}^{+0.04}_{-0.05}$ | 0.36
$\theta_{e}$ (°) | 158.6 | ${}^{+1.5}_{-1.8}$ | ${}^{+3}_{-3}$ | 158.1
$c_{-2}$ | 2.3 | ${}^{+0.8}_{-0.8}$ | ${}^{+1.7}_{-1.7}$ | 2.2
$\alpha$ | 1.18 | ${}^{+0.02}_{-0.03}$ | ${}^{+0.11}_{-0.18}$ | 1.20
pJaffe |  $\sigma$ (km s${}^{-}1$) | 296 | ${}^{+5}_{-4}$ | ${}^{+10}_{-12}$ | 290
e | 0.40 | ${}^{+0.06}_{-0.04}$ | ${}^{+0.12}_{-0.10}$ | 0.37
$\theta_{e}$ (°) | 152 | ${}^{+4}_{-3}$ | ${}^{+9}_{-7}$ | 154
$r_{trun}$ (″) | 10 | ${}^{+2}_{-2}$ | ${}^{+6}_{-4}$ | 10
gals |  $\sigma_{*}$ (km s${}^{-}1$) | 351 | ${}^{+12}_{-10}$ | ${}^{+20}_{-20}$ | 337
$r_{*}$ (″) | 9 | ${}^{+3}_{-2}$ | ${}^{+6}_{-4}$ | 9
pert |  $\gamma$ ($10^{-2}$) | 5.8 | ${}^{+1.0}_{-1.2}$ | ${}^{+2}_{-2}$ | 5.7
$\theta_{\gamma}$ (°) | 66 | ${}^{+3}_{-5}$ | ${}^{+7}_{-9}$ | 65
mpole3 |  $\epsilon$ ($10^{-2}$) | 1.6 | ${}^{+0.3}_{-0.3}$ | ${}^{+0.6}_{-0.6}$ | 1.4
$\theta_{\epsilon}$ (°) | -9 | ${}^{+2}_{-3}$ | ${}^{+6}_{-7}$ | -7
mpole4 |  $\epsilon$ ($10^{-2}$) | 1.1 | ${}^{+0.2}_{-0.2}$ | ${}^{+0.5}_{-0.4}$ | 1.1
$\theta_{\epsilon}$ (°) | 40 | ${}^{+3}_{-4}$ | ${}^{+8}_{-7}$ | 41
mpole5 |  $\epsilon$ ($10^{-2}$) | 1.39 | ${}^{+0.13}_{-0.2}$ | ${}^{+0.3}_{-0.3}$ | 1.25
$\theta_{\epsilon}$ (°) | 16.2 | ${}^{+2}_{-1.8}$ | ${}^{+6}_{-4}$ | 17.6
Figure 1: $\chi^{2}$ contours and profiles for the virial mass, concentration,
and inner slope parameters of the generalized NFW model (Eqn. 1). Red lines
and points mark the values for which the $\chi^{2}$ is minimized, orange lines
delimiter the 1$\sigma$ contours, and green lines the 2$\sigma$ contours. The
black (gray) shadow is the 1$\sigma$ (2$\sigma$) region.

The optimization on the image plane starting from the parameters which give
the minimum $\chi^{2}$ in the error computations is considered the best-fit
model and it is shown in the last column of Table 3. The central values of the
$\chi^{2}$ profiles do not exactly correspond to the best-fit parameters (see
columns 3 and 6 of Table 3) because the two results come from different
optimization procedures. Nonetheless, all the parameters of the best model fit
lie between the reported 1$\sigma$ confidence intervals except for $\sigma$
and $\sigma_{*}$ which lie between their 2$\sigma$ intervals.

The predicted image positions and critical curves of the four background
sources for this model are shown in Figure 2. This best-fit model has a
$\chi^{2}$ value of 59.91 in the image plane which corresponds to
$\chi^{2}_{\nu}=\chi^{2}/\nu=1.50$. The observational constraints on the
quasar contribute 11.76 to the total $\chi^{2}$ and the A, B, and C groups of
galaxies contribute 30.95, 2.92, and 8.16, respectively. Additionally, the
$\chi^{2}$ associated with the six parameter priors is 6.13. The model fits
all three delays well, with a $\chi^{2}$=1.33 for the delay constraints. The
model predicts a time delay between the leading image C and the fifth image E
of $\Delta t_{EC}$=2853.90 days or 7.81 yr, although this delay will be
virtually impossible to measure due to the faintness of image E. Additionally,
Table 4 gives the convergence $\kappa$ and shear $\gamma$ at the quasar
positions along with the total magnification.

Table 4: Convergence $\kappa$, Shear $\gamma$ and Magnification $\mu$ at the Quasar Image Positions Quasar Image | $\kappa$ | $\gamma$ | $\mu$
---|---|---|---
A | 0.728940 | 0.333242 |  $-$26.612171 |
B | 0.650724 | 0.232737 | 14.743285 |
C | 0.587686 | 0.225075 | 8.379131 |
D | 1.015191 | 0.489943 |  $-$4.169912 |
E | 5.669724 | 3.622670 | 0.115173 |

The goodness of the fit for the A group of galaxies is substantially worse
than for the rest of the background galaxies because the model predicts seven
images instead of the five observationally reported. This occurs because the
images A1.2, A2.2, and A3.2 lie very close to a critical curve. A very small
change in the source position can change the number of the images and that is
precisely what occurred in our model where one image has unfolded into three
(they are marked as ‘unfolded’ in Figure 2). The A group images are really
composed of extended arcs so the problem would be solved by increasing their
position uncertainties, although the model would be more realistic if we
modeled the galaxies as extended sources. Despite that, the reconstructed
image positions properly reproduce all the main structures of the system (see
Figure 3).

Figure 2: Critical curves and images positions for the four background
sources. The green dots mark the observed positions of the images. Circles
denote the quasar, A1, B1, and C1 images, squares represent the images of A2,
B2, and C2, and lastly the triangles mark the A3 images. Additionally, images
are color coded depending on their flux, from the brightest in red to the
dimmest in yellow.  Figure 3: The left panel shows the Hubble Space Telescope
Advanced Camera for Surveys (ACS)/WFC (GO-9744, PI: C. Kochanek) images of the
lens. The right panel shows the reconstructed mass distribution on a
logarithmic scale along with the images positions: quasar in white, and galaxy
group A in lime-green, B in blue, and C in turquoise. The orientation for both
figures is north up, west right.

##  5 Discussion and Conclusions

One of the main results from our model is the value of the inner slope of the
generalized NFW profile, $\alpha=1.18^{+0.02(+0.11)}_{-0.03(-0.18)}$ at the
68% (95%) confidence level. We must bear in mind that this value is the
asymptotic slope of the DM for $r\rightarrow 0$. However at this regime the
mass of the BCG becomes important and it is difficult to distinguish both
components because the lensing effect is sensitive to the total mass of the
system. A different parameterization of the BCG would alter the inner slope
value of the generalized NFW profile. If we wanted to determine the mass
distribution in the innermost regions of the cluster more precisely, we would
probably need a velocity dispersion profile in order to separate the two
components better.

Newman et al. (2013a, b) made such models for seven massive galaxy clusters
and found very shallow inner slopes of $0.50\pm
0.10\textnormal{(random)}^{+0.14}_{-0.13}\textnormal{(systematic)}$. On the
other hand, Smith et al. (2017) inferred a DM inner slope for the cluster
Abell 1201 of $\alpha=1.01\pm 0.12$ using the same procedure but with other
assumptions. Schaller et al. (2015) and He et al. (2020) tested this method
with the EAGLE hydrodynamical simulations and found that the results were
dependent on the assumptions about the mass-to-light ratio and the intrinsic
degeneracies of the profiles. The simulations suggest that the DM halo
distribution is compatible with the original NFW profile and our value. While
our model likely has degeneracies between the inner slope of the BCG (which we
hold fixed) and the DM halo, our result is compatible with these simulations
of massive galaxy clusters.

Our inferred virial mass and concentration parameter for the generalized NFW
cluster halo are in good agreement with those obtained with the lens model of
Oguri (2010). However they differ from the values reported by Ota et al.
(2006) based on X-ray observations, $M_{vir}=4.2^{+2.6}_{-1.5}\times
10^{14}\;h^{-1}M_{\sun}$ and $c_{-2}=6.1^{+1.5}_{-1.2}$. These differences can
be explained by the degeneracy between $M_{vir}$, $c_{-2}$ and $\alpha$ in
strong lensing models pointed out by Oguri (2010) and shown for our model in
Figure 1. The source of this degeneracy comes from the fact that the lensing
effect is mainly sensitive to the enclosed mass projected inside the region
delimited by the lensed images but the parameters $M_{vir}$ and $c_{-2}$ are
defined on the much larger scale of the virial radius and so are
extrapolations of our mass model. In this system, the outermost images are
$\sim$125 kpc from the center of the cluster whereas the parameters $M_{vir}$
and $c_{-2}$ are extrapolations to a radius $r_{vir}\approx$ 2 Mpc. It is
better to test the compatibility of the two results by comparing the mass
enclosed inside a smaller radius to which both methods are sensitive. Ota et
al. (2006) reported the projected mass inside 70 h-1kpc to be
$3.5^{+1.3}_{-0.8}\times 10^{13}\;h^{-1}M_{\sun}$, in good agreement with our
enclosed mass of $3.7\times 10^{13}\;h^{-1}M_{\sun}$.

We find an offset between the BCG and the DM halo of
$3.8^{+0.6(+1.4)}_{-0.7(-1.3)}$ kpc at the 68% (95%) confidence level.
Previous parametric models have found similar offsets for this system (Oguri
et al., 2004; Fohlmeister et al., 2007) although Oguri (2010) found an offset
compatible with zero. From fitting 10,000 galaxy clusters with strong lensing
measurements, Zitrin et al. (2012) found a mean separation of $18^{+37}_{-12}$
kpc where the uncertainty is the observed scatter. Kluge et al. (2020)
estimated the separation between the BCG and the intracluster light tracing
the DM distribution for local clusters, and found a mean offset of 36 kpc with
a sample scatter of 33 kpc. Our inferred offset is compatible with both
results.

We also find that the DM halo and the BCG are virtually aligned and with
similar ellipticities. On the other hand, the velocity dispersion scale
parameter for the cluster galaxy members differs from the BCG parameterization
even though they were started with values scaled to the BCG. This difference
suggests that BCGs are structurally different from the other cluster galaxy
members because they have more DM content than the rest. Another feature in
the BCG velocity dispersion is that the model fits the pseudo-Jaffe velocity
dispersion to a smaller value than the prior. That could be explained by the
fact that the observed velocity dispersion includes the effect of the DM halo
whereas the model velocity dispersion parameter represents only the mass of
the BCG.

In conclusion, we have made use of the recently measured time delay of image D
and previous observational data to model the lensing cluster SDSS J1004+4112.
The model parameters we obtained are broadly consistent with previous
estimates but they are now better constrained, in particular the inner slope
of the DM halo of the lensing cluster.

This research is based on observations made with the NASA/ESA Hubble Space
Telescope obtained from the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
under NASA contract NAS 5–26555. These observations are associated with
program GO-9744. J.A.M. and E.M. are supported by the Spanish Ministerio de
Ciencia e Innovación with the grants PID2020-118687GB-C32 and
PID2020-118687GB-C31. J.A.M. is also supported by the Generalitat Valenciana
with the project of excellence Prometeo/2020/085. C.S.K. is supported by NSF
grants AST-1908570 and AST-1814440.

##  Appendix A Model and Observed Velocity Dispersions

The glafic package parametrizes the pseudo-Jaffe profile with a velocity
dispersion which corresponds to the central velocity dispersion of a singular
isothermal sphere ($\sigma_{mod}$). On the other hand, the observed velocity
dispersion is the averaged velocity dispersion ($\sigma_{obs}$) along the line
of sight inside a finite size slit. The observed velocity dispersion is
produced by a combination of the BCG mass and the DM halo, however we assume
that the DM halo is negligible for a central velocity dispersion and we only
consider the dispersion velocity due to the BCG modeled as a pseudo-Jaffe
ellipsoid. In order to obtain their relationship we make use of the
expressions of the densities, masses, and velocities dispersion of Elíasdóttir
et al. (2007). Additionally, we transform the expressions into dimensionless
quantities plus a prefactor in order to simplify the expressions and we end
with the following definitions:

  * •

Surface density:

$\Sigma(R)=\frac{\sigma_{mod}^{2}}{2Gr_{trun}}\tilde{\Sigma}(x)\;\textnormal{
where
}\;\tilde{\Sigma}(x)=\frac{1}{\sqrt{f^{2}+x^{2}}}-\frac{1}{\sqrt{1+x^{2}}}$
(A1)

  * •

Volume density:

$\rho(r)=\frac{\sqrt{q}\left(1-f^{2}\right)\sigma_{mod}^{2}}{2\pi
Gr^{2}_{trun}}\tilde{\rho}(z)\;\textnormal{ where
}\;\tilde{\rho}(z)=\frac{1}{\left(f^{2}+z^{2}\right)\left(1+z^{2}\right)}$
(A2)

  * •

Enclosed projected mass:

$M_{2D}(R)=\frac{\pi
r_{trun}\sigma_{mod}^{2}}{Gq}\tilde{M}_{2D}(x)\;\textnormal{ where
}\;\tilde{M}_{2D}(x)=\sqrt{f^{2}+x^{2}}-\sqrt{1+x^{2}}+1-f$ (A3)

  * •

Enclosed mass:

$M_{3D}(r)=\frac{2r_{trun}\sigma_{mod}^{2}}{Gq}\tilde{M}_{3D}(z)\;\textnormal{
where }\;\tilde{M}_{3D}(z)=\arctan{\left(z\right)}-f\arctan{\left(z/f\right)}$
(A4)

  * •

Projected velocity dispersion:

$\sigma_{P}^{2}(R)=\frac{4\left(1-f^{2}\right)\sigma_{mod}^{2}}{\pi\sqrt{q}}\tilde{\sigma_{P}^{2}}(x)\;\textnormal{
where
}\;\tilde{\sigma_{P}^{2}}(x)=\frac{1}{\tilde{\Sigma}(x)}\int_{x}^{\infty}\frac{\tilde{M}_{3D}(z)\tilde{\rho}(z)}{z^{2}}\sqrt{z^{2}-x^{2}}dz$
(A5)

  * •

Averaged projected velocity dispersion inside R:

$\sigma_{obs}^{2}=\langle\sigma_{P}^{2}\rangle(R)=\frac{4\left(1-f^{2}\right)\sigma_{mod}^{2}}{\pi\sqrt{q}}\tilde{\langle\sigma_{P}^{2}\rangle}(x)\;\textnormal{
where
}\;\tilde{\langle\sigma_{P}^{2}\rangle}(x)=\frac{1}{\tilde{M}_{2D}(x)}\int_{0}^{x}\tilde{\sigma_{P}^{2}}(x^{\prime})\tilde{\Sigma}(x^{\prime})x^{\prime}dx^{\prime}$
(A6)

The dimensionless variables are $x\equiv R\sqrt{q}/r_{trun}$, $z\equiv
r\sqrt{q}/r_{trun}$, $f\equiv r_{core}/r_{trun}$, and $q=1-e$ is the axis
ratio. In our case, $f=0$ and the slit width employed to measure the velocity
dispersion was 0$\farcs$4 (Inada et al., 2008) so we used the average
dispersion inside R=0$\farcs$2\. Taking into account the central values and
errors for the ellipticity and the truncation radius, we obtain
$\sigma_{mod,prior}=325\pm 20$ km s${}^{-}1$ as the prior for the model
velocity dispersion.

## References

  * Chen et al. (2012) Chen, B., Dai, X., Kochanek, C. S., et al. 2012, ApJ, 755, 24
  * Elíasdóttir et al. (2007) Elíasdóttir, Á., Limousin, M., Richard, J., et al. 2007, arXiv:0710.5636
  * Fian et al. (2018) Fian, C., Guerras, E., Mediavilla, E., et al. 2018, ApJ, 859, 50
  * Fian et al. (2021) Fian, C., Mediavilla, E., Motta, V., et al. 2021, A&A, 653, A109
  * Fohlmeister et al. (2007) Fohlmeister, J., Kochanek, C. S., Falco, E. E., et al. 2007, ApJ, 662, 62
  * Fohlmeister et al. (2008) Fohlmeister, J., Kochanek, C. S., Falco, E. E., et al. 2008, ApJ, 676, 761
  * Gómez-Álvarez et al. (2006) Gómez-Álvarez, P., Mediavilla, E., Muñoz, J. A., et al. 2006, ApJ, 645, L5
  * He et al. (2020) He, Q., Li, H., Li, R., et al. 2020, MNRAS, 496, 4717
  * Inada et al. (2003) Inada, N., Oguri, M., Pindor, B., et al. 2003, Nature, 426, 810
  * Inada et al. (2005) Inada, N., Oguri, M., Keeton, C. R., et al. 2005, PASJ, 57, L7
  * Inada et al. (2006) Inada, N., Oguri, M., Morokuma, T., et al. 2006, ApJ, 653, L97
  * Inada et al. (2008) Inada, N., Oguri, M., Falco, E. E., et al. 2008, PASJ, 60, 27
  * Jackson (2011) Jackson, N. 2011, ApJ, 739, L28
  * Jing & Suto (2000) Jing, Y. P. & Suto, Y. 2000, ApJ, 529, L69
  * Kawano & Oguri (2006) Kawano, Y. & Oguri, M. 2006, PASJ, 58, 271
  * Keeton (2001) Keeton, C. R. 2001, astro-ph/0102340
  * Kluge et al. (2020) Kluge, M., Neureiter, B., Riffeser, A., et al. 2020, ApJS, 247, 43
  * Lamer et al. (2006) Lamer, G., Schwope, A., Wisotzki, L., et al. 2006, A&A, 454, 493
  * Liesenborgs et al. (2009) Liesenborgs, J., de Rijcke, S., Dejonghe, H., et al. 2009, MNRAS, 397, 341
  * McKean et al. (2021) McKean, J. P., Luichies, R., Drabent, A., et al. 2021, MNRAS, 505, L36
  * Mohammed et al. (2015) Mohammed, I., Saha, P., & Liesenborgs, J. 2015, PASJ, 67, 21
  * Motta et al. (2012) Motta, V., Mediavilla, E., Falco, E., et al. 2012, ApJ, 755, 82
  * Muñoz et al. (2022) Muñoz, J., Kochanek, C., Fohlmeister, J., et al. 2022 ApJ, in press, arXiv:2206.08597
  * Natarajan et al. (2009) Natarajan, P., Kneib, J.-P., Smail, I., et al. 2009, ApJ, 693, 970
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Newman et al. (2013a) Newman, A. B., Treu, T., Ellis, R. S., et al. 2013, ApJ, 765, 24
  * Newman et al. (2013b) Newman, A. B., Treu, T., Ellis, R. S., et al. 2013, ApJ, 765, 25
  * Oguri et al. (2004) Oguri, M., Inada, N., Keeton, C. R., et al. 2004, ApJ, 605, 78
  * Oguri (2010) Oguri, M. 2010, PASJ, 62, 1017
  * Ota et al. (2006) Ota, N., Inada, N., Oguri, M., et al. 2006, ApJ, 647, 215
  * Popović et al. (2020) Popović, L. Č., Afanasiev, V. L., Moiseev, A., et al. 2020, A&A, 634, A27
  * Press et al. (1992) Press, W. H., Teukolsky, S. A., Vetterling, W. T., et al. 1992, Cambridge: University Press, 2nd ed.
  * Richards et al. (2004) Richards, G. T., Keeton, C. R., Pindor, B., et al. 2004, ApJ, 610, 679
  * Ross et al. (2009) Ross, N. R., Assef, R. J., Kochanek, C. S., et al. 2009, ApJ, 702, 472
  * Saha et al. (2006) Saha, P., Read, J. I., & Williams, L. L. R. 2006, ApJ, 652, L5
  * Schaller et al. (2015) Schaller, M., Frenk, C. S., Bower, R. G., et al. 2015, MNRAS, 452, 343
  * Sharon et al. (2005) Sharon, K., Ofek, E. O., Smith, G. P., et al. 2005, ApJ, 629, L73
  * Smith et al. (2017) Smith, R. J., Lucey, J. R., & Edge, A. C. 2017, MNRAS, 471, 383
  * Tie & Kochanek (2018) Tie, S. S. & Kochanek, C. S. 2018, MNRAS, 473, 80
  * Williams & Saha (2004) Williams, L. L. R. & Saha, P. 2004, AJ, 128, 2631
  * Zitrin et al. (2012) Zitrin, A., Bartelmann, M., Umetsu, K., et al. 2012, MNRAS, 426, 2944


#  Wi-Fi Assist: Enhancing Vehicular Wi-Fi Connectivity with an
Infrastructure-driven Approach

Antonio Rodrigues  Carnegie Mellon UniversityPittsburghPAUSA Faculty of
Engineering, Univ. of PortoPorto, Portugal<EMAIL_ADDRESS>,  Peter
Steenkiste  Carnegie Mellon UniversityPittsburghPAUSA<EMAIL_ADDRESS>and  Ana
Aguiar  Faculty of Engineering, Univ. of PortoPortugal Instituto de
TelecomunicacoesPortoPortugal<EMAIL_ADDRESS>
###### Abstract.

Vehicles access the Internet via cellular networks, instead of Wi-Fi networks.
This choice has been mostly justified by the ubiquitous coverage of cellular
networks: Wi-Fi coverage has been shown to be inadequate in the past, even in
urban areas.

We argue that providing Internet connectivity to vehicles via Wi-Fi is worth a
revisit. Motivated by improvements in Wi-Fi network coverage in recent years,
we propose Wi-Fi Assist, an add-on to current Wi-Fi infrastructures which
differs from existing solutions in two key ways: (1) it is heavily
infrastructure-driven; and (2) defines an interface for low-latency
cooperation between different WLAN service sets, managed by different service
providers.

††copyright: none

##  1\. Introduction

In practice, vehicles access the Internet via cellular networks, instead of
Wi-Fi networks. Car manufacturers - e.g. BMW and Toyota - equip vehicles with
cellular radios to provide Internet access (bmw-assist, ; toyota-safety-
connect, ). Wi-Fi can deliver better throughput than cellular (3G
(Deshpande:2010:PCM:1879141.1879180, ), LTE (Deng:2014:WLB:2663716.2663727,
)), at lower energy and monetary costs, and reduce the increasing strain on
cellular network capacity via traffic offloading (6824752, ). Nonetheless, the
choice falls on cellular networks due to their ubiquitous coverage. Cellular
base stations cover large areas ($<1$ km to 10s of km) and their deployment is
carefully planned by mobile carriers. This keeps the number of handoffs low,
and allows for ‘make-before-break’ handoffs while under the service of the
same carrier. To a vehicle, this translates into continuous Internet
connectivity, with few interruptions. In contrast, the coverage of Wi-Fi APs
is narrower (10s of m), WLAN service sets are typically small - e.g. $<5$ APs
- and their deployment is often unplanned. As a result, in-car Internet
connectivity over Wi-Fi is highly susceptible to interruptions, even in urban
areas. Measurement studies in urban areas show that vehicular Wi-Fi sessions
may experience Internet connectivity disruptions up to 42% of non-overlapping
1 second segments, with a median length of 5 consecutive segments
(Deshpande:2010:PCM:1879141.1879180, ). Other studies report median periods
with no Internet connectivity up to $30$ seconds
(Bychkovsky:2006:MSV:1161089.1161097, ; Eriksson:2008:CVC:1409944.1409968, ).

In this short position paper, we argue that using Wi-Fi for vehicular Internet
is worth studying, for three reasons. First, the aforementioned measurement
studies were conducted in the late 2000s/early 2010s. Wi-Fi coverage in urban
areas is likely to have improved significantly by now. Second, these studies
only consider Internet connectivity via ‘open’ Wi-Fi hotspots
(Bychkovsky:2006:MSV:1161089.1161097, ; Eriksson:2008:CVC:1409944.1409968, )
or hotspots serviced by a particular provider
(Deshpande:2010:PCM:1879141.1879180, ). This limitation is likely to hide a
lot of connectivity opportunities. Finally, most research on Wi-Fi vehicular
connectivity focuses on client-centric approaches to reduce Wi-Fi connection
setup time, overlooking the benefits of infrastructure-side support. Solutions
such as IEEE 802.11k (4544755, ) and IEEE 802.11r (4573292, ) do include
infrastructure-side support, but are mostly tuned for ‘walking spped’
scenarios (Song:2017:WGT:3098822.3098846, ) and can only be applied within a
single WLAN service set. The infrastructure has a wider view over the whole
wireless system, both temporally and spatially (e.g., quality of Wi-Fi signal
overtime/at multiple locations). We believe that ISPs, which have some form of
control over most APs in urban scenarios, can enhance vehicular Internet
service by (1) preemptively configuring an appropriate set of APs in the
vehicles’ path; and (2) opening a set of cross-provider interfaces, allowing
configuration state to flow in-between APs managed by different administrative
domains. This can be enabled by the recently proposed concept of Network
Service Support (Panda:2016:ONI:2875951.2875953, ), which in turn relies on
Software Defined Networking (SDN), Network Function Virtualization (NFV) and
Cloudlets.

##  2\. Motivation

We obtain a preliminary profile of urban Wi-Fi coverage for vehicular
scenarios using data from the SenseMyCity project
(DBLP:journals/corr/RodriguesAB14, ). This dataset contains mobility traces
collected from smartphones of 100s of crowdsourcing volunteers in the city of
Porto, Portugal. More specifically, we use its Wi-Fi scan dataset, which
contains $\sim 57$ million entries of Wi-Fi scans, from a total 295.000
different APs. The data includes both pedestrian and vehicular traces: as a
conservative approach to isolate vehicular traces, we discard those with
median speeds lower than 5.6 m/s ($\sim 20$ km/h).

Figure 1. Map of Porto, Portugal, divided in $50\times 50$ meter cells.
Shading indicates number of different APs visible in the cell. Table 1. Data
summary

# of scans | 95300
---|---
# of visible APs | 80237
# of ESSIDs | 42843

# of APs per ESSID | 1-4 | 41988
---|---|---
5-9 | 684
10-49 | 150
50+ | 21
# of APs per auth. | Open | 21827
WPA2 | 48102
Other | 7599

Fig. 1 shows that Wi-Fi coverage is pervasive along the main roads: the
$50\times 50$ meter cells which overlap the thicker lines in the roadmap
report more than 20 different visible APs (more specifically, unique MAC
addresses). Moreover, additional data in Table 1 shows that most networks
(identified by ESSIDs) are composed by a small number of APs (1 to 4),
motivating the need for cooperation between networks managed by different WLAN
service sets. There are cases of ESSIDs which aggregate a large number of APs
- e.g. the ‘Fon’ network with $\sim 13$k different APs, or ‘MEO-WiFi’ with
$\sim 6$k APs111Data not shown in this document. \- however, we cannot derive
their degree of support for roaming optimizations from the data. Regarding
authentication, the most common type is ‘WPA2 Personal’. This is followed by
‘open’ APs, mostly part of large networks such as ‘Fon’ or ‘MEO-WiFi’, on
which users authenticate via a login/password combo. This shows there are
benefits in exploring ways of enabling APs protected by ‘WPA2-Personal’
authentication, which is effectively removing $>\sim 50$% of connectivity
opportunities for vehicles.

##  3\. Proposal : Wi-Fi Assist

We propose an add-on to current Wi-Fi infrastructures - Wi-Fi Assist\- which
aims at a general reduction of Wi-Fi connection setup time to any Wi-Fi
networks found by moving vehicles. If coupled with wide adoption by ISPs, this
could improve both Wi-Fi connectivity and throughput in vehicular scenarios,
previously shown to be at odds with each other due to DHCP delays in Wi-Fi
connection setup (Soroush:2011:CWM:2079296.2079300, ).

###  3.1. Overview

Despite its early stage, Wi-Fi Assist’s design differs from existing solutions
in two key ways: (1) it is heavily infrastructure-driven; and (2) defines an
interface for low-latency cooperation between different WLAN service sets,
managed by different service providers.

Figure 2. Overview of Wi-Fi Assist’s design

Infrastructure-driven design: Wi-Fi Assist aims at transferring work from the
client to the infrastructure side, so as to shrink non-connectivity periods as
much as possible. Non-connectivity periods have 3 main causes: (1) lack of Wi-
Fi signal coverage; (2) failure to connect to a Wi-Fi network; and (3) Wi-Fi
connection setup. While (1) is a ‘deal-breaker’, Fig 1 shows that it may be a
rare event. Wi-Fi Assist mitigates (2) and (3) by preemptively establishing
association/authentication state and encryption keys at Wi-Fi APs located
along the vehicle’s trajectory. Ideally, after a bootstrapping step, the
mobile client should only deal with transmitting and receiving data frames,
which are eventually handled by Wi-Fi APs within the vehicle’s range. To
accomplish this, Wi-Fi Assist servers keep track of a vehicle’s trajectory,
predict the possible future ramifications of the current trajectory and use
data aggregated from multiple participant APs to decide which are better
suited to handle a particular vehicle (steps 1, 2, 3 and 5 in Fig. 2).

Cross-provider interface: In order to expand the space of APs which can handle
a vehicle’s communication session, Wi-Fi Assist must be able to request the
establishment of association/authentication state in APs belonging to
different WLAN service sets, often managed by different providers (steps 4 and
5 in Fig. 2).

###  3.2. Challenges & Initial Design Options

Wi-Fi Assist’s realization involves the following challenges:

Latency minimization: Wi-Fi Assist adds new interfaces and system operations,
and thus new potential sources of latency, which must be minimized : (1)
communication with Wi-Fi Assist servers; (2) Wi-Fi Assist computations; and
(3) AP configuration operations. One option is to distribute and coordinate
Wi-Fi Assist servers across cloudlets (7807196, ;
Panda:2016:ONI:2875951.2875953, ) managed by the participating ISPs and
positioned near the configurable APs. The instantiation of Wi-Fi Assist
servers could be facilitated through NFV. Latency of configurations internal
to ISPs can be minimized by the pre-establishment of forwarding rules using
SDN capabilities.

Maintaining transport sessions: As seen in Fig. 2, a vehicle must maintain a
transport session while ‘hopping’ over different APs, included in different
WLAN service sets. An interesting option proposed by
(Croitoru:2015:TWM:2789770.2789786, ) is the use of Multipath TCP (MPTCP)
(RFC6824, ) for TCP transport, which allows for multiple sub-flows within a
the same TCP session, between different source IP:port pairs and a single
server.

Security: Wi-Fi Assist poses several security challenges along 3 axes: access
control, accountability and resource protection:

  * •

Access control: Credentials used by vehicles to access the scheduled Wi-Fi APs
along their path should not be abused by other entities. Furthermore, usage of
ISP configuration interfaces by Wi-Fi Assist should require strict
authentication policies, which must be fast to minimize latency. A token-based
authentication system (e.g. Kerberos) is an interesting solution.

  * •

Accountability: ISPs providing configuration interfaces must be able to
track/control traffic usage by Wi-Fi Assist instances, for billing purposes.

  * •

Resource protection: A malicious party - either a Wi-Fi Assist instance,
provider or vehicle - must not compromise the security of data traffic nor the
integrity of the network infrastructure.

## References

  * (1) BMW Assist.  www.bmwusa.com/bmwassist.  Last accessed: September, 2017.
  * (2) Toyota Safety Connect.  www.toyota.com/safety-connect.  Last accessed: September, 2017.
  * (3) IEEE Std 802.11k-2008 (Amendment 1: Radio Resource Measurement of Wireless LANs).  June 2008.
  * (4) IEEE Std 802.11r-2008 (Amendment 2: Fast Basic Service Set (BSS) Transition).  July 2008.
  * (5) J. G. Andrews, et al.  What Will 5G Be?  IEEE Journal on Selected Areas in Communications, 32(6):1065–1082, June 2014.
  * (6) V. Bychkovsky, et al.  A Measurement Study of Vehicular Internet Access Using in Situ Wi-Fi Networks.  In Proceedings of the 12th Annual International Conference on Mobile Computing and Networking, MobiCom ’06, pages 50–61, New York, NY, USA, 2006. ACM.
  * (7) A. Croitoru, et al.  Towards WiFi Mobility Without Fast Handover.  In Proceedings of the 12th USENIX Conference on Networked Systems Design and Implementation, NSDI’15, pages 219–234, Berkeley, CA, USA, 2015. USENIX Association.
  * (8) S. Deng, et al.  WiFi, LTE, or Both?: Measuring Multi-Homed Wireless Internet Performance.  In Proceedings of the 2014 Conference on Internet Measurement Conference, IMC ’14, pages 181–194, New York, NY, USA, 2014. ACM.
  * (9) P. Deshpande, et al.  Performance Comparison of 3G and Metro-scale WiFi for Vehicular Network Access.  In Proceedings of the 10th ACM SIGCOMM Conference on Internet Measurement, IMC ’10, pages 301–307, New York, NY, USA, 2010. ACM.
  * (10) J. Eriksson, et al.  Cabernet: Vehicular Content Delivery Using WiFi.  In Proceedings of the 14th ACM International Conference on Mobile Computing and Networking, MobiCom ’08, pages 199–210, New York, NY, USA, 2008. ACM.
  * (11) A. Ford, et al.  TCP Extensions for Multipath Operation with Multiple Addresses.  RFC 6824, RFC Editor, January 2013.  http://www.rfc-editor.org/rfc/rfc6824.txt.
  * (12) A. Panda, et al.  Open Network Interfaces for Carrier Networks.  SIGCOMM Comput. Commun. Rev., 46(1):5–11, jan 2016.
  * (13) J. G. P. Rodrigues, et al.  SenseMyCity: Crowdsourcing an Urban Sensor.  CoRR, abs/1412.2070, 2014.
  * (14) M. Satyanarayanan.  The Emergence of Edge Computing.  Computer, 50(1):30–39, Jan 2017.
  * (15) Z. Song, et al.  Wi-Fi Goes to Town: Rapid Picocell Switching for Wireless Transit Networks.  In Proceedings of the Conference of the ACM Special Interest Group on Data Communication, SIGCOMM ’17, pages 322–334, New York, NY, USA, 2017\. ACM.
  * (16) H. Soroush, et al.  Concurrent Wi-Fi for Mobile Users: Analysis and Measurements.  In Proceedings of the Seventh COnference on Emerging Networking EXperiments and Technologies, CoNEXT ’11, pages 4:1—-4:12, New York, NY, USA, 2011. ACM.