Bayesian Random Tomography of Particle Systems

Random tomography is a common problem in imaging science and refers to the task of reconstructing a three-dimensional volume from two-dimensional projection images acquired in unknown random directions. We present a Bayesian approach to random tomography. At the center of our approach is a meshless representation of the unknown volume as a mixture of spherical Gaussians. Each Gaussian can be interpreted as a particle such that the unknown volume is represented by a particle cloud. The particle representation allows us to speed up the computation of projection images and to represent a large variety of structures accurately and efficiently. We develop Markov chain Monte Carlo algorithms to infer the particle positions as well as the unknown orientations. Posterior sampling is challenging due to the high dimensionality and multimodality of the posterior distribution. We tackle these challenges by using Hamiltonian Monte Carlo and a global rotational sampling strategy. We test the approach on various simulated and real datasets.

for 2D point clouds where ξ n = (t n , s n , σ n ).
The Jeffreys' priors for the precision parameters are τ n ∼ 1/τ n where in case of the point cloud likelihood (Eq. S2) τ n = 1/σ 2 n . The conditional posteriors of these parameters are Gamma distributions: where χ 2 n = Mn m=1 g nm − α n − γ n k φ 2 u nm ; s n (P R n x k + t n ), σ 2 2 Efficient random number generators for Gamma-variates are offered by libraries such as numpy.random from the Python array processing library NumPy. The conditional posterior of the offset and scale is a two-dimensional Gaussian; to update these parameters, we generate a sample from the bivariate Normal distribution by calling numpy.random.multivariate normal. The shifts t n and magnification factors s n can be sampled with the Metropolis-Hastings algorithm.

FITTING POINT CLOUDS TO IMAGES WITH EXPECTATION MAXIMIZATION
Given a collection of pixels u i with associated intensities g i (i = 1, . . . , I), we aim to represent this information by a 2D point cloud {y m ; m = 1, . . . , M }. The model to relate the image intensities and pixels to a point cloud is a mixture of Gaussians (Eq. 12): where α is a suitable background parameter, γ is a scaling factor and σ the width of the Gaussian components. Using the approach described in subsection 3.2, we correct for the background by masking out the intensities that are greater than a suitable threshold and subtracting the threshold from the intensity such that α = 0. Without loss of generality we can set γ = 1. So we are trying to fit y m and σ such that the shifted intensities g i at pixels u i are reproduced as closely as possible with the Gaussian mixture model: This can be achieved with an Expectation Maximization (EM) algorithm that cycles over the following updates: • Soft assignment: Each pixel u i with associated intensity g i is assigned to a 2D point y m with probability p im : • Particle positions y m are the centers of mass of the assigned pixels weighted by the (positive) image intensity g i and the assignment probabilities p im computed in the previous step: where the denominator is the total image intensity represented by the m-th particle.
• The width of the Gaussians is estimated by For the entire series of 400 class averages provided by SIMPLE, we fitted 1000 particles to each image. The reasoning for choosing the same number of particles is that our model predicts the same total intensity which is identical to the number of points that represent the projection images. Figure S1 shows five representative class averages of the 80S ribosome and the their point cloud representation. 10 1 10 2 number of atoms per particle L/K 10 1 2 × 10 0 3 × 10 0 4 × 10 0 6 × 10 0 particle radius R [Å] Figure S2. Relation between the average number of atoms per particle (L/K) and the particle radius.

CHOICE OF PARAMETERS IN BOLTZMANN PRIOR
To enforce a reasonable packing of particle structures, we use a Boltzmann distribution (Eq. 17) with a repulsive pairwise distance potential (Eq. 18). The only parameters that need to be set are the particle radius R and the inverse temperature β of the Boltzmann ensemble.

Particle radius
Biomolecules pack in a way that is reminiscent of fluids (Liang and Dill, 2001). To design a prior distribution that favors particle configurations with similar packing characteristics, we analyzed a number of biomolecular structures at different degrees of coarse graining. The coarse-grained models were derived from atomic structures by running the DP-means algorithm (Kulis and Jordan, 2012) for different distance cutoffs λ = 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20Å, rather than using a computationally more demanding coarse graining approach proposed by us (Chen and Habeck, 2017).
The DP-means algorithm is a variant of the K-means clustering algorithm. DP-means assigns points to cluster centers, if the Euclidean distance of a point from a center does not exceed a distance cutoff. If no cluster exists that is close enough to a given point, a new cluster is created. This process is repeated multiple times until the iterations converge. In contrast to K-means, the number of clusters varies and depends on the cutoff value λ, which is the only input parameter: smaller cutoffs result in a larger number of clusters. The clusters centers form the particle positions; the distance cutoff λ corresponds to the particle diameter. We use the implementation of DP-means available at https://github.com/michaelhabeck/ DP-means.
The coarse-grained structures obtained with DP-means allow us to establish a relation between the number of particles K and the distance cutoff λ. In our implementation of Bayesian particle-based tomography, we work with a fixed number of particles and want to choose the particle radius R such that volume-exclusion observed in known biomolecular structures is enforced. Since different biomolecular structures pack similarly, we expect that the average number of atoms per particle can be related to the particle radius  Figure S3. Configurational temperature as a function of particle radius.
independent of the specific structure. If L denotes the total number of atoms in a biomolecular structure, then L/K is the average number of atoms per particle. We use half the DP-means cutoff as a proxy for the particle radius R = λ/2. Figure S2 shows that there is indeed a high correlation between L/K and R that can be captured with a linear relation between the logarithms of both quantities: log R ≈ −0.087 + 0.423 log(L/K) Mapping this back, we can predict a particle radius from the number of atoms L and the chosen number of particles K. For a given biomolecular system, we compute L based on the amino acid sequence.

Inverse temperature
The inverse temperature β is estimated from coarse grained structures using the configurational temperature formalism (Mechelke and Habeck, 2013). For each coarse-grained model that was produced to estimate the particle radii, we computed the configurational temperature. Figure S3 shows the relation between particle diameter and the configurational temperature, which can be fitted with a straight line in log space. The inverse temperature corresponds to the slope of the line and is set to β = 175.

SAMPLING PARTICLE POSITIONS AND PRECISIONS WITH FIXED ROTATIONS
Particle positions are sampled with HMC. To test the performance of HMC, we fitted coarse-grained models of GroEL/ES against 35 simulated projection images and 2D point clouds. For both types of data, HMC generates high-quality reconstructions.    Table S1. Resolution estimates and RMSD values for particle-based models of GroEL/ES using point clouds and projection images as input data.

Resolution assessment
In addition to the tests based on 2D point clouds derived from the projection images, we also used the first 50 class averages themselves as an input. We assessed the accuracy of the particle models inferred from both types of input data by computing the FSC correlating the high-resolution structure EMD-2660 with the initial model generated by our method. Figure S6 shows the FSCs and Table S2 lists the resolutions derived from these curves.  12.5 16.9 4000 9.3 12.5 8000 7.0 11.9 12000 5.9 10.6 Table S2. Resolution estimates for particle-based models of the 80S ribosome obtained with random tomography. The FSCs are computed by comparing the high-resolution reconstruction (EMD-2660) with the density map generated from the last 100 sampled particle configurations.

Computation times
Computation times for 100 steps of Gibbs sampling using 50 class averages as input. Tests were run on an intel i5 processor (1.6 GHz oct-core using a single thread only).  Table S3. Computation times for a 3D reconstruction from 50 class averages of the 80S ribosome using 100 Gibbs sampling steps.

Uncertainty Quantification
Bayesian methods and posterior sampling allow for uncertainty quantification. Since particle models can differ by a rigid transformation and a permutation of particle indices, we first need to superimpose the configurations generated by posterior sampling and establish a correspondence between particle positions across different samples. We do this by using the ICP (iterative closest point) method followed by linear assignment. After these steps, we can compute the standard deviation of each particle position which provides a structural error bar. Figure S7 shows a particle model where the error bar is encoded in the bead radius.   Figure S9. FSC curves comparing the high-resolution reconstruction EMD-5995 with Bayesian models using 100 and 2000 particles, respectively, as well as an initial structure obtained with RELION.

IMPACT OF BOLTZMANN PRIOR
To assess the impact of the Boltzmann prior (Eq. 17), we ran simulations with the beta-galactosidase images where the Boltzmann prior was switched on and off (β = 0). with prior without prior Figure S11. Impact of Boltzmann prior on particle packing. Shown are the radial distribution functions computed from sampled particle models of beta-galactosidase with (light blue line) and without Boltzmann prior (dark blue line).