1 Introduction

A fundamental question in neuroscience is how cortical circuit mechanisms drive perception and behavior. To tackle this question, the field has been collecting large-scale electrophysiology datasets under reproducible experimental settings (Siegle et al., 2021; Esmaeili et al., 2021; Urai et al., 2022; International Brain Laboratory et al., 2023). However, the field lacks data-grounded modeling approaches to generate and test hypotheses on the causal role of neuronal and circuit-level mechanisms. To leverage the high information density of contemporary recordings, we need both (1) modeling approaches that scale well with data, and (2) metrics to quantify when the models provide a plausible mechanism for the observed phenomena.

Biophysical simulations have been crucial for our understanding of single-cell mechanisms (Hodgkin, 1958), and have been used to describe interactions across cortical layers, columns, and areas (Markram et al., 2015; Billeh et al., 2020; Isbister et al., 2023; Chen et al., 2022; Rimehaug et al., 2023; Fraile et al., 2023; Spieler et al., 2023). A promising approach to constrain models to electrophysiological data lies in the optimization of the simulation parameters by gradient descent. These methods were successful in quantitatively classifying functional cell types (Pozzorini et al., 2015; Teeter et al., 2018), and modeling micro-circuit interactions (Pillow et al., 2008; Deny et al., 2017; Mahuas et al., 2020). To bridge the gap from single neurons or small retinal networks to cortical recordings in vivo, recent studies made substantial progress towards data-constrained recurrent neural network (RNN) models (Perich et al., 2020; Bellec et al., 2021; Arthur et al., 2023; Kim et al., 2023; Dinc et al., 2023; Sourmpis et al., 2023; Pals et al., 2024). In this line of work, neurons in the RNN are mapped one-to-one to recorded cells and optimized by gradient descent to predict recorded activity at large.

An important question is whether these data-constrained RNNs can reveal a truthful mechanism of neuronal activity and behavior. By construction, the RNNs can generate brain-like network activity, but how can we measure whether the reconstructed network faithfully represents the biophysical mechanism? To answer this question, we submit a range of RNN reconstruction methods to a difficult perturbation test : we measure the similarity of the network response to unseen perturbations in the RNN and the recorded biological circuit.

Optogenetics is a powerful tool to induce precise causal perturbations in vivo (Esmaeili et al., 2021; Guo et al., 2014). It involves the expression of light-sensitive ion channels (Aravanis et al., 2007), such as channelrhodopsins, in specific populations of neurons (e.g., excitatory/pyramidal or inhibitory/parvalbumin-expressing). In this paper, we use datasets combining dense electrophysiological recordings with optogenetic perturbations to evaluate RNN reconstruction methods. Since the neurons in our RNNs are mapped one-to-one to the recorded cells, we can model optogenetic perturbations targeting the same cell-types and areas as done in vivo. Yet, we observe that the similarity between the simulated and recorded perturbations varies greatly depending on the reconstruction methods.

Most prominently, we study two opposite types of RNN specifications. First, as a control model, we consider a traditional sigmoidal RNN (σRNN) which is arguably the most common choice for contemporary data-constrained RNNs (Perich et al., 2020; Arthur et al., 2023; Pals et al., 2024); and second, we develop a model with biologically informed inductive biases (bioRNN): (1) neuronal dynamics follow a simplified spiking neuron model, and (2) neurons associated with fast-spiking inhibitory cells have short-distance inhibitory projections (other neurons are excitatory with both local and long-range interareal connectivity). Following Neftci et al. (2019); Bellec et al. (2018b, 2021); Sourmpis et al. (2023), we adapt gradient descent techniques to optimize the bioRNN parameters of neurons and synapses to explain the recorded neural activity and behavior.

Strikingly, we find that the bioRNN is more robust to perturbations than the σRNN. This is nontrivial because it is in direct contradiction with other metrics often used in the field: the σRNN simulation achieves higher similarity with unseen recorded trials before perturbation, but lower than the bioRNN on perturbed trials. This contradiction is confirmed both on synthetic and in vivo datasets. To analyze this result, we submit a spectrum of intermediate bioRNN models to the same perturbation tests, and identify two bioRNN model features that are most important to improve robustness to perturbation: (1) Dale’s law (the cell type constrains the sign of the connections (Eccles, 1976)), and (2) local-only inhibition (inhibitory neurons do not project to other cortical areas). We also find that other biological inductive biases like spiking neuron dynamics and a sparsity prior improve the robustness to perturbation but to a lesser extent.

Furthermore, a perturbation-robust bioRNN enables the prediction of the causal effect of perturbations in the recorded circuit. This becomes particularly interesting with micro-perturbations (µ-perturbation) targeting dozens of neurons in a small time window. We show in silico that the causal effect of µ-perturbations can be well approximated by the RNN gradients, which has two important implications for experimental neuroscience: (1) in a close-loop experimental setup, we can use RNN gradients to target a µ-perturbation which optimally increase (or decrease) movement in a simulated mouse (this is demonstrated in silico); (2) our RNN reconstruction methodology enables the estimation of gradients of the recorded circuit. So fitting RNNs becomes a tool to measure “brain gradients” and potentially relate contemporary in vivo measurement to decades of theoretical results from machine learning where the gradient is a foundational concept (LeCun et al., 2015; Richards and Kording, 2023).

2 Results

2.1 Reconstructed Networks: biological inductive biases strengthen robustness to perturbations

Synthetic dataset for challenging causal inference

We build a toy synthetic dataset to formalize how we intend to reverse engineer the mechanism of a recorded circuit using optogenetic perturbations and RNN reconstruction methods. It also serves as the first dataset to evaluate our network reconstruction methods. This toy example represents a simplified version of large-scale cortical recordings from multiple brain areas during a low-dimensional instructed behavior (Steinmetz et al., 2019; Esmaeili et al., 2021; International Brain Laboratory et al., 2023), similarly to the vivo dataset of a GO/No-Go task Esmaeili et al. (2021) analyzed in the next section. Let’s consider two areas A and B which are either transiently active together (“hit trial” occurring with frequency p) or quiescent together (“miss trial” occurring with probability 1 − p). Since the two areas are active or in-active together, it is hard to infer if they are connected in a feedforward or recurrent fashion. In Method 4.1, we describe a theoretical example where it is impossible to decide between opposing mechanistic hypothesis (feedforward or recurrent) when recording only the macroscopic activations of areas A and B. In this case, performing optogenetic inactivation of one area is decisive to distinguish between the feedforward or recurrent hypothesis.

To generate artificial spike train recordings that capture this problem, we design two reference circuits (RefCircs) from which we can record the spike trains. Each RefCirc consist of two populations of 250 spiking neurons (80% are excitatory) representing areas A and B. To highlight the importance of optogenetic perturbations as in the Methods 4.1, the first circuit RefCirc1 is feedforward and the second Refcirc2 is recurrent: RefCirc1 (and not RefCirc2) has strictly zero feedback connections from B to A. Yet, the two RefCircs are almost identical without optogenetic perturbations: each neuron in RefCirc1 has been constructed to have an almost identical trial-averaged activity as the corresponding neuron in RefCirc2; and in response to a stimulus, the circuits display a similar bi-modal hit-or-miss response with a hit trial frequency p ≈ 50%. We consider that a trial is a hit if area A is active (averaged firing rate above 8Hz; Defining a “hit” trial based on area A is equivalent to saying that both areas need to be active during unperturbed trials with this dataset. But excluding B in this definition avoids that the hit rate is trivially impacted when manipulating the activity of B with optogenetic perturbations.). To simulate optogenetic inactivations of an area in the RefCircs, we inject a transient current into the inhibitory neurons, modeling the opening of light-sensitive ion channels (Symmetrically, an optogenetic activation is simulated as a positive current injected into excitatory cells). Figure S1 shows that optogenetic perturbations in area B reflect the presence or absence of feedback connections which differs in RefCirc 1 and 2. Methods 4.4 provides more details on the construction of the artificial circuits. Our perturbation test will be to compare the effect of optogenetic perturbations in the reconstructed RNN and in the reference circuits RefCirc1 and 2.

Network reconstruction methodology (synthetic dataset)

To reconstruct the recorded circuits with an RNN, we record activity from the spiking RefCircs, and optimize the parameters of an RNN to generate highly similar network activity. The whole reconstruction method is summarized graphically in panel A of Figure 1. In the simplest cases, the RNN is specified as a sigmoidal network model (Rosenblatt, 1960; Elman, 1990): σRNN1 and σRNN2 are optimized to reproduce the recording from RefCirc1 and RefCirc2 respectively. In this synthetic dataset, the reconstructed σRNNs have the same size as the RefCircs (500 neurons) and sigmoidal neurons are mapped one-to-one with RefCirc neurons (20% are mapped to inhibitory RefCirc neurons). They are initialized with all-to-all connectivity and are therefore blind to the structural difference of the RefCirc1 and 2 (feed-forward or recurrent). From each of the two RefCircs, we store 2,000 trials of simultaneous spike train recordings of all 500 neurons (step 1 in Fig. 1A). Half of the trials are used as training set and will be the basis for our data-driven RNN optimization. The second half of the recorded trials form the testing set and are used to evaluate the quality of the reconstruction before perturbations.

Figure 1:

We optimize the synaptic “weights” of the σRNN to minimize the difference between its activity and that of the RefCirc (step 2 in Fig 1A, see Methods). The optimization combines three loss functions defined mathematically in Methods 4.5: (i) the neuron-specific loss function ℒ_neuron is the mean-square error of the trial-averaged neural activity (e.g. the PSTH) between σRNN and RefCirc neurons. (ii) To account for fluctuations of the single-trial network activity, we use a trial-specific loss function ℒ_trial, which is the distance between the distribution of single trial population-averaged activity of σRNN and RefCirc (see Sourmpis et al. (2023)). (iii) Finally, we add a regularization loss function ℒ_reg to penalize unnecessary large weights.

We also developed a biologically informed RNN model (bioRNN) for which we have designed a successful optimization technique. The main differences between σRNNs and bioRNNs consist in the following biological inductive biases. Firstly, the bioRNN neuron model follows a simplified leaky integrate and fire dynamics (see Methods 4.2) yielding strictly binary spiking activity. Secondly, we constrain the recurrent weight matrix to describe cell-type specific connectivity constraints: following Dale’s law, neurons have either non-negative or non-positive outgoing connections; moreover, since cortical inhibitory neurons rarely project across areas, we assume that inhibitory neurons project only locally within the same area. Thirdly, we add a term to the regularization loss ℒ _reg to implement the prior knowledge that cross-area connections are more sparse than within an area. Adding these biological features into the model requires an adapted gradient descent algorithm and matrix initialization strategies (Methods 4.5). The reconstruction method with σRNNs and bioRNNs is otherwise identical: the models have the same size, and are optimized on the same data, for the same number of steps and using the same loss functions. The two models bioRNN1 and bioRNN2 are optimized to explain recordings from RefCirc1 and Refcirc2, respectively. Importantly, the structural difference between RefCirc1 (feedforward) and RefCirc2 (feedback) is assumed to be unknown during parameter optimization: at initialization, excitatory neurons in bioRNN1 or bioRNN2 project to any neuron in the network with transmission efficacies (aka as synaptic weights) initialized randomly.

After parameter optimization, we have four models, σRNN1, σRNN2, bioRNN1 and bioRNN2, that we call “reconstructed” models. To validate the reconstructed models, we verify that the network trajectories closely match the data on the test set in terms of (i) the “behavioral” hit-trial frequency, (ii) the peristimulus time histogram (PSTH) mean-square error of single neurons as evaluated by ℒ _neuron, and (iii) the distance between single-trial network dynamics as evaluated by ℒ _trial (see Suppl. Fig. S2 and Suppl. Table 1). At first sight, the σRNN displays a better performance when comparing the activity simulated by the RNN and recordings held out in the testing set: ℒ _trial is for instance lower with σRNN (see Suppl. Table 1). This is expected considering that the optimization of bioRNNs is less flexible and efficient because of the sign-constrained weight matrix and the imperfect surrogate gradient approximation through spiking activity. However, the two bioRNNs are drastically more robust when evaluating the models with perturbation tests.

Perturbation test

To test which of the reconstructed RNNs capture the causal mechanisms of the RefCircs, we simulate optogenetic activations and inactivations of area B (step 3 in Fig. 1A). We first compare the change of hit probability after perturbations in the reconstructed RNN and recorded in RefCirc in Figure 2. For the σRNN the activation or inactivation of area B changes drastically the peak firing rate in area A: all trials become a hit during inactivation of area B. This drastic increase of hit rate is not consistent the reference where the effect of the optogenetic inactivations is mild: the distribution of network responses remains bi-modal (hit versus miss) with only a moderate change of hit frequency for RefCirc2. For RefCirc1 we even expect by design because of the absence of feedback connections from B to A. In contrast, the bioRNN models capture these changes more accurately (see Fig. 1 and 2). Quantitative results are summarized in Fig. 2C, the error of hit probability changes is 7% with bioRNNs when averaged over all conditions (bioRNN1 and bioRNN2, with optogenetic inactivations and activations). The corresponding error is 48.5% on average for σRNNs. In this sense, we argue that the bioRNN provides a better prediction of the perturbed hit frequency than the σRNN. We also performed spike train recordings in the area that is not directly targeted by the light to compare the perturbed network dynamics in the fitted RNNs and the RefCirc. The perturbed dynamics are displayed in Fig. 2B. The quantity is a distance between the network dynamics (RNN versus reference) and is reported in Fig. 2D-E and Supplementary Table 1. Again, the perturbed dynamics of the bioRNN are more similar to those of the reference circuits , than with the σRNN (t-test p-value is 0.0003).

Figure 2:

To analyze which features of bioRNN explain this robustness to perturbation, we then derive a family of models where only one feature of the reconstruction is omitted. Namely, the “No Dale’s law” model does not have excitatory and inhibitory weight constraints, the “Non-local inhibition” model allows inhibitory neurons to project outside of their areas, the “No Spike” model replaces the spiking dynamics with a sigmoidal neuron model, and the “No Sparsity” model omits the cross-area sparsity penalty in ℒ _reg. Omitting all these features in bioRNN would be equivalent to using a σRNN. The accuracy metrics on the testing sets before perturbation are reported for all RNN variants on the Suppl. Fig. 2C and F. For reference, we also include the model “No TM” (trial-matching) which omits the loss function ℒ _trial during training.

The strongest effect measured with this analysis is that the Dale’s law and local inhibition explain most of the improved robustness of bioRNNs. This is visible in Fig. 2 as the perturbed trajectories of “No Dale’s law” and “Non-local inhibition” are most distant from the reference in Fig. 2B. This is confirmed numerically where both the hit-rate error and the distance of network dynamics increase the most when lifting these constraints (Fig. 2C-E and Supplementary Table 1. We explain this result as follows: the mono-synaptic effect of a cell stimulated by the light are always correct in bioRNN (according to Dale’s law, and inhibition locality), but often wrong in the alternative models (see Fig. 2A). For instance, a simple explanation may justify the failure of the “Non-local inhibition” model: the stimulation of inhibitory neurons in B induces (via the erroneous mono-synaptic inhibition) a reduction in the baseline activity in area A (see the green trace during inactivation in Fig. 2B). More generally for perturbation testing, we speculate that these features are measured to be important here because they are central to the biophysical nature of the perturbation considered: optogenetic perturbation targets specific cell types, and these features incorporate a biophysical connectivity priors which is hard to infer entirely from the unperturbed data.

In comparison, it is less clear how spiking is relevant to model optogenetic perturbations. Considering that train spiking networks is a long-standing open question (Bohte et al., 2002; Neftci et al., 2019; Wunderlich and Pehle, 2021), it is technologically nontrivial however that the spiking bioRNN matches the performance of the “No Spike” model on both perturbation test metrics (hit-frequency in Figure 2C and distances of network dynamics ℒ _trial Figure 2E). The spiking bioRNN is expected to achieve worse predictions because surrogate gradients Neftci et al. (2019) through the binarized spiking dynamics is an approximation. Indeed the σRNN models achieves lower loss functions on unperturbed data (see Suppl. Table 1), but the fitting does not generalize to perturbations. On perturbed data, it is the spiking bioRNN which achieves slightly better performance (spiking is better on all metrics, but without significant margins, t-test p-value is 0.31 for, see Suppl. Table 1 and Fig. 2E). Also, we speculate that simulating spikes will have a stronger predictive advantage on other perturbation test experiments where timing matters more than long-lasting area-wide activations.

So far, we have identified biological inductive biases of bioRNN that improve its robustness to perturbation. But the prediction error in Fig. 2B on the perturbed trajectory is still too large for our method to unequivocally reveal the characteristic differences between the feedforward RefCirc1 and the recurrent RefCirc2. To reveal this difference in the optimized bioRNN1 and bioRNN2 models, we find that adding a sparse connectivity prior across areas is necessary. This error of the “No Sparsity” model is visible in Fig. 2B (red curves): the PSTH traces for both RefCirc1 and 2 (in Fig. 2B) display the early increase in area A characteristic of feedback from area B to A, it should not exist for RefCirc1. This is corrected with the sparsity prior (bioRNN, blue curves). Altogether on this dataset, the model with constrained weights, spiking activity, and sparsity prior predicted best the perturbed trials. It was also able to capture this characteristic difference between feedforward and recurrent circuits RefCirc1 and RefCirc2. We will now investigate whether the same approach can predict optogenetic perturbations with data recorded in vivo.

2.2 Predicting perturbations on in vivo electrophysiology data

To test whether our reconstruction method with biological inductive biases can predict optogenetic perturbations in large-scale recordings, we used the dataset from Esmaeili et al. (2021). In this study, water-deprived mice were trained to perform a whisker tactile detection task. In 50% of the trials (Go trials), a whisker is deflected, followed by a 1-second delay, after which an auditory cue signals the mice can lick a water spout to receive a water reward. In the other 50% of trials (No-Go trials), no whisker deflection occurs, and licking after the auditory cue results in a penalty with an extended time-out period. While the mice performed the task, experimenters recorded 6, 182 units from 12 areas across 18 mice. Using this dataset, we focused on the 6 most relevant areas for executing this task (Esmaeili et al., 2021). From each area, we randomly selected 250 neurons (200 putative excitatory and 50 putative inhibitory), which correspond to 1500 neurons in total. These areas, all shown to be causally involved in the task (Esmaeili et al., 2021), include the primary and secondary whisker sensory cortex (wS1, wS2), the primary and secondary whisker motor cortex (wM1, wM2), the anterior lateral motor cortex (ALM), and the primary tongue-jaw motor cortex (tjM1). We fit the neuronal activity using the same reconstruction method as used for the synthetic dataset. In the model, we simulate the jaw movement of the mouse as a linear readout driven by the model’s neural activity. This readout is regressed with the real jaw movement extracted from video footage. The parameter optimization of the behavioral readout is performed jointly with fitting the synaptic weights to the neuronal recordings, see Methods 4.5. After training, our reconstructed model can generate neural activity with firing rate distribution, trial-averaged activity, single-trial network dynamics and behavioral outcome which are all consistent with the recordings (see Suppl. Fig S5). Before perturbations, we observe again that the σRNN model fits the testing set data better than the bioRNN model (see Table 2 and Fig. S5).

Table 2:

We then submit the reconstructed σRNNs and bioRNNs models to perturbation tests. In the dataset, the behavior of the animal is tracked during inactivation of the an area at a given time window (stimulus, delay, or choice periods). For each of the six areas and time windows, we extract the averaged hit frequency under optogenetic inactivation, and attempt to predict this perturbed behavior by inducing the same inactivations to the fitted RNNs. These perturbations are acute spatiotemporal optogenetic inactivations of each area during different time periods (see Figure 3B). As an example, we show the effect of an inactivation of wS1 during the whisker period in the model in Fig. 3. In panel C, we display the simulated trial of a fitted bioRNN with and without perturbations side by side. The two trials are simulated with the same random seed, and this example shows that an early perturbation in wS1 can change a lick decision from hit to miss in the model (Fig. 3C).

Figure 3:

We denote by Δp^D the in vivo change in lick probability across Go trials in response to optogenetic pertubations. The perturbations were performed in different periods for each area in Esmaeili et al. (2021) (stimulation, delay, or choice periods). For all areas and time windows, we measure the corresponding in the model. On average, the error change probability obtained with the σRNN model is which is significantly worse than the bioRNN model’s 16% (t-test p-value is 0.014, see Figure 3D). As in the synthetic dataset, we find this to be consistent over multiple bioRNN model variant, and we find that imposing Dale’s law and local inhibition best explain the improvement in perturbation-robustness. We also measure that the spiking bioRNN predicts the change in lick probability slightly better than the “No Spike” bioRNN model. Conversely, adding the sparsity prior does not seem to improve the perturbed hit-rate prediction on the real data as seen in the recurrent artifical dataset (RefCirc2) and not in the feedforward case (RefCirc1) as shown in Suppl. Fig. S4.

To further analyze the consistency of the perturbations in the model, we can compare the perturbation map showing changes in lick probability obtained from acute inactivation in the data and the model. The Suppl. Fig. S6 summarizes visually which area has a critical role at specific time points. The change of lick probability in area wS1, ALM and tjM1 are accurately predicted by the bioRNN. In contrast, our model tends to underestimate the causal effect induced by the inactivations of wS2, wM1 and wM2 (Suppl Fig S6). Overall, our model is consistent with a causal chain of interaction from wS1 to ALM and continuing to tjM1.

2.3 Applications for experimental electrophysiology

With future progress in recording technology and reconstruction methods, network reconstruction may soon predict the effect of optogenetic perturbation with even higher accuracy. In this section, we explore possible consequences and applications for experimental electro-physiology. We demonstrate in the following that (1) perturbation-robust bioRNNs enable us to estimate gradients of the recorded circuits, (2) which in turn enable us to target µ-perturbations in the recorded circuit and optimally increase (or decrease) induced movements in our simulated mouse. The so-called “recorded circuit” is a bioRNN trained on the in vivo dataset that we use as a proxy experimental preparation. Its mathematical underpinnings enables us for rigorous theoretical considerations and the design of forward-looking in silico experiments.

µ-perturbations measure brain gradients

We first prove a mathematical relationship between gradients in the recorded circuit and µ-perturbations. We define the integrated movement as where y_t is the movement of the jaw movement at time t generated by the model, and we denote ΔY ⁷ as the change of movement caused by the µ-perturbation. If the circuit has well-defined gradients (e.g. say a “No spike” bioRNN model trained on the in vivo recordings in the previous section), using a Taylor expansion, we find that:

where ℐ are the neuron and time indices selected for the optogenetic intervention. The error term ϵ is negligible when the current induced by the light is small. We first confirm this approximation with numerical visualization in Fig. 4A: we display movement perturbations ⟨ΔY ↯⟩ in the circuit with time windows of decreasing sizes (⟨·⟩ indicates a trial average). When the time window is small, and the perturbation is only applied to excitatory or inhibitory cells in Fig. 4A, one can appreciate visually the similarity with the binned gradient in Fig. 4B. Proceeding to a quantitative verification of equation (1), we now compare the effect of small perturbations targeting only 20 neurons on a singletrial. We use the gradient (see Fig. 4C) to predict the outcome of µ-perturbation as follows: for each trial, and each 100ms time window, we identify 20 neurons in the model with highest (or lowest) gradients. We then re-simulate the exact same trial with identical random seed, but induce a µ-perturbation on selected neurons (see rectangles in Figure 4). If we target neurons with strongly positive gradients, the perturbed jaw movements are strongly amplified ΔY ↯ > 0; conversely, if we target neurons with negative gradients the jaw movements are suppressed ΔY ↯ < 0. Although the equation (1) is only rigorously valid for models with well-defined gradients like the “No Spike” model, we also confirm in Fig. 4D that this numerical verification also holds in a spiking circuit model where the gradients are replaced with surrogate gradients (Neftci et al., 2019).

Figure 4:

An implication of equation (1) i s that the measurements ⟨ΔY ↯⟩ that can be recorded in vivo are estimates of the gradients in the recorded circuit. Yet, measuring detailed gradient maps (or perturbation maps) as displayed in Fig. 4 would be costly in vivo as it requires to average ΔY ↯ over dozens of trials for each spatio-temporal window. Instead, gradient calculation in a bioRNN model (that was fitted to the experimental preparation) is a rapid mathematical exercise. If the extracted model is valid, then the gradients in the bioRNN approximate (1) the effect of µ-perturbations ΔY ⁷ in the experimental preparation; (2) the gradient in the recorded circuit.

Targeting in vivo µ-perturbations with bioRNN gradients

Building on this theoretical finding, we build a speculative experimental setup where the bioRNN gradients are used to target a µ-perturbation and increase (or decrease) the movements Y in the experimental preparation in real time. We show a schematic of this speculative close-loop experiment in Fig. 5C extending contemporary read-write elecotrophysiology setups (Packer et al., 2015; Adesnik and Abdeladim, 2021; Grosenick et al., 2015; Papagiakoumou et al., 2020). We demonstrate in silico in Fig. 5A-B how this experiment could use bioRNN gradients to bias the simulated mouse movement Y. As a preparation step, and before applying perturbations, we assume that the bioRNN is well fitted to the recorded circuit and we collect a large databank ℬ of simulated trials from the fitted bioRNN. Then in real-time, we record the activity from the experimental preparation until the time t^∗ at which the stimulation will be delivered (Step 1 in Fig. 5A, t^∗ is 100ms before the decision period). Rapidly, we find the trial with the closest spike trains in the databank of simulated trials (Step 2) and use the corresponding gradient maps to target neurons with the highest gradient in the model (Step 3). The targeted stimulation is then delivered immediately at t^∗ to the experimental preparation (Step 4). When testing this in silico on our artificial experimental preparation, we show in Fig. 5C that this approach can bias the quantity of jaw movement Y driven by the circuit in a predictable way. The amount of movement is consistently augmented if we target neurons with the highest (or diminished if we target neurons with the lowest).

Figure 5:

3 Discussion

Finding the right level of detail to model recorded phenomena has sparked intensive debates in computational neuroscience. Our results show that the co-development of perturbation tests and biological modeling provides a quantitative measurement to address this problem. Perturbation tests can characterize two types of incomplete brain models: (1) physiologically detailed models intended to explain brain mechanisms but do not enable quantitative predictions; (2) machine learning models with high predictive power but capturing a wrong causal mechanism causing erroneous generalizations.

Our reconstruction method uses data-constrained RNNs with informed inductive biases, and is validated on perturbation tests. We achieve a reconstruction of the sensory-motor pathway in the mouse cortex during a sensory detection task from electrophysiology data. The model is optimized to explain electrophysiological recordings and generalizes better than standard models to in vivo optogenetic interventions. We found unambiguously that anatomically informed sign and connectivity constraints for dominant excitatory and inhibitiory cell types improve the model robustness to optogenetic perturbations. In hindsight, we conclude that adding biological constraints becomes beneficial when (1) they model the interaction between the circuit and the perturbation mechanism; (2) their implementation should not impair the efficiency of the optimization process. Broadly speaking this hindsight is also supported by other results elsewhere in neuroscience. For instance: biologically inspired topological maps improve convolutional models of the Monkey’s visual system Schrimpf et al. (2024), and detailed cell-type distribution and connectome improve vision model in the fly brain Lappalainen et al. (2023); Cowley et al. (2024). For future work, there is a dense knowledge of unexploited physiological data at the connectivity, laminar or cell-type level that could be added to improve a cortical model like ours (Harris et al., 2019; Liu et al., 2022; Udvary et al., 2022; Staiger and Petersen, 2021; Rimehaug et al., 2023). By submitting the extended models to the relevant perturbation tests, it becomes possible to measure quantitatively the goodness of their biological mechanism implementations. We do not rule out, that significant improvements on perturbation tests can also be achieved with other means (e.g. by training deep learning architectures Azabou et al. (2024); Pandarinath et al. (2018); Ye et al. (2023) on larger datasets to enable generalization, or with generic regularization techniques like low-rank connectivity Dubreuil et al. (2022); Valente et al. (2022)). However, in a similar way as the σRNN was apriori an excellent predictor on our initial test-set, any powerful brain model will likely have failure modes that can be well characterized and measured with an appropriate perturbation test. So perturbation tests could become a central component of an iterative loop to identify needed data collection or model improvements towards robust brain models.

To highlight the importance of perturbation-robust circuit models, we have discussed possible implications for experimental neuroscience in section 2.3. We have shown that network reconstructions generate detailed perturbation sensitivity maps from non-perturbed recordings only. These sensitivity maps correspond to the RNN gradients and can be used to make targeted micro-stimulation in vivo. As a result we could design a hypothetical close-loop setup combining read-write electrophysiology with a brain model to influence the brain activity or behavior, having potentially important practical and ethical consequences. More conceptually, we have shown theoretically that the gradients of a perturbation robust RNN are also consistent with the gradients of the true recorded circuits. In perspective with the foundational role of gradients in machine learning theory LeCun et al. (2015); Richards and Kording (2023), it enables the measurement of “brain gradients” and lays a computational link between in vivo experimental research and decades of theoretical results on artificial learning and cognition.

4 Methods

4.1 Mathematical toy model of the difficult causal inference between H1 and H2

Let’s consider two simplistic mathematical models that both depend on two binary random variables A and B which represent that putative area A is active as A = 1 and area B as B = 1. With this notation, we can construct two hypothetical causal mechanisms H1 (“a feedforward hypothesis”) and H2 (“a recurrent hypothesis”), which are radically different. The empirical frequency p(A, B) of the outcome does not allow us to differentiate whether the system was generated by a feedforward mechanism H1 or a recurrent mechanism H2. Schematically, we can represent the two mechanism hypotheses as follows:

For hypothesis H1: we assume that external inputs are driving the activity of area A such that A = 1 is active with probability p₀, and there are strong feed-forward connections from A to B causing systemically B = 1 as soon as A = 1. Alternatively, in H2, we assume that areas A and B receive independent external inputs with probability . Each of these two inputs is sufficient to cause A = 1 or B = 1, and the two areas are also strongly connected, so A = 1 always causes B = 1 and vice versa. Under these hypothetical mechanisms H1 and H2, one finds that the empirical probability table p(A, B) is identical: (“Hit trial”, both areas are active), p(A = 0, B = 0) = 1 − p₀ (“Miss trial”, the areas are quiescent) (To prove this: a and b to denote the binary external inputs into A and B, so we have: where we used that p(A = 1, B = 1|a, b) is 0 or 1, then using p(a = 1) = p(b = 1) = p1 and the independence between a and b we find: ). In both cases, the possibility that only one area is active is excluded by construction. So for any A and B p_H1(A, B) = p_H2(A, B) and in other words, even if we observe an infinite number of trials and compute any statistics of the binary activations A and B, discriminating the two possible causal interactions (H1 versus H2) is impossible.

A solution to discriminate between hypotheses H1 and H2 is to induce a causal perturbation. We can discriminate between our two hypotheses if we can impose a perturbation that forces the inactivation of area B in both mathematical models. In mathematical terms we refer to the do operator from causality theory. Under the feedforward mechanism H1 and inactivation of B, A is not affected p_H1 (A = 1 | do (B = 0)) = p₀. Under the recurrent hypothesis, H2, and inactivation of B, A is activated only by its external input such that p_H2 (A = 1 | do (B = 0)) = p1 ≠ p₀. So the measurement of the frequency of activation of area A under inactivation of B can discriminate between H1 and H2 which illustrates mathematically how a causal perturbation can be decisive to discriminate between those two hypothetical mechanisms.

4.2 Neuron and jaw movement model

We model neurons as leaky-integrate and fire (LIF) neurons. The output of every neuron j at time t is a binary outcome (spike if, no spike if) generated from its membrane voltage. The following equations give the dynamics of the membrane voltage :

where , and are the recurrent and input weight matrices. The timestep of the simulation δt is 2 ms when we simulate the real dataset and 1 ms otherwise. The superscript d denotes the synaptic delay; every synapse has one synaptic delay of either 2 or 3 ms. With, we define the integration speed of the membrane voltage, where τ_m = 30 ms for excitatory and τ_m = 10 ms for inhibitory neurons. The noise source is a Gaussian random variable with zero mean and standard deviation is typically initialized at 0.14). The input is a binary pulse signal with a duration of 10 ms. For the real dataset, we have two binary pulse input signals, one for the whisker deflection and one for the auditory cue. The spikes are sampled with a Bernoulli distribution , where v₀ is the temperature of the exponential function and v_thr,j is the effective membrane threshold. After each spike, the neuron receives a reset current with an amplitude of v_trh,j and enters an absolute refractory period of 4 ms, during which it cannot fire.

For networks fitted to the real dataset, we also simulate the jaw movement. The jaw movement trace y is controlled by a linear readout from the spiking activity of all excitatory neurons. Specifically, y is computed as, where b is a scaling parameter and is given by . Here, is the output weight matrix (linear readout) for the jaw, and τ_jaw = 5ms defines, which controls the integration velocity of the jaw trace.

4.3 Session-stitching and network structure

As in (Sourmpis et al., 2023), we simulate multi-area cortical neuronal activity fitted to electrophysiology neural recordings. Before we start the optimization, we define and fix each neuron’s area and cell type in the model by uniquely assigning them to a neuron from the recordings. For the real dataset from Esmaeili et al. (2021), the cell type is inferred from the cell’s action potential waveform (with fast-spiking neurons classified as inhibitory and regular-spiking neurons as excitatory). Most electrophysiology datasets include recordings from multiple sessions, and our method would typically require simultaneous recordings of all neurons. To address this challenge, similarly to (Sourmpis et al., 2023) we use the technique called “session-stitching” which allows neighboring modeled neurons to be mapped with neurons recorded across multiple sessions. This effectively creates a “collage” that integrates data from multiple sessions within our model. This approach has practical implications for our optimization process. Specifically, the trial-matching loss includes a term for each session, with the overall loss calculated as the average across all sessions (see 4.5).

For both the real and the synthetic datasets, we simulate each area with 250 LIF neurons and impose that each area has 200 excitatory neurons and 50 inhibitory. Respecting the observation that inhibitory neurons mostly project in the area that they belong to (Tamamaki and Tomioka, 2010; Markram et al., 2004), we don’t allow for across-area inhibitory connections. The “thalamic” input is available to every neuron of the circuit, and the “motor” output for the real dataset, i.e., jaw movement, is extracted with a trained linear readout from all the excitatory neurons of the network, see 4.2.

4.4 Reference circuits for hypotheses 1 and 2

To build a synthetic dataset that illustrates the difficulty of separating the feedforward (H1) and recurrent hypotheses (H2), we construct two reference spiking circuit models RefCirc1 and RefCirc2. The two networks consists of two areas A and B, and their activity follows the hard causal inference problem from method 4.1, making it hard to distinguish A1 and A2 when recording the co-activation of A and B. Moreover, to make the problem even harder, the two networks are constructed to make it almost impossible to distinguish between H1 and H2 with dense recordings: the two circuits are designed to have the same PSTH and single-trial network dynamics despite their structural difference, one is feedforward and the other is recurrent.

To do so, RefCirc1 and 2 are circuit models that start from random network initializations following the specifications described in Methods 4.2 and 4.3. The only difference is that we do not allow feedback connections from A to B in RefCirc1, the construction below is otherwise identical. The synaptic weights of the two circuits are optimized with the losses described in Methods 4.5 to fit the identical target statistics in all areas: the same PSTH activity for each neuron and the same distribution of single-trial network dynamics. The target statistics are chosen so the activity in RefCirc1 and 2 resemble kinematics and statistics from a primary and a secondary sensory area. The baseline firing rates of the neurons is dictated by the target PSTH distribution and it follows a log-normal distribution, with excitatory neurons having a mean of 2.9 Hz and a standard deviation of 1.25 Hz and inhibitory neurons having a mean of 4.47 Hz and a standard deviation of 1.31 Hz. The distribution of single-trial activity is given by the targeted single-trial dynamics: in RefCirc1 and 2, the areas A and B respond to input 50% of the time with a transient population average response following a double exponential kernel characterized by τ_rise = 5 ms and τ_fall = 20 ms. Mimicking a short signal propagation between areas, these transients have a 4 ms delay in area A and 12 ms delay in B (relative to the onset time of the stimulus). To impose a “behavioral” hit versus miss distribution that could emerge from a feedforward and recurrent hypothesis (see method 4.1), the targeted population-averaged response of each trial is either a double-exponential transient in both area A and B (”Hit trials”), or remains at a baseline level in both areas (”Miss trials”) in the remaining trials. At the end of the training, we verified that RefCirc1 and RefCirc2 generate very similar network activity in the absence of perturbation (see Figure S1). The circuits are then frozen and used to generate the synthetic dataset. We generate 2000 trials from these RefCircs, 1000 of which are used for the training set and 1000 for the testing set.

4.5 Optimization and loss function

The optimization method we use to fit our models is back-propagation through time (BPTT). To overcome the non-differentiability of the spiking function, we use surrogate gradients (Neftci et al., 2019). In particular, we use the piece-wise linear surrogate derivative from Bellec et al. (2018b). For the derivative calculations, we use and not exp. We use sample-and-measure loss functions that rely on summary statistics, as in (Bellec et al., 2021; Sourmpis et al., 2023), to fit the networks to the data. Our loss function has two main terms: one to fit the trial-averaged activity of every neuron (ℒ _neuron), and one to fit the single trial population average activity (ℒ _trial), = ℒ _neuron + ℒ _trial. The two terms of the loss function are reweighted with a parameter-free multi-task method (Défossez et al., 2023) that enables the gradients to have comparable scales.

As in Sourmpis et al. (2023): (1) To calculate the trial-averaged loss, we first filter the trial-averaged spiking activity using a rolling average window (f) of 8 ms. We then normalize it by the trial-averaged filtered data activity, (z^D are recorded spike trains)

where ⟨. ⟩ _t is the time average, and σ_t the standard deviation over time. The trial-averaged loss function is defined as:

where T is the number of time points in a trial and N is the number of neurons. For the real dataset, where we want to fit also the jaw movement, we have an additional term for the trial-averaged filtered and normalized jaw, , where y is the simulated jaw movement and y^D the recorded jaw movement.

(2)To calculate the trial-matching loss, we first filter the population-average activity of each area , using a rolling average window of 32 ms. We then normalize it by the population-averaged filtered activity of the same area from the recordings, , and concatenate all the areas that were simultaneously recorded, , where ⟨.⟩_k is the trial average, and σ_k the standard deviation over trials. The trial-matching loss is defined as:

where π is an assignment between pairs of K recorded and generated trials π : {1, … K} → {1, … K}. Note that the minimum over π is a combinatorial optimization that needs to be calculated for every evaluation of the loss function. For the real dataset, we consider the jaw movement as an additional area, and we concatenate it to the .

Based on this loss function, we optimize the following parameters: , and β for the RefCircs. For the RNNs, we optimize only the recurrent connectivity W ^rec,d, and the rest are fixed from the RefCircs. For the real dataset, additionally to the parameters optimized in the RefCircs, we also optimize the jaw’s linear readout and its scaling parameter b.

Implementing Dale’s law and local inhibition

In our network, the recurrent weights W ^rec are computed as the elementwise product of two matrices:, which encodes the strength of synaptic efficacies and is always positive, and, which has a fixed sign determined by the neurotransmitter type of the presynaptic neuron and :

To enforce Dale’s law during optimization, we set any negative values of to zero at each iteration as in Bellec et al. (2018a). Similarly, to constrain only local inhibitory connections during optimization, we zero out any changes in the synaptic efficacies of across-areas inhibitory connections at each iteration. In simplified models, Dale’s law or the local inhibition constraint can be disrupted by omitting this correction step.

The success of the network optimization highly depends on the initialization of the recurrent weight matrices. To initialize signed matrices we follow the theoretical Rajan and Abbott (2006) and practical insights Bellec et al. (2018b); Cornford et al. (2020) developed previously. After defining the constraints on the weight signs, the initialization amplitude for each target neuron is adjusted to a zero-summed input weights (the sum of incoming excitatory inputs is equal to the sum of inhibitory inputs). Then the weight amplitude is re-normalized by the modulus of its largest eigenvalue of W ^rec, so all the eigenvalues of this matrix W ^rec have modulus 1 or smaller.

Stopping criterion for the optimization

For the synthetic dataset, we train the models for 4000 gradient descent steps. For the real dataset, due to limited data and a noisy test set, we select the final model based on the optimization step that yields the best trial-type accuracy (closest to the trial-type accuracy from the data), derived from the jaw trace and whisker stimulus, along with the highest trial-matched Pearson correlation between the model and the recordings.

Sparsity regularization

There is a plethora of ways to enforce sparsity. In this work, we use weight regularization. In particular, we use the norm of the recurrent and input weights that promote a high level of sparsity (Xu et al., 2012). To avoid numerical instabilities, we apply this regularization only for synaptic weights above α and prune all synapses below α. (we set α = 1e⁻⁷). The regularized loss function becomes:

where W ^across,d are the connections from one area to the other.

For the synthetic dataset, we choose the level of across-area sparsity by performing a small grid search for λ₃. In particular, the sparsity level λ₃ is the maximum value λ₃ where the performance remains as good as without sparsity, see Suppl. Fig S3. For the real dataset, we use the same value λ₃ as the one we found for the full reconstruction method of bioRNN1.

4.6 Perturbation test of in silico optogenetics

In systems neuroscience, a method to test causal interactions between brain regions uses spatially and temporally precise optogenetic activations or inactivations (Esmaeili et al., 2021; Guo et al., 2014). Usually, inactivations refer to the strong activation of inhibitory neurons for cortical areas. These inhibitory neurons have strong intra-area connections that effectively “silence” their local-neighborhood (Helmstaedter et al., 2009).

Our model can simulate these perturbations and allow us to compare the causal mechanisms of two networks based on their responses to optogenetic perturbations. We implement activations and inactivations as a strong input current to all the neurons in one area’s excitatory or inhibitory population. For the RefCircs and reconstructed RNNs, we use a transient current that lasts 40 ms, from 20 ms before to 20 ms after the input stimulus. The strength of the current (light power) varies until there is an effect in the full reconstruction method bioRNN1. For the synthetic dataset in Figure 2 (except for panel D), we inject a current of into excitatory neurons for activations and into inhibitory neurons for inactivations. For the real dataset, we perform optogenetics inactivations in three different periods. As in Esmaeili et al. (2021), we silence the cortical circuit during the whisker presentation, the time between the whisker and auditory stimulus, or when the animal was licking for the reward. In particular, we use transient currents to the inhibitory neurons during (i.) 100 ms before and after the whisker presentation, (ii.) 100 ms after the whisker presentation till 100ms before the onset of the auditory cue, and (iii.) after the auditory cue till the end of our simulation. For cases (i.) and (ii.), we linearly decreased the strength of the current to avoid rebound excitation. The light power is chosen so that our model has the best results in reproducing the lick probability of the recordings. It is important to mention that the perturbation data are not used to optimize the network but to test whether the resulting network has the same causal interactions with the recordings.

For the RefCircs and bioRNNs, we evaluate the effect of the perturbations directly from the neural activity. We use the distance of network dynamics ℒ _trial to compare the two perturbed networks. For the real dataset, we compare the effect of the inactivations on the behavior; as behavior here, we mean whether the mouse/model licked. We classify the licking action using a multilayer perceptron with two hidden layers with 128 neurons each. The classifier is trained with the jaw movement of the real dataset, which was extracted from video filming using custom software, to predict the lick action, which was extracted from a piezo sensor placed in the spout. This classifier predicted lick correctly 94% of the time. We then used the same classifier on the jaw movement from the model to determine whether there was a “lick” or not.

Data availability statement

The code for this project is open sourced and published at https://github.com/Sourmpis/BiologicallyInformed. The dataset for the artificial dataset can be downloaded/generated on our code repository. The in vivo dataset was published openly for the previous publication Esmaeili et al. (2021). The dataset is accessible at: https://zenodo.org/records/4720013.

Supplemental information

Table 1:

Figure S1:

Figure S2:

Figure S3:

Figure S4:

Figure S5:

Figure S6:

Acknowledgements

We thank Alireza Modirshanechi, Shuqi Wang, and Tâm Nguyen for their valuable feedback on the manuscript. We are grateful to Vahid Esmaeili for collecting the dataset and ongoing support throughout this project. This research is supported by the Sinergia project CR-SII5 198612, the Swiss National Science Foundation (SNSF) project 200020 207426 awarded to WG, SNSF projects TMAG-3 209271 and 31003A 182010 awarded to CP, and the Vienna Science and Technology Fund (WWTF) project VRG24-018 awarded to GB.