A novel magnetic resonance imaging postprocessing technique for the assessment of intervertebral disc degeneration-Correlation with histological grading in a rabbit disc degeneration model.

INTRODUCTION: Estimation of intervertebral disc degeneration on magnetic resonance imaging (MRI) is challenging. Qualitative schemes used in clinical practice correlate poorly with pain and quantitative techniques have not entered widespread clinical use.
METHODS: As part of a prior study, 25 New Zealand white rabbits underwent annular puncture to induce disc degeneration in 50 noncontiguous lumbar discs. At 16 weeks, the animals underwent multi-echo T2 MRI scanning and were euthanized. The discs were stained and examined histologically. Quantitative T2 relaxation maps were prepared using the nonlinear least squares method. Decay Variance maps were created using a novel technique of aggregating the deviation in the intensity of each echo signal from the expected intensity based on the previous rate of decay.
RESULTS: Decay Variance maps showed a clear and well demarcated nucleus pulposus with a consistent rate of decay (low Decay Variance) in healthy discs that showed progressively more variable decay (higher Decay Variance) with increasing degeneration. Decay Variance maps required significantly less time to generate (1.0 ± 0.0 second) compared with traditional T2 relaxometry maps (5 (±0.9) to 1788.9 (±116) seconds). Histology scores correlated strongly with Decay Variance scores (r = 0.82, P < .01) and weakly with T2 signal intensity (r = 0.32, P < .01) and quantitative T2 relaxometry (r = 0.39, P < .01). Decay Variance had superior sensitivity and specificity for the detection of degenerate discs when compared to T2 signal intensity or Quantitative T2 mapping.
CONCLUSION: Our results show that using a multi-echo T2 MRI sequence, Decay Variance can quantitatively assess disc degeneration more accurately and with less image-processing time than quantitative T2 relaxometry in a rabbit disc puncture model. The technique is a viable candidate for quantitative assessment of disc degeneration on MRI scans. Further validation on human subjects is needed.

INTRODUCTION

Low back pain is a leading cause of years lived with disability (YLD) in both men and women worldwide, contributing approximately 57.6 million years to the total YLDs in 2016.1 In the United States of America, back pain is estimated to cost more than U$100 billion per year.2 The prevalence of chronic disabling back pain has tripled in 20 years by some measures.3

A major cause of low back pain is the degeneration of the Intervertebral Disc (IVD).4 The IVD consists of three parts in one functional unit: a gelatinous center called the nucleus pulposus (NP), an outer fibrocartilaginous ring called the annulus fibrosus (AF), and top and bottom caps called the cartilaginous endplates (EPs).5 Degeneration of the IVD is a complex pathological process characterized by a loss of hydration, changes in cell populations, a decrease in protein and polysaccharide concentrations and altered biomechanical properties.6

A number of new biological therapies are under active investigation in IVD degeneration including Growth and Differentiation Factor (GDF)‐5, GDF‐6, platelet‐rich plasma, mesenchymal stem cells, tocilizumab, and telomerase gene therapy.7, 8, 9, 10 The gold standard for assessing both degeneration and regeneration of the IVD is histology which requires the animal or sample to be sacrificed. The development of well‐validated nondestructive techniques with acceptable accuracy in estimating overall health of the IVD is essential in assessing longitudinal changes in the IVD in routine clinical practice, and also in evaluating the efficacy of various biological therapies aimed at regenerating the IVD. Novel postprocessing algorithms for magnetic resonance imaging (MRI) data are a promising option for nondestructive and objective assessment of IVD health.

The results from such MRI postprocessing algorithms must correlate well with the underlying disease processes and ideally be computationally efficient and practical to implement in routine clinical practice.

The extent of IVD degeneration, when visualized on a clinical MRI scan, is currently described using qualitative grades, especially those by Pfirrmann11 and Thompson.12 The field of quantitative MRI in the research of degenerative disc disease is advancing rapidly, including work using T2* and T1‐rho.13, 14

Recent contributions by Ellingson et al showed the predictive value of T2* quantitative relaxometry for sulfated glycosaminoglycans (S‐GAG) content within intervertebral discs.13 S‐GAG content has been demonstrated to correlate with IVD degeneration and likely has a physiological role in the normal function of the IVD in maintaining osmotic pressure.15

Paul et al demonstrated a strong correlation between T1‐rho and IVD degeneration in a cadaveric caprine model. The IVDs were kept in culture for 3 weeks and exposed to different concentrations of chondroitinase for inducing degeneration.14 T1‐rho is a difficult technique in that the estimation of molecular concentrations are susceptible to changes in temperature and pH, and this poses a significant challenge in its clinical implementation.16, 17 The clinical use of quantitative imaging techniques on MRI data requires significant time, not only to acquire images but to also perform postprocessing on the images; and therefore there are monetary barriers to its adoption. A review by The Association of University Radiologists' Radiology Research Alliance Quantitative Imaging Task Force found that there was “some antipathy among radiologists towards quantitative imaging” and that “with increasing emphasis on productivity, time spent performing quantitative imaging may be considered financially unrewarding”.18

Many authors in the field of intervertebral disc degeneration have also acknowledged the need for improved nondestructive imaging techniques for the assessment of animal and preclinical models of disc degeneration.19, 20

We hypothesized that the variability in the rate of signal decay over a multi‐echo MRI sequence encodes information about the tissue states that could help distinguish between a healthy and a degenerate IVD. We aimed to evaluate the difference between measured and expected signal decay for each echo in a multi‐echo T2 MRI sequence, and correlate the results with histologically graded IVD degeneration in a rabbit IVD model.

MATERIALS AND METHOD

Animals, MRI, and histology

Midline sagittal T2 weighted MRI data of 75 lumbar IVDs in 25 female New Zealand White Rabbits were collected as part of an earlier study.10 Degeneration was induced using 18 gauge needle‐puncture in 50 IVDs, and the remaining 25 served as nonpunctured controls. The degenerated IVDs were injected with either 10 μL of phosphate buffered saline or 1, 10, or 100 μg of GDF‐6.

Sixteen weeks following the needle‐puncture (12 weeks after injection), the rabbits were euthanized, and sagittal multi‐echo spin‐echo T2 MRI was performed using a 7‐Tesla BioSpec 70/30, (BRUKER, Billerica, Massachusetts). Unfixed samples were kept in a temperature controlled environment at 4°C, not fixed in formalin, and were MRI scanned within 5 hours of euthanasia. Data were collected in DICOM (Digital Imaging and Communications in Medicine) file format. Full DICOM header information containing acquisition parameters for a representative scan is provided in Supplementary Information‐1. Thereafter, the rabbit IVDs were fixed and prepared for H + E staining, and the degree of disc degeneration was scored by a pathologist as described by Miyazaki.10

Postprocessing

All postprocessing was performed using Matlab (v. R2017b; Mathworks Inc., Natick, Massachusetts). MRI data were analyzed in three ways: firstly, the raw T2 weighted signal intensity data at the first echo time (13 ms) were examined; thereafter, quantitative T2 maps were created using the pixel‐wise curve fitting operation and the nonlinear least squares method; and finally, Decay Variance maps were generated using the following method:

The Decay Variance of each pixel was calculated using the following definitions:

Signal Intensity (SI) recorded at the ith echo time i = 1,2,..,n

Signal Retention ratio (SR) at the (i + 1)th echo time i = 1,2,..,n−1

Signal Decay Change (SDC) at the (i + 2)th echo time i = 1,2,..,n−2

Thus the formula for the Decay Variance method for n acquisitions can be represented as:

Correction for the signal to noise ratio (which was taken to be inversely proportional to signal intensity at the first echo) was performed by dividing the Decay Variance by initial signal intensity on a pixel‐wise basis.

All three maps (raw T2, quantitative T2, Decay Variance) were standardized by taking the 85th percentile value as unity, and redistributing all pixel values between zero and one. Pixel values above the 85th percentile were deemed to have a value of unity. All calculations were performed in Matlab, with 32 GB of accessible RAM and 8 Intel Core i7‐7700K processor cores at 4.2 GHz. The processor time taken to generate each map for each technique was recorded using Matlab's inbuilt analytical tools.

A region of interest (ROI) was selected including the entire NP and AF but excluding the EPs for each of the 75 IVDs as shown in Figure 1.

Figure 1

Representative regions of interest used for the calculation of average quantitative T2 times, T2 weighted signal intensity and Decay Variance for each disc; in a healthy (left) and degenerate (right) intervertebral disc in a rabbit lumbar spine

For validation purposes, ROIs were created twice by two observers at least 1 week apart.

Statistical analyses were performed in IBM SPSS software (v. 20; IBM, Armonk, New York), with the level of significance at 5% (α = .05). The primary objective was to investigate the extent to which each technique correlated with the histological grading for IVD degeneration. The distribution of IVD degeneration scores was both ordinal categorical and non‐normally distributed, and therefore nonparametric Spearman's rho test was utilized.

For each of the three techniques, sensitivities and specificities for the presence of IVD degeneration using predefined thresholds of 90% sensitivity (for specificity calculation) and 90% specificity (for sensitivity calculations) were calculated, and receiver operating characteristic (ROC) curves were generated.

The research was performed with ethical approval from the University of California San Diego IACUC. Approval number S08258.

RESULTS

The Decay Variance technique generated images that showed well‐defined dark EP regions against a mid‐gray vertebral body, and well‐demarcated dark NP regions against a bright white AF in the healthy IVDs. With progressive IVD degeneration, the central dark NP became increasingly homogeneous with the AF (Figure 2).

Figure 2

Representative magnetic resonance imaging (MRI) of the rabbit lumbar spine. MRI scans were acquired 16 weeks following annular‐puncture in the intervertebral disc (IVD), and 12 weeks postinjection of phosphate buffered saline or Growth and Differentiation Factor‐6. T2 weighted image (A), inverse of T2 weighted image (B), quantitative T2 relaxometry map (C), and Decay Variance map (D) for two IVDs in one rabbit are shown. The lower IVD (white arrow) had the best possible histology score (score 4) indicating no degeneration, and the upper disc (black arrow) had the worst possible histology score (12). Despite minor difference on the T2 images (A), and modest differences in quantitative T2 relaxation times (C), on the Decay Variance map (D) the well‐demarcated dark nucleus pulposus seen in histologically healthy discs is absent in the upper disc. All windows were set programmatically between the 0th and 70th percentile of pixel values for that slice for each technique. Scale bar is 2 mm; color bar shows relaxation time for the T2 relaxometry map (C) in milliseconds

The average scores within the segmented ROI using the three MRI postprocessing techniques for different IVD grades are given in Table 1.

Table 1

Histology scores vs average scores from the three MRI techniques: 4 weeks following IVD needle‐puncture in rabbits, the injured IVDs were either treated with phosphate buffered saline (10 μL) or with 1, 10, or 100 μg of GDF‐6

Histological grade	Representative Decay Variance image	Sample size (n)	T2 score	Quantitative T2 score	Decay Variance score
4		24	0.59 (±0.08)	0.29 (±0.13)	0.49 (±0.07)
5‐11		24	0.64 (±0.08)	0.43 (±0.16)	0.62 (±0.07)
12		27	0.66 (±0.11)	0.43 (±0.16)	0.71 (±0.07)

Note: At 16 weeks postpuncture, MRI scans were performed in a 7 T MR scanner. Histological scores were compared to the average score in the region of interest in the IVD area generated by each postprocessing technique. Average scores from the three MRI techniques for IVDs in various states are presented here.

Abbreviations: GDF, Growth and Differentiation Factor; IVD, Intervertebral Disc; MRI, magnetic resonance imaging.

The three MRI postprocessing techniques all correlated with the histological scores for IVD degeneration (P < .01). T2 signal intensity was weakly correlated (r = 0.32, P < .01) with histological grade of IVD degeneration. Quantitative T2 relaxometry was weakly correlated (r = 0.39, P < .01) with histological grade of IVD degeneration. Decay Variance was strongly correlated (0.82) (P < .01) with histological grade of IVD degeneration (Table 2).

Table 2

The strength of correlation between each MRI postprocessing technique and the underlying histology score for intervertebral disc (IVD) degeneration: 4 weeks following the annular‐puncture in the lumbar IVDs in rabbits, the injured discs were either treated with phosphate buffered saline (10 μL) or with 1, 10, or 100 μg of GDF‐6

Technique	Correlation coefficient between histological grading of IVD and MRI value	P value
T2 weighted signal intensity	0.32	<.01
Quantitative T2 relaxometry	0.39	<.01
Decay Variance	0.82	<.01

Note: At 16 weeks postinjury, MRI scans of the lumbar spine were performed in a 7 T scanner. The average value in each disc area generated by each postprocessing technique was correlated with the severity of IVD degeneration captured using Masuda et al technique10 All the three techniques correlated significantly with the histology scores. The strength of correlation was greatest for the Decay Variance technique.

Abbreviations: GDF, Growth and Differentiation Factor; MRI, magnetic resonance imaging.

There were 25 IVDs without degeneration (ie, a histological IVD degeneration score of 4) and 50 IVDs with degeneration (ie, a histological IVD degeneration score ≥5).

Specificity with a minimum sensitivity of 90% was calculated for each technique (Table 3). Sensitivity with a minimum specificity of 90% was calculated for each technique (Table 4). ROC curves and tables were generated for every possible cutoff point for each of the three techniques. Complete sensitivity and specificity outcomes for every possible cut‐off are shown as ROC curves (Figure 3) and included in tabular form in the Supplementary Information section (SI 2).

Table 3

The best possible specificity for each technique with a minimum sensitivity of 90%

Technique	Sensitivity	Specificity
T2 signal intensity	90%	24%
Quantitative T2	90%	52%
Decay Variance	90%	92%

Note: For minimum 90% sensitivity.

Table 4

The best possible specificity for each technique with a minimum sensitivity of 90% (Note that due to the ratio of 50 positive to 25 negative discs no cutoff has a specificity of 90%, the next result was 92%)

Technique	Sensitivity	Specificity
T2 signal intensity	32%	92%
Quantitative T2	38%	92%
Decay Variance	92%	92%

Note: For minimum 90% sensitivity.

Figure 3

Receiver operating characteristics curves plotting sensitivity vs specificity for every possible cutoff point for each of the three techniques (T2 signal intensity at 13 ms, Quantitative T2 Relaxometry and Decay Variance) for 75 rabbit lumbar intervertebral discs (IVDs), 25 with no degeneration and 50 with degeneration induced by means of a needle puncture 16 weeks prior to magnetic resonance scanning. Decay Variance is superior at predicting the presence of IVD degeneration when compared with the other two techniques

The calculation time for the Decay Variance maps was less than 0.1% of the calculation time for the conventional quantitative T2 maps by non‐linear least squares curve fitting. The average calculation time for quantitative T2 maps was 1788.9 (±115.98) seconds and that for the Decay Variance maps was 1.02 seconds (±0.03) (Table 5).

Table 5

Calculation time (mean ± SD) to generate quantitative maps from magnetic resonance imaging data of rabbit lumbar spine using different techniques

Program	Method	Developer/Distributer	Calculation time (s)
Matlab Curvefit Toolbox	Nonlinear least squares	Mathworks USA	1788.8 (±116.0)
MR Map	Levenberg‐Marquardt	Deutsches Herzzentrum Berlin, German	37.3 (±3.5)
Segment	Maximum likelihood estimate	Lund Cardiac MRI Group/Medvisio Ab	5.0 (±0.9)
Decay Variance	Decay Variance	Authors of the present study	1.0 (±0.0)

Note: Quantitative T2 maps were calculated using nonlinear least squares in the Matlab curve fitting toolbox, the Levenberg Marquardt method using MRMap, and using a maximum likelihood estimate method using Segment (a commercially available vendor independent postprocessing tool). Calculations were performed with 32 GB of accessible RAM and 8 Intel Core i7‐7700K processor cores at 4.2 GHz.

Abbreviation: MRI, magnetic resonance imaging.

DISCUSSION

In this study, we developed a postprocessing tool for multi‐echo MRI scans that showed an excellent correlation with degenerative changes as quantified by histological grading in a rabbit lumbar IVD degeneration model. Decay Variance was able to separate degenerate from nondegenerate IVDs more reliably than quantitative T2 relaxometry or T2 signal intensity at 13 ms.

The calculation of quantitative T2 maps by nonlinear least squares methods assumes a single underlying exponential decay curve with a constant rate of signal decay.21 This assumption implies that each voxel contains only a single substance, which is valid only when performing an MRI on a chemically pure sample but will be invalid when applied to a biological material. Further, longitudinal relaxation cannot be “switched off” and will contribute to signal decay in the calculation of “T2” relaxation times. Reasonable acquisition times in living subjects limit the spatial resolution of MRI. The MRI scans in this study were acquired on a 7 T scanner with a voxel width of 180 μm. While this is a much finer spatial resolution than is generally available in clinical practice, it is a large area in histological terms (Figure 4).

Figure 4

Safranin O stained nucleus pulpous of a degenerate intervertebral disc. The micrograph (right) corresponds to the highlighted region (left) The entirety of the contents of this micrograph would fall within one pixel on a magnetic resonance image (180 × 180 μm)

The underlying decay curve for one voxel of a biological material (which likely contains a mixture of various cell types, proteins, lipids, water, and other molecules) is more likely to be a superposition of multiple discrete and varied decay curves of different weight. Thus, the degree to which the observed decay in signal intensity over time in a given pixel varies from a pure decay curve provides information, as a surrogate marker for the bio‐magnetic heterogeneity within the given pixel.

Calculation of quantitative T2 Maps by any accepted method is computationally intensive as it requires a curve fitting operation for each pixel. [21] Calculation of a Decay Variance map, by comparison, requires simple addition, subtraction, and division, and is, therefore, orders of magnitude less computationally intensive. Calculation of Decay Variance maps requires no special preparation or alteration to MRI image acquisition and can possibly be applied retrospectively to any multi‐echo T2 weighted sequence.

Low back pain is a condition afflicting 80% of adults in their lifetime, but the ability to identify causes on MRI using ordinal grading scales for IVD degeneration, such as described by Pfirrmann et al is limited.11, 12 In a study of 284 participants with low back pain, the authors reported no correlation between IVD degeneration quantified using Pfirrmann grade and pain or disability on standardized testing instruments, including the Roland‐Morris Disability Index, the Oswestry Disability Index, and The Short Form (SF)‐12 questionnaires.22

The need for better diagnostic techniques for low back pain assessment has been identified as a priority by both primary care practitioners23 and national funding bodies.24

This study has been conducted on an accepted animal model of IVD degeneration10; however, we do not have data yet on how this technique will perform on human lumbar MRIs and how the results may relate to low back pain.

In addition to potential intrinsic biological differences between humans and rabbits, this animal model further differs from human IVD degeneration principally in that most cases of IVD degeneration in the present rabbit model are severe. Most discs were either entirely in‐tact or very severely degenerate (degeneration scores of 11 or 12 out of 12), with relatively few cases of mild or moderate IVD degeneration.

Decay Variance has only been validated in T2 imaging and not in T2* imaging. The technique should be applied to and validated in T2* weighted multi‐echo imaging in the future. We have not tested Decay Variance in other tissues, and this is a direction for future research. A limitation of this study is that the imaging was performed at a single time point postmortem. Future studies in humans should consider repeated measures to assess change in Decay Variance over time.

These preliminary results show that Decay Variance is a novel and computationally efficient technique, which has a stronger correlation with histological defined severity of IVD degeneration than quantitative T2 weighted imaging.

Further validation on human subjects and at other magnetic field strengths is needed, but if the correlation with histological changes observed in the present study translates to a correlation with disease severity in humans, the implementation of this technique in clinical practice will allow for objective and accurate assessment of IVD health.

CONFLICT OF INTEREST

A.D.D. reports personal fees from Nuvasive Inc, other from Nuvasive Inc, other from Kunovus, other from Kunovus Technologies, outside the submitted work; In addition, A.D.D. and K.S. have a patent in discussion with the UNSW Tech Transfer Office related to this work. The other authors have no conflict of interest to declare.

AUTHOR CONTRIBUTIONS

K.S. developed the Decay Variance algorithm and code and performed the image postprocessing analysis. S.M. performed the histology. K.K. acquired the MRI scans. U.C., K.M., and A.D.D. contributed to study design. All authors contributed to manuscript preparation.

Appendix S1. Supporting information.

Click here for additional data file.

Appendix S2. Supporting information.

Click here for additional data file.

Appendix S3. Supporting information.

Click here for additional data file.