参数计算方法 技术领域

本发明涉及可示踪物质在活体器官组织间隙中 分布及细胞外间隙中 扩散的参数计算方法， 特别一种借助 MRI图像分析和计算可示踪物质在 活体器官中扩散参数的计算方法。

背景技术

药物代谢动力学研究中， 通常为了评价药物的扩散性能、 预测药物 在靶部位发挥的疗效， 需要在给药后的不同时间点处死、 解剖动物获得 组织样本或者采集健康志愿者的静脉血， 测量局部组织或其血样中的药 物浓度。 这些传统的研究方法不仅会对受试对象造成一 定的伤害， 而且 也仅能得到药物的宏观分布规律， 微观领域的测量参数无法获得。 此外， 进行离体测量前， 样本的保存和转移也需耗费一定的人力物力。

随着影像学技术的飞速发展， 放射性示踪法为实现药物的在体研究 和可视化测量提供了可能。但是， 碘剂等示踪剂需要依赖 X线、 CT等成 像设备， X线、 CT机和放射性同位素带来的电离辐射和电离损� ��会对受 试对象造成一定的危害。 .

磁共振成像 (magnetic resonance imaging, MRI)是近年来最常用的成 像检测技术， 用这种技术来观察人体或动物解剖结构与生理 功能具有实 时、 活体、 可视、 无创的优势。 MRI示踪剂的使用， 更加拓展了 MRI的 应用范围。 ' 目前， M I检查中主要应用两种可示踪物质： 一种是以钆喷酸葡胺 (Gd-DTPA )为代表的 T1阳性示踪剂， 另一种是以铁纳米微粒为代表的 T2阴性示踪剂。已有学者应用铁纳米微粒作为 MRI示踪剂对脑组织间液

1

替换页 （细则第 26条) ( interstitial fluid, ISF) 中代谢物的清除途径进行了初步研究， 证实了注 入脑内的铁纳米微粒经鼻粘膜处的淋巴最终进 入颈部淋巴结清除出脑。

然而， 研究结果也显示， 由于纳米颗粒铁磁性对梯度磁场的干扰， 图像产生变形， 并导致了大面积的信号缺失， 无法清楚显示示踪剂在活 体器官的准确观察和定量分析。

还有研究利用磁共振扩散加权成像 （diffusion weighted imaging, DWI)进行活体器官的测量。 这是一种用于测量水及体内其它分子的表 观扩散系数 （ apparent diffusion coefficient, ADC ) 和组织各向异性分数

(anisotropy fraction, FA) 的 MRI技术。 该方法是基于分子扩散变化导致 组织 MRI图像信号改变的原理，在某一方向上施加一 个扩散敏感线性梯度 场， 若组织内部某一位置在此方向上的扩散明显， 则采集到的磁共振信 号强度降低， 反之， 信号强度增加。 根据不同扩散敏感梯度下得到的不 同 MRI信号强度， 可计算得到 ADC值。 通过施加六个以上不同方向上的 扩散敏感梯度， 还可获得扩散张量特征参数， 如各向异性分数（fractional anisotropy, FA)扩散个各向异性与整个扩散的比值，代表了 扩散分子在扩 散主向量轴上的运动强度。

因此， 开发出在活体状态下可用于活体器官的扩散参 数指标定量监 测的方法， 对研究物质在该解剖部位的代谢， 以及未来利用该解剖部位 进行临床药物治疗， 都具有重大的理论和应用价值。

中国专利申请 CN2010102086244涉及一种脑组织间液与脑细胞外间 隙（extracellular space, ECS ) 的生理参数测量方法。 该技术利用小分子 示踪剂 Gd-DTPA的磁弛豫特性， 通过磁共振成像获取 Gd-DTPA在脑间质 内的三维空间分布和清除情况， 再应用经典扩散方程反演得到脑 ECS的 解剖和生理学参数。 但该专利中并没涉及如何进一步准确定量的分 析示 踪剂在脑 ECS中扩散的时空分布。 发明内容

为此， 本发明旨在提供可示踪物质在活体器官的扩散 参数的计算方 法， 该方法可以利用可示踪物质在活体器官中扩散 而引起的 MRI信号强 度变化， 计算得到可示踪物质在活体器官中的扩散性能 ， 以直接反映可 示踪物质在活体器官中因扩散而导致的分布变 化与清除过程。

进一步来说， 本发明利用 MRI示踪技术， 对获取的动态三维图像进 行区域分割和配准处理， 随后采用各向异性扩散方程， 最终得到活体器 官空间结构参数以及示踪物质在其中的扩散、 清除参数， 并可提供分子 的分布容积以及整体清除等药代动力学参数。

在可示踪物质在活体器官的扩散参数的计算方 法的示意性实施方式 中， 可示踪物质为 Gd-DTPA， 活体器官为脑 ECS。

本发明提供了一种可示踪物质在活体器官中的 扩散参数的计算方 法， 包括：

( 1 ) 在琼脂糖中的环境中测量可示踪物质不同浓度 （C)相应的 磁共振信号增量 （ASI);

(2) 根据步骤 （1 ) 测定的值确定磁共振信号增量 （ASI) 与可 示踪物质浓度（C)具有线性对应关系时的示踪� ��浓度上限， 并确实磁共 振信号增量 （ASI) 与可示踪物质浓度 （C ) 的线性对应关系；

(3 ) 将不同浓度的可示踪物质从注入点源（Q)注入� ��体器官中， 随后在磁共振设备中获得不同时刻相应的磁共 振图像；

(4) 对获得的磁共振图像中感兴趣的区域顺序进行 图像配准；

( 5 ) 根据经配准的图像对可示踪物质在活体器官扩 散进行参数 拟合， 它包括：

( 5.1 ) 确定可示踪物质的扩散方程，

( 5.2 ) 根据可示踪物质的浓度与 MRI信号增量的关系， 对步骤 (5.1)确定的方程进行换算，

(5.3) 对步骤 （5.2) 得到的方程求解， 以拟合出包括可示踪物质 的扩散系数 D*、 损失比率 和 /或局部体积比。的扩散参数。

在可示踪物质在活体器官中的扩散参数的计算 方法的再一种示意性 实施方式中， 其中：

步骤 （5.1) 确定的可示踪物质从注入点源 （Q) 扩散至空间内任意 一点 （Ρ) 的扩散方程为：

dc_ D_d .

其中： D为可示踪物质的自由扩散系数，

r表示计算点 （P) 到示踪剂注入点源 （Q) 的距离， C为可示踪物质的浓度， 由磁共振信号增量（ASI)计算得到， t为注入可示踪物质与磁共振图像采集之间的� �间间隔， α为细胞外间隙占活体器官全部容积的比例，

k'为可示踪物质在活体器官中的消除速率常数� �� λ _Ρ(5 为计算点（Ρ)到示踪剂注入点源（Q)之间� ��体器官的迂 曲度， 反映活体器官空间结构对分子扩散运动的阻碍 作 用， λ _Ρ(5 =λ (q,j) ,q和 j分别为示踪剂注入点源 （Q) 与 计算点 （P) 处连线分别与 x、 z轴的夹角； (5.2) 确定的转换后的方程为： d(ASI) _ D d ² (ASI) Q

c ¹ a 其中： k、 b为常数，

△SI为磁共振的信号强度值，

d表示单个像素对应的实际空间距离；

步骤 （5.3 ) 中对上述方程数学求解可以得到方程的解为:

且

其中： π和 e为数学常数，

ξ为磁共振图像上单个像素对应的实际空间距� ��。

在可示踪物质在活体器官中的扩散参数的计算 方法的另一种示意性 ；施方式中， 其中：

步骤 （5.1 ) 确定的可示踪物质在标准正交坐标系内的扩散 方程是：

• V C -÷—- .¾' C

i .¾( y, z )* a

其中： D为可示踪物质的自由扩散系数，

λ(χ,γ,ζ)为活体器官中某一点 ( _X ,y,z:>的迂曲度，

C为可示踪物质的浓度， 由磁共振信号增量 （ASI ) 计算得到， t为注入可示踪物质与磁共振图像采集之间的� �间间隔， α为细胞外间隙占活体器官全部容积的比例，

k'为可示踪物质在活体器官中的消除速率常数� ��

该扩散方程在各向异性的条件下可等价变换为 ：

C一 D d ² c ,' D C D d ² c < o〜

a 其中： D为可示踪物质的自由扩散系数，

、 、 λ _ζ 分别为活体器官在 x、 y、 z标准轴上的迂曲度， C为可示踪物质的浓度， 由磁共振信号增量 （ASI) 计算得到， t为注入可示踪物质与磁共振图像采集之间的� �间间隔， α为细胞外间隙占组织或器官全部容积的比例� ��

k'为可示踪物质在活体器官中的消除速率常数� ��

步骤 （5.2) 确定的转换后的方程为：

D —率 _+k ' _b

dt λ(χ,ν,ζ)' o a 其中： k、 b为常数，

Δ SI为磁共振的 MRI信号增量值，

d表示单个像素对应的实际空间距离；

步骤 （5.3) 对上述方程可转化为离散模型求解如下：

AS! ( + Δ/)― Δ57( = -- V ² (Μ/ ) + k― k ' AS1 ( ή + Α· ' b 其中： At为单位采样时间，

入（x,y.z)为空间内任意一点的迂曲度，

Δί表示任意方向 （i轴） 上的空间采样距离。

在可示踪物质在活体器官中的扩散参数的计算 方法的还一种示意性 实施方式中， 其中上述步骤 （5.3) 得到的方程可以采用如下最小二乘法 拟合： - k'(0Sr}+k'h-~ASJ ( + A j + AS/ (i) 其中^ ^')表示方差的最小值， 上述公式用矩阵形式表示为：

其中， 矩阵 A的第 （x,y，z，t，：）行的非零元为

h(x.y, zJ,:) '： Δ 7(/ + Αή- A ! (ί)

已知向量 b的第 （x,y，z, t，：） 行为：

h ( .v.、■., :.'·:)■- AS1 ( / + Δ/ ) - A.S7 ( I ) 替换变量 X为:

在可示踪物质在活体器官中的扩散参数的计算 方法的另一种示意性 实施方式中， 其中图像配准步骤包括： ·3

=έ!>( "^

其中： H _tn ， 〃ω和 / ⁽ /， _g )分别为两幅图像的边缘熵以及联合熵；

W和 分别为两幅图像中灰度值样本的个数；

P«， '】和 Ρ ⁽ ''， 表示灰度值出现的概率。

在可示踪物质在活体器官中的扩散参数的计算 方法的又一种示意性 实施方式中， 其中步骤 （4) 的配准步骤还包括如下的灰度校正：

I(x-. y ) = f _MP ( -...y .z)B(x.. y,∑) +n(x,y._z)

其中： 是未受噪声和非均匀场污染的真实磁共振信号 ；

为成像信号；

3为灰度非均匀场；

"为成像过程中引入的加性噪声。

在可示踪物质在活体器官中的扩散参数的计算 方法的又一种示意性 实施方式中， 其中的灰度校正进一步包括： 其中： （"P ^和 ("2， )是两个待校正区域内磁共振信号强度的均值� � 方差， 且 g和 ^分别为实测和校正后的磁共振图像信号强度� � 在可示踪物质在活体器官中的扩散参数的计算 方法的一种示意性实 施方式中， 可示踪物质为 Gd-DTPA， 动物组织器官的扩散参数包括 Gd-DTPA在脑细胞外间隙的有效扩散系数 D *、 迂曲度 A和损失比率 。

采用本发明的方法， 可以在活体状态下同时得到活体器官空间的反 应组织细胞间结构和生理特点的扩散参数和反 映组织液流动性以及示踪 分子在脑内经组织通道分布与清除过程的扩散 参数， 即该方法即可以计 算出包括有效扩散系数 *、消除速率常数 ：'、迂曲度 等微观生理参数， 也可以计算出如分布容积 Vd、 分布容积分数等宏观生理参数， 并在理论 上可以实现对疾病治疗给药方案的设计， 有助于指导临床用药。 例如， 可以对 Gd-DTPA在脑 ECS中的扩散与清除情况进行计算， 为脑部微循 环的研究提供有价值的分析数据。 附图说明

图 1显示了在 MRI设备中得到的不同浓度 Gd-DTPA的 MRI图像。 图 2是 Gd-DTPA在 0〜100mM浓度范围内， Gd-DTPA的浓度 C与 MRI信号强度增量 ASI的对应关系。

图 3显示了 Gd-DTPA在 0〜4mM浓度范围内， Gd-DTPA的浓度 C与 MRI信号强度增量 ASI的对应关系.

图 4显示了用腕线圈成像时， Gd-DTPA相关的信号强度增量 ASI与 Gd-DTPA的浓度 C处于 0〜1.5mM之间的线性关系。

图 5显示了用头线圈成像时， Gd-DTPA相关的信号强度增量 ASI与 Gd-DTPA的浓度 C处于 0〜1.5mM之间的线性关系。

图 6示意性地显示了处理以上得到的 MRI时间序列图像的方法。 图 7显示了在 MRI成像过程中大鼠的运动弓 I起大部分组织明显的形 变。

图 8用以说明本发明 MRI时间序列图像的配准框架。

图 9显示了使用 Matlab软件测量对减后的未增强和增强后的区域 。 图 10是图 9所示中增强区域 MRI信号测量后生成的测量结果图， 横坐标是信号强度值， 纵坐标是像素个数。

图 11 显示了 SD 大鼠脑内不同脑区注入 2μί浓度为 10 mM 的

Gd-DTPA示踪剂溶液， 在不同时间扫描后采集的图； 其中 （a) 的注入点 源位于尾状核区， （b) 的注入点源位于黑质纹状体， （c) 位于丘脑， （d) 位于枕叶皮层， （e)位于皮质下白质。

图 12是 Gd-DTPA注入大鼠骨骼肌间隙内在不同时间 MRI扫描后采 集的扩散图像。

图 13是示踪剂 Gd-DTPA在胶质瘤模型大鼠脑内的扩散图像， 其中 (a) 的 T1WI脑内高信号区域为 Gd-DTPA的扩散区域， 低信号区域为 瘤体部位， （b) 的 T2WI高信号区域为瘤体部位。

具体实施方式

为了对本发明的技术特征、 目的和效果有更加清楚的理解， 现对照 附图说明本发明的具体实施方式。

下面将结合一个具体的实施方式， 说明活体器官的生理参数的计算 方法。 在该实施例中， 示踪剂选用钆喷酸葡胺 （Gd-DTPA)， 即由拜耳先 灵药业生产的马根维显注射液。 测量琼脂糖环境中和大鼠脑部的的示踪 剂实际浓度下相应的 MRI信号强度。 成像动物为 250g〜300g SD雄性大 鼠，计算的活体器官为大鼠的大脑。大鼠脑部 的 MRI成像机选用 Siemens Magnetom trio a tim 3.0T磁共振机。

图 1显示了在 MRI设备中得到的不同浓度 Gd-DTPA的 MRI图像， 其中显示了在 0~100mM范围内，不同浓度的 Gd-DTPA注入 37°C的琼脂 糖后用在 8通道腕线圈中测量得到的 M I图像，成像序列 3D MP-RAGE, 参数设置如下： TE = 3.7 ms; TR = 1500 ms; 反转角 = 9°； TI = 900 ms _; FOV = 267mm; 体素 = 0.5 χ 0.5 χ 0.5 mm ³ ; 矩阵 = 512 x 96。

图 2是 Gd-DTPA在 0〜100mM浓度范围内， Gd-DTPA的浓度 C与 MRI信号强度增量 ASI的对应关系。 结果表明： 在 0〜100mM范围内， MRI信号强度增量 ASI随 Gd-DTPA浓度 C的增加呈先上升后下降的钟 形曲线， 在 4mM时达到顶点， 随后逐渐下降； 至 50mM时， 图像出现伪 影； 至 lOOmM时， ASI变为负值， 表明 lOOmM Gd-DTPA的 MRI信号强 度己低于纯水的信号强度。

图 3显示了 Gd-DTPA在 0〜4mM浓度范围内， Gd-DTPA的浓度 C与 M I信号强度增量 ASI的对应关系， Gd-DTPA的浓度 C在 0〜l mM范围 内 MRI信号强度增量 ASI随浓度 C上升较快， 且基本呈线性， 在 4mM 时达到顶点。 由此可以确定出 MRI信号强度增量 ASI的具有线性对应关 系时， 示踪剂浓度的上限值为 4mM。

图 4显示了用腕线圈成像时， Gd-DTPA相关的信号强度增量 ASI与 Gd-DTPA的浓度 C处于 0〜1.5mM之间的线性关系， 其中：

拟合的直线方程为： y=499.54x+194.76,

拟合优度为： R ² =0.9538。

图 5显示了用头线圈成像时， Gd-DTPA相关的信号强度增量 ASI与 Gd-DTPA的浓度 C处于 0~1.5mM之间的线性关系， 其中：

拟合的直线方程为： y=294.26x+197.9，

拟合优度为： R ² =0.8379。 由此可以看出， 应用不同体部线圈成像时， Gd-DTPA相关的信号强 度增量 ASI与浓度之间的线性关系存在差异， 成像时可能根据不同的活 体器官， 选择不同的体部线圈， 以保证拟合出较好的直线。

在本实施例中， 结合脑立体定位图谱， 依三维重建图像确定在大鼠 大脑微量注射 Gd-DTPA的注入点源 Q的位置与进针深度。 具体方法时： 备皮， 由大鼠两眼之间向后沿颅正中线剪幵约 4cm皮肤， 露出颅骨板， 然后将大鼠固定在立体定位仪上， 利用磁场兼容定位给药系统在注入点 源 Q位置处注药， 在进针点用电钻将颅骨板磨去， 将微量注射器降入颅 内， 在目标位置分三次将浓度为 10mM、 总量为 2 μ L的 Gd-DTPA注入， 进药时间间隔 5min。 进药完成留针 lOmin, 退出微量注射器， 缝合皮肤。 随后在 15min、 30min、 lh、 2h、 3h、 4h、 6h、 9h禾口 12h后进行扫描成像， 成像参数同前。

图' 6示意性地显示了处理以上得到的 M I时间序列图像的方法。 如 图所示， 首先进行非均匀场校正， 以消除序列图像之间的信号强度偏差； 然后利用图像配准技术实现增强后各组图像与 参考图像空间位置上的统 -; 最终通过信号强度增量分析获取不同时间点下 的示踪剂浓度分布， 进行扩散模型的求解。

由于 MR序列图像的采集时间比较长， 采集过程中被测对象体位不 可避免的会发生变化。 因此在图像对减之前， 需要利用图像配准技术将 时间序列对齐到同一坐标系下。 此外， 由于 MRI过程的不完善性， MRI 图像会因为强度不均匀场的影响呈现灰度不均 匀性， 即在同一均匀组织 的成像区域里， 图像灰度也会有平滑的变化。 如果不进行不均匀场的校 正， 即使两组图像经过配准后位置完全对齐， 但由于成像时所处不均匀 场的位置不同， 相应位置体素的灰度改变是不同的， 从而对减后得到的 灰度差并不能准确的反映示踪剂的浓度。 所以， 对 MR时间序列特点的 图像进行配准及灰度不均匀校正， 可以更好保证更好的精度验证， 以准 确求解扩散模型。

MR时间序列图像包含了示踪剂的扩散过程，各� ��时间点图像间的灰 度分布具有较大的差异， 给相似性测度带来困难。 此外， 图像间可能同. 时存在刚体变换以及弹性形变两种变换关系， 给 MR时间序列图像的配 准带来一定的影响。 比如图 7所示， 在成像过程中， 大鼠的运动引起大 部分组织明显的形变， 然而脑组织却只存在着刚体变换。

在本发明的一种示意性实施方式中， 采用结合感兴趣区域分割的图 像配准方法， 并选用互信息量做为相似性准则。 图 8示意性地显示了这 种算法的框架。

互信息量是一种基于信息论的图像相似性准则 ， 被广泛的应用多模 图像以及单模图像配准问题。 对于待配准的两幅图像 和 互信息量定 义为两幅图像边缘熵与联合熵之差：

I(f,g) = MiJ) + H(g)― E{f,g)

=: ) ( 1 )， τ ^ ' - p pijy

其中/ ， " ω和 w(/， _g )分别为两幅图像的边缘熵以及联合熵； w和 M 分别为两幅图像中灰度值样本的个数； p ^(/ )， 和 表示灰度值出现 的概率。 互信息是一种区域灰度统计测量， 并不要求两幅图像点与点之 间的灰度对应性， 因此对图像噪声、 灰度变化具有较高的鲁棒性。 尤其 是对于 MR时间序列图像， 由于存在着灰度不均匀性以及由示踪剂引起的 灰度变化， 互信息相比于其它基于图像灰度的相似性测度 更具有优势。

此外， 如图 8所示， 在图像配准的框架中采用了图像分割模块， 将互 信息相似性测度的计算限制在感兴趣的区域（ 比如大鼠的脑）， 从而在配

12

26 准过程中只需考虑刚体变换， 避免了周围组织的弹性形变给配准精度带 来的影响。 进一步地， 为实现快速、 准确的感兴趣区域提取， 采用基于 图论的交互式图像分割方法， 以有效处理因噪声干扰或缺少灰度对比度 造成起目标部分边缘丢失的情况， 并且而支持基于初始分割结果的快速 修改， 因此比较适合 MR时间序列图像复杂背景下感兴趣区域的提取� �� 以 提高图像配准的成功率。

针对 MR时间序列图像的信号强度变化， 在本发明一种示意性的方法 中， 图像间的修正方法为： 首先针对每一组图像进行非均匀场修正， 然 后再实现序列间信号强度一致性校正。根据最 常用的核磁图像组成模型， 成像得到的 MR信号可以表示为：

/( y ∑) ^ J _fR (x _: ...z)B(x._ y, z) + z) ( 2 ) 其中： 为成像信号；

是未受噪声和非均匀场污染的真实核磁信号； S为灰度非均匀场；

"为成像过程中引入的加性噪声。

这样， 通过估计出灰度非均匀场 S， 可以修正非均匀场修正问题， 从 而利用公式 （2) 恢复出真实的 MR信号。

针对此问题， 我们拟采用基于高频最大化的 （High-Frequency Maximization) 的校正方法。 从图像本身的信号强度分布入手， 以最大化 图像灰度分布的高频内容为目标， 通过迭代的过程获取非均匀场的分布。 这种方法不需要先验知识， 而且能够做到完全自动化， 因此比较适合大 数据量的 MR序列图像的处理。

非均匀场校正能够提升单组图像内信号强度分 布的一致性， 但由于 非均匀场的估计只考虑了单组图像的强度信息 ， 校正后可能会降低序列 间相同组织灰度分布的一致性。 尤其是存在信号强度漂移时， 问题可能 更明显。 因此在本发明一种示意性方法中， 还包括一种基于灰度统计特 征实现信号强度的一致性校正的方法。

给定待校正的两幅图像 和 ^，分别从均匀组织的成像区域内选取两 个位置相同的感兴趣区域 和 。 统计感兴趣区域内信号强度的均值和 方差， 得到 ("^)和 从而， 利用下面的线性变换实现两幅图像的 灰度一致性校正： g' =¾ +(.g -«ι) ^χ ·—

(3 )， 其中， g和 分别为原始和变换后的图像信号强度。

感兴趣区域即可以采用人工交互的方式确定， 也可以使用基于 MeanShift算法的区域聚类以及相应的区域匹配方 法自动提取。

如图 8所示， 将标准参考图像与待配准图像的感兴趣区域进 行图像分 割， 然后先根据公式 （1 )应用互信息量进行图像配准， 还可以再应用公 式 （3 )采用基于高频最大化的方法进行灰度校正， 得到配准后图像。 采 用上述方法可以根据 MRI的信号增量计算出的活体器官某处的示踪剂 浓 度， 根据得到的示踪剂浓度， 本发明的方法可以建立示踪剂在活体器官 的扩散模型并求解出反映活体器官特性的参数 。

为了反映 ECS的空间结构特征， 假设空间某一点 P到注入点源 Q空间 的迂回度 /仅与方向有关：

= λ(θ.

a (4 )， 其中， Θ和 φ分别为 PQ与 x、 z轴的夹角。

因此， 沿着方向 PQ扩散的方程可以表示为：

14

26 ( 5 )， 其中： D为示踪剂在琼脂糖中的自由扩散系数，

r表示计算点 P到所述示踪剂注入点源 Q的距离，

C为在局部脑组织中示踪剂的浓度， 由 MRI信号增量 （ASI) 计算得到，

t为注射示踪剂与 MRI图像采集之间的时间间隔， α为脑细胞外间隙占全脑容积的比例，

k'为示踪剂在脑内的消除速率常数，

λ _Ρ(3 为点 P到所述示踪剂注入点源 Q之间脑区的迂曲度。

可以看到， 与现有的扩散模型不同， 由方程（5 ) 建立的扩散模型在 各个方向上是可变的， 这样可以形成不规则的扩散面， 实现这一效果就 需要计算各个方向的迂回度。

由于示踪剂浓度跟它在 MRI图像中的信号增量可以近似为线性关系， 所以， MRI信号增量 ASr满足如下扩散方程：

(6 )， 其中， d表示单个像素对应的实际空间距离,

由此可以得到方程 （5 ) 的解为：

15

26 通过采集不同时刻直线 PQ上的所有像素点的值， 就可以拟合出 PQ方 向的扩散系数 ^{3 =} ^^和损失比率 以及局部体积比"。

尽管方程（5) 实现了不同方向上迂曲度是可以不同的， 但是同一方 向上的迂曲度是一样的， 而实际的迂曲度在同一方向是很可能不同的。 为了更准确地描述真实扩散过程， 需要确定迂曲度的空间可变性。 空间 中某一点 (x,y,z)的迂曲度为：

A = A( x , y , z ) (g) _; 或者写成极坐标形式：

A = Λ( Τ ₃ θ, φ ) ( 9 )， 此时扩散模型变为：

3C D 2 _C ^ Q__

di X( Xj ₃ 2 )^ (10)。

虽然方程 （10) 完全考虑了迂曲度的空间变化性， 但是要得到方程 (8) 的解通常是很困难的， 方程 （8) 是一个变系数的偏微分方程， 它 一般没有解析解。 为了求解方程（10)， 可以把该连续问题转化为如下的 离散问题：

Δ7(/ + Δ/)-Δ67( = ~~ - ~~ _r V ² ( /) + ^i -k'6S{t、 + k'b

Λ(Λ-, a (11)， 其中， 表示单位采样时间，

V ² (ASI)= \ ASJ(i + Δ) + &SJ {j - Ai) - 2 AS! (i)

A表示,轴方向上的空间采样距离，它取决于� ��像在,轴上的空间分辨 率。

对于离散方程 （11)， 可以采用最小二乘法拟合， 从而求解出空间迂 曲度 和损失比率 以及局部体积比。， 即需要求解如下优化问题：

16 26 D d ² ( ASl) Q

+ -/-/r' (厶 SI ) + /、-' Z, - Δ57 (ί + Δή + AS1 ( δν'

(12)。 将方程 （12) 写成矩阵形式， 得到:

其中， 矩阵 Α的第 (x,y, _Z ,t, _: )行的非零元为：

b (Λ-, y,∑,!：,：) = AS1 (/ + Δ/) - ASf (,)

(14)， 已知向量 b的第 (x,y，z,t,:)行为： b(x,y,z,†,^ ASI(( + At)― ASi(i)

(15), 替换变量 χ为： (]6)。

当时间采样率和空间采样率非常高时， 离散方程 （11) 可以等价于 连续方程 （10)。 通常情况下， 实际当中的采样率偏低， 这时得到的空间 迂曲度'和损失比率 以及局部体积比。不是严格意义上的概念， 而是"平 均"意义下的概念。

分子在 ECS中某点扩散， 可能具有各向异性的特点， 即沿着某个方向 扩散特别快， 沿着另一个方向被阻碍。 因此， 我们需要考虑分子各向异 性扩散的情形。 假设各向异性扩散张量的主方向跟坐标轴一致 ， 可以把 (2) 推广到各向异性的情况：

dc D & ² C D & ² C D d ² C Ο ' -

—k

(17),

SykovaE和 Nicholson C给出了方程 （3) 的解: (18)， ^中 ?■= - + ''. 在各向同性的介质中， A=A _S =' = ， 方程 （ ₁₈₎ 退化为 （ ₁₇₎ 。 可以看到方程（17) 中的迂曲度为常数， 无法反映出实际的 ECS的结 构。 为了反映 ECS的空间结构， 我们可以将离散方程（11)推广到各向异 性的情况，

c(t + - c(0 = ' , ^D 、 v' - li f)

' · (19)。 同样，用最小二乘法类似地求解出三个轴向的 空间迂曲度 和损失比 率 以及局部体积比。。 不同于各向同性情况， 此时的系数矩阵 Α和替换 变量 X变为：

(20),

D D D O 与各向同性的离散方程相比， 方程 （19) 的待定参数多了几乎 2倍， 因此， 利用方程 （19) 可以得到更高的拟合精度， 即目标函数方程 （12) 会更小， 从而能更好的预测实际的扩散浓度。

在求解脑组织间液流动性参数过程中， 需将大量的示踪剂浓度值与 其相应的位置参数与时间参数导入方程中进行 拟合。 因此， 采用本发明 的计算方法， 利用相应的计算机程序自动读取浓度分布图中 各像素点浓 度值与其相应位置参数与时间参数， 自动拟合给出组织间液流动性参数 值。

图 11显示了将 SD大鼠麻醉后放入磁共振成像仪腕线圈中， 用微量进 样器将 2 μΐ,浓度 10 mM的 Gd-DTPA溶液由不同的注入点源注射至脑内不 同脑区， 在不同时间扫描后采集的图像。 其中 （a) 的注入点源位于尾状 核区， （b ) 的注入点源位于黑质紋状体， （C ) 位于丘脑， （d ) 位于枕叶 皮层， （e) 位于皮质下白质。

采用本发明的模拟计算方法， 对 SD大鼠尾状核、 黑质纹状体、丘脑、 枕叶皮层以及大脑皮质下白质区域采集到的图 像进行了计算， 得到了如 表 1所示的结果， 可以计算出脑内不同区域的扩散参数； 表 2显示为大鼠 尾状核区扩散参数的各向异性；表 3为示踪剂 Gd-DTPA经尾状核注射后在 脑内的分布容积随时间变化的规律。

表 1 : SD大鼠不同脑区的扩散参数

表 2: 大鼠尾状核区域扩散参数的各向异性

19

26 表 3 : Gd-DTPA经尾状核注射后在脑内的分布容积随时� �变化的规律

可示踪物质在活体器官中的分布容积 Vd、 分布容积分数和大体清除 率的计算方法如下： 将配准和对减后的图像重建为轴位图像， 采用如图 所示的 Matlab程序定量测量各层图像中信号增强范围， 记录增强范围中 包含的体素数量和每个体素的信号强度值。 根据积分原理， 将各层增强 区域面积进行累加， 即得到各时刻示踪剂 Gd-DTPA在脑 ECS 内分布容 积 Vd ， 并将各时刻所得的 Vd除以大鼠大脑半球总体积获得 Gd-DTPA 在脑半球 ECS中的容积分数。 将每个体素的信号强度值相加得到增强范 围内的总信号强度。 宏观清除率的定义是前一个时间点的总信号强 度减 去后一个时间点的总信号强度， 再除以初始的总信号强度。

'图 9显示了使用 Matlab软件测量对减后的未增强和增强后的区域 。 图中圈中的区域是用增强后的图像减去增强前 的图像所得到的， 因此增 强区域显示为 "亮区"， 而其他区域为无增强区域， 因此显示为"暗区"。 图 10是图 9所示中增强区域 MRI信号测量后生成的测量结果图， 横坐 标是信号强度值， 纵坐标是像素个数。 应用上述方法测量出的尾状核区

20

替 （ 第 26条) 域的 Vd、 分布容积分数和宏观清除率的结果。 表 4所示为向八只 SD大 鼠尾状核区域注射 Gd-DTPA后， 在不同的时间内测量出的宏观清除率。

表 4: 八只 SD大鼠尾状核区域注射 Gd-DTPA后的宏观清除率

虽然上述的实施例是针对 SD大鼠的脑 ECS， 但本发明同样可用于其 它活体器官。 图 12是 Gd-DTPA注入大鼠骨骼肌间隙内在不同时间扫描� � 采集的扩散图像， 图中所示左侧股骨旁高信号区域为 Gd-DTPA注射后在 骨骼肌间隙内的扩散区域。

表 5显示了 Gd-DTPA经骨骼肌间隙内注入后的扩散参数。表 6和表 7分 别显示了 Gd-DTPA在骨骼肌间隙内的分布容积随时间变化� �规律和宏观 清除率。

表 5 : Gd-DTPA在骨骼肌间隙内的扩散参数

表 6: Gd-DTPA在大鼠骨骼肌间隙内的分布容积随时间� �化的规律

表 7: Gd-DTPA在大鼠骨骼肌间隙内的宏观清除率

图 13显示了示踪剂 Gd-DTPA在胶质瘤模型大鼠脑内的扩散图像， 其 中 （a) 的 T1WI脑内高信号区域显示为 Gd-DTPA的扩散区域， 低信号区 域为瘤体部位； （b) 的 T2WI高信号区域显示为瘤体部位。

表 8显示了 Gd-DTPA在胶质瘤大鼠脑内的扩散参数。

26 表 8: 胶质瘤大鼠脑内的扩散参数

采用与前面用于脑 ECS的方法相同方法， 可以了解 Gd-DTPA在胶 质瘤模型大鼠脑内的分布容积随时间变化的规 律。结果参见表 9和表 10。 表 9: Gd-DTPA在胶质瘤模型大鼠脑内的分布容积随时� �变化的规律

表 10: Gd-DTPA在胶质瘤模型大鼠脑内的宏观清除率

本发明提供的计算方法， 可以经过对感兴趣的活体器官的 MRI扫描 图像， 准确计算出可示踪物质在活体器官组织间隙分 布与细胞外间隙扩 散的参数， 有助于疾病的诊断及药代动力学等的研究。

上文所列出的一系列的详细说明仅仅是针对本 发明的可行性实施例 的具体说明， 它们并非用以限制本发明的保护范围， 凡未脱离本发明技 艺精神所作的等效实施例或变更均应包含在本 发明的保护范围之内。

24

26