Segmental differences of cervical spinal cord motion: advancing from confounders to a diagnostic tool

Increased cranio-caudal spinal cord motion is associated with clinical impairment in degenerative cervical myelopathy. However, whether spinal cord motion holds potential as a neuroimaging biomarker requires further validation. Different confounders (i.e. subject characteristics, methodological problems such as phase drift, etc.) on spinal cord motion readouts have to be considered. Twenty-two healthy subjects underwent phase contrast MRI, a subset of subjects (N = 9) had repeated scans. Parameters of interest included amplitude of velocity signal, maximum cranial respectively maximum caudal velocity, displacement (=area under curve of the velocity signal). The cervical spinal cord showed pulse synchronic oscillatory motions with significant differences in all readouts across cervical segments, with a maximum at C5. The Inter-rater reliability was excellent for all readouts. The test-retest reliability was excellent for all parameters at C2 to C6, but not for maximum cranial velocity at C6 and all readouts at C7. Spinal cord motion was correlated with spinal canal size, heart rate and body size. This is the first study to propose a standardized MRI measurement of spinal cord motion for further clinical implementation based on satisfactory phase drift correction and excellent reliability. Understanding the influence of confounders (e.g. structural conditions of the spine) is essential for introducing cord motion into the diagnostic work up.

Therefore, a reliable and sensitive assessment of cervical spinal cord motions along the cervical spine is required in order to foster its clinical application and potential introduction in the diagnostic work up of DCM.
The aim of this study was to investigate cranio-caudal spinal cord motions across the entire cervical spine in healthy volunteers, to evaluate its inter-rater and test-retest reliability for clinical application and its relation to anatomic and biometric conditions.

Methods
Population. This prospective study (2016)(2017)(2018) was approved by the local ethics committee (Kantonale Ethikkommission Zurich, KEK-ZH 2012-0343, BASEC Nr. PB_2016-00623). All methods were carried out in accordance with the relevant guidelines and regulations. Informed consent was obtained from all participants. 22 healthy subjects (age 63.3 ± 6.6 years) were recruited. Four subjects with incidental stenosis in structural imaging had to be excluded. Eighteen participants (9m, 9f; age 62.2 ± 6.5 years) were included in the final statistical analysis. Data on blood pressure before and after the MRI, age, body weight and body size were collected. Blood pressure, weight and body size was missing in one participant each.
Imaging. All subjects underwent a MRI scan (3T, MAGNETOM SkyraFit, Siemens Healthcare, Erlangen) including sagittal and axial standard clinical T2w sequences and axial 2D phase contrast imaging encoding cranio-caudal spinal cord motion in each cervical segment (Fig. 1a). Slice orientation in phase contrast imaging was adjusted perpendicular to the spinal cord. The phase contrast sequence was cardiac gated through peripheral pulse triggering using a finger clip. The venc (velocity encoding) value of the phase contrast sequence was set to 2 cm/s based on previous findings of cord motion [4][5][6][7]18 . During one cardiac cycle velocity signal was assessed within 20 time points and 128 measurements were averaged per segment. Total scanning time of the whole protocol was 23 min. Imaging analysis. Image analysis was performed using the OsirixTM free DICOM viewer (www.osirix.com).
Anatomical measurements. In sagittal T2w images the anterior-posterior diameter (ap) of the spinal canal at disc level, in axial images the anterior-posterior (ap) and right-left (rl) diameter, as well as spinal canal and spinal cord cross-sectional area (CSA) at disc level were measured in all cervical segments. To evaluate the apex of the cervical spine curvature tangential lines were drawn at the dorsal vertebral bodies with the point of intersection identifying the apex and additionally the angle was collected.
Velocity calculation. In phase contrast measurements, cranio-caudal spinal cord motion was analysed by a predefined ellipsoid shaped region of interest (30,52 mm 2 ) mid-centred into the spinal cord (Fig. 1b). According to the predefined velocity encoding (2 cm/s) of phase contrast imaging velocity signal within the MRI sequence ranges from −2 cm/s to 2 cm/s encoded in grey values from −4096 to 4096 in the Osirix DICOM viewer within 20 time points during one cardiac cycle. The mean of the measured greyscale values within the region of interest in each time point was exported with the extension "ROI export" of the Osirix program. The measured mean Steps in spinal cord motion evaluation. Spinal cord cranio-caudal motion was measured by axial MRI phase contrast imaging in each cervical segment (a; sagittal T2 image; red lines identifying the cervical segments; arrow illustrating motion). Motion was encoded into grey values of the MRI image within 20 timepoints of one heartcycyle. Velocity values were calculated by the mean of collected grey value within the predefined region of interest (b; axial phase contrast image; red circle illustrating predefined ROI) and the used velocity encoding in the sequence (2 cm/s). Due to a technically caused offset error (c; uncorrected motion signal during one heartcycle (0-100%); arrow indicating offset error) raw velocity values were misleading overestimated. Subtracting the mean of all 20 velocity values from each single value could sufficiently correct for the offset error (d; corrected motion signal during one heartcycle (0-100%)). Readouts of spinal cord motion used were amplitude of the motion signal, ranging from the maximum negative to the maximum positive value, maximum cranial velocity, maximum caudal velocity and the displacement (=area under curve of velocity signal). grey value was divided by 4096 and subsequently multiplied by 2 cm/s (velocity encoding of the phase contrast sequence) for calculation of the velocity.
Parameters of interest within one cardiac cycle included amplitude of the velocity signal ranging from the maximum negative to the maximum positive peak, maximum cranial respectively maximum caudal velocity, and the area under the curve (AUC) of the movement signal. AUC was calculated by stepwise summation of calculated squared areas in each time point (1/20 RR time multiplied by the measured velocity). Negative velocity values were transformed to a positive value for calculation of the area under curve.
Phase drift correction of spinal cord motion measurement. As phase contrast imaging is a relative measure of motion a so-called "phase drift" 10 leads to an offset error of the raw data ending in misleading velocity values (Fig. 1c). Therefore, phase drift conducts to an over-or underestimation in the velocity measurements and a correction for phase drift is needed. In our analysis phase drift correction was done by subtraction of the mean of all 20 velocity measurements within one cardiac cycle from the raw velocity value at each time point (Fig. 1d), as net motion of the spinal cord over one cardiac cycle is assumed to be 0 (start and end location of the spinal cord is expected to be at the same position -assuming mean velocity has to be 0).

Statistical analysis.
Statistical analysis was conducted using SPSS (IBM, Version 23). Metrics are reported as group mean ± SD. Between segment differences were calculated using the Friedmann test. For reliability analysis intraclass correlation coefficients (ICC) were calculated (two way mixed model, absolute agreement, average measures). For bivariate correlations Spearman Rho coefficients were calculated. Significance level of alpha was set to <0.05.
Inter-rater reliability. Inter-rater reliability was assessed by calculating intraclass correlation coefficients (two-way mixed model, absolute agreement, average measures). Two independent raters placed the ROI midcentred into the spinal cord as explained above (velocity calculation) in all cervical segments of each healthy participant (n = 18).
Test-retest reliability. For evaluation of test-retest reliability 9 of the 18 healthy volunteers were scanned a second time after 732+/− 164 days. Velocity measurements were done by the same rater at the first and second scan and intraclass correlation coefficient (two-way mixed model; absolute agreement; average measures) were calculated.
Spinal cord motion pattern. A characteristic, comparable motion pattern over all cervical segments could be observed (Fig. 2), with nearly no motion in the first half of the heart cycle (measured by peripheral pulse triggering), followed by a cranio-caudal and afterwards caudo-cranial oscillatory movement.
Phase drift correction. In our measurements a positive "phase drift" about 0.2-0.3 cm/s could be observed in all participants, which would have led to a misleading systematic overestimation of spinal cord motion (Supp. Fig. 1). By subtracting the mean velocity of all 20 time points from the raw value at each single time point in each segment, motion signal over one heart cycle could be sufficiently corrected (Fig. 2).
Interrater reliability. The inter-rater reliability of the 18 baseline measurements between two independent raters was with ICCs of almost 1 (p < 0.001) excellent in all cervical segments (Table 1).

Measured spinal cord motion values.
The mean amplitude of spinal cord motion ranged from 0.296 ± 0.088 cm/s in segment C2 to 0.445 ± 0.163 cm/s at level C5 (Table 3). Maximum cranial velocity ranged from 0.115 ± 0.041 cm/s at C2 to 0.186 ± 0.078 cm/s at C5. Maximum caudal velocity ranged from −0.169 ± 0.055 cm/s in segment C7 to −0.259 ± 0.094 cm/s at C5. For the displacement, the highest values could also be observed at C5 (0.054 ± 0.016 cm), with lowest values at C2 (0.036 ± 0.009 cm). Raw, uncorrected values are provided in Supplementary Material for comparison (Supp. Table 1).

Correlations of spinal cord motion to anatomic measures. Greater spinal cord displacement
at cervical level C5 was correlated with smaller anterior-posterior diameter in sagittal imaging (r = −0.540; p = 0.021) (Fig. 4a). Considering spinal cord motion across all cervical segments, higher displacement correlated

Amplitude
Max  Table 3. Segmental spinal cord motion readouts. Values of spinal cord motion were calculated for corrected velocity values in each cervical segment for amplitude (cm/s), maximum cranial (cm/s) and caudal velocity (cm/s) and for the displacement (cm). max = maximum, SD = standard deviation. www.nature.com/scientificreports www.nature.com/scientificreports/ Higher displacement at C2 correlated with smaller body size (r = −0.566; p = 0.018; Fig. 5b) and smaller mean arterial pressure (r = −0.487; p = 0.047; plot not shown). At C3 (r = −0.512, p = 0.036) and C5 (r = −0.523, p = 0.031) a correlation between higher displacement to lower systolic blood pressure after imaging was observed (plots not shown). Sex and body weight were not correlated with any readout.

Discussion
This study is the first measuring physiological spinal cord motions across all cervical spinal segments. Physiological, pulse-synchronic oscillations of the spinal cord have been observed and described in prior studies 2,5,18-22 , but mainly focused on individual stenotic and adjacent cervical segments in DCM patients. In our healthy study cohort, comparable motion patterns in all cervical segments could be observed with nearly no motion in the first half of the heart cycle (measured by peripheral pulse triggering), followed by biphasic cranio-caudal oscillatory movement (Fig. 2). This finding corresponds to studies observing analogue cervical spinal cord motions (first caudal and then cranial movements) using ECG triggering 2,5 .
All spinal cord motion readouts except amplitudes are highly influenced by phase drift of MRI measurements. Phase drift was reported to differ between scanners and also between different time points in the same scanner 10 . Phase drift correction was conducted by phantom measurements in earlier MRI studies 2,5,18 , but is time consuming and less applicable in a clinical setting. Also measurements of static tissue next to the spinal canal were used for calculation of the background velocity reflecting phase drift 18 . However, measuring static tissue was not reliable in our study between two different raters (data not shown). Vavasour et al. 6 adjusted velocity curves by adding or subtracting a uniform velocity from each cardiac time point to make the net motion over a cardiac cycle equal to zero but did not report how this value was obtained. In other studies assessing spinal cord motions no systematic correction for phase drift was mentioned 4,7 , while Chang et al. 4 did not need to correct as they used the half amplitude of the velocity signal, not altered by phase drift.
We propose phase drift correction by subtracting the the mean of all velocity values within one heart cycle from each single value, assuming oscillatory spinal cord motion starts and ends at the same position and net motion is equal to 0. The same assumption was reported before 6 , but no systematic calculation for correction was provided. Using this simplified calculation, sufficient correction could be shown in all measurements. Most  www.nature.com/scientificreports www.nature.com/scientificreports/ important, the amplitude of the velocity signal is not affected by the phase drift and may therefore be an attractive parameter for future studies on spinal cord motion as no correction is needed.
Most studies did not report on reliability of spinal cord motion measurements 2,[4][5][6]18 . In our study interrater reliability between two independent raters was excellent using a predefined ROI placed mid-centred into the spinal cord for all parameters at all levels. These results are in line with prior reported good inter-rater reliability using a smaller region of interest (19.68 mm 2 ), also mid-centred into the spinal cord 7 . Using a bigger ROI in our study further improved inter-rater reliability.
Our study is the first to report the test-retest reliability in phase contrast spinal cord motion measurements. Even though second measurements were performed in mean after two years (732+/− 164 days) in stable healthy controls test-retest reliability remained excellent in all segments except at C7, using the proposed phase drift correction. However, applying raw values without correction test-retest reliability was inferior, i.e. the displacement values were not reliable. Phase drift can differ between two measurements 10 , ending in different absolute raw values and low test-retest reliability, underlining the importance of systemic phase drift correction. Poor reliability at C7 might be mostly attributed to inferior signal to noise ratio in this segment as distance to the head coil in the MRI scanner is high 23 . Due to sufficient phase drift correction in the two different measurements, the proposed method should also be sufficient in different scanners with different phase drift 10 , but this was not evaluated in the present study.
Parameters of interest used in prior analyses were amplitude 4 , total displacement 6,7 , maximum displacement 6 , peak velocities 2,22 , mean velocity 6 and velocity difference between cervical segments 5 . However, due to different readouts comparability between measurements is limited. We introduce a standardized MRI evaluation, including a phase drift correction, ending up in four different parameters (amplitude, max cranial respectively max caudal velocity, displacement). Wolf et al. 7 did not correct for phase drift and therefore a overestimation of velocity measures cannot be excluded. Chang et al. 4 reported a mean value of cord motion of 5 ROIs at vertebral body level C2 to C6 of 0.28 cm/s for the half peak to peak amplitude, while no individual segmental values were reported. As significant differences could be shown between these segments in our study, calculating a mean appears not reasonable. Vavasour et al. 6 reported at C5 disc level a mean cord velocity of 0.144 cm/s and a displacement of 0.133 cm in controls, comparable to a mean velocity of 0.188 cm/s respectively a displacement of 0.176 cm in our uncorrected raw data. A correction for phase drift is mentioned in their study, but no detailed information was provided.
Impact of different confounders on spinal cord motion might explain discrepancies between reported velocity values along different studies. Magnetic field strength might influence accuracy of measurements. As initial MRI measurements were performed using lower field (0.5-1.5T) strength 2,5,18,22 , recent studies mostly measured at higher field (3T) MRIs 4,6,7 . Additionally, the influence of the chosen velocity encoding for the phase contrast MRI sequence on accuracy of motion measurements remains open. In our and some prior studies 5,18 an encoding of 2 cm/s, focusing on slow spinal cord motion was implemented, in contrast to 5 cm/s 2,7 and 6 cm/s 4 in other studies. As MRI digital imaging and communications in medicine (DICOM) data system encodes for 4096 grey values only, higher velocity encoding might end in lower accuracy of measurements of slow velocities like in spinal cord motions. Therefore, velocity encoding should be adapted as far as possible to expected velocity values.
Also the shape and the size of the region of interest will influence velocity values as they are collected as a mean value within the ROI. On the one hand it should cover a large part of the spinal cord, but partial volume effects at the spinal cord -CSF border should be avoided, and the size should be feasible to be used in stenotic segments with the cord being squeezed (i.e deformed). In contrast to our ellipsoid ROI of 30,52 mm 2 , prior studies used a rectangle ROI of 19.68 mm 2 7 , a manually drawn ROI around the spinal cord 6 , a round ROI covering the ap diameter of the spinal cord 4 , a circular ROI 2,5 with no information about the size and also differing sizes between patients in manually drawn ones 4,6 .
Influence of age on spinal cord motion, as measurements were performed in younger volunteers before 2,4,5,7,18 remains unknown. In our population no relation of spinal cord motion to age could be found, but all participants were older than 50 years comparable to only one prior study 6 .
In physiological conditions, the origin of spinal cord motions was mostly attributed to remote factors like intracranial CNS-and CSF-pulsation within the cardiac cycle 5 and breathing 24 . Breathing was not assessed in our study, but relation of higher spinal cord motion to lower RR time respectively higher heart rate could be observed. Interestingly studies assessing pulse wave velocities and blood flow in blood vessels showed also about 20% higher values with increasing heart rate [25][26][27][28] , underlining the impact of arterial pulsation on spinal cord motion. Additionally, focal influences on spinal cord motion due to expansion of segmental arteries could be shown in a canine model as oscillations diminished after transection of the local vascular support 21 . However, association of higher spinal cord motion to lower blood pressure in our measurements is contradicting previous findings in arterial blood vessels 29,30 and remains unexplained.
Intersegment comparisons in our study revealed significant differences between various cervical segments. Highest spinal cord motion could be observed at C5. No systematic evaluation of intersegment differences is provided in most prior studies 2,4,6,22 , or only in patients between compressed and non-compressed segments 4,5 . Wolf et al. 7 report a displacement ratio between C2 and C5 of 1 in healthy controls, reflecting no differences between motions in these segments, but uncorrected, potentially overestimated values were used for calculation. Max spinal cord motion was also reported at C2 and C5 in the latter study, but no statistical analysis about significant differences between these segments, as only differences to patients were evaluated.
Spinal cord motion was negatively correlated to spinal canal measures. These results are in line to the law of Hagen-Poiseuille 31 . The relation to anatomic conditions is comparable to reports in CSF dynamics before 32,33 . Specifically CSF velocity was higher in smaller spinal canal and smaller subarachnoidal space also with a maximum at C5 32,33 , but no relationship of spinal cord motion to anatomic conditions could be found in prior evaluations in healthy controls 6,7 . Additionally, a correlation between spinal cord and CSF motion at C2 was reported www.nature.com/scientificreports www.nature.com/scientificreports/ before 7 . As CSF-flow was frequently discussed to influence spinal cord motion 1,34-37 , our results underline the close interaction between CSF and spinal cord dynamics, with highest motion values at C5 either. Due to low interrater reliability of CSF measurements on level of stenosis in DCM patients 7 and low velocity encoding used in our study no CSF assessments were done.
Correlation of higher spinal cord motion to smaller body size in our population has to be attributed to smaller anatomic conditions in smaller people, as correlations of smaller body size to smaller anatomic measures, i.e. spinal canal could be observed (data not shown).
The configuration of the spine, i.e. the apex and the angle did not influence spinal cord dynamics in our population. Also no relation to weight and age could be observed. No prior evaluation on this could be found in literature.
Interestingly the C5 level exhibited the highest spinal cord motion and this level is also one of the most affected by spinal canal stenosis and in DCM 34 . Increased spinal cord motion, especially as shown in DCM patients 4-7 might reflect increased mechanical stress to the spinal cord conducting to damage and deterioration.
Different MRI techniques have been applied to evaluate the spinal cord complementary to standard clinical sequences, aiming for additional information about the cords tissue integrity. Most recently a multimodal MRI protocol in DCM was introduced 38,39 . The protocol including diffusion tensor imaging (DTI), magnetisation transfer ratio (MTR), T2* white matter/grey matter contrast ratio and T2 cross sectional area measurements at, below and above the level of the spinal stenosis shows promising ability to reveal spinal cord damage and deterioration over time, even in asymptomatic patients 40,41 . However, these advanced microstructural measurements require extensive postprocessing, and still await to be implemented in the daily clinical workup, while spinal cord motion measurements can be easily evaluated. Also, these MRI sequences are more time consuming compared to phase contrast measurements. In addition, increased spinal cord motion might indicate the risk of spinal cord damage due to increased mechanical stress, even if evident microstructural damage could not yet be observed. Therefore, it might be able to identify high risk patients and could help to decide for decompressive surgery.

Conclusion
This study provides promising data for a standardized evaluation and corrections for systemic bias in phase contrast imaging and further clinical implementation of spinal cord motion. The proposed protocol is applicable to develop normative values (common or per site) of cranio-caudal cervical spinal cord motions, to discern and monitor pathologic conditions and its association to DCM where spinal cord motion is influenced by spinal canal size, heart rate, body size and blood pressure.

Data Availability
The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.