Feasibility and resolution limits of opto-magnetic imaging of neural network activity in brain slices using color centers in diamond

We suggest a novel approach for wide-field imaging of the neural network dynamics of brain slices that uses highly sensitivity magnetometry based on nitrogen-vacancy (NV) centers in diamond. In-vitro recordings in brain slices is a proven method for the characterization of electrical neural activity and has strongly contributed to our understanding of the mechanisms that govern neural information processing. However, this traditional approach only acquires signals from a few positions, which severely limits its ability to characterize the dynamics of the underlying neural networks. We suggest to extend its scope using NV magnetometry-based imaging of the neural magnetic fields across the slice. Employing comprehensive computational simulations and theoretical analyses, we determine the spatiotemporal characteristics of the neural fields and the required key performance parameters of an NV magnetometry-based imaging setup. We investigate how the technical parameters determine the achievable spatial resolution for an optimal 2D reconstruction of neural currents from the measured field distributions. Finally, we compare the imaging of neural slice activity with that of a single planar pyramidal cell. Our results suggest that imaging of slice activity will be possible with the upcoming generation of NV magnetic field sensors, while single-shot imaging of planar cell activity remains challenging.

where the operator ∇ ′ acts on ′. Since there are no current sources in the extracellular medium (i.e., ∇ source = 0 in the extracellular space), the first term on the right hand side disappears. The extracellular potential is much smaller than the intracellular potential 2 (i. e. , σ i ϕ i ≫ σ e ϕ e ) and hence Eq. (S1) can be approximated as: Therefore, the extracellular potential is approximately independent of the external current distribution. If the nerve structure has a small cross-section A, we can assume that is only changing along the structure (this direction is parameterized using the longitudinal coordinate l) and the surface integral can be written in terms of a volume integral 2 : The numerical simulations of the neuron are based on core-conductor equations, which assume a constant current density across the cross-section of the cell, 2/23 where r i = 1/σ i A is the intracellular resistance per unit length and i m is the membrane current per unit length. Eq. (S4) allows us to write the extracellular potential in dependence on the membrane currents: (S5) The core-conductor model provides a good approximation as we are only interested in the transmembrane potential of the cell model 3 where I m n is the membrane current for compartment n, Δs n is its length, and a is its radius: where ρ n is the distance perpendicular to the center line of the compartment, h n is the longitudinal distance from the end of the compartment, and l n = h n + Δs n denotes the longitudinal distance from the start of the compartment 4,5 .
Figure S1: Model of a cylindrical neural process. The left side shows the volume-conductor model and the right side the equivalent core-conductor model. The volume conductor model distinguishes between the intracellular, membrane and extracellular regions. In the core-conductor model, the intracellular and extracellular regions are represented by axial resistances and they are connected by a membrane network 6 .
A similar solution can be derived for the extracellular magnetic field (starting from Eqs. 20 and 22 in 1 ) : The domain is divided into the intracellular, membrane and extracellular volumes, and the integral equations run over the complete domain. Since the extracellular and intracellular volumes are source free ohmic regions, the first integral is nonzero only at the membrane. We thus define source = as the biochemical current sources in the membrane. The contribution of the secondary currents in the second integral of Eq.
(S9) can be written for the intracellular and extracellular regions by considering the vector identities: (S10) For piecewise homogeneous regions, the volume integral on the right hand side is zero. Therefore, the second volume integral of Eq. (S9) can be converted into a surface integral and the magnetic field can be expressed as a sum of volume and surface integrals over the membrane 7 : (S11) The current densities are symmetric along the nerve structure and can be written in terms of their radial components J r and longitudinal components J l (r, l) = J r (r, l)� + J l (r, l) ̂ (S12) with � and ̂ being radial and axial unit vectors, respectively. The membrane conductivity is extremely small with respect to the conductivities of the inner and outer regions. We therefore can assume that the membrane currents have only a radial component 7 . As the membrane is very thin, the contribution of the membrane source to the magnetic field can be safely neglected and the volume integral in Eq. (S11) can be skipped.
At the inner and outer surfaces of the membrane, the radial components of the extracellular and intracellular current densities cannot contribute to the surface integrals since they are orthogonal to the surface normal n �.
Therefore, the total magnetic field is determined by the longitudinal (axial) components of currents J e and J i at the outer and inner surfaces, which are defined as: J e,i = σ e,i ∂Φ e,i ∂l (S13) It can be further assumed that the external current flow J e is distributed in a large volume, such that its longitudinal component at the outer membrane surface is small. In contrast, the internal current flow J i is confined to a small volume and its longitudinal component is strong and dominates Eq. (S11). Thus, Eq. (S11) can be approximated as: Discretizing the neuron into N compartments results in: where I axial n is the intracellular axial current for compartment n.
where ρ n the distance perpendicular to the line compartment, h n is the longitudinal distance from the end of the compartment, and l n = h n − Δs n the longitudinal distance from the start of the compartment.

S2 Comparison with the Numerical Solution for a Long Straight Axon
In the following, we validate equation (S17) by comparing the magnetic field obtained for a long straight axon with the fields calculated with the Finite-Element Method (FEM). For completeness, we also present the results obtained with two different analytical approaches, as described in 8 . The total extracellular magnetic field results from two different sources, namely the extracellular and intracellular currents. As discussed above, the magnetic field component due to the membrane currents J m can be safely neglected. As first step, the distributions of the electric potential in the intra-and extracellular spaces are calculated. The potential has to satisfy Laplace's equation in the volume conductor ∇. � , ∇ϕ e,i � = 0, (S18) whereby we set Neumann boundary conditions at the outer and inner boundaries of the cell membrane: . ( ∇ϕ e ) = −J n , and (S19) Vector n denotes an outwards pointing unit vector normal to the boundary. J n represents the current density at the boundaries. For the extracellular potential, we additionally set J n = 0 at the outer domain boundary.
with Δs n being the length and a the radius of compartment n. After determining ϕ e and ϕ i , the current densities are calculated using Ohm's law = , ∇ϕ e,i , and the Biot-Savart law is used to determine the contributions of the external and internal current flows to the total extracellular magnetic field.

6/23
The modelled geometry is depicted in Fig. S2A&B. The membrane potential and membrane currents were simulated using NEURON, and the results fed into the FEM calculations based on COMSOL Multiphysics (v5.3). We made use of the axisymmetric geometry around the Y axis to set up a 2D triangle mesh with a symmetry constraint at z=0 in order to describe the geometry in a computationally efficient way. The geometry corresponds to a cylinder of 40 mm height (in Y direction) and 60 mm diameter (corresponding to 30 mm in Z direction), with the nerve being placed in its center. The triangle density in the region close to the nerve membrane was strongly increased to account for the thin diameter of the axon. Quadratic finite elements were employed, and the relative tolerance for the conjugate gradient solver were required to be The axon has active Hodgkin-Huxley channel dynamics and a length of 3.2 mm. It is discretized in compartments with lengths smaller than the electrotonic space constant. Its diameter was chosen small (5 µm) in order to be in the range seen for hippocampus cells. The following parameters were used: agreement. From this, we can conclude that our forward modeling scheme based on Eq. (S17) is accurate when the assumption of a large homogeneous extracellular volume is valid.
In addition, the results obtained with two different analytical approaches, as described in 8 , are shown. Those approaches allow calculating the magnetic field from the simulated membrane potential, based on the representation of the axon as an infinite cylinder of radius a embedded in an isotropic unbounded outer region. The first analytical method, based on Ampere's law (B Ampere -purple curve in Fig. S2D), strongly overestimates the magnetic field caused by the internal current flow when the distance to the nerve is large, which is in accordance to the known limitations of this approach 8 . The second method that evaluates the Biot-Savart law in the spatial frequency domain (BS kspace -red dashed curve in Fig. S2D) performs better, but systematically underestimates the field strength at all distances. This effects only occurs for thin axon radii as tested here, while it is accurate for thicker axons as simulated in the original study 8 .

S3 Correction for the effects of tissue anisotropy and boundaries
The forward modeling scheme so far assumed that the extracellular volume V e can be well approximated as an infinite homogeneous conductor. However, the brain slice will be embedded in a saline bath, resulting in a conductivity boundary between brain and saline. In addition, an anisotropic electrical conductivity has been reported, e.g. for the guinea-pig hippocampus 9 . Finally, also the boundaries between the slice/saline and the 8/23 surrounding non-conductive air, plastic and diamond affects the extracellular current pathways. The effect of conductivity discontinuities on electric potential recording was explored for Multielectrode array in detail 10 . In the following, we investigate the impact of these conductivity discontinuities on the magnetic field created by a straight nerve fiber. For this, we employ numerical simulations based on the finite element method (FEM) to determine the magnetic field contribution by the extracellular currents (Eqs. S18 to S21), while determining the contribution of the internal currents directly from Eq. (S17). This allows us to derive multiplicative factors that can be applied to the results of the forward modeling scheme to approximately correct for the impact of the finite extracellular volume. We apply these factors to rescale the magnetic fields of the pyramidal neurons investigated in the main part of the study in dependence of their distance to the diamond surface. While the large number and the histological complexity of the simulated pyramidal neurons prevents a direct assessment of their extracellular currents via FEM, our strategy allows us to account for the major effects caused by the extracellular volume. the extracellular current density and the resulting magnetic field. The extracellular region was modelled as tetrahedral mesh, with an increased tetrahedral density close to the nerve fiber. Quadratic finite elements were employed, and the relative tolerance for the conjugate gradient solver were required to be <10 −3 .
Since the nerve is aligned along the Y direction, only the B X and B Z components of the magnetic field are strong at the diamond surface and are analyzed further. Figure S4 shows  Fig. S2C). Figure S5 shows the parts of the B X component caused by the internal and external current flow for the A+S scenario, revealing that the impact of the external currents is increased by around three orders of magnitude compared to an infinite extracellular volume. The nerve is far closer to the diamond surface than to the upper boundary of the saline bath, so that the major contribution comes from the diamond surface, which causes an increased current density in the extracellular space between the nerve cell and diamond. Interestingly, this effect gets stronger with increasing distance between the nerve and the surface, as this increases the overall amount of current in the space between the nerve cell and diamond. Figure S6 shows   We were further interested how the spatial distribution of B X and B Z at the diamond surface is affected by the finite extracellular volume conductor. Figure S7 exemplarily shows the results for a nerve placed at a distance of 150 µm from sensor plane for the S+A scenario. As expected, the spatial distribution of B Z is hardly affected by the external current flow. In contrast, B X gets "narrower" in X direction, while mostly keeping its shape in Y direction. . For B Z , the coefficient of determination between the magnetic field distributions of the internal and total current flow are close to 1 for all simulated distances between nerve and diamond surface (Fig. S8B). The impact of the finite extracellular volume conductor is larger for the spatial distribution of B X , with the coefficient of determination varying between 0.83 and 0.94. As visible in Fig. S7, this is mainly due to a narrowing of B X in the direction normal to the nerve fiber. This effect is neglected in the results shown for the pyramidal cells in the main part of the paper, which were only corrected using the scale factors.
As a consequence, the magnetic fields will be spatially slightly more confined than shown there, and also the achievable spatial resolution will be slightly underestimated. For comparison, we also tested the impact of the finite extracellular volume conductor on the electric potential at the diamond surface (Fig. S9). Interestingly, the electric potential depends more strongly on the volume conductor properties than the magnetic field. The tendency is that stronger electric potentials occur at the surface compared to the infinite homogenous volume conductor (denoted NB in Fig. S9A). However, this effect disappears when a thin saline "leakage" layer between the tissue compartment and the surface is introduced (scenario S+A+L). The relatively well conductive saline layer strongly reduces the potential differences at the surface. A similar result is seen for the spatial distribution of the electric potential at the surface (Fig. S9B), with the NB and S+A+L scenarios giving very similar results. Considering this similarity, we decided against applying correction factors for the electric potential in the main part of the study.
13/23 Figure S9: Effect of the finite extracellular volume conductor on the strength and spatial distribution of the electric potential at the diamond surface. (A) Electric potential at the diamond surface, sampled along a line in Y direction underneath the axon, for three different distances between the axon and the surface. The same four scenarios were tested as for the magnetic field. In addition, the results for an infinite homogeneous extracellular volume conductor are shown for comparison (termed NB -no boundaries). (B) Spatial distribution of the electric potential for a nerve placed at 200 µm distance from the diamond surface. When a thin saline "leakage" layer is introduced between the tissue and the diamond surface (scenario S+A+L), the field is almost identical to that obtained for an infinite homogeneous extracellular volume conductor.

S4 Comparison with Existing Studies: Single Axon case
We modeled a single axon to compare it with previously reported magnetometry measurements of axonal action potentials 12 . In that study, NV magnetometry measurements were performed for excised giant axons of two species, namely Myxicola infundibulum (worm) and Loligo pealeii (squid). In addition, recordings from axonal activity in intact worms were performed. Here, we compare our simulations based on the forward modeling scheme (Eq. S17) with their results for the excised worm axon and the intact worm. In their first experiment, the giant axon of the worm was excised and measured at a temperature of 21°C and a standoff distance of 0.3 mm to the nerve center. In their second experiment, the worm was kept at 10 °C and the measurements were taken at a distance 1.2 mm from the nerve center. To compare our numerical methods with the reported measurement results, we modeled single axons with σ i = 1.5/Ωm and diameters ranging between 200 µm and 400 µm in NEURON and estimated the resulting magnetic field based on the forward modeling scheme outlined above. The assumed temperatures were either 21°C (excised axon) or 10°C (axon of living worm). In the excised axon (Fig. S10A), the simulated peak magnetic field strength increases from around 1 to 3.5 nT with increasing diameter, which overlaps well with the measured field strengths. For the axonal activity of the living worm (Fig. S10B), the peak magnetic field strength decreases to one fourth of the one observed for the excised axon due to the increased distance. In addition, its duration is longer due to the lower ambient temperature that slows down the action potential. Our simulation results are in good agreement with the previously presented measurements 12 , which show that the peak field strength was around 0.4 nT for a diameter of 300 µm. It should be noted that we did not correct for the magnetic field of currents in the extracellular space, as the correction factors derived above are specific for the geometry of the envisioned recordings from the hippocampus slices. of the conductance of the synapses. These conductance changes are described by a double-exponential function g syn = g syn A −1 �e −t/τ rise − e −t/τ decay �, with constant τ rise = 1.5 ms controlling the rise time, constant τ decay = 2.5 ms controlling the decay time and g syn denoting the peak conductance 14 . Constant A 15/23 is used for normalization to set the maximum value of g syn to be g syn . These synaptic events inject currents of I syn = g syn (V m − V rev ) at the membrane positions of the synapses (a reverse potential of V rev = 0 mV is used). The simulations of the neural dynamics are conducted in NEURON, using the backward-Euler integration method to solve the cable equation with a fixed time step of 25 μs. For that, each section of a CA1 pyramidal cell (such as a dendritic branch) is subdivided into multiple cylindrical compartments to achieve compartment lengths smaller than the electrotonic space constant λ. Furthermore, each section is forced to have at least three compartments to be able to calculate axial currents from the membrane voltages as input to the above forward modelling scheme. The resulting models of CA1 pyramidal cells consist of 1427 compartments.
The stimulation of a CA1 pyramidal cell via the Schaffer collaterals is mimicked by synaptic events, which are generated at 40 randomly selected locations at the s. oriens and 40 random positions at the s. radiatum 15 .
The time points of the synaptic events are modelled by truncated random functions with a normal distribution as follows:

Generate an excitation time T exc
2. If T start ≤ T exc ≤ T stop , assign T exc as time point of the synaptic event, otherwise regenerate T exc until it remains in the desired time window.
(0, ) is a normal distribution with zero mean and standard deviation (σ = 0.25, unless indicated differently; corresponding to a standard deviation of 1.56 expressed in ms), and represents the time jitter in the excitatory inputs. The first wave of synaptic events is generated using T start = 0 and T stop = 25 ms, and a second wave is generated using T start = 25 and T stop = 50 ms in order to mimic repeated electric stimulation at 40 Hz. The synaptic strength is controlled by setting the peak conductance g syn . In order to evaluate an upper limit for the magnetic field strengths that can occur merely due to EPSPs without inducing action potentials, a value of 0.3 nS was chosen. Consistent spiking activity could be created by selecting a strength of 0.6 nS (or higher, if indicated). It is worth noting that the number of modelled excitatory synaptic inputs is far lower than occurring for real neurons. In order to reduce the computational complexity of the simulations, the effects of a high number of excitatory synapses is mimicked by increasing the synaptic strengths instead.

S6 Comparison with Existing Studies: Recordings from the hippocampus CA1 subfield with µ-SQUIDs
We modeled a hippocampus slice to compare it with previously reported magnetometry measurements from the hippocampus CA1 subfield of the guinea-pig based on a µ-SQUIDs and small pickup coils 16  However, we do not expect a too strong influence in this case, as the magnetic field component normal to the boundary between conductive tissue and non-conductive material (the diamond and dewar, respectively) is hardly affected (see results for B Z in Fig. S8).
Figure S11: Topography of the simulated neural magnetic fields due to activation of a longitudinal CA1 slice, recorded by µ-SQUIDs. The centers of the pickup coils were placed on a regular grid with a 4 mm spacing (as indicated by the black dots). The brain slice is indicated by the rectangle. Shown is the magnetic field component normal to the slice. The field topography (strong directly neighboring the slice, with opposite signs for left and right positions; weak directly above the slice) shows that the magnetic fields are generated by longitudinal currents flowing along the cell bodies.

S7 Computation of the True Volume Current Density
During the computations of the neural magnetic fields, the segments of neurons are modelled as currentcarrying wires. However, during the reconstruction of the current density from the simulated magnetic field measurements, the currents are assumed to be volume currents. In the following, we therefore derive a mapping between the axial currents and the corresponding volume current densities for comparison with the reconstructed current densities, y. In our modeling, we computed the total current within a voxel by integrating the axial current over the neuron segments in this voxel. The volume current density is then estimated by dividing by the volume of voxel where V i is the volume of voxel and L i n is the xy projection of segment vector n to voxel i.
For the CA1 slice, the volume current density was computed at 64x64 points with a total area of 1 mm 2 and an assumed thickness of 0.3 mm. For the planar cell, the voxel dimension was chosen to be 4x4x2 µm 3 for a total area of 0.09 mm 2 . To test the validity of our estimation, the magnetic fields of both the CA1 slice and the planar cell were determined from the volume current densities using the frequency domain approach described in 8 and compared with the magnetic field determined from the axial currents (Fig. S12). The 18/23 simulation parameters were kept same as for Fig. 2 and Fig. 3 of the main text. The volume currents, projected onto the measurement plane of the sensor, give consistent results compared to the actual axial sources, except small discrepancies that originate from discretization errors. Moreover, in the planar cell model, a thickness of 2 µm is chosen based on the mean radius of segments, resulting in higher fields around the thicker soma layer.   Fig. 7 in the main paper. Here, the neural source has a height of 50 µm (instead of 300 µm) in Z direction, starting at a distance of 50 µm from the diamond surface. Correspondingly, the total current induced by the source in Y direction is 1/6 of the current for the source in Fig. 7. A) Shown is the achievable system resolution (characterized by the full width at half maximum -FWHM -of the reconstructed neural axial current density) in dependence on the in-plane resolution of the sensor. The achieved system resolution is almost constant for higher in-plane resolutions, but drops off quickly for too low in-plane resolutions. (B) Peak signal-to-noise ratio (pSNR) for recordings from the CA1 subarea as a function of the in-plane resolution of the sensor (pixel size) and the signal-to-noise power of the source, η. As long as the in-plane resolution is chosen high enough to keep the best possible system resolution, also the pSNR of the reconstructed axial current densities remains constant. (F) The normalized axial current densities reconstructed from the corresponding magnetic fields in E. Each plot is normalized to the maximum of its noiseless reconstruction with the same spatial filter settings.