Method for determining a personalized cardiac model using a magnetic resonance imaging sequence

Abstract

A method is provided for determining a personalized cardiac model, including steps of (i) computing a velocity time profile of a blood flow across a selected area of the heart or the aorta during at least one cardiac cycle, using data acquired with a Magnetic Resonance Imaging (MRI) device; (ii) performing a segmentation of the velocity time profile so as to identify cardiac phases according to a predefined generic cardiac model; and (iii) computing normalized time location and/or duration of the cardiac phases within cardiac cycles so as to define a personalized cardiac model.

Claims

1. A method for determining a personalized cardiac model, comprising steps of: computing a velocity time profile of a blood flow across a selected area of the heart or the aorta during at least one cardiac cycle, using data acquired with a Magnetic Resonance Imaging (MRI) device; performing a segmentation of said velocity time profile so as to identify cardiac phases according to a predefined generic cardiac model; computing at least one of a time location or a duration of said cardiac phases within cardiac cycles so as to define a personalized cardiac model; acquiring and processing a Real-Time Phase Contrast (RTPC) signal sequence in a scan plane with a spatial orientation normal to a direction of flow of the blood, said RTPC sequence being based on a phase contrast MRI sequence where only the central k-space line is acquired; computing a 1D+t velocity map with one spatial axis and one time axis, by applying to the RTPC signal sequence unidimensional Fourier transforms along the frequency encoding direction; carrying out a Singular Value Decomposition (SVD) on the 1D+t velocity map; and with a threshold algorithm applied on the spatial weighting returned by the SVD, identifying in the 1D+t velocity map a set of regions of interest (ROI) of contiguous pixels for one or several of the first principal components issued from the SVD decomposition.

2. The method of claim 1, which further comprises a recording of an ElectroCardioGram (ECG) signal, said ECG signal being used for triggering the Magnetic Resonance Imaging (MRI) acquisitions.

3. The method of claim 1, which further comprises a step of identifying a region of the 1D+t velocity map in which the velocity time profile is representative of the cardiac phases.

4. The method of claim 1, which further comprises steps of: computing a spatial average value of the velocity time profile in the ROIs; identifying the ROI with the spatial average value having the highest power in the cardiac spectral band around the heart beat rate; and computing the velocity time profile by using the spatial average value of the ROI having the largest size along the spatial dimension axis among all ROIs having more than 40% of said highest power.

5. The method of claim 1, wherein the RTPC sequence uses a frequency encoding direction with an orientation chosen so as to project velocity fields from different sources to areas as distinct as possible along the frequency encoding axis.

6. The method of claim 5, wherein the acquisition and processing of the Real-Time Phase Contrast (RTPC) signal sequence is done in one of the following configurations: in a quasi axial scan plane cross-sectioning the ascending aorta with a frequency encoding direction going through the ascending and descending aorta cross-sections; in an intra-cardiac short-axis scan plane with a frequency encoding direction going through inferolateral and anteroseptal cardiac segments.

7. The method of claim 1, which further comprises steps of: computing a baseline of the velocity time profile representative of a null blood flow velocity; and for at least one cardiac cycle, (i) detecting a peak with the highest velocity corresponding to the systolic cardiac phase by fitting with a pre-defined curve shape, and (ii) identifying at least one of its time location or its duration within said cardiac cycle by locating the times of crossing of the baseline by the velocity time profile curve around said peak.

8. The method of claim 7, wherein the baseline is computed iteratively using the following steps: computing a baseline by applying a low-pass filter with a cut-off frequency of 0.35 Hz to the velocity time profile; and using histogram analysis and application of a threshold, generating a new velocity time profile by setting extreme or most distant values of the velocity time profile relative to the baseline to the baseline value.

9. The method of claim 7, which further comprises steps of, for at least one cardiac cycle: detecting the peak with the highest velocity among the still non-identified peaks by fitting with a pre-defined curve shape, and matching it with a cardiac phase using the generic cardiac model; and identifying at least one of its time location or its duration within said cardiac cycle by locating the times of crossing of the baseline by the velocity time profile curve around said peak.

10. The method of claim 7, which further comprises a step of computing a curve fitting at least one of the time location or the duration of the identified cardiac phases to the cardiac cycle duration, so as to generate the personalized cardiac model.

11. A medical imaging method, comprising steps of: determining a personalized cardiac model by using the method of claim 1; and acquiring medical imaging data of the heart using said personalized cardiac model for taking into account the cardiac phases.

12. The MRI medical imaging method according to claim 11, which comprises a step of triggering a MRI sequence using the personalized cardiac model.

13. The MRI medical imaging method according to claim 11, which comprises steps of: asynchronous acquisition of MRI data; and reconstruction of a temporal resolution Cine MRI sequence using the personalized cardiac model.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The methods according to embodiments of the present invention may be better understood with reference to the drawings, which are given for illustrative purposes only and are not meant to be limiting. Other aspects, goals and advantages of the invention shall be apparent from the descriptions given hereunder.

(2) FIG. 1 shows a flow chart of the method of the invention,

(3) FIG. 2 shows a velocity map corresponding to a scan plane located on the aorta,

(4) FIG. 3 shows (a) an intra-cardiac short-axis scan plane, and (b) a corresponding velocity map,

(5) FIG. 4 shows the spatial weights of the three first modes of the SVD of the velocity map of FIG. 2, in the aorta,

(6) FIG. 5 shows a small portion of velocity time profiles extracted from several possible ROIs in the velocity map of FIG. 2,

(7) FIG. 6 shows a velocity curve obtained in intra-cardiac short-axis scan plane,

(8) FIG. 7 shows portion of the velocity time profile obtained in aorta scan plane with overlaid R-wave detected on ECG, start of ejection and end of ejection,

(9) FIG. 8 shows portion of the velocity time profile obtained in intra-cardiac short-axis scan plane with overlaid R-wave detected on ECG and detection of: start of ejection, end of ejection, start of early filling, end of early filling, start of atrial filling and end of atrial filling,

(10) FIG. 9 shows measurements of end-systolic time for one person, and corresponding end-systolic time predicted by the model of the invention and models of the prior art.

DETAILED DESCRIPTION

(11) With reference to FIG. 1, we will now describe a mode of implementation of the invention which allows the construction of a personalized cardiac model adapted to each subject undergoing Cardiac Magnetic Resonance (CMR).

(12) The method of the invention comprises a first step 10 of selection of a slice or a scan plane to be used to calibrate the model.

(13) There for, a Cardiac MRI is performed. It can be done for instance using a 3 T Signa HDxt scanner from General Electric, Waukesha, Wis., with subjects in supine position and using an eight-element cardiac phased-array coil.

(14) Localizing scans are recorded, and a slice or a scan plane usable for the calibration is located.

(15) The choice of the scan plane depends on the kind of cardiac model which is to be used. For instance: a scan plane located at the aorta may be used to obtain a cardiac model with N=2 cardiac phases (systole/diastole); an intra-cardiac short axis scan plane may be used to obtain a cardiac model with up to N=7 cardiac phases (isovolumic contraction, ejection, isovolumic relaxation, early atrial filling, diastasis, atrial filling and atrioventricular delay).

(16) The method of the invention further comprises a step 11 of acquisition of a Real-Time Phase Contrast (RTPC) signal in the chosen slice, normal to the direction of velocity of the blood flow.

(17) The RTPC sequence is based on a bidimensional (2D) phase contrast MRI sequence, where only the central k-space line is acquired, similarly to the RACE sequence. One-directional through-slice interleaved velocity encoding is used. The velocity is measured successively in two opposite directions for each line of the k-space. The final value of the velocity is obtained by computing a difference between successive measurements in opposite directions.

(18) The acquisition may be done for instance with the following MRI scan parameters: Field Of View (FOV): 350 mm; Slice thickness: 8 mm; Flip angle: 15 degrees; Bandwidth: 62.5 KHz; Acquisition matrix (MP): 2561; Repetition time (TR): 6.6 ms Echo time (TE): 3.4 ms; Encoding velocity (Venc): 50 cm/s; Excitations repetitions: 1; k-Lines per segment: 1; Temporal resolution: 6.6 ms; Typical scan duration: 100 s.

(19) With these parameters, the temporal resolution for RTPC acquisitions is 6.6 ms.

(20) The orientation of the frequency encoding direction is chosen so as to project different velocity fields from different sources to separate areas on the frequency line, along the frequency encoding direction (at least as much as possible, some overlap may sometimes be unavoidable).

(21) The scan planes may comprise for instance: a quasi axial scan plane cross-sectioning the ascending aorta with a frequency direction going through the two aorta cross-sections (ascending and descending aorta); a scan plane in a so-called short axis with a frequency direction going through inferolateral and anteroseptal cardiac segments. The short axis corresponds to specific heart cross section orientations, in which both ventricles are visible.

(22) The RTPC acquisition is performed during several heartbeats (for instance 128 heartbeats) and during varying heart rates (for instance using natural heart rate variation during free-breathing).

(23) ElectroCardioGram (ECG) is recorded and used for triggering the MRI.

(24) Cardiac cycles for which an obvious error in QRS detection on the ECG occurs are discarded from the analysis (i.e. for instance cycles with instantaneous heart rate below 30 or above 120 beats per minute).

(25) The method of the invention further comprises a step 12 of computing of a velocity map.

(26) The RTPC data is transformed into a 1D+t image space by applying a unidimensional (1D) Fourier transform along the frequency encoding direction.

(27) The 1D+t image space is an image space with one coordinate being a dimension (1D, for instance in pixels), and one coordinate being the time (t).

(28) To limit noisy air regions and reduce processing time, only the central part corresponding to about 50% of the field of view (FOV) is kept, while checking that the organ of interest is included.

(29) For each of the MRI receiving coils, the 1D+t phase dataset is unwrapped using a fast bidimensional (2D) phase unwrapping algorithm. A velocity map is then computed by subtraction of data relative to consecutive echoes of the interleaved acquisitions. A temporal resolution equal to TR is obtained by computing the forward and the backward differences.

(30) Finally, a global 1D+t velocity map is computed, which corresponds to a magnitude-weighted sum of the respective coil's velocity maps.

(31) FIG. 2 shows an example of velocity map corresponding to a scan plane located on the aorta, which is suitable for determining a cardiac model with N=2 cardiac phases (systole/diastole).

(32) One axis corresponds to the time of acquisition. Three cardiac cycles (out of 128) are represented.

(33) The other axis corresponds to the spatial direction of the slice perpendicular to the direction of projection of the velocity field across the slice. So, one velocity profile is represented for each acquisition time.

(34) Positive velocity (in bright areas 20) and negative velocity (in dark areas 21) correspond respectively to systolic blood flow in the ascending and descending aorta.

(35) Similarly, FIG. 3b shows an example of velocity map corresponding to an intra-cardiac short-axis scan plane, which is suitable for determining a cardiac model with N=6 cardiac phases. One second of RTPC signal containing one cardiac cycle (out of 128) is represented.

(36) The intra-cardiac short-axis scan plane 30 from which the velocity map is computed is shown on FIG. 3a.

(37) The method of the invention further comprises a step 13 of computation of a velocity time profile.

(38) An automatic spatial segmentation of the 1D+t velocity map is performed in order to determine the area in which the velocity time profile is the most representative of the cardiac cycles.

(39) A Singular Value Decomposition (SVD) is carried out on the 1D+t velocity map.

(40) FIG. 4 shows the spatial weights of the three first modes of the SVD of the velocity map of FIG. 2, corresponding to the aorta. More precisely, it shows the spatial weight of the first mode 40, the spatial weight of the second mode 41 and the spatial weight of the third mode 42.

(41) The spatial weights of the modes are represented in function of the same spatial dimension (or direction) as the velocity map.

(42) A threshold algorithm using the spatial weighting returned by the SVD is applied to the velocity map in order to generate a set of spatial locations, corresponding to a set of regions of interest (ROI) of contiguous pixels, for each of the first three principal components issued from the SVD decomposition.

(43) For each ROI, the time course (of time profile) of the mean velocity in the ROI is computed. The ROI with the velocity time course having the highest power in the cardiac spectral band (using the mean heart rate+/0.03 Hz) is determined, and used to select all ROIs having more than 40% of this maximum. Among those selected ROIs, the ROI 31 having the largest size in pixels (or spatial unit) is chosen.

(44) That ROI 31 is shown on FIG. 3b for the intra-cardiac short-axis scan plane case.

(45) Finally, the velocity is computed for each acquisition time by cumulating or averaging the corresponding velocity values of the pixels of the selected ROI.

(46) As a result, a velocity time profile is obtained. In the current example it is 128 cardiac cycles long.

(47) FIG. 5 shows an example of velocity time profiles 50 extracted from several candidate ROIs of a velocity map corresponding to the aorta (FIG. 2). The selected velocity time profile 51 is the one automatically chosen by the algorithm, corresponding to the ROI having the largest size.

(48) In the aorta scan plane (corresponding to FIG. 2 and FIG. 5), the velocity time profile can be interpreted as the consequence of the ejected ascending aortic blood flow because, physiologically, ascending aorta is the location of highest variations of velocity in the chosen slice. The absolute values of velocity are not considered because they derive from the summation of all velocities along the phase encoding direction (as only the central k-space line is obtained). However, the duration of the highest peak of velocity (called S-wave) is considered as a measure of systole duration. So, only two phases (N=2) may be detected, namely the systole and the diastole.

(49) In the intra-cardiac short-axis scan plane location (corresponding to FIG. 3), the velocity time profile can be interpreted as the velocity of blood inside the heart during its phases, and so N=6 cardiac phases may be detected.

(50) FIG. 6 shows an example of velocity time profile obtained in intra-cardiac short-axis scan plane. Cardiac phases are shown, with the cardiac cycle (RR), the IsoVolumic Contraction (IVC), the Systolic Ejection (SE), the IsoVolumic Relaxation (IVR) and the Diastolic Diastasis (DD).

(51) The method of the invention further comprises a step 14 of detection of the cardiac phases, by processing the selected velocity curve 51.

(52) With reference to FIG. 7, systolic waves (S-waves) are automatically detected on the velocity curve 51 by detecting the peaks with the maximum velocity, with a shape fitted to a parabola (second order polynomial).

(53) The start time 71 and the end time 72 of forward systolic velocity, defined by the S-wave crossing a baseline 70, are recorded.

(54) The baseline 70 is computed iteratively as follows: a current baseline is computed by applying a low-pass filter with a cut-off frequency of 0.35 Hz, close to the heart beat rate, to the current velocity curve; then, using histogram analysis and application of a threshold, a new velocity curve is generated by setting extreme or most distant values of the current velocity curve relative to the current baseline to the baseline value.

(55) The new velocity curve is then used as current velocity curve in the next iteration; at each further iteration, the threshold for extreme values is lowered.

(56) The duration of systole is computed, as illustrated in FIG. 7, by measuring the delay between the ECG triggering (R-wave 73) and the end of forward systolic velocity 72. The online R-wave detection is post-processed in order to reposition the detection at the top of the QRS complex of ECG lead I.

(57) When systole duration is not plausible (for instance with values below 100 ms or above 500 ms), the corresponding cardiac cycle is discarded.

(58) The invention is not limited to N=2 cardiac phases.

(59) With reference to FIG. 8, when velocity containing distinct waves related to cardiac motion can be measured (using for instance intra-cardiac short-axis scan plane as shown in FIG. 3 and FIG. 6), the invention can also calibrate these waves in function of heart rate.

(60) Other waves (than systolic waves) are automatically detected on the velocity curve by iteratively detecting the peaks with the next maximum velocity, with a shape fitted to a parabola (second order polynomial).

(61) The start time and the end time of each peak, defined by the corresponding wave crossing the baseline 70, are recorded.

(62) The duration of the peaks and/or the gaps between the peaks (corresponding to the cardiac phases) is then computed by measuring the delay between the start time and the end time of these peaks or gaps.

(63) FIG. 8 illustrates shows portion of the velocity time profile 51 obtained in intra-cardiac short-axis scan plane, with overlaid R-wave 73 detected on ECG and detection of: start of ejection 71, end of ejection 72, start of early filling 82, end of early filling 83, start of atrial filling 80 and end of atrial filling 81 defined by baseline 70 crossings of second order polynomial fitting the peaks.

(64) In the intra-cardiac short-axis scan plane, six cardiac phases can be thus be computed (three waves: S, E, A and their gaps).

(65) The method of the invention further comprises a step 15 of definition of a personalized, Patient-Adaptive Cardiac Model.

(66) The cardiac phase durations are computed for each valid cardiac cycle as previously explained. Their respective duration (and/or their time location within the heart beat cycle) is then matched or fitted to the heart beat rate. Preferably, a first order (linear) law is used for the curve fitting.

(67) So, we obtain a personalized Patient-Adaptive Cardiac Model which consists in a linear model expressing wave peaks duration as a function of heart rate. As the model is derived from measurements made on the patient, it fits for each subject separately.

(68) FIG. 9 shows examples of measurements of end-systolic time for one volunteer in function of the heart rate. The measurements are made using an aortic scan plane.

(69) The plot shows the individual measurements 90, the values 91 predicted by the Patient-Adaptive Cardiac Model of the invention, and values 92 predicted by three models of the prior art (Weissler, Weissler global and Weissler shifted).

(70) As it can be seen, the Patient-Adaptive Cardiac Model of the invention predicts the end systole time more reliably than all other tested models.

(71) It should be noted that in the model of the invention both the calibration and the error in prediction are dependent upon the temporal resolution of the calibration acquisition. Of course the higher the temporal resolution, the better (in the described mode of realization we achieved 6.6 ms temporal resolution).

(72) We will now describe possible use or implementation of the adaptative cardiac model of the invention. Of course, these examples are in no way limitative.

(73) According to a mode of implementation, the adaptative cardiac model of the invention may be used to trigger a MRI sequence.

(74) In a first step, the model is calibrated for the patient as described previously.

(75) Then, the calibrated cardiac model is used as a drop-in replacement for the global cardiac model used in the adaptive black blood triggered sequence.

(76) So, advantageously, the trigger delay of the acquisition window is adapted to the specificity of the patient. Also, the duration of the acquisition window can be adapted because the length of the cardiac phase is known.

(77) As a result, the acquisition stays more consistently in the prescribed cardiac phase (systole or diastasis) and the image quality is improved (less mixing of cardiac phases).

(78) According to another mode of implementation, the adaptative cardiac model of the invention may be used for the reconstruction of a high temporal resolution Cine MRI sequence.

(79) In a first step, the model is calibrated for the patient as described previously.

(80) Then, the calibrated cardiac model is used as a drop-in replacement for the global cardiac model used in the reconstruction of free-breathing asynchronous balanced-ssfp sequence in the Cine-GRICS algorithm.

(81) As a result, the distribution of asynchronously acquired k-space lines into cardiac phase bins (images) is more correct regarding the real cardiac phases of the patient. The image quality is improved (less mixing of cardiac phases).

(82) Generally speaking, when using MRI imaging modalities, the cardiac model of the invention may be calibrated in a preliminary step with the MRI device, and then used in the measurement sequence. In recurrent exams however, the model may be calibrated once for a patient and then used in several measurement sequences.

(83) A model calibrated for a patient using a MRI device may of course be used with any other relevant imaging modality, such as computed tomography (CT) or untrasounds.

(84) While this invention has been described in conjunction with a number of embodiments, it is evident that many alternatives, modifications and variations would be or are apparent to those of ordinary skill in the applicable arts. Accordingly, it is intended to embrace all such alternatives, modifications, equivalents and variations that are within the spirit and scope of this invention.