Dense Super-Resolution Imaging of Molecular Orientation via Joint Sparse Basis Deconvolution and Spatial Pooling

TL;DR

This paper introduces a joint sparse deconvolution method for super-resolution microscopy that accurately estimates molecular position and orientation, revealing nanoscale DNA structure despite PSF overlap.

Contribution

It presents a novel algorithm combining basis decomposition and spatial pooling for robust multi-parameter molecular imaging in super-resolution microscopy.

Findings

Successfully imaged DNA with nanoscale spatial precision.

Revealed molecular orientational order and disorder.

Demonstrated robustness against PSF overlap and model mismatches.

Abstract

In single-molecule super-resolution microscopy, engineered point-spread functions (PSFs) are designed to efficiently encode new molecular properties, such as 3D orientation, into complex spatial features captured by a camera. To fully benefit from their optimality, algorithms must estimate multi-dimensional parameters such as molecular position and orientation in the presence of PSF overlap and model-experiment mismatches. Here, we present a novel joint sparse deconvolution algorithm based on the decomposition of fluorescence images into six basis images that characterize molecular orientation. The proposed algorithm exploits a group-sparsity structure across these basis images and applies a pooling strategy on corresponding spatial features for robust simultaneous estimates of the number, brightness, 2D position, and 3D orientation of fluorescent molecules. We demonstrate this method by imaging DNA transiently labeled with the intercalating dye YOYO-1. Imaging the position and orientation of each molecule reveals orientational order and disorder within DNA with nanoscale spatial precision.

Equations

$$I(\mathbf{u})$$

$$I(\mathbf{u})=$$

$$\displaystyle{h^{\prime}_{y}}^{j}(\bm{u}-\bm{d_{i}})\zeta_{y,i}^{j}\Big{)},$$

$$\hat{\mathbf{F}}=\underset{\mathbf{F}}{\mathrm{argmin}}\{\mathcal{L}(\mathbf{F})+\lambda \mathrm{R}(\mathbf{F})+{\mathcal{I}}_{\mathcal{C}}(\mathbf{F})\},$$

$$\underset{\mathbf{F}\in \mathcal{C}\cap \text{Supp}(\hat{\mathbf{F}})}{\min }\mathcal{L}(\mathbf{F}),$$

Full text

**Index Terms— ** Sparse deconvolution, group sparsity, single-molecule orientation, DNA intercalators

1 Introduction

Single-molecule localization microscopy (SMLM) relies on accurately and precisely estimating the position of many single emitters repeatedly blinking over time [1, 2, 3]. Besides their position, fluorescent molecules interacting with their nano-environment also report information such as orientation and emission spectra that are hidden from standard, diffraction-limited optical imaging systems. Augmenting the PSF, or the impulse response of the microscope, to uncover such fine details remains a vibrant research topic [4, 5]. In addition, biophysical measurements utilizing SMLM are more accurate when a physically-realistic image-formation model and recovery problem are utilized that include parameters such as molecular orientation and rotational mobility.

However, estimating multidimensional parameters such as molecular position, orientation, and brightness from noisy camera images poses a formidable challenge. For example, engineered PSFs have much larger footprints compared to the standard PSF (Fig. 1(b)), thereby causing frequent PSF overlaps [6] and lower pixel-wise signal-to-noise ratios, which can substantially reduce the detection rate of standard spot-detection methods. Further, a target structure must be densely sampled by blinking emitters in order to be resolved, and therefore, it is vital to develop an algorithmic framework to address the aforementioned challenges. Joint estimation of the molecular position and orientation for the standard PSF has been considered previously [7, 8, 9]. These methods, however, cannot be adapted to engineered PSFs for measuring jointly the position, 3D orientation, and rotational dynamics of fluorescent molecules. In addition, the problem of PSF overlap in the presence of multi-channel measurement errors has not been considered before. In this work, we tackle these challenges by 1) constructing a forward model that decomposes molecular parameters (i.e., brightness, position, and orientation) into six basis images corresponding to the six second moments of orientation dynamics [5]; 2) building a novel sparse basis deconvolution algorithm that exploits a group-sparsity norm to jointly recover all molecular parameters; and 3) applying spatial pooling across the six basis images to avoid false localizations due to experimental mismatches. We applied the proposed method to resolve DNA conformation via super-resolution imaging of DNA intercalating probes.

2 Image-formation model

A dipole, such as a single fluorescent molecule, is characterized by its position $\bm{r}=[x,y]^{T}$ and an orientation vector $\bm{\mu}=[\mu_{x},\mu_{y},\mu_{z}]^{T}=[\sin(\theta)\cos(\phi),\sin(\theta)\allowbreak\sin(\phi),\cos(\theta)]^{T}\in\mathbb{R}^{3}$, where $\theta$ and $\phi$ are the polar and the azimuthal angles in spherical coordinates, respectively (Fig. 1(a)). The image of a dipole, denoted by $I$, can be described in terms of an orientational second-moment vector $\bm{\mathcal{M}}=[\langle\mu_{x}^{2}\rangle,\langle\mu_{y}^{2}\rangle,\langle\mu_{z}^{2}\rangle,\langle\mu_{x}\mu_{y}\rangle,\allowbreak\langle\mu_{x}\mu_{z}\rangle,\allowbreak\langle\mu_{y}\mu_{z}\rangle]^{T}$, in which $\langle\cdot\rangle$ represents the temporal average over a camera frame or equivalently an ensemble average over the orientation domain [8, 10]. More precisely, I=s\Big{(}\sum_{j=1}^{6}\mathcal{M}^{j}\mathcal{B}^{j}\Big{)}. Here, $s$ denotes the total number of photons emitted by the dipole, and the $\mathcal{B}^{j}$’s represent the so-called basis images corresponding to a dipole exhibiting each orientational second-moment component. Note that here $I$ corresponds to two orthogonally-polarized images (i.e., x and y), concatenated together (Figs. 1(b), 2(c,d), 3(a inset)).

Consider the object domain $\Omega=\{\omega_{i}=(s_{i},\bm{\mathcal{M}}_{i},\bm{r}_{i})\ |\ s_{i}\allowbreak\geq 0,\|\bm{r}_{i}\|_{\infty}<r_{\text{max}},i\in\{1,2,\ldots,P\}\}$ as a collection of $P$ dipole-like molecules located within a square region of interest of length 2$r_{\text{max}}$. We represent the PSF of the microscope as $h(\bm{u};\omega):\Omega\rightarrow\mathcal{I}\times\mathbb{R}$ for which $\bm{u}$ denotes the coordinates in image space $\mathcal{I}\subset\mathbb{R}^{2}$. The noiseless image formed on the camera is the result of a weighted-convolution operation $I(\bm{u})=\sum_{i=1}^{P}s_{i}h(\bm{u};\omega_{i})*\delta(\bm{r}-\bm{r}_{i})$ ($*$ denotes the (linear) convolution operator and $\delta(\cdot)$ represents the Dirac delta function). We now exploit the fact that $h$ can be represented as the sum of six basis images:

[TABLE]

where $h^{j}(\cdot)$ corresponds to the basis image $\mathcal{B}^{j}$, and $\mathcal{M}_{i}^{j}$ represents the $j_{\text{th}}$ component of the second-order orientational moment of the $i_{\text{th}}$ emitter. As $s$ plays a scaling role here, we absorb it into $\bm{\mathcal{M}}$ and denote it by $\bm{\eta}=s\bm{\mathcal{M}}$. In the next section, we extend our model in Eq. (1) to explicitly take into account the continuous position of molecules, which proves to be crucial for the robustness of the recovery algorithm.

Let $D$ be a set of $N$ discrete grid points, where the distance between two adjacent points is given by $2\rho$. A position vector $\bm{r}$ can be uniquely mapped into $\bm{r}=\bm{d}+\bm{\delta}$ for some $\bm{d}\in D$ and $\bm{\delta}=[\delta_{x},\delta_{y}]^{T}$ with $\delta_{x},\delta_{y}\in[-\rho,\rho)$. Therefore, we can rewrite Eq. (1) as $I(\bm{u})=\sum_{i=1}^{N}\sum_{j=1}^{6}\eta^{j}_{i}h^{j}(\bm{u}-\bm{d_{i}-\delta_{i}})$. Note that $\bm{\eta}_{i}=0$ if no molecule can be mapped to a point in the neighbourhood of $\bm{d}_{i}$. For a fixed $i\in\{1,2,\ldots,N\}$ and $j\in\{1,2,\ldots,6\}$, we approximate $h$ via first-order Taylor expansion $h^{j}(\bm{u-d_{i}-\delta_{i}})\approx h^{j}(\bm{u}-\bm{d_{i}})-h^{\prime j}(\bm{u-d_{i}})^{T}\bm{\delta}_{i}$, in which $h^{\prime}(\cdot)$ denotes the derivative of $h$. Within the first-order approximation, we have

[TABLE]

where $\bm{\zeta}_{i}^{j}=\eta^{j}_{i}\bm{\delta}_{i}=[\zeta_{x,i}^{j},\zeta_{y,i}^{j}]^{T}$; ${h^{\prime}_{x}}^{j}$ and ${h^{\prime}_{y}}^{j}$ represent derivatives of $h$ along the $x$-axis and $y$-axis, respectively. Remarkably, the image model in Eq. (2) can be interpreted as the sum of a few PSFs weighted by object-domain parameters $\Omega$.

We also take into account camera pixelation by integrating the image $I$ over $m$ pixels, which effectively results in matrices $\bm{\Phi}^{j},\bm{Gx}^{j}$ and $\bm{Gy}^{j}$ for each basis with index $j$. Note that $\bm{\Phi}$ corresponds to $h$ in Eq. (2) whereas $\bm{Gx}$ and $\bm{Gy}$ correspond to $h^{\prime}_{x}$ and $h^{\prime}_{y}$, respectively. The final imaging model can be compactly represented as $\mathcal{A}\bm{F}=\sum_{j=1}^{6}\bm{A}^{j}\bm{f}^{j}$ with $\bm{A}^{j}=[\bm{\Phi}^{jT},\bm{Gx}^{jT},\bm{Gy}^{jT}]^{T}$ and $\bm{f}^{j}=[\bm{\eta}^{jT},\bm{\zeta}_{x}^{jT},\allowbreak\bm{\zeta}_{y}^{jT}]^{T}$. We further model the photon count in each pixel $i$ as an independent Poisson distribution g_{i}\sim\text{Pois}\Big{(}(\mathcal{A}\bm{F}+\bm{b})_{i}\Big{)} whose mean equals the sum of detected photons emitted by fluorescent molecules (i.e., $\mathcal{A}\bm{F}$) and background flux (i.e., $\bm{b}$).

Therefore, the nonlinear image-formation model for dipole emitters parameterized by brightness, 3D orientation, rotational dynamics, and position can be represented by a linear model utilizing a set of new parameters, accompanied by nonlinear constraints.

3 Recovery algorithm

Next, we focus on the regularized inverse problem for the derived model and show how we can handle nonlinear constraints to obtain efficient recovery in the presence of image overlap and experimental mismatches. The proposed algorithm consists of three main stages (Fig. 1(c)) as follows.

3.1 Joint sparse basis deconvolution

Our main observation is that the number of underlying molecules is much smaller than the number of ambient parameters. Taking advantage of the cascading structure of the proposed signal model, this prior knowledge can be encoded via a group/joint-sparsity norm across the six basis images, hence the name joint sparse basis deconvolution.

We formulate the recovery problem as a minimization:

[TABLE]

where $\mathcal{L}$ is the negative Poisson log-likelihood; $\mathrm{R}$ denotes the regularizer (described below); $\lambda$ is a penalty parameter; and $\mathcal{I}_{\mathcal{C}}$ denotes the indicator function of the set $\mathcal{C}$ representing the constraint set. The regularizer, $\mathrm{R}$, is a group-sparsity norm for brightness and position gradients across all six basis images defined by \mathrm{R}(\bm{F})=\sum_{i=1}^{N}\sqrt{\sum_{j=1}^{6}\Big{(}(\eta^{j}_{i})^{2}+(\zeta^{j}_{x,i})^{2}+(\zeta^{j}_{y,i})^{2}}\allowbreak\Big{)}. The set $\mathcal{C}=\cap_{i=1}^{6}C_{i}$ captures physical constraints on the signal model. In particular, $C_{j}=\{\bm{f}^{j}\ |\ \bm{\eta}^{j}\geq 0,\ -\rho\bm{\eta}^{j}\leq\bm{\zeta^{j}_{x}}\leq\rho\bm{\eta}^{j},\ -\rho\bm{\eta}^{j}\leq\bm{\zeta^{j}_{y}}\leq\rho\bm{\eta}^{j}\}$ for $j=\{1,2,3\};\ C_{j}=\{\bm{f}^{j}\ |\ \|\bm{\zeta^{j}_{x}}\|_{\infty}\leq\rho\|\bm{\eta^{j}}\|_{\infty},\ \|\bm{\zeta^{j}_{y}}\|_{\infty}\leq\rho\|\bm{\eta^{j}}\|_{\infty}\}$ for $j=\{4,5,6\}$. One may note that the sets $C_{j}$ ($j=\{4,5,6\}$) are non-convex. To preserve convexity, we neglect the first-order approximation in the last three basis images, implying that $\bm{f}^{j}=[\bm{\eta}^{jT},\bm{0}^{T},\bm{0}^{T}]^{T}$ for $j=\{4,5,6\}$. Henceforth we drop the corresponding constraint sets. This approximation is physically justified since the energy radiated by a dipole is mostly contained in the first three basis images [5, 10]. Consequently, it renders Eq. (3) as a convex program, and we develop a variant of FISTA to solve it [6, 11].

Deriving a closed-form proximal operator required for FISTA in Eq. (3) is not possible. We thus smooth the regularizer term with a differentiable function, e.g., its Moreau envelope. Let $w(\bm{F})=\lambda\mathrm{R}(\bm{F})$. The Moreau envelope of $w$, $\mathcal{E}_{\tau}(w)$, is continuous, and its gradient is given as $\nabla w_{\tau}(\bm{F})=\frac{1}{\tau}(\bm{F}-\text{prox}_{\tau w}(\bm{F}))$ in which $\text{prox}_{\tau w}(\bm{F})$ is the proximal operator of $\tau\lambda\mathrm{R}(\bm{F})$. The smoothing parameter $\tau$ controls the accuracy of the approximation. The modified optimization is now an instance of the gradient projection given as $\min_{\bm{F}}\mathcal{L}(\bm{F})+\mathcal{E}_{\tau}(\lambda\mathrm{R}(\bm{F}))+\mathcal{I}_{\mathcal{C}}(\bm{F})$, for which $\mathcal{I}_{\mathcal{C}}(\cdot)$ is the indicator of the set ${\mathcal{C}}=\cap_{i=1}^{3}C_{i}$. Another difficulty is that the projection operator associated with ${\mathcal{C}}$ does not admit a closed-form expression. To tackle this issue, we approximate each $C_{j}$ with a second-order cone as $\widehat{C}_{j}=\{\bm{f}^{j}\ |\ \|[\zeta^{j}_{x,i},\zeta^{j}_{y,i}]^{T}\|_{2}\leq\rho\eta_{i}^{j},\ i=\{1,2,\ldots,N\}\}$ ($j=\{1,2,3\}$), which quite remarkably admits a closed-form projection operator [6].

3.2 Spatial pooling for robust emitter identification and localization

An important challenge in sparse recovery is that of model mismatch, which has been shown to degrade the performance of grid-based recovery algorithms. However, for the problem of emitter localization, a model mismatch can cause false emitter localization, thereby introducing bias in the measurements (Fig. 2(a)). In our imaging system, the misalignment of two polarization channels leads to a mismatch; the estimated positions obtained separately from left and right channels may differ by an amount comparable to the localization precision. Notice, however, that the average of position estimates in two channels can actually be considered as the main parameter of interest. Using this insight and the fact that our model can provide continuous position estimates, we next show how to robustly identify the correct number of molecules.

As shown in [6], the recovered joint signal $\hat{\bm{F}}$ exhibits a specific structure in which position gradients converge to the true location of each emitter. To exploit this structure, we construct an operator called GradMap $\mathcal{G}:\mathbb{R}^{3N}\rightarrow\mathbb{R}^{N}$. Briefly, for each grid point or pixel, $\mathcal{G}$ computes how much the position gradients in neighbouring points converge to the grid point of interest. Put differently, $\mathcal{G}$ returns the “likelihood” that a true emitter belongs to each pixel. We apply $\mathcal{G}$ on each of first three deconvolved basis images to obtain $G^{j}=\mathcal{G}(\bm{\hat{f}}^{j})$ for $\bm{\hat{f}}^{j}=[\bm{\hat{\eta}}^{jT},\bm{\hat{\zeta}}_{x}^{jT},\allowbreak\bm{\hat{\zeta}}_{y}^{jT}]^{T}$ and $j=\{1,2,3\}$ (Fig. 2(b)). The final position estimates are obtained by pooling these three spatial maps, i.e., $G=\frac{1}{3}\sum_{i=1}^{3}G^{i}$ (Fig. 2(b)). The emitters’ initial positions and brightnesses correspond to the local maxima of $G$.

Let $\text{Supp}(\hat{\bm{F}})$ be the support set of estimated emitters’ positions obtained via spatial pooling. Localization is achieved via solving a constrained maximum-likelihood problem with an initial point obtained in the previous stage:

[TABLE]

which is a convex program. The molecular orientation parameters are estimated using the recovered second-moment vectors $\hat{\bm{\eta}}=\hat{s}\widehat{\bm{\mathcal{M}}}$.

4 Results

Fluorescent beads were utilized to validate our method’s detection and localization capabilities. A sparse layer of beads was imaged using the standard and Tri-spot PSFs with the system described in [5] (Fig. 2). Despite the substantial overlap and channel mismatch, as evidenced by the offset of the green dots from the peaks of the standard PSF (Fig. 2(c)), the algorithm is able to recover the correct position and number of emitters using the Tri-spot PSF (Fig. 2(d)).

$\lambda$-DNA (Thermo Scientific) was deposited onto coverslips using a molecular combing technique [12] and imaged using a reducing-oxidizing buffer [13] and Tri-spot PSF. The proposed method recovered the location, orientation, and rotational mobility $\gamma$ ($\gamma$ = 0 for an isotropic emitter and $\gamma$ = 1 for a fixed dipole) of many blinking YOYO-1 dyes transiently bound to DNA (Fig. 3(a,b)). The long-axis of the DNA strand was estimated by a simple least-squares polynomial fit, and the orientation of the dyes bound to the strand ($\Delta\phi$) was calculated relative to this fit. As shown in Fig. 3(b), molecules detected along a linear section of DNA are mostly oriented perpendicular to the DNA axis, which is consistent with the primary binding mode of YOYO-1 [14, 15]. The distribution of binding angles is likely broadened by secondary dye binding modes along the major and minor grooves of DNA, i.e., along its axis [14]. These modes have been reported at high concentrations of YOYO-1 and can be observed in regions of high localization density (Fig. 3(b)), where some dipoles are oriented parallel to the strand. The recovered rotational constraint yielded an average “wobble” cone angle of 91*∘*, assuming uniform rotation within a cone. Our 3D estimate of wobble angle is larger than previously observed for YOYO-1 [15] and SYTOX Orange [12] (another DNA intercalating dye), which were 2D orientation measurements. As YOYO-1 is a bis-intercalator, each molecule is effectively two dipole emitters with similar in-plane orientations but offset out-of-plane orientations; this superposition of dipoles likely causes a larger effective wobble angle in 3D and therefore an overestimate of the true rotational mobility of each bis-intercalator. Previous 2D studies of YOYO-1 were not sensitive to this out-of-plane offset.

The algorithm was also able to distinguish areas of relative organization and disorder. Regions where the average dye orientation was not perpendicular to the DNA axis indicate local fluctuations, or tangles, in the DNA strand (Fig. 3(e)) that are not observable via standard SMLM (Fig. 3(d)).

5 Conclusion

Engineered PSFs for orientation-sensitive super-resolution imaging pose major challenges, such as frequent overlapping PSFs and channel registration errors, for standard localization algorithms. We have presented a novel, robust algorithm for simultaneous recovery of the position and 3D orientation of fluorescent molecules using engineered PSFs. In contrast to methods based on defocus imaging, our algorithm can be applied to arbitrary orientation-sensitive PSFs and remedies mislocalizations due to PSF overlap and channel misalignment. We validated this method by imaging $\lambda$-DNA labeled with an intercalating dye, showing that the recovered molecular dipole orientation is primarily perpendicular to the DNA axis, which is consistent with previous observations. Interestingly, by measuring the full 3D orientation of YOYO-1, we observe a rotational constraint that is significantly smaller than that measured by 2D methods, suggesting that 3D orientation measurements may be necessary for revealing the true rotational dynamics of single molecules.