The correlation dimension: A useful objective measure of the transient visual evoked potential?

Abstract

Visual evoked potentials (VEPs) may be analyzed by examination of the morphology of their components, such as negative (N) and positive (P) peaks. However, methods that rely on component identification may be unreliable when dealing with responses of complex and variable morphology; therefore, objective methods are also useful. One potentially useful measure of the VEP is the correlation dimension. Its relevance to the visual system was investigated by examining its behavior when applied to the transient VEP in response to a range of chromatic contrasts (42%, two times psychophysical threshold, at psychophysical threshold) and to the visually unevoked response (zero contrast). Tests of nonlinearity (e.g., surrogate testing) were conducted. The correlation dimension was found to be negatively correlated with a stimulus property (chromatic contrast) and a known linear measure (the Fourier-derived VEP amplitude). It was also found to be related to visibility and perception of the stimulus such that the dimension reached a maximum for most of the participants at psychophysical threshold. The latter suggests that the correlation dimension may be useful as a diagnostic parameter to estimate psychophysical threshold and may find application in the objective screening and monitoring of congenital and acquired color vision deficiencies, with or without associated disease processes.

Introduction

Transient visual evoked potentials (VEPs) are often analyzed by examination of the morphology of their components, such as negative (N) and positive (P) peaks and their latencies (e.g., Campbell & Maffei, 1970; McCulloch & Skarf, 1991, 1994; Odom et al., 2004). However, methods that rely on subjective component identification may be unreliable when dealing with responses of complex and variable morphology (Boon, Suttle, & Henry, 2005; Kulikowski, 1977); therefore, objective methods are also useful. 

One commonly used objective method of VEP analysis is Fourier analysis, which decomposes the time series into a linear combination of sine and cosine waves of particular frequencies and amplitudes. It has been successfully used to characterize both steady-state (e.g., Norcia, Clarke, & Tyler, 1985; Victor & Mast, 1991) and transient VEPs (Boon et al., 2005). Fourier analysis is effective if the time series is from a linear dynamical system (with or without small added noise). A dynamical system is a mathematical model of a system that evolves in time. In this case, even if the time series is apparently irregular in the time domain, the linearity is manifest as sharp peaks in the frequency domain. On the other hand, if the time series is from a nonlinear dynamical system, with or without small added noise, or a random process (large noise) then the time series is generally irregular in both the time domain and the frequency domain. However, there are other methods of time series analysis, based on phase space reconstructions, which can reveal structure in time series from nonlinear dynamical systems when compared with time series from random processes. One of the most popular such methods of analysis involves correlation dimension estimates from phase space embeddings (defined later) using the Grassberger–Procaccia algorithm (Grassberger & Procaccia, 1983a, 1983b). In this paper, we will refer to phase space embeddings and correlation dimension estimates as nonlinear analysis to distinguish them from the intrinsically linear Fourier analysis. It should be noted that nonlinear analysis in this paper refers to the mathematical concept of nonlinear dynamical analysis and does not refer to nonlinearities in response functions (e.g., amplitude as a function of stimulus contrast), unless otherwise stated. 

Background concepts: Dimension, fractal dimension, correlation dimension

In mathematics, “dimension” measurements are used to quantify the space filling properties of a set. The space in which the set is embedded, which also has an associated dimension, may be an abstract mathematical space or a real space. There are many kinds of dimensions (Camastra, 2003). For example, the topological dimension of a point is 0, a line has dimension 1, a surface has dimension 2, a volume has dimension 3, and so on. Importantly, the topological dimension is always an integer. Now imagine a very wiggly line; intuitively we might think that it could fill up more space than a smooth line, but both lines have topological dimension of 1. In fact a wiggly line that is wiggly enough could indeed fill up more space than a smooth line, while still not filling up as much space as a surface. A more informative dimension measurement for such a line is a fractal dimension measurement that permits noninteger values. The line is then deemed to be a fractal if its fractal dimension is greater than its topological dimension (Mandelbrot, 1967). A fractal dimension measurement that has found widespread use in nonlinear dynamics time series analysis is the correlation dimension. 

In nonlinear dynamics, the time evolution of the system cannot be often obtained as a closed formula but it can be represented as a path (the phase space trajectory) in an abstract mathematical space called phase space. Using the time series for just one of the components of the system, it is possible to reconstruct a path that shares the same invariant properties (such as dimension) as the phase space trajectory. This process of reconstruction is called phase space embedding (see Figure 1 and 1). Phase space embedding and dimension measurements of reconstructed paths have been widely used to distinguish between irregular time series generated by noise and irregular time series generated by a nonlinear deterministic (one where the future state of the system is fully prescribed by the present state and the nature of its interactions) dynamical system (see, e.g., Henry, Lovell, & Camacho, 2001; Sprott, 2003). The key to this discrimination is that phase space trajectories of deterministic dynamical systems usually evolve towards a particular set of coordinates called an attractor and the dimension of the attractor is less than that of the full phase space. On the other hand, reconstructed phase space trajectories from noisy time series would not be expected to reduce to a low dimensional attractor but would be expected to fill the embedding space (subject to data limitations as described in 1). 

 

Figure 1

 

The top half of the animation shows a plot of a time series (amplitude of potential difference as a function of time). The time series traces out a sine wave in the time domain. To reconstruct the path of the phase space trajectory of the time series in two-dimensional phase space, a comb is created which has m = 2 teeth, where m is the dimension in which the reconstructed phase space trajectory is to be embedded (higher dimension embeddings would require higher numbers of teeth in the comb). The spacing between the teeth of the comb is equivalent to the time delay (τ). The points of the time series which touch the teeth of the comb define the (x, y) coordinates of the reconstructed phase space trajectory. When the animation is activated, the comb slides forward along the time series and a red dot traces out the path of the reconstructed coordinates in two-dimensional phase space (depicted in the bottom half of the animation). The reconstructed phase space trajectory is approximated as a broken line. It can be seen that the reconstructed phase space trajectory of a sine wave embedded in two dimensions is a one-dimensional object.

There are many different dimension measurements (examples of fractal dimensions include the Hausdorff dimension, information dimension, correlation dimension, e.g., Grassberger & Procaccia, 1983b) that may be performed on the reconstructed path and these different measurements may yield different values reflecting different types of correlations in the reconstruction. The simplest dimension measurement to perform is that of the correlation dimension (Grassberger & Procaccia, 1983a, 1983b), and this is the most widely used dimension measurement in applied time series analysis. The correlation dimension is an example of a fractal dimension measurement because it permits noninteger values. For future reference, we also note that the correlation dimension is related to the minimum number of variables needed to model the system's behavior in phase space. 

Further technical details on phase space trajectories, phase space embedding, and the correlation dimension have been provided in 1. 

Application of the correlation dimension to electrophysiological signals

It has previously been shown that the correlation dimension can provide a useful objective measure of physiological signals and may be used to supplement or complement Fourier (spectral) analysis. For example, when applied to EEGs recorded during sleep, the correlation dimension has been found to complement spectral analysis as a discriminator between stages I and II sleep (stage I was significantly higher in correlation dimension than stage II). On the other hand, spectral analysis could be used to discriminate between stage II and slow wave sleep (Fell, Röschke, Mann, & Schäffner, 1996). This ability to differentiate between electrophysiological signals arising from different electrophysiological states is also useful in the diagnosis of diseases of the nervous system 

An example of a disease for which the correlation dimension may find application is Parkinson's disease. In Parkinson's disease, the central nervous system degenerates (specifically, pigmented neurons of the substantia nigra and other brain stem dopaminergic groups are lost) (Mandir & Vaughan, 2000; Merck & Co. Inc., 2005) resulting in disruption of communication between nerve cells in the brain and causing impairments in movement. In Parkinson's disease, the correlation dimension of the EEG is more sensitive than spectral analysis at discriminating electrophysiological activity that is both task and disease dependent (Stam et al, 1995; Müller, Lutzenberger, Pulvermuller, Mohr, & Birbaumer, 2001). The correlation dimension also provides a means of differentiating EEGs arising from dementia due to Alzheimer's disease from that of vascular etiology (Jeong, Chae, Kim, & Han, 2001). This is useful as different dementias may produce similar disturbances in the EEG (Yagyu et al., 1997). 

The correlation dimension also shows promise as a diagnostic parameter for neural changes that occur in glaucoma (Schmeisser, 1993), a disease that results in progressive loss of the peripheral visual field. Glaucoma is associated with an impaired ability to respond to high flicker rates both psychophysically (Holopigian, Seiple, Mayron, Kory, & Lorenzo, 1990; Tyler, 1981) and electrophysiologically (Holopigian et al., 1990). This impaired ability is associated with the loss of magnocellular neurons (Holopigian et al., 1990; Schmeisser & Smith, 1989). Schmeisser (1993) found that flicker VEPs derived from individuals with glaucoma were lower in correlation dimension than from one healthy individual. 

As well as in the VEP, the correlation dimension has been found to be relevant to the visual system when applied to time series generated by eye movements, where it was found to be correlated with reading ability (Schmeisser, McDonough, Bond, Hislop, & Epstein, 2001). 

Time series derived from magnetoencephalogram (MEG) recordings have also been analyzed using the tools of nonlinear dynamical time series analysis and early results suggest that the correlation dimension may differentiate healthy brain states from those with schizophrenia (Kotini & Anninos, 2002) or brain tumours (Antoniou et al., 2004). 

In addition, the correlation dimension may be used as a parameter to describe the energy signature of different neural tasks. It has been found that the correlation dimension of the EEG (Mölle et al., 1996; Xu & Xu, 1988) and VEP (Arle & Simon, 1990; Schmeisser, 1993) can change significantly under different mental (Arle & Simon, 1990; Mölle et al., 1996; Stam, van Woerkom, & Pritchard, 1996; Price, Moore, Humphreys, & Wise, 1997; Tomberg, 1999; Xu & Xu, 1988), visual perceptual (Arle & Simon, 1990; Müller, Lutzenberger, Preissl, Pulvermüller, & Birbaumer, 2003), and motor tasks (Liu, Yang, Yao, Brown, & Yue, 2005). 

In the present study, the correlation dimension is examined to determine its relevance to the functioning of the visual system and its usefulness in analysis of the transient VEP by its application to the transient chromatic VEP. 

Methods

Design

It was hypothesized that the VEP time series is a manifestation of one dynamical variable (a variable that evolves over time) over a period of time to enable the VEP to be analyzed in the framework of nonlinear analysis (Schreiber, 1999). The time series was then tested for evidence of nonlinearity (correlation dimension and surrogate testing, as described below). The stimulus chromatic contrasts selected for analysis were 42%, two times psychophysical threshold (2T%), psychophysical threshold (T%), and zero (0%). A contrast level that is a multiple of psychophysical threshold (2T%) was useful in investigating the relationship between the correlation dimension and perception of the stimulus. 

Most of the transient chromatic VEPs analyzed in this study have been presented previously (Boon et al., 2005) where they were analyzed in terms of VEP component latency, amplitude, and their Fourier power spectra and compared with psychophysical chromatic contrast thresholds of each of these participants. 

Participants

Eleven adults participated in the study (mean age: 28 years, range 26–33, five males and six females). All of the participants had normal color vision (defined as a pass on screening with the Ishihara Pseudoisochromatic plate test), no ocular abnormalities, no history of amblyopia, and 6/6 or better Snellen acuity (both eyes), with refractive correction if required. Informed consent was obtained from each participant. 

Apparatus

Visual stimuli were generated using a VSG 2/5 card (Cambridge Research Systems, Rochester, England) and presented on a gamma corrected Sony CPD-G500 21-inch Trinitron color monitor. The frame rate of the monitor was 100 Hz. Luminance and chromaticity output over time were verified using a luminance meter (Minolta LS-110) and CRT tri-colorimeters (Minolta TV CA 2150 and Minolta CA-100) (the uncertainty of measurement was ±5% for luminance and ±0.010 for chromaticity). Ag/AgCl electrodes (9 mm) and a commercially available EEG paste were used for recording the VEPs. Single-channel transient VEPs were recorded using the Medelec Synergy averager (Radiometer Pacific, Sydney Australia), which generated time series of 1000 data points per analysis window. Recording was triggered by the external VSG 2/5 card. 

Stimuli

The stimuli were heterochromatic, subjectively isoluminant, obliquely oriented (45° or 135°) sine wave gratings of spatial frequency 1 cycle per degree (cpd) presented in square wave pattern onset–offset at a stimulus temporal frequency of 2 Hz and a duty cycle of 1:4 (stimulus on for 100 ms then off for 400 ms) on a background and surround of the same mean chromaticity ( x = 0.305, y = 0.310) and luminance (20 cd/m ²). The grating stimuli were squares with sides subtending an angle of 5°. At maximum contrast, the colors were of CIE chromaticity coordinates x = 0.380, y = 0.270 (magenta), and x = 0.230, y = 0.350 (cyan), the L-, M-, and S-cone contrasts were 0.16, −0.40, and 0.07, respectively, and the pooled root mean square (rms) cone contrast theoretically produced by the stimulus at maximum contrast was 0.25 (McKeefry, Parry, & Murray, 2003). In the present study, chromatic contrast level was expressed as a percentage of maximum available chromatic contrast. The spatial, temporal, and chromatic parameters of the stimuli were chosen to stimulate the L–M chromatic contrast system preferentially (Kulikowski, McKeefry, & Robson, 1997; McKeefry, Russell, Murray, & Kulikowski, 1996; Mullen, 1985; Murray, Parry, Carden, & Kulikowski, 1987; Rabin, Switkes, Crognale, Schneck, & Adams, 1994; Suttle & Harding, 1999) while minimizing the effects of chromatic aberration (Flitcroft, 1989). 

Isoluminance was determined using heterochromatic flicker photometry (mean of 11 determinations of minimal flicker of heterochromatic magenta-cyan stimuli with a spatial frequency of 1.0 cpd, field size 2.5°, and square wave phase reversed at a 15-Hz temporal frequency). The color ratio (L _color1 / (L _color1 + L _color2)) was varied by the participant. In 9 of the 11 participants, VEPs were recorded at their individual isoluminance color ratios. The remaining two participants were tested at the average isoluminance color ratio of a larger group of participants because their VEPs were originally recorded as part of a different experiment. The use of an average isoluminance value for those two participants may introduce luminance cues into their responses but this is not relevant to the question of whether VEPs may be described by the correlation dimension, so they are included in that part of our analysis. Data from these two participants, however, were not used to investigate correlations between the correlation dimension, Fourier power, and chromatic contrast. 

VEP procedure

The VEP recordings to be analyzed, using nonlinear techniques, were intended to be representative of those VEPs that would be measured using instrumentation commonly available in the clinical environment. As such, single-channel transient VEPs were averaged from at least 30 sweeps (Rabin et al., 1994) and limited to a recording of 1000 data points (limitation of the Medelec Synergy Averager). Recording sites were Oz (active site), Cz (reference site), and Fz (ground), according to the 10–20 system (Odom et al., 2004). Impedance was below 8 kΩ for all recordings. The signals were band-pass filtered (1–50 Hz) as most of the signal of the pattern VEP was expected to have a frequency of less than 30 Hz and to avoid aliasing of the signal (artificially low frequencies) (Odom et al., 2004). The sampling frequency was 1000 Hz, above the optimal sampling rate for the EEG of normal participants (512 Hz; Jing & Takigawa, 2000). The amplifier range was 2.5 mV. The effect of filtering on correlation dimension estimation has been examined by Lo and Principe (1989) who concluded that frequencies above 50 Hz should be filtered to remove unwanted noise frequencies as the EEG does not “have appreciable energy above 50 Hz.” Lo and Principe also noted that filtering of a signal, be it the EEG or white noise, results in a decrease in the estimated correlation dimension compared to an unfiltered signal thus filtered noise should also be used in any comparison between the correlation dimension of noise and signal. To remove electricity mains noise, notch filtering was activated at 50 Hz. Artifact rejection was set at 50 μV. The averaged VEP in response to gratings at each contrast was recorded at least twice to check repeatability. The criterion for repeatability in this experiment was the presence of an N–P complex with latencies within 10% of the longer latency on successive VEPs. In addition, the averaged sweeps were analyzed in their entirety; thus, they were not truncated in order to exclude baseline fluctuations in the electrophysiological signals nor to extract only the chromatic contrast components within the sweep. N-peak latency was recorded as the average of repeatable N-peak latencies. Otherwise, no N-peak latency was recorded. VEPs were not excluded from Fourier analysis or nonlinear dynamical analysis on the basis of lack of repeatability in the time domain. 

Psychophysical procedure

Psychophysical threshold was determined using a two alternative (between grating orientation; 45° or 135°) forced choice design, in which threshold was calculated as the mean of contrasts at 12 reversals in a 3-down 1-up staircase, with a 3-dB step size. On each trial, duration of presentation was 1 s (two duty cycles). The stimuli were presented centrally as isoluminance is known to vary with retinal location (Bilodeau & Faubert, 1999; Kulikowski, Robson, & McKeefry, 1996). 

Fourier analysis

Fourier power in the present study refers to the sum of the two highest peaks of the response present in the Fourier power spectrum of each VEP. This single parameter for Fourier power was used as it was found to be more closely correlated with chromatic contrast and psychophysical chromatic contrast threshold than other parameters in a previous study (Boon et al., 2005). In this method, at least two exponential windows are used to create two power spectra per VEP. The two power spectra are then compared and the peaks that are repeatable in frequency and relative height between the two power spectra were regarded as being more likely to be significant to the response. This method is based on the understanding that different kinds of windows produce different artifactual peaks or lobes in the frequency domain (Harvey & Cerna, 1993). Exponential windowing more clearly reveals the power spectra of transient responses than use of the default rectangular window (Harvey & Cerna, 1993). This feature of windowing was exploited such that if two different windows are used, then it is likely that the peaks that are due to the characteristics of the window will shift in frequency and relative height in the power spectra depending on the window hence will not be repeatable. On the other hand, the peaks that are repeatable are less likely to be characteristic of the window and more likely to be significant to the response. Thus, in the present study, the two window functions applied to the time series during Fourier analysis were exponential, with time constants (time for the exponential function to decay to one half of its initial value) of 512 ms and 256 ms. 

Rationale for method of correlation dimension estimation

Most applications of Grassberger and Procaccia's (1983a, 1983b) method of nonlinear analysis to electrophysiological time series have been to stationary time series such as an extended EEG or steady-state response. Its use in the assessment of transient responses has been less thoroughly explored despite the fact that its corollary, transient chaos, has been the subject of theoretical examination (Carroll, Pecora, & Rachford, 1987; Grebogi, Ott, & Yorke, 1986). However, preliminary investigations suggest that the fractal dimension of transient responses is a valid parameter (Arle & Simon, 1990; Dhamala, Lai, & Kostelich, 2001; Dhamala, Pagnoni, Wiesenfeld, & Berns, 2002; Jánosi & Tél, 1994). 

The fractal dimension of transient responses should be interpreted in a slightly different way than nontransient responses as the former do not represent the ultimate state to which the dynamical system will evolve. In other words, the phase space trajectory that is reconstructed from a transient time series may only reflect a limited subset of dynamical activity of a wider system and may not be ultimately characteristic of the underlying system as a whole (Jánosi & Tél, 1994). Thus, the correlation dimension of a transient will characterize the dynamics of a signal settling into, or just entering, the trajectory of an attractor (Jánosi & Tél, 1994) and may differ in correlation dimension from the underlying system (Arle & Simon, 1990; Jánosi, Flepp, & Tél, 1994). Jánosi and Tél (1994) noted that this does not preclude the measurement of the fractal dimension of transient time series. Arle and Simon (1990) reasoned that the underlying system will differ in some way when the EEG displays a series of spike transients compared to when it does not at baseline and found that this was reflected in the fractal dimension. Therefore, even though the spike transients within the EEG discussed in Arle and Simon's work are not precisely the same as transient evoked potentials, their reasoning that a difference in the underlying system during baseline and transient activity may be reflected in the fractal dimension also indicates that the correlation dimension is potentially useful in the analysis of visual evoked potentials. 

In addition, the application of Grassberger and Procaccia's (1983a, 1983b) algorithm to averaged transient electrophysiological data is not a common procedure. Previously, nonaveraged time series have been analyzed in preference to averaged time series as it has been reasoned that averaging will remove much of the variability that the nonlinear approach is designed to explore. However, in this study, the aim was not to determine the underlying structure of the visual system as manifest by the visual evoked potential but, instead, was to determine whether the correlation dimension may be a useful objective parameter of the transient visual evoked potential. Moreover, as the amplitude of the transient chromatic visual evoked potential is very small compared to the background fluctuations that may be related to other sources, such as the EEG (e.g., Friendly, 1993), it is necessary to average sweeps in order to access the transient VEP. 

Because of the limitations inherent in recording the transient VEP using clinically available instrumentation, such as sampling rate, number of data points, and band-pass filtering (Lo & Principe, 1989), it is acknowledged that the correlation dimension estimated in the present study will not necessarily be equivalent to the number of dynamical variables in the subset of electrical activity represented by the transient chromatic VEP. However, it may provide a reproducible, objective parameter that may be suitable for mapping of dynamical activity across stimulus contrasts, subject to the limitations outlined above. 

Correlation dimension analysis

The correlation dimension was assessed using Grassberger and Procaccia's (1983a, 1983b) algorithm for each VEP and their surrogate time series. 

Two control surrogate data sets were created for each VEP and were designed to be random and nondeterministic but with the same mean variance and power spectra as the original VEP data (equivalent to filtered (colored) noise (Theiler & Rapp, 1996)). As the software we used required the length of the time series to be analyzed to a power of two, each 1000 data point VEP time series was padded to a 1024 data series by the addition of 24 zeros. The fast Fourier transform function was applied to the padded VEP to determine the power spectrum. For the first half of the Fourier power spectrum, the surrogate transforms were given the original amplitudes at random phases. The second half of the Fourier power spectrum was taken to be the complex conjugates of the first half reflected about the midpoint to ensure that the surrogate data would be real. Then an inverse fast Fourier transform function was applied to generate the surrogate time series. 

The method of Grassberger and Procaccia's (1983a, 1983b) algorithm was then followed. This involved firstly converting the time series into evolving sets of coordinates in an abstract mathematical space called phase space following the procedure outlined in 11 for phase space embedding. In the present study, two delay times (τ = 4 ms and 6 ms) were used and the correlation dimension, D₂, was estimated for five embedding dimensions, m = 1, 2, 3, 4, and 5. The τ values were chosen based on an exploration of the reconstructed pathways in three-dimensional space arising from different values of τ. The final values of τ were selected so that (1) the reconstructed pathway would not be collapsed or aligned along the diagonal of three-dimensional phase space (Grassberger and Procaccia, 1983b) and (2) the smallest ripples in the VEP time series would not be excluded from the reconstructed phase space trajectory of the VEP. Apart from these considerations, the choice of one delay time over another should not greatly affect the estimated correlation dimension (Mazaraki, 1997; Takens, 1981). 

For each reconstructed phase space trajectory, the distances between all points in the trajectory were calculated and the logarithm of the smallest distance (represented by r _min) and the logarithm of the largest distance (represented by r _max) were computed. A series of ‘bins’ was then created to record the correlation sum, C( r), which is the normalized number of pairs of points with a separation distance of less than a specified distance r. The process of depositing counts of data into bins is analogous to recording counts of the occurrence of events within data in a frequency histogram. In this study, 64 bins (an arbitrary number) were used and the width of each bin was set to ( r _max − r _min) / 64. Thus, from first to last, the separation distances used in the analysis were defined by the series r ₁, r ₂,… r _n, where the radius r _n = r _min + n( r _max − r _min) / 64, where n = 1 to 64. 

In practice, D ₂ can be approximated by D ₂ ∼ log( C( r))/log( r) (Grassberger & Procaccia, 1983a, 1983b) as outlined in 1. Thus, D₂ was calculated as the slope, dlog(C(r))/dlog(r), of the linear scaling portion of the plot of log(C(r)) versus log(r). In this study, the slope of the scaling portion of the plot was calculated by determining the slopes of i consecutive local portions of the plots, with each portion of the plot consisting of k consecutive points (k = 6, 12), and finding the maximum slope out of all the calculated slopes for each k. The slope for each ith portion was determined to be the slope of the line of best fit (least squares method) drawn through the points that were part of the ith portion of the plot (for example, see Figure 2, i = slope number). 

Figure 2

 

Both plots are log(r) versus log(C(r)) for subject JV in response to 42% chromatic contrast for an embedding dimension of 2 and a delay time of 4. The left shows three consecutive slopes (in grey, numbered i = 1, 2, and 3 to the left of the slopes) determined from k = 6 consecutive points on the plot. The right shows three consecutive slopes determined from k = 12 consecutive points on the plot.

The rationale for using the maximum slope as an indicator of the correlation dimension was based on examination of many plots of log( C( r)) versus log( r) and the objective of estimating slope from a region of stable linear scaling. Such a region on the plot should contain successive local slopes of approximately equal magnitude and create a high plateau in a plot of local slope as a function of slope number (for example, see Figure 3). This is because progressively increasing slopes will occur over a range of small log( r) such that the scarcity of pairs of points separated by < r may result in any slope taken from this region being an inaccurate reflection of the phase-space filling properties of the trajectory. On the other hand, progressively decreasing slopes will occur as log( r) becomes so large that C( r) will contain almost all of the pairs of the points in the trajectory so will gradually form a plateau in the plot of log( C( r)) versus log( r) until every pair of points in the trajectory is included. Previous researchers have used the “middle-third” of the curve of log( C( r)) versus log( r) in order to avoid calculating slope from regions of too small and too large log( r), but this method can also result in error if it straddles two kinds of straight-line behaviors (noise and deterministic chaos) (Henry et al., 2001). On exploration of the data, the local maximum slope was found to lie within a range of local slopes of equal magnitude that lay between regions of progressively increasing and progressively decreasing slopes so the maximum slope was used as an indicator of D₂. In the present study, because the aim is to determine the relevance of the parameter to the functioning of the visual system by comparing across stimulus and perceptual conditions and not to estimate the absolute number of dynamical variables of the system, provided the method used is constant for all time series estimated, this method of D₂ estimation from the slope will permit comparison between the D₂'s of the VEPs. In total, four estimates of D₂ were made for each embedding dimension as two τ values and two k values were used. 

Figure 3

 

The plots on the left and right depict the slopes determined using k = 6 and 12 units, respectively, as a function of slope number (i). The steepest slope is taken as an estimate of the correlation dimension, D₂.

The embedding dimension, m, is the upper bound of the correlation dimension, D ₂, that may be estimated from the reconstructed data. Thus, as m is increased, D ₂ should also increase until D ₂ reaches a plateau value. A plateau could be achieved in two different ways; (1) a plateau representing the true dimension of the reconstructed path (a solid sphere would appear to be a two-dimensional object if embedded in a two-dimensional space but remains a three-dimensional object when embedded in any higher dimensional space), (2) a plateau that is an artefact of the finite data size (as the embedding dimension is increased, more and more data points are needed to fill the space). For a finite data set comprising N = 10 ³ points, a plateau value of D ₂ ≥ 2 log ₁₀ N = 6 should be discounted as an artefact of the finite data size (Eckmann & Ruelle, 1992), and so to avoid this ‘false plateau’ we have restricted the dimension of the embedding space to m ≤ 5. In the analysis below, D₂ was plotted as a function of m for m = 1, 2, 3, 4, and 5, and in cases where the plot appeared to plateau, the D₂ values for m = 5 were taken as estimates of the true dimension of the reconstructed path. Where D₂ did not appear to plateau with increasing embedding dimension, it was regarded as a result indistinguishable from uncorrelated noise. In the case where D₂ plateaued, a further test using surrogate data was employed to differentiate between a deterministic response and a correlated noise response (Theiler & Rapp, 1996). 

Nonlinearity of the VEPs containing chromatic contrast (42%, 2T%, and T%) was tested by comparing the estimated correlation dimensions of the VEPs and their surrogate data sets using the paired t-test. Nonlinearity of the 0% chromatic contrast responses, where there was no evoked potential in response to chromatic modulations, was similarly tested by comparison with surrogates of the 0% data. If both the experimental time series data and the surrogate data sets do not yield statistically significantly different correlation dimensions, then the possibility that the VEP is a linear system in combination with random noise may not be excluded. If the estimated correlation dimension of the data is found to be close to 5, the upper bound for reliable estimates of the fractal dimension, then the possibility that the data set represents random noise may not be excluded. If the slope of the plot of D ₂ as a function of m is not found to plateau but continues to increase with increasing m, then again the possibility that the data set represents random noise may not be ruled out. To determine whether the function reached a plateau, a plateau index was defined as D ₂(5) − D ₂(4), where values in parentheses are embedding dimensions of 4 and 5. An index of at least 0.3 (arbitrarily chosen) was taken to indicate no plateau. 

The correlation dimension at m = 5 was compared with the Fourier power of the VEP and the chromatic contrast of the stimulus to look for associations. Because one of the contrast levels selected is a multiple of psychophysical threshold, an association between the correlation dimension of the VEP and perception of the stimulus was also considered. 

Results

The group-averaged (in black) and individual (in grey) VEPs in response to 0%, T%, 2T%, and 42% chromatic contrast are presented in Figure 4. The N-peaks are indicated by the upward pointing arrows. The P-peaks of the N–P complex are indicated by the downward pointing arrows. The shoulder in the N–P complex is indicated by the diagonal arrows. It can be seen that the 42% response is an N–P complex with a shoulder. As the chromatic contrast of the stimulus is decreased, the interindividual similarity of VEP morphology between individuals is reduced and the N–P complex is less clearly defined in the group-averaged VEP, although present for all of the VEPs in response to 2T% and 42%, and six out of the nine individual VEPs in response to T%. There were no repeatable N–P complexes present in any of the VEPs in response to 0% chromatic contrast. 

Figure 4

 

VEPs of individual participants (grey lines, average of 70 to 80 sweeps) and group-averaged (black lines) VEPs in response to (a) 0%, (b) psychophysical threshold (T%), (c) twice psychophysical threshold (2T%), and (d) 42% chromatic contrast. It can be seen that the typical N–P complex is most distinct for the 42% chromatic contrast group-averaged VEPs. The upward pointing arrows indicate N-peaks, the downward pointing arrows indicate P-peaks, and the diagonal arrows indicate the shoulders of the N–P complexes.

Figure 5 shows a worked example of correlation dimension analysis of a VEP from participant BC in response to 42% chromatic contrast. It includes the VEP, a surrogate of the VEP, phase space reconstructions of the time series in three-dimensional phase space (animated), the respective plots of log C( r) versus log r, and m versus D ₂ for τ = 4, 6, and k= 6, 12. The fractal dimension for the dynamical system described by this VEP was eventually estimated at 2.41. It can be seen that, unlike the original VEP data, the surrogate data's corresponding D ₂ versus m plot does not plateau and that the mean correlation dimension at m= 5 is 4.23. 

 

 

 

Figure 5

 

An example of a VEP in response to 42% chromatic contrast from participant BC (left side, top row) and one surrogate (right side, top row), their reconstructed phase space trajectories embedded in three-dimensional phase space (row 2, clicking on the plots will take you to an animation of the reconstruction), their plot of log C( r) versus log r(row 3), the plot of m versus D ₂ for τ = 4, 6, and k= 6, 12 (row 4). Crosses indicate τ = 4, k= 6; circles indicate τ = 4, k= 12; filled diamonds indicate τ = 6, k = 6; filled triangles indicate τ = 6, k = 12 (row 4). Note that the plot of D ₂ as a function of m did not plateau for the surrogate data but did plateau for the VEP data.

The surrogate, chromatic contrast VEP (42%, 2T%, T%), and baseline 0% VEP data were compared in terms of their plateau index values and correlation dimension estimates for m = 5 (see Figure 6). Figure 6 shows box plots depicting these data for all of the VEPs in response to chromatic contrast, all the VEPs in response to 0% chromatic contrast and their respective surrogate data. An index of ≥0.3 (arbitrarily chosen) was taken as an indication that the function failed to plateau. The VEP and 0% data do plateau, as indicated by mean indices of 0.18 (SD 0.11) and 0.18 (SD 0.06), respectively. The surrogate data do not plateau, as indicated by mean indices of 0.95 (SD 0.34) and 0.83 (SD 0.33), respectively. A paired t-test comparison between the VEP, 0%, and their respective surrogate data sets, for both the parameters of plateau index and the correlation dimension (at m = 5), resulted in statistically significant differences between the VEP and surrogate data sets ( p < 0.0001) for both parameters. The correlation dimensions of the VEPs at m = 5 could therefore be regarded as the fractal dimension of the embedded time series, but this was not the case for surrogate data. 

Figure 6

 

Box plots comparing the plateau index values and correlation dimensions for m = 5 for the original VEPs and the electrophysiological signal in response at 0% chromatic contrast and their surrogate data. A plateau index value of ≤0.3 was taken to indicate that the slope does plateau. The median is indicated by the horizontal line intersecting each box. The upper and lower bounds of the boxes represent quartiles and the error bars represent 95% confidence limits for the data. The circles are outlier data.

Plateau index and correlation dimension estimates for the VEPs in response to 0% chromatic contrast were also compared to the 0% surrogates in order to determine whether the 0% chromatic contrast recordings were composed of stochastic or colored noise or whether they also contained a signal. Naturally, in the case of the 0% responses, the signal would not correspond to that of the chromatic contrast system but to other activity, most likely the resting EEG at the Oz location. 

The results agree with our understanding that noise sources of the VEP may not necessarily be stochastic. For example, noise may arise from the recording apparatus (Nunez, 1995; Sprott, 2003; Wong, 1996), alternative sources such as the EEG (e.g., delta waves, Moskowitz & Sokol, 1980; alpha waves, Hood & Greenstein, 2003), the pulse (Guy, ffytche, Brovelli, & Chumillas, 1999), or precortical spontaneous activity unrelated to visual function (Ghose & Freeman, 1992). Our results suggest, furthermore, that the variability of the 0% signal is unlikely to be due to stochastic noise. This is even despite the fact that the method of sweep averaging (which assumes that there are random deflections in the background signal that comprise of noise) works to realize the morphology of low amplitude transient evoked potentials (Friendly, 1993). 

The terms fractal and correlation dimension (correlation dimension is one type of fractal dimension as it permits noninteger values) may henceforth be used interchangeably for our VEP (42%, 2T%, T%, and 0%) data as it has been established that the correlation dimensions of the signals are unlikely to be artifacts of noise, unlike the correlation dimensions of the surrogate time series data. 

A strip plot of correlation dimensions of the electrophysiological recordings in response to different stimulus contrasts is presented in Figure 7. The correlation dimensions (mean (SD)) for 42%, 2T%, T%, and 0% are 2.6 (0.19), 2.8 (0.29), 3.1 (0.21), and 2.9 (0.26), respectively. Analysis of variance and paired comparisons (adjusted for multiple comparisons) indicate that the correlation dimensions of each of the stimulus groups are significantly different to each other except for 2T% and 0%. 

Figure 7

 

Strip plot of fractal dimensions of the electrophysiological recordings in response to different stimulus conditions (unfilled circles). The means are indicated by the larger, filled circles in red and the standard deviations by the error bars. Broken lines connect those points belonging to individuals. Note that the highest contrast stimulus had the lowest fractal dimension for all the stimuli for which chromatic contrast was present for all of the individuals.

The relationships between VEP correlation dimension, power, latency (first repeatable N-peak latency—if N-peak latency was not repeatable, it was excluded from analysis) and chromatic contrast are presented in Figure 8 ( N = 9). The data points that belong to each individual participant have been joined by lines. Separate lines have been used for each participant. In Figures 8A, 8B, and 8C, correlation dimension, Fourier power, and latency as a function of chromatic contrast of the stimulus are presented. In Figures 8D and 8E, points depicting fractal dimension as a function of Fourier power and latency are depicted with dotted lines joining the data points in order of the rank of the powers for Figure 8D and in order of the ranks of chromatic contrast for Figure 8E. 

Figure 8

 

VEP correlation dimension as a function of (a) chromatic contrast, (d) power, and (e) N-peak latency. (b and c) Power and latency as a function of chromatic contrast, respectively. The data points that belong to an individual participant have been joined by lines. Separate lines have been used for each participant. The lines join the data points in order of the rank of the chromatic contrasts in panel A, powers in panel D, and latencies in panel E. For the overall grouped data, there were significant relationships between the VEP correlation dimension and (a) chromatic contrast (Spearman's Rho: −0.88), (d) power (Spearman's Rho: −0.50), and (e) latency (Spearman's Rho: 0.41)

Rank correlation analysis (Spearman's Rho) of correlation dimension as a function of power and chromatic contrast was used for the data as a whole and for each individual. There was a significant negative relationship between the Fourier power and the correlation dimension of the VEPs overall for the group (Spearman's Rho: −0.502, df = 26, p = 0.007 (two-tailed)). However, this relationship was significant in only two out of the nine participants (Spearman's Rho). There was a significant positive relationship between the N-peak latency and the correlation dimension of the VEPs overall for the group (Spearman's Rho: 0.41, df = 23, p = 0.040 (two-tailed)) and was also strong for each of the participants (Spearman's Rho: 1.0). There was a significant negative relationship between the chromatic contrast of the stimulus and the correlation dimension of the VEP (Spearman's Rho: −0.889, df = 27, p < 0.0001 (two-tailed)), as shown in Figure 8A. The relationship was significant in seven out of the nine participants (Spearman's Rho). 

Furthermore, in order to determine whether there was agreement of the rankings of VEP correlation dimension (as a function of power, latency, and chromatic contrast) across participants, the nonparametric Kendall's W (coefficient of concordance) was calculated. For correlation dimension as a function of power, Kendall's W was 0.30 and the asymptotic significance was 0.09, indicating that there was little to no agreement in the ordering of correlation dimension as a function of power across participants. The nonsignificant asymptotic significance value indicates that there is a random distribution of ranks of the correlation dimension as a function of power. For correlation dimension as a function of latency, Kendall's W was 0.88 and the significance was 0.0004, indicating a high agreement in the ordering of correlation dimension as a function of latency. For correlation dimension as a function of chromatic contrast, Kendall's W was 0.790 and the asymptotic significance was 0.001, indicating that there was high agreement in the ordering of correlation dimension as a function of chromatic contrast across participants. The significant asymptotic significance value indicates that there is a nonrandom distribution of ranks of the correlation dimension as a function of chromatic contrast. 

Thus, the results of these analyses, considered together, indicate that there is a high negative correlation between chromatic contrast and VEP fractal dimension and a high positive correlation between N-peak latency and VEP fractal dimension and that both these relationships were evident in most of the individual participants. For the group, the VEP fractal dimension increases with decreasing power but this relationship was not consistent across participants. 

Discussion

These results show that the transient chromatic and 0% VEPs behave significantly differently from linear systems combined with random noise. In the framework of nonlinear analysis, this indicates that the chromatic and 0% VEPs recorded here are not purely linear, purely stochastic, or nonlinear time series combined with high amplitude stochastic noise. Instead they are nonlinear time series (which may comprise of a mix of both linear and nonlinear responses) and any stochastic noise that is present in the recorded VEPs is relatively low in amplitude. The evidence for nonlinearity of the response agrees with models of the visual system which describe visual function as a combination of linear and nonlinear activity (Carandini et al., 2005; Crognale, Switkes, & Adams, 1997). However, the fact that the 0% signal was not found to be stochastic indicates that a plateau in the plot of correlation dimension and embedding dimension is not necessarily characteristic of response to patterned visual stimuli. 

The fractal dimensions were negatively correlated with chromatic contrasts above zero, on both an individual and group basis to reach a maximum at psychophysical threshold. Stimulus chromatic contrast is a physical property but because one of the contrast levels selected is a multiple of psychophysical threshold, from the perceptual point of view this may indicate that there is also a negative correlation between the correlation dimension of the VEP and the visibility of the stimulus. These findings suggest that it may be possible to use the fractal dimension of the VEP as an indicator of psychophysical threshold, and to determine the presence of a transient response in a signal, in agreement with work by Arle and Simon (1990). However, the fractal dimensions of VEPs in response to different levels of stimulus chromatic contrast were not uniformly statistically significantly different in magnitude from the 0% (unevoked response). In fact, although 42% and T% fractal dimensions were statistically significantly different from the 0% fractal dimensions, the 2T% responses were not. This indicates that the fractal dimension may be useful for differentiating between similar percepts arising from different stimulus levels along a continuum such as chromatic contrast. However, it is possible that two entirely different perceptual states, such as when viewing a field containing chromatic contrast then viewing a uniform gray field, may be described by equal numbers of dynamical variables (even though the dynamical variables themselves may differ), hence similarity in correlation dimension between 2T% and 0% responses. Despite this limitation, our findings suggest that the correlation dimension may be used as an aid in the estimation of psychophysical threshold since in most of the participants the stimulus contrast with the highest correlation dimension was psychophysical chromatic contrast threshold. 

The lack of a significant difference between the 2T% and 0% fractal dimensions could indicate a lack of sensitivity of the fractal dimension but a closer examination of how fractal dimension changes across the stimulus conditions suggests otherwise. The fractal dimension peaks at psychophysical chromatic contrast threshold and then dips when there is no chromatic contrast present in the signal (the difference between T% and 0% correlation dimension is significant). This suggests that a change must occur within the signal at 0% to which the fractal dimension is sensitive. As the main change at 0% is that the visual system is no longer processing color contrast, this may involve a different population of cells activated by the stimulus and a perceptual change. Overall, this implies that the fractal dimension is related to stimulus contrast but may be more closely reflective of perceptual states or the cortical electrophysiological activity that leads to the percept. Other methods of visualizing activity at the level of the cortex have also found that cortical activity can be more closely related to the percept rather than the stimulus (Kreiman, Fried, & Koch, 2002). Kreiman et al. (2002) showed that during suppression of an incoming retinal image, cells which customarily signal for that stimulus in the middle temporal lobe do not fire even though the stimulus is present. On the other hand, work by Crognale et al. (1997) serves as a reminder that the source and type of cortical activity is also important when looking for correspondence between cortical electrophysiological and psychophysical responses. Crognale et al. found that it is possible for suprathreshold differences in amplitude, in the pattern-reversal chromatic VEP to be completely uncorrelated with suprathreshold visibility and suggested that this may occur if any increase in VEP amplitude is related to a natural resonance of the cortex (e.g., alpha waves) rather than activity within the visual system. Moreover, Crognale et al. noted that VEP latency was correlated with psychophysical measures and this is consistent with our finding of better correspondence between N-peak latency and both fractal dimension and chromatic contrast than the Fourier-derived VEP amplitudes. 

Interestingly, the mean fractal dimension of the VEPs in response to 42% chromatic contrast was significantly lower than that of the unevoked response (0%). This finding is similar in magnitude but opposite in direction to Arle and Simon's (1990) finding that the fractal dimension was significantly higher (by approximately 0.6) when VEPs were present than when they were absent by approximately 0.6. (In the present study, only the relative change in correlation dimension with stimulus condition, rather than the absolute correlation dimensions, are compared between studies as the estimate of absolute fractal dimension will vary depending on the specifics of the algorithm used, and this is an approach that others have taken; Anokhin, Lutzenberger, Nikolaev, & Birbaumer, 2000; Kirsch, Besthorn, Klein, Rindfleisch, & Olbrich, 2000.) This apparently opposite result, however, is still consistent with Arle and Simon's findings as they reasoned that a transient response should be distinguishable from a baseline response as it will have a different underlying dynamical structure to the baseline response. Their method did not make any presuppositions about the expected morphology of the transient nor the direction of the change (higher or lower fractal dimension). In addition, Arle and Simon noted that their technique could differentiate three transient responses in series, but not one or two, against an unevoked recording. Our results indicate that two transient responses in series may also be differentiated from an unevoked electrophysiological recording. 

Changes in fractal dimension as a function of a stimulus property have been observed in the EEG (Liu et al., 2005; Müller et al., 2003). However, it is acknowledged that comparisons with the VEP must be made with care, since the EEG (1) does not reflect only a subset of the activity, (2) may employ a different reference electrode location, with impact on the signal (Tomberg, Desmedt, & Ozaki, 1991), and (3) may involve different active electrode locations (EEGs recorded from different sites can show opposite movements in fractal dimension in response to the same fractal task; Mölle et al., 1996; Xu & Xu, 1988). However, there are similarities between the EEG and VEP. For example, both the VEP and the EEG may be composed of transients (Arle & Simon, 1990) and both may be recorded under different conditions of stimulation (e.g., Fell et al., 1996; Mölle et al., 1996). Moreover, the changes in dimension in both the EEG and VEP data will arise specifically from those electrodes overlying those parts of the brain required to process the stimulus and not from unrelated areas (Price et al., 1997; Tomberg, 1999). These are similarities that will allow comparisons between EEG and VEP fractal behavior. 

The fractal dimension of the EEG has previously been related to stimulus magnitude (Liu et al., 2005; Müller et al., 2003). For example, Müller et al. (2003) calculated the fractal dimension of EEGs from participants who were observing a moving stimulus at two different levels of complexity. For both conditions, the participants were asked to fixate a stationary fixation spot on a wall just behind the pendulums so the task was the same for the participants and only the properties of the stimulus had changed. Müller et al. found that the correlation dimension of the EEG at all electrode sites was higher for the more complex moving stimulus than for the simpler moving stimulus by <1.0 fractal dimension. Similarly, Liu et al. (2005) examined how the fractal dimension of the EEG changes when more than two levels of stimulus magnitude within a continuum are presented. Liu et al. used stimuli of graded magnitude and found that increasing the level of force applied in a hand-grip task resulted in an increase of the measured fractal dimension of the EEG derived from motor-related sites of the cortex. The change in fractal dimension, estimated using several algorithms, was in the range of 0.23 to 0.5. These findings are consistent with those of the present study, in that they too demonstrated a small but significant change in fractal dimension with a graded change in stimulus property. 

Why should fractal dimension of the VEP change with stimulus properties? The fractal dimension of the EEG has been hypothesized to reflect the number of coupled independent populations of neuronal assemblies involved in the generation of the EEG hence differences in fractal dimension may reflect differences in the number of populations of neuronal assemblies that are activated (Müller et al., 2003). As the VEP is also an electrophysiological signal that is created from a subset of the population of neurons that would contribute to the EEG, this hypothesis may also apply to the VEP. Based on this hypothesis, the number of independent neuronal assemblies may be fewer for VEPs recorded in response to highly visible chromatic stimuli (e.g., at 42%) than to stimuli that are only just visible (e.g., at T%). Another possibility is that the overall populations of neurons that mediate the processing of chromatic contrast and underpin the transient chromatic VEP are the same for both highly visible and just visible stimuli but that they may function more coherently (or more strongly thus will appear more coherent to our measurement apparatus) for highly visible stimuli. If two or more populations of cells act coherently, their activity would collapse into one activity which would appear as if they were only one independent neuronal assembly in this kind of an analysis. A possible dynamical picture is that when the stimulus is highly visible, regardless of how many independent neuronal assemblies there are, most of them will signal the presence of chromatic contrast. This largely coherent signal will be manifest as at least one population of neuronal assemblies. On the other hand, when the stimulus is only just visible, some neurons may signal the presence of chromatic contrast coherently (one independently active population) while others may either signal the absence of chromatic contrast or the presence of another feature of the stimulus (another independently active population) thereby increasing the number of independent, although coupled, signals in the transient chromatic VEP to a minimum of two. This hypothesis may also account for the relative decrease in fractal dimension for the recordings in response to 0% chromatic contrast compared to those VEPs recorded from psychophysical threshold. In theory, the complete absence of a chromatic contrast stimulus may cause most of the neurons to behave coherently as they signal the absence of chromatic contrast or, alternatively, the presence of another feature in the stimulus. 

The above interpretation is consistent with function of the cortical visual system as measured using cortical dye signaling. Grinvald, Arieli, Tsodyks, and Kenet (2003) visualized activity of cortical cell assemblies in response to optimal, nonoptimal, and no visual stimulation in the cat using voltage sensitive dyes. In the cat, it has been found that single neurons fire and behave as part of an assembly, even when spontaneously firing in the absence of stimulation. Grinvald et al. found that the relative amplitude of coherent activity in cortical cell assemblies in areas 17 and 18 in the absence of visual stimulation was slightly less but approximately equal to the amplitude of coherent activity in the presence of visual stimulation (Grinvald et al., 2003). In addition, the cells that were not time locked to the stimulus were also engaged in coherent on-going spontaneous activity that was almost equivalent in amplitude to the activity of the cell assemblies that were time locked to the stimulus. They also found that the more ideal the stimulus was for a particular cell assembly, the more likely it was that a neuron would fire coherently with other neurons in the same assembly. As the stimulus becomes less ideal for that particular cell assembly, it will fire less but it will still fire coherently with other cells in its assembly. When there was no stimulus at all, the neurons were found to behave the same way as if there were an optimal stimulus for that assembly, very coherent activity within an assembly but slightly reduced in amplitude. Grinvald et al. speculated that in the absence of activity, there is switching between internal cortical states (represented by cortical cell assemblies) which means that the firing rates of different cell assemblies at rest would cycle through maxima and minima at different times. Thus, although in combination the activity may appear stochastic in the time domain, it may be made up of multiple coherently functioning cell assemblies (Grinvald et al., 2003). 

Other experimental findings have also supported a hypothesis of multiple, coupled neuronal assemblies during sensory processing. Experimental work on the thalamus (MacLean, Watson, Aaron, & Yuste, 2005) and neocortex (Cossart, Aronov, & Yueste, 2003) and their effect on spontaneous cortical neuronal activity reveals that cell assemblies within the cortex function coherently according to a predefined structure and that outside influences, be they from the thalamus or inputs from sensory stimulation, will activate certain circuits, or pathways, within the structure. Interestingly, both MacLean et al. (2005) and Cossart et al. (2003) note that their experimental findings can be described by a dynamical attractor model, where an attractor is a “semi-stable point in the dynamics of a network that could implement an internal mental state … this model predicts that if a threshold of circuit component activation is achieved, the entire circuit ensemble becomes activated.” In the context of Maclean et al.'s and Cossart et al.'s comments, the pathways in the cortex that code for chromatic contrast are active whether or not there is chromatic contrast input, only the amplitude is different. In response to an optimal stimulus (such as the 42% contrast stimulus used here), the dynamical state of the cortex may switch to an internal dynamical state of high amplitude activity and very coherent activity (Cossart et al., 2003; MacLean et al., 2005). However, less optimal stimuli will result in less frequent, albeit still coherent, firing by that cell assembly or visual pathway. When there is no chromatic contrast stimulus, there is still visual input but it is not specific to the chromatic visual system so systems other than the chromatic pathways may be activated and act coherently. At T% chromatic contrast, the stimulus is less than optimal half the time and absent half the time so there would be neuronal cells that would normally be active at the 0% state which will be active in addition to cells that would normally be active at the 42% state but firing much less frequently. The presence of two very different patterns of coherent signaling could result in a minimum of at least 2 coherent variables dominating the VEP. The 42% dynamical state and the 0% dynamical state will dominate to a greater extent, respectively, resulting in a reduced fractal dimension at 42% and 0% compared to T%. Thus, the results we obtained for the 0% condition are consistent with previous work that has shown that spontaneous activity in the cortex is not stochastic but is instead precisely patterned (Grinvald et al., 2003) and may reflect internal cortical states. 

Fractal dimension of the VEP was less strongly related to Fourier power than to stimulus contrast or N-peak latency. This may be due to known interindividual variability of VEP amplitude (e.g., Rabin et al., 1994). As it is known that VEP amplitude may not be correlated with visual stimulus characteristics for some individuals (Boon et al., 2005; Jenkins, Douthwaite, & Peedle, 1985; Sokol and Moskowitz, 1985) and since the fractal dimension is moderately correlated with a stimulus characteristic such as chromatic contrast, fractal dimension might be expected to be weakly correlated with VEP amplitude. Fractal dimension and Fourier power are negatively correlated. This could arise because of the method used to calculate Fourier power in the present study, which was to sum the two highest powers of the power spectra derived from the two windows. This calculation favors signals that have narrow band power spectra where the power is concentrated in one or two frequencies at the most, rather than signals which produce broad band responses of relatively similar amplitudes over a wide range of frequencies. The Fourier power spectra show that when an optimally visible stimulus is presented at a specific stimulus temporal frequency, the power is concentrated at frequencies related to the stimulus temporal frequency while at psychophysical threshold, power is approximately of equal strength for those frequencies related to the stimulus temporal frequency and those that signal other nonchromatic features of the visual stimulus. Thus, the more complex response at psychophysical threshold would be expected to be associated with lower Fourier amplitudes in view of the fact that the signal would not be concentrated within one or two frequencies at psychophysical threshold. At 0% chromatic contrast, following the above reasoning, if more than one population of cells fires in their individual patterns but with a reduced firing rate compared to when there is a strong stimulus, this could result in multiple populations of cells firing low amplitude patterns, which would result in a broadband signal of small amplitude in the frequency domain. 

In the present study, it must be noted that only the chromatic visual system was explored and only the dimension of chromatic contrast was varied, such that all other stimulus parameters were held constant. Thus, the behavior of the correlation dimension of VEPs in response to stimuli of different temporal, spatial, and chromatic parameters is not known. However, the findings hold promise that the correlation dimension may also be a useful object parameter to describe evoked potentials designed to tap into other visual responses. In particular, it would be interesting to examine whether visual stimuli designed to activate more than one pathway result in VEPs of higher fractal dimension than stimuli designed to activate only one pathway. 

Conclusion

The present study has built on previous work (Arle & Simon, 1990; Schmeisser, 1993) that has explored the usefulness of the correlation dimension as an objective parameter of the VEP and as an indicator of visual system function (Schmeisser et al., 2001). The correlation dimension of the transient VEP was found to be relevant to the function of the visual system as it was negatively correlated with a stimulus property (chromatic contrast) and a known linear measure (the Fourier-derived VEP amplitude) and positively correlated with VEP response latency. It was also found to be related to visibility of the stimulus such that fractal dimension reached a maximum for most of the participants at psychophysical threshold. The latter suggests that the correlation dimension may be a useful objective diagnostic parameter in the estimation of psychophysical threshold and may find application in the objective screening and monitoring of congenital and acquired color vision deficiencies, with or without associated disease processes. Its ability to characterize evoked potentials other than the transient chromatic VEP and the flicker VEP is yet to be explored and may result in wider benefits on further investigation. 

Appendix A

The phase space trajectory and phase space embedding

Phase space embedding is a process whereby a time series is converted into a series of coordinates in a reconstructed phase space. The set of reconstructed coordinates defines a reconstructed phase space trajectory. If the original time series is one component of a deterministic dynamical system and certain conditions (discussed below) are fulfilled, then invariant properties such as the dimension of the phase space trajectory are faithfully captured in the reconstructed phase space trajectory (Takens, 1981). If the reconstructed phase space has m dimensions, then each reconstructed phase space coordinate X_i is an m component vector which is obtained from the time series y(t₁), y(t₂),… by the prescription X_i = (y(t_i), y(t_i+ τ), y(t_i+ 2τ),…, y(t_i + (m − 1)τ). Here m is called the embedding dimension, τ is a constant known as the delay time, and the index i denotes ordering in time. 

Figure 9 shows a single point in three-dimensional phase space created using the embedding procedure. The phase space embedding process can be visualized as a comb passing along the data. The number of teeth of the comb is equivalent to the embedding dimension and the distance between the teeth corresponds to the delay time. When the comb is positioned at the first point in the data, the points which the teeth of the comb touch define one component each of the single point in phase space represented by an m component vector. The comb is then slid along the data to create a series of phase space points that makes up the reconstructed phase space trajectory. 

Figure 9

 

The state of a dynamical system can be reconstructed from a single observable time series by a process of phase space embedding. The coordinates of a single point in three-dimensional phase space are shown here reconstructed from the values of three points in the original time series separated by a delay time ( τ). These three points can be visualized as points on the time series which touch the end of a comb with m = 3 teeth. The evolution of the dynamics of the system in phase space can be traced out by sliding the comb across the data.

The correlation dimension

Grassberger and Procaccia (1983a, 1983b) introduced a fractal dimension measurement called the correlation dimension, D₂, to characterize the phase space filling properties of attractors.¹ It is obtained by covering the set with boxes of a given size (r) and then computing the probability p_i(r) (equivalent to the relative frequency in sufficiently large data sets) of having a point of the set in the ith such box. The correlation dimension is defined by 

$D_2=\underset{r\to 0}{\lim }\sum _i\frac{\log \left({\sum }_i{p}_i\right(r{)}^2)}{\log r}.$

(A1)

The quantity ∑_ip_i(r)² is the probability of finding a pair of points in a box of size r. The so-called Grassberger–Procaccia algorithm (Grassberger & Procaccia, 1983a, 1983b) provides a computationally efficient way to implement this measurement. For small values of r, the probability of having a pair of points in a box of size r is equal to the probability of having a pair of points with separation distance less than r. For sufficiently large data sets, of number N, this latter probability is given by the correlation sum (also termed the correlation integral) 

$\begin{array}{ccc}C\left(r\right)& =& \underset{N\to \infty }{\lim }\frac{1}{N^2}\{\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{p}\mathrm{a}\mathrm{i}\mathrm{r}\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{s}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{s}\mathrm{e}\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}<r\}\hfill \\ & =& \underset{N\to \infty }{\lim }\frac{1}{N^2}\sum _{i,j=1}^N\theta (r-|{\mathrm{X}}_i-{\mathrm{X}}_j\left|\right)\hfill \end{array},$

(A2)

where θ is the Heaviside function (a discontinuous step function which has a value of either 0 or 1 and may be defined as θ(x) = 0 for x < 0, and as θ(x) = 1 for x > 0) which acts as a counter of the number of pairs of points with separation <r when combined with ∣X_i – X_j∣ (the separation distance between two points on the attractor, X_i and X_j). The multiplier 1/N² is included to normalize the count by the number of pairs of points on the attractor (without double counting). 

Grassberger and Procaccia (1983a, 1983b) established that for small values of the separation distance (r), the correlation sum grows like a power such that 

$C\left(r\right)\sim r^{\mathrm{D}2}.$

(A3)

(In Grassberger & Procaccia's, 1983a, 1983b, paper, r= l and D₂ = v.) The “∼” symbol in this expression is used to indicate that this is not an exact equality but is a scaling relation that is expected to be valid for sufficiently large N and small r. After taking logarithms of each side of the scaling relation and rearranging terms, we have the result 

$D_2\sim \frac{\text{log}\left(C\right(r\left)\right)}{\text{log}\left(r\right)}.$

(A4)

 

In practice, D ₂ is deduced from the slope of the straight line scaling region in a plot of log( C( r)) versus log( r). The scaling relation is a type of dimensional measurement that in this case describes how a reconstructed phase space trajectory fills a given phase space. Thus, the Grassberger–Procaccia algorithm (1983a, 1983b) has been widely used to measure the correlation dimension of reconstructed phase space trajectories from time series. In these applications, the time series is embedded at increasing embedding dimensions m until the measured D₂ reaches a plateau, as revealed in a plot of D₂ versus m. One of the most significant limitations in these applications is that the time series needs to be large enough to exhibit genuine scaling behavior at large embeddings. A rule of thumb (Eckmann & Ruelle, 1992) is that the upper limit on a reliable measurement of D₂ is given by D₂ = 2 log₁₀N. Thus, any plateau in the plot of D₂ versus m at values of D₂ higher than this bound should be discounted as artefacts of the finite data size. 

In summary, the correlation dimension is an objective parameter that describes the phase-space filling properties of phase space embeddings of time series data. It can be used to approximate the number of equations, or dynamical variables, required to model the activity represented by a single observable time series and in this sense is an indicator of the dynamical system's relative complexity. 

Acknowledgments

The first author was supported by the Research Training Scheme of the Commonwealth of Australia. 

Commercial relationships: none. 

Corresponding author: Mei Ying Boon. 

Email: [e-mail]. 

Address: School of Optometry and Vision Science, University of New South Wales, Sydney, 2052, Australia. 

Footnote

Footnotes

¹  In Grassberger and Procaccia's (1983a, 1983b) papers, they use the term “strange” attractors to describe the fractal attractors of nonlinear deterministic dynamical systems.