Here, we describe and evaluate a heuristic strategy that utilizes axon traces to detect and quantify boutons in 2PLSM images taken through a cranial window in vivo. This procedure consists of the following major steps: (1) tracing axons in 3D, (2) optimization of traces, (3) generation of axon intensity profiles, (4) detection of putative boutons based on the profiles, (5) normalization of intensity profiles and calculation of bouton weights, and (6) matching putative boutons across time-lapse images.

Axons can be traced automatically or manually with various tools (Acciai et al., 2016; Parekh and Ascoli, 2013). In this study, high density of labeled axons (Figure 1A) precluded the possibility of automated tracing. Therefore, axons were traced manually by using NCTracer software (Chothani et al., 2011; Gala et al., 2014), and traces were optimized as described in Materials and methods. Following optimization, two intensity profiles were generated for each axon by convolving specifically designed filters with the image at all trace node positions and scaling the results to unit means. While various filters can be used to generate axon intensity profiles, in this study, we settled on the following two: (1) a modified, multi-scale Laplacian of Gaussian filter (LoG_xy) to detect putative boutons and (2) a fixed size Gaussian filter (G) to provide an estimate of axon intensity in the regions devoid of boutons (see Materials and methods for details). In the following we will refer to such inter-bouton regions as axon shaft.

Figure 2A and B show that distinct putative boutons can be identified as peaks in the LoG_xy profile plotted against node positions along the trace, ${\displaystyle I^{Lo{G}_{xy}}\left(s_i\right)}$. This is because the LoG_xy filter is designed to sharpen boundaries between boutons by suppressing intensity in the regions immediately adjacent to boutons. In contrast, the G filter yields a smoother profile, ${\displaystyle I^G\left(s_i\right)}$, which is not very useful for resolving putative boutons that are in close proximity (arrows in Figure 2C), but is well suited for estimating shaft intensity. For these reasons, LoG_xy profiles were used to identify putative boutons, while G profiles were used to determine intensities of axon shafts.

Figure 2 with 2 supplements see all

Download asset Open asset

To automatically detect putative boutons in an LoG_xy profile we used an algorithm that is similar to the Backward-Stepwise Subset Selection method (Hastie et al., 2009). Here, a varying number of foreground peaks, ${\displaystyle N_f}$, and a constant number of background peaks, ${\displaystyle N_b}$, was fitted to ${\displaystyle I^{Lo{G}_{xy}}\left(s_i\right)}$ by minimizing the following objective function of peak positions, ${\displaystyle {\mu }_j^f}$and ${\displaystyle {\mu }_k^b}$, amplitudes, ${\displaystyle a_j^f}$and ${\displaystyle a_k^b}$, and widths, ${\displaystyle {\sigma }_j^f}$and ${\displaystyle {\sigma }_k^b}$ (see Materials and methods for details):

We note that, though profile background was fit with a sum of spatially distributed peak functions, constraints imposed on their widths (${\displaystyle {\sigma }_k^b\ge }$ 20 μm) ensure that the fit varies slowly along the axon (red line in Figure 2B).

Following the detection of putative boutons, intensity of bouton ${\displaystyle j}$ was defined as the sum of the ${\displaystyle j}$-th foreground peak amplitude and background intensity at the peak location:

Shaft intensity was estimated from the G profile by incorporating information about the positions of detected putative boutons, Figure 2C. To that end, we initialized the above described peak detection algorithm with the foreground and background peaks detected on the LoG_xy profile, but ran the algorithm on the G intensity profile. Shaft intensity for a given axon was defined as the fitted background intensity on the G profile, averaged over trace nodes,

In this expression, index G in the superscripts of peak amplitudes, ${\displaystyle a}$, positions, μ, and widths, σ, was added to emphasize that these quantities are calculated based on the G intensity profile. We note that background intensity, ${\displaystyle I_i^{Background}}$, is designed to vary smoothly along the axon (red line in Figure 2C) and be independent of bouton density, providing a robust estimate of axon shaft intensity.

Our goal is to use intensity profiles to extract structural information related to the physical sizes of boutons. This task is hindered by the facts that axon intensity depends strongly on expression levels of fluorescent molecules (Figure 1A) and microscopy conditions. Therefore, to measure unbiased structural information, intensity profiles must be properly normalized. One may consider using median (or mean) profile intensity or local shaft intensity for normalization. However, these types of normalizations can lead to errors. For example, if density of boutons varies across axons, normalization with the median (or mean) may bias boutons on higher bouton density axons towards lower intensity values. Also, if density of boutons is sufficiently large, normalization with local shaft intensity can lead to variability as the latter cannot be measured reliably between closely positioned boutons.

Our heuristic normalization approach is based on the idea that by convolving the LoG_xy or G filter with the image at a trace node position ${\displaystyle s_i}$ along axon ${\displaystyle a}$ in imaging session ${\displaystyle t}$, we obtain a quantity that is proportional to three factors: ${\displaystyle M_t}$, a factor related to imaging conditions [e.g. laser power, photomultiplier tube (PMT) voltage, and cranial window quality], ${\displaystyle {\rho }_{a,t}}$, volume density of fluorescent molecules, and ${\displaystyle A_{a,t}(s_i)}$, a structural factor which has been linked to axon cross-section area (drawn perpendicular to the xy projection of the axon centerline) convolved with the microscope point-spread function (Song et al., 2016) and profile filter:

In creating the LoG_xy and G profiles (Figure 2) we rescale ${\displaystyle {\tilde{I}}_{a,t}^{Lo{G}_{xy}}\left(s_i\right)}$ and ${\displaystyle {\tilde{I}}_{a,t}^G\left(s_i\right)}$ to unit means in order to minimize effects related to imaging conditions and expression levels, thus isolating structural information:

It may be tempting to use putative bouton intensity, ${\displaystyle I_j^{Bouton}}$ in Equation (3), which is detected based on ${\displaystyle I_{a,t}^{Lo{G}_{xy}}\left(s_i\right)}$ as proxy for bouton size. However, this may lead to bias as the denominators in Equation (6) depend strongly on bouton density. To address this issue, we use axon shaft intensity detected from G intensity profiles, ${\displaystyle I^{Shaft}}$ in Equation (4), for normalization. The resulting quantity, referred to as bouton weight, ${\displaystyle w_j^{Bouton}}$, conveys structural information, which is effectively independent of the above-mentioned bias,

Correlative light and electron microscopy (CLEM) was used to validate the described bouton detection procedure (Figure 3 and Figure 3—video 1). Four axon segments were selected for this analysis (Figure 3A). The axon segments were imaged in vivo with 2PLSM, the brain tissue was fixed shortly after, and subsequently imaged with EM (see Materials and methods for details). Putative boutons in the 2PLSM stack of images were detected and quantified as described above (Figure 3B). Axons in EM images were reconstructed in 3D (Figure 3C) and rendered in Blender software for further analysis.

Figure 3 with 1 supplement see all

Download asset Open asset

Putative boutons detected in 2PLSM images could be unambiguously matched with varicosities identified in EM (Figure 3—video 1, asterisks in Figure 4A–D, and Table 1). We used positions of the centers of matched boutons to register LM traces and EM centerlines with an optimal linear transformation (Hastie et al., 2009). Average distance between the registered traces (Gala et al., 2014) was small, 0.34 ± 0.17 μm (mean ± s.d.), confirming that the same set of axons was reconstructed in LM and EM. A small discrepancy between registered traces may be attributed to non-linear distortion of tissue in the EM experiment and the relatively low z-resolution of LM. Overall, 16 boutons en passant and one bouton terminaux could be identified in EM. The bouton terminaux was detected with our method, but was excluded from further analysis because it is not located on the axon centerline. With the omission of this bouton, the LM-based procedure detected 18 putative boutons with 0 false negatives and two false positives, both of which were small (1.1 and 2.2 in weight, #18 and #19 in Table 1).

Figure 4 with 1 supplement see all

Download asset Open asset

We tested the idea that 2PLSM intensity is correlated to the area of axon cross-section drawn perpendicular to the xy projection of axon centerline (Song et al., 2016). We refer to such cross-sections as z cross-sections (third row in Figure 4A–D) to distinguish them from normal cross-sections, which are drawn perpendicular to the axon centerline. For this, we calculated the normalized axon intensity profile, ${\displaystyle I_{a,t}^{Lo{G}_{xy}}\left(s_i\right)/I^{Shaft}}$, and compared it to EM z cross-section areas including and excluding mitochondrial z cross-section areas. Both z cross-section areas appear to be well correlated with normalized intensity profiles. However, it is clear from visual inspection that the normalized intensity profiles do not resolve small changes in axon cross-section, which is the result of limited resolution of 2PLSM.

To examine the extent to which LM-based measurements provide information about bouton size, we plotted bouton weight against bouton volumes including (Figure 4F) and excluding (Figure 4G) mitochondrial volume. Normal axon cross-sections were used to identify bouton boundaries and calculate bouton volume as described in Materials and methods (fourth row in Figure 4A–D and Figure 4E). Bouton #7 (red box in Figure 4B) was excluded from these analyses, because it is bounded on one side by the terminal bouton branch which biases weight measurement. The results show high degree of correlation in both cases (Pearson’s r = 0.93), supporting the idea that LM-based measurements can be used to quantify volumes of even very small varicosities (#3, 0.093 μm³). The results also provide support for the choice of filters, Equation (1), and the normalization procedure, Equation (7), showing that the proposed method could overcome axon-specific differences in the expression levels and bouton density, capturing meaningful fine-scale structural information.

It is important to emphasize that we did not attempt to find the best procedure for fitting the EM data. It is possible that some other combinations of filter types, normalizations, and parameters would lead to marginally better correlations. However, due to the small sample size, this would likely be a case of over-fitting. The described procedure was designed based on theoretical considerations, and the parameters were chosen based on the observed range of bouton sizes (see Materials and methods). Nonetheless, we explored alternative filtering and normalization strategies. Figure 4—figure supplement 1 shows profiles derived from raw voxel intensities, and profiles obtained with mean, Gaussian, and median filters of different but fixed sizes. Each of these profiles was normalized by its median value. The figure makes it clear that some amount of filtering benefits bouton detection, as jagged unfiltered profiles can lead to false positives. Similarly, relatively large filters of fixed size (e.g. 5 × 5 × 5 and R = 2) are not well suited for bouton detection as they often suppress intensities of small boutons and merge closely positioned boutons. When the filter size is carefully tuned (3 × 3 × 3 and R = 1), resulting profiles can look similar to those of multi-scale LoG_xy, although higher baselines in these profiles can make peak detection more challenging. Furthermore, we examined the combined effect of filtering and normalization on bouton measurement by comparing the amplitudes of profile peaks to corresponding bouton volumes. Multi-scale LoG_xy profile normalized with shaft intensity as described above leads to the highest correlation compared to the best considered alternatives (Figure 4—figure supplement 1E).

Next, we sought to evaluate the effects of various imaging conditions on bouton detection and measurement. To that end, a set of fluorescently labeled axons was imaged seven times in a span of 80 min with different laser power (LP) and PMT voltage (see inset in Figure 5A). In addition, in condition E, a thin layer of agarose was applied to the cranial window to mimic deterioration of window quality which often accompanies long-term imaging experiments. In the following, we assume that there is negligible structural plasticity of boutons throughout the duration of this short-term imaging experiment, and therefore differences in bouton measurements can be attributed to the effects of imaging conditions and measurement uncertainties. The inset in Figure 5B shows the maximum intensity projections of an axon segment imaged in all seven conditions. This inset illustrates that increasing (decreasing) LP and/or PMT voltage results in an overall increase (decrease) in intensity, and therefore proper normalization procedure must be used to minimize the effects of such changes on bouton measurements.

Figure 5

Download asset Open asset

Putative boutons on 16 axon segments were detected in all conditions independently (400 putative boutons per condition on average). Custom software was used to match the same putative boutons across conditions. Putative boutons detected in condition A were chosen to be the gold standard, and precision/recall in bouton detection in the remaining conditions B-G were evaluated as functions of bouton weight (Figure 5A and B). The results show that for 95% of putative boutons of weight ${\displaystyle w}$ > 2.0 both precision and recall equal one in all conditions. This number of unambiguously detected putative boutons goes up to 99% for ${\displaystyle w}$ > 2.5 and 100% for ${\displaystyle w}$ > 3.1. For reference, ${\displaystyle w}$ = 2.0 corresponds to the weights of the two smallest varicosities identified in EM (#3 and #11 in Table 1), and the weight of the next smallest bouton (#17) is ${\displaystyle w}$ = 2.8. Therefore, all but very small boutons can be detected with our method with high confidence.

To examine the effects of imaging conditions on bouton weight, we tested a subset of 290 putative boutons that were detected and matched across all seven experiments. Figure 5C shows bouton weights in conditions B-G plotted against the corresponding weights in condition A. Although axon intensities in conditions B and C are drastically different from A (see inset in Figure 5B), the normalization procedure could correct for these differences (blue and red lines in Figure 5C). However, small (≈ 7%) but significant bias was present in conditions D-G, which may have been caused by bleaching due to prolonged imaging. Figure 5C also reveals a considerable amount of variability in weight measurements. This variability, reflected in the R² coefficients, was similar across all comparisons, including the imaging experiments performed under the same conditions (e.g. A and G, cyan points). Therefore, bouton weight measurements are accompanied with uncertainty which is inherent to the LM-based methodology. This uncertainty cannot be eliminated entirely, and hence, it must be explicitly incorporated into models that derive biological information from bouton measurements.

To model the variability observed in Figure 5C, we examined bouton weight differences for all pairs of imaging conditions, ${\displaystyle \mathrm{\Delta }w=w_1-w_2}$. Because such changes clearly depend on bouton size, we looked at the statistics of ${\displaystyle \mathrm{\Delta }w}$ in various intervals of average bouton weight, ${\displaystyle ⟨w⟩=\left(w_1+w_2\right)/2}$. Figure 5D shows the cumulative distribution functions (CDFs) of bouton weight differences in eight intervals of ${\displaystyle ⟨w⟩}$. These distributions are not significantly different from Gaussian distributions, and their variances are roughly proportional to the mean bouton weights, ${\displaystyle \mathrm{var}\left(\mathrm{\Delta }w\right)=\alpha ⟨w⟩}$, Figure 5E. Therefore, all CDFs shown in Figure 5D could be standardized by rescaling the weight differences as ${\displaystyle \mathrm{\Delta }w/\sqrt{\alpha ⟨w⟩}}$, leading to distributions that are statistically indistinguishable from the Standard Normal CDF (all p > 0.05, one-sample KS test).

Based on this result, we propose a statistical model in which measured bouton weight, ${\displaystyle w}$, is the sum of the true weight, ${\displaystyle w_0}$, and a noise term, δ. The latter is randomly drawn from a Gaussian distribution with variance proportional to the measured weight, ${\displaystyle \mathrm{var}\left(\delta \right)=\frac{1}{2}\alpha w}$:

Equation (8) quantifies uncertainties in bouton weight measurements, making it possible to define bouton presence probabilistically. To that end, we impose a threshold on true bouton weight, ${\displaystyle w_{threshold}}$, and refer to putative boutons with ${\displaystyle w_0}$ > ${\displaystyle w_{threshold}}$, as LM (light microscopy) boutons. In the following we set ${\displaystyle w_{threshold}}$ = 2.0, which is motivated by several considerations. First, this value equals twice the average normalized axon shaft intensity. Therefore, by using ${\displaystyle w_{threshold}}$ = 2.0 we are only including peaks that are substantially larger than the axon shaft intensity. Second, ${\displaystyle w}$ = 2.0 corresponds to the weights of the two smallest varicosities detected in EM (#3 and #11 in Table 1). Finally, detection of putative boutons becomes unreliable for ${\displaystyle w}$ < 2.0 (Figure 5A and B). We note that as an alternative to choosing a single threshold value to define LM boutons, one could vary the threshold in a certain range (e.g. 1–3) and report results as functions of this parameter.

The probability that a putative bouton of measured weight ${\displaystyle w}$ belongs to the category of LM boutons, ${\displaystyle P\left(bouton|w\right)}$, can be calculated based on the noise model of Equation (8):

In this expression, ${\displaystyle erf}$ denotes the error function. Shaded regions in Figure 5F illustrate these probabilities for two putative boutons of measured weights ${\displaystyle w}$ = 1.5 and ${\displaystyle w}$ = 3.0. Note that even when ${\displaystyle w}$ is less than ${\displaystyle w_{threshold}}$, there is a non-zero probability that the detected peak is an LM bouton. For large ${\displaystyle w}$ (e.g. greater than 3), this probability approaches unity, and the LM bouton definition becomes virtually deterministic.

LM bouton definition can be used to calculate the probabilities of bouton addition, elimination, potentiation, and depression based on the measured weights in two imaging sessions (initial and final), ${\displaystyle w_i}$ and ${\displaystyle w_f}$:

We would like to clarify that LM boutons may or may not correspond to varicosities seen in EM (Table 1), and the latter may not always be associated with postsynaptic densities (PSDs) and, thus, functional synapses (Shepherd and Harris, 1998). Therefore, the relationship between an LM bouton and a synapse is not deterministic and is likely to contain false-positives and false-negatives. However, the number of such errors is expected to decrease with increasing ${\displaystyle w}$. For example, Table 1 shows that all LM boutons of ${\displaystyle w}$ > 2.2 correspond to EM varicosities, and all LM boutons of ${\displaystyle w}$ > 2.9 are varicosities associated with PSDs.

Next, we set out to determine if the above described bouton detection procedure can be used to identify significant structural changes in long-term, in vivo imaging experiments. To that end, axons of GFP labeled neurons were imaged in superficial layers of barrel cortex in five mice. Imaging was performed at 4 day intervals over a 24 day period (seven imaging sessions). On average, 968 bouton sites were detected on 20 axon segments in each animal and imaging session. These sites were tracked over the duration of the experiment (Figure 6—figure supplement 1) to quantify bouton addition, elimination, and weight changes as compared to the initial state of the circuit. The short-term imaging experiment of Figure 5 was used as control. Here, conditions B-G were compared to condition A, and the results were pooled.

Figure 6A shows the fraction of added LM boutons over time, as compared to the first imaging session. Specifically, we only consider significant bouton addition events, i.e. those for which the probability of addition according to Equation (10) is greater than 0.95. Analogous plots for the fractions of significant bouton eliminations and significant bouton weight changes (potentiation and depression combined) are shown in Figure 6B and C respectively. As was expected, the fractions of significant changes in the short-term imaging experiment (red points in Figure 6) are at the chance levels (dashed red lines). The latter was obtained with a bootstrap procedure in which the weights of individual boutons were independently shuffled across conditions. In contrast, in the long-term imaging experiment, the fractions of significant bouton additions, eliminations, and weight changes are significantly larger than chance already by the second imaging session (after 4 days), and these fractions continue to increase with time, consistent with the idea of gradual modification of the circuit. These results illustrate that the described methodology can be used to quantify circuit changes, despite the numerous challenges associated with long-term, in vivo imaging experiments.

Figure 6 with 1 supplement see all

Download asset Open asset

All experiments were performed according to the guidelines of the Swiss Federal Act on Animal Protection and Swiss Animal Protection Ordinance. The ethics committee of the University of Geneva and the Cantonal Veterinary Office of Geneva, Switzerland (approval code GE/61/17) approved all experiments.

For development of the detection methodology we used data from mice that had been repeatedly imaged as part of a larger optical micro-stimulation and behavioral study (not yet published). Only the methodology that is relevant for the present study will be mentioned. Adeno-associated viral (AAV) vectors encoding floxed GFP (AAV2/9.CAG.flex.eGFP.WPRE.bGH2 [Upen]) were co-injected with an AAV vector encoding Cre (AAV2/9.hSynapsin.hGHintron.GFP-Cre.WPRE.SV40 [Upen]) at a ratio of 0.15 × 10⁹: 1 (genome copies: genome copies). This produced an average of about 150 GFP-expressing cells per animal, mostly in L2/3 and L5 of the barrel cortex. We exclusively imaged the processes of these neurons in the upper layers of the cortex. Hence, the source of the fluorescent signal in the green channel is dominated by cytosolic GFP.

Following the AAV injection, we implanted a 3 mm diameter glass window over the barrel cortex as previously described (Holtmaat et al., 2009). Imaging experiments were started 6 weeks after the virus injection. Imaging was performed under anaesthesia (0.5–1.5% isoflurane). A custom-built alignment system was used to ensure identical positioning of the mouse’s cranial window at different time points, and to avoid rotations relative to the axial dimension of the objective. The system was based on aligning a beam reflected by the cranial window through two apertures, which uniquely determines a line in three spatial dimensions. The apertures and the laser position were constant, and the only free variables were the cranial window position and rotation (i.e., the x, y, and z-positions, as well as the row and pitch needed to be identical between sessions to satisfy the passage of light through the two apertures). We used a custom built 2PLSM (Janelia Research Campus, model Non-MIMMS in vivo microscope) to acquire in vivo anatomical images of the anesthetized animals through the cranial window. We used the ScanImage software (Janelia Research Campus, Vidrio Technologies) (Pologruto et al., 2003) to record the 2PLSM images, and additional custom software to align anatomical structures present on the central image stack by using a red-green overlay between current and past images.

For imaging we used a 20x water immersion objective (NA 0.95, XLUMPlanFI, Olympus, Japan). Images were acquired at a voxel volume of ≈ 0.13 × 0.26 × 0.8 μm³ (x × y × z). We binned the data in the x dimension to generate isometric voxels in x and y for all analyses. The voxel dwell time was 0.8 μs for the full resolution images (and thus twice as high after 2x binning of the x dimension). Each ROI comprised an image stack of ≈ 270 × 270 × 250 μm³. As the imaging light source we used a Ti:sapphire femtosecond pulsed laser (Chameleon ultra II, Coherent) that was lasing at 1010 nm. Emitted fluorescence light was split into two channels using a dichroic mirror (Semrock, FF735-Di01−35.5 × 49.0) and two bandpass filters (red channel, Semrock FF01-607/70; green channel, Semrock FF01-530/70). Each channel was equipped with a PMT (red channel, Hamamatsu R3896; green channel, Hamamatsu H10770PA-40SEL). Images were analyzed only in the green channel.

Electron microscopy was carried out on the region outlined in Figure 3A using a previously described method (Maco et al., 2014; Maco et al., 2013). An image of the blood vessel pattern below the cranial window was taken immediately after the last 2PLSM imaging session. The anesthetized animal was then chemically fixed by perfusing, via the heart, 10 ml of isotonic PBS immediately followed by 200 ml of 2.5% glutaraldehyde and 2% paraformaldehyde in phosphate buffer (0.1 M, pH 7.4). Two hours later, the brain was removed and 60-μm-vibratome sections cut from the imaged region, tangential to the cortical surface. The region containing the imaged axons was located in these sections by matching the pattern of blood vessels seen in the first 2–3 sections with the image of the brain’s surface taken prior to the perfusion. The 2-photon laser was then used to burn lines into the fixed section around the region of interest. The laser was tuned to λ = 850 nm with a power of ~300 mW at the back focal plane of the objective. A rectangle of ~25 μm x 25 μm was burnt around the ROI using the line scan mode (2 ms/line). This shape acts as a fiducial mark that can be seen after the tissue has been stained, and resin embedded, ready for election microscopy. This section was then washed in sodium cacodylate buffer (0.1 M, pH 7.4) 5 × 3 min, postfixed and stained in reduced osmium (1.5% potassium ferrocyanide with 1% osmium tetroxide in 0.1 M sodium cacodylate buffer) for 40 min, followed by another 40 min incubation with 1% osmium tetroxide in the same buffer, and finally for 40 min with 1% aqueous uranyl acetate. Next, the section was dehydrated in a series of increasing concentrations of alcohol, then infiltrated in 100% Durcupan resin, and flat embedded between two glass slides and placed in a 60°C oven for 24 hr.

A small piece of the section (2 mm x 2 mm), containing the laser marked region, was trimmed from the rest of the section and glued onto a slab of blank resin, keeping the region of interest closest to the surface. This block was then trimmed in the ultramicrotome, using glass knives, so that the laser marks were within 5 μm of its surface, and less than 2 μm from one edge. The block was then mounted on a 45° aluminum stub, using conductive carbon paint, so that this edge was uppermost. It was then sputter coated with a 30 nm thick layer of gold. The sample was then placed in the focused ion beam scanning electron microscope chamber (FIBSEM; Zeiss NVision 40, Zeiss SMT, Germany) and orientated so that the ion beam could mill parallel to the laser marks, and therefore parallel to the 2PLSM plane of focus. The sample was imaged with an electron beam of 1.7 kV and a probe current of 1.1 nA using the back-scattered electron detector (ESB). Images with a pixel size of 6 nm were collected with a dwell time of 12 μs per pixel. The ion beam removed 12 nm of resin after each image, using a current of 700 pA, with an energy of 30 kV.

The serial electron micrographs were aligned in FIJI (http://fiji.sc; Schindelin et al., 2012) and the structures of interest segmented in the TrakEM2 part of the software (Cardona et al., 2012). Models were then exported to the Blender software (Blender Foundation, Amsterdam; http://www.blender.org), and the NeuroMorph toolset (http://neuromorph.epfl.ch; Jorstad et al., 2015) was used to measure axon cross-section areas, bouton and mitochondrial volumes, and PSD surface areas. The beginning and end of each bouton were defined, using the centerline processing tool, as the points along the axon where its normal cross-section area changed by more than 30% within a distance of 0.2 μm. The boutons were then delineated and their volumes calculated in the measurement tool (fourth row in Figure 4A–D and Figure 4E).

Figure 2—figure supplement 1A shows three maximum intensity projection views of an axon segment, which was traced independently by five users. Inter-user trace variability, which is usually more pronounced along the z-dimension (optical axis), can hinder bouton detection and measurement. Therefore, it is essential to optimize the layout of manual traces, ensuring that they accurately follow axon centerlines in the underlying image.

In this study, we used the optimization algorithm described in (Chothani et al., 2011). In this version of the algorithm, positions of trace nodes ${\displaystyle r_k={\left(x_k,y_k,z_k\right)}^T}$ are updated to optimize a fitness function, which includes intensity integrated along the trace and a regularizing constraint on the positions of neighboring trace nodes:

Here, vectors ${\displaystyle l_m}$ denote the positions of voxel centers in the image stack, index ${\displaystyle k^{\prime }}$ enumerates the neighbors of node ${\displaystyle k}$, parameter λ denotes the average density of nodes in the trace (number of nodes per voxel), and Lagrange multiplier α > 0 controls the stiffness of the trace. The first term in this expression represents convolution of the image with the Laplacian of Gaussian filter of size R, and due to the fast decay of the Gaussian factor, summation in this term can be restricted to a small number of voxels in the vicinity of the trace. The following parameter values were used to produce the results of this study. R was set to three voxels (roughly equal to axon diameter in the images), λ was set to 0.5, and α equaled 0.001. The fitness function was maximized with Newton’s method as previously described (Gala et al., 2014). Trace node positions were synchronously updated at every iteration step, and optimization terminated when the relative change in ${\displaystyle \mathcal{ℱ}}$ fell below 10⁻⁶.

Results of this optimization procedure (Figure 2—figure supplement 1B) show marked improvement over manual traces. The optimized traces are smoother, closer to one another in z, and follow the underlying axon intensity in the image more accurately. Following optimization, all traces were subdivided to a higher density of nodes, λ = 4.

Parameters of filters used to generate the axon intensity profiles, Equation (1), were chosen in the following way. The xy size of the multi-scale LoG_xy filter was chosen to span bouton sizes observed in 2PLSM images, R_xy ∈ [1.5, 3.0] voxels (0.39 μm to 0.78 μm). Since the observed bouton size in the z dimension is dominated by the point spread function of the microscope, a single size was chosen for the z component of the filter, R_z = 2 voxels (1.6 μm). Upon convolution with the image, this multi-scale filter returns the maximum intensity calculated over the specified range of sizes. The Gaussian filter, G, was chosen to have a fixed size of R = 2 voxels in all three dimensions. This size was selected to be roughly equal to the typical axon shaft radius observed in the images (Figure 2—figure supplement 1A).

Figure 2—figure supplement 1C shows LoG_xy profiles for the five manual traces from Figure 2—figure supplement 1A. Small differences in the layout of manual traces can lead to significant variability in intensity profiles, which is undesirable. Trace optimization reduces this variability to an extent undetectable by visual inspection, Figure 2—figure supplement 1D.

Peak detection algorithm, Equation (2), was initialized with a number of foreground peaks, ${\displaystyle N_f=\lceil L/0.5\mu \mathrm{m}\rceil }$, which is much larger than the expected number of putative boutons. In this expression, ${\displaystyle \lceil \rceil }$ denotes the ceiling function and ${\displaystyle L}$ is the trace length in micrometers. The number of background peaks, ${\displaystyle N_b=\lceil L/25\mu \mathrm{m}\rceil }$, was determined based on the observed spatial scale of background variability. The peaks were initially distributed uniformly along the entire length of the trace. Both, Gaussian and Lorentzian peak functions were tested, but only the former was used in this study as it provided a better fit to intensity profiles as judged by the value of the objective function.

The objective function was minimized with the gradient descent method. Gradient steps were taken simultaneously along the ${\displaystyle {\mu }_j^f}$, ${\displaystyle {\mu }_k^b}$, ${\displaystyle a_j^f}$, ${\displaystyle a_k^b}$, ${\displaystyle {\sigma }_j^f}$, and ${\displaystyle {\sigma }_k^b}$ dimensions. At every gradient step, parameters that moved outside the bounds specified in Equation (2) were set to these bounds. Upon convergence, that is, when the relative change in the value of the objective function became less than 10⁻⁶, one small foreground peak was eliminated or a pair of closely positioned foreground peaks was merged, and the resulting set of peaks was re-optimized. This procedure was continued until there were no peaks left that passed the following heuristic criteria: (1) a single foreground peak is marked for elimination if the peak amplitude is less than 0.3, and (2) a pair of overlapping foreground peaks is marked for merger if distance between the peaks is less than 1.0 μm or overlap area of the peaks exceeds 50% of either area. Here, the threshold of 0.3 was set to exclude small profile peaks that result from fluctuations of intensity along the axon. This threshold was chosen to be significantly lower than the mean profile intensity, which is one according to Equation (6). Distance threshold of 1.0 μm (size of a small bouton, see e.g. Figure 3) and threshold on inter-peak overlap were introduced to merge peaks corresponding to the same varicosity.

Custom software was used to match putative boutons over time. Because most large boutons remains stable over weeks to months, they can serve as fiducial points to match the remaining boutons. Here, piecewise linear registration of normalized intensity profiles was performed based on few fiducial boutons marked by users across imaging sessions (Figure 6—figure supplement 1). This was followed by matching the remaining putative boutons with a greedy, distance-based algorithm. All results were then validated with visual inspection. For putative boutons detected in some, but not all imaging sessions the missing bouton weights were filled in with the normalized intensity profile values at the corresponding registered positions. This procedure was performed for axons traced by different users (Figure 2—figure supplement 2), axons repeatedly imaged under different conditions (Figure 5), and axons monitored in the long-term imaging experiment (Figure 6).

Small inaccuracies in layout of axon traces can introduce errors in bouton detection. Therefore, it is essential to verify that trace optimization can guarantee sufficiently high precision in bouton detection. For this test, putative boutons were detected automatically, both with and without trace optimization, on the same set of 16 axon segment traced manually by five users. In the absence of trace optimization, a total of 578 candidate bouton sites were identified and matched across all five user traces (Figure 2—figure supplement 2A). A candidate bouton site here is defined as a position on an axon where a putative bouton was detected based on at least one user trace. For the majority of candidate bouton sites, 348 (60%), a putative bouton was present in all five user traces, and thus, there was a full consensus in bouton detection. However, 141 (24%) candidate bouton sites had one conflict and 89 (16%) had two. A conflict is defined as a dissent from the majority, and given five traces, the maximum number of conflicts is two. In contrast, after trace optimization, 474 candidate bouton sites were detected, of which 409 (86%) had full consensus, 39 (8%) had one conflict, and 26 (6%) had two (Figure 2—figure supplement 2B). With trace optimization, it was possible to obtain complete agreement in all putative boutons of weights greater than 2.5. This is a marked improvement over bouton detection with no trace optimization where 17 putative boutons of weights greater than 2.5 had one or more conflicts.

Tracing inaccuracies also affect the weights of detected putative boutons. Figure 2—figure supplement 2C, shows that inaccurate manual traces can contribute to bias (slope difference from 1) and variance (deviation from the best fit line) in bouton weight, and that trace optimization significantly reduces such errors. To quantify the uncertainty in bouton weight we calculated root mean square (RMS) weight difference for putative boutons detected with no conflicts (Figure 2—figure supplement 2D). As was expected, trace optimization dramatically reduces this uncertainty in the entire range of bouton weights.

We would like to point out that some amount of uncertainty in bouton detection (Figure 2—figure supplement 2B) and measurement (Figure 2—figure supplement 2D) remains even after trace optimization. However, as shown in Figure 5A–C this uncertainty of the method is much smaller than variability originating from the experimental sources, and therefore, it is not a limiting factor in structural plasticity measurements.

Axons of fluorescently labeled neurons were traced by using the manual tracing module of NCTracer software (http://www.neurogeometry.org). Custom software written in MATLAB, BoutonAnalyzer (Gala et al., 2017; copy archived at https://github.com/elifesciences-publications/BoutonAnalyzer), was used to optimize the traces, generate intensity profiles, detect putative boutons, and match these boutons across multiple user traces and imaging sessions. BoutonAnalyzer enables the user to visually inspect the results of bouton detection and matching and edit them if necessary. The NeuroMorph Centerline Processing and NeuroMorph Measurement tools were used to measure axon cross-section areas, bouton and mitochondrial volumes, and PSD surface areas (Jorstad and Knott, 2017). A copy is archived at https://github.com/elifesciences-publications/NeuroMorph.

Data used in the CLEM experiment are available in the Dryad Digital Repository (https://dx.doi.org/10.5061/dryad.3q50t). This dataset includes a 2PLSM image stack in TIFF format (Figure 3A), optimized traces of four axon segments in SWC format (Figure 3A), and 3D EM reconstructions of these axons in Blender format (Figure 3C).

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

We are grateful to Jerome Blanc, Raeann Dalton, Rammy Dang, and Chi Zhang for their help in tracing neurons and annotating boutons. This work was supported by the NIH grant R01 NS091421 and AFOSR grant FA9550-15-1-0398 to AS, Swiss National Science Foundation (SNF) grants 331003A_153448 and 51NF40_158776 (National Centre for Competence in Research, SYNAPSY), and the International Foundation for Research on Paraplegia (chair Alain Rossier) to AH, and SNF grant CRSII3-154453 to AH and GK.

Animal experimentation: All experiments were performed according to the guidelines of the Swiss Federal Act on Animal Protection and Swiss Animal Protection Ordinance. The ethics committee of the University of Geneva and the Cantonal Veterinary Office of Geneva, Switzerland (approval code GE/61/17) approved all experiments.

© 2017, Gala et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.