Deep brain fluorescence imaging with minimally invasive ultra-thin optical fibers

A major open challenge in neuroscience is the ability to measure and perturb neural activity in vivo from well-defined neural sub-populations at cellular resolution anywhere in the brain. However, limitations posed by scattering and absorption prohibit non-invasive (surface) multiphoton approaches for deep (>2mm) structures, while Gradient Refreactive Index (GRIN) endoscopes are thick and cause significant damage upon insertion. Here, we demonstrate a novel microendoscope to image neural activity at arbitrary depths via an ultrathin multimode optical fiber (MMF) probe that is 5-10X thinner than commercially available microendoscopes. We demonstrate micron-scale resolution, multispectral and volumetric imaging. In contrast to previous approaches, we show that this method has an improved acquisition speed that is sufficient to capture rapid neuronal dynamics in-vivo in rodents expressing a genetically encoded calcium indicator. Our results emphasize the potential of this technology in neuroscience applications and open up possibilities for cellular resolution imaging in previously unreachable brain regions.


Introduction
The main potential advantage of MMF for biological endoscopy is their thin diameter (50-150 um) which induces less damage compared to thicker probes, such as GRIN lens 4,[9][10][11] . Fiber bending remains a major challenge to MMF imaging techniques that are based on Wavefront shaping (WFS) since fiber deformation changes how modes are coupled and invalidates precomputed transformations critical for image formation [12][13][14] . However, in some domains, such as neuroscience, major fiber bending may not always pose a direct problem since in many experimental designs fibers can be inserted into the brain along a straight trajectory to target specific brain regions. Furthermore, minor perturbations at the distal end of the fiber have only minor effects on the imaging capability if a proper fiber is selected 15 . Other challenges, on the other hand, need to be addressed to make this technology useful for neuroscience. First, acquisition speed needs to be sufficiently high to capture rapid neural firing events (ms scale), while maintaining sufficient spatial resolution to resolve cellular level details. Previous approaches 1,5-8 that use liquid crystal spatial light modulators (LC-SLM) have not been shown to be capable of this temporal resolution which would prohibit rapid sampling of neural signals. Second, the system needs to have high collection efficiency to capture small fluorescence changes evoked by calcium transients without causing photobleaching. To date, only fixed tissue has been successfully imaged 5,15 . Third, the system should image some distance away from the fiber tip since tissue in close proximity may not be functioning properly due to possible damage inflicted by fiber insertion.
Here, we focus on addressing these challenges and propose a novel microendoscope system that is based on a digital mirror device (DMD).

Generating phase modulation with a digital mirror device
The basic idea behind the Lee hologram 16 technique is to create a binary mask that creates a diffraction pattern with the desired phase modulation. Following the work of 16,17 , to create a 2D spatial phase mask (x,y), we first define a spatial carrier wave with frequency f and rotation ,over the mirror domain [x,y]. The mask is then defined as: (1) = cos( ) + sin( ) This mask generates three peaks on the Fourier plane since By blocking the first two terms (DC and one diffraction order) with an iris and allowing only one diffraction order to pass through, we create a 2D phase array that is up to a constant shift from the desired Φ( , ).In the limit, each phase Φ( , ) cannot be represented by a single mirror on the DMD, hence multiple mirrors are grouped to represent a single desired phase. We have experimented with varying sampling (number of mirrors per phase) and found that 10 is a good compromise between the size of each block and sufficient remaining mirrors to deliver the constant reference phase (see below).

Transmission matrix estimation and spot generation
The main derivation of the Transmission matrix method can be found in 1 . The goal of this procedure is to estimate the complex matrix K, which describes how the complex input field is mapped to the complex output field . The assumption is that this mapping can be approximated by a linear relationship: (where * denotes matrix multiplication). In our case, the input field can be described as a 2D array of wavelets (64x64), where each wavelet can be phase shifted independently using the Lee hologram presented on the DMD. In other words, The dimensionality of the output is determined by the resolution of the camera that images the tip of the fiber. In our case, dim ( ) = 128 128.
Since we determine the input, K can be recovered by However, cannot be observed directly because it is complex (i.e., only the intensity| | 2 can be measured. To recover the complex value of Eoutwe use an interferometric approach. We present both and a second reference beam . In our case, the second beam is actually the region surrounding , which we give a fixed constant phase value (i.e., 0). By phase shifting the entire input field by a certain values one can recover (up to an unknown fixed phase that is determined by the reference beam which surrounds our input field). Given a specific input field that has been phase shifted by and the reference beam r, we can write an expression for the observed intensity at pixel (i.e., the image is flattened), which describes the interference between the two complex waves (the reference and the input): Notice that the first term among the last three does not depend on , and the second is phase independent. If we define ≜ ��� , it can be shown that Where is a complex number and corresponds to a complex vector. Thus, a single presentation of an input field (4096x1) gives us (16384x1) constraints on K. To fully recover K (16384x4096), we need different inputs. The matrix derivation can be formulated as follows: (8) = ( ) where Z is a complex matrix calculated from the measured intensities, K is the unknown transmission matrix and E is a set of complex field inputs ( = [ 1 , 2 , … , 4096 ] of our choice. If we let = (the identify matrix), the solution is trivial (~). However, this is not efficient in terms of SNR because only a tiny fraction of the mirrors would be turned on for any given input. Any choice of E is reasonable as long as it is orthogonal. In practice we use the Hadamard-Walsh basis because it is balanced in terms of the mirrors that are set ON and OFF and also because it is easily invertible −1 = . We first map the Hadamard basis from the [-1,1] domain to the complex phase � * 0 , * , � (i.e., a complex field that we can generate with our Lee-hologram). We denote the complex phase basis as H. K is easily recovered by ~�. Under the assumption that K is unitary, * = , the inverse operation becomes: , where * denotes the complex conjugate operator. To generate a spot at location (x,y), K is inverted and multiplied by a delta function at the desired pixel location , where index is the flattened position of (x,y). Any pattern (with certain limitations) can be generated by multiplying it with −1 . For example, it is possible to generate larger excitation spots which are useful for rapid sub-sampling scanning or structured illumination to target specific structures. Only the phase angle of the above derivation is needed since only phase modulation is applied. In practice, the above calculations can be simplified and reduced to the following short Matlab code: Z=reshape(atan2(I 0 -I π/2 ,I π/2 -I π ), 128 2 , 64 2 )'; =reshape(atan2(Hsin(Z), Hcos(Z)), 64,64,128 2 ); Where I α represents a 3D volume (128x128x4096) of images obtained when projecting H with phase offset α and represents a 3D volume where each 2D slice contains the phases needed to form a single point. These calculations are sped up considerably with a custom C++ implementation that uses Cublas (Cuda optimized matrix library for GPU). A total of 12288 patterns need to be displayed (4096x3) on the DMD to measure K. With a reasonably fast CMOS camera (850 Hz) this operation takes ~14.5 seconds, while ~5 additional seconds are needed to perform the remaining matrix multiplications.

Multispectral imaging
To allow multispectral imaging we first recognized that the DMD will diffract the two excitation beams to different directions. These diffraction patterns can be steered by controlling the carrier wave frequency and rotation such that only one excitation wavelength passes through the iris with minimal interference from the other laser. Two calibrations (one for each laser source) are needed. Rapid transitioning between the two sources is then achieved by setting the DMD with the desired pattern.

Enhancement Factor Metric
To evaluate the quality of our generated spots we devised an automatic procedure that calculates the TM, generates several hundred spots and image each one of them under several neutral density filters. The multiple exposures of the same spot are then combined to form a high dynamic range image in the following way. First, for a given exposure under neutral density n, a dark image is subtracted and all pixels that are over-exposed or underexposed are removed . Then, the high dynamic range (HDR) image is calculated as where n corresponds to the neutral density value. The peak of the HDR image is found ( ), and a small neighborhood (5x5 pixels) around that region is removed. Then, the average signal inside the fiber core is calculated ( ). The enhancement factor is defined as / . We believe this approach is less susceptible to numerical instabilities from underexposure, which can artificially inflate enhancement values. This is typically the case since the peak is several orders of magnitude brighter than the background, which exceeds most camera's dynamic range.

Point Spread Function Estimation
To measure the (x,y) spread of the point spread function we imaged fluorescence microspheres and fitted each individually identified blob with a 2D Gaussian by fitting (a dark subtracted image) with the following function: The FWHM of the Guassian was defined as 2 √2 2 , The signal to noise ratio (SNR) was defined as A/N, where A represents the amplitude of the 2D Gaussian and N is the average intensity measured in a background region not containing microspheres.

Microspheres to background separability
We used the sensitivity index (d') from signal detection theory to quantify the Z-section with the sharpest focus. d' is defined as ′ ≜ − � 1 2 , where is the signal mean intensity, is the noise mean intensity, is the signal standard deviation and is the noise standard deviation.

Software, Data collection and processing
The software to operate the fiberscope was mostly written in Matlab (mathworks). Some modules were implemented in C++ and Compute Unified Device Architecture (CUDA) to speed up calculations. In an earlier version of the microscope we digitized PMT emissions using a DAQ capable of sampling at 500kHz. For each DMD flip, 10 samples were collected. In a newer version we developed a custom FPGA code running on a Xilinx chip that was triggered by a DMD flip, waited for 15sto allow mirror stabilization on the DMD and then generated TTLs at 8MS/s to drive the DAQ to sample emissions from the PMTs (total dwell time: 50 usecs at 20 kHz). This resulted in 180 samples / pixel, which were averaged to get the raw image data. We resampled movies such that each pixel was sampled at time t (to correct the small deviation caused by in-frame raster scan time). Unless otherwise stated, images were notch filtered to remove unwanted frequency bands (typically, 60Hz), then spatially smoothed with a Gaussian (=1.5-2m).For functional imaging we also temporally smoothed the sequence with a gaussian (=1.5-2 frames). In some plots, histogram equalization was performed on images (dynamic range stretching) to improve visualization of faint structures and to map a higher dynamic range to 8-bit images. Wild type mice (C56BL) from Jackson Labs were anesthetized with isoflurane (1-2%). A small borehole (~300-500um) was drilled and 50-100nl of AAV9-hSyn-GCaMP6s was pressure injected using a custom pump (25nl/minute) to deliver the GCaMP gene. We allowed 2-3 weeks of recovery and virus expression prior to imaging. Mice were then anesthetized again and a slightly larger craniotomy was performed to expose the brain. We used a thin wall hypodermic needle (28G) as reinforcement on the upper part of the fiber, leaving ~2mm of fiber exposed for insertion (the needle was not pushed in the brain). Animals were then raised on a custom movable translational stage in small increments until the fiber punctured through the dura and entered the brain. At the end of the imaging session, animals were euthanized with a lethal dose of pentobarbital.

Results
The microscope design (figure 1) is based on wave front shaping (WFS) 7,18 : the light wavefront is modulated before it is coupled at the proximal side such that when it exits at the distal tip, it creates a micron-scale excitation spot at a desired (x,y) position (figure 2) and distance (z) from the fiber tip (figure 3). Sections (or the entire volume) are reconstructed by scanning the excitation spot and measuring the fluorescence emissions collected by the same fiber. To generate wave front phase modulation we use a computergenerated hologram method 14,16,17 on a digital mirror device (DMD), that allows kilohertz phase modulation (methods). Random access scanning is possible since forming a spot at (x,y,z) corresponds to programming the DMD with the precomputed pattern. Structured illumination is also possible, allowing to dynamically changing the size of the PSF (methods). The 2D array of adjustable phases mix with a fixed reference beam in the fiber (figure 1) to generate the observed complex interference pattern at the distal tip (speckle pattern). The input-output transformation of the fiber is characterized using the transmission matrix (TM) algorithm 1 (methods) and the complex field is measured using a fast (three step) phase stepping algorithm. The entire procedure (termed "calibration") of measuring the TM and computing the holograms needed to scan the entire FOV (100x100 um) takes ~20 seconds using a fast CMOS camera and GPU optimized software (methods).
To assess the optical quality of the microendoscope we first developed a robust metric to measure the enhancement factor (EF), quantifying the strength of the excitation spot at the distal tip of the fiber (methods). We image generated spots with multiple neutral density filters and combine them to form a high dynamic range image (figure 4a) from which the EF can be calculated (methods). Typically obtained EF values were 305125(range: 50-550) and decreased as a function of radial distance (figure 4b). To further assess the optical quality we measured the point spread function (PSF) along the emission path by imaging 0.96um fluorescent microspheres ( figure 4c). The FWHM of a 2D Gaussian fitted to the blobs was 2.100.25 um (mean, std), with estimated SNR of 1.390.1(methods).    The ability to image away from the tip in the brain is important since cells nearest to the fiber are more likely to be damaged. To assess our ability to form excitation spots away from the fiber tip in a scattering medium we imaged 15um microspheres embedded in 2% intralipid agarose (figure 5). By shifting the calibration camera by a fixed offset and remeasuring the TM, it is possible to find the phases needed to raster scan at different planes away from the tip 19 , essentially acquiring a volume without the need to physically move the fiber during imaging. We quantified the separability between microspheres embedded in this scattering medium and the background using d' (methods) and found the highest d' at 100um, suggesting the sharpest focus was obtained at that Z-section (figure 5). Our approach can also be used to rapidly switch between laser sources (methods) to allow multispectral imaging. To demonstrate this capability we imaged two types of microspheres with different excitation and emission spectra (figure 6).  Fluorescent microspheres can be 10-1000 brighter than invivo biological fluorescence. To test whether our microendoscope has the sensitivity needed to capture enough photons from cells expressing GFP we set up an experiment to simultaneously image live cells on a glass slide with both the fiber from the top and a wide field epifluorescence microscope from the bottom. We imaged live baby hamster kidney (BHK) cells expressing GFP (figure 7) and found that thin fiberblast features could be easily recognized.

Figure 7. In-vitro imaging of BHK cells expressing GFP.
Comparison between full field epi imaging with a regular camera and the fiberscope.
The main innovation of our approach is the ability to rapidly scan the entire field of view in 3D to allow real-time measurements of calcium events evoked by spikes in-vivo. We can currently achieve full frame (100x100um) acquisition speeds of ~7-10 Hz, while smaller fields of view (~ 10x10um) can be imaged at speeds exceeding 200 Hz. To assess whether the microendoscope has collection efficiency sufficient to capture small fluorescence changes at such rates, we imaged a thin hippocampal neuronal tissue culture expressing the genetically encoded calcium indicator GCaMP6f 20 . Multiple neurons could be identified in the FOV ( figure 8). The fluorescence time course measured with the fiber was in good agreement with epi microscope measurements (mean correlation : 0.54, n=6). These results suggest our method can be applied to record activity in-vivo. To test this directly we targeted primary visual cortex in wild-type mice expressing GCaMP, delivered via viral-mediated transduction. We slowly raised animals on a motorized translational stage towards the fiber in steps of ~50um. Upon fiber insertion, several neurons could be identified as they entered the field of view (figure 9a). Since neural activity under anesthesia is low, clearer images were obtained by calculating the standard deviation over time (figure 9b). Neurons exhibited varying degree of activity and were not significantly photobleached over several minutes exposure at ~60nW (figure 9c).

Conclusion
In summary, we present a new microendoscope imaging system that is based on WFS using a DMD to control light in a single fiber. To the best of our knowledge, this is the first demonstration of using multimode fibers to image live cells and in-vivo neurons and represents a significant improvement over single pixel imaging fiber photometry 21 .
Although the EF in our imaging system may not be as high as one that can be obtained with LC-SLM based approaches, we demonstrate it is sufficient to acquire images which allow identification of fine features at the neuronal scale. A key improvement over previously proposed methods [22][23][24][25] , which is required for live imaging, is the ability to rapidly scan a large field of view, while maintaining sufficient spatial resolution to resolve cellular details and a high enough sensitivity to precisely measure small fluorescence fluctuations. The ability to rapidly switch between two light sources and to form structured light patterns may be useful to manipulate specific cells within the field of view using optogenetics. We expect this technology will open new avenues for addressing circuit level questions in deep brain regions that have thus far been unreachable with existing imaging modalities.