脑是中枢神经系统的主要部分，主要由大脑、小脑、脑干等构成[1]．大脑的结构错综复杂，有数千亿个神经元．脑中分布着很多由神经细胞集中而成的神经中枢，并有大量上、下行的神经纤维束通过，连接大脑、小脑和脊髓．小鼠脑模型的研究有助于理解临床试验的过程，对人脑的研究也有参考意义．临床上有许多用于大脑的成像技术，常用的X射线断层成像[2]需要特定的造影剂来显示血管，且有电离辐射．超声成像[3]对人体无辐射，运动灵活，因而被广泛应用，但其分辨率较差、对比度较低．正电子发射层析成像[4]作为灵敏度非常高的影像技术，也面临着成像分辨率低、使用成本高昂的问题．磁共振成像[5]使用成本较高，并且三维空间分辨率不高(约为1 mm)．与上述临床成像技术相比，光学成像技术具有高分辨率的特点．常用的光学成像方法有光学相干层析成像[6]、活体显微镜[7]、激光共聚焦扫描显微镜[8]和双光子显微镜[9]等．然而，由于颅骨对光有强烈的散射，光子很难穿透颅骨到达深处，因此这些光学成像模式在脑部成像上的使用也受到限制，通常采用磨薄颅骨、光透明和颅骨揭除等方法进行成像．为了探讨颅骨对光子传输的影响，本研究采用蒙特卡罗方法模拟了光子在小鼠大脑中的传播．分别研究了光子在未开颅和开颅下的小鼠脑组织中的传播，并分析了相应的光吸收分布、光能流率和光子轨迹．1 模型与方法1.1　脑组织特性脑组织是一种多层介质，每一层的结构和功能都不相同，光学特性也不相同．小鼠的大脑在结构上与人类的大脑相似．当光从外部照射到大脑时，光子须要通过头皮、颅骨、脑脊液、脑部(包含灰质和白质)，如图1所示．使用一个水平方向无穷大、垂直方向有限深的半无限模型来模拟光子在脑组织中的传输．表1为各层组织的光学参数[10]，表中：d为厚度；μa为吸收系数；μs为散射系数；g为各向异性系数；n为折射率．去除头皮和颅骨后，用生理盐水覆盖头部以维持大脑微环境．开颅后大脑各组织的光学参数如表2所示．10.13245/j.hust.238792.F001图1鼠脑多层介质模型10.13245/j.hust.238792.T001表1未开颅下四层脑组织的光学参数及其厚度脑组织d/cmμa/cm-1μs/cm-1gn头皮颅骨脑脊液脑部0.100.100.051.000.1910.1360.0260.186668681110.90.90.90.91.371.371.371.3710.13245/j.hust.238792.T002表2开颅下四层脑组织的光学参数及其厚度脑组织d/cmμa/cm-1μs/cm-1gn生理盐水生理盐水脑脊液脑部0.100.100.051.000.0260.0260.0260.1868881110.90.90.90.91.371.371.371.371.2　蒙特卡罗模拟蒙特卡罗方法是一种统计模拟方法[11]，也被称为随机抽样或统计检验方法．蒙特卡罗方法已被广泛应用于模拟光子在生物组织中的传输，其基本思路如下：运用光子吸收和散射的原理来跟踪光子通过混沌介质的过程；通过统计大量光子的模拟结果得到能量在组织内部的分布情况，其可行性被大量实验论证[12-15]．Wlison等运用蒙特卡罗模型预测不同吸收/散射比的均匀组织中光吸收分布，实验中使用植入组织中的光纤来收集光子，并对预测结果进行了实验论证[12]．Keijzer等运用蒙特卡罗方法模拟氩激光照射期间动脉组织中的有限宽度光分布[13]．Kanick等采用蒙特卡罗模型研究了光学特性对单光纤反射率测量中收集到的光子的传播影响[14]，蒙特卡罗模型估计的单光纤光子路径长度与实验测量的长度符合．Maeda等提出了一种用于蒙特卡罗模拟的九层皮肤组织模型，利用适当的光学参数和几何参数对传统三层模型中所使用的参数进行修正，得到与实测光谱符合较好的模拟光谱[15]．Wang等发表了蒙特卡罗程序MCML[16]，MCML程序实际上成为模拟光子在组织中传播的经典模型[17-18]，研究者对MCML程序进行了验证和改进．Boas等改进MCML中光子通量的记录方式，提出了基于体素化模型处理的3D蒙特卡罗模型[17]．Yun等基于MCML改进程序，模拟偏振光子在各项异性散射模型中的传播，偏振成像模拟和实验结果之间良好的一致性证实了蒙特卡罗模拟的有效性[18]．蒙特卡罗的模拟过程如图2所示．给定一个光子，设初始权重w为1，初始位置为(0，0，0)，初始方向余弦为(0，0，1)．设外部介质和生物组织表层的折射率分别为n0和n1，则镜面反射的反射系数Rsp为(n0-n1)2/(n0+n1)2，进入生物组织之后的光子权重减小为w=1-Rsp．光子开始在组织中传输，此时须要计算光子的步长s=(-lnξ)/μt，其中ξ为(0，1)之间均匀分布的随机数；μt=μa+μs，μa为吸收系数，μs为散射系数．设入射光子的当前位置为(x,y,z)，光子传播方向由其方向余弦(μx,μy,μz)确定，则下一位置坐标(x',y',z')为：10.13245/j.hust.238792.F002图2蒙特卡罗模拟流程图x'=x+μxs;y'=y+μys;z'=z+μzs. (1)设光子初始权重为w0，光子与组织之间的相互作用导致光子被吸收，权重减小．第n次相互作用后光子的权重wn=w0(μs/μt)n．光子在传播过程中发生散射，此时须要计算光子的散射角θ和方位角ψ．当各项异性因子g≠0时，光子的散射角为θ=arccos12g1+g2-1-g21-g+2gξ2．(2)由于散射有轴对称性，因此方位角ψ=2πξ．光子新的方向余弦(μx',μy',μz')可表示为：μx'=sin θ1-μz2(μxμzcos ψ-μysin ψ)+μxcos θ;μy'=sin θ1-μz2(μyμzcos ψ+μxsin ψ)+μycos θ;μz'=-sin θcos ψ1-μz2+μzcos θ.(3)当光子到达生物组织边界时，定义由光子的当前位置(x,y,z)沿其传播方向到生物组织边界的距离为约化步长s1，即s1= (z0-z)/μz(μz≤0);(z1-z)/μz(μz0), (4)式中z0，z1分别为组织上下边界在直角坐标系中的z坐标．若ss1，则光子未到达组织边界，将在组织中继续传播．若ss1，则光子到达组织边界，设此时光子的入射角为αi，透射角为αt，αi=arccos(|μz|)．到达边界的光子发生全反射的反射率R(αi)=12sin2(αi-αt)sin2(αi+αt)+tan2(αi-αt)tan2(αi-αt)．(5)然后判断光子是否发生全反射．当ξR(αi)时光子发生透射，此时根据每层的组织厚度判断光子是到达上下表面还是到达组织交界：若到达组织交界，则进入下层组织，进入下一步模拟；若到达上下表面，则逃逸出组织．当ξ≤R(αi)时，光子发生全内反射，光子方向余弦更新为(μx,μy，-μz)，此时剩余步长为s-s1．若剩余步长足够大，光子仍能到达其他边界，重复以上过程；若剩余步长不够大，光子继续留在组织中传输，发生吸收和散射．当光子的权重小于设定的阈值wth=0.000 1时结束跟踪，继续跟踪下一个光子，直至最后一个光子跟踪结束．本研究采用Wang等编写的蒙特卡罗模拟程序MCML[16]．对环境进行配置，包括配置光子数、二维栅格系统的格间距dr和dz、栅格数目、组织层数及各层组织的光学参数等．光子数设置为1×104个，二维栅格系统的格间距dr和dz分别设置为0.1 cm和0.1 cm，栅格数目分别设置为400和200，组织层数设置为4．将开颅和未开颅下小鼠大脑的光学参数代入，最终得到光能流率、光吸收能量密度等参量．2 结果与讨论2.1　光子在脑中传播的轨迹生物组织是一种高散射的随机介质，光子在不同特性组织中传播有很大的差异．图3(a)和图3(b)分别模拟开颅下和未开颅下脑组织中的光子轨迹，从图中可以看出开颅下光子在z方向的运动比未开颅下更深．由于头皮和颅骨对光子有较强的散射和吸收，而脑脊液层是一种低吸收和低散射的介质，因此光子可以在开颅下移动到更深的组织．在图3(a)中，光子运动的径向范围大于图3(b)，开颅下光子主要集中在径向-10~10 cm之间，未开颅下光子主要集中在径向-5~5 cm之间．开颅下获得的光子分布比未开颅下获得的光子分布更“胖”．10.13245/j.hust.238792.F003图3光子在鼠脑中的运动轨迹图2.2　光子在脑中传播的光吸收分布和光能流率分布图4显示了开颅和未开颅下光吸收能量密度A 和光能流率F的分布．从图4中可以看出：开颅后，光子传输显得较深，径向(即r方向)显得较宽，如图4(a)和图4(b)所示，光子可以不受阻碍地到达组织深处；脑脊液层和脑部之间有一条明确的分界线(图中白色箭头所示)，这是因为两层的吸收系数和散射系数相差较大；未开颅下，由于头皮和颅骨的强烈吸收和散射作用，光能量迅速衰减，只有少数光子能到达灰质层，如图4(b)所示；光能流率F也表现出和光吸收能量密度同样的变化规律，如图4(c)和4(d)所示．10.13245/j.hust.238792.F004图4光吸收能量密度、光能流率分布图2.3　光子在脑中传播的吸收曲线图5展示了光子在轴向和径向的传播情况．图5(a)展现的是光吸收能量密度A在深度方向的变化．当未开颅时，由于头皮和颅骨的强烈吸收和散10.13245/j.hust.238792.F005图5z和r方向上开颅和未开颅下的光吸收能量密度、光能流率射，lg {A}J/cm2从0.948 3(A=8.836 J/cm2)迅速下降到0.030 36(A=1.072 J/cm2)，剩余少量光子进入灰质．由于灰质的强吸收和散射，因此lg {A}J/cm2略有上升，从0.030 36(A=1.072 J/cm2)上升到0.054 31(A=1.133 J/cm2)，然后不断下降．当深度为1.3 cm 时，lg {A}J/cm2的值下降到-1.974(A=0.010 62 J/cm2)．开颅后，由于生理盐水和脑脊液的吸收和散射较弱，因此lg {A}J/cm2较低，为-0.107 7(A=0.780 3 J/cm2)，并缓慢上升，大多数光子可到达灰质．由于灰质的强吸收和散射，因此lg {A}J/cm2从-0.107 7(A=0.780 3 J/cm2)上升至0.832 7(A=6.804 J/cm2)，然后不断下降；当在1.3 cm深度时，lg {A}J/cm2的值下降到-1.797(A=0.015 97 J/cm2)．lg {F}J/cm2的变化与lg {A}J/cm2的变化相似，如图5(b)所示．在径向方向上，光子的传播范围定义为lg {A}J/cm2下降至其最大值一半的径向范围．如图5(c)所示，未开颅下，当光子扩散到径向距离约为-1.3 cm和1.3 cm时，lg {A}J/cm2下降到一半，即扩散宽度为2.6 cm．开颅下，当r为-2.6 cm和2.6 cm时，lg {A}J/cm2下降到一半，即扩散宽度为5.2 cm．可以看到在相同的条件下，开颅后光子在r方向的传播范围更广，是未开颅下的2倍．lg {F}J/cm2的变化与lg {A}J/cm2的变化类似，如图5(d)所示．3 结语本研究用蒙特卡罗方法分析了颅骨对光子传输的影响，并对光子轨迹、光吸收能量密度和光能流率进行了分析．开颅下光子轨迹运动范围更广，约为未开颅下的2倍．光吸收能量密度和光能流率的分布总体上是梭形的，其径向分布范围约为未开颅下的2倍．在深度方向，开颅后由于生理盐水和脑脊液的吸收和散射相对于颅骨较弱，到达脑部(灰质和白质)的剩余光子数是未开颅时的约6.35倍．