PROBABILISTIC BOOLEAN NETWORK MODELING FOR fMRI STUDY IN PARKINSON’S DISEASE by Zheng Ma B.A.Sc., The University of British Columbia, 2006 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Applied Science in The Faculty of Graduate Studies (Electrical and Computer Engineering) The University of British Columbia September, 2008 © Zheng Ma 2008 Abstract Recent research has suggested disrupted interactions between brain regions may contribute to some of the symptoms of motor disorders such as Parkinson’s Disease (PD). It is therefore important to develop models for inferring brain functional connectivity from data obtained through non-invasive imaging technologies, such as functional magnetic resonance imaging (fMRI). The complexity of brain activities as well as the dynamic nature of motor disorders require such models to be able to perform complex, large-scale, and dynamic system compu tation. Traditional models proposed in the literature such as structural equation modeling (SEM), multivariate autoregressive models (MAR), dynamic causal modeling (DCM), and dynamic Bayesian networks (DBNs) have all been suggested as suitable for fMRI data analy sis. However, they suffer from their own disadvantages such as high computational cost (e.g. DBNs), inability to deal with non-linear case (e.g. MAR), large sample size requirement (e.g. SEM), et., al. In this research, we propose applying Probabilistic Boolean Network (PBN) for modeling brain connectivity due to its solid stochastic properties, computational simplicity, robustness to uncertainty, and capability to deal with small-size data, typical for fIVIRI data sets. Applying the proposed PBN framework to real fMRI data recorded from PD subjects enables us to identify statistically significant abnormality in PD connectivity by comparing it with normal subjects. The PBN results also suggest a mechanism of evaluating the effectiveness of L-dopa, the principal treatment for PD. In addition to PBNs’ promis ing application in inferring brain connectivity, PBN modeling for brain ROTs also enables researchers to study dynamic activities of the system under stochastic conditions, gaining essential information regarding asymptotic behaviors of ROTs for potential therapeutic in tervention in PD. The results indicate significant difference in feature states between PD 11 patients and normal subjects. Hypothesizing the observed feature states for normal subject as the desired functional states, we further explore possible methods to manipulate the dy namic network behavior of PD patients in the favor of the desired states from the view of random perturbation as well as intervention. Results identified a target ROT with the best intervention performance, and that ROl is a potential candidate for therapeutic exercise. 111 Table of Contents Abstract . . . ii Table of Contents iv List of Tables . . vi List of Figures . . vii Acknowledgements ix 1 Introduction . 1 1.1 Background and Motivation 1 1.2 Background on functional Magnetic Resonance Imaging 3 1.3 Literature Survey . . . . 4 1.3.1 Structural Equation Modeling . . . 5 1.3.2 Multivariate Autoregressive Model . 7 1.3.3 Dynamic Causal Modeling 8 1.3.4 Dynamic Bayesian Networks 11 1.3.5 Group-Analysis Methods 11 1.4 Research Objectives and Methodology 13 1.5 Bibliography 17 2 PBN Connectivity Analysis 18 2.1 Introduction 18 iv 2.2 Methods 22 2.2.1 Preprocessing of fIVIRI data 22 2.2.2 Introduction of PBN 26 2.2.3 Learning PBNs 29 2.2.4 Group Analysis for Inter-Subject Variability 31 2.3 fMRI Case Study 32 2.3.1 fMRI Data Collection 32 2.3.2 Results and Discussions 34 2.4 Discussion and Conclusion 37 2.5 Bibliography 56 3 PBN Dynamic Analysis 59 3.1 Introduction 59 3.2 Methods 60 3.2.1 Theoretical Background of PBN 60 3.2.2 Dynamics of PBN 61 3.2.3 PBN Intervention 63 3.3 fIVIRI Case Study 65 3.3.1 fIVIRI Data Collection 65 3.3.2 Preprocessing of fIVIRI data 65 3.3.3 Results and Discussions 66 3.4 Conclusion 68 3.5 Bibliography 72 4 Conclusion 73 4.1 Conclusions 73 4.2 Future Work 74 V List of Tables 2.1 Truth Table for BN Example 40 2.2 Truth Table for PBN extended from BN example 43 2.3 Connections showing significantly difference due to the group factor, Left Hand 46 2.4 Performance of Medicine on PD, Left Hand 47 2.5 Connections showing significantly difference due to the frequency factor, Left Hand 52 2.6 Performance of Frequency on subjects, Left Hand 53 3.1 Intervention performance comparison based on Criterion II, Mmo (x(1), y) where = 20 71 vi List of Figures 1.1 fMRI Scan and Network Learning 15 1.2 Path Diagram of 4-ROT diagram 15 1.3 DBN example 16 2.1 Overall fMRI processing procedure 39 2.2 Boolean Functions for node x1, node x2, node x3, respectively 41 2.3 A basic building block for PBN 42 2.4 State Transition Diagram for PBN 43 2.5 An example of the learned influences from the fMRI data for a normal subject performing a High frequency squeezing using Left Hand 44 2.6 The complete experiment specification diagram 45 2.7 Examples of Connections in which L-dopa has normalizing effects 48 2.8 Projected Brain Connections where L-dopa has normalizing effects 49 2.9 Examples of Connections in which L-dopa deteriorates PD symptoms or dopa indicates “overshoot” effects 50 2.10 Projected Brain Connections where L-dopa either deteriorates PD symptoms or indicates “overshoot” effects 51 2.11 Multiple-Comparison Results of Influence Connections under Frequency factor 54 2.12 Projected Significant Brain Connections identified from Multiple-Comparison Method under Frequency factor 55 3.1 6-ROT case study 69 vii 3.2 3-ROTs state diagrams after modification to eliminate small transition prob abilities created by perturbation 70 3.3 Intervention performance comparison based on Criterion I, Hp (X(), y) for = 1...20, with y =1111’ 71 viii Acknowledgements I want to convey my greatest appreciation to my supervisor, Dr. Z.Jane Wang, whose patient guidance, persistent encouragement and deep insight in the research area have helped me immeasurably throughout the course of my thesis research. I am also grateful for Dr. Martin McKeown whose profound knowledge of Parkinson’s disease has provided me with precious comments and suggestions during my research. This thesis would never have been written without his assistance. I would like to thank my thesis examination committee members for their time and efforts in reading my thesis and providing valuable comments. My deepest gratitude goes to my family for their endless love, constant support, and immense encouragement over the years. This work was supported by NSERC agency funding. ix Chapter 1 Introduction 1.1 Background and Motivation Parkinson’s Disease, known as PD, is a degenerative disorder of the central nervous system that often impairs the sufferer’s motor skills and speech ability, as well as other functions [1]. PD’s primary symptoms include muscle rigidity, tremor, a slowing of physical move ment (bradykinesia) and, in extreme cases, a loss of physical movement (akinesia), all of which could be summarized as motor disorders. Motor disorders are caused by weak or disorganized functional connections in brain, whose primary task is to process informa tion received, and send message to body parts requiring actions through those connections. There are several factors contributing to the impairment of functional connections. In PD’s case, it is due to the insufficient formation and action of dopamine - a type of chemical compuound produced in the dopaminergic neurons of the brain, and acts as pathways for the basal ganglia to stimulate the motor cortex [2]. As a result, studying effective connec tions, represented as the neural influence that one brain region exerts over another, is of great importance for understanding the underlying pathological mechanism of PD as well as exploring future treatment. However, researchers and doctors were not able to learn the activity of effective connections in great depth until the introduction of non-invasive medical technologies such as functional magnetic resonance imaging (fMRI) and quantitative elec troencephalogram (EEG), both of which have enhanced understanding about human brain functioning in normal and disease states. As one of the most recently developed forms of neuroimaging, f1VIRI scans the brain ac 1 tivities of a subject performing certain tasks while minimizing his/her discomfort which may result in inaccurate measurements. A complet fMRI data learning framework is illustrated in flg.1.1. fMRI is characterized with two primary advantages. Firstly, it can noninvasively record brain signals without risks of radiation inherent in other scanning methods, such as CT scans. Secondly, fMRI can record on a high spatial resolution in a region of 3-6 millimeters, ideal for brain research [3]. However, it is often argued that fMRI may perform not as good as other imaging technologies like EEG in providing high temporal resolution. Such difference is not caused by the technique itself, but by the phenomena being measured since EEG measures electrical/neural activity while fMRI measures blood activity, which tends to have a longer response. Furthermore, the Blood-oxygen-level dependent or BOLD signals that fiVIRI records are believed to be closely related to the neural activity charac terizing effective connections. The oxygen consumption needed by active neurons leads to an increase in blood flow to regions of increased neural activity, occurring after a delay of approximately 1-5 seconds. Such increase rises to a peak over 4-5 seconds, before falling back to baseline (and typically undershooting slightly) [4]. This results iii local changes in the relative concentration of oxyhemoglobin and deoxyhemoglobin and changes in local cerebral blood volume in addition to this change in local cerebral blood flow. As an effort to intelligently extract and comprehensively understand massive informa tion embedded in fIVIRI data, researchers have been increasingly interested in mathematical networks that could facilitate systematic study of effective connections across different re gions of interest (ROIs). In this research, Probabilistic Boolean Networks (PBNs) [5] are employed to analyze fMRI data for the first time. PBNs aim to constructing logical interac tions among ROTs characterizing underlying connectivity pattern without knowing specific biological details. Therefore, PBNs allow researchers to study brain connectivity on the structure level, and such connectivity study currently is a hot issue for PD’s diagnosis as well as its treatment. Along with group analysis methods, we focus on detecting significant difference between PD patients and normal subjects in effective connections since those dif. 2 ference may shed light on further disease study. Another main goal of modeling PBN for brain ROIs is to identify potential therapy target through the asymptotic analysis, which could be facilitated by using the numerous theories and tools developed for markov chain. Meanwhile, we are also interested in suggesting way of changing the current trajectory of PD network by intentionally intervening identified target so that the system will be more likely to reach some desired states observed in normal subjects. To minimize the system disturbance brought by intervention, we plan to focus on transient intervention whose effects could be reversed by the network itself. The rest of this chapter is organized as follows: Section 1.2 presents a brief introduction of the fIVIRI signals. It is followed by a literature survey of traditional techniques for inferring brain connectivity, and current group analysis methods in Section 1.3. Finally, the research objective and methodology are described in Section 1.4. 1.2 Background on functional Magnetic Resonance Imaging Nerve cells do not have reserved energy necessary to maintain their functioning. Therefore, they obtain oxygen released by blood to power themselves up through a process called the hemodynamic response. Since the magnetic resonance (MR) signal of blood depends on the level of oxygenation, differential signals resulted from blood with different concerntration of oxygenated hemoglobin can be detected using an appropriate MR pulse sequence as blood- oxygen-level dependent (BOLD) contrast [4]. A high oxygen transfering rate from blood to active nerve cells results in the decrease of the concentration of oxygenated hemoglobin in blood, therefore decreasing intensities of BOLD signals since the blood magnetic susceptibil ity less closely matches the tissue magnetic susceptibility, conversly low oxygen transfering rate increases intensities of BOLD signals. By collecting data in an MRI scanner with mag netic susceptibility related parameters, researchers can estimate changes in BOLD contrast. 3 The direction of these changes (either positive or negative) depends on the relative changes in both cerebral blood flow (CBF) and oxygen consumption [6]. Increases in CBF that outstrip changes in oxygen consumption lead to stronger BOLD signal, while decreases in CBF that outstrip changes in oxygen consumption cause decreased BOLD signal intensity. Based on changes in BOLD contrast, researchers are able to identify which areas of brain are more activly involved in certain type of action through data of many repetitions of that action. Changes in BOLD signal have been long known to be well correlated with changes in blood flow. Furthermore, several studies during the past several decades have confirmed the fact that blood flow is tightly regulated spatially as well as temporally to sustain brain metabolism [4]. However, neuroscientists have still been seeking for a more direct relation ship between neural signals and BOLD signals. BOLD effects are measured using rapid volumetric acquisition of images with weighted contrast. Such images can be obtained with moderately good spatial and temporal reso lution. Images are usually taken every 14 seconds with the voxels in the resulting image typically representing cubes of tissue about 24 millimeters on each side in human brains [6]. Recent technical advancements, such as the use of high magnetic fields and advanced “multichannel” RF reception, have increased the spatial resolution up to the millimeter scale. Although responses to stimuli closely presented within one or two seconds can be distinguished from one another by the event-related fMRI method, the full time course of a BOLD response to a briefly presented stimulus lasts about 15 seconds for the robust positive response. 1.3 Literature Survey This section starts with providing a review on traditional techniques for inferring brain con nectivity. Brain connectivity is crucial to understanding how neurons and neural networks 4 function since neural activity are largely governed by connectivity. There are three types of brain connectivity patterns currently used in modeling brain functioning, with each of them defining the pathways between brain regions from a different perspective. Structural con nectivity, also known as anatomical connectivity, refers to a network of brain regions linked with structural (synaptic) connections, whose associated structural biophysical attributes are represented by parameters such as synaptic strength or effectiveness [7]. Structural connectivity is characteristized with relative stability over short time scales (seconds to minutes). In contrast, functional connectivity defines connections between brain regions by measuring deviations from regions’ statistical independence [8]. Unlike structural connec tivity, functional connectivity features high time-dependability.The last type of connectivity patterns is effective connectivity, which may be viewed as a combination of structural and functional connectivity. Effective connectivity describes networks of directional effects of one neural element over another, which can be inferred through systematic perturbations of the system or through time series analysis [8]. For each technique, a brief discussion regarding its theoretical basis is presented. In addition, this section also describes popular group-analysis methods necessary for meaningfully interpreting results from connectivity models. 1.3.1 Structural Equation Modeling Assuming that brain connectivity can be inferred from changes in the covariances of activi ties among various regions that are anatomically inter-connected, SEM employs the casual structure of the system to account for the observed pattern of correlation [9]. A simple version of “neural system” is defined in Fig.l.2, which is composed of four ROTs that are linked by anatomical connections. Its casual order is usually represented by a set of path equations in Eqn.l.l, in which the inter-ROT correlations are expressed as the sum of the compound anatomical connections. The path coefficients quatifying the anatomical con nections can be obtained through algebraic substitution. As a result, solving connectins 5 from path equations will incur high computation burden when dealing with high compli cated systems used to model brain functioning. Compared with path equation modeling, SEM (Eqn.1.2)represents the casual order of the system using a set of structural equations specifying the components of the variance of each ROT instead of the components of each correlation coefficients. SEM aims to inferring those anatomical connections as a solution to structural equations to minimize the differences between observed covariance relationships, and such solution is usually obtained mathematically from least-square or iterative mathods [10]. A. Path Equations: The following path equations indicate how the correlation can be decomposed to solve for the path coefficients: YAB = V, (1.1) /AC = x+(vw), 7AD = (vy) + (xz) + (vwz), YBC = w+(vx), 7BD = Y + (wz) + (vxz), 7CD = z + (wy) + (xvy), where ‘y’s are correlation coefficients between ROIs; low-case letters represent path coefficients to be solved, and residue influences are obmitted for simplicity. 6 B. Structural Equations: SEM represents the variance for each ROT as a function of weighted variance of other ROTs plus some residual influences: A = ‘i/A, (1.2) B = vA+’çbB, C xA+wB+b0, D = yB+zC+bD, where A, B, C, and D are known variances for each ROT, and ‘s are the residue influences not measurable. The resulted path coefficients quantifying anatomical connections can have different signs (positive or negative) to measure the directions of the covariance relationships between ROIs [11]. The success of SEM depends on large sample sizes, in most cases exceeding 200 individual cases, and most of its estimation functions require a multivariate normal distribution of the analyzed data. In addition, SEM’s covariance-based nature tends to discard useful information represented by other characteristics of the data, such as mean, kurtosis, or skewness. 1.3.2 Multivariate Autoregressive Model MAR models characterize the connections between ROTs in terms of the historical influ ences of ROIs upon each other [12]. Unlike SEM, which only captures the instantanous cormections, MAR models quantify ROl’s dependence based on data with time lags. As an extension of the well-known autoregressive (AR) process, a MAR (p) model extends the univariate case to represent a multi-dimensional time vector Yr-, as a pth-order weighted lin ear combination of previous measurements Y_ plus residue E (Eqn.1.3). For each weight matrix A(i), the magnitude of each coefficient in it can be interpreted as the strength of 7 either self-connections (diagonal entries) or inter-ROT connections at certain time lag. Y =Y1A(i) + E, (1.3) where Y,. is the nth measurement of a d-dimensional time series, while Y_ (i = 1... ,p) are previous measurements back to (n — p). Each A(i) is a d x d matrix of weight coefficients. The exclusive dependence of SEM upon instantaneous correlation tends to discard tem poral information useful for statistical inference. In contrast, MAR models derives temporal information from the auto as well as cross covariance functions that characterize the tem poral features of connections. In addition, the temporal persistence presented in time series makes MAR models to always have analytic solutions that can be obtained from linear algebra when dealing with self-connections. As a result, MAR models involves less com putational cost than SEM to represent the same effects. The application of MAR models in characterizing networks of cortical activity associated with event-related tasks enables people to quantify the dependence among all possible combinations of pairs of regions in the model. Therefore, connectivity networks can he compared across different event-related tasks. However, the application of MAR is restricted to modeling linear relationship, which is unsuitable for brain analysis since the dynamic of brain activity is often described as a non-linear system [13]. 1.3.3 Dynamic Causal Modeling DCM treats the brain as a deterministic nonlinear dynamic system, which is subject to external inputs, and produces output responses [14]. By feeding the system designed exper imental inputs, and measuring resulted output responses, DCM can estimate parameters used to quantify effective connections in terms of coupling among hidden neuronal activ ities in ROTs. As a departure from conventional connectivity-learning approaches, DCM features two primary distinctions. Firstly, DCM accommodates the non-linear and dynamic 8 nature of brain functioning, which is either neglected or inproperly modeled in conventional approaches such as MAR or SEM. Secondly, inputs in DCM for connectivity inference are known, and are designed for experimental purpose while conventional approaches assume the inputs to be unknown and stochastic. However, DCM shares with other analysis of effective connectivity a common property that the results are experiment-specific and can not be applied to general case. In other words, DCM is not an exploratory method, in stead, it is only used to test the hypothesis particular to the experiment [14]. In DCM, there are five state variables for each ROT with the first state variable representing the neu ronal activity, and the rest four representing the hemodynamic model. Therefore, the DCM could be considered as a combination of two sets of equations: neuronal state equations and hemodynamic state equations in Eqn.1.5. A. Neuronal State Equations: Consider a DCM system with m inputs and n outputs, with each ROT having one output. The m inputs are designed perturbations that are usually represented as design or stimulus functions in conventional analysis. The induced output for all n ROTs, defined as the change of the neuronal state variables with time, is represented as the neuro-physiological influences that are caused by neuronal activities, perturbations as well as model parameters in Eqn.1.4: 8Z/8t = F(Z, u, 0), (1.4) where Z = [zi,...,zj refer to neuronal state variables for n ROTs, F(.) is a nonlinear function, u refer to external perturbations, and 0 are model parameters whose posterior density is required for inference. By employing bilinear approximation to DCM, model parameters 0 can be reduced into three sets: the first set corresponds the extrinsic influence of inputs upon iieu ronal activities, the second set corresponds to the intrinsic connections that relate the response of one ROT to the state of others, the third set corresponds to changes in the 9 intrinsic coupling resulted from inputs. A mathematical demonstration is provided in Eqn.1.5: = (A+uB)Z+Cu, (1.5) A = 8F/0Z, B3 = C = aF/au, where Jacobian matrix A represents the first-order connectivity between ROTs without perturbations, corresponding to the second set of parameters. B3 refers to the change in coupling induced by the jth input, corresponding to the third set of parameters. The first set of parameters are defined in C. Those three set of parameters are used to construct the system architecture, and define inter-ROT connectivity. B. Hemodynamic State Equations: The second set of equations in DCM consists of biophysical state variables for hemodynamic response, including a vasodilatory signal, normalised flow, normalised venous volume, and normalised deoxyhemoglobin cotent. Hemodynamic state variables are determined only by the neuronal state variable z for each ROT i: zj results in an increase in the vasodilatory signal s, that is subject to auto-regulatory feedback. Inflow f responds in proportion to this signal with concomitant changes in blood volume v and deoxyhemoglobin content qj. This process is mathematically expressed in Eqn.1.6: z — Icjs3 — — 1), (1.6) 8f3/8t = Si, r8v/3t = fi — v 1/cr rdq/8t = f3E(.f1,Pi)/Pi — v qi/vi, 10 where E(f, p) = 1 — (1 — )1/f defines the oxygen extraction, and i, ‘y, r, c, and p are the biophysical parameters. 1.3.4 Dynamic Bayesian Networks DBNs were first proposed as graphical models to describe stochastic processes in [15, 16, 17]. Similar to regular BNs, DBNs represent the conditional independence/dependence among ROTs as directed acyclic graphs (DAGs). In addtion, featuring solid basis in statistical theory, and flexibility in handling both categorical and continuous data, DBNs are capable of modeling dynamic non-linear relationships among ROTs, making it particularly suitable for fMRI data. Since DBNs that describe a stochastic process in terms of all its past his tory can become intractably complicated, first-order Markov and stationary assumptions are usually applied [18]. Consier a discrete Gaussian stochastic process with X1, denoting its status variable at each time point. The process can he fully specified by the conditional distribution P(X/X_1)and the initial distribution P(Z1). Since P(X/X_1) does not change over time, the DBN grach could be defined as two time-slices: X1andX2,as illustrated in Fig 1.3. Despite being successfully applied for group analysis in fIVIRT study [18], the performance of DBNs is contrained by some of their limitations. First of all, DBN structure is represented by directed acyclic graphs (DAGs), not allowing cyclic connections between nodes. Secondly, DAGs must be converted to essential graphs (EGs) before con nectivity analysis since several different DAGs could represent the same set of conditional independence relationships, while there is only one EG uniquely encoding the set of con ditional independence relationships . Thirdly, the computational costs of large DBNs may also be prohibitive, making exploratory analysis less practical. 1.3.5 Group-Analysis Methods Group analysis involves two fundamental factors: the common features shared by individual subjects within the same group, and the subject-specific features that vary among the group 11 [18]. The common features usually represent significant properties of the underlying group, providing irreplaceable information for group analysis. They are of great interest in our case, where the primary goal is to identify significant difference between patients and normal subjects. However, the importance of feature distinctions among individual subjects can not be neglected since they may indicate the presence of outliers, potentially biasing the overall group result. This argument is supported by many studies that found that fMRI activation patterns may considerably differ among subjects within the same group. Therefore, it is important to maintain a delicate balance between similarity and diversity among groups [18]. Current group-analysis methods can be categorized into three broad categories: a virtual- typical-subject (VTS) approach, an individual-structure (IS) approach, and a common- structure (CS) approach [18j. The VTS approach assumes that each subject within the same group shares the same connectivity network while performing the same task. This approach is usually implemented by pooling together subjects from the same group, and treating them as a single subject with larger size, which tends to increase sensitivity. Another way to conduct VTS is to average group subjects to reconstruct a virtual network representing the entire group, so that the signal to noise ratio (SNR) can be significantly increased. However, VTS minimizes the effects of inter-subject variability, it is applicable in cases where the subjects are homogeneous within each group, and the inter-subject variability follows certain pattern such as Gaussian distribution. In contrast, the IS approach allows maximum inter-subject variability by learning individual networks separately, and conducting group analysis based on those networks. IS focuses on the ability of subjects to perform the same function or yield the same output yet with structurally different networks. The CS approach provides a trade-off between VTS and IS by imposing the same network structure on individual subjects, but allowing the parameters to vary across subjects. Therefore, the Cs approach address group similarity at the structural level and inter-subject variability at the parameter level. 12 1.4 Research Objectives and Methodology In this thesis, we examine the fIVIRI data collected from subjects performing three-frequency tasks. For comparison purpose, subjects were grouped into three categories: NP (normal subjects), PA (patients without medication), and PB (patients with medication). The objectives of this research are: 1. to identify significant difference between normal subjects and PD patients in terms of effective connectivity among ROTs, and to verify the usefulness of current L-dope medication for PD based on those connections. 2. to compare the network differences between normal subjects and PD patients in terms of their asymptotic behaviors, and study the feasibility of intervention in normalizing PD symptoms. To achieve these objectives, we propose to use PBNs for modeling brain connectivity by contructing a PBN for each subject, and learning the inter-ROl connections represented as the influence of one ROT upon another. Then we apply the Analysis of Variance (ANOVA) on the obtained connections with group factor and frequency factor to explore the abnormality caused by PD in the connectivity network. In addtion, by comparing estimated mean values of connections across groups, we can justify the effectiveness of L-dope medication on PD patients. Finally, we study the dynamic behaviors of group subjects in the context of first-order Markov Process. We generalize transition matrixs as well as the steady- state distributions within each groupto describe the dynamic features that characterize that group. The results also enables us to further explore the feasibility of ROl intervention as an effort to normalize the PD network. The rest of the thesis is organized as follows: chapter 2 presents the details regarding the application of PBNs on fMRI data. The results obtained from different groups indicate that PBNs help to identify statistically significant connections thta may provide promising directions for furture PD research. In chapter 3, we focus on studying the dynamic properties 13 of PBNs of the same fMRI data, and conducting intervention to reveal some ROTs that may become the potential key for studying the underlying machnism of PD. Lastly, we conclude the work along with suggestions for future work in chapter 4. 14 Background: Discovering Brain Connectivity from fMRI Data Figure 1.1: The figure depicts general procedures of fMRI networking learning. The exper irnent presented in this paper also follows this guidelim-. Figure 1.2: The path diagram of a neural system with four ROIs (A, B, C, D), and the anatomical connecions are indicated by arrows. v, x, y, z, w are the path coefficients representing the influences between structures through the anatomical connections. Scan the brain when the subject is performing a task, such as D Squeeze bulb Push button Multi-Channel Time Courses Learn the functional connectivity from the fMRI Data 15 r — — — — — — — — — — — — — DBN II I ®GG I I I I I I I I _________________ _________________ I I II I I I U U • I I oiI Yt+1 .I II Ib — — — — — N — — — — — — Figure 1.3: The figure represents a first-order Markov process including three sets of vari ables: input variables U, hidden state variables X, and observed variables Y. Arrows from X1 and U1 to X2 are associated with transition distribution. 16 1.5 Bibliography [1] J. Jankovic, “Parkinson’s disease: clinical features and diagnosis,” Neurosurg. Psychiatr, vol. 79, pp. 368—376, 2008. [2] T. Benes and F. Carisson, “The discovery of dopamine,” Trends in Pharmacological Sciences, vol. 22, pp. 46—47, 2001. [3] C. Weiller and etc., “Role of functional imaging in neurological disorders,” Journal of Magnetic Resonance Imaging, vol. 23, pp. 84—85, 2006. [4] S. Ogawa et al., “Oxygenation-sensitive contrast in magnetic resonance image of rodent brain at high magnetic fields,” Magn Reson Med, vol. 14, pp. 68—78, 1990. [5] I. Shmulevich et al., “Probabilistic boolean networks: A rule-based uncertainty model for gene regulatory networks,” Bioinformatics, vol. 18, pp. 261—274, 2002. [6] N. Logothetis, “Neurophysiological investigation of the basis of the fmri signa,” Nature, vol. 412, pp. 150—159, 2001. [7] B. Horwitz, “The elusive concept of brain connectivity,” Neuroimage, vol. 19, pp. 466—470, 2003. [8] K. Friston, “Functional and effective connectivity in neuroimaging: A synthesis,” Human Brain Mapping, vol. 2, pp. 56—78, 1994. [9] A. McIntosh, C. Grady, and etc., “Network analysis of cortical visual pathways mapped with pet,” J. Neurosci, vol. 14, pp. 665—666, 1994. [10] A. McIntosh, F. Gonzalez-Lima, and etc., “Structural modeling of functional visual pathways mapped with 2-deoxyglucose: effects of patterned light and footshock,” Brain Res, vol. 578, pp. 75—86, 1992. [11] M. Goncalves and D. Hull, “Connectivity analysis with structural equation modeling: An example of the effects of voxel selection,” Neurolmage, vol. 20, pp. 1455—1467, 2003. [12] L. Harrison, W. Penny, and K. Friston, “Multivariate autoregressive modeling of fmri time series,” Neurolmage. [13] K. Miller et al., “Nonlinear temporal dynamics of the cerebral blood flow response,” Hum Brain Mapp, vol. 13, pp. 1—12, 2001. [14] K. Friston, L. Harrison, and W. Penny, “Dynamic causal modeling,” Neurolmage, vol. 19, pp. 1273—1302, 2003. [15] M. Mitchell et al., “Classifying instantaneous cognitive states from fmri data,” American Medical Informatics Association Annual Symposium Proceedings, pp. 465—469, 2003. [16] L. Zhang et al., Modeling neu’ronal interactivity using dynamic Bayesian networks. Advances in Neural Information Processing Systems 18, MIT Press, Cambridge, MA, 2006. [17] C. Rajapakse and J. Zhou, “Learning effective brain connectivity with dynamic bayesian networks,” Neurolmage, 2007, accepted. [18] J. Li, Z. Wang, and M. McKeown, “Dynamic bayesian networks modelling of fmri: A com parison of group analysis methods,” Neurolmage, vol. 41, no. 2, pp. 398—407, 2008. 17 Chapter 2 Probabilistic Boolean Network Analysis of Brain Connectivity in . . 1Parkinsons Disease 2.1 Introduction The introduction of non-invasive medical technologies such as functional magnetic resonance imaging (fMRI) and quantitative electroencephalogram (EEG) has enhanced understanding about human brain functioning in normal and disease states. Understanding effective brain connectivity, the neural influence that one brain regioll exerts over another, is increasingly recognized as important for brain function, and its impairment may be associated with neuro-degenerative diseases such as Parkinson’s disease (PD). Inferring effective connections between brain regions of interest (ROTs) requires a math ematical model suitable for complex, large-scale, and dynamical system computation, and several such models have been proposed in the literature. Models such as structural equa tion modeling (SEM)[19, 20j, multivariate autoregressive models (MAR) [21], dynamic causal modeling (DCM) [22], and dynamic Bayesian networks (DBNs) [23] have all been suggested as suitable for fMRI data analyses, and each model has its particular advantages and disad 1A version of this chapter has been published/submitted for publication. Z. Ma, Z. J. Wang, and M. J. McKeown, “Probabilistic Boolean Network for Inferring Brain Connectivity Using fMRI Data,” in Proc. IEEE mt. Couf. Acoustics, Speech, and Signal Proc. (ICASSP), pp. 457-460, 2008. Z. Ma, Z. J. Wang, and M. J. McKeown, “Probabilistic Boolean Network Analysis of Brain Connectivity in Parkinsons Disease,” submitted to IEEE Journal of Selected Topics in Signal Processing, 2008. 18 vantages. SEM assumes that brain function can be monitored by examining the covariances of activities among various regions that are anatomically inter-connected [20] [19]. The sig nificance of condition-specific effects on brain connectivity is obtained by estimating the magnitude of each path connection in each experimental condition [19]. However, the suc cess of SEM depends on large sample sizes, and most of its estimation functions assume a multivariate normal distribution of the data [24]. Similarly, the application of IVIAR is restricted to modeling linear relationships, which may be unsuitable for brain analysis since the dynamic of brain activity may be described as a non-linear system [25]. A bilinear input-state-output model, DCM, has been proposed for inferring brain connectivity that may supplement linear modeling methods. DCM treats the brain as a deterministic nonlin ear system that is subject to external perturbations. It introduces perturbations into the system, and measures the response to obtain parameters quantifying effective connections in terms of coupling among hidden neural activities [22]. Very recently, dynamic Bayesian networks (DBNs) have been introduced to discover brain connectivity in fMRI data analysis [23, 26, 27] due to their solid basis in statistical theory, and their flexibility in handling both categorical and continuous data. Dynamic Bayesian Network (DBN) models represent ROl time courses as variable nodes connected by directed arcs representing their probabilistic dependence. The current fIVIRI modeling by DBNs usually takes on the form of a first-order Markov process, including dependent rela tionships between nodes within the same time frame as well as one-lag temporal dependence, both of which are useful for classification and predictive inference. Theoretically, DBNs are capable of dealing with stochastic and non-linear properties of multivariate data, making it particularly suitable for fMRI data. In addition, DBNs are capable of handling incomplete data as well as uncertainty, and can offer a framework for combining prior knowledge and data [28]. Despite the aforementioned advantages of DBN modeling, there are nevertheless some concerns that may restrict its current suitability for fIVIRI analysis. First of all, though 19 DBNs are theoretically able to model nonlinear interactions between node variables, current DBN modeling for fMRI normally assumes continuous variables and a linear (Gaussian) assumption due to the limited sample size typical of fMRI data sets. Directed acyclic graphs (DAGs), which represent the DBN structure, must be converted to essential graphs (EGs) before connectivity analysis since there is only one EG uniquely encoding the set of conditional independence relationships [28]. Finally, the computational costs of large DBNs may also be prohibitive, making exploratory analysis less practical. To address these concerns, we propose the use of Probabilistic Boolean Networks (PBNs) [29] for the analysis of fMRI data. PBNs are the stochastic extension of the well-known Boolean network framework proposed by Kauffman [30] and have proved successful in the bioinfornatics field, eg. gene expression analysis. PBNs combine the rule-based properties of Boolean networks and the probabilistic nature of DBNs. When fIVIRI experiments are performed in typical block fMRI designs, it is assumed that the hemodynamic response within a block is saturated, making the linearity assumption less tenable and also a binary representation more appropriate. Compared with DBNs, PBNs demonstrate computational simplicity [31] allowing a practical means to infer nonlinear relationships between brain ROTs. Like DBNs, PBNs describe the dynamics of the studied system in the probabilistic context of a Markov chain [29]. Therefore, independent PBNs and binary-valued DBNs can be converted between each other, representing the same joint distribution over their common nodes, as long as PBNs share the same initial distribution as DBNs. However, the mapping between PBNs and DBNs is not strictly one-to-one. Different PBN function sets can have the same probabilistic structure, making the PBN formalism a probabilistically redundant alternative for DBNs [32]. The conversion between DBNs and PBNs allows each one to benefit from the advanced tools originally designed for the other. For example, controlling the steady-state distribution by means of perturbation or intervention used in PBNs [32] could be potentially used in DBNs where no such methods has been introduced before. The sub-network projection originally developed for PBNs enables one to focus on 20 a small number of interesting nodes in DBNs. Meanwhile, the process of learning DBNs, consisting of model selection and parameter estimation, could also be adopted for variable and predictor selection in PBNs. Nevertheless there may be some features of PBNs that may affect their suitability for application to fMRI data analysis. For example, the necessary binarization preprocessing that must be first performed for PBNs may be less suitable for fIVIRI data obtained from event-related behavioral paradigms. We do not see the PBNs and DBNs as being in direct competition, and both of them have shown to provide good performance in biomedical applications [31] with the suitability of either method being application-dependent. The focus of both DBNs and PBNs on in teractions between nodes fits our interests in modeling brain functional connectivity, which, by definition, represents statistical dependencies among brain regions. This paper aims at introducing PBN modeling to fMRI analysis, and focuses on justifying its usefulness in inferring brain connectivity. Brain connectivity learning usually focuses on extracting connection features shared by subjects within a group and identifying significant feature difference between groups (e.g. subjects with disease vs normal subjects). Current group-analysis methods can be catego rized into three broad categories: a virtual-typical-subject (VTS) approach, an individual- structure (IS) approach, and a common-structure (CS) approach [28]. The VTS approach assumes each subject within the same group shares the same network structure. In con trast, the IS approach allows maximum variation between models by learning a network for each subject separately. Finally, the CS approach provides a trade-off between VTS and IS by requiring that each subject’s data is modeled with the same network structure, but the network parameters are allowed to vary. Since the IS has been demonstrated to provide more consistent performance in group analysis for PD patients in [28], we will adopt the IS approach for this paper. This paper is organized as follows, We describe the PBN framework for brain effective connectivity and group analysis in Section 2.2. A case study using fIVIRI data from normal 21 and PD subjects performing a motor task at three progressive levels of difficulty is discussed in Section 3.3. 2.2 Methods In this section, we present a PBN-based framework for inferring connectivity between ROTs. This framework not only quantifies the connections between brain regions, offering insight into mechanism of brain activity, but also enhances visualization. The PBN framework primarily consists of three components: the preprocessing of fiVIRI data, the inference of subject-specific PBNs, and finally group analysis of the individual PBNs. 2.2.1 Preprocessing of fMRI data The raw fMRI data first go through 3D motion correction and slice timing and then are further motion- corrected with Motion Corrected Independent Component Analysis (MCIMCA), a computationally expensive, but highly accurate method[33]. ROTs are then manually drawn on the high-resolution structural scans by a trained technician using a published atlas [34]. In this case 18 brain regions corresponding to the motor system were drawn. The resultant binary masks are then desampled to the fiVIRI resolution and then used to extract the voxels (and their associated time series) associated with each ROT. The proposed PBN analysis requires three important preprocessing steps: denoising, voxel selection, and binarization. Since a PBN is a discrete-valued probabilistic model, its input data is confined to binary vectors, in contrast to BN models that typically use continuous-valued data.The overall procedure combining those preprocessing stages is pro vided in Fig.2.1. fMRI Denoisirig: Most current denoising methods focus on building a reasonable model for serially correlated fMRI data. A traditional model for fMRI time series y of each voxel 22 is linear regression model [35][36]. y=X,3+c, (2.1) where X is the design matrix containing the task-related paradigm (and possibly other co variates), and 3 is the vector time-invariant regression coefficients, while Ct S the residual signal (the noise). Since fMRI residuals usually demonstrate an uneven spectrum distri bution with a large amount of power in the low frequency band, they cannot be naively modeled as white noise. There are two ways to solve this problem, either to whiten the noise by modeling the residual series Pt [35] [37], or to color it by convolving the Pt with a smoothing kernel [38]. For noise whitening approach, the residual series are often modeled as a AR(p) process. Pt = iPt-i + t, €t E N(0, 2) (2.2) Noise whitening attempts to remove the underlying autocorrelation by an iterative least square fit [39]. Similar to whitening, coloring replaces the unknown autocorrelation form with a known form by imposing a smoothing kernel, and is possibly more robust to bias caused by possible mis-specification of the form of residual autocorrelation. Other ap proaches, such as Independent Component Analysis (ICA) [40] and Canonical Correlation Analysis (CCA) [41] have also been proposed as efficient ways to extract the underlying components that are task-related from other noise components. Both methods are based on the assumption that each voxel time series x(t) (j = 1. . . n) is a linear combination of components s(t), (i = 1 n). xi(t) = aiisi(t) +... + ais(t), x2(t) = a2isi(t) + .. +a2s(t), x(t) = aisi(t) +... + as(t). 23 Each s(t) belongs to either ideal task-activated component, low-frequency drifts or other noise sources. ICA aims at obtaining statistically independent components s(t) measured by Gaussianity, while CCA focuses on finding components s(t) with maximum autocorrelation [41]. Though we investigated both ICA and CCA as denoising approaches, we note from empirical observation in our study that simply removing the low-frequency drift from each time series with a linear detrend operation resulted in the most consistent results. Voxel Selection: Given the fact that our study focuses on region-specific interactions, choosing appropriate voxels within each ROT to represent the overall behavior of that ROT may significantly affect the PBN analysis. There have been many ways proposed for voxel se lection including peak selection, eigenvariate selection, and average selection [42], suggesting that there exists no universal standard on how to choose a subset for optimal performance. Our approach for voxel selection combines several previous methods: we first use standard linear regression to obtain activation t-statistics for all voxels in the brain, and then find the one-tailed critical value c at a 95% level (uncorrected). Those voxels with t-statistics above the critical value are chosen as activated ones. To distinguish activated voxels from non-activated voxels, we further reconstruct the data by zeroing the time courses of all in-activated voxels. For each ROT, instead of simply taking the mean of all activated voxels within it we followed a previous approach [22] building a cubic with unit length centering at each activated voxel, and replacing its original time course with the mean of time courses of those voxels inside that cubic before extracting the overall time course that represents the activity of the entire ROT. Binarization: As our fMRI data are based on block designs where a task is done re peatedly within a block so that the hemodynamic response is saturated, we feel that the binarization process is appropriate in our analysis and the resulting binary data will be used as the inputs to PBN modeling. There have been several methods proposed to detect brain activation elevated by tasks, they are either based on a general linear model [43], or based on HMM [44] [45]. The binarization method proposed in [43] formulates the fMRI 24 time series ‘y as follows: W’y = Ky + W(X/3 + Hj + ), (2.3) where X accounts for the behavioral paradigm, H models confounding effects, and K is a discrete cosine transform matrix adjusted for low-frequency trends. The “pre-whitening” matrix W is generated from the estimated intrinsic autocorrelation of system noise c pre viously described in fMRI Denoisirig. The level of voxel activity for each scan due to the behavioral paradigm is determined by comparing the data, adjusted for known confounds, with a given threshold. If the adjusted data exceeds the threhold, then the level is activity is one, and is zero otherwise. R—_W7—KI WH,A=I(R>c*a), (2.4) where a2 is the variance of WE, c is the task-specific threshold, and A is the resulted binary sequence. We applied such approach upon our data, however, the subsequetial group analysis domenstrate poorly consistent results. Therefore, we adopt a feasible binarization method which can quantize real-valued fMRI data into binary values (on-state and off-state). General data-driven binarization methods, like the “big jump” method [46], aim at finding a critical threshold, such as the largest difference between two consecutive values in the sorted data. Since different ROIs reveal different levels of measurements, our binarization method is conducted separately on each ROT. Our preliminary investigations utilizing the “big jump” method were largely unsuccessful, since the largest difference often happens in the beginning of our sorted data, which makes the resulted binary data almost have no “off” states. As a result, we utilize another approach recently proposed [47], which focuses on capturing the first—order statistical trend in the fMR,I data. We suggest that it is reasonable to model the BOLD response as a hidden Markov model (HMM) process with parameters determined by the behavior paradigm. The key elements necessary to form a complete HMM include a set of observations Ot, which, in our case, are the time courses of the ROIs, a 25 two-state (1: off, 2: on) transition matrix A, and observation probability distribution B for each state (i = 1, 2), which is simplified as a Gaussian Distribution, and an initial state ir. Since we assume that the task starts from the rest condition, the initial state r could be set as 7t1 = 1, and 2 = 0. The mean and variance ci, of observation probability distribution for state i are determined by the structure of the paradigm p and the observations Ut. Take the off-state for example, I , , = 1 (O — (2.5) Poff tEpoff V Poff tEp0ff where Poff is the off-state fraction of the task paradigm. For the transition matrix A, we try a couple of values within an acceptable range, and especially pay attention to those results shared by different choices of A. We then employ Viterbi Algorithm to obtain the most likely state sequence, which is used to represent the fMRI voxel time course in the binary domain. To test whether our HMM—based binarization approach introduces extra correlation into the data, we replace the observed data with random Gaussian noises of the same power. The resulted binary sequences indicate no strong similarity with the experimental protocol, implying that the binarization does not force its ouput to be strong correlated with the protocol. 2.2.2 Introduction of PBN As an extension of a Boolean network (BN), a PBN models the interactions among different nodes of the system through logical rules with stochastic properties. Our introduction of PBN first starts with a brief review of a BN. A BN consists of a set of nodes (xi(t), ..., and a list of Boolean functions (f1t), ..., f(t)), where a Boolean function is defined as a combination of logical gates (AND, OR, NOT) used to regulate the connections among nodes. Assuming that each node of the studied system has a switch-like behavior, moving 26 between on and off states in the temporal dimension, and that nodes are not independent with each other, BN models the system as a highly-regulated binary network, where the activity of each node x(t) is governed by nodes known as parent nodes through one of those Boolean functions. Values of nodes are synchronously updated by corresponding functions based on the current state of parent nodes which the functions take as input variables. We construct a triple-node BN example represented by a truth table shown in Table.??. Functions targeting each node are illustrated in figure.2.2. For example, the value of xi(t+1) is governed by x1 (t), and x2 (t) through Boolean function: xi(t+ 1) = xi(t)flx2(t . (2.6) The dynamic of a BN is deterministic, implying that the temporal evolvement of the entire system completely depends on the initial state. Starting from the initial state, the system always follows a static transition mechanism regulated by those binary functions, and ultimately settles to a so-called attractor state, from which the system cannot move. By definition, all possible states which end up with the same attractor state are grouped as the basin of that attractor. For instance, we deduce from the simulation that, in Example of Table I, there are three disjoint basins of attraction: 000 forms a basin of attraction of its own, 001 and 010 are mutual attractors, which circulate between each other. The rest of the possible states are included in the basin of attractor 011. To avoid the potential risk caused by the deterministic rigidity of BN, Shmulevich et.,al proposed adding a probabilistic setting to BNs [29]. The basic idea of a PBN is to accom modate more than one Boolean function for each node, allowing the model to have more flexibility. Typically, for each node x,, there is a set of Boolean functions known as predictors, among which one is chosen to update the target node’s state according to its probability p at each transition. F = {f} ,j = 1 , m, Pr(x1(t+ 1) = f(t)) = p, (2.7) 27 where m is the number of predictors for x. As shown in Table.2.2, we could extend our previous BN example in the context of PBN by incorporating more functions with assigned probabilities. The PBN function set for node x1 is illustrated in Fig.2.3 The dynamic behavior of PBN is represented by a Markov chain whose transition matrix is defined by all possible Boolean functions as well as their probabilities. Of particular interest is the long- run (asymptotic) behavior of the PBN and if there exists a steady-state distribution the system will lead to regardless of its initial distribution. We study the asymptotic distribution of a n-state PBN in the strict context of a homogeneous Markov chain by first defining a stationary distribution ii-, so that the possibility for the system to be in state i at time any time r is always equal to ir, if the initial distribution is ir. (2.8) with p being the r-time transition probability from state i to state j. Whether a PBN system has a single or multiple stationary distributions depends on its reducibility. If a PBN is irreducible, implying that it has a communicating state space where any state can possibly be reached by any state, then it is assigned a unique stationary distribution. In contrast, a reducible PBN indicates that there exist irreducible sets of states so that no state outside them can be reached from inside. Each irreducible set of states has its own unique stationary distribution that can be extended to a stationary distribution for the overall system by setting the probability outside those sets to be zero. However, a unique stationary distribution does not guarantee the existence of a stead-state distribution. A PBN is said to have a steady-state distribution if and only if it is irreducible as well as aperiodic (ergodic). In other words, the stationary distribution is also the steady-state distribution if it satisfies the following condition: lirnp = ‘ir, (2.9) 28 Given the initial uniform joint distribution (iro = [1/8, ... 1/8]) in our example, we could obtain the asymptotic distribution ir = [0.15,0.04,0.15,O.11,0.1,0.16,0.1,0J9] which indi cates that this PBN system is irreducible. This is also supported by the state-transition diagram in figure 2.4 showing there exists no irreducible sets of states. In addition, changing the initial distribution always gives us the same asymptotic distribution, which verifies the existence of a steady-sate distribution since the state transition probability indicates that each state is periodic with period 1. 2.2.3 Learning PBNs Several learning approaches for PBN have been explored, like discussed in [48, 49]. Here we followed the approach specified in [29]. This learning process was implemented by the BN/PBN toolbox from Lhdesmki and Shmulevich. Our fIVIRI data is formatted into 18 x 130 matrics after pre-processing, where 18 is the number of pre-defined interested ROIs denoted as x (i = 1, 2, ...18), and 130 is the number of measurement steps. Because of the relatively large number of nodes compared with measurement steps, we limit the maximum number k of input variables of each predictor to be five, i.e (k < 5) to reduce system complexity as well as the search space for predictors. The optimal predictor for each possible input variable combination is obtained through greedy search. However, even if we limit the number of input variables, the computational cost is still quite high in our case since each Boolean function consisting of five inputs could be represented by a truth table with 2 entries, and the total number of Boolean functions to be searched is up to 232. For the same reason as to simplify modeled system, we also allow the maximum number of predictor functions each node can have to be four, and those predictor functions are ones with the highest Coefficient of Determination (COD) [50] [51] among optimal Boolean functions of all possible variable combinations. For each possible input variable combination 1 (1 = 1,2, ..., C), the COD of one of its Boolean functions f is used to measure the relative decrease in error by estimating 29 the target node x through that function instead of the best constant estimate: (2.10) where e is the error of the best constant estimate of x, and c(x, f) is the error of the estimate of x, by f. Given k,each possible variable combination has an optimal Boolean function with the highest COD among all candidates for that combination. However, we often observe that there are several combinations whose optimal Boolean functions share the minimal COD across all combinations. If the number of combinations with minimal COD is larger than 4, as previously mentioned, it implies that the current k is not able to precisely locate the best predictors. Therefore, we keep reducing the value of k to narrow our search range until we obtain no more than four predictors for each node. The probability p assigned to each predictor is determined as the weighted percentage of its COD across all predictors: p= (2.11) q=1 Upon obtaining predictors, we need to identify the impact of one node on another, which is measured in terms of influence value [29]. A large influence value means that the corre sponding node play a very important role in determining the target value. Consider node x with a function set F including all predictors f and corresponding probabilities p (j = 1,...,m).The influence of Xk upon x, i.e Ik(x), is defined as the sum of the influence of Xk upon each predictor J: 1k(X) = 1k(f)P, (2.12) where Ik(f)) is defined as: Ik(f) = Pr{ = i} = Pr (f(x) f(xk)), (2.13) 30 where f)(xk) is the same as f(x) except the value of Xk is toggled. This influence calculation is applied to all 18 ROTs, thus forming a 18x18 matrix G with each entry being the influence of x on x3. The influence matrix G enables us to obtain a complete connectivity network. This network can be demonstrated as a weighted directed graph, with each arrow representing the influence of one ROT to another. The value for each arrow is the magnitude of the influence. 2.2.4 Group Analysis for Inter-Subject Variability So far we have learned how to derive a PBN for each individual subject. However, to mean ingfully extrapolate PBN results from one subject to an entire population (e.g. PD patient group) first requires methods to meaningfully integrate results from individual subjects within a group and to rigorously compare PBNs across groups. To address this inter- subject variability issue, a common and critical challenging problem in many biomedical studies, we explore group analysis for each influence connection contained in G individually through the Analysis of Variance (ANOVA) [52]. ANOVA can test the effects of multiple factors on the data, although identifying sig nificantly different connections between PD sub.jects and normal subjects is our main goal. Nevertheless, given the fact that our experiment involves subjects squeezing a bulb at three different frequencies, we also consider the effects of frequency. As a result, we incorporate two factors in the ANOVA analysis: frequency and group (e.g. normal, PD on medication, and PD off medication), and consider these two factors independently. The ANOVA will return a p-value for each factor, and a small p-value below certain threshold (0.05 in our case) suggests a significant difference in the data as a result of the corresponding factor. The ANOVA test only offers general information about whether the means of data are sigiiificantly different. The multiple comparison method provides US with graphs mdicatiiig the estimated range of means for each level within the current factor under a predefined confidence level. Two means are significantly different if their intervals are disjoint, and 31 are not significantly different if their intervals overlap. In the present case, we have three different groups: Norm(normal), Poff (patients off medication), Pon (patients on medica tion). The existence of a large gap between the Poff bar and Norm bar may imply that PD causes a detectable abnormality in the studied connection. Furthermore, the multiple comparison method also helps us to evaluate the effects of L-dopa medication by comparing the differences between Norm and Poff with the difference between Pon and Poff [28J. For example, If (MNorm — i’ip.ff) shares the same sign as (M0 — ‘1poff), it suggests that the medicine alleviates the PD symptoms by pulling the mean interval of patients back close to that of normal subjects. 2.3 fMRI Case Study 2.3.1 fMRI Data Collection The study was approved by the University of British Columbia Ethics Board and all subjects gave informed written consent prior to participating. Ten volunteers with clinically diag nosed PD participated in the study (4 men, 6 women, mean age 66 years, 8 right-handed, 2 left-handed). All subjects had mild to moderate PD. We also recruited ten healthy, age- matched control subjects without active neurological disorders (3 men, 7 women, mean age 57.4 years, 9 right-handed, 1 left-handed). Exclusion criteria included atypical Parkinson ism, presence of other neurological or psychiatric conditions and use of antidepressants, sleeping tablets, or dopamine blocking agents. After completing the experiment in an off medication state, subjects were given the equivalent to their usual morning dose of L-dopa in immediate release form. They then repeated the same tasks post-medication following a 1 hour interval to allow L-dopa to reach peak dose. We employed an fMRI, ROT-based pre/post medication study of sinusoidal force produc tion at 3 frequencies (0.25, 0.5 and 0.75 Hz). Frequencies were comparable to prior tracking tasks and we confirmed that this range of frequencies could be accomplished by PD subjects 32 off-medication in a pilot study. Subjects lay on their back in the functional magnetic resonance scanner viewing a com puter screen via a projection-mirror system. All subjects used an in-house designed response device held in their left hand, which was a custom-built MR-compatible rubber squeeze-bulb connected to a pressure transducer outside the scanner room. They lay with their forearm resting down in a stable position, and were instructed to squeeze the bulb using an isometric hand grip and to keep their grip constant throughout the study. Each subject had their maximum voluntary contraction (MVC) measured at the start of the experiment and all subsequent movements were scaled to this, so they had to squeeze at 5-15 % of MVC to accomplish the task. Using the squeeze bulb, subjects were required to control the width of an “inflatable ring” (shown as a black horizontal bar on the screen) in order to keep the ring within an undulating pathway without scraping the sides. Applying greater pressure to the bilib increased the width of the bar, and releasing pressure from the bulb decreased the width of the bar. To avoid scraping the sides of the tunnel, the required pressure was between 5 and 15 % of MVC. The pathway used a block design with sinusoidal sections in three different frequencies (0.25, 0.5 and 0.75 Hz) in a pseudo-random order, and straight parts in between where the subjects had to keep a constant force of 10 % MVC. Each block lasted 19.85 seconds (exactly 10 TR intervals), alternating a sinusoid, constant force, sinusoid and so on for a total of 4 minutes. The complete experiment specification is shown in Figure.2.6. Functional MRI was conducted on a Philips Achieva 3.0 T scanner (Philips, Best, the Netherlands) equipped with a head-coil. We collected echo-planar (EPI) T2-weighted im ages with blood oxygenation level-dependent (BOLD) contrast. Scanning parameters were: repetition time 1985 ms, echo time 37 ms, flip angle 90, field of view (FOV) 240.00 mm, matrix size = 128 x 128, pixel size 1.9 x 1.9 mm. Thirty-six axial slices of 3mm thickness were collected in each volume, with a gap thickness of 1mm. We selected slices to cover the dorsal surface of the brain and include the cerebellum ventrally. A high resolution, 33 3-dimensional Ti-weighted image consisting of i70 axial slices was acquired of the whole brain to facilitate anatomical localization of activation for each subject. The following ROT’s were drawn separately in each hemisphere, based upon anatomical landmarks and guided by the Talairach atlas in [53]: primary motor cortex (Mi), supple mentary motor cortex (SMA), prefrontal cortex (PFC), caudate (CAU), putamen (PUT), globus pallidus (GLP), thalamus (THA), cerebellum (CER) and anterior cingulate cortex (ACC). 2.3.2 Results and Discussions First, to illustrate the learned PBN for individual subjects, an example of PBN for a normal subject performing the squeezing task with the left hand at the highest frequency is shown in Fig. 2.5. From this figure, we are able to visually identify important core nodes (ROTs) that play important roles for performing the corresponding fMRI task. For the purpose of backbone ROl identification, we employ a simple way here by detecting the ROTs that represent the largest topology change (the number of edges being deleted) if the node is removed from the network. It is noticed that ROTs RPUT, RTHA, RCER, L.M1, and R_PFC are identified as core regions. The R_CER, L_Mi regions represent an ipsilateral motor network that appears recruited at higher frequencies. The subcortical regions of RPUT, R.THA have been implicated in prior studies as being important for the scaling of movements [54] [55] [56]. Before we examine the group analysis results through the ANOVA test of individual connections, we carried out some preliminary work to justify the correctness of obtained influence table by comparing every entry with its corresponding partial correlation coefficient based on the assumption that high partial correlation coefficient usually implies strong influence. The test results confirm that our PBN-deduced influence connections are largely related to corresponding partial correlation coefficients, indicating the robustness of our PBN results. We further enforce our argument by randomly shuffling the temporal ordering 34 of original fMRI data, and counting the number of significant connections of each data set above certain threshold. The results show that randomly shuffled data results in far fewer significant connections, compared with the original, unpermuted original data. It is worth mentioning that a more rigorous assessment of statistical significance based on surrogate /permuted data will be further investigated in the future. We identify all influence connections whose values are significantly different due to ei ther group membership or frequency factors. Since the waveforms for all ROTs may likely be affected by the experimental protocol (the expected behavioral paradigm), it is natural to treat it as an input node in the network. Thus we incorporate the protocol into the PBN sys tem as an extra node (but zero all connections towards it) to represent the underlying model better and to reduce possible biased significance of learned connections between brain ROIs. Among all possible 361 connections, 16 are identified as statistically significantly different due to the group factor, as listed in Table.2.3. Another purpose of group analysis here was to determine the overall group effect of medication on PD patients and to identify the specific connections between R.OIs upon which the medication has significant impact. In order to do that, a feasible method is to compare the sign of the mean difference between normal subject and patients with the sign of the mean difference between patients with medication and pa tients without medication. If those two pairs share the same sign, it implies that L-dopa does provide effective treatment by normalizing the connection. Tahle.2.4 indicates that L-dopa medication has a normalizing effect on 10 connections among all 16 identified connections, including LPUT —* LM1, LPUT —* LSMA, LCAU —* RifHA, LCAU —* RGLP, LifHA —* RCER, LifHA —* LPFC, &CER -+ RCER, RM1 -* RM1, RSMA — RSMA, and LACC —* RCER, where influence magnitude in patients after medication approaches that of normal subjects (figure 2.7). Thus, the medication improves connec tivity in subcortical structures that are know to be dysfunctional in PD. However, there are several cases (L.PUT —* LM1, LCAU —* RGLP,et.al) where the medicine seems to “overshoot”( Fig.2.9(c)). These cortical areas may be effectively overdosed when the 35 medication is given to correct for the primary subcortical dopamine deficiency. In contrast, L-dopa medication seems to exacerbate the already abnormal connectivity in 2 connections including RPUT —* RPUT, RPUT —* RM1, and such phenomenon requires further clinical research. Those two connections are illustrated in Fig.2.9(a) & (b). Two cases of normalizing connections projected onto the brain anatomy are visualized in figure 2.8 while cases of deterioration and “overshoot” are visualized in figure 2.10. We also examined the effects of different frequencies by examining the frequency factor in the ANOVA test. In total, 16 connections are identified as having statistically significant modulation at various frequency levels (High (0.75Hz), Med (0.5Hz), and Low (0.25Hz)), as listed in Table 2.5. As listed in Table.2.6, a further multiple comparison test reveals that, among those 16 significant connections, L_THA —* R_SMA is particularly interesting because of its monotonic trend along frequency levels, as indicated in figure.2.11.(a). We also employ the multiple comparison method to find out whether different groups within the same frequency also indicate significant difference between normal subjects and PD subjects off of medication, and whether such difference is improved by medication. At the maximum frequency, we find that only RM1 — LPFC indicates mild improvement after medication. Meanwhile, we find three out of four significant connections at the middle frequency are due to the difference between normal subjects and PD subjects off medication, and all of them are greatly improved after medication, implying that medication has its maximal effect at the middle frequency, as shown in Fig.2.11.(b),(c),(d). At the lowest frequency of squeezing, no obvious medication effects are detected. This may reflect the fact at the middle frequency, compensatory activity within the normally-recruited network is maximal, as opposed to higher speeds were new, secondary networks are recruited in PD. Connections indicated in figure.2.11 are projected onto the brain anatomy for better visualization in figure.2.12. 36 2.4 Discussion and Conclusion We have presented a Probabilistic Boolean Network framework for inferring brain func tional connectivity and demonstrated its utility for examining brain ROl interactions in PD studies. The detected connectivity differences between the PD group and the control group suggest that brain connectivity between certain regions can be a sensitive marker for PD and may provide a better understanding of the underlying mechanisms of PD. The core regions detected, such as RCER, LM1 RPUT, and RTHA represent known regions such as an ipsilateral motor network or subcortical structures which are important for movements of increasing speed. There are nevertheless some limitations that may restrict the suitability of the proposed PBN approach. PBN modeling requires representing the fIVIRI voxel time course in the binary domain. Because of this binary assumption, the proposed technique should be re stricted to block designs, as it is assumed in block designs that the hemodynamic response is saturated within a block, and a binary representation is therefore appropriate. However, we do not think that this is necessarily a gross limitation. There has been a vigorous debate on the theoretical benefits of block-design vs event-related designs. For patient populations, the block design, because of its ability to reduce the total time within the scanner is often preferred. Statistical parametric maps are usually thresholded at a given level of signif icance which results in a binary mask of regions considered signficantly activated by the behavioural paradigm. Event-related designs are preferred if it is expected that the control and experimental task(s) result in heavily spatially overalapping activation patterns, and the relative amplitude of activation is important. In this case, the binarization procedure necessary for subsequent PBN analysis would be inappropriate. A potential risk of dealing exclusively in the binary domain is that with the relatively short time courses typical for fIVIRI data sets, random binarization patterns may appear signficantly correlated. This makes it more suitable to use the PBN to explore group behavior, as false connections are unlikely to be propagated through entire groups. For 37 real world applications, a desirable structure-learning method needs to account for the error rate of the claimed discovered connections. Thus, it is desirable for learning algorithms to control the error rate of the relationships discovered from a limited number of observations. However, as suggested by [31], the accuracy and reliability of PBNs are dependent upon not only model selection but also on inference schemes. The current exhausive-search learn ing method is time-consuming, which involves computation of up to 232 candidates in our case. Therefore, a more efficient algorithm will allow the asseessment of null data sets which will ultimately improve the accuracy of learned PBNs. Further work is needed to fully explore the domain over which PBNs are most suitable. Indirect connections, representing the influence of one ROl upon another through indirect roiltes need to be included for inference study since they show the total influence of one ROT and whether this influence is modified at any level of the system [57]. The current structure-learning algorithm in PBN could be modified to explicitly control the error rate (e.g. the false discovery rate) of the claimed connections in the network. The optimal way to combine results from individual PBN models also needs to be explored. Should the data from different subjects within a group be pooled, and thus enforce a common structure on each group, or should each individual structure model be computed and then the structures be combined in some way? [31] Despite these current shortcomings, the results presented here suggest that the PBN approach warrants further investigation as a complentary method for assessing functional connectivity between brain regions. 38 fMRI Raw Data Detrend T-rnapping F—- - Comparison of Comparison of t-value with t-value with Lshold threshold Voxel VoxelAdjustment Adjustment ROl Pre-processing • ROl Pre-processing ________ ForROll ________ ForROll8 Voxel Voxel Selection Selection ROl Mean ROl Mean Determination Determination HMM HMM Binarization Binarization Binary Vectorfor ROl I - BnaryVectorforROI 18 Probabilistic Boolean Network Figure 2.1: R,aw fMRI data are first detrended to remove noise, and standard linear regres sion methods are used to estimate a t-statistic associated with the likelihood of activation by the behavioral paradigm. The data are then divided into ROTs, and each ROT undergoes the same voxel-selection process and binarization procedure. The end result is a binary matrix with dimensions of 18 by 130 (the number of ROTs vs the number of time points in the experiment) — the input to the PBN model. 39 Table 2.1: Truth Table for BN Example x123 fi f2 f3 000 0 0 0 001 0 1 0 010 0 0 1 011 0 1 1 100 1 1 1 101 1 1 1 110 0 1 1 111 0 1 1 This BN system has three nodes (x1, x2, x3) with each having its corresponding binary function for value updating. Note that the number of input variables for each binary function varies. 40 Figure 2.2: This figure visualizes the logical function for the three-node BN system. 41 P1’ Figure 2.3: This figure demonstrates a network of predictors for our PBN example, which implies the stochastic properties of PBN’s dynamic behavior. 42 Table 2.2: Truth Table for PBN extended from BN example :1 :1 :2 :2 :3X123 ii J2 J1 J2 ii 000 0 1 0 1 0 001 0 1 1 0 1 010 0 1 0 1 0 011 0 1 1 1 1 100 1 1 1 0 1 101 1 1 1 0 1 110 0 1 1 0 1 111 0 1 0 1 0 p 0.6 0.4 0.5 0.5 1 This PBN expands previous BN by allowing x1 and x2 to have two predictors each. The bottom line indicates probability for each predictor 02 05 02 0.3 101 0.5 0.5 0. 0.2 100 001 000 0.2 0. 111 05 03 02 o 0: v.3 0 03 0.3 0 1010 Figure 2.4: The state transition diagram for PBN defined previously represents possible system states as vertices connected by directed arcs with transition probability, i. e.the probability for state ‘000’ to transfer to ‘100’ is 0.2, and et.al. 43 Figure 2.5: An example of the learned influences from the fMRI data for a normal subject performing a High frequency squeezing using Left Hand. The influnece diagram is a weighted directed graph with each arrow representing the influence of a ROl on another. The number next to the arrow is the magnitude of the influence. 44 Figure 2.6: The brain anatomy picture provides spatial locations of all 18 ROIs, among which RCER, LCER, LM1, RM1, RSMA, LSMA are of particular interest. Their fMRI time series are shown accordingly. The behavioral paradigm involves task performed with 3 frequencies. 45 Table 2.3: Connections showing significantly difference due to the group factor, Left Hand Connection P-value LPUT —* LM1 0.0156 LPUT —* LSMA 0.0335 RPUT -* RPUT 0.0011 RPUT —* RM1 0.0184 LCAU —* RifHA 0.0299 LCAU — RGLP 0.0311 LifHA - RCER 0.0417 LifHA — LPFC 0.0204 RCER -÷ RCER 0.0308 LLER —* RCAU 0.01 LM1 — LCAU 0.0117 R_M1 —* R_M1 0.0338 RSMA — RSMA 0.0269 RSMA —f LSMA 0.0218 RSMA —* RGLP 0.0187 L.ACC -÷ RCER 0.0022 Arrows —* denote influence connections with a time lag. The critical value here is .05. 46 Table 2.4: Performance of Medicine on PD, Left Hand Connection MNorm — Mp0ff M0 — Mp0ff LPUT —* L_Ml 0.0129 0.0863 LPUT —* LSMA -0.0356 -0.1251 RPUT — RPUT 0.1876 -0.2091 RPUT —* R_Ml 0.0745 -0.0277 LGAU —* RifHA 0.1132 0.0180 L.CAU — RGLP 0.0434 0.1066 LifHA —* RCER 0.0865 0.0049 LifHA -+ LPFG -0.0889 -0.0501 RCER —* RCER 0.0326 0.1886 LCER —* RCAU 0 0.0569 LM1 —* LCAU 0.0566 0 RM1 —+ Ru/h 0.1088 0.2294 RSMA —* RSMA 0.0397 0.1154 RSMA —* LSMA -0.0222 0.1042 RSMA —* RGLP -0.0052 0.0422 LACC —* RCER -0.1104 -0.1104 MNorm — Mp0ff denotes the difference of estimated means of influence coefficients between Norm and Poff. The same sign shared by MNorm — Mp0ff and MP0 — MPoff indicates that the medicine has positive effects on the connection, and vice versa. 47 CliCk Ot tflC 9 StP yot. t.eat to teat 2 -0.05 0 0.05 0.1 0.15 0.2 The OCIurCn means of gropsnana 2 goitcantly d,ffe,ent (b) LACC —* RCER Figure 2.7: The figure of each connection contains three bars with the top one being the estimated range of coefficient mean for Norm, the middle being the estimated range of coefficient mean for Poff, the bottom being the estimated range of coefficient xriean for Pon. A big gap between Norm and Poff indicates that PD cause obvious abnormality for this connection. Both figures show that the coefficient mean becomes closer to Norm after medication, which justifies current treatment. 2 3 1)4 -002 0 0.02 0.04 0.t 000 0.1 0.12 0 1 The column means of groups 1 and Pars significen try different (a) LifHA —* LPFC CliCk 00 the g.oup yOt te0ttotCt a 48 Figure 2.8: The brain anatomy picture provides specific locations of ROTs in each hemi sphere, with red circles denoting those connections where L-dopa has normalizing effects. 49 I I. Click o.. n.e go..p ye.. to test 3 S -o 00 0 0.05 0.1 0.15 0.2 No groups flats 0010 mn m eons signifloa,-.tly different from Group I (c) L.CAU -+ RGLP Figure 2.9: The figure of each connection contains three bars with the top one being the estimated range of coefficient mean for Norm, the middle being the estimated range of coefficient mean for Poff, the bottom being the estimated range of coefficient mean for Pon. A big gap between Norm and Poff indicates that PD cause obvious abnormality for this connection. Figures (a) & (b) show that the coefficient mean becomes further away from Norm after medication, which implies deterioration of PD symptoms. Figure (c) shows that the coefficient mean “overshoot” its Norm counterparts after medication. 2 2 0.1 0,2 03 0.4 05 06 0/ me column means of groups land S are significantly dIfferent (a) RPUT — RPUT Click o.. Ofl group you want to t000 04 -0.02 0 0.02 0.04 0.06 0.00 0.1 0.12 0,14 0. The colurt,n means of groups I arId Sore sigrliocan try diffCrCnt (b) RPUT -+ RM1 Click On One group you w000 to tCal 2 50 Figure 2.10: The brain anatomy picture provides specific locations of ROIs in each hemi sphere, with red circles denoting those connections where L-dopa either deteriorates PD symptoms or indicates “overshoot” effects. 51 Table 2.5: Connections showing significantly difference due to the frequency factor, Left Hand Connection P-value LP(JT —* RifHA 0.0123 LPUT -÷ RACC 0.046 RPUT —* RifHA 0.0328 RPUT —* RSMA 0.0498 RCAU — RCAU 0.0088 LifHA —* RSMA 0.0282 RCER — RCAU 0.0355 RCER —* RM1 0.0227 LCER —* LM1 0.0385 RM1 —* LSMA 0.0229 RSMA — LGLP 0.0361 LSMA —* R..M1 0.005 LPFC —* LifHA 0.0322 R.ACC - RGLP 0.0344 LACC —* LPFC 0.0237 LAGC — RACC 0.0347 Arrows —* denote influence connections with a time lag. The critical value here is .05. 52 Table 2.6: Performance of Frequency on subjects, Left Hand Connection MH — MM ML MM LPUT —* RifHA 0.0282 -0.0135 LPUT —* RACC -0.0604 -0.0229 RPUT — R,ff I-IA 0.0351 -0.0029 RPUT —* RSMA 0.0618 0.0173 RCAU - RCAU -0.007 0.0229 LifHA —* RSMA -0.0103 0.0610 RCER —* RCAU 0.0146 0.0008 &CER —* R_M1 -0.0789 -0.0382 L_CER —* LM1 -0.0253 0.0194 RM1 — L_SMA -0.0421 -0.0046 RSMA —* LGLP 0.0153 -0.0083 L_SMA —* R_M1 0.0711 -0.0567 LPFC —> LifHA -0.0585 -0.0639 RACC —* RGLP 0.0675 0.0208 LACC — LPFC -0.0068 0.0165 LACC — RACG -0.0009 0.0354 MH — MM denotes tne difference of estimated means of influence coefficients between high frequency and median frequency. The opposite sign shared by MH — MM and ML — MM indicates that the connection has monotonic trend along frequencies, and vice versa. 53 I I I I I -0.05 0 005 0.1 0.15 0.2 No groups have means signffcanlly different frsm np (b) L2”HA —* LM1 op I I I pa ph __________________________________________________ I I I I I I -0.04 -0.02 0 0.02 0.04 0.05 0.08 0.1 The means of grsops np and pa are significantly different (d)RSMA —> LfiFC Figure 2.11: Figure.(a) contains three bars with the top one being the estimated range of coefficient mean for high frequency, the middle being the estimated range of coefficient mean for median frequency, the bottom being the estimated range of coefficient mean for low frequency. It indicates that the coefficient mean decreases as the frequency increases. For figure.(b),(c),(d), each contains three bars with the top one being the estimated range of coefficient mean for Norni, the middle being the estimated range of coefficient mean for Poff, the bottom being the estimated range of coefficient mean for Fon. They illustrate that current medication imposes positive effects on those connections under median frequency. I) np pa ph -0.1 2 3 _1_ J .l J J -0.04 -0.02 0 0,02 0.04 0.06 0.08 01 The row means of groups land 3 are signffcantly different (a) LTHA -÷ RSMA op pa ph I I I I I -Ut -1105 0 0.05 Ot 0.15 0.2 0.25 The means of groups op and pa are signthcaotly different (c) RPFC — &ACC 0.25 II 54 Figure 2.12: The brain anatomy picture provides specific locations of ROTs in each hemi sphere, with red circles denoting those connections that either demenstrate monotonic trend along frequency levels or indicate significant difference between normal subjects and PD subjects off of medication within the same frequency range. 55 2.5 Bibliography [19] A. McIntosh, C. Grady, and etc., “Network analysis of cortical visual pathways mapped with pet,” J. Neurosci, vol. 14, pp. 665—666, 1994. [20] M. Goncalves and D. Hull, “Connectivity analysis with structural equation modeling: An example of the effects of voxel selection,” Neurolmage, vol. 20, pp. 1455—1467, 2003. [21] L. Harrison, W. Penny, and K. Friston, “Multivariate autoregressive modeling of fmri time series,” Neurolmage. [22] K. Friston, L. Harrison, and W. Penny, “Dynamic causal modeling,” Neurolmage, vol. 19, pp. 1273—1302, 2003. [23] L. Zhang et al., Modeling neuronal interactivity using dynamic Bayesian networks. Advances in Neural Information Processing Systems 18, MIT Press, Cambridge, MA, 2006. 1241 C. Wallenburg and J. Weber, “Structural equation modeling as a basis for theory development within logistics and supply chain management research,” Research Methodologies in Supply Chain Management,Physica- Verlag HD, pp. 171—186, 2005. [25] K. Miller et al., “Nonlinear temporal dynamics of the cerebral blood flow response,” Hum Brain Mapp, vol. 13, pp. 1—12, 2001. [26] M. Mitchell et al., “Classifying instantaneous cognitive states from fmri data,” American Medical Informatics Association Annual Symposium Proceedings, pp. 465—469, 2003. [27] C. Rajapakse and J. Zhou, “Learning effective brain connectivity with dynamic bayesian networks,” Neurolmage, 2007, accepted. [28] J. Li, Z. Wang, and M. McKeown, “Dynamic bayesian networks modelling of fmri: A com parison of group analysis methods,” Neurolmage, vol. 41, no. 2, pp. 398—407, 2008. [29] I. Shmulevich et al., “Probabilistic boolean networks: A rule-based uncertainty model for gene regulatory networks,” Bioinformatics, vol. 18, pp. 261—274, 2002. [30] S. Kauffman, “Metabolic stability and epigenesis in randomly constructed genetic nets,” J. Theoret. Biol, vol. 22, pp. 437—467, 1969. [31] P. Li et al., “Comparison of probabilistic boolean network and dynamic bayesian network approaches for inferring gene regulatory networks,” BMC Bioinformatics, vol. 8, pp. 118—129, 2007. [32] H. Lahdesmaki et al., “Relationship between probabilistic boolean networks and dynamic bayesian networks as models of gene regulatory networks,” Signal Processing, vol. 86, pp. 814—834, 2006. [33] R. Liao, J. Krolik, and M. McKeown, “An information-theoretic criterion for intrasubject alignment of fmri time series: motion corrected independent component analysis,” IEEE Transactions on Medical Imaging, vol. 24, pp. 29—44, 2005. [34] H. Damasio, “Human brain anatomy in computerized images,” Oxford University Press, Tech. Rep., 2005. 56 [35] E. B. et al., “Statistical methods of estimation and inference for functional mr image analysis,” in Magn Reson Med, vol. 35, 1996, PP. 261—277. [36] K. Worsley and K. Friston, “Analysis of fmri time-series revisited - again.” Neurolmage, vol. 2, pp. 173—181, 1995. [37] J. Locascio et al., “Time series analysis in the time domain and resampling methods for studies of functional magnetic resonance brain imaging,” Human Brain Mapping, vol. 5, pp. 168—193, 1998. [38] K. Friston et al., “To smooth or not to smooth? bias and efficiency in fmri time-series analysis,” Neurolmage, vol. 12, pp. 196—208, 2000. [39] E. Bulimore and etc., “Colored noise and computational inference in neurophysiological (fmri) time series analysis: resampling methods in time and wavelet domains,” Human Brain Map ping, vol. 12, pp. 61—78, 2001. [40] M. McKeown and J. Wang, “Ica denoising for event-related fmri studies,” 27th Annual Inter national Conference of the Engineering in Medicine and Biology Society, no. 2005. [41] M. Borga et al., “A canonical correlation approach to exploratory data analysis in fmri.” [42] C. Rorden and M. Brett, “Stereotaxic display of brain lesions,” Behav Neurol, vol. 2, PP. 191—200, 2000. [43] R. Patel, F. Bowman, and J. Ruling, “Determining hierarchical functional networks from auditory stimuli fmri,” Human Brain Mapping, pp. 462—470, 2007. [44] C. Hojen-Sorensen and L. Hansen, “Baysian modeling of fmri time series,” Neuin Advances in Neural Information Processing Systems, vol. 6, pp. 754—760, 2000. [45] S. Faisan et al., “Unsupervised learning and mapping of brain fmri signals based on hidden semi-markov event sequence models,” MICCA 12003, LNCS2879, pp. 75—82, 2003. [46] I. Shmulevich and W. Zhang, “Binary analysis and optimization-based normalization of gene expression data,” Bioinformatics, vol. 18, pp. 555—565, 2002. [47] R. Duan et al., “Activation detection on fmri time series using hidden niarkov model,” in The 2nd International IEEE EMBS Conference on Neural Engineering, 2005. [48] X. Zhou, X. Wang, and E. Dougherty, “Construction of genomic networks using mutual in formation clustering and reversible jump markov-chain-monte-carlo predictor design,” Signal Processing, vol. 4, pp. 745—761, 2003. [49] N. Friedman, K. Murphy, and S. Russell, “Learning the structure of dynamic probabilistic networks,” Proceedings of the Fourteenth Conference on Uncertainty in Artificial Intelligence, vol. 4, pp. 139—147, 1998. [50] E. Dougherty, S. Kim, and Y. Chen, “Coefficient of determination in nonlinear signal pro cessing,” Signal Process, vol. 20, pp. 2219—2235, 2000. [51] S. Kim et al., “A general nonlinear framework for the analysis of gene interaction via multi variate expression arrays,” J.Biomed. Opt, vol. 5, pp. 411—424, 2000. 57 [52] H. Lindman, “Analysis of variance in complex experimental designs,” San Francisco: W. H. Freeman Co., 2000. [53] J. Talairach and P. Tournoux, “Co-planar stereotaxic atlas of the human brain,” Thieme Medical Publisher, 1998. [54] F. Horak and M. Anderson, “Influence of globus pallidus on arm movements in monkeys. ii. effects of stimulation,” Journal of Neurophysiology, vol. 52, pp. 305—322, 1984. [55] K. Krakauer et al., “Differential cortical and subcortical activations in learning rotations and gains for reaching: a pet study,” Journal of Neurophysiology, vol. 91, pp. 924—933, 2004. [56] D. Vaillancourt et al., “Subthalamic nucleus and internal globus pallidus scale with the rate of change of force production in humans,” Neuroimage, vol. 23, pp. 175—186, 2004. [57] A. McIntosh, F. Gonzalez-Lima, and etc., “Structural modeling of functional visual pathways mapped with 2-deoxyglucose: effects of patterned light and footshock,” Brain Res, vol. 578, pp. 75—86, 1992. 58 Chapter 3 Dynamic Analysis of Probabilistic Boolean Network for fMRI Study in Parkinson’s Disease 2 3.1 Introduction Recent brain study has shown that the functionality of human brain is maintained by links between brain ROTs through efficient transcription regulations as well as sufficient formation and action of dopamine. As a result, the importance of single ROT perspectives has been greatly limited in obtaining insightful understanding of brain functioning. Researchers have been increasingly interested in systematic networks suitable for complex, large scale, and dynamical computation in an effort to intelligently extract and comprehensively understand massive information embedded in brain activity data from functional magnetic resonance imaging (fMRI). Several network models have been proposed in the literature to infer brain connectiv ity in fMRI, such as structural equation modeling (SEM) [58], multivariate autoregressive models (MAR) [59], dynamic causal modeling (DCM) [60], and dynamic Bayesian network (DBN) [61] , with each having its own advantages and disadvantages. BN approach does not require a prior structure based on anatomical connections, making it suitable for patholog 2A version of this chapter has been accepted for publication. Z. Ma, Z. J. Wang, “Dynamic Analysis of Probabilistic Boolean Network for fIVIRI Study in Parkinson’s Disease” in Proc. 30th IEEE Annual International Conference of the IEEE Engineering in Medicine and Biology Society. (EMBCSO), Vancouver, 2008. 59 ical conditions like PD, where functional connectivity may not be known before hand. DBN’s recent success for learning brain connectivity [61] has inspired us to explore the Probabilistic Boolean Networks (PBN) model [62], whose ancestor - Boolean Network - can be considered as a special case of DBN. Instead of focusing on quantitative biological details, PBN uses fundamental logical principles to regulate the studied system in the probabilistic context of markov chain, making it structurally simple yet dynamically complex [63]. Our previous work [64] has demonstrated the usefulness of PBN for brain connectivity modeffing and the effectiveness of L-dopa treatment for PD. However, another main goal of modelling PBN for brain ROTs is to identify potential therapy target through the asymptotic analysis, which could be facilitated by using the numerous theories and tools developed for markov chain. Meanwhile, we are also interested in suggesting way of changing the current trajectory of PD network by intentionally intervening identified target so that the system will be more likely to reach some desired states observed in normal subjects. To minimize the system disturbance brought by intervention, we plan to focus on transient intervention whose effects could be reversed by the network itself. The paper is organized as follows. In section.3.2,we first give a brief introduction of PBN framework in subsection 3.2.1, then followed by a detailed discussion of system dynamics during which we introduce the notion of random perturbation in subsection 3.2.2. Subsection 3.2.3 discusses a strategy for identifying the optimal intervention target based on the first- passage time theory. A case study using fMRI data from normal and PD subjects performing a squeezing motor task is discussed in Section 3.3. 3.2 Methods 3.2.1 Theoretical Background of PBN As an extension of well-known Boolean networks (BNs), a PBN models a system as a multivariate stochastic process under the binary domain, with the interactions among 60 system nodes being governed by logical rules. A BN consists of a set of node variables (xi(t), ..., x(t)), and a set of binary functions (f1t), ..., f(t)). Assuming that each node follows a switch-like behavior, moving between on and off states in the temporal dimension, and that nodes are not independent with each other, BN models the system as a highly- regulated binary network, where the activity of each node x(t) is determined by nodes known as parent nodes through one of those binary functions. Values of nodes are syn chronously updated by corresponding functions based on the current state of parent nodes which the functions take as input variables, as demonstrated by a triple-node BN example represented by a truth table in Table.2.1. Note that the dynamic of a BN is deterministic. The system always follows a static transition mechanism regulated by those binary functions, and ultimately settles to a so- called attractor state, from which the system cannot move. To avoid the potential risk caused by the deterministic rigidity of BN, Shmulevich proposed adding a probabilistic setting to BNs [621. The basic idea of a PBN is to accommodate more than one function for each node, allowing the model to have more flexibility. Typically, for each node x, there is a set of binary functions f known as predictors, among which one is chosen to update the target node’s state according to its probability p at each transition: F = {J} ,j = 1 ,m,P’r(x(t+ 1) = j(t)) =p, (3.1) where m is the number of predictors for x. As shown in Table.2.2, we could extend the previous BN example in the context of PBN by incorporating more functions with assigned probabilities. 3.2.2 Dynamics of PBN The dynamic behavior of PBN is represented by a Markov chain whose transition matrix is determined by the predictors as well as their probabilities. Of particular interest is the 61 asymptotic behavior of the PBN and if there exists a steady-state distribution the system will lead to regardless of its initial distribution. We study the asymptotic distribution of a n-state PBN in the strict context of a homogeneous Markov chain by first defining a stationary distribution ir, so that the possibility for the system to be in state i at any time r is always equal to ‘ir, if the initial distribution is ir. = (3.2) where p being the probability that the system transits from state i to state j within r times. Whether a PBN system has a single or multiple stationary distributions depends on its reducibility. If a PBN is irreducible, implying that it has a communicating state space where any state can possibly be reached by other states, then it is assigned a unique stationary distribution. However, a unique stationary distribution does not guarantee the existence of a steady-state distribution. A PBN has a steady-state distribution if and only if it is irreducible and aperiodic (ergodic). In other words, the stationary distribution is also the steady-state distribution if it satisfies: 1imp = rr. (3.3) Since the brain is always exposed to external stimuli, it is “biologically meaningful” to incorporate a random perturbation factor into the system. As proposed in [63j, this factor can be modeled as an iid vector -y e(l,O) where any element of ‘y could be one with probability p, implying that the corresponding ROl gets perturbed. A mathematical way of representing random perturbation is as follows: (x -y with probability 1-(1-p), (3.4) XI = 1. fk(x) with probability (1-p), (3.5) 62 where x is the original n-ROI system state vector, while x’ is the next-time state vector. If no perturbation happens, x’ will be determined by a PBN realization fk(x). Otherwise, x’ will be the result of modulo 2 addition of x and -y. The introduction of random pertur bation in PBN system makes the underlying Markov chain ergodic [63]. An ergodic Markov chain guarantees the existence of steady-state distribution so that the system’s asymptotic behavior can be learned regardless of the initial starting state. 3.2.3 PBN Intervention The system asymptotic study may help researchers to identify some states which are more easily reachable than others. A study regarding the sensitivity of PBN steady-state distri bution to perturbation indicates that stability of those states is more robust against random perturbation. Referred as phenotypes, those states may correspond to critical feature func tioning of the studied subjects. Therefore, significant differences between phenotypes of the PD subjects and the normal subject may shed light on the direction of potential therapy. This motivates us to develop an intervention strategy for letting the PD network to become more consistent with the normal network. In other words, we aim at intervening a certain target ROT so that the PD network will be more likely, more quickly to reach normal phe notypes. There are two types of criteria for locating the optimal target, and they could be used in conjunction [63]. The first one focuses on maximizing the probability (more likely) of reaching the desired states (normal phenotypes) while the second one minimizes the time (more quickly) spent for the network to reach the desired states. Though having different interpretation regarding the optimal intervention performance, both criteria use the same underlying algorithm - first passage time theory. We denote Fr(, y) as the probability that the system will visit state y for the first time at time r, starting from state x. It can be calculated in an iterated manner base on the system transition matrix T: Fr(,y) = T(x,z)Fr_i(z,y). (3.6) z(1,O)’—(y) 63 Criterion I measures the probability, Hi(x, y), that the system will reach y from x before time J?o: H,(x,y) = Fr(,y), (3.7) while Criterion II evaluates the mean first-passage time from x to y, M(x, y), as: M(x,y) = rFr(x,y). (3.8) Assume we have identified a PD phenotype x and a normal phenotype y, and we only consider single ROT perturbation for simplicity. Suppose the ith ROT is intervened by flipping its value while others remain the same, then we will get a new vector denoted as Given a sufficient large R0, Criterion I chooses the ROT resulting in the maximum HR0 as the optimal intervention target: = arg max H(x,y). (3.9) Criterion TI chooses the ROT with the minimum mean first—passage time as time target: = arg rnax M(x,y). (3.10) It is noted that phenotypes often form irreducible sub chains under PBN modeling when there is no random perturbation. The introduction of a small perturbation will create “bridges” between those subchain and the rest of the network. Therefore, we could extend our intervention strategy by slightly modifying the first-passage probability Fr(X, y) to deal with situations where the system needs to be transferred from a set of states X to a desired set of states Y. By employing a similar strategy, we can also identify intervention target(s) to help the system to stay away from certain undesired states. 64 3.3 fMRI Case Study 3.3.1 fMRI Data Collection By using a pressure-responsive bulb which was electronically connected to a computer, subjects including ten normal and ten PD patients were asked to squeeze the bulb by using left hand to control the height of a vertical bar to match a target bar moving up and down in a sinusoidal fashion. They lay down and squeezed a rubber bulb at two frequencies (0.00Hz and 0.75Hz) for 30-s blocks as instructed in a pseudo-random order. The patients took the experiment twice, once before medication and the other after medication. fMRI data of their brain activities during performing the task was collected with a Philips Achieva 3.0 T scanner. The following eighteen brain regions were selected as the ROTs in the study, both the left and right: primary motor cortex (Ml), supplementary motor cortex (SMA), lateral cerebellar hemisphere (CER), putamen (PUT), caudate (CAU), thalamus (THA), prefrontal cortex (PFC), anterior cingulate cortex (ACC), and globus pallidus (GLP). 3.3.2 Preprocessing of fMRI data The raw fMRI data first goes through motion correction and slice timing correction, and ROTs under study is drawn manually on each structural scan aligned with the functional data. The fiVIRI data further goes through denoising, voxel selection and binarization before being feed into PBN. Since those steps have been covered in great detail in [64], we only provide brief descriptions here. fMRI Denoising: Most current denoising methods focus on building a reasonable model for serially correlated fMRI data including recent approaches such as Independent Com ponent Analysis (ICA) and Canonical Correlation Analysis (CCA) [65]. Though we have included CCA as part of a complete denoising process in our study, for simplicity, here we simply remove the low-frequency drift from raw data through linear detrending as denoising. Voxel Selection: Choosing appropriate voxels to represent the overall behavior of each 65 ROT will significantly contribute to the success of our PBN analysis. There exists no uni versal standard on how to choose a subset for optimal performance. Our approach for voxel selection combines several previous methods: we first select activated voxels in each ROl based on t-value distribution. Then we apply spatial smoothing upon activated voxels by replacing its time course with the arithmetic mean of voxeLs contained in a cubic certering at that voxel. Finally we obtain ROT-wise time series by averaging activated voxels in each ROT. Binarization: Facing a challenge to develop a feasible binarization method of quantizing real-valued fMRT data into active (1) or inactive (0), we adopt the approach proposed in [66], which focuses on capturing the first-order statistical trend in the fMRT data because of its robustness against parameter change. It is reasonable to treat fMRT measurements, known as the blood oxygen level-depend (BOLD) response, like a hidden markov model (HMM) process because they directly reflects the tested event. 3.3.3 Results and Discussions We introduce a random perturbation into the PBN transition probability calculation to ensure that the system has a unique steady-state distribution. We set the perturbation probability to be 0.01 suggested in [63] to minimize its influence on asymptotic behavior of the system. We employ the “virtual-typical-subject” (VTS) group-analysis approach by pooling all subjects within the same group together and treating them as a single subject. Among 18 ROTs under study, both left and right CER, Ml, and SMA are of particular interest in PD study based on clinical knowledge. Therefore, though we learned the PBNs for 18 ROTs, instead of computing a 2’8-by- 218 transition matrix for the entire network, which could be extremely time-consuming, we focus on the above 6 ROTs and study the corresponding26-by-2 transition matrix by statistically incorporating effects of other ROTs. The steady-state distribution results, shown in Fig.3.1, indicate that the normal network is more likely to stay at state ‘111111’ (indexed by state 64) while the PD network is 66 more likely to be at state ‘000000’. Such a difference has been greatly alleviated for PD patients after medication, which reinforces the usefulness of current L-dopa in PD treatment. Furthermore, by computing the probability of each ROl being active, we observe that all interested 6 ROIs are more likely to be active in normal network than in PD network, as supported by recent clinical research.Once again, it is observed that generally ROT activation in PD subjects after medication starts approaching that of the age-matched controls, with an exception observed for right-SMA region. Since our fMRI data are left-hand squeezing experiment performed at high frequency (0.75Hz), we could further simplify our asymptotic analysis by focusing on left-CER, right Ml, and right-SMA, all of which are directly related to left-hand functioning. By comparing 3-ROT state diagrams for normal, and PD subjects shown in Figure 3.2, respectively, we observe that ‘000’ is more reachable in PD network while ‘111’ is more reachable in normal network, implying that we can consider ‘111’ as normal phenotype and ‘000’ as PD pheno type. We then look for the optimal intervention target, which provides the best first-passage time performance in reaching ‘111’ when starting from ‘000’. Figure 3.3 depicts how the first-passage time probability ((), y) converges as R0 increases, where y is the desired state ‘111’, and is obtained by intervening a single entry of ‘000’: ‘OOl’,’OlO’,’lOO’. It is shown that intervening right-SMA provides better results in reaching the desired state according to Criterion T. We also employ criterion II to estimate the mean first-passage time by setting a large “cut-off’ time value m0 so that m Mmo(x,y) = rFr(x,y), (3.11) as indicated in Table 3.1. The results show that intervening right-SMA yields the minimal mean first-passage time. Therefore, as suggested by both criteria, right-SMA could be a potential ROT candidate for future PD treatment. 67 3.4 Conclusion We have presented a Probabilistic Boolean Network (PBN)-based framework for fMRI data analysis in PD, aiming at identifying phenotypes as well as locating optimal intervention ROT target. The proposed approach opens up a new gate for PD treatment. Future work will focus on improving the group analysis in fMRI using PBNs. We will extend our approach to deal with the multiple ROI intervention case. In addition, Leave-One-Out cross validation will be used to test the robustness and consistency of the proposed framework. 68 0.35 0.3 0.25 0 (‘3 -o 0 a a) (‘3 tc 0.15 0 (‘3 a) 0.1 0.05 (a) 6-ROl steady-state distributions for normal, and PD subjects 1- 0.9 a) C) (‘3 0) C a) 0 0 4- 0 I 0.B 0.7 0.3 0.2 right-C ER left-CER left-Mi right-Mi * Normal • PD off-med • PD on-med right-SrvlA left-S MA (b) The probability of being active for interested ROIs Figure 3.1: 6-ROl case study 0.2 6-ROl state I I I __________________________ * * 4 * I . I I 69 001 .35 001 .33 0.14 000 . 010 .62 000 010 .36 .17 0.33 0.27 0 0.2 0. 011 .38 100 .44 011 .38 100 .42 0.26 0.28 0.26 o.i .21 0.15 101 •.16 111 .56 101 04 / 111 .54 0.19 110 .43 110 .32 (a)Normal Subject (b)PD Subject Figure 3.2: 3-ROTs state diagrams after modification to eliminate small transition probabil ities created by perturbation 70 0.9 I I I I I I0.8• I I I ê—right- Ml I IeftCER ,Ar AA,.0.7 I A’ c 06 A A’2 V’ 0)0.5 1/ — A 0.4 (15 / ‘ Aci) fr A(15 / ‘p /:0.2 0.1 fA I 0 I I I I I 0 2 4 6 8 10 12 14 16 18 20 RD Figure 3.3: Intervention performance comparison based on Criterion I, HR0 y) for R0 = 1...20, with y =11111. Table 3.1: Intervention performance comparison based on Criterion II, Mmo(3*, y) where 0=20. Intervention candidates ‘001’ I ‘010’ I ‘100’ 1 Mean first-passage time 5.2741 6.3232 6.9948 71 3.5 Bibliography [58] M. Goncalves and D. Hull, “Connectivity analysis with structural equation modeling: An example of the effects of voxel selection,” Neurolmage, vol. 20, pp. 1455—1467, 2003. [59] S. Kim et al., “A general nonlinear framework for the analysis of gene interaction via multi- variate expression arrays,” J.Biomed.Opt, vol. 5, pp. 411—424, 2000. [60) K. Friston, L. Harrison, and W. Penny, “Dynamic causal modeling,” Neurolmage, vol. 19, pp. 1273—1302, 2003. [61] J. Li, Z. Wang, and M. McKeown, “Dynamic beyesian networks (dbns) demonstrate impaired brain connectivity during performance of simultaneous movements in parkinson’s disease,” 3rd IEEE International Symposium on Biomedical Imaging: Nano to Macro, 2006. [62] I. Shmulevich et al., “Probabilistic boolean networks: A rule-based uncertainty model for gene regulatory networks,” Bioinformatics, vol. 18, pp. 261—274, 2002. [63] I. Shmulevich, E. Dougherty, and W. Zhang, “A tutorial on hidden Markov models and selected applications inspeech recognition,” Bioinformatics, vol. 18, pp. 1319—1331, 2002. [64] Z. Ma, Z. Wang, and M. McKeown, “Probabilistic booleari network for inferring brain con nectivity using fmri data.” [65] M. Borga et al., “A canonical correlation approach to exploratory data analysis in fmri.” [66] R. Duan et al., “Activation detection on fmri time series using hidden markov model,” in The 2nd International IEEE EMBS Conference on Neural Engineering, 2005. 72 Chapter 4 Conclusion 4.1 Conclusions In this research, we proposed using Probabilistic Boolean Networks (PBNs) to study brain connectivity based on fIVIRI data from PD patients and normal subjects performing motor tasks. In Chapter 2, we exploited the effectiveness of PBNs to identify significant connec tivity deviation possibly caused by PD. To fit the raw fMRI data into PBNs’ framework requires several pre-processing steps that have always been great challenges. In our exper iments, we combined some traditional pre-processing approaches that have proven feasible in previous study, with certain experiment-specific modification. The subsequent group analysis of learned networks takes into account both group factor and frequency factor to explore their influence on brain connectivity. The results indicate that PD subjects have significant connectivity difference, compared with normal subjects, and those difference is either consistent with clinical research or provides promising direction for future study. We also observed that the popular L-dopa medication does ease PD symptoms by normalizing some deviated connectivity. Furthermore, the group analysis results for frequency factor suggest that frequency can have profound effects on brain connectivity. In Chapter 3, we expanded our PBN modeling to study asymptotics of brain connec tivity under the context of markov chain. We adopted virtual-typical-subject approach by combining data within each group into a single subject to increase SNR as well as to have enough sample size for learning parameters. The simulation results based learned parame ters demonstrate higher level of task-induced activity on particular set of ROTs in normal 73 subjects than in PD patients. In addition, learned system transition networks allow us to identify certain states as phenotypes characterizing each group, upon which we employed First-time Passage Theory to locate optimal intervention target in normalizing PD networks. 4.2 Future Work As presented in this paper, PBNs privide an efficient way of studying brain connectivity at low computation cost while accommodating the dynamic properties of fIVIRI data. There are nevertheless some shortcomings that need further improvement for better performance. Firstly, more rigorous denoising and binarization methods must be derived before the pro posed analysis technique can be widely used in other applications. Secondly, the accuracy and reliability of PBNs are largely dependent upon both model selection and inference schemes. Therefore, more emphasis should be put upon how to optimize our model se lection, and more efficient inference algorithm should he explored to replace the current time-consuming exhausive-search learning method. Thirdly, influence connections should be expanded to include both direct and indirect cases in group analysis to reflect the overall system effects. Finally, Leave-one-out corss validation will be used to justify the robustness and consistency of the proposed network. 74
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Probabilistic Boolean network modeling for fMRI study...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Probabilistic Boolean network modeling for fMRI study in Parkinson's disease Ma, Zheng 2008
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Probabilistic Boolean network modeling for fMRI study in Parkinson's disease |
Creator |
Ma, Zheng |
Publisher | University of British Columbia |
Date Issued | 2008 |
Description | Recent research has suggested disrupted interactions between brain regions may contribute to some of the symptoms of motor disorders such as Parkinson’s Disease (PD). It is therefore important to develop models for inferring brain functional connectivity from data obtained through non-invasive imaging technologies, such as functional magnetic resonance imaging (fMRI). The complexity of brain activities as well as the dynamic nature of motor disorders require such models to be able to perform complex, large-scale, and dynamic system computation. Traditional models proposed in the literature such as structural equation modeling (SEM), multivariate autoregressive models (MAR), dynamic causal modeling (DCM), and dynamic Bayesian networks (DBNs) have all been suggested as suitable for fMRI data analysis. However, they suffer from their own disadvantages such as high computational cost (e.g. DBNs), inability to deal with non-linear case (e.g. MAR), large sample size requirement (e.g. SEM), et., al. In this research, we propose applying Probabilistic Boolean Network (PBN) for modeling brain connectivity due to its solid stochastic properties, computational simplicity, robustness to uncertainty, and capability to deal with small-size data, typical for fIVIRI data sets. Applying the proposed PBN framework to real fMRI data recorded from PD subjects enables us to identify statistically significant abnormality in PD connectivity by comparing it with normal subjects. The PBN results also suggest a mechanism of evaluating the effectiveness of L-dopa, the principal treatment for PD. In addition to PBNs’ promising application in inferring brain connectivity, PBN modeling for brain ROTs also enables researchers to study dynamic activities of the system under stochastic conditions, gaining essential information regarding asymptotic behaviors of ROTs for potential therapeutic intervention in PD. The results indicate significant difference in feature states between PD patients and normal subjects. Hypothesizing the observed feature states for normal subject as the desired functional states, we further explore possible methods to manipulate the dynamic network behavior of PD patients in the favor of the desired states from the view of random perturbation as well as intervention. Results identified a target ROT with the best intervention performance, and that ROl is a potential candidate for therapeutic exercise. |
Extent | 1677549 bytes |
Subject |
Brain connectivity Group analysis |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-02-04 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0066945 |
URI | http://hdl.handle.net/2429/4172 |
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 | 2008-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2008_fall_ma_zheng.pdf [ 1.6MB ]
- Metadata
- JSON: 24-1.0066945.json
- JSON-LD: 24-1.0066945-ld.json
- RDF/XML (Pretty): 24-1.0066945-rdf.xml
- RDF/JSON: 24-1.0066945-rdf.json
- Turtle: 24-1.0066945-turtle.txt
- N-Triples: 24-1.0066945-rdf-ntriples.txt
- Original Record: 24-1.0066945-source.json
- Full Text
- 24-1.0066945-fulltext.txt
- Citation
- 24-1.0066945.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}]}"
data-media="{[{embed.selectedMedia}]}"
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-0066945/manifest