WO2004051568A1

WO2004051568A1 - Brain connectivity mapping

Info

Publication number: WO2004051568A1
Application number: PCT/GB2003/001970
Authority: WO
Inventors: Timothy Edward John Behrens; Heidi Johansen-Berg; Stephen Mark Smith; Paul Mcmahan Matthews
Original assignee: Isis Innovation Limited
Priority date: 2002-11-29
Filing date: 2003-05-08
Publication date: 2004-06-17
Also published as: GB0227889D0; AU2003233884A1

Abstract

Computer-implemented techniques for mapping of the connectivity of the nervous system of a mammal use data derived from diffusion-weighted magnetic resonance imaging. Diffusion data (5) representing a probability density function for the diffusion parameters at respective voxels in a subject region of the nervous system is derived and used to derive a connectivity pattern data set (7) describing the anatomical connectivity of said nervous system to a seed voxel. Seed voxels are classified on the basis of their respective connectivity pattern data sets (7) according to a predetermined criterion such as the probability of connection of the seed voxels to portions of a target region. The techniques allow analysis of the anatomy of the nervous system because the unique connectivity pattern of a brain region determines the type of information available to it and therefore influences its function. The use of probabilistic techniques provides for robust and effective analysis and imaging.

Description

Brain Connectivity Mapping The present invention relates generally to mapping of the connectivity of the nervous system, and in particular the brain, of a mammal, in particular a human. It relates specifically to analysing the anatomical connectivity using imaging data, in particular derived from diffusion-weighted magnetic resonance imaging (DWI). Different aspects of the invention relate to computer-implemented techniques for deriving anatomical connectivity patterns and analysing the structure of the nervous system.

The unique connectivity pattern of a brain region determines the type of information available to it and therefore influences its function. It is desirable to define these connectivity patterns to enhance, not only our knowledge of human brain architecture, but also our ability to form hypotheses about regional brain function.

Although invasive tracer studies have produced a large body of evidence concerning connectivity patterns in non-human animals, our knowledge of brain connections in the human is much less complete. Methods for investigating human brain connectivity are largely limited to post-mortem studies with direct tracing of large tracts by dissection or histological studies of remote degeneration following focal injury. As such studies are limited by access to appropriate patients, this approach has led to an incomplete picture of human brain connections. It is therefore desirable to perform in vivo investigations of white matter pathways in the living human brain to extend our knowledge of human brain connectivity to complement functional imaging studies and better characterise patient groups in whom brain circuitry is thought to be abnormal. As an example, an important structure for brain connectivity studies is the thalamus. The thalamus is a deep grey matter structure through which nearly all incoming sensory information to the cortex is routed. It behaves as a central processor and relay station for sensory, motor and cognitive information. It is divided into cytoarchitectonically and neurochemically distinct nuclei that have widely different patterns of connectivity to cortical and other subcortical regions. The anatomy of cortico-thalamic and thalamo-cortical connections is well known for many non-human species. However, given the existence of inter-species differences in nucleic and cortical volumes and in thalamo-cortical connectivities, it may not be sensible to assume that connections established in the monkey brain will also be found in the human.

Boundaries between thalamic nuclei can be visualised histologically, but are not apparent using conventional magnetic resonance imaging (MRI). Some advances have recently been made in visualising major thalamic subdivisions using unconventional MRI sequences that selectively increase contrast-to-noise by using an inversion time which nulls signal from grey matter. However, the extent to which finer sub-divisions will be able to be defined is not clear and this method provides no information on potential functional relationships.

Recent developments in diffusion-weighted magnetic resonance imaging (DWI) techniques have enabled tracing of fibre tracts in the living human brain. DWI characterises the self-diffusion properties of water. In tissue with a high degree of directional organisation, the self-diffusion of protons is restricted more in some directions than in others. For example, in white matter the principal diffusion direction corresponds well to the major direction of connecting fibres passing through the voxel. In principle, by following the course of such principal diffusion directions it is possible to reconstruct fibre pathways in the brain and form a type of connectivity pattern representing the anatomical connectivity of the nervous system. However, until now, applications of DWI tractography have been limited to visualisation of major fibre bundles including the pyramidal and optic tracts, the corpus callosum and the superior longitudinal fasciculus. Known 'streamlimng' tract tracing algorithms can typically only progress when there is high certainty of fibre direction (i.e. when local diffusion anisotropy exceeds a specified threshold). This has therefore limited their ability to define pathways in or adjacent to grey matter.

According to the present invention, there is provided a computer- implemented method of analysing the structure of the nervous system of a mammal from imaging data representing local anatomical connectivity of a subject region of the nervous system, the method comprising: deriving from said imaging data, in respect of each of a plurality of first voxels in subject region of the nervous system, a connectivity pattern data set representing the anatomical connectivity between said first voxels and a at least one second voxel of said nervous system; and classifying the first voxels on the basis of their respective connectivity pattern data sets according to a predetermined criterion.

This method involves deriving the anatomical connectivity between a plurality of first voxels and at least one second voxel and classifying the first voxels on the basis of that connectivity pattern. Such a method provides for analysis of the structure of the nervous system. In particular, it allows one to differentiate voxels according to their connectivity to voxels elsewhere in the nervous system. As the function of the brain is dependent on the connectivity pattern to other parts of the nervous system, such classification has the capability of providing a powerful technique for structural analysis. This local connectivity information is used to define boundaries between sub-regions in the thalamus and provides information directly relevant to function. By considering not only diffusion measurements within a voxel, but also the mutual information available at a more global scale, it is possible to achieve thalamic segmentation at a finer spatial resolution than that of the DWI images. For example, applying the method to first voxels in the thalamus and taking the at least one second voxel in part of the cerebral cortex, it is possible to classify the thalamic voxels in a manner which closely corresponds to the human thalamic nuclei described from histo logical studies. On this basis, it is understood that this technique provides a reliable inference of the anatomy of voxels based on their connectivity pattern to other regions of the nervous system. In general, the technique may be applied to any first voxels and at least one second voxel which are inter-connected, but typically they will be separate regions of grey matter separated by white matter. The at least one voxel may be a single voxel or a region of the nervous system comprising plural voxels, for example a lobe, or even the entire subject region. Preferably, each connectivity pattern data set represents the probability of connection between the respective first voxel and the at least one second voxel.

Use of the probability of connection between voxels is a particularly effective way to represent the connectivity pattern. It allows the representation of complicated connectivity patterns, for example crossing or diverging paths. Another advantage is to provide effective recognition of patterns which are relatively weak, for example due to noise in the data set which noise is inevitable to some extent or due to the nature of the pattern being studied, eg where a connection path fans outwardly near grey matter or crosses other, stronger connection paths. It also accommodates situations where the boundaries between commonly-connected seed voxels are not sharp.

The predetermined criterion on which seed voxels are classified may be the probability of connection to a plurality of pre-selected regions each comprising at least one second voxel. For example, seed voxels may be classified according to the pre-selected region having the highest probability of connection or to each of the pre- selected regions for which probability exceeds a threshold.

The probability of connection is the preferred criterion for classifying the first voxels, because it is the most straightforward to implement. However, it is possible to apply other criteria. In particular, the seed voxels may be classified on the basis of any feature of the connectivity pattern, for example the shape, route or strength of connection paths, as well as the destination. The choice of criterion will depend on the anatomical nature of the structure it is desired to study. As studies advance, it is expected that the criteria will become more complicated.

There has been reported the use of local water diffusion properties of the nervous system to anatomically classify grey matter, for example the differentiation of thalamic regions with similar local diffusion properties. However, the present invention involves classification based, not on local properties, but on connectivity patterns extending across a plurality of voxels. This improves the ability to distinguish between areas with apparently similar local diffusion properties. The target portion will be separated from the seed voxels by a distance of at least the size of the seed voxels, preferably at least an order of magnitude greater than the size of the seed voxels.

The imaging data may be any imaging data that contains information about anatomical connectivity. For example, diffusion-weighted magnetic resonance imaging (DWI) characterises the self-diffusion properties of water. In tissue with a high degree of directional organisation, the self-diffusion of protons is restricted more in some directions than in others. For example, in white matter the principal diffusion direction corresponds well to the major fibre direction in the voxel. Therefore, the method is particularly applicable to diffusion-weighted magnetic resonance imaging data. To derive the connectivity pattern data set from the imaging data, such as diffusion-weighted magnetic resonance imaging data, any suitable method may be applied. For DWI imaging data, such methods are based on the known concept that paths in the nervous system are connected in the direction of diffusion. Most techniques use a maximum likelihood approach, for example by fitting a diffusion tensor to the imaging data. The tensors are used to trace connectivity paths. The first aspect of the present invention applies to connectivity patterns derived using such known techniques. However, particular advantage is achieved by using a connectivity pattern data set derived by the probabilistic method in accordance with the second aspect of the present invention which will now be described. According to the second aspect of the present invention, there is provided a computer-implemented method of deriving, from diffusion- weighted magnetic resonance imaging data of a subject region of the nervous system of a mammal, a connectivity pattern data set describing the anatomical connectivity of said nervous system, the method comprising: deriving, from said diffusion-weighted magnetic resonance imaging data, diffusion data representing a probability density function for the diffusion parameters at respective voxels in said subject region of the nervous system; and deriving from said diffusion data, as said connectivity pattern data set, a data set representing connectivity of a pre-selected seed voxel in said subject region of the nervous system to a target portion of said nervous system. This method provides a probabilistic approach in which uncertainty is accommodated by derivation, from the imaging data, of data representing a probability density function for the diffusion parameters, for example the principal diffusion direction and optionally also one or more of diffusivity, volume fraction of anisotropic diffusion or any other parameter in a local model of diffusion. This provides a considerable advantage over known techniques which find a point estimate of the diffusion parameters, using a maximum likelihood technique. The advantage arises because the imaging data in fact includes considerable uncertainty. This is due not only to noise, but significantly also to complexity of the anatomical structure. As an example of the latter, the structure may include paths which diverge, for example as paths in white matter approach grey matter, or which cross unrelated paths. Thus, many actual anatomical patterns cannot be properly represented by a point estimate such as principal diffusion direction from a diffusion tensor as in diffusion tensor imaging (DTI). As fibres approach grey matter, the diffusion anisotropy in the voxels becomes small and calculated principal diffusion direction vectors become increasingly variable due to noise in the image. This problem is so pronounced that the streamlining algorithms used to date have been forced to apply an arbitrary anisotropy threshold which forces the early termination of reconstructed fibres. For example, if a path passes through a voxel with low anisotropy, streamlining algorithms have too little information about the direction of progression of the fibre in an adjacent voxel to extend it further. This is even more of a problem when tracing fibres from grey matter, which typically has low anisotropy.

In contrast, the method in accordance with the second aspect of the present invention represents the uncertainty explicitly in a probability density function which is calculated from the imaging data. This provides a number of advantages, particularly when applied to a method of analysing structure in accordance with the first aspect of the present invention.

First, as an explicit representation of uncertainty is generated for a distribution of the possible paths, the relative likelihood of a pathway can be estimated even if the absolute probability is low. This reduces the confounds of complicated features such as crossing fibres. It also leads to a second advantage of robustness to noise. It can be difficult to track a non-probabilistic algorithm beyond a noisy voxel as it may initiate a meaningless change in path. However, with a probabilistic algorithm, paths which have taken errant routes tend to disperse quickly, so that voxels along these paths are classified with very low probability. In contrast, real paths tend to group together, giving a much higher relative probability of connection to the seed point for voxels on these paths.

It is to be noted that the second aspect of the invention involves deriving the diffusion data representing the probability density function from the imaging data. Therefore the probability density function represents the uncertainty in the diffusion parameters calculated from the imaging data itself. This provides significant advantages over the alternative possibility of associating a point estimate of the diffusion parameters with a predicted uncertainty not actually derived from the imaging data. Uncertainty in what prediction to apply leads to a lack of confidence about all such predictions. In contrast, the derived diffusion data is, by its deterministic nature, a reliable estimate of the uncertainty. Thus the derived diffusion data provides better results when used to derive connectivity patterns.

The derivation of diffusion data representing a probability density function (pdf) is preferably performed using a sampling technique, such as Markov Chain Monte Carlo sampling. However, any other suitable technique may be used. Generally, such techniques will derive a pdf for the parameters of a model of the anatomical structure. Any suitable model may be used, including a partial volume model or a diffusion tensor model (as in DTI). It may be a technique which derives a posterior pdf. Other suitable techniques include: global search; importance sampling; perfect sampling; variational Bayes; mixture sampling; Laplace approximation; ensemble learning rejector sampling; expectation maximisation; marginalised expectation maximisation; or empirical Bayes, but these techniques are merely examples.

Preferably, the derived connectivity pattern data represents the probability of connection of voxels of said target portion of said nervous system to the seed voxel. Such derivation of the probability of connection is particularly easy to calculate from the particular representation of the diffusion data as a local probability density function. For example, a probabilistic integration algorithm may be used to calculate the probability of a connected path existing.

As discussed in more detail below, the first and second aspects of the present invention have been applied together to analyse the structure of the brain for seed voxels in the thalamus as the first voxels which are classified, specifying regions of the cerebral cortex as the target region and second voxels. A method in accordance with the second aspect of the invention has been used to derive connectivity pattern data showing pathways, including major sensory and motor pathways, between the thalamus and the cortex, for example to the visual cortex and the optic tract.

Furthermore, the method in accordance with the first aspect of the present invention has been used to classify nuclear groups of the thalamus into clusters showing a strong correspondence with known locations of thalamic nuclei in humans and connections in non-human primates. On this basis, it is understood that the classified groups correspond to different thalamic nuclei or nuclear groups. The results demonstrate that it is possible to trace connectivity patterns from the thalamus to the grey matter of the human cerebral cortex in vivo using DWI. The approach is generalisable and is applicable to study other grey matter structures. Therefore, it is expected that the present invention will be applicable to analysis of the structure of any part of the nervous system, by appropriate selection of seed voxels, based on the connectivity pattern to any other regions of the nervous system, by appropriate selection of the target region.

The connectivity pattern data sets and the classified seed voxels derived by the present invention may be imaged by deriving image data for display on a display device.

According to further aspects of the present invention, there is provided a computer program, which may be encoded in a storage medium, for performing a method in accordance with the first or second aspect of the invention when executed on a computer. There is also provided a computer apparatus having means for performing the methods. The present invention provides a tool for study of the nervous system of a mammal, for both scientific and clinical research. It is expected that the non-invasive in vivo definition of connectivity patterns will complement functional imaging and has the potential to provide new insights into understanding disorders associated with developmental or regional alternations of brain connectivity. As a specific example, the present invention might be used to test the hypothesis that thalamic disfunction is a factor in schizophrenia, as well as aiding the interpretation of functional and structural differences in patients. In addition, impairments in cortico-cortical connectivity are found with learning difficulties and quantitative differences might explain variations in learning abilities and performance.

Another possible use is as a tool to plan and guide surgery. For example, specific thalamic nuclei are targets for stereotactic neurosurgery or deep brain stimulation in Parkinson's disease. At present, targeting uses stereotactic atlases (of another individual's brain) and direct electrical stimulation of potential target nuclei to localise sensor motor effects. Potentially, the present invention could be used to provide non-invasive identification of probable boundaries between major nuclei in the patient undergoing surgery and thereby improve targeting accuracy and outcomes.

An embodiment of the present invention will now be described by way of non-limitative example with reference to the accompanying drawings, in which:

Fig. 1 is a flow chart of a method in accordance with the present invention; Fig. 2 is a schematic view of a computer for implementing steps of the method of Fig. 1;

Figs. 3 A to 3D are images of a connectivity pattern data set produced using the method of Fig.1 for seed voxels in the thalamus;

Figs. 4A to 4D are images derived during a first performance of the method of Fig. 1 specifying seed voxels in the thalamus and specifying four regions of the cerebral cortex as the target portion, Fig. 4 A being an image of the regions of target portion, Fig. 4B being an image of the cerebral cortex illustrating histologically- identified regions and Figs. 4C and 4D being images of the classified seed voxels; Figs. 5 A to 5D are images derived during a second performance of the method of Fig. 1 specifying seed voxels in the thalamus and specifying seven regions of the cortex as the target portion, Fig. 5 A being an image of the regions of the target portion and Figs. 5B to D being images of the classified seed voxels; Fig. 6 is an image of classified seed voxels derived during a third performance of the method of Fig. 1; and

Figs. 7A to 7F are images are classified seed voxels derived during a fourth performance of the method of Fig. 1.

A method in accordance with the present invention will now be described with reference to Fig. 1.

In Step 1, DWI is performed to derive DWI imaging data 2. DWI is in itself a known type of MRI which is applied in a conventional manner to the present invention. For example, DWI is disclosed in: Basser, P.J., Mattiello, J. & LeBihan, D., "Estimation of the effective self-diffusion tensor from the NMR spin echo", J. Magn Reson. B 103, 247-254 (1994); Basser, P.J., Pajevic, S., Pierpaoli, C, Duda, J. & Aldroubi, A., "In vivo fiber tractography using DT-MRI data", Magn Reson. Med. 44, 625-632 (2000); Jones, D.K., Simmons, A., Williams, S.C. & Horsfield, M.A., "Non-invasive assessment of axonal fiber connectivity in the human brain via diffusion tensor MRI", Magn Reson. Med. 42, 37-41 (1999); and Conturo, T.E. et al., "Tracking neuronal fiber pathways in the living human brain", Proc. Natl. Acad. Sci. U.S.A. 96, 10422-10427 (1999).

Step 1 may be performed using conventional magnetic resonance imaging (MRI) equipment. DWI uses a specific spin echo sequence, for example as follows. Firstly, a pulse is applied to excite the spins. Then a large gradient is applied, followed by a pulse to provide a 180° flip. After that the same gradient is applied to re-phase the spins. The ratio, within each voxel, of the signal output in the presence of a diffusion gradient to the signal output in the absence of said gradient encodes information about the self diffusion properties in the direction of said gradient. Therefore the MRI is performed with and without said gradient to output the imaging data 2 as data representing said ratio for each voxel. The DWI in step 1 is performed over a large number of directions, as is conventional, typically at least 50 or 100 directions. Thus, the DWI imaging data 2 represents the measured ratio in all these directions.

The DWI is performed in Step 1 to image the subject region of the nervous system of a mammal, with particular advantage to a human. For example, the subject region may be the entire brain.

During Step 1, MRI may also be performed using the same MRI equipment to provide TI -weighted imaging data 3 for use in segmentation, skull stripping and registration of the DWI imaging data 2. The subsequent steps of the method are performed on the DWI imaging data

2 using a computer system programmed to implement the method. Thus, the subsequent steps of the method are implemented by a computer program executable on a computer system. When the code of the computer program is executed, the computer system is caused to execute the method. Thus, the programmed computer acts as an apparatus having means for performing the steps of the method. Any suitable computer system may be applied, for example a conventional PC 22 having a display 21, as illustrated schematically in Fig. 2. Optionally, the computer- implemented steps may be performed on pre-existing imaging data 2, in which case the DWI of step 1 will not have been performed contemporaneously with the remainder of the method.

In Step 3, diffusion data 5 representing a probability density function for the diffusion parameters is derived from the DWI imaging data 2. Due to factors such as noise in the diffusion measurement, subject motion and incomplete modelling, there is uncertainty associated with the diffusion parameters represented by the DWI imaging data 2. In the diffusion data 5, this uncertainty may be represented in the form of a probability density function (pdf); P(Ω|Y,M), where Ω are the model parameters, Y is the data and M is the model. This distribution is known as the full posterior pdf. The full posterior is formed using Bayes' theorem:

P(Ω I Y ) <* E(Y I Ω M)P(Ω) (1 ) and marginal posteriors on parameters of interest may be formed by integrating out uninteresting parameters from the full posterior.

R(Ω int|Y, ) = J (Ω I Y, M)d olher (2)

The diffusion data 5 is calculated for all the voxels in the subject region of the nervous system covered by the image data 2. The diffusion data 5 is estimated in Step 3 using Markov Chain Monte Carlo sampling techniques from a model of the image data 2. Any suitable model may be applied. The preferred model is a partial volume model, as follows.

The partial volume model assumes that the MR measurements are a sum of signal from within a white matter tract (which is modelled with wholly anisotropic diffusion), and signal from outside the tract (which is modelled with isotropic diffusion). Independent Gaussian noise may be assumed.

where n is the number of data collected at each voxel, y_{ is the i datum σ is the standard deviation of the Gaussian noise, and μ, is the predicted value of the z^'th datum.

" f_} exp(-b, r,^rR,AR/r,) (5)

where d is the diffusivity, b_{ and are the b-value and gradient direction associated with the i measurement respectively,/] and R_jAR_j ^T are the fraction of signal contributed by, and anisotropic diffusion tensor along, the 7^th fibre direction (θ_j.<p_j) respectively, and: A

In the work reported herein, the model has been limited to a single fibre, but alternatively a model including plural fibre directions could be used. The simple fibre model has 5 parameters: d, f_l5 θ_b φ_{ and S₀. Non-informative prior distributions on all parameters may be chosen.

In the sampling, samples are drawn from the parameters named above and the precision (1/σ²) using Monte Carlo Markov Chain sampling. This involves drawing samples of the fibre direction from this distribution at each voxel to develop a representation of the probability distribution function. The diffusion data 5 may represent the probability density function

P (θ,φ I Y). Based on the understanding that the diffusion direction conesponds to the connectivity direction of pathways in the voxel, then the diffusion data 5 may also be interpreted as the probability density function for connectivity in a particular direction (θ, φ). On this basis, Step 6 involves deriving a connectivity pattern data set 7 in respect of each of a plurality of pre-selected seed voxels representing the connectivity pattern for a pre-selected target region. Data 8 representing the seed voxels and data 9 representing the target region are input by the user specifying them using a segmentation technique. This may be performed using the TI -weighted imaging data 3 and conventional data manipulation tools. One possible segmentation technique, given by way of example, is as follows.

Firstly, using a segmentation tool, a probabilistic matter type segmentation and partial volume estimation on the TI -weighted imaging data 3 is formed. These results may optionally be thresholded to include only voxels estimated at greater than a threshold amount of grey matter, for example 35%, to mask the target regions. Using a skull-stripping tool, the DWI imaging data 2 and the TI -weighted imaging data 3 may be skull-stripped. Using a registration tool, affine registration between the subject region and the TI -weighted imaging data 3 may be performed to derive the transformation matrix between the two spaces.

If it is desired to image the connectivity pattern data set 7, then the target region may be specified as the entire subject region. Alternatively, the target region may be specified as a particular region, typically of grey matter, for which the connectivity to the seed voxels is of interest.

The connectivity pattern data set is calculated to represent the probability of connection of voxels of the target portion to the seed voxel. Thus, the connectivity pattern data set 7 may be calculated using a probabilistic integration algorithm which calculates the probability of a connecting path existing, as follows.

In principle, the connectivity pattern data set may be derived by solving the following path integral, based on the diffusion data 5.

We may define the probability of connection between two points A and B as P(3 some path through the data starting at A, ending at B). In general, the connecting path may have any length. We may now express the probability of connection with a specified path length n as a path integral:

P(3A → B\Y,n) = \ ,_A→BP(θ^ \n--P(θ_n,<t>_n\Y)dx_v..dx_n (6)

where A → B represents a possible path connecting A with B and x_; are the positions of the nodes of the connected path. P(x_{ | Xj-i) = P(Θ,Φ | Y) where JC; is the position of a particle at time i, and (θ,Φ) is the direction which takes x_Λ to J _;.

The path integral represented by equation (6) cannot be solved analytically, but may be solved numerically. The prefened technique of which is implemented in Step 6 is as follows.

The numerical solution is based on a process of simulated diffusion dependent on the local probabilities represented by the diffusion data 5 to draw independent samples from the spatial probability density function of a connective path between the seed voxel and all other voxels in the subject region.

Particles are notionally placed at the seed point in the data field and, at each time step, there is chosen a direction in which to move according to the local -Improbabilities. The spatial density function of x„|x₀ is therefore:

= P(3A → B\Y,n)

Hence the pdf of connection with path length n between points A and B is equal to the pdf of a single particle being at point B after n steps of simulated diffusion seeded at point A. The pdf of connection with any path length n ≤ m is the pdf of the particle being at point B after any number of steps of simulated diffusion up to m. The probability of connection between point A and voxel V is obtained by integrating R(3 A→B | Y),VB in voxel V.

The probabilistic nature of this algorithm automatically incorporates the estimated uncertainty in principle diffusion direction. Effectively, this process involves generating a pathway by producing one of the samples (Θ,Φ) from the distribution at a given point, moving a distance along (Θ,Φ) to a new point in the distribution. The process is repeated until the edge of the subject region is reached or, optionally, until the path begins to loop back on itself. This criterion of the path looping back on itself is implemented by use of a curvature threshold, which is preferably lenient, for example requiring that successive directions must be within 80° of each other, and by testing whether the path passes through the same area more than once. Crucially, the effect of this whole procedure is to integrate both the uncertainty in all local fibre directions to produce the best estimate of the spatial pdf for the global connectivity between the seed voxel and the target voxels. Each of the probable pathways generated represents a sample from this spatial pdf on global connectivity. Convergence of the solution is assumed after a predetermined number of samples have been drawn, for example 10,000 samples.

Optionally, in Steps 10 and 12, the connectivity pattern data set 7 is displayed using conventional display tools, for example, MEDx (Sensor Systems, Inc., VA, US); 3D Doctor MRIcro, etc. Image processing in Step 10 produces image data 11 which may be displayed in Step 12 on a display device, for example the display 21 of the computer 22. In Step 13, the connectivity pattern data set 7 in respect of each of the seed voxels is used to classify the seed voxels. In this case the seed voxels are the "first" voxels and the voxels of the target region are the "second" voxels in accordance with the first aspect of the invention. In particular, the criterion used for classification is the probability of connection of the seed voxels to pre-selected regions of the target portion. The user may input data 14 specifying the pre-selected regions. This may be achieved using the same tools as are used to enter the data 8 representing the seed voxels and the data 9 representing the target regions, as described above.

The classification may be based on the probability of connection of the seed voxels to the pre-selected regions of the target portion in various ways. For example, the seed voxels may be classified as being connected to the pre-selected region of the target portion which has the highest probability of connection. Alternatively, the seed voxels may be classified as being connected to each of the pre-selected regions of the target portion for which the probability of connection exceeds a pre-selected threshold. Indeed, many other criteria based on the connectivity pattern data set may be used, depending on the structure of interest. Such criteria will be based on the pattern of anatomical connectivity represented by the connectivity pattern data 7, for example the shape, route, destination or strength of the connectivity paths.

As an alternative in Step 13, the voxels of the target portion could be classified based on the same connectivity pattern data sets. For example, the voxels of the target portion could be classified based on the probability of connection to particular seed voxels, or on other criteria based on the pattern of anatomical connectivity. In this case, the voxels of the target region are the "first" voxels and the seed voxels are the "second" voxels.

Step 13 outputs data 15 representing the classification of the seed voxels. Optionally, in step 16, image processing may be performed to derive image data 17 representing an image of the seed voxels indicating their classification. The image processing in step 16 may be performed using conventional tools, for example those described for use in step 10 above. The resultant image data 17 may be displayed in Step 18 on a display device, for example the display 21 of the computer 22. To demonstrate the effectiveness of the method illustrated in Fig. 1, there will now be described an example of a specific imaging process which has been performed using this method. Reference is made to Figs. 3 to 6 which are images derived during this method. Figs. 3 to 6 are black-and-white representations of images which in reality are in colour. Whilst Figs. 3 to 6 are sufficient to demonstrate the effectiveness of the present invention, it should be borne in mind that the original colour images are much clearer.

DWI imaging data 2 was acquired in the DWI of step 1 by echo planar imaging on a GE 1.5 Tesla scanner with a standard quadrature head-coil (54 axes isotropically distributed with b-value of 1150 s/mm²). Six diffusion-weighted volumes were acquired with b-value 300 s/mm², and six volumes were acquired with no diffusion weighting. Each volume was acquired as 60 slices with slice thickness 2.3mm, field of view 220 x 220 mm², imaging matrix 96 x 96 zeropadded to 128 x 128, and in-slice pixel dimensions 1.2 mm². The subject region was the entire brain of a human. The data manipulation tools used were from the software library of the

Oxford Centre for Functional Magnetic Resonance Imaging of the Brain, cunently available at the URL ht ://www.fmrib.ox.ac.uk/fsl.

Steps 4 and 6 were carried out as described above, specifying all the voxels within the thalamus as seed voxels and specifying the cerebral cortex as the target region.

The resultant connectivity pattern data sets 7 were displayed in Steps 10 and 12. Figs. 3 A to 3D are images thus derived, in which voxels are colour coded according to the probability of connection to seed voxels in the thalamus. The thalamus itself is also colour coded at the centre of the respective images. As is clear from Figs. 3A to 3C, this technique provides the ability to trace pathways, including major sensory and motor pathways all the way from the thalamus to the cortex. For example, Fig. 3 A represents the connectivity pattern data set 7 for a seed voxel in the lateral geniculate nucleus (LGN) and clearly shows pathways to the visual cortex and optic tract. It can also be seen that the connectivity pattern data set 7 can define even relatively complex pathways. This is understood to be because the method does not assume any particular structure. For example, Fig. 3B is an image of the connectivity pattern data set 7 for a seed voxel within the optic tract and shows the generation of two separate branches of this pathway, one to the LGN and visual cortex and the other forming the brachium which connects to the superior colliculus. Step 13 of the method described above was performed specifying portions of the cortex as the portions of the target region. Classification of seed voxels in the thalamus was performed based on the criterion of the portion of the cortex having the highest probability of connection to the seed voxel. In a first example, large anatomically-defined regions of the cortex conesponding to known connection areas of the major thalamic nuclear groups in non-human primates were specified as the portions of the target region. Fig. 4A is an image of these portions and Fig. 4B is a diagram of the histologically-defined nuclei of the thalamus colour coded according to their major cortical connections on which the specification of the portions of the target region was based. Figs. 4C and 4D are images of the classified seed voxels of the thalamus.

It can be seen that the classification identifies commonly-connected clusters. On the basis of the strong conespondence between our connectivity-based clusters in the human thalamus, known locations of thalamic nuclei in humans, and connections in non-human primates, we understand that these correspond to different thalamic nuclei or nuclear groups. In the monkey, the mediodorsal nucleus (MD) projects primarily to prefrontal and cingulate cortices and is reciprocally connected to medial temporal regions including olfactory cortex and amygdala. We found that a large medial, dorsal region of the thalamus had highly probable connections to prefrontal and temporal regions, and suggest that this similar region includes MD. Studies in non-human primates have shown that the ventral posterior nucleus (VP) projects to primary and secondary somatosensory areas. As shown in Figs. 4C and 4D, a similar ventral posterior region was found with a strong probability of somatosensory connections. It is suggested that this conesponds to the human VP. In monkeys, the ventral lateral (VL) and ventral anterior (VA) nuclei project to primary (Ml) and premotor cortex (PMC). As shown in Figs. 4C and 4D, a lateral region, anterior to putative VP, was shown here to have a high probability of motor cortical connectivity. It is suggested that this region includes VL and VA. The lateral posterior nucleus (LP) and pulvina (Pu) project to posterior parietal cortex (PPC) and extrastriate areas. As shown in Figs. 4B and 4C, a posterior region was found that is connected to these areas. It is proposed that this region conesponds to the LP/Pu complex. Confidence in the connectivity parcellations was increased by the demonstration that the pattern of the emerging clusters was comparable between the left and right thalami. In a second example, the specified portions of the target region comprised further divisions of the cortical surface in the left hemisphere (for example, the primary motor cortex separated from the lateral and medial premotor cortices (PMC)). Fig. 5 A is an image of the portions of the target region. This was to test for sub-divisions within the major connectivity clusters, expected from cyoarchitectonic analyses. The classification resulted in further differentiation of thalamic subregions. Figs. 5B to 5D are images of the classified seed voxels.

As shown in Fig. 5B, within the large lateral clusters, smaller distinct regions connecting to somatosensory cortices, Ml, and PMC were distributed along a posterio-anterior axis. The known anatomy in the non-human primate, in which VP projects to somatosensory cortices, VLp projects to area 4, and VLa and VA project to lateral and medial area 6, suggests that these correspond to VP, VLp, and VLa/NA. On the cortical surface, the posterior parietal (PPC) and the occipital cortex were also separated. Connectivity-tracing generated further subdivision in the thalamus. As shown in Fig. 5C, it is proposed that the inferior/lateral area, with most highly probable connections to the occipital lobe, conesponds to the lateral geniculate nucleus (LGN). The posterior region that connects to PPC may include the lateral posterior (LP) nucleus and the anterior and lateral pulvinar. Medial to this, there are areas that connect to the occipital lobe and to the temporal cortex. These may include parts of the medial and inferior pulvinar. Finer structural assignments ultimately are limited by characteristics of the data, such as voxel size; the apparent connectivity of some of these more medial voxels may be influenced by the presence of the white matter tracts from MD and LGN that pass through this region en route to the visual and temporal cortices.

The probabilistic classifications of cortical connectivities demonstrate that connectivity-based boundaries between commonly-connected thalamic clusters are not sharp. In addition, the probabilistic representations demonstrate that within clusters showing a common most probable cortical target, pathways from many voxels show relatively high probabilities of connectivity to more than one cortical area. For example, Fig. 6 is an image of seed voxels classified as having a suprathreshold probability of connection to prefrontal cortex (pfc) and premotor cortex (pmc) and shows that in the cluster that we propose to include V A, some pathways reached premotor cortex and others reached prefrontal cortex consistent with the known distribution of cortical connections of this nucleus in the monkey. However, this could be due to multiple factors including partial volume effects (i.e. individual voxels will contain signal contributions from more than one nucleus), genuinely fuzzy histological borders between nuclei or between cortical target areas, limited conespondence between sulcal landmarks and cytoarchitectonic cortical boundaries, and noise in the MR signal.

The classification of thalamic voxels described thus far has been based only on probable connections to ipsilateral cerebral cortex. However, thalamic connections with subcortical structures also are generated. These provide further information regarding thalamic structure. There are a few voxels within the thalamus for which the probability of connection to cortical grey matter was very small. To better define regions that show lower probabilities of cortical connections, the criterion for classification may be based on a threshold probability. Figs 7A to 7E are images for classification based on the probability of reaching the cortex exceeded 10% (Figs. 7 A, 7B), 40% (Figs. 7C, 7D) and 80% (Figs. 7E, 7F). The classified clusters proposed to conespond to the major nuclei were preserved. As is clear from Figs. 7A to 7F, as the threshold was increased, an enlarging region emerged between the proposed lateral and medial nuclear groups that did not show suprathreshold connectivity probability to any cortical region. Pathways from this region were generated mainly to the basal ganglia or crossed the corpus callosum. It is proposed that this region includes parts of the internal medullary lamina, which has diffuse cortical connections, but projects predominantly to the striatum in the monkey.

Claims

1. A computer-implemented method of analysing the structure of the nervous - system of a mammal from imaging data representing local anatomical connectivity of a subject region of the nervous system, the method comprising: deriving from said imaging data, in respect of each of a plurality of first voxels in subject region of the nervous system, a connectivity pattern data set representing the anatomical connectivity between said first voxels and a at least one second voxel of said nervous system; and classifying the first voxels on the basis of their respective connectivity pattern data sets according to a predetermined criterion.

2. A method according to claim 1, wherein each connectivity pattern data set represents the probability of connection between the respective first voxel and the at least one second voxel.

3. A method according to claim 2, wherein said predetermined criterion is the probability of connection between the respective first voxel and a plurality of preselected regions each comprising at least one second voxel.

4. A method according to claim 3, said step of classifying comprises classifying the first voxels as being connected to the pre-selected region which has the highest probability of connection to the respective first voxels.

5. A method according to claim 3, said step of classifying comprises classifying the first voxels as being connected to each of the pre-selected regions for which the probability of connection to the respective first voxels exceeds a pre-selected threshold.

6. A method according to any one of claims 3 to 5, further comprising the step of receiving data from a user specifying the pre-selected regions.

7. A method according to claim 6, wherein said step of receiving data from a user specifying the pre-selected regions uses a segmentation technique.

8. A method according to any one of the preceding claims, wherein at least part of said target portion is positioned at a distance from said seed voxels of at least the size of said seed voxels.

9. A method according to any one of the preceding claims, further comprising deriving image data for display on a display device, representing an image of said seed voxels indicating the classification of the seed voxels.

10. A method according to any one of the preceding claims, wherein said imaging data representing local anatomical connectivity comprises diffusion- weighted magnetic resonance imaging data.

11. A method of analysing the structure of the nervous system of a mammal, comprising: performing imaging of a subject region of said nervous system of said mammal to obtain imaging data representing local anatomical connectivity; and performing a method according to any one of the preceding claims on said obtained imaging data.

12. A computer-implemented method of deriving, from diffusion- weighted magnetic resonance imaging data of a subject region of the nervous system of a mammal, a connectivity pattern data set describing the anatomical connectivity of said nervous system, the method comprising: deriving, from said diffusion-weighted magnetic resonance imaging data, diffusion data representing a probability density function for the diffusion parameters at respective voxels in said subject region of the nervous system; and deriving from said diffusion data, as said connectivity pattern data set, a data set representing connectivity of a pre-selected seed voxel in said subject region of the nervous system to a target portion of said nervous system.

13. A method according to claim 12, wherein said step of deriving diffusion data is performed using a sampling technique.

14. A method according to claim 13, wherein said sampling technique is a Markov Chain Monte Carlo sampling technique.

15. A method according to any one of claims 12 to 14, wherein said step of deriving diffusion data comprises deriving, as said diffusion data, probability density functions for parameters of a model of anatomical structure.

16. A method according to claim 15, wherein said model is a partial volume model.

17. A method according to any one of claims 12 to 16, wherein said step of deriving diffusion data is performed assuming independent Gaussian noise is present in said diffusion-weighted magnetic resonance imaging data.

18. A method according to any one of claims 12 to 17, wherein said connectivity pattern data set represents the probability of connection of voxels of said target portion of said nervous system to the seed voxel.

19. A method according to claim 18, wherein said step of deriving a data set as said connectivity pattern data set uses a probabalistic integration algorithm to calculate the probability of a connected path existing.

20. A method according to claim 19, wherein said step of deriving a data set as said connectivity pattern data set comprises solving the path integral, based on the diffusion data, for the probability of connection of the seed voxel to a voxel of said target portion integrated over possible paths.

21. A method according to claim 20, wherein said line integral is solved by sampling the diffusion data over simulated diffusion paths between the seed voxel and the voxel of the target portion.

22. A method according to any one of claims 12 to 21, wherein at least part of said target portion is positioned at a distance from said seed voxels of at least the size of said seed voxels.

23. A method according to any one of claims 12 to 22 further comprising deriving image data, for display on a display device, representing an image of said connectivity pattern data set.

24. A method according to claim 10, wherein said step of deriving a connectivity pattern data set comprises a method according to any one of claims 12 to 23.

25. A method of deriving a connectivity pattern data set describing the anatomical connectivity of the nervous system of a mammal, the method comprising: performing diffusion-weighted magnetic resonance imaging of a subject region of said nervous system of said mammal to derive diffusion- weighted magnetic resonance imaging data; and performing a method according to any one of claims 12 to 24 on said derived diffusion-weighted magnetic resonance imaging data.

26. A computer program executable on a computer system, the computer program being capable of causing the computer to execute a method according to any one of claims 1 to 10 or 12 to 24.

27. A storage medium storing a computer program according to claim 26.

28. A computer program storage medium readable by a computer system and encoding a computer program of instructions for executing a computer process for performing a method according to any one of claims 1 to 10 or 12 to 24.

29. A computer apparatus for analysing the structure of the nervous system of a mammal from diffusion-weighted magnetic resonance imaging data of a subject region of the nervous system, the apparatus comprising: means for deriving, from said diffusion-weighted magnetic resonance imaging data, in respect of each of a plurality of seed voxels in said subject region of the nervous system, a connectivity pattern data set representing the anatomical connectivity of said seed voxels to a target portion of said nervous system; and means for classifying the seed voxels on the basis of their respective connectivity pattern data sets according to a predetermined criterion.

30. Computer apparatus for deriving, from diffusion-weighted magnetic resonance imaging data of a subject region of the nervous system of a mammal, a connectivity pattern data set describing the anatomical connectivity of said nervous system, the apparatus comprising: means for deriving, from said diffusion-weighted magnetic resonance imaging data, diffusion data representing a probability density function for the diffusion parameters at respective voxels in said subject region of the nervous system, means for deriving from said diffusion data, as said connectivity pattern data set, a data set representing connectivity of a pre-selected seed voxel in said subject region of the nervous system to a target portion of said nervous system.