1. Introduction

1.1. Motivation

In the mid 2000s, optically-responsive proteins showed their use in in-vivo neuronal research. These proteins had the ability to indicate cell activity, as is the case for genetically encoded calcium indicators (GECI) [2], as well as to permit controlled neuronal excitation, as is the case for channel rhodopsins (ChR) [3]. The development of these molecular indicators and actuators opened an avenue of neuronal research which can be mediated entirely by light. Although this meant that the invasiveness of previous procedures could be greatly reduced, the ability to look at specific neural pathways remained limited to shallow regions of the brain no more than approximately 200 µm below the surface [4]. Highly scattering cortical tissue prevented deep probing by occluding and distorting light [5]. For this reason optically mediated research in deep nuclei, such as the basal ganglia, has been limited [1].

Since two photon microscopy offers cellular resolution at a maximum depth of approximately 1 mm [4,6], implants have been suggested to enable deeper probing [7–12]. To ensure cell vitality, these implants are required to have widths within a few hundred micrometers yet they must also maintain a high information bandwidth so that high temporal and spatial resolution images can be transmitted. For this reason, multimode fiber endoscopes have been developed for minimally invasive imaging [9–13]. Images, however, are not directly transmitted through the multimode fiber as is, but rather each image pixel is expressed as a superposition of fiber modes which, due to their differing propagation constants, interfere in a manner which gives a distorted image at the proximal end of the fiber. To mitigate this, phase shaping techniques have been used to account for modal dispersion and allow for high resolution imaging [10–13]. These techniques can offer diffraction limited spatial resolution yet they require scanning, thus making them considerably slow. When multimode fiber imaging is performed with acousto-optic deflectors (AOD) or digital micro-mirror devices (DMD) the temporal resolution is greatly increased to rates which can be suitable for functional imaging provided the spatial sampling rate is kept relatively low [13,14]. However, there remains another inherent problem in that multimode fiber endoscopes are highly sensitive to the fiber shape [15,16]. This results in a need for recalibration of the fiber endoscope every time its conformation changes.

Fiber bundles transmit an image through several cores with each core acting as a single pixel. As such, when imaging with incoherent light, they are impervious to bending. Furthermore, since the bundle acts simply as an image conduit, wide field microscopy can be performed making image acquisition fast. To minimize cross talk, however, the cores must be sufficiently far apart [17]. This intercore spacing dictates not only the highest spatial frequency detectable but also the number of sampling points for a given implant size which could be 10 to 50 times lower than for a multimode fiber of comparable size [18,19]. Low spatial resolution is the most common criticism of fiber bundles in the field of microendoscopy. As such, in this article, we present a method to overcome this low resolution problem by up-sampling fiber bundle images with compressive sensing (CS).

1.2. Compressive sensing

Compressive sensing is a numerical technique for reconstructing signals when they are sampled significantly below their Nyquist limit. It relies on the assumed sparsity of the acquired signal in a predefined representation basis [20]. The process can be thought of geometrically as is shown in Fig. 1, where a 2-pixel image is reconstructed from a single measurement. Let the point, labeled x, in Fig. 1(a) represent the true object in ℝ², in this case a 2-pixel image. In Fig. 1(b) the point, y, is measured in the basis vector shown by ψ₁. This gives us the information that the solution lies along the perpendicular orange line. This line is simply the null space of our incomplete measurement basis, Ψ, in this case a single vector in ℝ², offset from the origin by the measurement, y. In Fig. 1(c) the space undergoes a linear transformation, which is known to represent the original image in a sparse or quasi-sparse manner. This is shown by the Φ basis. The solution to the reconstruction problem, signified by the black point, x^, in Fig. 1(d). This is the intersection of the null space (orange line) and the smallest l₁-sphere in the representation basis (shaded diamond) which has an intersection with the null space. The l₁-norm is used as it seeks out sparse solutions while remaining computationally simple. Ideally, the l₀-quasi-norm should be used as it gives a count of non-zero elements. Using the l₀, however, gives a computationally difficult problem, therefore this constraint is typically relaxed [21]. For the general case where, instead of a two dimensional object reconstructed from a single measurement, an N dimensional object is reconstructed from M measurements (where M ≪ N), the reconstruction accuracy is dependent on the decay rate of the sorted coefficients in the representation basis. The problem above could be stated as

 

Fig. 1 Geometric depiction of compressive sensing reconstruction. a) Target in ℝ². b) Incomplete measurement of target in arbitrary basis. c) Space transformation to a sparsifying basis. d) Reconstruction by l₁ minimization

Download Full Size | PDF

2. Simulations

To apply CS to wide field fiber bundle imaging we must first know the measurement basis of our system. Since we are aiming for high resolution reconstructions of fiber bundle images, the basis is dictated entirely by the bundle structure. To ascertain this, we imaged the end face of a Fujikura FIGH-016-160S fiber bundle under uniform white light illumination. We then used a segmentation algorithm to determine the centers and diameters of each of the cores in the bundle. Making the assumption that a lensless bundle has a zero working distance, and that each core couples light in proportion to its area, we formed the theoretical M × N measurement matrix, Ψ. The number of measurements, M, was 1461 given by the number of cores and N was the number of pixels in the image to be reconstructed, which was varied to assess the performance of the reconstructions under different compression ratios. The whole-image discrete cosine transform (DCT) and total variation (TV) were used as representation bases, Φ, as their coefficients decay rapidly when imaging cells, thereby lending themselves to sparse reconstruction. Reconstructions were performed using the basis pursuit algorithm developed by Zhang et al. which solves Eq. (1) iteratively [22]. Results are summarized in Fig. 2.

 

Fig. 2 Simulated fiber bundle image reconstruction. a) Mean relative error at varying compression (N/M) ratios over 29 images of cortical and striatal cells. Shaded regions give the 95% confidence intervals. b) Sample red fluorescent protein (RFP) stained cortical neuron image. c) Image simulated through 160 µm diameter fiber. d)–g) Image reconstruction examples at a compression ratio of 95 by d) Fourier domain filtering, e) Interpolation, f) Compressive sensing with the DCT as a representation basis, and g) Compressive sensing with TV minimisation.

Download Full Size | PDF

Image reconstruction through CS was compared to two common techniques for removing pixelation during fiber bundle imaging, filtering, and linear interpolation. Whereas interpolation was found to perform slightly better than spatial filtering, both were outperformed by CS at all compression ratios. As the compression ratio increased, higher spatial frequencies, which cannot be recovered either by filtering or interpolation, became more vital in reconstruction, leading to a greater difference between CS and these techniques. We note, however, that anomalies in the dead space between the cores cannot be recovered with any of these methods as the zero-working-distance assumption dictates that each core accesses only the region of the image directly beneath itself, leaving the dead space between the cores unsampled. The merit of CS, however, lies in our prior knowledge of a sparse representation for our input images. To gauge the sparsity of the images in the DCT representation the errors were plotted alongside the best M-Term DCT approximation of the original image in Fig. 2(b). This error gives some indication of how sparsely the image is represented in Φ, as well as a minimum possible error at the given sampling rate. This was not done for TV-CS as the TV is not expressed as an invertible transform and does not lend itself to this analysis.

3. Image reconstruction

3.1. Theoretical measurement basis

For reconstruction of images through the fiber we first determined the theoretical measurement basis by computationally segmenting the end face of a Fujikura FIGH-10-350S fiber bundle. For this the same assumptions as those stated in the section above were used. To avoid being in direct contact with our sample we used a 500 µm diameter graded index rod lens at the end of the fiber. To determine resolvable features of each of the image reconstruction techniques we imaged a high resolution USAF 1951 resolution target under wide field transmission illumination. The results of the reconstructions are summarized in Fig. 3.

 

Fig. 3 Wide field image reconstruction with a theoretical measurement basis at a compression ratio of 10.46. a) Line traces shown in red in c), d), e), and f). b) Original image as seen through the fiber. Image reconstruction by c) filtered back projection, d) interpolation, e) DCT compressed sensing, and f) TV compressed sensing.

Download Full Size | PDF

In Fig. 3(e) a line trace with three peaks, corresponding to the three bars on the target, is clearly visible when using DCT-CS reconstruction. This is not the case for either interpolation or Fourier domain filtering. For TV-CS we see three levels however there are not three distinct peaks. This illustrates that the reconstruction quality relies on a selection of appropriate basis for sparse representation. Here the wave-like nature of the DCT works well to sparsely represent the similarly wave-like lines of the target. In the original image in Fig. 3, distortions are caused by core-to-core cross talk, fiber imperfections, and scattering particles, such as dust, on the fiber end face. These are generally minor imperfections, however; their effects are carried through the reconstructions. Hence, we cannot justify the assumption that the cores are completely homogeneous in the way in which they couple light. Furthermore, the introduction of a rod lens invalidates the assumption of zero working distance, meaning that a single point source at the object can, in fact, couple to multiple cores. We address this in the following subsection.

3.2. Observed measurement basis

To address the issues mentioned in the earlier section we determined the true fiber measurement basis. To do this, the segmentation algorithm was used to separate an image of the fiber face into 9,200 intensity values, one for each core in the bundle. Next, the distal end of the fiber was raster scanned across a point source in a 300×300 square grid, which covered the 115 µm field of view. After background subtraction and thresholding, the values were recorded in a sparse measurement matrix, Ψ, which mapped each of the 90,000 object points at the distal end of the fiber to 9,200 core intensities at the proximal end. To perform this calibration a 1 µm pinhole was used as a point source. This was placed on the stage S2, shown in Fig. 4(a) while the mirror was moved from position M(1) to M(2) to provide illumination. Figure 4(b) shows the image of the pinhole through the fiber when in focus. Once the calibration was completed the mirror was moved back to M(1) for fluorescence imaging. The fiber was used both for fluorescence excitation and collection. An example image of cells taken through the fiber with this method is shown in Fig. 4(c).

 

Fig. 4 Imaging set up and example fiber images. a) Set up arranged for fluorescent imaging. M(1) and M(2) are calibration and imaging mirror positions, respectively. ExF and EmF are excitation and emission filters. DC is a dichroic mirror. S1 is a 3-axis motorized stage whereas S2 and S3 are mechanical stages. O1 and O2 are objectives and TL is a tube lens. The flip mirror, FM, is used to verify the fiber images by imaging through O2 directly. b) An image taken during calibration of a point source through the fiber. c) Example of ex-vivo fluorescent cell imaging through the fiber with mouse motor thalamus neurons labelled by hSyn-mCherry-AAV₁ viral injection.

Download Full Size | PDF

After the measurement matrix Ψ was determined, it was preprocessed to account for the motor backlash error which occurs during scanning. Next the columns were normalized so that Ψ formed a unit-basis dictionary with which we can solve the problem stated in Eq. (2). This was solved by the L1/L2 model algorithm developed by Zhang et al. [22]. The normalization factors of the columns of Ψ were taken into account when performing the reconstruction. Using incoherent 532 nm excitation we acquired images of 1 µm fluorescent polystyrene beads. When reconstructing these images three compressed sensing methods were attempted and compared with results from filtering and interpolation. The difference in the three methods was in the selection of the representation basis. The first method used the DCT as described above, whereas the second method used the point (canonical) basis, and the third the total variation (TV). Both the point basis and total variation basis had been shown to be successful for compressed imaging of beads [23, 24]. This is because a sparse sample of beads often leaves much of the image dark making these bases suitably sparse representations. For TV reconstruction an algorithm developed by Li et al. was used [25]. The results are summarized in Fig. 5.

 

Fig. 5 Imaging 1 µm fluorescent beads. a) Image as seen through the fiber bundle. b) Image of beads as seen through the epifluorescence microscope (O2 in Fig. 4(a)). Reconstructions through c) filtering, d) interpolation, e) DCT compressed sensing, f) point compressed sensing, and g) TV compressed sensing. h) Line traces for each reconstruction method. The trace is taken across two beads shown by the location of the red line in b).

Download Full Size | PDF

Although there is a significant overlap, the line trace in 5(b) shows two distinct peaks for the two beads in the reference epifluorescence image. For the measurements on the distinction of the two beads, all the CS based methods showed two distinct features, whereas other methods did not. This reconstruction would not have been possible with the measurement matrix based on the assumptions made in the previous section. Since the spacing between the beads is very small, both beads would be detected within a single core or one bead would be in the intercore dead-space. The introduction of a rod lens and high resolution raster scanning ensures there are no blind spots, of size 1µm or greater, while also ensuring any two features 1µm or greater could be distinguished. In order to do this, knowledge of how light couples to and transmits through the fiber is necessary and this is exactly what is coded into the measurement matrix Ψ. The reconstructions with interpolation and filtering do not make use of this knowledge nor any predefined notions of sparsity and hence do not perform as well.

For reconstructing beads, which can be considered as point sources, the point basis is, theoretically, a very suitable representation basis. As can be seen in Fig. 5, the point basis works well to identify and localize the beads. This method, however, is very prone to type 1 errors (false positives) in the presence of noise, often resulting in an overestimation of the number of beads present in a sample. The relative unsuitability of this basis stems from the fiber structure itself. The Ψ matrix is dependent on the bundle arrangement and coupling, and, as such, it is completely constrained. In contrast, for the Φ matrix we are somewhat free to choose what represents the sample sparsely. There is, however, another constraint on Φ imposed by the already constrained Ψ. It has been shown that in order to ensure high probability of successful reconstruction for an S-sparse image, that is an image with S non-zero coefficients, the number of measurements, M, must scale with µ²Slog(N) [26]. Here, µ is a measure of coherence between bases and is defined as

Since Ψ is relatively similar to the point basis, µ has a relatively large value when using the point basis for Φ; therefore, with M and N held constant, the point basis representation would only work well for reconstructing images in very sparse cases. For cells, however, the situation is different. The DCT represents them sparsely while remaining largely incoherent with the fiber-constrained Ψ, thereby driving the value of µ down. We use this to our advantage when reconstructing images of cells. Finally, this technique was used to reconstruct images taken, through the fiber, of fluorescently labeled cholinergic interneurons. To prepare the slices for ex-vivo imaging, an anesthetized wild-type C57BL/6J male mouse was perfused intracardially using warm phosphate buffer 0.01 M (pH 7.4). This was followed by phosphate buffer perfusion containing 4% paraformaldehyde and 14% picric acid. The brain was removed after decapitation, and postfixed in the same solution for 24hrs at 8°C. The brain was then cryoprotected in a 50/50 mixture of fixative and 20% sucrose in 0.01 M phosphate buffered saline (PBS) for at least 24h or until the brain sunk in the solution. Next, 60 µm coronal sections were cut on a freezing microtome with a freezing stage (Yamato electrofreeze, MC-802A) and washed in chilled PBS (4×5min). Afterwards the sections were processed for immunocytochemistry to visualize cholinergic interneurons (ChIs) in the striatum. Sections were incubated for 4h at room temperature in 10% normal donkey serum (NDS) diluted in PBS. To identify cholinergic interneurons, the slices were incubated in antibody against choline-acetyltransferase (Millipore, California USA, 1:100 dilution) in antibody diluent (Triton0.3%, NaN3 Azide0.05%, PBS) in a dark room for at least 24h at 4°C on a shaker and washed in PBS (4 × 5 min). The brain sections were further incubated in a secondary antibody in a dark room for 12h at 4°C on a shaker (dilution 1:200). The secondary antibody used was Alexa fluor 488-donkey anti-goat (life technologies, Eugene, OR, USA; 1:200). The sections were washed again in PBS (4 × 5 min) and the samples were finally imaged in the fiber bundle setup and reconstructed. An example is shown in Fig. 6.

 

Fig. 6 Reconstruction of fluorescently labeled striatal interneurons. a) Image of an interneuron as seen through the fiber. b) Reconstructed image of the interneuron with DCT-CS.

Download Full Size | PDF

The reconstruction of the cell in Fig. 6(b) faintly brings into view some dendrites which are highlighted by the white arrows. The main issue hindering small features from coming into view is the low signal to background ratio which is further diminished due to the fiber autofluorescence [27]. Furthermore, there are reconstruction artifacts in the form of vertical and horizontal lines. These are characteristic of using the DCT as a sparsifying basis, which favors features aligned with the Cartesian grid. With the appropriate selection of representation basis, however, this problem could be mitigated. More complex bases have been suggested for reconstructing cells. For example, the curvelet transform, which can express images sparsely that are mostly uniform but have rare, yet abrupt, changes of intensity along edges, as is the case with cells and dendrites [28]. The artifacts, in such cases, would conform to the contours of the imaged features rather than cross over them as is the case for the DCT. The major draw backs of using such a representation basis is that it requires much more computational power and, therefore, reconstruction time. Furthermore it is expected that such a representation basis would have a higher coherence with the fiber-bundle-basis further increasing the convergence time and reducing the probability of accurate reconstruction. This remains an avenue to be explored.

4. Conclusion

We have exploited the strengths of fiber bundle micorendoscopy by making use of its flexibility with incoherent light, its ability to achieve high temporal resolution with wide field imaging, and its relatively small biological footprint. To overcome the major criticism of fiber bundles, which is that their spatial sampling is low, we have shown that it is possible to up sample fiber bundle images with the appropriate choice of basis. We foresee this as a viable technique for in-vivo wide field fluorescence microendoscopy.