一种基于乳腺DCE-MRI的药代动力学分析方法

一种基于乳腺dce-mri的药代动力学分析方法
技术领域
1.本发明属于医学图像处理领域，尤其涉及一种基于乳腺dce-mri的pk分析方法。


背景技术：

2.dce-mri可以用来无创性地评估乳腺肿瘤的结构特性和功能特性，它在乳腺肿瘤的精准诊疗中发挥着非常重要的作用。基于dce-mri的量化分析方法可以划分为两大类，即非参数分析与参数分析。特别地，参数分析需要构建示踪剂动力学模型来估计乳腺dce-mri的pk参数，包括容积转移常数k
trans
和速率常数k
ep
。这些参数可以用来评估血流动力学特性，它们反映了肿瘤的增长速度和生长程度。然而，采集到的图像中往往存在着噪声，这导致传统的示踪剂动力学模型估计的pk参数与真实值之间会存在偏差。除此之外，肿瘤的异质性阻碍了乳腺肿瘤的精准诊疗，异质性区域的pk分析是否可以提高模型在病理信息上的预测性能仍然是未知的。
3.为了便于读者理解，对本发明中出现的英文缩写作以下解释说明：
4.动态增强磁共振成像：英文全称为dynamic contrast-enhanced magnetic resonance imaging，本发明简称为dce-mri。
5.药代动力学：英文全称为pharmacokinetic，本发明简称为pk。
6.去卷积：英文全称为deconvolution，本发明简称为dec。
7.信噪比：英文全称为signal-to-noise ratio，本发明简称为snr。
8.均方根误差：英文全称为root mean square error，本发明简称为rmse。
9.达峰时间：英文全称为time to peak，本发明简称为ttp。
10.受试者操作特征：英文全称为receiver operating characteristic，本发明简称为roc。
11.roc曲线下的面积：英文全称为area under the roc curve，本发明简称为auc。


技术实现要素：

12.为了缓解图像中噪声对pk参数估计的影响，本发明在tofts模型的基础上，提出一种去卷积方法来估计乳腺dce-mri的pk参数。该方法通过将tofts模型中的卷积操作转换矩阵向量相乘的形式以增加抗噪性，从而可以更加准确地估计乳腺dce-mri的pk参数。进一步地，分析了异质性区域的pk参数，以更加精准地预测乳腺癌的病理信息，包括分子分型、ki-67与病理分级。
13.本发明提供一种基于乳腺dce-mri的药代动力学分析方法，包括以下步骤：
14.步骤1：收集乳腺癌患者的数据，包括影像数据与病理信息，并从中筛选出数据完整的乳腺癌患者；
15.步骤2：针对影像数据，采用模糊c均值聚类与马尔科夫随机场两种算法来进行肿瘤分割以获取肿瘤区域；
16.步骤3：合成不同snr下的浓度时间序列，使用去卷积方法估计合成浓度时间序列
的pk参数，并通过rmse指标来验证去卷积方法的有效性；
17.步骤4：采用去卷积方法来估计乳腺dce-mri的pk参数，以得到对应的pk参数图，进一步分析异质性区域的pk参数；
18.步骤5：从整个肿瘤区域与肿瘤子区域的pk参数图提取影像特征，构建诊断模型来预测分子分型、ki-67与病理分级。
19.作为优选，步骤3中，snr选择20、25、30和35，并为每个数值合成50个带有高斯噪声的浓度时间序列；
20.这些序列中k
trans
和k
ep
的值分别从生理范围区间0.05-0.55和0.75-2.75中随机选取且k
trans
＜k
ep
；通过公式(1)-(3)计算得到带有高斯噪声的浓度时间序列；
21.计算t时刻的血浆浓度c
p
(t)：
22.c
p
(t)＝d
·
(a1·
exp(-m1·
t) a2·
exp(-m2·
t))，t＝0，1，
…
，5
ꢀꢀꢀꢀꢀ
(1)
23.其中d＝0.2mmol/kg代表示踪剂的剂量，a1＝3.99kg/l与a2＝4.78kg/l代表两个振幅常数，m1＝0.1440min-1
与m2＝0.0111min-1
代表两个速率常数；
24.计算t时刻的组织的示踪剂浓度c
tissue
(t)：
[0025][0026]
其中k
trans
代表容积转运常数，k
ep
代表速率常数；
[0027]
计算t时刻真实的示踪剂浓度中添加零均值高斯噪声；
[0028]ctissue_s
(t)＝c
tissue
(t) ε(t)，t＝0，1，
…
，5
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0029]
ε的协方差矩阵为∑(0.9
|i-j|
)σ2，σ2与snr满足式(4)
[0030][0031]
得到50个带有高斯噪声的浓度时间序列的像素。
[0032]
作为优选，步骤3中使用去卷积方法从合成浓度时间序列中估计出对应的pk参数；
[0033]
对式(2)进行去卷积操作，将等式(2)重新表示为下面的向量矩阵形式：
[0034]ctissue
＝k
trans
·
h(k
ep
)c
p
[0035][0036]
这里，h(k
ep
)是一个toeplitz矩阵，l是时间序列的长度；
[0037]
求解式(6)最优化问题可以得到浓度时间序列的pk参数：
[0038][0039]
·
h(k
ep
)c
p
||2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0040]
s.t.0＜k
trans
＜k
ep
＜4
[0041]
设定小于4的k
ep
为合理的生理值；
[0042]
采用rmse指标来验证去卷积方法估计出的pk参数相对于真实pk参数的准确性：
[0043][0044]
k是去卷积方法估计的pk参数，k
*
是真实pk参数，n为合成的像素的总数目。
[0045]
作为优选，步骤4中，在进行乳腺dce-mri的pk参数估计之前，将信号强度时间序列转换为浓度时间序列，基于乳腺肿瘤中每个像素的信号强度，采用式(8)来计算t时刻的组织的示踪剂浓度c
tissue
(t)：
[0046][0047]
其中，s0是增强前的信号强度；基于式(5)与式(8)，通过求解式(6)，估计出对应的pk参数。
[0048]
本发明的实质性特点在于：针对现有的示踪剂动力学模型对噪声比较敏感的问题，本发明提出一种去卷积方法来增强抗噪性，以更加准确地估计乳腺dce-mri的pk参数。在此基础上，进一步分析了异质性区域的pk参数，以更加精准地预测分子分型、ki-67与病理分级，进而帮助医生为乳腺癌患者设计特定的、合适的治疗方案，以实现乳腺癌患者的个性化治疗。
附图说明
[0049]
图1为本发明采用去卷积方法来估计乳腺dce-mri的pk参数分析以预测病理信息的完整流程示意图。
[0050]
图2为本发明采用去卷积方法估计的合成浓度时间序列(snr＝30db)的pk参数的rmse的盒型图。
[0051]
图3为本发明采用去卷积方法得到的整个肿瘤区域与肿瘤子区域的k
trans
的伪彩图。
[0052]
图4为本发明提取pk参数的影像特征来预测分子分型、ki-67与病理分级的完整流程示意图。
[0053]
图5为本发明从k
trans
图中提取的影像特征在分子分型、ki-67与病理分级上的roc曲线。
具体实施方式
[0054]
下面通过具体实施例，并结合附图，对本发明的技术方案作进一步的具体说明。
[0055]
实施例1
[0056]
dce-mri数据处理：
[0057]
本发明采用的数据均来自某省三甲医院，采用一个3.0t的磁共振扫描仪来收集乳腺dce-mri数据，包括一个增强前的图像和五个增强后的图像，在这里，增强后的图像是在注射造影剂后采集得到。采用的详细扫描参数如下所示：重复时间＝4.51ms，回波时间＝1.61ms，采集矩阵＝448
×
448，时间分辨率＝60s。针对第三个增强后的图像，采用一个两级的肿瘤分割算法来获取肿瘤区域。在这里，两级的肿瘤分割方法是由模糊c均值与马尔科夫随机场两种算法组成，模糊c均值聚类算法用于确定初步的肿瘤区域，马尔科夫随机场算法
用于优化初步确定的肿瘤区域。
[0058]
浓度时间序列的合成：
[0059]
在snr为20、25、30与35下，分别合成了50个带有高斯噪声的浓度时间序列来验证去卷积方法的有效性，在这里，这些序列均有着真实pk参数。k
trans
和k
ep
的值分别从生理范围区间0.05-0.55和0.75-2.75中随机选取。在这里，pk参数满足下面的不等式关系：k
trans
＜k
ep
。一旦我们确定pk参数，带有高斯噪声的浓度时间序列可以通过等式(1)-(3)计算得到。采用下面的等式(1)计算t时刻的血浆浓度c
p
(t)：
[0060]cp
(t)＝d
·
(a1·
exp(-m1·
t) a2·
exp(-m2·
t))，t＝0，1，
…
，5
ꢀꢀꢀꢀꢀꢀ
(1)
[0061]
这里，d＝0.2mmol/kg代表示踪剂的剂量，a1＝3.99kg/l与a2＝4.78kg/l代表两个振幅常数，m1＝0.1440min-1
与m2＝0.0111min-1
代表两个速率常数。采用下面的标准tofts模型计算t时刻的组织的示踪剂浓度c
tissue
(t)：
[0062][0063]
这里，k
trans
代表容积转运常数，k
ep
代表速率常数。最后，采用下面的等式(3)向t时刻真实的示踪剂浓度中添加零均值高斯噪声
[104]
：
[0064]ctissue_s
(t)＝c
tissue
(t) ε(t)，t＝0,1，
…
，5
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0065]
在这里，ε的协方差矩阵为∑(0.9
|i-j|
)σ2，σ2与snr满足下面的等式(4)：
[0066][0067]
最终，在每个snr下，可以得到50个有着真实pk参数和带有高斯噪声的浓度时间序列的像素。
[0068]
去卷积方法的有效性验证：
[0069]
采用去卷积方法从合成浓度时间序列中估计出对应的pk参数。通过对等式(2)进行去卷积操作，可以将等式(2)重新表示为下面的向量矩阵形式：
[0070]ctissue
＝k
trans
·
h(k
ep
)c
p
[0071][0072]
这里，h(k
ep
)是一个toeplitz矩阵，l是时间序列的长度。通过求解下面的最优化问题可以得到浓度时间序列的pk参数：
[0073][0074]
·
h(k
ep
)c
p
||2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0075]
s.t.0＜k
trans
＜k
ep
＜4
[0076]
在这里，由于估计的pk参数存在超大的值(即非生理值)，因此我们将小于4的k
ep
定义为合理的生理值。
[0077]
采用rmse指标来验证去卷积方法估计出的pk参数相对于真实pk参数的准确性，它可以采用下面的等式(7)来计算得到：
[0078][0079]
在这里，k是去卷积方法估计的pk参数，k*是真实pk参数，n为合成的像素的总数目。
[0080]
除此之外，我们也采用tofts方法来估计不同snr下的浓度时间序列的pk参数。
[0081]
表1去卷积方法与tofts方法估计的pk参数在不同snr下的rmse
[0082][0083]
表1列出了去卷积方法与tofts方法估计不同snr下的浓度时间序列的pk参数的rmse。图2显示了去卷积方法与tofts方法估计出的snr＝30下的浓度时间序列的pk参数的rmse的盒型图。通过比较两种方法估计出的pk参数在不同snr下的rmse可以发现，与torts方法相比，去卷积估计出的的pk参数有着更小的rmse。该结果表明与tofts方法相比，去卷积方法可以估计出更加准确的pk参数。除此之外，随着snr的增加，两种方法估计的pk参数的rmse呈现下降的趋势。该结果表明随着snr的增加，两种方法可以估计出更加准确的pk参数。我们的结果也表明与k
ep
相比，去卷积方法估计出的k
trans
有着更小的rmse。tofts方法估计出的pk参数也可以得到类似的结果。这些结果表明与k
ep
相比，两种方法可以估计出更加准确的k
trans
。
[0084]
去卷积方法在乳腺dce-mri中的应用：
[0085]
将去卷积方法应用于乳腺dce-mri的pk参数估计中。在进行乳腺dce-mri的pk参数估计之前，需要将信号强度时间序列转换为浓度时间序列。以前的研究表明信号强度与示踪剂浓度是线性相关的。基于此，本发明采用相对信号增强来近似地反映组织的示踪剂浓度。基于乳腺肿瘤中每个像素的信号强度，可以采用下面的等式(8)来计算t时刻的组织的示踪剂浓度c
tissue
(t)：
[0086][0087]
这里，s0是增强前的信号强度。然后，基于等式(5)与等式(8)，通过求解等式(6)，可以估计出对应的pk参数。最终，我们可以得到乳腺肿瘤中每个像素的pk参数。除此之外，我们也采用tofts方法来估计乳腺dce-mri的pk参数。基于ttp的肿瘤异质性分析：
[0088]
采用一种基于ttp的肿瘤分区方法来将整个肿瘤划分为三个不重叠的肿瘤子区域。这里，ttp是指示踪剂浓度达到峰值所对应的时间。针对乳腺肿瘤中的每个像素，通过下
面的等式(9)来计算对应的ttp：
[0089][0090]
由于ttp＝1、2的像素较少，我们将其与ttp＝3的像素所在的区域定义为早期子区域，将ttp＝4的像素所在的区域与ttp＝5的像素所在的区域定义为中期子区域与晚期子区域。针对去卷积方法估计的pk参数，进一步提取每个肿瘤子区域的pk参数以分析异质性区域的pk参数。
[0091]
图3列出了tofts方法与去卷积方法估计的整个肿瘤区域与肿瘤子区域的k
trans
图。从图中可以看出，两种方法估计的k
trans
是非常接近的，几乎没有什么差别。除此之外，与其他的肿瘤子区域相比，早期子区域有着更高的k
trans
值。基于pk参数的多特征分析以预测病理信息：
[0092]
本发明中，病理信息是指分子分型、ki-67与病理分级。在这里，分子分型可以划分为四类，分别为lumihal a、luminal b、basal-like与her2分型，不同分子分型的患者有着不同的治疗反应与预后。ki-67可以划分为两类，分别为低表达与高表达，它反映了肿瘤增值活动的状态。肿瘤等级可以划分为三类，分别为i级、ii级与iii级，由于i级患者较少，本发明将i级与ii级均视为低肿瘤等级，将iii级视为高肿瘤等级。
[0093]
从pk参数图中提取104个影像特征来构建分类模型以预测分子分型、ki-67与病理分级。将所有数据集随机划分为一个训练集与一个测试集，首先采用五折交叉验证法在训练集上进行特征选择，然后将选择到的最优特征用来构建诊断模型，最后在测试集上来评估模型在分子分型、ki-67与病理分级上的诊断性能，其完整流程如图4所示。
[0094]
表2肿瘤子区域与整个肿瘤区域的k
trans
在病理信息上的诊断性能
[0095][0096][0097]
表2列出了肿瘤子区域与整个肿瘤区域的k
trans
在分子分型、ki-67与病理分级上的auc。图5显示了早期子区域与整个肿瘤区域的k
trans
的roc曲线。在整个肿瘤区域上，基于去卷积的k
trans
在luminal a、basal-like、her2肿瘤分型和ki-67上有着更高的auc，而基于tofts的k
trans
在luminal b肿瘤分型与病理分级上有着更高的auc。
[0098]
在三个肿瘤子区域中，早期子区域的k
trans
在分子分型、ki-67与病理分级上均得到了最高的auc，分别为0.8319、0.7032、0.7132、0.7490、0.8074与0.6950。进一步地，早期子
区域的k
trans
在分子分型、ki-67与病理分级上的诊断性能要优于整个肿瘤区域的k
trans
的诊断性能。该结果表明异质性区域的pk分析可以提高模型在分子分型、ki-67与病理分级上的诊断性能。