In-vitro Recordings of Neural Magnetic Activity From the Auditory Brainstem Using Color Centers in Diamond: A Simulation Study

Magnetometry based on nitrogen-vacancy (NV) centers in diamond is a novel technique capable of measuring magnetic fields with high sensitivity and high spatial resolution. With the further advancements of these sensors, they may open up novel approaches for the 2D imaging of neural signals in vitro. In the present study, we investigate the feasibility of NV-based imaging by numerically simulating the magnetic signal from the auditory pathway of a rodent brainstem slice (ventral cochlear nucleus, VCN, to the medial trapezoid body, MNTB) as stimulated by both electric and optic stimulation. The resulting signal from these two stimulation methods are evaluated and compared. A realistic pathway model was created based on published data of the neural morphologies and channel dynamics of the globular bushy cells in the VCN and their axonal projections to the principal cells in the MNTB. The pathway dynamics in response to optic and electric stimulation and the emitted magnetic fields were estimated using the cable equation. For simulating the optic stimulation, the light distribution in brain tissue was numerically estimated and used to model the optogenetic neural excitation based on a four state channelrhodopsin-2 (ChR2) model. The corresponding heating was also estimated, using the bio-heat equation and was found to be low (<2°C) even at excessively strong optic signals. A peak magnetic field strength of ∼0.5 and ∼0.1 nT was calculated from the auditory brainstem pathway after electrical and optical stimulation, respectively. By increasing the stimulating light intensity four-fold (far exceeding commonly used intensities) the peak magnetic signal strength only increased to 0.2 nT. Thus, while optogenetic stimulation would be favorable to avoid artefacts in the recordings, electric stimulation achieves higher peak fields. The present simulation study predicts that high-resolution magnetic imaging of the action potentials traveling along the auditory brainstem pathway will only be possible for next generation NV sensors. However, the existing sensors already have sufficient sensitivity to support the magnetic sensing of cumulated neural signals sampled from larger parts of the pathway, which might be a promising intermediate step toward further maturing this novel technology.

Introduction

Negatively charged nitrogen-vacancy (NV) color centers in diamond are a novel tool to sense weak magnetic fields with nano- to millimeter spatial resolution, reaching sensitivity levels of pT/Hz^1/2 and below even at low frequency ranges and under ambient temperatures (Taylor et al., 2008; Fescenko et al., 2020; Zhang et al., 2021). For this reason the technique has been suggested to be well suited for in vitro imaging of the neural activity of cultured nerve cells or of groups of nerves in brain slices (Hall et al., 2012). The sensing protocols are based on optically detected magnetic resonance (ODMR) in which external magnetic fields cause changes in the optically induced NV fluorescence. This method has recently been used to detect the magnetic field of living cells (Le Sage et al., 2013) and large axons (Barry et al., 2016), as well as to characterize the current flow in graphene structure (Tetienne et al., 2017).

Traditional in vitro electrophysiological techniques allow recordings of the transmembrane voltage or current which helps to interpret the underlying mechanism of neural activity. However, their insufficient spatial sampling has led to the development of multi-electrode arrays (MEA) (Obien et al., 2015) or voltage-sensitive dye imaging (VSDI) (Stepan et al., 2015). VSDI, however, require complex addition of dyes which may not be suitable for all systems and has poor signal to noise ratio (Homma et al., 2009), and the spatial resolution of MEA is limited by the electrode density in the array (Obien et al., 2015). As higher spatiotemporal resolution and coverage is desired, NV-based magnetometry might offer an interesting, stain-free, alternative for the spatially resolved wide-field imaging of neural magnetic fields. Using NV-based magnetometry any sample which produces electric currents can in principle be imaged with high spatiotemporal resolution, without the need for genetic manipulation or staining. This method is completely passive and a sample can simply be placed on to a diamond sensor without any alterations to it, and can furthermore offer multi-modal capabilities, such as being used to sense electric fields (Dolde et al., 2011) and temperature (Acosta et al., 2010). For more detailed comparisons of NV-based magnetometry and existing alternatives for imaging neural signals please refer to e.g., (Schirhagl et al., 2014) and (Hall et al., 2012). We previously investigated the NV-based magnetometry method for imaging magnetic fields emerging from pyramidal cells in the gray matter of rodent hippocampal slices (Karadas et al., 2018). Complementing these findings, our aim here is to explore the feasibility of using NV-based magnetometry to also measure fast action potentials traveling along white matter pathways. Using detailed simulations, we analyze the electromagnetic fields generated by synchronous neural activity of the auditory pathway from the ventral cochlear nucleus (VCN) to the medial nucleus of the trapezoid body (MNTB) in mouse brainstem slices (Figure 1). This pathway provides a fruitful test case as (i) it is constituted by a dense bundle of hundreds of large-diameter axons which arise from the globular bushy cells (GBCs) in the VCN (Ford et al., 2015), (ii) the pathway can accommodate high frequency (∼200–400 Hz) firing (Taschenberger and von Gersdorff, 2000; Lorteije et al., 2009), thus allowing more averaging trials to increase signal to noise ratio (SNR), and (iii) the pathway is well suited to image action potentials in isolation without strong contributions from other neural signals. It is also a highly interesting target for neuroscience research due to the large size of GBC synaptic terminals, the Calyx of Held, which makes contact with principal cells from the MNTB. This is the only mammalian synapse allowing for simultaneous patch-clamp recording of the post- and pre-synaptic activities (Barnes-Davies and Forsythe, 1995; Borst et al., 1995). Most of the knowledge of pre-synaptic calcium dynamics and short-term synaptic plasticity comes from this synapse. Neurons from the MNTB project to the lateral superior olive (LSO), which integrates auditory inputs from both ears and plays a pivotal role for interaural discrimination (Joris and Trussell, 2018). Interestingly, both pre- and post-synaptic neurons from the Calyx of Held express the calcium binding protein parvalbumin (PV) (Felmy and Schneggenburger, 2004). Thus, by using PV as a promoter for Cre-dependent expression of light sensitive cation channel channelrhodopsin-2 in these neuronal populations, it should be possible to stimulate VCN and MNTB neurons with light.

FIGURE 1

Figure 1. (A) Auditory brainstem pathway from VCN to MNTB and the calyx of Held synapse. MNTB cells project further to the lateral superior olive (LSO). (B) A single calyx of Held synapse connected to the soma of a MNTB principal cell. (C) The terminal part of a GBC axon and its different sections.

In the present study, we simulate the extracellular magnetic fields created by neural activity traveling along this pathway in an in vitro preparation of a mouse brainstem slice. Our aim is to characterize the spatial and temporal features of the neural magnetic signals that would be sensed by an underlying planar NV sensor in order to assess the feasibility of this new imaging approach. Specifically, we model the dynamic responses of the GBCs and the principal cells of the MNTB to external stimulation and the synaptic transmission between both in NEURON (Hines and Carnevale, 2001) using previously reported morphological and biophysical characteristics (Lorteije et al., 2009). In addition to simulating the responses to the well-established electric stimulation, we explore the use of optogenetic stimulation as a means to prevent interferences of the stimulation with the magnetic measurements (Gunaydin et al., 2010; Deisseroth, 2011). Our simulations allow estimations of the expected magnetic field strengths as well as the required sensitivity level of the imaging system.

Materials and Methods

Cell Models of the GBC-MNTB Pathway

The axons of the globular bushy cells project from the VCN to the principal neurons in the contralateral MNTB (Figure 1A), where they end in large axo-somatic excitatory glutamatergic synapses that are termed calyces of Held (Figure 1B). The principal neurons from the MNTB project mainly to the lateral superior olive (LSO). Most GBC axons form only one calyx that connects to a single principal cell. Thus, there are one-to-one connections between globular bushy cells and their corresponding MNTB principal cells (Borst and Soria van Hoeve, 2012). In the following, we describe the kinetic and geometric models employed in our simulations.

Globular Bushy Cells (GBC)

Each globular bushy cell is modeled by five different components: Soma, dendrites, axon hillock, initial segment, and myelinated axon (Xie and Manis, 2013; Manis and Campagnola, 2018). The soma is a sphere with a diameter of 20 μm. A primary dendrite of 50 μm length and 3 μm diameter connects to the soma, and two secondary dendrites branch from the distal part of the primary dendrite to the right and left with lengths of 50 μm and diameters of 2 μm. The morphological data of the axon is taken from Ford et al. (2015). The axon hillock is oriented 90° with respect to the primary dendrite and has 15 μm length and a maximal diameter of 2.5 μm and then decreases to a diameter of 1 μm for the initial segment. The myelinated axon has a tapered structure, so that its diameter increases from 3 μm at the initial segment to ∼3.3 μm when close to the synaptic terminals. Similarly, the node diameters increase from 1.5 to 2 μm close to the synaptic terminals. The internode length initially remains constant at ∼200 μm, but then systematically decreases on the last ∼800 μm of the axon and reaches ∼40 μm close to the synaptic terminals. The length of all nodes is 1 μm and constant. These systematic variations in the axon morphology increase the maximum depolarization of the nodes close to the calyx of Held and aids transmission of action potentials. The calyx itself is modeled as a half ring (diameter of 13 μm, thickness of 2 μm), unless otherwise stated. Please see the next section for a detailed description of the axon pathway.

The cell membrane of the soma, primary dendrites and initial segment contain inactivating Na⁺ channels, low-threshold K⁺ (K_LT) and high-threshold K⁺ (K_HT) channels, and hyperpolarization-activated cation (I_h) channels with the kinetics taken from Rothman and Manis (2003). The corresponding relations between the ionic currents and their kinetics, including the channel properties, are given in Supplementary Paragraph 1.

In order to model the cell membrane of the internodes of the axon, the passive myelinated axon model C from Richardson et al. (2000) was implemented in NEURON. The periaxonal space is not included in the model and it was further assumed that the myelin sheath membrane has constant conductivity (Supplementary Paragraph 2 contains detailed information). The nodes of Ranvier of the axon are modelled as containing inactivating Na⁺ and low-threshold K⁺ (KLT) channels with Hodgkin–Huxley-like kinetics as well as a leak conductance (Rothman and Manis, 2003). The last internode is followed by a longer node called heminode, which has the same channel types and densities as the other nodes of Ranvier (Ford et al., 2015). The heminode connects to the calyx of Held via a short internode named post-heminode, which we included in the modal as a capacitor without any conductance.

The membrane of the calyx of Held is modelled as having high-threshold K⁺ (KHT) and hyperpolarization-activated cation (h) channels and a voltage sensitive calcium conductance (Figure 1B; Supplementary Paragraph 1 contains details) (Borst and Sakmann, 1998).

Synaptic connections

The synaptic transmission model consists of two kinetics: Vesicle release and AMPA receptor. The vesicle release mechanism is modeled by 500 independent active release zones each with one readily releasable vesicle (Graham et al., 2001). The mechanism is activated by arrival of an action potential (AP) to the calyx of Held. It initiates fast and slow calcium transients for 1 and 2 ms, respectively. The fast calcium transient activates opening of the fast and slow gates. The fast gates go to open state quickly upon arrival of the AP and close rapidly. On the other hand, slow gates open partially and close slowly so that a potentiation can occur by a following AP. A release probability for all active zones is determined by the product of open states of slow and fast gates. Once a vesicle is released, a 1 ms square pulse of neurotransmitter T is produced that drives a corresponding AMPA receptor-mediated EPSC. The receptors are placed at the soma and partially at the axon initial segment of the principal cell. A 6-state gating model is used to simulate the AMPA receptor kinetics, composed of three closed states, fast and slow open states, and one desensitized state (Raman and Trussell, 1992). At the same time, the vesicles are replenished by a constant background replenishment from a large reserve pool and enhanced replenishment driven by the slow (residual) calcium transient. Further details on the kinetic model and parameters can be found in Supplementary Section 4.

MNTB Principal Cells

A principal neuron has a simple spherical perikaryon, one or two short and thin dendritic trees (Smith et al., 1998), and its axon provides inhibitory input to the LSO. Here, it is modeled by three compartments: Soma, dendrite, and axon. The soma is a sphere with a diameter of 20 μm. The dendrite consists of cylinders of 50 μm length and 3 μm diameter and makes 2 branches of 10 μm at the end. The axon is oriented 90° with respect to the dendrite and composed of 5 internodes and nodes. Each internode and node have constant diameters of 2 and 1 μm, respectively, and lengths of 100 and 1 μm. As for the axonal pathways of the GB cells (see next section), the axon pathway is generated by a spline interpolation of randomly generated five points on the way to the LSO. The model has Hodgkin–Huxley-like dynamics with a high density of Na⁺, low-threshold K⁺ and high-threshold K⁺ channels, similar to a previous model of the MNTB neuron (Wang et al., 1998). Whereas low-threshold K⁺ channels dampen excitability so that the calyceal inputs generates a single AP, high-threshold K⁺ channels guarantee rapid gating and shorten the AP (Brew and Forsythe, 1995). Details on the ion channel dynamics of the principal cells can be found in Supplementary Section 3.

Morphological Model of the GBC-MNTB Axonal Pathway

The white matter pathway from VCN to MNTB is reconstructed from the confocal images of transverse brainstem slices given in Ford et al. (2015) and Kronander et al. (2017). Here, the thickness of the brainstem slice is modelled as 350 μm, consisting of a bottom layer of 25 μm of inactive cells (affected by the cutting procedure), a 300 μm thick layer of active cells and a top layer of 25 μm of inactive cells (Figure 2B). For each globular bushy cell, the axonal path is determined by spline interpolation of eight points as follows:

FIGURE 2

Figure 2. (A) Four-state model, with two open and two closed states, of channelrhodopsin-2 (ChR2) as described in Nikolic et al. (2009). e12, e21, Kd1, Kd2, Kr are the rate constants (Supplementary Table 5) and Ka1 and Ka2 are the light dependent constants. (B) Schematic illustration of the simulated brainstem slice placed on the diamond sample. It is assumed that the neural cells in a distance of up to 25 μm to the diamond are dysfunctional due to the preparation. The changes in photoluminescence emitted by the NV-layers in the diamond, caused by the neural magnetic fields, can be recorded using an arrangement like an inverted microscope, with the camera replacing the photodetector (Wojciechowski et al., 2018). In optical stimulation, a light probe is placed on the MNTB (drawing is not to scale). (C) Three different illumination positions along at the pathway were simulated (VCN region containing the soma of the GB cells, the myelinated axons of the GB cells and the MNTB region) to test the feasibility of light excitation.

(i). The MNTB area is defined as an elliptical region (semi-major axis 200 μm and semi-minor axis 100 μm, active cell thickness 300 μm) which is 400 μm away from the brain midline (Weatherstone et al., 2017). A random position (x_c, y_c, z_c) is chosen in this area as the center location of the calyx of Held. This position is also the center of the soma of the corresponding principal cell.

(ii). Seven rectangular areas (each with 25 × 25 μm²) along the axonal pathway were manually predefined. A random position is picked from within each area, and the eighth point is added at (x_c, 0, 0) for the calyx of Held.

(iii). The axonal pathway is generated by spline interpolation of the eight points, and the complete pathway including the GBC soma is finally moved by y_c in y-direction and by z_c along z-direction.

(iv). The soma of the principal cell is placed at (x_c, y_c, z_c) and rotated around its z-axis by a random angle between -π/2 and π/2. The calyx of Held is also rotated randomly around its z-axis.

In order to model a brainstem stem slice with a thickness of 300 μm, in total 300 cells are created. In the literature, the total number of principal cells is given rather than the neural density, and the average number of cells per MNTB area is estimated to be around 2,500 in mice (Weatherstone et al., 2017). Considering 8–10 slices per brainstem and a total height of the VCN and MNTB area of around 3,000 μm, we assume to have around 300 cells in a 300 μm thick slice through the central part of the MNTB area as a worst case estimate. However, it is more likely that all the calyces present within 600 μm which results in ∼1,250 cells per slices (Kronander et al., 2017). We therefore report both the results for the worst case and the more likely case in the following.

The neuronal density also varies for different slices and species. Therefore, our results should be linearly rescaled to match those for a different neuron density.

Calculation of the Neural Magnetic Field and Electric Potential

It is known that the axial and transmembrane currents are the major sources of the magnetic field and electric potentials, respectively (Woosley et al., 1985). Therefore, we follow the two step approach given in Karadas et al. (2018) to determine the extracellular neural magnetic fields and electric potential: Initially, the membrane potentials V_mⁿ and transmembrane currents I_mⁿ of the simulated neuron are calculated using the NEURON software package (v7.4) (Hines and Carnevale, 2001) with a time resolution of 25 μs. Then, the extracellular magnetic fields B and the electric potential φ at the sensor surface are determined from the simulations using a forward modeling scheme for the simplified case of an unbounded homogenous volume conductor, as described in Karadas et al. (2018). Finally, the impact which the extracellular current flow in the non-simplified volume conductor exerts on the measured magnetic field is accounted for using a scaling factor as given in Karadas et al. (2018). This scaling factor approximates the effects of the conductivity anisotropy in the brain slice, the boundary between the brain slice and the underlying sensor, and the fluid surrounding the brain slice on the extracellular current flow.

Modeling of Optogenetic Stimulation

Optogenetic stimulation might be a good alternative to the well-established electric stimulation to avoid interferences of the stimulation with the magnetic measurements. To investigate the efficiency of optogenetic stimulation to generate APs across all axons of the pathway, we implemented the “light-neuron” model given in Foutz et al. (2012). The model simulates the neural dynamic in brain tissue generated by an optical probe which couples to light-sensitive neurons. The latter were modeled using the cable equation, extended by Channelrhodopsin-2 (ChR2) channels. In the following, we describe the ChR2 channel model and the light distribution in tissue created by an optical probe. In addition, the model used to estimate the thermal effects due to light absorption is outlined.

ChR2 Photo-Current Modeling

Channelrhodopsin-2 (ChR2) was modeled with two closed (C₁, C₂) and two open (O₁, O₂) channel states (Figure 2A; Nikolic et al., 2009; Foutz et al., 2012). The initial state of a channel is C₁. Under illumination, C₁transits into open state O₁ at a rate constant of K_a1. A channel in O₁ state exponentially decays to closed state C₂ at a rate K_d1, or instead transits to the second open state O₂ at a rate constant of e₁₂. State O₂ is more stable and has lower ion conductance with respect to O₁. The O₂ state can exponentially decay to C₂ at a rate K_d1 or return to O₁ at a rate constant of e₂₁. Channels being in C₂ state can go into O₂ state at an exponential increase rate K_a2,or they go to C₁ state by a slow thermal conversion with exponential decay rate K_r. The instantaneous rate of change of these states was defined by a set of rate equations (Grossman et al., 2011b):

$$\frac{d{O}_1}{dt}=K_{a1}{C}_1-\left(K_{d1}+e_{12}\right)O_1+e_{21}{O}_2(1)$$

$$\frac{d{O}_2}{dt}=K_{a2}{C}_2+e_{12}{O}_1-\left(K_{d2}+e_{21}\right)O_2(2)$$

$$\frac{d{C}_2}{dt}=K_{d2}{O}_2-\left(K_{a2}+K_r\right)C_2(3)$$

$$C_1+C_2+O_1+O_2=1(4)$$

In this state model, only the rate constants K_a1 and K_a2 that control the transition from the closed to the open states are dynamically changed by the light stimulation:

$${\text{K}}_{\text{ai}}=\{\begin{array}{cc}{\epsilon }_{\text{i}}\mathrm{\varphi }\left(1-{\text{e}}^{-\frac{\text{t}-{\text{t}}_{\text{ON}}}{\mathrm{\tau }}}\right)\mathit{\hspace{1em}}\hfill & \mathrm{during}\mathrm{illumination}\hfill \\ {\epsilon }_{\text{i}}\mathrm{\varphi }\left({\text{e}}^{-\frac{\text{t}-{\text{t}}_{\text{OFF}}}{\mathrm{\tau }}}-{\text{e}}^{-\frac{\text{t}-{\text{t}}_{\text{ON}}}{\mathrm{\tau }}}\right)\mathit{\hspace{1em}}\hfill & \mathrm{after}\mathrm{illumination}\hfill \end{array}(5)$$

Here, ε_i is the quantum efficiency of channelrhodopsin, ϕ is the photon flux per area during illumination. Once the irradiance is determined for each ChR2 molecule, ϕ is calculated by considering the cross-sectional area of the molecule (∼1.2 × 10^–8 μm²) (Hegemann et al., 2005) and the wavelength of the light (473 nm). τ is the time constant of the channel. t_ON and t_OFF are the time points at which illumination is switched on and off. Once O₁ or O₂ are reached, the photocurrent through the membrane for illumination ON is calculated as

$$I_{ChR2}=i_{max}\left(O_1+\mathrm{\gamma }{O}_2\right)\frac{\left(1-e^{-U/U_0}\right)}{U/U_1},(6)$$

where i_max = (V−E_ChR2)g_ChR2 is the maximum reachable ChR2 current depending on the reversal potential E_ChR2 (set to 0 mV) and the channel conductance per area (g_ChR2 = g₁ρ_ChR2). Constant γ (set to 0.05) is the conductance ratio between O₁ and O₂. The last term in (6) accounts for an inwardly rectifying current–voltage curve as shown experimentally in Grossman et al. (2011b), where U is the absolute transmembrane potential. U₀ and U₁ are empirical constants set to 40 and 15 mV, respectively.

After the illumination goes off, the ChR2 current shows an exponential decay which is explained by a fast (i_fast) and slow (i_slow) component as given by

$$i=i_{slow}{e}^{-{\mathrm{\Lambda }}_1\left(t-t_{OFF}\right)}+i_{fast}{e}^{-{\mathrm{\Lambda }}_2\left(t-t_{OFF}\right)},(7)$$

where Λ₁ and Λ₂ are decay factors. The complete set of equations defining the ChR2 dynamics and parameters are provided in Supplementary Section 5.

Light Model

There are four main factors defining the light distribution in tissue: (1) The light distribution at the source, (2) the geometric spread, (3) the scattering, and (4) the absorbance by the medium. The most established approaches to determine the light distribution in tissue are the Kubelka-Munk (KM) model and Monte Carlo (MC) simulations. The KM approach provides an analytical expression so that the irradiance can be determined at any point and easily incorporated as input into the NEURON solver. In the KM model, the light intensity is described by the transmitted and the backscattered components which are defined by two coupled differential equations (Mobley et al., 2014). The light beam in the forward direction decreases in intensity due to the absorption and scattering processes and gains intensity from the backscattering process from photons being scattered in the reverse direction. In a cylindrical coordinate system, the irradiance I at a point is defined by I(r,z) = T(r,z) I₀, where I₀ is the irradiance of the source and T(r,z) is the transmittance at the radial distance r and depth z. The transmittance factor is wavelength dependent and can be represented by three linear components (Foutz et al., 2012)

$$T(r,z)=G(r,z)C(z)M(r,z),(8)$$

where G describes the distribution of light emitted by the light source, C describes the geometrical spreading of the light, and M describes the scattering and the absorbance of the light according to the KM model. G(r,z) is generally defined by a Gaussian function which takes into the account expansion from the source. The geometric spread function C(z) is calculated by the law of conservation of radiant power. Therefore, as the light travels it will diverge and its intensity decreases. For the radius R(z) of the light cone at depth of z starting from R₀ emission, the spread function is:

$$C\left(z\right)=\frac{R_0}{R(z)}(9)$$

The solution for factor M is given by

$$M(x)=\frac{b}{asinh\left(b{\mathrm{\mu }}_{s}x\right)+bcosh\left(b{\mathrm{\mu }}_{s}x\right)}(10)$$

with

$$a=1+\frac{{\mathrm{\mu }}_a}{{\mathrm{\mu }}_s}\mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}b\mathit{\hspace{1em}}\sqrt{a^2-1},(11)$$

where the parameter $x=\sqrt{z^2+r^2}$ is the travelled distance, μ_s is the scatter coefficient per unit thickness (in [mm^–1]), and μ_a is the absorption coefficient per unit thickness (in [mm^–1]). The KM model relies on the simplifying assumptions of (1) isotropic scattering, which is valid in the diffusive regime for depths of multiple millimeters, and (2) isotropic illumination, thus neglecting the finite geometry and size of the illumination. In contrast, in optogenetics, scattering is highly anisotropic at the relevant distance scales and wavelengths and the light source typically has a comparable size to the tissue of interest (Yona et al., 2016). Here, we implemented the KM model using custom-written Python scripts.

For this reason, we additionally simulated the light distribution using the MC method, where the random walk of photon packets is traced through the tissue in three dimensions (Wang et al., 1995). In this method, photon packets are launched from a light source so that each packet has a random initial position and direction according to the geometry and numerical aperture of the source. Then, the light transport through the tissue is simulated as iterative process. Initially, a weight of 1 is given to each photon packet and it is moved by a random distance depending on the absorption and scattering coefficients. After each iteration, the photon packet loses weight by a factor of μ_a/(μ_a + μ_s), and this fraction is added to the total amount of light absorbed at the corresponding tissue position. The new trajectory is defined by an anisotropic scattering model using the Henyey–Greenstein (HG) phase function. This procedure is iterated for each photon packet until its weight decreases below a certain threshold. The details of the implementation are given in Jacques (2010). Here, we used the MC implementation that was made available by these authors.

In addition to these well-established techniques, another analytical approach based on the beam spread function (BSF) has been proposed (Yona et al., 2016). For validation of our light simulation, we also determined the light distribution predicted by BSF for same probe geometry and tissue parameters. The BSF simulations were done in MATLAB 2018a (The MathWorks, Natick, MA, United States) using the ScatterBrain (v0.8) software. The results of the validations are included in the Supplementary Material (Supplementary Figures 6, 7). Shortly, while the estimations obtained with the KM and BSF models agree well to each other and to existing experimental data, the MC simulations clearly overestimated the irradiance levels close to the light probe since there is almost zero scattering and absorption in the initial ballistic regime of photon propagation (for gray matter ∼140 μm) (Ntziachristos, 2010). For this reason, we chose the KM model for the rest of the study, unless indicated otherwise.

Optothermal Effect

Light absorption by tissue generates heat. The amount of absorbed energy is defined by the source term, which is the product of the light irradiance ϕ(mW/mm²) and the absorption coefficient μ_a. Here, we describe the resulting temperature increase by the well-established Pennes bioheat equation (Pennes, 1948), which we simplify by ignoring the convection term due to blood/fluid perfusion:

$$\mathrm{\rho }c\frac{\partial T}{\partial t}=\nabla \left(k\nabla T\right)+\mathrm{\varphi }{\mathrm{\mu }}_a.(12)$$

Here, the constant ρ denotes the mass density (ρ = 1.04 × 10³ kg/m³), c the specific heat capacitance (c = 3.6 J/gK) and k the thermal conductivity (k = 0.56 W/m K) of the brain tissue. The forward finite difference method is used to solve Eq. 12, using a spatial step width of Δx = 0.03 mm, and ensuring numerical stability by adjusting the temporal step as given in Stujenske et al. (2015) to

$$\mathrm{\Delta }t\le \frac{{\left(\mathrm{\Delta }x\right)}^2}{6k}\mathrm{\rho }c.(13)$$

To estimate the heating effect, we used the light irradiance estimates given by the Monte Carlo simulations, as they were higher than the values given by the KM model and were therefore selected as worst-case scenario.

Results

Validation of the Modelled Cell Dynamics

We confirmed that the model cells exhibited the expected dynamics using morphologically simplified models of a globular bushy cell connected to a principal cell (Figure 3A). First, we tested the responses of each individual cell to repetitive current pulses applied to the soma region. Consistent with electrophysiological data (Manis and Campagnola, 2018), the GB cell generates reliable firing for frequencies up to 200 Hz (Figure 3B). Similar results are obtained for the principal cell of the MNTB region (Brew and Forsythe, 1995; Macica et al., 2003; Figure 3C).

FIGURE 3

Figure 3. Validation of the modelled cell dynamics. (A) Morphologically simplified model of a GB cell connected to a principal cell. (B) GB cell response to current injections (0.25 ms, 20 nA) for two different stimulating frequencies (100 and 200 Hz) at the soma. (C) MNTB principal cell response to current injections (0.25 ms, 0.5 nA) for two different stimulating frequency (100 and 200 Hz) at the soma. (D) Synaptic transmission from the calyx of Held to the soma of the MNTB principal cell. An action potential and its negative time derivative are shown together with the resulting EPSC at the soma of the principal cell to indicate the synaptic delay. The soma of the MNTB principal cell was clamped to –70 mV to evoke the AP. (E) Trains of current pulses with frequencies of 100 and 200 Hz, respectively, are applied to the soma of the GB cell and the resulting EPSCs are determined at the soma of the principal cell. Both frequencies reliably generate APs at the calyx. (F) EPSCs normalized to the peak EPSC show initial facilitation followed by depression.

In order to verify the modelled synaptic transmission at the calyx, Figure 3D illustrates the synaptic currents under a voltage clamp with a holding level of −70 mV applied to the soma of the principal cell. The synaptic delay has a good match with experimental data obtained in rat brain slices for postnatal days 14 to 17 (Fedchyshyn and Wang, 2007). In addition, trains of current pulses with frequencies of 100 and 200 Hz, respectively, were applied to the GBC soma and the EPSCs were measured at the soma of the principal cell. Both cases generated highly reliable trains of APs at the calyx (Figure 3E). Since the calyx has 500 active transmission zones, it causes post-synaptic EPSCs larger than 1 nA. These EPSCs are strong enough to guarantee a suprathreshold EPSP for each pre-synaptic AP (Forsythe, 1994). For the chosen calcium concentration, the EPSC trains show an initial facilitation followed by a depression response (Figure 3F), in line with experimental data (Barnes-Davies and Forsythe, 1995). Continuing depression and facilitation, respectively, of the EPSCs could be observed when changing the extracellular calcium concentration correspondingly (data not shown) (Graham et al., 2001). In summary, the transmission model adequately represents paired-pulse facilitation, depression and AMPA receptor desensitization, as demonstrated in Barnes-Davies and Forsythe (1995).

Feasibility of Light Excitation of a Brainstem Slice

To study the feasibility of light stimulation as an alternative to the electrical stimulation, we systematically varied the intensity and duration of the light pulses to determine the irradiance threshold required to generate an action potential for three different parts of the pathway (as can be seen in Figure 2C), namely the soma of the GB cells, the calyx of Held and the soma of the MNTB principal cells (Figure 4A). The membrane potential at the 7th node of Ranvier, counted from the initial segment, was used to determine whether a cell fired or not. The irradiance threshold shows the typical dependence on pulse duration, as also commonly observed for electrical stimulation. The rheobase irradiance levels are ∼140, ∼24, and ∼1 mW/mm² for the GBC soma, the calyx and the MNTB principal cell soma, respectively. The corresponding chronaxies are ∼1.3, ∼2.0, and ∼3.5 ms. The calyx is more easily excitable than the GBC soma since it can accommodate more light excited currents due to its large surface area and higher calcium channel density. However, illumination of the calyx causes also independent firing of the MNTB principal cells (Figure 4B), as they have even lower thresholds. We also tested light stimulation of the myelinated GBC axons. However, action potentials could not be triggered even at very high (>10 W/mm²) illumination intensities. As the nodes of Ranvier have small surface areas, even fully saturated ChR2 currents are not enough to sufficiently depolarize the membrane.

FIGURE 4

Figure 4. (A) Dependencies of the irradiance threshold on the duration of the light pulse for three different parts (soma of GB cells, Calyx of Held, soma of MNTB principal cells) of the pathway. (B) Illumination of the calyx of Held activates light sensitive channels of the MNTB principal cells and causes firing independent of synaptic transmission. Since the irradiance threshold of MNTB cells is much lower, they already spike under far weaker illumination. (C) Firing frequencies in response to repetitive light stimulation (200 ms total duration, 3 ms pulses) for the three parts of the pathway. Repetition frequencies of the light pulses range from 10 to 60 Hz, and two different irradiance levels were tested. (D) Total stimulation duration is increased to 1 s, with the other parameters held identical to (C).

As long as magnetometer sensitivities are limited, it might be necessary to average the neural responses to many successive pulses in order to recover a sufficient SNR level. We therefore determined the firing frequencies of the modelled neurons in response to repetitive light stimuli, with varying the repetition frequencies from 10 to 100 Hz. The total duration of the pulse train was 200 ms, a single pulse had a duration of 3 ms, and illumination strength was chosen to be 1.5 and 3 times of the irradiance threshold (Figure 4C). The responses became quickly unstable for the lower illumination strength when increasing stimulation frequency. At the higher illumination strength, the responses became unstable for frequencies exceeding 15 Hz when stimulating the GBC soma. In contrast, stimulation of the calyx gave stable responses for frequencies up to 40 Hz. When increasing the duration of the pulse train to 1 s (Figure 4D), reliable stimulation of the GBC soma was not possible, while calyx stimulation was still stable up to a frequency of 30 Hz. The soma of the MNTB cells could accommodate even higher frequencies at both pulse train durations. However, as this activation does not result in action potentials traveling along the GBC axon pathway, frequencies exceeding 40 Hz might be of less practical relevance. On the other hand, varying between low and high frequency stimulation might give the possibility to differentiate between the relative contributions of the calyx and the soma of the principal cells to neural recordings. In general, these results suggest that optogenetic stimulation requires the use of lower stimulation frequencies than electric stimulation, as the slow closing mechanism of ChR2 channels causes a longer spiking latency (Grossman et al., 2011a). When firing in synchrony, the neural magnetic fields of the stimulated cells add up, increasing the signal recorded by the magnetometer. As this effect depends on the exact temporal alignment of the single fields, we investigated the timing of action potentials generated by 300 GB cells distributed in a 300 μm thick slice under three different lighting conditions (1 W/mm² light source irradiance and 3 ms duration; 4 W/mm² and 3 ms; 4 W/mm² and 5 ms). The light distribution in the tissue was simulated with the KM model. We selected a light irradiance at the source of 1 mW/mm² as this is at the upper limit for the available optogenetic solutions using “fiber-coupled” LED light sources, which can output a power of up to 30 mW (Prizmatix, 2019). This irradiance was increased four-fold to 4 mW/mm² as a test for an extreme case. The light probe (optical fiber) was simulated 25 μm away from the slice and centered above the MNTB area to excite the calyces. The diameter was set to 200 μm and the NA to 0.37. For stimulation with 1 W/mm² and 3 ms, 19% of the GB cells were antidromically activated in the complete slice. Mean firing times were 2.18 ms (±0.37 ms standard deviation, SD). Increasing the light irradiance at the source to 4 W/mm² (still at 3 ms duration) resulted in an activation of 79% of the cells with mean firing times of 1.50 ± 0.35 ms SD (Figure 5B). The spatial distribution of activated and inactivated neurons using different light stimulation methods can be seen in Supplementary Figure 9. Using longer pulses of 5 ms gave similar results (Figure 5C). These results indicate that optical stimulation is limited in its ability to activate all cells in high synchronicity, which will result in weaker and temporally more dispersed signals recorded by the magnetometer when compared to the signals caused by electrical stimulation.

FIGURE 5

Figure 5. (A) Spike times of the 300 cells in the brainstem slice under an illumination with an irradiance level of 1 W/mm² and 3 ms duration. Please see the main text for details on the light probe. Inactive cells are plotted as gray dots, activated as black dots. (B) The irradiance level at the source is increased to 4 W/mm², causing a clear increase in the activation rate. (C) The duration of the stimulus pulse is increased to 5 ms for an irradiance of 4 W/mm², still giving similar results as seen for the shorter pulse.

Spatial-Temporal Characteristics of the Neural Magnetic Activity of a Brainstem Slice

Our main interest was the estimation of the spatial-temporal characteristics of the neural magnetic fields that are caused by the stimulation of the auditory pathway in a transverse brainstem slice. An active cell layer thickness of 300 μm (Figure 6A) and a temperature of 34°C was chosen for the simulations. The distance from the diamond sensor to the active cell layer was set to 25 μm to account for the layer of inactive cells due to the cutting procedure. The magnetic field and the electric potential caused by the neural activity was calculated at a surface area with a size of 2 × 1 mm². The surface area was discretized into a 50 × 25 grid, resulting in an in-plane resolution of 40 μm. Considering a conduction velocity of 5 to 10 m/s and a depolarization time of 0.1 ms, such spatial sampling is considered sufficiently high to accurately capture the spatial features of the magnetic fields.

FIGURE 6

Figure 6. Simulation results showing the spatial distributions of the extracellular electric potential and neural magnetic fields of the auditory GBC-MNTB pathway in response to electrical and optical stimulation of the MNTB region. (A) The simulated pathway consists of 300 GB cells distributed in a 300 μm thick slice, connected to 300 MNTB principal cells. The stimulated MNTB region is indicated by the dashed ellipse. (B) Equivalent current dipole moment. (C) LFP and magnetic fields at t = 0.2 ms after the start of the electrical pulse. The field distributions are in accordance with an axial current flow from the axonal pathway to the MNTB region, caused by action potentials traveling antidromically along the pathway. (D) Magnetic fields at t = 1.88 ms after start of the light pulse. The spatial distributions around the MNTB region differ from those caused by electrical stimulation. Having lower light excitation thresholds, more MNTB principal cells than GB cells are activated, which makes the contribution of the former more dominant (see Supplementary Figure 3A). Still, the overall lower number of activated cells and the temporal dispersion of the activations decreases the peak field strength by a factor of approx. five.

Given that the calyces of Held and the soma of the principal cells had lowest irradiance thresholds, we simulated the electrical and optical stimulation of the MNTB region in the following. Our aim was to characterize the contributions of the neural activity at the different parts of the pathway to the overall magnetic fields that can be recorded at the surface of the brain slice. In the electrical stimulation case, the stimulus was chosen to have a duration of 0.1 ms. As idealized reference case, we used intracellular current injections that were placed at the calyces of Held and their amplitudes set to 5 nA, which is strong enough to generate synchronous events at all calyces of Held. In addition, to evaluate the feasibility of the generation of synchronous events by extracellular stimulation, a monopolar electrode was simulated at a depth of 200 μm of the slice and centered at the MNTB area to excite the calyces (details on the implementation can be found in Supplementary Section 6). In the optical stimulation case, a light probe with an NA of 0.37 and a diameter of 0.2 mm was simulated on the top of the slice at a distance of 25 μm. The source irradiance at the probe was set to 1 W/mm² (∼30 mW light power) and the stimulus had a duration of 3 ms.

Since neural magnetic fields strongly depend on the neural axial current densities (Swinney and Wikswo, 1980; Karadas et al., 2018), we also determined an equivalent current dipole (ECD) by summing the resulting neural currents (Figure 6B). The ECD is given by $\overrightarrow{\text{Q}}{\sum }_{\text{k}}\text{I}_{\text{i}}{}^{\text{k}}{\overrightarrow{\text{L}}}^{\text{k}}$, with ${\overrightarrow{L}}^k$ indicating the length of the k^th cylindrical compartment and $I_i^k$ the axial current flow within the compartment at the time frame i. The simulated ECD reached its peak strengths 0.20, 0.15, and 1.88 ms after start of the intracellular electrical, extracellular electrical and optical stimulation, respectively. Therefore, we determined the spatial distribution of the neural magnetic fields and local field potentials (LFPs) at these time points. The results for intracellular electrical stimulation are shown in Figure 6C. At the chosen time point, the GB cells are depolarized, and a strong axial current flow is started from the myelinated pathway at the left side of the midline to the calyces of Held. As a result, the LFP has a dipole distribution centered around the midline of the brain slice. Determined by the orientation of the axon pathway at this position, the neural source currents are predominantly oriented in x direction. The spatial pattern of the resulting magnetic field is in accordance to this direction, with the Bx component being almost zero. The magnitude of the By and Bz components reach peak values up to 120 pT at the distance of 25 μm for the worst-case estimate of 300 cells in the slice. Considering a more likely estimate of 1,250 cells in the slice, the peak signal reaches ∼500 pT (results obtained by linear rescaling of the field of the 300 cells). In general, axial currents in the myelinated axon pathway contribute most to the magnetic fields. For example, selecting a later time point at which only the soma of the MNTB principal cells are activated due to synaptic transmission (delay: ∼0.6–0.7 ms) results in almost four times lower peak fields (Supplementary Figure 4A). Likewise, the antidromic activation of the soma of the GB cells causes only weak fields (Supplementary Figure 3).

For extracellular rectangular stimuli with 100 μA and 100 μs, 82% of the GB cells are antidromically activated in the complete slice. Considering a titanium nitrite (TiN) electrode with 10 μm diameter, this level of current strength and duration is below the safe-charge injection limits to avoid damage to the electrode and cells (Multi Channels Systems MCS GmbH, 2018). Changing from intra- to extracellular stimulation only weakly affects the synchrony of the generated APs. For extracellular stimulation, the APs were mostly generated by activation of the 1st or 2nd node (counted from the calyx). The firing times of the (arbitrarily selected) the 6th node of the axon, again counted from the calyx, were 0.09 ms with a low standard deviation (SD) of ± 0.04 ms after stimulus onset, demonstrating the good synchrony of the APs. As a result, the magnetic peak fields seem to be only moderately reduced to ∼100 pT for our worst-case estimate of 300 cells in the slice (Supplementary Figure 12). Similar magnetic peak fields are obtained when an electrode is placed in the VCN region close to the soma of the GB cells, where ∼80% of the cells are activated with the same stimulation parameters via triggering APs at the ends of the axon initial segments (data not shown). Activation of 100% of the cells could be obtained when increasing the stimulus level to 200 μA, which, however, might result in damage of the electrode tip or the surrounding cells (Multi Channels Systems MCS GmbH, 2018) (data not shown).

In the optical stimulation case, the peak strengths of the fields decrease approximately fivefold (Figure 6D), in accordance to the weaker and temporally broadened ECD. This effect results from the temporal asynchrony of the generated action potentials and low number of excited cells (Figure 5A), which also causes a broadening of the LFP and the magnetic field of GB cells along x direction (i.e., along the main direction of the axon pathway; Supplementary Figure 4B). In addition, the magnitudes of the fields at the MNTB region are higher than those at the axon pathway and dominate the overall magnetic field in case of optical stimulation (Figure 6B). This is because the contribution of MNTB cells to the recorded magnetic fields is almost two times higher than that of the GB cells (Supplementary Figure 4B) due to their lower excitation thresholds.

Effect of Calyx Morphology

We aimed to ensure that the rather weak magnetic fields observed at early time points after stimulation were not biased by the chosen morphological models for the calyces of Held which were represented as half rings. Therefore, we tested three different variants for these models with increasing morphological complexity (Supplementary Figure 6A). The second variant modelled the calyx as an additional half ring perpendicular to the first one with two simple sticks of lengths of 6 μm, placed at random orientations. The third variant successively added additional sub-branches to reach better approximations of the real calyx morphology. However, the resulting magnetic field was only weakly affected by the model complexity (Supplementary Figures 6E,F).

Heating Effect of Light Excitation

We wanted to ensure that tissue heating due to optical absorption does not pose a problem when applying many light pulses at high repetition rates in order to be able to obtain averaged data with sufficient SNR. We thus estimated the light distribution in brain tissue and used this as input to the bioheat equation to estimate the temperature change. In order to obtain robust worst-case estimates, we used Monte Carlo simulations here as they predict higher irradiance levels than the KM model (see Supplementary Figures 7, 8). Two different numerical apertures (NA = 0.58 and NA = 0.37) and two probe diameters (200 and 500 μm) were tested to assess the impact of these parameters on the resulting light distribution. As NA had only little influence, it was kept fixed at 0.37 when evaluating the resulting heating.

We tested continuous and pulsed (3 ms pulse duration and 20 Hz repetition rate) stimulation with a light source irradiance of 1 W/mm² that illuminated a cubic tissue volume of 1 mm³. The initial temperature increase saturates after ∼10 s (Figure 7A), where it reaches 2.5–3.5°C for continuous stimulation (Figures 7A–C). For the pulsed stimulation, the effective source light irradiance is decreased to ∼64 mW/mm² and therefore also the heating effect is reduced by an order of magnitude. Maximal temperature changes for different stimulation frequencies were also tested (Figure 7D). For an irradiance level of 4 W/mm² (corresponding to a light output power of ∼120 mW), a temperature rise of 2°C is expected for 40 Hz stimulation. Considering that these choices of the stimulation frequency and maximal output power represent a worst-case scenario (see Figures 4, 5 and the corresponding results sections), we can safely conclude that the temperature rise will stay below 2°C. Very high light intensities, such as the ones used in this model, may cause other types of photochemical damage, however, such effects were beyond the scope of the present study.

FIGURE 7

Figure 7. (A) Maximal temperature changes in the tissue (z = 200 μm) in response to continuous and pulsed (3 ms duration and 20 Hz) stimulation with 30 mW output power for probes with 200 and 500 μm diameter (NA = 0.37). (B) Spatial profile of the temperature increase (t = 10 s) along the axial (z) direction. (C) Spatial profile of the temperature increases along radial direction for z = 200 μm. (D) Maximal temperature changes for different stimulation frequencies (source irradiance 4 W/mm², pulse duration 3 ms). An irradiance level of 4 W/mm² at the light probe corresponds to a light power of ∼120 mW from an optical fiber with a diameter of 200 μm. The temperature increases proportionally with the frequency.

Discussion

The main aim of this study was to evaluate the spatiotemporal characteristics of the neural magnetic fields of the auditory VCN-MNTB pathway in transverse slices of the mouse brainstem. For this purpose, we constructed a network model with realistic dynamics and morphology and simulated the magnetic fields in response to electrical and optical stimulations. Our results suggest that the peak magnetic field strengths of ∼0.5 and ∼0.1 nT should be expected from a 300 μm thick brainstem slice, after electrical and optical stimulations, respectively, when assuming that ∼1,250 MNTB cells are contained in the slice. That is, while optical stimulation would be beneficial to avoid artefacts in the signal that occur for electric stimulation, it also results in substantially weaker neural magnetic fields. Increasing the neural response to optical stimulation would likely require a higher light power than achievable with currently available commercial systems (>100 mW) to activate a sufficient amount of the cells at the same time, then reaching peak magnetic fields of ∼0.05 nT for 300 MNTB cells (Supplementary Figure 11), or ∼0.2 nT for the more optimistic assumption of 1,250 cells. While a light power of 100 mW can currently not be realized, it is worth noting that the temperature changes would still be limited to a physiologically acceptable range given that a pulsed stimulation mode is employed.

Characterizing the neural magnetic response to optical stimulation was one of the main aims of this simulation study. In the following, we will therefore discuss methodological factors that may affect the reliability of simulating optical stimulation. In addition, we give an overview of putative improvements of the efficacy of optogenetic stimulation. We will then relate the simulated neural response to the sensitivity levels that have been achieved for NV-based magnetic field recordings to give an estimate of the practical feasibility of the suggested recording approach. We conclude with a brief summary and outlook.

Simulating and Improving the Neural Magnetic Response to Optical Stimulation

Our estimations of the extracellular neural fields that can be triggered by optical stimulation were based on methods that rely in part on simplifying assumptions and uncertain parameters. We tested three different methods to estimate the light distribution in the brain slice. The Kubelka-Munk (KM) model is widely used in optogenetics, even though it relies on the simplifying assumption of isotropic scattering. In order to test for the impact of tissue anisotropy on the results, we therefore also performed simulations with a MC method (Liu et al., 2015). Finally, we tested a beam spread function (BSF) method that was proposed in a recent study to overcome the limitations of the KM model and especially to estimate the optical parameters from the measured attenuation curves (Yona et al., 2016). Generally, we found a good agreement between the results of the KM model and the BSF method, which both predict similar numbers of neuronal activations for the simulated stimulation protocols (Supplementary Figure 9). Their results were also in agreement with published experimental data on the light distribution and attenuation in brain gray matter (Aravanis et al., 2007; Gradinaru et al., 2009; Yizhar et al., 2011; Al-Juboori et al., 2013; Supplementary Figures 7, 8). We therefore utilized the KM model to study the cell responses to optical stimulation. However, our MC simulations gave deviating results, as apparent from a clear underestimation of the attenuation by the brain tissue (Supplementary Figure 7). The reasons for this discrepancy is likely that the MC does account for attenuations in the ballistic regime of photon propagation, which was our region of interest. Still, we decided to use the MC simulations to obtain robust worst-case estimates of the tissue heating.

Another important aspect of the modeling of optical stimulation was the selection of the density of the photo-sensitive channels in the neural membrane, as this has a direct influence on the irradiance threshold and illumination volume of the tissue (Foutz et al., 2012). Its absolute value is unknown and, in line with common practice (Grossman et al., 2011a), we assumed a constant ChR2 density of 1.3e10 channels/cm², based on results for the bacteriorhodopsin expression in an oocyte (Nagel et al., 1995). However, a higher density value was estimated (4.41e12 channels/cm²) from the experimental data of the mammalian nervous system (Arlow et al., 2013). When using this very high channel density, the strength of the equivalent current dipole of the optically triggered neural activity increases by a factor of three (Supplementary Figure 10), but only for the case of very high source irradiance (3.2 > W/mm², corresponding to an optical power of ∼100 mW). Therefore, we think that our finding of strongly reduced peak fields for optical versus electric stimulation holds robustly.

In a related manner, it was interesting to see that merely increasing the light power was not sufficient to achieve highly synchronous activations of all neurons in the brain slice, and that the maximally achievable repetition rate was lower for optical versus electric stimulation. These limitations were also caused by the comparatively slow dynamics and the low conductance of the light-sensitive ChR2 channels, which lowered the strength of the optically controlled membrane currents. Increasing channel density by increasing ChR2 expression might help to ameliorate these problems, but this approach can also change the native membrane environment (Lin, 2010). An alternative and favored approach is increasing the conductance and changing the channel kinetics by developing a ChR2 mutant to achieve (i) a lower threshold to light stimulation, (ii) a better efficiency at high frequencies, and (iii) red shifting the spectral response to reach deeper cells with minimal light scattering and absorption (Lin, 2010). With this aim, the E123T (ChETA) (Gunaydin et al., 2010) and ChD mutations (Lin et al., 2009) were developed for faster opening/closing rates. Channels bearing the abovementioned mutations achieve improved spiking efficiency at high frequencies (up to 200 Hz) at the cost of higher light thresholds (Grossman et al., 2011a). Another mutant (ChIEF) was developed to provide better light sensitivity with more stable photocurrent responses during continuous light or high frequency pulsed light (Lin et al., 2009; Lin, 2010). Therefore, choosing ChETA or ChIEF instead of ChR2 might help to increase synchrony of spiking or to reach higher spiking frequencies. However, we expect the achievable increases of the strength of the extracellular magnetic fields to be only modest, because the generation of synchronous APs throughout the slice would still require a light power which is currently not feasible in practice.

Recently, the expected temperature rise and its effect on the tissues during optical stimulation has been studied in rodent brain in vivo (Stujenske et al., 2015; Senova et al., 2017). These studies conducted optical stimulation of cortical regions in rats and found temperature rises in the range of 0.1–2.5°C, depending on the fiber and stimulation parameters (Senova et al., 2017). The authors did not report any change in neural activity due to heating for stimulation with up to 200 mW/mm² continuous power. Considering low frequency (∼30 Hz) pulsed stimulation, our simulation indicate temperature increases in a similar range with a light intensity of up to 4 W/mm² at the fiber tip (Figure 7D). This suggests that increasing light power to get higher activation should not be considered as a main concern for pulsed optical stimulation.

Relating the Neural Magnetic Fields to the Sensitivity of NV-Based Magnetic Field Recordings

In this section we discuss how the simulated magnetic fields can be detected by the use of current setups of NV diamond sensors. The aim was to predict the signal from an experimental setup similar to the one described in Webb et al. (2021). However, the modeled magnetic fields are independent on the sensor setup. In an experimental NV sensor setup, as seen in e.g., (Webb et al., 2021), the diamond sensor is illuminated from below with a laser with a wavelength of 532 nm (green), and the resulting fluoresce light at ∼650 nm (red) is measured which has an intensity that is directly proportional to the magnetic field. The light is collected through a 600 nm long pass filter to remove any potential interference from the blue light used for ontogenetic stimulation. In principal, the blue light could excite the NV centers. However, they have very low absorption at this wavelength, and the effect would be negligible. Furthermore, in Webb et al. (2021) the light was shielded from the brain slice with a layer of reflective aluminum, so that the interferences between the two different optical setups (i.e., optogenetic stimulation from above, and NV diamond sensor excitation from below) become negligible.

The sensitivity of diamond sensors are limited by the density of the NV ensemble in the diamond, the dephasing time $T_2^{*}$ and the on-resonance fluorescence contrast (Balasubramanian et al., 2008; Taylor et al., 2008). These parameters are determined by factors such as the isotopic composition of the diamond (Balasubramanian et al., 2009), material strain (Kehayias et al., 2019) and the growth method used (Acosta et al., 2009). For neural signals in the DC to low kHz range (millisecond timescale), continuous wave protocols are so far employed for sensing which have achieved volume-normalized sensitivities of up to ${\mathrm{\eta }}_{V}34nT{m}^{\frac{3}{2}}H{z}^{-\frac{1}{2}}$ using an isotopically engineered diamond with n_NV∼1ppm and ${\text{T}}_2^{*}\sim 0.5$ μs for sensing from a biological sample (Barry et al., 2016). When assuming a sampling rate of 10 kHz to reliably reconstruct the shape of the APs and a 5 μm-thick layer of NV centers in the diamond sample, this volume-normalized sensitivity corresponds to an area-normalized sensitivity of around 1,520 nT μm (applying Karadas et al., 2018, Eq. 9). Reasonable reconstructions of the 2D spatial distribution of neural magnetic fields in the range of 0.5 nT require area-normalized sensitivities of 10 nT μm or better (Karadas et al., 2018), so that several thousands of trials would need to be averaged to achieve sufficient reconstruction quality when using existing technology. However, sensitivity improvements of more than two orders of magnitude have been discussed for the next generation NV sensors (Barry et al., 2020), which would render the NV-based imaging of neural magnetic fields feasible.

Although sensitivities in the sub-pT/Hz^1/2 range have recently been reported for NV-based magnetic field sensing (Fescenko et al., 2020; Zhang et al., 2021), measurements of biological samples so far only reached picotesla ranges. This includes 15 pT/Hz^1/2 sensitivity for measurements on an invertebrate sample (Barry et al., 2016) and 50 pT/Hz^1/2 for mouse tissue (Webb et al., 2021), with a noise level of 12 pT being reached after 8 h of optical stimulation at 0.5 Hz. These sensitivities are already sufficient to measure the expected signals in the range of 0.5 nT and lower. They were reached when sensing the average magnetic field in the complete measurement area of the diamond (mm-scale) – essentially measuring the average of the magnetic fields as seen in Figure 6 – and thus sacrificing information about the spatial distribution of the fields. This is a promising intermediate step toward further maturing this novel technology. Obviously, the amplitude of the averaged signal depends on the sensor position relative to the sample. Even for close-by positions, the signal may be close to zero when the signal from areas with field in- and decreases are spatially averaged. This can be most clearly be seen for the Bz-field component in Figure 6 for the present model. In this respect, simulations of the expected spatial field distribution can be helpful to avoid such cancellation effects. Another possible method of getting rid of these cancellation effects and increasing the spatial resolution could be achieved by decreasing the sensor area, for example by using an NV diamond scanning probe, as done by Thiel et al. (2019).

Conclusion

In conclusion, our results serve as a detailed characterization of the electric and magnetic neural signals emitted from the auditory brainstem pathway of a mouse brain in response to optic or electric stimulation. Our results suggest that the APs traveling along the axonal pathway from GB cells to MNTB principal cells will dominate the imaged neural magnetic fields, while magnetic signals from the calyces of Held or the MNTB principal cells will be far weaker. The aim was to evaluate the feasibility of measuring the magnetic field from this neuronal pathway using existing NV diamond sensor technology, and the results showed that this pathway did not elicit sufficiently strong signals for spatially resolving the signals in these structures with either stimulation method. In contrast, the magnetic fields of pyramidal cells in hippocampal brain slices yield much stronger magnetic fields (Karadas et al., 2018) as they are stronger and have slower dynamics (thus allowing for lower sampling rates), which both will increase the SNR of the measurements. While avoiding artefacts in the recorded signals, the optical stimulation method consistently produced much lower peak magnetic fields (0.1 nT) than the electrical stimulation method (0.5 nT), even when the optical stimulation light irradiance at the source was increased to 4 W/mm², which exceeds what is practically reasonable.

Although NV sensor technology is not currently able to image neural magnetic signals to allow distinguishing between different neural sources along the investigated pathway, as is currently done using e.g., fluorescent imaging, the present results suggest that detecting the neural magnetic fields by spatial averaging the field over a whole mm-sized diamond (as done by Barry et al., 2016; Webb et al., 2021) would likely be possible. In order to reach maximal field strengths, the presented model and stimulation results are highly useful for determining good sample positioning for such methods. The simulation approach of neural magnetic fields in response to optical stimulation, as presented here, is flexible and can be easily extended to other brain regions that differ in light scattering and absorption properties. We also expect our approach to be useful to simulate the excitation of various neuronal types under different configurations of light excitation, and to optimize light distribution and illumination conditions.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

MK contributed to the study design, performed the simulations, analyzed and interpreted the data, and drafted a first version of the manuscript. CO contributed to the analysis and interpretation of the data, wrote sections of the manuscript, and led the manuscript revisions. AT contributed to the study design, the interpretation of the data, and the drafting of the first version of the manuscript. NW, J-FP, AH, and UA contributed to the study design. JW wrote a section of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This work was supported by the “EXMAD” project funded by the Innovation Fund Denmark and by an Interdisciplinary Synergy Grant “bioQ” (grant no. 28086) by the Novo Nordisk Foundation.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2021.643614/full#supplementary-material

References

Acosta, V. M., Bauch, E., Ledbetter, M. P., Santori, C., Fu, K.-M. C., Barclay, P. E., et al. (2009). Diamonds with a high density of nitrogen-vacancy centers for magnetometry applications. Phys. Rev. B 80:115202. doi: 10.1103/PhysRevB.80.115202

CrossRef Full Text | Google Scholar

Acosta, V. M., Bauch, E., Ledbetter, M. P., Waxman, A., Bouchard, L.-S., and Budker, D. (2010). Temperature dependence of the nitrogen-vacancy magnetic resonance in diamond. Phys. Rev. Lett. 104:70801. doi: 10.1103/PhysRevLett.104.070801

PubMed Abstract | CrossRef Full Text | Google Scholar

Al-Juboori, S. I., Dondzillo, A., Stubblefield, E. A., Felsen, G., Lei, T. C., and Klug, A. (2013). Light scattering properties vary across different regions of the adult mouse brain. PLoS One 8:e67626. doi: 10.1371/journal.pone.0067626

PubMed Abstract | CrossRef Full Text | Google Scholar

Aravanis, A. M., Wang, L., Zhang, F., Meltzer, L. A., Mogri, M. Z., Schneider, M. B., et al. (2007). An optical neural interface: in vivo control of rodent motor cortex with integrated fiberoptic and optogenetic technology. J. Neural Eng. 4, S143–S156. doi: 10.1088/1741-2560/4/3/S02

CrossRef Full Text | Google Scholar

Arlow, R. L., Foutz, T. J., and McIntyre, C. C. (2013). Theoretical principles underlying optical stimulation of myelinated axons expressing channelrhodopsin-2. Neuroscience 248, 541–551. doi: 10.1016/j.neuroscience.2013.06.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Balasubramanian, G., Chan, I. Y., Kolesov, R., Al-Hmoud, M., Tisler, J., Shin, C., et al. (2008). Nanoscale imaging magnetometry with diamond spins under ambient conditions. Nature 455, 648–651. doi: 10.1038/nature07278

PubMed Abstract | CrossRef Full Text | Google Scholar

Balasubramanian, G., Neumann, P., Twitchen, D., Markham, M., Kolesov, R., Mizuochi, N., et al. (2009). Ultralong spin coherence time in isotopically engineered diamond. Nat. Mater. 8, 383–387. doi: 10.1038/nmat2420

PubMed Abstract | CrossRef Full Text | Google Scholar

Barnes-Davies, M., and Forsythe, I. D. (1995). Pre- and postsynaptic glutamate receptors at a giant excitatory synapse in rat auditory brainstem slices. J. Physiol. 488, 387–406. doi: 10.1113/jphysiol.1995.sp020974

PubMed Abstract | CrossRef Full Text | Google Scholar

Barry, J. F., Schloss, J. M., Bauch, E., Turner, M. J., Hart, C. A., Pham, L. M., et al. (2020). Sensitivity optimization for NV-diamond magnetometry. Rev. Mod. Phys. 92:15004. doi: 10.1103/RevModPhys.92.015004

CrossRef Full Text | Google Scholar

Barry, J. F., Turner, M. J., Schloss, J. M., Glenn, D. R., Song, Y., Lukin, M. D., et al. (2016). Optical magnetic detection of single-neuron action potentials using quantum defects in diamond. Proc. Natl. Acad. Sci. U.S.A. 113:14133. doi: 10.1073/pnas.1601513113

PubMed Abstract | CrossRef Full Text | Google Scholar

Borst, J. G., Helmchen, F., and Sakmann, B. (1995). Pre- and postsynaptic whole-cell recordings in the medial nucleus of the trapezoid body of the rat. J. Physiol. 489, 825–840. doi: 10.1113/jphysiol.1995.sp021095

PubMed Abstract | CrossRef Full Text | Google Scholar

Borst, J. G. G., and Sakmann, B. (1998). Calcium current during a single action potential in a large presynaptic terminal of the rat brainstem. J. Physiol. 506, 143–157. doi: 10.1111/j.1469-7793.1998.143bx.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Borst, J. G. G., and Soria van Hoeve, J. (2012). The calyx of held synapse: from model synapse to auditory relay. Annu. Rev. Physiol. 74, 199–224. doi: 10.1146/annurev-physiol-020911-153236

PubMed Abstract | CrossRef Full Text | Google Scholar

Brew, H., and Forsythe, I. (1995). Two voltage-dependent K+ conductances with complementary functions in postsynaptic integration at a central auditory synapse. J. Neurosci. 15, 8011–8022. doi: 10.1523/JNEUROSCI.15-12-08011.1995

PubMed Abstract | CrossRef Full Text | Google Scholar

Deisseroth, K. (2011). Optogenetics. Nat. Methods 8, 26–29. doi: 10.1038/nmeth.f.324

PubMed Abstract | CrossRef Full Text | Google Scholar

Dolde, F., Fedder, H., Doherty, M. W., Nöbauer, T., Rempp, F., Balasubramanian, G., et al. (2011). Electric-field sensing using single diamond spins. Nat. Phys. 7, 459–463. doi: 10.1038/NPHYS1969

CrossRef Full Text | Google Scholar

Fedchyshyn, M. J., and Wang, L.-Y. (2007). Activity-dependent changes in temporal components of neurotransmission at the juvenile mouse calyx of Held synapse. J. Physiol. 581, 581–602. doi: 10.1113/jphysiol.2007.129833

PubMed Abstract | CrossRef Full Text | Google Scholar

Felmy, F., and Schneggenburger, R. (2004). Developmental expression of the Ca2+-binding proteins calretinin and parvalbumin at the calyx of Held of rats and mice. Eur. J. Neurosci. 20, 1473–1482. doi: 10.1111/j.1460-9568.2004.03604.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Fescenko, I., Jarmola, A., Savukov, I., Kehayias, P., Smits, J., Damron, J., et al. (2020). Diamond magnetometer enhanced by ferrite flux concentrators. Phys. Rev. Res. 2:23394. doi: 10.1103/PhysRevResearch.2.023394

PubMed Abstract | CrossRef Full Text | Google Scholar

Ford, M. C., Alexandrova, O., Cossell, L., Stange-Marten, A., Sinclair, J., Kopp-Scheinpflug, C., et al. (2015). Tuning of Ranvier node and internode properties in myelinated axons to adjust action potential timing. Nat. Commun. 6:8073. doi: 10.1038/ncomms9073

PubMed Abstract | CrossRef Full Text | Google Scholar

Forsythe, I. D. (1994). Direct patch recording from identified presynaptic terminals mediating glutamatergic EPSCs in the rat CNS, in vitro. J. Physiol. 479, 381–387. doi: 10.1113/jphysiol.1994.sp020303

PubMed Abstract | CrossRef Full Text | Google Scholar

Foutz, T. J., Arlow, R. L., and McIntyre, C. C. (2012). Theoretical principles underlying optical stimulation of a channelrhodopsin-2 positive pyramidal neuron. J. Neurophysiol. 107, 3235–3245. doi: 10.1152/jn.00501.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Gradinaru, V., Mogri, M., Thompson, K. R., Henderson, J. M., and Deisseroth, K. (2009). Optical deconstruction of parkinsonian neural circuitry. Science 324, 354–359. doi: 10.1126/science.1167093

PubMed Abstract | CrossRef Full Text | Google Scholar

Graham, B. P., Wong, A. Y. C., and Forsythe, I. D. (2001). A computational model of synaptic transmission at the calyx of Held. Neurocomputing 38–40, 37–42. doi: 10.1016/S0925-2312(01)00476-3

CrossRef Full Text | Google Scholar

Grossman, N., Nikolic, K., Grubb, M. S., Burrone, J., Toumazou, C., and Degenaar, P. (2011a). High-frequency limit of neural stimulation with ChR2. Proc. Annu. Int. Conf. IEEE Eng. Med. Biol. Soc. EMBS 4167–4170. doi: 10.1109/IEMBS.2011.6091034

PubMed Abstract | CrossRef Full Text | Google Scholar

Grossman, N., Nikolic, K., Toumazou, C., and Degenaar, P. (2011b). Modeling study of the light stimulation of a neuron cell with channelrhodopsin-2 mutants. IEEE Trans. Biomed. Eng. 58, 1742–1751. doi: 10.1109/TBME.2011.2114883

PubMed Abstract | CrossRef Full Text | Google Scholar

Gunaydin, L. A., Yizhar, O., Berndt, A., Sohal, V. S., Deisseroth, K., and Hegemann, P. (2010). Ultrafast optogenetic control. Nat. Neurosci. 13, 387–392. doi: 10.1038/nn.2495

PubMed Abstract | CrossRef Full Text | Google Scholar

Hall, L. T., Beart, G. C. G., Thomas, E. A., Simpson, D. A., McGuinness, L. P., Cole, J. H., et al. (2012). High spatial and temporal resolution wide-field imaging of neuron activity using quantum NV-diamond. Sci. Rep. 2:401. doi: 10.1038/srep00401

PubMed Abstract | CrossRef Full Text | Google Scholar

Hegemann, P., Ehlenbeck, S., and Gradmann, D. (2005). Multiple Photocycles of Channelrhodopsin. Biophys. J. 89, 3911–3918. doi: 10.1529/biophysj.105.069716

PubMed Abstract | CrossRef Full Text | Google Scholar

Hines, M. L., and Carnevale, N. T. (2001). NEURON: a tool for neuroscientists. Neuroscientist 7, 123–135. doi: 10.1177/107385840100700207

PubMed Abstract | CrossRef Full Text | Google Scholar

Homma, R., Baker, B. J., Jin, L., Garaschuk, O., Konnerth, A., Cohen, L. B., et al. (2009). “Wide-Field and two-photon imaging of brain activity with voltage and calcium-sensitive dyes,” in BT - Dynamic Brain Imaging: Multi-Modal Methods and In Vivo Applications, ed. F. Hyder (Totowa, NJ: Humana Press), 43–79. doi: 10.1007/978-1-59745-543-5_3

CrossRef Full Text | Google Scholar

Jacques, S. L. (2010). “Monte carlo modeling of light transport in tissue (steady state and time of flight),” in Optical-Thermal Response of Laser-Irradiated Tissue, eds A. Welch and M. van Gemert (Dordrecht: Springer). doi: 10.1007/978-90-481-8831-4_5

CrossRef Full Text | Google Scholar

Joris, P. X., and Trussell, L. O. (2018). The calyx of held: a hypothesis on the need for reliable timing in an intensity-difference encoder. Neuron 100, 534–549. doi: 10.1016/j.neuron.2018.10.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Karadas, M., Wojciechowski, A. M., Huck, A., Dalby, N. O., Andersen, U. L., and Thielscher, A. (2018). Feasibility and resolution limits of opto-magnetic imaging of neural network activity in brain slices using color centers in diamond. Sci. Rep. 8:4503. doi: 10.1038/s41598-018-22793-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Kehayias, P., Turner, M. J., Trubko, R., Schloss, J. M., Hart, C. A., Wesson, M., et al. (2019). Imaging crystal stress in diamond using ensembles of nitrogen-vacancy centers. Phys. Rev. B 100:174103. doi: 10.1103/PhysRevB.100.174103

CrossRef Full Text | Google Scholar

Kronander, E., Michalski, N., Lebrand, C., Hornung, J. P., and Schneggenburger, R. (2017). An organotypic slice culture to Study the formation of calyx of Held synapses in-vitro. PLoS One 12:e0175964. doi: 10.1371/journal.pone.0175964

PubMed Abstract | CrossRef Full Text | Google Scholar

Le Sage, D., Arai, K., Glenn, D. R., DeVience, S. J., Pham, L. M., Rahn-Lee, L., et al. (2013). Optical magnetic imaging of living cells. Nature 496, 486–489. doi: 10.1038/nature12072

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, J. Y. (2010). A user’s guide to channelrhodopsin variants: features, limitations and future developments. Exp. Physiol. 96, 19–25. doi: 10.1113/expphysiol.2009.051961

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, J. Y., Lin, M. Z., Steinbach, P., and Tsien, R. Y. (2009). Characterization of engineered channelrhodopsin variants with improved properties and kinetics. Biophys. J. 96, 1803–1814. doi: 10.1016/j.bpj.2008.11.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Jacques, S. L., Azimipour, M., Rogers, J. D., Pashaie, R., and Eliceiri, K. W. (2015). OptogenSIM: a 3D Monte Carlo simulation platform for light delivery design in optogenetics. Biomed. Opt. Express 6:4859. doi: 10.1364/BOE.6.004859

PubMed Abstract | CrossRef Full Text | Google Scholar

Lorteije, J. A. M., Rusu, S. I., Kushmerick, C., and Borst, J. G. G. (2009). Reliability and precision of the mouse calyx of held synapse. J. Neurosci. 29, 13770–13784. doi: 10.1523/JNEUROSCI.3285-09.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Macica, C. M., von Hehn, C. A. A., Wang, L.-Y., Ho, C.-S., Yokoyama, S., Joho, R. H., et al. (2003). Modulation of the kv3.1b potassium channel isoform adjusts the fidelity of the firing pattern of auditory neurons. J. Neurosci. 23, 1133–1141.

Google Scholar

Manis, P. B., and Campagnola, L. (2018). A biophysical modelling platform of the cochlear nucleus and other auditory circuits: from channels to networks. Hear. Res. 360, 76–91. doi: 10.1016/j.heares.2017.12.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Mobley, J., Vo-Dinh, T., and Tuchin, V. V. (2014). “Optical properties of tissue,” in Biomedical Photonics Handbook, Set, Vol. 3, ed. T. Vo-Dinh (Boca. Raton, FL: CRC Press), doi: 10.1201/b17290

CrossRef Full Text | Google Scholar

Multi Channels Systems MCS GmbH (2018). Microelectrode Array (MEA) Manual. Reutlingen: Multi Channels Systems MCS GmbH.

Google Scholar

Nagel, G., Möckel, B., Büldt, G., and Bamberg, E. (1995). Functional expression of bacteriorhodopsin in oocytes allows direct measurement of voltage dependence of light induced H + pumping. FEBS Lett. 377, 263–266. doi: 10.1016/0014-5793(95)01356-3

CrossRef Full Text | Google Scholar

Nikolic, K., Grossman, N., Grubb, M. S., Burrone, J., Toumazou, C., and Degenaar, P. (2009). Photocycles of channelrhodopsin-2. Photochem. Photobiol. 85, 400–411. doi: 10.1111/j.1751-1097.2008.00460.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ntziachristos, V. (2010). Going deeper than microscopy: the optical imaging frontier in biology. Nat. Methods 7, 603–614.

Google Scholar

Obien, M. E. J., Deligkaris, K., Bullmann, T., Bakkum, D. J., and Frey, U. (2015). Revealing neuronal function through microelectrode array recordings. Front. Neurosci. 8:423. doi: 10.3389/fnins.2014.00423

PubMed Abstract | CrossRef Full Text | Google Scholar

Pennes, H. H. (1948). Analysis of tissue and arterial blood temperatures in the resting human forearm. J. Appl. Physiol. 1, 93–122. doi: 10.1152/jappl.1948.1.2.93

PubMed Abstract | CrossRef Full Text | Google Scholar

Prizmatix (2019). Dual-Optogenetics-LED Fiber Coupled Dual LED for Optogenetics in Freely Moving Animals. Givat-Shmuel: Prizmatix

Google Scholar

Raman, I. M., and Trussell, L. O. (1992). The kinetics of the response to glutamate and kainate in neurons of the avian cochlear nucleus. Neuron 9, 173–186. doi: 10.1016/0896-6273(92)90232-3

CrossRef Full Text | Google Scholar

Richardson, A. G., McIntyre, C. C., and Grill, W. M. (2000). Modelling the effects of electric fields on nerve fibres: influence of the myelin sheath. Med. Biol. Eng. Comput. 38, 438–446. doi: 10.1007/BF02345014

PubMed Abstract | CrossRef Full Text | Google Scholar

Rothman, J. S., and Manis, P. B. (2003). The roles potassium currents play in regulating the electrical activity of ventral cochlear nucleus neurons. J. Neurophysiol. 89, 3097–3113. doi: 10.1152/jn.00127.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Schirhagl, R., Chang, K., Loretz, M., and Degen, C. L. (2014). Nitrogen-Vacancy centers in diamond: nanoscale sensors for physics and biology. Annu. Rev. Phys. Chem. 65, 83–105. doi: 10.1146/annurev-physchem-040513-103659

PubMed Abstract | CrossRef Full Text | Google Scholar

Senova, S., Scisniak, I., Chiang, C., Doignon, I., Palfi, S., Chaillet, A., et al. (2017). Experimental assessment of the safety and potential efficacy of high irradiance photostimulation of brain tissues. Sci. Rep. 7:43997. doi: 10.1038/srep43997

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, P. H., Joris, P. X., and Yin, T. C. (1998). Anatomy and physiology of principal cells of the medial nucleus of the trapezoid body (MNTB) of the cat. J. Neurophysiol. 79, 3127–3142. doi: 10.1152/jn.1998.79.6.3127

PubMed Abstract | CrossRef Full Text | Google Scholar

Stepan, J., Dine, J., and Eder, M. (2015). Functional optical probing of the hippocampal trisynaptic circuit in vitro: network dynamics, filter properties, and polysynaptic induction of CA1 LTP. Front. Neurosci. 9:160. doi: 10.3389/fnins.2015.00160

PubMed Abstract | CrossRef Full Text | Google Scholar

Stujenske, J. M., Spellman, T., and Gordon, J. A. (2015). Modeling the spatiotemporal dynamics of light and heat propagation for invivo optogenetics. Cell Rep. 12, 525–534. doi: 10.1016/j.celrep.2015.06.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Swinney, K. R., and Wikswo, J. P. (1980). A calculation of the magnetic field of a nerve action potential. Biophys. J. 32, 719–731. doi: 10.1016/S0006-3495(80)85012-0

CrossRef Full Text | Google Scholar

Taschenberger, H., and von Gersdorff, H. (2000). Fine-Tuning an auditory synapse for speed and fidelity: developmental changes in presynaptic waveform, epsc kinetics, and synaptic plasticity. J. Neurosci. 20, 9162–9173.

Google Scholar

Taylor, J. M., Cappellaro, P., Childress, L., Jiang, L., Budker, D., Hemmer, P. R., et al. (2008). High-sensitivity diamond magnetometer with nanoscale resolution. Nat. Phys. 4, 810–816. doi: 10.1038/nphys1075

CrossRef Full Text | Google Scholar

Tetienne, J.-P., Dontschuk, N., Broadway, D. A., Stacey, A., Simpson, D. A., and Hollenberg, L. C. L. (2017). Quantum imaging of current flow in graphene. Sci. Adv. 3:e1602429. doi: 10.1126/sciadv.1602429

PubMed Abstract | CrossRef Full Text | Google Scholar

Thiel, L., Wang, Z., Tschudin, M. A., Rohner, D., Gutiérrez-Lezama, I., Ubrig, N., et al. (2019). Probing magnetism in 2D materials at the nanoscale with single-spin microscopy. Science 364:973. doi: 10.1126/science.aav6926

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Jacques, S. L., and Zheng, L. (1995). MCML-Monte Carlo modeling of light transport in multi-layered tissues. Comput. Methods Programs Biomed. 47, 131–146. doi: 10.1016/0169-2607(95)01640-F

CrossRef Full Text | Google Scholar

Wang, L.-Y., Gan, L., Forsythe, I. D., and Kaczmarek, L. K. (1998). Contribution of the Kv3.1 potassium channel to high-frequency firing in mouse auditory neurones. J. Physiol. 509, 183–194. doi: 10.1111/j.1469-7793.1998.183bo.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Weatherstone, J. H., Kopp-Scheinpflug, C., Pilati, N., Wang, Y., Forsythe, I. D., Rubel, E. W., et al. (2017). Maintenance of neuronal size gradient in MNTB requires sound-evoked activity. J. Neurophysiol. 117, 756–766. doi: 10.1152/jn.00528.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Webb, J. L., Troise, L., Hansen, N. W., Olsson, C., Wojciechowski, A. M., Achard, J., et al. (2021). Detection of biological signals from a live mammalian muscle using an early stage diamond quantum sensor. Sci. Rep. 11:2412. doi: 10.1038/s41598-021-81828-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wojciechowski, A. M., Karadas, M., Huck, A., Osterkamp, C., Jankuhn, S., Meijer, J., et al. (2018). Contributed Review: camera-limits for wide-field magnetic resonance imaging with a nitrogen-vacancy spin sensor. Rev. Sci. Instrum. 89:031501. doi: 10.1063/1.5010282

CrossRef Full Text | Google Scholar

Woosley, J. K., Roth, B. J., and Wikswo, J. P. (1985). The magnetic field of a single axon: a volume conductor model. Math. Biosci. 76, 1–36. doi: 10.1016/0025-5564(85)90044-6

CrossRef Full Text | Google Scholar

Xie, R., and Manis, P. B. (2013). Target-Specific IPSC kinetics promote temporal processing in auditory parallel pathways. J. Neurosci. 33, 1598–1614. doi: 10.1523/JNEUROSCI.2541-12.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Yizhar, O., Fenno, L. E., Davidson, T. J., Mogri, M., and Deisseroth, K. (2011). Optogenetics in neural systems. Neuron 71, 9–34.

Google Scholar

Yona, G., Meitav, N., Kahn, I., and Shoham, S. (2016). Realistic numerical and analytical modeling of light scattering in brain tissue for optogenetic applications. eNeuro 3, 1–9. doi: 10.1523/ENEURO.0059-15.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, C., Shagieva, F., Widmann, M., Kuebler, M., Vorobyov, V., Kapitanova, P., et al. (2021). Diamond magnetometry and gradiometry with subpicotesla DC sensitivity. arXiv [preprint] arXiv:2012.15706.

Google Scholar