Correcting motion induced fluorescence artifacts in two-channel neural imaging.

Imaging neural activity in a behaving animal presents unique challenges in part because motion from an animal's movement creates artifacts in fluorescence intensity time-series that are difficult to distinguish from neural signals of interest. One approach to mitigating these artifacts is to image two channels simultaneously: one that captures an activity-dependent fluorophore, such as GCaMP, and another that captures an activity-independent fluorophore such as RFP. Because the activity-independent channel contains the same motion artifacts as the activity-dependent channel, but no neural signals, the two together can be used to identify and remove the artifacts. However, existing approaches for this correction, such as taking the ratio of the two channels, do not account for channel-independent noise in the measured fluorescence. Here, we present Two-channel Motion Artifact Correction (TMAC), a method which seeks to remove artifacts by specifying a generative model of the two channel fluorescence that incorporates motion artifact, neural activity, and noise. We use Bayesian inference to infer latent neural activity under this model, thus reducing the motion artifact present in the measured fluorescence traces. We further present a novel method for evaluating ground-truth performance of motion correction algorithms by comparing the decodability of behavior from two types of neural recordings; a recording that had both an activity-dependent fluorophore and an activity-independent fluorophore (GCaMP and RFP) and a recording where both fluorophores were activity-independent (GFP and RFP). A successful motion correction method should decode behavior from the first type of recording, but not the second. We use this metric to systematically compare five models for removing motion artifacts from fluorescent time traces. We decode locomotion from a GCaMP expressing animal 20x more accurately on average than from control when using TMAC inferred activity and outperforms all other methods of motion correction tested, the best of which were ~8x more accurate than control.

This is a PLOS Computational Biology Methods paper.

Introduction

Population fluorescent imaging of calcium-sensitive fluorescent indicators is a powerful approach for recording cellular neural dynamics [1,2]. Calcium imaging’s widespread adoption has benefited from extensive development of computational algorithms to find and segment neurons of interest [3], to extract and denoise calcium signals [4,5], and to infer their underlying voltage signals [6]. An important goal of systems neuroscience is to probe the neural basis of animal behavior [7,8] and calcium imaging has been used to measure neural activity in awake and behaving animals [9-11]. However, animal motion during calcium imaging poses unique challenges both for segmenting and tracking neurons and for accurately extracting calcium traces. Many computational approaches have been proposed to register, segment and track neurons in the presence of motion [9,12-21]. These approaches account for the gross movement of a neuron relative to its neighbors or within the field of view. By contrast, there have been relatively few efforts to account and correct for motion-related changes to the extracted fluorescence intensity time-series itself. Accounting for these latter effects of motion, however, is particularly important because motion related changes to fluorescence intensity can appear similar to behavior-related neural signals of interest. Here, we construct a method called Two-channel Motion Artifact Correction (TMAC) which uses the fluorescence traces from two-channel imaging to correct for shared motion artifacts between the two channels.

When a neuron deforms or moves relative to an imaging plane, its fluorescence intensity can change unrelated to neural activity (. While multiple factors could contribute to these motion induced changes, at least some of them arise from subtleties of image acquisition and segmentation. Even though the number of fluorescent molecules in a neuron is constant, the neuron’s shape, orientation and position are not. Accurately measuring fluorescence therefore requires carefully accounting for which voxels a neuron occupies, the intensity of each voxel, and the shape of the imaging hardware’s point-spread function used to calculate the voxel. For example, as a neuron changes its shape with respect to the point spread function, it may cover less voxels but be brighter in the voxels that remain. In real-world imaging conditions it can be challenging to detect subtle changes to the boundary of a neuron, especially in densely labeled neural populations, and similarly the point-spread function is not always sufficiently characterized. These challenges are particularly acute in recordings of moving animals, such as C. elegans [11,22,23], Hydra [24], Drosophila larvae [25,26], and zebrafish larvae [27,28] that all exhibit large head deformations during movement.

Behavior induces motion artifacts in optical recordings.

A) Diagram of a moving C. elegans undergoing motion and deformation during simultaneous measurements of GCaMP and RFP fluorescence. Deformations of the cell can conspire with imperfect segmentation to create motion artifacts. RFP intensity fluctuations reflect motion or noise. GCaMP intensity fluctuations reflect motion, calcium activity, and noise. B) Top, Animal body curvature. Bottom fluorescence of a neuron from a whole-brain recording of a moving worm expressing both GFP and RFP. Both fluorophores are activity-independent, yet we observe large highly correlated fluctuations in the two-channel fluorescence. C) Top, Body curvature. Bottom, fluorescence of a single neuron taken from a whole-brain recording in a moving worm expressing both GCaMP and RFP. The two channels still show correlated fluctuations despite the activity dependence of GCaMP. D) Histogram of the Pearson correlation coefficient squared between the red and the green channel for a dataset of whole-brain recordings from 10 moving GFP, RFP control worms. E) Same as D but in 9 GCaMP, RFP worms.

Because these artifacts arise from motion, they have the potential to confound interpretations of the neural correlates of behavior. For example, if a neuron appears to change fluorescence intensity during head bending, it may be hard to disambiguate whether that fluorescent transient reflects calcium activity or motion artifact. Two-channel imaging offers one strategy to account for motion-induced changes in fluorescence by measuring a calcium indicator like GCaMP in one channel and a calcium-independent indicator like RFP in another. Any intensity fluctuations of the calcium-independent channel must be either noise or motion artifact. In principle, this knowledge can be used to account for and correct artifacts in the calcium imaging channel. Crucially, whatever the source of these artifacts, as long as they affect the red and green channel equally, two-channel motion correction should be able to account for and remove them.

Two-channel calcium imaging has its origins in the use of FRET based calcium indicators that use a donor and receptor and reports calcium signals as the ratio between the two channels [29]. The ratiometric approach further proved useful in freely moving animals [10], because the ratio is less sensitive to motion-induced fluorescent changes common to both channels. Single-channel GCaMP based indicators have now surpassed FRET-based indicators in terms of popularity because of their speed and brightness [30]. However, the use of two-channel imaging has persisted in moving animals in part for its ability to remove motion artifacts with methods such as taking the ratio of the two channels [22,23,31-33], or by performing linear regression [11,34]. To date there has not been a systematic comparison of different mathematical approaches to account for motion induced changes to fluorescent intensity in two-channel imaging. TMAC explicitly models the statistical distributions of the noise, artifact, and cell activity and uses Bayesian inference to find activity uncontaminated by motion artifact. We use experimental data to compare TMAC with four other approaches to motion artifact removal and demonstrate that it outperforms previous methods.

Results

Motion artifacts

We inspected two-channel calcium imaging recordings of neurons in moving C. elegans from Hallinen et al., 2021 () for signs of motion-related calcium transients. We considered two types of recordings: those of control animals that express only the calcium-insensitive fluorophores GFP and RFP and those of calcium-sensitive animals that expressed the calcium indicator GCaMP in addition to RFP. In each case, GFP or GCaMP was imaged in the green channel and RFP was imaged in the red channel. In order to account for differences in intensity, we divide each channel by its time average and report fluorescence as fold change from the mean. In control animals, the green and red fluorescent intensity had large fluctuations relative to the mean ( and were highly correlated to each other (). Because these animals contained no neural-related signals and only motion, we conclude that motion artifacts contribute substantially to fluctuations in fluorescence. We also note the strong correlation between the GFP and RFP fluorescence, which suggests that the motion artifact is shared between the two channels. These shared fluctuations are what make motion correction using a second channel possible.

In GCaMP recordings, we also observed a correlation between the green and red channels (, although less so than in control recordings. The presence of activity-related signal in the green channel likely explains this difference. The correlation that remains between the red and the green channel suggests that the red channel could be used to correct for shared motion artifacts in the green channel.

Model

To account for fluctuations in fluorescence induced by motion artifact, we propose a method called Two-channel Motion Artifact Correction (TMAC) (). This method relies on a latent variable model of calcium fluorescence, which we use to infer the latent neural activity of each neuron. Intuitively, inferring the latent activity from the model amounts to subtracting the motion signals present in the red channel from the green channel while accounting for channel independent noise. We first process the fluorescence data by dividing each channel by its time average, such that the processed data has mean 1 and units of fold change from the mean. By specifying our model fluorescence as fold change from the mean the parameters do not depend strongly on experimental conditions such as cell expression, laser intensity, or image acquisition time, allowing us to compare across different neurons and fluorophores ().

TMAC infers activity and hyperparameters from synthetic data.

A) Diagram of the structure of TMAC. The green channel is modeled as the sum of the calcium activity, motion artifact, and independent gaussian noise. The red channel is modeled as the sum of the motion artifact and independent gaussian noise. The motion artifact is shared between the two channels. B) Top: fluorescence from a synthetic green and red channel. Middle and bottom: inferred activity and motion compared with the true activity and motion. C) Correlation squared between estimated activity from TMAC and true activity over many synthetic datasets. D) Violin plot of inferred and true parameter values when fitting TMAC.

We define the fluorescence from a neural activity-dependent green channel and a neural activity-independent red channel as follows. Let and denote vectors of the preprocessed red channel and green channel fluorescence from the same neuron. We assume that is the sum of a latent motion artifact time series and additive Gaussian white noise . Similarly, we assume the green channel measurements , arise as the sum of the same latent motion artifact , a latent time series of neural activity , and Gaussian white noise . To obtain a complete generative model, we assume Gaussian process priors over and , which specify that they evolve smoothly over time. Formally the model can be written as: Where 1 denotes a vector of 1s, denotes a multivariate Gaussian distribution, and σ and σ denote the variance of the additive Gaussian noise in the red and green channels, respectively. The latent time series of the motion artifact and neural activity are both governed by Gaussian process priors, with prior covariances and . The prior covariances and are each parameterized by a pair of hyperparameters, (σ, τ) and (σ, τ), where σ denotes the prior variance or amplitude, and τ denotes the temporal length scale, controlling the degree of smoothness. We use a stationary radial basis function (RBF) covariance function, such that the covariance between two time points depends only on the temporal separation between them:

Given this model structure, data and , parameters and , and hyperparameters θ (σ, τ, σ, τ, σ, σ), we fit the model in two steps. We first optimize the hyperparameters θ by maximizing the marginal likelihood p(, |θ). We then compute the maximum a posteriori (MAP) estimates for and given the optimized hyperparameters by maximizing the posterior p(,|,, ).

The marginal likelihood is obtained by integrating the joint distribution of data and latents over and , which can be computed in closed form:

This quantity is the likelihood function for the hyperparameters, which can be optimized numerically to provide a maximum (marginal) likelihood estimate:

By Bayes Rule the posterior is proportional to the product of the likelihood function for and and the prior : For this model, where prior and likelihood are both Gaussian, the posterior distribution is also Gaussian, with mean and covariance that can be expressed in closed form: where the covariance Λ and mean , which corresponds to the MAP estimates for and are given by: Note that our software implementation of TMAC relies on a Fourier-domain representation which diagonalizes the prior covariances, substantially reducing computation time [35].

Model validation

In order to demonstrate that our model correctly infers ground truth activity and variances, we used TMAC to generate a synthetic dataset ( top) with 5000 time points and a correlation time scale which roughly matches experimental datasets from , such that 6 time points approximately corresponds to one second. We then used TMAC to infer both the activity and motion artifact from this synthetic dataset ( middle and bottom). We find a good correspondence between the inferred activity from TMAC and the true activity (). The model also accurately infers the hyperparameters of the model ().

Decoding Behavior

We next sought to evaluate TMAC on experimentally acquired neural population calcium imaging recordings in moving worms from [11]. As in [11] the data was preprocessed to remove missing values by linearly interpolating over time, and photobleaching was corrected by separately dividing green and red channel fluorescence by an exponential fit to the fluorescence with the decay time constant shared between neurons. Datasets with a high number of missing values were excluded ().

The hyperparameters learned by TMAC provide a convenient method to characterize the recordings. We identified a putative high signal-to-noise neuron, by finding a neuron that the model associated with a high ratio of activity variance (σ) to motion and noise variances (σ) ( middle). We also identified a putative low signal-to-noise neuron by finding a neuron that the model associated with high ratio of motion artifact (σ) to activity and noise (σ) ( bottom). For completeness we also provide example fluorescence traces from neurons recorded in an immobilized worm (). In the case where TMAC estimates high signal-to-noise, TMAC predicts that a smoothed version of the green channel represents true calcium activity ( middle). In the case where TMAC estimates that most fluorescent fluctuations are due to motion, TMAC deviates from the measured fluorescent intensities and returns a flatter inferred activity ( bottom).

We next wanted to evaluate TMAC’s performance on removing motion artifacts. Since we lack access to the ground truth calcium activity, we considered the problem of decoding animal curvature from neural activity (). To decode the behavior of each animal we followed the same procedure from [11]. Namely, we trained a linear decoder with an L2 penalty on the weights determined using multifold cross validation. Behavior was predicted at each time point from the motion corrected activity, and all reported decoding accuracy scores are from data held out from the center 40% of the recording ().

Using decodability as a metric for motion correction performance is challenging because motion artifacts may also provide information about the target behavior. In some regimes, successfully removing motion artifacts could in principle reduce decodability rather than improve it. We therefore evaluated motion correction performance by inspecting the decodability of two types of recordings; one that had both activity-dependent and activity-independent fluorophores (GCaMP and RFP), and a control worm that had two activity-independent fluorophores (GFP and RFP), both from Hallinen et al., 2021. Although the GFP control worm contains no neural signal, the worm’s behavior can still be decoded, albeit poorly, by relying on behavior-related information in the motion artifact. We reasoned that a successful motion correction algorithm should reduce the decodability of motion artifacts while retaining the ability to decode from neural signals. We therefore define a new metric that evaluates motion correction performance as the ratio of its decodability of recordings from GCaMP animals to its decodability of recordings from GFP animals ().

Using this metric, we compare five models used for motion correction in two-channel imaging. The models are (): single channel GCaMP which is just the green channel fluorescence with no motion correction, the ratio model which is the green channel divided by the red channel, a linear regression model [11], an ICA based approach [36], and TMAC. The linear regression model finds the best linear fit of the red channel to the green channel and then subtracts off that best fit from the green channel. The ICA approach performs ICA using the red and green channel as inputs and returns the independent component least correlated to the red fluorescence.

The activity fit by TMAC is ~20x more effective at predicting curvature from a GCaMP expressing worm than a control worm (Figs and S2). We calculated the correlation between the activity inferred by TMAC and the red channel and they did not correlate strongly, suggesting that the model is successful at removing motion artifacts common to the two channels (). For completeness, we also compare each method on our synthetic dataset ().

Discussion

Motion artifacts are prevalent in recordings of behaving animals (). These artifacts can appear as signals of interest, reducing the interpretability of the data. Here we presented TMAC, a model which infers the latent neural activity without these artifacts by leveraging information in two-channel calcium imaging recordings. We demonstrate that TMAC substantially reduces decodable motion artifacts in experimental data () and outperforms four alternatives ().

The model we propose uses an additive interaction between the fluorescence transients caused by motion artifact and the fluorescence transients caused by neural activity. It is unclear in experimental data what the true interaction between artifact and activity is. However, even if the true interaction is not additive, TMAC will perform well under model mismatch because it can be thought of as a linear approximation of the true interaction between a and m. Consider the case where the activity and motion interact multiplicatively. Multiplicative interactions have linear components when the two variables have nonzero means. The linear component of the multiplicative interaction is clear by writing the equation for the green channel with no noise. We can be sure that the and have a mean (in this case 1), because fluorescent data is nonnegative. Furthermore, the term itself can be approximated by + as long as and are small. If the approximation of and as Gaussian is valid, we can be sure that the and higher order terms are small, because the standard deviations of and must be small relative to 1 to avoid negative values. For these reasons, while TMAC is a linear approximation to the true interaction between and , it is an approximation that is highly accurate. An interesting avenue for future work would be to consider asymmetric distributions which can retain both high variability and positivity and may better approximate the true activity distribution.

In this work we demonstrated TMAC’s ability to remove decodable motion artifacts from calcium induced fluorescence. However, TMAC will remove motion artifacts in any type of two-channel imaging, so long as one channel is activity-independent and both channels share the same motion artifact component. TMAC could therefore be applied to a wide range of two-channel imaging modalities including for voltage imaging [37], fiber photometry when using an isosbestic wavelength [38,39], or two-channel two-photon imaging [40].

TMAC inferred neural activity from an immobilized animal.

A) GCaMP and RFP fluorescence from a neuron that TMAC estimates to have a high ratio of activity variance (σ) to motion and noise variances (σ), recorded from an immobilized worm. B) GCaMP and RFP fluorescence from a different neuron in that same recording that TMAC estimates to have a high ratio of motion variance (σ) to activity and noise variances (σ). Because the worm is immobilized, the motion artifacts are still small even for the highest motion variance neurons. When there is low motion artifact, TMAC estimates the activity is similar to a smoothed version of the green channel.

(PDF)

Click here for additional data file.

Motion artifact correction for GCaMP and GFP expressing animals.

A) Decoding accuracy when decoding whole-body curvature from GFP recordings with different motion correction methods applied. These animals do not express an activity-dependent fluorophore so all decoding comes from motion artifacts. The mean and median of decoding of each method is listed. B) Decoding accuracy when decoding whole-body curvature from GCaMP expressing animals. The ratio of B to the median for each method in A is the value reported in . C) As in Fig 3C, but each decoding value has been divided by the mean (rather than median) decoding values from each metric in A.

Fig 3

TMAC reduces decodable motion artifacts in experimental data.

A) Top, animal body curvature over time. Middle, GCaMP and RFP fluorescence from a neuron that TMAC estimates to have high signal to noise, recorded from a behaving worm. Bottom, GCaMP and RFP fluorescence from a different neuron in that same recording that TMAC estimates to have large motion artifacts. B) Time trace of animal curvature and predicted behavior, decoded from activity inferred by TMAC in a GCaMP worm. Gray shaded regions were used to train the decoder, white region was held out and used to evaluate decoding performance. C) Ratio of decoding accuracy (ρ2) when decoding GCaMP divided by the median accuracy for a GFP worm across different models (). D) Histogram over all neurons of correlation squared between RFP and activity inferred by TMAC from a GFP worm. RFP vs GFP data the same as in . E) Same as F but in a GCaMP worm.

(PDF)

Click here for additional data file.

Accuracy of motion correction methods on synthetic data.

Each of the 5 methods for motion correction were tested on the synthetic dataset from . The reported value is the distribution of correlation squared between inferred activity and true activity over instantiations of neurons. This synthetic data was generated from TMAC itself so it is unsurprising that it outperforms other methods on this dataset. The linear regression method also performs well because, like TMAC, it assumes an additive interaction between motion and activity.

(PDF)

Click here for additional data file.

Data source for each figure.

(DOCX)

Click here for additional data file.

Models for motion artifact correction.

(DOCX)

Click here for additional data file.