Discovering sparse control strategies in C. elegans

Edward D. Lee1*, Xiaowen Chen2, Bryan C. Daniels3

1 Complexity Science Hub Vienna, Josefstaedterstrasse 39, Vienna 1080, Austria

2 Laboratoire de physique de l’École normale supérieure, CNRS, PSL Université,
Sorbonne Université, Université de Paris, 75005 Paris, France

3 School of Complex Adaptive Systems, Arizona State University, Tempe, AZ
85287, USA

E.D.L. designed the research. E.D.L. and X.C. wrote the code. All authors
analyzed the results and edited the manuscript.

*<EMAIL_ADDRESS>
## Abstract

Biological circuits such as neural or gene regulation networks use internal
states to map sensory input to an adaptive repertoire of behavior.
Characterizing this mapping is a major challenge for systems biology, and
though experiments that probe internal states are developing rapidly,
organismal complexity presents a fundamental obstacle given the many possible
ways internal states could map to behavior. Using C. elegans as an example, we
propose a protocol for systematic perturbation of neural states that limits
experimental complexity but still characterizes collective aspects of the
neural-behavioral map. We consider experimentally motivated small
perturbations — ones that are most likely to preserve natural dynamics and are
closer to internal control mechanisms — to neural states and their impact on
collective neural behavior. Then, we connect such perturbations to the local
information geometry of collective statistics, which can be fully
characterized using pairwise perturbations. Applying the protocol to a minimal
model of C. elegans neural activity, we find that collective neural statistics
are most sensitive to a few principal perturbative modes. Dominant eigenvalues
decay initially as a power law, unveiling a hierarchy that arises from
variation in individual neural activity and pairwise interactions. Highest-
ranking modes tend to be dominated by a few, “pivotal” neurons that account
for most of the system’s sensitivity, suggesting a sparse mechanism for
control of collective behavior.

## Author summary

The relationship between underlying biological circuitry and behavior is
complex and difficult to probe experimentally. Part of the problem is that,
for organisms of even modest size, there are an overwhelming number of
possible combinations of interacting circuit components. We develop a
theoretical framework to simplify this problem with experiments that change
the system minimally with small perturbations. In the realm of small
perturbations, not only does system behavior remain close to normal but the
range of possible perturbations is greatly reduced to only pairs of
experimental targets. We demonstrate such a perturbation using a minimal model
of neural activity in the C. elegans worm. We find that a few combinations of
“pivotal” neurons strongly affect the statistics of synchronous activity,
suggesting they may be important for neural control of behavior. Our work
suggests feasible, perturbative experiments to map how the physical components
of an organism control emergent collective activity.

## Introduction

Control of complex dynamical networks is a problem of major interest in
biology for finding drivers of diseased states [1, 2, 3] and determinants of
behavior [4, 5, 6, 7]. In neural systems, the question of which and how many
neurons correspond to controllability is largely open. While there is evidence
that single-neuron manipulation is sufficient to induce behavioral change in
some cases [8, 9], other research shows behavioral information to be encoded
amongst many neurons [10, 11, 4]. Identifying interesting neurons to probe
experimentally and choosing how to perturb them is difficult: in principle, a
thorough experiment would require a combinatorially large number of procedures
to test all the possible ways that neural targets could be modified.
Theoretical tools provide a way of winnowing down the number using control-
theoretic analysis of dynamical systems models [12, 13], structural analysis
[14, 5, 15], and network properties [16], amongst other approaches [17]. Here,
we develop a theoretical framework that explicitly considers perturbation
experiments to discover candidate control neurons. Our framework suggests a
way of leveraging recent developments in optogenetics [18, 19] that allow
fast, precise control of specific neruons and even the regulation of brain
states with the closed-loop optogenetic [20, 21, 22, 23] or optoelectronic
systems [24].

The question of control is in essence about how neural states map to behavior.
When the mapping is not deterministic, this is a problem of stochastic
encoding of the behavioral state $Z$ in the global neural state $Y$. Whether
from noise, experimental limitations, or other sources of uncertainty, there
is for a given neural state $Y$ a multiplicity of $Z$ characterized by the
conditional probability distribution $P(Z|Y)$. In practice, we cannot specify
a unique neural state $Y$ due to experimental uncertainty or limits to
measurement. Instead, we measure some distribution over configurations $P(Y)$
that we then map to some probability distribution of behavioral outcomes
$P(Z)$. By considering the realm of possible distinct neural configurations,
we obtain the stochastic mapping formulation in Figure 1, where at each point
we specify a distribution of neural activity, of behavior, and the
relationship between the two that accounts for imprecision and uncertainty.

Fig 1: (A) Local curvature of a manifold of probability distributions $P(Z)$
is captured by the Fisher information matrix (FIM). FIM eigenvectors $v_{i}$
represent axes of combinations of neural perturbations. Eigenvalues
$\lambda_{i}$ denote strength of curvature, or inverse sensitivity. (B) Set of
possible behavioral distributions $P(Z)$ as a high-dimensional manifold. (C)
Neural behavior represented as ensemble $P(Y)$ over states of activity $Y$.
Relationship between $P(Y)$ and behavioral distribution $P(Z)$ is given by
stochastic encoding $P(Z|Y)$ that maps neural activity to behavior.

A classic formulation of stochastic mapping in sensory encoding is that of
individual neurons that encode increasingly complex patterns of sensory
stimuli as information is funneled up a hierarchy [25, 26]. In this case, the
mapping between $Y$ and $Z$ is precise. For example, we take $Y$ to be the
firing rate of some neuron and $Z$ to be the similarity of the stimulus to
one’s grandmother. Since firing rate is a mean of a noisy measurement, it is
more appropriate to interpret the average firing rate to characterize a
distribution $P(Y)$, e.g. a Poisson distribution for the number of times a
neuron fires within a time frame. Then, varying the mean activity level of a
single neuron traces out a curve in the space of possible distributions
$P(Y)$. In the space of possible perceptions $P(Z)$, representing for example
a probability distribution over recognizing “grandma” or not and given by
$P(Z)=\sum_{Y}P(Z|Y)P(Y)$, we likewise trace out a curve as a single neural
parameter is varied.

Later paradigms generalized the focus on single neurons to a sparse set of
neurons that together span the range of possible sensory inputs. In the
formulation of sparse coding, where a basis set of neurons are extracted from
measuring neural response to a set of visual stimuli, leading to abstract
image features such as edges, corners, and contours [27, 28]. The idea of
sparse coding is that typical visual scenes can be represented with activity
in only a few neurons, with each neuron representing a separate high-level
feature. Then, it is the joint activity amongst these distributed components
that captures the whole image, an idea with multiple variations in the
literature [29, 30].

Despite the generalization away from single neurons in the study of sensory
encoding [31], an echo of it still rings in experimental studies of neural
control of behavior. Some experiments rely on the notion of individual (or
sets of similar) “driver” neurons that — when ablated, silenced, or otherwise
modified — substantially change behavioral outcomes [32, 33, 34]. While some
individual cells, like interneurons connecting distal parts of the body, are
essential for normal behavior, some higher-order behaviors depend on many
neurons in a way that is robust to the function of individual cells [35, 36,
37]. Indeed, the finding that some neural circuits encode sensory input in a
sparse, distributed way raises the question of whether neural control of
behavior also follows similar principles [26, 38, 4, 29].

Fig 2: Four possibilities for the sparsity of collective modes of sensitivity.
When collective activity shows sloppy structure, the local information
geometry is elongated, with both stiff sensitive directions and sloppy
insensitive ones. With dense collective encoding, multiple components each
matter equally. The starred combination, with sloppy, sparse combinations of
neurons, aligns with centralized control.

Here, we explore two potential sources of sparsity in neural control of
behavior. First, a desired output could be equally well produced by many
possible changes, or modes, at the neural level or could be much more easily
produced by only one particular change. Second, each mode can itself be sparse
or dense in the number of neurons involved [39, 40]. These possibilities lead
to the four hypotheses we draw in Figure 2, where control is sloppy or non-
sloppy with either sparse or dense collective modes of sensitivity.

In the context of distributions of neural activity, variation in sensitivity
of neural activity is described by the local dependence of the neural activity
distribution $P(Y)$ on parameters specifying its form. As pictured in the top
of Figure 1, “stiff” directions correspond to large changes in $P(Y)$ being
caused by small changes to neural parameters [41, 42]. In the “sloppy”
directions, changes to neural parameters are ineffectual, and only dramatic
perturbations cause a comparable change in $P(Y)$. The relative elongation of
the local information geometry disappears when encoding places no special
importance on any neural subgroup such that each plays a commensurate role, as
in some forms of population coding [43, 44]. In this way, the local
information geometry informs us about the first notion of sparsity mentioned
above. Additionally, the second notion of sparsity is captured by the degree
to which stiff directions involve changes to only a few neurons. Of particular
interest is the combination of sparse and sloppy structure, where collective
sensitivity would be concentrated in only a few perturbative sets that each
include only a few neurons. Such a possibility reflects centralized control,
which would seem to be optimal for simplifying the multiplicity and complexity
of control nodes. If real world examples display a range of structure
depending on organism and function, we need a flexible methodology for
inferring such variation from data.

We propose a perturbative, experimental approach to study stochastic encoding
between internal states and output behavior. Perturbing the statistics of
neural activity parameterizes motion in the space of all possible
distributions, the rate of which quantifies sensitivity to neural changes. Our
focus on small perturbations allows us to derive a local approximation that is
the simplest complete description of local sensitivity.

For concreteness, we focus on the model animal Caenorhabdits elegans. With
only 302 neurons and structural connectome mapped, C. elegans presents a
particularly interesting example with which to explore neural control. Current
experiments measure whole-brain activity in a freely moving worm and soon will
allow matching those neurons to the structural connectome [45, 46]. Combined
with tools for manipulating neural states [47], our proposed experimental
protocol could be used to search for neural centers important to collective
activity. As an example of the proposed method, we start with minimal models
matching the statistics of neural activity in C. elegans, and then compute the
local information geometry to extract combinations of neurons that, if
perturbed, would strongly affect collective statistics. Though we focus on the
neural activity of C. elegans, our procedure could be adapted to other
measures of internal states such as hormone concentration, gene expression, or
other experimentally accessible quantities, and to output behavior including
morphological descriptors such as orientation, velocity, extension, or a
stereotyped behavioral sequence [48]. In this sense, our work provides a
generalizable theoretical framework with specific measurements and predictions
that can be tested experimentally in C. elegans and other biological systems.

## A perturbative thought experiment

An exhaustive protocol for probing the C. elegans neural system would be to
test simultaneous perturbations to subsets of neurons in ways that are
sympathetic, antagonistic, and of varying magnitude. Unfortunately, the
combinatorial explosion in the number of possibilities makes such an
experiment impossible. In the limit of small perturbations, however, all
possible effects can be extracted from observing only the response of the
system to perturbations of pairs of neurons. We use this limit in part because
pairwise relationships between perturbations becomes a sufficient description
of system response. Additionally, it is often the case that collective
statistics of neural activity can be captured with pairwise neural statistics
without needing to invoke higher-order correlations [39, 49, 50], so we might
anticipate that the pairwise approximation could be a reasonable one even for
larger changes. We demonstrate how this simplification could be useful by
carrying out a thought experiment that mimics the perturbations in silico,
making use of a minimal model of neural activity.

To formalize the set of perturbations, we consider the discretized gradient of
calcium concentration $s_{i}$ of neuron $i$, which is either decreasing
$s_{i}=-1$, flat $s_{i}=0$, or increasing $s_{i}=1$ with probabilities $p(s)$
of observing the set $s\equiv\\{s_{i}\\}$. Such discretization has been
explored in previous work [10]. Then, a perturbation to underlying neural
activity is reflected in a modified probability distribution over discrete
states $\tilde{p}(s)=p(s)+\Delta(s)$. A unique measure of distinguishability,
the information divergence, between the original distribution and the modified
one reduces to

$\displaystyle D_{\rm KL}[p||\tilde{p}]$
$\displaystyle\approx\frac{1}{2}\sum_{s}\frac{\Delta(s)^{2}}{p(s)}=\frac{1}{2}\sum_{s}\sum_{ij}F_{ij}dv_{i}(s)dv_{j}(s).$
(1)

Then, the complete set of perturbations is given by the Hessian $F_{ij}$
describing the local sensitivity of $p(s)$ to perturbations along vectors
$dv$, which may be constrained by what is accessible in the experiment or
model. Because moving from one basis to another is a linear operation, the
particular form of the perturbations are not as important as the fact that the
set should span the possible set of perturbations. By exploiting this property
of analyticity, we can in principle reduce the complexity of the problem by
measuring a single set of experimental perturbations.

In a hypothetical experiment, we might effect perturbations by inserting
electrodes into immobilized worms — or, in a more elegant protocol using
optogenetic tools, by clamping membrane potential and effectively coupling
neurons to one another via an external circuit to control the strength and
timing of perturbation. Here, we use a perturbation that with small
probability copies the state of one neuron to another, chosen to resemble
increase or decrease in synaptic strength. Specifically, we select pairs,
identifying a “target” neuron in state $s_{t}$ and a “matcher” neuron in state
$s_{m}$, and clamp the matcher to the target’s state, $s_{m}=s_{t}$, with
small probability $\epsilon$ (see Appendix G for more details). Compared to
more common clamping and ablation techniques [23, 51], this proposed closed-
loop perturbation is not as drastic, mostly preserves natural dynamics, and is
closer to internal control mechanisms such as synaptic modulation that do not
require constant input from the outside [52, 53, 54].

## Maximum entropy (maxent) model example

We build a minimal model of the collective statistics of C. elegans on which
to run such a protocol in silico, systematically perturbing it and extracting
from its response the sensitive modes of collective neural statistics. The two
examples we consider are of $N=50$ neuron subsets from the anterior neural
network of the immobilized worm [4]. With a discretized time series
representing the three possible configurations for neuron $m$, we calculate
time-averaged probabilities of being in state $k$,

$\displaystyle r_{k}(s_{m})=\frac{1}{T}\sum_{t=1}^{T}\delta_{s_{m},k}(t),$ (2)

where the Kronecker delta function $\delta_{s_{m},k}$ indicates when the state
of the neuron $s_{m}$ is $k$ at time $t$ over the duration of the experiment
$T$. Similarly, the pairwise probabilities of agreement between any two
neurons is $\sum_{t=1}^{T}\delta_{s_{m},s_{t}}(t)/T$ and similar higher-order
correlations describe the probabilities of agreement between multiple neurons.
These higher-order correlations are not necessarily given by the lower-order
statistics, but only accounting for the pairwise correlations is sufficient to
capture the higher-order correlations [50] (see Figures A.6 and A.7).

An advantage of the maxent approach is that it captures the statistics of
neural activity in a way that reflects latent physical interactions while also
obeying a quantitative formulation of Occam’s Razor [55, 56]. The maxent
principle ensures that a model of the probability distribution $p(s)$ only
matches specified constraints but otherwise remains as structureless as
possible. This can be done by maximizing the model’s information entropy,
$S=-\sum_{s}p(s)\log p(s)$, here a sum over all $3^{N}$ possible
configurations, with the standard method of Langrangian multipliers [57, 58,
59, 39]. When the set of average individual neural activity
$\\{\left<r_{k}(s_{m})\right>\\}$ is constrained, the maxent procedure gives
the “independent model” where the probability of finding neuron $i$ in state
$k$,

$\displaystyle r_{k}(s_{i})$ $\displaystyle=e^{h_{i,k}}/Z_{i},$ (3)

is determined by the bias $h_{i,k}$, with normalization term
$Z_{i}=\sum_{k=-1}^{1}e^{h_{i,k}}$. A larger bias $h_{i,k}$ relative to
another state $k^{\prime}$ indicates that it is more likely to find the neuron
in state $k$ than $k^{\prime}$. This independent model fails to capture
correlations in neural activity.

A simple variation that does not assume independence involves constraining the
pairwise correlations. Given the finite-size of the data set, we constrain the
pairwise coincidence probabilities to obtain the pairwise maxent model with
distribution given the energy function $E(s)$,

$\displaystyle\begin{aligned} p(s)&=\left.\exp\left[-E(s)\right]\right/Z,\\\
Z&=\sum_{s}\exp\left[-E(s)\right],\\\
E(s)&=-\sum_{i<j}^{N}J_{ij}\delta_{s_{i},s_{j}}-\sum_{k=-1}^{1}\sum_{m=1}^{N}\delta_{k,s_{i}}h_{i,k}.\end{aligned}$
(4)

The log-probability of configuration $s$ varies with its energy such that
lower energy means larger probability. As with the independent model,
increasing neural bias $\\{h_{i,k}\\}$ pushes neurons to state $k$ and
increasing the magnitude of couplings $\\{J_{ij}\\}$ magnifies the tendency of
pairs to agree or to disagree when negative. Though the couplings are in
principle exactly specified by the pairwise correlations in the data, sampling
and experimental noise mean that there are many possible sets of couplings
that align within the statistical variation of the data sample. We focus on a
numerical solution that recovers topological structure in the coupling network
and has been shown to faithfully capture collective neural patterns [50] (see
Appendix A for solutions from an approximate gradient descent method). The
solution then defines a minimal interaction network that is distinct from the
pattern of pairwise correlations (see Figure A.6).

Fig 3: (A) Perturbation thought experiment consists of clamping matcher neuron
$m$ to the state of target neuron $t$ with some small probability $\epsilon$
when indicated by a random number generator (RNG). We picture electrodes
controlling membrane voltage, but optogenetic protocols may be more elegant.
(B) Perturbation corresponds to modifying fields and couplings in a pairwise
maxent model. (C) Principal eigenmatrix mapped to change in (D) biases
$\\{h_{i}\\}$, each direction indicated by a different color, and (E)
couplings $\\{J_{ij}\\}$ for $\epsilon=10^{-4}$. (F) Diffuse perturbation in
space of observables using replacement rule from Eq 5 corresponds to (G)
localized “natural” perturbation to one coupling (note the nonzero values in
the top left corner). Localized coupling perturbation is nontrivial and
nonintuitive in terms of the observable perturbation.

## Mapping perturbations from experiment to simulation

With the maxent model, we enact our in silico thought experiment. Continuing
with the formulation presented above, we connect a pair of target neuron $t$
and matcher neuron $m$ such that $s_{m}$ is clamped to match $s_{t}$ with
small probability $\epsilon\ll 1$. Then, the modified probability
$\tilde{r}_{k}$ that the matcher neuron is in state $k$ is the mixture

$\displaystyle\tilde{r}_{k}(s_{m})$
$\displaystyle\equiv(1-\epsilon)r_{k}(s_{m})+\epsilon\,r_{k}(s_{t}).$ (5)

As a result of the perturbation in Eq 5, the statistics describing neuron $m$
becomes more like those of neuron $t$. This procedure likewise modifies the
matcher neuron’s coincidence $p(s_{m}=s_{j})\rightarrow\tilde{p}(s_{m}=s_{j})$
with all other neurons indexed $j$ as if modifying synaptic weights [60]. If
we were to take sufficiently large $\epsilon$ in experiment, we would expect
such intervention not to remain localized to neuron $s_{m}$ but to alter the
neighbors and eventually the neighbors of neighbors and so on. If $\epsilon$
is small, then we expect a localized perturbation to be able to well-
approximate the desired change in statistics given in Eq 5. By taking the
limit where the perturbation is sufficiently small and assuming that it then
remain localized, we specify a perturbation that corresponds to a unique
change of parameters.

Such uniqueness implicitly depends on constraining ourselves to considering
perturbations that move us along the family of pairwise maxent models. This
constraint is not reflected in the replacement rule in Eq 5, which compatible
with an infinite variety of changes to higher-order correlations. In other
words, we do not necessarily need to restrict ourselves to considering only
the energy function defined in Eq 4, but we could allow for higher-order
interactions to appear under perturbation. This introduces ambiguity that can
only be resolved by choosing the form of the perturbations. Once we have
chosen to fix the structure of the energy function from the maxent principle,
however, we do not allow a perturbation to arbitrarily alter the form. This
assumption is consistent with the widespread observation that the pairwise
model generally captures well collective features of biological neural
networks of modest size [39, 61, 62, 63, 64, 65]. Thus, we use the pairwise
maxent model not only to specify a minimal, compressed representation of
neural statistics, but also to specify how the probability distribution
evolves under the replacement rule specified in Eq 5.

This replacement rule maps experimental perturbations to observed statistics
in contrast with starting with model parameters as is often the first instinct
for theorists. In the latter formulation, the fields and couplings would form
the basis for perturbations, also known as “natural” or “canonical”
perturbations [66]. This latter perspective adheres to the physical intuition
that the parameters from Eq 4 reflect latent physical interactions. When the
pairwise maxent model instead serves as an approximation of the underlying
physical interactions, a natural perturbation may be mediated by unknown
parameters instead of fields and couplings [49, 67], and this renders the
straightforward model perturbation difficult to translate into a literal
experiment. In comparison, defining perturbations in observable terms as in Eq
5 allows for easier experimental interpretation and for consistent comparison
of inferred models [60].

If it is collective neural activity that encodes behaviorally relevant
information [40, 39], then we can use collective properties as a proxy for
behavior. As one example measure of collective activity, we consider the
probability that the number of neurons in each of the three states, $\phi_{\rm
fine}(n_{1},n_{2},n_{3})$, is described when ordered by $n_{1}\geq n_{2}\geq
n_{3}$ such that $n_{1}$ corresponds to the size of the plurality and $n_{3}$
the smallest minority (i.e. there is no fixed correspondence to “rise,”
“fall,” and “flat”). This presents a measure of synchrony more fine-grained
than $\phi_{\rm coarse}(n_{1})$ that only considers the number of neurons in
the plurality $n$. Crucially, the pairwise maxent model also captures these
statistics. While neither statistic differentiates between the orientation of
the states, we limit ourselves here to simpler and less computationally
expensive statistics of collective neural activity [68, 69], noting that other
statistics may easily be incorporated in the way we outline.

Under perturbation of a pair of matcher and target neurons, the synchrony
distribution $\phi$ is mapped to an altered $\tilde{\phi}$ that depends on the
strength of perturbation $\epsilon$. To measure the distance between the two
distributions, we use the Kullback-Leibler (KL) divergence $D_{\rm KL}$ as in
Eq 1 [66]. In the limit of infinitesimal perturbation, $\epsilon\rightarrow
0$, KL divergence reduces to the Fisher information

$\displaystyle\begin{aligned}
F_{mt,m^{\prime}t^{\prime}}&\equiv\lim_{\epsilon\rightarrow
0}\frac{2}{\epsilon^{2}}D_{\rm
KL}\left[\phi\Big{|}\Big{|}\tilde{\phi}\right]=\lim_{\epsilon\rightarrow
0}\frac{2}{\epsilon^{2}}\sum_{n_{1}\geq n_{2}\geq
n_{3}}\phi\log\left(\frac{\phi}{\tilde{\phi}}\right)\equiv\lim_{\epsilon\rightarrow
0}\frac{2}{\epsilon^{2}}\left<\log\left(\frac{\phi}{\tilde{\phi}}\right)\right>\\\
&=\sum_{\theta_{i},\theta_{j}}\left<\frac{\partial^{2}\log\phi}{\partial\theta_{i}\theta_{j}}\mathcal{J}_{mt,m^{\prime}t^{\prime}}^{ij}\right>,\end{aligned}$
(6)

where the Jacobian $\mathcal{J}_{mt,m^{\prime}t^{\prime}}^{ij}$ transforms the
maxent parameters such as fields and couplings $\\{\theta_{i}\\}$, here
ordered along a single index $i$, into the vector of changes that correspond
to the observable perturbations specified in Eq 5. In contrast with Eq 1, we
consider in Eq 6 the impact of perturbation on collective synchrony in Eq 6
that accounts for the multiplicity of microscopic configurations belonging to
a single coarse-grained state [70].

The Fisher information matrix (FIM), whose entries are defined in Eq 6, is
spanned by eigenvectors that describe orthogonal perturbations of the
parameters $\\{\theta_{i}\\}$, where modes with large eigenvalues describe
perturbations to which collective synchrony is highly sensitive and small
eigenvalues represent ones to which the system is insensitive [60, 42]. As we
show in Eq 6, this basis describes the local curvature of a Riemannian
manifold generated by information distance. Importantly, the basis vectors can
involve a mix of antagonistic and sympathetic perturbations of neurons of
varying magnitude, one that would be nontrivial to extract a priori. Thus, the
FIM encodes how quickly coarse-grained configurations $\phi$ change as we
modify the system on a microscopic level by perturbing pairs of neurons at a
time, with perturbations that mimic changes to synaptic connectivity.

## Leading FIM eigenvalues show Zipfian decay

In Figure 4A, we show the rank-ordered eigenvalue spectrum of the FIM for
collective synchrony $\phi_{\rm fine}$ calculated with the pairwise model in
comparison with several null models: independent neurons, pairwise maxent
model with couplings randomly shuffled between all pairs of neurons (which
preserves the distribution of couplings but not the topology of the
interaction network), and canonical perturbations directly modifying each
coupling at a time by a fixed amount (Figure C.10). Across all cases, we
expect a hard cutoff at the dimensionality of the synchrony distribution
$\phi_{\rm fine}$ corresponding to rank $Z_{\rm max}=234$, which is much
smaller than a single dimension of the matrix $N^{2}-N={2{,}450}$ ($Z_{\rm
max}=33$ for the coarse-grained collective statistic $\phi_{\rm coarse}$ as we
show in Figure E.19). In contrast with the others, the independent neuron
model has short maximum cutoff well below the theoretical maximum, reflecting
the essential role of interactions in mediating pairwise perturbations. The
cutoff for the pairwise maxent model reveals the replacement rule in Eq 5
sufficient to very nearly span the full dimensionality of synchrony space.
This observation confirms that the pairwise model with the pairwise
replacement rule generates a useful basis to explore collective sensitivity.

The rank-ordered spectrum of eigenvalues $\lambda_{z}$ of rank $z$ initially
decays such that each successive level of perturbation returns
multiplicatively smaller response, following roughly on a very limited range
Zipf’s law, $\lambda_{z}\propto z^{-1}$. In contrast, a simple exponential
decay would indicate a sharp cutoff for sensitive modes beyond some rank. We
find that exponential decay alone does not describe the eigenvalue spectrum,
and instead find a reasonable fit using a power law with an exponential tail:

$\displaystyle\lambda_{z}$ $\displaystyle=Az^{-\alpha}e^{-z/\bar{z}};\qquad
z\leq Z_{\rm max}$ (7)

In Eq 7, we have parameters for vertical scaling $A$, exponent $\alpha$, tail
$\bar{z}$, and hard cutoff $Z_{\\-\rm max}$ beyond which there are no
additional eigenvalues or numerical precision errors may be important (see
Section C for more details about the fitting procedure). For random subsets of
$N=50$ neurons considered, we find that this model usually presents a
compelling fit and that the scale of the cutoff $\bar{z}\gtrsim 30$ indicates
a power-law-like regime for the dominant eigenvalues (Figures G.22 and G.23).

Such a scaling regime seems to be a feature of the system statistics,
consistently appearing across the null models. The shuffled null model agrees
closely with our thought experiment, in which eigenvalue decay tends to be
slowest, showing a power-law exponent closest to $\alpha=1$ (see Figure C.10).
Perhaps unsurprisingly, the spectra are also scaled to higher sensitivities
than the independent model. Though the latter is roughly commensurate in
overall sensitivities with canonical perturbations, we might expect
perturbations in the space of observable and model parameters to be scaled
differently from the transformation of variables (Appendix E). Thus, this
statistical hierarchy reflects the role of component disorder common across
the models whose overall collective sensitivity is magnified by interactions.

Fig 4: (A) Eigenvalue spectrum of FIM for pairwise maxent model, independent
model (indpt.), and shuffled couplings null, with observable perturbations
from Eq 5 averaged over $M$ Monte Carlo samples ($M=10$ for maxent model,
$M=10$ over random shuffles for coupling shuffled, $M=4$ for independent and
canonical). For comparison, we show response to canonical perturbation to
couplings. Insets on left show example eigenmatrices of rank 1 and rank 50.
Inset on right shows full eigenvalue spectrum. Error bars show standard error
of the mean. (B) Uniformity from Eq 8 averaged over samples for eigenmatrices.
Error bars represent standard deviation over rank-ordered plots. (C) Principal
eigenvalues averaged over random subsets of neurons as a function of subset
size rescaled by the average maximum eigenvalue for the maxent model.
Sensitivity tends to be much lower for the independent model but of similar
scale for shuffled couplings. Error bars represent standard deviation taken
over subsets for each permutation or model. Compare with Figure E.18.

## FIM eigenvectors reveal pivotal neurons

Inspecting the corresponding eigenvectors, we find some modes are dominated by
perturbations focused only on a few neurons. To better represent the
connection between pairwise perturbations and eigenvectors, we reshape
eigenvectors into eigenmatrices $V_{mt}$ such that the elements in a column
indicate how matcher neuron $m$ should imitate its target neighbors $t$ in
turn. When put into this representation, we often find vertical striations as
in the inset of Figure 4A. These striations are visible because they are
almost exclusively of the same sign, indicating that the mode describes
perturbations localized to a single neuron that tend to increase or decrease
its correlation with all neighbors simultaneously. In contrast, horizontal
striations that represent uniform perturbations across all the neighbors of a
particular neuron, a kind of global perturbation, tend to be sparser and
weaker on average. This emergent pattern contained in the block structure of
the FIM suggests that localized, uniform enhancement or suppression of
synaptic connections leading to a small set of pivotal neurons may serve as
effective mechanisms for modulating collective activity.

As a more direct analysis of such key neurons, we limit our analysis to
perturbations focused on a single matcher neuron at a time and put aside
perturbations combining different matcher neurons. When the FIM is ordered
first by matchers and then targets for each matcher, these perturbations
correspond to blocks along the diagonal of the FIM of dimension
$(N-1)\times(N-1)$ for fixed $m$ and variable $t$. From the principal
eigenvalues of the diagonal blocks, we find that the single, pivotal neurons
with the largest eigenvalues tend to coincide with the ones that manifest in
vertical striations of the full FIM (Appendix C). These striations correspond
to columns with high uniformity. As a measure of this, we define row and
column uniformity, respectively, as

$\displaystyle\begin{aligned}
U_{i}&\equiv\left(\sum_{j}v_{ij}\right)^{2};\qquad
V_{j}&\equiv\left(\sum_{i}v_{ij}\right)^{2},\end{aligned}$ (8)

for each normalized eigenmatrix $v_{ij}$. When we consider the subspace of
leading pivotal neurons and compare them with the subspace of randomly chosen
neurons, the principal eigenvalues of the former set saturate much more
quickly to the leading eigenvalue of the full matrix even when we have only
considered about half of the available neurons (Figure 4C), indicating that
there is a subset that in combination overwhelmingly determines collective
sensitivity. This concentration of sensitive modes in a few neurons is an
indication of centralized structure that emerges from neural heterogeneity.

The patterns that we note both in the eigenvalue spectrum and eigenvector
basis hold generally for random neuron subsets of $N=50$ sampled from each
experiment, about the maximum possible number that can be sampled reliably
[50]. In contrast with the typical notion that the identities of particular
neurons are essential, pivotal neurons fluctuate between subsets and random
samples as shown in Figure 5. While some neurons are more frequently
identified in spite of random fluctuations for the same subset of $N=50$
neurons, consistently identified pivotal neurons are less pronounced once we
consider different subsets. This means that the properties we find are not
tied to individual neurons, but are general features of the ensemble. Since
collective synchrony is a lower-dimensional statistic and in principle permits
exchange symmetries between neurons, this is not necessarily surprising. On
the other hand, we note that these collective properties are robust to
heterogeneity in network connectivity, bias, and interactions that would tend
to break such symmetries.

Fig 5: Fraction of times out of 10 Monte Carlo samples that neuron column
uniformity is above 99% and 99.9% percentile cutoffs. Distribution is column
uniformities over all neurons over all eigenmatrices for a single Monte Carlo
sample. See also Figure G.21.

## Discussion

Examples of how sensory and behavioral information is encoded in neural
activity suggest that, while information is sometimes localized to a few
neurons, it is at other times distributed amongst many. This may result from
information flow between collective to individual components [71] and because
the precise scale at which information is processed may vary with time,
function, and organism. We develop a perturbative approach that is sensitive
to the potential variety of involved scales. Using an information-geometric
perspective, we identify statistical aspects of the distribution of neural
behavior that are essential for preserving collective activity. Importantly,
we do not have to assume that collective properties will be sensitive to
individual neurons or certain combinations but can discover the appropriate
ones in a principled way. This becomes feasible to explore comprehensively
because the response of the system depends on only pairs of perturbations when
they are sufficiently small, a property of analyticity of the mapping from
activity to collective state.

As an in silico realization of the protocol, we use perturbations that mimic
internal neural coordination and calculate their impact on collective
synchrony. In a minimal model, we find that dominant perturbative modes do not
tend to be distributed amongst many neurons but are localized into pivotal
ones. The concentration of sensitivity in a few neurons is analogous to the
presence of driver neurons, modulation of which drives the system from one
configuration to another [72, 73, 74]. Interestingly, we find that pivotal
neuron identities fluctuate across subsets and Monte Carlo random samples that
introduce finite sample errors into the entries of the FIM. Instead of
specific pivotal neurons, it is the collective properties of localized
eigenvectors and scaling in sensitivity that are preserved.

These collective properties require neural heterogeneity. Since identical
neurons would imply uniform eigenvectors, variation in each neuron’s local
interaction network and bias is responsible for the emergence of pivotal
neurons [60]. By shuffling the inferred interaction network and the data, we
verify that these features seem to be a result of the distribution and
magnitude of interactions but not the specific topology. In this sense, the
neurons do not need to be labeled differently as if given some predetermined
role, but the collection of local interactions and biases is responsible for
the statistical properties of the FIM, echoing findings in theoretical models
of neural networks [75]. Though we do not address this question here, it would
be interesting to know how such properties emerge and if such ensemble
properties are encoded genetically, arise spontaneously, or generalize to
mobile worms.

One of major questions that we approach in our framework is how experimental
perturbations should be represented in the model. While the focus is often on
model parameters as representations of physical dials that are experimentally
accessible, direct correspondence is not obvious for a statistical model
inferred from data. As an example, when couplings inferred by a pairwise
maxent model were directly compared with physical contact between protein
residues, the model recovered only a subset of real interactions even while
recovering the statistical ensemble [55, 76]. Taking the pairwise maxent
model, it is clear exactly what increasing a coupling does to the energy
function, enhancing the tendency of a pair of neurons to coincide, but the
actual result is nontrivial modification of the entire distribution. For the
experimentalist, it seems more natural to consider the problem from the
perspective of the observable statistics that can be perturbed in a controlled
and precise manner. In the case of Boltzmann-type models, the relationship
between observables and model parameters can be made exact and results in a
simple form for small perturbations. Thus, our formulation is a theoretically
and experimentally tractable framework for predicting the effects of
perturbations.

To test these ideas requires innovations on top of previous experiments
(Section G). A central feature of our thought experiment is the small
perturbation. If the strength of the intervention ranges from perturbative to
discrete, the perturbative limit lies close to no intervention at all and
discrete interventions correspond to silencing or ablation. Discrete
interventions may be appropriate considering the logical structure of a
circuit but become combinatorially expensive to explore fully [77, 78]. While
avoiding these issues, small perturbations require a sufficient number of
measurements for the perturbed parameters to be measurable. By using the fact
that the Fisher information is proportional to the number of independent
samples and the Cramér-Rao bound, a lower bound on the number of experimental
observations $T$ required to distinguish a change in the mode for a
perturbation strength $\epsilon$ for Fisher information $F$ is on the order of

$\displaystyle T\sim(\epsilon^{2}F)^{-1}.$ (9)

Taking a relatively large perturbation of $\epsilon=10^{-2}$, we have $T\sim
10$ independent samples for an eigenvalue $\lambda\sim 10^{3}$, the order of
magnitude of the largest mode [79]. Linear growth in the number of required
samples, however, restricts us from measuring insensitive modes in a
reasonable amount of time (the experiments analyzed here last 8 minutes and
collect roughly 80 to 120 independent samples [50]). Though mapping the full
local information geometry of $N=50$ neurons would be difficult because the
number of pairwise perturbations exceeds $\sim 10^{6}$ (accounting two
distinct pairs of matcher and target neurons), a similar experiment with about
$N\sim 10$ neurons seems reasonable when coupled with techniques for
incomplete matrix estimation [80]. The feasibility of such an experiment stems
from our argument that perturbative experiments harness analyticity to vastly
simplify the range of possible perturbations. We exploit this property to
extract a basis for the local information geometry including multi-component
perturbations. One goal of such experimental intervention is then not to be
more precise but sufficiently varied to span the local basis, an idea that can
be generalized to other experimental systems besides the C. elegans model we
consider.

Until now, we have assumed that the stochastic encoding relating neural to
collective activity does not change, but this is not necessarily the case. To
consider this, we characterize the encoding with the fundamental notion of
channel capacity, the maximum rate at which information can conveyed by such a
mapping given by the mutual information
$I[Y;Z]=\sum_{Y,Z}P(Y,Z)\log\left(P(Y,Z)/P(Y)P(Z)\right)$ [81]. A targeted
perturbation to the set of states $Y$ corresponds to a new distribution
$\tilde{P}(Y)=P(Y)[1+a_{r}v_{r}(Y)]$ giving arbitrary weight $a_{r}$ to
perturbation vector $v_{r}$ of index $r$ [82]. The perturbation vector $v_{r}$
may refer to a projection from a single-neuron perturbation or eigenvector
considered above. Generally, such a perturbation may change not only the
distribution of neural statistics $P(Y)$ but also the mapping
$\tilde{P}(Z|Y)=P(Z|Y)[1+\Delta_{Z|Y}]$, which corresponds to an adaptive code
that changes in response to incoming statistics. Under small perturbations
$a_{r}\ll 1$ and small adaptation $\Delta_{Z|Y}$, the change in the mutual
information decomposes into

$\displaystyle\tilde{I}[Y;Z]-I[Y;Z]$
$\displaystyle=\underbrace{a_{r}\sum_{Y}P(Y)v_{r}(Y)\sum_{Z}d(Z|Y)}_{\rm
direct}+\underbrace{\sum_{Y}P(Y)\sum_{Z}\Delta_{Z|Y}d(Z|Y)}_{\rm adaptive},$
(10)

where we have defined $d(Z|Y)\equiv P(Z|Y)\log[P(Z|Y)/P(Z)]$ such that
$\sum_{Z}d(Z|Y)$ is the information gained about $Z$ from observing $Y$. Eq 10
contrasts the “direct” result of a perturbation in contrast with an “adaptive”
response by the system. Whereas the former term decomposes into the product of
the local information geometry of $P(Y)$ and the information content of the
mapping $Y\rightarrow Z$, the latter depends on system response that modifies
the stochastic mapping itself, such as rate limiting or rescaled response
[83]. Deviations from the predictions in our perturbative thought experiment
might represent the effects of such complications from the properties of the
stochastic encoding.

The information geometry of scientific theories more generally suggests that
they show sparse structure, where a few parameter dimensions strongly change
the qualitative characteristics of phenomena and most dimensions are
unimportant. This is a result of the logarithmically spaced eigenvalues of the
FIM, or “sloppiness” [84, 85, 42, 86]. Whether by nature or by design, such
quantitative reduction vastly simplifies the level of detail required for
approximate theories, allowing for accurate prediction even with large
uncertainty in most parameters [87, 41]. We find here a variation on this idea
in the statistics of neural activity in C. elegans. Neural activity seems non-
sloppy, with largest eigenvalues decaying slower like Zipf’s law in contrast
with regular logarithmic spacing when we consider experimentally realistic
perturbations. Thus, sensitivity, while concentrated into a few dominant
modes, has non-negligible smaller modes that decay in a self-similar way.
Might this be a feature of biological neural networks to reduce the
dimensionality of control parameters (if interestingly not as strongly as the
sloppy case) or to nest varying levels of control? The state-of-the-art today
with single neuron measurement in C. elegans and optical genetics is
approaching a point where experiments to test this hypothesis may become
feasible.

## Acknowledgments

We would like to thank the UNM Center for Advanced Research Computing,
supported in part by the National Science Foundation, for providing high
performance computing resources used in this work. We thank Matthew Fricke for
invaluable assistance with these resources. We thank the Santa Fe Institute
for providing computational resources used for this work. The computational
results presented have been achieved in part using the Vienna Scientific
Cluster (VSC). E.D.L. acknowledges funding from the Omega Miller Program, NSF
grant PHY1838420, and Medical University of Vienna. X.C. acknowledges
Francesco Randi for insightful discussion. B.C.D. was supported by the ASU–SFI
Center for Biosocial Complex Systems.

##  Appendix A Numerical solutions to inverse maxent problem

In principle, the maxent formulation presents a unique mapping from
statistical correlations to model parameters such that there is no parameter
fitting in the usual sense. The problem of determining the fields $h_{i,k}$
and couplings $J_{ij}$ that closely match mean activity and pairwise
correlations is known as the inverse maxent problem [88, 89]. In practice,
however, limitations to numerical precision and finite-sample noise mean that
fitting the parameters for even moderately sized systems is not without
ambiguity. With this ambiguity in mind, we present the different approaches we
use to solve the inverse problem for the pairwise maxent and independent
models.

The first method, which we focus on in the main text, allows us to identify a
sparse interaction network between neurons mirroring the sparse structural
connectivity of the neural connectome [50]. The parameters of the maxent model
are initialized at zero and are updated individually as to maximally increase
of the likelihood of the data at each step [90, 91]. This learning procedure
is stopped once the model reproduces the constrained observables within
experimental variability, estimated by bootstrapping random halves of the
data. If experimental variability is relatively large, this training procedure
can return a sparse interaction matrix $J_{ij}$. We show the results of such a
procedure in Figures A.6 and A.7.

The second method we use for comparison relies on the Monte Carlo Histogram
(MCH) approach [91], which is an approximate gradient descent algorithm. At
each step, the parameters are simultaneously updated by an amount that is
proportional to the difference between the corresponding observable calculated
from the model and the data. We obtain a solution that is within a norm error
threshold

$\displaystyle\sqrt{\sum_{i=1}^{N}(\left<s_{i}\right>_{\rm
data}-\left<s_{i}\right>)^{2}+\sum_{i<j}^{N}(\left<s_{i}s_{j}\right>_{\rm
data}-\left<s_{i}s_{j}\right>)^{2}}<1/2,$ (11)

where the cutoff is arbitrarily set to obtain a relatively fast and close fit
to the statistics of the data. Unlike the first method, MCH returns a dense
network of connections since all couplings are updated at every iteration til
convergence. While we find that the resulting model aligns well with the
features of the data, including collective synchrony, this procedure does not
make realistic assumptions about the topological structure of the underlying
physical network. For these solutions, we again find a signal for
distinguishing pivotal neurons in the uniformity of columns and rows, but it
is not as consistent as evident in the column and row uniformities plotted in
Figure C.14.

We also consider an independent neuron model. When solving the corresponding
maxent model, we penalize large fields by maximizing the log-likelihood
$\log\mathcal{L}_{i}$ for each spin $i$ along with a penalty such that the
total cost function $\mathcal{C}_{i}$ is defined as

$\displaystyle\mathcal{C}_{i}$
$\displaystyle\equiv-\log\mathcal{L}_{i}+\frac{1}{\sigma}\sum_{k=-1}^{1}h_{i,k}^{2}$
(12) because large fields make it especially costly to calculate the FIM
accurately. Indeed, a sparse cost function, one where the cost scales with the
absolute value of the fields, does not sufficiently penalize large fields,
rendering it infeasible for our current implementation of the FIM calculation.
With the constraint given in Eq 12, the goal is then to find the fields such
that $\displaystyle h_{i,k}^{*}$
$\displaystyle=\min_{h_{i,k}}\mathcal{C}_{i}.$ (13)

To determine the weight $\sigma$, we compute the cost function
$\mathcal{C}_{i}$ in Eq 12 for a range of $\sigma$. At small $\sigma$, the
quadratic penalty dominates and the function approaches a large constant. At
large $\sigma$, we recover the original maximum likelihood problem. We set
$\sigma=\sigma^{*}$ as is determined by the midpoint between these two
extremes for the cost function averaged over all spins $i$. Typically, this is
in the interval $\sigma^{*}\in[2,10]$. These steps return an independent model
of neurons, where the biases of the most extreme neurons have been tempered.

Fig A.6: Pairwise maxent model of C. elegans anterior neural activity for
experiment 170419 from reference [4]. (A) Pairwise correlations between subset
of $N=50$ neurons. Average individual neuron $r_{k=-1}(s_{i})$ shown along
diagonal. (B) Inferred biases $h_{i,k=-1}$ along diagonal and interaction
matrix of couplings $J_{ij}$ off the diagonal. (C) Collective synchrony,
probability that a plurality of $n$ neurons coincide $\phi_{\rm coarse}(n)$.
Inset compares pairwise correlations for data, pairwise maxent model,
independent model, and one example of shuffled couplings. See reference [50]
for more details about the model. See Figure A.7 for experiment 154139. Error
bars show one standard deviation over bootstrapped samples. Fig A.7: Maxent
model fit overview for experiment 154139. See Figure A.6 for experiment 170419
and explanation of panels.

##  Appendix B Calculation of Fisher information matrix (FIM)

To calculate the entries of the FIM, we rely on a Monte Carlo Markov Chain
(MCMC) sample of the distribution of neural states $p(s)$, which we then
coarse grain to approximate the distribution over collective synchrony $\phi$.
With the distribution, we then calculate its Kullback-Leibler divergence under
perturbation, which is straightforward to do by altering the statistical
correlation functions calculated on the MCMC sample according to Eq 5.
Relating this change to a change in the fields $\\{h_{i}\\}$ and couplings
$\\{J_{ij}\\}$ is a linear matrix problem in the perturbative limit. With the
corresponding change in the parameters, we calculate the change in energy
$\Delta E(\\{s_{i}\\})$ for each configuration $s$ sampled, the coarse
graining of which maps onto the corresponding effective energy of the system.
This algorithm is specified in more detail in the supplementary information of
reference [60].

Fig B.8: Comparison of (top) eigenvalue spectrum and (bottom) eigenmatrices
calculated from Monte Carlo samples of sizes $K=10^{5}$ and $K=10^{6}$.
Eigenvalue spectrum for $K=10^{5}$ is an average calculated over 10 different
Monte Carlo samples. Top two eigenvectors are taken from a single sample as an
example, and the two columns share similar, if not identical, features.

Overall, the calculation of the FIM is an expensive computational task where
we must obtain each entry of the $(N^{2}-N)\times(N^{2}-N)$ FIM which are
averages over $K$ MCMC samples for each of $|\phi|$ coarse-grained statistics.
To parallelize such a calculation, we combine custom code with the Metropolis
algorithm implemented in the ConIII Python package [89]. We relied on a number
of computational resources including local workstations and computing clusters
at the Santa Fe Institute, University of New Mexico, and Vienna Scientific
Computing (VSC).

In the main text, we show results from MCMC samples of size $K=10^{5}$ and
compare a few cases with a larger samples of size $K=10^{6}$ to verify that we
can obtain a good approximation of FIM properties with the smaller sample
having fixed $\epsilon=10^{-4}$.111At the very least, perturbation magnitude
should be smaller than the inverse of the MCMC sample size. We verify that our
estimates of the FIM entries converge within a relative tolerance of $0.1\%$
when compared with a larger perturbation of $\epsilon=2\times 10^{-4}$. The
example of the eigenvalue spectra in Figures 4 and E.18 involve an average
over 10 samples of size $K=10^{5}$ for the pairwise maxent model the shuffled
null model. See Figure B.8 for a comparison of spectral properties of the FIM
between the two sample sizes considered.

##  Appendix C Analyzing eigenmatrices

We compare vertical striations in comparison withrows by measuring comparing
the sum of row and column uniformities as defined in Eq 8. In Figures C.11,
C.12,C.13, and C.14, we show uniformity for 3 additional neuron subsamples for
the main approach discussed in the text, the independent model, the shuffled
null model, and the MCH solutions, respectively. For pairwise perturbations,
we find matrices strongly biased to large column sum norms compared to rows.
Given the nature of the pairwise perturbations that we consider, that means
that perturbations localized to a particular neuron (in contrast with
perturbations that would impact each of the neighbors in turn) would have a
dominant impact on the collective outcomes. This is a feature of localized
control because it means that turning off all local synaptic connections, or
turning them up, for a particular neuron is effective.

As a more direct measure of this, we can analyze the subspace of the FIM
corresponding to perturbation of each neuron at a time, fixing $m$ and
iterating over all $t$. As we show in Figure C.9, the principal eigenvalues
extracted from single neuron perturbations are correlated with the fraction of
MC samples in which the same neuron has unusually strong column uniformity.
Thus, two different measures of pivotal neurons, one based on single neuron
subspaces and another based on the structure of the entire system, align and
show that the strong collective tendencies that we find in the data do not
prevent individual neurons from playing important roles.

Fig C.9: Principal subspace eigenvalues vs. fraction above column uniformity
threshold of 99-percentile (compare with Figure 5). Subspace corresponds to
diagonal blocks of FIM for fixed matcher index $m$, which is used in Figure 4C
to order neurons. Average taken over ten MC samples with standard error of the
mean shown. This means that neurons often above threshold tend to have large
collective sensitivity, linking these two different measure of pivotal
strength. Points are randomly offset on x-axis for visibility.

To determine the scaling of the rank-ordered eigenvalue spectrum, we fit the
function detailed in the main text and reproduced here

$\displaystyle\lambda(z)$ $\displaystyle=Cz^{-\alpha}e^{-z/\bar{z}}$ (14)

for eigenvalue $\lambda(z)$ with rank $z$ by performing least-squares
minimization on the logarithmic differences. However, there no a priori
guarantee that the spectrum of eigenvalues is full rank and the tail of the
spectrum may be rife with numerical precision errors, and so we allow for the
possibility of a hard cutoff that should be applied. As a heuristic, we apply
a hard cutoff when the logarithmic slope falls below $-3$ and do not fit any
points for rank above the cutoff. While we anticipate that numerical precision
makes it difficult to estimate the smallest eigenvalues and thus the exact
value of the cutoff, we take the sharp tail to be an effective cutoff for when
the eigenvalues fall below a threshold of $\lambda<10^{-7}$, minuscule in
comparison with the typical principle eigenvalue in the range of $10^{2}$ to
$10^{4}$. As is visible in our fit in Figure 4, Eq 14 with a hard cutoff hews
very well to the averaged FIM spectrum. In contrast, a simple exponential
decay, having set $\alpha=0$, cannot capture the scaling form that we observe.
We show an overview of the truncated power law fit exponents and exponential
tails in Figure C.10.

We compare our findings for the pairwise maxent model with null models
including an independent model for neural activity and one where the couplings
are randomly assigned to a pair called “coupling shuffled.” While the
independent model gives qualitatively different results, we find that
shuffling the couplings amongst pairs preserves the qualitative features of
the FIM.

Fig C.10: (A, B) Power law exponent from fitting FIM eigenvalue spectra
comparing observable and canonical perturbations. (C, D) Scale of exponential
tail $\bar{z}$. Fig C.11: Uniformity of eigenmatrices for different neural
populations of $N=50$ neurons. We show averages over Monte Carlo samples.
Error bars show one standard deviation. (top row) Subsets from experiment
154139. (bottom row) Subsets from experiment 170419. Fig C.12: Uniformity of
the eigenmatrices for independent model. There is no dramatic separation
between column and row uniformity. See Figure C.11 for details. Fig C.13:
Uniformity of eigenmatrices after shuffling couplings. Permutation mostly
preserves separated column and row uniformities. See Figure C.11 for details.
Fig C.14: Uniformity of eigenmatrices for MCH solutions. See Figure C.11 for
details.

##  Appendix D Classes of perturbations

We distinguish between two principal classes of perturbations in the main text
denoted “observable” and “natural” or “canonical,” taking the nomenclature
used in Amari’s well-known textbook on information geometry [66]. Essentially,
these distinguish between perturbations defined in the space of correlations
compared with perturbations defined in the space of parameters. Since
perturbations in either representation can be transformed to one other by a
linear operation, there is no a priori reason to prefer one basis over
another. Theoretical treatments tend to consider the canonical picture because
of the underlying intuition that they correspond to physical forces.

In physics experiments, it is possible to access directly quantities such as
fields, couplings, and temperature to directly modify the Hamiltonian. With a
phenomenological model of the neural network, however, we do not expect that
perturbations protocols simply map to canonical perturbations. This means that
for the natural perturbations we define in Eq 5, there is a corresponding
representation in terms of the fields $\\{h_{i}\\}$ and couplings
$\\{J_{ij}\\}$, but it will generally be a complicated combination of many
changes across the system (see Figures 3 and D.15). In this sense, the
pairwise maxent model serves as a representation of our inference process
rather than a literal model of the neural network.

Fig D.15: Example of simple coupling perturbation mapped to space of
replacement rule. Compare with Figure 3C.

Besides the distinction between types of perturbations, we have a choice of
which protocols to consider within each space. In the main text, we discuss
individual neuron perturbations relative to neighbors as a way of approaching
closed-loop control. However, other experimental protocols might be of
interest such as probabilistic clamping to a fixed external reference frame
always increasing, decreasing, or staying flat. An infinite variety of such
variations are possible, and it may be the case that other representations
turn out to be more appropriate for capturing simple experimental
interventions. Though a complete basis can always, in principle, be
transformed from one set of perturbations to another, such transformations may
be difficult to perform accurately from limited and potentially noisy data
extracted from experiments.

For the sake of completeness, we additionally compute the FIM for the example
of clamping each neuron to an external reference frame—akin to clamping to an
imaginary neighbor that is always in the up, down, or flat state. For the case
of $K>2$ possible states of the neuron, clamping a neuron to a single state
involves specifying the relative chance in the probabilities that the neuron
is in the remaining $K-1$ states. We consider the case where clamping the
neuron to a particular state $y$ changes the probabilities of the remaining
two configurations in a way that traces out the geodesic towards $p_{y}\equiv
p(s_{m}=y)=1$ in the two-dimensional simplex (Figure E.17). Experimental
results may suggest more empirically grounded way of accounting for the
relative change in probabilities. Taking this formulation, the local
perturbation for the matcher neuron changes its bias as

$\displaystyle\begin{aligned}
\tilde{p}_{x}&=p_{x}-\epsilon(1-p_{y})\frac{\cos(\pi/2-\theta_{y})}{\cos(\theta_{y}-\pi/6)},\\\
\tilde{p}_{y}&=(1-\epsilon)p_{y}+\epsilon,\\\
\tilde{p}_{z}&=p_{z}-\epsilon(1-p_{y})\frac{\cos(\pi/6+\theta_{y})}{\cos(\theta_{y}-\pi/6)}.\end{aligned}$
(15)

In Figure D.16, we show a summary of results from such a protocol on the
pairwise maxent model. We find that eigenvalues decay faster with rank when
considering observable perturbations compared with natural perturbations for
both fine and coarse collective synchrony. Across the maxent models we
consider and the fine and coarse measure of synchrony, we consistently find
that the decay exponent is faster for natural perturbations.

Fig D.16: FIM eigenvalue spectra from field and bias perturbations for
experiment 170419. We compare spectra for $\phi_{\rm fine}$ with $\phi_{\rm
coarse}$. Panel letters indicate different neuron subsets of $N=50$ and
averages over Monte Carlo samples with error bars representing standard
deviations.

##  Appendix E Relating observable and canonical perturbations

In the main text, we primarily focus on perturbations to the observed
correlations, but we also find consistent differences between observable and
natural perturbations. As we show in Figure 4, the spectrum tends to decay
faster for canonical perturbations and they are smaller in general. To gain
some insight into these differences, we consider the two types of
perturbations for a highly simplified case of a binary neuron characterized by
a mean activity level.

Consider states $s$ with probabilities $p(s)=\exp[-E(s)]/Z$ coarse-grained
into a distribution $\phi(k)=\sum_{|s|=k}p(s)$, where the notation $|s|=k$
means that there are $k$ spins in the majority. Under the replacement rule,
the probability of neuron $i$ being in one configuration $p_{i}$ is modified
to become

$\displaystyle\begin{aligned} \tilde{p}_{i}&=(1-\epsilon)p_{i}+\epsilon\\\
&=p_{i}+(1-p)\epsilon.\end{aligned}$ (16)

Note that the derivative that returns the appropriate is of the form
$-d/d[\log(1-p_{i})]=(1-p_{i})d/dp_{i}$. This would correspond to our
observable perturbation, whereas the canonical perturbation would be with
respect to the Langrangian multipliers, the fields and couplings for the
pairwise maxent model.

The quantity of interest is the “score function,” or the expectation value of
the second derivative of the probability distribution that gives us the
entries of the FIM,

$\displaystyle-\left<\frac{\partial^{2}\log\phi}{\partial[\log(1-p_{i})]^{2}}\right>.$
(17)

The terms involving $\log\phi$ are of course the susceptibilities of the $k$th
correlation function with the spin $i$. Since $\phi$ is a coarse-grained
version of $p(s)$, we expect that the derivatives will reduces to linear
combinations of correlation functions.

To simplify notation, we take the term inside the brackets having defined
$y_{i}\equiv\log(1-p_{i})$ to obtain

$\displaystyle\begin{aligned} \frac{\partial^{2}\log\phi}{\partial
y_{i}^{2}}&=\frac{\partial}{\partial y_{i}}\left[\frac{\partial
h_{i}}{\partial y_{i}}\frac{\partial\log\phi}{\partial h_{i}}\right]\\\
&=\frac{\partial^{2}h_{i}}{\partial y_{i}^{2}}\frac{\partial\log\phi}{\partial
h_{i}}+\left[\frac{\partial h_{i}}{\partial
y_{i}}\right]^{2}\frac{\partial^{2}\log\phi}{\partial
h_{i}^{2}},\end{aligned}$ (18)

where we have introduced the field $h_{i}$ for neuron $i$. By taking the
derivative with respect to the field, we have incurred a term proportional to
the curvature in the change of variables as well as the square of the jacobian
bringing us from $y_{i}$ to $h_{i}$.

Now, we calculate the jacobian terms. Note that we are only considering the
domain where the relationship between $\log(1-p_{i})$ and $h_{i}$ is analytic,
otherwise we would have to deal with branch cuts. Given this caveat, we
differentiate

$\displaystyle\frac{\partial y_{i}}{\partial h_{i}}$
$\displaystyle=\frac{1}{p_{i}-1}\frac{\partial p_{i}}{\partial h_{i}}.$ (19)
Using the fact that $p_{i}=(\left<s_{i}\right>+1)/2$ and that the derivative
is the susceptibility, or the variance of spin $i$, we obtain
$\displaystyle=\frac{1}{\left<s_{i}\right>-1}\left(1-\left<s_{i}\right>^{2}\right)$
(20) $\displaystyle=-(1+\left<s_{i}\right>).$ (21)

Then, the second derivative is

$\displaystyle\frac{\partial^{2}y_{i}}{\partial h_{i}^{2}}$
$\displaystyle=-\left(1-\left<s_{i}\right>^{2}\right).$ (22)

Thus, the maxent model allows us to explicitly calculate the jacobian in terms
of linear response quantities that tell us how the observables change under a
small perturbation of the fields.

Now, let us go back to the score function in Eq 17. Note that because all the
above quantities are already ensemble averages, we can factor them out of the
brackets in Eq 17. This means that we have

$\displaystyle\frac{1}{\left<s_{i}\right>+1}\left[\frac{1}{\left<s_{i}\right>-1}\left<\frac{\partial\log\phi}{\partial
h_{i}}\right>-\frac{1}{\left<s_{i}\right>+1}\left<\frac{\partial^{2}\log\phi}{\partial
h_{i}^{2}}\right>\right].$ (23)

In contrast, if we were to do the same thing but take the derivative with
respect to the fields in Eq 18 instead of the natural parameter, we would
simply swap the $y_{i}$’s with the $h_{i}$’s. This means that we have the
reciprocal of the jacobians,

$\displaystyle(\left<s_{i}\right>+1)\left[(\left<s_{i}\right>-1)\left<\frac{\partial\log\phi}{\partial
y_{i}}\right>-(\left<s_{i}\right>+1)\left<\frac{\partial^{2}\log\phi}{\partial
y_{i}^{2}}\right>\right].$ (24)

In other words, changing variables gives us a different factor that depends on
the mean magnetization of the entries of the FIM. This means that we may
expect the overall eigenvalue spectrum to be, at the least, scaled
differently. From numerical calculations, we find that the spectrum for
observable perturbations to almost always be strictly greater than that for
natural perturbations, which may reflect a change in variables.

Fig E.17: Probability simplex for three-state neuron with probabilities of
being in each of three states denoted by $p_{x}$, $p_{y}$, and $p_{z}$ such
that $p_{x}+p_{y}+p_{z}=1$. Fig E.18: FIM overview for worm experiment 154139.
In main text, we show corresponding plots for experiment 170419 in Figure 4.
Fig E.19: FIM overview for worm experiment 170419 with $\phi_{\rm coarse}$
instead of $\phi_{\rm fine}$ as shown in main text in Figure 4.

Fig E.20: FIM overview for worm experiment 154139 with $\phi_{\rm coarse}$
instead of $\phi_{\rm fine}$ as shown in main text in Figure 4.

##  Appendix F Data sets and samples

We use the C. elegans worm data from reference [4] and the same inverse maxent
procedures detailed in reference [50]. We analyze the two immobilized worms
labeled with experiment numbers 170419 and 154139. For each worm, we consider
a limited number of different subsets of $N=50$ neurons (about the maximum
number of neurons that can be reliable measured from the data), where the
neurons have been randomly selected amongst those that have visited all three
states at least once.

##  Appendix G Experimental implementation

A realization of our thought experiment depends on developments in
simultaneous use of recording, perturbation, and computational analysis. The
experiment would require tracking time-averaged neural statistics before and
during perturbation with the ability to extract simultaneously the discrete
state of recorded neurons, a procedure that currently relies on post-
experimental analysis to handle noise and changing fluorescence [50].
Furthermore, single-neuron tracking is difficult and may be feasible with
immobilized worms, but presents a challenge to perform accurately with freely
moving ones. Furthermore, it is essential that the nature of the perturbation
on membrane voltage be precisely calibrated in order to replicate as closely
as possible theoretical clamping, which may require some characterization of
the particular experimental techniques used. This presents an abbreviated list
of some of the experimental difficulties that we foresee for implementing such
an experiment.

Though we assume that the perturbation randomly “flips” the neural
configuration of the focus neuron with some small probability at any given
moment in time, the perturbation specified in Eq 5 is compatible with
variations of the experiment that may be more natural. For example, it may be
more natural to force neurons up after a flat state but not directly from a
down state. More specifically, it may be important to keep the slope of
underlying neural activity from diverging. While such history-dependent
clamping is compatible with our formulation that relies only on time averages,
a rigorous comparison of techniques might reveal some protocols preferable
over others in terms of feasibility or correspondence to theoretical
predictions.

Finally, our proposed thought experiment generalizes beyond neural activity or
the C. elegans model system, where it may be possible to simultaneously
observe system configuration while imposing a perturbation. With optogenetic
tools, it may be interesting to consider other neural systems such as in vitro
cortical cultures or the neurons involved in Drosophila flight.

Fig G.21: Frequency of repeated pivotal neurons from experiment 154139. See
Figure 5 for more details. Fig G.22: Power law with truncated exponential fit
as dashed black lines to mean eigenvalue spectra for experiment 170419.
Sensitivity for $\phi_{\rm fine}$ in blue and $\phi_{\rm coarse}$ in orange.
Error bars show standard deviation over 10 Monte Carlo samples used to
calculate FIM over 4 different subsets of neurons shown in each panel. Fig
G.23: Power law with truncated exponential fit to mean eigenvalue spectra for
experiment 154139. See Figure G.22 for more details.

## References

  * 1. Goh KI, Cusick ME, Valle D, Childs B, Vidal M, Barabasi AL.  The Human Disease Network.  Proc Natl Acad Sci USA. 2007;104(21):8685–8690.  doi:10.1073/pnas.0701361104.
  * 2. Vidal M, Cusick ME, Barabási AL.  Interactome Networks and Human Disease.  Cell. 2011;144:986–998.
  * 3. Zhang B, Gaiteri C, Bodea LG, Wang Z, McElwee J, Podtelezhnikov AA, et al.  Integrated Systems Approach Identifies Genetic Nodes and Networks in Late-Onset Alzheimer’s Disease.  Cell. 2013;153(3):707–720.  doi:10.1016/j.cell.2013.03.030.
  * 4. Scholz M, Linder AN, Randi F, Sharma AK, Yu X, Shaevitz JW, et al.  Predicting Natural Behavior from Whole-Brain Neural Dynamics.  Neuroscience; 2018.
  * 5. Morone F, Makse HA.  Symmetry Group Factorization Reveals the Structure-Function Relation in the Neural Connectome of Caenorhabditis Elegans.  Nat Commun. 2019;10(1):4961.  doi:10.1038/s41467-019-12675-8.
  * 6. Liu YY, Slotine JJ, Barabási AL.  Controllability of Complex Networks.  Nature. 2011;473(7346):167–173.  doi:10.1038/nature10011.
  * 7. Tang Y, Gao H, Zou W, Kurths J.  Identifying Controlling Nodes in Neuronal Networks in Different Scales.  PLoS ONE. 2012;7(7):e41375.  doi:10.1371/journal.pone.0041375.
  * 8. Houweling AR, Brecht M.  Behavioural report of single neuron stimulation in somatosensory cortex.  Nature. 2008;451(7174):65–68.
  * 9. Huber D, Petreanu L, Ghitani N, Ranade S, Hromádka T, Mainen Z, et al.  Sparse optical microstimulation in barrel cortex drives learned behaviour in freely moving mice.  Nature. 2008;451(7174):61–64.
  * 10. Kato S, Kaplan HS, Schrödel T, Skora S, Lindsay TH, Yemini E, et al.  Global brain dynamics embed the motor command sequence of Caenorhabditis elegans.  Cell. 2015;163(3):656–669.
  * 11. Susoy V, Hung W, Witvliet D, Whitener JE, Wu M, Graham BJ, et al.  Natural sensory context drives diverse brain-wide activity during C. elegans mating.  bioRxiv. 2020;doi:10.1101/2020.09.09.289454.
  * 12. Yan G, Vértes PE, Towlson EK, Chew YL, Walker DS, Schafer WR, et al.  Network Control Principles Predict Neuron Function in the Caenorhabditis Elegans Connectome.  Nature. 2017;550(7677):519–523.  doi:10.1038/nature24056.
  * 13. Borriello E, Daniels BC.  The basis of easy controllability in Boolean networks.  in review. 2021; p. arXiv:2010.12075.
  * 14. Morone F, Leifer I, Makse HA.  Fibration Symmetries Uncover the Building Blocks of Biological Networks.  Proc Natl Acad Sci USA. 2020;117(15):8306–8314.  doi:10.1073/pnas.1914628117.
  * 15. Zañudo JGT, Yang G, Albert R, Levine H.  Structure-based control of complex networks with nonlinear dynamics.  Proceedings of the National Academy of Sciences of the United States of America. 2017;114(28):7234–7239.  doi:10.1073/pnas.1617387114.
  * 16. Del Ferraro G, Moreno A, Min B, Morone F, Pérez-Ramírez Ú, Pérez-Cervera L, et al.  Finding Influential Nodes for Integration in Brain Networks Using Optimal Percolation Theory.  Nat Commun. 2018;9(1):2274.  doi:10.1038/s41467-018-04718-3.
  * 17. Lynn CW, Bassett DS.  The Physics of Brain Network Structure, Function, and Control.  Nat Rev Phys. 2019;1(5):318–332.  doi:10.1038/s42254-019-0040-8.
  * 18. Boyden ES, Zhang F, Bamberg E, Nagel G, Deisseroth K.  Millisecond-timescale, genetically targeted optical control of neural activity.  Nature neuroscience. 2005;8(9):1263–1268.
  * 19. Mardinly AR, Oldenburg IA, Pégard NC, Sridharan S, Lyall EH, Chesnov K, et al.  Precise multimodal optical control of neural ensemble activity.  Nature neuroscience. 2018;21(6):881–893.
  * 20. Grosenick L, Marshel JH, Deisseroth K.  Closed-loop and activity-guided optogenetic control.  Neuron. 2015;86(1):106–139.
  * 21. Zrenner C, Belardinelli P, Müller-Dahlhaus F, Ziemann U.  Closed-loop neuroscience and non-invasive brain stimulation: a tale of two loops.  Frontiers in cellular neuroscience. 2016;10:92.
  * 22. Kim CK, Adhikari A, Deisseroth K.  Integration of Optogenetics with Complementary Methodologies in Systems Neuroscience.  Nat Rev Neurosci. 2017;18(4):222–235.  doi:10.1038/nrn.2017.15.
  * 23. Pokala N, Liu Q, Gordus A, Bargmann CI.  Inducible and Titratable Silencing of Caenorhabditis Elegans Neurons in Vivo with Histamine-Gated Chloride Channels.  Proc Natl Acad Sci USA. 2014;111(7):2770–2775.  doi:10.1073/pnas.1400615111.
  * 24. Gutruf P, Krishnamurthi V, Vázquez-Guardado A, Xie Z, Banks A, Su CJ, et al.  Fully implantable optoelectronic systems for battery-free, multimodal operation in neuroscience research.  Nature Electronics. 2018;1(12):652–660.
  * 25. Barlow HB.  Single Units and Sensation: A Neuron Doctrine for Perceptual Psychology?  Perception. 1972;1:371–394.
  * 26. Quiroga RQ, Kreiman G, Koch C, Fried I.  Sparse but Not ‘Grandmother-Cell’ Coding in the Medial Temporal Lobe.  Cell. 2007;12(3):87–91.
  * 27. Daniels BC, Krakauer DC, Flack JC.  Sparse Code of Conflict in a Primate Society.  Proc Natl Acad Sci USA. 2012;109(35):14259–14264.  doi:10.1073/pnas.1203021109.
  * 28. Olshausen BA, Field DJ.  Sparse Coding with an Overcomplete Basis Set: A Strategy Employed by V1?  Vision Research. 1997;37(23):3311–3325.  doi:10.1016/S0042-6989(97)00169-7.
  * 29. Spanne A, Jörntell H.  Questioning the Role of Sparse Coding in the Brain.  Trends in Neurosciences. 2015;38(7):417–427.  doi:10.1016/j.tins.2015.05.005.
  * 30. Beyeler M, Rounds EL, Carlson KD, Dutt N, Krichmar JL.  Neural Correlates of Sparse Coding and Dimensionality Reduction.  PLoS Comput Biol. 2019;15(6):e1006908.  doi:10.1371/journal.pcbi.1006908.
  * 31. Oizumi M, Okada M, Amari SI.  Information Loss Associated with Imperfect Observation and Mismatched Decoding.  Front Comput Neurosci. 2011;5.  doi:10.3389/fncom.2011.00009.
  * 32. Gray JM, Hill JJ, Bargmann CI.  A circuit for navigation in Caenorhabditis elegans.  Proceedings of the National Academy of Sciences. 2005;102(9):3184–3191.
  * 33. Ikeda M, Nakano S, Giles AC, Xu L, Costa WS, Gottschalk A, et al.  Context-Dependent Operation of Neural Circuits Underlies a Navigation Behavior in _Caenorhabditis_ _Elegans_.  Proc Natl Acad Sci USA. 2020;117(11):6178–6188.  doi:10.1073/pnas.1918528117.
  * 34. Ouellette MH, Desrochers MJ, Gheta I, Ramos R, Hendricks M.  A Gate-and-Switch Model for Head Orientation Behaviors in _Caenorhabditis_ _Elegans_.  eNeuro. 2018;5(6):ENEURO.0121–18.2018.  doi:10.1523/ENEURO.0121-18.2018.
  * 35. Chang AJ, Bargmann CI.  Hypoxia and the HIF-1 Transcriptional Pathway Reorganize a Neuronal Circuit for Oxygen-Dependent Behavior in Caenorhabditis Elegans.  Proc Natl Acad Sci USA. 2008;105(20):7321–7326.  doi:10.1073/pnas.0802164105.
  * 36. Lashley KS.  Brain Mechanisms and Intelligence: A Quantitative Study of Injuries to the Brain.  Behavior Research Fund Monographs. University of Chicago Press; 1929\.
  * 37. Damaslo H, Grabowski T, Frank R, Galaburda AM, Damasio AR.  The Return of Phineas Gage: Clues About the Brain from the Skull of a Famous Patient.  Science. 1994;264:5.
  * 38. Stephens GJ, Johnson-Kerner B, Bialek W, Ryu WS.  Dimensionality and Dynamics in the Behavior of C. Elegans.  PLoS Comput Biol. 2008;4(4):e1000028.  doi:10.1371/journal.pcbi.1000028.
  * 39. Schneidman E, Berry MJ, Segev R, Bialek W.  Weak Pairwise Correlations Imply Strongly Correlated Network States in a Neural Population.  Nature. 2006;440(7087):1007–1012.  doi:10.1038/nature04701.
  * 40. Hopfield JJ.  Neural Networks and Physical Systems with Emergent Collective Computational Abilities.  Proc Natl Acad Sci USA. 1982;79:2554–2558.
  * 41. Transtrum MK, Qiu P.  Model Reduction by Manifold Boundaries.  Phys Rev Lett. 2014;113(9):098701.  doi:10.1103/PhysRevLett.113.098701.
  * 42. Transtrum MK, Machta BB, Brown KS, Daniels BC, Myers CR, Sethna JP.  Perspective: Sloppiness and Emergent Theories in Physics, Biology, and Beyond.  J Chem Phys. 2015;143(1):010901.  doi:10.1063/1.4923066.
  * 43. Maunsell JH, Van Essen DC.  Functional Properties of Neurons in Middle Temporal Visual Area of the Macaque Monkey. I. Selectivity for Stimulus Direction, Speed, and Orientation.  Journal of Neurophysiology. 1983;49(5):1127–1147.  doi:10.1152/jn.1983.49.5.1127.
  * 44. Wu S, Amari Si, Nakahara H.  Population Coding and Decoding in a Neural Field: A Computational Study.  Neural Computation. 2002;14(5):999–1026.  doi:10.1162/089976602753633367.
  * 45. Nguyen JP, Shipley FB, Linder AN, Plummer GS, Liu M, Setru SU, et al.  Whole-brain calcium imaging with cellular resolution in freely behaving Caenorhabditis elegans.  Proceedings of the National Academy of Sciences. 2016;113(8):E1074–E1081.
  * 46. Venkatachalam V, Ji N, Wang X, Clark C, Mitchell JK, Klein M, et al.  Pan-Neuronal Imaging in Roaming _Caenorhabditis_ _Elegans_.  Proc Natl Acad Sci USA. 2016;113(8):E1082–E1088.  doi:10.1073/pnas.1507109113.
  * 47. Stirman JN, Crane MM, Husson SJ, Wabnig S, Schultheis C, Gottschalk A, et al.  Real-time multimodal optical control of neurons and muscles in freely behaving Caenorhabditis elegans.  Nature methods. 2011;8(2):153–158.
  * 48. Berman GJ, Bialek W, Shaevitz JW.  Predictability and Hierarchy in _Drosophila_ Behavior.  Proc Natl Acad Sci USA. 2016;113(42):11943–11948.  doi:10.1073/pnas.1607601113.
  * 49. Merchan L, Nemenman I.  On the Sufficiency of Pairwise Interactions in Maximum Entropy Models of Networks.  J Stat Phys. 2016;162(5):1294–1308.  doi:10.1007/s10955-016-1456-5.
  * 50. Chen X, Randi F, Leifer AM, Bialek W.  Searching for Collective Behavior in a Small Brain.  Phys Rev E. 2019;99(5):052418.  doi:10.1103/PhysRevE.99.052418.
  * 51. Xu S, Chisholm AD.  Highly Efficient Optogenetic Cell Ablation in C. Elegans Using Membrane-Targeted miniSOG.  Sci Rep. 2016;6(1):21271.  doi:10.1038/srep21271.
  * 52. Suk HJ, van Welie I, Kodandaramaiah SB, Allen B, Forest CR, Boyden ES.  Closed-loop real-time imaging enables fully automated cell-targeted patch-clamp neural recording in vivo.  Neuron. 2017;95(5):1037–1047.
  * 53. Newman JP, Fong Mf, Millard DC, Whitmire CJ, Stanley GB, Potter SM.  Optogenetic feedback control of neural activity.  Elife. 2015;4:e07192.
  * 54. Prinz AA, Abbott L, Marder E.  The dynamic clamp comes of age.  Trends in neurosciences. 2004;27(4):218–224.
  * 55. Morcos F, Pagnani A, Lunt B, Bertolino A, Marks DS, Sander C, et al.  Direct-Coupling Analysis of Residue Coevolution Captures Native Contacts across Many Protein Families.  Proc Natl Acad Sci USA. 2011;108(49):E1293–E1301.
  * 56. Volkov I, Banavar JR, Hubbell SP, Maritan A.  Inferring Species Interactions in Tropical Forests.  Proc Natl Acad Sci USA. 2009;106(33):13854–13859.  doi:10.1073/pnas.0903244106.
  * 57. Cover TM, Thomas JA.  Elements of Information Theory.  2nd ed. Hoboken: John Wiley & Sons; 2006.
  * 58. Jaynes ET.  Information Theory and Statistical Mechanics.  Phys Rev. 1957;106(4):620–630.  doi:10.1103/PhysRev.106.620.
  * 59. Lee ED, Broedersz CP, Bialek W.  Statistical Mechanics of the US Supreme Court.  J Stat Phys. 2015;160(2):275–301.  doi:10.1007/s10955-015-1253-6.
  * 60. Lee ED, Katz DM, Bommarito II MJ, Ginsparg PH.  Sensitivity of Collective Outcomes Identifies Pivotal Components.  J R Soc Interface. 2020;17(20190873).
  * 61. Tkacik G, Schneidman E, Berry II MJ, Bialek W.  Spin Glass Models for a Network of Real Neurons.  arXiv:09125409 [q-bio]. 2009;.
  * 62. Tkačik G, Marre O, Amodei D, Schneidman E, Bialek W, Berry MJ.  Searching for Collective Behavior in a Large Network of Sensory Neurons.  PLoS Comput Biol. 2014;10(1):e1003408.  doi:10.1371/journal.pcbi.1003408.
  * 63. Meshulam L, Gauthier JL, Brody CD, Tank DW, Bialek W.  Coarse Graining, Fixed Points, and Scaling in a Large Population of Neurons.  Phys Rev Lett. 2019;123(17):178103.  doi:10.1103/PhysRevLett.123.178103.
  * 64. Barton J, Cocco S.  Ising Models for Neural Activity Inferred via Selective Cluster Expansion: Structural and Coding Properties.  J Stat Mech. 2013;2013(03):P03002.  doi:10.1088/1742-5468/2013/03/P03002.
  * 65. Roudi Y, Nirenberg S, Latham PE.  Pairwise Maximum Entropy Models for Studying Large Biological Systems: When They Can Work and When They Can’t.  PLoS Comput Biol. 2009;5(5):e1000380.  doi:10.1371/journal.pcbi.1000380.
  * 66. Amari Si.  Information Geometry and Its Applications. vol. 194 of Applied Mathematical Sciences.  Springer Japan; 2016.
  * 67. Bialek W, Ranganathan R.  Rediscovering the Power of Pairwise Interactions.  arXiv:07124397 [q-bio]. 2007;.
  * 68. Tkačik G, Marre O, Mora T, Amodei D, Berry II MJ, Bialek W.  The Simplest Maximum Entropy Model for Collective Behavior in a Neural Network.  J Stat Mech. 2013;2013(03):P03011.  doi:10.1088/1742-5468/2013/03/P03011.
  * 69. Bialek WS.  Biophysics: Searching for Principles.  Princeton, NJ: Princeton University Press; 2012.
  * 70. As a more intuitive formulation of Eq 6, we could assign to each possible collective configuration a “synchrony energy” $\mathcal{E}$ such that its probability can be written $\phi\propto\exp[-\mathcal{E}]$. Such an effective energy represents a coarse-graining over the microscopic states corresponding to a collective configuration. Then, Eq 6 can also be written as the limiting quantity involving the change in energies $\Delta\mathcal{E}$ under such a perturbation, the variance $\lim_{\epsilon\rightarrow 0}\left(\left<\Delta\mathcal{E}^{2}\right>-\left<\Delta\mathcal{E}\right>^{2}\right)/\epsilon^{2}$ [60]. This is a measure of how differently the log-probability of each collective configuration changes. Thus, the Fisher information is proportional to the variance of the effective energy such that more sensitive directions of change are ones that maximally scatter the collective distribution.
  * 71. Daniels BC, Flack JC, Krakauer DC.  Dual Coding Theory Explains Biphasic Collective Computation in Neural Decision-Making.  Front Neurosci. 2017;11:313.  doi:10.3389/fnins.2017.00313.
  * 72. Mitra D.  Wmatrix and the Geometry of Model Equivalence and Reduction.  Proc Inst Electr Eng UK. 1969;116(6):1101.  doi:10.1049/piee.1969.0206.
  * 73. Roy P, Cela A, Hamam Y.  On the Relation of FIM and Controllability Gramian.  In: 2009 IEEE International Symposium on Industrial Embedded Systems. Lausanne, Switzerland: IEEE; 2009. p. 37–41.
  * 74. Liu J, Elia N.  Convergence of Fundamental Limitations in Information, Estimation, and Control.  In: Proceedings of the 45th IEEE Conference on Decision and Control. San Diego, CA, USA: IEEE; 2006. p. 5609–5614.
  * 75. Rajan K, Abbott LF.  Eigenvalue Spectra of Random Matrices for Neural Networks.  Phys Rev Lett. 2006;97(18):188104.  doi:10.1103/PhysRevLett.97.188104.
  * 76. Weigt M, White RA, Szurmant H, Hoch JA, Hwa T.  Identification of Direct Residue Contacts in Protein-Protein Interaction by Message Passing.  Proc Natl Acad Sci USA. 2009;106(1):67–72.  doi:10.1073/pnas.0805923106.
  * 77. Oliveri P, Tu Q, Davidson EH.  Global Regulatory Logic for Specification of an Embryonic Cell Lineage.  Proc Natl Acad Sci USA. 2008;105(16):5955–5962.  doi:10.1073/pnas.0711220105.
  * 78. Li E, Cui M, Peter IS, Davidson EH.  Encoding Regulatory State Boundaries in the Pregastrular Oral Ectoderm of the Sea Urchin Embryo.  Proc Natl Acad Sci USA. 2014;111(10):E906–E913.  doi:10.1073/pnas.1323105111.
  * 79. In calculating the Fisher information matrix, we use $\epsilon=10^{-4}$, but the results do not depend strongly at this range of perturbation size. Furthermore, this lower bound on the strength of perturbation is complemented from above by two additional bounds. First, it is unlikely that the system can relax back from a large perturbation. Second, a sufficiently large perturbation breaks the linear assumptions we make. These upper bounds might not separate with the lower bound for measurement, which would pose a problem for measurement, but this is a question that must be answered experimentally.
  * 80. Candès EJ, Recht B.  Exact Matrix Completion via Convex Optimization.  Found Comput Math. 2009;9(6):717–772.  doi:10.1007/s10208-009-9045-5.
  * 81. Shannon CE.  A Mathematical Theory of Communication.  Bell Syst Tech J. 1948;27:379–423, 623–656.
  * 82. In the case of a Boltzmann-type model, this perturbative term can be expressed as a change in the energy of state $Y$ relative to the average change in energy. For the pairwise maxent model, the energy is a sum over the set of all couplings, not necessarily only the ones associated with the matcher neurons and its neighbors. If the structure of the coupling network were locally tree-like such that local neighbors interacted more strongly with the matcher than ones further away, then the change in energy could be approximated as a sum only over the couplings connecting the matcher with its neighbors. This is not generally the case for the networks we consider in the data.
  * 83. Brenner N, Bialek W, de Ruyter van Steveninck R.  Adaptive Rescaling Maximizes Information Transmission.  Neuron. 2000;26(3):695–702.  doi:10.1016/S0896-6273(00)81205-2.
  * 84. Transtrum MK, Machta BB, Sethna JP.  Geometry of Nonlinear Least Squares with Applications to Sloppy Models and Optimization.  Phys Rev E. 2011;83(3):036701.  doi:10.1103/PhysRevE.83.036701.
  * 85. Transtrum MK, Machta BB, Sethna JP.  Why Are Nonlinear Fits to Data so Challenging?  Phys Rev Lett. 2010;104(6):060201.  doi:10.1103/PhysRevLett.104.060201.
  * 86. Machta BB, Chachra R, Transtrum MK, Sethna JP.  Parameter Space Compression Underlies Emergent Theories and Predictive Models.  Science. 2013;342(6158):604–607.  doi:10.1126/science.1238723.
  * 87. Daniels BC, Dobrzyński M, Fey D.  Parameter Estimation, Sloppiness, and Model Identifiability.  In: Munsky B, editor. Quantitative Biology: Theory, Computational Methods, and Models. Cambridge, Massachusetts: The MIT Press; 2018. p. 271–292.
  * 88. Nguyen HC, Zecchina R, Berg J.  Inverse Statistical Problems: From the Inverse Ising Problem to Data Science.  Advances in Physics. 2017;66(3):197–261.  doi:10.1080/00018732.2017.1341604.
  * 89. Lee ED, Daniels BC.  Convenient Interface to Inverse Ising (ConIII): A Python 3 Package for Solving Ising-Type Maximum Entropy Models.  JORS. 2019;7(1):3.  doi:10.5334/jors.217.
  * 90. Dudik M, Phillips SJ, Schapire RE.  Performance guarantees for regularized maximum entropy density estimation.  In: International Conference on Computational Learning Theory. Springer; 2004. p. 472–486.
  * 91. Broderick T, Dudik M, Tkacik G, Schapire RE, Bialek W.  Faster Solutions of the Inverse Pairwise Ising Problem.  arXiv:07122437 [cond-mat, q-bio]. 2007.