Locally Linear Embedding and fMRI feature selection in psychiatric classification
Abstract
Background:
Functional magnetic resonance imaging (fMRI) provides noninvasive measures of neuronal activity using an endogenous Blood OxygenationLevel Dependent (BOLD) contrast. This article introduces a nonlinear dimensionality reduction (Locally Linear Embedding) to extract informative measures of the underlying neuronal activity from BOLD timeseries. The method is validated using the LeaveOneOutCrossValidation (LOOCV) accuracy of classifying psychiatric diagnoses using restingstate and taskrelated fMRI.
Methods:
Locally Linear Embedding of BOLD timeseries (into each voxel’s respective tensor) was used to optimise feature selection. This uses Gauß’ Principle of Least Constraint to conserve quantities over both space and time. This conservation was assessed using LOOCV to greedily select time points in an incremental fashion on training data that was categorised in terms of psychiatric diagnoses.
Findings:
The embedded fMRI gave highly diagnostic performances (> 80%) on eleven publiclyavailable datasets containing healthy controls and patients with either Schizophrenia, AttentionDeficit Hyperactivity Disorder (ADHD), or Autism Spectrum Disorder (ASD). Furthermore, unlike the original fMRI data before or after using Principal Component Analysis (PCA) for artefact reduction, the embedded fMRI furnished significantly better than chance classification (defined as the majority class proportion) on ten of eleven datasets.
Interpretation:
Locally Linear Embedding appears to be a useful feature extraction procedure that retains important information about patterns of brain activity distinguishing among psychiatric cohorts.
I Introduction
Over a century ago, Charles Darwin alluded to an experimental paradigm that involved direct observation of the brain’s physical mechanisms (nervous matter) [14], where such observations [55] would serve as the physical basis for dichotomising speciesspecific behaviour. In the early 90s, an indirect and noninvasive measurement of mental activity over uniformlyspaced time points became possible through functional Magnetic Resonance Imaging (fMRI), which allows paramagnetic deoxyhemoglobin to act as an endogenous Blood OxygenationLevel Dependent (BOLD) contrast [49]. This article pursues Darwin’s proposed experimental paradigm by Locally Linear Embedding (LLE) [54] the BOLD timeseries to produce precise summaries of cerebral activity that may optimise the classification of different brain states, such as mental disorders.
Logothetis et al found that an increase in invasivelymeasured neural activity directly and monotonically reflects local BOLD signal increases and, for short stimulus presentations, there is a linear relationship between BOLD and neural responses [39]. This suggests the unobservable neural activity is spatiallylocalised in anatomical space. Locally Linear Embedding of fMRI data in space and time can, in principle [2, 20, 24] (see Appendix Locally Linear Embedding), summarise these local measurements of neuronal mass activity [40]–via the notion of analytic capacity [45]–to disclose information about global (i.e., wholebrain) activity patterns.
Ii Methods: Testing Locally Linear Embedding (LLE) with Cross Validation
The discriminatory power of Locally Linear Embedding was compared to the original fMRI, both before and after applying Principal Component Analysis (PCA) [31] for artefact reduction [39]. This comparison was performed on eleven datasets containing cohorts with different mental disorders, using a combination of Leave One Out Cross Validation (LOOCV) on the training set, with greedy feature selection based on Fisher discriminability [17, 18]. The purpose of using LOOCV and feature selection is to find the time points that discriminate patients from controls on the respective dataset. The feature selection step initially starts with an empty candidate set of time points and proceeds to select the time point with the highest discriminatory power [63]. Then, time points that improve discrimination in conjunction with those already in the candidate set are added to this set in incremental fashion; this process is terminated when there are no more time points that can be added to the candidate set to improve discriminability. Note that the selection of time points is based upon cross validation and does not induce any biased sampling.
To illustrate the patterns that best discriminate between groups, a paired twosample ttest between the patient and control groups is performed to both threshold and identify the statisticallysignificant differences (p 0.05 uncorrected) in space (at the time points identified by the greedy feature selection). Per Mill’s Methods of Induction (Method of difference), the functional differences (depicted by the statistically significantlydifferent regions) at the respective time point are therefore a necessary part of the cause of the phenomena that distinguish the subject groups [46], which in this case pertain to a neuropsychiatric disorder. Neuropsychiatric disorders are diagnosed using clinical assessments that include: evaluating the background demographics, collecting first and third party observations, and a structured psychiatric interview with the subject [53]. In detail:
Every subject’s fMRI timeseries is treated as a fourdimensional array with dimensional voxel waveforms for . Assume each subject scan is associated with a binaryvalued class label representing the diagnosis and that, for any subject scan, every voxel waveform is generated by a vector corresponding to a point on the manifold. Our approach to fMRIbased diagnosis involves two stages:
fMRI reconstruction takes the subject’s fMRI as input and outputs a reconstructed fMRI that is more informative than the original. Formally, this reconstruction is a mapping
In other words, Locally Linear Embedding reduces a timeseries of length T to a smaller number of spatial modes of dimensionality d; these modes contain all the information used for the subsequent step.
Classification builds a classifier that takes the subject’s reconstructed fMRI as input, and outputs a class label . The classifier is therefore a mapping .
The reconstructed, or reduced, fMRI data produced from step 1 is hereon referred to as . All reconstructions initially vectorise the fMRI data to produce a twodimensional array , and conclude by reshaping the resulting twodimensional reconstruction into a fourdimensional array .
Principal Component Analysis (PCA) [5, 31]
reconstructs by finding an orthogonal rotation that minimises the reconstruction cost
where is the mean over all voxel waveforms, and . To find the optimal , compute the righthand matrix for the singular value decomposition (SVD) of , which contains the rightsingular vectors of [26]. It follows that is the rotated matrix that minimises the reconstruction cost of the subject’s fMRI, where every column of is a principal component. It is assumed the first principal components (columns) of capture the “systematic structure", where the confounding factors are relegated to the remaining T  d principal components, which produces the twodimensional reconstructed fMRI .
Using PCA’s “systematic structure" for distinguishing humans with different neurological disorders has been met with caution [16], largely because PCA’s application to fMRI has some subjective components [1]. The main limitation behind the PCA reconstruction is that it assumes the lowerdimensional manifold is a linear subspace. We demonstrate that introducing the Cauchy stress tensor [9] on the Cartesian space with the Pythagorean distance metric enables threedimensional measurements over time, thereby revealing the local (group) action in the physical system.
Locally Linear Embedding (LLE) [54, 56]
reconstructs by constructing the Cauchy stress tensor [9] at every voxel for , which is achieved by minimising the reconstruction cost of its waveform in terms of its spatiallyadjacent neighbours:
(1) 
where the neighbourhood set for voxel is the complement of its spatial neighbours on the surface of the sphere with radius r, and are the reconstruction weights containing the spatiallyinvariant geometric properties of the Cauchy stress tensor at voxel i.
To determine the analytic capacity [45] at every voxel location i, LLE first subtracts the K spatial patterns of the voxels on the boundary of the sphere centred around voxel i to determine the separation distance from the origin of the tensor at the respective voxel. Then, it computes the local (symmetric) spatiotemporal covariance matrix:
(2) 
where is a nonnegative regularisation term to enforce positivedefiniteness (for this study, ). LLE calculates the reconstruction weights by finding the unique minimumnorm solution [50] to the constrained leastsquares problem defined by:
(3) 
where is a vector of ones and the element of can be thought of as the average height of a curve representing the mean transit time of the indicator [34] of voxel from voxel over the duration of the scan. Since represents the squared distance of the surface forces from voxel , the reconstruction weights are Lebesgue measures [36] summarising the analytic capacity, or spatiallyinvariant geometry [25], of the spacefilling curve [27], where the constraint ensures that the areas between the imaginary surface (acting as the origin that divides the body) and curves (defined by the stress vectors) are 1 in each of the directions. In practice the weights can be brittle [65] due to any number of reasons. Modified Locally Linear Embedding (MLLE) therefore computes the linearlyindependent (orthogonal) vectors^{1}^{1}1When using MLLE it is possible for , thus the optimal number of weight vectors for each voxel is determined by setting so that K  1 –i.e., is set to span as large of a basis as possible. After this step the desired dimensionality d is then input to the eigensolver. of using the eigendecomposition , thereby allowing the definition of multiple weight vectors for each voxel. Assuming the columns (eigenvectors) are sorted in descending order of their respective eigenvalues , MLLE uses the first columns to compute multiple local weight vectors for a single voxel:
where if (else ), , , , and . Since every measure’s invariant properties are determined in a squareintegrable space [58], LLE performs a global leastsquares optimisation based on Gauß’ Principle of Least Constraint [23] to calculate the vectors corresponding to points on the manifold:
(4) 
where is the dimensionality parameter selected by the user and is the twodimensional reconstructed fMRI. The global optimisation therefore calculates the points on the manifold [59, 57] that act as the fourdimensional orthogonal basis that best retains the geometry of the stress vectors. These represent the second order invariant properties of each voxel’s Cauchy stress tensor. A detailed explanation of this optimisation is provided below.
Define as the local sparse adjacency matrix, where:
The optimisation in Equation 4 can be written as a minimisation of the expected reconstruction cost, or error:
(5) 
where is the orthogonal projection for voxel , and is the sparse, symmetric and positivedefinite alignment matrix, and therefore admits the eigendecomposition [42]:
(6) 
where is the diagonal matrix containing the smallest eigenvalues of , and are the corresponding eigenvectors; Rayleigh’s variational principle [52, 12] enables calculation of these bottom eigenvectors. Each eigenvector represents a degree of freedom in space and time, where the eigenvector is the global unit vector that fills threedimensional space. The global unit vector is discarded to enforce the constraint that the manifolds have mean zero.
Note: To avoid degenerate solutions, LLE requires the manifolds to be centred around the origin in both space and time – i.e., and – and also have outer products with unit covariance – i.e., . Centring the manifolds about the origin ensures they are of the same scale, which is superficially similar to the common practice of signal, or count rate, normalisation [51]. The unit covariance constraint imposes the requirement that the reconstruction errors of the extracted manifolds are measured on the same scale.
Dataset  Partition  Specificity  Sensitivity  Precision  Accuracy  

Chance  Original  LLE  PCA  Chance  Original  LLE  PCA  Chance  Original  LLE  PCA  Chance  Original  LLE  PCA  
Beijing  Training  100%  73.3% 15.8%  80% 14.3%  66.7% 16.9%  0%  53.3% 17.9%  93.3% 8.9%  80% 14.3%  0%  66.7% 16.9%  82.4% 13.6%  70.6% 16.3%  50% 17.9%  63.3% 17.2%  86.7% 12.2%  73.3% 15.8% 
(Peking_3)  Holdout  100%  75%  87.5%  25%  0%  100%  50%  25%  0%  66.7%  66.7%  14%  67%  83.3%  75%  25% 
COBRE  Training  100%  57.7% 13.6%  92.3% 7.3%  96.2% 5.3%  0%  60% 13.4%  88% 8.9%  52% 13.7%  0%  57.7% 13.6%  91.7% 7.6%  92.9% 7.1%  51% 13.7%  54.9% 13.7%  90.2% 8.2%  74.5% 6.1% 
Holdout  100%  60%  100%  100%  0%  70%  100%  60%  0%  63.6%  100%  100%  50%  65%  100%  80%  
Mind Research  Training  100%  N/A  88.2% 11.5%  100%  0%  N/A  84.6% 12.9%  23.1% 15.1%  0%  N/A  84.6% 12.9%  100%  56.7% 17.7%  N/A  86.7% 12.2%  66.7% 16.9% 
Network (MRN)  Holdout  100%  N/A  67%  100%  0%  N/A  71%  14%  0%  N/A  63%  100%  56%  N/A  68.8%  62.5% 
Stanford  Training  100%  73.3% 15.8%  80% 14.3%  100%  0%  66.7% 16.9%  86.7% 12.2%  33.3% 16.9%  0%  71.4% 16.2%  81.3% 14%  100%  50% 17.9%  70% 16.4%  83.3% 13.3%  66.7% 16.9% 
Holdout  100%  80%  80%  100%  0%  40%  60%  0%  0%  67%  75%  0%  50%  60%  70%  50%  
University of  Training  100%  69.2% 20.2%  100%  100%  0%  42.9% 21.7%  71.4% 19.8%  14.3% 15.3%  0%  42.9% 21.7%  100%  100%  65% 20.9%  60% 21.5%  90% 13.1%  70% 20.1% 
Michigan (UM_2)  Holdout  100%  66.6%  100%  83.3%  0%  25%  50%  0%  0%  20%  100%  0%  60%  30%  80%  50% 
Dataset  Experimental  Pulse Sequence  Patient Disorder  Source  Training Data  Holdout Data  fMRI Dimensionality  Parameters  Age Range  

Design  total  patients  controls  total  patients  controls  ()  LLE (,)  PCA ()  [min,mean,max]  
Beijing (Peking_3)  restingstate  EPI  ADHD  ADHD200  30  15  15  12  4  8  (2,236)  236  [11, 13.24, 16]  
COBRE  restingstate  EPI  Schizophrenia  COBRE  51  25  26  20  10  10  (2,150)  150  [18, 37.45, 65]  
Mind Research Network (MRN)  blockdesign  EPI  Schizophrenia  MCIC  30  13  17  16  7  9  (2,177)  177  [18, 34.15, 60]  
Stanford  restingstate  SpiralIN/OUT EPI  Autism  ABIDE  30  15  15  10  5  5  (2,60)  60  [7.52, 9.95, 12.93]  
University of Michigan (UM_2)  restingstate  SpiralIN EPI  Autism  ABIDE  20  7  13  10  4  6  (2,222)  222  [13.1, 16.26, 28.8] 
Iia Feature & Parameter Selection
Hyperparameter Selection involves selecting the algorithm parameters that will be used to generate the reconstructed fMRI from LLE and PCA, respectively. Both PCA and LLE require one to specify the number of reconstruction dimensions , where LLE’s additional hyperparameter, , selects the neighbouring voxels whose coordinates are on the boundary of a sphere with radius that is centred around the respective voxel. For both PCA and LLE, the number of reconstruction dimensions is bounded by the number of time points. In the case of PCA, the optimal d is often chosen by the proportion of variance captured by the first d eigenvectors, expressed as . For LLE, however, there is no analogous interpretation because the d eigenvectors are uniformly spaced “time points" on the world line [47]. Thus crossvalidation procedures are implemented (see following subsection for details) on the training data to systematically select the best hyperparameters, iterating over , where is generated on a logscale from 1 to .
Sequential Forward Selection (SFS)^{2}^{2}2Implemented using MATLAB’s (R2013a) sequentialfs function [63]
is a nonparametric method for measurement (feature) selection that starts with an empty “candidate", or nearoptimal, set of “time points". The method first finds the “time point", defined as the second moment [6], with the highest classification accuracy and adds the corresponding image volume to the set. The method then finds an additional “time point" that strictly improves the classification accuracy in conjunction with the volume(s) whose “time points" are already in the “candidate" set, and terminates when no such volume can be found. Thus, SFS produces a nearoptimal set of “time points" (volumes) that distinguish the subject groups with high classification accuracy, where this set is determined by adding “time points" to the near optimal set in a onebyone fashion. It follows that the nearoptimal set of “time points" represent points in time during the scan that enable discrimination of patients from controls with high classification accuracy. Note that SFS is only used on the training data (see Section IIC).
IiB Classifier
A slight abuse of notation is introduced by redefining variable as the onedimensional representation of the c diagnostic volumes from the respective subject’s reconstructed fMRI produced in the previous step. Here, it is assumed there are c diagnostic time points and that is a random variable containing the collection of cohorts’ reconstructed fMRI.
Fisher’s Linear Discriminant (LDA)^{3}^{3}3Implemented using MATLAB’s (R2013a) classify function with the ‘diaglinear’ argument to estimate the positive diagonal covariance matrix [17, 18]
is a linear classification rule [26] that assumes both the patient and control class densities (at each voxel location) can be represented as multivariate Gaussians in threedimensional space, each with some intrinsic curvature [21, 22]. Each class density is expressed as
where it is assumed the classes have a common covariance matrix– i.e., . This assumption allows the log ratio between the posterior distribution of each class to form a decision boundary that lies between patients (class 1) and controls (class 0), written as , which is linear in . LDA therefore calculates the dimensional hyperplane that best discriminates, or separates, the diagnostic volumes that have been determined using SFS with LOOCV on the respective dataset, for each class’ reconstructed fMRI. It follows these diagnostic volumes contain spatial locations with sufficiently different patient and control class densities, where statisticallysignificantly different regions possess the requisite margin between the class densities such that they are perceptible [51].
IiC Evaluation Criteria
Each dataset was split such that the proportion of patients in the training partition was nearequal to the proportion of controls, and the holdout dataset contained at least 10 cohorts. In all but two cases, the aforementioned criteria could not be met because the datasets were unbalanced. For these situations, the training dataset was constructed such that it was a representative sample of the overall data, with the remaining cohorts being assigned to the holdout dataset. The holdout set therefore contained group proportions that could deviate from the training set. The cohorts used to define the training and holdout partition for each dataset are provided in Table 3 of the Supplementary Information.
Performance was compared to the original fMRI, both before and after applying PCA for artefact reduction [39], and chance, which is defined as the proportion of the majority class on the respective dataset. The performance of the reconstruction methods are evaluated using a combination of LeaveOneOut Cross Validation (LOOCV) and holdout data classification performance [26]. LOOCV is used on the training data to find both the best reconstruction parameter , and the diagnostic volumes that produce the highest accuracy for this . The holdout data classification accuracies use the parameters found from performing LOOCV on the training data. Note: the cohorts in the holdout set are never involved in determining the optimal hyperparameter , or the diagnostic volumes produced from this .


Iii Results
Performance & Visualisation
Tables I and II demonstrate that Locally Linear Embedded fMRI can distinguish various mental disorders from healthy controls with high discriminatory power ( 80%); results for six additional restingstate datasets are provided in Table 2 of the Supplementary Information. The datasets contained healthy controls and patients with either Schizophrenia, AttentionDeficit Hyperactivity Disorder (ADHD), or Autism Spectrum Disorder (ASD). Ten of these datasets contained cohorts in the restingstate, with the remaining containing schizophrenic patients and healthy controls performing the Sternberg Item Recognition Paradigm (SIRP) task [61]. For a description of the data sources and technical details, please consult the appendix.
Given that performance on the training partition uses Leave One Out Cross Validation (LOOCV) to individually predict each subject’s diagnosis, which is analogous to performing n Bernoulli trials [4] (where n is the number of cohorts), the training data’s performance metrics can be interpreted as the mean of successes over n binomiallydistributed observations. To calculate the error of these estimates, we follow Laplace’s approach of employing a normal distribution to estimate the error of binomiallydistributed observations [33]. Given that discrimination performance on the holdout partition is determined in a onetime fashion, variance estimates are not applicable.
The HarvardOxford Subcortical/Cortical and Cerebellum atlases [30] are used to identify the statisticallysignificant differences between patients and controls in the diagnostic volumes for the respective dataset, shown in Figures 1, 2, and 3 (Figures 1, 2, 3, 4, 5, and 6 for datasets in the Supplementary Information). Each dataset’s figure contains six different views of the statisticallysignificant differences at the time reflected by the respective time point(s); the coloured voxels at these time points denote statistically significant (p < 0.05 uncorrected) physical differences between the patient and control groups, where these groups include cohorts from both the training and holdout partitions. The proportion of significantly different voxels in each region, calculated by dividing the number of significantly different voxels by the total number of voxels in the respective region defined by the atlas, are provided in Tables 4 and 5 of the Supplementary Information.
Neurobiological Interpretation
Schizophrenia
David Ingvar & Göran Franzèn found that healthy controls exhibit increased flows in the prefrontal regions and decreased flow in the post central regions, and schizophrenic patients exhibited the reversed pattern, with low flows prefrontally and high flows postcentrally [28]. Furthermore, they noticed that a lower flow in the premotor and frontal regions was associated with symptoms of indifference, activity and autism, and a higher postcentral flow over the temporooccipitoparietal regions was associated with disturbed cognition [28]. Inspecting the statistical maps for the volumes in Figure 1, the significantly different areas seem to further substantiate Ingvar & Franzèn’s observations, as there are significant differences in the various temporooccipitalparietal regions in all of the volumes.
The Sternberg Item Recognition Paradigm (SIRP) task evaluates cohorts’ shortterm, or working, memory. Each of the seven tasks during the scan involves the subject memorising a set of objects, followed by presentation of a new object whose membership in this set is identified by a ‘yes‘ or ‘no‘. The tasks therefore evaluate the hypothesised information processing differences between schizophrenics and healthy controls [16]. It has been shown that the prefrontal and medial temporal regions are involved in encoding information, and it is believed the interactions between these regions are central to retrieval of stored information [60]. Figure 1, especially volume 136, illustrates significant differences in the areas associated with the prefrontal and medial temporal regions; this is further supported by the fact that the accuracy on the holdout data rose from 68.8% to 75% when using only volume 136. Thus, it is possible that schizophrenia indeed affects the physical mechanisms associated with retrieving stored information, as these mechanisms are central to the SIRP task. The reconstruction method therefore successfully reveals physical differences associated with task performance between patients and controls, which are different from the restingstate differences for the same subject groups.
AttentionDeficit Hyperactivity Disorder (ADHD)
Similar to the goals of Ingvar & Franzèn [28], previous work used PET scans to compare the regional cerebral blood flow of children with AttentionDeficit Hyperactivity Disorder (ADHD) to healthy controls, where it was found that the disorder was associated with hypoperfusion in the striatal and posterior periventricular regions [41]; these results provide biological evidence that is consistent with the canonical model for ADHD as a frontostriatal deficient disorder. Figure 2 shows significant differences in the various occipital, striatal, cerebellar and ventral regions of the brain.


Autism Spectrum Disorder (ASD)
Similar to Ingvar & Franzèn’s observations that schizophrenic patients and healthy controls had normal hemisphere flows [28], studies using PET to compare the regional cerebral blood flow of Autism Spectrum Disorder (ASD) patients to healthy controls observed normal metabolism and blood flow. Hypoperfusion in the temporal lobes, centred in the associative auditory and adjacent multimodal cortex [66], was observed in autistic children. Furthermore, this temporal hypoperfusion was individually identifiable in 75% of autistic children [66]. Figure 3 illustrates statistical differences in many areas of the temporal lobe.
In summary, Locally Linear Embedding appears to have conserved spatiotemporal patterns in restingstate fMRI (and taskrelated responses) that are consistent with literature on regionallyspecific abnormalities of cerebral activity in the psychiatric conditions used to assess classification performance.
Iv Discussion
The multidisciplinary nature of this work undoubtedly introduces difficulty when discussing its motivations, which is the deployment of this methodology in a clinical setting. Such a goal imposes some conditions. First, while the results support deployment, the methodology must be further evaluated by trained clinicians wellversed in the etymology of the disease under investigation. More importantly, however, experimental designs are compulsory when discovering biological markers for disease; that is, patients must be subject to the same stimulus presentation at the same time during the scan in order to homogenise comparisons. With respect to restingstate fMRI, it is felt that a consensus is required to glean neurobiological insight that can generalise, which makes further discussion more appropriate for future work.
This exposition focused on Locally Linear Embedding as a promising and effective form of dimensionality reduction as a preprocessing step for the analysis of fMRI timeseries. The approaches and aims of this form of preprocessing share a close relationship with other approaches in imaging neuroscience. These approaches include Independent Component Analysis, the use of Support Vector Machines (and regression) to classification problems (and prediction), and algorithms based on adaptive smoothing. In future work, it will be interesting to explore the formal connections between other approaches and assess their relative sensitivity in the context of the classification problems considered above.
One hundred and fourtyeight years after Darwin ascertained that mental activity invokes physical mechanisms in the brain [14], the brain of man and ant alike are among the marvellous collection of atoms in the world.
Acknowledgements
This work is dedicated to the late Sam Roweis (19722010), Donald MacCrimmon MacKay (19221987) and his son David John Cameron MacKay (19672016). Additionally, many thanks to colleagues Ruitong Huang & Dr. Csaba Szepesvàri, Joshua T. Vogelstein, Dr. Vincent D. Calhoun, Dr. Neil Lawrence, Dr. Bert Vogelstein, Dr. Michael Milham, Dr. Karl J. Friston, Dr. Klaas Enno Stephan, Dr. Nikos K. Logothetis, Dr. Christof Koch, and Dr. Yifan Hu of AT&T Labs’ Information Visualisation department for his assistance with GraphViz, which was used to depict the Cauchy Stress Tensor in three dimensions.
Locally Linear Embedding
fMRI contain T uniformlyspaced time points, where time point describes a threedimensional space comprised of voxels (volumetric elements); this is the global description. Each voxel is a volumetric measurement of the brain’s physical mechanisms in both space and time
– i.e., every voxel represents a measurement at spatial coordinates at time , which can also be expressed as a fourtuple . The BOLD contrast over T uniformlyspaced time points can be viewed as a random variable, or measurable function, [6] that is caused by local neuronal fluctuations. The terms global and local are only valid under a rigorous definition of threedimensional space. A topological space requires specification of the relationships between points in a set and their respective open sets, which are defined as the sets that generalise the concept of an open interval in the real line while also providing a rigorous definition for nearness of points in this space [7].


Locally Linear Embedding (LLE) [54] is a congruent transformation [13] that first extracts the local mathematical (geometric) structure of the data, and then performs a global optimisation that best conserves this local latent structure [3], or extension [38]. The application of LLE is challenging, primarily due to the implicit assumption that each data point and its neighbours lie on, or close to, a locally linear subspace [56]; earlier applications of LLE to fMRI discarded spatial properties of the data in the neighbourhood selection step [44], thereby failing to preserve the inherent spatial configuration. In contrast, specifying the topological space allows LLE to treat each data point’s open set as its neighbourhood, thereby forgoing the neighbourhood selection step entirely. It follows that the benefit of LLE depends on the suitability of the distance metric used in determining the neighbourhood of each data point. Figures 4 and 5 suggest that the physical reality represented by topological space is a collection of locallylinear subspaces whose properties are subject to change at any time.
LLE uses the local description to construct the Cauchy stress tensor [9] at every voxel, which proceeds as follows: Initially, the spatial pattern of each voxel i, given by , is subtracted from both itself and the spatial patterns of its neighbours , thereby allowing its “zero waveform" (represented as ) to serve as the imaginary plane, or deflexion axis [43, 15], that divides the spherical body; the subtracted spatial patterns, given by (vector space) , represent the distances (over time) from the respective voxel’s spatial location.
We then use the innerproduct for to compute the squared distances, which results in the realsymmetric positivedefinite spatiotemporal covariance matrix for every voxel i. Since the contact forces are inversely proportional to the squared differences represented in elements of [48], we use the MoorePenrose inverse [50] to solve for the minimumnorm solution () to the constrained least squares problem . Here, is the unitlength direction vector in each of the directions, and weights represent the stress vectors originating from voxel i. By the spectral theorem, there exists an orthogonal basis such that and [62]; thus the stress vectors represented by are shift, rotation, and translation invariant in the space defined on basis . The substantial overlap between spatiallyadjacent voxels’ tensors precludes independently calculating orthogonal bases , as this defines each tensor’s basis on separate vector spaces.
The stress vector for each voxel’s Cauchy stress tensor, , represents the spatiallyinvariant properties of the threedimensional forces applied over the duration of the scan in the respective subspace. LLE retains these spatiallylocal invariant properties by constructing the (global) adjacency matrix such that ; by construction contains every voxel’s stress vectors, which were determined using the Cauchy stress tensor. Given that the stress vectors’ invariant properties are determined using squared distances, LLE then computes the squared distances between every voxel in the (global) normalised Laplacian matrix [10] . This is expressed by , where [35]. By construction, is symmetric and positivedefinite, which means there exists an orthogonal basis such that [42], where represents the global unit vector in threedimensional space.^{4}^{4}4 is also the eigenvector of , which is discarded after the global optimisation. Since the Cauchy stress tensor is of second order, LLE uses Rayleigh’s variational principle [52, 12] to calculate the resonance frequencies that best preserve the geometry of the deflexion axis at every voxel’s stress tensor. These are given by the bottom (d+1) eigenvectors, each of which represent one degree of freedom in space and time. LLE therefore uses a global optimisation to embed the relative measurements of every voxel’s Cauchy stress tensor [9] in a global coordinate system that conserves quantities over both space and time [19], where this coordinate system contains a mechanical system [32] in static equilibrium [8]. Applying LLE to fMRI data can therefore be viewed as using Carl Friedrich Gauß’ Principle of Least Constraint [23] to determine the true motion of the mechanical system defined on the topology , where the Cauchy stress tensor [9] allows preservation of the intrinsic local Gaussian curvature [21, 22] in space and time.
References
 [1] Andersen, A. H., Gash, D. M., and Avison, M. J. Principal Component Analysis of the dynamic response measured by fMRI: a generalized linear systems framework. Magnetic Resonance Imaging, 17(6):795 – 815, 1999.
 [2] Archibald, R. C. Time as a fourth dimension. Bulletin of the American Mathematical Society 20, (8):409–412, 1914.
 [3] Bacon, F. Francisci de Verulamio, summi Angliae cancellarii, Instauratio magna. apud Joannem Billium, Typographum Regium, 1620.
 [4] Bernoulli, J. Ars Conjectandi. Impensis Thurnisiorum, Fraetum, 1713.
 [5] Bishop, C. M. Pattern Recognition and Machine Learning. Springer, 2006.
 [6] Billingsley, P. Probability and measure. Wiley, 1979.
 [7] Bourbaki, N. Éléments de mathématique. Topologie algébrique. Chapitres 1 à 4. SpringerVerlag Berlin Heidelberg, Paris, 2016.
 [8] Carathéodory, C. Untersuchungen über die grundlagen der thermodynamik. Mathematische Annalen, 67(3):355–386, 1909.
 [9] Cauchy, A. L. De la pression ou tension dans un corps solide. Exercices de Mathèmatiques 2 (1827), 42–56.
 [10] Chung, F. R. K. Spectral Graph Theory. American Mathematical Society, 1997.
 [11] Cormen, T. H., Leiserson, C. E., Rivest, R. L., and Stein, C. Introduction to Algorithms (2. ed.). MIT Press, 2001.
 [12] Courant, R. Variational methods for the solution of problems of equilibrium and vibrations. Bulletin of the American Mathematical Society 49, (1):1–23, 1943.
 [13] Coxeter, H. S. M. Introduction to geometry (2. ed.). Wiley, 1989.
 [14] Darwin, C. The Descent of Man. London: Murray, 1871.
 [15] Deeley, E. M. and MacKay, D. M. Multiplication and division by electronicanalogue methods Nature 163, 4147 (4 1949), 650–650
 [16] Demirci, O., Clark, V. P., Magnotta, V. A., Andreasen, N. C., Lauriello, J., Kiehl, K. A., Pearlson, G. D., and Calhoun, V. D. A review of challenges in the use of fMRI for disease classification. Brain Imaging Behav 2, 3 (Sep 2008), 147–226.
 [17] Fisher, R. A. The Use of Multiple Measurements in Taxonomic Problems. Annals of Eugenics 7, 2 (1936), 179–188.
 [18] Fisher, R. A. The Statistical Utilization of Multiple Measurements. Annals of Eugenics 8, 4 (1938), 376–386.
 [19] Friston, K. J. The freeenergy principle: a unified brain theory? Nature Reviews Neuroscience 11, 2 (2010), 127–138.
 [20] Garabedian, P. R., Schiffer, M. On existence theorems of potential theory and conformal mapping. Annals of Mathematics, 52(1):164–187, 1950.
 [21] Gauß, C. F. Disquisitiones generales circa superficies curvas. Typis Ditericianis, Germany, 1828.
 [22] Gauß, C. F. Theoria motus corporum coelestium in sectionibus conicis solem ambientium. Perthes et Besser, 1809.
 [23] Gauß, C. F. Ueber ein allgemeines grundgesetz der mechanik. Journal für die reine und angewandte Mathematik (1829), 232–235.
 [24] Gillespie, R. J. Fifty years of the vsepr model. Coordination Chemistry Reviews, 252(12–14):1315–1327, 7 2008.
 [25] Gregory, J. Geometriae pars universalis. Patavii: Typis heredum Pauli Frambotti, Italy, 1668.
 [26] Hastie, T., Tibshirani, R., and Friedman, J. The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second Edition, 2nd ed. 2009. corr. 3rd printing 5th printing. ed. Springer Series in Statistics. Springer, Sept. 2009.
 [27] Hilbert, D. Ueber die stetige Abbildung einer Line auf ein Flächenstück Mathematische Annalen 3 38 (1891), 459–460
 [28] Ingvar, D. H., and Franzén, G. Distribution of cerebral activity in chronic schizophrenia. The Lancet 304, 7895 (12 1974), 1484–1486.
 [29] Ingvar, D. H., and Lassen, N. A. Quantitative determination of regional cerebral bloodflow in man. The Lancet 278, 7206 (10 1961), 806–807.
 [30] Jenkinson, M., Beckmann, C. F., Behrens, T. E. J., Woolrich, M. W., and Smith, S. M. Fsl. Neuroimage 62, 2 (Aug 2012), 782–790.
 [31] Jolliffe, I. Principal component analysis. Springer Verlag, New York, 2002.
 [32] Koch, C., and Hepp, K. Quantum mechanics in the brain Nature 440, 7084 (2006), 611612.
 [33] Laplace, P. Théorie analytique des probabilités. Courcier, Paris, 1812.
 [34] Lassen, N. A., and Ingvar, D. H.; Potchen, E. J., McCready, V. R. (Eds.) Radioisotopic assessment of regional cerebral blood flow. Progress in Nuclear Medicine 1: NeuroNuclear Medicine. University Park Press, 1972. 376–409
 [35] Lawrence, N. D. A unifying probabilistic perspective for spectral dimensionality reduction: Insights and new models. Journal of Machine Learning Research 13 (June 2012), 1609–1638.
 [36] Lebesgue, H. Leçons sur l’intégration et la recherche des fonctions primitives: professées au Collège de France. Collection de monographies sur la théorie des fonctions. GauthierVillars, 1904.
 [37] Liddle, P. F., Friston, K. J., Frith, C. D., Hirsch, S. R., Jones, T., and Frackowiak, R. S. Patterns of cerebral blood flow in schizophrenia. The British Journal of Psychiatry 160, 2 (1992), 179–86.
 [38] Locke, J. An Essay Concerning Humane Understanding. Thomas Basset, 1690.
 [39] Logothetis, N., Pauls, J., Augath, M., Trinath, T., and Oeltermann, A.. Neurophysiological investigation of the basis of the fMRI signal. Nature 412, 6843 (7 2001), 150–157.
 [40] Logothetis, N. K. What we can do and what we cannot do with fMRI. Nature 453, 7197 (6 2008), 869–878.
 [41] Lou, H. C., Henriksen, L., and Bruhn, P. Focal cerebral dysfunction in developmental learning disabilities. The Lancet 335, 8680 (1 1990), 8–11.
 [42] MacKay, D. J. C. Information Theory, Inference, and Learning Algorithms Cambridge University Press, 2003
 [43] MacKay, D. M. A highspeed electronic function generator Nature 159, 4038 (3 1947), 406–407
 [44] Mannfolk, P., Wirestam, R., Nilsson, M., Stahlberg, F., and Olsrud, J. Dimensionality reduction of fMRI time series data using locally linear embedding. MAGMA, 23(56):327–338, Dec 2010.
 [45] Mel’nikov, M. S. Analytic capacity: discrete approach and curvature of measure Sbornik: Mathematics 186, (6):827846, 1995
 [46] Mill, J. S. A System of Logic, Ratiocinative and Inductive: Being a Connected View of the Principles of Evidence and the Methods of Scientific Investigation, Volume I. John W. Parker, 1843.
 [47] Minkowski, H. Raum und Zeit. Physikalische Zeitschrift 10 (1908), 75–88.
 [48] Newton, I. Philosophiae naturalis principia mathematica. J. Societatis Regiae ac Typis J. Streater, 1687.
 [49] Ogawa, S., Lee, T. M., Kay, A. R., and Tank, D. W. Brain magnetic resonance imaging with contrast dependent on blood oxygenation. Proceedings of the National Academy of Sciences 87, 24 (1990), 9868–9872.
 [50] Penrose, R. A generalized inverse for matrices. Mathematical Proceedings of the Cambridge Philosophical Society, 51:406–413, 7 1955.
 [51] Popham, M. G.; Potchen, E. J., McCready, V. R. (Eds.) Numerical Methods for the Detection of Abnormalities in Radionuclide Brain Scans. Progress in Nuclear Medicine 1: NeuroNuclear Medicine. University Park Press, 1972. 117141
 [52] Rayleigh, J. W. S. The theory of sound. MacMillan and Co. London, 1877.
 [53] Robins, L. N., and Helzer, J. E. Diagnosis and clinical assessment: The current state of psychiatric diagnosis. Annual Review of Psychology 37, 1 (1986), 409 – 432.
 [54] Roweis, S. T., and Saul, L. K. Nonlinear dimensionality reduction by locally linear embedding. Science 290, 5500 (2000), 2323–2326.
 [55] Roy, C. S., and Sherrington, C. S. On the regulation of the bloodsupply of the brain. Journal of Physiology 11, 12 (Jan 1890), 85–158.
 [56] Saul, L. K., and Roweis, S. T. Think globally, fit locally: unsupervised learning of low dimensional manifolds. Journal of Machine Learning Research 4 (Dec. 2003), 119–155.
 [57] Seung, H. S., and Lee, D. D. The manifold ways of perception. Science 290, 5500 (2000), 2268–2269.
 [58] Schmidt, E. Über die Auflösung linearer Gleichungen mit Unendlich vielen unbekannten Rendiconti del Circolo Matematico di Palermo 25 1 (1908), 53–77.
 [59] Sherrington, C. S. Man on His Nature. Gifford lectures, Edinburgh. University Press, 1951.
 [60] Simons, J. S., and Spiers, H. J. Prefrontal and medial temporal lobe interactions in longterm memory. Nature Reviews Neuroscience 4, 8 (08 2003), 637–648.
 [61] Sternberg, S. Highspeed scanning in human memory. Science 153, 3736 (1966), 652–654
 [62] Sylvester, J. J. A demonstration of the theorem that every homogeneous quadratic polynomial is reducible by real orthogonal substitutions to the form of a sum of positive and negative squares Philosophical Magazine Series 4 (4) 23 (1852), 138–142
 [63] Whitney, A. A direct method of nonparametric measurement selection. Computers, IEEE Transactions on C20, 9 (sept. 1971), 1100 – 1103.
 [64] Williams, J., Gibbon, M., First, M., Spitzer, R., Davies, M., Borus, J., Howes, M., Kane, J., Pope, H., and Rounsaville, B. The structured clinical interview for dsmiiir (scid). ii. multisite testretest reliability. Archives of general psychiatry 49, 8 (08 1992).
 [65] Zhang, Z. and Wang, J. MLLE: Modified Locally Linear Embedding Using Multiple Weights Advances in Neural Information Processing Systems 19 (2007): 1593
 [66] Zilbovicius, M., Boddaert, N., Belin, P., Poline, J., Remy, P., Mangin, J., Thivard, L., Barthélémy, C., and Samson, Y. Temporal lobe dysfunction in childhood autism: a PET study. American Journal of Psychiatry 157, 12 (2000), 1988–1993.