1. 引言

定量磁共振成像(magnetic resonance imaging, MRI)是基于特定的MRI序列和数学建模方法，实现对组织微观结构特征(以下简称微结构)的无创性量化分析，对揭示肿瘤的异质性具有重要价值[1]。扩散加权成像、弛豫成像等定量技术是评估组织微观结构的常用MRI方法[2]-[5]。基于同时改变扩散编码参数(如b值)以及弛豫编码参数(如回波时间TE)的混合扩散–弛豫MRI技术[6]，可以提供MRI扩散和弛豫之间的相关性和耦合性，能够分离和量化多个相关成分，为肿瘤的异质性评估提供新视角，在大脑[7]-[11]，前列腺[12]-[16]、胎盘[17] [18]等部位已有相关报道。实现快速鲁棒的扩散弛豫分析方法是虚拟病理特征预测的热点问题之一。胰腺癌(pancreatic cancer, PC)作为一种恶性程度极高的消化系统肿瘤，典型病理特征为丰富的纤维成分和压缩的血管，表现出有别于多数实体肿瘤如乏血供的影像特征和复杂微环境，实现无创性组织成分如上皮、间质和管腔的评估对创新药物开发、肿瘤诊断等具有重要意义[19]。

分离混合信号成分是一种经典的逆问题，不同的数据分析方法根据构建的不同模型假设，如体素分量的个数、多个分量的对比度变化函数的形式等进行区分。对于最通用的一些解混方法，包括非负矩阵分解[20]，独立成分分析[21]和低秩张量分解[22]等，这些方法仅适用于分量个数相对较少且假设的对比度变化函数相对简单的情况，可以将数据分解为有意义的分量，已应用于混合扩散–弛豫数据分析，但该方法具有一定的局限性，如可能将数据分解成缺乏生物学意义的成分，且分解不稳定。另外，由于混合扩散–弛豫空间需较长的TE，因此这些方法可能会忽略短T₂分量。混合扩散–弛豫的每个分量都与独立的单指数扩散和弛豫衰减相关联，连续体模型是分析混合扩散–弛豫数据的技术之一[23]，其逆问题即逆拉普拉斯变换(inverse Laplace transform, ILT)，通常采用最小范数最小二乘法(minimum-norm least-squares, MNLS)作为解决方案，但模型中出现的指数衰减特性会导致逆问题严重病态且高度不稳定。虽然通过施加一定的约束方法能获得合理的ILT解决方案，如非负性约束[24]、正则化约束[25] [26]、边际分布约束[25]等，但这些方法存在依赖于任意选择的正则化项等不足。

本研究采用一种数据驱动的方式克服以上限制，通过假设固定的分量个数，将不同b值和不同回波值对应的整体数据的体素信号视为高斯混合模型，利用期望最大化算法对混合信号中各个成分的权重进行迭代求解，获得反应组织上皮、间质和管腔成分的MRI参数图。并通过肿瘤的病理分析与MRI结果的比较，探究所建立模型的可行性。

2. 材料与方法

2.1. 研究对象

人源胰腺癌细胞(PANC-1)和人胰腺癌相关成纤维细胞(pancreatic cancer-associated fibroblasts, pCAFs)购买自上海葵赛生物科技有限公司。根据细胞密度和生长情况，每两到三天更换一次培养基，当细胞密度达到80%至90%时，进行传代以避免过度生长。待细胞量达到实验要求后，将PANC-1与pCAFs按数量1:1或1:5混合制成细胞悬液，用于建立肿瘤间质成分比例不同的胰腺癌裸鼠模型。用1 mL注射器吸取细胞悬液(0.3 ml)，注射到5周龄的BALB/c裸鼠皮下(裸鼠源于斯贝福(苏州)生物技术有限公司)，待肿瘤长径至0.8 cm左右时用于MRI实验。所有动物实验均由中国上海海军军医大学第一附属医院动物福利与伦理小组批准。

2.2. MRI成像

采用BioSpec 11.7-T磁共振成像仪对构建的裸鼠模型进行扫描，其具备11 cm的孔径宽度，740 mT/m最高梯度强度，6600 T/m/s梯度爬升率。扫描前，标记裸鼠，采用异氟烷吸入性麻醉药剂对裸鼠进行麻醉，气体流量调节为300 ml/min，浓度调节为3%~4%，维持约3分钟，等待裸鼠完全麻醉，进行MRI扫描。主要扫描序列及参数如下：

(1) 单次激发自旋回波–平面回波成像序列(spin-echo single-shot echo-planar imaging sequence, SE-EPI)进行多回波、多b值扩散加权成像(diffusion weighted imaging, DWI)。扫描参数：重复时间TR (3000 ms)，4个回波时间TE (分别是25、50、75、100 ms)，4个b值(分别是0、150、7500、1500 s/mm²)；10层，层厚/层间距：1 mm/1 mm，成像视野FOV：30 mm * 30 mm，采集矩阵：120 * 120，重复采集次数：3次。

(2) 采用多回波自旋回波序列(multi-echo spin-echo, MESE)实现T2_mapping成像。具体参数如下：重复时间(TR: 3600 ms)，30个回波时间(TE：7.5~225 ms，间隔7.5 ms)；15层，层厚/层间距：0.5 mm/0.5 mm，FOV：30 mm * 30 mm，采集矩阵：150 * 150，重复采集次数：1次。

MRI扫描过程中，连续使用异氟烷/氧气混合气体维持裸鼠的麻醉状态，浓度调整为1.0%~2.0%，持续监测裸鼠的呼吸，并使用加热垫维持裸鼠体温，以确保实验动物的生理状态稳定。

2.3. 病理分析

荷瘤鼠接受MRI检查后进行全肿瘤切除，并将肿瘤组织固定于福尔马林中，制作病理蜡块和切片。使用苏木精–伊红(hematoxylin-eosin, H&E)对肿瘤最大层面染色。为验证算法结果的准确性，我们对病理图像进行了不同组织成分的定量评估。采用Image J软件，利用StarDist工具[27]对HE染色区域进行细胞核自动分割，利用阈值法设置阳性与阴性细胞范围，并由病理医生确定各组织成分的面积，计算基于病理图像的上皮、间质和管腔成分所占比例[28]。

2.4. 扩散弛豫数据分析理论

将混合扩散–弛豫模型视为高斯混合模型，利用期望最大化(expectation-maximization, EM)算法[29]来估计模型参数，从而实现对组织微结构的表征。

2.4.1. 扩散–弛豫模型

对于一维扩散或弛豫模型，通常采用单指数衰减形式，其连续体模型的理想信号分别由下式给出

$S\left(TE\right)={\displaystyle \int F\left(T_2\right)*{\text{e}}^{-TE/T_2}\text{d}{T}_2}$ (1)

$S(b)={\displaystyle \int F\left(D\right)*{\text{e}}^{-bD}\text{d}D}$ (2)

我们将其中的分布函数 $F\left(T_2\right)$ 和 $F\left(D\right)$ 称为一维谱。

将一维方法进一步扩展形成高维方法。对于单个体素而言，其理想多维信号如下：

$S\left(b,TE\right)={\displaystyle \iint F\left(D,T_2\right)*{\text{e}}^{-bD}{\text{e}}^{-TE/T_2}}=KF$ (3)

其中 $K={\text{e}}^{-\text{bD}}{\text{e}}^{-\text{TE}/T_2}$ 。 $F\left(D,T_2\right)$ 是一个二维分布函数，表示扩散系数D和弛豫时间T₂的联合分布，我们称之为二维扩散–弛豫相关谱。

根据高斯混合模型的原理，我们假设整个体素信号有M个单高斯模型，每个模型有相应的概率分布及参数。但在实践中，概率分布是未知的，为隐变量，需要对其进行估算，因此这些体素隶属度指标是模型的潜在状态，若第n个体素属于第m个高斯模型，我们表示为 $Z_{nm}$ 。因此，完整的体素权重集为

$Z_n={\left\{Z_{nm}\right\}}_{m=1}^M,{{\displaystyle \sum }}_{m=1}^M{Z}_{nm}=1$ (4)

根据式(3)可得，单个体素相关谱可以由光谱分量的加权和给出：

$F\left(Z_n\right)={{\displaystyle \sum }}_{m=1}^M{Z}_{nm}{F}_m$ (5)

在考虑存在高斯噪声的情况下，单个体素的离散模型满足以下条件：

$S_n~N\left(KF\left(Z_n\right),{\sigma }_n^{2}I\right)$ (6)

其中令所有观测值都具有相同的方差 ${\sigma }_n^2$ ，选取b = 0 s/mm²、TE = 25 ms时的图像体积方差来估计协方差矩阵 ${\sigma }_n^{2}I$ 。从单体素衍生到整个图像(N个体素)，完整的数据集为 $D={\left\{S_n\right\}}_{n=1}^N$ ，则本研究建立的高斯混合模型所要推断的参数集为：

$\theta =\left\{F_1,F_2,\cdots ,F_M,{\left\{Z_{n1},Z_{n2},\cdots ,Z_{nM}\right\}}_{n=1}^N,{\left\{{\sigma }_n^2\right\}}_{n=1}^N\right\}$ (7)

为准确推断参数 $\theta$ ，我们引入了一种无监督学习策略——期望最大化算法，解决含有隐变量的混合模型参数估计问题。

2.4.2. 期望最大化算法

期望最大化算法的具体流程如图1所示。

Figure 1. The flow chart of the EM algorithm for solving GMM parameters

图1. 期望最大化算法求解GMM参数流程图

本研究中实现EM算法计算参数的具体流程如下：

(1) 初始化参数 ${\theta }^0$ ：

首先初始化光谱分量 $F_m$ ，权重系数 $Z_{nm}$ 。

(2) E步(计算隐变量的期望值)：

计算完整数据的对数似然函数，其中完整数据为 $\left\{D,Z_n\right\}$ ：

$\begin{array}{l}\log \left(L\left(\theta \right)\right)=\log P\left(D,Z_n|\theta \right)=\log P\left({\left\{S_n\right\}}_{n=1}^N,Z_n|F_1,F_2,\cdots ,F_M,{\left\{Z_{n1},Z_{n2},\cdots ,Z_{nM}\right\}}_{n=1}^N,{\left\{{\sigma }_n^2\right\}}_{n=1}^N\right)\\ =\underset{n=1}{\overset{N}{{\displaystyle \sum }}}\log N\left(S_n;KF\left(Z_n\right),{\sigma }_n^{2}I\right)\\ =\underset{n=1}{\overset{N}{{\displaystyle \sum }}}-\frac{1}{2}\log 2\pi {\sigma }_n^2-\frac{1}{2{\sigma }_n^2}{\Vert S_n-KF\left(Z_n\right)\Vert }_2^2\end{array}$ (8)

计算对数似然函数关于在给定观测数据和当前参数 ${\theta }^t$ 下对未观测数据Z的条件概率分布的期望，即Q函数。根据琴生(Jensen)不等式可得：

$Q\left(\theta ,{\theta }^{\left(t\right)}\right)=E_{z|D,{\theta }^{\left(t\right)}}\left[\log \left(L\left(\theta \right)\right)\right]=\underset{z}{{\displaystyle \sum }}P\left(Z_n|D,{\theta }^{\left(t\right)}\right)\log P\left(D,Z_n|\theta \right)$ (9)

其中 $\underset{z}{{\displaystyle \sum }}P\left(Z_n|D,{\theta }^{\left(t\right)}\right)=1$ 。

(3) M步(最大化Q函数)：

求使 $Q\left(\theta ,{\theta }^{\left(t\right)}\right)$ 极大化的θ，使得下一步迭代的 ${\theta }^{t+1}=\theta$ 。

${\theta }^{\left(t+1\right)}=\underset{\theta }{\arg \max }Q\left(\theta ,{\theta }^{\left(t\right)}\right)$ (10)

在最大化Q函数迭代过程中，可以得到对参数的迭代估计，如下：

${\left\{F_m\right\}}_{m=1}^M=\underset{F_m\ge 0}{\arg \min }{\Vert {{\displaystyle \sum }}_{n=1}^N\left(K{Z}_{nm}\right)F_m-{{\displaystyle \sum }}_{n=1}^N{S}_n-K\left(\underset{m\ne j}{{\displaystyle \sum }}{Z}_{nm}{F}_m\right)\Vert }_2^2$ (11)

${\left\{Z_{nm}\right\}}_{m=1}^M=\underset{ {{\displaystyle \sum }}_{m=1}^M{Z}_{nm}=1}{\arg \max }\log N\left(S_n;KF\left(Z_n\right),{\sigma }_n^{2}I\right)$ (12)

分别利用非负最小二乘法和内点法进行求解。

重复(2)和(3)步骤，直到函数收敛或达到最大迭代次数。我们设置收敛到容差为10^-3或完成5~10次迭代后，计算结束。

3. 结果

3.1. 数据处理

利用Python将所得的T2 mapping以及多回波多b值dicom格式图像批量转换成nii格式。对于扩散-弛豫数据，将总体拆分成每个b值和TE组合对应的图像，共16种对比度编码。如图2所示。

将4b × 4TE的数据转换成四维矩阵 $\left(x,y,z,j\right)=\left(120,120,10,16\right)$ ，其中x，y分别表示图像的维度，z表示扫描层数，j表示b值和TE的组合数。

Figure 2. DWI of pancreatic cancer with multiple b- and TE values

图2. 多b值多TE胰腺癌DWI图像

3.2. 初始值设定

根据病理图像，可假设肿瘤组织有3种不同成分：上皮、间质和管腔，则将M设置为3。在进行参数求解时，需要进行参数初始化。根据式(3)，在已知体素信号情况下，通过设置T₂与D值的初始值，可得初始化的K值。为减少误差，利用T2 mapping图像进行三高斯拟合，得到三个单高斯分布，对应不同组织成分特性，具体处理如下：

利用T2 mapping图像的整体灰度直方图，采用Origin软件对数据进行Gauss函数多峰拟合分析，可得三种成分高斯分布的均值，其作为不同组织成分的T₂初始值[39, 63, 90] ms。如图3所示。

Figure 3. Gaussian fitting results of T2 values

图3. T2值高斯拟合结果

对于上皮、间质和管腔三种不同成分的初始值D，根据经验[30]设置为[0.5, 1.0, 2.7] µm²/ms。

通过T₂、D值可得初始化的K值，根据式(3)，利用非负最小二乘法得到整个图像的相关谱分布F，将3个最大的独立分布作为初始光谱分量 $\left\{F_1,F_2,F_3\right\}$ 。利用式(5)，基于光谱分量，初始化权重系数 ${\left\{Z_{nm}\right\}}_{n=1}^N$ 。

3.3. 数据分析

基于Matlab2019b平台，分别将两例多b值多回波数据及T₂、D的初始值作为输入，集成到算法中，算法自动初始化光谱分量及权重系数，再利用EM算法不断迭代优化这两个参数值，直至收敛至稳定状态，获得胰腺癌不同组织成分相应权重图(图4)以及肿瘤区域均值(表1)。

Figure 4. Weight maps of different tissue components in two tumors

图4. 两例肿瘤对应不同成分的权重图

Table 1. Mean volume fractions of different tissue components in two tumors

表1. 两例肿瘤不同组织成分体积分数的均值

注：所得数值对应不同组织所占的体积分数(%)。ROI：感兴趣区域。

3.4. 病理结果比较

基于肿瘤最大截面切片HE染色图像，对两例肿瘤进行上皮、间质和管腔分割，分割过程如图5所示，其对应的比例见表2。

Figure 5. (a): Original pathology image; (b): Partial region of pathology; (c): Segmentation of cell nuclei to estimate cell distribution; (d): Tissue components are divided into epithelium (yellow = nuclei, gray = cytoplasm), stroma (blue), and lumen (green)

图5. (a)：原始病理图；(b)：病理部分区域；(c)：对细胞核进行分割以估计细胞分布；(d)：组织成分分为上皮(黄色 = 细胞核，灰色 = 细胞质)、间质(蓝色)、管腔(绿色)

Table 2. Pathological analysis of the percentage area occupied by different tissue components (%)

表2. 病理分析不同组织成分面积占比(%)

4. 讨论

磁共振混合扩散–弛豫的研究最初受限于图像质量和成像时间，主要用于脑部成像，并在病理结果预测方面显示出潜在价值。随着硬件和软件的进步，该技术的应用已扩展至体部器官，如前列腺[12]-[16]和胎盘[17] [18]中已有相关报道。Chatterjee等人[16]利用混合多维MRI技术，通过耦合每个组织成分的T2和ADC值，定量估计前列腺及肿瘤的组织成分，用于前列腺癌的诊断和侵袭性预测。Slator等人[17]采用ZEBRA [31]采集策略和最小振幅能量正则化的逆拉普拉斯变换方法，在单一重复时间(TR)内对多个TE和扩散编码进行采样，估计T2-ADC频谱，探索了胎盘的扩散和弛豫特性，并展示了可能反映多个组织微环境结构的二维峰。这些研究采用的分析技术均基于逆拉普拉斯变换原理，不同之处在于优化算法和正则化处理方式的差异。例如，在采用边缘分布[25]、空间正则化[6]等方面需要高信噪比的数据，且需要手动划分光谱域，这在图像存在非均匀区域时容易导致定义不清的问题。为解决这一问题，本研究采用数据驱动的方式，利用高斯混合模型对具有相似光谱的体素进行软聚类，自动划分多维MR图像，带来了潜在优势。通过对相似体素进行平均，减少了规范频谱估计中的噪声。此外，该方法完全为数据驱动，不易错过数据中的任何重要光谱分量，且能够持续捕获数据中的内在变化以防止过拟合。

针对微结构分析，目前存在一些先进的定量MRI技术，包括基于T2弛豫技术的腔内水成像(Luminal water imaging, LWI) [32]和基于扩散的技术如肿瘤细胞数中的血管、细胞外和受限扩散(VERDICT) [33]、限制谱成像(Restriction spectrum imaging, RSI) [34]和张量值扩散编码[35]，以及本文研究的混合弛豫–扩散技术。其中混合弛豫扩散技术包含多个方向，如混合多维MRI (Hybrid multi-dimensional MRI, HM-MRI) [12] [13]，时变扩散成像(Time-dependent diffusion) [36]，和扩散–弛豫相关光谱成像(Diffusion-relaxation correlation spectrum Imaging, DR-CSI) [6]等。本文主要采用其中的HM-MRI技术，其基于组织MRI信号中ADC和T2之间的相互依赖性，并利用一种软聚类模型及迭代求解的方式进行多维数据拟合，一方面该技术扫描时间相对较短、参数设置较为简便，另一方面，从结果来看可得到较好的临床转化可能。本研究中，利用了高场强MRI的优势，获得高质量图像和信噪比，使得数据分析技术在鲁棒性方面有很好的保证。

本研究利用动物模型和专用的MRI设备进行实验，并引入了高斯混合模型以及期望最大化算法，实现混合扩散弛豫的MRI参数自动化计算。尽管目前仅基于两例肿瘤组织进行分析，但本研究重点关注的是混合弛豫扩散模型算法可行性研究，仅需学习单一肿瘤的扩散–弛豫数据即可实现对应参数自动化估计。另外，从获得参数的结果来看，对两种不同间质成分的胰腺癌病理成分上已表现出明显差异，与本研究所建立的扩散弛豫分析模型计算结果有初步相似性。本研究也存在一定的局限性，需增加数据集的规模，进一步验证模型的稳定性和泛化能力，同时可对数据分析结果与病理结果之间进行相关性分析。

5. 结论

本研究成功建立了胰腺癌裸鼠模型，并实现高分辨MRI成像。通过无监督学习技术结合人工智能算法，实现了混合扩散–弛豫MRI的组织成分比例分析，为胰腺癌中包括上皮、间质和管腔三种病理成分比例的预测提供了新途径，有望用于胰腺癌异质性研究。未来工作将聚焦于算法的进一步优化，包括开发创新的初始化策略、改进收敛标准、并行化算法处理等。通过这些工作，我们期望为胰腺癌的早期诊断和治疗提供更可靠的科学依据，以改善患者的临床治疗效果，同时为基础研究和临床应用提供方法学参考。

NOTES

^*通讯作者。

参考文献