Topics:

• Visualisation: We will look at some microscope recordings of fluorescent calcium activity in the brain of larval zebrafish in ImageJ
• Classification: We will identify auditory responsive neurons based on their functional activity profile in Matlab

Fluorescent calcium imaging is a neuroimaging technique in animal research which involves using a calcium-sensitive fluorophore as a proxy for neuronal activation. In this dataset we are using genetically encoded GCaMP6s, which is based on the fluorophore GFP (green fluorescent protein). The modification for GCaMP6s is that the fluorophore is weakly fluorescent in the absence of Ca²⁺ ions, but changes configuration when it binds to Ca²⁺ ions and fluoresces much more strongly. Since Ca²⁺ ions flood into the neuronal cell body when it fires action potentials, this fluorescent marker of Ca²⁺ concentration gives a readout of neuronal activity levels.

The data you will work from in this tutorial comes from larval zebrafish, which are particularly well suited to the calcium imaging technique as they are optically transparent and have relatively small brains. We can therefore capture functional neuronal activity at single-neuron resolution throughout the brain, which is not possible in any other vertebrate. The imaging was performed on a custom light-sheet microscope which moves the excitation plane through the z-dimension of the brain while adjusting the acquisition focal plane with an electrically tuneable lens. By scanning through the brain like this, we acquire volumetric functional imaging of neural activity. We can also present stimuli to the fish to see the neural responses, in the data you will work with today they heard a sequence of white noise bursts.

Optional further reading for your interest:

Custom light-sheet microscope

Suite2p for segmentation

Auditory responses in larval zebrafish: Poulsen et al. , Wilde et al.

• [6 pts] Classification of extracted neurons

Data for this tutorial can be found here. If the image files are to large to download, the data without images can be found here.

Download the FIJI ImageJ application (open source). Open the calcium_imaging_raw.tif file in FIJI (you can use click and drag).

At the bottom of the window you can scroll through the frames of the acquisition. Re-shape the stack using Image→Hyperstacks→Stack to Hyperstack and use the below settings. Now you can scroll through depth and time independently.

1. Adjust the display brightness and contrast. Open the B&C window through Image→Adjust→Brightness/Contrast. The aim here is to choose levels that allow you to visualise fluctuations in fluorescence, without over-saturating or darkening the image.
2. Identify functional populations. Scroll through the images within the tiff, and find some neuronal populations with interesting activity. Use one of the selection tools to draw a circle around your chosen population (as shown in the example). Now extract the fluorescent activity over time within that region of interest using Image→Stacks→Plot Z-axis Profile. Select to plot over ‘time’. You can toggle on ‘live’ on the plot to see a blue line that tracks as you scroll through time on the image. Try selecting different populations in different parts of the brain. How similar or different is their activity?

Now we will look at the activity of some neurons that have been segmented and pre-processed for you using suite2p. This tool identifies regions of interest (ROIs) that represent individual neurons. It then calculates the average pixel value (fluorescence brightness) within that neuron at each time step to produce a raw activity trace of fluorescence. It also outputs the dimensions and some other statistics of each ROI.

In Matlab, load in the ‘fish_data.mat’ file. This contains: a mean image of one z-plane (‘mean_img’), statistics about the ROIs that are the segmented neurons (‘roi_stats’), the raw fluorescence traces (‘raw_fluor’: neurons x timepoints), the ΔF/F fluorescent traces (‘DeltaF’: neurons x timepoints), a theoretical trace of auditory responses to use as a regressor (‘aud_regressor’), and a smoothed motion-correction output to use as a motor regressor (‘motor_regressor’). The acquisition rate at a given z-plane is 2 frames per second.

The pre-processing involves taking the ΔF/F, which smooths and baseline-subtracts the fluorescent activity to account for bleaching over the course of the recording and differences in baseline fluorescent intensity between neurons. In this task, we will identify which of the segmented neurons respond to the auditory stimuli that were presented during the recording.

• HW4.1. [1 pt] Observe the locations of segmented ROIs on the mean image. First, let’s look at the mean image of the z-plane we are interested in, and where neurons were identified. Use the code provided to create this plot, and save it to add to your report. Note there are some gaps where neurons were not identified: what might contribute to a neuron failing to be segmented? Comment in your report.

% Plot the mean image from this z-plane
figure; hold on
imagesc(mean_img);
h=gca;
h.YDir='normal';
colormap('gray');
% Now add the locations of the identified ROIs (neurons)
for i=1:length(roi_stats)
    plot(roi_stats{1,i}.xpix,roi_stats{1,i}.ypix);
end

• HW4.2. [1 pt] Explore the raw fluorescence traces. Let’s see how the raw fluorescent activity traces look. Plot the mean raw fluorescence of all neurons over time and the raw activity traces of three example neurons using the code below. When we want to compare activity between neurons, or from the same neuron at different times in the recording, it is important that the baseline fluorescence is stable and similar between neurons. Look at the figures you have generated, are these criteria met? Include these figures and your assessment of the baseline activity in your report. Also include some suggestions of factors that could contribute to baseline fluorescence changing over time, and being different between neurons.

% plot the mean of raw fluorescence across all neurons
figure;
plot(mean(raw_fluor,1))
xticks(0:120:size(DeltaF,2))
xticklabels((0:120:size(DeltaF,2))/2)
xlabel('Time (seconds)')
ylabel('Mean raw fluorescence')
 
% plot raw fluorescence from some example neurons
ex_neurons=[1,2,3]; % choose a different 3 (or more) neurons amongst the 1216 available
figure; hold on
legend_labels={};
for i =1:length(ex_neurons)
    plot(raw_fluor(ex_neurons(i),:))
    legend_labels{i}=sprintf('neuron %g',ex_neurons(i));
end
xticks(0:120:size(DeltaF,2))
xticklabels((0:120:size(DeltaF,2))/2)
xlabel('Time (seconds)')
ylabel('Raw fluorescence')
legend(legend_labels)

• Try out baseline correcting and getting the change in fluorescence. Now we want to process the traces to achieve two things: get all neurons in the same frame of reference and achieve stable baseline values. Using the same example neurons you chose in the previous step, edit the code to plot (F – mean(F))/mean(F), where F is the raw fluorescence of a given neuron.

• HW4.3. [1 pt] Get the change in fluorescence using a moving baseline. As you will have seen, if the baseline is consistently changing then subtracting the mean does not flatten it. The approach we use to combat this is to use a moving baseline. This means at each time point the change in fluorescence is calculated using a smaller window of time on either side of that time point, rather than the whole recording. Write a function 'my_deltaf' to get the ΔF/F using a moving baseline (the window size can be 200 timepoints). Run it on your example neurons and plot the output. Include the plot and your function in your report. Comment on how successfully this approach has placed the traces in the same frame of reference and flattened the baseline.

• Create an exploratory rasterplot of the activity of all the neurons. For the rest of the analysis, we will use the traces that have been pre-processed for you (the DeltaF variable). We have over 1000 neurons segmented in this z-plane, let’s get an idea of what all their activity looks like over time. Use the code provided to generate a rasterplot. You will need to adjust the lower and upper limits of the colour bar to better observe the range of activity of the neurons. Save this plot such that the limits you chose are visible in the title and include it in your report.

% create a rasterplot of all the neurons in this plane
figure;
colormap('hot') % you may change the colourmap if you like
low_lim=; % adjust the value for the lower limit
up_lim=; % adjust the value for the upper limit
imagesc(DeltaF,[low_lim,up_lim]) 
colorbar()
ylabel('Neurons')
xticks(0:120:size(DeltaF,2))
xticklabels((0:120:size(DeltaF,2))/2)
xlabel('Time (seconds)')
title(sprintf('lower limit %g, upper limit %g',low_lim,up_lim))

• HW4.4. [0.5 pts] Run a linear regression to the auditory stimuli. Now we will use a theoretical trace (the ‘aud_regressor’ variable) to find neurons with an auditory-responsive activity profile. We will run a linear regression to this theoretical trace, and take neurons with a high r-squared value as being auditory responsive. Use the below code to run a linear regression on each of the neurons and generate an output variable r2 containing the r-squared value of the regression for each neuron. Then plot the distribution of r-squared values and choose a threshold cut-off to define auditory-responsive neurons. Save this plot and include it in your report.

% run a linear regression to the theoretical trace on each neuron
r2=[];
for i=1:size(DeltaF,1)
    mdl=fitlm(DeltaF(i,:),aud_regressor);
    r2(i)=mdl.Rsquared.Adjusted;
end
% now plot the distribution of r-squared values and choose a threshold
thresh=0; % set the threshold for auditory classification
figure; hold on
histogram(r2)
xline(thresh,'--','LineWidth',2)
xlabel('r^2 to auditory trace')
ylabel('Number of neurons')
title(sprintf('Auditory threshold = %g',thresh))

• HW4.5. [0.5 pts] Plot the mean activity trace of the auditory neurons. Now we can check whether the activity of neurons with r2 values above the threshold really looks auditory-responsive. Use the below code to generate a plot of the mean activity of your auditory-classified neurons, compared to the theoretical trace. Add in the mean activity of the non-auditory neurons, and include this plot in your report.

figure; hold on
% plot the theoretical trace for reference
plot(aud_regressor)
% plot of the mean activity of auditory-classified neurons
plot(mean(DeltaF(r2>thresh,:),1))
% add the mean activity of non-auditory neurons
%
xlabel('frames')
ylabel('DeltaF/F')
legend({'Auditory stimuli', 'Auditory neurons', 'Non-auditory neurons'})

• HW4.6. [1 pt] Plot the locations of the auditory neurons on the mean image. Now we can also check whether this seems like an anatomically cohesive cluster of neurons that have similar activity. Modify the code from earlier to plot the masks of only the auditory-responsive neurons on the mean image. Include this figure in your report. Based on your figures of the activity and location of your auditory-classified neurons, comment on whether your classification was successful.

• HW4.7. BONUS [1 pt] Find motor-responsive neurons. Typically, in this kind of experiment we need to control for other types of activity that are not the one we are interested in. For example, if the subject moves this can create artefacts in the image and will be represented by activity in neurons involved in motor circuitry. The ‘motor_regressor’ variable has been made from motion correction that was applied during the Suite2p segmentation. Following the same steps as we did above for the auditory neurons, identify neurons with motor-related activity. Include plots of the distribution of r2 values with annotated threshold and mean activity trace compared to the regressor in your report, and finally create a plot that shows the auditory and motor neuron populations on the mean image.