Trial matching: capturing variability with data-constrained spiking neural networks
Authors: Christos Sourmpis
Carl C. H. Petersen
Wulfram Gerstner
Guillaume Bellec
Venue: NeurIPS 2023 poster (2023)
Keywords: neuroscience
spiking networks
data-constrained modeling
electrophysiological recordings
optimal transport
trial variability
RNN
interpretable machine learning
TL;DR: We fit the parameters of a large recurrent spiking model to explain electrophysiological recordings. To capture the trial-to-trial variability in the data, a loss function based on optimal transport is defined.

Abstract: Simultaneous behavioral and electrophysiological recordings call for new methods to reveal the interactions between neural activity and behavior. A milestone would be an interpretable model of the co-variability of spiking activity and behavior across trials. Here, we model a mouse cortical sensory-motor pathway in a tactile detection task reported by licking with a large recurrent spiking neural network (RSNN), fitted to the recordings via gradient-based optimization. We focus specifically on the difficulty to match the trial-to-trial variability in the data. Our solution relies on optimal transport to define a distance between the distributions of generated and recorded trials. The technique is applied to artificial data and neural recordings covering six cortical areas. We find that the resulting RSNN can generate realistic cortical activity and predict jaw movements across the main modes of trial-to-trial variability. Our analysis also identifies an unexpected mode of variability in the data corresponding to task-irrelevant movements of the mouse.

Trial matching: capturing variability with data-constrained spiking neural networks
Christos Sourmpis, Carl C.H. Petersen, Wulfram Gerstner, Guillaume Bellec Brain Mind Institute,
School of Computer and Communication Sciences and School of Life Sciences, École Polytechnique Fédérale de Lausanne (EPFL). CH-1015 Lausanne, Switzerland, (f (firstname. (esstnamm@@ppf....
Abstract
Simultaneous behavioral and electrophysiological recordings call for new meth- ods to reveal the interactions between neural activity and behavior. A milestone would be an interpretable model of the co-variability of spiking activity and be- havior across trials. Here, we model a mouse cortical sensory- -otoor pathway in tactile detection task reported by licking with a large recurrent spiking neural network (RSNN), fitted to the recordings via madieent-based optimization. We fo- cus specifically on the difficulty to match the trial-to-trial variability the data. Our solution relies on optimal transport to define a distance between the distri- butions of generated and recorded trials. The technique applied to artificial data and neural recordings covering six cortical areas. We find that the resulting RSNN can generate realistic cortical activity and predict jaw movements across the main modes of trial-to-trial to-trial variability. Our analysis also identifies an unex- pected mode variability in the data corresponding to of the mouse.
movements
Introduction
Over the past decades, there has been remarkable advancement in neural recording technology. Today, we can simultaneously record hundreds, even thousands, of neurons with millisecond time precision. Coupled with behavior measurements, modern experiments enable stand how brain activity and behavior are intertwined [T] In these experiments, it is often observed that even well-trained animals respond the same stimuli with considerable variability. For exam- ple, mice trained on a simple tactile detection task occasionally miss the water reward (2), possibly because of satiation, lack of attention or neural noise. is also clear that there is additional uncon- trolled variability in the recorded neural activity 3 5 induced for instance by wide range of task-iree iirelvrreleaatt movements.
better under-
Our goal is to reconstruct a simulation of the sensory-motor circuitry driving the variability of neural activity and behavior. To understand the generated activity at circuit level, develop gener- ative model which biologically interpretable: all the spikes are generated by a recurrent spiking neural network (RSNN) with hard-biological constraints (i.e. the voltage and spiking dynamics are simulated with millisecond precision, neurons are either inhibitory or excitatory, spike transmission delay takes 4 ms]
First contribution, we make a significant advance in the simulation methods for data-constrained RSNNs. While most prior works
were limited to single recording sessions, model is
Our code is available at aft.hblcawnPEE....loppoaompppppp (code DOI hosted by zenodo
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

constrained to spike recordings from 28 sessions covering six cortical areas. The resulting spike- based model enables a data-constrained simulation of cortical sensory-motor pathway (from so- matosensory to motor cortices responsible for the whisker, jaw and tongue movements). As far as we know, our model is the first RSNN model constrained to multi-session recordings with automatic differentiation methods for spiking neural networks
TO
Second contribution, using this model we aim to pinpoint the circuitry that induces variability in be- havior (asking for instance what circuit triggers a loss of attention). Towards this goal, we identify an unsolved problem: "how do we enforce the generation of a realistic distribution of neural activity and behavior?" To do this, the model is fitted jointly to the recordings of spiking activity and move- ments to generate a realistic trial-to-trial co-variability between them. Our technical innovation is to define a supervised learning loss function to match the recorded and generated variability. Con- cretely the trial matching loss function is the distance between modeled and recorded distributions of neural activity and movements. It relies on recent advances in the field of optimal transport (T]] providing notions of distances between distributions. matching enables the recovery of the main modes of trial-to-trial variability which includes the neu- ral activity related to instructed behavior (e.g. miss versus hit trials) and uninstructed behavior like spontaneous movements.
our data-constrained RSNN, trial
Related work While there is a long tradition of data fitting using the leaky integrate and fire (LIF) model, spike response models [14] or generalized linear models (GLM) (6), most of these models were used to simulate single neuron dynamics [15] [6] or small networks with dozens of neurons recorded in the retina and other brain areas 6 7 17 A major drawback of those fitting algorithms
was the limitation to a single recording session. Beyond this, researchers have shown that FORCE
methods (18) could be used to fit up to 13 sessions with a large RSNN [17 19 20]. But in contrast
with back-propase proppaarpaattion through time (BPTT) in RSNNS 9 IO 2T.. FORCE is tied to the theory of recursive least squares making it harder to combine with deep learning technology or arbitrary loss functions, We know only one other study where BPTT is used to constrain RSNN to spike-train
recordings but this study was limited to single recording session.
Regarding generative models capturing trial-to-trial variability in neural data, many methods rely on trial-specific latent variables (22 23 24 25 26). This is often formalized by abstracting away the physical interpretation of these latent variables using deep neural networks (e.g. see LFADS [22
or spike-GAN (27 but our goal is here to model the interpretable mechanisms that can generate the recorded data. There are hypothetical implementations of latent variables in RSNNs, most com- monly, latent variables can be represented as the activity of mesoscopic populations of neurons (25]
or linear combinations of the neuronal activity [26 28 29). These two models assume respectively an implicit grouping of the neurons [25
low-rank connectivity matrix (26] 28 29] Here, we want to avoid making any structural hypothesis of this type a priori. We assume instead that the variability is sourced by unstructured noise (Gaussian current or Poisson inputs) and optimize the network parameters to transform it into a structured trial-to-trial variability (c.g. multi-modal distri- bution of hit versus miss trials). The optimization therefore decides what is the network mechanism that best explains the trial-to-trial variability observed in the data.
This hypothesis-free approach is made possible by the trial matching method presented here. This method is complementary to previous optimization methods for generative models in neuroscience. Many studies targeted solely trial-averaged statistics and ignored single-trial activity, for instance methods using the FORCE algorithm [17] 19
20 30 3TJ. RSNN methods using back-propagation
through time [8] and multiple techniques using (non-interpretable) deep generative models (32). There exist other objective functions which can constrain the trial-to-trial variability in the data,
namely: the maximum likelihood principle 6 15 spike-GANs 27 33 We illustrate however in
the discussion section why these two alternatives are not a straightiorward replacement for the trial matching loss function with our interpretable RSNN generator.
2 Large data-constrained Recurrent Spiking Neural Network (RSNN)
This paper aims to model the large-scale clectrophysiology recordings from (2 where they recorded 6.182 units from 12 areas across 18 mice whisker tactile detection task described in Figure I in 50% of the trials (the GO trials), a whisker is
P All animals in this dataset were trained to perform the
of which use 3,810 units from 6 cortical areas in this study

---

A
Data recording
B
Data-constrained RSNN
Um
ALM
Whisker
wM2
stimulus
ijM1
WM1
Auditory cue
WS1
<
wS2
Jaw movement
41)
(at)
Excitatory
Inhibitory
C
41)
33%
Data Model
False alarm
10% (FA)
Go
Lick
Hit
?%
No-Go
Lick
?%
-/
0
No Lick
Miss
17%
?%
-
0
No Lick
Correct rejection
(CR)
40% ?%
Time (s)
Time (s)
Figure 1: Modeling trial-varaability in electrophysiological recordings. A. During delayed whisker detection task, the mouse should report the sensation of a whisker stimulation by licking to obtain a water reward. Neural activity and behavior of the mouse are recorded simultancously. B. A recurrent spiking neural network (RSNN) of the sensorimotor pathway receives synaptic input modeling the sensory stimulation and produces the jaw movement as behavioral output. C. The stimuli and the licking action of the mouse organize the trials into four types (hit, miss, false alarm, and correct rejection). Our goal is to build model with realistic neural and behavioral variability. Panels A and are adapted from
deflected and after a 1 delay period an auditory cue indicates water availability if the mouse licks, whereas in the other 50% of trials (the No-Go trials), there is no whisker deflection and licking after the auditory cue is not rewarded. Throughout the paper we attempt to create data-constrained model of the six areas that we considered play major role in this behavioral task: the primary and secondary whisker somatosensory cortices (wS1, wS2), whisker motor cortices (wM1, wM2), the primary tongue-jaw motor cortex (tjMI) and the anterior lateral motor cortex (ALM), also known as ijM2 (see Figure MA and A).. While we focus on this dataset, the method described below aims to be broadly applicable to most contemporary large-scale electrophysiological recordings.
We built a spiking data-constrained model that simulates explicitly a cortical neural network at multiple scales. At the single-cell level, each neuron is either excitatory or inhibitory (the output weights have only positive or negative signs respectively), follows leaky-integrate and fire (LIF) dynamics, and transmits information in the form of spikes with synaptic delays ranging from 2 to 4 ms. At a cortical level, we model six brain areas of the sensory- consists of 250 recurrently connected neurons (200 excitatory and 50 inhibitory) shown in Figure BA, such that only excitatory neurons project to other areas. Since the jaw movement defines the behavioral output in this task, we also model how the tongue-jaw motor cortices (tjM1, ALM) drive the jaw movements.
pathway where each area
Mathematically, we model the spikes z/,, of the neuron j at time in the trial ber. The spiking dynamics are then driven by the integration of the somatic currents I into the membrane voltage v/J,- by integrating LIF dynamics with a discrete time step & = 2 ms. The jaw movement v is simulated with leaky integrator driven by the activity of tjM1 and ALM neurons, followed by an exponential non-linearity This can be summarized with the following equations, the trial index k is omitted for simplicity:
binary num-
of
ajut + (1-11,))) (1-a,))1 U +o-+5
(1)
15
Σ d.
Wt t-d
2
+
Σ
W x
(2)
d,
3

---

8' ajauD' + (1 ajow) Σ Wijnn 2/
(3)
y' = exp(i") +
(4)
where Wi. (jr W Wajaw, Uthraj and b are model parameters. The membrane time constants Tm
30 ms for excitatory and Tm = 10 ms for inhibitory neurons define aj = exp 61 and Tjaw
Tmm.f
50 ms define similarly Orjam which controls the velocity of integration of the membrane voltage and the jaw movement. To implement a soft threshold crossing condition, the spikes inside the
recurrent network are sampled with Bernoulli distribution zj Bio (" (((((")))))))))) )), where 20 is the temperature of the sigmoid (o). The spike trains x model the thalamic inputs simple Poisson neurons producing spikes randomly with firing probability of 5 Hz and increasing their firing rate when a whisker or auditory stimulation is present (see Appendix A). The last noise source Ej is an instantaneous Gaussian noise 5 of standard deviation Buthr V öt modeling random inputs from other areas (3 is model parameter that is kept constant over time).
Session stitching An important aspect of our fitting method is to leverage a dataset of electro- physiological recordings with many sessions. To constrain the neurons in the model to the data, we uniquely assign each neuron in the model to a single neuron from the recordings as illustrated in Figure A and AA. There are 3,810 recorded neurons in the dataset but we take a subset of 250 neurons per area to keep the same number of excitatory and inhibitory neurons in all areas (we were limited by the 290 regular spiking neurons recorded in wS1, and computation resources, mainly, the RAM of the GPU), we ignore the remaining data when fitting the model. This bijective mapping between neurons in the data and the model is fixed throughout the analysis and defines the area and cell types of the neurons in the model. The area is inferred from the location of the corresponding neuron in the dataset and the cell type is inferred from the action potential waveform of the cell (for simplicity, fast-spiking neurons are considered to be inhibitory and regular-spiking neurons as excitatory). Given this assignment, we denote =] as the spike train of neuron j in the dataset and 2j as the spike train of the corresponding neuron in the model; in general, an upper script DD always refers to the recorded data. With this notation the tensor of recorded data x)P only exists if neuron is recorded during session S (denoted j E S), and the activity of neurons not belonging in session S (j & S) is considered to be unknown. Consequently, two neurons i and j might be synaptically connected in the model although they correspond to neurons recorded in separate sessions. This choice is intended model network sizes beyond what can be recorded during single session. Our network is therefore a "collage" of multiple sessions stitched together as illustrated in Figure A and BA. This network is then constrained to the recorded data by optimizing the parameters to minimize the loss functions defined in the following section. Altogether, when modeling the dataset from Esmaeili and colleagues (2). the network consists of 1,500 neurons where each neuron is assigned to one neuron recorded in one of the 28 different recording sessions. Since multiple sessions are typically coming from different animals, we model "template mouse brain" which is not meant to reflect subject-to-subjec differences. In our view, this is coherent with the classical assumption that is made when recording this type of dataset: the features of the circuit being studied are shared among subjects that participate in the experiment. To study differences between sub- sub-groups of sub- jects (c.g. age group) one could separate the data into two session groups and compare two fitted models. This is beyond the scope of this paper and below we fit single model to the entire dataset.
3 Fitting single-trial variability with the trial matching loss function
We fit the network to the recordings with gradient descent and we rely on surrogate gradients to ex- tend propaggation to RSNNs 1o). At each iteration until convergence, we simulate batch of K = 150 statistically independent tralls. We measure some trial trial-average and single-trial statistics of the simulated and recorded activity, calculate a loss function, and minimize it with respect to all the trainable parameters of the model via gradient descent and automatic differentiation. This protocol sometimes referred to as a sample-and tion in GLMS where the network trajectory is clamped to the recorded data during optimization ( For the sound evaluation of the model, we separate the recorded trials of every session into a training set (75% of the total trials) and testing set (25%), that is never used in the fitting procedure. The
method [ as opposed to the likelihood optimiza-

---

split is done in a stratified manner so both sets have the same distribution of trial types. The full optimization lasts for approximately one to three days on GPU A100-SXM4440GB.
Trial-average loss We consider the trial-averaged activity over time of each neuron from every session Tneuron.j. sometimes referred also as a neuron's 's peristimulus time histogram (PSTH). This defined by Tneuron. (z) = * >x J, f where f is a rolling average filter with a window of 12 ms, and K is the number of trials in a batch of spike trains z. The statistics T (2P) are computed similarly on the K trials recorded during the session corresponding to neuron j. We denote the statistics T meuroon, after normalizing each neuron's 's trial-averaged activity, and we define the trial-averaged loss function as follows:
Lneuron IlTTneuron -Tneurnn (=P)
(5)
It is expected from 8 that minimizing this loss function alone generates realistic trial-averaged statistics like the average neuron firing rate.
Trial matching loss: fitting trial-to-trial variability Going beyond trial-averaged statistics, we now describe the trial matching loss function capture the main modes of trial-specific activity. From the previous neuroscience study (2). it appears that population activity in well-chosen areas is characteristic of the trial-specific variability. For instance, intense jaw movements are preceded by increased activity in the tongue-jaw motor cortices, and hit trials are characterized by secondary transient appearing in the sensory cortices a hundred milliseconds after a whisker stimulation. To de- fine statistics for a single trial & which can capture these features we denote the appulation-averaeen firing rate of an area A as (z)
JAnS| SrEAns(2j,A * f) where |A S| the number of neurons in area A recorded in session S. the smoothing filter f has a window size of 48 ms and the resulting signal is downsampled to avoid unnecessary redundancy. We write TA when each time bin is normalized to mean 0 and standard deviation 1 using the training data from session S. We, therefore, use 7 (z) as a feature vector to characterize the recorded area A during each trial & of session S, and we build richer feature vector including all recorded areas and the behavior with the concatenation T = (TAS, T A2 y/ * f). where Al and A2 are the two areas recorded in session S. and is the generated jaw trace from trial k. The same feature vector can be computed with the simulated jaw and simulated neural activity by sub-selecting the neurons from areas Al and A2 from session S to compare generated and recorded data.
1
T The challenging part is now define the distance between the recorded statistics trial, Px and the generated ones Terial. (z). Common choices of distances like the mean square error are not appro- priate to compare distributions. This is because the order of trials in batch of generated/recorded trials has no meaning priori: there is no reason for the random noise of the first generated trial to correspond to the first recorded trial - rather we want to compare unordered sets of trials and penalize if any generated trial is very far from any recorded trial.
Formalizing this mathematically we consider distance between distributions inspired by the opti- mal transport literature. Since the plain mean-squared error cannot be used, we use the mean-squared error of the optimal assignment between pairs of recorded and generated trials: we select randomly K' = min (, KP) generated and recorded trials (K and KP are respectively the number of gen- erated and recorded trials in one session), and this optimal assignment is formalized by the integer permutation x : {1....K' K'} + (1.... K"). Then using the feature vector T for any trial k we define the hard trial matching loss function as follows by the sum Licial = -ss CS with:
L = min Σ T (z) - - (R) (xP) 112
=min)
(6)
We compute this loss function identically on all the recorded sessions and take the summed gra- dients to update the parameters. Each evaluation of this loss function involves the computation of the optimal trial assignment æ which can be computed with the Hungarian algorithm (35) (see linear secear-summase sumaasignnment for implementation in scipy). This is not the only way to define distance between distributions statistics T In fact. this choice poses a potential problem because the optimization over æ a discrete optimization problem, so we have to assume that æ is constant with respect to the parameters when computing the loss gradients. We also tested alternative
5

---

A
B
Artificial Data
Without Trial matching
With hard Trial matching
With soft Trial matching
Artificial Recordings
A2
Miss like
A1
Firing rate (H2)
20
Session 1, 2,
Model Schematic
A1
A2
www
A
A2
A
Without Trial matching
With herd Trial matching
With soft Trial matching
Artificial Data
Miss like
Hit like
Excitatory
Trial Count
Inhibitory
50
All-to-all
25
E-E & E-I
25
25
90
25 50
25
Max firing rate of population average A2 (Hz)
Figure 2: Artificial Dataset. A. Session stitching: every neuron from the recordings is uniquely mapped to a neuron from our model. For example, an excitatory neuron from our model that be- longs in the putative Al is mapped constrain the connectivity so that only excitatory neurons can project across different brain regions. hit-like and a miss-like trial, while the second area (A2)
excitatory neuron "recorded" in A1. In network,
B. The first area (A1) responds equally responds only in hit-like trials. A model that does not use trial matching cannot capture the bimodal distribution of A2. C. Distribution of max firing rate of the population average of A2 from each trial. Only the trial matching algorithms retrieve the bimodal behavior of A2.
choices relying on relaxation of the hard assignment into a smooth and differentiable bi-stochastic matrix. This results in the soft trial matching loss function, which replaces the optimization over i (see the geomloss package for implementation in pytorch (13) In practice, to minimize both Ltrial (either the soft or hard version) and Lueuron simultaneously we optimize them in an additive fashion with a varameter-free multi-task method from deep learning which weightss the two loss functions to ensure that their gradients have comparable scales (see 36 for similar implementation).
by the Sinkhorn divergence 12 13
Simulation results
Validation using an artificial dataset We generated an artificial dataset with two distinct areas with 250 neurons each to showcase the effect of trial variability. In this dataset A1 (representing sensory area) always responds to a stimulus while A2 (representing a motor area) responds to the stimulus only 80% of the trials, see Figure This is a toy representation of the variability that is observed in the real data recorded in mice. so we construct the artificial data so that recording resembles a hit trial ("hit-like"") if the transient activity in A2 higher than 30 Hz (otherwise it's "miss-like" trial). From the results of our simulations Figure B-C, we can observe that the models that use trial matching (either soft trial matching or hard trial matching) can re- re-generate faithfully this bimodal response distribution ("hit-like" and "miss-like") in A2. In this dataset we saw little difference between the solutions soft and hard trial matching, any, soft trial matching reached its optimal performance with fewer iterations (see Figure A.2 As expected, when the model is only trained minimize the neuron loss for trial-averaged statistics, cannot generate stochastically this bimodal distribution and consistently generates instead noisy average response.

---

A
Model Schematic
ALM
WM2
wS1
wM1
ALM
Jaw movement
Session 1
tjM1
Session 2
wS2
wM1
Session
wS2
wM2
tjM1
Session 4
wS1
: Session 28
Excitatory Inhibitory
All-to-all E-E E-1
B 150
Tneuron. ((2)
Triat (2)
Data Model
Pearson correlation with testing set
1.0
100
Count
0.5
0.0
0
0.0
Firing Rate (Hz)
10'
10'
Data training
Model with triel matching
Modei without trail
training
Data
Model with trial
Model without trial
Figure 3: Large-scale electrophysiology recordings dataset. A. Session stitching: Every neuron from the recordings across sessions is uniquely mapped to neuron from model. For example, excitatory neuron from our model that belongs in the putative area tjMI is mapped to an excitatory neuron recorded from tjM1. In the pink box are the areas from which we decode the jaw trace. B. Baseline firing rate histogram, 200 ms before the whisker stimulus, from each neuron of our model and the recordings. Left: Pearson correlation of the PSTH, the violin plots represent the Pearson correlations across neurons. Right: trial-matched Pearson correlation of T the violin plots represent the distribution over 200 generated and recorded trial pairs.
Delayed whisker tactile detection dataset We then apply our modeling approach to the real large- scale electrophysiology recordings from 2 After optimization, we verify quantitatively that our model generates activity that is similar to the recordings in terms of trial-averaged statistics.
First, that the 1,500 neurons in the network exhibit realistic diversity of averaged firing rates: the distribution of neuron firing rates is log-normal and matches closely the distribution extracted from the data in Figure 3. Second, the single-neuron PSTHs of model the PSTHs from the recordings, This can be quantified by the Pearson trial-averaged correlation between the generated and held-out trials which we did not use for parameter fitting. We obtain an averaged Pearson correlation of 0.30 which is very close to the Pearson correlation obtained when comparing the training and testing sets 0.31 + 0.01. Figure 3 shows how the trial-averaged correlation is distributed over neurons. As expected, this trial averaged metric is not affected if we do use trial matching (0.30 + 0.01).
close match to
To quantify how the models capture the trial-to-trial variability, we then quantify how the distri- butions of neural activity and jaw movement consistent between data and model. So need to define the trial-matched Pearson correlation to compute a Pearson correlation between the dis- tribution of trial statistics TiS which unordered sets of trials. So compute the optimal assignment between trial pairs from the data and the recordings, and report the averaged Pear- son correlation overall trial pairs. Between the data and the model, we measure a trial rrial-macched Pearson correlation of 0.48 = 0.01, with performance ceiling at 0.52 + 0.01 obtained by compar- ing the training and testing set directly (see Figure c for details). For reference, the model without trial matching has lower trial- mrice- Pearson correlation 0.28 0.003.
Successful recovery of trial type distribution While the neuronal activity recorded, the behav- ioral response of the animal is also variable. When mice receive stimulation they perform correctly

---

A
12 Data о Model
B
Hit
Miss
Model
Data
0 Model without trial matching
wS1
wS2
wM1
10-
10
type
0.6
Firing Rate (Hz).
0.4
wM2
ALM
7.5
tjM1
Probability
0.0
Miss
Hit
CR
FA
A
<b>)
V
250
Active Hit
Quiet Hit
Active Hit
WS1
wS2
WM1
10-
10-
UMAP2
Firing Rate (Hz)
3-
Quiet Hit
wM2
ALM
tjM1
7.5
7.5
* Data
5.0
Model
2.5
UMAP 15
250 me
Figure 4: Emergent trial types. A. Trial type distribution from the recordings and the models. The whiskers show the 95% confidence interval of the trial type frequency. In this panel, trial types use the template matching method to avoid disadvantaging models without trial matching which have almost no variability in the jaw movement, see Figure A.3 B. Population and riel-averaged neuronal activity per area and per Hit and Miss trial type from the 400 simulated trials from the model against the averaged recordings from the testing set. C. Two-dimensional UMAP representation of the T of 400 simulated trials. The jaw movement is used in this representation. D. For the model, we separated the active hit and quiet hit trials based on their location the UMAP. For the data, separated the active hit and quiet hit trials based on the jaw movement in
with a 66% hit rate, while the absence of stimulus, mice still falsely lick with a 20% false alarm rate. Even in correct trials, the neural activity reflects variability which is correlated uninstructed jaw and tongue movements [
We evaluate the distribution of trial types (hit, miss, correct rejection, and false alarm) from fitted network model. Indeed, the 95% confidence intervals the estimated trial type frequencies always overlapping between the model and the data (see Figure AA). In this Figure, clas- sify the trial type with bearret-neighhor-lle
classifier using only the neural activity (see Figure A.3 In contrast, model without the trial matching would fail completely because it always pro- duces averaged trajectories instead of capturing the multi-modal variability possible classify trial types using jaw movement. To define equivalent trial types in the model, we rely on the presence or absence of the stimulation and classifier to identify lick action given the jaw movements. This classifier tron trained to predict lick action on the water dispenser given the recorded jaw movements (like in the data, the model can "move" the jaw without inducing lick action). After optimization with trial matching, since the jaw movement vy is contained in the fitted statistics T the distribution of jaw movement and the trial types are similar in the fitted model and the trial type distribution remains consistent. In Figure 41 we show population-averaged activity traces where the jaw is used to determine the trial type.
the data as seen in
Figure A With trial matching is
multi-layer percep-
Unsupervised discovery of modes of variability So we have analyzed whether variability among the main four trial types was expressed in the model, but the existence of these four trial types is not enforced explicitly in the loss function. Rather, the trial matching loss function aims to match the overall statistics of the distributions and it has discovered these four main modes of trial-

---

variability without explicit supervision. A consequence is that our model has possibly generated other modes of variability which are needed for the model to explain the full distribution of recorded data.
To summarize the full distribution of trials, we project the neural activity in 2D space in Figure EC. This representation recapitulates the modes of co-variability in all six areas and is not directly pos- sible to be generated from the dataset because trials from distinct sessions have non-overlapping recorded units. This representation becomes possible with a generative model like ours that can capture the distribution of trial-to-triall variability from incomplete recordings. Technically, we gen- erate 400 trials and we apply UMAP to the concatenation: (TwS1,k.... TALM,K) where each vector TA.:: includes now all neurons simulated in area A so the generated vectors contain data about all sessions and areas. To confirm that the generated representation of the distribution is consistent with the data we display template vectors for each trial condition e that are calculated directly on the recorded data trials. These templates are drawn with stars in Figure E and are computed as follows: the feature TR.. of this template vector computed by averaging the population activity of area A in all recorded trials from all sessions (see Appendix C for details), the averaged vectors are then concatenated and projected in the 2D UMAP space.
The consistency between the condition- condition-averaged templates and the full trial distribution is confirmed visually in Figure CC. The emerging distribution in this visualization is grouped in clusters. We observe that the template vectors representing the correct rejection, miss, and false alarm trials are located at the center of the corresponding cluster of generated trials. It also appears that the generated trial distribution provides richer information about the trial-to-trial variability that is hard to study using condition-averaged templates. For instance, we were curious why the generated hit trials are split into two clusters (see the two boxed clusters in Figure CC). Looking more closely at the difference between those two clusters of generated trials, it turned out that their difference can be explained by a simple feature: 85% of the generated hit trials on the left-hand cluster of panel AC have intense jaw movements during the delay period (max, lyy - y > 46 where is the standard deviation of ly' - y in the 200 ms before whisker stimulation). This criterion happens to be similar to the one used by experimentalists (2 to separate the hit trials in the recorded data, so we also refer to them the "active hit" and "quiet hit" trials and show the population activity in Figure 4D. This interesting correspondence highlights the potential for our method to study the structure of the trial-to-trial variability without a priori assumptions of the trial condition partitions. We conclude that our modeling approach can be used for a hypothesis-free identification of modes of trial-to-trial variability, even when they reflect task-irreeevant behavior.
"y'-
5 Discussion
We introduced generative modeling approach where data-constrained RSNN is fitted to multi- session electrophysiology data. The two major innovations of this paper are (1) the technical progress towards multi-session RSNN fitted with automatic differentiation, and (2) a trial match- ing loss function to match the trial-to-trial variability in recorded and generated data.
Interpretable mechanistic model of activity and behavior Our model has radically different objective in comparison with other deep-learning models: our RSNN is aimed to be biophysically in- terpretable. In FigureA.33| in the appendix, we illustrate a strong difference between our interpretable RSNN model and the more generic LFADS method which generates neural network activity with deep learning models. The most striking difference is that LFADS 2 relies on Gated recurrent units (GRU) which are abstract recurrent models for which there is no explicit mapping with real neurons. In our case, by construction there is a one-to-one mapping between an RSNN spiking unit and a recorded neuron. We further illustrated a consequence of our model in Figure A.5 we de- signed an extension of the artificial dataset where two putative areas are activated with a substantial delay of 20, or 200 ms. In real recordings, this would be only possible if unrecorded areas were maintaining some network activity during this long delay, but there is no plausible network hypoth- esis with two areas that could solve the task. When the delay between the areas is plausibly short (20 ms), the generated data from our RSNN model and LFADS are similarly accurate. For longer delays, the constraint in the RSNN model cannot allow for an implausible network solution by de- sign so our optimization collapses and does not provide any speculative circuit hypothesis for this artificial dataset unless we allow our network to have implausible long synaptic transmission delays,
9

---

see Figure A.5BB-C C. In contrast, LFADS can generate solutions independently of the delay duration, partly thanks to its bidirectional GRU units with efficient internal m gating which can handle very large delays. In the long term, we hope that the interpretability of our method will be able to capture biological mechanisms (e.g. predicting network structure, causal interactions between areas, and anatomical connectivity), but in this paper, we have focused on numerical and methodological questions which are getting us one step closer to this long-term objective.
debate is whether the brain computes with low- dimensional latent representations and how that is implemented in a neural circuit. Another reason for the LFADS to fit implausibly long delays, see Figure A.3 is inherent to the deep auto-encoder setting. By construction, LFADS generates trial-to-trial variability from low-dimensional latent rep- resentations, so the variability sourced by the latent variable which contains all the trial- sial-speciffic information. This is in stark contrast with our approach, where we see the emergence of structured manifolds in the trial-to-trial variability of the RSNN from the noise source of the network dynam- ics. In the UMAP representation of Figure FC. we visualize low-dimensional manifold of trial distribution although we did not enforce expliitlyy the presence of now-dimensional latent dynam- ics. Structure in the trial-to-trin toitial variability emerges because the RSNN is capable of transforming the unstructured noise sources (stochastic spikes and Gaussian input current) into a low-dimensional trial-to-trial variability typical variational auto-encoder setting would not achieve this.
Mechanisms of latent dynamies A
Note that it is also possible to add a random low-dimensional latent as source of low-dimensional variability in our RSNN as it is done with LFADS. In the FigurelA.4 we reproduce our results on the multi-session dataset from while assuming that all voltages vi have a trial-specific excitability offset Es,A using 5-dimensional gaussian noise 4h and one- seee-hidde-llaye perceptron Fe such that S. = (o..(1b)) The latent variable, therefore implements neuron-specific excitability offset coming from an unknown source. Unsurprisingly this latent noise model accelerates drastically the optimization, probably because fi,& is an efficient noise source for minimizing Leriat. More remarkably, the final fitting performance metrics are the same with and without the low-dimensional drive. It suggests that the extra assumption of low-dimensional input from unknown source is convenient but not necessary to generate realistic variability. Our first hypothesis that the recorded circuit alone produces the distribution of possible activity and behavior is sufficient. We speculate that providing this low-dimensional input to the model is sometimes counterproductive if the end goal is to identify the mechanism by which the circuit produces the low- -indeimesional dynamics.
Alternative loss functions to capture variability In our optimization, the trial-matching loss function constrains the network to produce realistic trial-to-trial variability. The main alternative methods to constrain the trial-to-trial variability would be likelihood- approaches 66 15 or
spike-GANs (3). These methods are appealing as they do not depend on the choice of trial statistics T TH Since these methods were never applied with multi-session data-connttaiieed
RSNN we explored how to extend them to our setting and compare the results. We tested these alternatives on the Artificial dataset in Figure A.2. The likelihood of the recorded spike trains (6]
15 cannot be defined with multiple sessions because cannot clamp neurons that are not recorded
(see S for details). The closest implementation that we could consider was to let the network simulate the data "freely" which requires, therefore, an optimal assignment between recorded and likelihood). With this loss function, we could not retrieve the bi-model hit versus miss trial type distribution unless it is optimized jointly with Cirial-
generated data, so it is a form of
We also tested the implementation of a spike-GAN GAN discriminator. In GANs the min-mac max optimiza- tion is notoriously hard to tune, and we were unable to train our generator with a generic spike-GAN discriminator from scratch (probably because the biological constraints of our generator affect the robustness of the optimization). In our hands, it only worked when the GAN discriminator was fed directly with the trial statistics T and the network was jointly fitted to the trial-averaged loss Lueuron- shows that GAN objective and the trial matching loss function hold similar role. We conclude that both of these clamping-fre free methods are promising to fit data contt-cont RSNNs. What differs between them, however, is that trial matching replaces the discriminator with the optimal assignment and the statistics T which are farameter-free, making them easy to use and numerically robust. It is conceivable for future work the best results are obtained by combining trial matching with other GAN-like generative methods.
10

---

Acknowledgments and Disclosure of Funding
This research was supported by the Swiss National Science Foundation (no. 31003A 182010, TMAG-3 7MM2-1.209271, 200020 20742), and Sinergia Project CRSII5_ 1916122. Many thanks to Lénaic Chizat, Vivien Seguy, James lsbister, Shuqi Wang and Vahid Esmaeili for their helpful discussions.
References
[1] Anne E Urai al. "Large-scale neural recordings call for new insights to link brain and behavior". In: Nature Neuroscience 25.1 (2022), pp. 11-19.
[2] Vahid Esmaeili et al. "Rapid suppression and sustained activation of distinct cortical regions for delayed sensory-triggered motor response". Neuron 109.13 (2021), pp. 2183-2201.
[3] Simon Musall et al. "Single-trial neural dynamics are dominated by richly varied move- ments". In: Nature Neuroscience 22.10 (2019), pp. 1677-1686.
[4] Carsen Stringer al. "Spontaneous behaviors drive multidimensional, brainwide activity". In: Science 364.6437 (2019), eaav7893.
[5] James B Isbister et al. "Clustering and control for adaptation uncovers time-warped spike time patterns in cortical networks in vivo". In: Scientific Reports 11.1 (2021), pp. 1-20.
[6] Jonathan W Pillow et al. "Spaaio-temporaa correlations and visual signalling in a complete neuronal population", In: Nature 454.7207 (2008), pp. 995-999.
[7] Gabriel Mahuas et "A new inference approach for training shallow and deep generalized linear models of noisy interacting neurons". In: Advances in Neural Information Processing Systems (2020), pp. 5070-5080.
[8] Guillaume Bellec et al. "Fitting summary statistics of neural data with a differentiable spiking network simulator". In: Advances in Neural Information Processing Systems (2021).
[9] Emre O Neftci, Hesham Mostafa, and Friedemann Zenke. "Surrogate gradient learning in spiking neural networks: Bringing the power of gradient- based optimization to spiking neural networks" In: JEEE Signal Processing Magazine 36.6 (2019). pp. 51-63.
[10] Guillaume Bellec et al. "Long short-term memory and learning-to-lue to-learn in networks of spik- ing neurons". In: Advances in Neural Information Processing Systems (2018).
[11] Gabriel Peyré, Marco Cuturi, al. "Computational optimal transport: With applications to data science". In: Foundations and Trends@ in Machine Learning 11.5-6 (2019), pp. 355- 607.
[12] Marco Cuturi. "Sinkhorn distances: Lightspeed computation of optimal transport". In: Ad- vances in Neural Information Processing Systems (2013).
[13] Jean Feydy et al. "Interpolating between Optimal Transport and MMD using Sinkhorn Di- vergences". In: The 22nd International Conference on Artificial Intelligence and Statistics. 2019, pp. 2681-2690.
[14] Wulfram Gerstner. "Time structure of the activity in neural network models". In: Physical Review E 51.1 (1995), 738.
[15] Christian Pozzorini et al. "Automated tigh-throughput characterization of single neurons by means of simplified spiking models". In: PLoS Computational Biology 11.6 (2015), e1004275.
[16] Corinne Teeter al. "Generalized leaky integrate-and-fire models classify multiple neuron types". In: Nature Communications 9.1 (2018), 709.
[17] Matthew G Perich et al. "Inferring brain-wide interactions using data-constrained recurrent
neural network models". In: BioRxiv (2020), pp. 2020-12.
[18] David Sussillo and Larry F Abbott. "Generating coherent patterns of activity from chaotic neural networks". In: Neuron 63.4 (2009), pp. 544 544-557.
[19] Benjamin J Arthur et "A scalable implementation of the recursive least-squares algorithm for training spiking neural networks". In: Frontiers in Neuroinformatics 17 (2023).
[20] Christopher M Kim et al. "Distributing task-related neural activity across cortical network through task-independent connections". In: Nature Communications 14.1 (2023), p. 2851,
[21] Dongsung Huh and Terrence Sejnowski. "Gradient descent for spiking neural networks". In: Advances in Neural Information Processing Systems (2018).

---

[22] Chethan Pandarinath et al. "Inferring single-trin trial neural population dynamics using sequential ritooeencoers.... In: Nature Methods 15.10 (2018), pp. 605-815.
[23] Ding Zhou and Xue-Xin Wei. "Learning identifiable and interpretable latent models of high- dimensional neural activity using pi- pi-V.E". In: Advances in Neural Information Processing Systems (2020), Pp. 7234-7247.
[24] Byron M Yu et al. "Gaussian ceeussian-process factor analysis for low-dimensional single-trial analysis of neural population activity". In: Advances Neural Information Processing Systems 21 (2008).
[25] Shuqi Wang et al. "Mesoscopic modeling of hidden spiking neurons". In: Advances in Neural Information Processing Systems 35 (2022), pp. 23566-23579)
[26]
Brian Depasquale et al. "The centrality of population-leve factors to network computation is demonstrated by a versatile approach for training spiking networks". In: Neuron (2023).
[27] Manuel Molano-Mazoo
al. "Synthesizing realistic neural population activity patterns using Generative Adversarial Networks". In: International Conference on Learning Representa-
tions (2018).
[28]
Adrian Valente, Jonathan W Pillow, and Srdjan Ostojic. "Extracting computational mecha- nisms from neural data using low-rank RNNs". In: Advances in Neural Information Process- ing Systems 35 (2022), pp. 24072-24086.
[29] Valentin Schmutz, Johanni Brea, and Wulfram Gerstner. "Convergence of redundancy-free spiking neural networks to rate networks". In: arXiv preprint arX1v:2303.05174 (2023).
[30] Kanaka Rajan, Christopher D Harvey, and David W Tank. "Recurrent network models of sequence generation and memory". In: Neuron 90.1 (2016), pp. 128-142.
[31] Aaron S Andalman al. "Neuronal dynamics regulating brain and behavioral state transi- tions". In: Cell 177.4 (2019), pp. 970-985.
[32]
Constantin-Klemen Lurz et al. "Generalization in data-driven models of primary visual cor- tex". In: International Conference on Learning Representations (2021).
[33] Poornima Ramesh, Mohamad Atayi, and Jakob Macke. "Adversarial training of neural encoding models population spike trains". In: NeurlPS Workshop on Real Neurons & Hidden Units: Future directions at the intersection of neuroxcience and artificial intelligence.
2019.
[34] Vahid Esmaeili al. "Learning-relats related congruent and incongruent changes of excitation and inhibition in distinct cortical areas". In: PLoS Biology 20.5 (2022), e3001667.
[35] Harold W Kuhn. "The Hungarian method for the assignment problem". In: Naval Research Logistics Quarterly 2.1-2 (1955), pp. 83-97.
[36] Alexandre Défossez et al. "High fidelity neural audio compression". In: arXiv preprint
arXiv:2210. 10.33330 (2022).
12


Reviews:
[{"review_id": "O3Izzpr837", "reviewer": "Reviewer_Ec9F", "summary": "This paper extends previous works on generative model of neural data using a recurrent spiking neural network.\nThe main contributions are design of new loss function based on optimal transport to match the trial-to-trial variability in the data, and the approach to extend RSNN generative model to leverage recordings from multiple sessions. \nThe proposed method is evaluated on both synthetic and real neural data, and authors showed that the method can be used for a hypothesis-free identification of modes of trial-to-trial variability.\n", "questions": "- Definition of $T_{trial}$ in trial matching, Equation 6, and template matching,  line 256, seems to differ.  Is subset of $T_{trial}$ (containing averaged response of area recorded in the particular session/trial) used to compute $L_{trial}$?\n- (typo)  $T_{trial}$ → $\\mathcal{T}_{trial}$  (line 222 and 244)", "limitations": "Yes", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- Potential usage of the method in answering neuroscience questions is appropriately illustrated from the real data results\n- Extension of existing RSNN to data constrained model is novel, and will be valuable step in increasing the usability of related methods in neuroscience.\n- Previous work is appropriately discussed and the contributions are clearly stated\n\n\n", "weaknesses": "- Although the authors mention in the discussion section the difference in the approach between the proposed method and existing methods, the (potential) advantage of the proposed method (biophysically interpretability) is not clearly illustrated. \n\n\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for reading our work in depth. We address below the main weakness and questions raised by the reviewer.  \n\n- *\"Although the authors mention in the discussion section the difference in the approach between the proposed method and existing methods, the (potential) advantage of the proposed method (biophysically interpretability) is not clearly illustrated.\"*:  \n\nTo illustrate more clearly the potential advantage of our method we prepared an extra simulation during the rebuttal period. In Figure 1 shown in the PDF attached to our general rebuttal, we generated an artificial dataset with two putative areas A1 and A2. Area A1 always responds acutely to a stimulus while area A2 responds 80% of the time with a long delay (100, 200 and 300ms). A deep learning method like LFADS [20] can fit and generate the activity perfectly in all cases. For our method, we see that the fitting is as good as LFADS when the delay between the two areas is within a plausible range, and collapses when the network assumption of two areas connected with short 4ms synaptic delays and realistic membrane time constants makes it impossible to produce the activity. This highlights that a successful fit with our model provides a plausible circuit hypothesis relying solely on the recorded neurons. \n\n-  *\"Definition of  $T_{\\mathrm{trial}} $  in trial matching, Equation 6, and template matching, line 256, seems to differ. Is subset of $T_{\\mathrm{trial}}$ (containing averaged response of area recorded in the particular session/trial) used to compute $\\mathcal{L}_{\\mathrm{trial}}$?\"*\n\nThank you for spotting this subtle difference. Indeed, $T_{trial}$ in equation (6) only considers the areas and neurons which were recorded in a specific session S, while in line 256, as used in the UMAP plot, we compute the feature $T_{A}$ from all neurons from area A of the simulated network. To clarify this in the main text we suggest distinguishing the two by adding a subscript S when $\\mathcal{T}_{trial}^S$ contains only the features recorded in session S and explaining the difference. \n\n- Concerning the typo reported by the reviewer this will be corrected."}, {"review_id": "oQWzfeYfXk", "reviewer": "Reviewer_Myz7", "summary": "This paper presents a framework for training a large recurrent spiking neural network on multi-session recordings by leveraging \nan optimal transport-based trial matching between the real data and generated data. This model is use to model the cortical sensory-motor pathway during a tactile detection task. ", "questions": "1. While the proposed model is promising, there are questions about how the model can be scaled to larger sets of recordings, as well as more heterogeneous recordings (example: groups of animals differentiated by age, state, disease propagation etc...) Would a single model which ignores the individual differences be adequate? \n2. Most neural recordings are currently structured by trials. To truly capture the full breadth of neural code complexity, it is necessary to study these neural dynamics in more unconstrained / complex / free-behavior settings, in which case the notion of trial no longer exists. Might the authors have ideas on how their approach can be adapted to such settings?", "limitations": "No limitations identified.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. In neuroscience, the challenge of training a model on multiple sessions is very relevant. It is expected to record data over multiple sessions and across multiple animals, hence having tools to jointly analyze the underlying dynamics of the neural population across different regions and recordings is important and critical. The use of trial averaging, as motivated by the authors, can be limiting and fails to capture the trial variability during complex behavior. The model presented in this work enables the analysis of the neural population dynamics. The design of the model is clearly well thought-out but remains very simple and is strongly supported by 1) good arguments 2) great visualizations 3) supporting empirical results.\n2. The empirical results support the effectiveness of the model at revealing the underlying modes of trial variability. The unsupervised discovery of new modes is also promising.\n3. The discussion section is very thorough. In particular the comparaisons with LFADS are interesting and insightful. This work challenges the idea that low-dimensional spaces are required to be able to interpret neural dynamics.", "weaknesses": "1. The main weakness of this work is that it was not directly compared to other baselines or test on other recordings during a different behavioral task.\n2. The model requires the selection of multiple hyperparameters. One in particular is the number of neurons (1500). The total number of neurons across all recordings is 4415 so there is a huge decrease in the number of units that are modeled. It is unclear how this choice can be made as the size of the dataset increases or for different datasets. The current model takes 3 days of training, are there expected computational limitations for the number of neurons that can be modeled? ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for the detailed review, and the positive feedback. Regarding the weaknesses of the paper raised by the reviewer: \n1. About the comparison to other baselines, in case this was overlooked by the reviewer, we would first like to highlight that we had compared trial matching with spike-GANs and likelihood-like approaches in the discussion and the appendix section B, see Figure A2. As far as we know, these algorithms were the only two alternatives capable of guiding the gradient-based optimization of a spiking activity generator like ours. We believe these were the best references to benchmark the efficiency of the trial-matching loss function.  \n \n    We agree that the application and comparison of our method to other tasks would be insightful. During the time available for the rebuttal, we prepared an extra simulation and comparison with the LFADS baseline. This new comparison highlights the difference between our biophysically interpretable models and typical deep learning methods. Figure 1 of the attached PDF shows our results. We generated an artificial dataset where area A1 always responds acutely to a stimulus while area A2 responds 80% of the time after a long delay (100, 200 and 300ms). Our methods achieve Pearson correlations similar to LFADS for the 100ms delay showing the efficiency of our numerical method. A striking difference emerges with the larger delays: when the assumption of our network model cannot provide a plausible circuit hypothesis, the trial-matched Pearson correlation saturates at a lower level. In contrast, the performance of LFADS is not affected, which is expected because the computation is not meant to model a plausible network computation.  \n \n2. In this dataset, we chose the number 1500 because, as written in line 119 in the paper, 1) we wanted to have an equal number of neurons per area, and 2) the area with the least number of recorded neurons is wS1 with 258 units (see ref [9]). So, we took 250 neurons per area. Also, the total number of units from the areas of interest is 2724, while the number of 4415 that we had written includes areas that we did not consider to be relevant in our simulation. We agree that this was not clear and should be clarified in the paper.  \n \n    We agree with the reviewer that it will be important to scale up our method to larger datasets. So far, we have not had to invest efforts to improve the simulation efficiency, but it seems that the bottleneck to scale for large datasets would be to overcome the memory limit of a GPU. We believe multiple implementation strategies could help. For instance, we currently use full backpropagation through time over 700 timesteps which could be split into shorter time segments using truncated backpropagation. We could also optimize the state representation to be more memory efficient when integrating the network dynamics. Overall, we think that there is a broad margin for improvement.  \n\nQuestions: \n1. As explained above, we believe our method can scale to larger networks. Regarding modeling cross-subject differences, it is commonly assumed when using a dataset with dozens of different rodent subjects that the circuit being studied is shared across subjects. We make the same assumption in our model. However, our method also offers possible ways to compare groups of animals differentiated by age (or other features mentioned by the reviewer), by fitting distinct models for subsets of subjects.  \n\n2. We believe that our method could easily be adapted to behavior-free recordings where the trial structure is not straightforward. The notion of trial is only used in our method to define a basis of time alignment to construct the feature vectors. If trials are not defined in the dataset, one could for instance generate \"pseudo-trials\" aligned to specific events (e.g. onset of jaw movements and/or peaks of activity in sensory areas) or consider other features. "}, {"review_id": "n2ABauqbZS", "reviewer": "Reviewer_Rsk5", "summary": "The authors present a method to fit recurrent spiking neural network models to experimental data, extending prior work to multi-session recordings that involve different subpopulations in each session. To capture single trial variability, they introduce a trial matching loss function. The technique is applied to neural recordings from six cortical areas in the sensory-motor pathway during a delayed whisker detection task. The model successfully captures trial-averaged firing rates of individual neurons and population averages. In addition to the expected matching results enforced by the loss function, the model reveals an unexpected mode of variability related to movements during the delay period, resulting in two distinct clusters within hit trials.", "questions": "-    Are all results based on held-out data? If so, how was the data held out?\n-    What are the new predictions that are not expected results from the training method? Aren't the four main trial types obvious due to the optimization objective including sensory neurons that discriminate between stimulus present/absent and jaw movements that discriminate between lick/no lick, resulting in true positive, false negative, false positive, and true negative categories?\n-    If you hold out all neurons from one or multiple entire areas, can your model accurately predict the activity in those areas?\n-    Since you use population averages to study trial variability, wouldn't a simple population-level approach using six units, each representing a cortical area, yield the same results?", "limitations": "The authors acknowledge that this paper represents only a step towards their long-term goal of developing biophysically interpretable mechanistic models.", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- Models trial variability beyond trial-averaged firing rates.\n- Considers a large-scale, cortex-wide multi-region dataset.", "weaknesses": "- Analysis of single trial variability still relies on population averages rather than individual trial averages.\n- Lacks a quantitative comparison to other methods.\n- Represents an incremental contribution, with the closest related work being Bellec et al., NeurIPS 2021.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for evaluating our work. We reply to each weakness and question below.  \nWeaknesses: \n\n- We are not sure that we understand correctly what the reviewer meant by \"individual trial averages\". In case the reviewer meant neuron specific individual trials, we tried to fit our model with neuron and time specific features in the form of the trial-matched likelihood as shown in the appendix section B, however it failed to fit the trial-to-trial variability. This most likely happens due to the noise of spiking data.  \n\n- Regarding comparison with other methods, in case this was overlooked by the reviewer, we made comparisons between trial matching versus spike-GANs and likelihood-like methods in the discussion and the appendix section B, see Figure A2.  \nDuring the rebuttal period, we have prepared an additional comparison to provide quantitative results between our method and LFADS, the results are provided in the PDF attached to the general rebuttal. If the model assumption can find a model hypothesis that can generate the recorded activity, it seems that our method is not less efficient than LFADS. \n \n- This paper was indeed an inspiring reference for our work, but here we explore and solve new problems (e.g. modeling trial variability and considering a cortex-wide multi-region dataset, the two strengths that are highlighted by the reviewer himself/herself). To be precise, the main methodological innovations of this paper which were not present in [3] are: \n    1) the dataset that was used in [3] was a single recorded session with a single recorded area. The method presented here is scaled to six areas and 28 recording sessions which raised multiple theoretical and practical problems. \n    2) The statistics used in ref [3] were all trial averaged and the method was incapable of producing a realistic distribution of trial-to-trial variability which is central to this paper. \n    3) There are many additional subtleties detailed in our paper/code that were not performed in ref [3]. To cite a few: we use excitatory and inhibitory neurons which are matched coherently to the excitatory/inhibitory cell types detected from the action potential waveform; the weight matrices are restricted to respect Dale’s law; the synaptic delays are restricted in a plausible range and the voltage dynamics follow LIF dynamics with plausible time constants. \n\nQuestions: \n- Yes, as stated in the paper lines 214, 216, 226 and in the caption Figure 4B all the reported results are always quantified on the held-out test trials. Those test trials consist of 25% of the trials which are randomly selected from each session. To clarify this, we suggest adding to the paper a sentence to describe the selection of the held-out test trials. \n \n- The split between the four trial types serves in our paper as an explicit way to quantify our success at generating the right distribution of trials. Further, our results go beyond showing \"just\" the emergence of trial types related to stimulation or movement: 1) the network dynamics are broadly consistent with all trial types and in all areas (see Figure 4B). 2) We have also reported in Figure 4 C-D and between lines 248 and 276 how trial-matching can capture other sources of trial variability like spontaneous motion which is independent of the trial condition structure. Independently of trial types, our method captures all the main modes of trial variability: it was a truly coincidental observation that the highlighted cluster in 4C represented active and quiet hit trials which happen to be trials that were studied in the initial dataset.  \n \n\n- Predicting activity in unrecorded areas is an interesting and open question. To some extent, we would like to argue that our model does it in the following sense: given a specific session of the same task where for example only areas wS1 and wS2 are recorded, our model appears to be capable to interpolate faithfully the remaining network activity. This is supported in part by Figure 4 B-C: if the recordings from wS1 and wS2 are sufficient to match the trials relevant to the simulated trials, the corresponding simulated trials will share a correct trial type and many other features reflected in other areas (see Figure 4D). \n \n - A population-level model would be a fundamentally different model. Given that we take the biophysical interpretability of the model very seriously in this paper, the relevant population level would ideally be a biophysical mean-field like process. Deriving precise and convincing mean-field models is theoretically difficult and -- despite encouraging progress [33] -- we are not aware of a data fitting method using accurate mean-field models with convincing biophysical interpretability.  \n\n    However, to demonstrate that a population model would work at the price of many advantages of our model, we tried to sketch out a sensible \"simple population level\" model. Each area is represented by one or few recurrent tanh units, and these units are mapped to the recorded neurons from the corresponding areas with a linear projection. The reviewer can see the results of our implementation in the General rebuttal. In short, a single tanh unit per area was not enough to yield a competitive generative model even on the artificial dataset. More importantly, this model compromises all the biophysically interpretable aspects of the model (e.g. spiking communication, cell types, Dale's law, and neuron-to-neuron mapping with the data) and we do not know what simple model could easily reconcile that."}, {"review_id": "P91RSu6siy", "reviewer": "Reviewer_Cbej", "summary": "This paper considers fitting high-fidelity recurrent spiking neural networks to multiple populations of neurons across sessions. The networks include LIF neurons for each fitted neuron and model within & across region interactions at fine timescales, and include a soft-thresholding mechanism for neural spiking. The model is fit to match the responses of the recorded neurons both across trials and to generate variability across trials that is consistent with the variability observed in the data. For the latter, the authors propose an optimal-transport based loss on the distributions of sampled recorded trials and trials simulated from the model. \n\nThe authors validate the proposed method in simulation, showing that the trial matching term in the objective is critical for capturing variability across trial types in neural data. Then, they apply it to a large multi-population neural dataset, showing the model captures trends across trial types in this dataset as well.", "questions": "* Session stitching. Can the authors test their method for session stitching in a simulation, where they stitch together neurons recorded in different sessions as is done in the real data? This would help guide when to trust the results of stitching together non-simultaneously recorded neurons.\n\n* When computing the trial target $T^\\prime_{trial}$ for neurons from different sessions, are the trials matched based on condition before summing neurons from different sessions? For example, `Hit` trials are combined from different sessions, but not `Hit` trials and `Miss` trials.", "limitations": "Yes", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The motivation and methods are clear and novel, and the resulting model appears to be significant for many neuroscience modeling tasks.", "weaknesses": "Fitting the model jointly to neurons that are not simultaneously recorded (and not necessarily from the same animal) seems to be conceptually at odds with the goal of capturing trial by trial variability. By definition, neurons from different sessions should be considered different for trial to trial variability. I would argue that a more accurate term for the loss is \"condition to condition\" variability. The authors show it is possible to capture responses across conditions using their objective. It would be great for the authors to expand on this point.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for reading our paper and providing useful feedback.  \nGiven the second question of the reviewer, we wonder whether there might be a misunderstanding on the nature of the trial-matching loss function. It appears that the reviewer understands that we only model \"condition-to-condition\" variability and not \"trial-to-trial\" variability. But in fact, the concept of trial condition (hit/miss/etc...) is only used for measuring the validity of the model after optimization, and the concept of condition does not appear in the definition of the trial-matching loss function or any other algorithmic step. The trial-matching loss function is only a distance across distributions which happens to capture the trial types among other things. This is the reason why we claim the non-trivial important result that our model can generate the correct distribution of simulated trial types (correct proportion in Figure 4A and correct co-variability of activity in the recorded areas in Figure 4B). \n \nBelow we respond to the weaknesses and questions raised by the reviewer:  \n\n- *“Fitting the model jointly to neurons that are not simultaneously recorded (and not necessarily from the same animal) seems to be conceptually at odds with the goal of capturing trial-by-trial variability.”*\n\nFollowing the tradition of rodent research, we consider that we study a system (the sensory-motor pathway) for which the features are shared among animals of the same species and seek to understand the common feature of the recorded circuits and not their differences. For this reason, the implicit statistical assumption of our model is that each session corresponds to a different electrode placement in the same circuit. Since the point by the reviewer is linked to some fundamental modeling choices, we are eager to clarify this in the paper. \n \n\n- *“By definition, neurons from different sessions should be considered different for trial to trial variability.”*\n\nWe are not sure if we correctly understand this remark. In the model, neurons from different sessions are indeed considered differently for the trial-to-trial variability. During training for the calculation of the loss function, we consider that neurons from one session are evolving while considering the neurons from the other sessions “hidden”. The loss is constrained on the neuron-specific trial averaged data, every neuron sees the rest as “hidden”, and on the session-specific trial-specific population averages, neurons from one session see the rest of neurons as “hidden”. \n\n- *“I would argue that a more accurate term for the loss is \"condition to condition\" variability. The authors show it is possible to capture responses across conditions using their objective. It would be great for the authors to expand on this point.”*\n\nAs explained above, we believe that it would be reductive to limit the effect of the algorithm to trial types of conditions. We are glad that the reviewer appreciates seeing the responses of neurons across conditions. If you think of specific additional material to clarify our algorithm, please let us know. \n\n- *“Session stitching. Can the authors test their method for session stitching in a simulation, where they stitch together neurons recorded in different sessions as is done in the real data? This would help guide when to trust the results of stitching together non-simultaneously recorded neurons.”*\n\nWe agree that a large-scale model recovery experiment would be great to quantify how closely the algorithm is capable of capturing the network structure despite the limited number of recorded neurons per session which is compensated with session stitching. At this point, we can mainly say that session stitching was crucial to leverage all the recorded data available and reconstruct a model of the entire sensory-motor pathway although we never had more than two areas recorded simultaneously.  This would not have been possible without session stitching. "}]

Author rebuttal:
We would like to thank all reviewers for reading our paper in depth and providing insightful feedback. All reviewers understood the paper and wrote that sections were "clear and novel" and "appropriately discussed and clearly stated". Three of the reviewers commented that the submitted manuscript was lacking baseline comparisons. This could be partly explained because our comparisons with the spike-GAN and the likelihood-based approach were in the discussion and in the appendix, and, in the revision, we will point to them more clearly.  We also ran extra simulations during the rebuttal to compare with other baselines and tackle important questions raised by the reviewers: 1) a new simulation will illustrate the difference between our "biophysically interpretable" model and other approaches, as requested by reviewer Ec9F and 2) we added a population-level model baseline as requested by reviewer Rsk5. The new simulations and associated results are displayed in the PDF attached to the present rebuttal. We are grateful to the reviewers for their comments which help improve the quality of the manuscript. 

**New simulation 1: comparison with uninterpretable generative models.** To illustrate the fundamental difference between biophysically interpretable methods against more abstract methods like LFADS, we show in the PDF attached to the rebuttal that our method may only find a functional circuit hypothesis when there exists a biologically plausible circuit hypothesis that can solve the task. We generated a variant of the previous artificial dataset but assumed that the putative motor area (A2) has a bump of activity that arrives with a variable delay from the putative sensory area (A1). When the delay is very large (i.e. 300 ms), a network with direct synaptic connectivity of 2-4ms delay would not be able to bridge this delay, see Figure 1 in attached pdf. So the optimization collapses for large delays and our method does not provide any network solution which explains the artificial data. In contrast, for a method like LFADS, since there is no explicit causal or biophysical interaction that is explicitly required between the two areas, the activity of both areas is well explained independently of the delay and even the input. For short delays, however, we found that our network is as good as LFADS for the metrics that we considered (LFADS was implemented with the tutorial https://github.com/snel-repo/lfads-tutorial, most of the hyper-parameters were the default). 

**New simulation 2: comparison with simpler population-level models.** To comment on the difference with population-level modeling, we tried to sketch out a sensible “simple population level” model with crude and imprecise approximations. Compromising many biophysically interpretable aspects of our model (e.g. spiking communication, cell types, Dale's law, and neuron-to-neuron mapping with the data), we have now implemented a spiking activity generator driven by a tanh-RNN with a time constant of 10ms. Each cortical area is represented by one or a few of such units, and the units are recurrently connected to model the interaction between areas. The model generates spikes as Poisson spike trains with a linear projection of the unit's activity to the corresponding area (note that true mean-field models are much more complex [33]).  
Testing this model on our artificial dataset, we noticed that one unit per area is not enough to model the two putative areas – it could not explain fully the trial-averaged statistics and failed to capture single trial statistics. It becomes closer to our results using 2 tanh units per area as summarized in Figure 2 in the attached PDF. On the real dataset, we expect that this population-based model might predict activity well given many more than 1 tanh-unit per area. In our mind the main failure of this type of population-level model is not the quantitative fit, rather, there is a tanh-unit that has no explicit mapping with any biophysical substrate. We would argue that this type of model has unfortunately no satisfactory meaning and blurs the line further between biophysically interpretable models and deep learning models like LSTM or transformers. We believe that our choice of neuron- and spike-level modeling is "a sweet spot" enabling simultaneous biophysical interpretability and, as we show, reasonably efficient optimization for parameter fitting. 

Meta-review: The reviewers were unanimous in their appraisal of this paper's contributions as above the bar for acceptance to NeurIPS, and I congratulate the authors on their detailed and comprehensive rebuttals, which were critical in leading some reviewers to raise their scores.  I'm pleased to report that this paper has been accepted to NeurIPS.  Congratulations!  Please revise the manuscript to address all reviewer comments and questions.
Decision: Accept (poster)