University of Amsterdam & VU University
Physics and Astronomy
Bachelor’s Thesis
Chaos and dimensionality
in wormspace
Author
Abel Jansma
Supervisor
Dr. Greg Stephens
Student ID
10433627
Second assessor
Dr. Jaap Kaandorp
July 7, 2015
”We’re physicists, we’re gonna do the simple thing first, even if it’s wrong.” G. Stephens
Abstract
In this research we looked at postures of C. elegans across di↵erent behaviours. A posture is represented as a linear combination of eigenworms: a set of stationary body shapes found through PCA/ICA of time series with temporal data of angles from individual worm segments over time. Almost all behaviour turns out to take place on a low (4- or 5-)dimensional manifold within wormspace: the shape space spanned by the eigenworms. We then look for the dynamical system that could underlie the dynamics on this manifold. As a first try, a least-squares linear model is learned from the data through multiple linear regression. Then, a visual Turing test is performed by comparing visualisations of real time series with simulated ones, and it is seen that some key dynamics are not captured by the linear model.
When learning a model for crawling behaviour, oscillations are seen but quickly decay and the worm ends up as a straight line in the origin of wormspace in about 20 seconds. Adding zero-mean Gaussian noise shows not to be the right way of driving the system and starts to dominate the dynamics. For the omega-turns, there is a clear contraction along the body, but it bears little similarity to observed turns. This still leaves a nonlinear system as a viable option, and to get a first indication of the type of system, estimates of it’s Lyapunov exponent and correlation dimension are made. Running an ensemble of simulations for di↵erent parameters we find a matrix of possible exponents with no clear plateaus. Most embeddings of the attractor seemed to have a positive largest Lyapunov exponent, but this could reasonably be attributed to stochastic elements in the data rather than actual dynamics.
Contents
1 Popular abstract 3 2 Introduction 4 3 Theory 7 3.1 Eigenworms . . . 7 3.2 The Model . . . 73.3 Lyapunov exponent from time series . . . 10
3.4 Correlation dimension . . . 12
4 Methods 13 5 Results 13 5.1 Linearising the worm . . . 13
5.2 Lyapunov exponents and Correlation Dimension . . . 19
5.2.1 The Lorenz system and the R¨ossler system . . . 19
5.2.2 Wormspace . . . 26
6 Discussion 31
1
Popular abstract
Vroeger werd wiskunde eigenlijk alleen op simpele problemen toegepast. Simpel betekent hier dat er maar een een kleine hoeveelheid dingen zijn die beschreven worden, en dat er weinig onderlinge interactie is. Een bowlingbal die valt onder zwaartekracht bijvoorbeeld, of een planeet die om een ster draait. Als je in dat laat-ste voorbeeld een tweede planeet toevoegd wordt het probleem meestal al te moeilijk om exact op te lossen.
Maar op een gegeven moment kwamen wis- en natuurkundigen er achter dat als je in plaats van een paar, een hele boel objecten met interacties bij elkaar gooit, dat je dan kan negeren wat elk individu doet, en dat het gedrag van het hele gezelschap soms weer eigenschappen heeft die goed beschreven kunnen worden met wiskunde. Met deze nieuwe manier van wiskunde toepassen konden ineens nieuwe, veel complexere systemen onderzocht worden. Dit begon in de natuurkunde, maar bleek toepasbaar op een steeds breder scala aan fenomenen. Al gauw vonden de technieken ook hun plek in de biologie. Met de komst van de computer konden deze systemen ineens ook gesimuleerd worden, en achter ogenschijnlijk zeer complexe fenomenen zoals vogelzw-ermen en netwerken van genen of hersencellen bleken relatief simpele regels verstopt te zitten.
In dit onderzoek proberen we iets te leren over de wiskunde die het kruipen van een heel klein (⇠1mm) wormpje beschrijft. In plaats van proberen te begrijpen wat er precies in zijn hersenen en spieren gebeurt, kiezen we voor een meer abstracte benadering en doen alsof het enige relevante de vorm van zijn lichaam is. Bovendien willen we niet ´e´en individu begrijpen, maar zijn we op zoek naar de algemene regels die bepalen hoe de vorm van een typische worm verandert tijdens bepaald gedrag.
Dan maken we nog een stap richting het abstracte: Op dezelfde manier dat elk muzikaal akkoord gezien kan worden als een optelling van een aantal noten, zo kan elke vorm die de worm aanneemt gezien worden als een optelling van een aantal worm-vormen die we de eigenwormen noemen. Echter, waar een piano twaalf ver-schillende fundamentele noten nodig heeft om alle akkoorden te maken, hebben wij maar vier of vijf fundamentele worm-vormen of eigenwormen nodig om alle wormen na te maken. Als we deze eigenwormen bij elkaar op tellen om een bepaalde vorm te
maken gebruiken we wel wat meer van de een dan de ander (alsof je bepaalde noten in een akkoord harder speelt); we geven ze als het ware een bepaalde weging, of volume mee. Omdat de wormen steeds van vorm veranderen, veranderen deze wegingen ook de hele tijd, en de manier waarop deze wegingen veranderen is het centrale fenomeen dat we bestuderen in dit onderzoek.
Er is een theorie die zegt dat leven zich op onstabiele punten verzamelt tussen chaos en orde, en dit wilde we graag uitzoeken voor de worm. Hier zijn een aantal tests voor die beschrijven hoe gevoelig het gedrag is voor een iets andere beginpositie, iets wat in de wiskunde gedefinieerd is als chaotisch. De tests leken te laten zien dat de worm inderdaad chaotisch is, maar er is ook onderzoek gedaan dat laat zien dat de tests niet perfect zijn voor korte waarnemingen met veel ruis, dus harde conclusies durven we niet te trekken.
2
Introduction
With the rise of the use of computers in bioinformatics, understanding the behaviour of an organism through a simulation of its neural network has become a more re-alistic endeavor, and a lot of progress has been made in C. elegans by for example the OpenWorm project (www.openworm.org). Since C. elegans has so extensively been studied, not only has it’s genome been sequenced, but also its entire neural connectivity network and many biochemical processes are known, and great results regarding insight in it’s behaviour have been booked by e.g. Gray et al. (2005).
However, it remains an open question if this is the right way to think about behaviour. Contrasting the neural network simulation's bottom up approach, argu-ments could also be made for a more abstract, top down approach.
First of all, a philosophical argument is proposed by Elsasser (1981) and Ulanow-icz & Kau↵man (2009). They base their arguing on the findings by Russel & White-head (1910) that show the equivalence between continuous mathematics and discrete operations in certain set theories. These sets had to satisfy a strong demand: they had to be perfectly homogeneous, which means that all elements of a set had to be fundamentally indistinguishible. This was the case for most sets of objects described by physics (e.g. electrons), which strongly justifies the use of deterministic laws and mathematics to explore the dynamics of certain phenomena.
that the true nature of the studied phenomena lies not at the underlying determin-istic, physical evolution, but rather at the higher up, stochastic behaviour of the system. Maybe, if interactions between heterogeneous classes can not be accurately and consistently described by deterministic laws, all the interesting information that we might want to learn about the system should be sought at a scale that is more complex or emergent.
Secondly, an economical argument can be made, the advantage of the top-down approach being the fact that all of the neural networks complexity could be just as well represented by a potentially much smaller set of more abstract observables, o↵ering the possibility of significant dimensionality reduction.
Here, following the work of Stephens et al. (2008), I try to construct a model for the undulatory locomotion in C. elegans using only temporal information about the posture. Changing the basis in which the system is described, the postures are decomposed into a linear combination of so called eigenworms (Stephens et al. 2008). Almost all (96 %) of the behaviour, both crawling and escaping, turns out to take place on a lower, 4- or 5-dimensional manifold in the space spanned by all eigenworms. Using these eigenworms as new coordinates, Stephens et al. (2008) were able to show that there are some very reproducable dynamics were visible in this low-dimensional wormspace, and their surprising results motivate us to investigate the dynamics of trajectories.
While not denying the biological relevance of, and insight to be gained from, sim-ulations of neural networks, this approach o↵ers the possibility of gaining information about the worm that is not apparent in purely biological simulations. If it turns out that the trajectories through wormspace are described by some dynamical system, the whole toolbox of the mathematical analysis of these systems can be applied to the worm and quantitative points can be made regarding complexity, stochasticity, determinism and chaos, all of which are not as easily defined in a purely biological sense.
Finally, what motivates this research is an exploration of the hypothesis that life evolves at, or to, the loosely defined edge of chaos, or more precicely, near critical points (Bialek & Mora 2011). This region, where complexity is at a maximum but stability at a minimum, marks the boundary between a frozen, ordered, crystallised state (in living systems often referred to as death) and a chaotic, seemingly noisy state (in living systems an uncontrollable, epileptic state of being). An analysis of the wormspace dynamics could give insight into where on the order-chaos scale the
locomotion of the worm could lie. Since this is a behaviour that should be reasonably ordered in order to e↵ectively transport the worm to a new position, but also has to be highly dynamic in the sense that rapid, drastic changes are necessary for an e↵ective escape response, a calculation of the lyapunov exponent could provide interesting results regarding the edge of chaos hypothesis.
The research started o↵ with constructing a linear model from the time series us-ing first order multiple linear regression over a part of the series. This approach was chosen to start o↵ experimenting with the simplest, most na¨ıve way to find equations of motion that characterise the time series. With this model and the original data in hand, both were visualised with the same custom visualisation software and a visual Turing test decided if the model passed as representative (a visual test turned out to be conclusive so no further quantitative analysis was needed).
When looking at larger ensembles of real time series, a relatively large variability in e.g. escape responses was seen, and the question arises if this variability should be attributed to inherent stochasticity or a chaotic dynamical system. To gain further insight into the possibly chaotic dynamics, attempts were made to come up with a value for the crawling behaviour’s Lyapunov exponent and correlation dimension, after verification on some known systems.
3
Theory
3.1
Eigenworms
The time series for the worms consisted of 20Hz measurements of the relative angle of 100 arbitrarily chosen and evenly spaced segments of the worm. The data cap-turing the worm’s postures thus consisted of a 100-dimensional array. This is not a very insightful representation of it’s behaviour, and a certain amount of meaningful dimensionality reduction would be nice. This is where a technique called Principal Component Analysis (PCA) can help us.
Given a set of n-dimensional datapoints, PCA can be used to find new coor-dinates that best capture the variance in the data. More precisely: by finding the eigenvectors of the data’s covariance matrix, it finds the (orthogonal) directions along which the data has the highest variance. The precise portion of the variance along a direction is then determined by it’s corresponding eigenvalue. In many situations, PCA is able to deduce important characteristic feautures from the data, and it is used a lot in e.g. computer vision, facial recognition and machine learning. For our purpose, when applied to the worm series, the found eigenvectors comprise 100 angles that form the ”spine” of a worm. These principal components of the time series, the eigenvectors of it’s covariance matrix, are referred to as eigenworms.
The eigenworms can then be used as new coordinates in the sense that every original datapoint can be expressed as a linear combination of eigenworms. However, if it turns out that there was little variance in some directions in the first place, some eigenworms can be neglected and we would be able to reduce the dimensionality of the shape space.
3.2
The Model
For the derivation of the least-squares solution to the multiple regression problem, we follow the reasoning from the book Pattern Recognition and Machine Learning (Bishop 2007). The time series D that we want to fit a model to is an N by 5 array of real-valued datapoints representing respective weights of 5 eigenworms (a1 5) over N observations:
D = 0 B B B B B B B B B B @ a1,1 a2,1 · · · a5,1 a1,2 a2,2 · · · a5,2 .. . ... . .. ... a1,N a2,N · · · a5,N 1 C C C C C C C C C C A
In assuming that this system is described by an order 1 linear model with Gaussian noise, we assume that the following holds for every timestep i and some parameter matrix W:
W· ai = ai+1+ ✏i+1
Where aj is the jth row of D, or the state of the system at t = j, and ✏j is a vector with random, zero-mean, Gaussian noise.
To construct the model, we want to find the W that maximises the likelihood function p(y| x, w, 2), defined as the probability of result y given current state x, system parameters w and variance 2.
We first look at coupling the time evolution of the first eigenworm a1 to the com-plete state of the system:
a1,t= W1· at 1+ ✏
Using the assumumptions about our noise distribution, we can then assume that the likelihood function is given by:
p(a1,t| at 1, W, 2) =N (a1,t| W1· at 1, 2)
Where N (y, | µ, 2) is the area below the value y of a Gaussian distribution function with mean µ and variance 2.
Generalising this reasoning to an N-dimensional vector y:
p(y| x, W, 2) = N Y n=1 N (yn,| µn, 2) = N Y n=1 1 (2⇡ 2)12 e 2 21 (yn µn) 2 !
Since the Log function is monotonically increasing we know that max(Log(p)) occurs at the same value for p as max(p), so we might just as well maximise the Log of the likelihood function.
Taking the natural log of the right hand side yields:
Log N Y n=1 1 (2⇡ 2)12 e 2 21 (yn µn) 2 ! = Log 1 (2⇡ 2)N2 exp ✓ 1 2 2 N X n=1 (yn µn)2 ◆! = 1 2 2 N X n=1 (yn µn)2+ Log ✓ 1 (2⇡ 2)N2 ◆ = 1 2 2 N X n=1 (yn µn)2 N 2 Log(2⇡ 2)
We can now substitute the values from our specific problem, take the gradient with respect to W and set it to zero to find the W that corresponds to the maximum of the likelihood function.
y = at x = at 1 µn= Wn· an,t 1 =) rW 1 2 2 N X n=1 (an,t Wn· an,t 1)2 N 2Log(2⇡ 2) ! = ( 12 N X n=1 (an,t Wn· an,t 1)an,t 1 = 0
the W for which these hold sets the least-squares solution for our linear fit to the time series, and can be easily done with e.g. Matlab.
3.3
Lyapunov exponent from time series
Given a dynamical system we can, by looking at the rate at which infinitesimally close initial conditions diverge, quantify the extent to which the system’s time evolution is sensitive to initial conditions. We define a dynamical system to be chaotic if arbitrarily close initial conditions diverge exponentially.
For some then, the following must hold:
| (t)| = e t| (0)| =) (t, (0)) = 1 tLog ⇣ | (t)| | (0)| ⌘
Where| (t)| denotes the distance between the two trajectories. Then, the Lyapunov exponent is defined as:
= lim
t!1 (0)lim!0 (t, (0))
One can imagine that in an n-dimensional system, the rates of divergence are di↵erent for di↵erent dimensions, giving rise to a spectrum of n distinct Lyapunov exponents. This is indeed the case and they can be seen as the rates at which the di↵erent axes of an infinitesimal sphere of initial conditions get stretched over time, resulting in an ellipsoid.
Here, however, we are only interested in the largest Lyapunov exponent since that is the one that determines to what extent a system is chaotic.
Rosenstein et al. (1993) propose an algorithm for finding the largest Lyapunov exponent based on the approach by Wolf et al. (1985). We start o↵ by looking at just one pair of nearby trajectories. The distance between the two points of a pair j as a function of timestep i is then proposed to be given by:
where Cj takes care of any normalisation needed due to non-unity dj(0). To get to a simpler equation involving lambda we take the logarithm.
Log(dj(i)) = Log(Cj) + j(i t)
For any pair of exponentially diverging points, we now end up with a straight line when plotting Log(dj(i)) vs. (i t), where the slope of the line is given by j. To get the largest Lyapunov exponent, we look at the average line:
y(i) = 1
thLog(dj(i)ij,
where we average over all pairs j, and define the largest Lyapunov exponent to be the slope of the least-squares linear fit to this line. Note that we neglected all directions but one, attributing all divergence to the largest exponent. This is justi-fied by assuming that even a small di↵erence in the largest two exponents quickly results in the largest of the two dominating the divergence. This is the case for all non-degenrate Lyapunov-spectra, but when the largest exponent value is actually degenerate, this assumption no longer holds.
A problem arises however when the system describing the dynamics is unknown, as is the case for the trajectories through wormspace. In that case, it’s impossible to artificially create two nearby trajectories and look at the way initial separations diverge. Furthermore, the attractor often has a complicated metric on it in it’s folded-up state and distances are not easily defined. Two topological theorems can help us overcome this difficulty.
First of all, we invoke Whitney’s Embedding Theorem that states that a generic map from an n-dimensional manifold to (at most) 2n+1-dimensional Euclidean space is an embedding, meaning that the image of the original manifold is completely un-folded inR2n+1, so no two points of the original manifold are mapped onto the same point in the image (i.e. the image is di↵eomorphic to the original manifold).
Now say we have a dynamical system with an attracting manifold, but don’t know the metric or topology of this attractor, it would be very nice to find such a map that embeds the attractor in Euclidean space. Finding this map is done using Takens’
theorem. It states that given a trajectory on the attractor, we can construct a set of new (Euclidean) coordinates, referred to as delay coordinates, and that this map is an embedding of the attractor. The delay coordinates are defined as follows:
Given a time series that contains one-dimensional data: {x(t)} = {x1, x2, x3, ..., xN}
We construct a vector-trajectory through delay-coordinate space as: xdelay(t) = (x(t), x(t ⌧ ), x(t 2⌧ ), ..., x(t m⌧ )),
which still has free parameters ⌧ , the delay time and m, the embedding dimension. This means that if we have a time series over a complicated, n-dimensional attract-ing manifold in our original shape space, we can simply construct a path in 2n+1 dimensional delay-coordinates and be assured of the fact that the map has preserved all relevant topological properties of the attractor (after all, it is a di↵eomorphism), but that the trajectory this time takes place in Euclidean space, and the rate of separation of nearest neighbours can be easily calculated using the Euclidean norm.
3.4
Correlation dimension
As we saw previously, the type of embedding to choose in order to perform the correct distance calculations for the Lyapunov exponent, is, amongst other things, determined by the attractor’s dimensionality. There are many ways to calculate the dimension of a fractal object, but for time series analysis, the correlation dimension often proves to be a good and easily calculated estimate. The algorithm we use to calculate the Lyapunov exponent (Rosenstein et al. 1993) needs the same distance calculations as the calculation of the correlation dimension, so we can get a feel for the dimensionality of the attractor without too much extra computation. The defi-nition of the correlation dimension makes use of the correlation integral:
C(✏) = lim N!1 1 N2 N X i6=j=1 ⇥(✏ |x(i) x(j)|,
where x(i) can be a vector of arbitrary dimension,|...| denotes the distance according to some metric, and ⇥(✏ ...) is the Heaviside step function. With this definition, we have a measure of the fraction of all pairs of points on the trajectory that lie within a ball of radius ✏ of each other. If we now look at how Log(C(✏)) scales with Log(✏), we define this proportionality to be the correlation dimension.
4
Methods
All the data that we used in this research was acquired by Stephens et al. (2008). The fundamental data they worked with were tracked videos of worms across di↵erent behaviours. Details about the camera tracking technique, worm preparation and laser heat-impulse can be found in the corresponding paper (Stephens et al. 2008). The videos were translated to time series for relative angles of 100 abritrarily chosen, equally spaced segments of the worm. They then applied PCA/ICA to these time series and were able to transform the data into the eigenworm basis. In this research we only worked with data already in the eigenworm-basis. Gaps in the data were filled by linear interpolation in Matlab, and representative series were chosen by looking at an ensemble of visualisations and choosing them by sight. All visualisations were created with Processing, a Java-based programming language.
5
Results
5.1
Linearising the worm
First, the visualisation software was created and tested on a selection of the data. Figures 17a and 17b redirect to videos of visualisations of the original time series comprising forward and backward crawling as well as full omega-turns. To trigger the escape response, the worms were hit on the head with a laser to give a short heat impulse at t = 10 seconds. Especially the crawling shows some noise and gaps in the data. To get rid of the gaps, linear interpolation was used to fill in intervals with no data.
Convinced that the visualisation accurately represented the original movements, we then started to construct a deterministic (i.e. noiseless) model based on the simplest kind of behaviour: undisturbed forward crawling. Given a time series, a least-squares solution is obtained for parameter-matrix W. New time series were
then generated by deterministically evolving the original initial state of a series with W. Results of this can be seen in the video that figure 17c links to.
It can immediately be seen that all states show damped oscillations around the straight worm that correspond to crawling behaviour. These oscillations quickly (⇠20 sec) die out and the worm ends up in the origin of wormspace without any curvature. An example of this can be seen in figure 1.
Figure 1: A wormspace representation of the typical evolution of a linearised crawler.
This is somewhat expected since a deterministic linear model can only show 3 types of behaviour: exponential decay, exponential rise and oscillatory movement. Furthermore, when we look at the determinant and trace of di↵erent W matrices (see figure 2) we find that both are larger than 0 and a little less that those of the identity. Since the frequency of observation was much higher than the natural frequency of the dynamical system, we’d expect W to be close to the identity, since change per time step is small, and consecutive values are mostly correlated with themselves, rather than other components of the system. It’s also reassuring to note that these properties of W predict a stable system, which is what we see.
Figure 2: Determinants and traces of W-matrices based on 10 di↵erent worms. The chosen trajectories were through 4D-wormspace.
From the work of Stephens et al. (2008) we can see that there is an a1 a2 limit cycle that corresponds to forward/backward crawling. This o↵ers the opportunity to change the coordinates to ones that are easier to interpret.
In figure 3, real and simulated time series are looked at in new coordinates:
R = q
a2
1+ a22, the radius of the limit cycle = tan 1(a2
a1
), the phase of the oscillation in radians !, the di↵erence between two consecutive values for
0 100 200 300 400 500 600 -10
0 10 20
30
Real time series: ?
0 100 200 300 400 500 600 -0.4 -0.2 0 0.2 0.4
!
0 100 200 300 400 500 600 -10 0 10 20 a3 0 100 200 300 400 500 600 2 4 6 8 10 12 R 0 100 200 300 400 500 600 -40 -30 -20 -10 0Linearised model: ?
0 100 200 300 400 500 600 -0.6 -0.4 -0.2 0 0.2 0.4!
0 100 200 300 400 500 600 -4 -2 0 2 4 a3 0 100 200 300 400 500 600 0 1 2 3 4 RFigure 3: Both a real (left) and a linearised (right) worm trajectory in the new coordinates. The heat impulse was given at frame 200.
On the left hand side of figure 3 we clearly see a forward crawl (growing , pos-itive !), followed by a reversal right at the moment the heat impulse is given, and finally, after a peak in a3, a return to the initial state. Comparable dynamics can be seen in the linearised series, but also in these coordinates, the linearised time series only shows oscillations and decays, and clear reversals are absent. There are a few anomalies seen in and !, but these turn out to correspond to peaks in a3 and moments in time where is not well defined, so these anomalies are the results of our choice of coordinate definitions rather than actual properties of the model.
Next, the possibility of noise being the source of the extra dynamics is explored. Figure 4 shows the di↵erence between a value in the time series and a value predicted by the model based on the previous value. At first glance, these show some resem-blance to a Gaussian distrubution, which motvates driving the system with some Gaussian noise. When generating the new time-series, this time a random noise-term is added at each frame, drawn from a zero-mean Gaussian distribution with the same variance as the errors of the original deterministic series. Results seen in the video from figure 17d.
Now, we turn to the omega-turns. There is a relatively large evolutionary pressure on escape-responses, so it is worth considering that at least this part of the behaviour should be insensitive to initial conditions and very repeatable, creating the possibility of it being accurately represented by a linear model. Learning a model from time series that start at frame 300 (to somewhat isolate the omega-turn from the reversal), which is 5 seconds after the heat impulse, results in the simulated behaviour as seen in the video from figure 17e.
Even though an ensemble of these simulations obviously doesn’t pass a visual Turing test when placed next to real time series, it does show some interesting prop-erties. In many worms, an obvious contraction takes place along the whole body, which is reminiscent of the omega-turn of the escape response, and a few even show very life-like turns. In most, however, the contraction is induced in more that one place in the body, resulting in the body wave having half the wavelength of the real turn. Again, in all worms the only visible dynamics are decay and oscillation.
5.2
Lyapunov exponents and Correlation Dimension
The results from the previous section showed that linear models weren’t able to capture all the dynamics from the original time series. This leads to the conclusion that any accurate description is probably going to lie in the realm of non-linear systems. This, however, is so unconstraining that it motivates anyone wanting to find a description to first explore some of the non-linear properties of the system. Two of the most telling features of a non-linear system are it’s largest Lyapunov exponent and (if the system has an attractor) the attractor’s correlation dimension. 5.2.1 The Lorenz system and the R¨ossler system
We first started to see if we could reproduce the results for the Lorenz and R¨ossler systems as found by Sprott (2004). Plots of trajectories over the attractors and the reconstructed attractors in delay-coordinates can be found in figures 5 and 6.
For two versions of the Lorenz system with embedding dimension 3, a lag time of 16 and a time step of 0.005, we got Lyapunov exponents of 1.53 and 0.95 where Sprott (2004) found respectively 1.50 and 0.9056. For the R¨ossler system with em-bedding dimension 3, lag time 10 and time step 0.05, we arrived at 0.076, where
Sprott (2004) found 0.0714. As shown in the plots for the Lorenz system the precise values are sensitive to what interval in time is chosen to fit to. It can be seen that it takes a while for the points to reach constant exponential divergence, and when comparting this settling time to the typical oscillation time of figure 9, we see that both have a simliar time scale and we choose to ignore the time steps before which not a single typical oscillation has occured (i.e. the first 100 frames).
Further exploring (⌧, m)-parameter space gives us figure 11. When looking at the Lyapunov exponents, relatively consistent values are seen in the top left corners and middle region of the matrices of respectively the Lorenz and the R¨ossler system. Los-ing stability when we go to the right side of the matrix makes sense, since there we increase the lag time and therefore start decorrelating the dimensions. Both systems also show sensitivity to increasing the embedding dimension, the Lorenz system even more so than the R¨ossler system. Both plots show decreasing values for the exponent when increasing the embedding dimension or the lag time.
The opposite happens with the correlation dimension. They show a larger stable area where the values are in good agreement with Sprott (2004), but start increas-ing with embeddincreas-ing dimension and lag time. Note that the bottom right corner the Lorenz systems has no values for the correlation dimension, since combining a large embedding dimension with a large lag time creates the need for a very long trajectory, which for computational reasons was not used in that simulation. Still, all plots show an area where the values show some plateauing, and comparing these values with those of Sprott (2004) convinced us that we were accurately calculating the exponents and dimensions, so we moved on to the wormspace trajectories.
Figure 5: A 3D-plot of a trajectory over both the real and the reconstructed Lorenz attractor, showing visually similar topological structure. A fourth order Runge-Kutta technique was used to integrate the system. The used parameters were:
= 16, = 4, ⇢ = 45.92.
Figure 6: A 3D-plot of a trajectory over both the real and the reconstructed R¨ossler attractor, showing visually similar topological structure. A fourth order Runge-Kutta technique was used to integrate the system. The used parameters were:
Figure 7: Linear fits to the average Log of separation of trajectories over the reconstructed Lorenz attractor. Used parameters:
Figure 8: Linear fits to the average Log of separation of trajectories over the reconstructed Lorenz attractor. Used parameters:
Figure 9: Time evolution of the x coordinate of the Lorenz system. Shows a typical time scale of about 100 time steps. The used parameters were:
= 16, = 4, ⇢ = 45.92.
Figure 10: Linear fit to the average Log of separation of trajectories over the reconstructed R¨ossler attractor. Used parameters:
Figure 11: Distribution of Lyapunov exponents for Lorenz- and R¨ossler system. a = 0.2, b = 0.2, c = 5.7, t = 0.05
5.2.2 Wormspace
Figures 12 and 13 show an example of an original and a reconstructed attractor based on a1. Before we make a final choice about what embedding dimension to choose and how large the time lag for the delayed coordinates should be, we first look at an ensemble of di↵erent calculations resulting in the matrix plots of figure 14. It can immediately be seen that across the four individuals, the results proved to be very consistent. Whether this proves that we accurately isolated a general principle be-hind the behaviour, or just means that noise and autocorrelation artefacts dominate the outcome of the calculation, is not yet clear. It can also be seen that there is little indication for plateaus within the explored parameter-space. Where figure 11 shows that for both systems a large region in parameter-space has a pretty constant correlation dimension, figure 14 shows that this is not the case for the worms. Since values for correlation dimension and Lyapunov exponents turn out to be so depen-dent on the embedding dimension and the lag time, we have to be very careful to make claims about the findings.
Tanaka et al. (1996) showed that in completely random time series with no tem-poral correlation at all, you can find spurious positive Lyapunov exponents that are completely due to the finite series length and their stochastic nature, and have noth-ing to do with underlynoth-ing dynamics. In fact, these random trajectories can’t even be said to contain any dynamics, so the Lyapunov exponent isn’t well defined for them. Since we don’t know the signal to noise ratio of our data, we should keep in mind that our results could be mostly attributed to noise. To get some context, the same calculations were done for time series that contained a random shu✏ing of the time steps, destroying all temporal correlation between time steps. The results of these calculations are seen in figure 15. Comparing figures 14 and 15, similar patterns can be seen in matrices for real and shu✏ed data, indicating that a lot of the structure in the matrix plot of the real data might be due to inherent noise. The matrix for the correlation dimension of the shu✏ed data shows that there is indeed no temporal correlation, since there is no dependence on lag time. Also, the correlation dimension grows with the embedding dimension, showing that there is no structure that at one point gets completely unfolded, similarly to how a trajectory over a chaotic attractor can become e↵ectively stochastic when the attractor is embedded in too small an embedding dimension.
In figure 16, comparing distributions of Lyapunov exponent values across all worms, we see that both the real time series as well as the random shu✏es (which are not supposed to contain any dynamics) show a definite tendency to have a pos-itive largest Lyapunov exponent. This tendency is much stronger in the real series, but still this is further reason to be careful to make claims about the actual largest Lyapunov exponent and could be a result of the e↵ect observed by Tanaka et al. (1996).
Figure 12: A trajectory through a1 a2 a3 wormspace corresponding to forward crawling.
Figure 13: An axample of a reconstructed wormspace attractor, based on a1, with embed-ding dimension 12 (first three shown) and lag time 10.
Lyapunov exponent - Worm 1 lag time 5 10 15 20 25 embedding dimension 5 10 15 20 0 0.1 0.2 0.3 0.4 0.5 Worm 2 5 10 15 20 25 5 10 15 20 0 0.1 0.2 0.3 Worm 3 5 10 15 20 25 5 10 15 20 0 0.1 0.2 0.3 0.4 Worm 4 5 10 15 20 25 5 10 15 20 0 0.1 0.2 0.3 0.4 0.5
Correlation dimension - Worm 1
5 10 15 20 25 5 10 15 20 2 3 4 5 6 7 8 Worm 2 5 10 15 20 25 5 10 15 20 2 3 4 5 6 Worm 3 5 10 15 20 25 5 10 15 20 2 3 4 5 6 7 Worm 4 5 10 15 20 25 5 10 15 20 2 3 4 5 6
Lyapunov exponent - Worm 1 lag time 5 10 15 20 25 embedding dimension 5 10 15 20 -0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 Worm 2 5 10 15 20 25 5 10 15 20 -0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 Worm 3 5 10 15 20 25 5 10 15 20 -0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 Worm 4 5 10 15 20 25 5 10 15 20 -0.005 0 0.005 0.01 0.015 0.02 0.025
Correlation dimension - Worm 1
5 10 15 20 25 5 10 15 20 4 6 8 10 12 14 Worm 2 5 10 15 20 25 5 10 15 20 4 6 8 10 12 Worm 3 5 10 15 20 25 5 10 15 20 4 6 8 10 12 14 Worm 4 5 10 15 20 25 5 10 15 20 4 6 8 10 12 14
Figure 16: Distribution of Lyapunov exponents across four worms when varying embedding dimension and lag time in both real and shu✏ed time series. Generated from data in figures 14 and 15.
6
Discussion
As for the linear models not accurately describing the behaviour, there is the ques-tion if we could’ve done better. We could, for example, have chosen a di↵erent noise-distribution since the errors showed to be not perfectly Gaussian. We could also, especially for the omega-turns, have been more careful in selecting the data that the model used to learn from. Plots of a3 throughout the escape response show a sharp peak right around the moment the head touches the tail. With this informa-tion we could’ve created a much cleaner dataset of e.g. just the first part of the turn. Segmenting the behaviours in a meaningful way is also appealing in the light of find-ings of Gray et al. (2005) that led to their conclusion that ”reversal and omega turn behaviors have distinct final motor pathways”. It is possible that although specific parts of behaviour could be reasonably well linearised, we’re looking at a multiple of di↵erent behaviours which makes it impossible.
For the Lyapunov exponents and the correlation dimension, we first got results that didn’t agree with historical values. This turned out to be a result of integrating the equations using Euler’s method, and switching to a more precise fourth order Runge-Kutta method fixed this. Both the correlation dimension and the largest Lyapunov exponent are found by generating a linear fit to a measured signal, and the finiteness of our data and the presence of noise makes finding the right interval to fit to ambiguous. When running a simulation with a specific set of parameters, we saw that locating the interval of interest often had to be done by sight and no obvious indicators were present. Selecting intervals by hand is not a big problem when running a small number of simulations, but when generating the matrix plots of figure 14, we had to choose a specific interval to fit to that was the same for every simulation. This is a very crude way of estimating the slope and will certainly result in less accurate results. Because of this it was also not possible to come up with a reasonable metric for certainty and errors of the results, since we could not be sure if the variability seen in the matrices comes from ill-fitting or fundamental e↵ects of changing parameter values.
The most ambiguity is with the attractor-reconstruction. Whitney’s and Takens’ theorems give an upper bound on the minimum embedding dimension, and delay-coordinate time lag should not be very important, but these statements assume an
infinitely long, completely deterministic time series over an attractor with known fractal dimension. We, however, only have a finite time series, with stochastic ele-ments (noise), over an attractor for which no definite dimensionality is known. This forced us to explore intervals in parameter space and come up with a number of Lyapunov-exponents from which not one could be chosen as representative of the system.
Increasing the lag time showed to destroy the stability of the values of the expo-nents and dimensions, which is as expected. When you increase the lag time between di↵erent delay-coordinates, you start decorrelating the values and start losing dy-namics in the noise. However, the same is not expected to happen for increasing embedding dimension. When increasing embedding dimension, at a certain point the attractor is completely unfolded and all calculations should yield the same result regardless of adding extra dimensions to the problem. This is not what we’re seeing in our results and could inspire further investigation.
Lastly we note that results across individual worms proved to be very repeatable, which begs the question if this is also the case for mutants with a di↵erent genome. Given an accurate way to characterise the dynamics of certain worm behaviours, it could be very interesting to see how changing the genotype a↵ects the dynamics.
7
Conclusion
We wanted to try to understand the behaviour of C. elegans. This left us with two approaches: bottom-up, from the neural network and biochemical processes up to-wards muscle activity, or top-down, from the postures of the organism as a whole down to the underlying dynamical system. Although the first approach is certainly easier to interpret and shows very significant biology of behaviour, there are good arguments to be made for the second approach. Elsasser (1981), inspired by Russel & Whitehead (1910), came up with his term Biotonic law, referring to a law of nature that, contrary to physical law, is applicable even to hetergeneous classes of objects. These laws are to be arrived at by looking at larger scales of organisation, and more abstract phenomena, like postures of the worm as a whole rather than the underlying biology. We thus set out to learn a predictive model for spine-curvature based on time series of postures. These postures turned out to be very well describable by
linear combinations of so called eigenworms, enabling us to look at the dynamics of trajectories through wormspace and construct a linear model for di↵erent behaviours.
For regular crawling, the model was able to capture the bodywaves propagating along the spine, but these oscillations quickly decayed and every worm ended up as a straight, motionless line. When trying to explain this by the need of a noise-term to drive the system we saw that added noise with the same variance as the errors of the model starts to completetely dominate the movement and adds no interesting dynamics.
For the omega-turn of the escape response, the same technique yielded clear con-tractions along the body, bearing similartity to the actual turning, but these also soon ended up as straight, motionless worms.
Since we then knew that a linear model was not sufficient to account for all dynamics, we set out to determine whether the dynamics came from a chaotic or non-chaotic system, and went on to calculate the Lyapunov exponent. This was first done as verification of the algorithm with the known Lorenz and R¨ossler systems.
Then, the algorithm was applied to wormspace trajectories corresponding to undisturbed crawling. We ran an ensemble of simulations with di↵erent embed-ding parameters and ended up with a matrix of Lyapunov exponents that showed no plateauing. From this we could not come up with a single, representative exponent. Looking at the the mostly positive distribution of exponents makes it tempting to conclude that the actual Lyapunov exponent will also probably be positive, making the worm chaotic. However, temporally uncorrelated shu✏ed data also showed this tendency towards positive exponents.
Research has been done on Lyapunov exponents of random time series that showed that given enough non-deterministic, random data, spurious positive Lyapunov ex-ponents can be found, and Tanaka et al. (1996) warned to be very carefeul with concluding anything from these positive exponents. Respecting this warning and given the fact that there’s no obvious plateauing visible in figure 14, it seems impos-sible to conclude that the crawling behaviour has a positive Lyapunov exponent and is therefore chaotic.
(a) Visualisation of original crawling series. Also available at http://bit.ly/1eLkIj4.
(b) Visualisation of original escape response series with the heat impulse at t = 10 seconds.. Also available at http://bit.ly/1GBrSMH.
(c) Visualisation of linearised crawling. Also available at http://bit.ly/1dkm6aP.
(d) Visualisation of linearised crawling with noise. Also available at http://bit.ly/1HknMNN.
References
Bialek, W. & Mora, T. 2011, Journal of Statistical Physics, 144, 268 Bishop, C. M. 2007, Pattern Recognition and Machine Learning (Springer) Elsasser, W. M. 1981, Progress in Theoretical Biology, 6
Gray, J. M., Hill, J. J., & Bargmann, C. I. 2005, Proceedings of the National Academy of Sciences of the United States of America, 102, 3184
Rosenstein, M., Collins, J., & De Luca, C. 1993, Physika D, 65, 117
Russel, B. & Whitehead, A. N. 1910, Principia Mathematica (Cambridge University Press)
Sprott, J. 2004, Common Chaotic systems, http://sprott.physics.wisc.edu/chaos/comchaos.htm Stephens, G. J., Johnson-Kerner, B., Bialek, W., & Ryu, W. S. 2008, PLoS Comput
Biol, 4, e1000028
Tanaka, T., Aihara, K., & Taki, M. 1996, Physical Review E, 54, 117
Ulanowicz, R. E. & Kau↵man, S. A. 2009, A Third Window - Natural life beyond Newton and Darwin (Templeton Foundation Press)
