Imaging and / or detecting pulsatile displacement using magnetic resonance imaging

Abstract

The present disclosure relates to imaging and / or detecting pulsatile displacement using MRI. A scan is performed on an imaging subject using an MRI device by causing the MRI device to apply an MRI pulse sequence with a 2D- or 3D- based readout to an imaging subject and acquire MR data of the imaging subject. The MRI pulse sequence comprises a plurality of repetitions TR. Each of the repetitions has displacement encoding gradient for sensitization to pulsatile displacement and has the 2D- or 3D- based readout applied either before or after the displacement encoding gradient. The displacement encoding gradient comprises a gradient lobe with a finite zero-order moment. The MR data is processed to obtain phase and / or magnitude data associated with the pulsatile displacement.

Description

[0001] IMAGING AND / OR DETECTING PULSATILE DISPLACEMENT USING MAGNETICRESONANCE IMAGING TECHNICAL FIELDEmbodiments disclosed herein relate to techniques (e.g., systems and methods) forimaging and / or detecting pulsatile displacement using magnetic resonance imaging (MRI). BACKGROUNDFluid and tissue motion may be present in the brain (human or animal) on various time-and length-scales and may be associated with brain health. For example, thecerebrospinal fluid (CSF) may play a key role in clearance of waste substances from thebrain. The pulsatile motion of CSF may be the result of a complex interplay between thedriving forces such as cardiac pulsation and respiratory motion and the anatomical and mechanical properties of the brain.Sensitive and / or robust methods for detecting, imaging, and / or characterizing the motionof CSF and / or surrounding brain tissue, on short and / or long timescales, may be usefulto advance the understanding of brain pulsation (healthy and pathological). MRI is intrinsically sensitive to motion due to the application of spatial B0gradients aspart of spatial encoding. For example, in phase-contrast MRI (PC-MRI), the motionsensitivity may be enhanced by the addition of extra bipolar gradients, which may enable detection and / or quantification of fluid flow. Problematically, PC-MRI may not be suitable for measuring low flow velocities (such asCSF flow), as doing so would require large gradient moments, which may increase theecho-time and overall scan duration. This means that existing PC-MRI techniques may suffer from low sampling efficiency and signal-to-noise ratio (SNR). SUMMARY In a first aspect, there is provided a method of imaging and / or detecting pulsatile displacement using MRI. The method comprises performing a scan on an imaging subject using an MRI device by causing the MRI device to: apply an MRI pulse sequencewith a 2D- or 3D- based readout to an imaging subject and acquire MR data of theimaging subject. The MRI pulse sequence comprises a plurality of repetitions (TR). Eachof the repetitions includes displacement encoding gradient for sensitization to pulsatiledisplacement and has the 2D- or 3D- readout applied either before or after thedisplacement encoding gradient. The displacement encoding gradient comprises agradient lobe with a finite zero-order moment. The method comprises processing the MRdata to obtain phase and / or magnitude data associated with the pulsatile displacement. The pulsatile displacement may be associated with physiological motion such as cardiac

[0002] PC933674WO - specificationmotion and / or respiratory motion. The phase and / or magnitude data may be processedfor display as MRI image(s). In some embodiments of the first aspect, the zero-order moment of the gradient lobe has a magnitude of at least 2π / γ∆x, where ∆x is voxel size of the MR data and γ is gyromagnetic ratio.In some embodiments of the first aspect, if the 2D- or 3D- based readout is applied beforethe displacement encoding gradient, the 2D- or 3D- based readout may acquire a freeinduction decay (FID) signal; and if the 2D- or 3D- based readout is applied after thedisplacement encoding gradient, the 2D- or 3D- based readout may acquire an echosignal. In some embodiments of the first aspect, the pulsatile displacement is pulsatile displacement of a biological tissue and / or a biological fluid. For example, the pulsatile displacement may be pulsatile displacement of at least one of: brain tissue (such as white matter) or CSF.In some embodiments of the first aspect, the 2D- or 3D- based readout is an echo planarimaging (EPI) based readout. In one example, the EPI based readout is a 3D EPI basedreadout and the MR data comprises volumetric MR data. For example, the 3D EPI basedreadout may be a segmented 3D EPI readout. In some embodiments of the first aspect, for each of the repetitions, the MRI pulsesequence lacks velocity encoding gradients before the 2D- or 3D- based readout.In some embodiments of the first aspect, the displacement encoding gradients of the plurality of repetitions comprise orthogonal gradients for sensitization to pulsatile displacements along orthogonal directions. In some embodiments of the first aspect, for each of the repetitions, the MRI pulse sequence applies RF spoiling.In some embodiments of the first aspect, for each of the repetitions, the MRI pulsesequence does not apply RF spoiling. In some embodiments of the first aspect, the MRI pulse sequence is a non-balancedsteady-state free precession (nb-SSFP) sequence.In some embodiments of the first aspect, the MRI pulse sequence is a gradient-echo (GRE) based sequence. In some embodiments of the first aspect, an amplitude of the pulsatile displacement has an order of magnitude of 10-3mm, an order of magnitude of 10-2mm, an order ofmagnitude of 10-1 mm, or an order of magnitude of 1 mm. In some embodiments of thefirst aspect, a temporal resolution of the pulsatile displacement is in an order of 1 s to 10-1s.

[0003] PC933674WO - specification In some embodiments of the first aspect, the imaging subject comprises a brain of a subject, and the acquisition of the MR data is based at least in part on prospectiveelectrocardiogram (ECG) or respiratory gating associated with the subject.In some embodiments of the first aspect, the imaging subject comprises a brain of a subject, and the processing of the MR data is based at least in part on retrospective ECG or respiratory gating associated with the subject. In some embodiments of the first aspect, the method further comprises converting phase data to displacement data associated with the pulsatile displacement based at least inpart on a lookup or mapping table, which has information for associating phase withdisplacement. The displacement data may be processed for display as MRI image(s). In some embodiments of the first aspect, the method further comprises converting phase data to displacement data associated with the pulsatile displacement based at least inpart on a model, such as a Bloch Simulation model or an extended phase graph (EPG)model, which can associate phase with displacement. The displacement data may beprocessed for display as MRI image(s).In a second aspect, there is provided a system for imaging and / or detecting pulsatiledisplacement using MRI. The system includes an MRI device and at least one processor.The MRI device is configured to: apply an MRI pulse sequence with a2D- or 3D- basedreadout to an imaging subject and acquire MR data of the imaging subject. The MRI pulse sequence comprises a plurality of repetitions (TR). Each of the repetitions includesdisplacement encoding gradient for sensitization to pulsatile displacement and with the2D- or 3D- based readout applied either before or after the displacement encodinggradient. The displacement encoding gradient comprises a gradient lobe with a finite zero-order moment. The at least one processor is configured to: process the MR data toobtain phase and / or magnitude data associated with the pulsatile displacement. Thepulsatile displacement may be associated with physiological motion such as cardiacmotion and / or respiratory motion. The phase and / or magnitude data may be processedfor display as MRI image(s). In some embodiments of the second aspect, the zero-order moment of the gradient lobe has a magnitude of at least 2π / γ∆x, where ∆x is voxel size of the MR data and γ is gyromagnetic ratio.In some embodiments of the second aspect, if the 2D- or 3D- based readout is appliedbefore the displacement encoding gradient, the 2D- or 3D- based readout may acquirea FID signal; and if the 2D- or 3D- based readout is applied after the displacementencoding gradient, the 2D- or 3D- based readout may acquire an echo signal.In some embodiments of the second aspect, the pulsatile displacement is pulsatile displacement of a biological tissue and / or a biological fluid. For example, the pulsatile displacement is pulsatile displacement of at least one of: brain tissue (such as white matter) or CSF.

[0004] PC933674WO - specificationIn some embodiments of the second aspect, the 2D- or 3D- based readout is a EPI basedreadout. In one example, the EPI based readout is a 3D EPI based readout and the MRdata comprises volumetric MR data. For example, the 3D EPI based readout may be asegmented 3D EPI readout. In some embodiments of the second aspect, for each of the repetitions, the MRI pulsesequence lacks velocity encoding gradients before the 2D- or 3D- based readout.In some embodiments of the second aspect, the displacement encoding gradients of theplurality of repetitions comprise orthogonal gradients for sensitization to pulsatile displacements along orthogonal directions. In some embodiments of the second aspect, for each of the repetitions, the MRI pulse sequence applies RF spoiling. In some embodiments of the second aspect, for each of the repetitions, the MRI pulse sequence does not apply RF spoiling. In some embodiments of the second aspect, the MRI pulse sequence is a nb-SSFP sequence. In some embodiments of the second aspect, the MRI pulse sequence is a GRE based sequence. In some embodiments of the second aspect, an amplitude of the pulsatile displacement has an order of magnitude of 10-3mm, an order of magnitude of 10-2mm, an order ofmagnitude of 10-1 mm, or an order of magnitude of 1 mm. In some embodiments of thesecond aspect, a temporal resolution of the pulsatile displacement is in an order of 1 sto 10-1s. In some embodiments of the second aspect, the imaging subject comprises a brain of a subject, and the MRI device is configured to acquire the MR data based at least in part on prospective ECG or respiratory gating associated with the subject. In some embodiments of the second aspect, the imaging subject comprises a brain of a subject, and the at least one processor is configured to process the MR data based at least in part on retrospective ECG or respiratory gating associated with the subject. In some embodiments of the second aspect, the at least one processor is configured to convert phase data to displacement data associated with the pulsatile displacementbased at least in part on a lookup or mapping table, which may contain information forassociating phase with displacement. The displacement data may be processed for display as MRI image(s). In some embodiments of the second aspect, the at least one processor is configured to convert phase data to displacement data associated with the pulsatile displacement

[0005] PC933674WO - specificationbased at least in part on a model, such as a Bloch Simulation model or an EPG model,which may be arranged to associate phase with displacement. The displacement datamay be processed for display as MRI image(s).In a third aspect, there is provided a computer-implemented method that comprisesreceiving MRI data associated with a movable substance. The MRI data is acquired by an MRI device using a plurality of MRI pulse sequences operable to encode motion. In some cases, the plurality of MRI pulse sequences may each correspond to a single cycle and may be combined to provide a single MRI pulse sequence with multiple repetitions (TR). The computer-implemented method further comprises processing the received MRI data based on a model to detect and / or quantify motion associated with the movable substance. The model is arranged to associate phase information of MRI data with motion. By processing the received MRI data using the model, motion, in particularpulsatile motion, associated with the movable substance may be detected and / orquantified. The pulsatile motion may be periodic or aperiodic. Optionally, the detection and / or quantification of motion comprises a determination of a direction of the motion. Optionally, the detection and / or quantification of motion comprises a determination of an amplitude or relative amplitude of the motion. Optionally, the detection and / or quantification of motion comprises a determination of a velocity or relative velocity of the motion. Optionally, each of the plurality of MRI pulse sequences respectively comprises a radiofrequency pulse for excitation and defines a repetition time ^^^^, and the model maybe defined based on where ^^^ is transverse coherence state of order ^^, + indicates a transverse coherencestate immediately after application of a radiofrequency pulse in the MRI pulse sequence,− indicates a transverse coherence state immediately before application of theradiofrequency pulse in the MRI pulse sequence, ^^^^ is the repetition time of the MRIpulse sequence, ^^2 is transverse relaxation time of the movable substance, ^^ is aninteger indicating number of TR-periods that the movable substance has experienced,and ^^^^ is a phase term that considers an effect of the motion associated with themovable substance on phase information of the received MRI data. In one example, adirection measure of motion of the movable substance may be provided by a sign of the phase term. Optionally, each of the plurality of MRI pulse sequences respectively comprises a readout sequence and gradient spoiling applied along a direction, and the phase term may be defined as: ^^^^^(^∙^^)where ^^(^^ ∙ ^^^^) at least includes a sum of of ^^ (^^ ∙ ^^^^) ∙ ∫^^ ^^⃗^(^^) ^^^^ and ^⃗^(^^ ∙ ^^^^) ∙∫^^ ^^⃗^(^^) ^^ ^^^^ , where ^^ is a position (displacement) measure of the movable substance, ^⃗^

[0006] PC933674WO - specificationis a velocity measure of the movable substance, ^⃗^ is the gradient spoiling applied in theMRI pulse sequence, and ^^ is time.Optionally, processing the received MRI data based on the model comprises obtaining (e.g., extracting) phase information from the received MRI data, and processing the phase information based on the model to determine a position (displacement) measureand optionally a velocity measure of the movable substance.Optionally, each of the plurality of MRI pulse sequences is based on a non-balanced steady-state free precession technique and comprises a readout sequence. It has been found that the use of non-balanced steady-state free precession technique can effectively encode motion and is in some cases particularly suited for use with the model. Optionally, the plurality of MRI pulse sequences comprise a plurality of radiofrequency pulses for excitation, and the plurality of radiofrequency pulses have a linear or quadraticincrement of phase ^^. The linear or quadratic increment of phase ^^ may enable the MRIdata acquired by the MRI device to more closely follow the Ernst equation, which facilitates subsequent processing of the MRI data based on the model. Optionally, the readout sequence comprises an echo planar imaging sequence. For example, the echo planar imaging sequence may comprise a 3-dimensional echo planar imaging sequence. Optionally, the non-balanced steady-state free precession technique comprises gradient spoiling applied along a single direction. Optionally, the gradient spoiling is applied before the readout sequence. In this case, the readout sequence may acquire an echo signal. Optionally, the gradient spoiling is applied after the readout sequence. In this case, the readout sequence may acquire a FID signal. Optionally, the plurality of MRI pulse sequences comprises: a first type of MRI pulse sequence comprising a readout sequence and gradient spoiling applied along a first direction; a second type of MRI pulse sequence comprising a readout sequence and gradient spoiling applied along a second direction perpendicular to the first direction; and a third type of MRI pulse sequence comprising a readout sequence and gradient spoilingapplied along a third direction perpendicular to each of the first direction and the seconddirection. This arrangement facilitates 3-dimensional motion encoding hence motion detection and / or quantification. Optionally, the movable substance is in a body of a subject. The subject may be human or animal. Optionally, the movable substance comprises brain tissue, such as white matter. Optionally, the movable substance comprises body fluid, such as cerebrospinal fluid.

[0007] PC933674WO - specification Optionally, the processing of the received MRI data is further based on physiological data of the subject, and the physiological data is obtained from the subject while the MRIdata is acquired. In other words, a retrospective gating technique is applied to organizethe received MRI data. This can help to improve accuracy of the detection and / or quantification. Optionally, the processing of the MRI data comprises temporally synchronizing or aligning the received MRI data with the physiological data of the subject, and processing the temporally synchronized or aligned MRI data based on the model to detect and / or quantify motion associated with the movable substance. Optionally, the MRI data is acquired from the subject by the MRI device based on physiological data of the subject. In other words, a prospective gating technique is applied. In this way, the MRI data acquired and hence received would have been synchronized with the one or more physiological cycles (e.g., one or more cardiac or respiratory cycles) of the subject. Optionally, the physiological data comprises data associated with one or more physiological cycles. Optionally, the physiological data comprises data associated with one or more cardiac cycles of the subject. Optionally, the physiological data comprises data associated with one or more respiratory cycles of the subject. Optionally, the MRI data acquired by the MRI device is non-contrast MRI data (i.e., the MRI data is acquired without using any extrinsic MRI contrast agent). Optionally, processing the received MRI data based on the model comprises: processing the received MRI data to reconstruct MRI images, and processing the MRI images based on the model to detect and / or quantify motion associated with the movable substance. Optionally, the computer-implemented method further comprises causing the MRI device to acquire the MRI data using the plurality of MRI pulse sequences.In a fourth aspect, there is provided a system comprising at least one processorconfigured to perform the computer-implemented method of the third aspect. Optionally,the system further comprises the MRI device arranged to acquire the MRI data using the plurality of MRI pulse sequences and operably connected with the processor to provide the MRI data to the at least one processor. In a fifth aspect, there is provided a carrier medium comprising instructions which, whenexecuted by at least one processor, cause the at least one processor to perform thecomputer-implemented method of the third aspect. For example, the carrier medium maybe a computer-readable storage medium. For example, the computer-readable storage medium may be a non-transitory computer-readable storage medium. In a sixth aspect, there is provided a data processing system comprising meansfor carrying out the steps or operations of the computer-implemented method of the thirdaspect.

[0008] PC933674WO - specification In a seventh aspect, there is provided a computer program product comprising instructions which, when the program is executed by a computer, cause the computer tocarry out the computer-implemented method of the third aspect.In an eighth aspect, there is provided a method for identifying brain disorder in a subject,comprising: (i) performing the method of the first aspect, wherein the imaging subjectcomprises a brain of a human or animal subject, and comparing phase and / or magnitudedata associated with the pulsatile displacement, or the pulsatile displacement, with reference data to identify presence of brain disorder in the human or animal subject; or (ii) performing the computer-implemented method of the third aspect, wherein the movable substance is in a human or animal subject, and comparing the measure of motion with reference motion measure to identify presence of brain disorder in the human or animal subject. In some examples, the reference data or the reference movement is determined based on biology of the human or animal subject such as age, gender, health status, size, weight, etc. Other features and aspects will become apparent by consideration of the detailed description and accompanying drawings. Any feature(s) described herein in relation to one aspect or embodiment may be combined with any other feature(s) described herein in relation to any other aspect or embodiment as appropriate and applicable.Term of degree such as “generally”, “about”, “substantially”, and the like, is used herein,depending on context, to account for one or more of: manufacture tolerance,degradation, trend, tendency, imperfect practical condition(s), etc. BRIEF DESCRIPTION OF THE DRAWINGS Embodiments will now be described, by way of example, with reference to the accompanying drawings, in which: Fig.1 is a block diagram illustrating a system in one embodiment; Fig.2 is a block diagram illustrating an MRI device in one embodiment; Fig.3 is a block diagram illustrating a data processing system in one embodiment; Fig. 4 is a flow diagram illustrating a method for detecting and / or quantifying motion based on magnetic resonance imaging in one embodiment; Fig.5 is a picture showing MRI images obtained from an experiment in one embodiment; Fig.6 is a picture showing MRI images obtained from the experiment in one embodiment; Fig.7 is a picture showing MRI images and graphs obtained from the experiment in one embodiment;

[0009] PC933674WO - specification Fig.8 shows a simplified MRI sequence and phase graph in one embodiment; Fig.9 shows simulation results obtained from a simulation in one embodiment; Fig.10 shows a simplified MRI sequence diagram in one embodiment; Fig.11A shows a simplified MRI sequence diagram in one embodiment; Fig.11B shows a simplified MRI sequence diagram in one embodiment; Fig.11C shows a simplified MRI sequence diagram in one embodiment; Fig.12A shows a simplified MRI sequence diagram in one embodiment; Fig.12B shows a simplified MRI sequence diagram in one embodiment; Fig.12C shows a simplified MRI sequence diagram in one embodiment; Fig.13 shows MRI images and graphs obtained from an experiment in one embodiment; Fig.14 shows a graph obtained from a simulation in one embodiment.Fig.15 shows a flow diagram illustrating a method for imaging and / or detecting pulsatiledisplacement using MRI in one embodiment;Fig.16 illustrates an MRI pulse sequence (simplified) and motion encoding effect in oneembodiment;Fig. 17 shows the simulation results of the magnitude and phase of the MR signal asfunction of time for a 1 Hz sinusoidal displacement function of increasing amplitude (forthe GRE and nb-SSFP versions of the 3D-EPI sequence in some embodiments);Fig.18 shows simulation results of the simulated nb-SSFP phase amplitude as a functionof displacement amplitude (from 0 to 1 mm) for grey matter (GM), white matter (WM),and CSF;Fig. 19 shows the simulation results of displacement sensitivity and relative phase-to-noise ratio for nb-SSFP FID and nb-SSFP echo in some embodiments;Fig. 20 shows the simulation results of CSF for sinusoidal displacement at ^^^ = 1^^^^ forthe EPG model in one embodiment versus Bloch simulations for GRE and nb-SSFP attwo selected displacement amplitudes;Fig. 21 shows in vivo measurement results of phase time-curves and correspondingfrequency spectra for single voxel ROIs located in the pons (graphs on the left) and the 4th ventricle (graphs on the right) respectively;

[0010] PC933674WO - specificationFigs.22A to 22I show in vivo measurement results of signal-pulsation for both nb-SSFPand GRE in single voxel ROIs after cardiac retrogating in selected areas and for differentdirections of the applied spoiler gradient: A- C) Repeatability of nb-SSFP phase-pulsationacross two scans in the same subject, wherein the same single voxel ROI location isused in both B and C; D-F) Comparison of nb-SSFP phase-pulsation values fromrepeated scans with varying spoiler gradient within the same subject and voxel ROI; G-I) Comparison of phase and magnitude pulsation in nb-SSFP versus GRE in subject 2;Fig.23 shows in vivo measurement results of phase-values at three different time pointswithin the cardiac cycle for the three different spoiler gradient directions;Fig.24 shows in vivo measurement results of displacement vectors at three different timepoints within the cardiac cycle; andFig.25 shows in vivo measurement results of phase-pulsation of the retrogated nb-SSFPsignal during the cardiac cycle measured in a representative single voxel ROI in the pons of different subjects. DETAILED DESCRIPTION Embodiments disclosed herein relate to imaging and / or detecting pulsatile displacement using MRI. Specifically, some embodiments disclosed herein provide MRI-based techniques for imaging and / or detecting displacement of fluids (e.g., CSF) and / or tissue, e.g., with sub-second time resolution. Fig.1 shows a system 100 in one embodiment. The system 100 includes an MRI device (MRI scanner) 102 and a data processing system 104 operably coupled with the MRI device 102 via one or more communication links 106. The one or more communication links 106 enable at least data communication between the MRI device 102 and the data processing system 104. The MRI device 102 may be a preclinical MRI device or a clinical MRI device. The data processing system 104 may be at least part of an operating console of the MRI device 102. In one example, the data processing system 104 is at least partly integrated with the MRI device 102. In another example, the data processing system 104 is separated from the MRI device 102. The one or more communication links 106 may include one or more wired communication links and / or one or more wireless communication links. The one or more communication links 106 may enable communication of MRI data acquiredby the MRI device 102 from the MRI device 102 to the data processing system 104 forprocessing. In one example, the one or more communication links 106 may enable communication of instructions or commands from the data processing system 104 to the MRI device 102 to affect its operation. Fig.2 shows an MRI device 200 in one embodiment. It should be noted that Fig.2 only illustrates the main components of the MRI device 200 to facilitate understanding of some embodiments. The MRI device 200 may be the MRI device 102 in Fig.1.PC933674WO - specification As shown in Fig. 2, the MRI device 200 includes a support structure 202, a magnet arrangement 204, a gradient coil arrangement 206, a radiofrequency (RF)coil arrangement 208, a power supply 210, a pulse sequence controller 212, and an RFtransmitter and receiver arrangement 214. The magnet arrangement 204 and the gradient coil arrangement 206 may be housed in a gantry of the MRI device 200. The support structure 202 may include a bed, a holder, or like support, for supporting asubject (e.g., animal, human) or object (e.g., phantom) to be imaged. The part of thesubject or object to be imaged may be referred to as imaging subject. The supportstructure 202 may be movable relative to the gantry, hence relative to the magnet arrangement 204 and / or the gradient coil arrangement 206, to facilitate imaging of the subject or object. The magnet arrangement 204 is arranged to provide a main, static magnetic field. The magnet arrangement 204 may include a cryomagnet. The strength of the main, static magnetic field may be, e.g., 1 Tesla, 3 Tesla, or 7 Tesla. In other examples, the main, static magnetic field may have a lower field strength or a higher field strength. The gradient coil arrangement 206 is arranged to generate directionally and / or temporally variable gradient magnetic fields in mutually perpendicular directions, e.g., inx-, y-, and z- directions. The gradient coil arrangement 206 may include three gradientcoils, each for generating gradient magnetic field in a respective direction.The RF coil arrangement 208 may include one or more RF coils. In one example, the RFcoil arrangement 208 includes an RF coil operable as both transmission and receivercoil. In another example, the RF coil arrangement 208 includes at least one transmissioncoil and at least one receiver coil. The RF coil arrangement 208 is operably connectedwith the RF transmitter and receiver arrangement 214. The RF coil arrangement 208 isarranged to receive RF signals from the RF transmitter and receiver arrangement 214 and to transmit RF pulses to a region of the subject or object to be imaged. The RFcoil arrangement 208 is further arranged to receive nuclear magnetic resonance (NMR)signals provided by the subject or object in response to the transmission of the RF pulses and to communicate the NMR signals to the RF transmitter and receiver arrangement214 for processing.The power supply 210 is arranged to provide power to the magnet arrangement 204 and / or the gradient coil arrangement 206 to operate them. The power supply 210 may independently provide power to the gradient coils of the gradient coil arrangement 206. The pulse sequence controller 212 is operably coupled with the power supply 210 and the RF transmitter and receiver arrangement 214, to control their operation hence to control the MRI pulse sequences used to image the subject or object. The pulse sequence controller 212 is arranged to control power provided to the magnet arrangement 204 and / or the gradient coil arrangement 206 (e.g., by controlling the intensity, timing, amplitude, direction, frequency, etc. of pulse current provided to the gradient coil arrangement 206) so as to control or affect the main magnetic field and / or the gradient magnetic fields. The pulse sequence controller 212 may individually orPC933674WO - specification independently control power provided to different gradient coils of the gradient coil arrangement 206. The pulse sequence controller 212 is arranged to control the RF signals transmitted from the RF transmitter and receiver arrangement 214 to the RFcoil arrangement 208. The pulse sequence controller 212 may further be operablycoupled with the support structure to control its movement relative to the gantry hence to control the part of the subject or object to be imaged. The RF transmitter and receiver arrangement 214 may include a transmitter and areceiver. The transmitter is arranged to provide RF signals to the RF coil arrangement208 based on the control of the pulse sequence controller 212. The receiver is arranged to process the received NMR signals to obtain MRI data. Fig.3 shows a data processing system 300 in one embodiment. The data processingsystem 300 is operable to process MRI data (also referred to as MR data) acquired byan MRI device. The data processing system 300 is operable to implement or facilitate implementation of the methods in some embodiments. The data processing system 300 may be the data processing system 104. The data processing system 300 may be at least partly integrated with an MRI device or it may be separated from the MRI device. The data processing system 300 includes components necessary to receive, store, and execute appropriate computer instructions, commands, and / or codes. In thisembodiment, the data processing system 300 includes at least one processor 302 andat least one memory 304. The at least one processor 302 may include one or more:CPU(s), MCU(s), GPU(s), logic circuit(s), Raspberry Pi chip(s), digital signal processor(s) (DSP), application-specific integrated circuit(s) (ASIC), field-programmable gate array(s) (FPGA), or any other digital or analog circuitry / circuitries configured to interpret and / or to execute program instructions and / or to process signals and / or information and / or data. The at least one memory 304 may include: one or more volatile memory (such as RAM, DRAM, SRAM, etc.), one or more non-volatile memory (such as ROM, PROM, EPROM, EEPROM, FRAM, MRAM, FLASH, SSD, NAND, NVDIMM, etc.), or any of their combinations. Appropriate computer instructions, commands, codes, models, information and / or data may be stored in the at least one memory 304. Computer instructions for executing or facilitating execution of the method embodiments may bestored in the at least one memory 304. The at least one processor 302 and the at leastone memory 304 may be integrated or separated but operably connected. Optionally, the data processing system 300 further includes at least one input device306. For example, the at least one input device 306 may include one or more of:keyboard, mouse, stylus, image scanner, microphone, tactile / touch input device (e.g., touch sensitive screen), image / video input device (e.g., camera), etc. Optionally, the data processing system 300 further includes at least one output devices 308. For example, the at least one output device 308 may include: display (e.g., monitor, screen, projector, etc.), speaker, headphone, earphone, printer, etc. The display may include a LCD display, a LED / OLED display, or other suitable display, which may or may not be touchsensitive. The data processing system 300 may further include at least one disk drive312 which may include one or more of: solid state drive, hard disk drive, optical drive, flash drive, magnetic tape drive, etc. A suitable operating system may be installed in thePC933674WO - specificationdata processing system 300, e.g., in the at least one disk drive 312 or in the at least onememory 304. The at least one memory 304 and the disk drive 312 may be operated by the at least one processor 302. Optionally, the data processing system 300 also includes at least one communication device 310 for establishing one or more communication links (not shown) with one or more other devices, such as servers, personal computers, terminals, tablets, phones, watches, or other computing devices, or MRI devices. The at least one communication device 310 may include one or more of: a modem, a Network Interface Card (NIC), an integrated network interface, a NFC transceiver, a ZigBee transceiver, a Wi-Fi transceiver, a Bluetooth® transceiver, a radio frequency transceiver,a cellular (e.g., 2G, 3G, 4G, 5G, beyond 5G, or the like) transceiver, an optical port, aninfrared port, a USB connection, or other wired or wireless communication interfaces. Transceiver may be implemented by one or more devices (integrated transmitter(s) and receiver(s), separate transmitter(s) and receiver(s), etc.). The communication link(s) may be wired or wireless for communicating commands, instructions, information and / or data.In one example, the at least one processor 302, the at least one memory 304 (optionallythe at least one input device 306, the at least one output device 308, the at least onecommunication device 310 and the at least one disk drive 312, if present) are connectedwith each other, directly or indirectly, through a bus, a Peripheral Component Interconnect (PCI), such as PCI Express, a Universal Serial Bus (USB), an optical bus,or other like bus structure. In one embodiment, at least some of these components maybe connected wirelessly, e.g., through a network, such as the Internet or a cloud computing network. It should be noted that the data processing system 300 is merely an example and that the data processing system 300 can have different configurations (e.g., include additional components, has fewer components, etc.) in other embodiments.It will also be appreciated that where the methods and systems disclosed herein areeither wholly implemented by computing system or partly implemented by computing system then any appropriate computing system architecture may be utilized. This may include stand-alone computers, network computers (cloud-based computers), or dedicated or non-dedicated hardware devices. Where the terms “computing system” and “computing device” are used, these terms are intended to include any appropriatearrangement of computer or processing hardware capable of implementing the functiondescribed. Fig.4 shows a method 400 for detecting and / or quantifying motion based on magnetic resonance imaging in one embodiment. As shown in Fig. 4, method 400 includes, in 402, receive MRI data associated with a movable substance. The MRI data is acquired by an MRI device, such as the MRI device 102, 200, using MRI pulse sequences operable to encode motion.Each of the MRI pulse sequences (each of the MRI pulse sequence corresponds to oneTR period) may be based on a non-balanced steady-state free precession technique andmay comprise a readout sequence for acquisition of a FID signal and / or an echo signal. The readout sequence may be arranged to image a volume of interest. For example, thePC933674WO - specification readout sequence may be an echo planar imaging sequence, such as a 3-dimensional echo planar imaging sequence. As another example, the readout sequence may be a FLASH-type imaging sequence, such as one with a golden angle readout to a volume- TR (total TR time for imaging the entire volume of interest). For example, the non- balanced steady-state free precession technique comprises gradient spoiling applied along a single direction. In one example, the gradient spoiling may be applied before, at, or near the beginning of the readout sequence. In another example, the gradient spoiling may be applied after, at, or near the end of the readout sequence. Each of the MRI pulse sequences may respectively include a combination of a readout sequence and a gradient spoiling applied along a direction (as a non-balanced steady-state free precession technique). In one example, the MRI pulse sequences include three types of MRI pulse sequences: a first type of MRI pulse sequence with a readout sequence and gradient spoiling applied along a first direction, a second type of MRI pulse sequence with a readout sequence and gradient spoiling applied along a second direction perpendicular to the first direction, and a third type of MRI pulse sequence with a readout sequence and gradient spoiling applied along a third direction perpendicular to each of the first direction and the second direction. Each of the three types of the MRI pulse sequence is for encoding motion in a respective direction. In one example, the first type of MRI pulse sequence is performed multiple times, then the second type of MRI pulse sequence is performed multiple times, and then the third type of MRI pulse sequence is performed multiple times. In one example, the MRI pulse sequences include multiple radiofrequencypulses for excitation (e.g., one radiofrequency pulse per TR period), and theradiofrequency pulses may have a linear or quadratic increment of phase ^^. As shown in Fig.4, method 400 also includes, in 404, process the MRI data based on a model, to detect and / or quantify motion associated with the movable substance. The model is arranged to associate phase information of MRI data with motion. The detection and / or quantification of motion may include the determination of one or more of: a direction of the motion, an amplitude of the motion, a relative amplitude of the motion, a velocity of the motion, and / or a relative velocity of the motion. In one embodiment, the model is associated with the MRI pulse sequences used to acquire the MRI data. For example, each of the MRI pulse sequences may respectively include a radiofrequency pulse for excitation and may define a repetition time ^^^^, and the model may be defined based on wherein ^^^ is transverse coherence state of order ^^, + indicates a transverse coherencestate immediately after application of a radiofrequency pulse in the MRI pulse sequence,− indicates a transverse coherence state immediately before application of theradiofrequency pulse in the MRI pulse sequence, ^^^^ is the repetition time of the MRIpulse sequence, ^^2 is transverse relaxation time of the movable substance, ^^ is aninteger indicating number of ^^^^-periods that the movable substance has experienced,and ^^^^ is a phase term that considers an effect of the motion associated with themovable substance on phase information of the received MRI data. In one embodiment,each of the MRI pulse sequences is respectively based on non-balanced steady-statePC933674WO - specification free precession technique and respectively comprises a readout sequence and gradient spoiling applied along a direction, and the phase term is defined as ^^^^^(^∙^^)where ^^(^^ ∙ ^^^^) at least includes a sum of of ^ ^⃗^ (^^ ∙ ^^^^) ∙ ∫^^ ^^⃗^(^^) ^^^^ and ^⃗^(^^ ∙ ^^^^) ∙ ^^ ^^^^ , where ^^ is a position (displacement) measure of the movable substance, ^⃗^is a velocity measure of the movable substance, ^⃗^ is the gradient spoiling applied in theMRI pulse sequence, and ^^ is time. In other embodiments, the phase term may bedefined using one or more equivalent equations. In one example, processing the received MRI data based on the model includes obtaining phase information from the received MRI data, and processing the phase information based on the model to determine a position (displacement) measure and / or a velocity measure of the movable substance. In one example, processing the received MRI data based on the model includes processing the received MRI data to reconstruct MRI images, and processing the MRI images based on the model to detect and / or quantify motion associated with the movable substance. In one embodiment of method 400, the motion may include pulsatile motion, which may be periodic or aperiodic, or slow motion, such as slow coherent motion. In one example,a slow motion may refer to a displacement in the order of 0.1 to 1 mm or a velocity ofsuch displacement over a duration of a cardiac cycle. In one embodiment of method 400,the movable substance may be present in a body of a subject, such as a human or an animal. The movable substance may include brain tissue, such as white matter, or body fluid, such as cerebrospinal fluid. The movable substance may be substance associated with the glymphatic system of the subject. In one embodiment of method 400, in which the movable substance is in the body of a subject, the processing of the received MRI data may be further based on physiological data of the subject, which is obtained from the subject while the MRI data is acquired. Examples of the physiological data include: data associated with one or more cardiac cycles of the subject, data associated with one or more respiratory cycles of the subject, etc. In one example, the processing of the MRI data in 404 includes temporally synchronizing or aligning the received MRI data with the physiological data of the subject and processing the temporally synchronized or aligned MRI data based on the model todetect and / or quantify motion associated with the movable substance. In other words,the processing of the MRI data in 404 may include the use of a retrospective gating technique. In one example, the MRI data may be acquired from the subject by the MRI device based on physiological data of the subject, i.e., using a prospective gating technique. By applying the prospective gating technique, the MRI data acquired by the MRI device is synchronized with the one or more physiological cycles (e.g., one or more cardiac orPC933674WO - specification respiratory cycles) of the subject. As a result, further temporal synchronization or alignment of the received MRI data with the physiological data in 404 may not be required. Preferably, in some embodiments, the MRI data acquired by the MRI device is non- contrast MRI data (i.e., the MRI data is obtained without using extrinsic MRI contrast agent). In method 400, operations 402 and 404 may be performed using a data processing system such as the data processing system 104, 300. In one embodiment, the method 400 may further include acquiring the MRI data using the MRI device. In this case, method 400 may be performed using both an MRI device, such as the MRI device 102, 200, and a data processing system, such as the data processing system 104, 300. In some embodiments, the method 400 may include further operations not specifically illustrated. For example, the method 400 may further include reconstructing MRI images based on the processing of the received MRI data and outputting the reconstructed MRI images for display. In this example, the operations of reconstructing MRI images and outputting them may be performed by the data processing system. For example, the method 400 may further include causing the MRI device to acquire the MRI data using the plurality of MRI pulse sequences. In this example, the operation of causing the MRI device to acquire the MRI data may be performed by the data processing system. Example studies related to method 400 have been performed. One example study relates to a method for characterization of CSF motion during cardiac cycle of a human subject. An object of this example study is to devise a fast and sensitive MRI method for measurement of CSF motion during the cardiac cycle. The method maybe considered as a specific example based at least partly on method 400 (in particularthe image acquisition part).There is a need to devise sensitive and / or robust MRI-based methods to characterizethe motion of CSF on both short and long timescales, e.g., to help advance the understanding of healthy and pathological CSF pulsation. Single-shot 3D volumetric imaging can be applied to study brain pulsation. 3D echo- planar imaging (EPI) based method, e.g., at 7T, may allow repeated acquisition of wholebrain images with sub-second temporal resolution, and it may have intrinsic sensitivity tomotion due to the non-zero zero- and first-order gradient-moments of the spoilergradients.3D EPI method may be suitable for probing CSF motion. In this example, an experiment is performed to evaluate the suitability of a 3D EPI based method for probing CSF motion. In the experiment, 3 healthy human subjects are scanned on a Siemens Magnetom Terra 7T MRI system (Siemens Healthineers, Erlangen, Germany) equipped with an 8TX / 32RX head coil (Nova Medical Inc, Wilmington, Delaware) using a gradient-echo 3D-EPI sequence. An example of the gradient-echo 3D-EPI sequence is disclosed in Stirnberg et al., “Segmented K-spacePC933674WO - specification blipped-controlled aliasing in parallel imaging for high spatiotemporal resolution EPI” (2021). In this experiment, the protocol is optimized for minimum volume-TR, with a 3D sagittally oriented matrix of 64 x 64 x 48, a partial Fourier (PF) factor of 6 / 8 along each of the phase -encoding direction and the partition encoding direction, and echo-spacing of 0.53 ms. The PF is placed at the end of the read-out to maximize TE within TR and maintain the first order gradient-moment at TE.3 x 2 CAIPIRINHA (controlled aliasing in parallelimaging results in higher acceleration) acceleration are applied, resulting in a TE of 6.7ms, a TR of 11.0 ms, and a volume-TR of 187 ms. In this experiment, a binomial water excitation (WE) pulse at flip angle (FA) = 10ois used for excitation, and 5003D volumes are acquired (imaged) with a total scan time of 93.5 seconds. During imaging, physiological data of the human subjects are detected using a respiration belt and a pulse oximeter and are recorded. The scan is repeated at different phase encoding (PE)directions. For one of the subjects, the same scan is repeated twice within the samesession. For anatomical reference, Magnetization-Prepared Rapid Acquisition Gradient Echo (MPRAGE) images at 0.75 mm isotropic resolution is acquired and brain-masked using software tool SPM12. Software tool FSL-Feat is used for motion correction and spatial smoothing of 3D-EPI, and the first 20 volumes are discarded to ensure the MRI data reflects a steady-state. Echo planar image correction (EPIC) is applied for EPI distortion correction. In this experiment, retrospective gating is applied: the peak slope of the ascending ridge of each cardiac pulse is defined as the start of the cardiac cycle,and for each of the 3D volume acquired (imaged), the time delay from the latest triggerpoint is found. A cine-video with 20 time bins is created by averaging across all image volumes within each time bin, resulting in around 20 images in each bin after discarding those falling outside a subject-specific predefined RR-interval (the time elapsed between two successive R waves of the QRS signal / wave on the electrocardiogram). In each voxel, the magnitude and phase oscillations are calculated as percentage / degrees change from the mean across the cardiac cycle in each voxel. For region of interest (ROI) analysis, single voxel ROIs in the 3D-EPI images are defined, with visual guidance from MPRAGE. Figs.5 to 7 show some results of this experiment. Fig. 5 shows mean 3D-EPI image and temporal standard deviation (STD) in both magnitude and phase for the four different acquisitions in one subject. From Fig.5, it can be seen that the temporal STD is similar across all scans (but some spatial differences can be observed). Residual geometric distortion is visible even after EPI distortion correction using EPIC. This may be due to relatively low tissue contrast available in this example. Fig. 6 shows the magnitude and phase images. The magnitude images show the difference in the lateral ventricle dynamics between the Head-Feat (HF) phase encoding (PE) direction and Posterior-Anterior (PA) phase encoding (PE) direction. The images show the sensitivity to directional flow in the MR signal magnitude. Signal fluctuations also observed in cortical subarachnoid space, both in magnitude and phase images.Phase oscillations in the order of ±2 degrees as observed in ventricles corresponds to aPC933674WO - specificationvelocity of ±1 cm / s, generally consistent with literature values disclosed in Matsumae etal., “Changing the Currently Held Concept of Cerebrospinal Fluid Dynamics Based on Shared Findings of Cerebrospinal Fluid Motion in the Cranial Cavity Using Various Types of Magnetic Resonance Imaging Techniques” (2019) and Yatsushiro et al., “Cardiac- driven Pulsatile Motion of Intracranial Cerebrospinal Fluid Visualized Based on a Correlation Mapping Technique” (2018). A global phase variation over the cardiac cycle is also consistently observed. Fig. 7 shows temporal evolution of magnitude and phase across the cardiac cycle for some example regions of interest. in Fig. 7, Graph A) illustrates the difference time- evolution that is typically observed between magnitude and phase; Graph B) shows the strong dependence in magnitude on the PE direction; Graph C) illustrates the variation across the three subjects (note that subject 1 had higher heart rate, 76 bpm, than the two other subjects (54 bpm and 60 bpm)); and Graph D) illustrates both the relatively low magnitude oscillations from a medium-sized artery and the high repeatability between scans.From example cine-videos obtained for one of the subjects at the two different PE-directions, significant fluctuations in both magnitude and phase across the cardiac cycle are observed in CSF and blood compartments. Especially in the lateral ventricles, there is a significant difference between the Head-to-Feat versus Posterior-to-Anterior phase encoding directions (also see graph B) in Fig.7). Significant CSF dynamics are detected even in the cortical subarachnoid space as illustrated in graph A) in Fig.7, where it can also be observed that magnitude and phase fluctuate differently. The three subjects demonstrate significant variation as illustrated for a ROI in the Right Lateral Ventricle in graph C) in Fig.7. The within-session repeatability is illustrated in graph D) in Fig.7 for a ROI placed in the pericallosal artery, which also illustrates the fairly low magnitude fluctuations for this medium-sized artery.Based on the results obtained in this example study, it can be observed that fast 3D-EPIis sensitive to fluid motion in the brain during the cardiac cycle. This sensitivity can partly be explained by the relatively high first-order gradient moment at TE along the PE direction, corresponding to an encoding velocity Vencof 1.0 m / s for the protocol settings in this example. Preliminary results from simulations using an extended phase graph model (e.g., the model described with reference to method 400) indicate that contributions from higher order pathways may also play a role, due to the long T2 of around 2.0 s for CSF. It is envisaged that an extended phase graph model can be used to quantify CSF motion with 3D-EPI. This example study demonstrates that dynamic 3D-EPI at low spatial resolution and short volume TR is very sensitive to CSF dynamics and can be optimized for characterization of CSF motion. This represents a valuable extension of the portfolio of fast MRI methods aimed at studying brain pulsation. In this example, brain 3D-EPI at 3 mm resolution and volume-TR = 187ms is acquired at 7T for 94 seconds and retrospectively sorted into 20 cardiac phases based on pulse oximeter. CSF dynamics is observed in all ventricles as well as in subarachnoid space, in addition to arterial pulsation. Both magnitude and phase pulsations are present in the 3 subjects scanned.PC933674WO - specification Highest sensitivity to motion is observed along the phase-encoding direction. This example study shows that fast 3D-EPI at 7T is very sensitive to the motion of CSF during the cardiac cycle and that it can be used to characterize CSF dynamics on individual patient level and aid understanding of brain diseases. Various MRI based motion encoding methods (sequence) exist. The inventors have found that short-TR / low flip angle sequence, such as spoiled gradient echo (GRE) sequence and non-balanced steady-state free precession (nb-SSFP), have intrinsicposition encoding through the spoiler gradient within each TR period, and the resultingdisplacement-dependent phase for each coherence state may build up coherently orincoherently over time. The inventors have found that in this field of art, sensitivity tomotion in nb-SSFP has been treated as a problem causing image artefacts and has notbeen considered to provide any useful image contrast. The inventors have devised,against this prejudice in this field of art, motion sensitivity of nb-SSFP may be exploitedto obtain useful motion-related information. The inventors have further discovered thatdue to the long lifetimes of the longitudinal states the gradient encoding becomes extremely sensitive to small periodic displacements, which can provide a useful signal to detect motion such as flow.Fig.8 shows an example nb-SSFP sequence (without readout) and phase graph in oneembodiment. The nb-SSFP sequence may be considered as part of the nb-SSFP technique in the method 400. As shown in Fig.8, the nb-SSFP sequence includes, ineach period, a radiofrequency pulse ⍺^^ (^^ = ^^, ^^ + 1, … ) followed by a spoiler gradientapplied in a direction. Depending on the spoiler gradient, the NMR signal may be an echo signal or an FID signal, which can be acquired to obtain related phase information.In one example, the flip angle ⍺ and / or the RF phase ^^ may be constant. In one example,the flip angle ⍺ and / or the RF phase ^^ may be varied (e.g. linearly or quadraticallyincrease for spoiled gradient echo FID imaging) to alter displacement and motion sensitivity. In one example, the radiofrequency pulses have a linear or quadratic increment of phase ^^.The inventors have devised that highly accelerated GRE 3D-EPI can be applied studybrain pulsation. In one example, to resolve the cardiac frequency sufficiently, the highly accelerated GRE 3D-EPI can provide a temporal resolution of below 0.2 seconds and a spatial resolution of around 3 mm isotropic. Such an example has been described with reference to Figs. 5 to 7. In that example, in addition to the expected magnitude fluctuations in the tissue, large high-frequency fluctuations in both magnitude and phase are observed, and these fluctuations are different for different phase-encoding directions.Based on the example described with reference to Figs. 5 to 7, the inventors havemodified the original 3D-EPI from a spoiled GRE to a non-balanced SSFP, and have focused on the phase of the signal. The inventors have devised that an extended phase graph model (e.g., the model described with reference to Fig.4) may be used to associate phase information of MRIdata with motion. In one embodiment, the model is defined by:^^^^^^ (^^ ∙ ^^^^) = ^^^^^^^^^⁄ ^^^^^^^(^∙^^)PC933674WO - specificationwhere ^^(^^ ∙ ^^^^) = ^ ^⃗^ (^^ ∙ ^^^^) ∙ ∫^^ ^^⃗^(^^) ^^^^ + ^⃗^(^^ ∙ ^^^^) ∙ ∫^^ ^^⃗^(^^) ^^ ^^^^ + …Meanings of the symbols have been defined in the model described with reference toFig. 4 so will not be repeated here. In this example, the phase term isdetermined by Taylor expansion. In this example, the phase term is included to account for the propagation of transverse states within each TR, to better model the phase-effect from pulsatile motion in the extended phase graph model. In this example, the phase term contains the phase value induced by the spoiler gradient and the current position x. Different pathways that contribute to the F-states have experienced different pairs ofdephase-rephase gradients at different time points and positions. The model in thisexample has been validated against Bloch simulations. The extended phase graph model in this above example has been used to simulate pulsation motion. In the simulation, the extended phase graph model is applied, where an oscillating displacement is introduced as an additional phase-term for all F-states during each TR. The value of this phase-term varies from one TR to the next TR depending on the details of the oscillation, and the resulting F0-state will be a complex summation of states with varying phase-term. For simplicity, as shown in Fig.9, a sinusoidal oscillation at 1 Hzand amplitudes of 10 microns, 0.1 mm, and 1 mm are used.Fig.9 shows example simulation results for the phase of the MR signal to illustrate the different regimes identified. In the graphs in Fig.9, solid lines represent CSF and dotted lines represent white matter. In the graphs in Fig.9, the upper row of graphs relates to the spoiled GRE sequence, the lower row of graphs relates to the nb-SSFP sequence (i.e., the GRE sequencewithout applying RF-spoiling). As shown in the left-most graph in the upper row with 10-micron displacement amplitude in the spoiled GRE sequence, there is a regime for CSF which provides extremely high sensitivity for motion. However, for larger displacement amplitudes as shown in the right-most graph in the upper row, the sequence is in a spoiled state with pseudo-random phase values over time. In between as shown in the middle graph in the upper row there is a transition regime. In the graphs of the upper row, no effect of displacement on white matter can be observed, due to much shorter T1 and T2 of white matter. As shown in the graphs in the lower row, a coherent regime extends to larger displacement amplitudes, and displacements in white matter also affects the phase in a coherent way.The inventors have devised that operation in a regime which is sensitive to both fluid(e.g., CSF) and tissue (e.g. white matter WM) motion and where the phase-oscillations scales linearly with the displacement amplitude may be preferred. To this end, theinventors have optimized the MRI imaging protocol (sequences).Fig.10 shows a simplified illustration of an optimized MRI imaging protocol (nb-SSFP sequences) in one embodiment. The optimized MRI imaging protocol shown can be used in the method 400 to obtain the MRI data, which can then be processed using the modelPC933674WO - specification to detect and / or quantify motion. Fig.10 only shows one sequence (period). As shown in Fig.10, the optimized MRI imaging protocol includes a (balanced, interleaved multi- shot) 3D-EPI imaging sequence and a spoiler gradient applied after the 3D-EPI imaging sequence to acquire the FID signal pathway. It should be noted that while Fig.10 showssimultaneously show three spoiler gradients, only one spoiler gradient is applied along asingle direction in each period. In other words, the optimized MRI imaging protocol in this example includes one type of sequence with spoiler gradient applied along one direction, one type of sequence with spoiler gradient applied along another direction, and one type of sequence with spoiler gradient applied along yet another direction (these directions are mutually perpendicular). In this example, the imaging protocol is for 7T. The imaging protocol (Fig.10) is based on nb-SSFP with option to modify spoiler gradient. The imaging protocol can be considered as nb-SSFP 3D-EPI, where nb-SSFP technique is applied along with a 3D- EPI readout. The nb-SSFP technique includes applying spoiler gradient along a singledirection at a time. For each sequence (period) of the protocol, the spoiler gradient isapplied only in one direction (although three are shown simultaneously overlaid on the 3D-EPI readout sequence in Fig.10). In this example, the spoiler gradient is applied at the end of the 3D-EPI readout sequence. However, in another example, the spoiler gradient is applied at the beginning of the 3D-EPI readout sequence. In this example,the imaging protocol further includes the following parameters: 60 x 60 x 40 matrix – 3.2mm isotropic voxels; sagittal orientation; phase encoding direction is anterior-posterior (A-P) direction; a non-selective water excitation flip angle of 20o, volume-TR = 148 ms (6.7 frames per second), 2π spoiler gradient along an axis, 500 volumes free breathing or 80 volumes during breath-hold, sequence perform at least 3 times with spoiler gradient applied along three different axes (perpendicular to each other); total scan time of 1-4 minutes. Figs.11A to 11C show MRI pulse sequences in one embodiment, which can be used to in the method 400 to obtain the MRI data, which can then be processed using the model to detect and / or quantify motion. The MRI pulse sequences in Figs.11A to 11C are to acquire FID signals (which contain phase information). Each of Figs.11A to 11C only show one sequence or period. In this example, the imaging (readout) sequence is balanced (e.g. one shot of an interleaved multi-shot 3D-EPI readout) such that the spoiler gradients define the amount and direction of dephasing at the end of TR following signal readout. In one example, the MRI pulse sequences used to obtain the MRI data include multiple applications of the sequence in Fig. 11A (to encode motion in one direction), followed by multiple applications of the sequence in Fig. 11B (to encode motion in another direction), and then followed by multiple applications of the sequence in Fig.11A (to encode motion in yet another direction). Figs.12A to 12C show MRI pulse sequences in one embodiment, which can be used to in the method 400 to obtain the MRI data, which can then be processed using the model to detect and / or quantify motion. The MRI pulse sequences in Figs.12A to 12C are to acquire echo signals (which contain phase information). Each of Figs.12A to 12C only show one sequence or period. In this example, the imaging (readout) sequence itself is balanced (e.g. one shot of an interleaved multi-shot 3D-EPI readout) such that the spoilerPC933674WO - specification gradients define the amount and direction of dephasing at the beginning of TR before signal readout. In one example, the MRI pulse sequences used to obtain the MRI data include multiple applications of the sequence in Fig. 12A (to encode motion in one direction), followed by multiple applications of the sequence in Fig. 12B (to encode motion in another direction), and then followed by multiple applications of the sequence in Fig.12A (to encode motion in yet another direction). Fig.13 shows the MRI data obtained using the imaging protocol of Fig.10. The MRI data shown is MRI phase data obtained after cleaning up (unwrapping and voxel-wise subtraction of mean). The three images on the left show spoiler gradient applied along different directions (Head-Feat (H-F) direction, Anterior-Posterior (A-P) direction, Left- Right (L-R) direction). Based on the MRI data, it has been found that CSF and tissue pulsation can be observed in real-time and that the three different acquisitions show different signal behaviour. The three graphs on the right correspond to a one-minute plot of the phase values in a single voxel placed in the pons (right in front of the brainstem) and in left lateral ventricle. As shown in these graphs, a high SNR (not noise, but actual phase fluctuations in a single voxel) is obtained, fluctuations in both tissue and CSF can be observed, fluctuations associated with CSF are on a much longer timescale than the cardiac cycle. In one example, the MRI data obtained is further processed by retrogating against a pulse oximeter measuring and creating a CINE video with 20 bins across the cardiac cycle. In addition, the different displacement encoding acquisitions are combined into displacement vectors. It has been found that CSF is pushed down from the lateralventricles via the cerebral aqueduct and down in the spinal canal during systole, only toreturn in the opposite direction during diastole. In one example, scaling factors from the extended phase graph model simulations are used to estimate the displacement amplitudes in a few different regions. Fig.14 shows a graph indicating the scaling factors (to convert between peak phase obtained in the MRI data and displacement associated with the motion). In other words, by using the phase information in the MRI data, and based on the extended phase graph model, motion can be detected and quantified. Example values obtained from single subject, individual voxel ROI on the pons include a peak phase of 0.18 rad and an estimated displacement of 0.17 mm (obtained based on the scaling factors). The estimated displacement of 0.17 mm appears to agree well with an example literature value range of 0.15mm to 0.18 mm.Through some of the simulations and experiments above, it has been found that 3D-EPInon-balanced SSFP FID (among other protocols) can be used to effectively imagepulsatile motion (e.g., displacement and / or velocity) in both CSF and tissue. The imagingof the pulsatile motion may be performed in real-time. It has also been found that the devised extended phase graph model with a phase term incorporated agrees withexperimental observations. The method with the devised model can be used to detectand quantify motion, and it is highly sensitive to small displacements. It is envisaged thatPC933674WO - specification the method can be used for characterization of brain pulsation, e.g., in the clinic, due to short scan time and high SNR. Some embodiments provide a method for imaging pulsatile displacement of tissue and fluid in the human brain and central nervous system with the use of MRI. In one example, the method uses a non-balanced steady-state free precession sequence with spoiler gradient applied along a single physical direction at a time, combined with a highly accelerated 3D echo-planar imaging read-out to acquire repeated images over the whole brain with a time-resolution of around 5 frames per second. Multiple runs with the spoiler gradient applied along different directions are required to get directional information. Total acquisition time ranges from 30 seconds to 5 minutes. Resulting magnitude and phase images are then sorted based on additional physiological recordings to create 4D time-curves in 3D within chosen periodic interval, such as the cardiac or respiratory cycle. Based on simulations using extended phase graph model (the model described with reference to Fig. 4), the displacement amplitudes can be quantified in physical units. Compared to some existing methods, this method has better combined spatial coverage and temporal resolution as well as higher sensitivity to small displacements. The method in this example is sensitive to displacements down to around 10 microns at a spatial resolution of around 3 mm. As pulsatile motion of brain tissue and neurofluids play important roles in brain health, the method in some embodiments can be used for basic neurological research and forclinical diagnostic purposes. The simplicity and robustness of the method may make itparticularly attractive for use in patients.Fig. 15 shows a method 1500 related to detecting and / or quantifying motion based onMRI in one embodiment. Specifically, in this embodiment, the method 1500 is for imagingand / or detecting pulsatile displacement using MRI. For example, the pulsatiledisplacement may be associated with physiological motion such as cardiac motion and / or respiratory motion. For example, the pulsatile displacement may be pulsatile displacement of a biological tissue and / or a biological fluid. For example, the pulsatile displacement may be pulsatile displacement of at least one of: brain tissue (such as white matter) or CSF.The method 1500 includes, in 1502, apply an MRI pulse sequence with a 2D- or 3D-based readout to an imaging subject. The MRI pulse sequence includes multiplerepetitions (TR) each including displacement encoding gradient for sensitization topulsatile displacement and each having the 2D- or 3D- based readout applied eitherbefore or after the displacement encoding gradient. The displacement encoding gradient includes a gradient lobe with a finite zero-order moment. In some implementations, the zero-order moment of the gradient lobe may have a magnitude of at least 2π / γ∆x, where ∆x is voxel size of the MR data and γ is gyromagnetic ratio. In some implementations, ifthe 2D- or 3D- based readout is applied before the displacement encoding gradient, the2D- or 3D- based readout may acquire a FID signal; and if the 2D- or 3D- based readoutis applied after the displacement encoding gradient, the 2D- or 3D- based readout mayacquire an echo signal. Examples of the 2D- or 3D- based readout include: MREG-typePC933674WO - specification 3D readout (single-shot, non-Cartesian per volume, repeated over multiple TRs), 2D-spiral readout, 3D-spiral readout, etc. In some implementations, the 2D- or 3D- basedreadout is a 3D EPI based readout such as a segmented 3D EPI readout, and the MR data is volumetric MR data. In some implementations, for each of the repetitions (TR),the MRI pulse sequence lacks velocity encoding gradients before the 2D- or 3D- basedreadout. In some implementations, the displacement encoding gradients of the plurality of repetitions (TR) may include orthogonal gradients for sensitization to pulsatile displacements along orthogonal directions. For example, the repetitions (TR) may have different displacement encoding gradients for sensitization to pulsatile displacement along different directions. In some implementations, the MRI pulse sequence may apply RF spoiling, e.g., for each of the repetitions. In some implementations, the MRI pulse sequence may not apply any RF spoiling for any of the repetitions. In some implementations, the MRI pulse sequence is a non-balanced steady-state free precession (nb-SSFP) sequence. In some implementations, the MRI pulse sequence is a gradient-echo (GRE) based sequence, e.g., a spoiled GRE sequence. The method 1500 includes, in 1504, acquire MR data of the imaging subject. In some implementations, the imaging subject may include a body part, such as a brain, of a subject, such as a human or animal, and the acquisition of the MR data may be based at least in part on prospective ECG or respiratory gating associated with the subject. The method 1500 includes, in 1506, process the MR data to obtain phase and / or magnitude data associated with the pulsatile displacement. In some implementations, the imaging subject may include a body part, such as a brain, of a subject, such as a human or animal, and the processing of the MR data may be based at least in part on retrospective ECG or respiratory gating associated with the subject. The method 1500 may be used to image and / or detect pulsatile displacement which may have an order of magnitude of 10-3mm, an order of magnitude of 10-2mm, an order ofmagnitude of 10-1 mm, or an order of magnitude of 1 mm. The method 1500 may be usedto image and / or detect pulsatile displacement which may have a temporal resolution inan order of 1 s to 10-1 s.In some embodiments, the method 1500 may include operations not specificallyillustrated. For example, the method 1500 may include converting phase data todisplacement data associated with the pulsatile displacement based at least in part on alookup or mapping table which includes information for associating phase withdisplacement. For example, the method 1500 may include converting phase data todisplacement data associated with the pulsatile displacement based at least in part on amodel such as a Bloch Simulation model or an EPG model. The model can associatephase with displacement. The phase, magnitude, and displacement data may beprocessed for display as MRI image(s).In some implementations, the method 1500 may be performed using the system 100. Insome implementations, the imaging aspect of the method 1500 may be performed usingthe MRI device 200. In some implementations, the data processing aspect of the method1500 may be performed using the data processing system 300.PC933674WO - specificationTurning now to MR theory associated with the method 1500, an extended phase graphmodel for low-amplitude, pulsatile motion in one embodiment is provided.In general, the EPG model may allow efficient simulation of the formation of various MRcontrasts, e.g., short-TR sequences such as steady-state free precession with or withoutRF- and gradient-spoiling. In one embodiment, this formalism is extended to incorporatethe effect of low amplitude, pulsatile motion. For a standard EPG model for a non-balanced steady state free precession (nb-SSPF) sequence, where the propagation of the transverse state ^^^^ (^^ − 1) immediately afterRF-pulse ^^ − 1 to the transverse state ^^ ^^^^ (^^) immediately before the next RF-pulse ^^can be represented as: where ^^ is the order of the state, ^^^ represents the FID, and ^^^^ represents the echo thatrefocuses towards TR.An associated MRI pulse sequence in one embodiment is illustrated (in simplified form)in Fig.16.The MRI pulse sequence in Fig. 16 can be considered as a short-TR sequence withspoiler gradients. In Fig.16, the ^^^ state created by RF-pulse number ^^ − 1 may acquirea net phase change if the isochromat position is changed along the direction of thespoiler gradient between the two gradient lobes within TR-periods ^^ − 1 and ^^. The totalphase of ^^^states may be determined by the motion experienced by all individual pathways contributing to the given state.In some implementations, for each TR, a 2D- or 3D- based readout sequence, such asa 3D EPI based readout, may be applied before the gradient to readout the FID signal.In some implementations, for each TR, a2D- or 3D- based readout sequence, such as3D EPI based readout, may be applied after the gradient to readout the echo signal.The above formalism may apply when the zeroth order gradient moment along all threeorthogonal axes is finite (not zero) and identical for all TR-periods. This may require rewinding of the spatial encoding gradients and the addition of a spoiler gradient. Thespoiler gradient may result in the lifting of all transverse states by one index ^^. Successfulseparation of the transverse states may be achieved if the k-spaces of the states do notoverlap, or if the spoiler gradient results in at least 2^^ phase-shift across a single voxel.This may provide the following requirement for the spoiler gradient moment ^^^^: where Δ^^ is the voxel size and ^^ is the gyromagnetic ratio.The EPG model, in its standard form, does not include the effect of motion or flow. In thisembodiment, the effect of motion or flow is included in the EPG model by the addition ofa phase-term ^^^^(^) in equation (1), e.g.:PC933674WO - specificationwhere ^^(^^) may be expressed as: where ^^^^(^^^^) is the jth-order gradient moment vector as integrated over one TR: and^^^^(^)^^^is the jth-order time-derivative of the position vector ^^ evaluated at time ^^ = ^^ ⋅^^^^. The number of required terms in the Taylor expansion may depend on the choice ofMR protocol and / or the expected motion regime.The gradient moments for a 3D-EPI sequence applied in the example and thecorresponding phase-terms up to 4thorder for a sinusoidal pulsatile displacement of 1 mm amplitude and 1 second period have been calculated. The calculation results areshown in Table 1. Specifically, Table 1 shows the gradient moments at the end of eachTR for 3D-EPI sequence with corresponding maximum contribution to phase term inequation (4) for sinusoidal displacement function with amplitude ^^^ = 1 mm at ^^ = 1 Hz.Maximum velocity becomes 2^^^^^and maximum acceleration becomes 4^^^^^^and soforth. Smaller displacements lead to smaller values. For the regime of motion in thisexample, two terms may be sufficient.Table 1 – Calculated gradient moments in one exampleGradient Moment Order Value ^^^^[^^^^^^]Zeroth Order (position) 7.3e − 6 2.0First Order (velocity) 1.1e − 7 0.2Second Order (acceleration) 1.1e − 9 0.01Third Order (jerk) 1.2e − 11 0.001For this MR protocol, in a regime of smooth motion, the Taylor expansion may be limitedto the two first terms: ^^(^^) ≃ ^^[^^^^(^^^^) ⋅ [^^(^^) − ^^(0)] + ^^^^(^^^^) ⋅ ^^(^^)] (6)where ^^^^ and ^^^^ are the zero and first order gradient moment vectors, and ^^ is thevelocity vector.The FID signal at TE may be given by the zero-order pathway ^^^^, propagated in thesame way from RF-pulse ^^ : where ^^(^^^^, ^^) ≃ ^^[^^^^(^^^^) ⋅ [^^(^^) − ^^(0)] + ^^^^(^^^^) ⋅ ^^(^^)] (8)In this example a left-rotating coordinate system is assumed. The small change in ^^ and^^ from TE to TR is not taken into account in this formulation.PC933674WO - specificationSimulations are performed based on the EPG model for low-amplitude, pulsatile motion(EPG model with phase term) in the above embodiment.In the simulations, the EPG model is implemented in Matlab R2023B (MathWorks) andis used for GRE and nb-SSFP simulations. An overview of the input parameters for thesimulations is shown in Table 2.Table 2 – Input parameters used for simulations of a EPG model for low-amplitude,pulsatile motion in one embodimentParameter Name Symbol Unit ValuesRepetition Time TR ms 5,10,20,40,100Echo Time TE ms 0Flip Angle FA deg 5, 7, 10, 15, 20, 30, 40, 60, 90Gradient Spoiler M0 Rad 2^^Spatial Resolution Δ^^ mm 1, 3, 50th Gradient Moment M0 ^^^^ / ^^ ⋅ ^^^^ 4.7, 7.8, 23.51th Gradient Moment M1 ^^^^ / ^^ ⋅ ^^^^^ 47, 78, 235RF-spoil quadratic increment deg 0, 50Longitudinal Relaxation Time T1 sec 1.6, 4.0Transverse Relaxation Time T2 ms 50, 1000Displacement Amplitude ^^^ mm 0.001-5Displacement Frequency ^^^ Hz 0.5, 1, 4The T1 and T2 relaxation times used in the simulations generally correspond to those ofwhite matter (WM), grey matter (GM), and cerebrospinal fluid (CSF) at 7T. In thesimulations, for simplicity, TE is set to 0 (thus the effect of ^^2∗ is not accounted for in thesimulations). In this example, to simulate cardiac and respiratory driven pulsatile motion, a sinusoidaldisplacement function is assumed along a single physical direction ^^ :Δ^^(^^) = ^^^sin (2^^^^^^^) (9)where ^^^ is the amplitude of the pulsatile displacement, ^^^ is the frequency and ^^ is time(^^ = 0 at first RF pulse). Both ^^^ and ^^^ are varied in the simulations (Table 2). Thevelocity function ^^^(^^) along direction ^^ may be defined as the time-derivative of Δ^^(^^):^^^(^^) = 2^^^^^^^^cos (2^^^^^^^) (10)In the simulations, Δ^^(^^) and ^^^(^^) are only evaluated at times ^^ = ^^^^^^ and areassumed to be constant between RF pulses.PC933674WO - specification In relation to the simulation parameter settings of the simulations, due to the large totalparameter space (large amount of parameters), not all combinations of input parametersare explored in detail in the simulations. In this example, the combination of TR = 10 ms,FA = 20 degrees, and voxel size Δ^^ = 3 mm is used as default sequence settings incases where the effect of different displacement functions are explored, as well as for comparison of GRE versus nb-SSFP and EPG versus isochromat averaging models. Correspondingly, in this example, the amplitude and frequency of the displacementfunction are set fixed at ^^^ = 1 mm and ^^^ = 1 Hz, respectively, when the values of TR,FA and Δ^^ are varied. In all cases of the simulations, the minimum required 0th orderspoiler gradient moment ^^0 is applied according to equation (2) while the 1st ordermoment M1 is estimated assuming the spoiler gradient is applied at the very end of theTR-period and with Gmax = 25mT / m.In this example, after each EPG simulation, the phase of both the FID (^^^^state) and the echo (^^^^^ state) signals are extracted for the last full displacement period (1 / ^^^) of thesimulation. From this data, the peak amplitude ^^^ of the phase variation is measured. Inaddition, the time-shift Δ^^ between the displacement function and the phase function isestimated using cross-correlation, and the smoothness of the phase function is estimated using the root-mean squared error (RMSE) between the actual phase function and a best-fitted sinusoidal function. The displacement sensitivity is defined as ^^^ / ^^^and is measured in units of [Rad / mm].The relative phase-to-noise ratio is also calculated, by multiplying the sensitivity with themean magnitude of the FID and echo signals respectively.For default sequence parameters, the EPG model in the above embodiment is comparedto corresponding Bloch simulations of a single 3.0 mm voxel in isocenter, with 301 isochromats positioned along z separated by 0.01 mm. Specifically, in this example, atimestep of 5μs is used, and a simple rectangular RF-pulse together with the actualgradient shape of the spoiler gradient in the 3D-EPI sequence are applied. Thesimulation runs for at least 10 seconds (= 1000 TR intervals) to allow T1 effects tostabilize. For calculation of the effective B-field and rotation matrix for each time step,the core code from Okell Thomas W. Bloch Simulator for MRI pulse sequences: Initialrelease, Zenodo, 2023 is used.In vivo imaging experiments are also performed.For the in vivo imaging experiments, four healthy volunteers (aged 23-28 years, 2 maleand 2 female) are scanned on a Siemens Magnetom Terra 7T MR system (SiemensHealthineers, Erlangen, Germany) that is equipped with a 8TX / 32RX head coil (Nova Medical Inc, Wilmington, MA) using a custom 3D-EPI sequence with rewinding of spatial- encoding gradients. Details pertaining to the custom 3D-EPI sequence can be found inR Stirnberg and T Stöcker, “Segmented K-space blipped-controlled aliasing in parallelimaging for high spatiotemporal resolution EPI”, Magnetic resonance in medicine vol.85,3 (2021): 1540-1551.PC933674WO - specificationThe MR sequence has the option to vary both the RF-spoiling between the standard 50degrees quadratic RF-phase increment (= spoiled GRE) and no RF-spoiling (= non-balanced steady-state free precession). Both are experimented. Furthermore, thegradient spoiler moments that are added to the rewinding moments can be individuallyadjusted on all three axes. The basic protocol is optimized for minimum volume TR, witha 3D sagittally oriented matrix of 60×60×40 with phase-encoding along the posterior-anterior direction and nominal isotropic voxel size of 3.2 mm. Partial Fourier = 6 / 8 isapplied along both the phase-encoding and the 3D direction, and the echo-spacing is 0.53 ms. In this example, the Partial-Fourier-skipped phase-encoding lines are placed at the beginning of the EPI read-out to minimize TE within TR and to limit the first ordergradient-moment at TE. 3x2 blipped-CAIPI acceleration is applied, resulting in TE / TR =4.05 / 10.6 ms and volume TR = 148.4 ms. A non-selective binomial 1-2-1 water excitationRF-pulse at nominal FA = 20 degrees is used for excitation. 500 volumes are acquiredwith a total scan time of 74.2 seconds. The reference lines required for online parallel imaging reconstruction (GRAPPA) are acquired at the beginning of each scan using afast spoiled GRE readout. The scan is repeated for different values of gradient spoilermoments and directions, as well as with RF-spoiling enabled and disabled, while otherparameters are kept unchanged. For each subject, the same scan is repeated twicewithin the same session. Physiological data are logged using a respiration belt and apulse oximeter mounted on the index finger of the subject. For anatomical reference, anMPRAGE with universal pulses at 1 mm isotropic resolution is acquired.In respect of image processing for the in vivo experiments, FSL Feat is used for motioncorrection of each 3D magnitude image volume, the first 30 volumes are deleted tostabilize the T1 effects. No spatial smoothing or temporal high-pass filters are applied inthis example.The phase images are pre-processed. First, real and imaginary 3D images are calculatedfrom the reconstructed magnitude and phase images. Then, transforms based on themotion estimates from FSL Feat are applied to the real and imaginary images separately. Each motion-corrected complex image volume is subsequently multiplied with thecomplex conjugate of the mean across all 3D volumes and, from the resulting Hermitianinner product, new net phase images are calculated. In addition, for the phase imagesto be used in retrospective cardiac gating, a high-pass filter with cross-over frequency at0.35 Hz is applied in each voxel.For retrospective gating based on the pulse oximeter data, the peak slope of theascending ridge of each cardiac pulse is identified and used as the reference time point.For each acquired 3D volume, the time delay from the latest cardiac reference time pointis found and the image is assigned to a corresponding time bin. The number of time binsis set to 30, and the RR-interval is estimated based on the average cardiac rate of thesubject during the acquisition. A cine-video is created by averaging across all imagevolumes within each time bin. In effect, around 10 volumes are assigned to each timebin after discarding those falling outside the defined RR-interval. In each voxel, the magnitude and phase oscillations are calculated as percentage / degrees change from the mean across the cardiac cycle in each voxel. ThePC933674WO - specificationabove described processing is performed using in-house code in Matlab R2023B. EPICis used for EPI distortion correction. The overlay videos are created with FSLeyes.For ROI analysis, single voxel ROIs in the 3D-EPI images are extracted based onlocation guided by the MPRAGE image.In some implementations, quantification may be possible by using (e.g., inverting) theEPG-model in the above embodiment based on the measured phase data. In thisimplementation example, however, a semi-quantitative mapping approach is appliedusing a lookup or mapping table. In this example, the measured phase values in eachvoxel and time point are substituted by its corresponding displacement value based onthe lookup or mapping table, using the sensitivity-relation obtained from forward EPG-simulations. Since the value may depend strongly (may be influenced by) on therelaxation parameters, SPM12 is used to create a WM-GM-CSF segmentation map fromthe MPRAGE, and in each voxel, the sensitivity-value for the correct tissue-type is usedaccording to the segmentation map, where partial volume effect is taken into account.The approach in this example assumes a linear displacement-to-phase response, whichmay be an approximation in practice, particularly in CSF, but may provide a usefulestimate of the real displacement amplitudes. Subsequently, the peak-to-peakdisplacement in each voxel is calculated from the largest distance between any pair ofdisplacement vectors during the cardiac cycle. The results of the simulations are now presented.Fig.17 shows example simulation results for the magnitude and phase of the MR signalas function of time for a 1 Hz sinusoidal displacement function of increasing amplitude^^^, for the GRE and nb-SSFP versions of the 3D-EPI sequence in some embodiments.The simulation in this example only considers CSF and WM. In Fig.17, for GRE (the first case), only very small displacements produce a smoothlyvarying phase, and only in CSF. In this regime, GRE may be extremely sensitive to smalldisplacements and may be able to detect motion below 0.01 mm in CSF. In Fig. 17, fornb-SSFP (FID / S1) (FID signal is readout) and nb-SSFP (ECHO / S2) (echo signal isreadout), smooth magnitude and phase oscillations can be observed for both CSF andWM and for displacement amplitudes up to around 1 mm. For both GRE and nb-SSFP,large displacement amplitudes result in random fluctuations in magnitude and phase.In this example simulation, for CSF, a sinusoidal phase oscillation is observed for very small displacement amplitudes in both GRE and nb-SSFP while apparently random phase values are observed at large displacements. In an intermediate regime, the phase is still periodic but with an irregular shape. The difference between GRE and nb-SSFP is that GRE is much more sensitive than nb-SSFP to CSF displacements and enters theirregular and random phase regimes for smaller displacement amplitudes than nb-SSFP.In this example simulation, for WM relaxation times, no effect on the MR signal is foundfor GRE and, in nb-SSFP, similar phase-curves as for CSF are observed. The doublepeak observed for the magnitude fluctuations in nb-SSFP can be understood by considering that reduced magnitude can result from the summation of individualPC933674WO - specification coherence pathways with different net phase values, which may occur during periods of rapid displacement changes (peak absolute velocity), regardless of the polarity (directionof the velocity). Hence, in this example, the magnitude function has two peaks perdisplacement period.Fig.18 shows the EPG simulated nb-SSFP phase amplitude ^^^ for increasing pulsationamplitude ^^^ at ^^^ = 1 Hz for CSF, WM, and GM in the FID signal of nb-SSFP (TR =10.6 ms, FA = 20 deg, Δ^^ = 3.2 mm). It is found that the displacement sensitivity is abouttwice as high in CSF compared to WM and GM due to the long T1 and T2. The functional relation is strictly monotonic and can be used for or as a look-up table.In Fig. 18, the simulated nb-SSFP phase ^^^ is plotted as a function of displacementamplitude ^^^ up to 1 mm for GM, WM, and CSF. It can be seen that in this range, thenb-SSFP phase is a smooth and strictly monotonically increasing function of displacement amplitude and vice versa. Relaxation times may have a significant impacton the displacement sensitivity, as the phase for CSF is about 2 times higher comparedto the phase for WM and GM tissue, for the same displacement amplitudes.The phase displacement sensitivity and relative phase-to-noise ratio of the FID signaland the echo signal are compared. Figure 19 shows the displacement sensitivity andrelative phase-to-noise ratio for nb-SSFP FID and nb-SSFP echo as a function ofdisplacement amplitude ^^^ at ^^^ = 1 Hz and with TR = 10 ms, FA = 20 degrees, and voxelsize Δ^^ = 3 mm.As shown in the Fig.19, the sensitivity is higher for the nb-SSFP echo than for nb-SSFPFID. However, the relative phase-to-noise ratio is almost identical, as the higher signalmagnitude in FID may compensate the lower phase-displacement sensitivity. As a result,in this example, there is no significant gain from collecting (readout) the echo signalinstead of the FID signal (see, e.g., Fig. 16).Fig.20 shows the comparison of the simulation results from the EPG model in the aboveembodiment and the isochromat averaging simulations. Specifically, Fig. 20 showssimulation results in CSF for sinusoidal displacement at ^^^ = 1^^^^ for the proposed EPGmodel versus Bloch simulations for GRE and nb-SSFP at two selected displacement amplitudes ^^^(0.1mm and 1.0mm). Good agreement is observed in both magnitude and phase of the MR signal. In fact, the agreement is so good that the two curves overlap almost completely, even for irregular and quickly changing phase and magnitudefunctions. In this example, the standard deviation of the difference between the two typesof simulations for the cases shown in the figure are 1.0% for the magnitude and about0.012-0.013 rad for the phase. It should be noted that in this example the EPG modelsimulations are at least 2-3 order of magnitude faster than corresponding Blochsimulations, and therefore it may allow exploration of a much larger set of inputparameters within reasonable computation times. The results of the in vivo measurements are now presented.PC933674WO - specificationFig. 21 shows example nb-SSFP phase time-curves and corresponding frequencyspectra for single voxel ROIs located in the pons (graphs on the left) and the 4thventricle(graphs on the right) from one subject respectively. The shaded area in the frequencyspectra represent the fundamental cardiac and respiration frequency ranges. As shownin Fig. 21, for the 4th ventricle CSF voxel, phase variation is present both in the cardiac,respiratory, and ultra-low frequency range < 0.2 Hz. For the WM voxel in the pons, thephase variations are generally more than a factor of 2 lower, with less clear contributionfrom the respiratory and ultra-low frequency ranges.Figs. 22A to 22I show several single-voxel phase values during the cardiac cycle afterretrograting processing. Specifically, Figs.22A to 22I show in vivo measurement resultsof signal-pulsation for both nb-SSFP and GRE in single voxel ROIs after cardiacretrogating in selected areas and for different directions of the applied spoiler gradient: A to C) Repeatability of nb-SSFP phase-pulsation across two scans in the same subject.The same single voxel ROI location is used in both B and C; D to F) Comparison of nb-SSFP phase-pulsation values from repeated scans with varying spoiler gradient withinthe same subject and voxel ROI; G to I) Comparison of phase and magnitude pulsationin nb-SSFP versus GRE in subject 2.It can be seen that repeatability is good for all three example cases shown, with an average root-mean-square error between two measurements of 0.010 rad for pons and0.019 rad for the lateral ventricle CSF. Based on the graph of Fig. 18, this maycorrespond to a displacement measurement accuracy of around 10μm in both cases.Note that the opposite direction of phase changes in the left versus right lateral ventricles for the image acquisition with spoiler gradient applied along the LR-axis, as shown in Fig.22B. For comparison, the phase change in the same single voxel ROI as in Fig.22B but from the acquisition with spoiler gradient applied along the AP-direction is shown inFig. 22C. Here, the phase changes are in the same direction for both ventricles. Theseobservations match closely the expected symmetry of the CSF motion, and they demonstrate that the sequence is sensitive to motion along the axis of the applied spoiler gradient and with polarity information intact. For tissue motion and CSF in the cerebral aqueduct, where fluid motion can be assumed fairly unidirectional and coherent, the scaling of phase values is linear with the value of the zeroth order spoiler gradient moment M0.In the 4th ventricle, on the other hand, phase amplitudes are generally larger, moreirregular and do not scale with M0. This indicates that the motion is not coherent within the voxel volume and / or that the displacement amplitude is large compared to theexpected smooth range as discussed below (different phase regimes).Figs. 22G to 22I provide comparisons between nb-SSFP FID and GRE for the samelocation single voxel ROIs in Subject 2. It can be seen that for the pons ROI, both thephase and magnitude of the GRE signal show small and non-smooth variations through the cardiac cycle. This is in overall agreement with EPG simulations. For the nb-SSFP magnitude, a tendency towards a double peak behaviour can be observed. For the cerebral aqueduct ROI, the signal fluctuations are generally higher, but show the samePC933674WO - specification overall pattern. For the 4thventricle ROI the signal fluctuations are even higher and more complicated.Fig.23 shows phase-values at three different time points within the cardiac cycle for thethree different spoiler gradient directions. Specifically, a single sagittal slice of nb-SSFPphase is shown for Subject 1 at three different time points during the cardiac cycle andfor each of the three different physical directions of the spoiler gradient.Fig.24 shows displacement vectors at three different time points within the cardiac cycle.Note that the length of the vectors is amplified. As an example, the peak displacementin the brainstem is only around 0.1 mm. More specifically, Fig. 24 shows the estimateddisplacement vector calculated using the function ^^(^^^) obtained from EPG simulations.Together, these results illustrate the rich information content of the brain pulsationimaging method disclosed herein.Fig.25 shows the nb-SSFP phase variation through the cardiac cycle for a representativesingle voxel ROI in the pons for all four subjects. It can be seen that the overall pattern of the phase variation is very similar among the four subjects, who are all young healthy volunteers. Table 3 shows some estimated peak-to-peak displacement values for representativesingle voxel tissue ROIs in Subject 1 together with corresponding values reported in theliterature in this example. It can be seen that the measurements are within the range ofthe literature values.Table 3 – Values for the estimated peak-to-peak displacement in representative singlevoxel ROIs in Subject 1 compared to values found in the literature Voxel Location Displacement [mm] Literature values [mm]Pons 0.19 0.15-0.18Brainstem 0.15 0.19Thalamus Left 0.11 0.05-0.13Thalamus Right 0.09 0.05-0.13Caudate Left 0.07 0.03-0.08Caudate Right 0.09 0.03-0.08Frontal Lobe WM 0.04 0.04-0.09Occipital Lobe WM 0.04 0.03-0.04In some embodiments disclosed herein, for the first time, the intrinsic motion-sensitivityof nb-SSFP has been harnessed as an active and useful contrast in MRI. Compared toexisting MRI methods for detecting motion, the method in some embodiments may bePC933674WO - specification particularly suited for imaging small amplitude pulsatile motion such as the cardiac and respiratory pulsation that may be found in brain tissue and fluid. Based on the simulationsand the in-vivo measurements described, the method in some embodiments may beapplied for quantitative measurements of pulsatile motion. While the phase-to-displacement relation may be non-linear and may involve both varying displacementsensitivity and response delays, it appears to be smooth and strictly monotonicallyincreasing for applicable displacement ranges and tissue types. This may qualify it for alookup or mapping table approach. Further quantification may require solving a complexinverse problem with more variables, including relaxation rates and sequenceparameters. The computational efficiency of the EPG extension model in someembodiments may be useful for future quantification applications.The motion sensitivity of MRI may be known and it may be both a problem which causesimage artefacts and a feature that can be used to detect and measure motion. It shouldbe noted, however, that difference exists between velocity sensitivity versusdisplacement-sensitivity, e.g., depending on whether the 0th or 1st order gradient momentis the dominating contributor to the phase in Equation (4).One classic approach for flow imaging may be to apply a set of bipolar gradients beforeread-out in a GRE sequence. Coherent motion of the spin isochromat during theapplication of these gradients may result in a phase shift proportional to the velocity alongthe direction of the applied gradient. The proportionality coefficient may be determinedby the duration and amplitude of the bipolar gradients. This phase-contrast MRI methodmay be useful for imaging and / or quantification of blood flow, as it may be best suited forrelatively high flow velocities (0.1-2 m / s). With reference to the theory discussed above,classic phase-contrast is only sensitive to velocity and not the position, as the 0th ordergradient moment M0 = 0 while the 1st order gradient moment is high. Similarly, balancedsteady-state free precession (b-SSFP) has been utilized for blood flow imaging. B-SSFPdiffers from nb-SSFP in that there are no spoiler gradients within each TR to separate the different coherence pathways, but rather perfect rewinding so that all pathways are coherently added. In the same way as for the GRE-based phase-contrast MRI, the 0thorder gradient moment is zero by design, which means that the phase is not sensitive tothe position, but only to the velocity through the second term of the Taylor expansion.Some embodiments disclosed herein focus on the motion sensitivity of non-balanced SSFP. Due to the spoiler gradient, the 0thorder moment within each TR-period is notzero, and the phase at the end of the TR-period can be sensitive to the position of thespin isochromat during that TR. The total phase of the FID and echo signals may bedetermined by the complex addition of pathways with different history and thereforedifferent phase values, if the position of the spin isochromat is changing over time. Theinventors have found that steady motion in terms of constant flow velocity will not resultin a measurable phase effect in nb-SSFP. This is due to the inversion of phase valuesthrough the refocusing of pathways, so that ^^^and ^^^states may have opposite phasevalues which may effectively cancel each other out. nb-SSFP may therefore not becategorized as a velocity-encoding method, but rather a displacement-encoding method.It is noted that in some cases even coherent displacement may affect the magnitude ofboth the GRE and the nb-SSFP signal, as long as it is pulsatile.PC933674WO - specificationAs illustrated in Fig. 17, in some embodiments, the phase-effect on the nb-SSFP signalcan be separated in three regimes: (1) a smooth regime for very small displacements,where the phase may closely follow the displacement function; (2) a transition regime,where the phase-effect may still resemble the overall shape of the displacement function,but with some deviation; and (3) a quasi-random regime, where there are seeminglyrandom jumps in both phase and magnitude from one TR to the next. Taken together with the lack of sensitivity to constant velocity flow, nb-SSFP may become rather insensitive to arterial blood flow. Any systematic phase fluctuations observed in the datamay likely originate from low amplitude tissue or fluid pulsatile displacement.The displacement amplitude at which the transition from one phase-regime to the next occurs may depend on a number of parameters, including the 0thorder gradient momentM0, TR, FA, and relaxation parameters, for example. Some embodiments disclosedherein focus on the combination of short TR and low M0 (e.g., corresponding to 2π dephasing at 3mm voxel size) to enable the smooth range to extend up to around^^^ = 1 mm. It may be possible to reach similar “smooth range” for higher M0 (and therebyat higher spatial resolution) if longer TR is used. As the volume-TR should remainreasonably short and each volume should not require too many excitations (TR-intervals), this may require slab-selective thin-slab coverage or higher undersampling.Various examples of simulations versus in-vivo measurements are provided herein. Ingeneral, the in-vivo measurements may reflect the findings from EPG simulations. Both GRE and nb-SSFP data may demonstrate phase and magnitude fluctuations associated with cardiac and respiratory induced tissue and fluid motion. For nb-SSFP, significant and reproducible temporal phase fluctuations have been observed during the cardiac cycle that reflect expected displacement patterns in both tissue and CSF. For GRE,random signal fluctuations of relatively low amplitude in tissue and some moresystematic signal fluctuations in CSF have been observed. Overall, these observationsof the in-vivo experiments are in agreement with the simulation results.However, there are also some differences. For example, no clear double-peak behaviouris observed for the nb-SSFP magnitude in the in-vivo data (as predicted from simulations), although there is a tendency towards such a double-peak in the pons andcerebral aqueduct ROIs as shown in Figs. 22A to 22I. This discrepancy may be due toother effects that could affect the signal magnitude, especially in CSF, e.g., non-coherentmotion such as turbulent flow. Non-rigid motion may also lead to an underestimation ofthe peak displacement values within each voxel. The comparison of phase-values for different amplitudes of the spoiler gradient revealsthat in the pons where rigid coherent motion of low amplitude is expected, the phasevalues scale as expected. A similar observation exists for the cerebral aqueduct. Thismay indicate that the motion in the aqueduct is dominated by laminar flow. In the 4thventricle on the other hand, no systemic scaling of the phase with the value of M0 isobserved. A possible explanation for this is that the displacements may be large andbeyond the smooth phase regime.PC933674WO - specification Repeated measurements in the same subject during the same scan-session as shownin Figs.22A to 22I allow an estimation of the measurement precision. It has been foundthat the average root-mean-squared-difference between the curves in Figs.22A to 22Cmay represent a displacement of around 10 µm after applying the semi-quantitative map-ping approach based on the graph in Fig.18. Such high measurement precision may allow the detection of longitudinal variation and changes in the brain pulsation patterns,e.g., as result of natural variation or pathological processes.The accuracy of the displacement measurements is more difficult to evaluate. Thecomparison between the estimated peak-to-peak displacement values to other valuesreported in the literature in Table 3 suggests that the quantification approach applied inthe above embodiments is generally sound, and that the EPG model may provide acorrect theoretical framework for the method and image contrast.Some embodiments may be only based on the signal phase and therefore may onlydetect displacements that are coherent on the length scale of one voxel. However, sincethe magnitude of the nb-SSFP signal may carry information about both coherent andincoherent motion, an extension of the approach in some other embodiments could beenvisioned where combined analysis of the phase and magnitude may allowquantification of both coherent and incoherent motion.Another component of measurement accuracy relates to the timing of the cardiacretrogating relative to the arterial pulsation in the brain. The retrograting applied may bebased on a pulse oximeter on the index finger, where T=0 is defined as the peak slope of the ascending ridge of the pulse-curve. This time point should in theory be after the arrival of the arterial pulse wave to the brain, since it is known that the arterial pulse wave arrives later in the index finger than in the brain. This decay can vary between individuals. However, as illustrated such variation is not significant among the four subjects in the experiments. In some embodiments, the additional phase-term in the EPG model included to describe the effect of time-dependent position and velocity may have other applications. EPG is in general an order of magnitude faster than Bloch-simulations, and may allow simulations of a much wider set of parameters than Bloch simulations. In addition, the EPG model may offer more intuitive understanding of the dynamics behind the effectscaused by motion. The formalism may be general, and may find further applications inother short-TR / low FA type sequences, for example MR fingerprinting.Based on the importance of brain pulsation in CSF circulation and brain health, themethod disclosed in some embodiments may have several potential clinical applications.In the embodiments, it is somewhat unexpected that the spoiled GRE sequence can beused to detect coherent displacements with extremely small amplitudes. Based on the simulations, pulsatile motion with amplitudes well below 10 µm could potentially bedetected and quantified using such sequence. This may only hold in the case of long T1and T2, such as in CSF or other fluid compartments. This effect, which is believed to bePC933674WO - specificationnovel, may have novel applications in various fields, e.g., because it can be used tomeasure displacements within a full 3D volume. Embodiments disclosed herein have provided a MRI based method for imaging and / ordetecting pulsatile displacement, e.g., based on non-balanced steady state freeprecession or spoiled GRE, combined with a 2D- or 3D- readout, such as a highlyaccelerated 3D-EPI readout. The method may be used to measure pulsatile displacement of both fluids and tissues, e.g., with sub second time resolution. Gradient spoilers may be used for controlled sensitization to displacements along a chosen direction. In some embodiments, by repeating the measurement at least three times with orthogonal spoiler gradients and optionally applying cardiac retro-gating, full 3D motion vectors may be estimated in each voxel throughout the cardiac cycle.It will be appreciated by one skilled in the art that variations and / or modifications may bemade to the described and / or illustrated embodiments to provide other embodiments. The described and / or illustrated embodiments should therefore be considered in all respects as illustrative, not restrictive. For example, in some embodiments, the MRI device can be any imaging device operable to perform, at least, magnetic resonance imaging. For example, in some embodiments, the MR data processing may be performed at least partly online (e.g., in real time as an imaging subject is being imaged by the MRI scanner). For example, in some embodiments, the MR data processing may be performed offline (i.e., after the imagingsubject has been imaged by the MRI scanner). For example, in some embodiments, thespecific parameters of the MRI imaging protocol / sequence (e.g., the pulse timing, thenumber of pulses, the number of gradient, the timing of gradient, etc.) may be modifiedor optimized depending on applications and may be different from those describedherein.PC933674WO - specification

Claims

CLAIMS1. A method of imaging and / or detecting pulsatile displacement using magneticresonance imaging (MRI), the method comprising: performing a scan on an imaging subject using an MRI device by causing the MRI device to apply an MRI pulse sequence with a 2D- or 3D- based readout to an imagingsubject, the MRI pulse sequence comprises a plurality of repetitions (TR), each of the repetitions having displacement encoding gradient for sensitization to pulsatiledisplacement and the 2D- or 3D- based readout applied either before or after thedisplacement encoding gradient, the displacement encoding gradient comprises agradient lobe with a finite zero-order moment, andacquire MR data of the imaging subject; andprocessing the MR data to obtain phase and / or magnitude data associated with the pulsatile displacement.

2. The method of claim 1, wherein the pulsatile displacement is pulsatile displacementof at least one of: brain tissue, such as white matter, or cerebrospinal fluid (CSF).

3. The method of claim 1 or 2, wherein the zero-order moment of the gradient lobe hasa magnitude of at least 2π / γ∆x, where ∆x is voxel size of the MR data and γ is gyromagnetic ratio.

4. The method of any one of claims 1 to 3, wherein the 2D- or 3D- based readout is a3D echo planar imaging (EPI) based readout and the MR data comprises volumetricMR data.

5. The method of claim 4, wherein the 3D EPI based readout is a segmented 3D EPIreadout.

6. The method of any one of claims 1 to 5, wherein for each of the repetitions, the MRIpulse sequence lacks velocity encoding gradients before the 2D- or 3D- basedreadout.

7. The method of any one of claims 1 to 6, wherein the displacement encoding gradientsof the plurality of repetitions comprise orthogonal gradients for sensitization to pulsatile displacements along orthogonal directions.

8. The method of any one of claims 1 to 7, wherein for each of the repetitions, the MRIpulse sequence applies RF spoiling.

9. The method of any one of claims 1 to 7, wherein for each of the repetitions, the MRIpulse sequence does not apply RF spoiling.

10. The method of any one of claims 1 to 7, wherein the MRI pulse sequence is a non-balanced SSFP sequence.PC933674WO - specification11. The method of any one of claims 1 to 7, wherein the MRI pulse sequence is agradient-echo (GRE) based sequence.

12. The method of any one of claims 1 to 11, wherein an amplitude of the pulsatiledisplacement has an order of magnitude of 10-2mm.

13. The method of any one of claims 1 to 11, wherein an amplitude of the pulsatiledisplacement has an order of magnitude of 10-1mm.

14. The method of any one of claims 1 to 11, wherein an amplitude of the pulsatiledisplacement has an order of magnitude of 1 mm.

15. The method of any one of claims 1 to 14, wherein a temporal resolution of thepulsatile displacement is in an order of 1 s to 10-1 s.

16. The method of any one of claims 1 to 7, wherein the imaging subject comprises abrain of a subject, and the acquisition of the MR data is based at least in part on prospective electrocardiogram (ECG) or respiratory gating associated with thesubject.

17. The method of any one of claims 1 to 7, wherein the imaging subject comprises abrain of a subject, and the processing of the MR data is based at least in part on retrospective ECG or respiratory gating associated with the subject.

18. The method of any one of claims 1 to 17, further comprising:converting phase data to displacement data associated with the pulsatile displacement based at least in part on a lookup or mapping table.

19. The method of any one of claims 1 to 17, further comprising:converting phase data to displacement data associated with the pulsatile displacement based at least in part on a model such as a Bloch Simulation model or an extended phase graph (EPG) model.

20. A system for imaging and / or detecting pulsatile displacement using MRI, comprising:an MRI device configured to: apply an MRI pulse sequence with a2D- or 3D- based readout to an imagingsubject, the MRI pulse sequence comprises a plurality of repetitions (TR), each of the repetitions having displacement encoding gradient for sensitization to pulsatiledisplacement and the 2D- or 3D- based readout applied either before or after thedisplacement encoding gradient, the displacement encoding gradient comprises a gradient lobe with a finite zero-order moment; and acquire MR data of the imaging subject, and at least one processor configured to: process the MR data to obtain phase and / or magnitude data associated with the pulsatile displacement.PC933674WO - specification

Images