- Email: [email protected]
Fixed point analysis of nystagmus
Journal of Neuroscience Methods 161 (2007) 134–141
Fixed point analysis of nystagmus Maria Theodorou ∗ , Richard A. Clement Institute of Child Health, University College London, 30 Guilford Street, London WC1N 1EH, UK Received 16 February 2006; received in revised form 3 September 2006; accepted 29 September 2006
Abstract Motor disorders frequently contain a rhythmic component, but the associated oscillations are not usually precisely periodic. This lack of strict periodicity can make it difficult to identify the effects of experimental manipulations on the oscillation. In this report, we describe the application of a numerical technique for identifying fixed points of a nonlinear map to the recovery of underlying periodicities of the eye movement disorder of nystagmus. The technique is illustrated by application to two different types of nystagmus. In addition we use a local analysis of the behaviour at the fixed points to distinguish between different bifurcations in the two examples with changes in gaze angle. We conclude that the technique reveals consistent effects of experimental manipulations, which may be useful for quantitative characterisation of experimental and therapeutic manipulations of motor disorders. © 2006 Elsevier B.V. All rights reserved. Keywords: Nystagmus; Oscillations; Deterministic; Nonlinear dynamics; Fourier analysis; Wavelet analysis; Phase space analysis
1. Introduction Motor disorders frequently lead to the development of oscillations. In order to evaluate any treatment of the disorders quantitative techniques for describing the oscillations are required. An example of such an oscillation is the eye movement disorder of early-onset nystagmus. Early-onset nystagmus is an involuntary oscillation of the eyes which develops at birth or shortly afterwards, often within the first 6 months and persists throughout life. It is characterised by predominantly bilateral, conjugate, horizontal eye movements. It may be idiopathic, associated with a sensory abnormality, or more rarely, a neurological lesion (Harris, 1997). The main difficulty in analysing earlyonset nystagmus waveforms is that they are very variable, so that it has been difficult to obtain consistent results (Reccia et al., 1990; Clement et al., 2002; Optican and Miura, 2004). This is a common problem amongst motor disorders, so that effective techniques for analysing eye movement recordings should have widespread applicability. Eye movement recordings from subjects with early-onset nystagmus show a variety of different waveforms which can be described by combinations of jerk and pendular waveforms
∗
Corresponding author. Tel.: +44 207 905 2346. E-mail address: [email protected] (M. Theodorou).
0165-0270/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jneumeth.2006.09.025
(Dell’Osso and Daroff, 1975). A jerk waveform is comprised of an increasing velocity exponential drift of the eye away from the target direction followed by a quick movement which returns the eye to the direction of the target. A typical pendular waveform consists of an approximately sinusoidal oscillation. The oscillatory movement of early-onset nystagmus can be characterized by the Fourier components of the eye position signal (Dell’Osso et al., 1974; Dickinson and Abadi, 1985; Reccia et al., 1989, 1990; Abadi and Worfolk, 1991; Worfolk and Abadi, 1991; Clement et al., 2002). Typically there is a peak in the power spectrum at the fundamental frequency of the oscillation, together with harmonic components with smaller amplitudes. The main peak provides a good measure of the average frequency and amplitude of the early-onset nystagmus waveform. Fourier analysis has revealed that even pure pendular nystagmus, which is often described as sinusoidal, has harmonic components (Reccia et al., 1990; Abadi and Worfolk, 1991). A recent alternative to Fourier analysis is wavelet analysis (Miura et al., 2003), which allows localization in both the time and frequency domain. However, a much longer segment of data is required to get consistent results which can lead to changes in the waveform being obscured by attention changes in the subject (Dell’Osso, 2004). A more promising approach has arisen through the application of dynamical systems theory. The deterministic components of the nystagmus oscillations can be identified using a tech-
M. Theodorou, R.A. Clement / Journal of Neuroscience Methods 161 (2007) 134–141
nique developed in nonlinear dynamics. In this framework the behaviour of the neural system can be described by a map from one state to the next. The state of the neural system can be represented geometrically by a point; successive states form a sequence of points, which is referred to as an orbit of the system. If the system shows stable behaviour, such that the state of the system is unchanged after each application of the map, then the corresponding orbit is referred to as a period 1 orbit. If the system shows stable behaviour, in which it cycles between two different states, then the corresponding orbit is referred to as a period 2 orbit. Neural systems are known to have unstable periodic orbits and the irregular behaviour of the system derives from the state of the system wandering between periodic orbits (So et al., 1998; Harrison et al., 2004). In this paper we describe the use of a technique for the recovery of the periodic orbits of the oculomotor system during nystagmus, and show that the technique can be used to reveal changes in the nystagmus with gaze angle despite the variability of early-onset nystagmus. 2. Methodology
135
The simplest type of behaviour of a system is that it settles into an equilibrium state. When the system is in equilibrium: dx1 = f1 (x1 , . . . , xn ) = 0 dt dx2 = f2 (x1 , . . . , xn ) = 0 dt .. .
(2)
dxn = fn (x1 , . . . , xn ) = 0 dt Close to an equilibrium point x* , the behaviour of the system is approximately linear so if x is a state in a neighbourhood of the equilibrium, then the behaviour of the system can be described by a set of linear equations: d(x1 − x1∗ ) = a11 (x1 − x1∗ ) + a12 (x2 − x2∗ )+· · · + a1n (xn − xn∗ ) dt d(x2 − x2∗ ) = a21 (x1 − x1∗ ) + a22 (x2 − x2∗ ) + · · · + a2n (xn − xn∗ ) dt .. . d(xn − xn∗ ) = an1 (x1 − x1∗ ) + an2 (x2 − x2∗ ) + · · · + ann (xn − xn∗ ) dt (3)
2.1. Linearisation Our approach begins with the assumption that the oculomotor system can be modelled by a nonlinear dynamical system. Following this approach, the behaviour of the eye movement system can be described by a set of n differential equations:
where ai1 , ai2 , . . ., ain are constants. The quantity aij describes how the jth variable affects the ith variable, so if the jth variable does not affect the ith variable then aij = 0 (Stark et al., 2003). The value of each constant aij corresponds to the partial derivative of fi with respect to xj evaluated at the equilibrium state x* . By applying this equivalence the system Eq. (3) can be re-written as:
∂f1 (x1 , . . . , xn ) ∂f1 (x1 , . . . , xn ) ∂f1 (x1 , . . . , xn ) d(x1 − x1∗ ) = (x1 − x1∗ ) + (x2 − x2∗ ) + · · · + (xn − xn∗ ) dt ∂x1 ∂x2 ∂xn ∂f2 (x1 , . . . , xn ) ∂f2 (x1 , . . . , xn ) ∂f2 (x1 , . . . , xn ) d(x2 − x2∗ ) (x1 − x1∗ ) + (x2 − x2∗ ) + · · · + (xn − xn∗ ) = ∂x1 ∂x2 ∂xn dt .. . d(xn − xn∗ ) ∂fn (x1 , . . . , xn ) ∂fn (x1 , . . . , xn ) ∂fn (x1 , . . . , xn ) (x1 − x1∗ ) + (x2 − x2∗ ) + · · · + (xn − xn∗ ) = dt ∂x1 ∂x2 ∂xn
dx1 = f1 (x1 , x2 , . . . , xn ) dt dx2 = f2 (x1 , x2 , . . . , xn ) dt .. .
(1)
(4)
The array of partial derivatives can be written as a matrix with i rows and j columns which is referred to as the Jacobian. If F is the set of nonlinear functions f1 –fn , then the Jacobian of F may be denoted by F. In summary, the process of linearisation replaces the original nonlinear differential equations which are valid for all states of the system, with linear equations which are only valid for small changes around the equilibrium.
dxn = fn (x1 , x2 , . . . , xn ) dt
2.2. Finding fixed points
where the variables of x1 , . . ., xn constitute the state of the system and the functions f1 , . . ., fn determine the dynamics of the system. Several eye movement models have now been formulated in terms of such systems of differential equations (Gancarz and Grossberg, 1998; Broomhead et al., 2000; Laptev et al., 2006).
In practice it is not usually possible to measure all the state variables x1 , . . ., xn . However, geometric techniques are available which enable estimation of the state of the system. Within the geometric framework, the values of the n state variables are treated as the coordinates of a point in an n-dimensional state
136
M. Theodorou, R.A. Clement / Journal of Neuroscience Methods 161 (2007) 134–141
space, and the successive states adopted by a system over time appear as a curve in this state space. The ‘method of delays’ can be used to reconstruct the curve from the time series (y[1], y[2], . . ., y[n]) obtained from the measurements of a single variable of the system. The method involves taking successive measurements and using these d measurements as the co-ordinates of a point in a d-dimensional delay space. Provided d is large enough, the curve in delay space corresponds to the actual curve in state space. If the nth set of d successive samples are labelled Y[n] = (y[n], y[n + 1], . . ., y[n + d − 1]) then the discrete map from one state to the next can be described by the difference equation: Y [n + 1] = M(Y [n]) where M is a set of nonlinear maps m1 –md . The geometric interpretation of an equilibrium state Y* of the map M is that the coordinates of the point Y* do not change under the map, so it is natural to refer to the equilibrium state as a fixed point. If Y* is a fixed point, the behaviour of the system in the neighbourhood of Y* can be described by a linearisation of the map M: Y [n + 1] = Y ∗ + ∇M(Y [n] − Y ∗ )
Y [n + 1] = (I − ∇M)Y ∗ + ∇MY [n]
Y ∗ = (I − ∇M)−1 (Y [n + 1] − ∇MY [n]) (M(Y [n]) − ∇MY [n])
y[n + 2] − y∗ = a1 (y[n + 1] − y∗ ) + a2 (y[n] − y∗ )
(8)
y[n + 1] − y∗ = a1 (y[n] − y∗ ) + a2 (y[n − 1] − y∗ )
(9)
y[n] − y∗ = a1 (y[n − 1] − y∗ ) + a2 (y[n − 2] − y∗ )
(10)
y[n + 2] − y[n + 1] = a1 (y[n + 1] − y[n])
where I is the d × d identity matrix. This latter equation can be re-arranged to obtain an expression for the fixed point Y* :
= (1 − ∇M)
where a1 and a2 are constants corresponding to the partial derivatives. The Jacobian can be estimated numerically by using finite differences. In a two-dimensional system, there are 2 partial derivatives which can be solved with 2 linear equations: The linear map has the form:
Subtracting Eq. (9) from (8), and (9) from (10) gives two equations which can be solved for the two unknown constants.
which can be re-written as:
−1
Close to the equilibrium state, the Jacobian may be represented by: ⎤ ⎡ ∂m1 (y[n], y[n − 1]) ∂m1 (y[n], y[n − 1]) a 1 a2 ⎦= ⎣ ∂y[n] ∂y[n − 1] 1 0 1 0 (7)
+ a2 (y[n] − y[n − 1])
(11)
y[n] − y[n + 1] = a1 (y[n − 1] − y[n]) (5)
From the definition of a fixed point it follows that the ddimensional delay vector corresponding to a fixed point will consist of d repeats of the measurement value at the fixed point. Any such vector will lie on the long diagonal in delay space. So the estimates of the fixed points can be further refined by dropping any of the transformed data points which lie beyond a threshold distance from the diagonal. The measurement values associated with the fixed points can then be identified by histogramming the frequency of transformed data points against measurement value. 2.3. Example one: two-dimensional systems In the case of a two-dimensional system, at least two successive measurements, {y[n], y[n + 1]}, are required to reconstruct the state changes. With d = 2, the equation for the discrete map from one delay space vector to the next has the form: y[n + 1] m1 (y[n], y[n − 1]) m1 (y[n], y[n − 1]) = = y[n] m2 (y[n], y[n − 1]) y[n] (6)
+ a2 (y[n − 2] − y[n + 1])
(12)
In early-onset nystagmus a variable which is accessible to measurement is horizontal eye position (Abadi et al., 1997). Individual cycles of nystagmus can be identified by thresholding the velocity of the eye movement at a value chosen to be approximately midway through the fast phase of the nystagmus. The thresholding procedure can be used to convert the eye position recording into a time series of the durations of successive cycles, as illustrated in Fig. 1. We carried out this transformation on our data and then applied the fixed point technique for a two-dimensional system. In this case the identified fixed points correspond to the cycle length of the periodic orbits of the system. When histogramming the transformed data points we used bin widths of 25 ms. To find examples of cycles close to the periodic orbit we selected cycles which had a duration which were within +/−12.5 ms of the midpoint of the bin of the peak of the histogram. 2.4. Example two: three-dimensional systems In the case of a three-dimensional system where at least three successive measurements, {y[n], y[n + 1], y[n + 2]}, are required,
M. Theodorou, R.A. Clement / Journal of Neuroscience Methods 161 (2007) 134–141
137
where a1 , a2 and a3 are constants corresponding to the partial derivatives. The linear map has the form: y[n + 3] − y∗ = a1 (y[n + 2] − y∗ ) + a2 (y[n + 1] − y∗ ) + a3 (y[n] − y∗ )
(15)
y[n + 2] − y∗ = a1 (y[n + 1] − y∗ ) + a2 (y[n] − y∗ ) + a3 (y[n − 1] − y∗ )
(16)
y[n + 1] − y∗ = a1 (y[n] − y∗ ) + a2 (y[n − 1] − y∗ ) + a3 (y[n − 2] − y∗ )
(17)
y[n] − y∗ = a1 (y[n − 1] − y∗ ) + a2 (y[n − 2] − y∗ ) + a3 (y[n − 3] − y∗ )
(18)
Subtracting Eq. (17) from each of the Eqs. (15), (16) and (18) gives you three equations: y[n + 3] − y[n + 1] = a1 (y[n + 2] − y[n]) + a2 (y[n + 1] − y[n − 1]) + a3 (y[n] − y[n − 2]) (19) y[n + 2] − y[n + 1] = a1 (y[n + 1] − y[n]) + a2 (y[n] − y[n − 1]) + a3 (y[n − 1] − y[n − 2]) (20) Fig. 1. Application of the fixed point technique to finding periodic orbits within recordings of early-onset nystagmus. (A) 1 s of the eye position trace recorded from subject 1 at a sampling rate of 1000 Hz. (B) The corresponding velocity trace, with the velocity threshold crossings used to identify the interval lengths identified by dots. (C) Time series of interval lengths between successive threshold crossings obtained from 60 s of eye position measurements. (D) Histogram of estimated fixed point intervals obtained with bin widths of 25 ms. The peak of the histogram identifies the interval of the periodic orbit, which lies between 250 and 275 ms in this example.
expanding the nonlinear function S gives: ⎡ ⎤ ⎡ ⎤ y[n + 2] m1 (y[n + 1], y[n], y[n − 1]) ⎢ ⎥ ⎢ ⎥ ⎣ y[n + 1] ⎦ = ⎣ m2 (y[n + 1], y[n], y[n − 1]) ⎦ y[n] m3 (y[n + 1], y[n], y[n − 1]) ⎡ ⎤ m1 (y[n + 1], y[n], y[n − 1]) ⎢ ⎥ y[n + 1] =⎣ ⎦ y[n]
(13)
Close to the equilibrium state the Jacobian may be represented by: ⎡ ∂m (y[n + 1], y[n], y[n − 1]) ∂m (y[n + 1], y[n], y[n − 1]) 1
⎢ ⎢ ⎣
2
∂y[n + 1] 1 0
y[n] 0 1
y[n] − y[n + 1] = a1 (y[n − 1] − y[n]) + a2 (y[n − 2] − y[n − 1]) + a3 (y[n − 3] − y[n − 2]) (21) which can be solved for the three unknown constants. In this paper we use the three-dimensional implementation of the fixed point technique to find the fixed points of eye movement systems with nystagmus. The state space behaviour of the system is reconstructed using delay vectors consisting of three measurements of eye position separated by 20 ms intervals. In addition, we go on to characterise the linear stability of the systems close to the fixed points, by applying techniques which have previously been used to show that the behaviour of an eye movement system producing jerk nystagmus corresponds to a three-dimensional linear system with a fixed point at the delay vector produced by steady fixation of the target direction. The details of these local linear techniques have been fully described
∂m3 (y[n + 1], y[n], y[n − 1]) ⎤ ⎡ a1 ⎥ ⎢ y[n − 1] ⎥ ≡ ⎣1 ⎦ 0 0 0
a2 0 1
⎤ a3 ⎥ 0 ⎦ 0
(14)
138
M. Theodorou, R.A. Clement / Journal of Neuroscience Methods 161 (2007) 134–141
in previous publications (Abadi et al., 1997; Akman et al., 2006) and will only be briefly described here. The dimensionality of the dynamics of the eye movements system close to a fixed point was established by application of a local principal components method. For a fixed point at the origin of the delay space, the method involves the following steps. First, a seven-dimensional delay space reconstruction of the behaviour of the system is formed. Second, the delay vectors within nested spheres centred on the origin are selected. Third, for each set of vectors singular value decomposition is used to calculate orthogonal filters (singular vectors) which are ordered according to the proportion of the variance (singular values) of the eye position recording that they pass. Finally, the singular values which scale linearly with the radius of the spheres are found. In the case of jerk nystagmus, there were three such vectors which implied that the delay vectors spanned a linear space of dimension 3. Subsequently a local linear model was fitted to the movement of the delay vectors in the neighbourhood of the fixed point. The construction involved choosing the radius of a sphere around the origin of the delay space within the range in which the singular values scale linearly with radius, and then calculating the linear map from a selection of the differences between pairs of the delay vectors within the sphere, to the differences between the vectors 20 sample intervals later. Finally, the eigenvalues of the linear maps are calculated. In our application 500 delay vectors were selected and the procedure was repeated 100 times with different selections, in order to obtain an estimate of the variance associated with the eigenvalues.
ity of 6/18 Snellen in each eye and a small left head turn. Subject 2 was a 37-year-old male with best corrected visual acuity of 6/12 Snellen in each eye and a left head turn. 3. Results Figs. 2 and 3 illustrate the results of applying the fixed point technique to the eye movements of two subjects with earlyonset nystagmus. Subject 1, whose eye movements are shown in Figs. 1 and 2, has jerk nystagmus in which the eyes make alternate movements of a slow drift away from the target followed by a fast movement back. Although the waveform has the same form in all gaze directions, the peak-to-peak amplitudes of the cycles are least in the straight ahead position. Analysis of variance showed there was a significant trend with gaze angle (F(d.f.=2.24) = 37.88, P < 0.001). The dynamics of the movement correspond to those found previously for jerk nystagmus (Abadi et al., 1997). Eigenvalue 1 is larger than unity and describes unstable behaviour in that the state of the system moves exponentially away from the fixed point along the direction of the corresponding eigenvector. In jerk nystagmus this eigenvalue characterises the increasing velocity slow phase component of the movement. Eigenvalue 3 is less than unity and describes stable behaviour in which the state of the system converges rapidly towards the fixed point. In jerk nystagmus this characterises the stable fast phase movement back to the target. The second eigen-
2.5. Eye movement recording The stimulus for eye movements consisted of a row of 3 red LEDs which were spaced 15 degrees apart (0◦ , −15◦ (left) and +15◦ (right)). Each LED subtended 0.35◦ at the eye. The subject was asked to look at the illuminated LED with both eyes open. The stimulus sequence consisted of cycles in which first the LED in the primary position was illuminated (0◦ ), followed by the LED at −15◦ , +15◦ and 0◦ again. Each LED was illuminated for 1 min. The head was stabilised in the primary position with a chin rest. The eye movements were recorded using an IRIS 6500 infrared limbal tracker (Skalar Medical, Delft, The Netherlands). The analogue output was filtered by a 100 Hz low-pass filter, digitised to 12-bit resolution, and then sampled at 1 ms intervals. The system was linear over a range of ±25◦ (horizontal), with a resolution of 0.03◦ . Calibration was carried out using the 15◦ saccades (i.e. right and left), manually adjusting the graphs so that the foveating part of the waveform was set at approximately 0◦ and ±15◦ , respectively. 2.6. Subjects Eye movement recordings were made from two adult subjects with early-onset nystagmus. Both subjects had bilateral conjugate horizontal eye movements. There were no other ocular or neurological abnormalities detectable on clinical examination. Subject 1 was a 58-year-old male with best corrected visual acu-
Fig. 2. Examples of periodic orbits found in the nystagmus of subject 1. In this and the subsequent figure, examples of the eye movements that occur when the state of the oculomotor system is following a periodic orbit are shown on the left and a reconstruction of the periodic orbit is shown on the right. The examples for which the beginning and end of the orbit are closest are shown in black and the reconstruction of these orbits that are shown on the right. The examples are shown for three different gaze positions, 15◦ left, straight ahead and 15◦ right.
M. Theodorou, R.A. Clement / Journal of Neuroscience Methods 161 (2007) 134–141
139
Table 1 Amplitude and the real parts of the eigenvalues (mean (standard deviation)) of the waveform of subject 1 in the 3 horizontal positions of gaze Target position
Amplitude (◦ )
Eigenvalue 1
Eigenvalue 2
Eigenvalue 3
Right Straight Left
9.13 (1.87) 3.56 (0.65) 6.00 (0.82)
1.591 (0.095) 1.446 (0.037) 1.331 (0.027)
0.922 (0.085) 0.989 (0.029) 0.885 (0.031)
0.242 (0.060) 0.116 (0.019) 0.035 (0.020)
which are approached at the extremes of the oscillations, and their associated eigenvalues are given in Table 2. In all cases where there was a saccade toward the fixed point, an unstable and a stable eigenvalue were found, which appear to correspond to the unstable and stable eigenvalues found in jerk nystagmus. However, instead of a neutral eigenvalue, a pair of complex conjugate eigenvalues were found, which describe a decreasing amplitude oscillation about the fixed point. The reason the eigenvalues for the left fixed point do not appear to have a complex conjugate pair of eigenvalues is because of the computational procedure used to calculate the eigenvalues. In order to estimate the distribution of the eigenvalues, a bootstrap technique was used in which a number of different selections of delay vectors and the eigenvalues were computed for each set. For each set a pair of complex conjugate values was always found, but in the case of the fixed point the real parts of the pair were sometimes the largest or the second largest eigenvalues. Averaging over many examples of these two cases obliterates the evidence of the complex conjugate pairs. Physiologically we interpret this change in the relative amplitude of the eigenvalues as signifying a shift in the mechanism underlying the nystagmus away from a jerk mechanism. Fig. 3. Examples of periodic orbits found in the nystagmus of subject 2.
4. Discussion value is approximately unity, which is neutral with respect to stability. The direction of the eigenvector associated with this eigenvalue corresponds to the line of positions in phase space in which the end points of the fast phase of nystagmus lie. The three eigenvalues of the single fixed point in subject 1 are given in Table 1. The eye movements of subject 2 are shown in Fig. 3. This subject has a combination of pendular and jerk waveforms. In this case the form of the waveform changes, becoming bidirectional on right gaze. The peak-to-peak amplitudes of the oscillations are also least on right gaze. Analysis of variance showed that there was a significant trend with gaze angle (F(d.f.=2.60) = 155.41, P < 0.001). Two fixed points were found,
The fixed point technique has previously been applied to the analysis of epileptiform discharges in hippocampal slices (So et al., 1997; Slutzky et al., 2001; Harrison et al., 2004), and so far has not had as wider application to neuroscience as might be expected, with the only other application being to nystagmus (Clement et al., 2002). One possible reason for the limited take up is that the original description (So et al., 1997) contained some necessary mathematical justifications which may have obscured the simplicity of the technique. In adapting the technique for motor disorders we have been able to simplify some of the technical procedures given in the original description. For instance, the addition of a tensor of random numbers to the equation for estimating the fixed point is necessary to be
Table 2 Amplitude and the real parts of the eigenvalues (mean (standard deviation)) of the waveform of subject 2 in the 3 horizontal positions of gaze Target position
Amplitude (◦ )
Eigenvalue 1
Eigenvalue 2
Eigenvalue 3
Eigenvalue 4
Right
3.10 (0.31)
1.795 (0.084) 1.701 (0.076)
0.247 (0.017) 0.370 (0.033)
0.090 (0.041) 0.370 (0.033)
0.090 (0.041) 0.162 (0.079)
Straight
6.51 (0.74)
1.627 (0.014) 1.659 (0.060)
0.668 (0.003) 0.128 (0.049)
0.131 (0.001) 0.126 (0.048)
0.131 (0.001) 0.097 (0.062)
Left
7.21 (0.97)
1.212 (0.182) 1.243 (0.041)
1.104 (0.094) 0.611 (0.021)
0.946 (0.149) 0.611 (0.021)
0.572 (0.308) 0.043 (0.015)
140
M. Theodorou, R.A. Clement / Journal of Neuroscience Methods 161 (2007) 134–141
sure of avoiding spurious fixed points (So et al., 1997) and was used in the previous application to nystagmus, (Clement et al., 2002). However, the use of the tensor to exclude spurious fixed points is not needed with eye movement data because each cycle of such data has to be examined to exclude any artefacts due to attentional changes in fixation or blinks, and anomalous cycles would be picked up at this stage of the analysis. Another possible reason is that the technique has been described as a language for neuronal dynamics in which the behaviour of the system is described by a hierarchy of periodic orbits (So et al., 1998). But finding higher order periodic orbits becomes increasingly difficult, because there are relatively fewer examples, and the evaluation of the eigenvalues becomes sensitive to the examples used. The difficulty of obtaining reliable estimates of the eigenvalues simply, does not support the use of the fixed point technique as a method for building a quantitative model of a neural system. In the context of nystagmus, the major advantage of the technique is that it enables quantitative descriptions of the effects of experimental manipulations to be obtained despite the variable nature of the waveform (Reccia et al., 1990; Clement et al., 2002; Optican and Miura, 2004). After using the technique to recover examples of the underlying periodic orbits in two cases of nystagmus, we have been able to show that there are significant changes in the peak-to-peak amplitudes of the example cycles with changes in gaze angle. This finding suggests that the technique could be a useful additional technique for analysing other oscillatory motor phenomena, such as the tremor that arises in Parkinson’s disease. The most immediate application of the technique is likely to be in providing reliable outcome measures when treatments are being tested. Furthermore, a recent study of drug therapy for nystagmus (Shery et al., 2006) showed that in some cases the effect of treatment was greater on pendular than on jerk waveforms. By enabling quantitative estimates of the relative instabilities of the jerk and pendular components of the waveform, it may be possible to use the technique to predict which patients will benefit most from the treatment. Application of the fixed point technique to two different forms of early-onset nystagmus has revealed previously unsuspected differences in the dynamics of their eye movements. In the jerk form of nystagmus (subject 1), the orbits of the system repeatedly approach a single unstable fixed point. In the pendular form of nystagmus (subject 2), the orbits of the system cycle between two unstable fixed points. A previous analysis of a model for the mechanism of nystagmus has shown that it is an example of a symmetrical system with a pair of homoclinic orbits, that is orbits which begin and end at the same unstable fixed point (Akman et al., 2005). With appropriate parameter changes, the pair of homoclinic orbits can join at the fixed point to form a single cycle by a process referred to as a gluing bifurcation. In the case of jerk nystagmus, the dynamics correspond to having left and right beating jerk waveforms which join to make a bi-directional jerk waveform. But as shown by the results of subject 2 the fixed point technique also reveals a type of bifurcation which was not expected. Subject 2 has a pair of heteroclinic orbits, that is orbits which join two different fixed points during a complete cycle. On left gaze the orbit contains a left beat-
ing jerk component and on right gaze the orbit contains a right beating jerk component. As is common in pendular nystagmus with foveating saccades, as the gaze angle is shifted from left to right, a bidirectional waveform appears with both left and right beating jerk components. In this case a simultaneous gluing at both fixed points appears to occur. Although the possibility of such a bifurcation remains to be confirmed by numerical simulation and analysis, the provisional finding suggests that the fixed point approach may be a powerful technique for revealing novel bifurcations in biological systems which have not previously been discovered in physics or engineering applications. Acknowledgements We thank Ozgur Akman for helpful discussions. This work was supported by grants from The Ormsby Charitable Trust and The British Eye Research Foundation. References Abadi RV, Broomhead D, Clement RA, Whittle JP, Worfolk R. Dynamical systems analysis: a new method of analysing congenital nystagmus waveforms. Exp Brain Res 1997;117:355–61. Abadi RV, Worfolk R. Harmonic analysis of congenital nystagmus waveforms. Clin Vision Sci 1991;6:385–8. Akman OE, Broomhead DS, Abadi RV, Clement RA. Eye movement instabilities and nystagmus can be predicted by a nonlinear dynamics model of the saccadic system. J Math Biol 2005;51:661–94. Akman OE, Broomhead DS, Clement RA, Abadi RV. Nonlinear time series analysis of jerk congenital nystagmus. J Comput Neurosci 2006;21(2):153–70. Broomhead DS, Clement RA, Muldoon MR, Whittle JP, Scallan C, Abadi RV. Modelling of congenital nystagmus waveforms produced by saccadic system abnormalities. Biol Cybern 2000;82:391–9. Clement RA, Whittle JP, Muldoon MR, Abadi RV, Broomhead D, Akman O. Characterisation of congenital nystagmus waveforms in terms of periodic orbits. Vision Res 2002;42:2123–30. Dell’Osso LF. Tenotomy and congenital nystagmus: a failure to answer the wrong question. Vision Res 2004;44:3091–4. Dell’Osso LF, Daroff RB. Congenital nystagmus waveforms and foveation strategy. Doc Ophthalmol 1975;39:155–82. Dell’Osso LF, Flynn JT, Daroff RB. Hereditary congenital nystagmus. An intrafamilial study. Arch Ophthalmol 1974;92:366–74. Dickinson CM, Abadi RV. The influence of nystagmoid oscillation on contrast sensitivity in normal observers. Vision Res 1985;25:1089–96. Gancarz G, Grossberg G. A neural model of the saccade generator in the reticular formation. Neural Netw 1998;11:1159–74. Harris C. Nystagmus and eye movement disorders. In: Taylor D, editor. Paediatric ophthalmology. 2nd ed. Blackwell; 1997. p. 869–96. Harrison PK, Tattersall JE, Clement RA. Periodic orbit analysis reveals subtle effects of atropine on epileptiform activity in the guinea-pig hippocampal slice. Neurosci Lett 2004;357:183–6. Laptev D, Akman OE, Clement RA. Stability of the saccadic oculomotor system. Biol Cybern 2006;95:281–7. Miura K, Hertle RW, FitzGibbon EJ, Optican LM. Effects of tenotomy surgery on congenital nystagmus waveforms in adult patients. Part I. Wavelet spectral analysis. Vision Res 2003;43:2345–56. Optican LM, Miura K. Tenotomy and congenital nystagmus: a null result is not a failure, for “It is not the answer that enlightens, but the question”. Vision Res 2004;44:3095–8. Reccia R, Roberti G, Russo P. Spectral analysis of pendular waveforms in congenital nystagmus. Ophthalmic Res 1989;21:83–92. Reccia R, Roberti G, Russo P. Computer analysis of ENG spectral features from patients with congenital nystagmus. J Biomed Eng 1990;12:39–45.
M. Theodorou, R.A. Clement / Journal of Neuroscience Methods 161 (2007) 134–141 Shery T, Proudlock FA, Sarvananthan N, McLean RJ, Gottlob I. The effects of gabapentin and memantine in acquired nystagmus: a retrospective study. BJO 2006;124:839–43. Slutzky MW, Cvitanovic P, Moul DJ. Deterministic Chaos and noise in three in vitro hippocampal models of epilepsy. Ann Biomed Eng 2001;29:607–18. So P, Francis JT, Netoff TI, Gluckman BJ, Schiff SJ. Periodic orbits: a new language for neuronal dynamics. Biophys J 1998;74:2776–85.
141
So P, Ott E, Sauer SJ, Gluckman BJ, Grebogi C, Schiff SJ. Extracting unstable periodic orbits from chaotic time series data. Phys Rev E 1997;55:5398–417. Stark J, Brewer D, Barenco M, Tomescu D, Callard R, Hubank M. Reconstructing gene networks: what are the limits? Biochem Soc Trans 2003;31: 1519–25. Worfolk R, Abadi RV. Quick phase programming and saccadic re-orientation in congenital nystagmus. Vision Res 1991;31(10):1819–30.
Fixed point analysis of nystagmus Maria Theodorou ∗ , Richard A. Clement Institute of Child Health, University College London, 30 Guilford Street, London WC1N 1EH, UK Received 16 February 2006; received in revised form 3 September 2006; accepted 29 September 2006
Abstract Motor disorders frequently contain a rhythmic component, but the associated oscillations are not usually precisely periodic. This lack of strict periodicity can make it difficult to identify the effects of experimental manipulations on the oscillation. In this report, we describe the application of a numerical technique for identifying fixed points of a nonlinear map to the recovery of underlying periodicities of the eye movement disorder of nystagmus. The technique is illustrated by application to two different types of nystagmus. In addition we use a local analysis of the behaviour at the fixed points to distinguish between different bifurcations in the two examples with changes in gaze angle. We conclude that the technique reveals consistent effects of experimental manipulations, which may be useful for quantitative characterisation of experimental and therapeutic manipulations of motor disorders. © 2006 Elsevier B.V. All rights reserved. Keywords: Nystagmus; Oscillations; Deterministic; Nonlinear dynamics; Fourier analysis; Wavelet analysis; Phase space analysis
1. Introduction Motor disorders frequently lead to the development of oscillations. In order to evaluate any treatment of the disorders quantitative techniques for describing the oscillations are required. An example of such an oscillation is the eye movement disorder of early-onset nystagmus. Early-onset nystagmus is an involuntary oscillation of the eyes which develops at birth or shortly afterwards, often within the first 6 months and persists throughout life. It is characterised by predominantly bilateral, conjugate, horizontal eye movements. It may be idiopathic, associated with a sensory abnormality, or more rarely, a neurological lesion (Harris, 1997). The main difficulty in analysing earlyonset nystagmus waveforms is that they are very variable, so that it has been difficult to obtain consistent results (Reccia et al., 1990; Clement et al., 2002; Optican and Miura, 2004). This is a common problem amongst motor disorders, so that effective techniques for analysing eye movement recordings should have widespread applicability. Eye movement recordings from subjects with early-onset nystagmus show a variety of different waveforms which can be described by combinations of jerk and pendular waveforms
∗
Corresponding author. Tel.: +44 207 905 2346. E-mail address: [email protected] (M. Theodorou).
0165-0270/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jneumeth.2006.09.025
(Dell’Osso and Daroff, 1975). A jerk waveform is comprised of an increasing velocity exponential drift of the eye away from the target direction followed by a quick movement which returns the eye to the direction of the target. A typical pendular waveform consists of an approximately sinusoidal oscillation. The oscillatory movement of early-onset nystagmus can be characterized by the Fourier components of the eye position signal (Dell’Osso et al., 1974; Dickinson and Abadi, 1985; Reccia et al., 1989, 1990; Abadi and Worfolk, 1991; Worfolk and Abadi, 1991; Clement et al., 2002). Typically there is a peak in the power spectrum at the fundamental frequency of the oscillation, together with harmonic components with smaller amplitudes. The main peak provides a good measure of the average frequency and amplitude of the early-onset nystagmus waveform. Fourier analysis has revealed that even pure pendular nystagmus, which is often described as sinusoidal, has harmonic components (Reccia et al., 1990; Abadi and Worfolk, 1991). A recent alternative to Fourier analysis is wavelet analysis (Miura et al., 2003), which allows localization in both the time and frequency domain. However, a much longer segment of data is required to get consistent results which can lead to changes in the waveform being obscured by attention changes in the subject (Dell’Osso, 2004). A more promising approach has arisen through the application of dynamical systems theory. The deterministic components of the nystagmus oscillations can be identified using a tech-
M. Theodorou, R.A. Clement / Journal of Neuroscience Methods 161 (2007) 134–141
nique developed in nonlinear dynamics. In this framework the behaviour of the neural system can be described by a map from one state to the next. The state of the neural system can be represented geometrically by a point; successive states form a sequence of points, which is referred to as an orbit of the system. If the system shows stable behaviour, such that the state of the system is unchanged after each application of the map, then the corresponding orbit is referred to as a period 1 orbit. If the system shows stable behaviour, in which it cycles between two different states, then the corresponding orbit is referred to as a period 2 orbit. Neural systems are known to have unstable periodic orbits and the irregular behaviour of the system derives from the state of the system wandering between periodic orbits (So et al., 1998; Harrison et al., 2004). In this paper we describe the use of a technique for the recovery of the periodic orbits of the oculomotor system during nystagmus, and show that the technique can be used to reveal changes in the nystagmus with gaze angle despite the variability of early-onset nystagmus. 2. Methodology
135
The simplest type of behaviour of a system is that it settles into an equilibrium state. When the system is in equilibrium: dx1 = f1 (x1 , . . . , xn ) = 0 dt dx2 = f2 (x1 , . . . , xn ) = 0 dt .. .
(2)
dxn = fn (x1 , . . . , xn ) = 0 dt Close to an equilibrium point x* , the behaviour of the system is approximately linear so if x is a state in a neighbourhood of the equilibrium, then the behaviour of the system can be described by a set of linear equations: d(x1 − x1∗ ) = a11 (x1 − x1∗ ) + a12 (x2 − x2∗ )+· · · + a1n (xn − xn∗ ) dt d(x2 − x2∗ ) = a21 (x1 − x1∗ ) + a22 (x2 − x2∗ ) + · · · + a2n (xn − xn∗ ) dt .. . d(xn − xn∗ ) = an1 (x1 − x1∗ ) + an2 (x2 − x2∗ ) + · · · + ann (xn − xn∗ ) dt (3)
2.1. Linearisation Our approach begins with the assumption that the oculomotor system can be modelled by a nonlinear dynamical system. Following this approach, the behaviour of the eye movement system can be described by a set of n differential equations:
where ai1 , ai2 , . . ., ain are constants. The quantity aij describes how the jth variable affects the ith variable, so if the jth variable does not affect the ith variable then aij = 0 (Stark et al., 2003). The value of each constant aij corresponds to the partial derivative of fi with respect to xj evaluated at the equilibrium state x* . By applying this equivalence the system Eq. (3) can be re-written as:
∂f1 (x1 , . . . , xn ) ∂f1 (x1 , . . . , xn ) ∂f1 (x1 , . . . , xn ) d(x1 − x1∗ ) = (x1 − x1∗ ) + (x2 − x2∗ ) + · · · + (xn − xn∗ ) dt ∂x1 ∂x2 ∂xn ∂f2 (x1 , . . . , xn ) ∂f2 (x1 , . . . , xn ) ∂f2 (x1 , . . . , xn ) d(x2 − x2∗ ) (x1 − x1∗ ) + (x2 − x2∗ ) + · · · + (xn − xn∗ ) = ∂x1 ∂x2 ∂xn dt .. . d(xn − xn∗ ) ∂fn (x1 , . . . , xn ) ∂fn (x1 , . . . , xn ) ∂fn (x1 , . . . , xn ) (x1 − x1∗ ) + (x2 − x2∗ ) + · · · + (xn − xn∗ ) = dt ∂x1 ∂x2 ∂xn
dx1 = f1 (x1 , x2 , . . . , xn ) dt dx2 = f2 (x1 , x2 , . . . , xn ) dt .. .
(1)
(4)
The array of partial derivatives can be written as a matrix with i rows and j columns which is referred to as the Jacobian. If F is the set of nonlinear functions f1 –fn , then the Jacobian of F may be denoted by F. In summary, the process of linearisation replaces the original nonlinear differential equations which are valid for all states of the system, with linear equations which are only valid for small changes around the equilibrium.
dxn = fn (x1 , x2 , . . . , xn ) dt
2.2. Finding fixed points
where the variables of x1 , . . ., xn constitute the state of the system and the functions f1 , . . ., fn determine the dynamics of the system. Several eye movement models have now been formulated in terms of such systems of differential equations (Gancarz and Grossberg, 1998; Broomhead et al., 2000; Laptev et al., 2006).
In practice it is not usually possible to measure all the state variables x1 , . . ., xn . However, geometric techniques are available which enable estimation of the state of the system. Within the geometric framework, the values of the n state variables are treated as the coordinates of a point in an n-dimensional state
136
M. Theodorou, R.A. Clement / Journal of Neuroscience Methods 161 (2007) 134–141
space, and the successive states adopted by a system over time appear as a curve in this state space. The ‘method of delays’ can be used to reconstruct the curve from the time series (y[1], y[2], . . ., y[n]) obtained from the measurements of a single variable of the system. The method involves taking successive measurements and using these d measurements as the co-ordinates of a point in a d-dimensional delay space. Provided d is large enough, the curve in delay space corresponds to the actual curve in state space. If the nth set of d successive samples are labelled Y[n] = (y[n], y[n + 1], . . ., y[n + d − 1]) then the discrete map from one state to the next can be described by the difference equation: Y [n + 1] = M(Y [n]) where M is a set of nonlinear maps m1 –md . The geometric interpretation of an equilibrium state Y* of the map M is that the coordinates of the point Y* do not change under the map, so it is natural to refer to the equilibrium state as a fixed point. If Y* is a fixed point, the behaviour of the system in the neighbourhood of Y* can be described by a linearisation of the map M: Y [n + 1] = Y ∗ + ∇M(Y [n] − Y ∗ )
Y [n + 1] = (I − ∇M)Y ∗ + ∇MY [n]
Y ∗ = (I − ∇M)−1 (Y [n + 1] − ∇MY [n]) (M(Y [n]) − ∇MY [n])
y[n + 2] − y∗ = a1 (y[n + 1] − y∗ ) + a2 (y[n] − y∗ )
(8)
y[n + 1] − y∗ = a1 (y[n] − y∗ ) + a2 (y[n − 1] − y∗ )
(9)
y[n] − y∗ = a1 (y[n − 1] − y∗ ) + a2 (y[n − 2] − y∗ )
(10)
y[n + 2] − y[n + 1] = a1 (y[n + 1] − y[n])
where I is the d × d identity matrix. This latter equation can be re-arranged to obtain an expression for the fixed point Y* :
= (1 − ∇M)
where a1 and a2 are constants corresponding to the partial derivatives. The Jacobian can be estimated numerically by using finite differences. In a two-dimensional system, there are 2 partial derivatives which can be solved with 2 linear equations: The linear map has the form:
Subtracting Eq. (9) from (8), and (9) from (10) gives two equations which can be solved for the two unknown constants.
which can be re-written as:
−1
Close to the equilibrium state, the Jacobian may be represented by: ⎤ ⎡ ∂m1 (y[n], y[n − 1]) ∂m1 (y[n], y[n − 1]) a 1 a2 ⎦= ⎣ ∂y[n] ∂y[n − 1] 1 0 1 0 (7)
+ a2 (y[n] − y[n − 1])
(11)
y[n] − y[n + 1] = a1 (y[n − 1] − y[n]) (5)
From the definition of a fixed point it follows that the ddimensional delay vector corresponding to a fixed point will consist of d repeats of the measurement value at the fixed point. Any such vector will lie on the long diagonal in delay space. So the estimates of the fixed points can be further refined by dropping any of the transformed data points which lie beyond a threshold distance from the diagonal. The measurement values associated with the fixed points can then be identified by histogramming the frequency of transformed data points against measurement value. 2.3. Example one: two-dimensional systems In the case of a two-dimensional system, at least two successive measurements, {y[n], y[n + 1]}, are required to reconstruct the state changes. With d = 2, the equation for the discrete map from one delay space vector to the next has the form: y[n + 1] m1 (y[n], y[n − 1]) m1 (y[n], y[n − 1]) = = y[n] m2 (y[n], y[n − 1]) y[n] (6)
+ a2 (y[n − 2] − y[n + 1])
(12)
In early-onset nystagmus a variable which is accessible to measurement is horizontal eye position (Abadi et al., 1997). Individual cycles of nystagmus can be identified by thresholding the velocity of the eye movement at a value chosen to be approximately midway through the fast phase of the nystagmus. The thresholding procedure can be used to convert the eye position recording into a time series of the durations of successive cycles, as illustrated in Fig. 1. We carried out this transformation on our data and then applied the fixed point technique for a two-dimensional system. In this case the identified fixed points correspond to the cycle length of the periodic orbits of the system. When histogramming the transformed data points we used bin widths of 25 ms. To find examples of cycles close to the periodic orbit we selected cycles which had a duration which were within +/−12.5 ms of the midpoint of the bin of the peak of the histogram. 2.4. Example two: three-dimensional systems In the case of a three-dimensional system where at least three successive measurements, {y[n], y[n + 1], y[n + 2]}, are required,
M. Theodorou, R.A. Clement / Journal of Neuroscience Methods 161 (2007) 134–141
137
where a1 , a2 and a3 are constants corresponding to the partial derivatives. The linear map has the form: y[n + 3] − y∗ = a1 (y[n + 2] − y∗ ) + a2 (y[n + 1] − y∗ ) + a3 (y[n] − y∗ )
(15)
y[n + 2] − y∗ = a1 (y[n + 1] − y∗ ) + a2 (y[n] − y∗ ) + a3 (y[n − 1] − y∗ )
(16)
y[n + 1] − y∗ = a1 (y[n] − y∗ ) + a2 (y[n − 1] − y∗ ) + a3 (y[n − 2] − y∗ )
(17)
y[n] − y∗ = a1 (y[n − 1] − y∗ ) + a2 (y[n − 2] − y∗ ) + a3 (y[n − 3] − y∗ )
(18)
Subtracting Eq. (17) from each of the Eqs. (15), (16) and (18) gives you three equations: y[n + 3] − y[n + 1] = a1 (y[n + 2] − y[n]) + a2 (y[n + 1] − y[n − 1]) + a3 (y[n] − y[n − 2]) (19) y[n + 2] − y[n + 1] = a1 (y[n + 1] − y[n]) + a2 (y[n] − y[n − 1]) + a3 (y[n − 1] − y[n − 2]) (20) Fig. 1. Application of the fixed point technique to finding periodic orbits within recordings of early-onset nystagmus. (A) 1 s of the eye position trace recorded from subject 1 at a sampling rate of 1000 Hz. (B) The corresponding velocity trace, with the velocity threshold crossings used to identify the interval lengths identified by dots. (C) Time series of interval lengths between successive threshold crossings obtained from 60 s of eye position measurements. (D) Histogram of estimated fixed point intervals obtained with bin widths of 25 ms. The peak of the histogram identifies the interval of the periodic orbit, which lies between 250 and 275 ms in this example.
expanding the nonlinear function S gives: ⎡ ⎤ ⎡ ⎤ y[n + 2] m1 (y[n + 1], y[n], y[n − 1]) ⎢ ⎥ ⎢ ⎥ ⎣ y[n + 1] ⎦ = ⎣ m2 (y[n + 1], y[n], y[n − 1]) ⎦ y[n] m3 (y[n + 1], y[n], y[n − 1]) ⎡ ⎤ m1 (y[n + 1], y[n], y[n − 1]) ⎢ ⎥ y[n + 1] =⎣ ⎦ y[n]
(13)
Close to the equilibrium state the Jacobian may be represented by: ⎡ ∂m (y[n + 1], y[n], y[n − 1]) ∂m (y[n + 1], y[n], y[n − 1]) 1
⎢ ⎢ ⎣
2
∂y[n + 1] 1 0
y[n] 0 1
y[n] − y[n + 1] = a1 (y[n − 1] − y[n]) + a2 (y[n − 2] − y[n − 1]) + a3 (y[n − 3] − y[n − 2]) (21) which can be solved for the three unknown constants. In this paper we use the three-dimensional implementation of the fixed point technique to find the fixed points of eye movement systems with nystagmus. The state space behaviour of the system is reconstructed using delay vectors consisting of three measurements of eye position separated by 20 ms intervals. In addition, we go on to characterise the linear stability of the systems close to the fixed points, by applying techniques which have previously been used to show that the behaviour of an eye movement system producing jerk nystagmus corresponds to a three-dimensional linear system with a fixed point at the delay vector produced by steady fixation of the target direction. The details of these local linear techniques have been fully described
∂m3 (y[n + 1], y[n], y[n − 1]) ⎤ ⎡ a1 ⎥ ⎢ y[n − 1] ⎥ ≡ ⎣1 ⎦ 0 0 0
a2 0 1
⎤ a3 ⎥ 0 ⎦ 0
(14)
138
M. Theodorou, R.A. Clement / Journal of Neuroscience Methods 161 (2007) 134–141
in previous publications (Abadi et al., 1997; Akman et al., 2006) and will only be briefly described here. The dimensionality of the dynamics of the eye movements system close to a fixed point was established by application of a local principal components method. For a fixed point at the origin of the delay space, the method involves the following steps. First, a seven-dimensional delay space reconstruction of the behaviour of the system is formed. Second, the delay vectors within nested spheres centred on the origin are selected. Third, for each set of vectors singular value decomposition is used to calculate orthogonal filters (singular vectors) which are ordered according to the proportion of the variance (singular values) of the eye position recording that they pass. Finally, the singular values which scale linearly with the radius of the spheres are found. In the case of jerk nystagmus, there were three such vectors which implied that the delay vectors spanned a linear space of dimension 3. Subsequently a local linear model was fitted to the movement of the delay vectors in the neighbourhood of the fixed point. The construction involved choosing the radius of a sphere around the origin of the delay space within the range in which the singular values scale linearly with radius, and then calculating the linear map from a selection of the differences between pairs of the delay vectors within the sphere, to the differences between the vectors 20 sample intervals later. Finally, the eigenvalues of the linear maps are calculated. In our application 500 delay vectors were selected and the procedure was repeated 100 times with different selections, in order to obtain an estimate of the variance associated with the eigenvalues.
ity of 6/18 Snellen in each eye and a small left head turn. Subject 2 was a 37-year-old male with best corrected visual acuity of 6/12 Snellen in each eye and a left head turn. 3. Results Figs. 2 and 3 illustrate the results of applying the fixed point technique to the eye movements of two subjects with earlyonset nystagmus. Subject 1, whose eye movements are shown in Figs. 1 and 2, has jerk nystagmus in which the eyes make alternate movements of a slow drift away from the target followed by a fast movement back. Although the waveform has the same form in all gaze directions, the peak-to-peak amplitudes of the cycles are least in the straight ahead position. Analysis of variance showed there was a significant trend with gaze angle (F(d.f.=2.24) = 37.88, P < 0.001). The dynamics of the movement correspond to those found previously for jerk nystagmus (Abadi et al., 1997). Eigenvalue 1 is larger than unity and describes unstable behaviour in that the state of the system moves exponentially away from the fixed point along the direction of the corresponding eigenvector. In jerk nystagmus this eigenvalue characterises the increasing velocity slow phase component of the movement. Eigenvalue 3 is less than unity and describes stable behaviour in which the state of the system converges rapidly towards the fixed point. In jerk nystagmus this characterises the stable fast phase movement back to the target. The second eigen-
2.5. Eye movement recording The stimulus for eye movements consisted of a row of 3 red LEDs which were spaced 15 degrees apart (0◦ , −15◦ (left) and +15◦ (right)). Each LED subtended 0.35◦ at the eye. The subject was asked to look at the illuminated LED with both eyes open. The stimulus sequence consisted of cycles in which first the LED in the primary position was illuminated (0◦ ), followed by the LED at −15◦ , +15◦ and 0◦ again. Each LED was illuminated for 1 min. The head was stabilised in the primary position with a chin rest. The eye movements were recorded using an IRIS 6500 infrared limbal tracker (Skalar Medical, Delft, The Netherlands). The analogue output was filtered by a 100 Hz low-pass filter, digitised to 12-bit resolution, and then sampled at 1 ms intervals. The system was linear over a range of ±25◦ (horizontal), with a resolution of 0.03◦ . Calibration was carried out using the 15◦ saccades (i.e. right and left), manually adjusting the graphs so that the foveating part of the waveform was set at approximately 0◦ and ±15◦ , respectively. 2.6. Subjects Eye movement recordings were made from two adult subjects with early-onset nystagmus. Both subjects had bilateral conjugate horizontal eye movements. There were no other ocular or neurological abnormalities detectable on clinical examination. Subject 1 was a 58-year-old male with best corrected visual acu-
Fig. 2. Examples of periodic orbits found in the nystagmus of subject 1. In this and the subsequent figure, examples of the eye movements that occur when the state of the oculomotor system is following a periodic orbit are shown on the left and a reconstruction of the periodic orbit is shown on the right. The examples for which the beginning and end of the orbit are closest are shown in black and the reconstruction of these orbits that are shown on the right. The examples are shown for three different gaze positions, 15◦ left, straight ahead and 15◦ right.
M. Theodorou, R.A. Clement / Journal of Neuroscience Methods 161 (2007) 134–141
139
Table 1 Amplitude and the real parts of the eigenvalues (mean (standard deviation)) of the waveform of subject 1 in the 3 horizontal positions of gaze Target position
Amplitude (◦ )
Eigenvalue 1
Eigenvalue 2
Eigenvalue 3
Right Straight Left
9.13 (1.87) 3.56 (0.65) 6.00 (0.82)
1.591 (0.095) 1.446 (0.037) 1.331 (0.027)
0.922 (0.085) 0.989 (0.029) 0.885 (0.031)
0.242 (0.060) 0.116 (0.019) 0.035 (0.020)
which are approached at the extremes of the oscillations, and their associated eigenvalues are given in Table 2. In all cases where there was a saccade toward the fixed point, an unstable and a stable eigenvalue were found, which appear to correspond to the unstable and stable eigenvalues found in jerk nystagmus. However, instead of a neutral eigenvalue, a pair of complex conjugate eigenvalues were found, which describe a decreasing amplitude oscillation about the fixed point. The reason the eigenvalues for the left fixed point do not appear to have a complex conjugate pair of eigenvalues is because of the computational procedure used to calculate the eigenvalues. In order to estimate the distribution of the eigenvalues, a bootstrap technique was used in which a number of different selections of delay vectors and the eigenvalues were computed for each set. For each set a pair of complex conjugate values was always found, but in the case of the fixed point the real parts of the pair were sometimes the largest or the second largest eigenvalues. Averaging over many examples of these two cases obliterates the evidence of the complex conjugate pairs. Physiologically we interpret this change in the relative amplitude of the eigenvalues as signifying a shift in the mechanism underlying the nystagmus away from a jerk mechanism. Fig. 3. Examples of periodic orbits found in the nystagmus of subject 2.
4. Discussion value is approximately unity, which is neutral with respect to stability. The direction of the eigenvector associated with this eigenvalue corresponds to the line of positions in phase space in which the end points of the fast phase of nystagmus lie. The three eigenvalues of the single fixed point in subject 1 are given in Table 1. The eye movements of subject 2 are shown in Fig. 3. This subject has a combination of pendular and jerk waveforms. In this case the form of the waveform changes, becoming bidirectional on right gaze. The peak-to-peak amplitudes of the oscillations are also least on right gaze. Analysis of variance showed that there was a significant trend with gaze angle (F(d.f.=2.60) = 155.41, P < 0.001). Two fixed points were found,
The fixed point technique has previously been applied to the analysis of epileptiform discharges in hippocampal slices (So et al., 1997; Slutzky et al., 2001; Harrison et al., 2004), and so far has not had as wider application to neuroscience as might be expected, with the only other application being to nystagmus (Clement et al., 2002). One possible reason for the limited take up is that the original description (So et al., 1997) contained some necessary mathematical justifications which may have obscured the simplicity of the technique. In adapting the technique for motor disorders we have been able to simplify some of the technical procedures given in the original description. For instance, the addition of a tensor of random numbers to the equation for estimating the fixed point is necessary to be
Table 2 Amplitude and the real parts of the eigenvalues (mean (standard deviation)) of the waveform of subject 2 in the 3 horizontal positions of gaze Target position
Amplitude (◦ )
Eigenvalue 1
Eigenvalue 2
Eigenvalue 3
Eigenvalue 4
Right
3.10 (0.31)
1.795 (0.084) 1.701 (0.076)
0.247 (0.017) 0.370 (0.033)
0.090 (0.041) 0.370 (0.033)
0.090 (0.041) 0.162 (0.079)
Straight
6.51 (0.74)
1.627 (0.014) 1.659 (0.060)
0.668 (0.003) 0.128 (0.049)
0.131 (0.001) 0.126 (0.048)
0.131 (0.001) 0.097 (0.062)
Left
7.21 (0.97)
1.212 (0.182) 1.243 (0.041)
1.104 (0.094) 0.611 (0.021)
0.946 (0.149) 0.611 (0.021)
0.572 (0.308) 0.043 (0.015)
140
M. Theodorou, R.A. Clement / Journal of Neuroscience Methods 161 (2007) 134–141
sure of avoiding spurious fixed points (So et al., 1997) and was used in the previous application to nystagmus, (Clement et al., 2002). However, the use of the tensor to exclude spurious fixed points is not needed with eye movement data because each cycle of such data has to be examined to exclude any artefacts due to attentional changes in fixation or blinks, and anomalous cycles would be picked up at this stage of the analysis. Another possible reason is that the technique has been described as a language for neuronal dynamics in which the behaviour of the system is described by a hierarchy of periodic orbits (So et al., 1998). But finding higher order periodic orbits becomes increasingly difficult, because there are relatively fewer examples, and the evaluation of the eigenvalues becomes sensitive to the examples used. The difficulty of obtaining reliable estimates of the eigenvalues simply, does not support the use of the fixed point technique as a method for building a quantitative model of a neural system. In the context of nystagmus, the major advantage of the technique is that it enables quantitative descriptions of the effects of experimental manipulations to be obtained despite the variable nature of the waveform (Reccia et al., 1990; Clement et al., 2002; Optican and Miura, 2004). After using the technique to recover examples of the underlying periodic orbits in two cases of nystagmus, we have been able to show that there are significant changes in the peak-to-peak amplitudes of the example cycles with changes in gaze angle. This finding suggests that the technique could be a useful additional technique for analysing other oscillatory motor phenomena, such as the tremor that arises in Parkinson’s disease. The most immediate application of the technique is likely to be in providing reliable outcome measures when treatments are being tested. Furthermore, a recent study of drug therapy for nystagmus (Shery et al., 2006) showed that in some cases the effect of treatment was greater on pendular than on jerk waveforms. By enabling quantitative estimates of the relative instabilities of the jerk and pendular components of the waveform, it may be possible to use the technique to predict which patients will benefit most from the treatment. Application of the fixed point technique to two different forms of early-onset nystagmus has revealed previously unsuspected differences in the dynamics of their eye movements. In the jerk form of nystagmus (subject 1), the orbits of the system repeatedly approach a single unstable fixed point. In the pendular form of nystagmus (subject 2), the orbits of the system cycle between two unstable fixed points. A previous analysis of a model for the mechanism of nystagmus has shown that it is an example of a symmetrical system with a pair of homoclinic orbits, that is orbits which begin and end at the same unstable fixed point (Akman et al., 2005). With appropriate parameter changes, the pair of homoclinic orbits can join at the fixed point to form a single cycle by a process referred to as a gluing bifurcation. In the case of jerk nystagmus, the dynamics correspond to having left and right beating jerk waveforms which join to make a bi-directional jerk waveform. But as shown by the results of subject 2 the fixed point technique also reveals a type of bifurcation which was not expected. Subject 2 has a pair of heteroclinic orbits, that is orbits which join two different fixed points during a complete cycle. On left gaze the orbit contains a left beat-
ing jerk component and on right gaze the orbit contains a right beating jerk component. As is common in pendular nystagmus with foveating saccades, as the gaze angle is shifted from left to right, a bidirectional waveform appears with both left and right beating jerk components. In this case a simultaneous gluing at both fixed points appears to occur. Although the possibility of such a bifurcation remains to be confirmed by numerical simulation and analysis, the provisional finding suggests that the fixed point approach may be a powerful technique for revealing novel bifurcations in biological systems which have not previously been discovered in physics or engineering applications. Acknowledgements We thank Ozgur Akman for helpful discussions. This work was supported by grants from The Ormsby Charitable Trust and The British Eye Research Foundation. References Abadi RV, Broomhead D, Clement RA, Whittle JP, Worfolk R. Dynamical systems analysis: a new method of analysing congenital nystagmus waveforms. Exp Brain Res 1997;117:355–61. Abadi RV, Worfolk R. Harmonic analysis of congenital nystagmus waveforms. Clin Vision Sci 1991;6:385–8. Akman OE, Broomhead DS, Abadi RV, Clement RA. Eye movement instabilities and nystagmus can be predicted by a nonlinear dynamics model of the saccadic system. J Math Biol 2005;51:661–94. Akman OE, Broomhead DS, Clement RA, Abadi RV. Nonlinear time series analysis of jerk congenital nystagmus. J Comput Neurosci 2006;21(2):153–70. Broomhead DS, Clement RA, Muldoon MR, Whittle JP, Scallan C, Abadi RV. Modelling of congenital nystagmus waveforms produced by saccadic system abnormalities. Biol Cybern 2000;82:391–9. Clement RA, Whittle JP, Muldoon MR, Abadi RV, Broomhead D, Akman O. Characterisation of congenital nystagmus waveforms in terms of periodic orbits. Vision Res 2002;42:2123–30. Dell’Osso LF. Tenotomy and congenital nystagmus: a failure to answer the wrong question. Vision Res 2004;44:3091–4. Dell’Osso LF, Daroff RB. Congenital nystagmus waveforms and foveation strategy. Doc Ophthalmol 1975;39:155–82. Dell’Osso LF, Flynn JT, Daroff RB. Hereditary congenital nystagmus. An intrafamilial study. Arch Ophthalmol 1974;92:366–74. Dickinson CM, Abadi RV. The influence of nystagmoid oscillation on contrast sensitivity in normal observers. Vision Res 1985;25:1089–96. Gancarz G, Grossberg G. A neural model of the saccade generator in the reticular formation. Neural Netw 1998;11:1159–74. Harris C. Nystagmus and eye movement disorders. In: Taylor D, editor. Paediatric ophthalmology. 2nd ed. Blackwell; 1997. p. 869–96. Harrison PK, Tattersall JE, Clement RA. Periodic orbit analysis reveals subtle effects of atropine on epileptiform activity in the guinea-pig hippocampal slice. Neurosci Lett 2004;357:183–6. Laptev D, Akman OE, Clement RA. Stability of the saccadic oculomotor system. Biol Cybern 2006;95:281–7. Miura K, Hertle RW, FitzGibbon EJ, Optican LM. Effects of tenotomy surgery on congenital nystagmus waveforms in adult patients. Part I. Wavelet spectral analysis. Vision Res 2003;43:2345–56. Optican LM, Miura K. Tenotomy and congenital nystagmus: a null result is not a failure, for “It is not the answer that enlightens, but the question”. Vision Res 2004;44:3095–8. Reccia R, Roberti G, Russo P. Spectral analysis of pendular waveforms in congenital nystagmus. Ophthalmic Res 1989;21:83–92. Reccia R, Roberti G, Russo P. Computer analysis of ENG spectral features from patients with congenital nystagmus. J Biomed Eng 1990;12:39–45.
M. Theodorou, R.A. Clement / Journal of Neuroscience Methods 161 (2007) 134–141 Shery T, Proudlock FA, Sarvananthan N, McLean RJ, Gottlob I. The effects of gabapentin and memantine in acquired nystagmus: a retrospective study. BJO 2006;124:839–43. Slutzky MW, Cvitanovic P, Moul DJ. Deterministic Chaos and noise in three in vitro hippocampal models of epilepsy. Ann Biomed Eng 2001;29:607–18. So P, Francis JT, Netoff TI, Gluckman BJ, Schiff SJ. Periodic orbits: a new language for neuronal dynamics. Biophys J 1998;74:2776–85.
141
So P, Ott E, Sauer SJ, Gluckman BJ, Grebogi C, Schiff SJ. Extracting unstable periodic orbits from chaotic time series data. Phys Rev E 1997;55:5398–417. Stark J, Brewer D, Barenco M, Tomescu D, Callard R, Hubank M. Reconstructing gene networks: what are the limits? Biochem Soc Trans 2003;31: 1519–25. Worfolk R, Abadi RV. Quick phase programming and saccadic re-orientation in congenital nystagmus. Vision Res 1991;31(10):1819–30.