%======================================
\chapter{Introduction}
\label{ch:intro}
%======================================
%======================================
\section{Overview}
%======================================
\label{s:overview}

Advances in neuroimaging technology such as functional magnetic resonance imaging (fMRI), positron emission tomography (PET), and
%single-photon emission computerized tomography (SPECT),
diffusion tensor imaging (DTI), have revolutionized our understanding of the human brain. The development of these neuroimaging techniques
has given researchers an opportunity for in-depth investigations into both the anatomy and the function of the human brain. %and investigate biological mechanisms of different neurological and psychiatric disorders.
Neuroimaging plays an increasingly important role in the diagnosis, therapeutic monitoring, and scientific investigation of different neurologic
and psychiatric disorders.
Neuroimaging studies yield massive data sets that pose challenges for statistical analyses due, in part, to the intricate brain neurophysiology.
The functional organization of the brain is based on information processing properties of specialized cells (neurons) and the complex and
malleable modes of their interconnectivity.
Our main objective was to study aspects of the complex spatial relationships present in functional neuroimaging data and to develop new
statistical methods that provide accurate inference in the presence of associated patterns and correlations.

The dissertation is organized as follows:
The remainder of Chapter 1 provides background info on the human brain, various neuroimaging technologies, and standard methods of statistical
analyses and also outlines our proposed research objectives.
Chapter 2 presents a novel statistical approach to model the complex spatio-temporal structure of neuroimaging data.
Chapter 3 presents methods to evaluate the level of connectivity within functionally defined neural networks, and
Chapter 4 presents a novel prediction framework for neuroimaging data.

% in Section 1.5 by proposing
%a set of questions of interest and subsequently describe our contributions to the statistical
%methodology of functional neuroimaging which serve to better answer these questions.
%I present the current state of the art before describing, in three chapters, novel methods
%to determine functional connectivity of the human brain and subsequently determine and
%visualize functional brain networks elicited during particular experimental conditions

%\section{Background Information}
%==================================================================
\section{An Introduction to Organization of the Human Brain}
\label{ss:brain}
%==================================================================
{\it Basic Neuroanatomy.}
The brain is the center of the human nervous system and constantly receives, interprets,
and directs sensory information and organizes motor and other responses throughout the body. The human brain has intricate structure and
function and has been the target of scientific inquiry for centuries.
The largest and uppermost section of the brain, %(anterior portion of the forebrain),
operating at the highest functional level, is called the cerebrum.
% The higher region of the brain that controls intelligence, learning, perception, and emotion.
The cerebrum consists of two hemispheres that communicate with each other by a thick bundle of fibers called the corpus callosum.
The right hemisphere controls the left side of the body and vice-versa.
Each hemisphere is covered by a thick layer (1.5mm to 5mm) of gray matter called the cerebral cortex.
\begin{figure}[!h]
\begin{center}
\begin{tabular}{cc}
  %\hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  \includegraphics[scale=0.4]{Figures/lobes.eps} & \hspace{-2em} \includegraphics[scale=0.475]{Figures/brain3.ps} \\  %axes.eps
  (a) & (b) \\
  %\hline
\end{tabular}
\end{center}
\caption[(a) Principal fissures and lobes of the cerebrum viewed laterally. (b) The axes of the central nervous system]
{\label{fig:brain_regions} (a)~Principal fissures and lobes of the cerebrum viewed laterally, from {\it Anatomy of the Human Body},
the 1918 edition, available on-line at {\tt www.bartleby.com/107} (Gray, 1918).
(b)~The axes of the central nervous system.}
\end{figure}
%
The cerebral cortex is heavily folded, resulting in ridges called gyri and grooves called sulci.
A cerebral hemisphere can be divided into four lobes: frontal, parietal, temporal, and occipital (Figure \ref{fig:brain_regions} (a)), each of which specializes in different functions.
%Subcortical structures include various ventricles and nuclei, the thalamus, hypothalamus, cerebellum, and the parts that make up the brainstem.
The cerebral cortex is connected to numerous subcortical structures such as the thalamus, the basal ganglia, hypothalamus, cerebellum, and the  brainstem.
%The main regions of the brain discussed later in the proposal are outlined in Figure
%\ref{brain_regions} (b).

Neuroscientists use several terms to describe orientation, some of which are illustrated in Figure \ref{fig:brain_regions} (b).
The following terms are used to describe direction in the brain:
%ventral- toward the belly (front), dorsal Toward the back
%Rostral Toward the nose
%Caudal Toward the tail
%Superior Toward the top of the head/body
%Lateral Away from the middle
%Medial Toward the middle
% Three directional planes exist in the brain: rostral/caudal, dorsal/ventral, and medial/lateral.
rostral/anterior means front or head end; caudal/posterior means hind or tail end; dorsal means back or top side; ventral means belly or bottom side; lateral means away from the midline, medial means toward the midline.
We typically view neuroimaging scans in three orthogonal planes referred to as axial (transverse or horizontal), coronal, and sagittal
(see Figure \ref{fig:brain_planes}).
%The coronal plane, horizontal plane and sagittal plane are shown in the figure A.1. The coronal plane is also called the frontal plane. It cuts the brain into the anterior and posterior parts. The horizontal (transverse or axial) plane cuts the brain into top (dorsal) and bottom (ventral) parts. The sagittal plane divides the right and left side of the brain into parts. The midsagittal plane would divide the right and left sides of the brain into two equal parts.
%Since the brain is three dimensional, any point on or inside the brain can be localized on three ``axes'' or ``planes'' - the $x$, $y$ and $z$ axes or planes.


%{\it Neuron.}
The functional cellular unit of the central nervous system is a nerve cell or {\it neuron}.
%The neuron consists of a cell body (soma), dendrites attached to the body, and an axon.
Brain tissue is divided into two types: gray matter and white matter.
% These names are derived from their appearance to the naked eye.
Grey matter contains neural cell bodies, in contrast to white matter, which mostly contains myelinated axons
(myelin is a collection of lipid fats and proteins that sheaths the axons). The color difference arises mainly from the white color of myelin.
%Groups of nerve cell bodies in other parts of the brain are called nuclei, or columns if they occur in long rows.
%Accumulations of nerve cell bodies outside the CNS (central nervous system) are called ganglia.''

\vspace{2ex}
{\it Brain Atlases and Coordinate Systems.}
When group studies using functional imaging data are performed, the individual brain images are usually transformed into a common coordinate space
to accomodate the variation.
The Talairach space and the Montreal Neurological Institute (MNI) space are the two most widely used atlas spaces in the neuroscience community.
%The Talairach coordinate system has become the standard reference for reporting the brain locations in scientific publication, even when the data have been spatially transformed into different brain templates (e.g., MNI space).
The Talairach coordinate system is based on a stereotaxic atlas of the human cerebral cortex published by Talairach and Tournoux
\citep{talairach:88}. It is related to anatomical landmarks % It is based on three brain structures (AC, PC and interhemispheric fissure)
%The yaxis is defined by the line connecting the most superior of AC and
%the most inferior of the posterior commissure (PC); the x-axis is
%defined by the line that passes through the AC point and
%orthogonal to the AC-PC line; whereas the z-axis is the line that
%passes through the interhemispheric fissure and the AC point.
and specifies locations relative to their distance from the midpoint of a brain white matter structure called the anterior commissure (AC) - a small but easy to spot region. Each location is described by three numbers which describe the distance in millimeters from the AC: $x$ is the left/right dimension, $y$ is the posterior/anterior dimension, and $z$ is the ventral/dorsal dimension. In this atlas the axial slices are referred to by their $z$ coordinate, coronal by their $y$ coordinate, and sagittal by their $x$ coordinate. The atlas is based on one brain, a brain of a 60-year-old French woman.
%It is important to stress that normalization strives to retain the unique features of each individual brain, and therefore Talairach coordinates are only approximate when comparing locations to other individuals.
%
\begin{figure}[!h]
\includegraphics[scale=0.4]{Figures/planes.eps}
\makebox[2ex]{}
\includegraphics[scale=0.45]{Figures/planes_ch2_labels.eps}
\caption[Planes of brain sections]{Planes of Brain Sections.}
\label{fig:brain_planes}
\end{figure}

The MNI defined a new standard brain by using a large series of MRI scans of healthy normal controls \citep{evans:93}. These atlases differ in shape and size. For example, the Talairach brain is considerably smaller than the average brain by up to 10 millimeters in each dimension. Software packages are available to approximately convert one coordinate system to the other.

%\subsection{A Brief History of Study of the Human Brain Function (?)}
%=========================================================
\section{Functional Neuroimaging}
\label{s:FNI}
%=========================================================
Neuroimaging is traditionally divided into {\it structural} and {\it functional} imaging. Structural imaging maps the brain anatomy and
includes computed tomography (CT) and MRI. Functional techniques seek to examine the physiological properties of the brain,
either at rest or during task-induced activation. %, and include PET, single photon emission computed tomography (SPECT), MRI spectroscopy and fMRI,
% and encephalographic brain mapping.
A variety of methods allow mapping of human brain functioning.
PET and fMRI measure localized changes in cerebral blood flow (referred to as activations) that accompany neuronal activity with
relatively high spatial resolution ($\leq$ 1mm), but with temporal resolution (2-20 sec) limited by the much slower rate of brain
regional blood flow and blood oxygenation.
In contrast, techniques such as electroencephalography (EEG) and magnetoencephalography (MEG) map
the underlying electrical activity of the brain cortex. These methods allow high temporal resolution of neural processes, but have poor
spatial resolution (over 1 cm).
%Temporal resolution determines our ability to separate brain events in time, while spatial resolution
%determines our ability to distinguish changes in an image across spatial locations.
% Spatial resolution: ability of the system to distinguish two small high-contrast objects separated by only a small distance.

While each modality is interesting in its own right, in this dissertation we focus and introduce new statistical methods applied to fMRI and PET data.
%Among various imaging modalities, functional fMRI and PET have enjoyed the most technical advancement.
In the following two sections, we give a brief description of the main principles on which these two neuroimaging techniques are based.
In this dissertation we concentrate on how they are applied to {\it human} brain neuroimaging.

%---------------------------------------------------------------
\subsection{Functional Magnetic Resonance Imaging (fMRI)}
\label{ss:fMRI}
%---------------------------------------------------------------
As opposed to PET, which requires the use of radioactive markers, % which limits the frequency of exposures for a single participant,
fMRI is a relatively safe and non-invasive technique for generating maps of and studying brain activity.
%It provides an indirect measurement of brain activity through changes in blood vasculature that accompany neural activity,
%either during a specific task or during spontaneous electrical activity. <-- add if there is space
%Unlike MRI data, which provide a static picture of the structure of the brain,
fMRI data consist of a 3-D sequence of individual magnetic resonance
(MR) images - a ``movie'' of a subject's brain activity. To understand the principles of fMRI, we need to understand how individual MR images are acquired.
Below, we present a brief overview of the principles of MRI.
% here, and refer the interested readers to \cite{higgins:96} and \cite{haacke:99} for more information. %the background information on the basic principles of the structural MRI and MR physics.

%----------------------------------------
\subsubsection{Basic Principles of MRI.}
``MRI is based on ability of organic tissue, placed in a magnetic field, to absorb energy from radio waves, and then re-emit it, in proportion to mobile hydrogen ion concentration'' \citep{higgins:96}.
The basis of this process is an interaction between radio waves and atomic nuclei called {\it nuclear magnetic resonance} (NMR).
%MR signal is a product of a phenomenon called {\it nuclear magnetic resonance} (NMR) - an interaction between radio waves and atomic nuclei (Bloch et. al., 1946, Purcell et. al., 1946). In the presence of a magnetic field, certain elements absorb energy from radiation at the radio frequency (RF) end of the electromagnetic spectrum. When the external source of radiation is turned off, energy is released again and can be picked up by a receiver, so the absorbed energy can be detected.
%Since much of the human brain is water, hydrogen nuclei are abundant, which makes them ideal candidates for brain imaging studies.
To obtain an MR image, a subject is placed in a field of a large electromagnet (generally from 1.5 to 4.0 Tesla).
The magnetic field aligns the magnetization of hydrogen atoms in the brain.
(Since much of the human brain is water, hydrogen nuclei are abundant, which makes them ideal candidates for brain imaging studies.)
A hydrogen nucleus whose spin is oriented
parallel to the applied magnetic field is said to be relaxed or in the low energy state, while a nucleus whose spin is
oriented against the magnetic field is is said to be in an excited or high energy state.
The scanner injects a pulse of radio frequency (RF) to excite the nuclei and raise them out of their low energy states.
Upon removal of the RF pulse, the hydrogen nuclei return to their original aligned position (i.e., lower energy states),
and in doing so, emit RF energy measured by the scanner. If the frequency of the input radio waves equals the natural resonance
frequency (the Larmor frequency) of the element in the magnetic field, the tissue will absorb the radiation and when radio
waves are turned off, it will emit a return signal of the same frequency, but lesser amplitude.
If the frequency does not match the natural resonance of the element, the element will neither resonate nor send a signal.
% So, if radio signal is injected at the Larmor frequency for hydrogen nuclei, the nuclei that absorbed the RF will dissipate
% this energy and return to their lower energy state when the RF pulse terminates.
The emission of this RF signal is the MRI signal that is recorded in fMRI studies. Devices (coils) that can generate additional
small local magnetic field gradients are used to slightly shift the Larmor frequency across brain areas and the scanner to recognize
the location of the signal by its frequency.

%{\it Signal weighting (T1, T2, PD) and sequences parameters : TR, TE.}
By adding RF or gradient pulses and carefully choosing their timing, it is possible to highlight different characteristics of the imaged tissue.
In general, MRI maps the water distribution in the brain. However, the useful contrast in MR images comes not only from spatial variation in
the density of water, but also from differences in nuclear magnetic properties known as relaxation.
They are characterized by distinct rates or {\it relaxation times}, used in MRI to distinguish between tissue types.
Three relaxation times are of primary interest in MRI: $T_1$, $T_2$ and $T_2^*$.
%$T_1$ relaxation is a measure of the time taken to realign with the external magnetic field. It can be thought of as the longitudinal component,
%since the protons are `tilting' back to their `upright' position.
% (will indicate how quickly the spinning nuclei will emit their absorbed RF into the surrounding tissue).
%a time constant for the return of the magnetization to its equilibrium position aligned along the static magnetic field of the scanner whenever it is disturbed,
%$T_1$ time refers to interval where 63\% of longitudinal magnetization is recovered.
%$T_2$ corresponds to the decay of magnetization perpendicular to the main magnetic field (transverse magnetization),
%which occurs as the protons gradually fall out of phase in their precession.
%$T_2$ time refers to interval where only 37\% of original transverse magnetization is present.
%$T_2$ and $T_2^*$ are time constants associated with loss of signal one the magnetization has been sampled.
$T_1$ effects measure recovery of longitudinal magnetization (parallel to the main magnetic field), while $T_2$ refers to decay of
transverse magnetization (perpendicular to the main magnetic field). $T_1$ time refers to interval where 63\% of longitudinal magnetization is recovered, and $T_2$ time refers to the interval where only 37\% of original
transverse magnetization is present.
%$T_1$ and $T_2$ vary for different tissues.
%In short, $T_1$ effects measure recovery of longitudinal magnetization and $T_2$ refers to decay of transverse magnetization.
%Relaxation times are strictly MR-based parameters that describe the re-growth of longitudinal magnetization ($T_1$) and
%loss of transverse magnetization ($T_2$) after a RF pulse flips magnetization out of alignment with the externally applied static magnetic field.
$T_2$ decay results from inherent, random variations in the precession of individual protons.
However, large-scale inhomogeneities in the magnetic field can also produce differences in precession frequency, resulting in signal reduction.
When $T_2$ dephasing (which occurs as the protons gradually fall out of phase in their precession) is due to to one or more localisable sources,
it is referred to as as $T_2^*$. So, $T_2$ has only one reason for dephasing and is intrinsic to tissue type (e.g. different $T_2$ of cerebral spinal
fluid and fat), while %Pure T2 dephasing is intrinsic to sample (e.g. different T2 of CSF and fat).
$T_2^*$ dephasing includes true $T_2$ as well as field inhomogeneity and tissue susceptibility.
(The magnetic field inhomogeneity can result from the microscopic magnetic field gradients resulting from
an increase in blood oxygenation.
%deoxyhemoglobin in the brain.
More detail will be given in the next section.)
%$T_2^*$ decay is always faster than $T_2$ decay.
%$T_2^*$ is usually much smaller than $T_2$.
%Due to these artifacts, Larmor frequency varies between locations.
%T2* leads to rapid loss of signal
%$T_2^*$ is the most relevant relaxation time for understanding contrasts in fMRI images.
fMRI is a $T_2^*$ image.

The raw data obtained from an MRI scanner are collected in the frequency domain. The inverse Fourier transform is then used to
transfer the data into image space, where data analysis is performed.

%-------------------------------------
\subsubsection {Principles of fMRI}
%--- define BOLD
fMRI relies on hemodynamic correlates of neural activity, not the activity itself.
%The fMRI technique that has been used most frequently in cognitive neuroscience research is the BOLD (Blood Oxygen Level Dependent) method.
It is most commonly performed using blood oxygenation level dependent (BOLD) contrast
\citep{ogawa:90}. % ,ogawa:92}.
The activation signal measured with BOLD fMRI is indirectly measuring changes in the concentration of paramagnetic deoxyhemoglobin
(hemoglobin without bound oxygen molecules) which arise from an increase in blood oxygen utilization in the vicinity
of neuronal firing.
%It is based on the increase in blood flow and blood oxygenation (collectively known as {\it hemodynamics}) to the
%local vasculature, which is known to be closely linked to neural activity. % (although the relationship is complex).
%Images are constructed based on blood-oxygenation-level-dependent (BOLD) contrast:
%a change in the signal strength of brain water protons produced by the paramagnetic effects of venous blood deoxy-Hb.

%The blood oxygen level dependent (BOLD) response is one of the most widely used functional imaging contrasts for brain
%activation studies.

%===== fMRI technology =======
fMRI is possible because of two basic phenomena. The first is the fact that regional changes in brain activity are associated with local changes in blood flow and blood oxygenation in the brain (collectively known as {\it hemodynamics}).
%When a certain area of the brain is active, the blood flow to that area increases over a wider volume and to a greater extend than needed.
%%it pulls more oxygenated blood to that area then is actually needed.
%This results in a surfeit of oxygenated blood associated with increases in neuronal activity.
%
The second phenomenon is the ability of paramagnetic agents to produce contrast in the MR signal.
(Paramagnetism is the ability of an otherwise nonmagnetic material to exhibit magnetic properties in the presence of magnetic field.)
%Deoxygenated blood has magnetic properties (because oxygen neutralizes the effect of iron in blood). The ratio of oxygenated blood to deoxygenated blood can be picked up as signal by the magnetic field of fMRI.
%We now give a brief description of the two phenomena.

{\it The physiological basis of BOLD fMRI.}
When a subject performs a behavioral task, for example a repetitive finger-tapping task, neuronal networks involving multiple brain regions are activated. These collections of neurons require large amounts of energy as adenosine triphosphate (ATP) to sustain their metabolic activity. Since the brain does not store its own energy, increased blood flow is required to deliver the necessary glucose and oxygen (which is bound to hemoglobin) to meet this metabolic demand.
There is a disproportionate increase in regional cerebral blood flow (CBF) relative to the cerebral oxygen metabolic rate
%(CMRO$_2$)
in response to neuronal activation (e.g., 30\% vs. 5\%).
%, i.e. the blood flow to the `active' area increases over a wider volume and to a greater extend than needed.
% CMRO2 increases less than CBF during activation (5\% versus 30\%)
This results in a decreased oxygen extraction fraction (more oxygen supply than consumption) and hence local `hyperoxygenation'
-- a surfeit of oxygenated blood. %This mismatch is the basis of BOLD fMRI.
%
%This leads to less deoxygenated hemoglobin in blood and increase in BOLD fMRI signal.
% Hence, there is relatively more oxygenated hemoglobin (oxy-Hb) than deoxygeneted Hb (deoxy-Hb) in more active areas.

{\it Biophysics of BOLD signal.} The reduced oxygen extraction leads to an increase in
the ratio of oxygenated hemoglobin (oxy-Hb) to deoxygenated Hb (deoxy-Hb) in the region of neural activation. Deoxy-Hb is slightly paramagnetic relative to brain tissue while oxy-Hb is dimagnetic. Because of its paramagnetic properties, deoxy-Hb has the effect of suppressing the MR signal, while oxy-Hb does not. Presence of deoxy-Hb causes microscopic magnetic field inhomogeneities (gradients) which leads to distractive interference
from signal within the tissue element called {\it voxel} (a unit of 3D volume).
This process tends to shorten the $T_2^*$ relaxation time.
As oxygen extraction falls with increased local blood flow in a region of greater neuronal activity, the $T_2^*$ becomes longer and the MRI signal
intensity increases relative to the baseline state.
%--- name: BOLD ------
%When the content of deoxy-Hb changes in the blood, the relaxation process of water protons is modified and these changes can be seen in MRI.
The image intensity that varies with deoxy-Hb content has been termed Blood Oxygenation Level Dependent (BOLD) and was suggested for potential use in functional study of the brain by \cite{ogawa:90}.

%In summary, during the execution of a mental task, the areas that are active are refreshed with oxygenated blood and their local magnetic susceptibility and the measured MRI signal are changing. Therefore, a series of MR images, if properly acquired, can can be used to study changes in blood oxygenation, which can be used to infer brain activity.

%Paramagnetic deoxy-Hb is washed out and signal increases.

%Activation of brain neurons results in localized change in blood flow and oxygenation. These hemodynamic changes are delayed by
%several seconds following the electrical activity of the brain.

{\it Basic parameters of image acquisition.}
There are many imaging design parameters (fixed at the beginning of the study) which control the spatial and temporal resolution of the data
and image quality in general. The most important intrinsic parameters, which affect the signal measured at each voxel, are
the repetition time (TR) and echo time (TE).  TR is the time, in milliseconds (ms), between successive applications of RF pulses to a particular
volume of tissue.  It is impossible to measure the signal immediately after the RF is applied, due to hardware limitations. The short
waiting time (also measured in ms's) during which the peak signal is obtained is called TE.
% EPI

The most common imaging sequence used in fMRI is the fast method of echo planar imaging (EPI), which allows collection of
whole brain data in a few seconds or less.
%Echo planar imaging (EPI) technique is a technique that provides a mean of acquiring MRI images in
%much faster than with the ordinary equipment.(A TR of less than 500 milliseconds is considered short, while a TE of more
%  than 15 milliseconds is considered long.)

%{\it Advantages and Limitations of fMRI.}
%The main advantages of fMRI over other brain imaging methods are its non-invasive nature
%and its good spatial resolution typically of order $3 \times 3 \times 5$ mm$^3$, corresponding to
%image dimension of order of $64 \times 64 \times 30$ (which can be samples in approximately 2s).
%
%The main disadvantage of the method is the rather poor temporal resolution (TR)
%of a few seconds (generally ranging from 0.5 s to 3 s), which
%is due mainly to the poor temporal resolution of the BOLD effect, not the MR technique itself.
%Namely, the temporal resolution of the fMRI is dependent on the hemodynamic response of the brain and
%how frequently the response is sampled. It is limited by the fact
%that the blood flow response typically lags the actual electrical signal by 1-2 s, peaks in 4-12 s, and does not
%track activity on a ms-by-ms basis.
%Since the blood flow response is influenced by activity levels over some time interval (a few
%hundred ms or more) it is less temporally specific than the neuronal activity with which it is
%associated.
%However, in contrast to PET, fMRI has better temporal and spatial resolution, although its biochemical sensitivity and selectivity
%is not as good as PET.
%
%% {\it Challenges in analysis of fMRI data.}
%Because BOLD contrast is fairly small, with activation inducing
%signal increase ranging from 5 to 10\%, a sensitive method is needed to detect
%the activation sites. Moreover, changes in fMRI brain images are not only due to
%activation, but also to other variables, such as pulsatile blood flow, imaging artifacts
%or patient motion, which all pose challenges for the analysis.

%Despite these problems, two objectives are currently pursued
%using fMRI: 1) Mapping the response to a stimulus onto an anatomical image, i.e.,
%identifying regions for which there is a specific correlation between a time-dependent
%sensorimotor task and a measured change in physiology; 2) Characterizing the
%transfer function of the brain.

%---- This is GOOD, but don't have space --------
%FMRI studies yield large amounts of noisy data sets that contain complex spatio-temporal correlations, driven by
%highly sophisticated neurophysiology and aspects of typical experimental
%designs. These complex dependence structures pose analytical challenges for statistical modeling.
%Statistics plays a crucial role in understanding the nature of the data and obtaining relevant
%results that can be used and interpreted by neuroscientists.

%`A major factor in the success in of activation studies with BOLD fMRI should be attributed to the development of dedicated statistical data processing (Friston et al 1995), which has provided easy access to the methodology in spite of the fact that the evoked fMRI signal is fairly weak relative to noise.'

%{\it Design of fMRI studies.}
%
%In addition to PET and fMRI, there are a number of other kinds of brain imaging techniques,
%each with distinctive advantages and limitations.

\subsection{Positron Emission Tomography (PET) Imaging}
PET is a nuclear imaging technique for mapping brain function, or other molecular processes in the body.
It measures emissions related to positron decay from radioactively labeled chemicals that have been injected into the bloodstream and generates
images of the distribution of the radio chemicals throughout the brain and body.

PET is also not an exact measure of brain function. It depends upon certain assumption about what happens when an area of the brain becomes active: 1) cerebral metabolism requires glucose metabolism, which requires oxygen from blood flow, i.e. in parts of the brain that are more active there is more blood flow; 2) cerebral blood flow varies locally with corresponding local variations in neural activity.

In PET, images of the brain's activity are generated by tracking the brain's use of a radioactively tagged compound, such as glucose, oxygen or a particular drug. A machine called a cyclotron (used to accelerate charged particles to high energies) is used to tag the specific substance with small amounts of radioactivity. Before scanning, the labeled compound (called a radiotracer) is injected into the subject's bloodstream (or inhaled). When the tracer gets into the bloodstream, it goes to areas of the brain that use it, while the subject is engaged in some type of mental activity. So, oxygen and glucose accumulate in brain areas that are metabolically active. When the radioactive material decays, it gives off a neutron and a positron. When a positron hits an electron, both are destroyed and two gamma rays are released approximately 180$^\circ$ opposed to each other. Sensors in the PET scanner detect and record the brain area where the gamma rays are emitted. By reconstructing the sites of the positron-electron collisions, the location of active regions can be imaged. A computer uses the data gathered by the sensors to construct multicolored 2- or 3-dimensional images that show where the compound is in the brain. Different colors or degrees of brightness on a PET image represent different levels of tissue or organ function. This method hence provides a functional view of the brain.

%{\it Advantages and Limitations of PET.}
%PET provides an image of brain activity. Disadvantages: Expensive to use. Radioactive material used, low temporal resolution.

PET has very high biochemical sensitivity and selectivity (which allow probing the neurochemical processes at the molecular level), but its temporal and spatial resolution are inferior to that of fMRI. The mean free path of the positrons in brain tissue limits the spatial resolution of PET scanning to about 4 mm. However, PET images can be superimposed on subject's MRI images, providing detailed information about specific brain areas involved in a wide variety of functions. Spatial resolution of PET data depends upon several other factors: the size and type of the crystal used in the scanner scintillator to detect the gamma radiation emission, %(the smaller the size the better the resolution?),
the energy of the positron emitted etc. Temporal resolution depends mainly upon half-life of the isotope. Safety regulations require to wait 5 half-lives between injecting the radioactive tracers.

%``A unique feature of PET is its biochemical selectivity, in that
%specific radiotracers can be designed that
%bind selectively to the molecular targets, such as receptors,
%transporters, or enzymes that are involved in the synthesis or
%metabolism of neurotransmitters. PET thus has been extensively
%used for biochemical and pharmacological imaging of
%the human brain, although it also has been widely used to
%identify anatomical correlates of cognitive process''.

{\it Experimental Design}.
In a typical PET experiment, a brain scan is taken during a control task (e.g., resting with eyes closed).
The control task is compared to brain scan taken while the subject is exposed to the experimental treatment or performing the experimental task.
To determine the brain activity that can be attributed to the experimental condition, the difference between the PET scans is calculated.

%==================================================
\section{Analysis of Functional Neuroimaging Data}
%==================================================
\label{FNI_analysis}
Typical objectives in the analysis of functional neuroimaging data are sought through the following:
(i) activation studies attempt to localize particular brain regions that are source(s) of task-related neural processing,
% essentially looking at different patterns of activity between different experimental tasks,
(ii) connectivity studies seek to identify what brain areas show similar patterns of activity
over time, yielding distributed networks of brain function,
and (iii) prediction studies try to use neuroimaging scans to predict future behavioral, psychological, or cognitive outcomes, or
to classify subgroups of individuals (e.g. for diagnostic purposes).
%In this proposal each of these objectives will be addresses.
%In Sections \ref{ss:activation_studies} -- \ref{ss:prediction} we give a brief summary of each of these three research areas. % highlights

Data from a functional neuroimaging study consist of a series of 3-D images, typically obtained while the subject
performs a certain cognitive, behavioral or emotional task, or while at rest.
%The tasks (stimuli) are carefully designed to ....
In an fMRI study, typically hundreds of such 3-D images are obtained, taken approximately 2-4 seconds apart. In a PET study, the number of scans
is significantly smaller (the maximum number of scans per subject is limited due to the total isotope dose allowed), and the brain images are
obtained much less frequently. Each 3-D image comprises of a large number ($>$ 100,000 in an fMRI study; usually smaller for PET) of voxels.
%In case of fMRI, the image intensity from each voxel represents the spatial distribution of the nuclear spin density in that area.
%In case of PET, the voxel intensity measures the regional cerebral blood flow (rCBF) as a proxy for brain activity.
In addition, the experiment may be repeated for the same subject, as well as for multiple subjects.
%Figure \ref{f:study_design} shows a simplified illustration of an interventional (repeated-measures) study, for a single subject.
Because of the neurophysiology of the network organization of the brain, spatial correlations are very likely. Also, temporal correlations (both within and between scanning sessions) are present due to repeated scanning.

%To remove uninteresting variability from the data; increase SNR, to prepare the data for statistical analysis
%A fair amount of preprocessing is required before any statistical analysis can be conducted.
To reduce artifacts caused by the scanner and subject movement, and to map individual brains to a standard brain atlas
(when performing group analysis and making population inference), the functional scans undergo several preprocessing steps
before statistical analysis.

\subsection{Preprocessing Pipeline}
\label{ss:preprocessing}
The main steps involved in fMRI preprocessing are: slice timing correction (shifting each voxel's time course so that we can assume they were measured simultaneously), realignment (to correct for motion; usually a rigid body transformation using 6 variable parameters is used), coregistration of structural and functional images (to map the results obtained from functional data to a high resolution structural image for presentation purposes), normalization (important in group analysis, attempts to register each subject's brain to a standard template brain (see Section \ref{ss:brain}), it is commonly done using non-linear transformations to match local features), and spatial smoothing (typically done by convolving the functional images with a Gaussian kernel). Spatial smoothing is done for several reasons: it may improve inter-subject registration, it ensures that the assumptions of the random field theory (often used for multiple testing correction) are satisfied, and to increase signal to noise ratio. However, spatial smoothing causes a loss of acquired data and may introduce artificial spatial correlations between nearby voxels.
Preprocessing steps and the order in which they are performed are important since they affect both the spatial and temporal correlation structure of the data.

%Preprocessing steps for the PET data are similar. % ask Tim or Dana for help; maybe write a paragraph
We refer to \cite{woods:98a}a and \cite{woods:98b}b for the details on the preprocessing steps for PET data.

\subsection{Statistical Modeling for Activation Studies}
\label{ss:activation_studies}
The general objective of activation studies is to identify brain locations that are involved in the
neural processing associated with tasks that subjects perform while in the scanner and possibly to compare
these neural processing traits between tasks or between subgroups of individuals.
Commonly, a two-stage statistical modeling procedure is used, consisting of single subject
analysis at stage 1 and a group level model at stage 2.
We present the modeling procedure in context of fMRI data, and PET analyses proceed using a similar framework.
To set the notation, assume that the data consist of a brain volume of $V$ voxels (indexed by $v=1\ldots,V$), repeatedly measured at
$T$ time points (indexed by $t=1,\ldots,T$). In addition, assume that the experiment is repeated for $K$ subjects
(indexed by $i=1,\ldots,K$).
\subsubsection{Single-Subject Analysis}
%{\it Modeling the measured fMRI data.}
%In this section we describe a typical single-subject statistical analysis of task fMRI data.
%For single-subject analysis of PET data, which includes coregistration of significant PET activation to every
%individual's MRI image and transformation into Talairach space, we refer to ?? for more details.
The first stage of the model characterizes distributed brain activity associated with various experimental
conditions, separately for each individual.
A single subject model for fMRI activation in a single voxel can be expressed as
\begin{equation}\label{eq:fmri_model}
y_{iv}(t)=\sum_{l=1}^L x_{ivl}(t) B_{ivl} + \sum_{m=1}^M h_{ivm}(t) \nu_{ivm} + \ve_{iv}(t),
\end{equation}
for $v=1,\ldots,V$, $i=1,\ldots,K$ and $t=1,\ldots,T$. Here $y_{iv}(t)$ represents the measured fMRI signal at time $t$
for voxel $v$ and subject $i$; $x_{vil}$ represents the task-related BOLD response (the signal of interest) corresponding
to $l$-th condition at time $t$;
and $h_{ivm}$ represents the contribution of the nuisance parameters (such as high-pass filtering
parameters - to remove the drift component typically present due to scanner instabilities, periodic fluctuations
due to heart rate and respiration, and head motion) at time $t$. The terms $B_{ivl}$ and $\nu_{ivm}$ represent the
unknown amplitude of $x_{ivl}$ and $h_{ivm}$ and $\ve_{iv}(t)$ represents the noise process.
Typically, the relationship between stimuli and BOLD response %(a complex nonlinear function of the results of neuronal and vascular %changes)
is modeled using a linear time invariant system where the stimulus acts as the input and the hemodynamic response function (HRF) as the impulse response function.
In a linear system framework the signal $x(t)$ is modeled as the convolution of a stimulus function $s(t)$ and the
hemodynamic response $f(t)$, i.e. $x(t)=(s * f)(t)=\int_0^\infty s(u)f(t-u)\, \text{d}u$.  %(Figure~\ref{fig:BOLD} gives a graphical illustration.).
%\begin{figure}\label{fig:BOLD}
%\begin{center}
%\includegraphics[scale=0.7]{Figures/BOLD_response.eps}
%\caption{Neural signal convolved with  a canonical HRF to produce the predicted BOLD signal.}
%\end{center}
%\end{figure}
%, where $h$ is either assumed to take a canonical form or is modeled as a set of linear basis functions.

In summary, \eqref{eq:fmri_model} can, for most fMRI experiments, be represented in matrix form as
%
\begin{equation}\label{eq:fmri_model_matrixF}
\mathbf{y}_i(v)=\mathbf{X}_{iv}\mathbf{B}_i(v)+\mathbf{H}_{iv}\boldsymbol{\nu}_i(v)+\boldsymbol{\ve}_i(v) ,
\end{equation}
%
where $\boldsymbol{\ve}_i(v)$ is usually assumed to follow an AR(2) process.
Here we arranged the serial BOLD fMRI responses for subject $i$ in the vector $\mathbf{y}_i(v)=(y_{iv1},\ldots,y_{ivT})'$
(at a given voxel location).
$\mathbf{X}_{iv}$ is a $T\times L$ matrix with columns corresponding to predicted BOLD response for each condition $l$, and
$\mathbf{B}_i(v)=(B_{iv1},\ldots,\be_{ivL})'$. If, for example, the nuisance matrix
represents a drift, modeled using a $p$th order polynomial function, then $\boldsymbol{\nu}_i(v)=(\nu_{vi1},\ldots,\nu_{vip})'$.
$\mathbf{H}_{iv}$ is a $T\times p$ matrix with columns corresponding to the polynomial functions.
%The model \eqref{eq:fmri_model_matrixF} can further be combined across voxels, and can be written as
%%
%\begin{equation}\label{eq:fmri_model_matrixF_1}
%\mathbf{Y}_i=\mathbf{X}_i \mathbf{B}_i + \mathbf{H}_i \mathbf{N}_i +\mathbf{E}_i.
%\end{equation}
%$\mathbf{Y}_i$ is a $T\times V$ matrix in which each column is a time series corresponding to a single brain voxel and
%each row is the collection of voxels that make up an image at the specific time. Matrices $\mathbf{X}_i$ and $\mathbf{H}_i$
%are the common design matrices used for each voxel, and $\mathbf{B}_i=(\bb_{1i},\ldots,\bb_{Vj})$,
%$\mathbf{N}_i =(\boldsymbol{\nu}_{1i},\ldots,\boldsymbol{\nu}_{Vj})$, $\mathbf{E}_i=(\boldsymbol{\ve}_{1i},\ldots,\boldsymbol{\ve}_{Vj})$.
%
%Typically, the vectorized variance of $\mathbf{E}$ is assumed to be separable in time and space. The spatial
%covariance is also often assumed to be negligible (compared to the temporal covariance) and therefore ignored.
%% (NOTE: this is almost verbatim from Martin, so need to rephrase).
%%
%Model \eqref{eq:fmri_model_matrixF_1} provides a full spatio-temporal model of brain activity, but is not considered
%feasible due to extreme computational demands and, instead, model \eqref{eq:fmri_model_matrixF} is applied to each voxel separately.

\vspace{1ex}
The {\it general linear model} (GLM)\citep{worsley:95, friston:02} is the most popular way to analyze
the single subject fMRI data. The GLM assumes that the
HRF is known apriori. If, in addition, the stimulus function $v_l(t)$ is known and equivalent to experimental paradigm (for example:
a vector of zero's and one's, representing times when the signal is `off' and `on', respectively), then
\eqref{eq:fmri_model_matrixF} becomes a standard multiple regression model with known signal components and
unknown amplitudes. HRF is usually modeled using a canonical HRF, typically a gamma function or a difference between
two gamma functions. The BOLD response is summarized in the design matrix $\mathbf{X}_{iv}$, containing separate
columns for each of the $l$ predictors. The parameters $\mathbf{B}_i(v)$ and $\boldsymbol{\nu}_i(v)$ represent subject specific
effects corresponding to $\mathbf{X}_{iv}$ and $\mathbf{H}_{iv}$, respectively.
The vector $\boldsymbol{\ve}_i(v)$ denotes random errors representing characteristics of the measurement process that
are unaccounted for by the model, and we assume that the errors follow a zero-mean multivariate normal distribution
$\boldsymbol{\ve}_i(v) \sim \text{N}(0, \si^2_i(v) \bI)$, where $\si^2_i(v)$ is the error variance at
voxel $v$.
% Temporal correlations - if space
% Serial measurements of brain activity exhibit temporal correlations. Model above does not explicitly address those correlations,
% but few commonly use strategies exist....

After fitting the model, a contrast vector $\mathbf{c}$ can be defined to test for an effect $\th_i=\mathbf{c}' \mathbf{B}_i$.
For example, we can test for the increased brain activity during the `on' signal versus during the `off' signal.
Hypothesis testing can be performed in a usual way: individual model parameters can be tested using a $t$-test and
subset of parameters using a partial $F$-test.

%A comprehensive summary of the analysis of fMRI data from the acquisition of the raw data to its role in
%localizing brain activity is given by \cite{lindquist:08}.
%its role in localizing brain activity
%\cite{lindquist:08} gives a very good summary of the analysis of fMRI data from the acquisition of the raw data to
%its role in localizing brain activity, stressing the crucial role of statistics in understanding the nature of the data and
%obtaining relevant results.
%
%\begin{figure}\label{fig:hrf}
%\begin{center}
%\includegraphics[scale=0.6]{Figures/hrf.eps}
%\caption[Typical BOLD response to an impulse stimulation]{Typical BOLD response to an impulse stimulation.
%The response peaks approximately 5 sec after stimulation, and is followed by
%an undershoot (at high magnetic files, an initial undershoot can sometimes be observed).
%}
%\end{center}
%\end{figure}
%
%Various approaches for single-subject analyses of PET data can be found in the literature.
%We refer to [refs] for details.

\subsubsection{Multi-Subject (Group Level) Analysis}
Typically, multi-subject analysis of fMRI data from an activation study is performed using hierarchical models, which provide a
framework for performing mixed-effects analysis.
%In studies investigating task-related changes in measured brain activity,
%two-stage modeling is commonly used, where the first stage involves subject-specific models and the second stage specifies group (or population) parameters.

The stage 2 model combines the subject-specific effects $\mathbf{B}_i(v)$ from equation \eqref{eq:fmri_model_matrixF} to estimate
the associated group level or population parameters. The second level model can then be written as
\begin{equation}\label{eq:second_level}
\mathbf{B}_i(v)=\mathbf{W}_{iv}\bb(v)+\mathbf{e}_i(v)
\end{equation}
where $\mathbf{W}_{iv}$ is the second-stage design matrix that often consists of subject-specific characteristics and $\bb(v)$ contains the group
level parameters representing the effects related to different sessions or the effects associated with different
subpopulations (e.g. treatment groups). Often, the interest is in making inferences about the differences between sessions/groups.
A common approach is to obtain $\th_i$ from the single subject analysis and
then model it in terms of the corresponding group level parameter.

%The first-level model can be written as
%\begin{equation*}
%\mathbf{y}=\mathbf{X} \mathbf{B} +\bve,
%\end{equation*}
%where $\mathbf{y}=(\mathbf{y}_1',\ldots,\mathbf{y}_K')'$, the design matrix $\bX=\text{diag}(\bX_1,\ldots, \bX_K)$,
%$\bB=(\bB_1',\ldots, \bB_K')'$, $\bve=(\bve_1',\ldots,\bve_K')'$ and $\text{Var}(\bve)=\text{diag}(\bV_1',\ldots,\bV_K')$.
%The second level can then be written as:
%$$
%\bB=\bX_{\text{gr}} \bb_{\text{gr}}+ \bve_{\text{gr}},
%$$
%where $\bX_{\text{gr}}$ is the second level (group) design matrix, $\bb_{\text{gr}}$ is the vector of second-level parameters, and
%$\bve_{\text{gr}} \sim N(0,\bI \si_{\text{gr}}^2)$.

%{\it Spatial modeling}
\subsection{Data-driven Descriptive Analysis Methods}
\label{ss:data_driven}
Data-driven approaches to the analysis of functional neuroimaging data employ pattern-based exploratory methods to
``characterize'' the nature of the signal present in the data.
A lot of recent work regarding the analysis of brain imaging data has focused on examining functional connectivity of the human brain.
 %Functional connectivity refers to the association or dependency of activation between regions. One way to examine functional connectivity is to use temporal correlations between spatially remote neurophysiological events. This represents a ''model free'' characterization of brain connectivity (Lee et al., 2003). A second, potentially more powerful method, uses effective connectivity to measure the interaction of one brain region with another, mediated by anatomical connections between them as opposed to the direct effect of task or activation in a common driving region (Friston et. al., 1997).
%Functional connectivity studies seek to identify the ares of brain that share similar temporal task-related neural responses.
Functional connectivity has been defined as the temporal correlation between spatially remote neurophysiologic events.
%A related term is {\it effective connectivity}, defined as measure of the interaction of one brain region with another, mediated by anatomical connections between them.
For example, Independent Component Analysis (ICA) and cluster analysis attempt to divide the brain into different functional
networks involved in performing certain tasks or characterizing a ``default mode'' of the resting-state brain.
\subsubsection{Clustering}
Cluster analysis uses PET or fMRI data to help identify dissociable networks or clusters, each consisting of voxels that show correlated
patterns of measured brain activity.  These clustering solutions do not define the underlying neuroanatomical connections,
but instead identify functional associations between voxels.  Ideally, intracluster voxels should exhibit high functional
(or spatial) autocorrelation, validating that the neural responses within clusters are functionally related.

Clustering procedures classify the $V$ voxels in an image into $G$ groups, with each cluster consisting of $V_g$ voxels, where $g=1,\ldots, G$ and
$V=\sum_{g=1}^G V_g$. Most methods classify voxels by measuring the distance (dissimilarity) between the activity time courses for every pair of voxels and combining voxels with small distances. Two typical dissimilarity measures are Euclidean and Mahalonobis distances, which are easy to compute and naturally interpretable.
%For example, one popular measure of distance between
%the activity in voxels $i$ and $j$ is the Euclidean distance given by
%$
%d\bigl(\bT(v_i),\bT(v_j)\bigr)=\bigl[\bigl(\bT(v_i)-\bT(v_j)\bigr)' \bigl(\bT(v_i)-\bT(v_j)\bigr)\bigr]^{1/2}.
%$

Numerous investigators have proposed the use of clustering methods for neuroimaging data including
\cite{balslev:02}, \cite{baumgartner:00}, \cite{bowman:04a}, and \cite{cordes:02}.
% ,fadili:00,fadili:01,filzmoser:99,goutte:99,goutte:01,poline:94,simon:04,sommer:03,stanberry:03}.
Clustering algorithms generally fall into one of two categories, namely hierarchical clustering algorithms or partitioning algorithms.
Hierarchical clustering performs a series of nested merges or divisions, ranging from a solution with one large cluster containing all voxels to a solution in which every voxel represents a separate cluster.  Among the available hierarchical clustering procedures are Ward's, beta-flexible, centroid linkage, variable linkage, median linkage, single linkage, and complete linkage algorithms \citep{bowman:04a}.
Partitioning algorithms, such as K-means and fuzzy K-means, specify the number of clusters at the outset and sequentially reallocate voxels to
clusters until obtaining a final clustering solution \citep{fadili:01}.
Descriptions of particular clustering algorithms are available in \cite{rencher:02}, and \cite{bowman:04b} highlight several of these algorithms for neuroimaging applications.

\subsubsection{Independent Component Analysis (ICA)}
\label{ss:ICA}
Among the data-driven techniques, ICA has become a very popular method, successfully employed to decompose functional neuroimaging data into
sets of spatial maps and associated time-courses. ICA is an application of a blind source separation that attempts to decompose the data set
into components that are as statistically independent from each other as possible \citep{common:94,herault:86}.
%As a method for recovering underlying signals, or independent components (ICs) from linear data mixtures, it can be applied to fMRI time-series to %spatially localize and temporally characterize the sources of blood-oxygen-level-dependent (BOLD) activation.
It is possible to pursue either the temporal %\citep{biswal:99}
or the spatial %\citep{mckeown:98}
independence of the target components.
% In this paper, we will concentrate on spatial ICA.
One method of performing ICA minimizes the mutual information between components \citep{bell:95,mckeown:97,mckeown:98,calhoun:00}.
The fixed-point algorithm \citep{hyvarinen:99} pursues the same goal of minimizing the mutual information, but uses the concept of normalized
differential entropy or negentropy \citep{common:94}. One interpretation of negentropy is as a measure of non-normality, so maximizing
the negentropy finds directions of maximal-non-normality in the data.
%\cite{beckmann:04} propose a probabilistic ICA (PICA) model aimed at solving the problem of overfitting in classical ICA applied to fMRI data,
%by including a Gaussian noise term in the classical ICA decomposition.
%
%Several methods for extending the single-session probabilistic ICA model to higher dimensions and making group inferences
%from fMRI data have been proposed \citep{calhoun:01,beckmann:05}.
%In our analysis, we consider the group PICA approach proposed by \cite{beckmann:05} to our experimental fMRI data.

\vspace{1ex}
\noindent{\it Classical ICA.}
Let $\bY$ denote a $S \times V$ matrix of the observed voxel time courses, where $S$ is the number of scans, and $V$ is the number of voxels included in analysis; $\mathbf{C}$ is a $N \times V$ random matrix whose rows $C_i$ are to be filled with the unknown realizations of the $N$ spatial components (images, reshaped in 1-D vectors), and $\mathbf{A}$ is a $S \times N$ mixing matrix, whose columns contain the associated time-courses of the N components. The spatial ICA problem for fMRI time series can be formulated as an estimation of the following model for the data:
\begin{equation}\label{ICA1}
\bY =\bA \bC
\end{equation}
There are no assumptions about the mixing matrix $\bA$, while the constraint on the spatial processes $C_i$ is that they are
(ideally) mutually statistically independent. The amount of statistical dependence within a fixed number of spatial components can
be quantified by means of their mutual information, an important function in information theory \citep{common:94}.
Thus, the ICA decomposition of $\bY$ can be defined (up to a multiplicative constant and to the sign) as an invertible
transformation: $\bC = \bW \bY$, where the matrix $\bW$ (so-called unmixing matrix) is determined such that the mutual
information of the target components $C_i$ is minimized (i.e., such that $C_i$'s are ``as independent as possible'').
Matrix $\bA$ is the pseudoinverse of $\bW$. In \cite{mckeown:98}, the sources are estimated by iteratively optimizing the
unmixing matrix $\bW$ so that $\bC = \bW \bY$ contains mutually independent rows, using the `infomax' algorithm.

\vspace{1ex}
\noindent{\it Probabilistic ICA.}
\cite{beckmann:04} propose a probabilistic ICA (PICA) model aimed at solving the problem of overfitting in classical ICA applied to fMRI data,
by including a Gaussian noise term in the classical ICA decomposition.
%The model of equation \eqref{ICA1}  does not include the noise model.
The PICA model, which extends model \eqref{ICA1}, is formulated as a generative linear latent variables model.
The model is ``characterized by assuming that the $S$-variate vector of observations is generated from a set of $q$ statistically
independent non-Gaussian sources via a linear instantaneous mixing process, corrupted by additive noise
$\boldsymbol{\eta}(t)$''\citep{beckmann:04}:
$$
\mathbf{y}_v=\mathbf{A} \mathbf{c}_v + \boldsymbol{\mu} + \boldsymbol{\eta}_v, \quad \forall v\in \{1,\ldots,V\}
$$
where $\mathbf{y}_v$ denotes the $S$-dimensional column vector of individual measurements at voxel location $v$, $\mathbf{c}_v$ denotes the $q$-dimensional column vector  of non-Gaussian source signals contained in the data and $\boldsymbol{\eta}_v$ denotes Gaussian noise $\boldsymbol{\eta}_v \sim N(0,\si^2\boldsymbol{\Si}_v)$.
The number of source processes is assumed to be smaller than the number of observations in time ($q < S$).
The model is closely related to Factor Analysis \citep{bartholomew:87}.
%Summarize briefly ML estimation of the mixing matrix and sigma2, model order selection, estimation of the unmixing matrix  and inference] In brief,
The individual component maps are obtained using a modified fixed-point iteration scheme to
optimize for non-Gaussian source estimates via maximizing the negentropy. The maps are then separately
transformed to spatial $Z$-scores. %In contrast to raw IC estimates, which only encode the estimated signal,
These $Z$-score maps depend on the amount of variability explained by the entire decomposition at each voxel
location relative to the residual noise. % similar to statistical parametric maps from a GLM analysis.
Next, Gaussian/Gamma mixture models are fitted to the individual $Z$-maps in order to infer voxel
locations that are significantly modulated by the associated time-course.
For more details on PICA, we refer the reader to \cite{beckmann:04}.

\vspace{1ex}
\noindent{\it Group ICA}.
Several methods for extending the single-session probabilistic ICA model to higher dimensions and making group inferences
%from fMRI data
have been proposed \citep{calhoun:01,beckmann:05,guo:08a}.
In \cite{beckmann:05}, the single-session PICA has been extended to higher dimensions allowing for a model-free analysis of multi-subject or multi-session fMRI data. For group ICA, either tensorial ICA (where the data is decomposed into spatial maps, time courses and subject/session modes), or a simpler temporal concatenation approach (a single 2-D ICA is run on the concatenated data matrix obtained by stacking all 2-D data matrices of every single data set on top of each other) are possible. The latter approach is recommended to use when the interest is in finding a common spatial pattern, but we cannot assume that the associated temporal response is consistent between subjects, such as in the analysis of data acquired without stimulation (e.g. resting-state).
%We hence apply this approach for IC analysis of our fMRI data set.
The technique is derived from the Parallel Factor Analysis \citep{harshman:70,harshman:94}.
%[Briefly summarize temporal and tensorial approach if enough space?].
\footnote{The research by Beckmann and Smith described above has been implemented as MELODIC (Multivariate Exploratory Linear Optimized Decomposition into Independent Components - a stand alone C++ program). It is freely available as part of FSL (FMRIB's Software Library).}
%
\subsection{Prediction and Classification}
\label{ss:prediction}
In recent years there has been growing interest in the use of neuroimaging data as a tool for classification
(e.g. classification of mental and behavioral disorders) and prediction (e.g., predicting an early onset of Alzheimer's disease).
%Typical objectives of prediction and classification methods applied to functional neuroimaging data include:
%(i) predicting behaviors (many attempts have been made, e.g. \cite{pugh:96,hoeft:07,giessing:07),
%%(ii) clinical diagnoses, or (iv) outcomes of treatments.
%(ii) predicting future patterns of brain activity, and predicting health outcomes from brain activity data (e.g, patients clinical response to
%treatment).
%
Typical prediction objectives in functional neuroimaging studies include: (i) predicting human experiences and behaviors from brain imaging data, (ii) predicting brain activity patterns, e.g. following a treatment intervention, using pre-treatment data and patient characteristics, and (iii) predicting treatment response and relapse. %, e.g. in substance use disorder research. %(prediction of treatment effectiveness)
Many attempts have been made to address (i), e.g. \cite{pugh:96}, \cite{hoeft:07}, \cite{giessing:07}.
There have also been attempts to address (ii). For example, \cite{guo:08} introduced a Bayesian hierarchical
model for PET and fMRI data to forecast brain activity in schizophrenic patients following a
specific treatment.
%This work represents an important first step in attempting to help treatment decisions through use of functional neuroimaging data.
%and provides a foundation on which future research, including that on predicting clinical symptom responses, can build.
%Some other related work includes [refs].
%Predicting patient response to treatments has also been recently addressed
Functional neuroimaging has recently begun to show promise as a clinical tool in the prediction of treatment response.
For example, \cite{evans:06} demonstated functional neuroimaging to be promising as a clinical tool in the prediction of treatment
response in psychiatric disorders such as major depression and obsessive-compulsive disorder. The findings they reviewed
suggest that treatment outcome may be predicted by patterns of pre-treatment brain activity in psychiatric patients.
They conclude, however, that the actual clinical utility of such tests remains to be shown.
%
%Also, there has been growing interest in predicting stimuli directly from functional data (refs, LaConte)
% (``this would allow for the possibility to infer information from the scans about the subject's though process and use
% brain activation patterns to characterize subjective human experience.''- Lindquist)
% and for monitoring of disease progression.
%
%One of the promising directions in neuroimaging is in increasing its clinical applicability in terms of diagnostic
% purposes and to help predict future patient's outcomes
%This field is very new and limited methods were used for this purpose (only standard methods e.g. logistic classifiers and SVC.
%Standard methods typically involve feature selection or a data reduction step and from that point each feature is combined as an
% independent piece of information to contribute to classification
%Idea of predicting is still relatively new and there haven't been any attempts to perform formal statistical testing
%Typical candidates for statistical methods are logistic and SV classifiers (they consider each voxel separately)
%

%- short section: some of the examples: discriminant analysis (Fisher); SVM (give ref's: what it is, give specific appl. in NI);
% logistic regression; Bayes classifiers.

%The most popular classification methods in neuroimaging use linear classifiers, which include
%Logistic Regression (known as Multinomial Regression if there are more than 2 classes), linear Support Vector Machines (SVM),
%as well as the classification decisions of Gaussian Naive Bayes (GNB) and Fisher's Linear
%Discriminant Analysis (LDA). SVMs arise from the Statistical Learning Theory of Vapnik (Vapnik, 1995) and possess
%several unique properties appropriate for real world applications, including fMRI.

Popular classification methods in functional neuroimaging include Fisher's Linear (and Quadratic) Discriminant
Analysis (LDA) \citep{fisher:36}, logistic regression (LR), %(known as multinomial regression if there are more than 2 classes),
Support-Vector Machines (SVM) \citep{vapnik:95}, Neural Network classifiers (NN) \citep{hertz:91}, Partial Least Squares (PLS) analysis
\citep{wold:66,mcintosh:96}, and Bayesian classifiers (e.g. naive Bayes classifier).

%Recent studies have demonstrated that spatial patterns of fMRI BOLD activity distribution over the brain may be used to classify
%different groups or mental states. These studies are based on the application of advanced pattern recognition approaches
%and multivariate statistical classifiers.
%Pattern-based classifiers (which include SVM, NN, LDA) are being applied to the analysis of functional neuroimaging data with
%increasing frequency \citep{laconte:05,pessoa:07}.

Pattern recognition methods (which include SVM, NN, LDA) have been successfully applied in functional neuroimaging studies
\citep{laconte:05,pessoa:07}. These methods can be used to infer cognitive states (so-called brain decoding).
Using such approaches, it is possible to predict the mental state of a subject or a stimulus class by analyzing the spatial
distribution of neural responses.
%In addition it is possible to identify the regions of the brain containing the information that underlies the classification.
The SVM is one of the most popular methods used to carry out this type of analysis.
SVMs arise from the Statistical Learning Theory of Vapnik \citep{vapnik:95} and their
formulation  was motivated to deal with small sample sizes and high dimensional inputs, which match the
situation involved for temporally predictive modeling of fMRI data.
%The goal of pattern-based classification of functional neuroimaging data is to link individual brain activation
%patterns to the experimental conditions experienced during the scans.
%Recent studies have demonstrated that spatial patterns of fMRI BOLD activity distribution over the brain may
%be used to classify different groups or mental states.
\cite{cox:03} used statistical pattern recognition algorithms including LDA and SVM
to separate brain activation maps from an fMRI experiment in which participants viewed images of objects (baskets, birds, butterflies,
chairs, cows, etc.). The objects belonged to various categories, both of similar and differing forms.
A notable finding of their study was that it was possible, with small amounts of fMRI data,
to determine what the participant was viewing at levels well above chance.

Hanson et al. (2004) implemented a variety of feed-forward, neural network
architectures with both linear and nonlinear decision rules.
PLS, introduced to neuroimaging community by \cite{mcintosh:96}, is a useful method for forming prediction equations when there are a large number of explanatory variables,
particularly when the random error variance is large. \cite{giessing:07} applied PLS analysis to study behavioral effects of nicotine in nonsmokers
and found that neural data under placebo can be used to predict individual behavioral effects of nicotine.
\cite{mitchell:04} present case studies in which they have successfully trained three classifiers:
a Gaussian Naive Bayes classifier, k-nearest neighbor, and linear SVMs, to distinguish
cognitive states such as whether the human subject is looking at a picture
or a sentence or whether the word the subject is viewing is a word describing food, people, buildings, etc.
% brain states prediction <-- Pessoa ref, important to mention
%Classifiers such as LDA and SVM are pattern-based and multivariate. The correlation classifier is pattern-based, but not multivariate

%========= patter recognition ================
%
%Pattern-based classifiers are being applied to the analysis of functional neuroimaging data with increasing frequency.
%The goal of pattern-based classification of functional neuroimaging data is to link individual brain activation
%patterns to the experimental conditions experienced during the scans.
% These analyses ask a different question than either inferential or exploratory approaches.
%They combine the quantitative advantages of inferential statistics while retaining a pattern-based analysis of the scans.
%
%Pattern recognition methods (including SVM, NN, LDA) have been successfully applied in functional neuroimaging studies
% \citep{laconte:05,pessoa:07}. These methods can be used to infer cognitive states, so-called brain decoding.
% Using such approaches, it is possible to predict the mental state of a subject or a stimulus class by analyzing the spatial
% distribution of neural responses. In addition it is possible to identify the regions of the brain containing the information
%that underlies the classification. The Support Vector Machine (SVM) is one of the most popular methods used to carry out this type of analysis.
%
%Pattern-based classification algorithms have been applied recently to the analysis of functional neuroimaging data in a
%way that links patterns of brain activity to the experimental conditions experienced during the scans
% (Carlson, Schrater & He, 2003; Cox & Savoy, 2003; Hanson, Matsuka & Haxby, 2004; Haxby, Gobbini, Furey, Ishai, Schouten, & Pietrini, 2001;
%Haynes and Rees, 2005; Mitchell, Hutchinson, Niculescu, Pereira, Wang, Just, & Newman, 2004; O'Toole, Jiang, Abdi, & Haxby, 2005;
%Spiridon & Kanwisher, 2002; Yukikasu & Tong, 2005).

%------ conclusion --------
%Most studies use linear classifiers, while others employ nonlinear classifiers. (see refs in Prediction of Successful Memory Encoding from
%fMRI Data)
Methodology for statistical prediction based on functional neuroimaging data represents an important area for
future research, and preliminary work in this area provides a promising outlook for the potential usefulness of functional neuroimaging
data in a clinical setting.
%
%Clearly, functional neuroimaging is playing an important role in substance use disorder research.
For the functional neuroimaging modality, however, to be useful in defining diagnostic categories or monitoring treatment
success, it needs to be shown that the technology has the ability to define clinically relevant information on a single-subject basis.

\section{Motivating Examples}
\label{motivating_examples}
\subsection{An fMRI Study on Inhibitory Control in Cocaine Addicts}
A common characteristic of drug addiction is an impairment in the ability to exert inhibitory control over
drug-related behaviors, in spite of adverse consequences. This inhibitory
control deficit can be elicited using response inhibition tasks outside of the drug seeking context.
The inhibitory control task, referred to as the {\it STOP-signal} task, was designed to evaluate the ability to cancel
a prepotent motor response \citep{aron:06}. The sample included 12 cocaine addicts who received an extensive outpatient
behavioral treatment program and 15 healthy controls. The cocaine addicts were scanned while
performing an inhibitory control task in two separate sessions, before and after treatment,
and similarly control subjects had baseline and follow-up scans.

The aim of our analysis is to identify brain locations where there are inhibitory control-related increases in brain activity
following treatment for cocaine addiction. In our analysis we take into account that there are both spatial correlations
present between neighboring voxels in the brain, as well as temporal correlations between repeated brain scanning sessions.
%
\subsection{A PET Study on Working Memory in Schizophrenia Patients}
Schizophrenia is a devastating mental illness that may cause severe disturbances in cognition, social behavior, and emotion.
One important brain function that often goes awry in patients with schizophrenia is working memory, characterized by the ability to manipulate information in short-term memory.
Our data are from a parametric study of working memory in an executive function task in individuals with schizophrenia.
There were 16 subjects in the study, each having a total of 8 PET scans obtained using a blood flow tracer.
The aim is to examine the functional connectivity related to working memory tasks and to evaluate the validity of the neural
processing networks identified by functional clustering. % or by ICA.
%
\subsection{An fMRI Resting-state Study of Depression}
This data consist of fMRI resting-state scans of seven women with a history of major depression. Functional data were acquired while the subjects were in the scanner with their eyes open and while looking at a visual fixation cross.
A group independent component analysis (ICA) was applied to investigate the resting-state functional connectivity of a group of depressed patients. There is substantial interest in determining associations between the brain activity characteristics of different regions while subjects are in a resting state. For example, a set of regions has been consistently identified in these investigations and has been labeled as the {\it default mode network} (DMN). Our goal is to introduce a measure for the degree of functional autocorrelation within neural processing networks identified by ICA and to evaluate the statistical significance of the observed associations.

\subsection{A PET Study of Alzheimer's Disease}
Alzheimer's Disease Neuroimaging Initiative (ADNI) \\
({\tt http://www.loni.ucla.edu/ADNI/}) is a large national project with the goal to develop biomarkers of Alzheimer's Disease (AD) in elderly subjects, to define the rate of progress of mild cognitive impairment and Alzheimer's disease, and to provide a large database which will improve design of treatment trials. We analyzed PET data from 73 healthy controls and 73 Alzheimer's disease patients, obtained from the ADNI database.
Our goal is to predict subject's follow-up (6 month) brain activity, based on the baseline activity.


%%========================================
%\section{Previous Work/Current Methods}
%%=======================================
%%====================================
%\subsection{Spatio-temporal Modeling}
%%====================================
%Spatial-temporal models arise when data are collected across time as well as space.
%
%%==================================
%\subsection{Functional Connectivity}
%%==================================
%%%\subsection{Resting-state Data}
%\subsection{Data-driven Techniques}
%\subsubsection{Clustering}
%\subsubsection{Independent Component Analysis (ICA)}

%==================================
\section{Proposed Research}
\label{proposed_research}
Functional neuroimaging studies yield large data sets characterized by complex dependence structures
driven by highly sophisticated neurophysiology and neuroanatomy, and aspects of the experimental designs.
These complex dependence structures pose analytical challenges for statistical modeling.
To date, there is a paucity of methods that incorporate spatial considerations. % spatial correlations present in functional neuroimaging data.
Most analyses are performed in a univariate setting, where each voxel is modeled independently of the others
(``massive univariate approach'').
Ignoring correlations affects the precision of estimates of model parameters and consequently may lead to inaccurate
statistical tests. Spatio-temporal modeling may mitigate these shortcomings by incorporating more physiologically
plausible assumptions and by borrowing strength across related measures of neural activity.
%==================================
\subsection{Simultaneous Spatio-temporal Modeling of fMRI data}
fMRI data sets contain temporal correlations from repeated scanning (within and between scanning sessions)
and complex spatial correlations. Typical analyses investigating task-related changes in measured brain
activity use a two-stage procedure in which the first stage involves
subject-specific models % relating neural processing to experimental tasks
and the second-stage specifies group (or population) level parameters.
Customarily, the first-level accounts for temporal correlations between
the serial scans acquired during one scanning session.
Despite accounting for these correlations, fMRI studies often include multiple sessions
%(e.g., before and after treatment),
and temporal dependencies may persist between the corresponding estimates of mean neural activity.
Further, spatial correlations between brain activity measurements in different locations are often
unaccounted for in statistical modeling and estimation.
\cite{bowman:05} proposed an extended two-stage model for the estimation and testing of localized activity
in which the second stage accounts for spatial dependencies between voxels within the same neural
processing cluster (defined by a data-driven cluster analysis). This model, however, did not account
for repeated measures type associations between the multiple experimental effects for each subject.
We propose a two-stage, spatio-temporal, autoregressive model which simultaneously
accounts for spatial dependencies between voxels within the same anatomical
region and for temporal dependencies between a subject's estimates from multiple
sessions.
%The model can be written in matrix form as follows:
%\begin{equation*}
%\bB_{igs} = \bb_{gs} +\rho_{gs} \mathbf{W}^* (\bB_{igs} -\bb_{gs}) + \xi_{gs} \mathbf{M}^*
%(\bB_{igs} -\bb_{gs}) +\boldsymbol{e}_{igs}\, ,
%\end{equation*}
%where $\boldsymbol{e}_{igs} \sim \text{MVN} (\mathbf{0},\phi_{gs}^2 \mathbf{I})$;
%add if not too long:
% The approach can also be applied to account for correlations between estimates of mean neural
% activity from multiple experimental conditions/stimuli from a single session, e.g.
% related to a particular cognitive domain.
We develop an algorithm that leverages the special structure of our covariance model, enabling relatively
fast and efficient estimation.
%
%We propose an extended two-stage, spatio-temporal, simultaneous autoregressive model which
%accounts for both spatial dependencies between voxels within the same anatomical region and for
%temporal dependencies between the multiple experimental effects/sessions for a subject. We develop
%an algorithm that leverages the special structure of our covariance model, enabling relatively fast
%and efficient estimation.

\subsection{Functional Autocorrelation within Neural Processing Networks}
Data-driven statistical approaches, such as cluster analysis or independent component analysis, applied to in vivo functional
neuroimaging data help to identify neural processing networks that exhibit similar task-related or resting-state patterns of activity. Ideally, the measured brain activity for voxels within such networks should exhibit high autocorrelation.  An important limitation of the existing approaches is that they do not quantify or statistically test the strength or nature of the within network relatedness between voxels.  To extend the results given by such data-driven analyses, we propose the use of Moran's $I$ statistic to measure the degree of functional autocorrelation within identified neural processing networks and to evaluate the statistical significance of the observed associations.
We adapt the conventional definition of Moran's $I$, for applicability to neuroimaging analyses, by defining the global autocorrelation index using network-based neighborhoods.  Also, we compute network-specific contributions to the overall autocorrelation.

\subsection{A Novel Spatial Prediction Model}
Neuroimaging has been used in clinical practice for over 30 years, but it is still perceived as rarely offering clinicians
much help in direct patient management. One of the important directions in neuroimaging is increasing its clinical applicability
in terms of diagnostic purposes and to help predict future patients' outcomes.

We propose a novel Bayesian hierarchical framework for predicting follow-up neural activity based on the baseline functional neuroimaging data that attempts to overcome some shortcomings of the presently used modeling methods (such as small sample sizes, heterogeneity, as well as recording of brain activity under resting-state conditions) by borrowing strength from the spatial correlations present in the data. The spatial correlations are incorporated in the model in two ways: the short-range correlations between neighboring voxels are incorporated through a multivariate conditional autoregressive (CAR)  prior of the spatial parameters, while the long-term correlations between anatomical brain regions are incorporated through the covariance matrix of the random effect parameters.

