Neural anatomy and optical microscopy (NAOMi) simulation for evaluating calcium imaging methods.

New-method: We provide a new framework for evaluating two-photon microscopy methods via in silico Neural Anatomy and Optical Microscopy (NAOMi) simulation. Our computationally efficient model generates large anatomical volumes of mouse cortex, simulates neural activity, and incorporates optical propagation and scanning to create realistic calcium imaging datasets.

Comparison with existing methods: NAOMi is a first-of-its kind framework for assessing optical imaging modalities. Existing methods are either anatomical simulations or do not address functional imaging. Thus there is no competing method for simulating realistic functional optical microscopy data.

4. Materials and Methods Our TPM simulator is designed to permit testing of many different aspects of the calcium imaging process. To achieve this flexibility, our simulator is divided into five distinct modules, each focused on a portion of either the tissue or scanning simulation ( Fig. 1 ). The five modules are: 1) the neuron module responsible for generating single neurons, 2) the volume module responsible for assembling the neurons into a tissue volume that includes neuropil and vasculature, 3) the activity module that generates the temporal calcium traces for each neuron and the neuropil, 4) the optics module that simulates the point-spread function and occlusion due to the optical mask, and 5) the imaging module that simulates the TPM noise model and object motion. 4.1.

Neuron model

The first module creates simple, yet plausible models of neurons that can be placed throughout a volume and scanned in simulation. We model the neural shape via a probability distribution over smooth deformation of a sphere, followed by a nonlinearity. This model allows for fast sampling of unique neurons, meaning that each simulated volume will contain a completely new set of neurons. Additionally, we provide for each neuron a nucleus modeled as a shrunken and smoothed version of the soma shape. This model captures the relationship observed in detailed field emission scanning electron microscopes (FESEM) [ 101 ]. Finally, we simulate for each neuron a number of dendrites, one of which is created thicker and at a downward orientation, as to model the apical dendrites. The model of the smoothly deformed cell body is an isotropic Gaussian process [ 102 ] defined over a sphere. To sample from this distribution and create the cell body, we sample uniformly over a sphere [ 103 ], sample i.i.d. a Gaussian random variable for each point, and smooth the points according to the process covariance. Denoting the sample points p i ∈ ℝ 3 , the height (distance from center of the sphere) can be sampled from r i ~ N 0 , K , K i , j = e − d p i , p j / l , where l is the length-scale that controls the smoothness of the cell body, and d (·, ·) is the geodesic distance between any two points. For the unit sphere (radius one), this distance is the arc length along the great circle connecting the two points d p i , p j = 2 sin − 1 p i − p j 2 2 . When unconstrained, the radial height of this function can, at times, exceed the maximum and minimum realistic deformations r max = max i | r i | and r min = min i | r i |. We thus rescale the radii values as r i = r max − r min r i − min i r i max i r i − min i r i + r min The resulting points r i p i form the points for a mesh that define the interior of the cell body. To account for the non-spherical shape found in pyramidal neurons, we can modify the radii values by making the base radius at each point dependent on a function of its location on the sphere. Specifically, we use the equation for a tear-drop that is defined parametrically by the azimuth and elevation angles ϕ , and θ as (1) cos ϕ sin θ sin m θ 2 , sin ϕ sin m θ 2 , cos θ . The final step in creating the cell body is to create the nucleus, which is accomplished by shrinking and smoothing the cell wall shape as defined by r i p i as r i ∗ = r 50 % r i − min i r i r 5 % + r i + r min − Δ r m i n Dendrites are added to each neuron via a stochastic growing process [ 59 ]. The process generates start and end points for each dendrite, and iteratively grows the dendrite through the volume while avoiding any obstacles (i.e. other cell’s somas, dendrites, or blood vessels). Apical and basal dendrite endpoints are separately set within the volume and the grown dendrites are dilated to widths consistent with measured anatomy ( Tab. 1 , Sup. Fig. 7 ). 4.2.

Volume generation

To create the tissue volume, we initialize an empty volume and begin by placing blood vessels throughout the volume. For computational feasibility, the volume is modeled as a 3-D grid of points with sub-micron sampling (we typically use 0.5 µ m distances). The blood vessels are grown in three parts: surface vasculature, diving arterioles, and capillaries. Surface vasculature is grown by connecting nodes randomly placed upon the surface of the volume. The connected paths are smoothly varied and dilated. Diving arterioles are set at endpoints of surface vasculature and connected to the bottom of the simulated volume. Capillaries are connected to the diving arterioles and pseudo-randomly placed within the volume in a space-filling fashion. Vessel diameters, concentration, branching frequency, and orientation were compared and fit to two-photon microscopy data of mouse vasculature (Schaffer-Nishimura lab, unpublished data). Once initialized, the volume is then filled with the neuron somas. We sequentially place the neurons randomly throughout the empty space in the volume, with a minimum distance that allows cell bodies, but not nuclei, to overlap. The random placement can be modified to encourage neurons to be more spread out, or more cluttered. When an overlap occurs, the overlapping region is given to the latest cell to be placed. This allows our volume to contain touching cell bodies. Once all the cell somas are placed, dendrites are grown for each neuron sequentially, such as to avoid location conflicts with other cells. Apical dendrites are grown in the same fashion, only thicker, axially oriented, and having fewer transversal deviations. Separate apical dendrites corresponding to neurons in deeper cell layers are grown in a similar fashion from the bottom of the volume to the top. As a final step, axons fill up the remaining empty space, up to the typical 0.7 fill fraction of layer 2/3 in mouse V1. The same dendrite growing algorithm [ 59 ] is used to create millions of short axon segments throughout the entire volume. To obtain the global correlated background components, axon segments are locally grouped together. The axon groups are assigned to individual cells by minimizing their centroid distance to cell bodies, and then all remaining axons groups are randomly assigned. 4.3.

Time-trace generation

To simulate temporal activity, we provide a number of options to generate time-traces for each neuron. We provide both statistical models that generate stereotypical activity as well as more detailed [Ca 2+ ] dynamics model. The statistical model provides a simple way to input basic behaviors of various fluorescent proteins (i.e. rise-time and decay). The [Ca 2+ ] dynamics model simulates the molecular kinetics over time, and provides a way to test the time-trace assumptions made in calcium imaging analysis algorithms. 4.3.1.

Spike-time generation

We provide two methods to generate spike trains to drive the fluorescence activity simulation. The first method creates independent activity for each neuron, including bursting behavior. The second model simulates a Hawkes process which accounts both for self-excitation, driving busting behavior, as well as inter-neuron spiking correlations [ 104 ]. To generate independent spike trains, we model each neuron as a bursting neuron, where bursts occur at independent, exponential intervals (2) P Δ t b u r s t = λ b u r s t e − Δ t b u r s t / λ b u r s t for ∆ t burst > 0. The rate of bursting λ burst is chosen differently for each neuron. The rates can be given to the simulator, or the simulator can automatically draw burst rates from a Gamma distribution with a provided mean rate and parameter α = 1. For each burst, the number of spikes are chosen as (3) N b u r s t = 1 + Poisson λ N where the parameter λ N controls the length of the bursting. The inter-spike times between spikes in a burst were modeled as uniformly random between 5 ms and 7 ms. Alternative distributions of spiking activities can easily be implemented by passing different λ burst , α , or λ N values to the simulator, or by direct modification to the code to implement different distributions that better reflect activity in other cortical areas. For the Hawkes model simulation, we first generate a connectivity matrix that encodes how each neuron’s firing excites other neurons. We model this connectivity with a Watts-Strogatz small-world network model [ 64 ]. To correlate the processes to the network activity, we allow for all neurons to influence the background processes, while not allowing many return connections. To stabilize the point-process, we normalize the resulting connectivity matrix to have maximum eigenvalue magnitude of 0.98. We then run the Hawkes process using Lewis’ method [ 105 ], with an exponential distribution over the neuron’s base firing rates and a higher base firing rate for the background components. Finally, we bin the resulting continuous-time spike events into 1 ms bins to create the discretized spikes that are then fed into the calcium dynamics simulation. 4.3.2. AR-p dynamics For each cell, we generate a baseline fluorescence, β i = 1 + z where z ~ N 0 , σ 2 is a Gaussian random variable. The variance σ 2 controls the distribution of baselines, and we set a default value to σ 2 = 0.04. The next step is to simulate the spike or “event” times for each neuron. As most neurons are sparsely active, we draw the firing rate of each neuron as λ i ~ Gamma ( a , θ ). The parameter θ gives the average inter-spike distance in time and should be set according to the temporal sampling rate set in the simulation. The parameter a is the shape parameter and modulates the distribution of the firing rates. We find that a = 1 (where the Gamma distribution collapses to an exponential distribution) yields realistic distributions of neuron activity levels. The actual event times are then sampled for each neuron according to a Poisson process with rate λ i . To model the different calcium levels at each event (e.g. due to multiple spikes or to adaptation [ 106 ]), we sample the overall concentration as coming from a unit log-normal distribution (i.e. an exponentiated Normal distribution N 0 , 1 ). Once the spike times are obtained, an auto-regressive model with p degrees of freedom (AR- p ) is used to simulate the calcium and fluorescence impulse response. As a difference equation, AR- p models can be written as y n = ∑ i = 1 p a i y n − i + b x n , where the a i ’s are the AR coefficients and b is a scalar multiple of the input. The impulse response can be obtained by solving the inverse Laplace transform y n = h n ∗ x n , h n = L − 1 b z p z p − ∑ i = 1 p a i z p − i . Standard linear systems theory shows that h [ n ] will be composed of the exponentiated roots of the characteristic polynomial z p − ∑ i = 1 p a i z p − i and therefore will be an exponentially decaying function. Higher order polynomials can result in a rise time as well. For this work we find that an AR-2 model ( p = 2) sufficiently models the rise and fall of observed GCAMP responses. The filter h is convolved with the spike-time vector to create the temporal activity per neuron. 4.3.3. [ Ca 2+ ] dynamics The fluorescence of a cell is dependent on the number of calcium ions bound to the indicator. if we denote [Ca 2+ ] as the amount of free calcium in the cell and [B] as the number of proteins in the cell, we can use the binding/unbinding dynamics, coupled with the entry/exit dynamics of [Ca 2+ ] in the cell to determine the fluorescence level at any given time. Specifically, we use the nonlinear diffusion of [Ca 2+ ] d Ca 2 + d t = − γ 1 + κ s + B K d Ca 2 + + K d 2 − 1 Ca 2 + − Ca 2 + r e s t , where [Ca 2+ ] rest represents the baseline free [Ca 2+ ] , γ is the [Ca 2+ ] diffusion constant, κ s is the endogenous [Ca 2+ ] binding ratio, and K d is the protein binding affinity constant [ 69 , 70 ]. As γ is a function of the volume-to-surface area, we use a different γ value for dendrite dynamics as for dynamics in the soma [ 107 ]. While this model permits simulation of the [Ca 2+ ] concentration over time, the model does not include the on/off time constants τ on and τ off that describe how long it takes for the bound proteins to become active. We can model this effect, as in [ 70 ], by convolving with a double-exponential function (4) h t = A 1 − e − t / τ o n e − t / τ o f f where the amplitude A and the time constants τ on and τ off can be fit to the particular protein kinetics. The final step in simulating the fluorescence time-traces is to convert the calcium concentrations to fluorescence levels. For this task, we use the Hill equation (5) Δ F / F = 1 1 + K D / Ca 2 + n H where the parameters K D and n H have been measured in the literature (specifically [ 72 ], Table 1 ), and the absolute florescence is (6) F = F 0 Δ F / F + F 0 where the baseline fluorescence F 0 can be tuned to the protein statistics. 4.4.

Optics simulation

The optics module consists of modeling the shape and intensity of the point-spread function (PSF) within the scanned tissue. For computational purposes, we assume the shape of the PSF is constant across the scanned volume and only the amplitude is modulated. We estimate the PSF within the scanned tissue by propagating a specified field through the simulated tissue across the field of view. We describe the scalar field at the front aperture of the objective lens as a Gaussian with a circular aperture and spherical phase: (7) U 0 ρ = exp − ρ 2 ρ e 2 − i k ρ 2 2 f ρ ⩽ ρ 0 0 ρ > ρ 0 where U 0 is the scalar field, ρ = x 2 + y 2 is the polar position, k is the wavenumber, ρ 0 is the radius of the objective lens, ρ e is the radius of the excitation beam, and f is the focal length of the objective lens. The wavefront is multiplied by any additional specified aberrations due to the microscope or the sample: (8) U ρ , θ = U 0 ρ exp i k ∑ i = 1 ∞ a i Z i ρ , θ where θ is the polar angle, a i are the Zernike coefficients, and Z i are the Zernike polynomials. By default, only spherical aberration approximating the contribution of the refractive index mismatch of the sample and astigmatism approximating the contribution of offset scanning galvanometers are included. The field U ( ρ , θ ) is propagated through the sample to the focal plane along a 2D grid of positions within a simulated refractive index volume δn . The volume δn is generated from the simulated vasculature and a 3D Gaussian process with a weight distribution approximating the refractive index distribution of mouse cortical tissue (see Supp Fig. 11 ) [ 108 ]: (9) δ n x , y , z = n d i f f V x , y , z + G P x , y , z where V is the vasculature and GP is the smooth Gaussian Process representing the optical properties in the non-vasculature areas. The vasculature provides the bulk of the long range refractive index shifts in the simulation, while the Gaussian process approximates the local shifts. The Fresnel diffraction integral is used to estimate the field throughout the volume, and the split-step beam propagation method [ 109 ] is used to apply the effects of inhomogeneity within the volume. The simulated phase-difference volume is summed into optical phase masks corresponding to each propagation step: (10) ϕ x i , y i , z i = k ∫ z i z i + 1 δ n x , y , z d z . This quantity is multiplied after each optical propgation step as (11) U i + 1 = e i k z i λ z i + 1 − z i ∫ ∞ ∞ ∫ − ∞ ∞ e − i ϕ i U i e i k 2 z i + 1 − z i x i + 1 − x i 2 + y i + 1 − y i 2 d x i d y i , where ϕ i = ϕ ( x i , y i , z i ) is the optical phase mask and U i = U ( x i , y i , z i ) is the scalar field at each position. The resultant 3D field generated by the propagation is then used to calculate the two-photon PSF: (12) P S F x , y , z = U x , y , z 4 . The aberrations caused by the phase differences approximate the effects of wavefront distortions caused by refractive index inhomogeneity within the imaged sample [ 110 , 111 ]. The two-photon PSFs at each location across the field of view are averaged to obtain the PSF to be scanned through the simulation, and the summed intensity of the PSFs across the field are used to generate an intensity scaling mask for scanning. For runtime considerations, the PSF near the focal plane is sampled at the resolution of the volume while the out of focus PSF and scaling mask is sampled at a reduced resolution. For alternative optical setups, we adjust the input field U 0 accordingly. For a low numerical aperture excitation beam, ρ e is reduced, and for a Bessel beam excitation U 0 is replaced with an excitation ring. See supplementary information for more details. An additional optical mask is also calculated by estimating the reduction in signal from absorption of the collected light by the vasculature. The collected light at each scanned position is reduced by a collection cone corresponding to the simulated collection objective numerical aperture: (13) r z = tan sin − 1 N A / n z where r ( z ) is the collected cone radius as a function of depth, and: (14) C x , y = 10 A ∫ 0 d r z ∗ V x , y , z d z where C is fraction of light collected, d is the tissue depth, and A is the adjusted light absorbance of light emitted from GFP normalized by the arterial blood absorbance factor. This absorbance mask is multiplied to the optical excitation mask to give the combined spatial signal scaling mask. 4.5. Scanning in silico The final module takes the generated volume, the generated PSF and time-traces, and generates the TPM output frames. The first step here is to use the time-traces and fluorescence distribution for each neuron to “color in” the corresponding volume with the current fluorescence level for that neuron. Similarly, the background level is set by repeating this process with the neuropil. The PSF is then convolved with the current volume, and the result is masked with the optical path mask to create an initial image. To simulate motion in the movie, we select a portion of this initial frame to treat as the entire image. The starting position (upper left corner) for the with-motion frame is moved according to a small ± 0.5 µ m jitter with occasional larger jumps (up to 2–3 µ m). Options to include per-line motion and shearing are also implemented by choosing different sub-sections of each row as the with-motion frame is extracted from the larger motionless frame. This frame represents the fluorescence level at each point in the sampled image. To obtain the actual electrical signals sampled by the TPM device, we apply a noise model that simulates the number of photons incident on the array (modeled as Poisson) followed by an electrical noise model that is Gaussian, with increasing mean and variance with larger numbers of incident photons. If λ is the true florescence for a pixel, x is the number of incident photons, and y is the measured electrical signal, the noise model can be expressed as x ~ Poisson λ y ~ log N μ 0 + μ x , σ 0 2 + σ 2 x where µ 0 and σ 0 3 are the baseline noise mean and variance (with no photons), and µ and σ 2 are the parameters controlling how the measurement mean and variance grow with increased incident photons. As a final step, we simulate the analog-to-digital accumulators’ property where photons arriving in one pixel’s accumulation time can cause an analog shape that bleeds through to the accumulation for the next pixel (Fig. 16). We simulate this effect by noting that if a photon arrives early in the sample period, then the analog PMT response g ( t ) is completely inside of the sample period and no bleed-through occurs. On the other hand, if the photon arrives within ∆ of the end of the sampling period, where ∆ is the temporal extent of g ( t ) (Fig. 16), then the tail end of g ( t ) that continues beyond the end of the period is integrated into the next sample. The probability of a given bleed-through level for one photon can thus be quantified as (15) p b = 1 − Δ T f o r b = 0 Δ T 1 g T − τ b f o r 0 < b ⩽ ∫ 0 Δ g t d t , where τ ( b ) represents the delay τ that is needed to result in a given bleed-through b . Since the relationship between b and τ , b = ∫ T − τ Δ g t d t , is monotonically increasing when g ( t ) ⩾ for all t , τ ( b ) is a well defined function. Since photon arrivals are approximately independent, the bleed-through probability distribution for multiple photons is the convolution of the distribution for a single photon. The resulting statistical model then takes a random fraction (uniformly chosen between zero and 50%) of each pixel with probability 0.2, and adds that amount to the next pixel, (16) p b = 1 − P b l e e d f o r b = 0 p b l e e d b m a x f o r 0 < b ⩽ b m a x Vasculature scanning To compare the axial spread of the PSF in deep volumes, we simulated datasets comparable to deep z-stacks of blood vessels [ 112 ]. The real z-stacks were taken in mouse neocortex at 830 nm excitation wavelength and a numerical aperture of 1.0 with 1 µ m spatial sampling at depths up to 450 µ m. For our analyses we generated local vasculature z-stacks (60 µ m) at a depth of 330 µ m by using the scanning module to iteratively scan a uniformly labeled simulated blood vessel volume. To adjust for the scattering at the new imaging wavelength, we used the equation provided in [ 108 ] along with a quadratic relationship between scattering length and refractive index. These values reasonably approximated values provided in [ 113 ] for nearby wavelengths. Additionally, we generated simulated datasets with randomly placed point sources (0.5 µ m in size) that were scanned in a local z-stack. We considered the axial and lateral spread of point sources (point source volumes) and line sources (capillaries in blood vessel volume) around 330 µ m in depth for this analysis. For point sources we isolated a 7×7 pixel subregion centered around a point source, background subtracted the isolated values and estimated the axial full-width half-max (FWHM) of each source. For line sources capillaries were isolated into 30 × 30 pixel subregions, rotated, and a 5 µ m central region of the brightest slice was used to estimate the lateral spread (FWHM). The central 5 µ m subvolume was background subtracted and used to estimate the axial spread (FWHM). The axial PSF was estimated as the difference between the axial spread and half the lateral spread, an estimate which assumes the width of the lateral PSF and the line sources are roughly comparable. Hemoglobin absorption To calculate the absorption due to hemoglobin, we assumed default concentrations of 150mg/ml Hemoglobin (Hb), 64500 g/mol Hb, and 2.9 (abs/ µ m)/(mol/L) in units of abs/ µ m. The absorbance was then calculated using Scott Prahl’s Hb curve [ 114 ] and eGFP emission spectrum [ 115 ].

Estimation of per-trace noise variance

To estimate the noise variance for each time-trace, we begin with the basic per-pixel noise model (17) y = Φ f + d + ϵ , where the noise is heteroskedastic in that the variance is proportional to the mean (18) ϵ ~ N 0 , diag Φ f + d . The least-squares estimate of the activations under the imperfect spatial profiles Φ ^ is (19) f ^ = arg min f y − Φ ^ f 2 2 (20) = Φ T Φ ^ − 1 Φ ^ T y (21) = Φ ^ T Φ ^ − 1 Φ ^ T Φ f + d + ϵ The covariance of this estimate is then (22) Cov f ^ = Cov Φ ^ T Φ ^ − 1 Φ ^ T Φ f + d + ϵ (23) = Cov Φ ^ T Φ ^ − 1 Φ ^ T ϵ . (24) = E Φ ^ T Φ ^ − 1 ϵ ϵ T Φ ^ Φ ^ T Φ ^ − 1 . (25) = Φ ^ T Φ ^ − 1 Φ ^ T E ϵ ϵ T Φ ^ Φ ^ T Φ ^ − 1 . (26) = Φ ^ T Φ ^ − 1 Φ ^ T diag Φ f + d Φ ^ Φ ^ T Φ ^ − 1 .

Calcium analysis algorithms

All analyses were computed with the 2017 versions of CNMF [ 52 ] and Suite2p [ 18 ]. Parameters were adjusted manually for each algorithm to optimize the output. For CNMF we used fr = 30, tsub = 5, patch size = [ 40 , 40 ], overlap = [ 8 , 8 ], K = 7, tau = 6, p = 0, and nu_bg = 1. For Suite2P we used diameter = 12, DeleteBin = 1, sig = 0.5, nSVDforROI = 1000, NavgFramesSVD = 5000, signalExtraction = ‘surround’, innerNeuropil = 1, outerNeuropil = Inf, minNeuropilPixels = 400, ratioNeuropil = 5, imageRate = 30, sensorTau = 0.5, maxNeurop = 1, sensorTau = 0.5, and redmax = 1. For PCA/ICA we used fr = 30, ssub = 2, tsub = 10, nPCs = 1000, smwidth = 3, thresh = 2, arealims = 10, mu = 0.5, dt = tsub/fr, deconvtau = 0, spike_th = 2, norm = 1. While figures displayed typical outputs (i.e., denoised traces from CNMF), all quantitative comparisons were computed using the raw DF/F traces returned by each algorithm. Local correlation calculation To calculate the local correlations, V1 two-photon recordings and simulations were motion corrected using correlation-based rigid motion correction [ 116 ]. A 1500-frame subsection of each dataset over a 250 × 250 pixel area was extracted. For each pixel, the Pearson correlation between its fluorescence activity and that of each of the neighboring pixels in a 51 × 51 pixel square neighborhood (up to 25 pixels away in each direction) were calculated. The results were averaged over all pixels (to create the mean images) and histograms were created to depict the spread for correlations along the fast-scan direction. Computation of the ideal spatial profiles The ground-truth template for each cell was obtained from NAOMi by lighting up that individual cell and zeroing out all other cells, and then scanning the volume at the simulated imaging depth plane. We further refined this profile to the ‘visible’ pixels by removing any pixel whose estimated noise variance was larger than the signal. Such a ‘binary mask’ is common, for example, in increasing intelligibility for speech data by applying the binary mask in the time-frequency domain.

Calculation of auxiliary time-traces

To calculate the auxiliary, noisy “ground truth” time traces, we used Profile assisted least squares (PALS), which uses the ground truth ideal profiles in a least-squares estimation. We considered the movie frames y t for t = 1… T and the calculated ground-truth spatial profiles X = [ x 1 , …, x N ]. The noisy time trace estimates were then calculated via the least-squares estimation procedure at each time-step t (27) s ^ t = arg min S y t − X s 2 2 .

Fitting

GCaMP parameters Fitting

GCaMP6f and GCaMP6s parameters was accomplished by using either constrained optimization (interior-point implementation in MATALB) or BayesOpt [ 117 ]. Data used to fit these models were downloaded from [ 118 ] and are available on the Collaborative Research in Computational Neurosciences (CRCNS) web portal at https://crcns.org/data-sets/methods/cai-1/?searchterm=chen%20gcamp . Baseline estimation Baseline fluorescence estimation is critical for assessing how much of the total acquired signal originates from each individual cell, which in turn permits more accurate ∆F/F estimates. More accurate estimates both reduce both reduce the expected noise level, and also produce time traces that are more meaningful in that they directly relate to the actual intracellular calcium levels and spike rates. This is especially true in preperations with variable expression levels. Accurate baseline estimation is challenging with no clear present solution. Analyzing current methods for baseline estimation, however, can provide an approximate error range of this measure, which is vital for establishing error bars in the scaling of spike rate estimates and other factors for detemining data quality for ensuing analyses. We analyzed ∆F/F distributions and estimates per pixel (pixelwise estimation) to obtain a range for the values that this ratio can take ( Sup. Fig. 38 ). To estimate the baseline fluorescence level for each neuron we first calculated its spatial profile S ( k ) and the modal image M ( k ) for each pixel k. We assumed M ( k ) = c * S ( k ) + B ( k ) , where B ( k ) is the contribution to the fluorescence from all other sources and the baseline fluorescence level is given by c * S ( k ) . When the neuron was active, the pixel values were approximated as F ( k ) = ( dF/F ) * c S ( k ) + B ( k ), in the absence of background activity, which we assumed for this estimate. We estimated c * S ( k ) by plotting a scatter-plot of S ( k ) against M ( k ) and noting that for a given S value, there exist a range of possible M values with a minimum value that increases as S increases. These minimum values correspond to the case where the contribution by B is minimal and approximately uniform, corresponding to fluorescence from the tails of the PSF. By fitting the line of these minimum values, we used this slope to approximate the value of c . Care may be taken to estimate this minimum B contribution by examining the values of pixels immediately outside the spatial profile. We used PCA to locally estimate S ( k ) , which better approximates these values immediately outside of the neuron, improving estimation accuracy. For each pixel k , we note that B(k) must be positive. Therefore, for each pixel k , c