A pharmacokinetic analysis method based on breast DCE-MRI

Abstract

The present application proposes a deconvolution method to accurately estimate the corresponding pharmacokinetic parameters from breast DCE-MRI, analyzes the pharmacokinetic parameters of the heterogeneous region, and explores their clinical value in predicting pathological information. The steps include: collecting data of breast cancer patients; tumor segmentation is performed on the image to obtain the tumor region; the concentration time series with noise is synthesized, the deconvolution method is used to estimate the pharmacokinetic parameters of the synthesized concentration time series, and the root mean square error is used to verify the effectiveness of the deconvolution method; the deconvolution method is used to estimate the pharmacokinetic parameters of breast DCE-MRI, and the pharmacokinetic parameters of the heterogeneous region are further analyzed; the image features of the pharmacokinetic parameter map are extracted, and a classification model is constructed to predict the pathological information of breast cancer, such as molecular typing, Ki-67 and pathological grading.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of medical image processing, and particularly relates to a PK analysis method based on breast DCE-MRI. BACKGROUND

[0002] DCE-MRI can be used to non-invasively evaluate the structural and functional characteristics of breast tumors, and it plays a very important role in the precise diagnosis and treatment of breast tumors. The quantitative analysis method based on DCE-MRI can be divided into two categories, namely non-parametric analysis and parametric analysis. In particular, parametric analysis needs to construct a tracer kinetic model to estimate the PK parameters of breast DCE-MRI, including volume transfer constant K trans and rate constant K ep . These parameters can be used to evaluate the hemodynamic characteristics, which reflect the growth rate and degree of tumor growth. However, there is often noise in the collected images, which leads to a deviation between the PK parameters estimated by the traditional tracer kinetic model and the true value. In addition, the heterogeneity of tumors hinders the precise diagnosis and treatment of breast tumors, and whether the PK analysis of heterogeneous regions can improve the predictive performance of the model on pathological information is still unknown.

[0003] In order to facilitate the reader's understanding, the English abbreviations appearing in the present application are explained as follows:

[0004] Dynamic contrast-enhanced magnetic resonance imaging: English full name is dynamic contrast-enhanced magnetic resonance imaging, and the present application is referred to as DCE-MRI.

[0005] Pharmacokinetics: English full name is pharmacokinetic, and the present application is referred to as PK.

[0006] Deconvolution: English full name is deconvolution, and the present application is referred to as DEC.

[0007] Signal-to-noise ratio: English full name is signal-to-noise ratio, and the present application is referred to as SNR.

[0008] Root mean square error: English full name is root mean square error, and the present application is referred to as RMSE.

[0009] Time to peak: English full name is time to peak, and the present application is referred to as TTP.

[0010] Receiver operating characteristic: English full name is receiver operating characteristic, and the present application is referred to as ROC.

[0011] Area under the ROC curve: the full name in English is area under the ROC curve, and the invention is abbreviated as AUC. SUMMARY

[0012] In order to alleviate the influence of noise in the image on the estimation of PK parameters, the present application proposes a deconvolution method to estimate the PK parameters of breast DCE-MRI on the basis of the Tofts model. The method converts the convolution operation in the Tofts model into matrix vector multiplication to increase the noise resistance, so that the PK parameters of breast DCE-MRI can be more accurately estimated. Further, the PK parameters of the heterogeneous region are analyzed to more accurately predict the pathological information of breast cancer, including molecular typing, Ki-67 and pathological grading.

[0013] The present application provides a pharmacokinetic analysis method based on breast DCE-MRI, comprising the following steps:

[0014] Step 1: Collect the data of breast cancer patients, including image data and pathological information, and select breast cancer patients with complete data from them;

[0015] Step 2: For image data, fuzzy C-means clustering and Markov random field algorithms are used for tumor segmentation to obtain tumor regions;

[0016] Step 3: Synthesize concentration time series under different SNRs, use deconvolution method to estimate the PK parameters of the synthesized concentration time series, and verify the effectiveness of the deconvolution method through RMSE index;

[0017] Step 4: Use the deconvolution method to estimate the PK parameters of breast DCE-MRI to obtain the corresponding PK parameter map, and further analyze the PK parameters of the heterogeneous region;

[0018] Step 5: Extract image features from the PK parameter maps of the whole tumor region and the tumor sub-region, construct a diagnostic model to predict molecular typing, Ki-67 and pathological grading.

[0019] As preferred, in step 3, SNR is selected as 20, 25, 30 and 35, and 50 concentration time series with Gaussian noise are synthesized for each value;

[0020] The values of K trans and K ep are randomly selected from the physiological range of 0.05-0.55 and 0.75-2.75 respectively, and K trans <K ep ; the concentration time series with Gaussian noise is calculated by formulas (1)-(3);

[0021] The plasma concentration C at time t is calculated p (t) :

[0022] C p (t) = D - (a1- exp(-m1-t) + a2- exp(-m2-t)), t = 0, 1,..., 5 (1)

[0023] where D = 0.2 mmol / kg represents the dose of tracer, a1= 3.99 kg / L and a2= 4.78 kg / L represent two amplitude constants, m1= 0.1440 min -1 and m2= 0.0111 min -1 represent two rate constants;

[0024] The tracer concentration C in the tissue at time t is calculated tissue (t) :

[0025]

[0026] where K trans represents the volume transport constant, K ep represents the rate constant;

[0027] The real tracer concentration at time t is calculated adding zero-mean Gaussian noise;

[0028] C tissue_s (t) = C tissue (t) + ε(t), t = 0, 1,..., 5 (3)

[0029] The covariance matrix of ε is∑(0.9 |i-j| )σ 2 , σ 2 and SNR satisfy equation (4)

[0030]

[0031] 50 pixels with concentration time series with Gaussian noise are obtained.

[0032] As a preferred, the deconvolution method is used in step 3 to estimate the corresponding PK parameters from the synthetic concentration time series;

[0033] The deconvolution operation is performed on equation (2) to re-express equation (2) in the following vector matrix form:

[0034] C tissue = K trans · H(K ep ) C p

[0035]

[0036] Here, H(K ep ) is a Toeplitz matrix, and L is the length of the time series;

[0037] Solving the optimization problem of formula (6) can obtain the PK parameters of the concentration time series:

[0038]

[0039] H(K ep )C p ||2 (6)

[0040] s.t.0<K trans <K ep <4

[0041] K ep <4 is set as a reasonable physiological value;

[0042] The RMSE index is used to verify the accuracy of the PK parameters estimated by the deconvolution method relative to the true PK parameters:

[0043]

[0044] K is the PK parameter estimated by the deconvolution method, K * is the true PK parameter, and N is the total number of synthesized pixels.

[0045] As preferred, in step 4, before the PK parameter estimation of the breast DCE-MRI, the signal intensity time series is converted into the concentration time series, and based on the signal intensity of each pixel in the breast tumor, formula (8) is used to calculate the tracer concentration C tissue (t) of the tissue at time t:

[0046]

[0047] Wherein, S0 is the signal intensity before enhancement; based on formula (5) and formula (8), by solving formula (6), the corresponding PK parameter is estimated.

[0048] The substantial feature of the present application is that: in view of the problem that the existing tracer kinetic model is more sensitive to noise, the present application proposes a deconvolution method to enhance the anti-noise performance, so as to more accurately estimate the PK parameters of the breast DCE-MRI. On this basis, the PK parameters of the heterogeneous region are further analyzed to more accurately predict the molecular typing, Ki-67 and pathological grading, and then help doctors to design specific and appropriate treatment plans for breast cancer patients, so as to realize the personalized treatment of breast cancer patients. BRIEF DESCRIPTION OF DRAWINGS

[0049] Figure 1 The complete flowchart of the present application using deconvolution method to estimate PK parameters of breast DCE-MRI for predicting pathological information.

[0050] Figure 2 The boxplot of RMSE of PK parameters of synthetic concentration time series (SNR=30dB) estimated by deconvolution method of the present application.

[0051] Figure 3 The pseudo-color maps of K trans and K trans of the whole tumor region and sub-regions obtained by deconvolution method of the present application.

[0052] Figure 4 The complete flowchart of the present application extracting image features of PK parameters for predicting molecular subtypes, Ki-67 and pathological grade.

[0053] Figure 5 The ROC curves of image features extracted from K trans in the K trans maps on molecular subtypes, Ki-67 and pathological grade of the present application. DETAILED DESCRIPTION

[0054] The technical solutions of the present application will be further specifically explained below by means of specific embodiments and in combination with the drawings.

[0055] Embodiment 1

[0056] DCE-MRI data processing:

[0057] The data used in the present application are all from a third-grade hospital in a province, and a 3.0T magnetic resonance scanner is used to collect breast DCE-MRI data, including one image before enhancement and five images after enhancement, where the images after enhancement are obtained after injection of contrast agent. The detailed scanning parameters used are as follows: repetition time = 4.51 ms, echo time = 1.61 ms, acquisition matrix = 448 x 448, and time resolution = 60 s. For the third image after enhancement, a two-stage tumor segmentation algorithm is used to obtain the tumor region. Here, the two-stage tumor segmentation method is composed of fuzzy C-means and Markov random field algorithms, the fuzzy C-means clustering algorithm is used to determine the preliminary tumor region, and the Markov random field algorithm is used to optimize the preliminary determined tumor region.

[0058] Synthesis of concentration time series:

[0059] Under SNR of 20, 25, 30 and 35, 50 concentration time series with Gaussian noise are synthesized respectively to verify the effectiveness of the deconvolution method, where the sequences all have true PK parameters. K trans and K epvalues are randomly selected from the physiological range intervals 0.05-0.55 and 0.75-2.75, respectively. Here, the PK parameters satisfy the following inequality relationship: K trans < K ep Once we determine the PK parameters, the concentration time series with Gaussian noise can be calculated by equations (1)-(3). The plasma concentration C p (t) at time t is calculated by the following equation (1):

[0060] C p (t) = D-(a1-exp(-m1-t) + a2-exp(-m2-t)), t = 0, 1,..., 5 (1)

[0061] Here, D = 0.2 mmol / kg represents the dose of the tracer, a1= 3.99 kg / L and a2= 4.78 kg / L represent two amplitude constants, m1= 0.1440 min -1 and m2= 0.0111 min -1 represent two rate constants. The tracer concentration C tissue (t) in the tissue at time t is calculated by the following standard Tofts model:

[0062]

[0063] Here, K trans represents the volume transfer constant, K ep represents the rate constant. Finally, the following equation (3) is used to add zero-mean Gaussian noise to the real tracer concentration at time t:

[104] :

[0064] C tissue_s (t) = C tissue (t) + ε(t), t = 0, 1,..., 5 (3)

[0065] Here, the covariance matrix of ε is∑(0.9 |i-j| )σ 2 , σ 2 and SNR satisfy the following equation (4):

[0066]

[0067] Finally, at each SNR, 50 pixels with real PK parameters and concentration time series with Gaussian noise can be obtained.

[0068] Verification of the effectiveness of the deconvolution method:

[0069] The deconvolution method is used to estimate the corresponding PK parameters from the synthetic concentration time series. By deconvolving equation (2), equation (2) can be re-expressed in the following vector-matrix form:

[0070] C tissue = K trans · H(K ep ) C p

[0071]

[0072] Here, H(K ep ) is a Toeplitz matrix and L is the length of the time series. The PK parameters of the concentration time series can be obtained by solving the following optimization problem:

[0073]

[0074] · H(K ep ) C p ||2 (6)

[0075] s.t. 0 < K trans < K ep < 4

[0076] Here, K ep < 4 is defined as the reasonable physiological value because of the existence of the overlarge value (i.e. non-physiological value) in the estimated PK parameters.

[0077] The RMSE index is used to verify the accuracy of the PK parameters estimated by the deconvolution method relative to the true PK parameters, which can be calculated using the following equation (7):

[0078]

[0079] Here, K is the PK parameter estimated by the deconvolution method, K* is the true PK parameter, and N is the total number of pixels in the synthetic image.

[0080] In addition, the Tofts method is also used to estimate the PK parameters of the concentration time series under different SNRs.

[0081] Table 1 RMSE of the PK parameters estimated by the deconvolution method and the Tofts method under different SNRs

[0082]

[0083] Table 1 lists the RMSE of the PK parameters of the concentration time series estimated by the deconvolution method and the Tofts method under different SNRs. Figure 2A box plot of the RMSE of the PK parameters of the concentration time series estimated by the deconvolution method and the Tofts method at SNR = 30 is shown. By comparing the RMSE of the PK parameters estimated by the two methods at different SNRs, it can be found that the PK parameters estimated by the deconvolution method have smaller RMSE compared with the Tofts method. This result shows that the deconvolution method can estimate more accurate PK parameters compared with the Tofts method. In addition, the RMSE of the PK parameters estimated by the two methods presents a downward trend with the increase of SNR. This result shows that the two methods can estimate more accurate PK parameters with the increase of SNR. Our results also show that K ep transit estimated by the deconvolution method has smaller RMSE compared with K trans transit estimated by the Tofts method. These results show that the two methods can estimate more accurate K ep transit compared with K trans transit.

[0084] Application of the deconvolution method in breast DCE-MRI:

[0085] The deconvolution method is applied to the PK parameter estimation of breast DCE-MRI. Before the PK parameter estimation of breast DCE-MRI, the signal intensity time series needs to be converted into the concentration time series. Previous studies have shown that the signal intensity is linearly related to the tracer concentration. Based on this, the relative signal enhancement is used to approximately reflect the tracer concentration of the tissue in the present application. Based on the signal intensity of each pixel in the breast tumor, the tracer concentration of the tissue at t time C tissue (t) can be calculated using the following equation (8):

[0086]

[0087] Here, S0 is the signal intensity before enhancement. Then, based on equation (5) and equation (8), the corresponding PK parameters can be estimated by solving equation (6). Finally, we can get the PK parameters of each pixel in the breast tumor. In addition, we also use the Tofts method to estimate the PK parameters of breast DCE-MRI.

[0088] A TTP-based tumor partitioning method is used to divide the entire tumor into three non-overlapping tumor sub-regions. Here, TTP refers to the time when the tracer concentration reaches the peak value. For each pixel in the breast tumor, the corresponding TTP is calculated by the following equation (9):

[0089]

[0090] Since there are relatively few pixels with TTP=1 and 2, we define the region containing these pixels, along with the region containing the pixel with TTP=3, as the early sub-region. We define the region containing the pixel with TTP=4 and the region containing the pixel with TTP=5 as the mid-term and late-term sub-regions, respectively. For the PK parameters estimated by the deconvolution method, we further extract the PK parameters for each tumor sub-region to analyze the PK parameters of heterogeneous regions.

[0091] Figure 3 The Tofts method and the deconvolution method are listed for estimating the K values ​​of the entire tumor region and tumor subregions. trans The figure shows that K is estimated by both methods. trans They are very similar, with almost no difference. In addition, the early-stage subregion has a higher Kc compared to other tumor subregions. trans Value. Multi-feature analysis based on PK parameters to predict pathological information:

[0092] In this invention, pathological information refers to molecular subtyping, Ki-67, and pathological grading. Molecular subtyping can be divided into four categories: lumihal A, luminal B, basal-like, and HER2 subtyping. Patients with different molecular subtypings have different treatment responses and prognoses. Ki-67 can be divided into two categories: low expression and high expression, reflecting the state of tumor proliferation activity. Tumor grading can be divided into three categories: grade I, grade II, and grade III. Since grade I patients are relatively few, this invention considers grades I and II as low-grade tumors and grade III as high-grade tumors.

[0093] 104 image features were extracted from the PK parameter map to construct a classification model for predicting molecular typing, Ki-67, and pathological grading. All datasets were randomly divided into a training set and a test set. First, five-fold cross-validation was used to select features on the training set. Then, the optimal selected features were used to construct the diagnostic model. Finally, the model's diagnostic performance on the test set for molecular typing, Ki-67, and pathological grading was evaluated. The complete process is as follows: Figure 4 As shown.

[0094] Table 2. K values ​​for tumor subregions and the entire tumor region trans Diagnostic performance in pathological information

[0095]

[0096]

[0097] Table 2 lists the K values ​​for tumor subregions and the entire tumor region. trans AUC in molecular typing, Ki-67 and pathological grading. Figure 5The ROC curves of early sub-region and whole tumor region K trans The ROC curves of early sub-region and whole tumor region K trans The AUCs of early sub-region K trans The AUCs of early sub-region K

[0098] The ROC curves of early sub-region and whole tumor region K trans The AUCs of early sub-region K trans The AUCs of early sub-region K trans The AUCs of early sub-region K

[0098] The ROC curves of early sub-region and whole tumor region K trans The AUCs of early sub-region K trans The AUCs of early sub-region K trans The AUCs of early sub-region K

[0098] The ROC curves of early sub-region and whole tumor region K trans The AUCs of early sub-region K trans The AUCs of early sub-region K trans The AUCs of early sub-region K

[0098] The ROC curves of early sub-region and whole tumor region K trans The AUCs of early sub-region K trans The AUCs of early sub-region K trans The AUCs of early sub-region K

[0098] The ROC curves of early sub-region and whole tumor region K trans The AUCs of early sub-region K trans The AUCs of early sub-region K trans The AUCs of early sub-region K

[0098] The ROC curves of early sub-region and whole tumor region K trans The AUCs of early sub-region K trans The AUCs of early sub-region K trans The AUCs of early sub-region K

[0098] The ROC curves of early sub-region and whole tumor region K trans The AUCs of early sub-region K trans The AUCs of early sub-region K trans The AUCs of early sub-region K

[0098] The ROC curves of early sub-region and whole tumor region K trans The AUCs of early sub-region K

Claims

1. A method of pharmacokinetic analysis based on breast DCE-MRI, characterized in that, The method comprises the following steps: Step 1: Collecting data of breast cancer patients, including image data and pathological information, and screening breast cancer patients with complete data from the data; Step 2: For image data, two algorithms of fuzzy C-means clustering and Markov random field are used for tumor segmentation to obtain tumor regions; Step 3: Synthesizing concentration time series under different SNRs, using deconvolution method to estimate the PK parameters of the synthesized concentration time series, and verifying the effectiveness of the deconvolution method by RMSE index; Step 4: Using the deconvolution method to estimate the PK parameters of breast DCE-MRI to obtain the corresponding PK parameter maps, and further analyzing the PK parameters of the heterogeneous regions; Step 5: Extracting image features from the PK parameter maps of the whole tumor region and the tumor sub-region, constructing a diagnostic model to predict molecular subtypes, Ki-67 and pathological grades; In step 3, SNR is selected as 20, 25, 30 and 35, and 50 concentration time series with Gaussian noise are synthesized for each value; The values of K trans and K ep are randomly selected from the physiological range intervals 0.05-0.55 and 0.75-2.75, respectively, and K trans < K ep ; the concentration time series with Gaussian noise is calculated by formulas (1)-(3); Ct= Ct0e"t / τ p (t): C p (t) = D - (a1- exp(-m1- t) + a2- exp(-m2- t)), t = 0, 1,..., 5 (1) where D = 0.2 mmol / kg denotes the dose of the tracer, a1 = 3.99 kg / L and a2 = 4.78 kg / L represent two amplitude constants, m1 = 0.1440 min -1 and m2 = 0.0111 min -1 represent two rate constants; calculating the tracer concentration C of the tissue at time t tissue (t): where K trans represents the volume transport constant, K ep represents the rate constant; Calculate the real tracer concentration at time t and add zero-mean Gaussian noise; C tissue_s (t) = C tissue (t) + ε(t), t = 0, 1,..., 5 (3) The covariance matrix of ε is ∑(0.9 |i-j| )σ 2 , σ 2 and SNR satisfy equation (4) Get 50 concentration time series with Gaussian noise; Perform deconvolution operation on equation (2) to represent equation (2) as the following vector matrix form: C tissue = K trans · H(K ep ) C p Here, H(K ep ) is a Toeplitz matrix, and L is the length of the time series. Solving the optimization problem of equation (6) can obtain the PK parameters of the concentration time series: Setting K less than 4 ep is a reasonable physiological value; The RMSE index is used to verify the accuracy of the PK parameters estimated by the deconvolution method relative to the true PK parameters: K is the PK parameter estimated by the deconvolution method, K * is the true PK parameter, N is the total number of synthetic pixels.

2. The method of pharmacokinetic analysis based on breast DCE-MRI according to claim 1, characterized in that, In step 4, before performing PK parameter estimation for breast DCE-MRI, the signal intensity time series is converted to a concentration time series, and the tracer concentration of the tissue at time t is calculated based on the signal intensity of each pixel in the breast tumor using equation (8) tissue (t): Where S0 is the signal intensity before enhancement; Based on equation (5) and equation (8), the corresponding PK parameters are estimated by solving equation (6).

Images