11institutetext: School of Mathematics and Statistics, Xi’an Jiaotong
University, Xi’an, China
11email<EMAIL_ADDRESS>
22institutetext: Tencent Jarvis Lab, Shenzhen, China

# Residual Moment Loss for
Medical Image Segmentation

Quanziang Wang 11 Renzhen Wang 11 Yuexiang Li 22 Kai Ma 22 Yefeng Zheng 22
Deyu Meng 11

###### Abstract

Location information is proven to benefit the deep learning models on
capturing the manifold structure of target objects, and accordingly boosts the
accuracy of medical image segmentation. However, most existing methods encode
the location information in an implicit way, _e.g.,_ the distance transform
maps, which describe the relative distance from each pixel to the contour
boundary, for the network to learn. These implicit approaches do not fully
exploit the position information (_i.e.,_ absolute location) of targets. In
this paper, we propose a novel loss function, namely residual moment (RM)
loss, to explicitly embed the location information of segmentation targets
during the training of deep learning networks. Particularly, motivated by
image moments, the segmentation prediction map and ground-truth map are
weighted by coordinate information. Then our RM loss encourages the networks
to maintain the consistency between the two weighted maps, which promotes the
segmentation networks to easily locate the targets and extract manifold-
structure-related features. We validate the proposed RM loss by conducting
extensive experiments on two publicly available datasets, _i.e.,_ 2D optic cup
and disk segmentation and 3D left atrial segmentation. The experimental
results demonstrate the effectiveness of our RM loss, which significantly
boosts the accuracy of segmentation networks.

###### Keywords:

Location information Image moment Medical image segmentation.

††footnotetext: Q. Wang and R. Wang contribute equally.

##  1 Introduction

Image segmentation plays an essential role in the field of medical image
analysis, such as disease diagnosis and surgery planning. Witnessing the
success of convolutional neural networks (CNNs) in computer vision, an
increasing number of researchers developed deep-learning-based frameworks for
medical image segmentation [15, 12]. However, most of the existing methods
simply brought the strategy adopted in computer vision to process medical
images, _i.e.,_ formulating the segmentation task as pixel-wise classification
(_e.g.,_ using the cross entropy/Dice loss to supervise the training of
segmentation networks). Such a formulation ignores the prior-knowledge—the
absolute position information of target objects provides useful information
for medical image segmentation, since the to-be-segmented targets often lie in
a more prominent low-dimensional manifold [17]. For example, the same organs
of different patients commonly have similar topology structures and positions.
Instead, the objects of the same category (_e.g.,_ dog) may appear in any
position of the natural images with varying sizes. Thus, we uncover a natural
question—how to fully exploit the location information to improve the accuracy
of medical image segmentation.

In the recent literature, there are studies embedding location information
into deep CNNs for performance improvement, which can be roughly separated
into two groups—explicit and implicit embedding. The former one explicitly
encoded the geometric position of pixels as features to enhance the capacity
of models. For example, Liu _et al._ [9] proposed CoordConv, which let
convolution access to its own input coordinates through concatenating extra
coordinate channels. Wang _et al._ [18] verified that CoordConv can enable a
spatially variant prediction in instance segmentation tasks. Hu _et al._ [4]
exploited spatial distances as geometry prior to promote the representation
learning of models. Similarly, to learn the local or global geometric
dependency, relative geometric position was incorporated into a self-attention
module in [1, 11]. For medical images, recent studies proposed to use the
encoded position information as an additional channel of network input to
assist CNNs better dealing with retinal nerve fiber layer defect detection [2]
and organ segmentation [13]. The latter implicit embedding methods usually
generated the distance transform maps (e.g., the distance of each pixel to the
closest boundary contour) from the ground truth, and designed corresponding
loss functions [10, 8, 7, 20] for optimization. These methods implicitly
encoded the relative position information rather than the absolute one.

All the above methods verified the benefit yielded by embedding the position
information into the deep CNNs to the segmentation performance. In fact, the
image moment, which simply embeds the object position with pixel values by
element-wise multiplication, is a potential direction to further boost the
segmentation accuracy of deep learning networks. Due to the excellent property
of affine invariant [5], image moments have been widely used as a feature
extractor in existing works [3, 6, 21]. Nevertheless, to our best knowledge,
no previous study tries to implement the image moment as a loss function for
neural networks to directly optimize. In this paper, we propose a novel image-
moment-based loss function, termed residual moment (RM) loss. Such a loss
function can explicitly embed the coordinate information into the training of
segmentation network by calculating the mean squared error of residual map
between the prediction map and ground-truth. We mathematically prove that the
RM loss is equivalent to a higher order center moment w.r.t. the square of
residual map. In essence, this indicates that our RM loss is a specific
weighted loss function, where the weights are location-aware and encode the
spatial relation of pixels. In a word, our RM loss can promote the
segmentation networks to easily locate the target objects and capture their
manifold structure.

In summary, our contributions are mainly three-fold. First, residual moment
loss is proposed to make segmentation network capture more position
information during training by incorporating coordinate information
explicitly. Second, we theoretically prove that the RM loss is equivalent to
the higher order center moment of the square of residual map implying that the
RM loss is essentially a weighted location-aware loss function. Third, the
loss function can be equipped with any existing segmentation networks without
changing their architecture and without extra inference burden during testing.
We execute comprehensive experiments with both 2D and 3D architectures to
verify its effectiveness.

##  2 Method

As aforementioned, our residual moment loss works as an auxiliary loss to help
the segmentation networks explicitly capture the location information during
training. In the following subsection, we first revisit image moment (Sec.
2.1), and then mathematically formulate the proposed residual moment loss
(Sec. 2.2). Finally, we introduce the way to implement the proposed residual
moment loss function in a deep learning network (Sec. 2.3). Without loss of
generality, we introduce our methods based on 2D images.

###  2.1 Revisiting of Image Moment

Many previous works regard the image moments as statistics to extract some
affinity invariant features for the subsequent processing. Mathematically, the
two-dimensional $(p+q)$th order moments in statistics are defined as

$m_{pq}(f)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}i^{p}j^{q}f(i,j)didj,\quad
p,q=0,1,2,\cdots,$ (1)

where $f(i,j)$ is a bivariate density distribution function of the random
variables $i$ and $j$. The center moments are reformulated by replacing $i$
and $j$ with $i-\bar{i}$ and $j-\bar{j}$ in Eq. 1, respectively, where
$\bar{i}=m_{10}/m_{00}$ and $\bar{j}=m_{01}/m_{00}$. Due to the sound
theoretical properties [5], center moments are widely-used in practice.

Since the random variables $i$ and $j$ are discrete for image processing, the
central moments of an $H\times W$ image are rewritten as:

$m_{pq}(f)=\sum_{i=1}^{H}\sum_{j=1}^{W}(i-\bar{i})^{p}(j-\bar{j})^{q}f(i,j),\quad
p,q=0,1,2,\cdots,$ (2)

where $i$ and $j$ represent the pixel position, $f(i,j)$ is the pixel value at
position $(i,j)$ of image $f$ and $(\bar{i},\bar{j})$ is the coordinate of
image’s center calculated via $\bar{i}=H/2,\quad\bar{j}=W/2$. For simplicity,
we marked the $\tilde{i}=i-\bar{i}$ and $\tilde{j}=j-\bar{j}$ in the
following. Referring to Eq. 2, image central moments embed the information of
object positions and pixel values by element-wise multiplication.

###  2.2 Residual Moment Loss

In this section, we present the proposed residual moment loss in details. For
notation simplicity, we first denote $\mu_{pq}(f)$ as:

$\mu_{pq}(f)=[\tilde{i}^{p}\tilde{j}^{q}f(i,j)]_{H\times W},\quad
p,q=0,1,2,\cdots.$ (3)

Here, $\mu_{pq}(f)$ is an $H\times W$ matrix with the same size as the input
image $f$. Formally, the proposed residual moment loss is defined by the mean
squared error (MSE) between $\mu_{pq}$ of the segmentation prediction
$\hat{y}$ and the ground-truth $y$, which can be formulated as:

$\textit{l}_{RM}^{(pq)}(\hat{y},y)=\sum_{i=1}^{H}\sum_{j=1}^{W}(\mu_{pq}(\hat{y})-\mu_{pq}(y))^{2}.$
(4)

Eq. 4 reveals that the segmentation prediction explicitly interacts with the
ground-truth not only in terms of object location but also the information of
pixel values. Such a flexible information flow enables the segmentation
networks to accurately locate the target objects and therefore fully exploit
useful information on the manifold structure. Beyond that, we mathematically
prove that the proposed loss function is equivalent to the higher order center
moment w.r.t. the residual map between the prediction map and the ground-truth
map.

###### Theorem 2.1

Denote the residual map between the prediction map and the ground-truth map as
$R=(\hat{y}-y)$. Then the loss function Eq. 4 is equivalent to the $(2p+2q)$th
order center moment of the element-wise square of residual map:

$\textit{l}_{RM}^{(pq)}(\hat{y},y)=m_{2p,2q}(R^{2}).$ (5)

Henceforth, we name $\textit{l}_{RM}^{(pq)}(\hat{y},y)$ as Residual Moment
Loss.

###### Proof

Taking the definition of $\mu_{pq}(f)$ into the Eq. 4, we have

$\displaystyle\textit{l}_{RM}^{(pq)}(\hat{y},y)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{H}\sum_{j=1}^{W}(\mu_{pq}(\hat{y})-\mu_{pq}(y))^{2}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{H}\sum_{j=1}^{W}\tilde{i}^{2p}\tilde{j}^{2q}(\hat{y}(i,j)-y(i,j))^{2}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{H}\sum_{j=1}^{W}\tilde{i}^{2p}\tilde{j}^{2q}R^{2}(i,j)=m_{2p,2q}(R^{2}).$

The proof is then completed. $\hfill\blacksquare$

From the proof, we can see that the proposed RM loss actually is the weighted
MSE loss, which equals the MSE loss as $p=q=0$. Note that the weights of our
RM loss are the coordinates of pixels, which implicitly model the spatial
relation between pixels and accordingly make our RM loss location-aware. This
is also the underlying reason why the proposed RM loss can help networks to
capture the position information naturally. In addition, from the Theorem 2.1,
the Hu invariant moments [5] can also be implemented with the proposed RM
loss.

As previously mentioned, our RM loss is an auxiliary loss function assisting
the segmentation network. Assuming the commonly used segmentation loss
function as $l_{Seg}(\hat{y},y)$, then, the full objective function can be
defined as:

$\textit{l}_{total}(\hat{y},y)=l_{Seg}(\hat{y},y)+\alpha
l_{RM}^{(pq)}(\hat{y},y),$ (6)

where $l_{Seg}(\hat{y},y)$ can be the widely-used cross entropy loss $l_{CE}$
or Dice loss $l_{Dice}$. In our following experiment, both losses, _i.e.,_
$l_{Seg}=l_{CE}+l_{Dice}$, are involved for a better segmentation performance.
The hyperparameter $\alpha>0$ is the weight of $l_{RM}^{(pq)}$. And the choice
of order $(p,q)$ is very flexible and we can further use a combination of
different orders (verified in Sec. 3.2).

Figure 1: The pipeline of the residual moment loss function. We use the
coordinate information to calculate the $\mu_{pq}$. The residual moment loss
is the mean squared error (MSE) between the $\mu_{pq}$ of predicted
segmentation and ground-truth. And the Segmentation Loss means the common used
loss functions, such as cross entropy and dice loss.

###  2.3 Coordinate Information

The whole pipeline of the residual moment loss is shown in Fig. 1. Referring
to the proof in the previous section, the proposed RM loss adopts the
coordinate information as the weights for pixels. However, it is wasteful and
unnecessary to train the neural network to automatically learn the coordinate
information. Hence, we construct the coordinate matrices ($C_{I}$ and $C_{J}$
with the size of $H\times W$) for the coordinate information embedding, which
can be represented as:

$C_{I}(i,j)=i,\quad C_{J}(i,j)=j,$ (7)

where $i=1,2,\dots,H$ and $j=1,2,\dots,W$. By Eq. 2, the center coordinate
matrices are $\tilde{C}_{I}=C_{I}-H/2$ and $\tilde{C}_{J}=C_{J}-W/2$. To
prevent numerical calculation from being too large, we normalize the
coordinate matrices through dividing $\tilde{C}_{I}$ and $\tilde{C}_{J}$ by H
and W, respectively. Therefore, the $\mu_{pq}$ can be reformulated as:

$\mu_{pq}(f)=\tilde{C}_{I}^{p}\odot\tilde{C}_{J}^{q}\odot f,\quad
p,q=0,1,2,\dots,$ (8)

where the $\odot$ is the element-wise production. Then we can calculate the RM
loss by Eq. 4 conveniently. The $\mu_{pq}$ and our RM loss can be extended to
3D by simply switching the coordinate matrices to 3D. The 3D $(p+q+r)$th RM
loss is also equivalent to the $(2p+2q+2r)$th center moment of the residual
volume.

##  3 Experiments

To verify the effectiveness and generalization of our method, we apply the
residual moment loss to 2D and 3D neural networks with various public
datasets. All our experiments are implemented with the PyTorch (1.3.0)
framework [14].

###  3.1 Implementation Details

We evaluate our method on two datasets, which are the 2D optic cup and disk
segmentation dataset: Drishti-GS (GS) [16] and the left atrial (LA) 3D
gadolinium-enhanced magnetic resonance imaging
(MRI)††http://atriaseg2018.cardiacatlas.org/ [19].

For the GS dataset, we use 50 images for training and 50 for testing. These
images are resized to $512\times 512$ for computational efficiency. We apply
the 2D U-Net [15] as the baseline and set the weight of RM loss in Eq. 6 as
$\alpha=1$. The model is trained by the stochastic gradient descent (SGD) with
the learning rate of 0.001 for 2000 iterations.

The LA dataset contains 100 training cases and 54 testing cases. Following the
experimental setting in [10], we randomly selected 16 cases for training and
20 cases for testing for a fair comparison with [10]. We cropped all cases
centering at the heart region to alleviate the problem of unbalanced
categories. All cases are normalized by subtracting the mean and dividing by
the standard deviation. We use 3D V-Net [12] as the baseline which is
optimized by SGD with the 0.01 learning rate. To prevent overfitting, we use
dropout during training but turn it off in the inference stage. Besides, the
RM loss weight is $\alpha=0.01$ in Eq. 6.

The evaluation metrics employed in our experiments include two region-based
metrics (_i.e._ , Dice and Jaccard), and two boundary-based metrics (_i.e._ ,
Average Surface Distance (ASD) and 95% Hausdorff Distance (95HD)). The region-
based metrics and boundary-based metrics can reflect the performance of
segmentation and shape similarity, respectively.

###  3.2 Drishti-GS Segmentation Results

Figure 2: The optic cup and disk results of the GS dataset. The green and blue
lines represent the contours of ground-truth and the segmentation results,
respectively. The first row is the results of the baseline and the second row
is our method, while the first two columns are the cup results and the last
two columns are the disk results.

The quantitative results are listed in Table 1. We can see that: 1) Our
residual moment loss with moments of different orders consistently outperforms
the baseline model. 2) Compared to the MSE loss (_i.e._ , the zeroth order RM
loss $l_{RM}^{(0,0)}$ with no explicit position information introduced), the
proposed RM loss achieves a superior result. 3) The RM loss with a combination
of three orders (_i.e._ , (2, 0), (0, 2), and (2, 2)) performs the best on all
the metrics. This can be rationally explained as the coordinate information of
x and y axes is asymmetric, which leads to the superposition of multiple RM
losses of different orders that can introduce various coordinate information.
From Fig. 2, it can be seen intuitively that our method leads to lower false
positive results than the baseline, which is mainly due to its strong capacity
of capturing the location information of the cup and disk.

Table 1: Optical cup and disk segmentation accuracy with mean (standard
deviation) on the GS dataset.

Dataset | Method |  Dice (%) $\uparrow$ |  Jaccard (%) $\uparrow$ |  95HD $\downarrow$ |  ASD $\downarrow$
---|---|---|---|---|---
GS cup | U-Net |  80.9 (3.50) |  69.1 (4.66) |  13.9 (1.93) |  6.80 (1.16)
MSE Loss ($l_{RM}^{(0,0)}$) |  81.3 (2.67) |  69.5 (3.46) |  15.9 (2.40) |  7.42 (0.78)
| RM Loss $l_{RM}^{(1,0)}$ |  81.3 (2.17) |  69.7 (2.75) |  14.5 (1.69) |  6.81 (0.50)
| RM Loss $l_{RM}^{(1,0)}+l_{RM}^{(0,1)}$ |  82.3 (3.94) |  71.0 (5.41) |  15.2 (2.84) |  7.09 (0.93)
| RM Loss $l_{RM}^{(2,0)}$ |  83.1 (4.07) |  72.5 (5.46) |  14.4 (2.48) |  6.68 (0.78)
| RM Loss $l_{RM}^{(2,0)}+l_{RM}^{(0,2)}$ |  84.0 (2.74) |  73.5 (3.72) |  14.1 (4.12) |  6.72 (1.63)
| RM Loss $l_{RM}^{(2,0)}+l_{RM}^{(0,2)}+l_{RM}^{(2,2)}$ |  84.4 (1.84) |  73.7 (1.35) |  13.2 (3.28) |  6.46 (1.56)
GS disk | U-Net |  92.2 (2.18) |  86.6 (3.28) |  19.7 (15.3) |  7.25 (5.32)
MSE Loss ($l_{RM}^{(0,0)}$) |  91.6 (3.46) |  85.9 (4.96) |  20.2 (16.5) |  7.02 (6.12)
| RM Loss $l_{RM}^{(1,0)}$ |  93.0 (2.18) |  88.0 (3.25) |  18.3 (18.01) |  7.14 (6.81)
| RM Loss $l_{RM}^{(1,0)}+l_{RM}^{(0,1)}$ |  93.4 (1.29) |  88.5 (2.10) |  18.0 (16.1) |  7.42 (6.94)
| RM Loss $l_{RM}^{(2,0)}$ |  93.3 (1.55) |  88.5 (2.04) |  13.0 (2.86) |  4.89 (1.48)
| RM Loss $l_{RM}^{(2,0)}+l_{RM}^{(0,2)}$ |  93.0 (1.93) |  87.8 (2.90) |  12.3 (2.67) |  4.60 (1.37)
| RM Loss $l_{RM}^{(2,0)}+l_{RM}^{(0,2)}+l_{RM}^{(2,2)}$ |  94.5 (0.64) |  90.1 (0.93) |  10.7 (2.24) |  3.80 (0.73)
|  |  |  |  |

###  3.3 Left Atrial MRI Segmentation Results

Our method can be easily extended to 3D. Table 2 and Fig. 3 show the
quantitative and qualitative results of the LA dataset. Here, we compare with
some implicit position embedding methods, which are boundary loss [8],
Hausdorff distance loss [7] and signed distance function loss [20]. The
experimental results of these three comparison methods are directly taken from
Ma _et al._ [10] and our experimental settings followed their work.

As listed in Table 2, it can be easily observed that: 1) For Dice and Jaccard,
our residual moment loss achieves the best results among all the comparison
methods; 2) Our method obtains an obvious improvement on 95HD and ASD compared
with the baseline, but fails to outperform implicit position encoding methods
for the two boundary-based metrics. This is probably because the loss
functions of these methods are specifically designed to minimize the distance-
based metrics. As shown in Fig. 3, it is seen that the residual moment loss
can significantly improve the ability to extract location information to
reduce the misclassification around the objects, which verifies that the
proposed loss function is effective in extracting location information.

Figure 3: Visualization of the left atrial segmentation results on the LA MRI
dataset. The first column is ground-truth. The second and third columns are
the results of baseline and our method, respectively. Table 2: Left atrial
segmentation accuracy with mean (standard deviation) on the LA MRI dataset.

Method |  Dice (%) $\uparrow$ |  Jaccard (%) $\uparrow$ |  95HD $\downarrow$ |  ASD $\downarrow$
---|---|---|---|---
V-Net |  85.4 (0.69) |  75.2 (1.00) |  26.0 (6.28) |  7.44 (1.81)
Boundary Loss [8] |  85.0 (5.64) |  74.2 (7.87) |  20.8 (15.0) |  5.43 (3.43)
Hausdorff Distance Loss [7] |  85.5 (4.96) |  75.0 (7.30) |  15.9 (13.3) |  4.46 (3.68)
Signed Distance Function Loss [20] |  84.2 (8.48) |  73.5 (11.0) |  13.5 (11.2) |  3.24 (3.10)
RM Loss $l_{RM}^{(2,0,0)}+l_{RM}^{(0,2,0)}+l_{RM}^{(0,0,2)}$ |  85.9 (0.76) |  75.7 (1.01) |  26.0 (5.67) |  7.13 (1.50)
RM Loss $l_{RM}^{(2,0,0)}+l_{RM}^{(0,2,0)}+l_{RM}^{(0,0,2)}+l_{RM}^{(2,2,2)}$ |  86.5 (0.51) |  76.6 (0.75) |  21.4 (3.30) |  6.16 (1.01)
|  |  |  |

##  4 Conclusion

In this work, we proposed a novel residual moment loss function to extract
location information in medical image segmentation. Motivated by image
moments, we explicitly encoded the coordinate information of pixels (or
voxels) to the RM loss, which is simple but can capture the target location
effectively. In addition, our method is also easy to optimize with high
computational efficiency. The experimental results demonstrated that our
method could be adapted to various data types and network architectures. The
method can also be easily embedded into other network training strategies and
used in different practical problems, which will be investigated in our future
research. Source code will be publicly available.

## References

  * [1] Bello, I., Zoph, B., Vaswani, A., Shlens, J., Le, Q.V.: Attention augmented convolutional networks. In: Proceedings of the IEEE/CVF International Conference on Computer Vision. pp. 3286–3295 (2019)
  * [2] Ding, F., Yang, G., Ding, D., Cheng, G.: Retinal nerve fiber layer defect detection with position guidance. In: International Conference on Medical Image Computing and Computer Assisted Intervention. pp. 745–754. Springer (2020)
  * [3] Flusser, J., Suk, T.: Rotation moment invariants for recognition of symmetric objects. IEEE Transactions on Image Processing 15(12), 3784–3790 (2006)
  * [4] Hu, H., Zhang, Z., Xie, Z., Lin, S.: Local relation networks for image recognition. In: Proceedings of the IEEE/CVF International Conference on Computer Vision. pp. 3464–3473 (2019)
  * [5] Hu, M.K.: Visual pattern recognition by moment invariants. IRE Transactions on Information Theory 8(2), 179–187 (1962)
  * [6] Iscan, Z., Yüksel, A., Dokur, Z., Korürek, M., Ölmez, T.: Medical image segmentation with transform and moment based features and incremental supervised neural network. Digital Signal Processing 19(5), 890–901 (2009)
  * [7] Karimi, D., Salcudean, S.E.: Reducing the Hausdorff distance in medical image segmentation with convolutional neural networks. IEEE Transactions on Medical Imaging 39(2), 499–513 (2019)
  * [8] Kervadec, H., Bouchtiba, J., Desrosiers, C., Granger, E., Dolz, J., Ayed, I.B.: Boundary loss for highly unbalanced segmentation. Medical Image Analysis 67, 101851 (2021)
  * [9] Liu, R., Lehman, J., Molino, P., Such, F.P., Frank, E., Sergeev, A., Yosinski, J.: An intriguing failing of convolutional neural networks and the coordconv solution. In: Proceedings of the 32nd International Conference on Neural Information Processing Systems. pp. 9628–9639 (2018)
  * [10] Ma, J., Wei, Z., Zhang, Y., Wang, Y., Lv, R., Zhu, C., Chen, G., Liu, J., Peng, C., Wang, L., et al.: How distance transform maps boost segmentation CNNs: an empirical study. In: Medical Imaging with Deep Learning. pp. 479–492. PMLR (2020)
  * [11] Ma, X., Fu, S.: Position-aware recalibration module: Learning from feature semantics and feature position. In: Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intelligence (2020)
  * [12] Milletari, F., Navab, N., Ahmadi, S.A.: V-Net: Fully convolutional neural networks for volumetric medical image segmentation. In: Fourth International Conference on 3D Vision. pp. 565–571. IEEE (2016)
  * [13] Murase, R., Suganuma, M., Okatani, T.: How can CNNs use image position for segmentation? arXiv preprint arXiv:2005.03463 (2020)
  * [14] Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., Lerer, A.: Automatic differentiation in pytorch (2017)
  * [15] Ronneberger, O., Fischer, P., Brox, T.: U-Net: Convolutional networks for biomedical image segmentation. In: International Conference on Medical Image Computing and Computer Assisted Intervention. pp. 234–241. Springer (2015)
  * [16] Sivaswamy, J., Krishnadas, S., Joshi, G.D., Jain, M., Tabish, A.U.S.: Drishti-GS: Retinal image dataset for optic nerve head (ONH) segmentation. In: IEEE 11th International Symposium on Biomedical Imaging. pp. 53–56. IEEE (2014)
  * [17] Wang, R., Cao, S., Ma, K., Meng, D., Zheng, Y.: Pairwise semantic segmentation via conjugate fully convolutional network. In: International Conference on Medical Image Computing and Computer Assisted Intervention. pp. 157–165. Springer (2019)
  * [18] Wang, X., Kong, T., Shen, C., Jiang, Y., Li, L.: SOLO: Segmenting objects by locations. In: European Conference on Computer Vision. pp. 649–665. Springer (2020)
  * [19] Xiong, Z., Xia, Q., Hu, Z., Huang, N., Bian, C., Zheng, Y., Vesal, S., Ravikumar, N., Maier, A., Yang, X., et al.: A global benchmark of algorithms for segmenting the left atrium from late gadolinium-enhanced cardiac magnetic resonance imaging. Medical Image Analysis 67, 101832 (2021)
  * [20] Xue, Y., Tang, H., Qiao, Z., Gong, G., Yin, Y., Qian, Z., Huang, C., Fan, W., Huang, X.: Shape-aware organ segmentation by predicting signed distance maps. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 34, pp. 12565–12572 (2020)
  * [21] Zhang, Y., Wang, S., Sun, P., Phillips, P.: Pathological brain detection based on wavelet entropy and hu moment invariants. Bio-Medical Materials and Engineering 26(s1), S1283–S1290 (2015)


# Using Bayesian Evidence Synthesis Methods to Incorporate Real World Evidence
in Surrogate Endpoint Evaluation

Lorna Wheaton  Biostatistics Research Group, Department of Health Sciences,
University of Leicester, University Road, Leicester LE1 7RH, UK  Anastasios
Papanikos  Biostatistics Research Group, Department of Health Sciences,
University of Leicester, University Road, Leicester LE1 7RH, UK
GlaxoSmithKline R&D Centre, GlaxoSmithKline, Stevenage, UK  Anne Thomas
Leicester Cancer Research Centre, University of Leicester, Leicester, UK
Sylwia Bujkiewicz

## Abstract

Objective: Traditionally validation of surrogate endpoints has been carried
out using RCT data. However, RCT data may be too limited to validate surrogate
endpoints. In this paper, we sought to improve validation of surrogate
endpoints with the inclusion of real world evidence (RWE). Study Design and
Setting: We use data from comparative RWE (cRWE) and single arm RWE (sRWE), to
supplement RCT evidence for evaluation of progression free survival (PFS) as a
surrogate endpoint to overall survival (OS) in metastatic colorectal cancer
(mCRC). Treatment effect estimates from RCTs, cRWE and matched sRWE, comparing
anti-angiogenic treatments with chemotherapy, were used to inform surrogacy
patterns and predictions of the treatment effect on OS from the treatment
effect on PFS.

Results: Seven RCTs, four cRWE studies and three matched sRWE studies were
identified. The addition of RWE to RCTs reduced the uncertainty around the
estimates of the parameters for the surrogate relationship. Addition of RWE to
RCTs also improved the accuracy and precision of predictions of the treatment
effect on OS obtained using data on the observed effect on PFS.

Conclusion: The addition of RWE to RCT data improved the precision of the
parameters describing the surrogate relationship between treatment effects on
PFS and OS and the predicted clinical benefit.

Keywords— Surrogate Endpoint, Real world evidence, Bayesian meta-analysis,
Colorectal cancer

##  1 Introduction

Surrogate endpoints are often used when it takes too long, is too expensive or
too difficult to observe the effect on the final clinical outcome of interest
[1]. However, before surrogate endpoints can be used, for example, for
regulatory approvals, they should be validated [2, 3]. Surrogate endpoints can
be validated based on three levels of association; (1) biological
plausibility, (2) individual-level surrogacy and (3) trial-level surrogacy
[4]. However, identifying and validating potential surrogate endpoints can be
difficult when data for such analysis are limited. Traditionally, surrogate
endpoint evaluation has been carried out using only data from randomised
controlled trials (RCTs). Shortages of RCT data are becoming more common as
precision medicine evolves and treatments become more effective and are
targeted to specific patient populations, leading to smaller cohorts of
patients where it takes longer to observe the treatment effect on the final
outcome with reasonable precision. This is due to fewer events (such as
deaths) recorded in patients receiving targeted therapies and thus high
uncertainty around the effectiveness estimates and, as a consequence, around
the estimates of association between the treatment effects on surrogate
endpoint and final outcome. It is therefore possible that a putative surrogate
endpoint cannot be validated and treatments may not be granted conditional
approval based on treatment effects on the questionable surrogate endpoint.
However, in recent years there has been increased interest in the use of real
world evidence (RWE) at all stages of drug development [5, 6, 7, 8, 9, 10, 11,
12]. This is because RWE can increase the evidence base for decision-making,
increase follow-up times and is often more generalisable to the target
population. Addition of RWE could improve validation of surrogate endpoints
which could not be validated using RCT data alone.

In this paper, we explored how RWE can be used to strengthen the evidence base
for surrogate endpoint evaluation. We made use of comparative RWE (cRWE) and
single arm RWE (sRWE) to supplement RCT data on the effectiveness of anti-
angiogenic therapies in metastatic colorectal cancer (mCRC). We then
investigated the impact of the addition of RWE on the estimates of the
surrogate relationship between treatment effects on PFS and OS.

The remainder of this paper is structured as follows. Data sources and the
statistical methods are described in Section 2. The results are presented in
Section 3 which is followed by discussion and conclusions in Sections 4 and 5.

##  2 Methods

###  2.1 Data Sources

####  2.1.1 Randomised Controlled Trials

Data were obtained from a prior literature review conducted by Ciani et al
[13], which included treatment effect estimates from 11 RCTs in mCRC that
assessed anti-angiogenic treatments such as Bevacizumab, combined with various
chemotherapies, such as FOLFOX. In the review, Ciani et al [13], defined OS as
time from randomisation to time of death and PFS as time from randomisation to
tumour progression or death from any cause.

For this paper, we extracted the treatment effects (logHRs) on PFS and OS.
RCTs were only included in the analysis where the control arm was chemotherapy
and the treatment arm was an anti-angiogenic treatment plus chemotherapy.

####  2.1.2 Comparative Real World Evidence

We carried out a literature review to identify cRWE evaluating anti-angiogenic
treatments for mCRC. The following combinations of terms were used to search
for studies published between 2000 and 2020 in the PubMed database: (1)
“metastatic colorectal cancer”, (2) “cohort”, “cohort study”, “retrospective”
or “prospective”, (3) “PFS”, (4) “OS” and (5) “antiangiogenic” or
“bevacizumab”.

Abstracts, titles and, where necessary, full articles were screened and
studies which were not relevant were removed. LogHRs on PFS and OS and their
corresponding standard errors were extracted from the remaining studies. To
account for potential bias, cRWE studies were only included if they reported
treatment effects adjusted for baseline characteristics or potential
confounders.

####  2.1.3 Single Arm Real World Evidence

Papanikos et al [14] identified 16 single-arm observational studies evaluating
anti-angiogenic treatments alone. To obtain relative treatment effects
required for surrogate endpoint evaluation, we carried out a literature review
of single-arm studies of chemotherapy alone, which subsequently were used as
control arms. The following combinations of terms were used to search for
single-arm studies on chemotherapy published between 2000 and 2020 in the
PubMed database: (1) “metastatic colorectal cancer”, (2) “chemotherapy”, (3)
“cohort”, “cohort study”, “retrospective” or “prospective”, (4) “progression”
or “PFS” and (5) “overall survival” or “OS”.

The following terms could not be contained in the title or abstract of the
studies; “antiangiogenic”, “bevacizumab”, “cetuximab”, “aflibercept”,
“randomised trial”, “randomized trial”, or “phase”. These terms were excluded
to prevent RCTs and cRWE being returned. Any additional studies found outside
the database search were also included.

###  2.2 Matching Single Arm Studies

Unlike RCTs and cRWE, treatment effects cannot be extracted from sRWE as
single-arm studies do not make comparisons. To obtain relative treatment
effects, sRWE studies were matched using aggregate level data according to the
method proposed by Schmitz et al [15]. The distance $\Delta_{tot}$ between any
two single-arm studies $j$ and $k$ was determined as the weighted average of
differences in covariates.

$\Delta_{tot}[j,k]=\frac{\sum_{i=1}^{n}w_{i}\cdot\Delta_{i}[j,k]}{\sum_{i=1}^{n}w_{i}}$
(1)

where $n$ is the number of covariates, $w_{i}$ refers to the weights given to
the covariates and $\Delta_{i}[j,k]$ is the normalised difference between
studies $j$ and $k$ in covariate $i$. This distance takes a value between $0$
and $1$, where smaller values indicate more similar studies. Distance measures
between treatment arms for RCTs and cRWE were also calculated. Since there is
no consensus on a threshold for similarity, the maximum distance metric
between arms of RCTs was used as the threshold. Where multiple matches were
possible, matches with the smallest distance measure were used. The weight of
each of the covariates was decided based on rankings from a consensus
statement [16].

###  2.3 Obtaining Treatment Effects for Matched Single Arm Studies

WebPlotDigitizer was used to extract data from Kaplan-Meier curves from each
arm of the matched sRWE studies. Kaplan-Meier curves from RCTs and cRWE were
also digitized to compare digitized and reported logHRs. Data from risk
tables, reporting the number of patients at risk in each arm at regular time
intervals, was also extracted to improve approximated IPD [17]. Where risk
tables were not reported, the number of patients and total number of events in
each arm were used.

Data extracted from Kaplan-Meier curves and risk tables (where available) were
used in Stata to reconstruct IPD using the ipdfc command by Wei and Royston
[18]. The Cox proportional hazards model with a single covariate for treatment
arm was used to analyse the reconstructed IPD to obtain logHRs on PFS and OS
for matched sRWE studies, cRWE and RCTs.

###  2.4 Statistical Analysis

The standard model for surrogate endpoint evaluation by Daniels and Hughes
[19], denoted here as D&H, and bivariate random-effects meta-analysis (BRMA)
using the product normal formulation (PNF) were used as alternative methods to
model the correlated treatment effects (logHRs) on the surrogate endpoint
(PFS) and final outcome (OS) using a Bayesian framework. The models were
applied to RCT data alone, RCTs and cRWE and RCTs, cRWE and matched sRWE.
Sensitivity analyses to vague prior distributions were conducted for both
models.

####  2.4.1 Daniels and Hughes Model

Daniels and Hughes proposed that the observed treatment effects measured on
the surrogate endpoint ($Y_{1i}$) and final outcome ($Y_{2i}$) come from a
bivariate normal distribution and estimate the underlying true effects on the
surrogate and final outcomes ($\delta_{1i}$ and $\delta_{2i}$) from each study
$i$ with corresponding within-study standard deviations ($\sigma_{1i}$ and
$\sigma_{2i}$) and within-study correlation ($\rho_{wi}$);

$\left(\begin{array}[]{c}Y_{1i}\\\ Y_{2i}\end{array}\right)\sim
N\left(\left(\begin{array}[]{l}\delta_{1i}\\\
\delta_{2i}\end{array}\right),\left(\begin{array}[]{cc}\sigma_{1i}^{2}&\sigma_{1i}\sigma_{2i}\rho_{wi}\\\
\sigma_{1i}\sigma_{2i}\rho_{wi}&\sigma_{2i}^{2}\end{array}\right)\right)$ (2)

The true effects measured on the surrogate endpoint ($\delta_{1i}$) are
assumed to be independent in each study (fixed effects). It is also assumed
there is a linear relationship between the true treatment effects on the final
outcome and the surrogate endpoint.

$\delta_{2i}\mid\delta_{1i}\sim
N\left(\lambda_{0}+\lambda_{1}\delta_{1i},\psi_{2}^{2}\right)$ (3)

Daniels and Hughes referred to a surrogate relationship as perfect when the
following conditions were met: (a) $\lambda_{0}=0$, (b) $\lambda_{1}\neq 0$
and (c) $\psi_{2}^{2}=0$. These conditions state that (a) no treatment effect
on the surrogate endpoint implies no treatment effect on the final outcome,
(b) the slope is not zero implying an association between treatment effects on
the surrogate and final outcomes and (c) treatment effects on the final
outcome can be perfectly predicted by treatment effects on the surrogate
endpoint.

To implement this model in a Bayesian framework, non-informative prior
distributions were placed on the fixed effects
$\mathrm{\delta}_{1i}\sim\mathrm{N}\left(0,10^{4}\right)$ and regression
parameters $\mathrm{\lambda}_{0,1}\sim\mathrm{N}\left(0,10^{4}\right)$. To
ensure a non-informative prior distribution on the conditional variance, a
uniform prior distribution was placed on the conditional standard deviation
$\mathrm{\psi_{2}}\sim\mathrm{Unif}\left(0,2\right)$. A minimally informative
prior distribution $\mathrm{\rho_{wi}}\sim\mathrm{Unif}\left(0,1\right)$ was
placed on the within-study correlation.

####  2.4.2 Bivariate Random-Effects Meta-Analysis (Product Normal
Formulation)

The D&H model does not estimate correlation or the study-level $R^{2}$ which
are often used to assess the strength of a surrogate relationship [20, 21,
22]. For example, the German Institute for Quality and Efficiency in
Healthcare (IQWIG) defined an acceptable surrogate endpoint by setting a lower
bound for the confidence interval on the correlation coefficient to be 0.85
[23]. To be able to estimate the between-studies correlation and the study-
level $R^{2}$, we used BRMA PNF as an alternative method to model the
surrogacy pattern. BRMA PNF has the same within-study model as the D&H model
(2) but the between-studies model assumes exchangeability of the correlated
true (random) treatment effects on both outcomes. In the PNF, the bivariate
normal distribution is represented as a sequence of univariate conditional
distributions (4):

$\left\\{\begin{array}[]{l}\delta_{1i}\sim
N\left(\eta_{1},\psi_{1}^{2}\right)\\\ \delta_{2i}\mid\delta_{1i}\sim
N\left(\eta_{2i},\psi_{2}^{2}\right)\\\
\eta_{2i}=\lambda_{0}+\lambda_{1}\delta_{1i}\end{array}\right.$ (4)

where $\delta_{1i}$ and $\delta_{2i}$ are the true effects in the population
which are correlated, assumed exchangeable and normally distributed.

The parameters of the BRMA PNF model can be represented in terms of the
parameters of the bivariate normal distribution using the following formulae
[24].

$\psi_{1}^{2}=\tau_{1}^{2},\psi_{2}^{2}=\tau_{2}^{2}-\lambda_{1}^{2}\tau_{1}^{2},\lambda_{1}=\frac{\tau_{2}}{\tau_{1}}\rho,$
(5)

where $\tau_{1}$ and $\tau_{2}$ are the between-studies heterogeneity
parameters and $\rho$ is the between-studies correlation. To implement this
model in a Bayesian framework, vague prior distributions were placed on the
between-studies parameters
$\mathrm{\tau}_{1,2}\sim\mathrm{Unif}\left(0,2\right)$,
$\mathrm{\rho}\sim\mathrm{Unif}\left(-1,1\right)$ and the intercept
$\mathrm{\lambda_{0}}\sim\mathrm{N}\left(0,10^{4}\right)$ implying prior
distributions on $\lambda_{1}$, $\psi_{1}^{2}$ and $\psi_{2}^{2}$ through
rearranging the relationships (5).

In addition to the surrogacy criteria from the D&H model, a perfect surrogate
relationship is defined in the BRMA PNF model when $\rho=\pm 1$ [25]. This
implies a perfect association between treatment effects on the surrogate
endpoint and final outcome. The study-level $R^{2}$ in this random-effects
model is equal to $\rho^{2}$ [22].

####  2.4.3 Bias Adjustment

To account for systematic differences in treatment effects between RCTs, cRWE
and sRWE, the BRMA PNF model was extended to allow for bias adjustment. The
between-studies model (4) remains the same for all studies and the within-
study model for RCTs remains the same as in (2). However, the within-study
models for cRWE and matched sRWE include bias terms, $\alpha_{j}$ and
$\beta_{j}$ respectively, for the surrogate and final outcomes ($j=1,2$):

$\left(\begin{array}[]{c}Y_{1i}\\\ Y_{2i}\end{array}\right)\sim
N\left(\left(\begin{array}[]{l}\delta_{1i}+\alpha_{1}\\\
\delta_{2i}+\alpha_{2}\end{array}\right),\left(\begin{array}[]{cc}\sigma_{1i}^{2}&\sigma_{1i}\sigma_{2i}\rho_{wi}\\\
\sigma_{1i}\sigma_{2i}\rho_{wi}&\sigma_{2i}^{2}\end{array}\right)\right)$ (6)

$\left(\begin{array}[]{c}Y_{1i}\\\ Y_{2i}\end{array}\right)\sim
N\left(\left(\begin{array}[]{l}\delta_{1i}+\beta_{1}\\\
\delta_{2i}+\beta_{2}\end{array}\right),\left(\begin{array}[]{cc}\sigma_{1i}^{2}&\sigma_{1i}\sigma_{2i}\rho_{wi}\\\
\sigma_{1i}\sigma_{2i}\rho_{wi}&\sigma_{2i}^{2}\end{array}\right)\right)$ (7)

Additional non-informative prior distributions were placed on the bias terms,
$\mathrm{\alpha}_{1,2}\sim\mathrm{N}\left(0,10^{4}\right)$ and
$\mathrm{\beta}_{1,2}\sim\mathrm{N}\left(0,10^{4}\right)$.

####  2.4.4 Cross-validation

To assess whether addition of RWE improves accuracy or precision of
predictions, a “take-one-out” cross-validation procedure was conducted for the
D&H and BRMA PNF models. For each study $i~{}(i=1,...,N)$, the treatment
effect on the final outcome, $Y_{2i}$, was removed and assumed missing at
random. The treatment effect on the final outcome was predicted from the
treatment effect on the surrogate endpoint, conditional on data on both
outcomes from all other studies in the meta-analysis. The mean predicted
effect is equal to the mean predicted true effect from MCMC simulation. The
variance of the predicted effect is
$\sigma_{2i}^{2}+\operatorname{var}\left(\hat{\delta}_{2i}\mid
Y_{1i},\sigma_{1i},Y_{1(-i)},Y_{2(-i)}\right)$ where $Y_{1,2(-i)}$ are the
observed treatment effects on the surrogate and final outcomes for the
remaining studies not omitted in the $ith$ iteration [19, 3]. For a valid
surrogate, the 95% predicted interval (constructed using the variance) will
contain the observed treatment effect in at least 95% of studies.

###  2.5 Software and Computing

All models were implemented using WinBUGS [26] where estimates were obtained
using MCMC simulation with 150000 iterations (including 50000 burn-in).
Convergence was checked via visual assessment of history, density and
autocorrelation plots. Posterior estimates are presented as means
(approximately normal posterior) or medians (skewed posterior) with 95%
credible intervals (CrI). R was used for data manipulation, to execute WinBUGS
code using the R2WinBUGS package [27] and to produce figures using the ggplot2
package.

##  3 Results

###  3.1 Summary of Data

Of the 11 RCTs obtained from the prior literature review, 4 were excluded for
not investigating the effect of anti-angiogenics in combination with
chemotherapy against chemotherapy alone. Overall 7 RCTs were included in the
analysis.

The database search of PubMed returned 166 publications for cRWE studies and
145 publications for sRWE studies on the chemotherapy arm. After screening
titles, abstracts and, where appropriate, full articles, 7 cRWE studies
comparing bevacizumab against chemotherapy remained and 8 sRWE studies of
chemotherapy alone remained. Of the 7 cRWE studies, 4 adjusted for covariates.

Five covariates were reported in all sRWE studies and of these 5, sex was the
only covariate not recommended for reporting in the consensus statement.
Following the consensus statement ranking, sex was given a weight of 1 and all
other covariates a weight of 2. The covariates selected for matching were:
treatment line (weight=2, current treatment line normalised between range
1-3), age (weight=2, median age normalised between range 18-100), performance
score (weight=2, ECOG/WHO score normalised between range 0-3), tumour location
(weight=2, proportion with colon tumour compared to rectum tumour) and sex
(weight=1, proportion of females).

Table 1 shows the distance measures between the sRWE studies. A maximum
distance measure of 0.035 was applied as this was close to the maximum
distance measure from RCTs (0.032). This resulted in an exploration of 9%
(n=11) of possible matches. In Table 1, possible matches are shaded and final
matches (lowest distance measures) are shown in bold. Overall, three matched
sRWE studies were included in the analysis.

Table 1: Distance Metric between Single Arm Observational Studies

|  | Dong
---
2015
| Matsumoto
---
2007
| Catalano
---
2009
| Fuse
---
2008
| Suenaga
---
2008
| Fuse
---
2007
| Hochster
---
2003
| Yoshino
---
2007
Bendell 2012 (1) | 0.029 | 0.273 | 0.048 | 0.134 | 0.077 | 0.126 | 0.043 | 0.057
Bendell 2012 (2) | 0.034 | 0.273 | 0.043 | 0.139 | 0.074 | 0.131 | 0.041 | 0.062
Hurwitz 2014 | 0.151 | 0.156 | 0.158 | 0.211 | 0.171 | 0.248 | 0.169 | 0.179
Van Cutsem 2009 (1) | 0.051 | 0.234 | 0.092 | 0.093 | 0.040 | 0.085 | 0.073 | 0.018
Van Cutsem 2009 (2) | 0.034 | 0.232 | 0.085 | 0.082 | 0.044 | 0.074 | 0.083 | 0.018
Van Cutsem 2009 (3) | 0.047 | 0.221 | 0.096 | 0.072 | 0.045 | 0.064 | 0.094 | 0.013
Van Cutsem 2009 (4) | 0.051 | 0.224 | 0.095 | 0.079 | 0.045 | 0.071 | 0.087 | 0.017
Bennouna 2017 (1) | 0.058 | 0.227 | 0.069 | 0.110 | 0.073 | 0.102 | 0.089 | 0.041
Bennouna 2017 (2) | 0.208 | 0.115 | 0.208 | 0.175 | 0.204 | 0.212 | 0.242 | 0.191
Buchler 2014 (1) | 0.045 | 0.233 | 0.079 | 0.094 | 0.070 | 0.086 | 0.080 | 0.025
Buchler 2014 (2) | 0.052 | 0.213 | 0.088 | 0.082 | 0.050 | 0.072 | 0.091 | 0.014
Ocvirk 2011 (1) | 0.043 | 0.260 | 0.072 | 0.100 | 0.066 | 0.094 | 0.066 | 0.065
Ocvirk 2011 (2) | 0.051 | 0.226 | 0.097 | 0.080 | 0.047 | 0.062 | 0.096 | 0.044
Moriwaki 2012 (1) | 0.204 | 0.062 | 0.233 | 0.121 | 0.160 | 0.148 | 0.243 | 0.165
Moriwaki 2012 (2) | 0.184 | 0.099 | 0.228 | 0.132 | 0.129 | 0.169 | 0.211 | 0.150
Kotaka 2016 | 0.068 | 0.200 | 0.095 | 0.072 | 0.033 | 0.064 | 0.104 | 0.034

Figure 1 presents the data from the included studies as a scatter plot of the
observed treatment effects on PFS and OS. The plot shows a possible strong
positive relationship between treatment effects on PFS and OS, suggesting a
potential valid surrogate relationship.

Figure 1: Scatterplot of observed logHR on OS against observed logHR on PFS
for RCTs (blue circle), cRWE (red square) and matched sRWE (green triangle).
Black line shows the linear relationship between logHR on PFS and logHR on OS
obtained from D&H model conducted using all sources of evidence. Figure 2:
Forest plots of HRs from RCTs, cRWE and sRWE. Graph shows observed HRs and 95%
confidence intervals on PFS (black lines - left) and OS (black lines - right)
and predicted HRs and 95% predicted intervals on OS obtained from cross-
validation using D&H model (grey lines - right).

###  3.2 Daniels and Hughes Model

Table 2 shows results from the D&H model. There is evidence of a surrogate
relationship regardless of the type of evidence used as the intercept and
conditional variance are close to zero and the slope is non-zero. Figure 2
illustrates this relationship, highlighting that studies with larger treatment
effects on PFS generally have larger effects on OS.

Table 2: Surrogacy criteria obtained from D&H model applied to data from RCTs,
comparative RWE (cRWE) and matched single-arm RWE (sRWE). Last two rows give
cross-validation results from D&H model where absolute discrepancy is the
absolute difference between the observed logHR and the predicted logHR on OS,
while $w_{\hat{Y}_{2j}}/w_{Y_{2j}}$ is the ratio of the width of the 95%
predicted interval of the logHR on OS to the width of the 95% confidence
interval of the observed estimate of the logHR on OS.

| RCTs | RCTs & cRWE | RCTs, cRWE & sRWE
---|---|---|---
$\lambda_{0}$ | 0.10 (-0.13, 0.34) | 0.051 (-0.13, 0.23) | -0.015 (-0.20, 0.15)
$\lambda_{1}$ | 0.71 (0.13, 1.30) | 0.69 (0.20, 1.17) | 0.55 (0.099, 0.97)
$\psi_{2}^{2}$ | 0.0091 (0.0001, 0.11) | 0.010 (0.0002, 0.051) | 0.010 (0.0002, 0.054)
Absolute Discrepancy, Median (Range) | 0.091 (0.0022, 0.27) | 0.070 (0.0079, 0.30) | 0.076 (0.0019, 0.73)
$w_{\hat{Y}_{2j}}/w_{Y_{2j}}$, Median (Range) | 2.90 (1.76, 3.58) | 2.09 (1.33, 2.40) | 1.79 (1.22, 2.42)

Addition of cRWE to RCTs improved the precision of all three estimates for the
surrogate relationship while having minimal impact on the point estimates. For
example, using RCT data alone resulted in the conditional variance 0.0091 (95%
CrI: 0.0001, 0.11), addition of cRWE gave a conditional variance of 0.010 (95%
CrI: 0.0002, 0.051). Thus, addition of cRWE reduced uncertainty by 54% in
terms of the credible interval width. Addition of sRWE to RCTs and cRWE
further improved precision of the intercept and slope estimates and only
slightly decreased precision of the conditional variance.

###  3.3 Bivariate Random-Effects Meta-Analysis (Product Normal Formulation)

Table 3 shows there was weaker evidence for a surrogate relationship using the
BRMA PNF compared to the D&H model. Although correlation was 0.75 (when using
all evidence) the credible interval was wide (95% CrI: -0.081, 0.98).
Furthermore, the credible interval for the slope contained zero, suggesting no
relationship between treatment effects on the surrogate and final endpoint.
Such differences in results between the models could be a result of the random
effects assumption of the BRMA PNF model. When assuming normal random effects
is appropriate, greater borrowing of information can lead to more precise
estimates. However, when this assumption is violated, the model can lead to
over-shrinkage of the true effects, thus potentially reducing the between-
studies correlation [28]. This can lead to bias and reduce precision of
estimates for surrogacy parameters.

Despite differing results to the D&H model, addition of RWE generally improved
the precision of estimates obtained from the BRMA PNF model. The correlation
obtained using RCT data only was 0.75 (95% CrI: -0.24, 0.99) whereas
correlation obtained using all sources of evidence was 0.75 (95% CrI: -0.081,
0.98). Thus, addition of RWE reduced uncertainty by 14% in terms of the width
of the credible interval while the point estimate remained the same.

Table 3: Surrogacy criteria obtained from BRMA PNF model applied to data from
RCTs, comparative RWE (cRWE) and matched single-arm RWE (sRWE). Last two rows
give cross-validation results from BRMA PNF model where absolute discrepancy
is the absolute difference between the observed logHR and the predicted logHR
on OS, while $w_{\hat{Y}_{2j}}/w_{Y_{2j}}$ is the ratio of the width of the
95% predicted interval of the logHR on OS to the width of the 95% confidence
interval of the observed estimate of the logHR on OS.

|  | RCTs
---
| RCTs & cRWE
---
| RCTs, cRWE & sRWE
---
| Bias Adjusted
---
RCTs, cRWE & sRWE
$d_{1}$ | -0.36 (-0.62, -0.13) | -0.34 (-0.49, -0.21) | -0.37 (-0.51, -0.23) | -0.36 (-0.57, -0.16)
$d_{2}$ | -0.14 (-0.35, 0.046) | -0.18 (-0.31, -0.052) | -0.22 (-0.34, -0.11) | -0.14 (-0.29, -0.0020)
$\rho$ | 0.75 (-0.24, 0.99) | 0.66 (-0.30, 0.97) | 0.75 (-0.081, 0.98) | 0.69 (-0.32, 0.99)
$\tau_{1}$ | 0.27 (0.11, 0.59) | 0.19 (0.079, 0.36) | 0.22 (0.11, 0.38) | 0.23 (0.10, 0.43)
$\tau_{2}$ | 0.21 (0.075, 0.48) | 0.16 (0.072, 0.31) | 0.17 (0.086, 0.30) | 0.15 (0.048, 0.28)
$\lambda_{0}$ | 0.054 (-0.21, 0.33) | 0.0068 (-0.26, 0.28) | -0.012 (-0.25, 0.22) | 0.011 (-0.22, 0.26)
$\lambda_{1}$ | 0.54 (-0.16, 1.31) | 0.54 (-0.25, 1.40) | 0.56 (-0.053, 1.19) | 0.42 (-0.17, 1.14)
$\psi_{2}^{2}$ | 0.013 (0.0009, 0.10) | 0.011 (0.0014, 0.048) | 0.011 (0.0010, 0.050) | 0.0077 (0.0003, 0.047)
$R^{2}$ | 0.57 (0.0045, 0.98) | 0.43 (0.0023, 0.94) | 0.56 (0.0061, 0.97) | 0.49 (0.0022, 0.97)
$\alpha_{1}$ | - | - | - | 0.041 (-0.32, 0.38)
$\alpha_{2}$ | - | - | - | -0.12 (-0.36, 0.12)
$\beta_{1}$ | - | - | - | -0.13 (-0.52, 0.29)
$\beta_{2}$ | - | - | - | -0.28 (-0.59, 0.039)
| Absolute Discrepancy,
---
Median (Range)
0.16 (0.017, 0.26) | 0.16 (0.0046, 0.23) | 0.12 (0.026, 0.57) | 0.17 (0.0092, 0.66)
|  $w_{\hat{Y}_{2j}}/w_{Y_{2j}}$,
---
Median (Range)
2.74 (1.61, 3.30) | 2.09 (1.15, 2.49) | 1.78 (1.14, 2.56) | 1.66 (1.14, 3.07)

###  3.4 Bias Adjustment

Table 3 shows the results of the bias adjusted BRMA PNF model. When adjusting
for bias, there was no improvement in precision and the 95% credible interval
for the slope estimate still contained zero, indicating weak evidence for
surrogacy. In addition, when using bias adjustment, correlation was lower
(0.69; 95% CrI: -0.32, 0.99) than when not using bias adjustment (0.75; 95%
CrI: -0.081, 0.98), providing less evidence for a valid surrogate
relationship. The lower correlation could be explained by the bias adjustment
reducing between-studies heterogeneity for treatment effects on the final
outcome.

###  3.5 Cross-validation

The last two rows of Table 2 show the results of cross-validation using the
D&H model. Median absolute discrepancy between predicted and observed
treatment effects on OS decreased with the addition of cRWE and slightly
increased with the further addition of sRWE. This increase is likely explained
by matched Suenaga 2008 and Kotaka 2016 studies which had a very low relative
treatment effect on OS in comparison to the treatment effect on PFS (Figures 1
and 2). Table 2 also shows that the length of the predicted interval, relative
to the observed confidence interval, fell with the addition of cRWE (2.90 to
2.09) and sRWE (2.09 to 1.79). Overall, the precision of the predicted
intervals improved with the addition of RWE. Cross-validation for the BRMA PNF
model showed similar results (Table 3).

##  4 Discussion

When existing clinical trial data are limited, surrogate endpoint validation
may fail. As a result, new therapies may not receive conditional marketing
authorisation or, if approved, they may still fail at the health technology
assessment (HTA) decision-making stage, by HTA agencies such as the National
Institute for Health and Care Excellence (NICE) [2, 3]. In this paper, we
provide an approach for using RWE to strengthen the evidence base for
surrogate endpoint validation.

There are several limitations of this research. Inclusion of sRWE relied on
digitizing Kaplan-Meier curves. However, such curves are not always published
and thus, potentially useful studies would not be included. One method to
overcome this issue is to extract median survival times for the surrogate and
final outcomes and use an exponential hazard assumption, as proposed by
Schmitz et al [15], to obtain treatment effects. In a preliminary analysis we
applied this method to the mCRC dataset; however, the exponential hazard
assumption did not provide a good approximation.

Matching of sRWE was based on study level covariates and thus was prone to
bias as patients were not randomised or matched at the individual-level and,
therefore, the assumption of exchangeability may have been violated [29]. This
bias was further exacerbated by only matching on 5 covariates, when 10
additional characteristics were recommended by the consensus statement and
identified as risk factors. However, these variables were not reported in the
included studies.

While the bias adjustment was used to account for potential systematic
differences in treatment effects between data sources, all sources of evidence
contributed the same weight to the model. This suggests that RCTs, cRWE and
matched sRWE were of equivalent quality. However, it is widely acknowledged
that RCTs are above RWE in the hierarchy of evidence [30]. Further
methodological research is carried out to allow for accounting for such
differences in quality.

##  5 Conclusions

RWE can be used to improve the precision of estimates for surrogate endpoint
validation relative to using RCT data alone. The addition of RWE to RCT data
also allows for more precise predictions to be made of the treatment effects
on the final clinical outcome based on the treatment effect measured on the
surrogate endpoint.

## Acknowledgements

LW was funded by a National Institute for Health Research (NIHR) Pre-doctoral
Fellowship [NIHR301013]. The views expressed are those of the authors and not
necessarily those of the NIHR or the Department of Health and Social Care. The
funding bodies played no role in the design of the study and collection,
analysis, and interpretation of data and in writing the manuscript. SB was
funded by Medical Research Council [grant no. MR/R025223/1].

## Author Contributions

SB conceptualised the study and provided supervision to LW and TP. LW and TP
curated the data. LW carried out all formal analyses, methodological
development and visualisation, and provided the original draft of the
manuscript. AT provided clinical advice. All authors contributed to the
manuscript revisions.

## References

  * [1] T. Burzykowski, G. Molenberghs, and M. Buyse.  The evaluation of surrogate endpoints.  Springer, 2006.
  * [2] Dalia Dawoud, Huseyin Naci, Oriana Ciani, and Sylwia Bujkiewicz.  Raising the bar for using surrogate endpoints in drug regulation and health technology assessment.  bmj, 374, 2021.
  * [3] S Bujkiewicz, F Achana, T Papanikos, RD Riley, and KR Abrams.  Nice dsu technical support document 20: multivariate meta-analysis of summary data for combining treatment effects on correlated outcomes and evaluating surrogate endpoints, 2019.
  * [4] Julian Elston and Rod S Taylor.  Use of surrogate outcomes in cost-effectiveness models: A review of united kingdom health technology assessment reports1.  International journal of technology assessment in health care, 25(1):6–13, 2009.
  * [5] Brett K Beaulieu-Jones, Samuel G Finlayson, William Yuan, Russ B Altman, Isaac S Kohane, Vinay Prasad, and Kun-Hsing Yu.  Examining the use of real-world evidence in the regulatory process.  Clinical Pharmacology & Therapeutics, 107(4):843–852, 2020.
  * [6] Robert Flynn, Kelly Plueschke, Chantal Quinten, Valerie Strassmann, Ruben G Duijnhoven, Maria Gordillo-Marañon, Marcia Rueckbeil, Catherine Cohet, and Xavier Kurz.  Marketing authorisation applications made to the european medicines agency in 2018–2019: what was the contribution of real-world evidence?  Clinical Pharmacology & Therapeutics, 2021.
  * [7] Marc L Berger, Harold Sox, Richard J Willke, Diana L Brixner, Hans-Georg Eichler, Wim Goettsch, David Madigan, Amr Makady, Sebastian Schneeweiss, Rosanna Tarricone, et al.  Good practices for real-world data studies of treatment and/or comparative effectiveness: recommendations from the joint ispor-ispe special task force on real-world evidence in health care decision making.  Value in Health, 20(8):1003–1008, 2017.
  * [8] Amr Makady, Ard van Veelen, Páll Jonsson, Owen Moseley, Anne D’Andon, Anthonius de Boer, Hans Hillege, Olaf Klungel, and Wim Goettsch.  Using real-world data in health technology assessment (hta) practice: a comparative study of five hta agencies.  Pharmacoeconomics, 36(3):359–368, 2018.
  * [9] Pablo E Verde and Christian Ohmann.  Combining randomized and non-randomized evidence in clinical research: a review of methods and applications.  Research synthesis methods, 6(1):45–62, 2015.
  * [10] Orestis Efthimiou, Dimitris Mavridis, Thomas PA Debray, Myrto Samara, Mark Belger, George CM Siontis, Stefan Leucht, Georgia Salanti, and GetReal Work Package 4.  Combining randomized and non-randomized evidence in network meta-analysis.  Statistics in medicine, 36(8):1210–1226, 2017.
  * [11] Nicky J Welton, Phillippo David M, Rhiannon Owen, Hayley J Jones, Sofia Dias, Sylwia Bujkiewicz, Anthony E Ades, and Keith R Abrams.  Chte2020 sources and synthesis of evidence; update to evidence synthesis methods.  DSU Report, March 2020.
  * [12] Susanne Schmitz, Roisin Adams, and Cathal Walsh.  Incorporating data from various trial designs into a mixed treatment comparison model.  Statistics in medicine, 32(17):2935–2949, 2013.
  * [13] Oriana Ciani, Marc Buyse, Ruth Garside, Jaime Peters, Everardo D Saad, Ken Stein, and Rod S Taylor.  Meta-analyses of randomized controlled trials show suboptimal validity of surrogate outcomes for overall survival in advanced colorectal cancer.  Journal of clinical epidemiology, 68(7):833–842, 2015.
  * [14] Anastasios Papanikos.  Bayesian evidence synthesis for surrogate endpoints in the era of precision medicine.  PhD thesis, University of Leicester, 2021.
  * [15] Susanne Schmitz, Áine Maguire, James Morris, Kai Ruggeri, Elisa Haller, Isla Kuhn, Joy Leahy, Natalia Homer, Ayesha Khan, Jack Bowden, et al.  The use of single armed observational data to closing the gap in otherwise disconnected evidence networks: a network meta-analysis in multiple myeloma.  BMC medical research methodology, 18(1):1–18, 2018.
  * [16] Kaitlyn KH Goey, Halfdan Sørbye, Bengt Glimelius, Richard A Adams, Thierry André, Dirk Arnold, Jordan D Berlin, György Bodoky, Aimery de Gramont, Eduardo Díaz-Rubio, et al.  Consensus statement on essential patient characteristics in systemic treatment trials for metastatic colorectal cancer: Supported by the arcad group.  European Journal of Cancer, 100:35–45, 2018.
  * [17] Jayne F Tierney, Lesley A Stewart, Davina Ghersi, Sarah Burdett, and Matthew R Sydes.  Practical methods for incorporating summary time-to-event data into meta-analysis.  Trials, 8(1):1–16, 2007.
  * [18] Yinghui Wei and Patrick Royston.  Reconstructing time-to-event data from published kaplan–meier curves.  The Stata Journal, 17(4):786–802, 2017.
  * [19] Michael J Daniels and Michael D Hughes.  Meta-analysis for the evaluation of potential surrogate markers.  Statistics in medicine, 16(17):1965–1982, 1997.
  * [20] M. Buyse, G. Molenberghs, T. Burzykowski, D. Renard, and H. Geys.  The validation of surrogate endpoints in meta-analyses of randomized experiments.  Biostatistics, 1(1):49–67, 2000.
  * [21] T. Burzykowski, G. Molenberghs, M. Buyse, H. Geys, and D. Renard.  Validation of surrogate end points in multiple randomized clinical trials with failure time end points.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 50(4):405–422, 2001.
  * [22] L. A. Renfro, Q. Shi, D. J. Sargent, and B. P. Carlin.  Bayesian adjusted r2 for the meta-analytic evaluation of surrogate time-to-event endpoints in clinical trials.  Statistics in Medicine, 31(8):743–761, 2012.
  * [23] Institute for Quality and Efficiency in Health Care.  Validity of surrogate endpoints in oncology. executive summary of rapid report a10-05, version 1.1., https://www.ncbi.nlm.nih.gov/books/nbk198799/, Published 2011. Accessed December 3, 2021.
  * [24] Sylwia Bujkiewicz, John R Thompson, Richard D Riley, and Keith R Abrams.  Bayesian meta-analytical methods to incorporate multiple surrogate endpoints in drug development process.  Statistics in medicine, 35(7):1063–1089, 2016.
  * [25] Sylwia Bujkiewicz, Dan Jackson, John R Thompson, Rebecca M Turner, Nicolas Städler, Keith R Abrams, and Ian R White.  Bivariate network meta-analysis for surrogate endpoint evaluation.  Statistics in medicine, 38(18):3322–3341, 2019.
  * [26] David J Lunn, Andrew Thomas, Nicky Best, and David Spiegelhalter.  Winbugs-a bayesian modelling framework: concepts, structure, and extensibility.  Statistics and computing, 10(4):325–337, 2000.
  * [27] Sibylle Sturtz, Uwe Ligges, and Andrew E Gelman.  R2winbugs: a package for running winbugs from r.  2005\.
  * [28] Sylwia Bujkiewicz, John R Thompson, Enti Spata, and Keith R Abrams.  Uncertainty in the bayesian meta-analysis of normally distributed surrogate endpoints.  Statistical methods in medical research, 26(5):2287–2318, 2017\.
  * [29] Joy Leahy, Howard Thom, Jeroen P Jansen, Emma Gray, Aisling O’Leary, Arthur White, and Cathal Walsh.  Incorporating single-arm evidence into a network meta-analysis using aggregate level matching: assessing the impact.  Statistics in medicine, 38(14):2505–2523, 2019.
  * [30] Brian Brighton, Mohit Bhandari, Paul Tornetta, David T Felson, et al.  Hierarchy of evidence: from case reports to randomized controlled trials.  Clinical Orthopaedics and Related Research (1976-2007), 413:19–24, 2003.


# An approximate expectation-maximization for two-dimensional multi-target
detection

Shay Kreymer, Amit Singer, and Tamir Bendory S. Kreymer and T. Bendory are
with the School of Electrical Engineering of Tel Aviv University, Tel Aviv,
Israel, e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>A. Singer
is with the Department of Mathematics and PACM, Princeton University,
Princeton, NJ, USA, e-mail<EMAIL_ADDRESS>S.K. is supported by the
Yitzhak and Chaya Weinstein Research Institute for Signal Processing. A.S. is
partly supported by AFOSR Award FA9550-20-1-0266, the Simons Foundation Math+X
Investigator Award, NSF BIGDATA Award IIS1837992, NSF Award DMS-2009753, and
NIH/NIGMS Award R01GM136780- 01. T.B. is partially supported by the BSF grant
no. 2020159, the NSF-BSF grant no. 2019752, and the ISF grant no. 1924/21.

###### Abstract

We consider the two-dimensional multi-target detection (MTD) problem of
estimating a target image from a noisy measurement that contains multiple
copies of the image, each randomly rotated and translated. The MTD model
serves as a mathematical abstraction of the structure reconstruction problem
in single-particle cryo-electron microscopy, the chief motivation of this
study. We focus on high noise regimes, where accurate detection of image
occurrences within a measurement is impossible. To estimate the image, we
develop an expectation-maximization framework that aims to maximize an
approximation of the likelihood function. We demonstrate image recovery in
highly noisy environments, and show that our framework outperforms the
previously studied autocorrelation analysis in a wide range of parameters.

###### Index Terms:

Expectation-maximization, multi-target detection, cryo-electron microscopy.

##  I Introduction

We study the multi-target detection (MTD) problem of estimating a target image
$f:\mathbb{R}^{2}\rightarrow\mathbb{R}$ from a noisy measurement that contains
multiple copies of the image, each randomly rotated and translated [1, 2, 3,
4, 5, 6, 7]. We consider a measurement $M\in\mathbb{R}^{N\times N}$ of the
form

$M[\vec{\ell}]=\sum_{i=1}^{p}F_{\phi_{i}}[\vec{\ell}-\vec{\ell}_{i}]+\varepsilon[\vec{\ell}],$
(1)

where $F_{\phi_{i}}[\vec{\ell}]:=f_{\phi_{i}}(\vec{\ell}/n)$ is a discrete
copy of $f$, rotated by angle $\phi_{i}$ about the origin; $n$ is the known
radius of the image in pixels;
$\\{\phi_{i}\\}_{i=1}^{p}\sim\operatorname{Unif}[0,2\pi)$ are uniformly
distributed rotations;
$\\{\vec{\ell}_{i}\\}_{i=1}^{p}\in\\{n+1,\ldots,N-n\\}^{2}$ are identically
distributed random translations (though not independent, and the underlying
distribution is not assumed to be known); and $\varepsilon[\vec{\ell}]$ is
i.i.d. Gaussian noise with zero mean and variance $\sigma^{2}$. The rotations,
translations, and the number of occurrences of $f$ in $M$, denoted by $p$, are
unknown. Importantly, since the rotations are unknown, it is possible to
reconstruct the target image only up to a rotation.

Following [3, 4, 5, 8], we assume that the image $f$ is supported on the unit
disk $\\{\vec{x}\in\mathbb{R}^{2}:|\vec{x}|\leq 1\\}$ and has a finite
expansion in the basis of Dirichlet Laplacian eigenfunctions. In particular,
the image $f$ can be expanded as

$f(r,\theta)=\sum_{(\nu,q):\lambda_{\nu,q}\leq\lambda}\alpha_{\nu,q}\psi_{\nu,q}(r,\theta),\quad
r\leq 1,$ (2)

in polar coordinates $(r,\theta)$, where
$\psi_{\nu,q}(r,\theta)=J_{\nu}\left(\lambda_{\nu,q}r\right)e^{i\nu\theta}$,
$\nu\in\mathbb{Z}_{\geq 0}$, $J_{\nu}$ is the $\nu$-th order Bessel function
of the first kind, $\lambda_{\nu,q}>0$ is the $q$-th positive root of
$J_{\nu}$, $\lambda$ is called the bandlimit frequency, and $\alpha$ is the
vector of expansion coefficients. Hereafter, by estimating the image we mean
estimating the vector of coefficients $\alpha$. Notably, the basis of
Dirichlet Laplacian eigenfunctions is steerable: rotating $f$ is equivalent to
modulating the expansion coefficients $\alpha_{\nu,q}$. Specifically, the
expansion of the rotated image $f_{\phi}(r,\theta):=f(r,\theta+\phi)$ is given
by
$f_{\phi}(r,\theta)=\sum_{(\nu,q):\lambda_{\nu,q}\leq\lambda}\alpha_{\nu,q}\psi_{\nu,q}(r,\theta)e^{i\nu\phi}$.
For real valued images, $\alpha_{-\nu,q}=\alpha_{\nu,q}^{*}$.

We focus on the well-separated case of the 2-D MTD problem, which was
introduced in [3, 4]. In this case, each translation in the measurement $M$
(1) is separated by at least a full image diameter from its neighbors.
Specifically, we assume

$|\vec{\ell}_{i_{1}}-\vec{\ell}_{i_{2}}|>4n,\quad\text{ for all }i_{1}\neq
i_{2}.$ (3)

Figure 1 presents an example of a measurement $M$ at different signal-to-noise
ratios (SNRs). We define
$\text{SNR}:=\frac{\|F\|_{\text{F}}^{2}}{A\sigma^{2}}$, where $A$ is the area
in pixels of $F$.

(a) No noise

(b) $\text{SNR}=50$

(c) $\text{SNR}=2$

Figure 1: Three measurements at different SNRs: (a) no noise; (b)
$\text{SNR}=50$; (c) $\text{SNR}=2$. Each measurement contains multiple
rotated versions of the target image. We focus on the low SNR regime (e.g.,
panel (c)) in which the locations and rotations of the image occurrences
cannot be detected reliably.

The MTD model serves as a mathematical abstraction of the cryo-electron
microscopy (cryo-EM) technology for macromolecular structure determination [9,
10, 11]. In a cryo-EM experiment, individual copies of the target biomolecule
are dispersed at unknown 2-D locations and 3-D orientations in a thin layer of
vitreous ice, from which 2-D tomographic projection images are produced by an
electron microscope [12]. It is necessary to keep the electron dose low in
order to minimize irreversible structural damage. Consequently, the projection
images are considerably noisy. In the current data processing pipeline of
cryo-EM [13, 14, 15, 16], the 2-D projections are first detected and extracted
from the micrograph, and later rotationally and translationally aligned to
reconstruct the 3-D molecular structure. This approach fails for small
molecules, which are difficult to detect and align [7, 9, 13, 17].

The MTD model was devised in [7] in order to study the recovery of small
molecular structures using cryo-EM, below the current detection limit [18]. In
[3, 4, 5], an autocorrelation analysis technique was devised for the 2-D MTD
problem (1). Autocorrelation analysis is a special case of the method of
moments, and it consists of finding an image that best matches the empirical
autocorrelations of the measurement, thus bypassing detecting the locations
and rotations of individual image occurrences. In this work, we propose to
replace autocorrelation analysis by expectation-maximization (EM): a classical
iterative algorithm to compute the maximum likelihood estimator [19].
Similarly to autocorrelation analysis, the EM algorithm estimates the target
image $F$ directly by marginalizing over the translations and rotations.

Previous works demonstrated that EM outperforms autocorrelation analysis in
terms of estimation accuracy for the 1-D MTD problem [2], as well as for the
closely related multireference alignment model [20, 21]. EM is also the most
popular computational framework for reconstructing molecular structures using
cryo-EM, see for example [15, 16]. Moreover, a recent paper [22] shows that
likelihood optimization in the low SNR regime reduces to a sequence of least
squares optimization problems that match the moments of the image estimate to
the observable moments one by one, and by that suggests that EM has the
potential to surpass the estimation accuracy achieved by autocorrelation
analysis.

At each iteration, EM assigns probabilities to all possible rotations and
translations (see Section II). Unfortunately, for the MTD model (1), the
number of possible translations grows quickly with the measurement size, and
therefore direct application of the EM algorithm to the MTD problem is
computationally intractable, even for very small measurements. Thus, following
[2], we suggest mitigating the computational burden by developing an EM
algorithm that maximizes an approximation of the likelihood function. In the
approximate EM scheme, the number of possible translations is linear in
$N^{2}$ (the size of the measurement), making the algorithm tractable.

The main contribution of this paper is in developing an approximate EM
framework for the 2-D MTD problem; see Section II. In Section III, we
demonstrate successful reconstructions in noisy regimes. We also conduct
extensive numerical experiments that demonstrate significant improvement in
estimation accuracy compared to autocorrelation analysis. Section IV concludes
the paper and introduces future work. In particular, we discuss potential
implications for the ongoing effort to estimate small molecular structures
using cryo-EM [7].

##  II Approximate expectation-maximization

Given a measurement $M$ that follows the MTD model (1), the maximum marginal
likelihood estimator (MMLE) for the vector of coefficients $\alpha$, that
represents the target image $f$ (2), is the maximizer of the likelihood
function $p(M|\alpha)$. The translations and rotations of the target images
within the measurement are treated as nuisance variables. The EM algorithm
estimates the MMLE by iteratively applying the expectation (E) and
maximization (M) steps [19]. Specifically, given the current estimate
$\alpha_{k}$, the E-step computes the expected log-likelihood function, where
the expectation is taken over all admissible configurations of translations
and rotations. The estimate is then updated in the M-step by maximizing the
function with respect to $\alpha$. Unfortunately, for the MTD model, the
number of possible translations grows quickly with the problem size, rendering
direct application of EM computationally intractable. Hence, based on [2], we
suggest to apply an approximate EM, in which the number of possible
translations is linear in $N^{2}$.

The approximate EM begins with partitioning the measurement $M$ into
$N_{d}=(N/L)^{2}$ non-overlapping patches; each patch is of size $L\times L$,
where $L=2n+1$ is the diameter of the target image $F$. In EM terminology, the
patches are called the observed data. The separation condition (3) implies
that each patch $M_{m}$ can contain either no target image, a full target
image, or part of a rotated image; overall there are $(2L)^{2}$ possibilities,
excluding rotations. In particular, each patch can be modeled by

$M_{m}=CT_{\vec{\ell}_{m}}ZF^{{{L}}}_{\phi_{m}}+\varepsilon_{m},\quad\varepsilon_{m}\sim\mathcal{N}(0,\sigma^{2}I_{L\times
L}),$ (4)

where $F^{L}$ is a square image of size $L\times L$ enclosing the disk-shaped
image $F$ of radius $n$, the operator $Z$ zero-pads $L$ entries to the right
and to the bottom of a rotated copy of $F$, and $T_{\vec{\ell}_{m}}$
circularly shifts the zero-padded image by
$\vec{\ell}_{m}=({\ell_{m}}_{x},{\ell_{m}}_{y})\in\mathbb{L}:=\\{0,1,\ldots,2L-1\\}^{2}$
positions, that is,

$(T_{\vec{\ell}_{m}}ZF^{{{L}}}_{\phi_{m}})\left[i,j\right]=\\\
(ZF^{{{L}}}_{\phi_{m}})\left[(i+{\ell_{m}}_{x})\bmod
2L,(j+{\ell_{m}}_{y})\bmod 2L\right].$ (5)

The operator $C$ then crops the image to size $L\times L$, i.e.,
$CT_{\vec{\ell}_{m}}ZF^{{{L}}}_{\phi_{m}}=T_{\vec{\ell}_{m}}ZF^{{{L}}}_{\phi_{m}}[0:L-1,0:L-1]$,
and the result is further corrupted by additive white Gaussian noise. In
addition, since the EM algorithm assigns probabilities to rotations (in the
expectation step), we need to discretize the search space of rotations:

$\phi_{m}\in\Phi:=\left\\{k\frac{2\pi}{K}\right\\},\quad k=0,\ldots,K-1,$ (6)

where $K$ is a parameter chosen by the user. Higher $K$ provides higher
accuracy at the cost of running time (see Figure 4). The rotations,
translations, and patches are referred to, in EM terminology, as the complete
data. We assume that in each patch
$p(\vec{\ell}_{m},\phi_{m}|\alpha_{k})=p(\vec{\ell}_{m})p(\phi_{m})$, namely,
$\vec{\ell}_{m}$ and $\phi_{m}$ are independent of $\alpha_{k}$ and of each
other. Since the rotations are drawn from a uniform distribution over the set
$\Phi$ from (6), we can write
$p(\vec{\ell}_{m},\phi_{m}|\alpha_{k})=\rho[\vec{\ell}_{m}]/K$, where
$\rho[\vec{\ell}]$ is the distribution of 2-D translations in the patch (which
should be estimated simultaneously with $\alpha$).

In the E-step, our algorithm calculates the expected log-likelihood function
of the model, given explicitly by

$Q(\alpha,\rho|\alpha_{k},\rho_{k})=\mathbb{E}\left[\log\mathfrak{L}|M_{0},\ldots,M_{N_{d}-1};\alpha_{k},\rho_{k}\right],$
(7)

where $\mathfrak{L}$ is the likelihood function, defined as

$\mathfrak{L}:=p(M_{0},\ldots,M_{N_{d}-1},\vec{\ell}_{0},\ldots,\vec{\ell}_{N_{d}-1},\phi_{0},\ldots,\phi_{N_{d}-1};\alpha,\rho)\\\
\approx\prod_{m=0}^{N_{d}-1}p(M_{m},\vec{\ell}_{m},\phi_{m};\alpha,\rho),$ (8)

where we neglect statistical dependencies between patches.

Bayes’ rule dictates

$p(\vec{\ell}_{m},\phi_{m}|M_{m},\alpha_{k})\\\
=\frac{p(M_{m}|\vec{\ell}_{m},\phi_{m},\alpha_{k})p(\vec{\ell}_{m},\phi_{m}|\alpha_{k})}{\sum_{\vec{\ell^{\prime}}\in\mathbb{L}}\sum_{\phi^{\prime}\in\Phi}p(M_{m}|\vec{\ell^{\prime}},\phi^{\prime},\alpha_{k})p(\vec{\ell^{\prime}},\phi^{\prime}|\alpha_{k})},$
(9)

which is the normalized likelihood function

$p(M_{m}|\vec{\ell}_{m},\phi_{m},\alpha)\propto\exp\left(-\frac{\|M_{m}-CT_{\vec{{{\ell}}_{m}}}ZF_{\phi_{m}}\|_{\text{F}}^{2}}{2\sigma^{2}}\right),$
(10)

with the normalization
$\sum_{\vec{\ell}\in\mathbb{L}}\sum_{\phi\in\Phi}p(M_{m}|\vec{\ell}_{m},\phi_{m},\alpha)=1$,
weighted by the prior distribution
$p(\vec{\ell}_{m},\phi_{m}|\alpha_{k})=\rho[\vec{\ell}_{m}]/K$.

Utilizing the approximation in (8), we can write the expected log-likelihood
function, up to a constant, as:

$Q(\alpha,\rho|\alpha_{k},\rho_{k})=\sum_{m=0}^{N_{d}-1}\sum_{\vec{\ell}\in\mathbb{L}}\sum_{\phi\in\Phi}p(M_{m}|\vec{\ell},\phi,\alpha_{k})\rho_{k}[\vec{\ell}]\\\
\times\left(\log p(M_{m}|\vec{\ell},\phi,\alpha)+\log\rho[\vec{\ell}]\right).$
(11)

The M-step updates the image estimate $\alpha$ and $\rho$ by maximizing
$Q(\alpha,\rho|\alpha_{k},\rho_{k})$ under the constraint that $\rho$ lies on
the simplex $\Delta_{4L^{2}}$:

$\arg\max_{\alpha,\rho}Q(\alpha,\rho|\alpha_{k},\rho_{k})\quad\text{s.t.}\quad\rho\in\Delta_{4L^{2}}.$
(12)

The constrained maximization of (12) can be achieved with the unconstrained
maximization of the Lagrangian

$\mathcal{L}(\alpha,\rho,\eta)=Q(\alpha,\rho|\alpha_{k},\rho_{k})+\eta\left(1-\sum_{\vec{\ell}\in\mathbb{L}}\rho[\vec{\ell}]\right),$
(13)

where $\eta$ is the Lagrange multiplier. As we will see next, the constraint
of (12) is automatically satisfied at the maximum of the Lagrangian.

Since $Q(\alpha,\rho|\alpha_{k},\rho_{k})$ is additively separable for
$\alpha$ and $\rho$, we maximize $\mathcal{L}(\alpha,\rho,\eta)$ with respect
to $\alpha$ and $\rho$ separately. At the maximum of
$\mathcal{L}(\alpha,\rho,\eta)$, we have

$0=\frac{\partial\mathcal{L}}{\partial{(\alpha)_{\nu,q}}}=\sum_{m=0}^{N_{d}-1}\sum_{\vec{\ell}\in\mathbb{L}}\sum_{\phi\in\Phi}p(M_{m}|\vec{\ell},\phi,\alpha_{k})\rho_{k}[\vec{\ell}]\\\
\times\frac{\partial\log
p(M_{m}|\vec{\ell},\phi,\alpha)}{\partial{(\alpha)_{\nu,q}}},$ (14)

resulting in a set of linear equations which is solved to update $\alpha$. In
order to update $\rho$, we maximize $\mathcal{L}(\alpha,\rho,\eta)$ with
respect to $\rho$:

$0=\frac{\partial\mathcal{L}}{\mathcal{\rho}[\vec{\ell}]}=\sum_{m=0}^{N_{d}-1}\sum_{\phi\in\Phi}p(M_{m}|\vec{\ell},\phi,\alpha_{k})\rho_{k}[\vec{\ell}]\frac{1}{\rho[\vec{\ell}]}-\eta,$
(15)

for $\vec{\ell}\in\mathbb{L}$. We thus obtain the update rule for $\rho$ as

$\rho[\vec{\ell}]=\frac{1}{\eta}\sum_{m=0}^{N_{d}-1}\sum_{\phi\in\Phi}p(M_{m}|\vec{\ell},\phi,\alpha_{k})\rho_{k}[\vec{\ell}],$
(16)

and $\eta=N_{d}$ from the normalization
$\sum_{\vec{\ell}\in\mathbb{L}}\rho[\vec{\ell}]=1$. The approximate EM
algorithm is summarized in Algorithm 1.

Input: measurement $M$; noise variance $\sigma^{2}$; initial guesses
$\alpha_{0}$ and $\rho_{0}$; parameter $K$ (6); stopping parameter $\epsilon$

Output: an estimate of $\alpha$ and $\rho$

1 set $k\rightarrow 0$;

2 calculate $Q_{0}$ according to (11) and set $Q_{-1}\rightarrow-\infty$;

3 while  _$Q_{k}-Q_{k-1} >\epsilon$_ do

4 calculate $p(M_{m}|\vec{\ell},\phi,\alpha_{k})$ according to (10);

5 update $\alpha_{k+1}$ by solving (14);

6 update $\rho_{k+1}$ according to (16);

7 calculate $Q_{k+1}$ according to (11);

8 set $k\rightarrow k+1$;

Algorithm 1 Approximate EM for 2-D MTD Figure 2: The mean estimation error of
recovering the target image $F$, as a function of the SNR, by approximate EM
(Algorithm 1) and autocorrelation analysis. Evidently, the approximate EM
algorithm outperforms autocorrelation analysis.

##  III Numerical experiments

In this section, we present numerical results for the approximate EM described
in Section II. As a baseline, we compare the results against autocorrelation
analysis with the first three autocorrelations based on the framework (and
code) of [3, 4, 5]. To take the in-plane rotation symmetry into account, we
measure the estimation error by
$\min_{\phi\in[0,2\pi)}\frac{\|\alpha^{*}-\alpha_{\phi}\|_{2}}{\|\alpha^{*}\|_{2}},$
where $\alpha^{*}$ is the true vector of expansion coefficients, and
$\alpha_{\phi}$ is the vector of coefficients of the estimated image, rotated
by $\phi$. In all experiments, the measurements were generated according to
(1) with density $\gamma:=p\frac{\pi n^{2}}{N^{2}}=0.04$. The rotations were
drawn from a uniform distribution on $[0,2\pi)$, while the search space of the
EM was discretized with a parameter $K$ according to (6). The target images
are of diameter $L=5\text{ pixels}$. Each entry of the target images was drawn
i.i.d. from a uniform distribution on $[0,1]$. Then, each image was normalized
such that $\|F\|_{\text{F}}=10$, and expanded using its first $10$
coefficients as in (2). The initializations of the EM and autocorrelation
analysis iterations were drawn from the same distribution as the ground truth
images, and $\gamma_{\text{init}}=0.03$. If the algorithms were initialized
from several random points, we calculated the error of the image estimate
whose likelihood is maximal (for approximate EM), or whose objective function
is minimal (for autocorrelation analysis). Figures 2, 3 and 4 present the mean
error over $40$ trials. The code to reproduce all experiments is publicly
available at https://github.com/krshay/MTD-2D-EM.

###  III-A Recovery error as a function of the SNR

Figure 2 presents recovery error as a function of the SNR. The measurements
are of size $N^{2}=2500^{2}\text{ pixels}$, and we use $K=8$ possible
rotations (6), and $5$ random initial guesses for $\alpha$. To save
computation time (see Section IV), we initialized the approximate EM algorithm
using the estimate achieved by autocorrelation analysis. We achieve a
significant improvement in recovery accuracy using approximate EM, even though
the search space of rotations is coarsely sampled.

In an additional numerical experiment with $\text{SNR}=2$ (low SNR, see Figure
1), $N^{2}=10000^{2}\text{ pixels}$, and $K=16$, the relative error of
approximate EM was 0.017, whereas the error of autocorrelation analysis was
$0.073$.

Figure 3: The mean estimation error of recovering the target image $F$ (left)
and running time (right), as a function of the measurement size $N^{2}$ by
approximate EM and autocorrelation analysis. For the estimation error, the
black dashed lines illustrate a slope of $-1/2$, as predicted by the law of
large numbers. For the computation time, it illustrates a slope of $1$,
implying a linear increase in computation time, since the number of patches
$N_{d}=N^{2}/L^{2}$ grows linearly in $N^{2}$.

Figure 4: The mean estimation error of recovering the target image $F$ (left)
and corresponding running time (right), as a function of the size of the
search space of rotations $K$, by approximate EM. The recovery error and
running time using autocorrelation analysis are marked by dashed red lines.
For the running time, the black dashed line illustrates a slope of $1$, which
implies a linear increase in computation time, as the number of computations
per patch depends linearly on $K$.

###  III-B Recovery error as a function of the measurement size

Figure 3 presents recovery error and running time as a function of the
measurement size $N^{2}$, with $\text{SNR}=5$, $K=16$, and $5$ random initial
guesses for $\alpha$. Using approximate EM, the error decays as
$1/\sqrt{N^{2}}$. The same trend is visible also for autocorrelation analysis
for sufficiently large measurements. We achieve a significant improvement in
recovery accuracy using approximate EM. However, the computation time of our
method is greater, and grows linearly with the measurement size. In
particular, autocorrelation analysis is faster since it requires a single pass
over the measurement, while the computational complexity of matching an image
to the observed autocorrelations scale with $L^{2}\ll N^{2}$.

###  III-C Recovery error as a function of discretization of rotations

Figure 4 presents recovery error and running time as a function of $K$, the
size of the search space of rotations, for measurements with $\text{SNR}=5$,
$N^{2}=1500^{2}\text{ pixels}$, and one initial guess for $\alpha$.
Remarkably, even when the EM searches over only $4$ rotations (recall that the
rotations are drawn from a continuous distribution), the obtained estimation
error is similar to the estimation error of autocorrelation analysis, which
takes all possible (infinitely many) rotations into account. As expected, the
computation time grows linearly with the parameter $K$. This implies that one
can save running time by coarsely sampling the search space of rotations
without severely degrading the estimation quality.

##  IV Conclusion

This paper is motivated by the effort of reconstructing small 3-D molecular
structures using cryo-EM, below the current detection limit [7]. The main
contribution of this paper is in introducing an approximate EM scheme for the
2-D MTD problem, and comparing it numerically to autocorrelation analysis. The
numerical experiments show an improvement in estimation accuracy, but at the
cost of computational time. As Figure 4 shows, the parameter $K$ provides an
accuracy-running time trade-off. A possible improvement is to increase the
resolution of the search space as the iterations progress [2], or to design a
branch-and-bound algorithm that accelerates EM by quickly and inexpensively
ruling out large regions of the search space that have very low probability to
contain the optimum of the objective function [16].

Our ultimate goal is to develop an approximate EM scheme for recovering small
molecular structures using cryo-EM [7]. In order to achieve a computationally
efficient algorithm for the 3-D case of cryo-EM, parallel processing and
randomized algorithms, such as stochastic or online EM [23, 24, 25, 26, 27],
must be utilized. Further accelerations can be achieved by applying the EM
algorithm on a lower dimensional representation of the data [28], and
initializing the EM iterations with efficient computational techniques, such
as autocorrelation analysis or stochastic gradient descent [16]. Moreover,
adding a prior on the target image is expected to improve robustness at the
cost of possible model bias. Another research direction is replacing EM by
more intricate techniques that aim to approximate the posterior distribution,
such as variational inference [29] or variational auto-encoders [30].

## References

  * [1] T. Bendory, N. Boumal, W. Leeb, E. Levin, and A. Singer, “Multi-target detection with application to cryo-electron microscopy,” _Inverse Problems_ , vol. 35, no. 10, p. 104003, 2019.
  * [2] T.-Y. Lan, T. Bendory, N. Boumal, and A. Singer, “Multi-target detection with an arbitrary spacing distribution,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 1589–1601, 2020.
  * [3] N. F. Marshall, T.-Y. Lan, T. Bendory, and A. Singer, “Image recovery from rotational and translational invariants,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 5780–5784.
  * [4] T. Bendory, T.-Y. Lan, N. F. Marshall, I. Rukshin, and A. Singer, “Multi-target detection with rotations,” _arXiv preprint arXiv:2101.07709_ , 2021.
  * [5] S. Kreymer and T. Bendory, “Two-dimensional multi-target detection: An autocorrelation analysis approach,” _IEEE Transactions on Signal Processing_ , vol. 70, pp. 835–849, 2022.
  * [6] Y. Shalit, R. Weber, A. Abas, S. Kreymer, and T. Bendory, “Generalized autocorrelation analysis for multi-target detection,” _arXiv preprint arXiv:2109.11813_ , 2021.
  * [7] T. Bendory, N. Boumal, W. Leeb, E. Levin, and A. Singer, “Toward single particle reconstruction without particle picking: breaking the detection limit,” _arXiv preprint arXiv:1810.00226_ , 2018.
  * [8] Z. Zhao and A. Singer, “Fourier–Bessel rotational invariant eigenimages,” _JOSA A_ , vol. 30, no. 5, pp. 871–877, 2013.
  * [9] R. Henderson, “The potential and limitations of neutrons, electrons and X-rays for atomic resolution microscopy of unstained biological molecules,” _Quarterly Reviews of Biophysics_ , vol. 28, no. 2, pp. 171–193, 1995.
  * [10] E. Nogales, “The development of cryo-EM into a mainstream structural biology technique,” _Nature Methods_ , vol. 13, no. 1, pp. 24–27, 2016.
  * [11] X.-C. Bai, G. McMullan, and S. H. Scheres, “How cryo-EM is revolutionizing structural biology,” _Trends in Biochemical Sciences_ , vol. 40, no. 1, pp. 49–57, 2015.
  * [12] J. Frank, _Three-dimensional electron microscopy of macromolecular assemblies: visualization of biological molecules in their native state_. Oxford University Press, 2006.
  * [13] T. Bendory, A. Bartesaghi, and A. Singer, “Single-particle cryo-electron microscopy: Mathematical theory, computational challenges, and opportunities,” _IEEE Signal Processing Magazine_ , vol. 37, no. 2, pp. 58–76, 2020.
  * [14] A. Singer and F. J. Sigworth, “Computational methods for single-particle electron cryomicroscopy,” _Annual Review of Biomedical Data Science_ , vol. 3, pp. 163–190, 2020.
  * [15] S. H. Scheres, “RELION: implementation of a Bayesian approach to cryo-EM structure determination,” _Journal of Structural Biology_ , vol. 180, no. 3, pp. 519–530, 2012.
  * [16] A. Punjani, J. L. Rubinstein, D. J. Fleet, and M. A. Brubaker, “cryoSPARC: algorithms for rapid unsupervised cryo-EM structure determination,” _Nature Methods_ , vol. 14, no. 3, pp. 290–296, 2017.
  * [17] C. Aguerrebere, M. Delbracio, A. Bartesaghi, and G. Sapiro, “Fundamental limits in multi-image alignment,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 21, pp. 5707–5722, 2016.
  * [18] E. D’Imprima and W. Kühlbrandt, “Current limitations to high-resolution structure determination by single-particle cryoEM,” _Quarterly Reviews of Biophysics_ , vol. 54, 2021.
  * [19] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm,” _Journal of the Royal Statistical Society: Series B (Methodological)_ , vol. 39, no. 1, pp. 1–22, 1977\.
  * [20] T. Bendory, N. Boumal, C. Ma, Z. Zhao, and A. Singer, “Bispectrum inversion with application to multireference alignment,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 4, pp. 1037–1050, 2017.
  * [21] E. Abbe, T. Bendory, W. Leeb, J. M. Pereira, N. Sharon, and A. Singer, “Multireference alignment is easier with an aperiodic translation distribution,” _IEEE Transactions on Information Theory_ , vol. 65, no. 6, pp. 3565–3584, 2018.
  * [22] A. Katsevich and A. Bandeira, “Likelihood maximization and moment matching in low SNR gaussian mixture models,” _arXiv preprint arXiv:2006.15202_ , 2020\.
  * [23] S. F. Nielsen, “The stochastic EM algorithm: estimation and asymptotic results,” _Bernoulli_ , pp. 457–489, 2000.
  * [24] J. Chen, J. Zhu, Y. W. Teh, and T. Zhang, “Stochastic expectation maximization with variance reduction,” in _NeurIPS_ , 2018, pp. 7978–7988.
  * [25] P. Liang and D. Klein, “Online EM for unsupervised models,” in _Proceedings of Human Language Technologies: The 2009 Annual Conference of the North American Chapter of the Association for Computational Linguistics_ , 2009, pp. 611–619.
  * [26] O. Cappé and E. Moulines, “On-line expectation–maximization algorithm for latent data models,” _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , vol. 71, no. 3, pp. 593–613, 2009.
  * [27] O. Cappé, “Online EM algorithm for hidden Markov models,” _Journal of Computational and Graphical Statistics_ , vol. 20, no. 3, pp. 728–749, 2011.
  * [28] N. C. Dvornek, F. J. Sigworth, and H. D. Tagare, “SubspaceEM: A fast maximum-a-posteriori algorithm for cryo-EM single particle reconstruction,” _Journal of Structural Biology_ , vol. 190, no. 2, pp. 200–214, 2015.
  * [29] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe, “Variational inference: A review for statisticians,” _Journal of the American Statistical Association_ , vol. 112, no. 518, pp. 859–877, 2017.
  * [30] D. Rosenbaum, M. Garnelo, M. Zielinski, C. Beattie, E. Clancy, A. Huber, P. Kohli, A. W. Senior, J. Jumper, C. Doersch _et al._ , “Inferring a continuous distribution of atom coordinates from cryo-EM images using VAEs,” _arXiv preprint arXiv:2106.14108_ , 2021.