Fractal analysis design for distinguishing subject characteristics on motor control of neck pain patients

Literature reflects that there is considerable amount of research has been conducted concerning motor control and postural control of human body. Research findings also reveal the importance of stability and variability of movement in relation to spine, in particular. However, the evaluation method on stability and variability of spinal movement is limited on the linear methods in revealing the significance of movement. The knowledge on the nonlinearity and the dynamic features of spinal movement is still inadequate in describing the motor control and musculoskeletal characteristics.


Introduction
Literature reflects that there is considerable amount of research has been conducted concerning motor control and postural control of human body. Research findings also reveal the importance of stability and variability of movement in relation to spine, in particular. However, the evaluation method on stability and variability of spinal movement is limited on the linear methods in revealing the significance of movement. The knowledge on the nonlinearity and the dynamic features of spinal movement is still inadequate in describing the motor control and musculoskeletal characteristics.
The aim of this chapter is to explore into area of understanding how human body maintain posture in a dynamic manner under the context of non-analytical method of spinal motor control and the kinematic resultant of musculoskeletal system. Based on the research outcomes, contribution can also be made into the application domain of evaluating motor control characteristics as reflected from different subject profiles.

Literature review
In general, there are two types of postural equilibrium, namely, static and dynamic, to maintain a state of balance in which all forces acting on the body (Kandel, 2000). Static equilibrium is to allow the body rests in an intended position. Dynamic equilibrium is to maintain a balance while progressing through an intended movement. In the case of spinal stability, it is defined as the ability to maintain intervertebral and global torso equilibrium (Granata, 2006). Kinematic variability of postural control has been studied for postural stability in many researches. Small biomechanical or neuromotor disturbances are continuously perturbing the system, causing kinematic variance. Therefore, stability can be observed when the measured kinematics appears to be attracted toward the posture of static equilibrium. Measurements on motor feedback can provide substantial information about the postural changes and kinematic variability. However, the dynamic features of motor control particularly the association between adjacent measurement time points as well as the serial ordering of kinematic outcomes are not well considered.

www.intechopen.com
Kinematic variability usually include the measurements of center-of-pressure (COP) trajectory traveled within a period of time, root mean square and ellipse area of COP. These measurements are commonly assumed to represent stability when compared healthy controls to pathological cases (Radebold, 2001). From the perspective of human balancing performance, these measurements on kinematic variability are the output from the disturbance in amplitude observed. They represent the total disturbance distance, average amplitude of output variance and the size of involuntary movement (Prieto et al., 1996) Therefore, the dynamic features on spatio-temporal dimensions of input disturbances should be taken into account when investigation of variability and stability are considered (England & Granata, 2005;Li et al., 2005) A number of attempts have been made to analyze more specifically the dynamic features of biological time series. Sets of analysis methods have been designed to reveal the hidden fractal properties of time series. The methods differ in many points, but they, in general, attempt to assess, in multiple temporal intervals of various lengths, the dispersion or the displacement of variables (Delignières et al., 2003;Rosenstein et al., 1993) One of the most classical method is the Hurst rescaled range analysis (HRRA). It is a common method to extract fractal features in economics (Peters 1994), geophysics (Nickolaenko et al., 2000), biology (Hoop et al., 1995) and motor control (Yamada 1995). Hurst's approach is based on the assessment of the range of displacements of the locally integrated time series for each interval. A second method, detrended fluctuation analysis (DFA) has been developed especially for the analysis of biological time series, in particular heartbeat time series (Plastisa & Gal, 2008) and gait stride interval series (Hausdorf et al., 1996;Dingwell & Cusumano, 2010), which are often highly nonstationary. Collins and De Luca (1993) have proposed the stabilogram diffusion analysis (SDA) in the analysis of center-of-pressure (COP) trajectories during unperturbed stance. The idea of SDA is rooted from the study of Brownian motion by Einstein (1905). The Einstein's relation has been generalized by Mandelbrot and van Ness (1968) to account for a family of stochastic processes that they termed fractional Brownian motion (fBm). In HRRS, DFA and SDA, they have the same interpretation of results and seem very similar. However, HRRA and DFA include integration which distinguishes from SDA. They avoid bi-linear characteristics of SDA and show linear plot.

Stability measures
A number of attempts have been made to analyze more specifically the dynamic features of biological time series. Sets of analysis methods have been designed to reveal the hidden fractal properties of time series. The methods differ in many points, but they, in general, attempt to assess, in multiple temporal intervals of various lengths, the dispersion or the displacement of variables (Delignières et al., 2003;Rosenstein et al., 1993)

Hurst rescaled range analysis
One of the most classical method is the Hurst rescaled range analysis (HRRA). It is a common method to extract fractal features in economics (Peters, 1994), geophysics (Nickolaenko et al. 2000), biology (Hoop et al., 1993) and motor control (Yamada, 1995). Hurst's approach is based on the assessment of the range of displacements of the locally integrated time series for each interval. A time series of N numbers as x(t) is divided into non-overlapping intervals of length n. The integrated series X(t, n) within each interval is calculated as in the following equation. The difference between maximum and minimum of the integrated data X(t, n) for each interval is defined as the range R n .
The range is then divided for normalization by the local standard deviation S n for the original series. This calculation is repeated over all possible interval lengths. In practice, the shortest length is suggested to be around 10 data points, and the largest is N/2. Finally, the resealed ranges R/S are averaged for each interval length n. The R/S is related to n by a power law: where a is a constant and n is the interval length Because H R/S cannot exceed 1.0, the particular application of the method to biological bounded series is not possible if H > 0.5 meaning a persistent series.

Detrended fluctuation analysis
A second method, detrended fluctuation analysis (DFA) has been developed especially for the analysis of biological time series, in particular heartbeat time series (Plastisa & Gal, 2008) and gait stride interval series (Hausdorf et al., 1996;Dingwell & Cusumano, 2010), which are often highly non-stationary. A time series with N numbers as x(t) is integrated for each t the accumulated departure from the mean of the whole series.
This integration step maps the original time series to a self-similar process. Next, the vertical characteristic scale of the integrated time series is to be measured. The integrated series is divided into non-overlapping intervals of length n. In each interval of length n, a least squares line is fit to the data, representing the trend in that interval. The integrated time series X(t) is then detrended by subtracting the local trend X n (t) given by the regression. For a given interval length n, the characteristic size of fluctuation for this integrated and detrended time series is calculated.
The computation is repeated over all possible interval lengths to provide a relationship between F(n) and interval length n. In practice, the shortest length is suggested to be around 10 data points, and the largest is N/2. Typically, F(n) increases with interval length n. The F(n) is related to n by a power law.
() Fn a n   where a is a constant and α is the scaling exponent The scaling exponent α is expressed as the slope of the log-log plot of F(n) as a function of n. It can vary between 0.0 and 1.5. In the original formulation, the DFA is conceived to be applied directly on raw data. A value α = 0.5 indicates the data are uncorrelated random process, i.e., white noise. α = 1.0 represents pink or 1/f noise, which lies at the boundaries between stationarity (α < 1.0) and non-stationarity (α > 1.0). When α is 1.5, the original series is a Brownian motion. Higher values are mathematically obtainable, up to 2.0, for persistent series. Nevertheless, the reliability of such high exponents is still unknown. The relationship between exponent H and α is given below.

Stabilogram diffusion analysis
Collins and De Luca (1993) have proposed the stabilogram diffusion analysis (SDA) in the analysis of center-of-pressure (COP) trajectories during unperturbed stance. The idea of www.intechopen.com SDA is rooted from the study of Brownian motion by Einstein. It is shown that the mean square displacement is related to the time interval Δt .  The Einstein's relation has been generalized by Mandelbrot and van Ness (1968) to account for a family of stochastic processes that they termed fractional Brownian motion (fBm). The scaling law is proposed. This scaling exponent H can be calculated from the slope fitting to a log-log plot of mean square displacement versus time span. An important feature of fBm is the presence of correlations between past and future movement. Such motion exhibits long-memory processes, and each value is dependent upon the global history of the data series. The correlations indicate that, on average, the fluctuations as recorded from one time span are statistically similar to the fluctuations on other time span. The interpretation of the scaling exponent is that the value H = 0 indicates the data are white noise which can be interpreted as a random signal with a flat power spectral density. A value of H = 0.5 represents the classical Brownian motion which is the random movement. In the case when H > 0.5, there exists persistence phenomenon which means the direction of future movement tends to be positively correlated with the direction of current movement.
While H < 0.5, it is associated with anti-persistent phenomenon in which the direction of future movement direction tends to be negatively correlated with the direction of current movement.
The plot also demonstrates a transitional behavior which is distinguished by a transition point named as critical point representing a particular value of time interval Δt c . The critical point is defined as the end of a short-term region of the plot. Rougier (1999) has proposed an automatic determination method for identifying the transition point. The logarithm of the resultant curvature-diffusion plot and a theoretical straight line representing the logarithm of a pure stochastic process with the initial point at the former plot are determined. If the spinal curvature variability behaves in a persistent manner, its distance from the pure stochastic process line will increase. However, the distance from the pure stochastic process line will decrease if the spinal curvature variability behaves in an anti-persistent manner. The point with maximum vertical distance between the spinal curvature variability and the pure stochastic process line indicates the transition between the persistent and anti-persistent processes and is therefore identified as the transition point.
In the aforementioned experiments conducted by Collins and De Luca (1993), the COP trajectory has been found to have critical time approximately at Δt = 1s. Similar results can also be obtained by subjects attempted to actively balance an inverted pendulum (Treffner & Kelso, 1995).
The plot has been interpreted by postulating an open-loop control mechanism. It is also suggested that the short-term region can refer to exploratory processes by considering persistence as information gathering, (Riley et al., 1997). Subsequent researches have shown systematic evolutions of the parameters of the model, for instance, H exponents for shortterm process, coordinates of the point of inflexion, etc., under the manipulation of factors such as vision, leaning, or haptic touch (Riley et al., 1997;Rougier et al., 2001).

Spinal curvature and kinematic measurement
In predicting and quantifying the spinal curvature, researches have been done to develop various biomechanical models. The spine is a multi-joint structure with non-linear geometry, rather than a rigid body. This fact has to be considered in biomechanical models of the spinal model. Imaging techniques, such as X-ray, computer tomography (CT), or magnetic resonance imaging (MRI) have been the most widely used techniques for obtaining spinal geometry. However, these methods are invasive and not recommended for all subjects for prolonged and repeated exposure. Besides, these methods are costly, inappropriate for field measurements and require licensed technicians for operation. Therefore, non-invasive method would be an advantage for evaluating the spinal geometry, based upon easily obtained measures of trunk posture.
There are several non-invasive methods for spinal curvature evaluation described in literature. These include skin markers, external marker photography or videography, electromagnetic devices, electrogoniometer, accelerometer, flexible tape measurement, ultrasonic digitizer, etc. In compared to image-based methods (CT, MRI, etc.), the benefits of non-invasive methods of spinal curvature prediction can be well distinguished, provided that accurate prediction can be obtained. The criteria of the method would be based on tool or device that is applicable to work environment as well as laboratory settings. The method should also allow quick measurements that permit screenings of subject groups.

Research method
This chapter draws upon previous empirical research and theoretical background to establish a modeling and analytical framework on non-analytical method, which results in key parameters that describe the dynamic features during postural analysis of cervical spine. A fractal analysis model is developed for postural analysis based on extraction of spinal curvature. The research demonstrates how the model can be employed and to evaluates clinical application with neck pain patients. The methodology was attempted to be applied into potential practical contribution by means of the dynamic features concluded after computation. The aspect trying to touch upon here was to differentiate various subject characteristics. The hypothesis was that subjects who carried different characteristics would exhibit different dynamic features in control.
The objective of this experiment was to try distinguishing the dynamic features of patients and normal subjects by means of the fractal analysis with respect to cervical spine control. Experimental setup was prepared for capturing cervical movements during cervical spine retraction training. Neck pain patients and normal subjects were invited to conduct the data acquisition. The normal subjects were reported to have no history of neck pain lasting for more than three days during the last one year time. The patients suffered from mechanical neck pain including myofascial neck pain and degenerative changes. The neck pain happened for at least six weeks. There was no radiating symptom noted.
The acquisition sequence was divided into stages based on self-adopted posture, before, during and after training. Positional data along cervical spine was first captured. Angular data on curvature was then extracted corresponding to the flexion and extension of spine in different regions. Data were then put forward to the computation of fractal analysis designed. In distinguishing the subject characteristics between a patient and normal subject, hypothesis was made on the dynamic features. This included critical time, diffusion coefficient and Hurst exponent. Evaluation on these parameters demonstrated that the two types of subject, who carried different characteristics, were exhibiting different dynamic features. www.intechopen.com

Experiment
The acquisition sequence was divided into four stages, named "S -self-adopted posture", "I -upright posture before training", "T -training stage" and "F -upright posture after training". During the "S" phase, subjects were requested to keep a sitting posture that they adopted in daily living. The capture was lasted for 10 seconds in each trial. A total of 10 trials were collected with 15 seconds rest in between consecutive trials. In the "I" phase, subject were requested to keep an upright sitting posture as steady as possible for 10 seconds in each trial. Again, a total of 10 trials were collected with 15 second rest in between consecutive trials. After that, a neck retraction exercise was taught. After the subjects had managed to conduct the exercise by themselves, the "T" phase was started to capture their performance during the exercise. The capture was lasted for 10 seconds in each trial and there were 10 trials of data acquired. In between each consecutive trial, there were 15 seconds of rest. At the final stage, the subjects were requested to keep an upright sitting posture as steady as possible for 10 seconds in each trial. A total of 10 trials were conducted with 15 seconds rest in between consecutive trials.

Optical motion capture
An optical motion capture system was setup for the purpose of data acquisition of subject movement. The system was optical camera digital system manufactured by the Vicon (UK). Each camera was attached in the front a visible red LED ring-lights which emitted narrow band visible red light with fixed frequency customized by the capturing software. The frequency that we used for this experiment was set to be 100 Hz, meaning that for each second, the system was able to capture 100 frames of image as appeared within the view of camera. In order to capture the movement of subject, spherical markers, with 3M (USA) Scotchlite optical reflective material taped on surface, were used to attach on skin proximal to bone prominences. The trajectories of markers were determined by the system from consecutive image frame sequences captured. The three dimensional coordinates of the markers in space were then calculated from various camera angles to result in X, Y and Z coordinates. The accuracy of marker trajectories in three dimensional space was within 0.3 mm as determined by control setup of a dummy subject.
A total of 6 markers were attached to the skin proximal to bony prominences on the subject for the measurement. The marker placements were illustrated in the following

Fig. 2. Placement of markers
For each marker, the X, Y and Z coordinates were captured according to the physical space location in the standard unit of millimeter. On the other hand, a dummy subject was also setup as a control to evaluate the digital noise characteristics of the system relative to the subject data captured. On the dummy, the same number of markers was attached according to the various positions located over the subject body.
The data acquisition was controlled by an utility software from the Motion Analysis Corporation called EVaRT, which has an interface for visualization of data in real-time. The data were exported in a file format named TRC. It was a plain text format with each row corresponding to the number of frames during the data acquisition, while the columns recorded the individual position according to the axis of each marker.

Development of SDA model for spinal curvature
An attempt was made to develop the stabilogram diffusion analysis (SDA) model for spinal curvature. The objectives were to develop a methodology suitable and valid for analyzing spinal curvature through SDA model. The result was then determined if similar control mechanism could be found when analyzing spinal curvature as compared to center-ofpressure (COP) appeared in previous section. First of all, an experimental environment was setup for capturing subject data as described in previous section. The subject was given a set of instructions to perform during the experiment. Data was then capture through the markers attached to the skin of subject. With the numerical data, a number of steps were developed to determine the spinal curvature. On the other hand, the programming procedures for computation of SDA were implemented. The curvature data were then used for the application of SDA computation. Fractional Brownian motion (fBm) characteristics were then extracted, together with a number of graphical representations of data. Finally, these dynamic features were compared. The aim was to determine if the purposed method of developing SDA model for spinal curvature was valid in distinguishing subject characteristics. With the captured data of markers in X, Y and Z coordinates, equations were derived for determination on angular displacement of the head relative to shoulder. First of all, the local coordination systems for head and shoulders were found by the corresponding vector spaces of markers.  The sampling rate of data in this case was 100 Hz and the cut-off frequency was 4 Hz. Zero phase digital filtering was applied from a built-in MATLAB function named filtfilt( ).
The inclined angle was calculated from the filtered raw data of each marker position. The inclined angle, spinal curvature, could be illustrated into a time plot according to the time dimension.

Curvature diffusion plot
The data was then studied as a one-dimensional random walk with the assumption of stochastic process. First of all, a time interval t was defined as the time difference between data points. Given two adjacent intervals, the shortest length was 1/100 s, and the largest was 10 seconds. The squared difference between any two data values was calculated according to various data intervals. By grouping the squared differences together according to each data interval, the mean square angle <ΔX 2 > was calculated from the average over the number of entries making up the group. For example,

Diffusion coefficient
In posturographic investigation, it would be impractical to have subjects perform for extended periods of time. Physiological factors such as fatigue would tend to interfere the results. In this study, 10 trials were collected for a subject, the analysis was therefore averaged sets of results derived from these trials. Specifically, stabilogram-diffusion plot was computed by averaging ten curves to obtain a resultant plot for a particular subject.
Trajectory along time could be quantified by a nonfinite integer or factional space dimension in providing a quantitative measurement of evenness. For the one dimensional random walk with stepwise displacement X, a diffusion coefficient D was calculated from the slopes of the resultant linear-linear plots of mean square angle versus time interval curves, i.e., <ΔX 2 > vs t.
where <ΔX 2 > was the mean square angle, which was the arithmetic mean of ΔX 2 ; D denoted the level of stochastic activity and was the half slope of curvature-diffusion plot

Critical time
The short-term diffusion coefficients D s was computed from the slops of the lines fitted to the short-term region. The critical point (Δt c , <ΔX 2 > c ) is defined by the intersection of the lines fitted to the end of short-term region of the plot. The computation was derived from the difference between a logarithmic plot of a pure stochastic process and the logarithmic plot of curvature-diffusion curve. The critical time was defined by the time interval of maximum difference. (Rougier, 1999)

Hurst exponent
On the other hand, scaling exponents H s was computed from the resultant log-log plot of the above curvature-diffusion curve. In all the above cases, the slopes were determined by utilizing the method of least squares to fit straight lines through defined portions of the plots. Similarly, the short-term region was defined for scaling exponents. The level of correlation between past and future increments was indicated by scaling regime of H. When H was equal to 0.5, it denoted that the process is totally random and was a classical Brownian motion. When H was greater than 0.5, the stochastic process was positively correlated and exhibited persistence behavior. In this case, a Brownian particle moving in a particular direction will tend to continue in the same direction. Conversely, H smaller than 0.5 denoted the stochastic process was negatively correlated and exhibited anti-persistence behavior. In such a case, an increasing trend in the past implied a decreasing trend in the future.

Results analysis
We hypothesized that the spinal curvature exhibited dynamic features and open-loop control behavior as modeled by SDA using fractional Brownian motion method. Moreover, different dynamic characteristics could be distinguished between normal and pathological subjects.
To verify the hypothesis, the first dynamic feature was taken as the critical time point. First of all, the critical time of a patient would have a larger value than a normal subject. From the experimental results, the critical time of patient was consistently having larger value than normal subject in all the four phases of capture. between patient and normal subject were consistent with the hypothesis across the four phases. Exceptions were only found at the short-term diffusion coefficient of phase "S" and "F".

Acknowledgement
The author would like to acknowledge and extend the gratitude to Pricilla and Miko who have contributed their knowledge in discussions and made these research findings possible. www.intechopen.com