Analysis of motor performance inParkinson’s disease through LTIdynamical systemsbyAhmad AshooriA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Electrical and Computer Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)December 2014c© Ahmad Ashoori 2014AbstractMotor performance in Parkinson’s disease is marked by a variety of impairments, and is indicativeof faulty feedback mechanisms in the brain that arise due to a lack of dopamine. Quantifying motorperformance as well as investigating the correlation of motor behaviour and clinical measures hasbeen sought as in important avenue to study disease pathophysiology. We propose the use of lineardynamical systems to finely characterize motor performance, and for possible interpretations ofneurological phenomena, such as compensatory mechanisms that can be cast in the framework ofcontrol theory. While coarse measurements of motor performance exist (e.g., maximum velocity,mean velocity, root mean square error), they often cannot distinguish between qualitatively differenttypes of movement. (Consider for example, the wide variety of step response behaviors possiblein second order systems which may have the same settling time, but a wide range of dampingratios.) In contrast, LTI models provide a dynamical mapping of the sensorimotor transformationrequired in tracking tasks. Model parameters of linear dynamical systems (such as damping ratio,and natural frequency) can more readily describe qualitatively similar motor phenomenon thanexisting coarse measurements.From time-series data of various manual tracking experiments, we first perform a set of prepro-cessing techniques dependent on the particular experimental setup. The pre-processing techniqueshave been carefully evaluated to appropriately address possible noise, nonlinearities, and otherdata irregularities. We then apply existing tools for system identification (e.g., ARX, ARMAX,subspace methods) to identify numerically robust second-order LTI system models. The choiceof system identification method (as well as choice of pre-processing techniques) is critical in de-termining the variability in the resulting model parameters. We identify parameters of the LTImodels which are amenable to comparison across subjects (and across groups), and evaluate theirstatistical significance.The main contributions of this thesis are:1. Selection of pre-processing and system identification techniques for three experimental setupsfor manual tracking experiments2. Control theoretic framework for tremor as a compensatory mechanism that is advantageousiiAbstractin some tracking tasks3. Sensor-based assessment of rigidity via decay rate4. Assessment of cognitive inflexibility through mode detection and delayiiiPrefaceDifferent parts of the work presented in this thesis has been published in or submitted to severalinternational conferences and scientific journals.Results from Chapter 3 were published in:[i] N. Baradaran, S. N. Tan, A. Liu, A. Ashoori, S. J. Palmer, Z. J. Wang, M. Oishi, andM. J. McKeown, “Parkinson’s disease rigidity: Relation to brain connectivity and motorperformance,” Frontiers in Movement Disorders, vol. 4, no. 67, June 2013.Results from Chapter 5 were published in:[ii] A. Ashoori, M. J. McKeown, M. Oishi, “Switched manual pursuit tracking tasks to measuremotor performance in Parkinson’s disease,” in IET Journal of Control Theory & Applications,vol. 5, no. 17, pp. 1970–1977, April 2011.[iii] M. Oishi, N. Matni, A. Ashoori, and M. J. McKeown, “Switching restrictions for stabilitydespite switching delay: application to switched tracking tasks in Parkinson’s disease,” inJournal of Nonlinear Systems and Applications, Special issue on hybrid systems, pp. 16–25,February 2011.[iv] M. Oishi, A. Ashoori, and M. J. McKeown, “Mode detection in switched pursuit trackingtasks: Hybrid estimation to measure performance in Parkinson’s disease,” in Proc. IEEEConference on Decision and Control, Atlanta, GA, December 15–17, 2010, pp. 2124–2130.Results from Chapter 4 were presented at:[v] N. Niksirat, A. Ashoori, N. Baradaran, M. Oishi, and M. J. McKeown,“The Updated Warten-berg Pendulum Test as an Objective Quantitative Measure of Rigidity in Parkinson’s disease,”Canadian Physiotherapy Association’s Congress and the Orthopaedic Division’s Symposium,Montreal, Canada, May 23–26, 2013.In papers [ii], [iv], [v], I pre-processed the raw data, identified LTI model parameters, andinvestigated their significance across groups as well as their correlation with clinical measures. InivPrefaceall papers, I helped interpret significant trends in the context of known phenomena in Parkinson’sdisease. I was the primary author for paper [ii]. For the papers published in medical journals,I provided text relevant to the technical analysis. My supervisors, Professor Oishi and ProfessorMcKeown, provided technical and editing feedback on my contribution to all of the above papers.In paper [i], I did the system identification and computed the model parameters. N. Baradaraninvestigated their correlation with fMRI data and I investigated the correlation of model parameterswith clinical measures. For paper [iii], I identified LTI models and computed Lyapunov functionsto analyze stability. Professor Oishi investigated the switching stability and completed the bulk ofwriting.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Parkinson’s Disease . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Problem Formulation: Manual Tracking Tasks . . . . . . . . . . . . . . . . . . . . . 41.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.4 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91.5 Contributions of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121.6 Organization of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 Pre- and Post-Processing Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . 132.1 Data Pre-Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.1.1 Low-Pass Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.1.2 De-Saturation Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.1.3 Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.1.4 Downsampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.1.5 Time-Lag Removal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.2 Post-Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20viTable of Contents2.2.1 Non-Gaussian Data for Statistical Analysis . . . . . . . . . . . . . . . . . . . 203 Insights Into Compensatory Mechanisms of Tremor . . . . . . . . . . . . . . . . . 223.1 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.2 Pre-Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.3 Model Creation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.3.1 Post-Processing Whitening Filter . . . . . . . . . . . . . . . . . . . . . . . . 273.3.2 Bootstrapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.4.1 Models Quality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.4.2 Correlation with Clinical Measures . . . . . . . . . . . . . . . . . . . . . . . 293.4.3 Group Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 324 Kinematics of Rigidity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.1 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.2 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.3 Model Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 404.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 474.4.1 Variations Across Groups . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 494.4.2 Symmetry of Oscillations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 514.4.3 Medication Effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 524.4.4 Correlation with Rigidity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 534.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 545 Assessment of Cognitive Inflexibility Through Mode Detection and Delay . . 595.1 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 605.2 Model Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 615.3 Model Identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 645.4 Detecting Mode Changes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 655.4.1 Multiple Model Adaptive Estimation . . . . . . . . . . . . . . . . . . . . . . 655.4.2 Implementing the MMAE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 675.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 685.5.1 Switching Detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 705.5.2 Delay in Switching Detection . . . . . . . . . . . . . . . . . . . . . . . . . . . 725.5.3 Estimation in PD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73viiTable of Contents6 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 786.1 Summary of Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 796.2 Future Research Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 796.2.1 The Effect of Galvanic Vestibular Stimulation . . . . . . . . . . . . . . . . . 796.2.2 Reward Mechanisms in Forward Models . . . . . . . . . . . . . . . . . . . . . 826.2.3 Telemonitoring Parkinson’s Disease . . . . . . . . . . . . . . . . . . . . . . . 84Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86AppendixA Estimator dynamics in switching task . . . . . . . . . . . . . . . . . . . . . . . . . . 98viiiList of Tables2.1 Order of data pre-processing techniques for various experimental setups. . . . . . . . 133.1 Demographic characteristics of PD patients. . . . . . . . . . . . . . . . . . . . . . . . 233.2 Time lags (in seconds) removed by the correlation-based time-lag removal algorithmfor different PD subjects off medication. Note that the sampling time is 0.0167 seconds. 263.3 Time lags (in seconds) removed by the correlation-based time-lag removal algorithmfor different normal subjects. Note that the sampling time is 0.0167 seconds. . . . . 263.4 Accuracy of the models in different levels of noise for normal subjects. . . . . . . . . 283.5 Accuracy of the models in different levels of noise for PD subjects off-medication. . . 293.6 Correlation between model parameters and total tremor score. . . . . . . . . . . . . 304.1 Detailed characteristics of PD subjects. For rigidity, ‘LH’ denoted left hand, ‘RH’denotes right hand, ‘LL’ denotes left leg, ‘RL’ denotes right leg, and ‘N’ denotes neck. 374.2 The weights used for errors in normal LS and improved LS methods. The errors arein the format of [Q*(error before t seconds) + R*(error after t seconds)] . . . . . . . 474.3 Detailed dynamic parameters of Normal subjects. . . . . . . . . . . . . . . . . . . . . 494.4 Detailed dynamic parameters for the left arm of PD subjects off medication. Theentries shown as ζωn are the ones that no effective decay is used for their modelingand therefore the effective decay rate is just the multiplication of damping ratio andnatural frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 494.5 Detailed dynamic parameters for the right arm of PD subjects off medication. Theentries shown as ζωn are the ones that no effective decay is used for their modelingand therefore the effective decay rate is just the multiplication of damping ratio andnatural frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 504.6 Detailed dynamic parameters for the left arm of PD subjects on medication. Theentries shown as ζωn are the ones that no effective decay is used for their modelingand therefore the effective decay rate is just the multiplication of damping ratio andnatural frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50ixList of Tables4.7 Detailed dynamic parameters for the right arm of PD subjects on medication. Theentries shown as ζωn are the ones that no effective decay is used for their modelingand therefore the effective decay rate is just the multiplication of damping ratio andnatural frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 514.8 Static and dynamic model parameters for normal and PD subjects. . . . . . . . . . . 514.9 Static and dynamic differential model parameters for normal and PD subjects.Parameters used here are the difference of the ones for dominant hand and non-dominant hand. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 555.1 Dynamics for subjects shown in Figures 5.6, 5.7, and 5.8, with Qq = diag(αq, 0). . . 685.2 Switching detection and delays between ‘Better’ and ‘Worse’ modes in the ‘Normal-Better-Worse’ sequence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 725.3 Mean damping ratio and natural frequency in estimator dynamics for normal andPD subjects in ‘Better’ mode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 756.1 Static and dynamic model parameters for PD subjects off medication in ‘Better’ and‘Worse’ tasks. All of statistical significances are computed using paired ttest. . . . . 82A.1 Damping ratio and natural frequency in estimator dynamics for all normal and PDsubjects in ‘Better’ mode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99xList of Figures1.1 A schematic of the brain with dopamine pathways [133]. . . . . . . . . . . . . . . . . 21.2 A manual tracking task using a squeeze bulb. . . . . . . . . . . . . . . . . . . . . . . 41.3 Tracking performance of a PD subject before and after medication. . . . . . . . . . . 52.1 Pre-processing block diagram. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.2 Output signal for a typical PD subject performing the squeeze-bulb experiment be-fore (solid cyan) and after (solid blue) pre-processing. The input signal the subjectis asked to track is also shown in dashed red. . . . . . . . . . . . . . . . . . . . . . . 152.3 Inattentive subject conducting a tracking experiment. . . . . . . . . . . . . . . . . . 162.4 A low-pass filter with 8Hz cutoff frequency removes high-frequency noise from theraw output signal (cyan solid line), resulting in the pre-processed output signal (bluesolid line). The input signal (red dashed line) that the subject is trying to track isshown for completeness. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.5 The de-saturation filter (2.2) primarily effects large magnitude errors that occurduring actuator saturation. The raw output signal is shown in cyan solid line, thepre-processed output signal is shown in blue solid line and the input signal is shownin red dashed line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.6 Output signal for a typical PD subject performing the squeeze-bulb experiment (solidblue) that has been low-pass filtered as well as predicted output by the model (solidblack). The desired input signal the subject is asked to track is also included (dashedred). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.7 Output signal for a typical PD subject performing the squeeze-bulb experiment be-fore (solid cyan) and after pre-processing (low-pass filtering and correlation-basedtime-lag removal). The solid blue line represents the output shifted based on cor-relation analysis, and the solid black line represents the predicted output by themodel. For completeness, the desired input signal the subject is asked to track isalso included (dashed red). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21xiList of Figures3.1 Experimental setup. Subjects are instructed to keep the horizontal bar within avertically-scrolling, 0.75 Hz sinusoidal pathway by squeezing and releasing the bulb.The scrolling pathway was either unaltered, or perturbed with a Gaussian noise equalto 25% or 50% of the target signal amplitude. . . . . . . . . . . . . . . . . . . . . . . 243.2 Typical PD subject’s tracking performance in different tasks. The target signal isunaltered in the top picture, perturbed with a Gaussian noise equal to 25% and 50%of the target signal amplitude respectively in the middle and the bottom picture. Thedashed red line represents target signal, the solid cyan line represents squeeze bulbraw output signal, and the solid blue line is the preprocessed output. As expected,as noise amplitude increases, tracking performance deteriorates . . . . . . . . . . . . 253.3 The correlation of total tremor with damping ratio (p = 0.0009, robustfit). Notethat the logarithm of the parameters are used since the original distribution of themwas not Gaussian. The damping ratio of the normal subjects are also providedfor completeness. No tremor score is associated with normal subjects and they areprovided at a low tremor position for visualization purposes. . . . . . . . . . . . . . . 303.4 The correlation between total tremor and damping ratio over RMS error. The posi-tive correlation was statistically significant (p = 0.0822e-3, robustfit). Note that thelogarithm of the parameters are used since the original distribution of them was notGaussian. The damping ratio over RMS for the normal subjects are also providedfor completeness. No tremor score is associated with normal subjects and they areprovided at a low tremor position for visualization purposes. . . . . . . . . . . . . . . 313.5 Damping ratio for normal subjects and PD subjects (p = 0.6462, ANOVA). Notethat the logarithm of the parameters are used since the original distribution of themwas not Gaussian. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.6 Damping ratio over RMS error for normal subjects and PD subjects (p = 0.4913,ANOVA). Note that the logarithm of the parameters are used since the originaldistribution of them was not Gaussian. . . . . . . . . . . . . . . . . . . . . . . . . . . 334.1 Experimental setup. Subjects were instructed to relax to the best of their abilitieswithout interfering with the pendulum motion. Then, the arms of the subject werereleased to fall without informing the subjects. . . . . . . . . . . . . . . . . . . . . . 364.2 A pendulum. CG is the center of gravity. . . . . . . . . . . . . . . . . . . . . . . . . 384.3 Typical PD subject’s arm swing angular velocity. . . . . . . . . . . . . . . . . . . . . 41xiiList of Figures4.4 Typical Normal subject’s arm swing angular velocity. The dashed blue one is exper-imental, the dotted purple one is estimated by LS method, and the dotted black oneis estimated by LS after hypothesizing new damping ratio and re-run LS weightingthe beginning part more. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.5 Typical off-medication PD subject’s arm swing angular velocity. The dashed blueone is experimental, the dotted purple one is estimated by LS method, and the dottedblack one is estimated by LS after hypothesizing new damping ratio and re-run LSweighting the beginning part more. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 454.6 Typical on-medication PD subject’s arm swing angular velocity. The dashed blueone is experimental, the dotted purple one is estimated by LS method, and the dottedblack one is estimated by LS after hypothesizing new damping ratio and re-run LSweighting the beginning part more. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 464.7 An off-medication PD subject’s arm swing angular velocity with nonlinear behavior.The dashed blue one is experimental, the dotted purple one is estimated by weightedLS method for linear model, and the dotted black one is estimated by LS after usingthe model with effective decay rate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 484.8 Decay rate for normal subjects and PD subjects on and off medication. . . . . . . . 524.9 Effective decay rate for normal subjects and PD subjects on and off medication. . . 534.10 Differential decay rate for normal subjects and PD subjects on and off medication.Parameters used here are the difference of the ones for dominant hand and non-dominant hand. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 544.11 The decay rates for non-dominant arm. . . . . . . . . . . . . . . . . . . . . . . . . . 564.12 The decay rates for dominant arm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 574.13 The correlation of effective decay rate and rigidity. . . . . . . . . . . . . . . . . . . . 585.1 Experimental setup. The target trajectory is u(t) = sin(f1t) + sin(f2t). Users areinstructed to keep the cursor y(t) level with the target. The error u(t)−y(t) is scaledby 0.3 in ‘Better’ mode, by 2.0 in ‘Worse’ mode, and unchanged in ‘Normal’ mode. . 605.2 Typical normal subject’s tracking performance in different tasks. The solid blue linerepresents target position, and the dashed red line represents cursor position. Ascan be seen, when the manipulated error shown to them increases, they track thetarget worse. This figure will be replaced with a 3-mode hybrid system in Section 5.5. 62xiiiList of Figures5.3 Hybrid system H that represents the switched pursuit tracking task. Referenceinput trajectory u represents the target position, and output trajectory y representsthe cursor position. Modes ‘Better’, ‘Normal’, and ‘Worse’ refer to the attenuated,unchanged, or exaggerated error reflected in the manipulation of the cursor position. 635.4 Block diagram of model identification process. . . . . . . . . . . . . . . . . . . . . . . 645.5 Schematic of the MMAE algorithm. Kalman filters (KF) associated with each modegenerate a state estimate xˆq[k] that is used to create a residual for that mode.Residuals rq[k], Sq[k] are weighted and compared in the posterior probability evalu-ator (PPE), which generates a likelihood function Wq[k] for each mode. The mostlikely mode µˆ[k] is then selected. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 665.6 Typical normal subject. The solid blue line represents the target position, and thedashed red line represents the cursor position. The subject detects both switches,from ‘Normal’ to ‘Better’ as well as from ‘Better’ to ‘Worse’. The switching delayis very small in comparison to the time scale of the entire ‘Normal-Better-Worse’sequence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 695.7 Typical PD subject off-medication. The solid blue line represents the target position,and the dashed red line represents the cursor position. The subject is unable to detectboth switches, from ‘Normal’ to ‘Better’ as well as from ‘Better’ to ‘Worse’. We namethis type, a non-switching off medication PD subject. . . . . . . . . . . . . . . . . . 705.8 Typical PD subject on-medication. the solid blue line represents target position, andthe dashed red line represents cursor position. This subject does detect the switchat t = 50s from ‘Better’ to ‘Worse’. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 715.9 Typical PD subject off-medication. The solid blue line represents the target position,and the dashed red line represents the cursor position. The subject is able to detectthe switch from ‘Better’ to ‘Worse’. We name this type, a switching off medicationPD subject. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 735.10 The delay in switching detection for different groups in ‘Normal-Better-Worse’ block.Each time step is 0.03 seconds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 745.11 Damping ratio of estimator dynamics in the ‘Better’ task in the ‘Normal-Better-Worse’ sequence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 765.12 Natural frequency of estimator dynamics in the ‘Better’ task in the ‘Normal-Better-Worse’ sequence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 776.1 The change in damping ratio going from the GVS off state to the GVS on state forthe better task. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81xivList of Figures6.2 Pursuit tracking block diagram. ‘r’ is the reference trajectory. ‘u’ is the controlsignal generated by the brain. ‘e’ is the tracking error. ‘yactual’ is the actual cursorposition on the computer screen. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83xvAcknowledgmentsThis thesis would not have been possible without the help and support of many individuals. I wouldlike to express my sincere appreciation to Professors Meeko Oishi and Martin J. McKeown for theirsupervision and mentorship. Their enthusiasm to explore new ideas, expertise in the field, andtheir availability to provide me with feedback, have played a key role in conducting this research.I would not be where I am today if I had not their continuous support and giving me the insightof a big picture thinking. And for all that I am forever grateful.I have had wonderful collaborations with Nazanin Baradaran, Pouria Talebifard, Diana Kim,and Negin Niksirat—collaborations that I hope continue beyond our time at UBC. I would like tothank them for diligently pursuing our research paths and providing bright feedbacks. There is noway I would have possibly known how to navigate myself through this type of research had you notshared your insight.I have had the pleasure of being part of a great research group. Nikolai Matni, Neda Eskandari,Shahab Kaynama, and Pouyan Taghipour have all made my stay at UBC a great experience. Ithank Shahab in particular for being a great resource (both before and after his departure fromUBC). He has been both an amazing friend and a great brain to discuss ideas with. The ControlSystems Reading Group meetings were very supportive and inspiring: Special thanks to ProfessorsBhushan Gopaluni, Farrokh Sassani, Ryozo Nagamune, and the rest of the team, particularlyDaniel, Ehsan, Devin, Marius, Jin-Oh, and Klaske.Many scholars at UBC, and University of Tehran have had eternal influence on my understand-ing of controls and mathematics. I would like to particularly thank Professors Parviz JabehdarMaralani and Behzad Moshiri (University of Tehran) whose fascinating knowledge of the field,being so considerate and thoughtful about their students’ success, and consistent support of myacademic endeavors throughout the years have helped shape what I am today.I could not have done this work without the gracious help of the patients and their families whoenthusiastically participated in this research. You shall always be an inspiration to me.Last but not least, I thank my family—Dad, Mom, and Amir, Fatemeh—for their unconditionallove, motivating encouragements, and moral support. Thank you for teaching me how challengesin education, and life makes me well-rounded. The energy that I acquired on my trips backxviAcknowledgmentshome would never be forgotten because of your sincere care and devotion. My wonderful parents,to whom I earnestly and genuinely dedicate this thesis, have been an incredible source of love,support, and motivation. I am who I am today because of your selfless sacrifices and having faithin me throughout my whole life.This work has been financially supported by NSERC Discovery Grants #327387 and MichaelSmith Health Foundation for Research Team Start-up Grant, and the Institute for Computing,Information and Cognitive Systems (ICICS).xviiTo Mom and Dad–When my parents did not have my hand, they had my back.xviiiChapter 1Introduction1.1 Parkinson’s DiseaseParkinson’s disease (PD) is the second most common neuro-degenerative disorder (after Alzheimersdisease), affecting more than one million people in North America [73]. It is clinically characterizedby motor symptoms of rigidity, tremor, and bradykinesia (slowness of movement), and is sometimesaccompanied by cognitive impairments, such as deficits in working memory and executive functions[132]. The performance of sequential movements, coordination, and balance may also deteriorate[60]. In the late stages of the disease, mobility is severely impaired and in combination with severebradykinesia, can result in impaired balance and falls. Symptoms can be mitigated through L-dopa (a dopamine precursor) oral intake in early stages of the disease, and through more invasivemeasures (e.g., deep brain stimulation) in later stages of the disease [28].The motor cortex of the brain, which ultimately drives the muscles via the spinal cord andperipheral nerves, is modulated by two major systems: the basal ganglia and cerebellum [67, 68]. InPD, the unexplained degeneration of dopaminergic neurons of the basal ganglia nuclei – specifically,substantia nigra pars compacta (SNc) neurons – triggers functional changes affecting the wholebasal ganglia network [93, 126, 143] (Figure 1.1). Thus, SNc projections to striatum is impairedwhich leads to Parkinsonism symptoms [29, 70, 144]. The basal ganglia is most active when asubject performs an internally guided (IG) action that is recalled from memory [57]. On the otherhand, externally guided (EG) actions are normally associated with activity in the cerebellum wheresensorimotor computational processes are important [68]. Nevertheless, this association is not astrict dichotomy and the cerebellum can function during internally-cued movements and the basalganglia and supplementary motor area (SMA) are active in externally-cued movements [68].So while Parkinson’s disease is considered a classic basal ganglia disease, connectivity betweenthese two systems implies that both systems are affected. These two systems are critical in feed-forward and feedback processes that enable dextrous motor control, adaptation to changing envi-ronments, and effective executive function.11.1. Parkinson’s DiseaseFigure 1.1: A schematic of the brain with dopamine pathways [133].Motor tasks in Parkinson’s diseaseThe relationship between altered pathways in the brain and characteristics of motor performancein PD is not clearly understood, and may yield important insights into the disease.Consider, for example, the heterogeneity in progression of PD [44]. Some PD subjects havesignificant tremor that appears at an earlier stage of the disease, while others present with significantrigidity. Tremor-predominant PD subjects (PDT) tend to have a better prognosis and more slowdisease progression than those with a more Akinetic-rigid (PDAR) presentation [61]. Unlike theclinical signs of rigidity and bradykinesia, there is no known correlation between the severity ofPD tremor and the severity of dopamine depletion [36], but the mechanisms through which tremorcould be protective are unknown. A growing body of literature suggests compensatory changes atthe macroscopic, systems level. Structures outside the BG can potentially provide compensationfor the neurodegeneration seen in PD, as overactivity of the output regions of the BG has beenshown even in the absence of motor symptoms [15].In addition, the underlying mechanisms of rigidity are poorly understood, and no direct rela-tionship exists between dopamine deficiency and rigidity [110]. Researchers have suggested thatimproper commands from one or several descending spinal pathways defect the short reflex path-ways at the spinal level [34]. Others explain its cause by increased muscle tone, and excessive andcontinuous contraction of muscles. Clinical observations also suggest that reinforcing movements(e.g., voluntary actions of the contralateral limb), known as the Froment’s manoeuvre, increaserigidity [60]. Measures of motor performance that describe the underlying dynamics, as opposed to21.1. Parkinson’s Diseasemere gross measures (e.g., maximum speed, mean speed, distance traveled), could be particularlyuseful in elucidating various theories of rigidity, as well as its role in PD.Hence motor tasks which are designed to elucidate different pathways in the brain may provideinsight into the potential existence of disease subtypes, as well as a possible measure of diseaseseverity. Further, rigorous and detailed assessment of motor tasks may also provide a means ofcharacterizing phenomenon such as rigidity that are often assessed in clinical settings that canobscure the subject’s true motor capabilities.Lastly, we note that it may be possible to identify cognitive effects in carefully designed motortasks. The Wisconsin Card Sorting Task (WCST) [5, 103] is a common measure of cognitiveinflexibility. In the WCST, the player receives a reward for correctly sorting cards according to acriterion that can be learned as the game is played. However, at some point during the task, therules for sorting change (without enunciation to the subject), and the player must learn the newrules in order to maximize their reward. PD subjects often have difficulty re-learning the rules,and continue to play by the old rules even though they may incur a loss for doing so [103]. Inthe context of a motor task, cognitive inflexibility may be evaluated in measuring the delay in asubject’s response to a sudden change in the underlying dynamics of the task.Assessing Parkinson’s diseaseAssessment and inference of mechanisms in the brain is complicated by the fact that assessmentof disease severity in Parkinson’s disease is not simple and straightforward. Many qualitativeand quantitative assessments have been conducted to evaluate disease severity in PD, spanningcognitive [74, 79, 98, 118] and motor [19] domains. However, qualitative approaches such as neu-roanatomy analysis [118] or attention control [23, 149] are not capable of describing disease severityor mechanisms in the brain to a decent extent.Given the inherent heterogeneity of the disease, assessments of motor performance that arehighly sensitive to individual subjects are important. How motor tasks are accomplished mayallow insight into specific compensatory mechanisms in the brain. The development of a suite ofmotor tasks, well-suited in particular to elucidate some of the faulty feedback mechanisms as wellas compensatory mechanisms in PD, has the potential to improve the assessment of disease stateand to provide insight into disease progression and subtypes. Further, the low-tech nature of suchtasks could enable remote health assessment in rural communities. For some subjects, performingmotor tasks at home could alleviate performance anxiety that may affect motor performance in aclinical setting. Further, motor tasks that are small-scale, manual tracking tasks can be performedconcurrently with fMRI or EEG in clinical settings, allowing for direct correlations between observedmotor performance and brain data.31.2. Problem Formulation: Manual Tracking TasksFigure 1.2: A manual tracking task using a squeeze bulb.1.2 Problem Formulation: Manual Tracking TasksManual tracking tasks, in which subjects are asked to track a target or moving object with a hand,or a cursor under manual control (Figure 1.2), are an increasingly popular means to assess motorperformance, as they can be executed during brain mapping procedures (e.g. functional MagneticResonance Imaging (fMRI), electroencephalography (EEG), positron emission tomography (PET))in which excessive head motion must be prevented. Pursuit tracking tasks involve multiple brainfunctions, including planning, execution, motor control, sensory feedback, and saccadic and smoothpursuit eye movements, so it is difficult to discern with certainty the particular structures of thebrain involved in a given task. However, two particular regions of the brain, the basal ganglia andthe cerebellum, can be emphasized by introducing visual errors in the cursor position [79], tactilechanges in the input device [63], and by using target trajectories with sudden and unexpectedchanges [7, 32].While the exact mechanism by which these two brain regions influence motor output has notbeen fully explained [68], the basal ganglia seems to emphasize the importance of movement se-lection [10], while the cerebellum is considered essential for anticipating sensations that will beexperienced by motor actions (indeed, one model of the cerebellum is as a Kalman filter [88, 105]).However, the presumed dichotomy is by no means clear cut.Common measures of performance such as root mean square (RMS) error, absolute mean error,maximum error, lag time, peak velocity, and mean velocity [6, 48, 63, 65, 79, 141], convey agross measure of tracking performance, but do not convey subtleties of how that performance isachieved [74, 101]. Such distinctions may be critical in exploring potentially different strategiesthat Parkinson’s disease subjects employ to maintain overall tracking performance, and may havediagnostic and discriminatory importance.We use linear dynamical systems to model manual tracking tasks [8, 95, 131] and to charac-terize motor performance. Oftentimes these tasks are sinusoidal nature, however, step responseparameters are useful in capturing some of the more qualitative aspects of tracking, such as thewell-known, clinically observed overshoot in subjects on L-dopa (and undershoot in subjects with-41.2. Problem Formulation: Manual Tracking TasksFigure 1.3: Tracking performance of a PD subject before and after medication.out L-dopa) (Figure 1.3). We have found (as have other researchers) that second-order modelscapture important distinguishing features of the tracking data. Further, model parameters that aredependent on pole locations can be compared across subjects (and across groups).We presume a black-box formulation with input that is the trajectory the subject should track,and output that is the actual position of the cursor under subject control. We presume thatthe black-box represents either a linear system or a linear switched system, depending on theexperimental setup.u(t)−→ x˙(t) = Aq(t)x(t) +Bq(t)u(t)y(t) = Cq(t)x(t) +Dq(t)u(t)−→ y(t) (1.1)We solve a standard system identification problem in Matlab using established time-domain meth-ods (e.g., ARX, ARMA, ARMAX) available in the System Identification Toolbox [80, 81] to obtainthe system matrices (Aq, Bq, Cq, Dq) for a given mode q ∈ Q, the set of discrete modes. These meth-ods are based on minimizing the error between the predicted outputs and measured outputs usingpast values of input/output measurement. Pre-processing of the data is required to remove high-frequency noise, low-frequency bias, and saturation effects, to average multiple trials of the sameexperiment under the same conditions, to resample data as appropriate to the system identificationmethods used, and to remove any data incongruities (e.g., if the subject was not paying atten-tion or did not complete the experiment). Details on the pre-processing used for each particularexperimental setup is described in Chapter 2.Data are gathered at the Pacific Parkinson’s Research Centre / University of British Columbia51.3. MethodsMovement Disorders Clinic, where more than 1100 Parkinson’s disease subjects are seen annually,under the supervision of neurologist Prof. Martin J. McKeown. Data are gathered according tostandard protocols and with approval from University of British Columbia’s Ethics Board. In allexperimental setups, subjects practice the task during a training session until they are familiarwith the task requirements (that is, until their tracking errors converge to a constant value). Weonly consider data gathered after the training period is completed. The experiments are completedby subjects with clinically definite mild to moderate Parkinson’s Disease (measured clinically asHoehn & Yahr stage 1-3), both off and on medication, as well as by “normal” age-matched subjectswithout Parkinson’s disease or other confounding ailments.Among different types of models (e.g., linear, nonlinear, different orders, neural networks [66],etc) linear time-invariant (LTI) dynamical systems have proven by researchers to be a valid ab-straction of the human response in manual tracking, and can be correlated with brain data instudies involving control and PD subjects [2, 8]. LTI model parameters can provide insights intobrain function in disease state, by quantifying dynamic motor performance. In addition, stabilityanalysis is straightforward in linear models by just assessing the poles of the transfer function [37].On the other hand, other types of models, particularly nonlinear models need much more data totrain and lead to higher computational cost in many analyses. They are hard to be compared witheach other in terms of speed, stability, etc. Lastly, we note that while many real processes in natureare nonlinear, they can often be considered to be approximately linear over a small operating range(a small frequency domain in our experiments).It should be noted that the linear models that we build only describe the behavior of the systemwithin the frequency range of the input signals (less than 1 Hz). Also, if the experiment (in thiscase the input-output tracking data) is inherently nonlinear, the linear model cannot capture all ofthe dynamics. This limitation has been discussed in manual tracking of sinusoidal tasks by otherresearchers [2, 48] and they have validated the use of linear models to be sufficient. Linear modelsare sensitive to outliers. This limitation is mitigated by numerous amount of data gathered foreach experiment. Linear models, in general, tend to have higher mean square errors compared tononlinear models [3]. However, the simplicity, scalability, and comparability of linear models haveled to their extensive use wherever approximation could be tolerated.1.3 MethodsWe use Matlab’s system identification toolbox [80] to create LTI dynamic models. We use ARMAX[50, 148] and continuous time state space models as the main tools. For the squeeze-bulb and61.3. Methodspendulum experiments, continuous time models of the following form are identified.x˙(t) = Ax(t) +Bu(t)y(t) = Cx(t) +Du(t) (1.2)For the glass rod experiment, discrete-time models are identified via ARMAX [50, 148] whichoccurs by solvingy[k] = Hu[k] + e[k] (1.3)for the linear operator H that relates the input data u to the output data y. The error e[k] iszero-mean white noise, presuming both the input and the output have also been detrended, suchthat E(u) = 0 and E(y) = 0, where E(·) indicates expected value. This is essentially a “black-box” method [26, 82], meaning that no structure in H (aside from linearity and time-invariance) ispresumed. The model H is expressed as the ARMAX modely[k] + a1y[k − 1] + a0y[k − 2] = b2u[k] + b1u[k − 1] + b0u[k − 2] (1.4)Hence the resulting discrete-time model (shown below in canonical form)x[k + 1] =[−a1 1−a0 0]x[k] +[b1 − a1b2b2 − a0b2]u[k]y[k] =[1 0]x[k] + b2u[k](1.5)is identified. Let us denote these matrices by:A =[−a1 1−a0 0], B =[b1 − a1b2b2 − a0b2], C =[1 0], D = b2 (1.6)Note that the discrete identified matrices (A,B,C,D) are related to the continuous-time form(via a zero-order hold with sampling interval T ), that we mentioned in equation (1.2), viaA = eAT , B =∫ T0eAτBdτ, C = C, D = D (1.7)Solution to the ARX / ARMA / ARMAX identification problem (for an LTI system as describedin equation (1.1) for a given mode) is computationally efficient, and algorithms have been imple-mented in the System Identification Toolbox in Matlab [80]. We then analyze the eigenvalues of Amatrix to obtain damping ratio, natural frequency, and other model parameters in all experiments.71.3. MethodsThe resulting models are validated via standard measures: a) error coherence, b) the ability ofthe model to match data it was created from (e.g., Akaike Information Criterion or Final PredictionError), and c) the ability of the model to predict the output of an independent data set (e.g., thefitness score). Since we have already established the desired model order, AIC / FPE are lessrelevant for our data sets and overfitting less of a concern. Hence our method is to check the errorcoherence for significant outliers (to assure validity of the model), and to then compute the fitnessscore.Fitness score =1−√√√√ 1NN∑i=1(y[i]− yˆ[i])2√√√√ 1NN∑i=1(y[i]− y¯)2× 100 (1.8)Essentially, system identification methods try to account for as much of the relationship as possiblebetween the input and the output in the operator H, and what cannot be accounted for is lumpedinto the white noise term e. For systems in which second order is simply not adequate to capturethe input-output behavior, the noise term will be very large. As it is possible to create modelswhich have high fitness score but do not intuitively match the data well, we also follow this analysiswith visual checks, and correct models which are poor matches. For example, we can initiate theARMAX identification from model parameters are closer to a visual match.In general, tradeoffs will exist between accuracy and model order, meaning that with increasedmodel order, higher model accuracy will result. This could result in overfitting, in which theresulting model is very good at predicting the output from the data set it was generated in, butpoor at predicting the output from an independent data set (that it was not generated from).However, because we are limiting our models to second order, overfitting is not a concern. (Wehave evaluated higher order models for completeness, but our analysis is similar to previous work,in that we also found that second order models handle this tradeoff in the best manner.)In addition, because we have inputs of limited spectrum in some tasks, we note the limitationsof the resulting model for these data sets – the model is accurate only for other input signals inthat same spectrum. This limitation has also been discussed in similar analyses of manual trackingof sinusoidal tasks ([2, 48]).Other time-domain approaches to solving LTI models are grey-box methods [80, 128], andsubspace-based methods [43, 139] subspace identification methods (e.g., for direct identificationof a state-space model). Frequency domain methods estimate the discrete-time transfer functionG(z) = b2z2+b1z+b0z2+a1z+z0 via iterative estimation techniques. We primarily use the ARMAX approachbecause the error coherence seems to be more consistent than for subspace or transfer function81.4. Related Workmethods (which are more likely to violate bounds for model validity for our particular data sets).However, we do make use of grey-box system identification in one experimental setup, for whicha specific coordinate system is required for the state-space model in order to analyze switchingbetween models (Chapter5).1.4 Related WorkLTI models of manual trackingResearchers explored linear dynamical system models of manual tracking tasks as early as the late1970s [48], with application to human-in-the-loop automation [31, 98, 129]. Abdel-Malek et al. [1, 2]used ARMA models to characterize tracking in response to sinusoidal input trajectories in operatordynamics as well as to quantify tracking performance with Parkinson’s disease. Interestingly, thesestudies all found second-order LTI models sufficient for capturing human response to sinusoidaltracking tasks. Another small study (5 subjects) found that with a more complex tracking task,involving 5 harmonically related frequencies, the resultant LTI models varied with different thetracking tasks [31]. In all of these scenarios, the limited spectrum of the tracking trajectory affectsthe spectrum over which the identified model is valid [81], an inherent disadvantage of sinusoidaltracking. All of these studies all used a black-box approach that ignores any theoretical constructionof motor control processes.In [98], a pursuit tracking task is done by the subjects in which the input is the summationof two sinusoidal signals with different frequencies. This would generate a relatively unpredictablewaveform for the subjects to track and hence allows the causality assumption for the models.Then, assuming the human motor system as a black-box, they created transfer function modelsand analyzed their spectral properties as well as time series properties. They found that time-series analysis (using eigenvalues) outperforms Fourier analysis (in frequency domain) in describinghuman tracking behaviour.Abdelmalek and his colleague in [2] chose sum of five sinusoidal signals as well as random ternaryinputs of various bandwidths for pursuit tracking task. They used time-series models which first wasintroduced by Shinners in [125]. Shinners had developed autoregressive moving average (ARMA)models for human motor behaviour in compensatory tracking experiments. Among different cri-teria for choosing model order (e.g., final prediction error, Akaike’s information criterion, residualvariance, entropy) they chose entropy as it does not need model creation prior to calculation. Cre-ating models up to 5th order ones using least square method, they found that second order modelsare capable to characterize those pursuit tracking task to a decent extent.91.4. Related WorkState estimation and forward modelsOther researchers have focused on computational modeling of specific motor control processes [133,163] in the brain and in the spine.The cerebellum as well as parts of the spine are hypothesizedto create internal models [146], such as a Kalman filter [88, 105], that is used as feedforward inmotor control processes. Feedback occurs when sensory information is compared with predictionsobtained from the internal model [119], and is thought to occur in the basal ganglia [68], althoughsuch a strict dichotomy is likely an oversimplification of actual brain processes [120]. Experimentshave been designed to evaluate the relative role of feedforward and feedback mechanisms, by, forexample, changing the visual field to create mismatches between sensory information and internalmodel predictions [74, 79]. Researchers have hypothesized optimal control [78, 135] as a frameworkfor how movements are generated, and reinforcement learning [33, 59] as a framework for how newmovements are acquired. Much of the work in these areas focuses on specific tasks (e.g., reaching ortapping) and is based on measurements of EMG, saccadal eye movements, and simple assessmentsof motor performance (e.g., time to reach).State estimation via a Kalman filter has been proposed as a theoretical framework to explorefeedforward and feedback processes during normal motor performance, and may therefore providea more direct measure of motor dysfunction in PD [71, 86, 105, 118, 134, 146].For example, Miall and his colleagues investigated the role of state estimation by the cerebellum[86]. They hypothesized that the cerebellum combines afferent sensory inputs to the brain andefferent motor commands from it. They had already figured that transcranial magnetic stimulation(TMS) over the cerebellum will lead to faster responses to visual stimuli in healthy subjects [89].Their results are suggestive of the idea that cerebellum plays an important role in updating stateestimates.In [9], a visuomotor pointing task is studied to investigate hand movement performance inresponse to cumulative small changes in the visual feedback. Subjects were asked to point toa visual target in discrete blocks via a cursor whose position was perturbed with noise in threedifferent levels. Observer analysis for this task explores the ability of subjects to suppress the noiseand compares that to that of Kalman filters. Their results suggest that the subjects exponentiallyweighted recent errors with a weight that was dependent to the level of the noise. These weightswere observed in both Kalman filters and a modified delta learning rule that appeared to modelhuman behavior with a relatively good fit. Furthermore, they hypothesized that subjects performbetter if the perturbation is introduced as small incremental steps as opposed to a big change. Thisfinding is consistent with previous studies [57] and is well suited with exponential weight of errors.Many studies indicate that discrepancies between actual and expected sensory consequences of101.4. Related Workmotor actions may be the incentive for motor learning and eventually enhancing motor performance[105, 118, 146]. However, the specific brain structures engaged in developing and updating offeedforward and feedback processes are still investigated by many researchers [25, 88, 118]. In [118],Shadmehr and his colleagues try to describe brain computations to form a motor command. Theyclaim that imagining the movement is interpreted as forming the cost of that task. Minimizationof this cost is done through a forward model that potentially is located in the cerebellum. On thesame note, basal ganglia is responsible for learning costs and rewards associated with certain tasksand thus uses “optimal control” as its main tool. Parietal cortex is engaged in state estimation andprovides sensory feedback for the cerebellum to predict sensory consequences of motor commands.In [86], the effect of varying TMS in a reach-to-target task is tested. Miall and his colleaguesshowed the error (the difference of the final finger position and the target) is not correlated withthe TMS pulse-train duration. Furthermore, the estimated time of the movement is likely to be incorrelation with physiological processes in the cerebellum. In other words, they draw upon the ideathat the cerebellum builds a “forward model” that predicts motor commands according to sensoryinputs [88]. Their results support the idea of the powerful role of the cerebellum in updating stateestimates and its dysfunction will lead to motor errors.Deep brain stimulation (DBS)We also note the use of control theory and system identification to model and control neuronal pat-terns in the β-band (8-40 Hz) related to pathological oscillations that are associated with Parkin-son’s disease.Within the controls community, work has been done on modeling of low-level neuronal patterns[27] associated with Parkinsonism symptoms. Other work has focused on low-level neurologicalprocesses in the basal ganglia [113], a part of the brain which is directly affected by the loss ofdopamine in Parkinson’s disease. Researchers have begun to study the effects of deep brain stim-ulation [114], in which electric pulses are applied in an open-loop fashion. Deep brain stimulationappears to alter pathological oscillation in the beta band (8-30 Hz) associated with Parkinson’sdisease [111]. Some work has been done on how to use feedback to modify the application of deepbrain stimulation [30, 46, 58, 107]. Models of how the stimulation affects the brain are difficult toformulate [13], and accurate measurements of how brain waves are altered through the stimulationare difficult to obtain [12]. Other work has focused on the motor effects of applying deep brainstimulation [99, 109, 122].111.5. Contributions of the Thesis1.5 Contributions of the ThesisThis thesis focuses on the use of control theoretic measures to help elucidate mechanisms in thebrain in Parkinson’s disease. The main contributions are:• Selection of pre-processing and system identification techniques for three experimental setups formanual tracking experiments• Control theoretic framework for tremor as a compensatory mechanism that is advantageous insome tracking tasks• Sensor-based assessment of rigidity via decay rate• Assessment of cognitive inflexibility through mode detection and delay (via hybrid estimation)The above contributions have resulted in three journal publications and two refereed conferencepublication, as described in the Preface.1.6 Organization of the ThesisThe remainder of the thesis is organized as follows. Chapter 2 provides necessary pre-processingsteps to prepare input-output signals for system identification. In Chapter 3, existence of com-pensatory mechanisms in the brain is investigated. Specifically, the relationship between tremorand motor performance is analyzed to assess the effect of tremor as a compensatory mechanism.An objective measure for rigidity through a simple pendulum test is provided in Chapter 4. InChapter 5, we use a novel test to evaluate both motor performance and cognitive inflexibility –manual pursuit tracking with sudden changes in the tracking dynamics. Chapter 6 summarizes thethesis repeating its main contributions, and provides directions for future research.Finally, supplementary materials are provided in Appendix A, which contains the model pa-rameters presented in Chapter 5.12Chapter 2Pre- and Post-Processing TechniquesA key step in modeling manual tracking data, especially data from subjects with disease, is ap-propriate pre-processing of the input and output signals: biological data may be corrupted bymovement artifacts, sensor noise, actuator nonlinearities, data gaps, and human error, for exam-ple. Further, given that our models are solely second-order, we downsample the original data toensure that we can capture phenomenon at frequencies of interest. We have found that for theexperimental setups that we focus on, careful inspection of the data as well as systematic pre-processing is vital to identify models that have high accuracy and good predictive ability. As inmany biomedical applications [47, 51, 62, 145], pre-processing techniques must be tailored to theparticular experimental setup and data gathering process.We preprocess the data by employing filters for detrending, downsampling, removal of satura-tion effects, and removal of time-lags [124], depending on the experimental setup and the desiredidentification problem. Some of these filters are applied to the output, while others are applied toboth the input and the output. Because some of these pre-processing steps are nonlinear transfor-mations of the original data, the ordering of the various filters does matter (Figure 2.1). Table 2indicates the order of pre-processing steps for various experimental setups.The main goal of the pre-processing is to remove those features which are artifacts, or whichare problematic for LTI system identification algorithms (or cause significant problems in resultingSqueeze Bulb Pendulum Glass RodData pre-processing Experiment Experiment Experiment(Ch. 3) (Ch. 4) (Ch. 5)Low-pass filter 5 2 2De-saturation filter – – 3Downsampling 4 – 4Correlation-based shift 3 – –Interpolation 1 1 –Averaging Trials 2 – 1Table 2.1: Order of data pre-processing techniques for various experimental setups.13Chapter 2. Pre- and Post-Processing TechniquesRawDataSignalConditioningTime-lagRemoval ResamplingPre-processedDataFigure 2.1: Pre-processing block diagram.model accuracies or predictive capabilities) and are not indicative of important movement charac-teristics. These pre-processing techniques exploit the fact that we have a reasonable idea a prioriof what the resulting model should look like – namely, that it is likely a LTI second-order system.Figure 2.2 shows the effect of pre-processing on glass-rod experimental data, which significantlyimproves the resulting model without altering important features of the input-output relationship(such as overshoot).However, there are inevitably some phenomenon (e.g., nonlinearities inherent to the trackingtask) which are not directly addressed by these pre-processing techniques. While in most cases alinear model appears to be a reasonable approximation to the input-output relationship, in someinstances, neither the raw nor filtered input-output data is well represented by an LTI model.Nonlinear system identification and comparative analysis of nonlinear systems is beyond the scopeof this work, but may be an interesting avenue for future investigations.Careful screening of the data is required to assure that the pre-processing techniques are appliedto valid data. Cases such as that shown in Figure 2.3 (in which the subject likely fell asleep briefly)were rare in the data sets examined, and hence have been treated as outliers. In cases in which theexperiment was conducted multiple times, and the other experiments showed valid data, we simplyignore the clearly invalid data. While the uncertainty associated with resulting model parametersincreases with fewer data, such an approach still allows us to make full use of the valid data.After identifying the LTI models and computing LTI model parameters (damping ratio, naturalfrequency, decay rate, settling time, rise time, peak time, overshoot) we post-process the resultingparameters to aid in analysis across groups. We use log transformations as well as linear regressions,depending on the particular experimental setup and the similarity of the parameter distributionsto a Gaussian distribution.142.1. Data Pre-Processing0 5 10 15 20 25−150−100−50050100150Time [seconds]Target and cursor position [mm] Pre−processed outputraw outputdesired inputFigure 2.2: Output signal for a typical PD subject performing the squeeze-bulb experiment before(solid cyan) and after (solid blue) pre-processing. The input signal the subject is askedto track is also shown in dashed red.2.1 Data Pre-Processing2.1.1 Low-Pass FilterA variety of high frequency phenomena appear in the output signal. As with may data acquisitionsystems, transducers may exhibit to high-frequency noise.To attenuate unwanted power in high frequencies, we employ a butterworth low-pass filter[69]. The cutoff frequency is chosen to be 8 Hz for data sets from squeeze bulb (Chapters 3),pendulum (Chapter 4), and glass rod (Chapter 5) experimental setups. This filter has a fairly sharpattenuation without excessive distortion in the magnitude of nearby frequencies. We evaluatedbutterworth, Bessel, and elliptic low-pass filters, and found the fourth order butterworth filterto be the most effective on all data sets (Figure 2.4). Equation (2.1) shows the filter frequencyresponse. Parkinson’s tremor usually occurs at 2–6 Hz, and the input signals occur at around0.25–1 Hz.G(ω) =√11 + ω8 (2.1)152.1. Data Pre-Processing0 5 10 15 20 25 30406080100120140160180200220Target and cursor position [mm]Time [seconds]Figure 2.3: Inattentive subject conducting a tracking experiment.2.1.2 De-Saturation FilterIn the squeeze bulb and joystick experimental setups, saturation would occur infrequently whenthe subjects’ input signal reached the physical constraints of the device (for example, pushing thejoystick position to its upper or lower extremum, as in Figure 2.5). Such nonlinearity cannot becaptured through an LTI model with white noise. Hence we created a de-saturating filter thatreduces large magnitude errors that are created during periods of saturation, while having minimaleffects on smaller magnitude errors (since dde tan−1(e)|e=0 = 1).e[k] = y[k]− u[k]esat[k] = tan−1(e[k])ynew[k] = u[k] + esat[k](2.2)2.1.3 InterpolationIn rare instances, small segments of data may not be recorded (e.g., sensor malfunction). Thesegaps are relatively small (about 4–5 time samples) compared to the length of recorded signals.Hence, rather than segmenting the data around the data gaps, we use cubic interpolation to fillin reasonable values for all sampling points within the data gap. A third order polynomial was162.1. Data Pre-Processing0 5 10 15 20−80−60−40−20020406080Target and cursor position [mm]Time [seconds] desired inputraw outputfiltered outputFigure 2.4: A low-pass filter with 8Hz cutoff frequency removes high-frequency noise from the rawoutput signal (cyan solid line), resulting in the pre-processed output signal (blue solidline). The input signal (red dashed line) that the subject is trying to track is shown forcompleteness.selected for interpolation to preserve data smoothness.2.1.4 DownsamplingData was gathered at a very high frequency (128 Hz) in pendulum, and 60 Hz in glass rod, andsqueeze bulb experimental setups. However, with merely second-order models, we are unable toaccurately capture behavior at a significantly lower frequency than the sampling frequency. (Withhigher order models, accurate model identification would not be a problem, but with second-ordermodels, only the state from at most two time-steps prior can recursively affect the current state). Weemploy the standard rule of thumb [81] and select a sampling rate 6-10 times above the bandwidthof the system (which we approximate by the highest input frequency). No downsampling wasdone for the pendulum experiment. For the squeeze-bulb, and glass rod experimental setups, thedownsampled sampling frequency was 15, and 30 Hz, respectively.172.1. Data Pre-Processing0 5 10 15 20 25 30−400−300−200−1000100200300400Target and cursor position [mm]Time [seconds] desired inputraw ouputde−saturated outputFigure 2.5: The de-saturation filter (2.2) primarily effects large magnitude errors that occur duringactuator saturation. The raw output signal is shown in cyan solid line, the pre-processedoutput signal is shown in blue solid line and the input signal is shown in red dashedline.2.1.5 Time-Lag RemovalIn some data sets, visual inspection of the input-output data suggests a time-varying lag timebetween the input and the output. While small constant lags can be accommodated through azero in the LTI model, large or time-varying lags cannot and more rigorous techniques to handlenon-stationary data is needed.We deploy a time lag removal technique based on correlation analysis of input and output. Inour experiments, the time-lag is not remarkably varying with time but a large delay sometimesexists between the input and output signals. Hence, shifting the output (or the input) signal wouldefficiently help model identification.Correlation-based time-lag removalCross-correlation is a measure of similarity of two waveforms as a function of a time-lag applied toone of them [96]. In this technique, we form the cross correlation of the input and output signals182.1. Data Pre-Processingaccording to the following equation:R(τ) =∫ ∞−∞u(t)y(t+ τ) (2.3)Assuming that the input and output are similar in a tracking experiment, the maximum (orminimum if the signals are negatively correlated) of the cross-correlation function indicates thepoint in time where the signals are best aligned and hence is a good candidate to be the delay ofthe system (τ).τdelay = argmaxτ R(τ) (2.4)Depending on whether the highest correlation provided in equation (2.4) occurs at a positiveor negative time, the input or output signal is shifted along the temporal axis. This method is notcomputationally expensive, which makes batch processing of signals easier.We compared the model fidelity of correlation-based shifted data to that of the original data(without shifting) on a sample subject in the squeeze bulb experiment. Figure 2.6 shows thepredicted output by the model built from non-shifted data while Figure 2.7 shows the predictedoutput by the model built from correlation-based shifted data. The accuracy of the model built fromnon-shifted data is 65.35%, while the accuracy of the model built from correlation-based shifteddata is 72.65%. The better predicted data is also visually observed at peaks and valleys. Note thatthe difference of these two models prediction capability is more evident in more complicated tasks,e.g., squeeze bulb task when the desired input is perturbed with noise.We note that there are other methods for estimating the delay such as methods based on Fouriertransform. The delay shows its effect on the phase of the signal as ωτ in each frequency ω, whilenot changing the amplitude of the signal in the frequency domain. This method did not turn outto be very efficient as it led to different delays in different frequencies.Lastly, careful inspection of different experiments revealed that among different experimentsetups, the squeeze bulb experiment is mainly concerned with severe time-lags. There is no inputfor the pendulum experiment as the arm is swinging passively due to gravitational force. So, thetime delay between input and output is not defined for it. The glass rod experiment which isconducted with joystick is not suffering from significant delays between input and output.192.2. Post-Processing0 5 10 15 20 25−150−100−50050100150Time [seconds]Target and cursor position [mm] output−−not shiftedoutput−−predicted by modelinput−−desiredFigure 2.6: Output signal for a typical PD subject performing the squeeze-bulb experiment (solidblue) that has been low-pass filtered as well as predicted output by the model (solidblack). The desired input signal the subject is asked to track is also included (dashedred).2.2 Post-Processing2.2.1 Non-Gaussian Data for Statistical AnalysisOnce appropriate input-output data was prepared for system identification, models are built. Themodel parameters are used for statistical analysis and potential correlation with clinical measures.However, due to heterogeneous and distributed nature of biomedical data these parameters mightrequire to be transformed prior to analysis. Most of statistical tests rely on the assumption thatthe data has a Gaussian distribution which rarely is the case for biomedical data. Therefore, trans-formations could be used to make the distribution Gaussian. The most common transformation islog-ratio transformation that is used in this thesis.202.2. Post-Processing0 5 10 15 20 25−150−100−50050100150Time [seconds]Target and cursor position [mm] output by Correlation−based shiftoutput−−predicted by modelinput−−desiredFigure 2.7: Output signal for a typical PD subject performing the squeeze-bulb experiment before(solid cyan) and after pre-processing (low-pass filtering and correlation-based time-lagremoval). The solid blue line represents the output shifted based on correlation analysis,and the solid black line represents the predicted output by the model. For completeness,the desired input signal the subject is asked to track is also included (dashed red).21Chapter 3Insights Into CompensatoryMechanisms of TremorA strange symptom of Parkinson’s disease, is the considerable heterogeneity in its progression[44]. There is no correlation between the severity of PD tremor and the severity of dopaminedepletion [36]. But, PET scans have shown that resting activity of rostral, medial, and intermediatecerebellum is increased in PD tremor [39]. Nonetheless, the relationship of tremor and motorperformance is incompletely known.Assessment and biological interpretation of compensatory mechanisms in the brain is a com-plicated and delicate task. In spite of many studies and empirical observations, the structural andfunctional effects of PD tremor has been difficult to elucidate [36]. Although PD is considered abasal ganglia disease, we want to explore structures outside the BG that can potentially providecompensation for the neurodegeneration seen in this disease. To this end, in this chapter we willuse linear dynamical systems to investigate the relationship of motor performance and tremor intracking tasks.Ten PD subjects (off medication) and ten age-matched healthy subjects performed a manualpursuit tracking task involving squeezing a rubber bulb to guide a horizontal bar displayed on acomputer screen through a vertically scrolling sinusoidal pathway. While subjects were instructedto maintain a sinusoidal force output, the sinusoidal pathway was corrupted with different noiselevels: shown either unaltered or perturbed, with a Gaussian noise equal to 25% or 50% of thetarget signal amplitude. After pre-processing the input-output signals, second order continuous-time models are identified to explore potential correlation of dynamical parameters with clinicalmeasures.3.1 Experimental SetupTen PD subjects (off L-dopa medication) with clinically diagnosed, mild to moderate PD (Hoehnand Yahr stage 2–3) [54], were recruited for this study at the Pacific Parkinson’s Research Centre(PPRC)/Movement Disorders Clinic [95]. All subjects provided written informed consent, and all223.1. Experimental SetupTable 3.1: Demographic characteristics of PD patients.Demographic Characteristics PDGender 4 male, 6 femaleHandedness 8 RH, 2LHH&Y Stage 2− 3Mean symptom duration 5.8± 3 yearsMean age 66± 8 yearsAverage daily dose of L-dopa 685± 231 mgMean UPDRS 26± 8research was approved by the University of British Columbia Clinical Research Ethics Board. Allsubjects’ clinical Unified Parkinson’s Disease Rating Scale (UPDRS) score, including tremor scores,were assessed by a physician trained in the UPDRS (see Table 3.1 for demographic details).The manual tracking data were collected during a functional Magnetic Resonance Image (fMRI)study, where subjects were instructed to use their right hand to squeeze a rubber bulb in responseto visual stimuli displayed on the computer screen (See Figure 3.1 [95]). The squeeze bulb wasconnected to a pressure transducer (Honeywell, Inc., Plymouth, MN, model PPT0100AWN2VA)outside the scanner room via low-compliance water-filled tubing. Using the squeeze bulb, patientswere required to control the width of a horizontal bar on the screen to keep it within a sinusoidal,vertically scrolling pathway. Increasing pressure on the bulb led to increasing width in the bar andvice versa. The scrolling pathway was shown either unaltered, or perturbed with a Gaussian noiseequal to 25% or 50% of the pathway width amplitude. Subjects were instructed to maintain aspure a sinusoidal output as possible in all conditions, including when the path had added noise.Maximum Voluntary Contraction (MVC) for each subject was assessed by asking them to squeezethe bulb with their maximum force at the beginning of a 30-min training session while the pressurewas measured. First, subjects were asked to squeeze 5–15% of their MVC, so that they couldperform the task for longer periods. Then, a sequence of three separate tracking tasks (100 secondseach) was performed, with a short delay (20 seconds) between each task to mark its end. Allsqueezing tasks were performed at the frequency of 0.75 Hz.PD subjects performed the task after an overnight withdrawal of their anti-Parkinsonian med-ication (a minimum of 12 hours since their last dose of L-dopa, and 18 hours since their last doseof dopamine agonists). All subjects practiced the task in a 30-min training session beforehand toget familiarized with the task requirements.233.2. Pre-ProcessingFigure 3.1: Experimental setup. Subjects are instructed to keep the horizontal bar within avertically-scrolling, 0.75 Hz sinusoidal pathway by squeezing and releasing the bulb.The scrolling pathway was either unaltered, or perturbed with a Gaussian noise equalto 25% or 50% of the target signal amplitude.3.2 Pre-ProcessingInspection of the raw signals for a typical PD subject immediately revealed multiple transients thatcould affect model accuracy, so we further explored suitable preprocessing strategies.We used the same steps mentioned in Section 2.1 during the pre-processing part: signal condi-tioning, time-lag removal, and re-sampling.• Signal conditioning: We low-pass filtered the data at 8Hz according to Chapter 2 [69]. Inaddition, since each identical experiment had been carried out three times, we took the meanof the input and output of the three experiments.• Time-lag removal We deployed a correlation-based time-lag removal technique mentionedin Chapter 2.• Resampling: The original sampling rate for our data was 60 Hz and the frequency of thesqueezing tasks was 0.75 Hz. Downsampling the data to 15 Hz turned out to improve modelaccuracies without losing their prediction capability.All of the mentioned preprocessing steps are provided as a block diagram in Figure 2.1. Pre-processed and raw signals are provided in Figure 3.2.243.2. Pre-Processing0 5 10 15−20002000 5 10 15−2000200Target and cursor position [mm]0 5 10 15−2000200Time [seconds]Figure 3.2: Typical PD subject’s tracking performance in different tasks. The target signal isunaltered in the top picture, perturbed with a Gaussian noise equal to 25% and 50%of the target signal amplitude respectively in the middle and the bottom picture. Thedashed red line represents target signal, the solid cyan line represents squeeze bulb rawoutput signal, and the solid blue line is the preprocessed output. As expected, as noiseamplitude increases, tracking performance deterioratesThe time lags removed by the correlation-based time-lag removal algorithm for different PDsubjects off medication and also normal subjects are provided in tables 3.2, and 3.3 respectively.The average delay for normal subjects was 0.10 ± 0.06 seconds while that of PD subjects was0.09 ± 0.05 seconds which showed no statistical significance between two groups (p = 0.4546,ANOVA). Note that the sampling time for the squeeze bulb data is 0.0167 seconds and these delaysare demonstrating about 10 sample shifts. The correlation of model parameters and time-lags areinvestigated in Section 3.4.2.253.3. Model CreationTable 3.2: Time lags (in seconds) removed by the correlation-based time-lag removal algorithm fordifferent PD subjects off medication. Note that the sampling time is 0.0167 seconds.No noise:25% of the noise:50% of thenoise signal amplitude signal amplitude1 0.0668 0.0835 0.03342 0.2171 0.0167 0.28393 0.2672 0.0835 0.20044 0.1503 0.1169 0.16705 0.0501 0.1169 0.03346 0.0334 0.0501 0.01677 0.1503 0.1169 0.06688 0.1503 0.1002 0.11699 0.1169 0.1336 0.066810 0.1336 0.0668 0.0835Table 3.3: Time lags (in seconds) removed by the correlation-based time-lag removal algorithm fordifferent normal subjects. Note that the sampling time is 0.0167 seconds.No noise:25% of the noise:50% of thenoise signal amplitude signal amplitude1 0.0167 0.0167 0.03342 0.1336 0.0668 0.05013 0.1837 0.2338 0.21714 0.1169 0.1002 0.15035 0.1336 0.1503 0.16706 0.1169 0.1002 0.13367 0.1503 0.1002 0.03348 0.1670 0.1169 0.01679 0.1169 0.0167 0.066810 0.0668 0.2171 0.01673.3 Model CreationWe build upon previous work in manual pursuit tracking [2, 98], to model the dynamics as a second-order LTI system. However, in contrast to their work, here we utilize bootstrapping methods(averaging hundreds of models which is described in Section 3.3.2 of this chapter) to create robustestimates of model parameters.Consider a continuous-time linear system that is represented in the canonical observer form:x˙(t) =[a0 1a1 0]x(t) +[b0b1]u(t)yˆ(t) =[1 0]x(t) + d0u(t)(3.1)For each set of input-output data (sampled sequences of u(t), y(t)), constant parameters a0, a1, b0, b1, d0in equation (3.1) are estimated by Matlab’s System Identification toolbox [80]. The parameters263.3. Model Creationwere numerically evaluated via Matlab’s pem (Prediction Error Minimization) function which isasymptotically unbiased and efficient for Gaussian distributed prediction errors. The function es-sentially tries to minimize the variance of the modeling error, i.e., the difference of model predictedoutput and the actual output given in equation (3.2) [81].w(t) = y(t)− yˆ(t) (3.2)Without loss of generality, we chose the observer canonical form to ensure the same staterepresentation x for each task. The observer canonical form simplifies conversion to transfer functionand allows us to easily assess the properties of the system as well. Note that different forms of statespace representations are convertible to each other through linear transformations [37].While we investigated other modeling frameworks (not shown) we chose the aforementionedformulation because: 1) Previous work has verified that second-order models are a reasonablemodel for manual pursuit tracking [2, 8, 98]. 2) The computational cost of averaging hundreds ofmodels makes smaller order models desirable.3.3.1 Post-Processing Whitening FilterIdeally, w(t) mentioned in equation (3.2) should be a white Gaussian noise whose temporallyadjacent samples are uncorrelated from each other. If this is not the case, then biased estimates inthe model parameters may exist. If model residuals are not temporally uncorrelated, a “whiteningfilter” [53] can be constructed through an autoregressive(AR) model fit to w(t):N∑i=0aiw(t− i) = e(t) (3.3)where e(t) is a white Gaussian noise. This filter is then applied to u(t), y(t), and w(t) inequation (3.4).y(t) = Gu(t) + w(t) (3.4)Then, the identification process is done for the filtered input and output data. The modelingerror of this model identification, i.e. e(t), is guaranteed to be white. This process was doneonce for every subject and task and an average whitening filter was constructed and applied to allinput-output sequences.273.4. ResultsTable 3.4: Accuracy of the models in different levels of noise for normal subjects.No noise:25% of the noise:50% of thenoise signal amplitude signal amplitudeBuilt fromraw 63.67± 7.12 55.93± 8.06 32.12± 12.13data (%)Built frompre-processed 79.53± 2.95 76.22± 4.67 70.39± 5.73data (%)3.3.2 BootstrappingOne model was built for each subject and task. Under model assumptions, permuting white,Gaussian residuals will not affect model estimates, and thus we created surrogate data throughequation (3.5).ynew(t) = yˆ(t) + permuted version of w(t) (3.5)Then new set of input-output data (sampled sequences of u(t), ynew(t)) were used for identification.This process was done 500 times to have an average model that is robust to numerical errors.The following model parameters were considered: damping ratio, natural frequency, decay rate,peak time, rise time, settling time, steady state output to the reference, and RMS error. Differenttremor scores for PD subjects gathered are: right and left rest tremor, right and left posturetremor, total tremor, right side tremor, left side tremor, as well as differentiated tremor (right sidetremor minus left side tremor). Since the common measure of tracking performance, namely RMSerror, was not significantly correlated with clinical PD measures, we investigated the correlation ofdynamical model parameters with tremor scores. We wanted to determine if a linear relationshipexisted between model parameters and clinical scores.3.4 Results3.4.1 Models QualityOverall, model accuracy was 73.42%± 5.76%. Model accuracies for different levels of noise beforeand after pre-processing for normal and PD subjects are provided in Tables 3.4 and 3.5, respectively.The numbers in the bottom row clearly demonstrate the enhanced performance of models trainedon the pre-processed biomedical data. We considered these accuracies sufficient to interrogate thederived models for motor performance. In addition to improvements in accuracy, the bootstrapping283.4. ResultsTable 3.5: Accuracy of the models in different levels of noise for PD subjects off-medication.No noise:25% of the noise:50% of thenoise signal amplitude signal amplitudeBuilt fromraw 51.23± 13.72 44.70± 12.85 14.21± 14.06data (%)Built frompre-processed 74.31± 4.87 71.27± 6.54 68.83± 9.83data (%)approach ensures robustness to outliers or transients that may have occurred in the original data.3.4.2 Correlation with Clinical MeasuresWe used Matlab’s robustfit function to investigate the correlation between model parametersand clinical scores. This command is an iteratively reweighted least-square regression methodthat de-weights the outliers present in the data points at each iteration. As we found that thetremor scores were non-gaussian, we first applied the logarithm as a normalizing transformation(See Section 2.2.1).The correlation of log(total tremor) and log(right side tremor) with log(damping ratio) wasstatistically significant for PD subjects off medication (p = 0.0009, and p = 0.0023 respectively).As this was a positive correlation (see Figure 3.3), more tremor leads to faster response. Notethat the model parameters are the difference of C-level and A-level (e.g., ζ = ζClevel − ζAlevel).The same pattern of significant correlations was found for B-level minus A-level parameters (e.g.,ζ = ζBlevel−ζAlevel). The absolute parameters for different noise levels did not significantly correlatewith clinical measures.While gross tracking measure (RMS error) did not significantly correlate with any of the clinicalmeasures, the ratio of damping ratio over RMS error did. Let us define a profit function for trackingthat incorporates both static properties (RMS error) and dynamic characteristics (damping ratio):P (u, y) = damping ratioRMS error (3.6)The higher the profit function in equation (3.6), the better the tracking performance. Interestingly,the correlation between tracking profit and total tremor and right side tremor was statisticallysignificant for PD patients(p = 0.0822e-3, and p = 0.0019 respectively). This positive correlationis depicted in Figure 3.4. A summary of significant correlations are provided in Table 3.6.To ensure that these correlations were due to biological phenomenon and not merely as a result293.4. Results0.8 1 1.2 1.4 1.6 1.8 2 2.2−2−10123Log (total tremor)Log (damping ratio) y = −3.5387 + 2.8708xR2 = 0.7913PD subjectsNormal subjectsFigure 3.3: The correlation of total tremor with damping ratio (p = 0.0009, robustfit). Note thatthe logarithm of the parameters are used since the original distribution of them was notGaussian. The damping ratio of the normal subjects are also provided for completeness.No tremor score is associated with normal subjects and they are provided at a lowtremor position for visualization purposes.Table 3.6: Correlation between model parameters and total tremor score.Clevel −Alevel Blevel −AlevelDamping ratio p = 0.0009 p = 0.0087R2 = 0.7913 R2 = 0.5596Damping ratio/ p = 0.0822e− 3 p = 0.0019RMS error R2 = 0.6973 R2 = 0.5442of time-lag removal, we investigated the correlation between lag shift amounts and damping ratios.There was no significant correlation detected (p = 0.8572, robustfit).3.4.3 Group AnalysisNormal subjects had a mean differential damping ratio of 0.96 ± 0.74, while PD subjects off-medication had a lower mean differential damping ratio of 0.70 ± 1.61. This difference was notstatistically significant (p = 0.6462, ANOVA) but shows faster dynamics for normal subjects. Also,as can be seen in Figure 3.5 PD subjects’ damping ratios have higher standard deviation (p =0.0319, VARTEST), which shows nonhomogeneous effect of the disease. Note that the damping303.4. Results0.8 1 1.2 1.4 1.6 1.8 2 2.2−0.8−0.6−0.4−0.200.20.40.60.81Log (total tremor)Log (damping ratio/rms error) y = −1.3904 + 1.0604xR2 = 0.6973PD subjectsNormal subjectsFigure 3.4: The correlation between total tremor and damping ratio over RMS error. The posi-tive correlation was statistically significant (p = 0.0822e-3, robustfit). Note that thelogarithm of the parameters are used since the original distribution of them was notGaussian. The damping ratio over RMS for the normal subjects are also provided forcompleteness. No tremor score is associated with normal subjects and they are providedat a low tremor position for visualization purposes.ratios are the difference of C-level and A-level (i.e., ζ = ζClevel − ζAlevel). The same pattern (higherdifferential damping ratios for normal subjects) was observed for B-level minus A-level dampingratios (i.e., ζ = ζBlevel − ζAlevel) (p = 0.6412, ANOVA). The absolute parameters for different noiselevels did not significantly differ across groups either.As we found that the damping ratios were non-gaussian, we first applied the logarithm as anormalizing transformation (See Section 2.2.1).For the profit function mentioned in equation (3.6), the following pattern was observed: Normalsubjects had a mean differential profit of 0.29±0.20. PD subjects off-medication had a lower meandifferential profit of 0.17± 0.49 which is lower than that of normal subjects (p = 0.4913, ANOVA).Also, PD subjects’ profit functions have higher standard deviation (p = 0.0148, VARTEST). Thisdifference is depicted in Figure 3.6. Note that the profit values are the difference of C-level and A-level in this case as well. The same pattern (higher profit fuction for normal subjects) was observedfor B-level minus A-level damping ratios (i.e., P = PBlevel − PAlevel) (p = 0.6412, ANOVA).313.5. DiscussionNormal PD off−medication−2−1.5−1−0.500.511.522.53Log(Damping Ratio)Figure 3.5: Damping ratio for normal subjects and PD subjects (p = 0.6462, ANOVA). Note thatthe logarithm of the parameters are used since the original distribution of them wasnot Gaussian.3.5 DiscussionOur results provide a novel insight into the pathophysiology of PD. The classic model of PD hasemphasized the “sole” role of the BG in modulating motor function through segregated striato-thalamo-cortical (STC) circuits, leading to bradykinesia and rigidity in the dopamine-deficientstate [4]. However it is increasingly evident that the classic model of the BG is insufficient toexplain the more benign prognosis of PDT compared with PDAR and why PD tremor is less re-sponsive to dopaminergic medication. Nevertheless, the motor system is also influenced by thecerebello-thalamo-cortical (CTC) circuit, in addition to the BG. Although typically investigatedindependently, these two circuits may dynamically influence one another. We speculate that otherbrain regions beside the BG may be crucial for PD tremor, which are highly influenced by CTCcircuit as a compensatory mechanism in PDT patients.Generally, normal subjects had higher damping ratios compared to PD subjects, which showsfaster response for the normal subjects. Also, the profit function was also generally higher thanthat of PD subjects supporting the faster tracking with less error. In addition, some PD subjects’damping ratios are even higher than that of normal subjects which could be related to super-normal activity of tremor-dominant patients that have been reported in [52, 76]. Both tremor andbradykinesia are symptoms of Parkinson’s disease that degrade movement. However, tremor could323.5. DiscussionNormal PD off−medication−0.6−0.4−0.200.20.40.60.81Log(damping ratio / rms error)Figure 3.6: Damping ratio over RMS error for normal subjects and PD subjects (p = 0.4913,ANOVA). Note that the logarithm of the parameters are used since the original distri-bution of them was not Gaussian.be beneficial in certain aspects [24]. Unlike tremor, bradykinesia did not significantly correlatewith any of the model parameters. These results support the commodious effect of tremor in motortracking tasks.The existence of compensatory mechanisms in PD may at least partially account for the find-ing that imaging measures of pathological disease progression, such as F-dopa Positron EmissionTomography (PET) and Single-Photon Emission Computed Tomography (SPECT), do not appearto correlate with clinical measures of symptom severity, such as the Unified Parkinson’s DiseaseRating Scale (UPDRS) [90, 142].Despite extensive clinical and basic studies of these empirical observations, the pathophysiologyof PD tremor has been difficult to decipher [36]. Lesion experiments in monkeys have revealedthat parkinsonian-like tremor can only be induced when there is damage to the cerebellum orits connections, in addition to impairments in nigrostriatal dopaminergic pathways [106]. PETscans have shown that resting activity of rostral, medial, and intermediate cerebellum (vermis andparavermis) is increased in PD tremor [39]. The cerebellum, traditionally associated with puremotor control, is now considered to be essential for the development of “forward models”, suchas predicting the sensory consequences of motor actions [16, 87, 88]. Activation of these areasare potentially responsible for better tracking performance, and thus we speculate that augmentedcerebellar activity, in order to compensate for diseased BG circuits, results in enhanced motor333.5. Discussiontracking performance, but leads to the epi-phenomenon of tremor.34Chapter 4Kinematics of RigidityOne of the earliest manifestations of PD is decreased arm swing [75]. Parkinsonian rigidity is acommon example of hypertonicity where in spite of subject’s effort to relax, there is a persistentsteady state of muscle contraction [104]. Increase in the muscle tone leads to resistance to passivemovement or stretch through the range of motion.Clinically, rigidity is one component of the Unified Parkinson’s Disease Rating Scale (UPDRS)[104], where stiffness of a manipulated joint is ranked on a subjective scale of 0–4, and thereforemay be susceptible to inter-rater variability. Traditionally, Wartenberg pendulum test [22] hasbeen used to measure rigidity in which the subjects lying supine or in a sitting position such thattheir limb could hang freely. Then, their limb is raised to a horizontal position against gravityand subsequently released. The number of oscillations is then recorded as a measure for rigidity.Wartenberg test has resulted in marked reduction of parkinsonian leg swing’s maximum velocitycompared to spastic stroke patients [22].In this chapter, we greatly extend the traditional test using motion sensor quantitative ac-celerometers and the physical properties of subject’s arm to fit linear dynamic models during thewhole period of the arm swing. We show that parameters of the fitted models sensitively detectchanges in the passive qualities of the limb. In clinical settings where massive instruments cannotbe used, this simple test will enable the clinicians to quantitatively measure rigidity as opposed tothe qualitative approaches used in UPDRS.The main advantages of using the modified pendulum test for objective quantification of rigidityare: 1) simplicity; not using huge biomechanical devices, 2) sensitivity; using the whole period ofswing as opposed to first swing excursion usually analyzed 3) reliability, safety, and reproducibility,4) accuracy; early detection and diagnosis of Parkinson’s disease have been reported using the armswing magnitude and asymmetry rather than the leg [14, 75].4.1 Experimental SetupThe experimental setup included a goniometer and two sensor units, electromyography sensing unit(Kinetisense, CleveMed, Cleveland, OH) and the Urias air splint (length of 8 inch and maximum354.1. Experimental SetupFigure 4.1: Experimental setup. Subjects were instructed to relax to the best of their abilitieswithout interfering with the pendulum motion. Then, the arms of the subject werereleased to fall without informing the subjects.circumference of 12.5 inches ) which prevented the flexion of the elbow. The goniometer was used tomeasure the arm joint angle. One of the surface EMG electrodes was attached to the anterior belly ofthe biceps brachii, while the other on the lateral head of triceps brachii (See Figure 4.1. Each sensorunit included three orthogonally positioned accelerometers with gyroscope-measuring units along 3axis (x, y, z). They measure the linear acceleration along each direction (e.g., ay = g = 9.8m/s2)and angular velocity (in deg/sec), respectively. For each muscle, two EMG electrodes where placedone inch apart. The motion sensor data sampling at 128 Hz was synchronized with user unitsampling EMG at 2048 Hz. The command module with an embedded 2.4 GHz radio providedreal-time wireless data transmission to a computer. Each subject was set up with 2 kinetisensesensor units on the skin using medical tape positioned on ventrolateral aspect of wrist and on theanterior surface of the shoulder joint.Ten PD subjects (on and off L-dopa medication) with clinically diagnosed, mild to moderate PD(Hoehn and Yahr stage 2–3) [54], and ten healthy, age-matched subjects without active neurologi-cal disorders were recruited for this study. Subjects performed the task at the Pacific Parkinson’sResearch Centre (PPRC)/Movement Disorders Clinic [95]. They provided written informed con-sent, and all research was approved by the University of British Columbia Clinical Research Ethics364.2. Problem FormulationTable 4.1: Detailed characteristics of PD subjects. For rigidity, ‘LH’ denoted left hand, ‘RH’ de-notes right hand, ‘LL’ denotes left leg, ‘RL’ denotes right leg, and ‘N’ denotes neck.Subj. Gender Handed Mass Arm Arm UPDRS RigidityNo. (lb) length radius LH,RH,LL,RL,N1 M R 175 22.5” 1.45” 18 0, 0, 0, 1, 12 F R 120 20” 1.48” 22 0, 1, 0, 1, 03 M L 197 22” 1.75” 21 1, 2, 1, 2, 24 F R 130 22” 1.32” 32 2, 1, 2, 2, 25 M R 162 20.5” 1.48” 35 2, 2, 2, 2, 26 F R 175 19.5” 1.42” 32 1, 1, 1, 2, 27 M R 165 20.5” 1.34” 22 1, 1, 1, 1, 18 M R 190 23” 1.67” 24 1, 1, 1, 1, 29 M R 155 21.5” 1.53” 21 1, 1, 1, 1, 110 M R 166 22.5” 1.57” 10 1, 1, 0, 1, 0Board. PD subjects’ clinical Unified Parkinson’s Disease Rating Scale (UPDRS) score, includingrigidity scores, were assessed by a physician trained in the UPDRS (see Table 4.1 for demographicdetails). Note that the arm radius in the table is the average of wrist, elbow, and shoulder radiuses.Subjects were in the seated position on an armless chair without head or back support. Sinceparkinsonian rigidity is typically more pronounced in extension than flexion, the EMG of triceps wasmonitored. Because of the air splint, the arm swing in the anteriorposterior plane is ensured. Twomotion sensors were attached to the wrist and shoulder joints. The subject’s arm was lifted whereit reached a horizontal status, then the arm was released by the examiner, permitting the arm toswing freely. The final arm position was in vertical position at approximately 0◦ after undergoing apendulum-like swinging motion. Subjects were instructed to be relaxed to the best of their abilitieswithout interfering with the pendulum motion. Five seconds later, the data collection was started.Then, the arm of the subject was released to fall without informing the subjects. Each experimentwas performed twice for each hand. The EMG data were observed during the experiment in orderto discard data sets with large activity of EMG and therefore, ensuring the involuntary swingingmotion of the arm.PD subjects performed the task after an overnight withdrawal of their anti-Parkinsonian med-ication (a minimum of 12 hours since their last dose of L-dopa, and 18 hours since their last doseof dopamine agonists) and again 45 minutes following the administration their prescribed dose ofLevodopa.374.2. Problem FormulationLh2rmgCGFigure 4.2: A pendulum. CG is the center of gravity.4.2 Problem FormulationAs we know, for a pendulum of length ‘L’ shown in Figure 4.2 that has small swings (sin(θ) = θ),in the presence of friction we can write:Jθ¨ = −bθ˙ −mghθ (4.1)in which ‘J ’ is pendulum moment of inertia, ‘b’ is the damping coefficient, ‘m’ is pendulum mass,and ‘h’ is the distance of the center of gravity from the pivot that will be ‘L/2’ for a uniformpendulum. For a cylinder with radius ‘r’ that has a uniform mass, ‘J ’ could be calculated via:J = 13mL2 + 14mr2 (4.2)The mechanical energy of the system will diminish in time. Hence, the angular position of thependulum would be underdamped and defined by:θ(t) = θmaxe−b2J t cos(ωdt+ φ) (4.3)in whichωd =√mghJ −b24J2 (4.4)384.2. Problem Formulationand ‘φ’ is a constant phase to take into account the initial position. The maximum swing is shownas θmax.Now, we could consider the arm swing (without elbow flexion) as a uniform mass cylinder swing.We have experimental data for angular velocity θ˙(t) (measured in deg/sec). The angular positionθ(t) could be written as:θ(t) = θ0e−b2J t cos(ωdt+ φ) (4.5)in which θ0 cosφ is the initial position.Theoretical θ˙(t) could be calculated by taking the derivative of equation (4.5):θ˙(t) = θ0e−b2J t(− b2J cos(ωdt+ φ)− ωd sin(ωdt+ φ)) (4.6)In equation (4.6), ‘J ’ is known (according to equation (4.2)) while ‘θ0’, ‘ωd’, ‘b’, and ‘φ’ are notknown. Among ‘ωd’, and ‘b’ only one of them is considered unknown as they are dependent to eachother according to equation (4.4)–Note that ‘m’, ‘g’ are known and ‘h = L2 ’. Among ‘φ’, ‘θ0’, and‘b’ only one of them is unknown as they are dependent to each other according to the followingequations:θ(0) = 90 =⇒ θ0 = 90cos(φ) (4.7)θ˙(0) = 0 =⇒ φ = tan−1( −b2Jωd ) (4.8)Hence, equation (4.6) could be rewritten as:θ˙(t) = 90cos(φ)e− b2J t(− b2J cos(ωdt+ φ)− ωd sin(ωdt+ φ)) (4.9)in which ‘φ’ and ‘ωd’ are functions of ‘b’. So, the only unknown parameters here is ‘b’. Notethat θ(0) = 90 (and not pi2 ) as the angular velocity is measured in deg/sec.Now, by having the experimental θ˙(t), we can estimate ‘b’ via least squares (LS) optimizationwhich is explained in the sequel (See Section 4.3). Having ‘b’ estimated, we can calculate modelaccuracy as well as damping ratio and natural frequency for each subject. Recall from equation(4.1) thatωn =√mghJ =√gL213L2 + 14r2(4.10)394.3. Model Estimationζ = b2Jωn (4.11)ωn is known for each subject and does not depend on ‘b’ as could be seen in equation (4.10).However, it would be different for different subjects as they have different arm lengths, and radiuses.We could also calculate ‘J ’ for each subject according to equation (4.2). Note that the mass of thearm is considered roughly as 0.1 of the body mass for each subject. We measured arm length andradius for each subject as well. Thus, by estimating ‘b’ via LS optimization we can compute modelparameters.4.3 Model EstimationAmong all of the preprocessing steps mentioned in Chapter 2 only signal conditioning is used here.There is no input to this system as the arm is passively swinging and hence no input-output timelag is defined. Nonlinear effects are also mitigated to a decent extent using the air splint. Wehad angular velocity in x, y, z directions. The sensor on the wrist was mounted in a way thatthe direction along the arm was ‘y’, the one perpendicular to the palm was ‘z’, and the third onewas ‘x’. Since the subjects’ arm might have not necessarily swung in the xy plane (the plane thatis perpendicular to the body plane), singular value decomposition (SVD–also known as principalcomponent analysis [64]) was done to find linearly uncorrelated directions. The first principalcomponent was chosen as the angular velocity in the new ‘x’ direction. A sample angular velocityof the arm swing (raw output) is depicted in Figure 4.3.We carefully inspected each of the trials to ensure that there are no large activities of EMGcausing voluntary resistance to swings. We also selected the initial time point for the output signalmanually. We, essentially did the modeling with an output signal starting from the first nonzerovelocity logged by the sensor. To avoid artifacts, we ensured the nonzero values last for at least500 milliseconds. The sampling rate for this data set was 128 Hz. We detrended the data to makeit zero-mean for system identification. We then low-pass filtered the data at 8 Hz to get rid of thehigh-frequency content.A simple pendulum swing is governed by a second order equation stated in Section 4.2 (seeequation (4.1)). Considering θ(t) and θ˙(t) as state variables, the following continuous-time linear404.3. Model Estimation0 1 2 3 4 5 6 7 8−5−4−3−2−1012Time (seconds)Angular Velocity(rad/sec)Figure 4.3: Typical PD subject’s arm swing angular velocity.system can describe the system:x˙(t) =[0 1−a0 −a1]x(t)yˆ(t) =[0 1]x(t)(4.12)Note that a0 does not depend on the output data in this framework. It only depends on physicalcharacteristics of the arm and can be computed via the following equation:a0 = mghJ (4.13)The parameter a1, on the other hand, has a direct relationship with the damping term b viaequation (4.14).a1 = bJ (4.14)The relationship between a0, a1 and damping ratio and natural frequency are as follows:ωn = √a0 (4.15)414.3. Model Estimationζ = a12√a0 (4.16)For each trajectory of output data (sampled sequence of y(t)), constant parameter a1 in equation(4.12) has to be estimated by least squares optimization. Without loss of generality, we chose thisstate space form to have the position and the velocity as state vector x. We chose the C matrixas [0 1] instead of [1 0] to have the angular velocity as the output. This makes comparisonswith the experimental data easier as angular velocity is measured by the sensor. This canonicalform simplifies conversion to transfer function and allows us to easily assess the properties of thesystem as well. Note that different forms of state space representations are convertible to eachother through linear transformations [37].To find b, the following optimization problem has to be solved to minimize the difference ofmodel predicted output and the actual output:minb∈U⊂R+N1∑k=0Q( 90cos(φ)e− b2J t[k](− b2J cos(ωdt[k] + φ)− ωd sin(ωdt[k] + φ))− θ˙[k])2+N∑k=N1+1R( 90cos(φ)e− b2J t[k](− b2J cos(ωdt[k] + φ)− ωd sin(ωdt[k] + φ))− θ˙[k])2(4.17)in which ‘φ’, ‘ωd’ are also functions of ‘b’ as mentioned in equations (4.4), (4.8) and ‘N ’ is thenumber of samples. ‘Q’ and ‘R’ are the weighting factors for the beginning part of the swing (thefirst ‘N1’ samples that start from arm drop until the arm velocity becomes zero for the first time)and the leftover of the swing (the N −N1 samples that start right after the first time that the armreaches zero velocity to the end of the swing). As the frequency of damped oscillations (ωd) hasto be real, a constraint on b is imposed for this optimization problem that is provided in equation(4.18).0 ≤ b ≤ 1√4Jmgh (4.18)The aforementioned optimization problem is not convex and hence numerical methods are de-veloped to deal with that. We used Matlab’s fmincon (Constrained Minimization) function whichis efficient for non-convex optimization problems. This command inherently uses SQP (sequentialquadratic programming) to solve the problem. Since this method (like most other numerical meth-ods) is vulnerable to get trapped in local minima, we tried different initial guesses and did theoptimization starting from them. The one with higher accuracy and visually captured more data424.3. Model Estimationcharacteristics was chosen.One issue arose in implementing the least squares algorithm. The algorithm tries to minimizethe error between the predicted output and the experimental output for the “entire” period ofswing. This may lead to models that have a good prediction of output, but fail to estimate thedamping ratio and/or decay rate correctly. In other words, in an attempt to minimize the errorfor the whole period, the errors in the beginning of the swing (the first peak and valley) may besacrificed because of the number of swings at the end with lower amplitudes. This will severelyaffect the estimation of damping ratio and therefore decay rate. To mitigate this effect, we ran theLS algorithm first. Then for those that did not match well, we improved the models with emphasison the data in the beginning and middle of the experiment via the following algorithm:1. Estimate ‘b’ via LS algorithm and then calculate the damping ratio through equation (4.11).2. Hypothesize an improved damping ratio in order to better match data. Note that the naturalfrequency only depends on physical properties of the arm. Convert improved damping ratio[and natural frequency] to damping term ‘b’ in (4.11).3. Initialize least squares algorithm with the improved system parameters considering more‘weight’ on the beginning and the middle of the swing.4. Simulate and form the output for the newly estimated system, and the original system.Compute corresponding fitness scores and visually observe the qualitative fit.5. Repeat steps 2–4 as necessary, and select the “best” model.This issue mainly affected the PD subjects off medication and some of PD subjects on medi-cation. Normal subjects in general turned out to have models that fairly captured the data withgood accuracy and good visual fit even with equal weighting of the errors for the entire swing. Asample of the experimental output and model predicted output (both via normal LS and improvedLS) for normal and PD subjects is provided in Figures 4.4, 4.5, and 4.6. Note that for visualpurposes since the angular velocity was measured by the sensor, the angular velocity (instead ofthe angle) is shown in the figures.Accuracy of the models can be assessed by a “fitness score” which depends on the differencebetween the actual output y and modeled output yˆ, normalized by the variations in the outputfrom its mean value y¯ [81] (See equation (1.8) in Section 1.3.Note that not necessarily in every case the improved LS led to a better fit nor better capturingof the system dynamics. Between the output formed by improved LS method and the one formed byLS method, we selected the one which better captures system dynamics (and has a good accuracy).434.3. Model Estimation0 1 2 3 4 5 6 7 8 9−200−150−100−50050100150200250angular velocity (deg/sec)Time (sec)A typical normal subject experimentaloriginal modelupdated modelFigure 4.4: Typical Normal subject’s arm swing angular velocity. The dashed blue one is experi-mental, the dotted purple one is estimated by LS method, and the dotted black one isestimated by LS after hypothesizing new damping ratio and re-run LS weighting thebeginning part more.The improvement in the capability of the model for capturing system dynamics, did not significantlychange after improved LS for most of normal and PD subjects on medication. However, for PDsubjects off medication the improved LS led to a model which captures system dynamics muchbetter than LS. The weights for improved LS are chosen in a way to have more emphasis on theerrors of beginning of the swing. The weights for different subjects are chosen by trial and errorand the ones that led to a good fit are provided in Table 4.2.Even using the weighted LS algorithm did not lead to good models for more than half of PDsubjects off medication and 2 PD subjects on medication. The reason is the nonlinear behaviorof these subjects in the beginning of the swing causing the swing profile to have a huge overshootin the beginning. The swings will damp fairly quickly after that overshoot. Such behavior is notretractable by linear second order models. If more weighting is allocated to the beginning of theswing (larger Q) to capture that overshoot, some oscillations will appear in the model output atthe end of the swing profile while the experimental data damp so quickly. On the other hand, if lessweighting is used for the beginning of the swing (larger R) to capture the dynamics of the swing,that huge overshoot is underestimated by the model and the fidelity of the model parameters such444.3. Model Estimation0 1 2 3 4 5 6−150−100−50050100150200250angular velocity (deg/sec)Time (sec)A typical PD subject off medication experimentaloriginal modelupdated modelFigure 4.5: Typical off-medication PD subject’s arm swing angular velocity. The dashed blue oneis experimental, the dotted purple one is estimated by LS method, and the dotted blackone is estimated by LS after hypothesizing new damping ratio and re-run LS weightingthe beginning part more.as damping ratio (or decay rate) is questionable.To capture the dynamics of this behavior, yet not using complicated nonlinear models, weassumed the swing profile to be a multiplication of an exponential term and a sinusoidal term(same as linear models). However, we assumed the decay rate of the envelope to be independentof the damped frequency of the swings. This way the angle is defined by the following equation:θ(t) = θ0e−at cos(ωdt+ φ) (4.19)in which ‘a’ is the “effective decay rate”. Note that in the linear model the angle trajectory wasgoverned by equation (4.5) in which the decay rate was b2J and ωd was also defined by equation(4.4). The difference of this model and the linear model is that the envelope of the exponentialterm (e−at) does not depend on ωd. An interpretation of such modeling is that the weighting of Qand R that we just mentioned is incorporated in an adaptive manner, i.e., a time-dependent weightis considered for optimization that considers more weight for the beginning of the swing and astime passes, this weighting becomes less and less.454.3. Model Estimation0 1 2 3 4 5 6 7 8−200−150−100−50050100150Angular Velocity (deg/sec)Time (sec)A typical PD subject on medication experimentaloriginal modelupdated model Figure 4.6: Typical on-medication PD subject’s arm swing angular velocity. The dashed blue oneis experimental, the dotted purple one is estimated by LS method, and the dotted blackone is estimated by LS after hypothesizing new damping ratio and re-run LS weightingthe beginning part more.Hence the angular velocity profile mentioned in equation (4.9) is rewritten as:θ˙(t) = 90cos(φ)e−at(− a cos(ωdt+ φ)− ωd sin(ωdt+ φ)) (4.20)and hence the optimization problem mentioned in equation (4.17) is modified as follows:minb,a∈U⊂R+N1∑k=0Q( 90cos(φ)e−at[k](− a cos(ωdt[k] + φ)− ωd sin(ωdt[k] + φ))− θ˙[k])2+N∑k=N1+1R( 90cos(φ)e−at[k](− a cos(ωdt[k] + φ)− ωd sin(ωdt[k] + φ))− θ˙[k])2(4.21)We still keep Q,R weighting factors to get the best possible trajectory. The initial guess for ‘a’was chosen as the decay rate of the envelope of the swing profile which is obtained by fitting anexponential curve into the peaks of the absolute value of angular velocity. The initial guess for b, asbefore, was chosen within the boundaries defined in equation (4.18). A sample of this method that464.4. ResultsTable 4.2: The weights used for errors in normal LS and improved LS methods. The errors are inthe format of [Q*(error before t seconds) + R*(error after t seconds)]Subject The error used in improved LS The error used in improved LSfor left arm for right armN1 equal weights (Q=R) equal weights (Q=R)N2 Q=10, R=1, t=1.5 equal weights (Q=R)N3 Q=10, R=1, t=1 Q=10, R=1, t=1N4 Q=10, R=1, t=1.5 equal weights (Q=R)N5 equal weights (Q=R) equal weights (Q=R)N6 equal weights (Q=R) equal weights (Q=R)N7 Q=10, R=1, t=1 equal weights (Q=R)N8 equal weights (Q=R) equal weights (Q=R)N9 equal weights (Q=R) equal weights (Q=R)N10 equal weights (Q=R) Q=10, R=1, t=1PD1 off med equal weights (Q=R) Q=10, R=1, t=1.5PD2 off med Q=10, R=1, t=1.5 Q=10, R=1, t=1PD3 off med Q=10, R=1, t=1 Q=10, R=1, t=1PD4 off med Q=10, R=1, t=1 Q=10, R=1, t=1PD5 off med Q=10, R=1, t=1.5 Q=10, R=1, t=1.5PD6 off med Q=10, R=1, t=1 Q=10, R=1, t=1PD7 off med equal weights (Q=R) equal weights (Q=R)PD8 off med Q=10, R=1, t=1.5 equal weights (Q=R)PD9 off med Q=10, R=1, t=1.5 Q=10, R=1, t=1.5PD10 off med Q=10, R=1, t=1.5 equal weights (Q=R)PD1 on med No-data No-dataPD2 on med equal weights (Q=R) equal weights (Q=R)PD3 on med equal weights (Q=R) Q=10, R=1, t=1.5PD4 on med equal weights (Q=R) Q=10, R=1, t=1.5PD5 on med equal weights (Q=R) Q=10, R=1, t=1PD6 on med equal weights (Q=R) Q=10, R=1, t=1PD7 on med equal weights (Q=R) equal weights (Q=R)PD8 on med Q=10, R=1, t=1 Q=10, R=1, t=1PD9 on med equal weights (Q=R) equal weights (Q=R)PD10 on med Q=10, R=1, t=1.5 Q=10, R=1, t=1.5led to a better model for one of the PD subjects with nonlinear behavior is depicted in Figure 4.7.Note that for normal subjects and most of PD subjects on medication (and some of PD subjectsoff medication) the effective decay rate is simply equal to the decay rate (ζωn) as the swing profileis captured by a linear model.4.4 ResultsOverall, model accuracy was 72.1% ± 10.5%. Model accuracies were higher for normal subjectsas more oscillations were accomplished. We believe these accuracies are sufficient to explore thederived models’ parameters for motor performance. Note that in some cases a model with a lower fit474.4. Results0 0.5 1 1.5 2 2.5 3 3.5 4 4.5−300−250−200−150−100−50050100A PD subject off medication with nonlinear behaviorTime (sec)Agular velocity (deg/sec) experimentaloriginal modelmodel with effective decayFigure 4.7: An off-medication PD subject’s arm swing angular velocity with nonlinear behavior.The dashed blue one is experimental, the dotted purple one is estimated by weightedLS method for linear model, and the dotted black one is estimated by LS after usingthe model with effective decay rate.was chosen as the dynamics of the data was better captured visually. We considered the followingmodel parameters: natural frequency (ωn), decay rate (ζωn), effective decay rate (a), decay ratecorrection (a− ζωn), damping term (b), peak time, rise time, settling time, steady state output tothe reference, and a static parameter which is maximum velocity of the swing. We did not considerthe damping ratio (ζ) as the nonlinear behavior of PD subjects and having huge overshoots mightinterfere with our results. Thus decay rate would be a better candidate to investigate how fast thearm swing decays. The parameters for subjects are provided in Tables 4.3, 4.4, 4.5, 4.6, and 4.7.Note that for PD subjects on medication the data for the first subject was corrupted. So, for theanalysis related for medication effect, we only considered 9 subjects. For group comparisons, thepopulation of PD subjects off medication had 20 elements (10 for left and 10 for right side) whilethe population of PD subjects on medication had 18 elements (9 for left and 9 for right side).Different clinical measures for PD subjects gathered are: UPDRS, right and left rest tremor,right and left posture tremor, total tremor, right side tremor, left side tremor, bradykinesia, andrigidity. We investigated the variations of model parameters across different groups. We alsoinvestigated the correlation of dynamical model parameters with rigidity and UPDRS scores.484.4. ResultsTable 4.3: Detailed dynamic parameters of Normal subjects.Subj. Natural Damping Ratio Damping Ratio Accuracy AccuracyNo. Frequency Left Arm Right Arm Left Arm(%) Right Arm(%)1 5.2440 0.0913 0.0837 81.23 79.522 5.4997 0.1198 0.1444 79.24 77.473 5.3722 0.1825 0.1206 54.59 54.044 5.5881 0.1646 0.1167 45.78 69.445 5.3726 0.1563 0.2664 81.06 68.916 5.5096 0.0989 0.1543 69.79 69.987 5.5825 0.0701 0.0951 58.52 62.058 5.1837 0.2419 0.3953 62.95 67.519 5.5307 0.1285 0.1395 62.51 62.9110 5.5064 0.4577 0.3956 62.61 54.90Table 4.4: Detailed dynamic parameters for the left arm of PD subjects off medication. The entriesshown as ζωn are the ones that no effective decay is used for their modeling and thereforethe effective decay rate is just the multiplication of damping ratio and natural frequency.Subj. Natural Damping Ratio Effective AccuracyNo. Frequency Left Arm Decay Left Arm(%)1 5.0663 0.1510 ζωn 79.232 5.3710 0.0763 ζωn 74.853 5.1194 0.3567 ζωn 64.034 7.1585 0.1808 2.6803 65.345 9.0987 0.1574 3.8271 71.926 5.4397 0.5326 1.3608 56.117 5.3076 0.3784 ζωn 86.168 6.5012 0.3501 0.9515 57.799 7.9020 0.3621 1.3875 53.4610 5.0650 0.2912 1.5980 77.144.4.1 Variations Across GroupsThe following Model parameters for all subjects were statistically significant across groups: decayrate (ζωn, p = 0.0029, ANOVA), effective decay rate (a, p = 0.0185, ANOVA). Note that equation(4.11) could be rewritten as:ζωn = b2J (4.22)Hence, any significance or correlation deciphered for decay rate holds for the damping term b aswell. Normal subjects had a mean decay rate of 0.95±0.59, while PD subjects off-medication had amean decay rate of 1.68±0.70 (See Figure 4.8). The mean decay rate for PD subjects on-medicationwas 1.11± 0.70. Similar pattern was observed for the effective decay rate (See Figure 4.9).While the aforementioned graphs include values from both left and right arm, we note that asignificant difference was also seen across groups for non-dominant hand decay rate and effective494.4. ResultsTable 4.5: Detailed dynamic parameters for the right arm of PD subjects off medication. Theentries shown as ζωn are the ones that no effective decay is used for their modeling andtherefore the effective decay rate is just the multiplication of damping ratio and naturalfrequency.Subj. Natural Damping Ratio Effective AccuracyNo. Frequency Right Arm Decay Right Arm(%)1 5.0663 0.2414 ζωn 56.662 5.3710 0.1582 ζωn 74.763 5.1194 0.5237 ζωn 67.044 4.7531 0.4756 ζωn 77.825 7.1540 0.1545 ζωn 74.036 5.9123 0.2843 0.8136 58.447 5.3076 0.1603 ζωn 79.778 6.5012 0.3194 2.1134 58.599 5.1810 0.3462 1.7199 69.4510 5.0650 0.3808 2.0737 75.63Table 4.6: Detailed dynamic parameters for the left arm of PD subjects on medication. The entriesshown as ζωn are the ones that no effective decay is used for their modeling and thereforethe effective decay rate is just the multiplication of damping ratio and natural frequency.Subj. Natural Damping Ratio Effective AccuracyNo. Frequency Left Arm Decay Left Arm(%)1 No− data No− data No− data No− data2 5.3710 0.0706 ζωn 79.323 5.4175 0.0969 1.1033 59.554 5.1246 0.5045 ζωn 72.915 5.3056 0.2737 ζωn 56.456 5.4397 0.2563 ζωn 85.657 5.3076 0.2450 ζωn 88.138 6.0074 0.3336 2.2168 58.109 5.1810 0.2209 ζωn 83.4410 5.0650 0.3805 ζωn 62.23decay rate (p = 0.0035 and p = 0.0398, ANOVA respectively) but no significance across groups forthe dominant hand was observed.Table 4.8 shows the parameters for both normal and PD subjects. Note as shown in Figures 4.8and 4.9, PD subjects off medication have the highest decay rates and effective decay rates amongthe three groups while normals have the lowest. Normals have the fastest velocity of swing amongthe three groups, while PD subjects off-medication are the slowest. These patterns are consistentwith previous studies [22, 138, 147], demonstrating more damped response in PD subjects.504.4. ResultsTable 4.7: Detailed dynamic parameters for the right arm of PD subjects on medication. Theentries shown as ζωn are the ones that no effective decay is used for their modeling andtherefore the effective decay rate is just the multiplication of damping ratio and naturalfrequency.Subj. Natural Damping Ratio Effective AccuracyNo. Frequency Right Arm Decay Right Arm(%)1 No− data No− data No− data No− data2 6.0536 0.0318 ζωn 59.653 5.1194 0.4268 2.8203 79.414 4.4521 0.0175 1.3225 56.305 5.3056 0.2172 ζωn 72.426 4.2821 0.2154 0.9322 69.867 5.3076 0.1193 ζωn 53.768 6.0074 0.0294 2.0204 58.449 5.1810 0.2030 ζωn 78.7910 5.0650 0.2193 ζωn 68.35Table 4.8: Static and dynamic model parameters for normal and PD subjects.Normal PD PD Statisticalon-med. off-med. significanceDecay rate 0.95± 0.59 1.11± 0.70 1.68± 0.70 p = 0.0029,ANOVAEffective 0.95± 0.59 1.35± 0.69 1.62± 0.83 p = 0.0010,decay rate ANOVAMax. Velocity 6.00± 1.48 4.92± 1.18 4.78± 1.60 p = 0.0185,[rad/sec] ANOVA4.4.2 Symmetry of OscillationsWe further investigated the symmetry of oscillations in dominant and non-dominant hand for eachsubject. To this end, we formed new differential parameters as the difference of dominant andnon-dominant hand for each group and subject such as:(ζωn)diff = (ζωn)non-dominant − (ζωn)dominant (4.23)Then, we explored the significance of these parameters across groups. Interestingly, normalshad more symmetric oscillations compared to PD subjects. PD subjects on medication had moresymmetric oscillations compared to PD subjects off medication (See Figure 4.10).Normal subjects had a mean differential decay rate of −0.12 ± 0.31, while PD subjects off-medication had a mean differential decay rate of 0.34 ± 0.76. PD subjects on medication had amean differential decay rate of −0.76± 1.06 and these differences were statistically significant (p =514.4. ResultsNormal PD on−medication PD off−medication00.511.522.53Decay RateFigure 4.8: Decay rate for normal subjects and PD subjects on and off medication.0.0148, ANOVA). Similar pattern was seen for maximum velocity (See Table 4.9). No significantdifference was seen between normal and PD subjects on medication. These results are consistentwith previous studies stating asymmetry in Parkinson’s disease [75].4.4.3 Medication EffectAs can be seen in Table 4.8, model parameters of PD subjects on-medication are closer to thatof normals. Interestingly, the values have been more homogenized (having less variance) aftermedication. This could be related to the homogenizing effect of the medication.We used paired ttest to explore the statistical significance of medication effect on each subject.We formed new differential parameters as the difference of on and off medication state for eachsubject, e.g.,(ζωn)med-diff = (ζωn)on-med − (ζωn)off-med (4.24)Consistent with results in the Section 4.4.1, we found improvements in decay rate. This dif-ferential decay rate derived in equation (4.24) was significantly different from zero (p = 0.0084,ttest) showing medication decreases decay rate. No significant difference for effective decay rate524.4. ResultsNormal PD on−medication PD off−medication0.511.522.533.54Effective Decay RateFigure 4.9: Effective decay rate for normal subjects and PD subjects on and off medication.nor maximum velocity were seen between before and after medication state. We note that thedifferential decay rate and the differential decay rate correction for the dominant hand was signifi-cantly different from zero (p = 0.0225 for ζωn and p = 0.0329 for a− ζωn), but those parameterswere not significantly different from zero for the non-dominant hand (p = 0.0665 and p = 0.8534respectively). This could be due to the fact that medication does not symmetrically affect bothhands. The visualization of decay rate for individual subjects have been depicted in Figures 4.11,and 4.12.4.4.4 Correlation with RigidityWe used Matlab’s robustfit function to investigate the correlation between model parametersand clinical scores. This command is an iteratively reweighted least-square regression method thatde-weights the outliers present in the data points at each iteration. The correlation of UPDRS witheffective decay rate (average value of left and right arms) was statistically significant for PD subjectsoff medication (p = 0.0024). We also note that the effective decay rate for the non-dominant handindividually increased with an increase in UPDRS but not statistically significant (p = 0.1061).534.5. DiscussionNormal PD on−medication PD off−medication−2−1.5−1−0.500.51Differential Decay RateFigure 4.10: Differential decay rate for normal subjects and PD subjects on and off medication.Parameters used here are the difference of the ones for dominant hand and non-dominant hand.Rigidity, as a component of UPDRS, also significantly correlated with effective decay rate. Therelationship of total rigidity and effective decay rate (average of left and right arm) is providedin Figure 4.13 (p = 0.0007). Non-dominant side rigidity almost significantly correlated with non-dominant arm effective decay rate (p = 0.0560). No significance was seen for the dominant side(p = 0.7926). All of these were positive correlations that are consistent with previous correlationsmeaning more severe state of the disease leads to more damped response.4.5 DiscussionOur results are consistent with prior studies that made speculations about the correlation of damp-ing coefficient (a.k.a decay rate) with rigidity [117, 138]. However, those works only address the firstswing excursion or at most two swings to approximately measure damping term. The advantageof this work is in consideration of the whole period of arm swing to model its dynamical system.In addition, we have incorporated the physical property of the arm such as radius and length to544.5. DiscussionTable 4.9: Static and dynamic differential model parameters for normal and PD subjects. Pa-rameters used here are the difference of the ones for dominant hand and non-dominanthand.Normal PD PD Statisticalon-med. off-med. significanceDifferential −0.12± 0.31 −0.76± 1.06 0.34± 0.76 p = 0.0148,Decay rate ANOVADifferential 0.49± 0.44 0.76± 0.59 1.52± 1.32 p = 0.0398,Max. Velocity ANOVAhave patient-specific models with high fidelity. This will provide a novel insight into an objectivemeasure for quantifying rigidity.In addition, the improvements in decay rate were observed after investigating the significanceof medication effect on each subject. Thus, the effect of medication was manifested by reductionin decay rate (and effective decay rate in some PD subjects). Arm swing belongs to the motorfunction domain and the presence of motor asymmetry increases the accuracy of PD diagnosis[75]. The asymmetry of the arm swing oscillations between the dominant and non-dominant handwas improved for all of the model parameters after medication. Our models were able to depictthe positive effect of medication on arm swing asymmetry among PD subjects. Thus, arm swingasymmetry can be used for early and differential diagnosis of PD and our models had high sensitivityfor detecting changes in rigidity due to administration of levodopa.The positive correlation of total UPDRS and quantified rigidity (both non-dominant side, andtotal score) showed that more severe state of the disease leads to more damped response and the armswing oscillations would decay faster. Therefore, the dynamic parameters of this test are correlatedwith clinical assessment of rigidity. This is a novel contribution to the sparse literature on objectivemethods to quantify rigidity among Parkinson’s disease patients. The pendulum movement doesnot impose any constant velocity and the increase in muscle tone is more pronounced in extension.In contrast to lead-pipe rigidity, which refers to the relative uniformity of rigidity that occurs inall direction, cogwheel rigidity is not diagnostic and may be encountered in any condition withrigidity and tremor [21]. Since the lead-pipe characteristic of rigidity is due to co-contraction ofagonist and antagonist muscles, using this test and ensuring that the subjects have no neck orback support would lead to lead-pipe rigidity with disappearance of cogwheel phenomenon [21].However none of the previous studies have reassured whether they are quantifying the cogwheel orlead-pipe rigidity. Since reduction of arm swing magnitude as well as arm swing asymmetry havebeen reported as early signs for PD diagnosis, we believe that this simple diagnostic test providesan objective quantitative measure of rigidity and an early marker for the disease.554.5. DiscussionNormal PD off med PD on med00.511.522.533.544.555.5Decay Rate−−non−dominant armFigure 4.11: The decay rates for non-dominant arm.Wartenberg pendulum test has shown to differentiate spasticity from rigidity quantitatively [22].Clinically, it is normally taught that spasticity is a velocity dependent resistance to movement, whilerigidity is velocity independent. PD subjects, who are normally considered to have “pure” rigidity,might also have velocity-dependent changes in movement of their arms. Prior studies have suggestedthat the damping ratio, which was referred to as maximal swing excursion, may be a marker ofspasticity in the upper limb and damping coefficient may be a marker of rigidity limb [77, 138].Moreover, by measuring angular position of elbow joint and passive torque, another study showedthat damping coefficient (decay rate) is a good measure for rigidity [117]. Therefore, in the future,by applying pendulum test to the arm using this model, it can be deduced whether the hypertonicityin rigidity is differentiated from spasticity where increase in muscle tone is velocity dependent andmore pronounced in extension than flexion. This study used both flexion and extension cycle. Ithad been shown that the viscous component of stiffness had a better correlation with UPDRS thanthe elastic one [117]. This is in alignment with our finding with regards to the viscous term ‘b’,which was significantly higher among PD subjects off-medication.This is a simple test based on sensors that are mounted on the wrist and no huge biomechanicaldevices are needed. It is safe, reliable, and repeatable. Unlike the common assessment of rigiditybased on UPDRS [42], which is subjective, this is an objective measure. This is another application564.5. DiscussionNormal PD off med PD on med00.511.522.533.544.555.5Decay Rate−−dominant armFigure 4.12: The decay rates for dominant arm.of second-order LTI models to assess motor performance in Parkinson’s disease.574.5. Discussion0.5 1 1.5 2 2.5 3024681012Effective Decay Rate−−averagedTotal Rigidityy = −1.0206 + 4.1331xR2 = 0.6120Figure 4.13: The correlation of effective decay rate and rigidity.58Chapter 5Assessment of Cognitive InflexibilityThrough Mode Detection and DelayIn addition to motor performance limitations, a key cognitive deficit seen in PD is that of “cognitiveinflexibility,” often assessed by the Wisconsin Card Sorting Test (WCST) [5, 103] and other tests. Inthis chapter, we explore the idea of sudden change in manual pursuit tracking because continuous-time models can provide insight into compensatory mechanisms in PD.We use a novel test to evaluate both motor performance and cognitive inflexibility – manualpursuit tracking with sudden changes in the tracking dynamics. We evaluate the performance ofnormal and PD subjects by evaluating 1) the dynamical parameters describing tracking dynamicswithin a task, and 2) detection of switching between tasks, as well as the corresponding delay indetection. Our approach builds on an input-output view of tracking tasks to create dynamicalmodels in which the human is essentially a black-box [2, 31, 98, 146]. In previous work [8], theauthors demonstrated statistically significant sensitivity in the parameters of second-order lineartime-invariant (LTI) dynamical models that had been fitted to manual pursuit tracking data. Al-though a standard measure of tracking performance, RMS error, was not significantly differentbetween normal and PD subjects, differences in damping ratio and natural frequency reflected themotor effects of medication in PD. This chapter extends those results by evaluating the ability ofsubjects to modify their performance in response to sudden changes in task dynamics.There are three possible “modes” of pursuit tracking in our experiment, based on whether thevisual feedback of the actual tracking error is attenuated, exaggerated, or unchanged. Each mode ismodeled as a second-order continuous system with white noise. Hence with sudden mode changes,the problem of pursuit tracking can be modeled as one of hybrid estimation and control. Thehuman must regulate an error trajectory whose dynamics may suddenly change. Since the modechange is not explicitly announced to the subject, the subject must also estimate the current modein order to maximize tracking performance. We draw upon multiple model adaptive estimation(MMAE) [11, 17, 55, 108] to identify a) whether the subject can adapt to the sudden change indynamics, and b) the time it takes to do so – that is, to identify a) the current mode, and b) thedelay in detecting the correct mode after a mode change. The perceived unpredictability of the595.1. Experimental Setupf2f1ErrorTargetCursor‘Glass rod’Figure 5.1: Experimental setup. The target trajectory is u(t) = sin(f1t) + sin(f2t). Users areinstructed to keep the cursor y(t) level with the target. The error u(t)− y(t) is scaledby 0.3 in ‘Better’ mode, by 2.0 in ‘Worse’ mode, and unchanged in ‘Normal’ mode.continuous-time tracking target is handled through the use of a Kalman filter, as in [9, 105, 146].The novelty of this work relates to a) the application of hybrid estimation techniques to amanual pursuit tracking task with multiple modes, b) a statistically significant assessment of thepresence and delay of mode detection in PD compared to normal subjects, and c) a statisticalassessment of the damping ratio and natural frequency of tracking performance in non-switchingand switching PD subjects.5.1 Experimental SetupFourteen PD subjects (on and off L-dopa medication) with clinically diagnosed, mild to moderatePD and ten healthy, age-matched subjects without active neurological disorders were recruited forthis study at the Pacific Parkinson’s Research Centre at the University of British Columbia atVancouver, Canada, after first providing informed consent (a full description of the experimentalsetup can be found in [8]). Subjects were asked to perform a tracking task by using a joystick inresponse to visual stimuli displayed on a computer screen, as shown in Figure 5.1. A horizontal“glass rod” connecting two boxes (each 60mm × 45mm) was shown on the display, where the boxon the left (Target) oscillated in the vertical direction at a linear combination of two frequencies (f1and f2), thus giving it a smooth but fairly complex appearing motion. Subjects were instructed tomove the box on the right (Cursor) by using the joystick so that the glass rod remained horizontalat all times. All subjects practiced for 5 – 10 minutes, during which time f1 and f2 were determinedfor each subject’s hand to maintain an error rate between 60 – 70% of the time. The individuallydetermined frequencies, f1 and f2, were then held constant throughout the rest of the study.PD subjects performed the task once after an overnight withdrawal (minimum of 12 hours sincetheir last dose of L-dopa, minimum of 18 hours since the last dose of dopamine agonists) of their605.2. Model Structureanti-Parkinson drugs and again one hour after admission of L-dopa.Part 1: Over a single 90-second interval, a sequence of three separate tracking tasks was per-formed, with a short delay (5–10 seconds) between each task to mark its end. In each task, thevisual feedback of the actual tracking errors was either amplified, attenuated or unaltered (but didnot switch between the three options). In the ‘Normal’ task, the vertical distance between thetarget and cursor displayed on the monitor reflected the true error generated by the subject. In the‘Better’ task, this distance was artificially reduced on the computer screen to 30% of the true error.In practice, the attenuation essentially made the tracking error better than expected. Finally, inthe ‘Worse’ task, the distance between the target and the cursor was artificially doubled, makingthe tracking error worse than expected. Subjects performed eight sets of the 90-second intervals(i.e., a total of 8× 3 tasks). The work upon which this chapter is built ([8]) focused solely on thispart of the experiment.Part 2: The same sequences of three different tasks was again performed over a single 90-secondinterval, but without a delay between tasks. The subject was not provided with any additionalsignal that might indicate that the task had changed. In effect, two unenunciated mode switchesoccurred in every 90-second interval. With a 10-second pause at the start of each interval, thefirst task lasted 20 seconds, and the remaining two tasks each lasted 30 seconds. This pattern wasrepeated eight times, resulting in a total of 8× 2 mode switches. A total of 4 sequences were eachtested twice: ‘Normal-Better-Worse’, ‘Worse-Normal-Better’, ‘Better-Worse-Normal’, and ‘Better-Normal-Worse’.The focus of previous work [8] was solely on the first part of the experiment. The models wereunderdamped (having complex conjugate eigenvalues) with distinction of natural frequency anddamping ratio significance across groups and tasks [8]. Different behaviors of a normal subject (asan example) is provided in Figure 5.2 in order to give an insight of the tracking errors they makein different tasks. As can be seen, as the manipulated error shown to them increases, they trackthe target better. The errors are relatively larger in corresponding tasks for PD patients [8].5.2 Model StructureFor modeling purposes, we consider the motor system as a black box, with the target trajectory asthe input, and the cursor trajectory as the output. Each subject’s tracking dynamics are modeledas discrete-time stochastic hybrid system H = (Q,X , R, f,Σ,U ,D, h) with discrete modes Q ={‘Better’, ‘Normal’, ‘Worse’} corresponding to each tracking task, continuous state x ∈ X ⊆ R2,transition function R as shown in Figure 5.3, f the set of LTI discrete-time dynamics fq : Q ×X ×U ×D → X indexed by the current mode q ∈ Q, discrete inputs Σ = {σBetter, σNormal, σWorse},615.2. Model Structure0 5 10 15 20 25 30−5000500Better0 5 10 15 20 25 30−5000500NormalTarget and cursor position [mm]0 5 10 15 20 25 30−5000500WorseTime [seconds]Figure 5.2: Typical normal subject’s tracking performance in different tasks. The solid blue linerepresents target position, and the dashed red line represents cursor position. As canbe seen, when the manipulated error shown to them increases, they track the targetworse. This figure will be replaced with a 3-mode hybrid system in Section 5.5.continuous input u ∈ U ⊆ R that represents the target position, continuous disturbance thatconsists of two mode-dependent i.i.d gaussian processes [wq, vq]T ∈ D ⊆ R2 with covariances Qqand Rq, respectively, and output map h : Q×X ×U ×D → Y ⊆ R such that y ∈ Y ⊆ R representsthe actual cursor position. The reset map is the identity. To model the dynamics in each mode, webuild upon previous work in manual pursuit tracking [2, 98], which characterized human trackingperformance to a sinusoidal input as a second-order LTI system. We correspondingly model thedynamics in each mode (e.g., for each of ‘Better’, ‘Normal’, and ‘Worse’ tracking tasks) as asecond-order LTI system in canonical observer form. In contrast to previous work [8], white noiseis explicitly included in the dynamics. Consider a discrete-time Markov jump linear system withthree modes: ‘Better’, ‘Normal’, ‘Worse’. In mode q the dynamics are given byx[k + 1] = Aqx[k] +Bqu[k] + Fqwq[k], Qq = E(wq[k]wq[k]T )y[k] = Cx[k] +Dqu[k] + vq[k], Rq = E(vq[k]vq[k]T )(5.1)625.2. Model StructureBetterNormalWorseσWorsey[k] = x1[k] + DNormalu[k] + vNormal[k]y[k] = x1[k] + DWorseu[k] + vWorse[k]y[k] = x1[k] + DBetteru[k] + vBetter[k]x[k + 1] = ABetterx[k] + BBetteru[k] + FBetterwBetter[k]x[k + 1] = ANormalx[k] + BNormalu[k] + FNormalwNormal[k]x[k + 1] = AWorsex[k] + BWorseu[k] + FworsewWorse[k]σBetterσBetterσNormalσNormalσWorseFigure 5.3: Hybrid system H that represents the switched pursuit tracking task. Reference inputtrajectory u represents the target position, and output trajectory y represents the cursorposition. Modes ‘Better’, ‘Normal’, and ‘Worse’ refer to the attenuated, unchanged, orexaggerated error reflected in the manipulation of the cursor position.For a given task, u[k] is the target position (the trajectory subjects are asked to track) and y[k] isthe actual cursor position at time k. State space matrices Aq, Bq, Fq, Dq and noise variances Qq, Rqcan be extracted from these two sets of values for each mode [80, 128], and C =[1 0]in allmodes. Hence with 10 Normal subjects and 14 PD subjects, and 3 tasks per subject, 10×3+14×3 =72 stochastic LTI models must be identified.The second order stochastic LTI formulation was chosen amongst different modeling frameworksbecause it has several advantages. 1) Previous work has established the Kalman filter as a reasonablemodel of the cerebellum [88]. 2) Previous work has established second-order linear systems aseffective models of manual pursuit tracking [2, 98, 129]. 3) Using higher-order LTI systems ornonlinear systems with extended Kalman filters improves model accuracy (measured accordingto accuracy in reproducing an independent input-output data set) only marginally with morecomplicated model formulations. 4) The computational cost of switching becomes much higherwith higher-order LTI systems or nonlinear systems.635.3. Model Identification+TargetPositionCursorPositionLTI SystemOutputGeneratorARX-VarianceyˆyBlack-boxSystemIdentificationGrey-boxSystemIdentificationu Grey-boxmodelmodelFigure 5.4: Block diagram of model identification process.5.3 Model IdentificationWe use Matlab’s System Identification toolbox [80] to derive stochastic LTI modes in observercanonical form, to ensure the same state representation x for each task and subject. A grey-boxmodel identification scheme (shown in Figure 5.4) was used since it handles both process noise andmeasurement noise [128], required to compute the Kalman gains necessary for state estimation.Deterministic second-order ARX models calculated in [8] via black-box LTI system identificationwere used to first estimate the measurement noise variance Rq, by computing the difference betweenpredicted and actual outputs. Then grey-box identification was used to determine unknown systemmatrices Aq, Bq, Fq, Dq and process noise variance Qq.Among the preprocessing steps mentioned in Chapter 2 the following ones are used:• Signal conditioning: We low-pass filtered the data at 8Hz according to Chapter 2 [69]. Inaddition, since each task had been carried out eight times in part 1 of the experiment, wemerged their input-output data for system identification.• Dealing with apparatus saturation The experiment was done using a joystick. Dueto apparatus saturation at peaks and/or valleys of the input signal (See Section 2.1.2), wedecided to de-saturate the input-output error according to equation (2.2).Four sets of input-output data were used to create a single grey-box model, and the remainingfour sets of input-output data were used to validate the model [81]. Validation was accomplishedby using Matlab’s compare function (as in [8]), which determines how well the model’s predictedoutput represents the actual output in an independent data set. Across all subjects, the averagegrey-box model accuracy was 80.05% ± 10.27%. Numerical issues due to the limited spectrum ofthe input signal are mitigated by the amount of data gathered.The limited frequency content of the input signal precludes identifiability, however we are onlyinterested in behavior of the model at frequencies in the immediate neighborhood of the frequenciesin the input signal [2, 98], and the process noises are relatively small. With all models restricted645.4. Detecting Mode Changesto be second-order (and hence limited in the parameter space), numerical identification techniqueswere robust to several initializations, showing local robustness under the iterative methods used. Inaddition, reasonable loss function values as well as prediction errors were obtained for each model.Hence we believe that validity of comparisons across groups and tasks based on the numericallyidentified models is maintained.5.4 Detecting Mode ChangesWe consider the problem in which a sequence of input-output data is provided, along with dynamics(5.1) that correspond to a known set of modes, to determine which mode mode is most likely. Thisis essentially the problem all subjects face when subjected to the merged sequence experiment, inwhich the current task is not enunciated and suddenly changes at seemingly arbitrary intervals. Byreverse engineering the “best” mode sequence for an input-output data set, we attempt to recreateeach subjects’ assessment of the correct task sequence. That is, given the hybrid system withdynamics (5.1), and a known input/output sequence, we wish to determine the switching sequencewhich maximizes the likelihood that the estimated mode is the actual mode.5.4.1 Multiple Model Adaptive EstimationMultiple Model Adaptive Estimation (MMAE) identifies the most likely mode by comparing resid-uals amongst different modes simultaneously, and selecting the one which is most likely to generatethe current data. Residuals associated with each mode are created through a Kalman filter, whichcreates an optimal state estimate xˆ. For a standard Kalman filter in mode q [20, 35],Prediction:xˆq[k|k − 1] = Aqxˆq[k − 1] +Bqu[k − 1]Pq[k|k − 1] = AqPq[k − 1]ATq + FqQq[k − 1]F TqUpdate:rq[k] = y[k]− yˆq[k]Sq[k] = CPq[k|k − 1]CT +Rq[k]Kq[k] = Pq[k|k − 1]CTS−1q [k]xˆq[k] = xˆq[k|k − 1] +Kq[k]rq[k]Pq[k] = (I −Kq[k]C)Pq[k|k − 1](5.2)655.4. Detecting Mode ChangesTargetPositionCursorPosition][kyHumanSubject][kuKF‘Better’KF‘Normal’KF‘Worse’][kw ][kv][ˆ kxBetter][krWorse][krNormal][kS Normal][kSWorse][ˆ kxNormal][ˆ kxWorse][kS Better][krBetterPosteriorProbabilityEvaluator(PPE)MAX][kWNormal][kWWorse][kWBetterFigure 5.5: Schematic of the MMAE algorithm. Kalman filters (KF) associated with each modegenerate a state estimate xˆq[k] that is used to create a residual for that mode. Residualsrq[k], Sq[k] are weighted and compared in the posterior probability evaluator (PPE),which generates a likelihood function Wq[k] for each mode. The most likely mode µˆ[k]is then selected.with error covariance Pq[k], residual rq[k], and residual covariance Sq[k], the predicted output isyˆq[k] = Cxˆq[k] +Dqu[k]. (5.3)The optimal Kalman gain is Kq[k].MMAE algorithms vary in the calculation of the modes’ likelihood function [17, 45, 56, 84, 91,97]; we use the MMAE described in [45], that models discrete events as a Markov process withconstant transition probabilities. While various methods of likelihood update were evaluated, wechose this method for its successful detection of the mode switches that best matched the actualmode changes in normal subjects.A schematic of the MMAE algorithm is shown in Figure 5.5. Residuals for each mode areweighted adaptively to create a likelihood function for each mode, with the residual associated withmode q equal to the difference between the actual output and predicted output using dynamics665.4. Detecting Mode Changesassociated with mode q. A posterior probability evaluator (PPE) generates the likelihood functionΛq[k] = 1√2piSq[k]· exp{− r2q [k]2Sq[k]}(5.4)that determines the probability Vq[k] that mode q is the true mode µ[k]. Since the true mode doesnot change very frequently, we useVq[k] = Wq[k − 1]Wq[k] = Λq [k]Vq [k]∑3j=1 Λj [k]Vj [k](5.5)with final probability Wq[k], and Wq[0] = 13 (since there are 3 modes of equal likelihood). Finally,the mode estimateµˆ[k] = argmaxq Wq[k] (5.6)is the mode with the highest likelihood Wq[k].5.4.2 Implementing the MMAENormally, noise covariances are determined by heuristic, ad hoc approaches, which leads to theclassical “tuning of the filter” problem [83]. For the LTI system in equation (5.1), process noisevariances estimated by the system identification algorithms are relatively large since Matlab’s grey-box identification algorithm assumes that the system is driven by a known measurement noise.In an attempt to reproduce the output with a measurement noise of known covariance, a largeprocess noise is associated to the state equation. This creates residual covariances in equation (5.2)that are incapable of representing how well each model represents the actual output. However, thevalues for process noise covariance are unreliable due to their high variance. Hence, we initializethe Kalman filters with process noise covariances approximately five orders of magnitude lowerthan those calculated by the identification algorithm, in order to accommodate higher magnitudemeasurement noise.Lastly, we note that since all dynamics are represented in observer canonical form, withoutprocess and measurement noise, it would be impossible to uniquely reconstruct the state trajectoryfor the autonomous system. That is, the modes would be indistinguishable since for all mode pairsp, q ∈ Q, the observability rank test fails: rank([Op Oq]) ≤ 2n [140]. However, for any transitionfrom mode p to mode q, the variance of the continuous state will change, and hence so will thevariance of the output.675.5. ResultsTable 5.1: Dynamics for subjects shown in Figures 5.6, 5.7, and 5.8, with Qq = diag(αq, 0).Type Fig Matrices for mode predictionABetter =[ 1.998 1−0.999 0], BBetter =[ 0.002−0.003],DBetter = 0.006, αBetter = 2.2e7, RBetter = 1.1e4Normal 5.6 ANormal =[ 1.994 1−0.994 0], BNormal =[ 0.005−0.006],DNormal = 0.004, αNormal = 7.1e7, RNormal = 2.8e3AWorse =[ 1.917 1−0.925 0], BWorse =[ 0.039−0.031],DWorse = 0.035, αWorse = 1.5e9, RWorse = 3.9e3ABetter =[ 1.959 1−0.960 0], BBetter =[ 0.018−0.017],DBetter = 0.018, αBetter = 2.8e8, RBetter = 3.7e3PD 5.7 ANormal =[ 1.967 1−0.969 0], BNormal =[ −0.0440.046],off med DNormal = 0.060, αNormal = 1.2e6, RNormal = 2.6e3AWorse =[ 1.874 1−0.874 0], BWorse =[ 0.115−0.114],DWorse = 0.134, αWorse = 4.3e7, RWorse = 1.9e3ABetter =[ 1.991 1−0.993 0]BBetter =[ −0.0150.017],DBetter = −0.001, αBetter = 7.1e7, RBetter = 5.2e3PD 5.8 ANormal =[ 1.991 1−0.992 0], BNormal =[ −0.0040.005],on med DNormal = −0.005, αNormal = 2.9e8, RNormal = 1.8e3AWorse =[ 1.907 1−0.914 0], BWorse =[ 0.056−0.050],DWorse = 0.055, αWorse = 2.2e8, RWorse = 1.4e35.5 ResultsAmong the four switching sequences (‘Normal-Better-Worse’, ‘Worse-Normal-Better’, ‘Better-Worse-Normal’, ‘Better-Normal-Worse’), switching between ‘Better’ and ‘Worse’ tasks was the most ob-vious across all subjects. While switching between ‘Better’ and ‘Normal’ or between ‘Normal’ and‘Worse’ was also evident in some subjects, we focus on switching that occurred from ‘Better’ to‘Worse’ modes. Typical tracking performance for the ‘Normal-Better-Worse’ sequence is shown fora normal subject in Figure 5.6, for a PD subject off medication in Figure 5.7, and for a PD subjecton medication in Figure 5.8. The top part of each of these figures shows the target position and685.5. Results0 10 20 30 40 50 60 70 80−400−300−200−1000100200300400Normal Better WorseTarget and cursor position [mm]0 10 20 30 40 50 60 70 80024Time [seconds]NormalBetter WorseModeFigure 5.6: Typical normal subject. The solid blue line represents the target position, and thedashed red line represents the cursor position. The subject detects both switches, from‘Normal’ to ‘Better’ as well as from ‘Better’ to ‘Worse’. The switching delay is verysmall in comparison to the time scale of the entire ‘Normal-Better-Worse’ sequence.cursor position (e.g., input u and output y), and the bottom part of each of these figures shows theswitching sequence estimated according to equation (5.6). When switching occurs, large trackingerrors immediately result in the new mode. This can be seen at t = 20 and t = 50 in Figures 5.6, 5.7,and 5.8. System matrices used to estimate the current mode for the particular subjects in thesefigures are listed in Table 5.1.In the top portion of each figure, the input-output trajectory shows the expected poorer perfor-mance in the ‘Worse’ task as compared to the ‘Normal’ task, and the expected improved trackingperformance in the ‘Better’ task. The bottom portion of the figure shows the most likely modesequence generated by the MMAE algorithm. The normal subject shown in Figure 5.6, like all othernormal subjects in the study, successfully detected the task change between ‘Better’ and ‘Worse’.Figure 5.7 depicts tracking performance of a typical PD subject off medication who does notdetect the task change, and Figure 5.8 depicts the tracking performance of a typical PD subjecton medication who does detect the task change. When a switch occurs, large tracking errors lead695.5. Results0 10 20 30 40 50 60 70 80−400−300−200−1000100200300400Normal Better WorseTarget and cursor position [mm]0 10 20 30 40 50 60 70 80024Time [seconds]Better Better BetterModeFigure 5.7: Typical PD subject off-medication. The solid blue line represents the target position,and the dashed red line represents the cursor position. The subject is unable to detectboth switches, from ‘Normal’ to ‘Better’ as well as from ‘Better’ to ‘Worse’. We namethis type, a non-switching off medication PD subject.to a decrease in the likelihood of the current mode (5.4), prompting a change in the mode withthe maximum likelihood. This is most noticeable in Figures 5.6, 5.7, and 5.8 at times t = 20s andt = 50s.Analysis of the performance characteristics between tracking tasks is examined in detail in [8].Across groups, subjects showed difficulties in tracking the target at extrema. While normal subjectshad lower RMS error on average than PD subjects, the difference was not statistically significant.5.5.1 Switching DetectionSince switching between ‘Better’ and ‘Worse’ tasks was the most evident among the four switchingsequences (‘Normal-Better-Worse’, ‘Worse-Normal-Better’, ‘Better-Worse-Normal’, ‘Better-Normal-Worse’) across all subjects, it would be the main focus of this section. As described in Table 5.2,the sudden change between ‘Better’ and ‘Worse’ tasks was detected by all 10 normal subjects inthe ‘Normal-Better-Worse’ sequence. But only 5 out of 14 PD subjects off medication detected the705.5. Results0 10 20 30 40 50 60 70 80−400−300−200−1000100200300400Normal Better WorseTarget and cursor position [mm]0 10 20 30 40 50 60 70 80024Time [seconds]Better Better WorseModeFigure 5.8: Typical PD subject on-medication. the solid blue line represents target position, andthe dashed red line represents cursor position. This subject does detect the switch att = 50s from ‘Better’ to ‘Worse’.switch, while 7 out of same 14 PD subjects on medication detected the switch. Not unexpectedly,proportionally far fewer PD subjects detected the switch than did normal subjects. Further, theeffect of L-dopa is to improve PD subjects’ ability to detect a mode change, and hence make theirperformance slightly more similar to that of normal subjects.The transition from ‘Better’ to ‘Worse’ was detected by only 4 out of 14 PD subjects offmedication, and 9 PD subjects on medication in the ‘Better-Worse-Normal’ sequence. However,all 10 normal subjects detected the switch in this sequence as well. There are different number ofin PD subjects who detected the switch off and on medication. This may partly be due to fatigue[150], since more subjects improve after medication in the sequence in which the transition from‘Better’ to ‘Worse’ occurs earlier in the switching sequence.These results show difficulty for PD subjects in adapting to a sudden change, which are con-sistent with the previous studies [5, 74, 79, 103]. This difficulty may be attributed to the basalganglia, which plays a major role in regulating motor commands for switching from one task toanother [103]. We also note that mode switches between ‘Normal’ and ‘Better’ tasks and between715.5. ResultsTable 5.2: Switching detection and delays between ‘Better’ and ‘Worse’ modes in the ‘Normal-Better-Worse’ sequence.Normal PD PD Statisticaloff med on med SignificanceNumber of subjectswho detected the 10/10 5/14 7/14mode changeMean delay inswitching detection 3.77 5.0 4.14 p = 0.07,[sampling units] ANOVAVariance of delay in p = 0.88,switching detection 0.83 0.70 0.89 VARTEST‘Worse’ and ‘Normal’ tasks were not detected by many normal subjects or PD subjects. Maybebecause these tasks are not as distinctly different as ‘Better’ and ‘Worse’.To determine whether the failure to detect switching was due to biological phenomenon ormerely to poorly tuned switching algorithms, we also implemented higher-order LTI models (notshown). Although these models had higher accuracies, and hence were more likely to allow fordifferentiation between modes in the MMAE algorithm, no additional mode changes were detected.Hence we conclude that the failure to detect switching is a feature of the underlying biology.5.5.2 Delay in Switching DetectionLatencies for the subjects capable of detecting mode change between ‘Better’ and ‘Worse’ in the‘Normal-Better-Worse’ sequence approached statistical significance across three groups (p = 0.07,ANOVA). The mean delay was 3.77 ± 0.83 time steps (each time step is 0.03 seconds) for normalsubjects, while for PD subjects on medication the mean delay was 4.14± 0.89 time steps, and forPD subjects off medication the mean delay was 5.0± 0.70 time steps (see Table 5.2). As expected,L-dopa decreased the amount of time required for PD patients to detect a mode change.A similar pattern appeared in the ‘Better-Worse-Normal’ sequence, although not statisticallysignificant (p = 0.3351, ANOVA). For normal subjects, the mean delay was 4.0± 1.41 time steps,while for PD subjects on medication the mean delay was 4.14±1.34 time steps, and for PD subjectsoff medication the mean delay was 6.5±3.69 time steps. So, the delays have been more homogenized(having less variance) after medication. The delay variance was significant across groups in thisblock (p = 0.05, variance test), which could be related to the homogenizing effect of the medication.While previous studies have shown that reaction time is delayed in PD subjects during switchingexperiments [41, 72], a key result of this work is that many PD subjects simply do not seem to725.5. Results0 10 20 30 40 50 60 70 80−400−300−200−1000100200300400Normal Better WorseTarget and cursor position [mm]0 10 20 30 40 50 60 70 80024Time [seconds]Better Better WorseModeFigure 5.9: Typical PD subject off-medication. The solid blue line represents the target position,and the dashed red line represents the cursor position. The subject is able to detectthe switch from ‘Better’ to ‘Worse’. We name this type, a switching off medication PDsubject.detect the switch from ‘Better’ to ‘Worse’ tasks.5.5.3 Estimation in PDThe pattern observed in delays made us have a closer look into ‘normal’ and ‘PD subjects offmedication’ groups. For PD subjects off medication in the ‘Normal-Better-Worse’ sequence, themean delay was 5.0 ± 0.70 time steps, while it was 3.77 ± 0.83 time steps for normal subjects(See Section 5.5.2). This difference was statistically significant across these ‘two’ groups (p =0.017, ANOVA). Figure 5.10 provides information about latencies for normals and PD subjects offmedication who detected the switching. Figure 5.9 depicts the tracking performance of a switchingPD subject off medication (who does detect the task change).We further investigated the differences in tracking behavior between those PD subjects whodetected the task change and those who did not. The estimator dynamics in mode q can be735.5. ResultsNormal PD off−medication33.544.555.56Time stepsFigure 5.10: The delay in switching detection for different groups in ‘Normal-Better-Worse’ block.Each time step is 0.03 seconds.rewritten from equation (5.2) asxˆq[k] = (Aq −Kq[k]C)xˆq[k|k − 1] + (Bq −Kq[k]Dq)u[k − 1] +Kq[k]y[k]. (5.7)In steady-state, estimator dynamics have a damping ratio ζq and natural frequency ωn,q as deter-mined by the eigenvalues λ1, λ2 of (Aq −KqC), with Kq = limk→∞Kq[k] [20].ζq = −12 · λ1+λ2√λ1λ2ωn,q = +√λ1λ2(5.8)Table 5.3 shows the mean damping ratio and natural frequency for both normal and PD subjectsduring the ‘Better’ task of the ‘Normal-Better-Worse’ sequence (values for all subjects × modesare in Appendix A). Note as shown in Figures 5.11 and 5.12, PD subjects have damping ratiosand natural frequencies that are lower than the damping ratios and natural frequencies of nor-mal subjects. These patterns are statistically significant across groups (damping ratio p = 0.015,ANOVA; natural frequency p = 0.012, ANOVA). Our results are consistent with the results ofnoise-free models reported in [8] and other studies [72], demonstrating impaired reaction time inPD subjects.As expected, estimator dynamics are less damped in non-switching PD subjects (damping ratio0.0044±0.0030) than in normal subjects (damping ratio 0.0117±0.0100), indicating a more reactiveresponse in PD subjects who do not detect the task change, as compared to normal subjects.Switching PD subjects’ mean damping ratio (0.0110 ± 0.0070) falls between these two groups.745.5. ResultsTable 5.3: Mean damping ratio and natural frequency in estimator dynamics for normal and PDsubjects in ‘Better’ mode.Normal PD PD Statistical(switching) switching non-switching significanceDamping ratio 0.0171 0.0116 0.0044 p = 0.015,ANOVANatural frequency 38.90 38.76 38.60 p = 0.012,[rad/s] ANOVAEssentially, those PD subjects who are more similar to normal subjects in the ability of detectingthe sudden task change, also have estimator dynamics that are “closer” to that of normal subjects(as compared to PD subjects who cannot detect the task change).A similar pattern was also observed for natural frequency across groups (See Table 5.3). Thenatural frequency of estimator dynamics was 38.90±0.25 rad/sec for normal subjects, while switch-ing and non-switching PD subjects had natural frequencies of 38.76±0.17 rad/sec and 38.60±0.06rad/sec, respectively.In both Figures 5.11 and 5.12, considerable overlap exists between groups in both estimatordamping ratio, and estimator natural frequency. While the overall pattern is statistically significant,the overlap means that these values cannot be used diagnostically to discriminate PD subjectsfrom normals, or within PD subtypes. Lastly, no clear distinction between groups was seen in the‘Normal’ or ‘Worse’ tasks.Previous work has attempted to develop a standardized approach in the clinical evaluation ofthe motor performance in PD patients [61, 115], possibly allowing the disease to be stratified intodifferent subtypes with different rates of progression. Our demonstration of impaired switchingand altered dynamics in a subset of PD subjects suggests that the proposed theoretical frameworkprovides a comprehensive way to examine the motor deficits that may be seen in PD.755.5. ResultsNormal switching PD non−switching PD00.0050.010.0150.020.0250.03Damping ratioFigure 5.11: Damping ratio of estimator dynamics in the ‘Better’ task in the ‘Normal-Better-Worse’sequence.765.5. ResultsNormal switching PD non−switching PD38.538.638.738.838.93939.139.2Natural frequencyFigure 5.12: Natural frequency of estimator dynamics in the ‘Better’ task in the ‘Normal-Better-Worse’ sequence.77Chapter 6Conclusions and Future WorkWe considered the problem of motor performance assessment and inference of mechanisms in thebrain in Parkinson’s disease. We proposed linear dynamical systems to finely characterize motorperformance, and for possible interpretations of neurological phenomena, such as compensatorymechanisms. While gross measurements of motor performance such as root mean square error,delay in tracking tasks, and maximum velocity cannot distinguish between qualitatively differenttypes of movement, dynamical model parameters (e.g., damping ratio and natural frequency) canmore readily describe those differences.We used an input-output view of tracking tasks to create second order LTI dynamical modelsin which the human is essentially a black-box. LTI model parameters can provide insights intobrain function in disease state, by quantifying dynamic motor performance. Prior to model iden-tification, a set of preprocessing techniques dependent on the particular experimental setup wasperformed. The pre-processing techniques have been carefully evaluated to appropriately addresspossible nonlinearities, non-stationarity, and other data irregularities.Manual tracking experiments were utilized for quantification of motor performance and biolog-ical interpretation of compensatory mechanisms in the brain. We demonstrated that the classicmodel of the BG is insufficient to explain the more benign prognosis of PDT compared with PDARand the motor system is also influenced by the cerebello-thalamo-cortical (CTC) circuit, in addi-tion to the BG. PD tremor led to better motor performance with higher damping ratios and thiscorrelation could be associated to compensatory mechanism. We speculated that augmented cere-bellar activity, in order to compensate for diseased BG circuits, results in enhanced motor trackingperformance.We used a simple test based on sensors that are mounted on the wrist and no huge biomechanicaldevices to measure rigidity through arm swing. It is safe, reliable, safety, and repeatable. Rigidityis commonly assessed based on UPDRS, which is a subjective tool used by physicians. We showedthat LTI models could be exploited based on modified pendulum test to come up with an objectivemeasure for rigidity.To address the issue of cognitive “inflexibility”, we explored the idea of sudden change inmanual pursuit tracking. Since the mode change is not explicitly announced to the subject, the786.1. Summary of Contributionssubject must also estimate the current mode in order to maximize tracking performance. Weevaluated the performance of normal and PD subjects by evaluating 1) the dynamical parametersdescribing tracking dynamics within a task, and 2) detection of switching between tasks, as well asthe corresponding delay in detection. We demonstrated impaired switching and altered dynamics ina subset of PD subjects that suggests the proposed theoretical framework provides a comprehensiveway to examine the motor deficits in PD.6.1 Summary of ContributionsWe summarize our contributions as follows:• Selection of pre-processing and system identification techniques for three experimental setupsincluding manual tracking experiments• Control theoretic framework for tremor as a compensatory mechanism that is advantageousin some tracking tasks• Sensor-based assessment of rigidity via decay rate• The application of hybrid estimation techniques to a manual pursuit tracking task with mul-tiple modes to assess cognitive inflexibility6.2 Future Research DirectionsThere are several directions for future work and further investigations.6.2.1 The Effect of Galvanic Vestibular StimulationHypothesisWe assessed motor performance as well as cognitive inflexibility in Chapter 5 and we explored theeffect of Levodopa on tracking performance as well as the ability in switching between differenttasks. While Levodopa is known to be the “gold standard” treatment for PD, it might lead to someside effects such as dyskinesia and psychiatric disorders [123]. Thus, non-invasive brain stimulationtechniques are of great interest to researchers because of their tolerability and safety; for PD aswell as other neurological disorders.Non-invasive brain stimulation techniques rely on electromagnetic basis to transcranially modifybrain activity [121]. Transcranial direct current stimulation (tDCS) and transcranial magnetic796.2. Future Research Directionsstimulation (TMS) are the two most famous methods. In tDCS, an electrical current is appliedto a region of scalp to depolarize neurons via decreasing GABA transmission, and thus improvespontaneous firing [85]. TMS involves applying a transient magnetic field to the brain througha coil, which causes reversible shifts in cortical excitability in focalized regions [83]. tDCS andTMS can particularly target the cortical-basal ganglia circuitry by primary motor cortex (M1)stimulation as well as prefrontal loops [116].Recent studies have shown that noisy galvanic vestibular stimulation (GVS) could be beneficialfor PD patients [100, 102]. GVS is a technique known to alter the firing rates of vestibular afferentnerves. It is fairly inexpensive compared to TMS and does not have the low focality drawbackfrom which tDCS is suffering [85, 92]. GVS affects cortical-vestibular regions without inducing anyelectrical activation in peripheral vestibular afferents [137]. A promising avenue of future study isto investigate the effect of noisy GVS on motor performance in a visuomotor tracking task in PD.Preliminary analysisSome pursuit tracking data have been collected at the Pacific Parkinson’s Research Centre at theUniversity of British Columbia from normal and PD subjects. Twelve PD subjects (10 males, 2females; mean age 61.4 ± 6.5 years; 11 right-handed, 1 left-handed) off L-dopa medication wererecruited in this study. Twelve healthy age-matched control subjects (4 males, 8 females; mean age58.3± 9.0 years; all right-handed) without neurological disorders participated in the study.Similar to Chapter 5, a sequences of three different tasks was performed over a single 90-secondinterval without a delay between tasks (See Section 5.1). The subject was not provided with anyadditional signal that might indicate that the task had changed. In effect, two unenunciated modeswitches occurred in every 90-second interval. Each task lasted 30 seconds. Unlike the three tasksin Chapter 5, only two different tasks were tested; The ‘Better’ task, in which the distance betweenthe target and the cursor was artificially reduced on the computer screen to 30% of the true error.In practice, the attenuation essentially made the tracking error better than expected. The ‘Worse’task, in which this distance was artificially doubled, making the tracking error worse than expected.A total of 4 sequences were each tested twice: ‘Worse-Better-Worse’, ‘Better-Worse-Better’,when the GVS was on or off. The sequence at which the GVS was on or off was not enunciatedto the subjects and did not have any particular order. Each sequence was followed by a break(30 s) to preclude a hysteretic effect carrying over to the next sequence. Subjects were allowed topractice tracking the target and using the joystick prior to the experiment as needed in at leastone practice sequence. This experiment would be a means to evaluate the effect of GVS on bothmotor performance and cognitive inflexibility – manual pursuit tracking with sudden changes inthe tracking dynamics.806.2. Future Research Directions1 1.2 1.4 1.6 1.8 2−0.100.10.20.30.40.50.60.70.8arrows pointing from GVSoff to GVSonDamping ratio−−Better taskFigure 6.1: The change in damping ratio going from the GVS off state to the GVS on state for thebetter task.Some initial work on these data have been done and black box model have been identified foreach subject and each task. We investigated the effect of GVS on model dynamic parameters suchas damping ratio and natural frequency as well as some static parameters such as RMS error andmaximum velocity. While RMS error, was not statistically significant between GVS on versus offstate, differences in damping ratio and natural frequency reflected the motor effects of GVS in PD.In the “Better” task, the damping ratio of PD patients with GVS on was significantly higherthan that of PD patients with GVS off (p = 0.0072, ttest). These improvements are depicted inFigure 6.1. No significant change in damping ratio was observed for normal subjects. The maximumvelocity of PD patients with GVS on was also higher than that of PD patients with GVS off (p =0.0309, ttest). Interestingly, similar effect of GVS was seen in normal subjects (p = 0.0270, ttest).A summary of all significant changes in parameters are provided in Table 6.1. The values that arelisted as N.S in the table had no statistical significance.These results are consistent with previous studies showing the positive effect of GVS on motortracking data [100, 102]. Galvanic vestibular stimulation has caused the tracking response tobe faster and with a higher damping ratio. These tests would help us to quantify the effect ofnon-invasive brain stimulations to some degree and have a better understanding of their effect onmotor performance. Yet, GVS has cognitive effects as well [127] and more analysis has to be doneto evaluate cognitive inflexibility through hybrid switching algorithms. The switching algorithm816.2. Future Research DirectionsTable 6.1: Static and dynamic model parameters for PD subjects off medication in ‘Better’ and‘Worse’ tasks. All of statistical significances are computed using paired ttest.Better WorseDamping Ratio p = 0.0072 p = 0.0085Decay Rate N.S p = 0.0363Settling Time N.S p = 0.0434Max. Velocity p = 0.0309 p = 0.0335based on MMAE has led to high covariance matrices for Kalman filter observer and we are tryingto come up with a more robust algorithm for detecting the switches if they occur.6.2.2 Reward Mechanisms in Forward ModelsAt a lower level, it would be good to come up with some models, which represents the feedforwardand feedback processes in the brain as well as the tracking dynamics. A model predictive control(MPC) approach seems to be a valid choice according to reward/punishment mechanisms [118] thathappens in the brain. Shadmehr has investigated some qualitative measures of this approach [118].Some other researchers (e.g., [18]) demonstrated that PD patients (off-med) have a very specificdeficit in learning from reward, but are normal at learning from punishment. In contrast, whenthey are placed on dopamine agonists this learning sensitivity reverses: they have impairments inlearning from punishment, but are normal at learning from reward [18]. So, the reward/punishmentmechanism seems to be very important in forming a motor command and modeling this would helpus better understand the brain impairments in this disease.In addition, “apathy” is a disorder of diminished motivation and interest, resulting in a decreasedgoal-oriented behavior [38]. It is a frequent symptom early on in the disease [40], and has asignificant impact on quality of life in PD [94]. We want to use control theory application inanalyzing apathy effect in tracking performance of PD subjects. That would help us quantifyingmotor performance that may provide a more direct measure of motor dysfunction in PD. Hence,it would be useful in exploring feedforward and feedback processes in the brain. We essentiallywant to reverse engineer what happens in the brain that control (energy) signal ‘u’ is generatedaccording to the actual tracking error ‘e’ visually seen by the subject (See Figure 6.2). This couldbe an indicator of effort in subjects assuming that all other parameters are fixed. The goal (reward)and effort (energy) to reach that goal are best formulated in the MPC framework. So, using MPCformulation for tracking data, we potentially are expecting less effort in PD patients with apathy.Let us assume that ‘e’ is the error that subjects see as the difference of what output theygenerate and what the desired trajectory is. Moreover, let the control signal ‘u’ be the effort that is826.2. Future Research DirectionsReferencetrajectoryARX model Ay = Buu+ - Optimization actualyerFigure 6.2: Pursuit tracking block diagram. ‘r’ is the reference trajectory. ‘u’ is the control signalgenerated by the brain. ‘e’ is the tracking error. ‘yactual’ is the actual cursor positionon the computer screen.transformed into a force applied to the joystick generating the output y that we see on the screen.Assume that this transformation is done through an ARX model. So, we can write:y = Gu+ f (6.1)in which, y is a vector of predicted outputs, G is a matrix, whose rows are formed from theARX model obtained for this plant, u is a vector of current and future inputs, and f is a vector ofprevious outputs. Considering Fi(z−1) as a polynomial with respect to the backward shift operatorz−1, we can write:y[k + 1|k]y[k + 2|k]...y[k +Np|k]=g0 0 . . . 0g1 g0 . . . 0... ... ... ...gNp−1 gNp−2 . . . g0×u[k]u[k + 1]...u[k +Np − 1]+F1(z−1)F2(z−1)...FNp(z−1)y[k] (6.2)The optimization problem in MPC is to minimize the cost function (6.3) subject to stateevolution constraints mentioned in equation (6.1) over ‘u’ in a prediction horizon (Np) and coming836.2. Future Research Directionsup with vector ‘u’.J =Np∑i=1((y[k + i|k]− r[k + i])TQ(y[k + i|k]− r[k + i]) + u[i]TRu[i]) (6.3)In the cost function defined above, y[k] is the output (that is a function of states, which in-cludes cursor position and its velocity for our system), u[k] is the control effort, and Q,R are theweighting matrices associated to tracking error and control effort. Np is the horizon, until whichthe minimization is performed at each step (i.e., prediction horizon). The optimization problem isthen defined as:minu∈U⊂RJ s.t. y[k + i|k] = Gi(z)u[k] + Fi(z−1) ∗ y[k], i = 1, . . . , Np (6.4)However, only the first element of the vector ‘u’ is applied to the system and the optimizationproblem will be solved for a shifted horizon in the next step.Now, if the transfer function introduced in equation (6.1) is not known, it could be identifiedthrough recursive least square (RLS) method. This process is done online and we artificially assumeeach u applied to the orange block in Figure 6.2 would yield output y (reverse engineering). Byidentifying the parameters of this ARX model, optimal u could be calculated in each step. Thismay seem a bit unreasonable to generate optimal control signals with the identified parametersthat have not converged yet. However, RLS is a fast identification method and converges in lessthan 5–6 time steps and hence, it would yield control signals using identified parameters, which arevery close to true parameters.We would like to find the control signals ‘u’ at each time step over each task experiment andcompare them across groups. Previously, second order LTI models turned out to work well for thewhole system with feedback (i.e., the dashed rectangle in Figure 6.2). We could assume secondorder LTI model for the open loop plant (i.e., the orange box shown in Figure 6.2) as well.The novelty of this work is to ‘quantify’ the energy required to form a motor command inpursuit tracking tasks, while previous work (e.g., [118]) have qualitatively investigated that. Thiswould help to quantitatively assess apathy, which is related to reward/punishment mechanism toform motor commands in the brain.6.2.3 Telemonitoring Parkinson’s DiseaseUnified Parkinson’s disease rating scale (UPDRS) is often used to monitor disease severity, whichis a subjective test and requires patient’s physical presence at the clinic. Thus, monitoring thesymptoms is a costly and time-consuming task which is susceptible to potential errors. A great846.2. Future Research Directionsavenue of interest is to come up with a portable device or a remote test that can monitor diseaseseverity from afar [49, 136]. Most of these tests are based on speech and voice assessment andempirical tests that are subjective as well. Hence, developing a safe, reliable, repeatable, andobjective measure for telemonitoring the disease progression that can be self-administered quicklyseems necessary.Ideally, it would have been good to have a simple (not using huge biomechanical devices) andportable device that can measure disease severity based on motor performance. In Chapter 4we came up with a modified pendulum test that objectively measures rigidity as a component ofUPDRS based on a motor performance parameter (i.e., damping ratio). Ultimately, if tremor,bradykinesia and other components of UPDRS could be assessed through a statistical mappingof motor parameters, it’s worth to monitor the disease via those objective measures. Thus, “acombination of motor behavioral parameters” could be a good candidate to help elucidate diseaseseverity from afar and with simple tests. Some studies have been conducted to quantify bradykinesiaor tremor [112, 130] which include mounting miniature gyroscopes for hours on the forearm. Despitetheir relatively good accuracy, the period of time wearing these accessories is irritating for PDpatients.Finding this objective measure would be of great value for patients that are reluctant to makefrequent visits to the clinic. It could also be a means to future research areas that are attemptingto come up with novel treatments of Parkinson’s disease. Because those studies usually requirehigh-frequency, and remote samples from a large populations.85Bibliography[1] A. Abdel-Malek, C. Markham, P. Marmarelis, and V. Marmarelis. Quantifying deficienciesassociated with parkinson’s disease by use of time-series analysis. Electroencephalography andClinical Neurophysiology, 69(1):24–33, 1988.[2] A. Abdel-Malek and V. Z. Marmarelis. Modeling of task-dependent characteristics of hu-man operator dynamics pursuit manual tracking. IEEE Transactions on Systems, Man, andCybernetics, 18(1):163–172, January/February 1988.[3] J.-C. Aguero, C.R. Rojas, and G.C. Goodwin. Fundamental limitations on the accuracy ofmimo linear models obtained by pem for systems operating in open loop. In Decision andControl, 2009 held jointly with the 2009 28th Chinese Control Conference. CDC/CCC 2009.Proceedings of the 48th IEEE Conference on, pages 482–487, Dec 2009.[4] R. L. Albin, A. B. Young, and J. B. Penney. The functional anatomy of basal gangliadisorders. Trends in Neuroscience, 12(10):366–375, 1989.[5] A. Alevriadou, Z. Katsarou, S. Bostantjopoulou, G. Kiosseoglou, and G. Mentenopoulos.Wisconsin Card Sorting Test variables in relation to motor symptoms in Parkinson’s disease.Percept Mot Skills, 89(3):824–830, 1999.[6] D. P. Allen, J. R. Playfer, N. M. Aly, P. Duffey, A. Heald, S. L. Smith, and D. M. Halliday. Onthe use of low-cost computer peripherals for the assessment of motor dysfunction in Parkin-son’s disease - quantification of bradykinesia using target tracking tasks. IEEE Transactionson Neural Systems and Rehabilitation Engineering, 15(2):286–294, June 2007.[7] A. Ashoori, M.J. McKeown, and M.M.K. Oishi. Switched manual pursuit tracking to measuremotor performance in parkinson’s disease. Control Theory Applications, IET, 5(17):1970–1977, April 2011.[8] W. L. Au, N. Lei, Meeko M. K. Oishi, and M. J. McKeown. L-Dopa induces under-damped visually guided motor responses in parkinson’s disease. Experimental Brain Research,202(3):553–559, 2010.[9] R. J. Baddeley, H. A. Ingram, and R. C. Miall. System identification applied to a visuo-motor task: Near-optimal human performance in a noisy changing task. The Journal ofNeuroscience, 23(7):3066–3075, April 2003.[10] I. Bar-Gad and H. Bergman. Stepping out of the box: Information processing in the neuralnetworks of the basal ganglia. Current Opinion in Neurobiology, 11(6):689–695, 2001.86Bibliography[11] Y. Bar-Shalom, X. Rong Li, and Thiagalingam Kirubarajan. Estimation with Applicationsto Tracking and Navigation. Wiley, New York, 2001.[12] I Basu, D Graupe, D Tuninetti, P Shukla, KV Slavin, LV Metman, and DM Corcos. Patho-logical tremor prediction using surface electromyogram and acceleration: potential use in’on-off’ demand driven deep brain stimulator design. Journal of Neural Engineering, 10(3),June 2013.[13] Ishita Basu, Daniela Tuninetti, Daniel Graupe, and Konstantin V Slavin. A model for sim-ulating local field potential in the thalamus of essential tremor patient during deep brainstimulation. BMC Neuroscience, 13(Suppl 1):P40, 2012.[14] G. Becker, A. Muller, S. Braune, T. Buttner, R. Benecke, W. Greulich, W. Klein, G. Mark,J. Rieke, and R. Thumler. Early diagnosis of parkinson’s disease. Neurology, 249:40–48, 2002.[15] E. Bezard, A. R. Crossman, C. E. Gross, and J. M. Brotchie. Structures outside the basal gan-glia may compensate for dopamine loss in the presymptomatic stages of parkinson’s disease.The Journal of the Federation of American Societies for Experimental Biology, 15(6):1091–1094, 2001.[16] S. J. Blakemore, C. D. Frith, and D. M. Wolpert. The cerebellum is involved in predictingthe sensory consequences of action. Neuroreport, 12(9):1879–1884, 2001.[17] H. A. P. Blom and Y. Bar-Shalom. The interacting multiple model algorithm for systemswith Markovian switching coefficients. IEEE Transactions on Automatic Control, 33(8):780– 783, Aug. 1988.[18] N. Bodi, S. Keri, H. Nagy, A. Moustafa, C. E. Myers, N. Daw, G. Dibo, A. Takats, D. Bereczki,and M. A. Gluck. Reward-learning and the novelty-seeking personality: a between- andwithin-subjects study of the effects of dopamine agonists on young Parkinson’s patients.Brain, 132.[19] B. R. Brewer, S. Pradhan, G. Carvell, and A. Delitto. Application of modified regressiontechniques to a quantitative assessment for the motor signs of Parkinson’s disease. IEEETransactions on Neural Systems and Rehabilitation Engineering, 17(6):568–575, December2009.[20] W. Brogan. Modern Control Theory. Prentice Hall, Upper Saddle River, NJ, 3rd edition,1991.[21] E. Broussolle, P. Krack, S. Thobois, J. Xie-Brustolin, P. Pollak, and C. G. Goetz. Contribu-tion of jules froment to the study of parkinsonian rigidity. Movement Disorders, 22(7):909 –914, 2007.[22] R. A. Brown, D. A. Lawson, G. C. Leslie, A. MacArthur, W. J. MacLennan, M. E. McMurdo,W. J. Mutch, and N. J. Part. Does the wartenberg pendulum test differentiate quantitatively87Bibliographybetween spasticity and rigidity? a study in elderly stroke and parkinsonian patients. Journalof Neurology, Neurosurgery, and Psychiatry, 51(9):1178–1186, 1988.[23] R. G. Brown and C. D. Marsden. Internal versus external cues and the control of attentionin Parkinson’s disease. Brain, 111(2):323–345, 1988.[24] D. J. Burn, E. N. Rowan, L. M. Allan, S. Molloy, J. T. O’Brien, and I. G. McKeith. Motorsubtype and cognitive decline in parkinson’s disease, parkinson’s disease with dementia, anddementia with lewy bodies. Journal of Neurology, Neurosurgery, and Psychiatry, 77(5):585–589, 2006.[25] H. Chen-Harris, W. M. Joiner, V. Ethier, , D. S. Zee, and R. Shadmehr. Adaptive control ofsaccades via internal feedback. J. Neurosci., 28(11):2804 – 2813, 2008.[26] Ju-Yeop Choi, B.H. Cho, H.F. VanLandingham, Hyung-Soo Mok, and Joong-Ho Song. Systemidentification of power converters based on a black-box approach. Circuits and Systems I:Fundamental Theory and Applications, IEEE Transactions on, 45(11):1148–1158, 1998.[27] M. Cloutier, R. Middleton, and P. Wellstead. Feedback motif for the pathogenesis of parkin-son’s disease. IET Systems Biology, 6(3):86–93, 2012.[28] R Cools, RA Barker, BJ Sahakian, and TW Robbins. L-Dopa medication remediates cognitiveinflexibility, but increases impulsivity in patients with parkinson’s disease. Neuropsychologia,41:1431–1441, 2003.[29] P. Damier, E. C. Hirsch, Y. Agid, and A. M. Graybiel. The substantia nigra of the humanbrain: Ii. patterns of loss of dopamine-containing neurons in parkinson’s disease. Brain,122(8):1437–1448, 1999.[30] C.M. Davidson, A.M. De Paor, and M.M. Lowery. Insights from control theory into deepbrain stimulation for relief from parkinson’s disease. In ELEKTRO, 2012, pages 2–7, 2012.[31] P. R. Davidson, R. D. Jones, J. H. Andreae, and H. R. Sirisena. Simulating closed-andopen-loop voluntary movement: A nonlinear control-systems approach. IEEE Transactionson Biomedical Engineering, 49(11):1242–1252, November 2002.[32] B. Day, J. Dick, and C. Marsden. Patients with parkinson’s disease can employ a predictivemotor strategy. Journal of Neurology, Neurosurgery, and Psychiatry, 47(12):1299–1306, 1984.[33] Peter Dayan and L. F. Abbott. Theoretical Neuroscience: Computational and mathematicalmodeling of neural systems. MIT Press, Cambridge, MA, 2005.[34] P. J. Delwaide. Parkinsonian rigidity. Functional Neurology, 16:147–156, 2001.[35] S. Denve, J. R. Duhamel, and A. Pouget. Optimal sensorimotor integration in recurrentcortical networks: A neural implementation of kalman filters. The Journal of Neuroscience,27(21):5744–5756, 2007.88Bibliography[36] G. Deuschl, J. Raethjen, R. Baron, M. Lindemann, H. Wilms, and P. Krack. The pathophys-iology of parkinsonian tremor: a review. Neurology, 247(5):33–48, 2000.[37] R. C. Dorf and R. H. Bishiop. Modern Control Systems. Addison-Wesley Longman PublishingCo., Inc., Boston, MA, 1994.[38] Rosa L. Drijgers, Kathy Dujardin, Jennifer S.A.M. Reijnders, Luc Defebvre, and Albert F.G.Leentjens. Validation of diagnostic criteria for apathy in Parkinson’s disease. Parkinsonismand Related Disorders, 16(10):656 – 660, 2010.[39] H. Duffau, N. Tzourio, D. Caparros-Lefebvre, F. Parker, and B. Mazoyer. Tremor and vol-untary repetitive movement in parkinson’s disease: comparison before and after l-dopa withpositron emission tomography. Experimental Brain Research, 107(3):453–462, 1996.[40] Kathy Dujardin, Pascal Sockeel, David Devos, Marie Delliaux, Pierre Krystkowiak, AlainDeste, and Luc Defebvre. Characteristics of apathy in Parkinson’s disease. Movement Dis-orders, 22(6):778 – 784, 2007.[41] E.V. Evarts, H. Teravainen, and D.B. Calne. Reaction time in parkinson’s disease. Brain,104:167–186, 1981.[42] S. Fahn, R. L. Elton, and UPDRS Development Committee. In S. Fahn, C. D. Marsden,D. B. Calne, and M. Goldstein, editors, Unified Parkinsons Disease Rating Scale, RecentDevelopments in Parkinson’s Disease, pages 153–163. Florham Park, NJ: Macmillan, 1987.[43] Wouter Favoreel, Bart De Moor, and Peter Van Overschee. Subspace state space systemidentification for industrial processes. Journal of Process Control, 10(2–3):149–155, 2000.[44] J. M. Fearnley and A. J. Lees. Ageing and parkinson’s disease: substantia nigra regionalselectivity. Brain, 114:2283–2301, 1991.[45] S. Fekri, M. Athans, and A. Pascoal. Issues, progress and new results in robust adaptivecontrol. International Journal of Adaptive Control and Signal Processing, 20(10):519–579,August 2006.[46] X.-j. Feng, B. Greenwald, H. Rabitz, E. Shea-Brown, and R. Kosut. Toward closed-loopoptimization of deep brain stimulation for parkinson’s disease: concepts and lessons from acomputational model. Journal of Neural Engineering, 4(2):L14, 2007.[47] E.B.L. Filho, N.M.M. Rodrigues, E.A.B. da Silva, M.B. de Carvalho, S.M.M. de Faria, andV.M.M. da Silva. On ecg signal compression with 1-d multiscale recurrent patterns allied topreprocessing techniques. Biomedical Engineering, IEEE Transactions on, 56(3):896 –900,March 2009.[48] K. Flowers. Some frequency response characteristics of parkinsonism on pursuit tracking.Brain, 101(1):19–34, 1978.89Bibliography[49] J. Ghika, A.W. Wiegner, J.J. Fang, L. Davies, R.R. Young, and J.H. Growdon. Portablesystem for quantifying motor abnormalities in parkinson’s disease. Biomedical Engineering,IEEE Transactions on, 40(3):276–283, 1993.[50] L. Guo and D. Huang. Least-squares identification for armax models without the positivereal condition. Automatic Control, IEEE Transactions on, 34(10):1094–1098, 1989.[51] P.S. Hammon and V.R. de Sa. Preprocessing and meta-classification for brain-computerinterfaces. Biomedical Engineering, IEEE Transactions on, 54(3):518–525, March 2007.[52] Rick C. Helmich, Mark Hallett, Gunther Deuschl, Ivan Toni, and Bastiaan R. Bloem. Cerebralcauses and consequences of parkinsonian resting tremor: a tale of two circuits? Brain,135(11):3206–3226, 2012.[53] J. Herault, N. Guyader, and A. Guerin-Dugue. Scene variability and perception constancy inthe visual system: a model of pre-processing before data analysis and learning. In MachineLearning for Signal Processing, 2009. MLSP 2009. IEEE International Workshop on, pages1 –12, Sept. 2009.[54] M. M. Hoehn and M. D. Yahr. Parkinsonism : onset, progression, and mortality. Neurology,17(5):427–442, 1967.[55] Michael W. Hofbaur and Brian C. Williams. Mode estimation of probabilistic hybrid systems.In C. J. Tomlin and M. Greenstreet, editors, Hybrid Systems: Computation and Control,volume LNCS 2289, pages 253–266. Springer-Verlag, London, UK, 2002.[56] I. Hwang, H. Balakrishnan, and C. Tomlin. State estimation for hybrid systems: Applicationsto aircraft tracking. IEE Proceedings on Control Theory and Applications, 153(5):556–566,September 2006.[57] H. A. Ingram, P. van Donkelaar, J. Cole, J. L. Vercher, G. M. Gauthier, and R. C. Miall. Therole of proprioception and attention in a visuomotor adaptation task. Experimental BrainResearch, 132(1):114–126, 2000.[58] Daniel Graupe Ishita Basu, Daniela Tuninetti and Konstantin V. Slavin. Adaptive controlof deep brain stimulator for essential tremor: Entropy-based tremor prediction using surface-emg. In IEEE Engineering in Medicine and Biology Conference, pages 7711–7714, Boston,MA, 2011.[59] Makoto Ito and Kenji Doya. Multiple representations and algorithms for reinforcement learn-ing in the cortico-basal ganglia circuit. Current Opinion in Neurobiology, 21(3):368–373, 2011.[60] J. Jankovic. Parkinson’s disease: clinical features and diagnosis. Journal of Neurology,Neurosurgery and Psychiatry, 79(4):368 – 376, 2008.90Bibliography[61] J. Jankovic, M. McDermott, J. Carter, S. Gauthier, C. Goetz, L. Golbe, S. Huber, W. Koller,C. Olanow, I. Shoulson, M. Stern, C. Tanner, W. Weiner, and Parkinson Study Group. Vari-able expression of Parkinson’s disease: A base-line analysis of the dat atop cohort. Neurology,40(10):1529 – 1534, 1990.[62] Peifeng Ji, Ee-Leng Tan, Woon-Seng Gan, and Jun Yang. A comparative analysis of pre-processing methods for the parametric loudspeaker based on the khokhlov-zabolotskaya-kuznetsov equation for speech reproduction. Audio, Speech, and Language Processing, IEEETransactions on, 19(4):937–946, May 2011.[63] M. T. V. Johnson, A. N. Kipnis, J. D. Coltz, A. Gupta, P. Silverstein, F. Zwiebel, and T. J.Ebner. Effects of levodopa and viscosity on the velocity and accuracy of visually guidedtracking in Parkinson’s disease. Brain, 119:801–813, 1996.[64] I. T. Jolliffe. Principal Component Analysis. Springer Verlag, New York, 2002.[65] R. Jones, M. Donaldson, and B. Sharman. A technique for removal of the visuopercep-tual its application to parkinson’s disease. IEEE Transactions on Biomedical Engineering,43(10):1001–1010, 1996.[66] C.-F. Juang, C.-T. Chiou, and Chun-Lung Lai. Hierarchical singleton-type recurrent neu-ral fuzzy networks for noisy speech recognition. Neural Networks, IEEE Transactions on,18(3):833–843, May 2007.[67] M. Jueptner, I. H. Jenkins, D. J. Brooks, R. S. Frackowiak, and R. E. Passingham. The sensoryguidance of movement: a comparison of the cerebellum and basal ganglia. Experimental BrainResearch, 112(3):462–474, 1996.[68] M. Jueptner and C. Weiller. A review of differences between basal ganglia and cerebellarcontrol of movements as revealed by functional imaging studies. Brain, 121(8):1437–1449,1998.[69] B.S. Kim and S.K. Yoo. Motion artifact reduction in photoplethysmography using indepen-dent component analysis. Biomedical Engineering, IEEE Transactions on, 53(3):566–568,March 2006.[70] Alexxai V. Kravitz, Benjamin S. Freeze, Philip R. L. Parker, Kenneth Kay, Myo T. Thwin,Karl Deisseroth, and Anatol C. Kreitzer. Regulation of parkinsonian motor behaviours byoptogenetic control of basal ganglia circuitry. Nature, 466:622–626, 2010.[71] Arthur D Kuo. An optimal state estimation model of sensory integration in human posturalbalance. Journal of Neural Engineering, 2(3):235–249, 2005.[72] Yasar Kutukcu, William J. Marks, Douglas S. Goodin, and Michael J. Aminoff. Simple andchoice reaction time in Parkinson’s disease. Brain Research, 815(2):367–372, 1999.91Bibliography[73] A. E. Lang and A. M. Lozano. Parkinson’s disease. first of two parts. The New EnglandJournal of Medicine, 339(15):1044–1053, 1998.[74] Sarah Lemieux, Mehrdad Ghassemi, Mandar Jog, Roderick Edwards, and Christian Duval.The influence of levodopa-induced dyskinesias on manual tracking in patients with parkinson’sdisease. Experimental Brain Research, 176(3):465–475, January 2007.[75] Michael D. Lewek, Roxanne Poole, Julia Johnson, Omar Halawa, and Xuemei Huang. Armswing magnitude and asymmetry during gait in the early stages of parkinson’s disease. Gaitand Posture, 31(2):256–260, 2010.[76] S. J. G. Lewis, T. Foltynie, A. D. Blackwell, T. W. Robbins, A. M. Owen, and R. A. Barker.Heterogeneity of parkinson’s disease in the early clinical stages using a data driven approach.Journal of Neurology, Neurosurgery, and Psychiatry, 76(3):343–348, 2005.[77] C. C. Lin, M. S. Ju, and C. W. Lin. The pendulum test for evaluating spasticity of the elbowjoint. Archives of Physical Medicine and Rehabilitation, 84(1):69 – 74, 2003.[78] Dan Liu and Emanuel Todorov. Evidence for the flexible sensorimotor strategies predictedby optimal feedback control. Journal of Neuroscience, 27(35):9354–9368, August 2007.[79] Xuguang Liu, Robert Osterbauer, Tipu Z. Aziz, R. Christopher Miall, and John F. Stein.Increased response to visual feedback of drug-induced dyskinetic movements in advancedParkinson’s disease. Neuroscience Letters, 304:25–28, 2001.[80] L. Ljung. System Identification Toolbox for use with MATLAB. The MathWorks Inc., SouthNatick, Mass., 1986.[81] L. Ljung. System Identification - Theory For the User. Prentice Hall, Upper Saddle River,NJ, 1999.[82] L. Ljung and Zhen-Dong Yuan. Asymptotic properties of black-box identification of transferfunctions. Automatic Control, IEEE Transactions on, 30(6):514–530, 1985.[83] J.A.R. Macias and A.G. Exposito. Self-tuning of kalman filters for harmonic computation.Power Delivery, IEEE Transactions on, 21(1):501 – 503, Jan. 2006.[84] P. S. Maybeck. Stochastic Models, Estimation and Control, volume 2. Academic Press, NewYork, 1982.[85] L. F. Medeiros, De Souza I. C., L. P. Vidor, A. De Souza, A. Deitos, M. S. Voltz, F. Fregni,W. Caumo, and I. L. S. Torres. Neurobiological effects of transcranial direct current stimu-lation: A review. Front Psychiatry, 3(110), 2012.[86] R. Miall and Dominic King. State estimation in the cerebellum. The Cerebellum, 7:572–576,2008. 10.1007/s12311-008-0072-6.92Bibliography[87] R. C. Miall and E. W. Jenkinson. Functional imaging of changes in cerebellar activity relatedto learning during a novel eye-hand tracking task. Experimental Brain Research, 166(2):170–183, 2005.[88] R. C. Miall, D. J. Weir, D. M. Wolpert, and J. F. Stein. Is the cerebellum a Smith predictor?Journal of Motor Behavior, 25(3):203–216, 1993.[89] R. Chris Miall, Lars O. D Christensen, Owen Cain, and James Stanley. Disruption of stateestimation in the human lateral cerebellum. PLoS Biol, 5(11):2733–2744, 11 2007.[90] P. K. Morrish, G. V. Sawle, and D. J. Brooks. An [18f]dopa-pet and clinical study of the rateof progression in parkinson’s disease. Brain, 119(2):585–591, 1996.[91] K. Narendra and C. Xiang. Adaptive control of discrete-time systems using multiple models.IEEE Transactions on Automatic Control, 45(9):1669–1686, September 2000.[92] M. A. Nitsche, S. Doemkes, T. Karakse, A. Antal, D. Liebetanz, N. Lang, F. Tergau, andW. Paulus. Shaping the effects of transcranial direct current stimulation of the human motorcortex. Journal of Neurophysiology, 97(4):3109–3117, 2007.[93] Robert L. Nussbaum and Christopher E. Ellis. Alzheimer’s disease and parkinson’s disease.New England Journal of Medicine, 348(14):1356–1364, 2003.[94] Miyako Oguru, Hisao Tachibana, Kazuo Toda, Bungo Okuda, and Nobuyuki Oka. Apathyand depression in parkinson disease. Journal of Geriatric Psychiatry and Neurology, 23(1):35– 41, 2010.[95] M.M.K. Oishi, P. Talebifard, and M.J. McKeown. Assessing manual pursuit trackingin parkinson’s disease via linear dynamical systems. Annals of Biomedical Engineering,39(8):2263–2273, 2011.[96] A. V. Oppenheim and R. W. Schafer. Discrete-time signal processing. Prentice Hall, 1989.[97] Charles D. Ormsby, John F. Raquet, and Peter S. Maybeck. A new generalized residualmultiple model adaptive estimator of parameters and states. Mathematical and ComputerModeling, 43(9-10):1092–1113, 2006.[98] Frank Osafo-Charles, Gyan C. Agarwal, William D. O’Neill, and Gerald L. Gottlieb. Appli-cation of time-series modeling to human operator dynamics. IEEE Transactions on Systems,Man and Cybernetics, 10(12):849–860, December 1980.[99] D. Page and M. Jahanshahi. Deep brain stimulation of the subthalamic nucleus improves setshifting but does not affect dual task performance in Parkinson’s disease. IEEE Transactionson Neural Systems and Rehabilitation Engineering, 15(2):198–206, June 2007.[100] S. Pal, S. M. Rosengren, and J. G. Colebatch. Stochastic galvanic vestibular stimulationproduces a small reduction in sway in parkinson’s disease. Journal of Vestibular Research:equilibrium and orientation, 19(3):137 – 142, 2009.93Bibliography[101] S.J. Palmer, J. Li, Z.J. Wang, and M.J. McKeown. Joint amplitude and connectivity com-pensatory mechanisms in parkinson’s disease. Neuroscience, 166(4):1110–1118, 2010.[102] W. Pan, R. Soma, S. Kwak, and Y. Yamamoto. Improvement of motor functions bynoisy vestibular stimulation in central neurodegenerative disorders. Journal of Neurology,255(11):1657 – 1661, 2008.[103] A. M. Paolo, B. N. Axelrod, A. I. Tro¨ster, K. T. Blackwell, and W. C. Koller. Utility of aWisconsin Card Sorting Test short form in persons with Alzheimer’s and Parkinson’s disease.Journal of Clinical and Experimental Neuropsychology, 18(6):892–897, 1996.[104] S. K. Patrick, A. A. Denington, M. J. Gauthier, D. M. Gillard, and A. Prochazka. Quantifi-cation of the updrs rigidity scale. IEEE Transactions on Neural Systems and RehabilitationEngineering, 9(1):31 – 41, 2001.[105] M. G. Paulin, L. F. Hoffman, and C. Assad. A model of cerebellar computations for dynamicalstate estimation. Autonomous Robots, 11(3):279–284, 2001.[106] J. C. Pchadre, L. Larochelle, and L. J. Poirier. Parkinsonian akinesia, rigidity and tremor inthe monkey. histopathological and neuropharmacological study. Journal of the NeurologicalSciences, 28(2):147–157, 1976.[107] G. Pedoto, S. Santaniello, G. Fiengo, L. Glielmo, M. Hallett, Ping Zhuang, and S.V. Sarma.Towards automated navigation of deep brain stimulating electrodes: Analyzing neuronalactivity near the target. In IEEE International Conference on Control Applications, pages782–787, 2012.[108] A. Plotnik and S. Rock. Improved estimation of target velocity using multiple model esti-mation and a dynamic bayesian network for a robotic tracker of ocean animals. In RoboticsResearch, pages 402–415. Springer-Verlag, 2007.[109] Janey Prodoehl, Daniel M. Cross, John C. Rothwell, Leo Verhagen Metman, Roy A. E.Bakay, and David E. Vaillancourt. Effects of STN DBS on memory guided force control inParkinson’s disease. IEEE Transactions on Neural Systems and Rehabilitation Engineering,15(2):155–165, June 2007.[110] M. C. Rodriguez-Oroz, M. Jahanshahi, P. Krack, I. Litvan, R. Macias, E. Bezard, and J. A.Obeso. Initial clinical manifestations of parkinson’s disease: features and pathophysiologicalmechanisms. Lancet Neurology, 8(12):1128–1139, 2009.[111] B. Rosin, M. Slovik, R. Mitelman, M. Rivlin-Etzion, S.N. Haber, Z Israel, E Vaadia, andH. Bergman. Closed-loop deep brain stimulation is superior in ameliorating parkinsonism.Neuron, 72(2):370–384, 2011.[112] A. Salarian, H. Russmann, C. Wider, P.R. Burkhard, F. J G Vingerhoets, and K. Aminian.Quantification of tremor and bradykinesia in parkinson’s disease using a novel ambulatorymonitoring system. Biomedical Engineering, IEEE Transactions on, 54(2):313–322, 2007.94Bibliography[113] S. Santaniello, G. Fiengo, L. Glielmo, and W. M. Grill. Basal ganglia modeling in healthyand Parkinson’s disease state – isolated neurons activity. In Proceedings of the AmericanControl Conference, pages 4089 – 4094, New York, July 2007.[114] S. Santaniello, G. Fiengo, L. Glielmo, and W.M. Grill. Closed-loop control of deep brainstimulation: a simulation study. IEEE Transactions on Neural Systems and RehabilitationEngineering, 19(1):15–24, 2011.[115] M. C. Schiess, H. Zheng, V. M. Soukup, J. G. Bonnen, and H. J. W. Nauta. Parkinson’s diseasesubtypes: clinical classification and ventricular cerebrospinal fluid analysis. Parkinsonism andRelated Disorders, 6(2):69 – 76, 2000.[116] R. Schulz, C. Gerloff, and F. C. Hummel. Non-invasive brain stimulation in neurologicaldiseases. Neuropharmacology, 64:579 – 587, 2012.[117] B. Sepehri, A. Esteki, E. Ebrahimi-Takamjani, G. A. Shahidi, F. Khamseh, and M. Moin-odin. Quantification of rigidity in parkinson’s disease. Annals of Biomedical Engineering,35(12):2196–2203, 2007.[118] R. Shadmehr and J. W. Krakauer. A computational neuroanatomy for motor control. Ex-perimental Brain Research, 185(3):359–381, 2008.[119] R Shadmehr, MA Smith, and JW Krakauer. Error correction, sensory prediction, and adap-tation in motor control. Annual Reviews of Neuroscience, 33:89–108, 2010.[120] Reza Shadmehr and Sandro Mussa-Ivaldi. Biological learning and control: How the brainbuilds representations, predicts events, and makes decisions. MIT Press, Cambridge, MA,2012.[121] M. M. Shafi, M. B. Westover, M. D. Fox, and A. Pascual-Leone. Exploration and modu-lation of brain network interactions with noninvasive brain stimulation in combination withneuroimaging. European Journal of Neuroscience, 35(6):805 – 825, 2012.[122] M. B. Shapiro, D. E. Vaillancourt, M. M. Sturman, L. V. Metman, R. A. Bakay, and D. M.Corcos. Effects of stn dbs on rigidity in parkinson’s disease. IEEE Trans. Neural Syst. Rehabil.Eng., 15(2):173 – 181, 2007.[123] K. M. Shaw, A. J. Lees, and G. M. Stern. The impact of treatment with levodopa onparkinson’s disease. Quarterly Journal of Medicine, 49(195):283 – 293, 1980.[124] Kaiquan Shen, Ke Yu, A. Bandla, Yu Sun, N. Thakor, and Xiaoping Li. Multiple time-lag canonical correlation analysis for removing muscular artifacts in eeg. In Engineering inMedicine and Biology Society (EMBC), 2013 35th Annual International Conference of theIEEE, pages 6792–6795, July 2013.[125] Stanley M. Shinners. Modeling of human operator performance utilizing time series analysis.Systems, Man and Cybernetics, IEEE Transactions on, SMC-4(5):446–458, 1974.95Bibliography[126] Joshua M. Shulman, Philip L. De Jager, and Mel B. Feany. Parkinson’s disease: Geneticsand pathogenesis. Annual Review of Pathology: Mechanisms of Disease, 6(1):193–222, 2011.[127] P. F. Smith, L. H. Geddes, J. H. Baek, C. L. Darlington, and Y. Zheng. Modulation ofmemory by vestibular lesions and galvanic vestibular stimulation. Frontiers in Neurology,1(141), 2010.[128] B. Sohlberg. Hybrid grey box modelling of a pickling process. Control Engineering Practice,13(9):1093 – 1102, 2005. Modelling and Control of Biomedical Systems.[129] S. Soveyni and R. B. Gillespie. Cancellation of biodynamic feedthrough in vehicle controltasks. IEEE Transactions on Control Systems Technology, 15(6):1018–1029, November 2007.[130] S. Spieker, A. Boose, S. Breit, and J. Dichgans. Long-term measurement of tremor. MovementDisorders, 13(3):81–84, 1998.[131] J. K. R. Stevenson, M. Oishi, S. Farajian, E. Cretu, E. Ty, and M. J. McKeown. Responseto sensory uncertainty in parkinson’s disease: A marker of cerebellar dysfunction? EuropeanJournal for Neuroscience, 83(2):298–305, February 2011.[132] C. M. Tanner and Goldman S. M. Epidemiology of parkinson’s disease. Neurol. Clin.,14(2):317–335, 1996.[133] Denmark The Lundbeck Institute. The dopamine pathways in ocd, available online at http ://www.cnsforum.com/imagebank/item/Neuro path DA OCD/default.aspx. 2001.[134] Chung Tin and Chi-Sang Poon. Internal models in sensorimotor integration: perspectivesfrom adaptive control theory. Journal of Neural Engineering, 2(3):147–163, 2005.[135] Emanuel Todorov. Optimality principles in sensorimotor control. NATURE NEURO-SCIENCE, 7(9):907–915, SEPTEMBER 2004.[136] A. Tsanas, M.A. Little, P.E. McSharry, and L.O. Ramig. Accurate telemonitoring of parkin-son’s disease progression by noninvasive speech tests. Biomedical Engineering, IEEE Trans-actions on, 57(4):884–893, 2010.[137] K. S. Utz, V. Dimova, K. Oppenlnder, and G. Kerkhoff. Electrified minds: transcranialdirect current stimulation (tdcs) and galvanic vestibular stimulation (gvs) as methods ofnon-invasive brain stimulation in neuropsychology–a review of current data and future impli-cations. Neuropsychologia, 48(10):2789–2810, 2010.[138] MariaS Valle, Antonino Casabona, Rosaria Sgarlata, Rosaria Garozzo, Maria Vinci, andMatteo Cioni. The pendulum test as a tool to evaluate passive knee stiffness and viscosity ofpatients with rheumatoid arthritis. BMC Musculoskeletal Disorders, 7(1):1–12, 2006.[139] Mats Viberg. Subspace-based methods for the identification of linear. Automatica,31(12):1835–1851, 1995.96[140] R. Vidal, A. Chiuso, S. Soatto, and S. Sastry. Observability of linear hybrid systems. InR. Alur and G. Pappas, editors, Hybrid Systems: Computation and Control, volume LNCS2293, pages 526–539. Springer-Verlag, 2003.[141] R. W. Watson, R. D. Jones, and N. B. Sharman. Two-dimensional tracking tasks for quan-tification of sensory-motor dysfunction and their application to Parkinson’s disease. Medicaland Biological Engineering and Computing, 35:141–145, 1997.[142] A. L. Whone, R. L. Watts, A. J. Stoessl, M. Davis, S. Reske, C. Nahmias, A. E. Lang,O. Rascol, M. J. Ribeiro, P. Remy, W. H. Poewe, R. A. Hauser, D. J. Brooks, and REAL-PET Study Group. Slower progression of parkinson’s disease with ropinirole versus levodopa:The real-pet study. Brain, 54(1):93–101, 2003.[143] T. Wichmann and M.R. DeLong. In P. Riederer, H. Reichmann, M.B.H. Youdim, and M. Ger-lach, editors, Parkinson’s Disease and Related Disorders, volume 70 of Journal of NeuralTransmission. Supplementa, pages 21–25. Springer Vienna, 2006.[144] Thomas Wichmann, Mahlon R. DeLong, Jorge Guridi, and Jose A. Obeso. Milestones inresearch on the pathophysiology of parkinson’s disease. Movement Disorders, 26:1032–1041,2011.[145] D. Widjaja, S. Vandeput, J. Taelman, M.A.K.A. Braeken, R.A. Otte, B.R.H. Van den Bergh,and S. Van Huffel. Accurate r peak detection and advanced preprocessing of normal ecg forheart rate variability analysis. In Computing in Cardiology, 2010, pages 533–536, Sept. 2010.[146] Daniel M. Wolpert, R. Chris Miall, and Mitsuo Kawato. Internal models in the cerebellum.Trends in Cognitive Sciences, 2(9):338–347, 1998.[147] Ruiping Xia, M. Radovic, A.J. Threlkeld, and Zhi-Hong Mao. System identification andmodeling approach to characterizing rigidity in parkinson’s disease: Neural and non-neuralcontributions. In Bioinformatics and Biomedical Engineering (iCBBE), 2010 4th Interna-tional Conference on, pages 1–4, 2010.[148] Hong-Tzer Yang, Chao-Ming Huang, and Ching-Lien Huang. Identification of armax modelfor short term load forecasting: an evolutionary programming approach. Power Systems,IEEE Transactions on, 11(1):403–408, 1996.[149] Einat Yehene, Nachshon Meiran, and Nachum Soroker. Basal ganglia play a unique rolein task switching within the frontal-subcortical circuits: Evidence from patients with focallesions. J. Cognitive Neuroscience, 20(6):1079–1093, 2008.[150] A. Zenzola, G. Masi, M. Mari, G. Defazio, P. Livrea, and P. Lamberti. Fatigue in parkinson’sdisease. Neurological Sciences, 225(3):225–226, 2003.97Appendix AEstimator dynamics in switching taskParameters for estimator dynamics for all subjects and tasks are in Table A.1.98Appendix A. Estimator dynamics in switching taskTable A.1: Damping ratio and natural frequency in estimator dynamics for all normal and PDsubjects in ‘Better’ mode.Subject Better Normal Worsenumber (ζ, ωn) (ζ, ωn) (ζ, ωn)Normal #1 (0.0011, 38.48) (0.0030, 38.55) (0.0439, 39.69)Normal #2 (0.0029, 38.56) (0.0265, 39.12) (0.0120, 38.84)Normal #3 (0.3914, 49.29) (1.937, 85.17) (0.0420, 39.57)Normal #4 (0.0064, 38.65) (0.0064, 38.69) (0.0020, 38.57)Normal #5 (0.0309, 39.22) (0.0132, 38.85) (0.0054, 38.65)Normal #6 (0.0187, 38.94) (0.1627, 42.62) (0.0433, 39.51)Normal #7 (0.0171, 38.90) (0.0015, 38.56) (0.0194, 39.02)Normal #8 (0.0302, 39.20) (0.0406, 39.52) (0.0741, 40.36)Normal #9 (1.330, 81.49) (0.0481, 39.64) (1.758, 74.20)Normal #10 (0.0139, 38.82) (0.0032, 38.55) (0.0321, 39.29)PD #1 (0.0017, 38.51) (1.787, 71.69) (2.070, 96.62)PD #2 (0.0038, 38.59) (0.0152, 38.84) (0.2691, 45.02)PD #3 (0.0036, 38.57) (0.0487, 39.64) (0.4433, 51.64)PD #4 (0.0030, 38.55) (0.0264, 39.14) (1.804, 71.43)PD #5 (0.0232, 39.03) (0.0202, 38.84) (0.0755, 40.21)PD #6 (0.0031, 38.69) (0.0265, 39.11) (0.0376, 39.39)PD #7 (0.0039, 38.62) (0.0101, 38.85) (0.0187, 39.08)PD #8 (0.0092, 39.70) (0.0145, 38.83) (0.0028, 38.63)PD #9 (0.0107, 38.73) (0.0160, 38.87) (0.0196, 38.91)PD #10 (0.0040, 38.61) (0.0044, 38.62) (0.0508, 39.81)PD #11 (0.0153, 38.83) (0.0070, 38.65) (0.0057, 38.63)PD #12 (0.0134, 38.79) (0.0015, 38.52) (0.0199, 38.97)PD #13 (0.0072, 38.64) (0.0032, 38.59) (0.0305, 39.23)PD #14 (0.0011, 38.53) (0.0087, 38.69) (0.0082, 38.70)99
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Analysis of motor performance in Parkinson's disease...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Analysis of motor performance in Parkinson's disease through LTI dynamical systems Ashoori, Ahmad 2014
pdf
Page Metadata
Item Metadata
Title | Analysis of motor performance in Parkinson's disease through LTI dynamical systems |
Creator |
Ashoori, Ahmad |
Publisher | University of British Columbia |
Date Issued | 2014 |
Description | Motor performance in Parkinson's disease is marked by a variety of impairments, and is indicative of faulty feedback mechanisms in the brain that arise due to a lack of dopamine. Quantifying motor performance as well as investigating the correlation of motor behaviour and clinical measures has been sought as in important avenue to study disease pathophysiology. We propose the use of linear dynamical systems to finely characterize motor performance, and for possible interpretations of neurological phenomena, such as compensatory mechanisms that can be cast in the framework of control theory. While coarse measurements of motor performance exist (e.g., maximum velocity, mean velocity, root mean square error), they often cannot distinguish between qualitatively different types of movement. (Consider for example, the wide variety of step response behaviors possible in second order systems which may have the same settling time, but a wide range of damping ratios.) In contrast, LTI models provide a dynamical mapping of the sensorimotor transformation required in tracking tasks. Model parameters of linear dynamical systems (such as damping ratio, and natural frequency) can more readily describe qualitatively similar motor phenomenon than existing coarse measurements. From time-series data of various manual tracking experiments, we first perform a set of preprocessing techniques dependent on the particular experimental setup. The pre-processing techniques have been carefully evaluated to appropriately address possible noise, nonlinearities, and other data irregularities. We then apply existing tools for system identification (e.g., ARX, ARMAX, subspace methods) to identify numerically robust second-order LTI system models. The choice of system identification method (as well as choice of pre-processing techniques) is critical in determining the variability in the resulting model parameters. We identify parameters of the LTI models which are amenable to comparison across subjects (and across groups), and evaluate their statistical significance. The main contributions of this thesis are: - Selection of pre-processing and system identification techniques for three experimental setups for manual tracking experiments - Control theoretic framework for tremor as a compensatory mechanism that is advantageous in some tracking tasks - Sensor-based assessment of rigidity via decay rate - Assessment of cognitive inflexibility through mode detection and delay |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2014-12-15 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0135634 |
URI | http://hdl.handle.net/2429/51518 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2015-02 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2015_february_ashoori_ahmad.pdf [ 2.05MB ]
- Metadata
- JSON: 24-1.0135634.json
- JSON-LD: 24-1.0135634-ld.json
- RDF/XML (Pretty): 24-1.0135634-rdf.xml
- RDF/JSON: 24-1.0135634-rdf.json
- Turtle: 24-1.0135634-turtle.txt
- N-Triples: 24-1.0135634-rdf-ntriples.txt
- Original Record: 24-1.0135634-source.json
- Full Text
- 24-1.0135634-fulltext.txt
- Citation
- 24-1.0135634.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0135634/manifest