Systems and methods for improved tractography images

This application is a Broadening Reissue of U.S. Pat. No. 10,702,156 (previously U.S. patent application Ser. No. 15/291,959) titled “SYSTEMS AND METHODS FOR IMPROVED TRACTOGRAPHY IMAGES” issued on Jul. 7, 2020, which is incorporated herein by reference in its entirety.

Multiple sclerosis (MS) is a disease of the central nervous system. MS is an inflammatory, demyelinating disease that affects more than 2 million people worldwide. MS can primarily affect the white matter (WM) and grey matter in the brain and the spinal cord. MS can result in neuronal and axonal degeneration, which can be observed as brain lesions or plaques, grey matter atrophy, and diffuse abnormalities. MS manifest itself in a number of different forms, including: clinically isolated syndrome, relapsing-remitting syndrome, secondary progressive syndrome, primary progressive syndrome, and progressive relapsing syndrome. Damage caused to the brain by MS or other neurodegenerative diseases can make it difficult to perform tractography to determine brain connections, which can limit the usefulness of tractography in clinical evaluation of patients with neurodegenerative diseases.

The following description of the drawings and detailed description are exemplary and explanatory and are intended to provide further explanation of the invention as claimed. Other objects, advantages, and novel features will be readily apparent to those skilled in the art from the following brief description of the drawings and detailed description of the invention.

The various concepts introduced above and discussed in greater detail below may be implemented in any of numerous ways, as the described concepts are not limited to any particular manner of implementation. Examples of specific implementations and applications are provided primarily for illustrative purposes.

The present disclosure discusses systems and methods for identifying biomarkers that can help with the diagnosis, prognosis, and treatment choices of patients with brain diseases, including demyelinating or neurodegenerative diseases, stroke or brain trauma. Patients can be evaluated using clinical scales, but the scales can be based on patient questionnaires making them an unreliable measure between patients and over the course of the patient's disease progression. Diffusion based magnetic resonance imaging (dMRI) techniques often use parameters such as fractional anisotropy (FA), mean diffusivity (MD), and radial diffusivity (RD) that are based on simple models that can fail in the presence of neurodegeneration, such as demyelination. The present disclosure discusses systems and methods that enhance dMRI images and enable tractography to be performed on images of a damaged nervous system. The damaged tracks identified by the present system can be used as a biomarker for the assessment of patients. In some implementations, the biomarkers are converted into clinical scales that can be used to compare patients to one another or a patient over time.

illustrates a block diagram of an example systemfor the assessment of neurological damage. The systemincludes an imaging system. The imaging systemprovides imaging data to the assessment engine. The assessment engineincludes a databasefor the storage of anatomical image recordsand diffusion image records. The assessment engineincludes a segmentation engine, an enhancement engine, a tractography engine, and a scoring engine. In some implementations, the assessment enginecan be coupled to a monitor or other system for the display of the results generated by the assessment engine.

The systemincludes the imaging systemthat provides imaging data to the assessment engine. The imaging systemcan be one or more magnetic resonance imaging (MRI) systems. The imaging systemcan be configured to acquire imaging data using different imaging acquisition modalities. The imaging systemcan be configured to capture and generate both anatomical image records and diffusion image records. For example, the imaging systemcan acquire T1, T2, high-angular resolution diffusion images (HARDI), function MRI (fMRI), magnetization-prepared rapid gradient-echo (MPRAGE), fluid-attenuated inversion recovery (FLAIR), diffusion tensor imaging (DTI), diffusion spectrum imaging (DSI), and optical coherence tomography (OCT), spectroscopy or any combination thereof. In other implementations, a first imaging systemcan capture and generate the anatomical image records and a second imaging systemcan capture and generate the diffusion image records. In some implementations, the imaging systemprovides the imaging data directly to the assessment enginethrough a direct data or network connection. In other implementations, the imaging systemcan provide the imaging data to the assessment enginethrough an intermediary device. For example, the imaging systemcan first provide the imaging data to an intermediary device such as a networked server, cloud based storage, or other computer, and the assessment enginecan retrieve the imaging data from the intermediary device prior to the analysis of the imaging data by the assessment engine.

The systemalso includes the assessment engine. The assessment engineand its components are described in greater detail below. As an overview, the assessment enginecan receive neurological imaging data (in the form of anatomical image records and diffusion image records), enhance the imaging data, and generate a score to quantify the level of neurological damage experienced by the patient. For example, and using MS as an example disease, the assessment enginecan generate a score that relates to the amount of white matter tracts affected by MS lesions. The systemcan also be used with patient suffering (or believe to be suffering) from Alzheimer's disease, Parkinson's disease, dementia, any other brain disease that presents in the MM with lesional change of gray scale intensity in the corresponding Mill image. In some implementations, the level of damage is converted to a clinical score, such as the Expanded Disability Status Scale (EDSS), Sloan score, multiple sclerosis functional composite (MSFC), paced auditory serial addition test (PASAT), brief repeatable battery-Neuropsychology (BRB-N) test, selective reminder test (SRT), symbol digit modality test (SDMT), spatial recall test, and the word list generation (WLG) test. The lesions change the diffusion properties of the tissue that affect the fiber tracking parameters and result in incorrect streamline reconstruction compared to the underlying anatomy. For example, when a fiber enters a brain lesion the traditional techniques have difficulty in determining whether the fiber continues through the lesion or stops. In some implementations, the techniques have difficulty because the values of the tracking are abnormal in the lesion area the fiber can continue in a random direction introducing an error that will propagate through the fiber tracking algorithm. Because these techniques cannot determine whether the track continues through or stops in the lesion, the techniques cannot accurately calculate fiber track damage caused by the lesion. The assessment enginecan improve the analysis of the areas containing lesions and plaques to generate improved tractography through the damaged areas so that the impairment can be assessed.

In some implementations, the assessment engineis a stand-alone device and in other implementations the assessment engineis a component of another device. As a stand-alone device, the assessment enginecan include special purpose logic circuitry (e.g., a field programmable gate array (FPGA) or an application specific integrated circuit (ASIC)), a microprocessor, or a combination thereof. The assessment enginecan be coupled with a computer or imaging systemvia a wired or wireless network connection or other wired or wireless connections, such as, but not limited to, a universal serial bus (USB) connection, FireWire connection, eSATA connection, or Thunderbolt connection. When provided with imaging data, the stand-alone assessment enginecan return a score for the patient or other result as described herein. In other implementations, the assessment enginecan be implemented as a component of the imaging systemor other system, such as a desktop computer, and one or more components of the assessment enginecan be implemented as components of the other system.

The assessment engineincludes the database. The databasecan be stored on a computer readable medium such as, but not limited to, a magnetic disk hard drive, random-access memory (RAM), electrically-erasable ROM (EEPROM), erasable-programmable ROM (EPROM), flash memory, optical media, or any other suitable medium for storing the anatomical image records, the diffusion image records, and processor executable instructions.

The assessment enginestores anatomical image records, diffusion image records, and other types of MRI images within the database. The anatomical image recordscan include T1 and T2 images from the imaging system. The diffusion image recordscan include HARDI and other diffusion-weighted images from the imaging system.

The assessment enginealso includes the segmentation engine. As an overview, the segmentation enginesegments the anatomical image records. In some implementations, the segmentation enginecan segment different portions of the anatomical image records. For example, the segmentation enginecan identify portions of the anatomical image recordthat correspond to the white matter, the grey matter, the skull, other anatomical structures, or any combination thereof. In some implementations, the segmentation enginecan also identify lesions within the anatomical image recordsand generate segments that contain the lesions. In some implementations, the segmentation enginesegments the anatomical image recordthrough the use of edge detection algorithms that identify boundary regions between different tissue types imaged in the anatomical image record. The segmentation enginecan also use threshold-based methods, histogram-based methods, or graph-based methods to segment the data.illustrates an example output of the segmentation engine. In some implementations, the segmentation engineidentifies and segments the brain lesions in a fully automated manner. In other implementations, the segmentation enginecan identify and segment the brain lesions in a semi-automated manner. For example, a user may identify the center (or other portion) of a brain lesion and then the segmentation enginecan identify the boundary of the lesion. In some implementations, the segmentation enginecan store the boundary coordinates of the lesion or an identification of the voxels corresponding to the lesion as a file in the database. In some implementations, the coordinates of multiple lesions are stored in separate files or as separate vectors within a single file.

illustrates an input anatomical imageand a segmented anatomical image. The anatomical imageis a T2 image that includes hyperintensities corresponding to MS lesions. The anatomical imageis the output of the segmentation engine. The segmentation engineidentified the MS lesions and generated segments(a)-(f), which each include one of the MS lesions. As described above, in some implementations, the output of the segmentation engineis an array of vectors that define the boarder of the segments(a)-(f) or an indication of the voxels contained in the segments(a)-(f).

Referring to, the assessment enginealso includes the enhancement engine. As an overview, the enhancement engineenhances the diffusion image recordsuch that the tractography enginecan track or estimate the correct fiber architecture even in the presence of lesions and other neurological damage. The enhancement engineretrieves one or more of the diffusion image recordsfrom the database. In some implementations, the diffusion image recordsare initially modeled as fiber orientation distribution functions (FOF) or spherical harmonics (SH) glyphs. The enhancement engineconverts the glyphs into an amplitude image by sampling the glyphs from a plurality of directions to generate a voxel from each glyph that defines the intensity of the diffusion along the plurality of directions. In some implementations, the glyphs are sampled along between about 100 and about 500, about 100 and about 400, about 100 and about 300, or about 100 and about 200 directions.

In some implementations, the enhancement enginemasks the amplitude image with a white matter mask. In the below described convolution step, the white matter mask can prevent the ventricles and the exterior of the brain from being included in the convolution.

The enhancement engineis also configured to generate a kernel. The enhancement enginecan convolve the masked, amplitude image with the generated kernel. The kernel can solve the diffusion equation for diffusion MRI images. The kernel represents the Brownian motion kernel on the coupled space of positions and orientations:

provides that the total integral over positions and orientations is 1. The 2D kernel is given by:

Where the short notation is:

In some implementations, to avoid numerical errors:

The diffusion parameters D, D, and the stopping time t allow the adaptation of the kernels to different patients. In some implementations, t determines the relevance of the neighborhood; Ddetermines the kernel width; and the quotient D/Dmodels the bending of the fibers along which diffusion takes place.

The enhancement enginecan then convolve the kernel with the masked, amplitude image (U):

Where pis the kernel at position y and orientation n; Δy′ is the discrete volume measure; Δn′ is the discrete surface measure; P is the set of lattice position neighboring y; and T is the set of tessellation vectors. The rotated and translated correlation kernel is:p(y′,n′)=p(R(y′−y),Rn)

Where Ris any rotation mapping onto n′.

In some implementations, the enhancement enginetunes the kernel to the patient or to a class of patient (e.g., to a type of manifestation of the disease or a specific disease), and in other implementations, the enhancement engineuses the same kernel parameters for each patient. The enhancement enginemay tune the kernel by setting the parameters D, D, t, or any combination thereof. For example, Dcan be between about 1 and about 2, Dcan be between about 0.001 and about 0.05, and t can be between about 1 and about 20.

In some implementations, the enhancement engineis configured to calculate voxel scattering coefficient (VSC). In some implementations, the VSC is calculated before the above-described kernel is applied and the kernel can be updated responsive to the VSC. The VSC can be an estimate of the number and distribution (or clustering) of the patient's lesions. For example, if the VSC is high, the voxels are very spread (maybe forming small and separated lesion groups). If the VSC is low, the voxels are concentrated, suggesting a relatively larger group of lesions or a number of clustered small lesion. The VSC can be calculated by calculating a matrix of covariances in the x, y, and z direction among the positions of all the voxels. The enhancement enginecan then calculate the determinate of the matrix. The VSC can be the summer of all the voxels divided by the determinant.

illustrate a tractography image and an ODF glyph image of crossing phantom nerve bundles with and without, processing with the enhancement engine, respectively.illustrates a tractography imagewith a first nerve bundle(a) and a second nerve bundle(b) crossing at about 60 degrees. The corresponding glyph imageillustrates the glyphspresent at the intersectionof the first nerve bundle(a) and second nerve bundle(b) in the tractography image. As described below, the tractography engineidentifies tracks responsive to the glyphs. In the tractography image(which was not processed by the enhancement engine), the tractography engineidentifies a plurality of tracks. As illustrated in the tractography image, a bundle of fibers(e.g. a track) starts in the lower right-hand corner of the tractography imageand projects towards the upper left-hand corner of the tractography image. The crossing of the first nerve bundle(a) and the second nerve bundle(b) create a conflictive point at the intersection, and, as illustrated, a number of the tracksincorrectly terminates at the intersection.

illustrates a tractography imageof the first nerve bundle(a) and the second nerve bundle(b) crossing at about 60 degrees after processing with the enhancement engine. The corresponding glyph imageillustrates the glyphspresent at the intersectionof the first nerve bundle(a) and second nerve bundle(b) in the tractography image. As illustrated in the tractography image, a bundle of fibersstarts in the lower right-hand corner of the tractography imageand projects towards the upper left-hand corner of the tractography image. In contrast to the tractography image, the trackscontinue through the intersectionand continue projecting toward the upper left-hand corner of the tractography image.

Referring to, the assessment enginealso includes a tractography engine. As discussed briefly above, the tractography engineanalyzes the enhanced image output by the enhancement engineand estimates neural tracks in the enhanced image. The tractography enginemodels tracks (e.g., bundles of axons) that connect different regions of the brain. The tractography enginemodels the tracks based on the diffusion of water molecules within the brain and the enhanced image to generate a mask estimating where fibers would be if there were no lesions or damage. In some implementations, the tractography engineuses streamline tractography to represent the underlying neural fibers. In general, the tractography enginecan use fiber-orientation descriptive models to find paths of minimal hindrance to water diffusion using local voxel-wise orientation information. The tractography enginecan model the tracks by calculating an orientation estimate at each voxel within the enhanced image. In some implementations, the tractography enginecan generate tracks responsive to seed locations. The seed locations can serve as a beginning location of interest, and the tractography enginecan identify tracks that begin at or near the seed locations. The seed locations can be provided to the tractography engineby a user of the systemor the seed locations can be automatically determined. For example, to determine the tracks of the optic radiations, the seed location can be the thalamus, which can be identified using a gray and white matter mask. In some implementations, inclusion locations can also be used by the tractography engineto refine the generated tracks. For example, tracks starting from the seed location that do not pass through inclusion locations can be discarded by the tractography engine. Using the visual cortex an example inclusion location and continuing the above example of determining the tracks of the optic radiations, tracks that don't pass through the visual cortex can be discarded. The inclusion locations can be areas along the path of the tracks or a termination location of the tracks.

In some implementations, the tractography enginecan filter identified fibers and remove spurious fibers from the tracks. The tractography enginecan remove fibers from the tracks based on fiber length. For example, fibers that are shorter or longer than a predefined range may be excluded from the track. In another example, fibers that deviate in length from the track's average length a predetermined amount can be excluded.

In another example, the tractography enginecan filter the identified fibers using a centroid. The centroid of each fiber can be calculated to determine a position (x, y, z) for the fiber. Each of the coordinates can be compared with the distribution of all coordinates in the bundle. For example, the x position of a fiber centroid is compared with the x position distributions for all the centroids. A thresholds Tcan be defined for the x coordinate. If x differs more than Tstandard deviations from the mean x, the fiber can be discarded. The process can then be repeated for they and z coordinates with their respective thresholds Tand T.

In another example, the tractography enginecan filter the identified fibers using coherence. The tractography enginecan estimate each fibers coherence by computing the increments in x, y, and z between each pair of adjacent point: incX, incY, incZ along the fiber's path. The tractography enginecan then compare the distribution of incX with the distribution of incX for all the fibers. The two distributions can be compared with a Kolmogorov Smirnov test to obtain a p-value. If the p-value is lower than a given threshold, the fiber can be discarded. The process can be repeated for they and z coordinates.

In some implementations, one or more of the filtering methods can be used in combination with one another. For example, first the fibers can be filtered using the coherence method. The remaining fibers can be further filtered using the fiber length method, and the remaining fibers can be further filtered using the centroid method. In some implementations, the filtering process can be repeated a predetermined number of times.

illustrates an example tractography imagegenerated by the tractography engine. As described above, the tractography enginedetermines or is provided with a seed location from which the tracks are to start. To identify the tracks starting in the thalamus and ending in the visual cortex, in the tractography image, the thalamusis selected as the seed location. In order to increase segmentation, the region of interest can be divided into multiple portions (e.g., into 4 parts). The center of mass can be estimated for each and the exterior-posterior portion is preserved where the seeding is to be done. The minimum cube that fits both seed and include region can be found and then all the cortical regions except the regions of interest are subtracted. The generated mask can be set to 1 outside the cube and inside the other cortical regions so that all fibers trespassing that mask will be excluded. The visual cortexis selected as an inclusion location. Responsive to the identification of the seed location and the inclusion location, the tractography engineidentifies the trackstraveling from the thalamus to the visual cortex.

Referring to, the assessment enginealso includes the scoring engine. The scoring enginecan generate one or more clinical scores for the images. To determine the amount of track damage the scoring enginecan register the segmented anatomical image with the tractography image. For the volumetric scores, the tractography is converted into a mask and registered to the anatomical image. The score can correspond to the volume or number of damaged tracks. The volumetric score can correspond to the percent of volume or the volume in the anatomical space. For the fiber-based scores, the lesions mask is registered to the diffusion space and then the score is computed. The score can correspond to the volume or number of damaged fibers (or tracks).

In some implementations, the lesion mask is registered with the tractography image by finding the lesions in the anatomical space (T1 or T2). The reference image (T1 or T2) is registered into the respective subject image (T1 or T2), and the same transformation is applied to the lesion mask so that now the lesion mask is in the same space as the subject. In some implementations, the tractography mask image can be registered into the lesion mask (in T1 or T2 space) in order to have more precision due to T1 space higher resolution.

Once the images are registered, the damage to the tracks can be converted to damage scores and clinical scores. In some implementations, damaged tracks are those that pass through one of the lesion segments. The damage score can be the percentage of tracks in a bundle that are damaged. For example, if 20% of the tracks from the thalamus of the visual cortex pass through a lesion then the score can be 20%. In another implementation, the score can be based on the volume of a lesion inside a bundle. The volume can be determined by the number of lesion voxels within the bundle divided by the volume of the bundle. In some implementations, a total score can be generated for the patient by combining the scores from different pathways. For example, the score from the corpus callosum, optical radiation, and corticospinal tract can be combined to generate a single score. In some implementations, the score assessing the level of damage can be converted into a clinical score, such as a Sloan score, EDSS score, or other above described clinical scores.

illustrates a block diagram of an example methodfor calculating a clinical score or a damage score (or more simply, a score). The methodincludes receiving anatomical and diffusion-weighted (DW) data (step). The methodalso includes generating one or more segments from the anatomical image data (step). The methodincludes enhancing the DWI data (step), and generating a tractography (step). The generated tractography is registered with the anatomical image data (step). A score is calculated responsive to the registered tractography and anatomical image data (step).

As set forth above, the methodincludes receiving anatomical image data and DWI data (step). Also referring to, the anatomical image data and the DWI data can be received by the assessment enginefrom an imaging systemor from an intermediary device, such as networked or cloud based storage. In some implementations, the received anatomical data can include 3D structural T1-Magnetization-Prepared-Rapid-Gradient-Echo (MPRAGE) data with a voxel size of about 0.9×0.9×0.9 mm. The anatomical data can also include 3D Structural Fluid Attenuated Inversion Recovery (FLAIR) data with a voxel size of about 0.9×0.9×0.9 mm. The DWI data can include High-Angular Resolution Diffusion Imaging (HARDI), with a voxel size of about 2×2×2 mm, b-value of 1500 s/mm, and about 70 gradients. The resolution of the data received by the assessment enginecan be of greater or less than the above described resolution. The received image data can be stored by the assessment engineinto the database.

The methodcan also include generating one or more segments from the anatomical data (step). In some implementations, each of the segments can include a lesion. The segment can identify the two-dimensional or three-dimensional area, boarder location, or volume within the anatomical images that is occupied by lesions or plaques. In some implementations, the assessment engineis configured to display the anatomical images to a user. The user can interact with the anatomical images and outline the lesions or plaques to generate the boarders of the segments. In other implementations, the segments are automatically generated using, for example, the FreeSurfer© software (made available by FreeSurfer, of Cambridge, Mass.). The location of each of the segments can be stored into the databasein association with the anatomical image. For example, an identifier of the volumes contained within each of the segments can be stored. In other implementations, the segmentation process is semi-automated and the user identifies the lesion (or a lesion of interest) and the assessment engine can then identify the boarder of the lesion.

The methodcan also include enhancing the DWI data (step). As described above in relation toand, in some implementations, unenhanced DWI data in patients with plaques can result in tractography that incorrectly terminates at intersections or other conflict points such as lesions. The enhancement engineof the assessment enginecan enhance the DWI data such that the tractography enginecan properly generate tracks that pass through intersections, lesions, and other conflict points. In some implementations, the enhancement engineenhances the DW data by convolving the DW data with the above described kernel. The enhancement enginecan tune the kernel to the patient's data prior to convolving the DW data with the kernel. For example, the enhancement enginecan tune the D, D, and t parameters of the kernel. In some implementations, different regions of the DW data can be convolved with different kernel settings. In some implementations, Dis between about 0.7 and about 5, between 0.8 and about 3, or between about 0.9 and about 1.2; Dis between about 0.02 and about 0.08, between about 0.03 and about 0.07, or between about 0.04 and about 0.06; and t is between about 1 and about 20, between about 1.2 and about 10, between about 1.3 and about 5, or between about 1.3 and about 1.6.

The methodalso includes generating a tractography (step). Also referring to, the tractography enginereceives the enhanced image from the enhancement engine. In some implementations, the tractography enginealso receives one or more seed locations and one or more inclusion locations. Given the seed locations and the inclusion locations, the tractography enginecan determine which tracks begin in the seed locations and then pass through or terminate in the inclusion locations. In other implementations, the tractography enginecan generate the tractography based on a template. The template can be generated by the tractography enginebased on tractographies from healthy individuals (e.g., individuals that do not have brain lesions or other damage). In some implementations, the tractography enginecan generate the template tractography by selecting and merging together tractographies form health individuals that are similar to the patient in race, weight, height, sex, age, or any combination thereof.

The methodalso includes registering the tractography with the anatomical image (step). In order to compare the tractography results with anatomical information (eg. the lesion mask) the tractography and anatomical image are overlapped into the same space. In some implementations, the anatomical image and segments are registered with the above-described template and in other implementations the anatomical image is registered with the tractography made for the patient as described in relation to step. In some implementations, the anatomical image is registered to the enhanced DWI using non-linear registration.

In some implementations, to maintain the tractography, the anatomical image can be moved to the same space as the tractography (e.g., the b0 space). In some implementations, the anatomical image can be a T1 image. The registration can include moving a specific anatomical image in T1 space (e.g., the skull-stripped T1 image itself) to an image in the b0 space. The methodcan include generating a transformation matrix that is employed to move the desired anatomical image (e.g., lesion mask) to the b0 space. Once this is done tractography and lesions are in the same space and can be compared. In some implementations, this method can be used in cases to evaluate the damage as a percentage of fibers that traverse a lesion because the spatial position of each fiber in the tractography is needed. In some implementations, a predetermined amount of the fiber must be damaged before the fiber is classified as damaged.

In another implementation, the tractography is converted into an image for registration with the lesion maps. In these implementations, the resulting image can be to the T1 space. The image can be generated by creating an empty (all zeros) 3D image. The 3D image can have the same dimensions as the image from which the tractography came from. The tractography is composed by a series of fibers, each composed by an ordered series of points (e.g., a x, y, and z coordinate). For each fiber and each point in the fiber, a “1” is added to the voxel value corresponding to that space location in the newly generated 3D image. The process is repeated for all of the points in all of the fibers. Once repeated for each fiber, each voxel has a value that indicates the number of fibers that have pass through the voxel in the tractography. The resulting image is in the b0 space and can be registered to other spaces

The methodalso includes calculating a score responsive to the registered tractography and anatomical image data (step). The clinical score can be responsive to a damage score calculated by the scoring engine, which can be converted into a clinical score using a linear correlation or a neural network. In some implementations, the scoring engine can calculate the damage score as the number of lesioned voxels within a brain track (e.g., the number of voxels within the lesion segments) divided by the total number of voxels within the brain track. In another implementation, the damage score can be calculated as the damage density divided by the bundle density. The damage density can be the sum of the voxels in the registered tractography that are located within one of the voxel segments and the bundle density can be the total sum of the voxels in the tractography.