2.1.

Interface Simulation

We calculate the fluorescence sensitivity profile for a single fiber inserted into the tissue using Monte Carlo simulations of 10,000,000 photon packets launched from the fiber as described in Refs. 15 and 16. The tissue is modeled as a $1\text{-}{\mathrm{mm}}^3$ volume broken into isotropic $5\text{-}\mu \mathrm{m}$ voxels, with the tip of the fiber positioned in the center of the volume. Photon packets are emitted from the fiber with an angle determined by a Gaussian distribution reflecting the numerical aperture of the fiber (0.377). The packet is modeled moving through voxels and at each step can be either scattered or fractionally absorbed. Absorption is dominated by hemoglobin in the blood, and as a result we calculate the absorption coefficient based on a 3% blood volume fraction¹⁷ with a 15-g/DL hemoglobin concentration¹⁸ and a 70% oxygenation fraction. Given the extinction coefficient of hemoglobin,¹⁹ our absorption coefficients are ${\mu }_a=0.337{\mathrm{mm}}^{-1}$ for 490-nm light, ${\mu }_a=0.343{\mathrm{mm}}^{-1}$ for 512-nm light. Scattering is described well in the literature, and we use a scattering coefficient of ${\mu }_s=20{\mathrm{mm}}^{-1}$ with anisotropy $g=0.9$.^20,21 By averaging the paths of all photon packets, we construct a 3-D photon distribution throughout the tissue.

The sensitivity profile for a bundle of fibers is obtained by taking the linear superposition of the profile for a single fiber. Based on our prior histology documenting the distribution of bundles of hundreds to thousands of fibers,¹⁵ we estimate the distribution of fibers in the tissue as a bivariate normal distribution in $xy$ space. To explore the interface properties under various conditions, we vary the fiber density by changing the number of fibers for a fixed standard deviation of the spatial distribution of the fibers. Unless otherwise noted, we use a standard deviation of $\sigma =150\mu \mathrm{m}$ for $xy$ splay, consistent with the observed distribution of fibers at a depth of 2 mm. The fiber depth distribution will vary based on how the fibers are cut prior to insertion; for our simulations, we model this variability as a normal distribution of depths with a standard deviation of $\sigma =15\mu \mathrm{m}$. The angle of each fiber is assumed to be uniform and parallel to the $z$ axis.

By superimposing the sensitivity profile on the tip of each fiber in the bundle based on the splay distribution,^22,23 we can evaluate the interface properties of the full bundle. The bundle’s ability to deliver excitation light to a given neuron is equal to the sum of all superimposed sensitivity profiles, whereas the bundle’s ability to collect fluorescent emissions from a given neuron is measured in terms of the few fibers that receive the majority of the emitted light from that neuron. The combination of the excitation and emission forms a mixing matrix

2.2.

Neural Population Simulation

Equipped with this model of the optical interface properties for a bundle of microfibers, we can now evaluate the ability to measure signals from a simulated neural population. We simulate a three-dimensional (3-D) volume of tissue, of sufficient size to ensure that it includes all regions of non-negligible fluorescence sensitivity (usually around $1.2{\mathrm{mm}}^3$), with a uniform distribution of neurons. We vary the density of neurons to understand the impact of signal density, but use densities that are consistent with relevant subpopulations of interest.^24,25 Unless explicitly stated, simulations use a density of 250,000 neurons per ${\mathrm{mm}}^3$. We assume all neurons express the relevant genetic probe.

Neural activity is simulated as independent spike events (Bernoulli processes).²⁶ The spike probability at each time step is calculated based on the spiking frequency, which unless otherwise stated is 0.4 Hz. Spikes are then convolved with a GCaMP6 waveform; unless otherwise stated, simulations use the GCaMP6f waveform with a rise time ($t_{\text{peak}}=0.14\mathrm{s}$) and an exponential decay ($t_{\frac{1}{2}}=0.32\mathrm{s}$).²⁷ Fluorescence was modeled at 100 time steps per second, and then downsampled based on the simulated frame rate for the recording setup.

The fluorescence traces for each neuron are then combined via the mixing matrix described in the previous section to produce fluorescence traces for each fiber. The fluorescence collected by a fiber can be written

2.3.

Source Separation

Given the linear mixing process inherent in the fluorescence signal collected by each fiber, the data are well suited to blind source separation techniques to estimate the underlying neural fluorescence signals. To improve the performance of the source separation, we developed an approach that incorporates our knowledge of the mixing process and associated fluorescent indicator dynamics. Similar to existing techniques that use deconvolution to estimate precise spike timing,²⁸ we initially inverse filter the recorded signal to remove the fluorescent indicator dynamics. This step requires an estimate of the exponential decay associated with the calcium response and indicator time constant and can be omitted for indicators where these dynamics are not known or are not stereotyped

Next, we apply the non-negative-independent component analysis algorithm²⁹ to perform blind source separation under the constraints of non-negativity. The algorithm first whitens the data, and then applies a series of orthonormal rotations to reduce the error between $Z$ and the reconstruction of $Z$ from the rectified (non-negative) components of $S$.

Having separated the independent signals, we can then apply the waveform filter ($\overrightarrow{w}$) to produce fluorescence traces that correspond with the identified independent components.

As there are many more neurons than fibers, separating all neural activity is an underdetermined problem. This technique is limited to extracting as many neural traces as there are fibers.

2.4.

Evaluating Source Separation

The evaluation of the source separation must consider both that the output of the source separation is a small subset of the underlying neural signals and the underdetermined nature of the blind source separation. To this end, we calculate the correlation coefficients of each separated signal (${\widehat{\overrightarrow{x}}}_i$) with each true, simulated fluorescence trace (${\overrightarrow{x}}_i$). The quality of the separation is then evaluated based on how many of the separated signals have a sufficiently high correlation with a true trace. In Sec. 3, we evaluate threshold choices and how they impact the number of extracted “true” signals.

To understand the source separation in the context of an experimental question, we further evaluate the separated signals on a spike detection task. Using the approach of detecting threshold crossings as a proxy for neural spikes, we can compare threshold crossing events in the separated signal with those in the matched underlying neural trace. This forms a binary classification task that can be measured in terms of a receiver operating characteristic plot and the accompanying area under the curve (AUC).³⁰ The AUC represents how well threshold crossings on the separated signal can be used to approximate threshold crossings on the true neural fluorescence; an AUC of 0.5 indicates poor performance (unable to separate), whereas a value of 1 indicates perfect performance.

As a control, we compare correlations between the separated signals with a set of randomly generated traces.

2.5.

Code

The code to generate the fiber sensitivity profile and to run all described models are publicly available at https://codeocean.com/2018/11/12/fiber-source-separation/code. The MATLAB modeling and source separation implementation is available on Github https://github.com/nathanntg/fiber-source-separation.

3.1.

Relevant Experimental Technique

The relevant experimental technique being modeled, including the optical interface, implant methods, and motivating histology, is fully described in Ref. 15. In brief, we use commercially available flexible endoscopes (Schott 1534180), which are manufactured by dissolving an acid soluble glass between the individual fibers in coherent imaging bundles. This produces an endoscope with thousands of dissociated fibers that come together in polished imaging surfaces at both ends [Fig. 1(b)]; each fiber has a core (diameter: $5.1\mu \mathrm{m}$), a cladding (thickness: $1\mu \mathrm{m}$), and the remnants of the acid soluble glass (thickness: $0.4\mu \mathrm{m}$). By cutting the bundle in half, the dissociated fibers are exposed [Fig. 1(c)].

Fig. 1

Illustration of the described method, using a bundle of optical microfibers as a multichannel, deep brain recording interface. (a) A standard fluorescent imaging configuration is used to interface with the fiber bundle. The polished imaging surface is mounted below an objective. Excitation light from an LED or other light source passes through a filter cube and is coupled into the fiber bundle; fluorescent emissions return through the objective, pass through the emission arm of the filter cube and are recorded by a camera. The fiber bundle consists of dissociated fibers, each with a diameter of $8\mu \mathrm{m}$. The individual fibers are implanted into the target brain region and secured to the skull. (b) The polished imaging surface as seen by the camera. (c) A bundle of 18,000 fibers prior to implant with gray silicone sheathing cut away. (d) Histology at the tip of a bundle implanted to 2.95 mm in a zebra finch, showing 530 optical microfibers (annotated in green). The fibers displace substantially less tissue (orange circle is cross section of displaced tissue) than existing optical techniques for recording from such a large volume of tissue.

After cutting away a fraction of the fibers in the bundle (based on the desired implant count), the remaining fibers are anchored together 4 to 5 mm from the tip using a bead of light-cured acrylic (Flow-It ALC, Pentron Clinical). The bead is secured to a digital manipulator in a stereotaxic rig. An anesthetized animal is mounted in the rig, and a craniotomy and durotonomy are opened above the target brain region. The fibers can then be slowly lowered into the tissue; the fibers have enough rigidity to enter the tissue without bending. Histology from implants shows that each fiber follows a path of least resistance, spreading through the target brain region [Fig. 1(d)].

The other end of the fiber (a ferrule and polished imaging surface) can be mounted below an objective in either a traditional fluorescence microscope or a purpose built optical configuration to provide excitation light and collect emission light via a CMOS sensor [Fig. 1(a)].

3.2.

Fiber Profile Model

In Fig. 2, we show a two-dimensional (2-D) slice of this 3-D photon distribution and compare it with a similarly generated profile of fluence in water (i.e., without negligible absorption and scattering³¹). The profile in water and tissue is similar, reflecting the fact that the relevant length scales are below the mean free path of light in brain tissue. That is, even in the tissue, photons undergo negligible scattering on the length scales relevant to the interface. In addition, we compare the two computationally generated profiles with an image of a single fiber in a fluorescein solution. Blue light is shone through the fiber, and the illuminated fluorescein is imaged through a 500- to 550-nm emission filter, revealing a distribution consistent with the model.

Fig. 2

(a) A single-fiber tip in a fluorescein solution with blue light (446 to 486 nm) emitted from the fiber. Due to the near total internal reflection, the light exits the tip of the fiber. The image, captured by a fluorescence microscope with a GFP filter, reveals the fluence excitation profile for the cut fiber. Grayscale image recolored to match modeled fluorescent profiles. (b) The linear fluence excitation profile of a single fiber, calculated via a Monte Carlo simulation of photon packets propagating through water. (c) The linear fluence excitation profile from a single fiber, calculated via a Monte Carlo simulation of photon packets propagating through the tissue. At each step, a photon packet may scatter or be absorbed based on tissue properties estimated for 490-nm light. Note that incorporating tissue scattering properties does not dramatically alter the profile, as the scale is below the mean free path of light in the brain. The dotted contours show the region that will receive $2.5\mathrm{mW}/{\mathrm{mm}}^2$ or greater excitation, given coupling 1.25, 2.5, or $5\mu \mathrm{W}$ of light into the fiber.

3.3.

Interface Model

Calculating the fluorescence sensitivity profiles for various distributions of optical microfibers, consistent with prior histology, we are able to model and evaluate the number and relative brightness of neurons contributing to the detected fluorescence signal. Figure 3 shows the modeled distribution of fibers and the relative brightness of individual neurons under different sets of conditions: two fiber densities (varying the number of fibers with the same standard deviation of splay) and two neural densities (reflecting different potential subpopulations). As the neural population is assumed to be uniformly distributed, increasing the density of neurons results in an increase in the number of neural traces being collected by the implant. Increasing the number of fibers has a more pronounced effect, as it increases the magnitude and uniformity of the excitation power delivered to the region and achieves a more dense sampling of fluorescence from the target neural population.

Fig. 3

The neural population we expect to interface with through a bundle of splayed optical microfibers. The spatial distribution of fibers (small pink dots) is based on the bivariate normal distribution seen in histology slices,¹⁵ and is shown in relation to the neurons (large circles) that are contributing fluorescence to the collected signal above a 1% threshold. Neurons are colored based on their fluorescence signal contribution to the fiber normalized by the maximum fluorescence signal that would be recorded when a cell is immediately under a fiber. Columns show two distinct fiber counts (with same splay parameter, $\sigma =150\mu \mathrm{m}$); rows show two distinct cell densities.

This relationship between the number of fibers and the number of neurons contributing fluorescence is explored further in Fig. 4, with each point based on the average of five randomly generated distributions of fibers. As the number of fibers increases, we initially observe a rapid increase in the number of neurons whose fluorescence contribution to a given fiber is above illustrative thresholds, due both to the increased sampling and to the more uniformly distributed excitation power. These effects eventually saturate as the distribution of fibers is kept fixed and all neurons within the field-of-view of the fixed distribution of fibers are contributing fluorescence signals.

Fig. 4

By increasing the number of fibers, while holding the splay diameter constant ($\sigma =150\mu \mathrm{m}$), the bundle can interface with more neurons. The visibility of a neuron is calculated based on their fluorescent contribution to a fiber normalized by the maximum fluorescence signal that would be recorded when a cell is immediately under a fiber. Each point is the average of five randomly generated fiber implant distributions. (a) The solid lines plot how many neurons are clearly visible (given different thresholds, represented as different color lines) to at least one fiber in the bundle. The dotted lines plot how many neurons are clearly visible to two or more fibers. (b) How many neurons are visible to a single fiber in the bundle. As the number of fibers increases, the excitation power increases and more neurons become visible.

These results demonstrate that the number of neurons contributing to the fluorescence signal increases more rapidly than the number of fibers. Under low densities of fibers, each fiber collects from a distinct set of neurons and the low-excitation power means that the signal will be limited to the 2 to 3 neurons closest to the fiber aperture. But as the fiber density increases, the paradigm shifts to one where many neurons are contributing fluorescence to multiple nearby fibers [(Fig. 4(b)]. As multiple neural traces are present in the signal from a single fiber, it becomes beneficial to be able to apply source separation techniques.

For source separation to be effective, there needs to be an overlap in the sensitivity profiles of fibers and, as a result, neurons contributing fluorescence to multiple fibers. The dashed lines in Fig. 4 show how many neurons contribute fluorescence to two or more fibers. As the number of fibers increases, this value also increases, suggesting that source separation becomes more feasible.

The breakdown of the fluorescence signals recording by an average fiber is shown in Fig. 5, depicting the relative contribution of the neurons that most strongly interface with the fiber. Each fiber collects the fluorescence from many neurons. Figure 5 shows that the brightest cell is on average 50% to 100% brighter than the next brightest cell contributing fluorescence to a given fiber. As the density of fibers increases (and, as a result, increases the uniformity of the excitation power), the drop off in brightness from subsequent neurons in the sorted list of neurons becomes less pronounced. We can look at this same breakdown from the other perspective: how many fibers capture the fluorescence from one neuron, as shown in Fig. 6. With a low number of fibers, the overwhelming majority of fluorescence from a neuron reaches a single fiber. As the density of fibers increases though, this shifts such that multiple fibers capture the fluorescence from a single neuron, permitting the application of source separation approaches.

Fig. 5

The contribution of individual neurons to the signal recorded through a single fiber. Values are normalized to the first (brightest) neuron visible to the fiber. Increasing the number of fibers does not substantially change the distribution, but delivers more light, which increases the brightness of the neurons further from the tip of the fiber. Each fiber captures fluorescence from a number of neurons, with an exponential drop off in the relative contribution of neurons. Each plot represents the average distribution of neural contributions for each fiber, averaged across fibers, and across five randomly sampled distributions of fibers. The four subfigures show increasing numbers of fibers, with a constant amount of splay ($\sigma =150\mu \mathrm{m}$).

Fig. 6

Each neuron primarily contributes fluorescence to a single fiber. Each box represents the strength of the interface between a given neuron and a fiber (round trip fluorescence), sorted in descending order of signal contributed to each fiber and normalized so that the strongest interface (the first fiber) is 100%. Plots are based on 5 sampled distributions of fibers, evaluating the 50 brightest neurons (most visible to a single fiber, in terms of round trip fluence). (a)–(d) show increasing number of fibers in a bundle, with a constant amount of splay ($\sigma =150\mu \mathrm{m}$).

3.4.

Source Separation

Given a higher density of fibers, we can apply the blind source separation technique described in Sec. 2 to approximate individual neural components that may be contributing to the recorded fluorescence. Figure 7 shows a toy model, created to be illustrative of the linear mixing of hypothetical neural traces (based on a Bernoulli spiking process convolved with the GCaMP6f waveform) and the subsequent results of the source separation technique.

Fig. 7

Simplified model used to exemplify the source separation process. (a) A simple 2-D configuration of fibers and neurons to demonstrate the source separation technique. Three fibers are shown in gray, each with contours representing their sensitivity profile. Three neurons (blue, red, and yellow) are positioned near the fibers while many neurons are distributed further from the fibers (creating background fluorescence). (b) Traces for the three neurons of interest generated as a random Bernoulli process convolved with a GCaMP6f waveform. (c) The signals recorded from the three fibers, representing a linear mixing of all the neurons based on the relative sensitivity. (d) The three outputs of the source separation process, sorted to most closely match the original neural traces.

In order to assess the performance of the source separation, the separated signals must be compared to all underlying neural activity to see if the extracted signals correspond with specific neural traces. We match extracted signals with underlying neural traces based on correlations; extracted signals having a sufficiently high correlation to the fluorescence of an underlying neuron ($r^2\ge 0.6$) are considered an accurate match. Figure 8 shows the percentage of extracted traces that are accurate matches with and without source separation, under a number of different model parameters (fiber counts and neural firing rate). Bundles with a larger number of fibers have a larger number of accurate matches, but calculating the percentage accuracy (normalizing by the number of fibers) better enables comparing performance across fiber counts. First, we assess accurate matches between the signals recorded via the fibers (without any source separation) and the underlying neural traces; as suggested by the earlier modeling, many fibers are dominated by 1 or 2 neurons and produce accurate matches without additional processing. By applying non-negative-independent component analysis, the number of accurate matches increases across all parameters. The increase is significant in all cases (paired $t$-test, $p$ $\le 0.01$). The benefit of source separation becomes more pronounced with an increasing number of fibers, consistent with our expectations given the greater overlap between neighboring fibers resulting from both the number of fibers and the increased uniformity of excitation power. As a control, we compare the fiber signals with newly-generated neural traces to confirm that the identified accurate matches are not simply a probabilistic result of the large neural population being modeled; with the newly-generated data, none of the traces accurately match (not pictured), indicating that matches are not false positives.

Fig. 8

Percentage of extracted signals that accurately match an underlying neuron, based on a correlation of $r^2\ge 0.6$. The first column (blue) shows the accuracy of the raw fiber intensity values; as some fibers are dominated by 1 to 2 neurons, there is a high percentage of traces that closely match underlying neurons. The second column (red) demonstrates the benefit of blind source separation, increasing the number of accurately extracted signals. These results are consistent and robust across various model parameters, including number of fibers and neural firing rates. In all cases, the source separation produces a significant increase in number of accurately matched traces (paired $t$-test, $p$ $\le 0.01$). Error bar shows std. dev. over 5 iterations.

Because this sort of evaluation may be sensitive to the $r^2$ threshold, we also evaluated the utility of the matched traces for different thresholds. A receiver operating characteristic curve (Fig. 9) evaluates the performance of using accurately matched, separated signals to detect action potentials (via threshold crossing) for a range of $r^2$ thresholds. As the $r^2$ accuracy threshold increases, the matched separated signals provide a more useful input for threshold detection, at the expense of reducing the number of separated signals. These results motivated the selected threshold of 0.6. As a control, the same threshold crossing detection was repeated comparing the separated signals to unrelated random traces. The AUC was between 0.5 to 0.51 for all $r^2$ thresholds (not pictured), consistent with the lack of a relationship between the ground truth and separated signals.

Fig. 9

By increasing the $r^2$ threshold used to match separated signals with underlying neural activity, we observe an increase in accuracy in spike detection based on threshold crossing for spike detection, but a decrease in the number of matched neuronal signals. By altering the spike detection threshold, we can generate a receiver operating characteristic curve measuring the accuracy of the spike detection algorithm. Using different $r^2$ thresholds to match separated signals with the simulated neural fluorescence traces, we can generate different ROC curves. By lowering the $r^2$ threshold, we match more traces ($N$), but sacrifice accuracy (measured as AUC). Average of five simulations based on 100 fibers with a splay of $\sigma =125\mu \mathrm{m}$ in a region with a cell density of 250,000 neurons per ${\mathrm{mm}}^3$.