Chaotic heteroclinic networks as models of switching behavior in biological systems

Biological systems are invariably highly complex, involving very large numbers of components or substances interacting in complicated ways. Two examples are cortical and metabolic networks. Without major simplification, analytical studies of detailed biological models seem hopelessly out of reach. Not all biological models can be dimension reduced, but when a system has a finite number of dominant states that are semi-stable, phenomenological models focusing on transitions between these states can be constructed and analyzed. One way to reduce to this more tractable situation is to sufficiently constrain the circumstances or scope of study, thereby limiting the set of relevant behaviors. When they exist, dominant states may be identified through empirical observation. They can also be deduced from methods of data-driven modeling such as PCA, HMM, Koopman approximation, or machine learning [9, 23, 46, 2, 18, 44].

This paper is about dynamical models of systems with a finite number of semi-stable “attractor states" around which the system stabilizes briefly before transitioning to another state. We are interested in situations where trajectories are not periodic but seemingly random; which state will occur next has the appearance of being unpredictable, as do transition times. Switching dynamics of this kind have been observed in many biological settings (see Discussion). The attractor states in question typically represent fundamental behaviors or processes [18, 11, 48, 41, 37], and the branching is often associated with decisions. Our motivating examples are primarily from biology, but we consider here general dynamical systems with the characteristics above.

Continuous-time Markov chains (the states of which correspond to the attractor states above) are the simplest natural models that come to mind for describing switching dynamics, but the Markovian assumption is strong. Biological events often have some degree of history dependence. What an animal will do next, for example, may depend on what it did last and how long it has been in its present state. It would be useful to have tractable models that permit a broader set of dynamical behaviors.

We propose that chaotic heteroclinic networks are excellent phenomenological models of switching dynamics. Given a homoclinic loop or heteroclinic cycle, it has been known since the days of Poincaré that transversal intersections of stable and unstable manifolds give rise to complicated behavior[32, 34, 35]. The qualitative picture is clarified by Smale, who showed that this complicated behavior included the existence of horseshoes [42]. Here we take these ideas one step further. We demonstrate that by manipulating the eigenvalues of the saddle fixed points in a heteroclinic network together with the geometry of stable and unstable manifolds and their intersections, one can impose quantitative control on the switching dynamics, including branching probabilities, dwell times, and other more detailed characteristics. These new tools have enabled us to design large classes of chaotic dynamical systems and to tailor transition statistics in heteroclinic networks in order to match observed behavior in biological modeling.

As proof of concept, we constructed chaotic heteroclinic networks to reproduce experimental C. elegans data. C. elegans are a well-studied model organism in biology due to their simple anatomy, simple nervous system, few observed behaviors, and the relative ease of performing experimental measurements of behavioral and neural activity [18, 30, 36, 40, 47]. It has been shown that such activity can be represented in low-dimensional PCA space where semi-stable states correspond to stereotyped behaviors (or fictive behaviors) with characterizable probabilities of transition between these states [18, 31, 24]. C. elegans transition between different behaviors seemingly at random; their tendency to reside in various dynamical and behavioral states can be modified with experimental conditions, such as oxygen levels, genetic strain, and developmental stage [31]. The stochastic switching dynamics of C. elegans, modulated by experimental conditions and corresponding to different locations in neural activity space, is an ideal system to demonstrate how heteroclinic networks can generate stochastic switching behavior observed in biological systems.

Researchers have modeled C. elegans’ activity with a variety of different paradigms. Markov models can accurately capture the observed switching statistics but do not shed light on how transitions are generated [39, 14, 24]; it has also been observed that target states vary with dwell times, pointing to the need for more detailed modeling [31]. Previous dynamical systems models of C. elegans activity used control inputs as the mechanism for inducing state changes [13, 29], but state changes in C. elegans appear to occur mostly spontaneously and not in direct response to external stimuli. The model we propose has the capability to address all of these issues.

Seemingly random switching from one pattern of behavior to another with no apparent trigger is a motif seen throughout biology and neuroscience. We propose that chaotic heteroclinic networks with quantitative control of transition properties may be suitable phenomenological models for this type of behavior.

This section contains the theory part of the paper. We begin by reviewing some basic ideas connected with heteroclinic dynamics. Then we specialize to 2D, discussing first the case of flows before perturbing their time-$t$-maps to produce branching behavior in chaotic heteroclinic networks, the main objects of interest in this paper. This is followed by examples of quantitative control on the switching dynamics, techniques that will be used to build phenomenological models of C. elegan behavior.

Consider an ordinary differential equation (ODE)

$\displaystyle\dot{x}=f(x),\quad x\in\mathds{R}^{d},$

(1)

where $d\geq 2$ is an arbitrary integer and $f$ is assumed to be smooth. The flow associated with Eq. 1 is denoted by $\varphi_{t}$, so that if $x(t)$ is the solution of the ODE with initial condition $x(0)=x_{0}\in\mathds{R}^{d}$, then $\varphi_{t}(x_{0}):=x(t)$ is the trajectory of the flow starting from $x_{0}$. Let $p$ be an equilibrium point of $\varphi_{t}$ i.e. $f(p)=0$. We say $p$ is a hyperbolic fixed point if none of the eigenvalues $\lambda$ of $Df(p)$ lies on the imaginary axis. If the real parts of $\lambda$ are $<0$ (resp. $>0$) for all $\lambda$, then $p$ is a sink (resp. a source); if some are $>0$ and some are $<0$, then $p$ is a saddle fixed point. Through each hyperbolic fixed point of saddle type are stable and unstable manifolds, defined to be

	$\displaystyle W^{s}(p)$	$\displaystyle=\{x\in\mathds{R}^{d}:\varphi_{t}(x)\rightarrow p\ \ \text{as}\ \ t\rightarrow\infty\},$		(2)
	$\displaystyle W^{u}(p)$	$\displaystyle=\{x\in\mathds{R}^{d}:\varphi_{t}(x)\rightarrow p\ \ \text{as}\ \ t\rightarrow-\infty\}$		(3)

respectively. A point $x$ is called a homoclinic point of $p$ if $x\in W^{u}(p)\cap W^{s}(p)$. If $p_{1}$ and $p_{2}$ are distinct hyperbolic fixed points of saddle type, then $x$ is a heteroclinic point from $p_{1}$ to $p_{2}$ if $x\in W^{u}(p_{1})\cap W^{s}(p_{2})$, i.e., $\varphi_{t}(x)\rightarrow p_{1}$ as $t\to-\infty$ and $\varphi_{t}(x)\rightarrow p_{2}$ as $t\to\infty$. A point $x$ is called a tranverse homoclinic point of $p$ if $W^{u}(p)$ meets $W^{s}(p)$ transversally at $x$. Transverse heteroclinic points are defined similarly.

In discrete time, i.e., for dynamical systems generated by iterating a map $F:\mathds{R}^{d}\to\mathds{R}^{d}$ which we assume to be smooth and invertible, the corresponding objects are defined analogously: e.g. a fixed point $p\in\mathds{R}^{d}$ is called a hyperbolic fixed point of saddle type if the spectrum $\Sigma$ of $DF(p)$ is such that $\Sigma\cap\{|z|=1\}=\emptyset$ and $\Sigma\cap\{|z|<1\},\ \Sigma\cap\{|z|>1\}\neq\emptyset$. Stable and unstable manifolds are defined as above with $n\in\mathbb{Z}$ in the place of $t$, and a point $x\in W^{u}(p_{1})\cap W^{s}(p_{2})$ is called a heteroclinic point from $p_{1}$ to $p_{2}$.

The existence of transverse homoclinic points is well known to lead to complicated dynamics, an observation that goes back to Poincaré. The qualitative picture is developed by Smale, who showed that homoclinic points are accompanied by the existence of horseshoes, a hallmark of chaos. For a more systematic treatment, see Refs. [32, 33, 26, 16].

The continuous-time picture in $\mathbb{R}^{2}$ is simple because for saddle fixed points, $W^{u}$ and $W^{s}$ are 1D, and if $x$ is a heteroclinic point from $p_{1}$ to $p_{2}$, then the orbit of $x$, $\mathcal{O}(x):=\{\varphi_{t}(x),t\in\mathbb{R}\}$, coincides with one of the branches of $W^{u}(p_{1})$ and one of the branches of $W^{s}(p_{2})$. That is, $\mathcal{O}(x)$ is a saddle connection from $p_{1}$ to $p_{2}$.

A sequence of saddle equilibrium points, $p_{1},p_{2},\cdots,p_{n}$, $p_{n+1}=p_{1}$ is said to form a heteroclinic cycle if there is a saddle connection $\Gamma_{i}$ from $p_{i}$ to $p_{i+1}$, $i=1,2,\cdots,n$. Heteroclinic cycles are prevalent in population dynamics [17, 1] and occur naturally in low codimension bifurcation theory [16]. Fig. 1(a) shows an example of a 3-cycle. Orbits starting from initial conditions sufficiently near $p_{1}$ in the domain bounded by $\Gamma:=\cup\Gamma_{i}$ will shadow the piecewise smooth curve $\Gamma$ for some time. Additionally, if the heteroclinic cycle is stable, then $\varphi_{t}(x)$ will be attracted to $\Gamma$ as $t\to\infty$, in much the same way that nearby orbits are attracted to limit cycles — except that as $t$ increases, the orbit spends a larger and larger fraction of time near the three fixed points.

A sufficient condition for the stability of a heteroclinic cycle connecting saddle fixed points $p_{1},\cdots,p_{n}$ is

$\displaystyle\frac{-\lambda_{s}^{(i)}}{\lambda_{u}^{(i)}}>1\qquad\mbox{for all }i$

(4)

where $\lambda_{s}^{(i)}$ is the eigenvalue in the stable direction of $Df(p_{i})$ and $\lambda_{u}^{(i)}$ is the eigenvalue in the unstable direction. Condition (4) implies that compression is greater than expansion at each $p_{i}$. This condition can in fact be relaxed to requiring only that the product of the ratios $-\lambda_{s}^{(i)}/\lambda_{u}^{(i)}$ be $>1$. See e.g. Refs. [22, 21, 27].

More general than heteroclinic cycles is the idea of a heteroclinic network, loosely defined to be the union of a finite collection of saddle fixed points joined by saddle connections. An example of such a network is shown in Fig. 1(b). For further discussion and examples, see e.g. Refs. [20, 26, 15, 12].

It is easy to see that in a heteroclinic network on $\mathbb{R}^{2}$, the phase space is divided into cells bounded by heteroclinic cycles. Imposing stability conditions on each cycle as above, the dynamical picture can be described as follows: Each orbit is trapped in a cell; it is attracted to the heteroclinic cycle that bounds the cell, lingering at each one of the saddle fixed points which act as metastable states before moving to the next. The set of metastable states visited depends on initial conditions, but they are visited in a cyclical order without exception. Dwell times at each metastable state, hence the time to complete each cycle, increase without bound as time goes to infinity.

Heteroclinic networks for 2D flows as described above are already models of switching behavior, but the dynamics they describe are very special. In particular, branching cannot occur, and the dynamics are not chaotic.

A number of authors have introduced random noise to heteroclinic networks of the type above to induce branching behavior. Ref. [37], for example, constructs dynamical systems containing such networks to model orbits moving between metastable cognitive states in the brain. Ref. [6] used small noise perturbations near saddle fixed points to model decision making. There is a fair amount of rigorous theory on randomly perturbed heteroclinic networks [3, 5, 4, 7]. We postpone further discussion of these results as our approach is orthogonal: The seemingly random switching behavior in our model comes not from the use of random noise but from chaotic dynamics as we now describe.

Chaotic behavior in our model comes from the nontrivial intersection of stable and unstable manifolds in heteroclinic networks. As discussed above, flows in 2D cannot support chaotic dynamics. The lowest phase dimension for which chaotic heteroclinic dynamics can occur is 3D for continuous time and 2D for discrete time. For simplicity, we will work with the latter, though many of the ideas are not confined to two phase dimensions.

Specifically, we start from a flow $\varphi_{t}$ on $\mathbb{R}^{2}$ with a heteroclinic network, fix a small number $t_{0}>0$, and let $F_{0}$ be the time-$t_{0}$-map $\varphi_{t_{0}}$. It is an easy fact that the saddle fixed points and their stable and unstable manifolds for $F_{0}$ are identical to those of $\varphi_{t}$. We describe below a procedure for modifying $F_{0}$ along a saddle connection to obtain a smooth map $F:\mathbb{R}^{2}\to\mathbb{R}^{2}$ with nontrivial intersection of stable and unstable manifolds. The qualitative ideas behind this procedure are standard in the dynamical systems literature [32], but the quantitative control to follow are, to our knowledge, novel.

Let $p_{1}$ and $p_{2}$ be saddle fixed points. We assume there is a heteroclinic orbit $\Gamma_{1}$ of $\varphi_{t}$ connecting $p_{1}$ to $p_{2}$. For simplicity, we assume, as in Fig. 2(a), that $p_{1},p_{2}$ and $\Gamma_{1}$ lie on a horizontal line, while $W^{s}(p_{1})$, $W^{u}(p_{2})$ are vertical lines. For the map $F_{0}$, a fundamental domain of $W^{u}(p_{1})$ on $\Gamma_{1}$ is a segment $S\subset\Gamma_{1}$ with the property that the left endpoint of $S$ is mapped to the right endpoint under $F_{0}$, so that $F_{0}(S)\cap S=\emptyset$ modulo end points and $\Gamma_{1}=\cup_{n\in\mathbb{Z}}\ F_{0}^{n}(S)$.

Now fix a ball $B$ centered at the midpoint of $S$ with radius half the length of $S$, and let $F$ be such that

	$\displaystyle F$	$\displaystyle=F_{0}\qquad$	$\displaystyle\text{outside of }B,$
	$\displaystyle F$	$\displaystyle=F_{0}+g\qquad$	$\displaystyle\text{on }B$

where for all $x\in B$, $F_{0}(x)+g(x)$ is assumed to lie directly above or below $F_{0}(x)$. That is, we perturb $F_{0}$ into $F$ by pushing the images of points in $B$ up or down. See Fig. 2(a), where $F(S)$ has the shape of a biased sinewave. Since $F=F_{0}$ outside of $B$, the $F^{n}$-images of $F(S)$ for $n=1,2,3,\cdots$ are stretched vertically, compressed horizontally and pressed against $W^{u}(p_{2})$ as $n$ increases.

Let $W_{F}^{u}(p_{1})$ denote the unstable manifold of $F$ through $p_{1}$ (to distinguish it from the corresponding object for $F_{0}$), and let $W_{F}^{s}(p_{2})$ denote the stable manifold of $F$ through $p_{2}$. Using shorthands such as $S=[x_{1},x_{2}]$ to refer to the horizontal segment from $x_{1}$ to $x_{2}$, and letting $F(x_{2})=F_{0}(x_{2})=x_{3}$ (see Fig. 2(a)), it is easy to check that

	$\displaystyle\mbox{right branch of }W_{F}^{u}(p_{1})$	$\displaystyle=$	$\displaystyle[p_{1},x_{1}]\ \bigcup\ \cup_{n\geq 0}F^{n}(S)\ ,$
	$\displaystyle\mbox{left branch of }W_{F}^{s}(p_{2})$	$\displaystyle=$	$\displaystyle[x_{2},p_{2}]\ \bigcup\ \cup_{n\geq 0}F^{-n}([x_{2},x_{3}])\ .$

Thus, if $g$ is nontrivial, meaning if $F$ maps some points in $S$ above $\Gamma_{1}$ and some points below, then $W_{F}^{u}(p_{1})$ and $W_{F}^{s}(p_{2})$ intersect nontrivially and branching behavior is created: Points in $B$ that are mapped to locations above $\Gamma_{1}$ will eventually follow the up-branch of $W_{F}^{u}(p_{2})$, while those mapped to locations below $\Gamma_{1}$ will follow the down-branch of $W_{F}^{u}(p_{2})$.

The chaotic behavior or sensitive dependence on initial condition that ensues stems from the fact that for a point $x$ near $p_{1}$, it is hard to predict whether its trajectory under $F$ will follow the up or down-branch of $W_{F}^{u}(p_{2})$. For $x\in\Gamma_{1}$ strictly to the right of $p_{1}$, $F^{n}(x)$ will be in $S$ for some $n>0$. In the scenario depicted in Fig. 2(a), if it lies on the left two-thirds of $S$, it will follow the up-branch, and if it lies on the right third, it will go down. But even though its future is predetermined, the precise location of $F^{n}(x)$ is hard to predict when $x$ is very closer to $p_{1}$ and $n$ is large. For $x\not\in\Gamma_{1}$ near $p_{1}$ and on the right side of $W^{s}(p_{1})$, due to the vertical compression and horizontal expansion near $p_{1}$, there will be an $n>0$ such that $F^{n}(x)$ is very closer to $S$. The same reasoning as above then applies, with the precise location of $x$ playing a role in the “decision" at $p_{2}$.

Fig. 2(b) elucidates the branching properties that result from the perturbation described in Fig. 2(a). Fixing a disk around $p_{2}$, we define dwell time near $p_{2}$ to be the duration of time a trajectory spends in the disk. Observe that points that are farther from $\Gamma_{1}$ will have shorter dwell times than those closer to $\Gamma_{1}$. In this example, the distribution of dwell times for trajectories that eventually follow the up- and down-branches are different due to the up-down bias in the map $g$. A histogram showing the heights of $\{F(z_{i})\}$ for a large collection of points $\{z_{i}\}$ evenly spaced on $S$ is also shown. The higher probability to be above $\Gamma_{1}$ (hence to follow the up-branch of $W^{u}_{F}(p_{2})$ is evident, as is the fact that among those $F(z_{i})$ above $\Gamma_{1}$, the center of mass is farther from $\Gamma_{1}$, resulting in shorter dwell times.

Consider a line segment $S^{\prime}$ slightly above $S$: It is easy to see that the height distribution of $F(z_{i})$ and dwell times (not shown) will be even more biased than starting from $S$. For $S^{\prime}$ higher still, all trajectories will follow the up-branch. A similar analysis can be carried out for all points near $S$.

Two other examples of $g$ are shown in Fig. 2(c). The top and bottom panels show transition statistics for trajectories perturbed with a sawtooth perturbation versus a square wave perturbation. (One can approximate these functions with smooth maps without affecting too much the distributions of interest.) The sawtooth distribution results in the heights of $F(z_{i})$ being uniformly distributed, whereas the square wave distribution pushes some points far from $\Gamma_{1}$ and leaves others close, resulting in distinct clusters. We illustrate also, in the right panels (c.1-c.4), the effect of the unstable eigenvalues on dwell times. These plots confirm that the shapes of dwell time distributions are determined by the height distributions of $F(z_{i})$ as well as the eigenvalues at the saddle fixed point: weaker expansion at $p_{2}$ pushes points away more slowly, thereby increasing dwell times.

C. elegans, nematodes that grow to about $1$mm in length, are a well-studied model organism in neuroscience [47, 10, 19, 45, 18]. Whole-brain calcium imaging techniques record the activity of most C. elegans’ neurons while real or fictive behavioral states are simultaneously observed [30, 18, 38]. Because their movements are mediated by waves of contractions of bands of muscles that run the full length of the body, C. elegans exhibit only a limited set of behaviors, predominantly forward crawling, turns, reversals, and quiescence [18, 31, 24]. The low-dimensional neural dynamics of C. elegans cluster into different regions in PCA space that correspond to distinct behaviors [31, 18]. Experimentalists have fit Markov models to the behavioral sequences observed, treating transitions between behavioral states as well as the dwell times in each state as stochastic [31, 24]. Variations of behavior patterns in different experimental conditions have also been documented [24, 31].

In the last two sections of this paper we build dynamical systems in the form of chaotic heteroclinic networks (as described in Section II) to simulate C. elegans data. Results reported from two different labs will be used to challenge the model. In this section we focus on results from Ref. [31].

Ref. [31] considered primarily four types of behaviors, abbreviated as forward, reversal, turn, and quiescence. Treating these four behaviors as discrete states, we construct a chaotic heteroclinic network (Fig. 3(a)) containing three saddle fixed points that correspond to the forward, quiescence, and reversal states, and model the “turn" state by a heteroclinic orbit connecting the reversal to the forward fixed point: turns in the data seem always to be preceded by reversal and followed by forward movements. As discussed in Section II, we first build a (nonchaotic) heteroclinc network (in the form of a 2D flow) with the properties above. The chaotic network is obtained by applying perturbations to the flow-map along relevant heteroclinic orbits. Construction details are provided in the Appendix. Dwell times are defined to be time spent in a specified region around each one of the saddle fixed points or a portion of the heteroclinic orbit in the case of turns.

For illustration, we pick one of the modifications of the flow-map, namely the one illustrated in Fig. 3(b), and elaborate on the considerations that went into our choice of the perturbation function, $g$, in the region shown. Experimental data show that following reversal, the turn state occurs with very high probability ($0.95$ according to the Markov matrix in (d)). Assuming that the quiescence state is not dominant (which is often the case), most orbits arriving at reversal are from the upper branch of the stable manifold $W^{s}(p)$ where $p$ is the fixed point corresponding to reversal. These two facts together suggest that $g$, the function that determines how the upper branch of $W^{s}(p)$ meets the unstable manifold from the fixed point corresponding to the forward state, should have a strong bias. Accordingly, we have chosen a perturbation that pushes most of the points to the right of $W^{s}(p)$ and eventually into the turn state. Ref. [31] did not include data on how to constrain the shape of $g$, giving only bias probabilities. In general, the shape of $g$ can be deduced from measured dwell time distributions. The $g$ shown in Fig. 3(b), e.g., produces a bimodal distribution of displacements (or deviations from $W^{s}(p)$) one step later; see Section II for more detail. Such a distribution (Fig. 3(b)) implies that a large fraction of the trajectories is pushed quite far from $W^{s}(p)$ and will have shorter dwell time at $p$, while a smaller fraction stays very close to it and will have longer dwell times.

To show how a typical trajectory in the constructed dynamical system may look, we create a behavioral state time series by labeling the system at each point in time with either one of the four identified states or as “in-between", referring to times not in any of the four sets. Transitions between states are assumed to be relatively fast, making the “in-between" category negligible. In the example shown in Fig. 3(c), we have cut out the “in-between" stretches showing only the transition dynamics between the four states.

A target Markov matrix containing transition statistics estimated from experimental data published in Ref. [31] is shown in Fig. 3(d). These statistics were used to guide our model construction. The corresponding statistics collected from our model are shown in Fig. 3(e). Similarity in the matrix entries shows that perturbation functions can be tailored to produce simulated data that closely match the target statistics.

We digress here to point out that unlike Markov models, the dynamical models we propose have the capability to encode memory of past events. The following are two examples to illustrate this flexibility: (1) From the discussion above, we see that most orbits from forward to reversal leave reversal on the outside of the heteroclinic orbit connecting reversal to forward. If we set the strength of the contraction during the time in the turn state to be weak, most orbits will remain on the outside of the loop at forward, making them more likely to switch back to reversal after a short visit to forward. A strong contraction along the heteroclinic orbit will pinch all orbits close to the stable manifold of the fixed point at forward, resulting in less correlation between past and future behaviors. (2) We have prescribed a function, $g$, for how the unstable manifold from forward meets the upper branch of $W^{s}(p)$. There is a corresponding function, $\hat{g}$, to connect the unstable manifold from quiescence to the lower branch of $W^{s}(p)$. By choosing $g$ and $\hat{g}$ to be different, one can cause dwell times and transition probabilities at reversal to depend on the previous step, i.e., whether the orbit came from forward or quiescence. Similar constructions can lead to dependencies on the steps before.

Another set of results from Ref. [31] that we would like to replicate concerns C. elegans’ transition statistics and dwell times under different experimental conditions: prelethargus vs lethargus, and $10\%$ vs $21\%$ oxygen. The first are developmental stages: It is known that during lethargus periods, C. elegans are more likely to transition to quiescent behavior and to remain in a quiescent state for longer stretches of time with brief periods of activity [31]. Nematodes in the prelethargus stage, in contrast, are much more active. Oxygen levels in the environment are also known to affect activity level: C. elegans in an atmospheric 21% oxygen environment are more likely to remain active than C. elegans in a 10% oxygen environment, which is amenable to sleep. This condition occurs naturally when worms socially aggregate, creating a preferred low oxygen environment [31]. Ref. [31] contains detailed information on dwell times at and transition probabilities from the forward state under different developmental and environmental conditions.

Fig. 4(a) reprints, with permission from AAAS, data from Ref. [31]. Each panel reports on data recorded under a different set of conditions. The heights of the bars in each histogram (ignoring colors for now) show the numbers of transitions out of the forward state within $0$-$3$, $3$-$30$ and $>30$ seconds after entering this state; they can be seen as a probability distribution of dwell times. Each bar in the histogram is further decomposed into a red and a blue portion representing the number of transitions to reversal (red) or quiescence (blue). The red and blue lines superimposed on the histogram indicate the fractions of transitions of these two types. In the middle two columns, we present the corresponding statistics collected from our heteroclinic network.

In the top row we fit transition statistics for C. elegans in a prelethargus state with a $10\%$ oxygen environment. In this developmental stage C. elegans are more active. It is observed experimentally that under these conditions, about half of the transitions out of forward occur in the first $3$ seconds. Moreover, when dwell time is $<3$ sec., the C. elegans are much more likely to switch from forward to reversal; the longer the dwell time in the forward state, the more likely the C. elegans will switch to quiescence. Overall, transitions from forward to reversal are much more prominent than transitions into quiescence, though there is a nonnegligible fraction of the latter due to the low oxygen condition. These data are reproduced in our dynamical systems model using the perturbation function $g$ shown in the rightmost panel. Comparing the middle two panels to the left, we see the strong quantitative resemblance of model outputs to experimental data.

Next we compare the top two rows, both of which are for the prelethargus state but at different oxygen levels. Transition statistics are qualitatively similar, except that at $21\%$ oxygen, transition probabilities to quiescence are a little bit lower as high oxygen levels promote activity in C. elegans. Our model is able to capture these relatively subtle distinctions using the perturbations shown.

Differences between the prelethargus and lethargus data (at 10% oxygen) are far more substantial. In contrast to the top row, dwell times in the forward state in the lethargus case (bottom row) are much longer, nearly half of them more than 30 seconds. As before, the longer the time to transition, the more likely it will go into quiescence. This translates into a very significant difference between the prelethargus and lethargus data: In the lethargus case, a good majority of the transitions are to quiescence. Fig. 4 shows excellent agreement between data and model in both cases.

We finish by observing that even though only four types of behaviors are considered, a $4$-state Markov chain is not sufficient for capturing the type of results in Fig. 4; our dynamical model with its tunable parameters offers a simple way to recreate these and other behavioral characteristics.

In other studies [24, 18], researchers have made more refined classifications of C. elegans behavioral states than those used to fit the model in the previous section. Fig. 5(a) is reproduced from Ref. [24] with permission. It shows a 2D projection (from high dimensional data) of trajectories segmented into states with characteristic dynamics shown in different colors for two different worms. Here the authors distinguished between 8 discrete behavioral states, among them sustained reversals (blue), dorsal and ventral turns (green and yellow), forward crawling (red and crimson) and the transition from forward back to reversal (orange and brown). As part of their study, they proposed an 8-state Markov model, the transition probabilities of which are shown in Fig. 5(c) (reproduced with permission). Our aim in this section is to demonstrate that a chaotic heteroclinic network model can be built to reproduce the more complex picture of C. elegans behavior reported in this paper.

As in the last section, we start by building a 2D flow, representing each of the eight states by either a saddle fixed point or a heteroclinic orbit. Connections between nodes are based on Figure 5(a),(c); a diagram is shown in Figure 5(b). The color code is as in Figure 5(a) and (c); one of the states (purple) appears twice as trajectory segments in that state are found in two distinct sets of circumstances. To keep the dynamics simple, we have opted to capture only the more significant transitions, omitting several that occur with very low probability. Modifications are then made on the time-$t$-map of the flow to create nontrivial intersections of stable and unstable manifolds resulting in a chaotic heteroclinic network. The perturbation functions used are chosen with the aim of reproducing the transition probabilities in Figure 5(c).

In addition to transition probabilities, Ref. [24], Fig. 6(b), provides some guidance on the duration of time spent in each state. Though numerical values are not provided, and it is impossible to deduce exact dwell times from these figures, they do offer clues on dwell times. We have incorporated these clues into our dynamical model through the choice of unstable eigenvalues at saddle fixed points and by adjusting the speeds of travel along heteroclinic orbits. By specifying different eigenvalues and perturbation functions for each fixed point, as well as specifying different speeds along the corridors connecting fixed points, we are able to tune the model to replicate both transition statistics and dwell times in the 8-state C. elegans behavioral data.

Model outputs of the resulting dynamical system are collected, and the likelihood of its trajectories visiting state $j$ immediately after state $i$ is tabulated and shown in Fig. 5(d). Though not perfect, these dynamical transition probabilities — generated by a purely deterministic dynamical system — show good resemblance to the stochastic transition probabilities in Fig. 5(c). Finally, we present in Fig. 5(f) three example trajectories from our dynamical network, in analogy with the three trajectories from Ref. [24] (Fig. 5(e)). The variable itineraries and patterns of transitions in the three trajectories shown exemplify the possibilities and unpredictable outcomes typical in chaotic dynamical systems.

Stochastic switching in biology. The paradigm of having a finite number of dominant states and switching between them in a seemingly random way in the absence of external stimuli is relevant beyond C. elegans. It has been observed in biological and neural activity in many organisms, and Markov models have been used to capture these behaviors. For example, mice have been found to have different behavioral modules, the switching dynamics between which can be modeled with a hidden Markov model [48]. Larval zebrafish have been observed to transition between discrete behavior states following a Markov model [41]. Fruit fly behavior is known to evolve along a low-dimensional attractor with periodic and nonperiodic behavioral sequences [8]. Other transitions are likely to be caused or facilitated by biological events, but when a mechanistic understanding is out of reach due to high complexities in the underlying biology, many authors have, as a first attempt, idealized such transitions as random. For example, activity in the ventromedial prefrontal cortex (vmPFC) of humans appears to transition stochastically through a series of discrete states that correspond to affective experiences; these state transitions have been modeled with an HMM [11]. The neural activity in the gustatory cortex appears to transition between discrete states stochastically, a fact some researchers have attributed to noise-induced variability [28]. On the phenomenological level, these and many other examples of switching dynamics have a great deal in common with those considered in this paper.

Dynamical vs. stochastic models. Seemingly random switching has been observed in dynamical models. For example, models of spiking neurons are known to support multiple semi-stable states visited by the system in ways that resemble random transitions [25]. It is also known that random switching in (purely deterministic) dynamical systems can be generated by the addition of random noise. Ref. [37], for example, posits that transient cognitive dynamics result from metastable cognitive states which, when linked together, produce stable heteroclinic channels. Sequential decision-making can be modeled with stable heteroclinic channels with noise added to the system. See also Refs. [4, 6, 5]. Differential dwell times near fixed points and random switching between heteroclinic cycles generated by random noise have been studied by a number of authors; see e.g. Refs. [43, 3].

As an alternative to Markov models or dynamical models perturbed by random noise, we propose in this paper the use of chaotic heteroclinic networks to model seemingly random switching behavior. In our heteroclinic network models, dynamical complexity alone — without the need for stochastic manipulations — produces the variability and unpredictable behavior observed. One argument in favor of dynamical models could be that they are more natural, in the sense that physical and biological phenomena are not really stochastic in nature. One could also argue that dynamical systems models are more flexible: In Sections III and IV, we have demonstrated that chaotic heteroclinic networks can be constructed to emulate C. elegans data in two influential papers. Our models accurately reproduce not only basic transition properties between semi-stable states, such as dwell times and branching biases (Figs. 3, 4 and 5), but more subtle behaviors, including exit distributions and dependence of transition probabilities on dwell time (Fig. 4). We have also explained how one could, through network design, influence the degree to which memory of prior events will be retained.

This is not to suggest that stochastic models cannot possess similar capabilities; they can, but the models will have to be somewhat more complicated than $n$-state Markov chains where $n$ is the number of behavioral states. We have shown that dynamical models that are relatively easy to construct are capable of recreating a range of behaviors, with excellent quantitative control achieved through tunable parameters. Hence we propose them as a viable alternative to purely stochastic models.

Remarks on global network construction. As discussed, we start by constructing a 2D flow which we perturb to obtain the desired switching behavior. Among the many ways to construct such a flow, we have found the following to be easy to tune and to analyze: First we lay down all the saddle fixed points in $\mathbb{R}^{2}$, then we connect them with stable and unstable manifolds for the 2D flow. Each stable (unstable) manifold is assumed to comprise a vertical and a horizontal segment with corners rounded (as in Fig. 5(b)). In a neighborhood of this system of “train tracks" of stable and unstable manifolds, we build the vector field to be contractive to trap all nearby orbits. Where stable and unstable manifolds cross in $\mathbb{R}^{2}$, “overpasses" can be introduced. This geometry does not impose constraints on the architecture of the network.

Extension to a larger class of dynamical models. For clarity of exposition, we have limited ourselves in this paper to heteroclinic networks that are perturbations of time-$t$-maps of 2D flows. There are many ways to enlarge this class of models. An obvious extension is to allow higher dimensions, so that the stable and unstable manifolds of the saddle fixed points can have dimension greater than one. In the rest of the Discussion, however, we would like to focus on a different extension, to a broad class of dynamical systems we call generalized heterclinic networks:

Continuing to consider switching behavior among a finite collection $\mathcal{S}$ of identifiable, semi-stable states, an extension that we think is natural — both from the perspective of modeling and from the dynamical systems theory point of view — is to consider internal dynamics within each of the states in $\mathcal{S}$. Specfically, the saddle fixed points in classical heteroclinic networks can be replaced by hyperbolic invariant sets such as hyperbolic periodic orbits, normally hyperbolic invariant tori supporting quasi-periodic behavior, or Smale’s horseshoes [42]. In the case of a periodic cycle, for example, trajectories approaching the cycle will be attracted to it, stay with it for some time, and then veer away after going around the cycle a seemingly random number of times. Stable and unstable manifolds of hyperbolic invariant sets are well defined [32] and their intersections will define the switching dynamics.