1. Introduction

Fluorophores absorb and emit polarized light in an orientation-dependent manner. Fluorescence polarization microscopy (FPM) techniques can be used to calculate fluorophore orientation from fluorescence intensity measurements. Widefield FPM can enable diffraction-limited orientation mapping across fluorescence microscopy images, which is particularly useful for characterizing molecular organization within living cells. One of the most common orientation-mapping techniques is fluorescence detected linear dichroism (FDLD) [1], which leverages the orientation-dependence of the intensity (I) of a fluorophore excited by plane polarized light:

FDLD is most commonly employed by using simple polarization optics to rotate ${\textbf E}$ by angle ${\upalpha }$ around a fixed optical axis that, as is standard in epifluorescence microscopy, is perpendicular to the sample plane (i.e. the x-y plane) and parallel to the z-axis. In this context:

 

Fig. 1. (A) VALiD microscope setup showing excitation laser beam (green), ${\textbf E}$ (orange), ${\bf {\upmu} }$ (red, superimposed on the chemical structure of an example fluorophore), rotating half wave plate (HWP), sliding lens (SL), and objective back focal plane (BFP). (B) $\textrm{I}$ as a function of ${\upbeta }$ calculated using Eq. (1) with ${\upalpha } = 45^\circ $, $\upphi = 60^\circ $, and ${\uptheta } = 45^\circ $ and illustrated for five distinct ${\upbeta }$ values. (C) Fluorophore azimuthal angle (${\varphi }$) and polar angle (${\uptheta }$) shown in cartesian coordinate system. (D) Bisected cone of $\textrm{E}$ orientations achievable when ${\upalpha } = 45^\circ $. (E) Five distinct fluorophore orientations and (F) corresponding intensity curves calculated using conventional FDLD (left) and VALiD (right).

Download Full Size | PDF

FDLD is a technologically accessible and versatile technique that has been used to image the orientation of fluorescent ensembles [1–10] and single fluorophores [11]. One limitation of FDLD is that it cannot be used to image fluorophores that are perpendicular to the sample plane (when ${\uptheta } = 0^\circ $, ${\Psi } = 90^\circ $ and, therefore, $\textrm{I} = 0$ for all ${\upalpha }$). Furthermore, and critically for this work, measurements of ${\upvarphi }$ are strictly two-dimensional because they are confined to the x-y plane. Any effort to measure three-dimensional orientation is thwarted by a fundamental limitation of this implementation of FDLD: All measurement vectors are co-planar with the x-y plane, which causes measurement degeneracy. In other words, any given ${\bf {\upmu} }$ cannot be distinguished from its reflection across the sample plane. To use an example, an orientation measurement of ${\upvarphi } = 45^\circ $, ${\uptheta } = 30^\circ $ is indistinguishable from ${\upvarphi } = 45^\circ $, ${\uptheta } = 120^\circ $ because the two orientations produce the same intensity at all ${\upalpha }$ values. This can generally be expressed mathematically using Eq. (2):

Several previously presented techniques that use the orientation dependence of emission and/or excitation for the measurement of individual fluorophores’ orientations are non-degenerate [11,15–26]. However, these techniques generally require fluorophores to be spatially and rotationally static and spatially isolated from each other. (While some techniques allow for limited rotational mobility [24–26], they still require individual fluorophores to be spatially isolated). These conditions are not always possible; there are many situations in which samples are densely labeled with fluorophores that exhibit localized organization. There are many applications in which it would be useful to measure the average 3D orientation of semi-ordered spatiotemporal ensembles of fluorophores. However, orientation mapping of densely labeled samples has thus far only been accomplished using, at best, four-fold degenerate techniques.

For example, we recently developed molecular force microscopy (MFM). MFM employs FDLD and DNA mechanotechnology [27] tension sensors to measure the orientation of traction forces at the cell-substrate interface with unprecedented spatial and temporal resolution [2,28]. The four-fold degeneracy of FDLD limits the power of this technique to provide biophysical insight by preventing direct mapping of force vectors; each measurement represents two distinct force orientations with equal likelihood, each of which has a distinct biological meaning. This limitation prevents useful forms of inquiry (such as force-balance analysis) that would be illuminating considering the increasingly apparent importance of force orientation in molecular biophysics research [29–32] and the increasing use of techniques like MFM in the development of active force generating nanomaterials [33,34].

Herein, we present the first non-degenerate technique, to the best of our knowledge, for ensemble fluorophore orientation mapping on densely labeled samples. This technique is a variant of FDLD that leverages a variable incidence polarized excitation beam and, as such, we call it variable incidence angle linear dichroism (VALiD).

VALiD eliminates measurement degeneracy by rotating ${\textbf E}$ in three-dimensions. This rotation can only be achieved by tilting the excitation beam with respect to the sample plane. One possible VALiD microscope setup is illustrated in Fig. 1(A): a rotatable half wave plate (HWP) is used to vary ${\upalpha }$ of a linearly polarized, collimated excitation laser beam. A sliding lens is then used to laterally shift the position of the excitation beam’s focal point on the objective lens’s back focal plane (BFP), which tilts the laser beam as it leaves the objective such that it enters the sample at angle ${\upbeta }$. While there are many ways to implement VALiD within this framework, we start by considering the simplest implementation (which we call configuration 1): ${\upalpha }$ is fixed at a constant value (${\upalpha } = {{\upalpha }_0}$) and multiple images are acquired at different distinct ${\upbeta }$ values. Setting ${{\upalpha }_0} = 90^\circ $ makes ${\textbf E}$ invariant with changing ${\upbeta }$ and setting ${{\upalpha }_0} = 0^\circ $ re-introduces degeneracy by restricting ${\textbf E}$ to y-z plane (thus making the y-z plane a new plane of symmetry). As such, an intermediate value of ${{\upalpha }_0}$ (such as $45^\circ $) must be used.

The expected intensity of a single static fluorophore imaged with ${{\upalpha }_0} = 45^\circ $ and varying ${\upbeta }$ angles is shown in Fig. 1(B). Note that varying ${\upbeta }$ rotates ${\textbf E}$ through a bisected cone with half angle ${{\upalpha }_0}$ centered on the x-axis (Fig. 1(D)). The measurement vectors in this cone are not coplanar, ensuring that each orientation of ${\bf {\upmu} }$ will produce a unique set of intensity measurements (Fig. 1(E),(F)). Theoretical VALiD curves for five static fluorophores with the same ${\upvarphi }$ but different ${\uptheta }$ calculated using ${{\upalpha }_0} = 45^\circ $ and ${\upbeta } \in [{ - 90^\circ ,90^\circ } ]$ are shown in Fig. 1(F) (right). For reference, we also show intensity curves for the same orientations calculated according to conventional FDLD, where ${\upbeta } = 0{\; }$ and ${\upalpha } \in [{0^\circ ,180^\circ } ]$ (Fig. 1(F) (left)). Note two key distinctions: first, in VALiD, all five curves are unique while curves generated from FDLD are degenerate; second, VALiD generates a curve for the ${\uptheta } = 0^\circ $ case while FDLD does not. In practice, three-dimensional fluorophore orientation can be calculated from experimental intensity measurements via least-squares residual fitting to Eq. (1).

2. Mathematical modeling and computational methods

In this work, we will use computational methods to predict the experimental conditions under which VALiD can be correctly implemented. We first derive idealized equations for the fluorescence intensity of axially symmetric distributions of fluorophores. We then perform Monte Carlo simulations of VALiD measurements to model the effect of noise on orientation measurements and understand the effects of various experimental conditions and parameters on measurement quality. Throughout this text, we use non-bold type for scalars (e.g. ${\upphi}, {\uptheta}, \textrm{A}$) and bold type for vectors (e.g. ${\bf {\upmu} }$, ${\textbf E}$). Vectors are described conventionally as 3 × 1 matrices. Unit vectors are defined with a hat; $|{\textbf v} |= \hat{{\textbf v}}$ . We use the standard x-y-z coordinate system to describe the microscope coordinate system and a local k-l-m coordinate system to describe fluorophore ensembles. Within this local coordinate system, we denote polar angle as ${\updelta }$ and azimuthal angle as ${\upgamma }$. These quantities are matched to specific vectors using subscripts (e.g. ${{\updelta }_{\upmu }}$ and ${{\updelta }_\textrm{E}}$). Finally, we will use rotation matrices; ${\textrm{R}_{\textbf v}}({\uptau } )$ is a 3 × 3 matrix that, when multiplied by a vector, rotates the vector by ${\uptau }$ degrees around the axis parallel to ${\textbf v}$. We can define the orientation of ${\textbf E}$ and ${\bf {\upmu} }$ with rotation matrices:

2.1 Fluorescence intensity of a conical ensemble

To enable the measurement of ensemble average orientation, we first derive equations for the expected fluorescence intensity of an ensemble of fluorophores that are uniformly distributed around a “symmetry axis”. The orientation of the symmetry axis is defined by the unit vector $\hat{{\textbf m}}$ (Fig. 2(A)), which has polar angle ${\Theta }$ and azimuthal angle ${\Phi }$ such that

 

Fig. 2. (A) Global coordinate system showing orientation of symmetry axis $\hat{{\textbf m}}$ denoted with polar angle ${\Theta }$ and azimuthal angle ${\Phi }$. (B,C) Local coordinate system centered on $\hat{{\textbf m}}$ showing orientation of (B) ${\bf {\upmu} }$ with local polar angle ${{\updelta }_{\upmu }}$ and local azimuthal angle ${{\upgamma }_{\upmu }}$ and (C) E with local polar angle ${{\updelta }_\textrm{E}}$ and local azimuthal angle ${{\upgamma }_\textrm{E}}$. (D) Depiction of ensembles of fluorophores distributed uniformly around symmetry axis $\hat{{\textbf m}}$ via conical and Watson distributions. These ensembles can respectively be exactly and approximately decomposed into a dipole and a disk. (E) Intensity curves of disk (blue) and dipole (orange) components of ensemble and dashed lines showing weighted averages for 15$^\circ $ increments of ${{\updelta }_0}$ with ${\Theta } = {\Phi } = {{\upalpha }_0} = 45 \circ $, and $\textrm{A} = 1$. (F) Surface plot showing the error (multiplied by 100 so that it loosely represents the percent error) of Eq. (12) when used to represent the Watson distribution (with $\textrm{A} = 1$). Because the error is low (<5%) for all combinations of ${\Delta }$ and ${{\updelta }_\textrm{E}}$, we conclude that Eq. (12) is a reasonable approximation for the intensity of Watson-distributed fluorescent ensembles.

Download Full Size | PDF

Assumption 1: We first assume that ${{\updelta }_{\upmu }}$ is a fixed constant such that fluorophore orientation is restricted to a cone centered on $\hat{{\textbf m}}$ (Fig. 2(D)). In this scenario, the ensemble average ${{\updelta }_{\upmu }}$ value, which we denote in this work as ${\Delta }$, is the same as the half angle of the cone (Fig. 2(D)). This assumption, which simplifies the forthcoming mathematics, is a reasonable approximation of many physically-relevant scenarios. For example, the commonly used fluorophore DiI has two lipid tails that readily insert into plasma membranes such that the transition dipole moment can rotate freely parallel to the plasma membrane. In this scenario, $\hat{{\textbf m}}$ is normal to the plasma membrane and ${\Delta } \approx 75^\circ $ [14,35].

Assumption 2: Second, we assume that each fluorophore’s rotational time constant is much shorter than the fluorescence lifetime. This assumption simplifies mathematics by removing the orientation dependence of fluorescence emission on the measured ensemble intensity. Put another way, if a fluorophore’s orientation within the cone is randomized between excitation and emission, then the effects of the orientation-dependence of emission will be equivalent for all orientations of ${\textbf E}$. It follows that, under this assumption, changes in fluorescence intensity will result strictly from changes in the orientation of ${\textbf E}$.

We also make three minor assumptions: first, we assume that photobleaching is negligible – which is reasonable for densely labeled samples that are excited with low power. Second, we assume that other time-dependent photophysical effects such as blinking are negligible, which is reasonable for ensembles of fluorophores where these effects can be averaged out and lumped in with the orientation-independent $\textrm{Q}$ term. Finally, we assume and that the transition dipole moment for emission is parallel to the transition dipole moment for excitation (${\bf {\upmu} }$), which is reasonable for most fluorophores. Unlike with assumptions 1 and 2, we do not study the effect of these assumptions within this work. However, we note that photobleaching corrections can be applied to densely labeled samples that do undergo meaningful photobleaching. We also recommend avoiding the use of fluorophores with non-parallel emission and excitation transition dipole moments for future implementations of this technique.

Using these assumptions, the intensity of an ensemble of fluorophores can be expressed in a relatively simple manner:

Next, we consider the local k-l-m coordinate system centered on $\hat{{\textbf m}}$ (Fig. 2(B),(C)). The k-l-m unit basis vectors can be obtained by rotating the unit basis vectors of the global x-y-z coordinate system:

2.2 Fluorescence intensity of a Watson ensemble

We next tested the necessity of assumption 1 by testing whether Eq. (12) is reasonable for the description of other axially-symmetric distributions. Specifically, we evaluated the use of Eq. (12) to represent fluorophores with orientations distributed according to the Watson distribution. The Watson distribution is an angular analog of the Gaussian distribution (Fig. 2(D)) and is defined by the probability density function (${\textrm{P}_{\textrm{Watson}}}$):

While we will not go into more sophisticated axially-symmetric distributions in this work, we also note that the intensity of any arbitrary axially-symmetric distribution can be calculated by treating the distribution as a summation of many cones and integrating over Eq. (5).

In addition, we note that when ${\Delta } = {\tan ^{ - 1}}\sqrt 2 \approx 54.7^\circ$ (also known as the “magic angle”) intensity ceases to depend on ${\textbf E}$’s orientation because $\textrm{dI}/\textrm{d}\Psi = 0$. We therefore predict that VALiD loses accuracy as ${\Delta }$ approaches the magic angle.

2.2 Microscope configurations

Below, we present three microscope configurations that could be used to implement VALiD. These three configurations exhibit qualitative differences in the set of measurement vectors that they can produce. We use a Monte Carlo simulation method – described below – to determine how well each configuration can be used to measure properties of simple axially symmetric fluorescent ensembles.

Configuration 1, which was described above, is the simplest implementation of VALiD, wherein ${\upalpha }$ is fixed at ${\upalpha } = {{\upalpha }_0}$ and $\textrm{N}$ images are acquired at ${\upbeta }$ values equally spaced between $- 60^\circ $ and $60^\circ $. We used a $60^\circ $ cutoff because higher incidence angles will result in a large proportion of excitation light being reflected at the sample-coverslip interface and may result in distortions to the polarization state and/or spatial profile of the excitation light.

Configuration 2 is similar to config. 1; $\textrm{N} - 1$ images are acquired at varying ${\upbeta }$ values with ${\upalpha } = {{\upalpha }_0}$. In addition, a single image is taken with ${\upalpha } = 90^\circ $ and ${\upbeta } = 0^\circ $. In practice, this configuration can be implemented by adding an appropriately oriented polarization modulator that can be manually inserted into and removed from the path of the excitation beam. Config. 2 is intended to overcome some of the specific limitations of config. 1 by adding an additional measurement vector that is not restricted to the ${\textbf E}$ cone.

Configuration 3 allows both ${\upbeta }$ and ${\upalpha }$ to vary freely and can be implemented using an electronically controlled polarization modulator (e.g. a HWP mounted in a rotational actuator or Pockels cell). This configuration is expected to maximize the spread of the vectors to equally cover the entire unit sphere. To simulate ${\textbf E}$ in this configuration, we generated $4\textrm{N} - 3$ vectors equally distributed on the unit sphere using the SpiralSampleSphere method [36] and then discarded all vectors except for those with polar angles between $30^\circ $ and $90^\circ $. As such, this approach mimics the situation where the maximum ${\upbeta }$ angle is $60^\circ $.

While these configurations are all based on a though-the-objective TIRF-style microscope wherein the excitation beam is tilted within a single plane of incidence, VALiD could also be implemented using alternative instruments such as azimuthal scanning fluorescence microscopes [37].

2.3 Monte Carlo simulation method

We used a Monte Carlo simulation method to evaluate each configuration over a range of experimental parameters. To summarize this method, we generate simulated intensity data for a given fluorescent ensemble and set of measurement vectors, add noise to the intensity data, and then use least-squares residual fitting to measure the best-fit orientation of the ensemble’s symmetry axis (Fig. 3(A)). We then repeat this process many times to measure the average accuracy of a given set of measurement vectors for a given fluorescent ensemble. This process is then repeated for a representative set of symmetry axis orientations to measure the overall performance of a set of measurement vectors. Finally, we repeated this entire process for different microscope configurations and different types of ensembles. The Monte Carlo algorithm is described in greater detail below. Least-squares residual fitting is performed using a two-layer lookup table approach. Steps 1-3 of the algorithm describe initialization of the lookup table. Steps 4-6 describe the process used for a single measurement simulation, and steps 7-9 describe various repetitions used to evaluate a wide parameter space.

 

Fig. 3. (A) Depiction of a single iteration of the Monte Carlo algorithm used for error estimation: ${\textrm{I}_{\textrm{ideal}}}$ is calculated with Eq. (12) for a given set of ${\Theta }$, ${\Phi }$, ${\Delta }$, and $\textrm{A}$ values (the curve shown was calculated with ${\Theta } = 30^\circ $, ${\Phi } = 40^\circ $, ${\Delta } = 0^\circ $, and $\textrm{A} = 1$). Noise is then added to the calculated ${\textrm{I}_{\textrm{ideal}}}$ values, and then the true orientation is fit to the noisy data. Finally, the angular error ($\varepsilon $) of the best-fit orientation is calculated. (B-M) Monte Carlo simulation outputs collected with ${\Delta } = 0^\circ $, ${{\upalpha }_0} = 30^\circ $ (when applicable), and $\textrm{A}$ and ${\Delta }$ treated as known parameters. (B-D) Representative heatmaps of orientation measurements for a given true orientation denoted by x (${\Theta } = 40^\circ $, ${\Phi } = 45^\circ $) with $\textrm{A} \times \textrm{N} = 700$ for (B) configuration 1 with $\textrm{N} = 4$ (C) configuration 2 with $\textrm{N} = 4$, and (D) configuration 3 with $\textrm{N} = 7$. The x may appear to deviate from the maximum likelihood measurement due to artifacts associated with data display. (E-M) Surface plots showing the (E-G) average error ($\left\langle {\upvarepsilon } \right\rangle $), (H-J) error anisotropy (${\upsigma }$), and (K-M) error ellipses as a function of orientation with $\textrm{A} \times \textrm{N} = 1,000$ for (E, H, K) configuration 1 with $\textrm{N} = 5$, (F, I, L) configuration 2 with $\textrm{N} = 5$, and (G, J, M) configuration 3 with $\textrm{N} = 10$. All surface plots have black dots showing intersection of the unit sphere with the z-axis (top dot), x-axis (bottom left dot), and y-axis (bottom right dot) and the plots in (B-J) have a black circle showing the intersection of the unit sphere with the x-y plane.

Download Full Size | PDF

We evaluated the accuracy of configs. 1-3 under different fitting conditions; ${\Theta }$ and ${\Phi }$ were always treated as unknown fit parameters, but $\textrm{A}$ and ${\Delta }$ were optional additional unknown fit parameters. The situation where $\textrm{A}$ is known prior to acquisition is not experimentally practical. As such, these simulations were primarily run to isolate fundamental characteristics of the orientation measurement process. When $\textrm{A}$ was an additional unknown parameter, each set of ${\textrm{I}_{\textrm{ideal}}}$ lookup table values (steps 1-2) or ${\textrm{I}_{\textrm{noise}}}$ values (step 4) was normalized to the average of the set. When ${\Delta }$ was an additional unknown parameter, step 5 was repeated 31 times (one for each $3^\circ $ increment of ${\Delta }$ from $0^\circ $ to $90^\circ $) and the best of the 31 best-fits was selected. Note that ${\Delta }$ could be estimated prior to VALiD acquisition using a variety of techniques or via existing reports in the research literature. For instance, fluorescence anisotropy techniques have been used to estimate ${\Delta }$ for DiI in supported lipid bilayers as ∼75° [14,35].

To maintain a consistent SNR when comparing sets of measurement vectors with different $\textrm{N}$, we defined $\textrm{A} \times \textrm{N}$ as a proxy for the “total photon count”. This definition imposes a tradeoff with increasing $\textrm{N}$ between the benefit obtained by using more measurement vectors and the cost of decreasing the SNR for each individual intensity measurement.

2.4 Monte Carlo simulation assessment metrics

In addition to quantifying the average error for each orientation, we also quantified the error anisotropy. In this context, we define error anisotropy as the tendency of noise to bias orientation measurements preferentially towards some orientations. We agree with Chandler et al. that the ideal technique should be able to “reconstruct the orientation with a small and nearly uniform uncertainty” for all orientations [16]. Substantial error anisotropy is likely to result in systematic errors and “sink” orientations that will be measured more commonly than others. For example, estimation of ${\uptheta }$ using FDLD with Eq. (2) results in fluorescent ensembles with low ${\uptheta }$ angles being systematically over-estimated [2]. In this case, ${\uptheta }$ overestimation occurs because noise generally increases the amplitude of a best-fit sinusoid, particularly when the amplitude of the true sinusoid is very small. In other cases, high error anisotropy may have other causes. Regardless of the underlying cause, our Monte Carlo simulation method allows us to quantify error anisotropy in a context-independent manner.

Error anisotropy can be roughly visualized by analyzing the distribution on a unit sphere of 30,000 measured orientations for a given true orientation of $\hat{{\textbf m}}$. When error is anisotropic, measurements will be concentrated along a given direction (e.g. around a given $\upgamma $ angle) on the unit hemisphere and the density map will resemble an elliptical Gaussian. When error is isotropic, measurements will be uniformly distributed around the true orientation and the density map will resemble an isotropic Gaussian. Figures 3(B) and 3(D) show examples of anisotropic and isotropic errors, respectively. Error anisotropy ($\upsigma $) can be quantified by fitting the distribution of measured orientations to an uncertainty ellipse. This calculation is performed by assembling the k and l components of all measured ${\hat{{\textbf m}}^{\ast }}$ values into two “component vectors”, (denoted ${\textbf K}$ and ${\textbf L}$, respectively, where each is a $\textrm{n} \times 1$ vector) calculating the covariance matrix (C) of these component vectors:

After calculating $\left\langle {\upvarepsilon } \right\rangle $ and ${\upsigma }$ for each of 360 representative orientations, we then calculated five assessment metrics that together reflect the performance of a given set of measurement vectors for a given pair of $\textrm{A}$ and ${\Delta }$ parameters. Those metrics are the average and maximum of the 360 calculated $\left\langle {\upvarepsilon } \right\rangle $ values (${\left\langle {\upvarepsilon } \right\rangle _{\textrm{mean}}}$ and ${\left\langle {\upvarepsilon } \right\rangle _{\textrm{max}}}$), the average and maximum of the 360 calculated ${\upsigma }$ values (${{\upsigma }_{\textrm{mean}}}$ and ${{\upsigma }_{\textrm{max}}}$), and finally an “overall quality” metric $\textrm{q}$, where

3. Results

3.1 Initial configuration assessment

We first ran representative Monte Carlo simulations for configs. 1, 2, and 3 to obtain a qualitative understanding for the contexts in which each configuration could be reasonably used. For our initial assessment, we set ${\Delta } = 0$, $\; {\upalpha _0} = 30^\circ $ (for configurations 1 and 2), ${\Theta } = 40^\circ $, and ${\Phi } = 45^\circ $. We also set $\textrm{A} \times \textrm{N} = 700$, which is an experimentally relevant value that would be expected for an ensemble consisting of a few dozen fluorophores. (This value is also somewhat low, ensuring that the effect of noise is meaningful). We treated both $\textrm{A}$ and ${\Delta }$ as known parameters. A quick interpretation of Eq. (12) suggests that at least four measurement vectors must be used to achieve Nyquist sampling. Therefore, for configurations 1 and 2 we used $\textrm{N} = 4$ measurement vectors. For config. 3, we used $\textrm{N} = 7$ measurement vectors to obtain more uniform sampling. We ran $\textrm{n} = 30,000$ iterations under these conditions and plotted errormaps showing the probability of measuring any given orientation when the true orientation is fixed (Fig. 3(B)–(D)). As expected, we see a single node of measured orientations centered around the true orientation, thus confirming that all three configurations can achieve nondegenerate orientation measurement. (As a control, we ran the same process with $\textrm{N} = 4$ for FDLD and observed two nodes separated by a $180^\circ $ rotation around the z-axis, not shown). Noise results in deviations from the true orientation such that $\left\langle {\upvarepsilon } \right\rangle $ is somewhat consistent between the three configurations ($11.0^\circ $ for C1, $7.9^\circ $ for C2, and $8.2^\circ $ for C3). However, the orientational nature of these deviations vary; for config. 3, error is highly isotropic ($\upsigma = 0.19$), producing measurement errors that are equally probable in all directions (Fig. 3(D)). For configs. 1 and 2, error is highly anisotropic ($\upsigma = 0.85$ and $0.64$ respectively – Fig. 3(B), (C)). These results suggest that spreading the measurement vectors across a wide variety of orientations shields against error anisotropy. These results are also in line with recent theoretical analyses [38].

We next varied $\hat{{\textbf m}}$ across a representative set of 360 orientations and measured how ${\upsigma }$ and $\left\langle {\upvarepsilon } \right\rangle $ vary as a function of orientation for all three configurations. For these simulations we used $\textrm{A} \times \textrm{N} = 1,000$ and ${\Delta } = 0^\circ $ and treated $\textrm{A}$ and ${\Delta }$ as known parameters. For configs. 1 and 2 we used $\textrm{N} = 5$ and ${{\upalpha }_0} = 30^\circ $, and for config. 3 we used $\textrm{N} = 10$. Unit sphere colormaps are shown in Fig. 3(E)–(J). To show the nature of error anisotropy across the entire unit sphere, we also plotted error ellipses for 100 representative orientations where the orientation of each ellipse shows the primary orientation of measurement error and the lengths of the major and minor axes are proportional to $\sqrt {{{\uplambda }_1}} $ and $\sqrt {{{\uplambda }_2}} $, respectively (Fig. 3(K)–(M)).

These results provide a more detailed depiction of measurement error; $\left\langle {\upvarepsilon } \right\rangle $ depends heavily on orientation for config. 1 (Fig. 3(E)), ranging from $5.4^\circ $ to ${\left\langle {\upvarepsilon } \right\rangle _{\textrm{max}}} = 25.0^\circ $ (with ${\left\langle {\upvarepsilon } \right\rangle _{\textrm{mean}}} = 11.9^\circ $). Generally, orientations that are proximal to the half-cone of measurement vectors are measured much more accurately, while orientations that are further from the half-cone are measured less accurately because they are poorly excited. Config. 1 also produces highly anisotropic error for nearly all $\hat{{\textbf m}}$ orientations, with ${{\upsigma }_{\textrm{mean}}} = 0.84$ and ${{\upsigma }_{\textrm{max}}} = 0.96$ (Fig. 3(H), (K)). Config. 2 generally is associated with lower error, with ${\left\langle {\upvarepsilon } \right\rangle _{\textrm{max}}} = 13.1^\circ $ and ${\left\langle {\upvarepsilon } \right\rangle _{\textrm{mean}}} = 7.0^\circ $ (Fig. 3(F)). Config. 2 also produces lower error anisotropy, with ${{\upsigma }_{\textrm{mean}}} = 0.60$ and ${{\upsigma }_{\textrm{max}}} = 0.96$, than config. 1 (Fig. 3(I), (L)). Measurement of many $\hat{{\textbf m}}$ orientations is improved with config. 2, but some orientations – particularly, those that are close to the x-y plane – exhibit relatively high error and measurement anisotropy. These results suggest that the additional measurement vector parallel to the y-axis is sufficient to mitigate some, but not all, of the shortcomings associated with config. 1. In contrast, config. 3 exhibits low error (${\left\langle {\upvarepsilon } \right\rangle _{\textrm{mean}}} = 7.4^\circ $ and ${\left\langle {\upvarepsilon } \right\rangle _{\textrm{max}}} = 9.1^\circ $) and relatively low error anisotropy (${{\upsigma }_{\textrm{mean}}} = 0.44$ and ${{\upsigma }_{\textrm{max}}} = 0.75$) across all orientations (Fig. 3(G), (J), M).

The improved performance of config. 3 ($\textrm{q} = 5.10$) over config. 1 ($\textrm{q} = 0.83$) and config. 2 ($\textrm{q} = 3.82$) further demonstrates the benefits of spreading the set of measurement vectors out to cover a wider range of orientations. However, the improved performance comes at the cost of increasing complexity of the microscope system. It is worth noting that these simulations were performed with a relatively low SNR (only a few hundred photons per intensity measurement), and the issues with configs. 1 and 2 may not be as prominent at the higher SNR values that will often be accessible for densely labeled samples.

3.2 Optimization of VALiD parameters

We next sought to determine optimized $\textrm{N}$ and ${{\upalpha }_0}$ values for each of the three configurations. As stated above, increasing $\textrm{N}$ adds more unique measurement vectors, but also reduces the SNR of each measurement (when $\textrm{A} \times \textrm{N}$ is held constant) and increases the amount of time required to obtain a full set of measurements. To determine optimal parameters, we used the Monte Carlo strategy described above to calculate the assessment metrics as a function of $\textrm{N}$ for configs. 1, 2, and 3 with $\textrm{A} \times \textrm{N}$ held constant and ${{\upalpha }_0} = 30^\circ $ (when applicable). We also performed simulations for configs. 1 and 2 with ${{\upalpha }_0}$ ranging from $0^\circ $ to $90^\circ $ with $\textrm{N} = 5$. As a control, we ran simulations of FDLD with varying $\textrm{N}$, in which there are $\textrm{N}$ distinct ${\upalpha }$ values ranging from 0 to $({\textrm{N} - 1} )\times 180^\circ $ and ${\upbeta } = 0^\circ $. For these simulations we selected $\textrm{A} \times \textrm{N} = 3,000$, set ${\Delta } = 0^\circ $, and ran simulations with different sets of unknown parameters. Optimal N and ${\upalpha _0}$ values were defined as the values that maximized $\textrm{q}$. Our results are shown in Fig. 4. Assessment metrics at optimized conditions are summarized in Table 1.

 

Fig. 4. Results of Monte Carlo simulation-based optimization. We refer readers to section 3.2 for a more detailed description of the results shown here. Each subplot contains three curves showing five assessment metrics (section 2.4) as a function of the number of images ($\textrm{N}$) or the measurement vector cone angle (${{\upalpha }_0}$). Blue curves show results when both $\textrm{A}$ and ${\Delta }$ were known, green curves show results when ${\Delta }$ was known and $\textrm{A}$ was unknown, and red curves show results when both $\textrm{A}$ and ${\Delta }$ were unknown. Simulations were performed with $\textrm{A} \times \textrm{N} = 3,000$ and ${\Delta } = 0$ with all integer values of N ranging from 3 to 20 (left four columns) and 5$^\circ $ increments of ${{\upalpha }_0}$ ranging from $0^\circ $ to $90^\circ $. When testing for optimal $\textrm{N}$, ${{\upalpha }_0} = 30^\circ $ was used. When testing for optimal ${{\upalpha }_0}$, $\textrm{N} = 5$ was used. To summarize the results, FDLD always fails due to measurement degeneracy, config. 3 always performs better than configs. 1 or 2, and config. 2 always performs better than config. 1. Configs. 1, 2, and 3 were found to have $\textrm{N} = 4$, 5, and 6 optimal, respectively, regardless of the number of free fit parameters. Configurations 1 or 2 were generally found to have an optimum of ${{\upalpha }_0} = 30^\circ $, although in some cases the ${{\upalpha }_0} = 25^\circ $ or ${{\upalpha }_0} = 35^\circ $ exhibited an up to 1% increase in quality ($\textrm{q}$) over the ${{\upalpha }_0} = 30^\circ $ condition.

Download Full Size | PDF

Table 1. Optimal experimental parameters and assessment metrics under optimal conditions.

View Table

We observed failure, which we define as occurring when ${\left\langle {\upvarepsilon } \right\rangle _{\textrm{max}}} \ge 45^\circ $, for FDLD in all scenarios due to the four-fold degenerate nature of FDLD measurements. This finding was expected, thus partially validating our simulation approach. Config. 1 performed well only when $\textrm{A}$ and ${\Delta }$ were known, and the optimal condition was $\textrm{N} = 4$ and ${{\upalpha }_0} = 30^\circ $. Config. 1 performed poorly when ${\rm A}$ was treated as an unknown parameter and failed when ${\Delta }$ was unknown. This failure was likely due to specific degeneracies that arise when constraining ${\textbf E}$ vectors to a half-cone. One such degeneracy is between the case where $\hat{{\textbf m}}$ is parallel to the y-axis and the case where ${\Delta } = 54.7^\circ $, both of which produce intensity values that are invariant with respect to ${\upbeta }$. When $\textrm{A}$ and ${\Delta }$ were known, performance was relatively insensitive to ${{\upalpha }_0}$ over a range of values from ∼$20^\circ $ to $45^\circ $. For config. 1, error anisotropy was always high. Taken together, these results suggest that config. 1 can be used with caution when $\textrm{N} = 4$ and ${{\upalpha }_0} \approx 30^\circ $ when ${\Delta }$ is known, but that config. 1 should not be used when ${\Delta }$ is unknown. While prior knowledge of $\textrm{A}$ helps to improve accuracy, it does not improve error anisotropy and so results obtained using config. 1 should also be assumed to exhibit anisotropic error. Furthermore, due to inherent variability in fluorescent samples, accurate estimation of $\textrm{A}$ is impractical without additional acquisitions with non-polarized light or with ${\textbf E}$ vectors that are not restricted to the cone of measurement vectors – as is exactly the case with config. 2.

Config. 2 demonstrated improved performance over config. 1. The optimal $\textrm{N}$ value increased to 5, with ${\left\langle {\upvarepsilon } \right\rangle _{\textrm{mean}}}$, ${\left\langle {\upvarepsilon } \right\rangle _{\textrm{max}}}$, ${{\upsigma }_{\textrm{mean}}}$, and ${{\upsigma }_{\textrm{max}}}$ all increasing monotonically with $\textrm{N}$ when $\textrm{N} > 5$. This unit increase in $\textrm{N}$ means that the optimal acquisition for config. 2 is the same as the optimal acquisition for config. 1, but with one additional measurement vector (with ${\upalpha } = 90^\circ $ and ${\upbeta } = 0^\circ $). Unlike config. 1, config. 2 performed well when $\textrm{A}$ was added as an additional unknown parameter. This finding suggests that config. 2 is inherently well-suited for simultaneous estimation of $\textrm{A}$ and orientation, which is attractive for practical purposes. However, although config. 2 exhibits slightly lower error anisotropy than config. 1, it still exhibits substantial error anisotropy for all $\textrm{N}$ and, as such, this configuration should also be used with a degree of caution. Config. 2 essentially failed (with ${\langle{\upvarepsilon }\rangle_{\textrm{max}}}$ near $45^\circ $) when ${\Delta }$ and $\textrm{A}$ were both unknown. Accordingly, our results suggest that config. 2 is well suited towards orientation measurement when ${\Delta }$ is known, but not when ${\Delta }$ is unknown. Optimization of ${{\upalpha }_0}$ showed that ${{\upalpha }_0} = 25^\circ $ is optimal when $\textrm{A}$ and ${\Delta }$ are known and ${{\upalpha }_0} = 35^\circ $ is optimal when $\textrm{A}$ is unknown and ${\Delta }$ is known. However, as with config. 1 we found that performance was relatively insensitive to ${{\upalpha }_0}$ over a range of values from ∼$20^\circ $ to ∼$45^\circ $, and in both cases the optimal $\textrm{q}$ was <1% higher than $\textrm{q}$ when ${{\upalpha }_0} = 30^\circ $. As such, we can reasonably conclude that ${{\upalpha }_0} = 30^\circ $ is optimal for both config. 1 and config. 2. It is not immediately clear why ${{\upalpha }_0} = 30^\circ $ is optimal, but we note that $30^\circ $ is the angle that minimizes the average angular difference between all orientations and the cone of measurement vectors (shown in Fig. 1(D)).

We observed the best performance with config. 3, for which there was another unit increment in the optimal number of measurement vectors to $\textrm{N} = 6$ (the specific set of [${\upalpha },{\upbeta }]$ values was [$152^\circ , - 34^\circ ]$, [$20^\circ , - 12^\circ $], [$115^\circ ,0^\circ $], [$158^\circ ,6^\circ $], [$58^\circ ,34^\circ $], [$118^\circ ,59^\circ $]). Unlike with configs. 1 and 2, config. 3 did not fail when both ${\Delta }$ and $\textrm{A}$ were treated as unknown parameters. Config. 3 was also unique in that increasing $\textrm{N}$ beyond the optimum of 6 generally caused ${{\upsigma }_{\textrm{mean}}}$ and ${{\upsigma }_{\textrm{max}}}$ to decrease. While increasing $\textrm{N}$ came with the cost of increased ${\left\langle {\upvarepsilon } \right\rangle _{\textrm{mean}}}$ and ${\left\langle {\upvarepsilon } \right\rangle _{\textrm{max}}}$, the increase in error was modest. For example, when $\textrm{A}$ was unknown and ${\Delta }$ was known, increasing $\textrm{N}$ from 6 to 20 caused ${{\upsigma }_{\textrm{max}}}$ to decrease from 0.84 to 0.73, while ${\left\langle {\upvarepsilon } \right\rangle _{\textrm{max}}}$ increased from $4.3^\circ $ to $6.0^\circ $. However, the benefit of increasing $\textrm{N}$ was less pronounced when both ${\Delta }$ and $\textrm{A}$ were unknown; in the cases where ${\Delta }$ was known, ${{\upsigma }_{\textrm{max}}}$ trended toward ∼0.7, but when ${\Delta }$ was unknown, ${{\upsigma }_{\textrm{max}}}$ trended towards ∼0.94. The similarity in the two ${{\upsigma }_{\textrm{mean}}}$ vs. $\textrm{N}$ curves when $\textrm{A}$ was unknown suggests that only a small subset of orientations exhibited highly anisotropic error when ${\Delta }$ was unknown, potentially due to specific degeneracies that arise when ${\Delta }$ is added as an additional free parameter. Together, these results suggest that 1) config. 3 can be used to obtain accurate orientation measurements, 2) reasonably isotropic orientation measurements can generally be obtained when ${\Delta }$ is known, and 3) the error anisotropy can be decreased – to an extent – at the cost of slightly reduced accuracy by increasing $\textrm{N}$. These results also suggest that anisotropic error is a fundamental issue with VALiD that will be difficult to fully address for practical scenarios when $\textrm{A}$ is unknown. However, we note that when the error is very low the fact that it is highly anisotropic has little practical relevance.

3.3 VALiD accuracy as a function of total photon count

To assess the manner in which the quality of VALiD-based orientation measurement scales with SNR, we measured ${\left\langle {\upvarepsilon } \right\rangle _{\textrm{max}}}$ and ${{\upsigma }_{\textrm{max}}}$ with $\textrm{A} \times \textrm{N}$ values ranging from 10² to 10⁵ for all three configurations using optimized conditions ($\textrm{N} = 4$, ${{\upalpha }_0} = 30^\circ $ for config. 1, $\textrm{N} = 5$, ${{\upalpha }_0} = 30^\circ $ for config. 2, and $\textrm{N} = 6$ for config. 3) and ${\Delta } = 15^\circ $. We ran these simulations under two fitting conditions: with $\textrm{A}$ unknown and ${\Delta }$ known, or with both $\textrm{A}$ and ${\Delta }$ unknown (Fig. 5(A)). The results display principles discussed above: All configs. can accurately measure orientation when ${\Delta }$ is known and SNR is high; config. 2 is much more accurate than config. 1; and config. 3 is substantially more accurate than config. 2; only config. 3 can measure orientation when ${\Delta }$ is unknown (although at extremely high SNR config. 2 appears to perform satisfactorily) and ${{\upsigma }_{\textrm{max}}}$ is generally close to 1, except when ${\Delta }$ is known for config. 3.

 

Fig. 5. VALiD performance as a function of $\textrm{A} \times \textrm{N}$ and ${\Delta }$. (A) ${\langle{\upvarepsilon }\rangle_{\textrm{max}}}$ (top) and ${{\upsigma }_{\textrm{max}}}$ (bottom) as a function of $\textrm{A} \times \textrm{N}$ with ${\Delta } = 15^\circ $ for config. 1 with $\textrm{N} = 4$ and ${{\upalpha }_0} = 30^\circ $ (red, squares), Config. 2 with $\textrm{N} = 5$ and ${{\upalpha }_0} = 30^\circ $ (green, circles), and config. 3 with $\textrm{N} = 6$ (blue, triangles) with $\textrm{A}$ treated as an unknown fit parameter and ${\Delta }$ known (solid curves) or ${\Delta }$ unknown (dashed curves). (B) ${\langle{\upvarepsilon }\rangle_{\textrm{max}}}$ as a function of $\textrm{A} \times \textrm{N}$ and ${\Delta }$ for config. 3 with $\textrm{N} = 6,$ $\textrm{A}$ unknown, and ${\Delta }$ known. Black curves show contours where error is $30^\circ $, $10^\circ $, $3^\circ $, and $1^\circ $. Simulations were performed with ${\Delta }$ varying in 2.5$^\circ $ increments from $0^\circ $ to $90^\circ $, and with 40 $\textrm{A} \times \textrm{N}$ values distributed logarithmically from 10² to 10⁵.

Download Full Size | PDF

If we define a cutoff photon count ${({\textrm{A} \times \textrm{N}} )_{\textrm{cutoff}}}$ as the $\textrm{A} \times \textrm{N}$ value that results in ${\langle {\upvarepsilon }\rangle _{\textrm{max}}} \approx 10^\circ$, we observe ${({\textrm{A} \times \textrm{N}} )_{\textrm{cutoff}}} = 18,900$, 4,200, and 2,000 for configs. 1, 2, and 3 respectively when ${\Delta }$ is known. When ${\Delta }$ is unknown ${({\textrm{A} \times \textrm{N}} )_{\textrm{cutoff}}} = 3,200$ for config. 3 and ${({\textrm{A} \times \textrm{N}} )_{\textrm{cutoff}}} > {10^5}$ for configs. 1 and 2.

For config. 3 with ${\Delta }$ unknown, ${\langle{\upvarepsilon }\rangle_{\textrm{max}}}$ converges to ${\sim} 6.1^\circ $ (rather than $0^\circ $, as expected) with increasing $\textrm{A} \times \textrm{N}$ values. We expect that this convergence is an artifact of our lookup table-based fitting algorithm, rather than an outcome of specific degeneracies. As such, future implementations of VALiD may be improved by alternate fitting methods. For example, it may be possible to derive analytical methods for fast best-fit calculations that are stable against noise. Alternatively, following the lookup table method with gradient descent optimization may help the fitting algorithm converge to the global optimum.

We next simulated ${\left\langle {\upvarepsilon } \right\rangle _{\textrm{max}}}$ as a function of $\textrm{A} \times \textrm{N}$ and ${\Delta }$ with config. 3, $\textrm{N} = 6$, $\textrm{A}$ unknown, and ${\Delta }$ known (Fig. 5(B)). For computational tractability, we reduced the number of Monte Carlo iterations per orientation to $n = 1,000$ and reduced the number of simulated orientations per condition to 120. Our results show that error is always negatively correlated with $\textrm{A} \times \textrm{N}$. As predicted, we also see that error increases monotonically as ${\Delta }$ approaches the magic angle of ∼$54.7^\circ $. Through parameter fitting to the $10^\circ $ contour line, we find that the relationship:

3.4 Accounting for non-conical distributions and orientation-dependence of emission

We close out our simulation-based analysis by evaluating the effect that lifting assumptions 1 and 2 has on orientation estimation. Specifically, we generated maps of ${\langle{\upvarepsilon }\rangle_{\textrm{max}}}$ as a function of ${\Delta }$ and $\textrm{A} \times \textrm{N}$ as described above for Fig. 5(B) under three conditions that we call “static cone”, “mobile Watson”, and “static Watson”. For these simulations, we used config. 3 with $\textrm{N} = 6$, ${\Delta }$ known, and $\textrm{A}$ unknown, and calculated ${\textrm{I}_{\textrm{ideal}}}$ using condition-specific methods as described below. We assess error when reference curves are calculated either with assumptions 1 and 2 (Eq. (12)) or with the same condition specific techniques used to calculate ${\textrm{I}_{\textrm{ideal}}}$. By running simulations in this way, we assess problems that may arise if Eq. (12) is used in situations where assumptions 1 and 2 do not hold true, as well as the extent to which these problems can be corrected by using appropriate calculation methods for orientation fitting.

The static cone condition relaxes assumption 2, which states that fluorophore orientation randomizes within the probability distribution between excitation and emission. Instead, the opposite situation, wherein fluorophore orientation does not change at all between excitation and emission, is made. Under this condition, Eq. (12) no longer holds true because the orientation-dependence of emission will also modulate intensity. The mobile Watson condition relaxes assumption 1 and assumes that fluorophores are distributed according to the Watson distribution, rather than a conical distribution, as described in section 2.2. The static Watson condition relaxes both assumptions 1 and 2 by assuming that fluorophores are both distributed according to the Watson distribution and do not re-orient between excitation and emission.

To develop a means of calculating ${\textrm{I}_{\textrm{ideal}}}$ for the static cone condition, we first calculated the expected collection efficiency (the fraction of emitted photons collected by the microscope) for a static fluorophore as a function of ${\uptheta }$. For these calculations, we accounted for 1) the anisotropic nature of fluorescence emission, wherein the intensity of fluorescence emission is strongest in directions perpendicular to ${\bf {\upmu} }$ and decreases in a ${\sin ^2}{\; }$fashion to 0 in directions parallel to ${\bf {\upmu} }$, 2) the limited collection angle of a microscope objective, and 3) partial reflection of emitted light at the sample-surface interface (the polarization of all emitted rays is perpendicular to the ray and parallel to the plane made by the ray and ${\bf {\upmu} }$). We also assumed that fluorophores were far enough from the sample-coverslip interface to prevent near-field effects on emission. Running these calculations we found that the collection efficiency of a fluorophore ${\Gamma }$ was well described by the equation:

 

Fig. 6. Surface plots showing ${\langle{\upvarepsilon }\rangle_{\textrm{max}}}$ as a function of ${\Delta }$ and $\textrm{A} \times \textrm{N}$ when the underlying ensemble can be represented by (A,B) a cone of fluorophores that do not reorient between excitation and emission, (C,D) a Watson distribution wherein fluorophores re-orient between excitation and emission, and (E,F) a Watson distribution wherein fluorophores do not re-orient between excitation and emission. All simulations were performed with config. 3, $N = 6$, ${\Delta }$ known, and $\textrm{A}$ unknown. Simulations were performed with reference intensity values (A, C, E) calculated using the static cone assumption (Eq. (12)), or (B, D, F) corrected in a manner that is specific to the true underlying distribution. Black contour lines are specific to each surface plot, while red dashed contour lines are the same contour lines that are present in Fig. 5(B) (static cone).

Download Full Size | PDF

We found that for the static cone condition, ${\langle{\upvarepsilon }\rangle_{\textrm{max}}}$ increased by a few degrees for all ${\Delta }$ and $\textrm{A} \times \textrm{N}$ when the lookup table fit values were calculated assuming a mobile cone (using Eq. (12)). This increase in error can be visualized by an outward shift of the contour lines in Fig. 6(A). Furthermore, error appeared to converge to non-zero values, which suggests that improper application of assumption 2 can result in small systematic errors. Systematic errors were eliminated and ${\langle{\upvarepsilon }\rangle_{\textrm{max}}}$ reduced back to the baseline level when reference curves were properly calculated using the static cone model (Fig. 6(B)).

Similarly, the mobile Watson condition demonstrated increased ${\langle{\upvarepsilon }\rangle_{\textrm{max}}}$ and introduced systematic errors when lookup table values were calculated using a mobile cone assumption (Fig. 6(C)). The increase in error was more pronounced than in the static cone case and exhibited a non-monotonic dependence on $|{{\Delta } - 54.7^\circ } |$. Again, error reduced to baseline when lookup table values were calculated using Eq. (16), which corresponds to the mobile Watson condition (Fig. 6(D)). The same trends were observed for the static Watson condition, with the most pronounced pre-correction increases in error of all three conditions (Fig. 6(E),(F)). Interestingly, in the Watson distribution cases contour lines converged to ${\Delta } \approx 60^\circ $, rather than the magic angle, when corrected reference values were used (Fig. 6(D),(F)).

Together, these results show that failure to properly parameterize the fluorescent ensemble under investigation may slightly reduce the accuracy of the ensemble. However, under all three of the conditions tested systematic errors caused by improper parameterization were generally minimal ($< 5^\circ $); the $3^\circ $ contour line substantially deviated from the static cone condition (Fig. 5(B)), but the $10^\circ $ contour line did not. As such, we conclude that the effects of assumptions 1 and 2 on performance are likely minimal. While efforts should generally be taken to use accurate models for calculations of ${\textrm{I}_{\textrm{ideal}}}$, in many cases the amount of information about a fluorescent ensemble is limited. In such cases, Eq. (12) can be used as a reasonable approximation for simple axially symmetric fluorescent ensembles.

4. Discussion and conclusions

Here we present an initial theoretical and computational exploration of a proposed technique, which we call VALiD, to perform non-degenerate measurement of fluorescent ensemble orientation. We present a mathematical model that can be used to fit the orientation of fluorescent ensembles to sets of intensity measurements acquired using the VALiD approach. We present three potential microscope configurations and optimal experimental parameters for each. We also present Monte Carlo simulation-based estimates for the SNR level necessary to implement VALiD with high accuracy.

We chose to use this Monte Carlo simulation-based approach over analytical approaches for several reasons. First, the Monte Carlo simulation algorithm described is simple and does not rely heavily on analytical formulations. This simplicity could, in the future, enable extension of the method to simulate more complex fluorophore orientation distributions that are not easy to describe analytically. Second, the least-squares residuals process used to calculate orientation from simulated intensity measurements resembles the process that would commonly be used to process experimental data, thus allowing for an unbiased means of assessing the predicted accuracy of VALiD under different conditions. Third, this approach enables the detection of systematic errors arising from measurement degeneracies between orientations that are not similar to each other. Nonetheless, future studies would likely benefit greatly in terms of computational speed and fundamental insight from an analytical treatment. For example, analysis of axially-symmetric distributions in terms of their second moment matrices [24,39] could enable rapid matrix-based (rather than lookup table-based) orientation calculation. Similarly, Cramer-Rao lower bound estimates may yield fundamental insight into the maximal expected accuracy of VALiD under various conditions and prove crucial to optimization efforts [16,38].

The three tested configurations are only a few examples of microscope setups that can be used to implement VALiD. For example, variants of VALiD could also be implemented using light-sheet, structured illumination, and variable azimuthal incidence angle (rather than variable inclination) microscope setups. Future efforts to implement VALiD would benefit from paired error-estimation simulations, similar to those presented here, that consider specific aspects of the microscope system and account for effects such as ${\upbeta }$- and ${\upalpha }$-dependent reflection of the excitation laser at the sample-surface interface. Such simulations could help experimenters optimize their experimental apparatus to maximize performance.

Practically, it will be necessary to correct for changes in the shape of the illumination profile (i.e. by normalizing images to control images of randomly oriented fluorophores taken at corresponding ${\upalpha }$ and ${\upbeta }$ values). Normalization should also account for the z-dependence of excitation beam intensity; the lateral component of the beam’s propagation with nonzero ${\upbeta }$ will result in a z-dependence of $|{\textbf E} |$, particularly when the beam has a non-flat (e.g. Gaussian) profile (although flat field illumination [40] could mitigate some of these issues). It may also be necessary to account for slight changes in the ${\textbf E}$ orientation of excitation light passing through the coverslip-sample interface due to reflection at the interface, which is polarization-dependent and can be accounted for using the Fresnel equations.

Note that a setup similar to the one we simulated is used for though-the-objective TIRF microscopy, meaning that VALiD can be implemented on many TIRF microscopes with slight modifications. Variable inclination microscopy techniques have been increasingly used in recent years, including: variable-angle TIRF [41–43] and scanning angle interference microscopy [44], both of which enable enhanced z-axis localization resolution; highly-inclined thin illumination microscopy (HiLo) [45], which enhances signal-to-noise by restricting illumination to a thin slice; variable-angle near-TIRF [46], which enables refractive index image mapping; and three-dimensional optical polarization tomography, which enables orientation mapping of single static fluorophores [17]. Accordingly, the proposed design is experimentally feasible.

While we have considered four types of ensembles (mobile/static, cone/Watson) and found that the static cone assumption is somewhat reasonable for all of them, physical scenarios will likely exhibit intermediate behavior such as partial fluorophore re-orientation between excitation and emission or behave as wobbly cones. In some cases, optimization routines that incorporate additional parameters may help to account for some of these properties. Constrained optimization routines (e.g. treating ${\Delta }$ as unknown but restricted to the range between $0^\circ $ and $30^\circ $) and/or combination with emission-splitting optical components may offer improvements in accuracy and yield additional information about the type of fluorescent ensemble that is being investigated. We have not explored the analysis of non-axially symmetric orientation distributions, axially symmetric distributions with less simple (e.g. multi-modal) shapes, or axially-symmetric distributions with separate static and temporal components to their orientational heterogeneity [47–49]. Future theoretical investigations, and potentially combination with techniques that enhance orientational resolution (such as two-photon excitation microscopy [50,51]) may enable resolution of arbitrarily complex orientation distributions. Combination of VALiD with computational microscopy [52] or structured illumination [53] approaches may enable the reconstruction of super-resolution images from VALiD acquisitions. However, fluorophores of certain orientations will likely produce systematic biases in the observed spatial illumination pattern, particularly when they are out-of-focus [54]. Accurately reconstructing such images with high spatial resolution may require joint spatio-angular deconvolution procedures [55,56].

In conclusion, VALiD overcomes a fundamental limitation of existing fluorophore ensemble orientation mapping techniques and is therefore a promising technique for fluorescence microscopy orientation mapping applications.