A flow-based latent state generative model of neural population responses to natural images
Authors: Mohammad Bashiri
Edgar Y. Walker
Konstantin-Klemens Lurz
Akshay Kumar Jagadish
Taliah Muhammad
Zhiwei Ding
Zhuokun Ding
Andreas S. Tolias
Fabian H. Sinz
Venue: NeurIPS 2021 Spotlight (2021)
Keywords: mouse visual cortex
neural system identification
latent variable models
normalizing flow
generative models
noise correlations

Abstract: We present a joint deep neural system identification model for two major sources of neural variability: stimulus-driven and stimulus-conditioned fluctuations. To this end, we combine (1) state-of-the-art deep networks for stimulus-driven activity and (2) a flexible, normalizing flow-based generative model to capture the stimulus-conditioned variability including noise correlations. This allows us to train the model end-to-end without the need for sophisticated probabilistic approximations associated with many latent state models for stimulus-conditioned fluctuations. We train the model on the responses of thousands of neurons from multiple areas of the mouse visual cortex to natural images. We show that our model outperforms previous state-of-the-art models in predicting the distribution of neural population responses to novel stimuli, including shared stimulus-conditioned variability. Furthermore, it successfully learns known latent factors of the population responses that are related to behavioral variables such as pupil dilation, and other factors that vary systematically with brain area or retinotopic location. Overall, our model accurately accounts for two critical sources of neural variability while avoiding several complexities associated with many existing latent state models. It thus provides a useful tool for uncovering the interplay between different factors that contribute to variability in neural activity.

A flow-based latent state generative model of neural population responses to natural images
Mohammad Bashiri,'' 1." Edgar Y. Walker, 1," Konstantin-Klemens Lurz,'
Akshay Kumar Jagadish, 122 Taliah Muhammad," Zhiwei Ding, 3 Zhuokun Ding,
Andreas S. Tolias, Fabian H. Sinz5...
Institute for Bioinformatics and Medical Informatics, University of Tübingen, Germany
Max Planck Institute for Biological Cybernetics, Tübingen, Germany
Department for Neuroscience, Baylor College of Medicine, Houston, TX, USA 4 Center for Neuroscience and Artificial Intelligence, Baylor College of Medicine, Houston, TX, USA 5 Department of Computer Science, University Göttingen, Germany
equal contribution, sinzzcs. uni goiereinnnngggeetine....eeee
Abstract
We present a joint deep neural system identification model for two major sources neural variability: stimulus-driven and stimulus-condittoond fluctuations. To this end, combine (I) state-of-the-art deep networks for stimulus-drive activity and (2) flexible, normalizing flow-based generative model to capture the stimulus- conditioned variability including noise correlations. This allows model end-to-end without the need for sophisticated probabilistic approximations associated with many latent state models for dimmulu--onditioons fluctuations. We train the model on the responses of thousands of neurons from multiple areas of the mouse visual cortex to natural images. We show that our model outperforms previous state-of- ofethe-rrt models in predicting the distribution of neural population responses novel stimuli, including shared thermore, it successfully learns known latent factors of the population responses that are related to behavioral variables such as pupil dilation, and other factors that vary systematically with brain area or retinotopic location. Overall, our model accurately accounts for two critical sources neural variability while avoiding several complexities associated with many existing latent state models. It thus provides a useful tool for uncovering the interplay between different factors that contribute to variability in neural activity.
train the
variability. Fur-
Introduction
Characterizing the activity of sensory neurons is major goal of neural system identification. While neural responses in the visual cortex vary with visual stimuli, they also exhibit variability to the
repeated presentations of identical stimuli [14. This ttiuulus-onnditions variability has significant
and sophisticated correlations among neurons commonly referred as noise correlations 6
and exhibits dependency on various factors such as the stimulus (70), the behavioral task [10] 11]
attention 12414), and the general brain state [15] Understanding the nature this correlated variability and its functional implication in the processing of sensory stimuli requires models that account for both cimmulus-driven and shared stimulus conuuuuu-cooiitioneed variability. The goal is thus to model the dimulus-conditionns response distribution p(r x) of population activity E R" over neurons responding to an arbitrary sensory stimulus x. However, models that account for stimulus- driven and cinullu--onniitioned correlated variability have been developed largely independently.
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

In the recent decade. we have seen significant progress in modeling stimulus-driven activity, largely driven by the use of deep neural networks (DNNs) [7 22.
Typically, the expected response of the neurons conditioned on the stimulus is captured as a function of the stimulus via a deep network fx(x) = E[r|x] with learnable parameters 0. These models can therefore predict how population responses depend on an arbitrary stimulus, and could even be used to derive stimuli that would yield desirable responses [23] [4]]. Typically, these networks are trained using Poisson-loss, assuming that the population activity r is distributed around the stimulus-conditiones mean fx(x) with an independent Poisson distribution. Therefore, existing state-of-the-art networks commonly ignore correlations among neural responses, and impose strong assumptions about the form of the marginal distribution (i.e. Poisson) for each neuron. As sensory populations are known to exhibit noise correlations and deviate from Poisson distributions 4 25 26], this conditional independence assumption might limit the ability of these models to accurately capture p(r|x).
On the other hand, many of the existing models for stimulus-conditioned variability capture the variations in the population activity by specifically modeling the responses to repeated presentations of an identical stimulus. Many of these approaches employ statistical techniques such as maximum- entropy or copula distributions to reduce the number of parameters needed to fit the target distribution (27-29] A popular approach has been to describe the stimulus-conditionee variability in terms of a typically lower-dimensional shared latent state z: p(r[x) = s p(r X z)p (z x) dz [16] 25 26 B5 Among these hierarchical generative models that can capture more sophisticated relationships between the stimulus and noise correlations, as well as deviations from Poisson, such as 351. While these approaches present powerful methods to capture stimulus-conditioned variability, they often fit p(r)x) separately for each unique stimulus and require responses repeated presentations of the stimulus [16 25 26
over-dispersion 225 26 32
34
29 35 This limits their ability to yield predictions to a novel stimulus without requiring some stimulus stimulus-specifet parameters to be learned. Furthermore, the increased complexity of the distribution usually requires substantially more involved probabilistic machinery to make latent state inference and parameter fitting feasible. Consequently, most latent state models for neural data either ignore stimulus-driven variability altogether 30 31 34]. or employ a very simple model of stimulus driven variations [16 25 26 32)
Here, we propose a new model that closes the gap between these two approaches by combining DNN- based models of stimulus-driven activity with a latent state model that accounts for shared stimulus- conditioned variability. While DNNs can be trained effectively via madient-based optimization, the challenge is to avoid the complex probabilistic machinery associated with existing latent state models, particularly those that require stimulus-specins of identical stimuli. To this end, we combine normalizing flows [36|4]] with Gaussian Factor Analysis (FF) models (42), where the stimulus-dependencce occurs through DNN that learns to shift the mean of the FA distribution based on the stimulus. FA models make use of multivariate Gaussian distributions with a particular low-rank structure of the covariance matrix. While the use of FA in capturing shared variability greatly simplifies inference and learning, it is not directly applicable to neural responses because neural responses are not diissiibndddttuueed, particularly for low firing rates. To circumvent this problem, variance-stabilizing transformations. such as the square-root function, have been used in the past to make the responses more naussian-ddisributeed [6] However, there may be other transformations that capture the response distribution more accurately. Furthermore, since the transformation for one neuron may not be applicable to other neurons, ideally it would be learned for each neuron separately. To achieve this flexibility, we allow our model to learn neuron-speciff transformations with a marginal normalizing flow.
parameters be learned over repeated presentations
Normalizing flow models are density estimators that use a series of diffeomorphisms to transform the source density underlying the data into a simple distribution typicalll an isotropic Gaussian of the same dimension. These transformations are usually chosen to have log-determinants, and typically act on the entire variable vector to capture any statistical dependencies between the dimensions. Here, we replace the isotropic Gaussian with an FA model to capture dependencies among dimensions and only use diffeomorphisms that act on each dimension separately, i.e. apply flow-based transformations on the marginals only. While this choice places certain restrictions on the complex dependencies between neurons that may be captured (refer to section 4 Discussion for details), it has two important advantages: (1) The generative model is easy to train while combining state-of-the-an deep networks with flexible latent state models, and (2) the use of marginal flows allows for easy mechanism to compute conditional distributions of one neuron given responses of other neurons that would not be easy to obtain with non- morginaaa flow models.
2

---

latent model
transformed responses "V"
neural responses 'r"
Flow-based Factor Analysis model. Schematic of the flow- based model relating all relevant variables in the study. b: Schematic of the sub-network used by the image-computable model map cortical positions into receptive field positions. Refer to section Section[[Mein Methodd for the details.
Figure 1:
latent states "
covariance
Cc
T
neuro
Vansform
T
Falx) moan
neuron 1
neuron
image "x'
(e)
b
readout position network
fr (as)
inage-computable model
cortex
receplive field posttions "6"
positions
In summary, we make the following contributions. We (I) combine state-of-the-an DNN-based models with flow-based latent state models jointly account for stimulus-driven and shared stimulus- conditioned variability in neural population activity. Our model can predict the distribution of neural responses to unseen stimuli, without the need for repeated presentations to learn commulus-conditionn variability. We (2) apply our method on the activity of thousands of neurons in response to natural images, recorded via two-photon Calcium imaging from multiple areas of the mouse visual cortex. We demonstrate that our model outperforms current state-of-the-ard methods in capturing the distribution of responses. Finally, (3) show that model infers latent state structures with meaningful relations to behavioral variables such as pupil dilation as well as other functional and anatomical properties of visual sensory neurons.
2 Methods
2.1 Models
Flow-based Factor Analysis model (FlowFA) For a given stimulus x and population response
E R", where n is the number of neurons, we define our normalizing flow-based Factor Analysis (FlowFA) model of the cimulus-conditionn population activity p(r[x) as
p(r|x,0, Φ) (T); CC + IdeV
(I)
FlowFA has two major parts: (1) A flow model To with learnable parameters à that transforms the population responses such that the transformed responses v (r) are well modelled by a (2) Gaussian Factor Analysis (FA) model A ((v; (x),, CC + v)) (Fig. 1p) Here, (Nv; Σ) denotes Gaussian distribution over with mean p and covariance E. According the FA model, the random variable v is generated via v = fo + Cz + S where z E R* is a low-dimensional latent state with k & n and an isotropic Gaussian prior NN [x) whose samples map via the factor loading matrix C Re*k The effect of the stimulus
the responses is captured by the mean of the FA
distribution that depends on the stimulus, modeled as a deep network fg (x) E RR with learnable parameters 9 (Fig. bb.. We further include neuron-specifie, independent noise NN00, v where E R is a diagonal covariance matrix.
Since the flow model is trainable change of variables, it introduces the absolute determinant (r)) the Jacobian V of To with respect to r into Eq. (T. The transform itself diffcomorphism, i.c. invertible differentiable mapping Tş RR" + R" allowing to evaluate the exact likelihood of each data point and easily draw samples from the model. Therefore, the model serves as a fully generative model from which samples of the stimulus-conditiones population responses can easily be generated for an arbitrary stimulus.
In the model formulation presented here, we choose T to act on each single dimension separately, ie. To (r) = [To, (r.),.., To. (r.))).. This choice results in a diagonal Jacobian which not only substantially simplifies the form of determinant to detV,Ta(r) = Th ăr. us easily compute conditionals and marginals (see appendix A for the details). This would not generally be possible for diffcomorphisms with non-diagonal Jacobian.
ST& but also allows

---

Zero-Inflated Flow-based Factor Analysis model (ZIFFA) For two-photon Calcium imaging, sigmificant portion of inferred neural activity is zero, resulting in a sharp peak at zero the response distribution (i.e. zero-inflated distribution) 43 This zero-inflation is potentially a problem for the FlowFA model since the model would attempt to generate the peak at zero by mapping large proportion of the Gaussian probability mass onto the "zero" responses, resulting in a poor fit to the response distribution. To avoid this, we extend FlowFA by modeling the zero responses with a separate peak (similar to Wei et al. (43) and applying the FlowFA model to capture only the positive responses, We refer to this model as Zero-Inflated Flow-based Factor Analysis (ZIFFA). More specifically, ZIFFA is mixture model that models neural responses below and above threshold value p with two separate, non-overlapping distributions. To capture the peak at zero, the responses below the threshold (i.e. "zero" responses) modeled by uniform distribution, while FlowFA is used to capture responses above the threshold:
p(r(x)
II (B:rsppp
1 - (x)
gi (xx) NCTs(E)): C+ C.C. + + V+)+)T((++)) (i:re>p)
II
P
(2) where (x) is the probability of the response being above the threshold p modeled, jointly with the mean of the FA, as a function of the stimulus via a DNN fo with learnable parameters A. r+ and fo,+(x) the sub-vectors, and C+ and V are the uub-matrices corresponding to responses above the threshold, and 0, C, v are the same as defined in Eq. ① Refer to appendix B for the derivation.
Control models We compare the FA-based models against two control models used for neural system identification that assume independence among neurons with specific forms marginal distributions inspired by existing work: (1) Poisson [18] [2] and (2) Zero-inflated Gamma (ZIG) [43). To capture continuous neural responses measured with Calcium imaging, we relax the discrete Poisson distribution into continuous distribution by assuming f + where ~ Poisson(A) and < Uniform(0, 1). This yields the likelihood function
Ppaiss (r(x)) = II )e(x)lnde-
[rs]!
where A(x) = fx) is the predicted firing rate of the neurons to input image x modeled as a DNN fa with learnable parameters 8. The ZIG distribution is mixture of a uniform and gamma distribution separated at the value p with no overlap [43]:
II (x()
exp (- Ti ))
Pri (r x)=
(4)
[x)))
1x(x))
where (x( is the scale parameter of the gamma distribution, and (x) is same as in Eq. (2. To formulate ZIG as mnage-computable model, (x(x) and q (x) are jointly modeled using DNN fo with learnable parameters f. Similar to Wei et al. [43]. we let the shape parameter N2 be neuron- specific, but independent of the input. Importantly, we used the same value for p in both ZIG and ZIFFA models.
Note that when the covariance matrix of the FA-based models is diagonal (i.e. O-dimensional latent state), these models assume independence among neurons and their performance directly comparable to the control models.
2.2 Model components
Deep convolutional neural network fo We capture the stimulus-driven changes in the neuronal response distribution using deep convolutional neural network fx()) with the same architecture as used by Lurz et al. [22] Briefly, the network consists of two parts: (1) shared four-layer core network, where each layer consists of a standard or depth-separable 44 convolution operation resulting in 64 feature channels, followed by batch normalization and ELU nonlinearity, and (2) a neuron-specific readout mechanism (referred to as "Gaussian readout") that learns the position of the neuron's receptive field (RF) and computes a weighted sum of the features at this position along the channel dimension (Fig. [p). In contrast to Lurz et al. 22 where the RF positions ò in image space were obtained by applying a shared affine transformation on the experimentally measured cortical positions A of the neurons, here allow this mapping to take on non- linear form to allow flips

---

in the representation of the visual field as function of cortical position (Fig. [p). to model vortex-to-visual space mappings for multiple brain areas, as the retinotopy of some areas are mirrored with respect to each other. During training, we apply 1.1 regularization to the readout feature weights and 1.2 regularization on the Laplace- Laplace-filtered filtered weights of the first convolution layer.
This is crucial
Normalizing flow T' We construct the marginal flow model T. = affine o exp oaffine o ELU o affine o E affine log caffine from a set of monotonic functions {affine, ELU, exp). of which only the affine transformation has learnable parameters. We restricted all the affine transformation layers to have positive scale, and additionally restricted the first affine layer to have a positive offset. For each neuron indexed by i. we leam separate marginal transformation To. We compare the flow transformation against two common fixed transformations: square-root [16] 30] and Anscombe [45]. These two transformations can be expressed by the general form u = exp(a log + b) + c) which is a series of affine, log, affine, and exp transformations, with a = 0.5, b O, and c = 0 for square-root, g. and c (2g for Anscombe. We specifically chose the components of To such
and = 0.5, b- 3 that these common fixed transformations exist as special cases, ensuring that the flow transformations are strictly more flexible than any choice of fixed transformations commonly found in the literature. For ZIFFA, we adjusted the formulation of the marginal flow To such that the predicted neuronal responses remain above p. the boundary between the uniform and the FlowFA components of the mixture model, by replacing the first affine transformation in Ta with a layer that only shifts by -p.
2.3 Neural and behavioral data
We recorded the response of neurons in mouse visual cortices (layer L2/3) to gray-scale natural images using wide-field two-photon microscope [46] (see appendix C for details). In this study, we used two scans from two mice spanning three visual areas: primary visual cortex (V1) and lateromedial area (LM) in scan 1: V1 and posteromedial area (PM) in scan 2. A total of 2,867 VI neurons and 907 LMM neurons were recorded in scan 1; 5,029 V1 neurons and 3,343 PM neurons were recorded in scan 2. Among these, we used 1,000 V1 and 907 LM neurons from scan 1, and 1,000 V1 and 1,000 PM neurons from scan 2. For both scans, neurons were randomly selected the area contained more than 1,000 neurons. We also recorded behavioral variables such as pupil dilation, simultaneously. The natural image stimuli were sampled from ImageNet [47), cropped to fit a monitor with 16:9 aspect ratio, and presented to the mice at resolution of 0.53 ppd (pixels per degree of visual angle). A total of 6,000 images were shown in each scan, of which 1, 1,000 images consist of 100 unique images each repeated 10 times to allow for estimate of the neural response variability. We used the repeated images for testing, and split the remaining images into 4,500 training and 500 validation images.
2.4 Model fitting and evaluation
Fitting We trained all models and-to-end via gradient- based optimization to maximize the log- likelihood obtained from Eqs. Q (2) (3) or (4) for the corresponding model, optimizing over all learnable parameters. To ensure that V. the diagonal covariance matrix, stays positive-valued, we re-parameterized V e" and optimized instead. To find the best image-computable DNN models, we used Bayesian optimization (8 to find hyper-parameters that maximized the final log-likelihood of the trained model. prpee-pprraeters include the learning rate and regularization coefficient on the readout weights. The log-like liooeloooood used for scheduling learning rate, early stopping, and finding was computed on the validation set. Additional details about training can be found in appendix D. The code can be found at httptttsss::sgigiee
2022 com/ sinzlab/ www.annlllabbabbbiiiaa
Evaluation We compared the FA-based models (ZIFFA, FlowFA, and FA with fixed transformations) to the control models based on likelihood and leave-neuron-out prediction correlation on the test set. For the former, we computed the likelihood of the responses in bits per neuron per image under each model, based on Eqs. 1 (3), and ( accordingly. For the correlation measure, we computed the Pearson correlation between the predicted and the measured responses of each neuron on the test set. For the FA-based models that may capture the statistical dependency (i.e. covariance) between neurons, we predicted the response of a given neuron conditioned on the responses of all other neurons recorded simultaneously response of all other neurons r, we estimated the response of a neuron r; to the image by computing the posterior mean of the neuron's 's response F(r: |x, r.d. We refer to this measure as conditional correlation (see appendix E for details).
the trial. More specifically, given an image and the
5

---

a
idiettiyy
sqrt
anscombe
example
example
b
neuron
neuron 1
10
Response
C
Model transformation
identity
0,001
0.002
0.004
0.016
0:051
0.527
0.004
scrt
0.003
0.001
0.002
0.027
0.013
0.082
0.010
0.001
0.007
0.056
0.013
0.005
0.052
anscombe
0.002
0.003
0.001
0.023
0.075
0.394
0.009
0.009
6.041
0.431
0.010
0.031
0.686
flow
0.001
0.001
0.001
Identity
sart
anscombe
example
example
example
example
examples
example 6
example
example
examples g
example 10
Simulated transformation
FlowFA model recovers the underlying transformation. Simulated responses for 2 neurons under various transformations. Across all transformations, transformed responses were sampled from Gaussian distributions with differing means (indicated by the color of the samples) but identical covariance. The covariance between the two neurons is shown in black text. b: Transforma- tions learned by the flow model are shown in black, overlaid on the ground-truth transformations. c: Performance of models with fixed or learned (flow) transformations (rows) trained on responses simulated with variety of transformations (columns). Cases where the simulating and trained transformations are the same are indicated by black outlines. Performance is measured as the KL divergence between the modeled and ground-truth distributions, where 0 would correspond to perfect fit.
Figure 2:
Results
3.1 Model performance
FlowFA model faithfully recovers invertible transformations on synthetic data We first used synthetic data illustrate that our FlowFA model with a learnable transformation can adequately learn and recover wide variety of transformations resulting in different response distributions. To this end, we sampled 5,000 data points for 100 neurons from models with different ground-truth transformations (see appendix F for details on data generation). The invertible transformations (Example 1-10) had the general form (a log(y + b) + c) with differing values of a, b, and e (Fig. [D). We trained FA-based models with either fixed (FixedFA) or learnable flow-based (FlowFA) transformation. As expected, the models with a fixed transformation performed well if the data was generated with a similar transformation, but the performance suffered when the transformations differed (Fig.2k first three rows). In contrast, the FlowFA model was able to ficxibly learn every underlying transformation (Fig. (2)) and effectively captured all distributions across all simulations (Fig.) c. last row).
Flow-based models capture cortical response distribution well After demonstrating that the flow- based model can effectively fit wide range of distributions, used to capture distributions the mouse visual cortex population responses to natural images, recorded in two different two-photon seans from two mice (scan and scan 2, refer to section 2.3 for details). We trained the FA-based models (ZIFFA, FlowFA, and FixedFA) for different values of latent dimensions & E {0, 1, 2, 3, 10). We measured the model performance by computing the log-likelihood as well as the conditional correlations (see section 2.4).

---

Scan 1
Scan 2
b
Scan
8.93
11.35
7.5
10
8.88
11.32
Log-llkellhodd
5,0
8.83
11.29
2.5
9.0
6.0000
11.40
00000
5
Scan
8.7
000
0
11.15
-2.
8.4
10
.0.90
# of latent dimensions
Scan 1
Scan 2
Response
Conditional Correlation
0.5
0.4
[Pearson]
0.4
0.3
8610
FA with learned transformation
FA with fixed transformation
Models with independence assumption
0.3
0.2
ZIFFA
sqrt anscombe identity
Poisson
0.2
0.1
0
FlowFA
ZIG
2 10
1
2 3 10
8 of latent dimensions
Figure 3: Comparison of models trained on the mouse visual cortical population responses natural images. log-likelihonn computed for models trained scan 1 (left panel) and scan 2 (right panel). Values for both individual (lighter shade) and average (darker shade) performance of a model trained under various random seeds are shown. Gray block provides zoomed-in view of the ZIFFA, FlowFA, and Prro-Inflated-aameee the flow-based models (ZIFFA in green, average across neurons in light green; FlowFA in pink, average across neurons in light pink) shown in comparison to fixed transformations. c: Conditional correlation. Format is similar to
(ZIG) models. b: Neuron-specific transformations learned by
The ZIFFA model outperformed all other models across all numbers of latent dimensions l: in terms of log- ligelikel (Fig. 3l). Furthermore, with increasing latent dimensions, the conditional correlation of the ZIFFA model improved significantly beyond the control models (Fig. 3p.) Interestingly, observed that the ZIFFA model exhibited slightly lower correlation performance compared to models with fixed transformations, reflecting that fitting models likelihood does necessarily yield optimal correlation. Importantly, the flow-based models outperformed all FixedFA models in terms of likelihood, which is corroborated by the fact that the learned transformation markedly differs from all fixed transformations and from one neuron to the other (Fig. 3b). Overall, the results suggest that the ZIFFA model is able capture the (marginal) neural response distributions more accurately than other models (Fig. S2) while at the same time it learns and takes advantage of the statistical dependencies between neurons.
3.2 Uncovering biological insights from the trained model
Here, explore the utility of our model in uncovering potential biological insights. All analyses were performed on the trained ZIFFA model with 3 latent dimensions.
Model-based visual area identification Several visual areas in mice show retinotopies that are "flipped" with respect each other 49). Intuitively, this means that if point moves along the cortical surface, as it crosses the boundary between two "mirrored" areas, its counterpart in visual space would reverse its movement direction. As described section 2.2] our model equipped with a component network that predicts the RF location å of each neuron in visual space function of its cortical location A (Fig. [])). This network can be used to infer distinct visual cortical areas by detecting where the retinotopy "flips" with respect to the cortical position. To detect this flip we looked the sign of the determinant of the Jacobian of the RF positions with respect to cortical positions det a6 aA* The sign can detect changes in the direction because (1) the sign of a determinant flips if one of the column row vectors of the Jacobian matrix flips and (2) the determinant is invariant under rotation. When we compare distinct areas identified via the model to the experimentally identified areas, find very good match (Fig. 44, left right panels). To assess the quality of the learned
7

---

Scan
Scan
Model-based
Model-based area identification
Cortex positions
RF positions
area identification
Cortex positions
RF positions
Cortex
Image
Cortex
Image
b
Latent
Latent
Latent
Latent
Latent
Latent
x10
Cortex
x
C
Latent
Latent 2
Latent
Latent
Latent
Latent
v1
V1
Image
LM
PM
Image
Image
d
Pupil dilation
Jadii Aww th Mallit (WMNA (MM "A (YM9
R=00000
0.00
y
Labent
R² =,50
Latent
0.02
0.09
Latent
0
250
500
750
1000
250
500
750
1000
Test trials (ardered time)
Figure 4: Analysis of the ZIFFA model with 3-dimensional latent state (k 3). a: Model-based area identification from responses of visual sensory neurons to natural images. Left panel (Cortex positions): cortical position of the recorded neurons color-coded by experimentally identified areas (green: V1; blue: LM; orange: PM). Middle panel (RF positions): learned receptive field position for each neuron function of cortical positions color-coded by experimentally identified areas. Right panel MModel-based area identification): visual areas identified via the model by computing the determinant of the relative changes in RF position with respect to changes in cortical position; blue color shows negative determinant (i.e. mirrored visual field representation) and red color shows positive determinant (i.e. mor-mirroeed visual field representation). b-c: Distribution of the latent-to- neuron weights across cortical positions (b) and receptive field positions (c). Pupil dilation (black) and the inferred latent states (red) across trials from the test set. R² values are computed between the inferred latent state and the pupil dilation.
mapping, we quantified how well our model can identify distinct visual brain areas via the sign of the determinant. Across models initialized and trained with different random seeds, the sign correctly classifies distinct brain areas with an accuracy of 84% + 3.4% (SEM) and 75% + 7.7% (SEM). Because the experimental methods to determine area assignment that we use as ground truth can be quite coarse, the actual accuracy could even be higher. This suggests that our model could in principle

---

allow neuroscientists to identify distinct visual areas from responses to natural images alone, without the need for an extra experiment for area identification.
Inferred latent states and their functional and anatomical implications We next explored the latent states and how they relate to anatomy or behavior. For any particular trial, the FA-based models allow us to infer the most probable latent state z (MAP estimate), where the effect of each latent dimension on the neural population is captured by the factor loading matrix C. However, as formulated in Eq. ① and (2) interpreting the inferred latent states can be difficult because the latent dimensions can be arbitrarily permuted and rotated (with corresponding changes in C) without affecting the fit of the model. To facilitate interpretability of the inferred latent states, we follow a similar procedure used by Yu et al. (30) extract orthonormalized latent states which are uniquely ordered by the amount of response variability each latent dimension accounts for (see appendix G for detailed explanation).
The orthonormalized latent states inferred from the ZIFFA model showed strong correlations with behavioral variables such as pupil dilation (Fig. 4J), as expected from previous works that use pupil dilation as a proxy for arousal and attention (50H54). Interestingly, pupil dilation correlated most strongly with the second latent dimension in both scans with It² values of 0.53 (p- < 0.001, two-tailed test for significance of correlation (55))) and 0.63 (p- < 0.001) for scan I and scan 2, respectively, comparable to values previously reported [56] To our surprise, this observation was consistent across models initialized and trained with different random seeds (Fig. S4b). To further quantify how well the latent states can jointly predict the pupil dilation. we regressed the pupil dilation against the latent states (Fig. S4a). The resulting R² values were 0.56 (p < 0.001) and 0.76 (p < 0.001) for scan 1 and scan 2, respectively. The high correlation between the latent states and the known surrogates of global brain state such as pupil dilation suggests that the latent model is able to learn meaningful dependencies and common factors in neural population.
Next, we explored whether the effect of the orthonormalized latent states on the neurons is related to their cortical or RF positions. To this end, we plotted the sign and magnitude of the weight mapping from the latent state to each neuron on the cortical position (Fig. 44)) or the RF positions of the neurons (Fig. 4b). We observed that the effect of some latent dimensions vary systematically across brain areas where the latent dimension has generally opposite effect on different areas (Fig, Ep: dimension 2 for both scans). In addition, some latent dimensions seemed to vary as a function of RF positionssretiiotopyy where a differential effect of the latent dimension is observed for both areas dimension for both scans). Interestingly, the first dimension which accounts for most of the shared variability in neural responses (refer to section G for more details) seemed to have a global effect that does not vary across different visual areas. These observations illustrate that our model can be useful tool for uncovering the functional and structural implications of the behavioral or internal processes associated with the inferred latent states.
While the result of the analyses we present here are promising, we would like to point out that all analyses are preliminary, and conclusive biological interpretations would require additional rigorous experiments and analyses.
Discussion
Getting the best of both worlds Two major components of the variability in the activity of cortical neurons are the variability due to stimulus and the variability due to unobserved or internal processes, such as behavioral tasks or general brain states, that affect population of neurons similar ways giving rise to correlated variability among neurons. Here, we presented a model that combines state- of-the-art DNN-based models to predict stimulus-driven changes in neural activity with a simple, yet flexible, flow-based factor analysis model to account for correlated neural activity. This formulation allows us to evaluate the exact likelihood of neural responses, easily sample responses, and efficiently compute conditional and marginal distributions of subsets of neurons. By fitting this model to the activity of thousands of neurons from multiple areas of mouse visual cortex in response to natural images, we obtained ofate-of-the-ar art performance in capturing neural response distribution while additionally yielding latent states that exhibit meaningful relations to anatomy and functional properties of visual sensory neurons.
Modeling zero-inflated response distribution Flow models use diffeomorphisms to map one distribution into another. However, diffeomorphisms cannot transform single peak at 0 0-typically
9

---

observed in neural responses recorded via Calcium imaging into a smooth distribution such as Gaussian used in our model. The ZIFFA model avoids this problem by only transforming the positive part of the response with a diffeomorphism while explicitly capturing the peak at 0 via uniform distribution as found in ZIG. Importantly, ZIFFA preserves all properties of the FlowFA model, while capturing the marginal distributions more accurately (Fig. S2), achieving higher likelihood (Fig. B and learning more consistent and less step-like transformations (Fig. S3).
Dependency of noise correlation on the stimulus The presented flow-based models learn a nonlin- ear transformation between a simple distribution (Gaussian FA) and the neural response distribution. While the learned covariance structure on the "transformed" neural responses captured by the FA model does not vary with the stimulus and the stimulus is only used to shift the mean of the FA model, this is not true for samples from the FA model transformed back into "neural response space" because the nonlinear flow transformation can introduce changes in the covariance as the mean varies (Fig- (b)). This mean-dependent change in the covariance potentially allows the model to capture changes in the covariance structure based on stimulus through the nonlinear transformation. A possible extension of our model is explicit dependence of the FA's covariance matrix the stimulus, which would allow the model to capture more complex dependencies between the stimulus and covariance structure.
Comparison to related methods Our approach in capturing stimulus-conditioned variability is related to many existing approaches, or can be seen as a generalization thereof, while being compu- tationally casier to handle at the same time. Recently, Keeley et al. [35] captured the bial-by-trial fluctuations by modeling the stimulus- specciff and trial-specific latents via Factor Analysis (FA) models much like in our model. Importantly, while we capture the dependence of the stimulus-specific latents on the stimulus explicitly via a trained DNN, they inferred it from repeated presentations of the stimulus. Furthermore, the final Poisson distribution used to map from the latents to the distribution of neurons can be captured in our model via the flow-based transformation (e.g. inverse Anscombe) that maps Gaussian-distribuce distribution. Moreover, the use of FA in combination with the marginal flow makes our approach related to copula-based distribution approximation and related approaches [28 29 577. However, by explicitly limiting the stimulus dependence to occur via the shift in the mean of the FA model along with flow-based transformation of responses, we avoid the reliance on the repeated presentations of the stimuli [29] or highly constrained forms of the marginal distribution (8)..
latents into a continuous approximation of Poisson
Limitations and future extensions As discussed above, our flow-based approach generalizes several existing methods to capture
variability of neural responses while being computationally more tractable. This allows to train our models end-to-end directly the likelihood via common gradient-based optimization algorithms. Within this general framework, we transformations,
presented specific case where we learned neuron-specific mapping responses into a FA model whose mean varies with the stimulus. As noted earlier, for each stimulus. this approach closely parallels Gaussian copula and thus shares much of the same limitations. Also, the fact that stimulus- deepedded changes in the covariance structure only occur through the learned transformation implies that the model can only capture changes in the covariance structure that varies with the mean (a limitation shared with many of the existing models). That being said, we believe that our general approach of flow-based modeling of neural response distributions allows for several generalizations that would overcome these limitations. Examples include an explicit dependence of the FA's covariance matrix on the stimulus, as well as the usage of richer, potentially stimulus-dependent, learnable transformations.
Broader impact Accurate models of neural variability such as the one presented here can lead to deeper scientific insights and understanding of how brains perceive and compute with sensory information, and can eventually also provide insights into how neurological and psychological disorders may disturb these functions, In particular, a more accurate model that relates internal brain states, stimulus-driven responses, and anatomical features such retinotopy or memberships to certain brain areas might provide deeper insights into the computational principles of cortex. Naturally, our model requires data from animal experiments to be trained. However, we used existing datasets with very general protocols that can be used in several analyses to make efficient scientific use of data from animal experiments. Furthermore, models such the one presented here do help to reduce the amount of animal experiments as faithful models allow us to explore the functional principles of neural populations in silico.
10

---

Acknowledgments and Disclosure of Funding
We thank all reviewers for their constructive and thoughtful feedback. Furthermore, we thank Alexander Ecker, Konstantin Willeke, Amne Nix, Christoph-Benjamin Blessing, and Emmanouil Froudarakis for comments and discussions. MB is supported by the International Max Planck Research School for Intelligent Systems. KKL is funded by the German Federal Ministry of Education and Research through the Tübingen AI Center (FKZ: 01IS18039A). FHS is supported by the Carl-Zeiss--tiftune and acknowledges the support of the DFG Cluster of Excellence "Machine Learning New Perspectives for Science", EXC 2064/1, project number 390727645. This work was supported by an AWS Machine Learning research award to FHS. Also supported by the Intelligence Advanced Research Projects Activity (IARPA) via Department of Interior/Interior Business Center (Dol/IBC) contract number DI6PC00003, The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. Disclaimer: The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of IARPA, DoI/IBC, the U.S. Government.
References
[1] F Dean. The variability of discharge of simple cells in the cat striate cortex. Exp. Brain Res., 44(4):437-444 440, 1981
[2] D J Tolhurst, J Movshon, and Dean. The statistical reliability of signals in single neurons in cat and monkey visual cortex. Vision Res., 23(8):775-785, 1983.
[3] George J Tomko and Donald R Crapper. Neuronal variability: non-stationary responses to identical visual stimuli, Brain research, 79(3):405-418, 1974.
[4] Michael Shadlen and William T Newsome. The variable discharge of cortical neurons: implications for connectivity, computation, and information coding. Journal of neuroscience, 18(10):3 1810033880-3 1998.
[5] David Tolhurst, Anthony Movshon, and Andrew Dean. The statistical reliability of signals in single neurons in cat and monkey visual cortex. Vision research, 23(8):775-78) 1983.
[6] Marlene R Cohen and Adam Kohn. Measuring and interpreting neuronal correlations. Nature neuroscience, 14(7):811, 2011.
[7] Krešimir Josie, Eric Shea-Brown, Brent Doiron, and Jaime de la Rocha. Stimulus-dependent correlations and population codes. Neural computation, 2740022774-2804, 2009.
[8] Adrián Ponce-Alvarez, Alexander Thiele, Thomas D Albright, Gene R Stoner, and Gustavo Deco Stimulus-dependent variability and noise correlations in cortical mt neurons. Proceedings of the National Academy of Sciences, 110333)13162-11117, 2013.
[9] Mihály Bányai, Andreea Lazar, Liane Klein, Johanna Klon-Lipok, Marcell Stippinger, Wolf Singer, and Gergö Orbán. Stimulus complexity shapes response correlations in primary visual cortex. Proceedings of the National Academy of Sciences, 116()):723-22332 2019.
[10] Marlene R Cohen and William T Newsome. Context-dependent changes in functional circuitry in visual area mt. Neuron, 60(1):162-173, 2008.
[11] Ralf M Haefner, Pietro Berkes, and József Fiser. Perceptual decision decision-making as probabilistic inference by neural sampling. Neuron, 90(3):649-666 2016
[12] Marlene Cohen and John HR Maunsell. Attention improves performance primarily by reducing interneuronal correlations. Nature newroscience, 12(12):1594, 2009.
[13] Jude F Mitchell, Kristy A Sundberg, and John H Reynolds. Spatial attention decorrelates intrinsic activity fluctuations in macaque area v4. Neuron, (3(6) 888, 2009.
[14] Farran Briggs, George Mangun, and W Martin Usrey. Attention enhances synaptic efficacy and the signal-to-noise ratio in neural circuits. Nature, 476--499 2013.
11

---

[15] Cristopher M Niell and Michael P Stryker. Modulation of visual responses by behavioral state mouse visual cortex. Neuron, 65(4):472-479, 2010.
[16] Alexander S Ecker, Philipp Berens, R James Cotton, Manivannan Subramaniyan, George H Denfield, Cathryn R Cadwell, Stelios M Smirnakis, Matthias Bethge, and Andreas S Tolias. State dependence of noise correlations in macaque primary visual cortex. Neuron, 82(1): 235-248 2014.
[17] David A Klindt, Alexander S Ecker, Thomas Euler, and Matthias Bethge. Neural system identification for large populations separating "what" and "where". Adv. Neural Inf. Process. Syst., November 2017.
[18) Eleanor Batty, Josh Merel, Nora Brackbill, Alexander Heitman, Alexander Sher, Alan Litke, J Chichilnisky, and Liam Paninski. Multilayer recurrent network models of primate retinal ganglion cell responses. November 2016.
[19] Lane T McIntosh, Niru Maheswaranathan, Aran Nayebi, Surya Ganguli, and Stephen A Baccus. Deep learning models of the retinal response to natural scenes. In Advances in neural information processing systems, volume 29, pages 1369-1377, February 2016.
[20] Fabian Sinz. Alexander S Ecker, Paul Fahey, Edgar Walker, Erick Cobos, Emmanouil Froudarakis, Dimitri Yatsenko, Zachary Pitkow, Jacob Reimer, and Andreas Tolias. Stim- ulus domain transfer in recurrent models for large scale cortical population prediction video. In Advances in Neural Information Processing Systems 31, pages 7199-7210, 2018.
[21] Santiago A Cadena, George H Denfield, Edgar Y Walker, Leon A Gatys, Andreas S Tolias, Matthias Bethge, and Alexander S Ecker. Deep convolutional models improve predictions of macaque v1 responses to natural images. PLoS computational biology,
2019.
[22] Konstantin- Konstantin-Klemens Lurz, Mohammad Bashiri, Konstantin Friedrich Willeke, Akshay Kumar Jagadish, Erie Wang, Edgar Y Walker, Santiago Cadena, Taliah Muhammad, Eric Cobos, Andreas Tolias, et al. Generalization in data-driven models of primary visual cortex. bioRxiv, 2020.
[23] Edgar Y Walker, Fabian H Sinz, Erick Cobos, Taliah Muhammad, Emmanouil Froudarakis, Paul G Fahey, Alexander S Ecker, Jacob Reimer, Xaq Pitkow, and Andreas S Tolias. Inception loops discover what excites neurons most using deep predictive models. Nature neuroscience, 22(12):000-2265, 2019.
[24] Pouya Bashivan, Kohitij Kar, and James J DiCarlo. Neural population control via deep image synthesis. Science, :33av4463)) May 2019.
[25] Adam S Charles, Mijung Park, J Patrick Weller, Gregory D Horwitz, and Jonathan W Pillow. Dethroming the fano factor: a flexible, model-based approach to partitioning neural variability. Neural computation, 30(4):1012-1045 2018.
[26] Robbe LT Goris, J Anthony Movshon, and Eero P Simoncelli. Partitioning neuronal variability. Nature neuroscience, 17(6):858-86) 2014.
[27] Cian O' Donnell, J Tiago Gonçalves, Nick Whiteley, Carlos Petter-Cailliia,, and Terrence J Sejnowski. The population tracking model: A simple, scalable statistical model for neural population data. Neural Comput., 29(1):50-93, January 2017.
[28] Pietro Berkes, Frank Wood, and Jonathan Pillow. Characterizing neural dependencies with copula models. lattps:/ppiloowal pdf Accessed: 2021-5-22.
NNPS
[29) Oleksandr Sorochynskyi. Stéphane Deny, Olivier Marre, and Ulisse Ferrari. Predicting syn- chronous firing of large neural populations from sequential recordings. PLoS Comput. Biol., 17 100se000 January 2021.
[30] Byron M Yu, John P Cunningham, Gopal Santhanam, Stephen I Ryu, Krishna V Shenoy, and Maneesh Sahani. Gaussian- factor analysis for low-dimensional single-trial analysis of neural population activity. J. Neurophysiol. 102(1):614-66 July 2009.
12

---

[31] Jakob H Macke, Lars Buesing, John P Cunningham, Byron M Yu, Krishna V Shenoy, and Maneesh Sahani. Empirical models of spiking in neural populations. In J Shawe-Taylos 7emel, P I. Bartlett, F Pereira, and K Q Weinberger, editors, Advances in Neural Information Processing Systems 24, pages 1350-1358. Curran Associates, Inc., 2011.
[32] Evan W Archer, Urs Koster, Jonathan W Pillow, and Jakob H Macke. of neural population activity in sensory cortical circuits. In Advances in Neural Information Processing Systems 27: 28th Conference on Neural Information Processing Systems (NIPS 2014), pages 343-351, 2015.
models
[33] Yuan Zhao and Il Memming Park. Variational latent gaussian process for recovering Single-Trial dynamics from population spike trains. Neural Comput., 1995):12931331, May 2017.
[34] Anqi Wu, Nicholas A Roy, Stephen Keeley, and Jonathan W Pillow. Gaussian process based nonlinear latent structure discovery in multivariate spike train data. Adv. Neural Inf. Process. Syst., 30:3496-3505, December 2017.
[35] Stephen L Keeley, Mikio C Aoi, Yiyi Yu, Spencer L Smith, and Jonathan W Pillow. Identifying signal and noise structure in neural population activity with gaussian process factor models. July 2020.
[36] E G Tabak. A family of non-parametric density estimation algorithms. https: h math nyu. httd:///aaabakkuullaaaaaatto///s Turner---uurrrepp pdf 2000. Accessed: 2021-5-25.
[37] Esteban G Tabak and Eric Vanden-Eijndee Density estimation by dual ascent of the log-
likelihood. CMS Books Math./Ouvrages Math. SMC, 8(1):217-233, March 2010.
[38] Oren Rippel and Ryan Prescott Adams. High-Dimensional probability estimation with deep density models. February 2013.
[39] J JP Agnelli, M Cadeiras, E G Tabak, C V Turner, and E Vanden-Eijnden. Clustering and classification through normalizing flows feature space. Multiscale Model. Simul., 8(5): 1784-1802, January 2010.
[40] Dinh, Sohl-Dickstein, and Bengio. Density estimation using real NVP. Technical report, 2017.
[41] Danilo Jimenez Rezende and Shakir Mohamed. Variational inference with normalizing flows. 1505. 050700 2015.
[42] Sam Roweis and Zoubin Ghahramani. A unifying review of linear gaussian models. Neural Comput., 11(2):305-345, 1999.
[43] Xue-Xin Wei, Ding Zhou, Andres Grosmark, Zaki Ajabi, Fraser Sparks, Pengcheng Zhou, Mark Brandon, Attila Losonezy, and Liam Paninski. A zero-inflated gamma model for deconvolved calcium imaging traces. June 2020.
[44] François Chollet. Xception:
Deep learning with depthwise separable convolutions. In Proceed- ings of the IEEE conference on computer vision and pattern recognition, pages 1251-1258, 2017.
[45] Shaul K Bar-Lev and Peter Enis. On the classical choice of variance stabilizing transformations and application for poisson variate. Biometrika, 75(4):803-804, 1988.
[46] Nicholas James Sofroniew, Daniel Flickinger, Jonathan King, and Karel Svoboda. A large field of view two-photon mesoscope with subcellular resolution for in vivo imaging. Elife, 5:e1 14424 2016.
[47] Olga Russakovsky, Jia Deng. Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma. Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, Alexander C. Berg, and Li Fei-Fei. ImageNet Large Scale Visual Recognition Challenge. International Journal of Computer Vision (IJCV), 115(3):111-25 252, 2015. doi: 10. 10007/s1163--15-0000
[48] Facebook. Adaptive experimentation platform, 2019. URL https: // //p..xxde/ dev/
13

---

[49] Marina E Garrett, Ian Naubaus, James H Marshel, and Edward M Callaway. Topography and areal organization of mouse visual cortex. Journal of Neuroscience, 34(37):12587-12661 2014.
[50] Jacob Reimer, Emmanouil Froudarakis, Cathryn R R Cadwell, Dimitri Yatsenko, George H H Denfield, and Andreas S S Tolias. Pupil fluctuations track fast switching of cortical states during quiet wakefulness. Neuron, (112))355-322) 2014.
[51] Martin Vinck, Renata Batista-Brito, Ulf Knoblich, and Jessica A Cardin. Arousal and locomotion make distinct contributions to cortical activity patterns and visual encoding. Neuron, 86(3): 740-754, May 2015.
[52] Matthew J McGinley, Martin Vinck, Jacob Reimer, Renata Batista-Brito, Edward Zagha, Cathryn R Cadwell, Andreas S Tolias, Jessica A Cardin, and David A McCormick. Waking state: Rapid variations modulate neural and behavioral responses. Neuron, 87(6):1143-1161 September 2015.
[53] Jacob Reimer, Matthew J McGinley, Yang Liu, Charles Rodenkirch, Qi Wang, David Mc- Cormick, and Andreas S Tolias Pupil fluctuations track rapid changes in adrenergic and cholinergic activity in cortex. Nat. Commun., 7:1 7:13289, November 2016.
[54] Siddhartha Joshi, Yin Li, Rishi M. Kalwani, and Joshua I. Gold, Relationships between Pupil Diameter and Neuronal Activity in the Locus Coeruleus, Colliculi, and Cingulate Cortex. Neuron, (9)))221-234, 2016. ISSN 0896-6273. doi: 00.0016f.neuron.211....00
[55] Student. Probable error of a correlation coefficient. Biometrika, pages 302-310, 1908.
[56] Carsen Stringer, Marius Pachitariu. Nicholas Steinmetz, Charu Bai Reddy, Matteo Carandini, and Kenneth D Harris. Spontancous behaviors drive multidimensional, brainwide activity. Science, 364(6437), 2019.
[57] Jakob H Macke, Philipp Berens, Alexander S Ecker, Andreas Tolias, and Matthias Bethge. Generating spike trains with specified correlation coefficients. Neural Comput., 21(2):397-423 February 2009.
[58] Emmanouil Froudarakis, Uri Cohen, Maria Diamantaki, Edgar Y Walker, Jacob Reimer, Philipp Berens, Haim Sompolinsky, and Andreas S Tolias. Object manifold geometry across the mouse cortical visual hierarchy. August 2020.
[59] Eftychios A Pnevmatikakis, Daniel Soudry, Yuanjun Gao, Timothy Machado, Josh Merel, David Pfau, Thomas Reardon, Yu Mu, Clay Lacefield, Weijian Yang, et al. Simultancous denoising, deconvolution, and demixing of calcium imaging data. Neuron, 89(2):285-299 2016.
[60] Kingma and Ba. Adam: A method for stochastic optimization. International Conference on Learning Representations, pages 1-13, 2014.
[61] Lutz Prechelt. Early stopping - but when? In Grégoire Montavon, Geneviève Orr, and Klaus-Robert Müller, editors, Neural Networks: Tricks of the Trade: Second Edition, pages 53-67. Springer Berlin Heidelberg, Berlin, Heidelberg, 2012.
[62] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang. Zachary Vito,o Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high- performance deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché- Buc, E. Fox, and Garnett, editors, Advances in Neural Information Processing Systems 32, pages 8024-8035. Curran Associates, Inc., 2019. URL http:: hww.:/ppppp......ppppppp.....ppp www.:new.rrrccccoop..eee cc/paprr/ -ibrary
9015- an hugpppacii...tttiiiii.rerrrrrrrr pdf
[63] Charles R. Harris, K. Jarrod Millman, Stéfan van der Walt, Ralf Gommers, Pauli Vir- tanen, David Counaaeaa,, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith,
14

---

Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Hal- dane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchand Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant. Array programming with NumPy. Nature, 58577825)3357-362 September 2020. doi: 10.1008//4115 020-6-020-24999
URLhttps //doi org/ 10. 1038/ #4 e415445844 586 020-26999000-2222
15


Reviews:
[{"review_id": "yehioafD92i", "reviewer": "Reviewer_BNcm", "summary": "The authors present a Factor Analysis + flow model for modeling stimulus dependent neural responses. The FA + flow architecture allows the model to capture stimulus dependence, trial to trial correlated variability, and some features of the conditional and marginal response distributions. Validation is performed on synthetic data and real data from mouse visual cortex using 2 photon calcium imaging. The networks and inferred latent states are then compared to anatomical information and behavioral states.", "questions": "", "limitations": "", "rating": 8, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is clearly written and the results and figures are clear. The FlowFA model is a concise way of including both stimulus dependent and stimulus conditioned (correlated) variability into a model. The exploration of potential uses of this type of methods in Fig 4 is also compelling (although see the technical concern). I also think this model bridges the gap between shallow (Factor Analysis) and deep (unconstrained flow models) in a way that provides both interpretability and more accurate modeling which could serve as useful techniques in future models.\n\nI have one technical concern about how some the latent state results are interpreted, some other clarifying questions/comments about the method and results, and a suggestion for the discussion.\n\n- I think it would benefit a broad neuroscience audience to spend a little more time describing flow-based models in Section 2.1. Contrasting them with standard feed forward models which are probably more familiar might help, i.e., what mappings can you do with flow based models (like inverse) that is not available with regular models. You mention exact likelihood computation and sampling, but expanding on this would help.\n- Providing a bit more information about the neural data in Section 2.4 would help. How is the response variable being preprocessed? I would assume is it florescence $\\Delta F / F$, but this isn't specified. Do you truncate it to $>=0$?\n- In Fig 3b, what latent dimensionality was used? Do you get qualitatively similar transformations as a function of latent dimensionality?\n- Fig 4a, in the Cortex Position panels, it would be helpful to explicitly state what color indicates in the Cortex positions and RF positions panels in the figure caption. I'm assuming experimentally measured area assignment.\n- For 4a, does the inclusion of correlated variability in the model improve the inferred estimates of $\\delta$ compared to a model fit with $C=0$, i.e., don't take correlations into account but share parameters in the same way? This analysis isn't strictly necessary to perform, but it would strengthen the result.\n\nTechnical concern (and my main reservation for a higher rating)\n- For Fig 4 b and c, $C$ is only determined up to a rotation in the latent space. Both the likelihood and equation for $\\mathbf v$ are unchanged by taking $C\\rightarrow CR$ for some rotation $R$. $CRR^TC^T=CC^T$ for the likelihood and the distribution of $\\mathbf v$ is unchanged since $\\mathbf z$ is isotropic. I believe this complicates the interpretation of Fig 4 b and c since it is not clear how sign relationships will be preserved under arbitrary rotations. There are schemes for defining a particular rotation [1] or a different measure of $C$ which is invariant to rotations could be shown, such as the leverage scores (see [2] for a definition).\n- Due to the above point, it might make more sense to report the $R^2$ from regressing pupil dilation onto all latent states since the particular latent state basis is not meaningful.\n\n1. Kaiser, H. F. (1959). Computer program for varimax rotation in factor analysis. Educational and psychological measurement, 19(3), 413-420.\n2. Mahoney, Michael W., and Petros Drineas. \"CUR matrix decompositions for improved data analysis.\" Proceedings of the National Academy of Sciences 106.3 (2009): 697-702.\n\n- For the results in Fig 4d, specify how $p(z|r)$ is computed in the FlowFA model (maybe earlier in section 2.1). Line 258 implies a MAP estimate is used, but it is not clear.\n- For the sentence on lines 127-129, are only the nonlinearities shared but affine parameters different per-neuron? or all all nonlinearities and affine parameters shared? If the later, did you consider a model where each neuron has an independent set of parameters? Was overfitting a problem? Do you get roughly the same transformation for every neuron anyway?\n\n\nThis is somewhat tangential to the main thrust of the paper, but the sentence that starts on line 297 in the Discussion contrasts healthy, functional brains and brains with neurological and psychological disorders. People with neurological or psychological disorders can have healthy, functional brains. For example, rates of depression in graduate students in the US have greatly increased during the pandemic, but the cause probably won't be found in the brains of the grad students.\nhttps://www.nature.com/articles/d41586-020-02439-6\n\nI would rephrase this sentence to lessen the contrast. Suggestion:\n\"Accurate models of neural variability such as the one presented here can lead to deeper scientific insights and understanding of how brains perceive and compute with sensory information, and can eventually also provide insights into neurological and psychological disorders.\"\n\nMinor comments\n- and this cannot -> and thus cannot", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors provide some discussion of the limitations. I think the limitations of diffeomorphisms could be more clear. Is another way of stating it that the diffeomorphisms cannot map a distribution with 1 peak (Gaussian) to a distribution with potentially multiple peaks (neural response)?", "review_text": "", "extra_scores": {"ethical_concerns": "The authors discuss some ethical concerns and contributions of improved modeling (including this method) to the use of animal models.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "uU38WFnQQEC", "reviewer": "Reviewer_Two8", "summary": "The authors propose an algorithm based on a combination of a DNN-based model (that account for stimulus-variability) and a flow-based model (accounting for inherent neural variability) to model the behavior of neural responses. In this paper, the algorithm, called FA-flow (compared with ZIG, Poisson and ZA-Fixed) is trained on pairs of image/neural recording from bi-photon microscope on 2 mices. The authors claim is that their algorithm perform similarly to current state of the model neural response while additionally yielding latent states that exhibit meaningful relations to anatomical and functional properties.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "## General comments :\nThis article is in general interesting, and bring new Machine Learning concepts (i.e. the normalizing Flow) into computational Neuroscience. The main claim of the article, that is adding a flow on top of a DNN-based model allow to model both types of neural variability (i.e. stimuli variability and inherent neural variability) is not properly demonstrated. Nevertheless, the experiments conducted by the authors demonstrate the FA-Flow algorithm produce a relatively good fit of neural responses, but some of the conclusion made by the authors are not support by quantitative results (i.e. RFs mappings) or suffer from low statistical significance (see below for more details). In addition, some of the conclusion made by the authors are complicated to grasp (especially the conclusion from the Fig 4b and c).  My score : 5. \n\n## More specifics comments :\nSection Methods : This section is clear and well written. Minor question : The authors are using 1000 V1 neurons to train the model on scan 1, and they report a total number of V1 neurons of 2867. For scan 2 the selection process is described (this is random) but  not for scan 1. Could please describe the neuron’s selection process for V1 neurons on scan 1 ?\n\n### Section Results : \n#### Section 3.1 : \nI think a crucial experiment is missing in this section. One of the main claim of the authors is that flexible latent model better capture intrinsic correlated neural activity, and this claim is not demonstrated for their FA-Flow model. On my understanding the authors have all the datas to conduct such an experiment : some of the images are repeated 10 times (see line 172), thus one can see on these specific images how the FA-Flow model is fitting the intrinsic neural variability compared to the other models. Such an experiments would support the main claim of the article.\n\nFigure 3 : Line 236, the authors state : « with Increasing latent dimension the FlowFA model … approaches the ZIG model in terms of likelihood’.  I do not see this improvement of log-likelihood with the increase of latent dimensions. As this is far from being obvious, the authors should conduct statistical test to support their claim.\n\nFigure 4a : We see qualitatively that the FAFlow algorithm allows deducing a relatively good mapping between the RF position in the image space and neurons position in the cortex. But without further quantification of the quality of the mapping, it is hard to draw a solid conclusion. Is it possible to quantify the quality of the mapping ?\n\nFigure 4b and c : I have trouble understanding the claim made by the authors : ’the effect of some latent states vary systematically accros brain area, as a function of RF location' (line 255). Could you please further explain how you are drawing this conclusion. Also a legend showing the meaning of the colors would help to better understand these 2 subfigures.\n\nFigure 4d : Without further comparison it is complicated to say that ’some of the inferred latent state correlated strongly with pupil dilatation' (see line 259). For example in scan 1, the latent 3 accounts for only 41% of the variability (r2=0.41) which is usually not high enough to consider a good fit. Further statistical tests need to be conducted to be able to conclude on the significance of the prediction (and ideally a comparison with other r2 given in the literature).\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors have not addressed the societal impact of their work", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4h"}, "extra_text": {}}, {"review_id": "Wa-GHOEB_Z", "reviewer": "Reviewer_mWE7", "summary": "In the paper \"A Flow-based latent state generative model of neural population responses to natural images\", the authors propose a model that can take into account stimulus driven stimulus-conditioned neural activity. The model is based on normalizing flows for the marginals and on a factor analysis model for modelling covariances and is validated on simulated data as well as data recorded from mouse visual cortex.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Originality:\nTo the best of my knowledge, the factor analysis approach combined with normalizing flows for the marginals is novel. That said, many non-cited models have previously been proposed, explicitly taking stimulus driven and stimulus-conditioned neural activity into account. For instance, Keeley, Aoi, Yu, Smith and Pillow, bioRxiv 2020 explicitly model signal and noise components in a Gaussian process factor analysis framework. Sorochynskyi, Deny, Marre and Ferrari, PLoS Comput Biol 2021 used a maximum entropy model and explicitly take noise correlations into account as well.\n\nQuality:\nThe flow transformations are compared against two fixed transformations: square-root and Anscombe. The alternatives being fixed, it is not surprising that the flow model outperforms these fixed transformations. A comparison to a more flexible transformation such as a kernel density estimator would have been more insightful. Figure 5 suggests that even the ZIG model does a better job of modelling the marginals of the analyzed dataset.\nSimilarly, the comparisons to the control models are not fair. Contrary to the proposed model, both the independent Poisson and the ZIG model assume independence (line 137).\nTractable maximum entropy models (e.g. O'Donnel et al. Neural Comput 2017) or the dichotomized Gaussian model with appropriate marginals (Macke et al., Neural Comput 2009) would have been better contenders.\n\nClarity:\nThe clarity of the manuscript is fine but could be improved. In line 87, we learn that neural responses are represented as real numbers but it is not clear what these numbers mean until much later in the manuscript. Sometimes the model is called FAFlow (e.g. Fig 5), other times FlowFA (e.g. line 88).\n\nSignificance:\nModels for taking both stimulus driven and stimulus-conditioned activity into account are important tools in systems neuroscience. Given the lacking comparisons, the value of the proposed model is unclear.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Limitations of the work with respect to the ZIP marginals are discussed to some extend in the Discussion but important points are left out. The FlowFA model assumes a stimulus-independent Gaussian dependence structure which is transformed by marginal stimulus-dependent nonlinearities. How flexible is the resulting dependence structure? Is the flow transformation of line 129 limiting? Such questions are not addressed in the manuscript. Societal implications with respect to animal experiments are briefly mentioned in the Discussion.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "8vvraOFtrJB", "reviewer": "Reviewer_QQF7", "summary": "This paper proposes a model of neural responses by capturing both stimulus dependent mean as well as variance. They transform the firing rate of each neuron using a reversible flow based model and model the transformed responses as a multi-variate gaussian. The mean of the multivariate gaussian is given by a neural network transformation of the input stimulus, and the low rank + diagonal covariance is independent of stimulus. They show the utility of this approach on synthetic as well as calcium imaging recordings from mouse visual cortex.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is well written, performs relevant baseline comparisons and shows the utility on real neuroscience data. \n\n* I think the overall idea is quite novel. It might be helpful to describe the approach as learning one transformation of stimulus and another transformation of responses, and mapping them both to a common space where geometry is simpler. This general idea is quite different than most models in computational neuroscience that either learn encoding (stimulus-> response), or decoding (response-> stimulus) models. \n\n* My high-level comment is that the ideas can be simplified by removing the connection to flow models at many places. The paper uses (1) only per-neuron invertible transformations and (2) uses same transformation for each neuron, thereby not using the full power of current flow-based models.\n\n* In terms of an ideal model, I think the main limitation is that the correlated noise in responses does not have a clear dependence on stimulus. That is, the degree of covariance between responses could be high for some stimuli and low for others in neuroscience data, but the model does not capture that. This limitation should be highlighted better in the text. \n\n* It seems that the transformation learned by the flow in Figure 3b is highly nonlinear and step-like (pink line). Would it be totally fine to lose the intertibility property of flow models then? All components of the model are not necessary for this real data case.\n\n* Similarly, the ZIG model outperforms FAFLOW model for zero latents. Losing the invertibility of FLOW model for a transformation that looks more like ZIG (Figure 5) might be better overall (even for >0 factors).\n\n* The conditional correlation plots in Figure 3a are not explained in the main text, and they seem to be the strongest results. \n* Moreover, why does flow model perform better than ZIG for 0 latents (where cells are modeled independently) for conditional correlation metric, but not for likelihood metric? \n\n* The interpretation of latents (Figure 4) could be cofounded by experimental factors. The spatial smoothness of latents could be explained by spatially smooth expression levels of calcium indicator, or spatially smooth recording noise. Either control for these confounds, or explain them as possible interpretation.\n\n* The det(Jacobian) explanation in Figure 4 and text is not clear. Please explain them in non-mathematical terms to highlight the relevance/importance. \n\n* Color bars missing in Figure 4", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "no societal impact.", "review_text": "", "extra_scores": {"ethical_concerns": "no ethical cconcern.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6 hours"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper introduces a system identification model for stimulus-driven and stimulus-conditioned fluctuations. They propose an architecture based on factor analysis and a flow model, and validate their performance on synthetic and mouse visual calcium imaging data. The paper addresses an important question in neuroscience, it is well written, introduces a new ML approach to computational neuroscience, and the authors have done a great deal to address all the concerns of the reviewers and convince them of the points they intended to make. Given this, I recommend this paper for acceptance and ask the authors to make all the modifications they promised in the final manuscript. 

Decision: Accept (Spotlight)