Automated customization of large-scale spiking network models to neuronal population activity

Abstract

Understanding brain function is facilitated by constructing computational models that accurately reproduce aspects of brain activity. Networks of spiking neurons capture the underlying biophysics of neuronal circuits, yet their activity's dependence on model parameters is notoriously complex. As a result, heuristic methods have been used to configure spiking network models, which can lead to an inability to discover activity regimes complex enough to match large-scale neuronal recordings. Here we propose an automatic procedure, Spiking Network Optimization using Population Statistics (SNOPS), to customize spiking network models that reproduce the population-wide covariability of large-scale neuronal recordings. We first confirmed that SNOPS accurately recovers simulated neural activity statistics. Then, we applied SNOPS to recordings in macaque visual and prefrontal cortices and discovered previously unknown limitations of spiking network models. Taken together, SNOPS can guide the development of network models, thereby enabling deeper insight into how networks of neurons give rise to brain function.

Fig. 1 |. Framework for automated customization of a spiking network model to neuronal recordings.

a, A SNN has a complicated dependency between its parameters and spiking output. For example, different parameter sets correspond to each of four previously identified activity regimes of a CBN: asynchronous irregular, synchronous regular, synchronous irregular and asynchronous regular. In this work, the SNN has eight parameters, including those that govern the connection strength between neurons as well as the timescale of synaptic decay (Supplementary Table 1). b, Our customization framework matches activity statistics of spike trains produced by the network model to those of neuronal recordings. It uses a guided searching algorithm to iteratively update the model parameters. The activity statistics are defined by the user and can include single-neuron, pairwise and population activity statistics.

Fig. 2 |. Activity statistics for comparing the activity of a spiking network model to neuronal recordings.

a, The three types of activity statistics based on single neurons (for example, fr and ff), pairs of neurons (for example, spike count correlation) and a population of neurons (for example, percent shared variance, number of dimensions and eigenspectrum of shared variance). The units of fr are spikes per second, those of ff are spike count and those of the eigenspectrum of shared variance are (spike count). All other activity statistics are unitless. These activity statistics are all based on spike counts within a 200 ms spike count bin (left), which can be represented in a population activity space (center). Each dot represents the activity across the neuronal population within a given time window. b, The activity statistics based on population recordings in macaque visual area V4 (dashed lines) were challenging to reproduce by the four parameter regimes of a CBN (colored symbols, cf. Fig. 1a, mean across five network instantiations of network connectivity graphs and initial membrane potentials corresponding to the same network parameter set). None of the four activity regimes accurately reproduced the activity statistics of the V4 population recordings (dashed lines). The V4 activity statistics are shown as the mean ± 1 s.d. across 19 recording sessions (Methods). The spike counts for V4 were computed using a fixed 200 ms time window preceding the onset of each stimulus presentation. All activity statistics were based on randomly subsampling 50 neurons from each CBN or V4 dataset.

Fig. 3 |. Customizing a spiking network model using BO with GPs.

The BO algorithm attempts to find a parameter set $\mathbf{\theta }$ for a spiking network model such that its activity statistics match those of neuronal recordings. a, Spike trains are recorded from the brain and their activity statistics, ${\mathbf{s}}_{\text{true}}$, are computed. This step is performed only once, since the same recorded activity is used for comparison on all iterations. b, On the $k$-th iteration, spike trains are generated from the network model using parameter set ${\mathbf{\theta }}_k$, proposed by the previous iteration. c, The activity statistics of the spike trains generated from the network model, $\mathbf{s}\left({\mathbf{\theta }}_k\right)$, are computed. The cost for ${\mathbf{\theta }}_k$ depends on how far each of those activity statistics is from the corresponding activity statistics of the neuronal recordings, ${\mathbf{s}}_{\text{true}}$. The intensification procedure occurs between b and c. The feasibility of ${\mathbf{\theta }}_k$ is determined here using a short simulation. d, A GP (solid line) is used to approximate the true, unknown cost function, $c\left({\mathbf{\theta }}_k\right)$ (dashed red line). We seek to find the minimum of this true, unknown cost function (denoted by ${\mathbf{\theta }}^{\star }$). Each iteration of the BO provides one evaluation of the cost at a particular setting of the model parameters (black dots). The cost at the current iteration is labeled $\hat{c}\left({\mathbf{\theta }}_k\right)$ and the other black dots represent the costs evaluated during previous iterations. The GP provides an uncertainty of our estimate of the cost function (gray shading). For illustrative purposes, we show here a single model parameter being optimized, whereas our algorithm typically optimizes multiple model parameters simultaneously. A separate GP is used to approximate the feasibility of ${\mathbf{\theta }}_k$ (not shown). e, An acquisition function is defined based on the two GPs in d to determine the next parameter set, ${\mathbf{\theta }}_{k+1}$, to evaluate. The acquisition function implements an exploration–exploitation tradeoff, where areas of low predicted cost and high uncertainty are desirable.

Fig. 4 |. Accurate customization of a CBN model to simulated spike trains using SNOPS.

a, A CBN was used to generate spike trains with randomly chosen parameter sets ${\mathbf{\theta }}_{\text{true}}$ (Methods). SNOPS (or other optimization algorithms) was then used to customize the parameters, ${\mathit{\theta }}_k$, of a separate CBN to match the ‘ground truth’ activity statistics, ${\mathbf{s}}_{\text{true}}$, of the generated spike trains. b, For a given amount of computer running time (Methods), SNOPS (blue) finds parameters with lower cost than accelerated random search (red) and random search (green). The vertical axis represents the lowest log(cost) up to the given running time and hence decreases monotonically. The solid lines and shading represent the mean ± 1 s.d. across 40 customization runs. Note that 10⁵ s equals 1.2 days or 27.8 hours. c, For a representative customization run, SNOPS (blue) identified model parameters whose activity statistics were closer to the ground truth (dashed lines) than accelerated random search (red) and random search (green). The error bars on the ground truth represent one s.d. across five network instantiations corresponding to the same ground truth parameter set. The circles represent the mean across five network instantiations corresponding to the network parameter set identified by each optimization algorithm. d, Across all 40 customization runs, SNOPS accurately reproduced the ground truth activity statistics (all points lie near the diagonal). Each dot represents the results from one SNOPS customization run to a randomly generated ground truth dataset. For visual clarity, only the first (that is, most dominant) mode of es is plotted in the rightmost image.

Fig. 5 |. Reproducing activity statistics of macaque V4 and PFC recordings with the CBN and SBN.

a, Left: stylized representation of the CBN. Right: activity statistics of the CBN (circles, mean across five network instantiations corresponding to the same identified parameter set) after being customized using SNOPS to the same V4 dataset as in Fig. 2b. The dashed line and shading represent the mean ± 1 s.d. across 19 sessions. b, Left: stylized representation of the SBN. The SBN is different from the CBN in that the connection probability depends on the distance between neurons. Right: activity statistics of the SBN (triangles, mean across five network instantiations corresponding to the same identified parameter set) after being customized using SNOPS to the same V4 datasets as in a. c, The SBN more accurately reproduced activity statistics than the CBN across 16 datasets, comprising four task conditions with recordings in two brain areas (V4 and PFC) in each of the two monkeys. The arrow indicates the example V4 dataset shown in a and b.

Fig. 6 |. Revealing the inflexibility of CBN relative to SBN with tradeoff cost.

a, Activity statistics of the CBN (circles, mean across five network instantiations corresponding to the same identified parameter set) and SBN (triangles, mean across five network instantiations corresponding to the same identified parameter set) after being customized using SNOPS to one V4 activity statistic (dashed line) at a time. The same V4 dataset as Fig. 2b was used. b, A high tradeoff cost represents the case where customizing the network to reproduce two activity statistics simultaneously yields a higher average cost of the two statistics than customizing each statistic individually (top). By contrast, a low tradeoff cost represents the case where the cost of customizing two activity statistics simultaneously yields a similar cost to customizing each statistic individually (bottom). c, Tradeoff costs between pairs of statistics for the CBN (top) and SBN (bottom) on the same V4 dataset as Fig. 2b. We observed similar effects when customizing these models to PFC recordings (Supplementary Fig. 9). d, Customizing the CBN and SBN to different numbers of activity statistics included in the cost function simultaneously, on the same V4 dataset as Fig. 2b. Each dot represents one particular subset of activity statistics (for example, highlighted dot indicates the average cost of $r_{\text{sc}}$ and ff when including only those two activity statistics in the cost function). The cost of each dot was computed over five network instantiations corresponding to the same identified parameter set of that dot. Each bar indicates the average cost across all subsets of the corresponding number of activity statistics.