Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Probabilistic Boolean network modeling for fMRI study in Parkinson's disease Ma, Zheng 2008

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

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

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 List of Tables  .  iv  vi  .  List of Figures  .  vii  .  Acknowledgements  ix  Introduction  1  1  2  .  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  1.3.2  Multivariate Autoregressive Model  1.3.3  Dynamic Causal Modeling  1.3.4  Dynamic Bayesian Networks  11  1.3.5  Group-Analysis Methods  11  .  .  .  .  5 7 8  1.4  Research Objectives and Methodology  13  1.5  Bibliography  17  PBN Connectivity Analysis  18  2.1  18  Introduction  iv  2.2  2.3  3  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  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  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  4  Methods  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  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  2.5  Connections showing significantly difference due to the frequency factor, Left  47  Hand  52  2.6  Performance of Frequency on subjects, Left Hand  53  3.1  Intervention performance comparison based on Criterion II, Mmo (x( ), y) where 1 =  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 x , node x 1 , node x 2 , respectively 3  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  3.1  Method under Frequency factor  55  6-ROT case study  69  vii  3.2  3-ROTs state diagrams after modification to eliminate small transition prob abilities created by perturbation  3.3  70  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 bloodoxygen-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,  /AC  =  x+(vw),  7AD  =  (vy) + (xz) + (vwz),  YBC  =  w+(vx),  7BD  =  Y + (wz) + (vxz),  7CD  =  z + (wy) + (xvy),  (1.1)  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,  B  =  vA+’çbB,  C D  (1.2)  , 0 xA+wB+b =  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  =  Y A 1 (i) + E,  where Y,. is the nth measurement of a d-dimensional time series, while Y_ (i previous measurements back to (n  —  (1.3)  =  1... ,p) are  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  where Z  =  =  F(Z, u, 0),  (1.4)  [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,  A  =  8F/0Z,  3 B  =  C  =  (1.5)  aF/au,  where Jacobian matrix A represents the first-order connectivity between ROTs without perturbations, corresponding to the second set of parameters. B 3 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  —  3 Icjs  — —  8f / 3 8t  =  Si,  r8v/3t  =  fi  rdq/8t  =  1 E 3 f , (.f Pi)/Pi  —  1),  (1.6)  v 1/cr  —  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 X , 1 denoting its status variable at each time point. The process can he fully specified by the conditional distribution P(X/X_ ). Since P(X/X_ 1 ) and the initial distribution P(Z 1 ) 1 does not change over time, the DBN grach could be defined as two time-slices: , 2 a 1 X as ndX 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 virtualtypical-subject (VTS) approach, an individual-structure (IS) approach, and a commonstructure (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 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  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.  15  r I I I I I I I I  —  — — — — — — — — — — — —  DBN  I  ®GG  I  I  I I  I  I  I  I I II I b  oi 1 Yt+  .  U  U  •  I II —  — — — —  N  — — — —  — —  I  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 X 1 and U 2 are associated with transition distribution. 1 to X  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 Parkinsons Disease 1 .  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 A version of this chapter has been published/submitted for publication. Z. Ma, Z. J. Wang, and M. 1 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 individualstructure (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. =  iPt-i +  Pt  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 components s(t), (i  =  1  =  1.  .  .  n) is a linear combination of  n). xi(t)  =  aiisi(t) +... + ais(t),  x ( 2 t)  =  a2isi(t) +  x(t)  =  aisi(t) +... + as(t).  ..  s(t), 2 +a  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—_W — 7 KI—WH,A=I(R>c*a),  (2.4)  where a 2 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  1 7t  =  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, 1 Poff I  ,  , =  tEpoff  V  Poff  (O  (2.5)  —  ff 0 tEp  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 (f (t), 1  ...,  ...,  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 x 1 (t), and x 2 (t) through Boolean function:  xi(t+ 1)  =  xi(t)flx ( 2 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(x (t + 1) 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 x 1 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 longrun (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-,  time r is always equal to  so that the possibility for the system to be in state i at time any 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 obtain the asymptotic distribution  ir  =  =  [1/8,  ...  1/8]) in our example, we could  [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 one of its Boolean functions  =  1,2,  ...,  C), the COD of  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, estimate of x, by  f) is the error of the  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: (2.11)  p= 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  f and corresponding probabilities p (j  =  upon x, i.e Ik(x), is defined as the sum of the influence of  Xk  x with a function set F including all predictors 1,...,m).The influence of upon each predictor  Xk  J: 1k(X)  =  (2.12)  1k(f)P,  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 x . The influence matrix G enables us to obtain a complete 3 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 (M 0  —  ‘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 2.3.1  fMRI Case Study 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, agematched 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 Lif HA  —*  —*  LM1, LPUT  RCER, LifHA  RSMA, and LACC  —*  —*  —*  LSMA, LCAU  LPFC, &CER  -+  —*  Rif HA, LCAU  RCER, RM1  -*  —*  RGLP,  RM1, RSMA  —  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  T-rnapping  Detrend  F—-  -  Comparison of t-value with  Comparison of t-value with threshold  Lshold  Voxel  Voxel Adjustment  Adjustment ROl Pre-processing ForROll  ROl Pre-processing ForROll8  •  Voxel Selection  Voxel Selection ROl Mean Determination  ROl Mean Determination  HMM Binarization  HMM 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 3 2 1 x  fi  000 001  0 0 0 0 1 1  010 011 100  f 2 0 1 0 1 1 1  f 3 0 0 1 1 1 1  101 110 1 1 0 111 1 1 0 This BN system has three nodes (x , x 1 , x 2 ) with each having its corresponding binary 3 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 ii  :1 J2  :2 J1  :2 J2  :3 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  3 2 1 X  0.6 0.4 0.5 0.5 p 1 This PBN expands previous BN by allowing x 1 and x 2 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  001 000  0.2  111 02  03  v.3 0.3 010  100  0.  0  o  05  0:  03 01  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 LPUT —* LSMA RPUT -* RPUT RPUT —* RM1 LCAU —* RifHA LCAU — RGLP LifHA - RCER LifHA — LPFC RCER -÷ RCER LLER —* RCAU LM1 — LCAU R_M1 —* R_M1  0.0156  —*  Arrows  —*  RSMA — RSMA RSMA —f LSMA RSMA —* RGLP L.ACC -÷ RCER denote influence connections with a time  0.0335 0.0011 0.0184 0.0299 0.0311 0.0417 0.0204 0.0308 0.01 0.0117 0.0338 0.0269 0.0218 0.0187 0.0022 lag. The critical value here is .05.  46  Table 2.4: Performance of Medicine on PD, Left Hand Connection  LPUT —* LPUT —* RPUT — RPUT —* LGAU —* L.CAU — LifHA LifHA  MNorm  L_Ml LSMA RPUT R_Ml RifHA RGLP  —  Mp f 0 f  0.0129 -0.0356 0.1876 0.0745 0.1132 0.0434  ff 0 M Mp 0 0.0863 -0.1251 -0.2091 -0.0277 0.0180 0.1066 —  RCER 0.0865 0.0049 -+ 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 ff denotes the difference of estimated means of influence coefficients between 0 MNorm — Mp Norm and Poff. The same sign shared by MNorm Mp ff and MP 0 0 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  3  1)4  -002  0 0.02 0.04 The column means of groups 1  (a) LifHA  0.t 000 0.1 and Pars significen try different  —*  0.12  0 1  LPFC  CliCk 00 the g.oup yOt te0ttotCt  2  a  -0.05  0 0.05 The OCIurCn means of gropsnana 2  (b) LACC  —*  0.1  0.15 goitcantly d,ffe,ent  0.2  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.  48  I  I.  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  to test  Click o.. n.e go..p ye..  2  3  0.1  03 0.4 0,2 05 me column means of groups land S are significantly dIfferent  (a) RPUT  06  0/  RPUT  —  Click o.. Ofl group you want to t000  2  04  -0.02  0 0.02 0.04 0.06 0.12 0.00 0.1 The colurt,n means of groups I arId Sore sigrliocan try diffCrCnt  (b) RPUT  -+  0,14  0.  RM1  Click On One group you w000 to tCal  2  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.  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 LP(JT —* Rif HA  Arrows  —*  LPUT -÷ RACC RPUT —* RifHA RPUT —* RSMA RCAU — RCAU LifHA —* RSMA RCER — RCAU RCER —* RM1 LCER —* LM1 RM1 —* LSMA RSMA — LGLP LSMA —* R..M1 LPFC —* LifHA R.ACC - RGLP LACC —* LPFC LAGC — RACC denote influence connections with a time  P-value 0.0123 0.046 0.0328 0.0498 0.0088 0.0282 0.0355 0.0227 0.0385 0.0229 0.0361 0.005 0.0322 0.0344 0.0237 0.0347 lag. The critical value here is .05.  52  Table 2.6: Performance of Frequency on subjects, Left Hand Connection LPUT —* RifHA LPUT —* RACC RPUT — R,ff I-IA RPUT —* RSMA RCAU - RCAU Lif HA —* RSMA RCER —* RCAU &CER —* R_M1 L_CER —* LM1 RM1 — L_SMA RSMA —* LGLP L_SMA —* R_M1 LPFC —> LifHA RACC —* RGLP LACC  LPFC  MH MM 0.0282 -0.0604 0.0351 0.0618 -0.007 -0.0103 0.0146 -0.0789 -0.0253 -0.0421 0.0153 0.0711 -0.0585 0.0675 —  ML MM -0.0135 -0.0229 -0.0029 0.0173 0.0229 0.0610 0.0008 -0.0382 0.0194 -0.0046 -0.0083 -0.0567 -0.0639 0.0208 0.0165  -0.0068 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)  np  2  pa  3  ph  _1_  -0.04  J  .l  J  J  -0.02 0 0,02 0.04 0.06 0.08 01 The row means of groups land 3 are signffcantly different  (a) LTHA  I  I  I  I  I  -0.05  0  005  0.1  0.15  No groups  RSMA  -÷  op  II  pa  0.2  0.25  have means signffcanlly different frsm np  (b) L2”HA  op  —*  I  LM1  I  I  I  I  pa  ph  ph I  -Ut  -0.1  I  I  I  I  -1105 0 0.05 Ot 0.15 0.2 0.25 The means of groups op and pa are signthcaotly different  (c) RPFC  —  &ACC  I  -0.04 -0.02  I  I  I  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.  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? analysis,” Neurolmage, vol. 12, pp. 196—208, 2000.  bias and efficiency in fmri time-series  [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 A version of this chapter has been accepted for publication. Z. Ma, Z. J. Wang, “Dynamic Analysis 2 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 firstpassage 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 3.2.1  Methods 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 1 (f ( t),  ...,  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 highlyregulated 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 socalled 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 r  is always equal to  ‘ir,  ir,  so that the possibility for the system to be in state i at any time  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)  with probability (1-p),  (3.5)  XI =  1. fk(x)  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:  T(x,z)Fr_i(z,y).  Fr(,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)  =  (3.7)  Fr(,y),  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 R , Criterion I chooses the ROT resulting in the maximum HR 0 0 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 3.3.1  fMRI Case Study 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 8 2’ by-  218  transition matrix for the entire network,  which could be extremely time-consuming, we focus on the above 6 ROTs and study the corresponding 2 -by-2 transition matrix by statistically incorporating effects of other ROTs. 6 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 state ‘111’, and  ((), y) converges as R 0 increases, where y is the desired  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 m 0 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  0.2  a a) (‘3  tc 0.15 0 (‘3  a)  0.1 0.05  6-ROl state  (a) 6-ROl steady-state distributions for normal, and PD subjects 1-  I  I  0.9 a) C) (‘3 0) C  a)  I  *  *  Normal  •  PD off-med  •  PD on-med  0.B  *  4  0.7  0  *  0 40  I  I .  0.3 0.2 right-C ER  left-CER  left-Mi  I  I  right-Mi  right-SrvlA  left-S MA  (b) The probability of being active for interested ROIs Figure 3.1: 6-ROl case study  69  001  .35  001  .33 0.14  000  011  010  .  .17 0.27 0.33  .62  000 0 0.  .38  100  .44  011  0.26  •.16  .36  100  .42  0.2  .38 0.28  0.26 101  010  .21  o.i 111  0.15  .56  101  04  111  /  .54  0.19 110  .43  (a)Normal Subject  110  .32  (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  ê—right- Ml I IeftCER  0.8• I 0.7  I  I  I  I  I  ,Ar AA,. A’  c 06 2 0)0.5  1  V’  A  /  A  —  0.4 (15 ci)  fr  (15  /  A’  ‘  /  A A  ‘p  /:  0.2 fA  0.1 0  I 0  2  I  I  I  I  I  4  6  8  10 RD  12  14  16  18  Figure 3.3: Intervention performance comparison based on Criterion I, HR 0 = 1...20, with y =11111.  20  0 y) for R  Table 3.1: Intervention performance comparison based on Criterion II, Mmo(3*, y) where =20. 0 m Intervention candidates ‘001’ I ‘010’ I ‘100’ Mean first-passage time 5.2741 6.3232 6.9948  1  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 multivariate 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  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0066945/manifest

Comment

Related Items