# Exploring a strongly non–Markovian animal behavior

Vasyl Alba,1,2 Gordon J. Berman,1,3 William Bialek,1,4 and Joshua W. Shaevitz1
1Joseph Henry Laboratories of Physics and Lewis–Sigler Institute for
Integrative Genomics, Princeton University, Princeton NJ 08544
2Department of Engineering Science and Applied Mathematics and NSF–Simons
Center for Quantitative Biology, Northwestern University, Evanston IL 60208
3Departments of Biology and Physics, Emory University, Atlanta GA 30322
4Initiative for the Theoretical Sciences, The Graduate Center, City University
of New York, 365 Fifth Avenue, New York NY 10016

###### Abstract

A freely walking fly visits roughly 100 stereotyped states in a strongly non-
Markovian sequence. To explore these dynamics, we develop a generalization of
the information bottleneck method, compressing the large number of behavioral
states into a more compact description that maximally preserves the
correlations between successive states. Surprisingly, preserving these short
time correlations with a compression into just two states captures the long
ranged correlations seen in the raw data. Having reduced the behavior to a
binary sequence, we describe the distribution of these sequences by an Ising
model with pairwise interactions, which is the maximum entropy model that
matches the two-point correlations. Matching the correlation function at
longer and longer times drives the resulting model toward the Ising model with
inverse square interactions and near zero magnetic field. The emergence of
this statistical physics problem from the analysis real data on animal
behavior is unexpected.

In the twentieth century there were two very distinct approaches to the
characterization of animal behavior. Ethologists focused on behavior in its
natural context, often describing in qualitative terms phenomena of great
complexity frisch_74 ; gould_82 . Psychophysicists, in contrast, brought
behavior into the laboratory, constraining subjects to choosing from a small
set of alternatives, quantifying the probabilities of different choices as a
function of sensory inputs and task constraints lawson+uhlenbeck_50 ;
green+swets_66 . Recently there has been considerable interest in quantifying
unconstrained and more naturalistic behaviors. Examples include exploring the
variability of eye movement trajectories in primates osborne+al_05 , and the
postural dynamics of freely moving C elegans stephens+al_08 ; ahamed+al_20 ,
walking flies branson+al_09 ; berman+al_14 ; berman+al_16 , and mice
wiltschko+al_15 ; marshall+al_20 . The combination of high resolution video
imaging and efficient AI tools mathis+al_18 ; pereira+al_19  is making these
approaches more generally applicable, and these data have focused attention on
the wide range of time scales in the behavior of single organisms, from
milliseconds to a lifetime brown+bivort_17 ; datta+al_19 ; pereira+al_20 ;
mathis+mathis_20 ; bialek_20 .

The behavior of a fly walking freely in a featureless arena can be described
as a sequence of transitions among discrete, stereotyped states berman+al_14 .
Some of these states correspond to actions with a simple verbal description,
such as grooming particular body parts, while other states are not so simple;
nonetheless all flies of a single species revisit the same $\sim$100 states,
and one can recognize the same states in closely related species. Although the
description of behavior as a sequence of states often is accompanied an
analysis of transitions from one state to the next, and the (possibly
implicit) hypothesis that these transitions are independent of one another,
fly behavior dramatically violates this Markovian assumption berman+al_16 . We
would like to find a simple phenomenological description of behavioral
sequences that captures this non–Markovian structure, ideally leading to some
insight into the internal states that make it possible.

The non–Markovian character of the fly’s trajectory through state space can be
seen in several ways. In particular, if we look at the probability that the
fly moves from one state to another after a time $\tau$, then in a Markov
model all of these transition matrices are powers of the fundamental
transition matrix describing the model, and the eigenvalues of this matrix
will decay at a constant rate. This predicts that essentially all memory of
the initial state will be lost after $\sim 30$ transitions. Instead, memory
persists out to $\sim 1000$ transitions, and the eigenvalues of the transition
matrix decay more and more slowly as we look across longer time scales
berman+al_16 .

Figure 1: Probability $P_{c}(\tau)$ that a behavioral state is revisited after
$\tau$ transitions. We normalize so that
${\tilde{P}}_{c}(\tau)=P_{c}(\tau)/P_{c}(\tau\rightarrow\infty)-1$ decays to
zero. Means and standard errors across 59 flies (red) compared with
predictions from a Markov model that matches observed transition probabilities
(green). In both cases we analyze individual flies and average
${\tilde{P}}_{c}$ at the end.

We can see the same effect in the mutual information between states separated
by a time $\tau$, or in the probability that the fly returns to the same state
after $\tau$ transitions, $P_{c}(\tau)$, shown in Fig 1. In both of these
representations we detect correlations in behavioral state extending over
thousands of transitions, far beyond what is predicted by a Markov model
matched to the measured transition probabilities at $\tau=1$. More subtly, the
decay of mutual information or return probability is gradual, and does not
seem to be characterized by discrete time scales, echoing the slowing of
eigenvalue decays.

The problem that we face in characterizing the non–Markovian structure of
behavior is that with more than one hundred states, there are $\sim 10^{4}$
possible transitions between pairs of states. If we observe a single fly for
one hour, we see roughly this number of transitions berman+al_14 ;
berman+al_16 . If we make much longer observations on single flies we
encounter non–stationarity, while if we merge data from many flies we may
obscure individual differences. To make progress we need to simplify our
description of the behavioral states, but we need to do this in a way that
preserves the long memory seen in the full description.

Concretely, we want to map the behavioral state at each moment in time into
some reduced or clustered description, $x_{t}\rightarrow z_{t}$, with the hope
that we can preserve the long memories that we see in the state sequence; in
general this mapping can be probabilistic, described by $P(z_{t}|x_{t})$. The
‘‘information bottleneck’’ problem tishby+al_99  is the search for a
compression that maintains information about some other variable, which we
could take to be a behavioral state in the future berman+al_16 , but here we
want a more symmetric formulation (Fig 2).

Figure 2: Mapping of behavioral states $x$ into a compressed representation
$z$. We compress the state at each moment in time independently, with a
mapping $x\rightarrow z$ that is in general probabilistic, $P(z|x)$. The
compressed representation captures an information $I(z_{t};x_{t})$ about the
original state, and the compressed variables share information
$I(z_{t};z_{t+\tau})$ across time $\tau$. We choose the compression to
maximize $I(z_{t};z_{t+\tau})$ at fixed $I(z_{t};x_{t})$, as in Eq (1).

We choose compressions that preserve temporal correlations by maximizing the
mutual information between compressed variables at different times,
$I(z_{t};z_{t+\tau})$. To control the complexity of our description, we limit
the information that we capture about the original variables,
$I(z_{t};x_{t})$. In principle the mapping $x_{t}\rightarrow z_{t}$ could be
probabilistic, so we solve the variational problem

$\max_{P(z|x)}\left[I(z_{t};z_{t+\tau})-TI(z_{t};x_{t})\right],$ (1)

where the Lagrange multiplier $T$ imposes the constraint on $I(z_{t};x_{t})$.
As in the original bottleneck problem, $T$ plays the role of a temperature,
such that as $T\rightarrow 0$ the mapping $x\rightarrow z$ becomes
deterministic.

Following the same steps as in Ref tishby+al_99 , we find that the solution to
the maximization problem in Eq (1) obeys the self–consistent equation

$\displaystyle P(z|x)$ $\displaystyle=$ $\displaystyle\frac{P(z)}{{\cal
Z}(x;T)}\exp\left[-\frac{1}{T}F(x,z)\right],$ (2) $\displaystyle
F(x_{t},z_{t})$ $\displaystyle=$ $\displaystyle
D_{KL}\big{[}P(z_{t+\tau}|x_{t})||P(z_{t+\tau}|z_{t})\big{]}$ (3)
$\displaystyle\,\,\,\,\,\,\,\,\,\,+D_{KL}\big{[}P(z_{t+\tau}|x_{t})||P(z_{t+\tau}|z_{t})\big{]},$

where $D_{KL}[P||Q]$ is the Kullback–Leibler divergence between the
distributions $P$ and $Q$, and ${\cal Z}(x;T)$ is a normalization constant.
Again following Ref tishby+al_99 , we can turn a self–consistent equation into
an iterative algorithm. We have explored solutions with different cardinality
for $z_{t}$, and find that they form a hierarchical clustering scheme, as
expected from earlier work berman+al_16 . Our focus here is on the most severe
compression, in which $z_{t}$ has just two states, and the limit $T\rightarrow
0$, where the mapping $x\rightarrow z$ becomes deterministic. Further, we
search for compressions that preserve information from one state to the next,
corresponding to $\tau=1$ in Eq (1).

Figure 3: Behavioral states and their compressed representation. (a) The
original $N_{s}=123$ states and their boundaries, grouped informally into
clusters berman+al_14 . (b) The optimal compression into two states, with
$\sigma=+1$ shown in red and $\sigma=-1$ shown in grey.

When we compress into just two states we turn behavior into a binary sequence,
or equivalently a one–dimensional chain of Ising spins $z_{t}=\sigma_{t}=\pm
1$; the structure of this mapping is illustrated in Fig 3. Each point in the
two–dimensional behavioral space represents a short temporal sequence from the
original video recording; the probability distribution in this plane has many
well resolved peaks, and the dynamics consist of sojourns in the neighborhood
of one peak followed by a quick jump to another, allowing us to define the
$N_{s}=123$ discrete states and their boundaries berman+al_14 . States which
are neighbors in these two dimensions are similar by construction, and can be
(informally) grouped into the clusters shown in Fig 3a–anterior movements,
posterior movements, locomotion gaits, etc. In Fig 3b we have the results of
the optimal mapping $x\rightarrow z$, which groups together a range of
relatively rapid movements into one state $\sigma=+1$ while walking and idling
are mapped to $\sigma=-1$ difference .

Once we have reduced behavioral trajectories to a chain of Ising spins, we can
measure temporal correlations with the usual spin–spin correlation function
Pc_C ,

$C(\tau)=\langle\sigma_{t}\sigma_{t+\tau}\rangle-\langle\sigma_{t}\rangle^{2}.$
(4)

Perhaps surprisingly, although we asked for a compression that preserves
information only across one behavioral transition, Fig 4 shows that the
resulting binary or Ising variables have measurable correlations out to
$\sim$1000 transitions. While the compression into two behavioral ‘‘states’’
obviously throws away a lot of detail (by construction), long time
correlations are preserved.

Figure 4: Connected correlation function for the binary description of
behavior, Eq (4). Data points and error bars (red) as in Fig 1, compared with
a Markov model over the two states (green). Solid line is a single exponential
decay with correlation time $\tau_{c}=7.9\pm 0.9$.

Having reduced behavior to binary sequences, we would like to understand the
full probability distribution for these sequences, $P(\\{\sigma_{t}\\})$,
going beyond the Markov approximation. We want our approximation to match the
observed bias between the two behavioral states, measured by the expectation
value $\langle\sigma_{t}\rangle$. We also want to capture the long–ranged
correlations in Fig 4, so we insist that the correlation function $C(\tau)$
computed from the model $P(\\{\sigma_{t}\\})$ match the correlation function
that we observe experimentally. The minimally structured, or maximum entropy
model that matches $\langle\sigma_{t}\rangle$ and $C(\tau)$ is

$P(\\{\sigma_{t}\\})={\frac{1}{Z}}\exp\left[h\sum_{t}\sigma_{t}+{\frac{1}{2}}\sum_{t,t^{\prime}}\sigma_{t}J(t-t^{\prime})\sigma_{t^{\prime}}\right],$
(5)

where the interactions $J(\tau)$ must be tuned to match the correlation
function $C(\tau)$, and the field $h$ must be tuned to match the asymmetry
$\langle\sigma_{t}\rangle$ jaynes_57 ; bialek_12 . We solve this inverse
problem following the same strategy as in previous work tkacik+al_14 , using
Monte Carlo to estimate $C(\tau)$ in the model and adjusting $J(\tau)$ in
proportion to the difference between this estimate and the measured values.

Figure 5: Interactions $J(\tau)$ needed to reproduce the correlation function
$C(\tau)$ for $\tau<\tau_{\rm max}$, with $\tau_{\rm max}=100$ (red), $400$
(green), $1000$ (blue points with errors). Error bars are the standard
deviations across models that match $C(\tau)$ from randomly chosen halves of
the flies in the experimental population. Black dashed line is
$J(\tau)=1/\tau^{2}$, for comparison.

Even without detailed calculation, we know that generating long–ranged
correlations in one dimension (here, time) requires long–ranged interactions,
and this is what we find (Fig 5). Specifically, if we try to match the
correlations $C(\tau)$ for $|\tau|\leq\tau_{\rm max}$, then
$J(|\tau|>\tau_{\rm max})=0$, but as we increase $\tau_{\rm max}$ we
‘‘uncover’’ interactions that couple states with longer and longer
separations. We see no end to this, within the limits of our data. More
subtly, although the asymmetry in the two clusters of behavioral states is
large, corresponding to a ‘‘magnetization’’ $\langle\sigma_{t}\rangle=0.424\pm
0.048$, the magnetic field that we find is very small, $h=(4.7\pm 0.7)\times
10^{-3}$ when we match correlations out to $\tau_{\rm max}=10^{3}$.
Intuitively, long–ranged correlations imply that sequences of states are
moving collectively, and hence a small intrinsic bias is amplified.

Over more than two decades in $\tau$, the interactions are very nearly
$J(\tau)=J_{0}/\tau^{2}$. The Ising model with such ‘‘inverse square’’
interactions has a fascinating history, and played an important role in the
development of scaling ideas that presaged the full development of the
renormalization group anderson_84 . It is quite startling to see this model
emerge from the analysis of data on animal behavior.

Perhaps the most important prediction of our model is the probability of the
behavioral state at time $t$ given the states at all other times $t^{\prime}$,

$\displaystyle P\left(\sigma_{t}|\\{\sigma_{t^{\prime}\neq t}\\}\right)$
$\displaystyle=$ $\displaystyle{{\exp\left[h_{\rm
eff}\left(\\{\sigma_{t^{\prime}\neq
t}\\}\right)\sigma_{t}\right]}\over{2\cosh[h_{\rm
eff}\left(\\{\sigma_{t^{\prime}\neq t}\\}\right)]}}$ (6) $\displaystyle h_{\rm
eff}\left(\\{\sigma_{t^{\prime}\neq t}\\}\right)$ $\displaystyle=$
$\displaystyle h+\sum_{t^{\prime}\neq t}J(t-t^{\prime})\sigma_{t^{\prime}}.$
(7)

We can test this prediction by walking through the data and collecting all the
moments in time where $h_{\rm eff}$ falls into some small range, and averaging
the behavioral states over these times testh ; we should find
$\langle\sigma\rangle=\tanh(h_{\rm eff})$. To avoid any dangers of
over–fitting overfit , we assume that the system is described exactly by the
inverse square Ising model, $J(\tau)=J_{0}/\tau^{2}$, and fit the one
parameter $J_{0}$. Results are shown in Fig 6a.

The agreement between theory and experiment that we see in Fig 6a is very
good, and the best fit value of $J_{0}$ is consistent with the estimate of
$J(\tau=1)$ in Fig 5. The largest deviations between theory and experiment are
at extreme values of $h_{\rm eff}$, where behavior is even more nearly
deterministic ($\langle\sigma\rangle\rightarrow\pm 1$) than predicted by the
theory. The distribution of $h_{\rm eff}$ is strongly bimodal, as if the fly
were alternating between strong internal biases, and we see this in the time
course of $h_{\rm eff}$ in Fig 6b.

Figure 6: (a) Mean behavioral state at one moment in time vs the effective
field determined by states at other times, from Eq (7) with
$J(\tau)=J_{0}/\tau^{2}$, $J_{0}=0.54\pm 0.01$. Points are means across 59
flies, and error bars are standard deviations across random halves of this
population; line is $\langle\sigma\rangle=\tanh(h_{\rm eff})$. (b) A short
segment of $h_{\rm eff}$ vs time for one fly.

Equation 6 has the form of a behavioral state responding at each moment in
time to a bias $h_{\rm eff}$, but these biases are determined by behavioral
states at other times through Eq (7). Importantly, the full model in Eq (5)
can be rewritten in a form where behavioral states respond independently at
each moment in time to a completely internal, fluctuating bias $\phi_{t}$
bialek_20 ,

$P\left(\\{\sigma_{t}\\}\right)=\int{\cal
D}\phi\,P\left(\\{\phi_{t}\\}\right)\prod_{t}{{\exp\left[\sigma_{t}(h+\phi_{t})\right]}\over{2\cosh(h+\phi_{t})}}.$
(8)

The probability distribution of biases is given by

$\displaystyle P\left(\\{\phi_{t}\\}\right)$ $\displaystyle=$
$\displaystyle{1\over Z}\exp\left[-S\left(\\{\phi_{t}\\}\right)\right]$ (9)
$\displaystyle S\left(\\{\phi_{t}\\}\right)$ $\displaystyle=$
$\displaystyle{1\over
2}\sum_{tt^{\prime}}\phi_{t}K(t-t^{\prime})\phi_{t^{\prime}}-\sum_{t}\ln\cosh(h+\phi_{t}),$

and the kernel $K(\tau)$ is the inverse of the interaction $J(\tau)$,

$\sum_{t^{\prime\prime}}K(t-t^{\prime\prime})J(t^{\prime\prime}-t^{\prime})=\delta_{tt^{\prime}}.$
(11)

For models in which $J(\tau)\sim 1/\tau^{2}$, for large $\tau$ we have
$K(\tau)\sim-1/\tau^{2}$, so that the internal biases also must have
interactions that extend over long spans of time; there must be yet more
internal degrees of freedom which carry memory across these spans. To generate
the $\sim 1/\tau^{2}$ behavior exactly requires that the dimensionality of
these internal dynamics be effectively infinite bialek_20 . In a different
language, if we wanted to generate $\sim 1/\tau^{2}$ behavior exactly using a
hidden Markov model, the number of hidden states would have to be effectively
infinite as well. While not suggesting a specific mechanism, the explicit
$\sim 1/\tau^{2}$ interactions provide a simple description of the
phenomenology, the minimum required to match the observed correlations
$C(\tau)$.

To summarize, we have tamed some of the complexity of the behavior in walking
flies by compressing the behavioral states into binary variables which
preserve correlations from one state to the next. The resulting binary
sequence nonetheless captures long ranged correlations in the behavior, with
memory detectable out to thousands of transitions. The simplest model that
describes this behavior turns out to be an Ising model with almost perfectly
inverse square interactions, and this model gives excellent quantitative
predictions for the behavioral state at one moment in time given the
surrounding sequence of states. Although very simple, this description points
to an effectively infinite number of hidden states.

###### Acknowledgements.

This work was supported in part by the National Science Foundation through the
Center for the Physics of Biological Function (PHY–1734030), the Center for
the Science of Information (CCF–0939370), the NSF–Simons Center for
Quantitative Biology (DMS–1764421), and Grant PHY–1607612; by the National
Institutes of Health (NS104889); and by the Simons Foundation.

## References

  * (1) K von Frisch, Decoding the language of the bee. Science 185, 663–668 (1974).
  * (2) JL Gould, Ethology: The Mechanisms and Evolution of Behavior (WW Norton, New York, 1982).
  * (3) JL Lawson and GE Uhlenbeck, Threshold Signals. MIT Radiation Laboratory Series 24. (McGraw–Hill, New York, 1950).
  * (4) DM Green and JA Swets, Signal Detection Theory and Psychophysics (Wiley, New York, 1966).
  * (5) LC Osborne, SG Lisberger, and W Bialek, A sensory source for motor variation. Nature 437, 412–416 (2005).
  * (6) GJ Stephens, B Johnson–Kerner, W Bialek, and WS Ryu, Dimensionality and dynamics in the behavior of C. elegans. PLoS Comput Biol 4, e1000028 (2008)
  * (7) T Ahamed, AC Costa, and GJ Stephens, Capturing the continuous complexity of behaviour in C elegans. Nat Phys in press (2020).
https://doi.org/10.1038/s41567-020-01036-8.

  * (8) K Branson, AA Robie, J Bender, P Perona, and MH Dickinson, High–throughput ethomics in large groups of Drosophila. Nat Methods 6, 451–457 (2009).
  * (9) GJ Berman, DM Choi, W Bialek, and JW Shaevitz, Mapping the stereotyped behaviour of freely moving fruit flies. J R Soc Interface 11, 20146072 (2014).
  * (10) GJ Berman, W Bialek, and JW Shaevitz, Predictability and hierarchy in Drosophila behavior. Proc Natl Acad Sci (USA) 113, 11943–11948 (2016).
  * (11) AB Wiltschko, MJ Johnson, G Iurilli, RE Peterson, JM Katon, SL Pashkovski, VE Abraira, RP Adams, and SR Datta, Mapping sub–second structure in mouse behavior. Neuron 88, 1–15 (2015).
  * (12) JD Marshall, DE Aldarondo, TW Dunn, WL Wang, GJ Berman, and BP Ölveczky, Continuous whole-body 3D kinematic recordings across the rodent behavioral repertoire. Neuron in press (2020).
https://doi.org/10.1016/j.neuron.2020.11.016.

  * (13) A Mathis, P Mamidanna, KM Cury, T Abe, VN Murthy, MW Mathis, and M Bethge, DeepLabCut: Markerless pose estimation of user-defined body parts with deep learning. Nat Neurosci 21, 1281–1289 (2018).
  * (14) T Pereira, D Aldarondo, L Willmore, M Kislin, SS Wang, M Murthy, and JW Shaevitz, Fast animal pose estimation using deep neural networks. Nat Methods 16, 117–125 (2019).
  * (15) AEX Brown and B de Bivort, Ethology as a physical science. Nat Phys 14, 653–657 (2017).
  * (16) SR Datta, DJ Anderson, K Branson, P Perona, and A Leifer, Computational neuroethology: A call to action. Neuron 104, 11–24 (2019).
  * (17) TD Pereira, JW Shaevitz, and M Murthy, Quantifying behavior to understand the brain. Nat Neurosci 23, 1537–1549 (2020).
  * (18) MW Mathis and A Mathis, Deep learning tools for the measurement of animal behavior in neuroscience. Current Opin Neuro 60, 1–11 (2020).
  * (19) W Bialek, What do we mean by the dimensionality of behavior? arXiv:2008.09574 [q–bio.NC] (2020).
  * (20) N Tishby, FC Pereira, and W Bialek, The information bottleneck method. In Proceedings of the 37th Annual Allerton Conference on Communication, Control and Computing, B Hajek and RS Sreenivas, eds, pp 368–377 (University of Illinois, 1999); arXiv:physics/0004057 (2000).
  * (21) This result is different from that in Ref berman+al_16 , which used the (asymmetric) information bottleneck to find compressed representations that preserve the maximum information about the future state, $I(z_{t};x_{t+\tau})$. In the asymmetric case the optimal binary representation corresponds almost perfectly with moving vs idle.
  * (22) In comparing Figs 1 and 4, note that for binary variables we have ${\tilde{P}}_{c}(\tau)=C(\tau)/[1+\langle\sigma_{t}\rangle^{2}]$.
  * (23) ET Jaynes, Information theory and statistical mechanics. Phys Rev 106, 620–630 (1957).
  * (24) W Bialek, Biophysics: Searching for Principles. (Princeton University Press, Princeton NJ, 2012).
  * (25) G Tkačik, O Marre, D Amodei, E Schneidman, W Bialek, and MJ Berry II, Searching for collective behavior in a large network of sensory neurons. PLoS Comput Biol 10, e1003408 (2014).
  * (26) PW Anderson, Basic Notions of Condensed Matter Physics (Benjamin/Cummings, Menlo Park CA, 1984).
  * (27) For a parallel discussion of states in a network of neurons: L Meshulam, JL Gauthier, CD Brody, DW Tank, and W Bialek, Collective behavior of place and non-place neurons in the hippocampal network. Neuron 96, 1178–1191 (2017).
  * (28) Although we have a large amount of data, the model we are considering has 1000 parameters if we try to match correlations out to $\tau_{\rm max}=10^{3}$. Precisely because correlations have a long tail, it is not clear how many effectively independent samples we really have. Thus, to test our theoretical approach, we take the simple inverse square model seriously.