孫愛扣,陳珍平,3?,匡藍(lán)珺,高克坤,劉程偉,于 濤

(1. 南華大學(xué) 核科學(xué)技術(shù)學(xué)院;2. 南華大學(xué) 先進(jìn)核能技術(shù)設(shè)計與安全教育部重點實驗室：湖南衡陽421001;3. 蘇州大學(xué) 放射醫(yī)學(xué)與輻射防護(hù)國家重點實驗室, 江蘇蘇州 215123)

惡性腫瘤現(xiàn)已成為導(dǎo)致人類死亡的頭號殺手,且發(fā)病率正在逐年增加。2022年,國家癌癥中心發(fā)布我國癌癥患者為406.4萬人[1]。目前,癌癥的治療方式主要有手術(shù)治療、化學(xué)治療及放射治療。國內(nèi)將近50%～70%腫瘤患者在治療過程中接受放射治療。其中,硼中子俘獲治療(boron neutron capture therapy,BNCT)是一種嶄新的二元靶向放射治療手段,主要原理是通過將具有選擇性的含硼藥物注射入人體內(nèi),讓10B富集于腫瘤區(qū),用熱中子束照射10B富集區(qū)域,由于10B具有異常大的熱中子吸收截面,大量的10B會發(fā)生10B(n,α)7Li反應(yīng),釋放出α,7Li 2個具有高傳能線密度(linear energy transfer,LET)的粒子,這兩種粒子在約12～13 μm的范圍內(nèi)沉積能量,使沉積能量區(qū)域的細(xì)胞DNA產(chǎn)生難以修復(fù)的損傷,最終殺死腫瘤細(xì)胞,治療效果優(yōu)于常規(guī)治療手段[2]。

BNCT的治療包括計算和分析患者體內(nèi)的照射劑量分布,用以確定中子射束照射時間及方向角度的最優(yōu)化配置,并同時遵守危及正常組織和器官的劑量限制。BNCT中組織劑量基本來自于：(1)硼劑量。由10B(n,α)7Li反應(yīng)所獲得,且由于硼對熱中子的吸收截面非常大(3 840 b),釋放能量多,因而硼劑量對總劑量的貢獻(xiàn)很大,為主要劑量;(2)熱中子劑量。熱中子與人體中的14N原子發(fā)生核反應(yīng)生成反沖核14C和質(zhì)子,在人體細(xì)胞中沉積能量所產(chǎn)生的的劑量;(3)超熱中子及快中子劑量。統(tǒng)計的是除熱中子外所有中子由彈性散射而沉積的能量;(4)光子劑量。主要有伴隨入射中子束產(chǎn)生的γ射線、硼中子俘獲反應(yīng)中產(chǎn)生的射線及熱中子與人體內(nèi)的H原子發(fā)生俘獲反應(yīng)生成的γ射線3種來源。上述4種劑量單元的計算非常復(fù)雜,目前尚無精確劑量計算的相關(guān)經(jīng)驗公式,現(xiàn)階段BNCT的治療計劃系統(tǒng)中主要采用蒙特卡羅模擬計算患者或體模的個體化模型劑量[3]。目前,國際上針對BNCT的治療計劃系統(tǒng)(treatment planning system,TPS)日漸成熟,開發(fā)了多種治療軟件,主要有：美國的NCTPlan和SERA、日本的JCDS、北京應(yīng)用物理與計算數(shù)學(xué)研究所研制的MCDB及中硼醫(yī)療公司研發(fā)的NeuMANTA等[4]。

本文基于成熟的蒙特卡羅方法研制了自主化BNCT蒙特卡羅劑量計算程序Magic。針對國際基準(zhǔn)題Snyder修正頭部模型,在此基礎(chǔ)上添加一個假設(shè)的腫瘤模型,建立1 mm精細(xì)體素模型來逼近解析模型,并在程序中采用網(wǎng)格中心點方法填充體素模型材料,同時對中子源、Kerma因子及S(α,β)熱散射等參數(shù)設(shè)置,使用MCNP和Magic程序進(jìn)行模型深度-劑量率曲線模擬對比,對結(jié)果進(jìn)行分析研究,初步驗證了Magic程序正確性,并以Magic為基礎(chǔ),研究10B藥代動力學(xué)模型下的硼濃度對BNCT劑量學(xué)影響規(guī)律,可為BNCT 劑量計算基準(zhǔn)化和臨床實踐應(yīng)用提供參考。

1 蒙特卡羅劑量計算程序Magic開發(fā)

目前,國際上BNCT治療計劃系統(tǒng)的劑量計算引擎主要使用MCNP作為核心計算傳輸工具,然而MCNP的禁止轉(zhuǎn)讓限制了其在中國的使用。劑量計算是整個BNCT治療計劃系統(tǒng)的核心部分,MCNP的限制將影響國內(nèi)自主化BNCT治療計劃系統(tǒng)的發(fā)展,同時MCNP程序體系龐大,主要用于核能、核工程及核技術(shù)等各領(lǐng)域的相關(guān)理論計算研究,并不適用于特殊研究需求的發(fā)展。鑒于此,本文基于蒙特卡羅粒子輸運方法自主開發(fā)了BNCT治療計劃系統(tǒng)專用計算引擎—蒙特卡羅劑量計算程序Magic。Magic主要采用C++開發(fā),借助于C++豐富的標(biāo)準(zhǔn)庫和第三方庫,覆蓋了多個可擴(kuò)展功能,便于未來的BNCT治療計劃系統(tǒng)的發(fā)展和維護(hù)。

Magic程序由7塊功能模塊組合構(gòu)成,包括幾何模塊(Geometry Module)、源項模塊(Source Module)、數(shù)據(jù)庫模塊(Database Module)、粒子輸運模塊(Particle Transport Module)、計數(shù)模塊(Tallies Module)、輸出模塊(Output Module)和輔助模塊(Auxiliary Module),圖 1為Magic開發(fā)框架。幾何模塊主要采用面方程描述方法構(gòu)造復(fù)雜的劑量計算模型,基于組合幾何(CSG)方法進(jìn)行幾何描述,同時支持基于Lattice重復(fù)結(jié)構(gòu)的超精細(xì)毫米量級體素(voxel)建模,并對相應(yīng)組織區(qū)域從數(shù)據(jù)庫模塊中填充組織材料的核素成分及比例;源項模塊旨在設(shè)置粒子的初始條件和源的屬性,如源的空間分布、能量分布、角度分布及粒子類型等信息,并對相應(yīng)源參數(shù)支持多種概率分布抽樣(離散、連續(xù)、直方圖及混合),滿足模擬BNCT復(fù)雜中子源的要求;數(shù)據(jù)庫模塊儲存了基于ENDF/B-VIII.0評價數(shù)據(jù)庫,利用具有層次結(jié)構(gòu),可處理大數(shù)據(jù)集等優(yōu)點的HDF5數(shù)據(jù)格式來存儲粒子輸運中物理相關(guān)的各反應(yīng)截面數(shù)據(jù)[5],同時根據(jù)ICRU46文和ICRU63文內(nèi)容[6-7],創(chuàng)建了包括106種人體組織成分材料庫和相應(yīng)的Kerma因子庫,并提供了支持修改和編輯相關(guān)數(shù)據(jù)的功能;粒子輸運模塊主要處理程序模擬粒子從產(chǎn)生到消失的過程,中子與物質(zhì)相互作用主要涉及到彈性散射、非彈性散射及吸收,其中,吸收包含BNCT所關(guān)注劑量相關(guān)反應(yīng)(n,α),(n,p),(n,γ)等,在此過程中由中子產(chǎn)生的次級光子將被儲存在粒子庫中,為后續(xù)的光子輸運做好準(zhǔn)備,光子與物質(zhì)相互作用主要考慮光電效應(yīng)、電子對效應(yīng)、相干散射及康普頓散射;粒子輸運完成后,調(diào)用計數(shù)模塊,根據(jù)用戶的需求,設(shè)置不同的統(tǒng)計參數(shù),包括感興趣區(qū)域、能量區(qū)間和粒子類型等,以獲取相應(yīng)的注量,結(jié)合數(shù)據(jù)庫模塊中的Kerma因子庫,得到劑量;最后通過輸出模塊對BNCT治療中關(guān)注的各種劑量,包括硼劑量、熱中子劑量、超熱中子劑量、快中子劑量、光子劑量及總相對生物效應(yīng)劑量進(jìn)行統(tǒng)計,并提供可視化處理;輔助模塊為上述模塊提供相關(guān)函數(shù)庫及異常捕獲等。

Magic程序以粒子輸運模塊為核心,與上游的幾何模塊、源項模塊及截面模塊相互連接,同時支持下游的計數(shù)模塊和輸出模塊。各模塊之間相互獨立,且模塊間通過接口進(jìn)行消息傳遞實現(xiàn)各模塊功能,利于程序維護(hù)和擴(kuò)展。

圖1 Magic開發(fā)框架

2 模型及方法描述

2.1 Snyder修正頭部模型

目前,國際上對于BNCT劑量理論計算普遍采用的基準(zhǔn)問題為Snyder修正頭部模型,由3個橢球面來界定不同組織的邊界,將頭部由里到外分為腦組織、顱骨及頭皮組織3部分。整個頭部又被空氣包圍,通常稱為解析模型[3]。各分界面滿足的3個橢圓面解析方程可表示為

腦組織與顱骨邊界：

(1)

顱骨與頭皮組織邊界：

(2)

頭皮組織與空氣邊界：

(3)

但在實際情況中,由于患者腫瘤形狀、位置及大小各不相同,無法使用解析方程去描述模型,因此在國際臨床中,普遍使用連續(xù)均勻的立方體網(wǎng)格模型來逼近模型結(jié)構(gòu),通常稱為體素模型。目前Snyder修正體素模型網(wǎng)格尺寸通常取4,8,16 mm,當(dāng)構(gòu)建的體素網(wǎng)格越小,越能更好地趨近解析模型。故本文利用Magic構(gòu)建1 mm×1 mm×1 mm精細(xì)Snyder體素模型,共分224(x方向)×224(y方向)×224(z方向)=11 239 424個體素網(wǎng)格,驗證正確性的同時,也對Magic構(gòu)建千萬網(wǎng)格數(shù)的承載力提出考驗。模擬計算的同時,在腦組織中加入半徑為1.5 cm的球面作為腫瘤邊界,模擬實際臨床中腦部膠質(zhì)腫瘤的情形。本文所模擬的含腫瘤Snyder修正模型如圖2所示。

含腫瘤Snyder修正模型中腦組織、顱骨及頭皮組織的元素成分和密度由Magic數(shù)據(jù)庫獲得,腫瘤詳細(xì)的組織成分參考文獻(xiàn)[8]。目前臨床試驗要求,腫瘤細(xì)胞的10B濃度必須大于正常組織細(xì)胞的2.5倍,故本文設(shè)置頭皮組織及腦內(nèi)正常組織硼濃度為10-5,腫瘤組織內(nèi)的硼濃度為3×10-5,顱骨內(nèi)不含硼[9-10]。對解析模型中材料,可精確填充描述,但對網(wǎng)格模型,會因網(wǎng)格在兩種物質(zhì)交界面處難以處理,傳統(tǒng)處理方式通過得到每個立方體網(wǎng)格內(nèi)部不同材料精確的體積比填充,當(dāng)網(wǎng)格數(shù)量較多時,材料填充極其耗時且會帶來一定偏差。為降低模擬問題的復(fù)雜度,在Magic程序模型構(gòu)造中使用網(wǎng)格中心點的方法,根據(jù)每個網(wǎng)格的中心點位置確定材料密度,能大大降低材料數(shù)目,可由原來的286種材料降低至4種材料[11]。

(a)Tumor modified Snyder analytical model

2.2 相對生物學(xué)效應(yīng)劑量

由于不同種類的電離輻射會導(dǎo)致生物體產(chǎn)生不同的生物效應(yīng),BNCT涉及多個不同輻射劑量組分的混合場照射,因此,對BNCT治療須衡量每種劑量組分的相對生物學(xué)效應(yīng)劑量,評估其對正常組織細(xì)胞和腫瘤細(xì)胞的殺傷效應(yīng)。用實驗測量的相對生物學(xué)效應(yīng)(relative biological effectiveness,RBE)來表征熱中子劑量、超熱中子及快中子劑量與光子劑量的生物學(xué)效應(yīng)高低。對于硼劑量的評估,須考慮α粒子和7Li離子的RBE及10B在組織和細(xì)胞中的微觀幾何位置對生物學(xué)效應(yīng)的綜合影響,被稱為復(fù)合生物學(xué)效應(yīng)(compound biological effectiveness,CBE)[12]。將上述4種劑量成分分別乘上各自對應(yīng)的加權(quán)因子得到總相對生物學(xué)效應(yīng)劑量H,表示為

H=WBDB+WγDγ+WnDn+WPDP

(4)

其中：WB為硼劑量的CBE因子;Wγ,Wn及WP分別為光子、快中子及氮俘獲產(chǎn)生熱中子劑量的RBE因子;DB,Dγ,Dn,DP分別為硼劑量、光子劑量、快中子劑量及氮俘獲產(chǎn)生的熱中子劑量。為使本文計算結(jié)果與國際上BNCT劑量結(jié)果具有可比性,選取腫瘤及健康組織不同劑量成分的RBE和CBE因子[13-16],如表1所列。

表1 腫瘤及健康組織不同劑量成分的RBE和CBE因子[13-16]

2.3 中子源設(shè)置

合適的中子源是BNCT治療腫瘤成功的關(guān)鍵要素之一,理想情況下,需選擇不同中子源來治療不同深度處的腫瘤。BNCT劑量計算中,要對入射源項粒子束的類型和形狀、粒子源的種類、粒子束能量及方向等參數(shù)進(jìn)行詳細(xì)描述。本文采取國際上針對Snyder修正模型所提出的寬譜超熱混合中子束,其中能量小于0.5 eV的熱中子占10%,0.5 eV～10 keV的超熱中子占89%,10 keV～2 MeV的快中子占1%,總體中子源能譜服從1/E(E為中子能量)分布,均勻分布在5 cm的圓盤面上,源強度為1010cm-2·s-1。

2.4 其他相關(guān)參數(shù)的設(shè)置

由于BNCT劑量計算時涉及中子-光子耦合輸運過程,劑量可通過注量率與Kerma因子得到,表示為

(5)

其中：φ(j,E)為第j個體素網(wǎng)格,能量為E的歸一化注量率;V(j)為第j個體素網(wǎng)格的體積;K(j,E)為j個體素網(wǎng)格,能量為E的中子/光子劑量轉(zhuǎn)換因子,即Kerma因子。本文采用ICRU63給出的相對于ICRU46成年人腦的Kerma因子,但該Kerma因子所對應(yīng)的最低能量為0.025 3 eV,而對于BNCT的劑量來說,低于0.025 3 eV能量的中子對熱中子劑量的貢獻(xiàn)非常重要,故本文采取雙對數(shù)插值,將低于0.025 3 eV對應(yīng)的Kerma值外推到0.000 1 eV對應(yīng)的數(shù)據(jù),本文中光子的Kerma因子基于Seltzer[17]所計算的質(zhì)能吸收系數(shù)得到。

低能中子對BNCT的劑量貢獻(xiàn)很大,但當(dāng)中子能量降低到幾個電子伏時,散射靶核的熱運動會對碰撞產(chǎn)生強烈的影響,從而對出射中子的能量及出射角度都會造成影響,最終會導(dǎo)致結(jié)果產(chǎn)生偏差。MCNP對于這種情況提供了S(α,β)熱散射模型,直接調(diào)用相應(yīng)核素的S(α,β)熱散射模型數(shù)據(jù)即可[18],本文Magic采取相同的S(α,β)熱散射模型處理Snyder材料中H的熱中子散射。

3 數(shù)值驗證

3.1 MCNP與Magic計算結(jié)果比較

由于MCNP程序通?？勺鳛楦黝惷商乜_計算的“標(biāo)準(zhǔn)”,常常用于標(biāo)定其他蒙特卡羅程序的準(zhǔn)確性測試,尤其是MCNP已成熟運用在BNCT的治療計劃中。本文利用MCNP程序計算模擬Snyder解析模型的各劑量率作為標(biāo)準(zhǔn),使用Magic程序采用上文參數(shù)設(shè)置模擬同一解析模型,比較二者計算結(jié)果的差異,并在此基礎(chǔ)上繼續(xù)利用Magic程序?qū)? mm千萬量級體素網(wǎng)格的Snyder模型進(jìn)行模擬,將模擬結(jié)果與MCNP解析模型結(jié)果進(jìn)行比較。本文重點在于探究Magic程序的正確性,對模型沿z軸統(tǒng)計相關(guān)劑量,其中每個統(tǒng)計點的計數(shù)網(wǎng)格大小設(shè)置為1.6 cm×1.6 cm×0.4 cm。在實際BNCT治療中有部分光子由中子源伴隨產(chǎn)生的,這部分光子是無法避免的,本文假設(shè)中子源無伴隨光子污染,只討論中子與體內(nèi)的核素發(fā)生反應(yīng)所產(chǎn)生的次級光子。

為減小蒙特卡羅程序的深穿透所帶來的結(jié)果影響,兩個蒙特卡羅程序輸運模擬的粒子數(shù)均為1×108,確保樣本數(shù)量充足。圖3-圖7分別為硼劑量率、熱中子劑量率、超熱中子及快中子劑量率、次級光子劑量率和總的相對生物效應(yīng)劑量率及相對偏差隨深度的變化關(guān)系。相對偏差是以MCNP解析模型結(jié)果為標(biāo)準(zhǔn)值,Magic的解析/體素模型的結(jié)果與之比較得到的偏差百分比。

(a)Dose rate vs. depth

(b)Relative deviation vs. depth

(a)Dose rate vs. depth

(b)Relative deviation vs. depth

(a)Dose rate vs. depth

(b)Relative deviation vs. depth

(a)Dose rate vs. depth

(b)Relative deviation vs. depth

(a)Dose rate vs. depth

3.2 結(jié)果分析

由圖3(a)可見,深度為0.5～1.5 cm,14～17 cm處為顱骨,由于顱骨不含硼,故沒有硼劑量率,這兩處的相對偏差為0;深度為4.5～7 cm時,分布曲線上的陡然變化是由于該處為腫瘤區(qū)域,硼含量較高,導(dǎo)致硼劑量率顯著增加,有助于在臨床中更精確地確定腫瘤的位置。由圖3(b)可見,相對偏差保持在BNCT臨床允許的5%以內(nèi)[19-20]。由圖3(a)、圖4(a)和圖6(a)可見,分布曲線趨勢一致,這是由于次級光子來自于中子反應(yīng),中子劑量越多的地方產(chǎn)生的光子也越多,在頭部的淺表組織產(chǎn)生的損傷最大;由圖3(b)、圖4(b)和圖6(b)可見,相對偏差也均保持在5%以內(nèi)。由圖5(a)可見,因中子源發(fā)射大部分為超熱中子,一開始就發(fā)生中子慢化到熱中子,隨著深度的加深,劑量率逐漸降低,深度小于12 cm時,相對偏差都基本保持在5%以內(nèi);深度大于12 cm時,為深穿透問題,蒙特卡羅程序計算時到達(dá)深處的粒子數(shù)較少,造成統(tǒng)計結(jié)果相對偏差較大。最后,將圖3-圖6的劑量率按照相應(yīng)的權(quán)重因子整合,得到總相對生物效應(yīng)劑量率,如圖7(a)所示。由圖7(b)可見,相對偏差均保持在5%以內(nèi)。總體上,由圖3-圖7可見,Magic的解析模型與體素模型的模擬計算結(jié)果和MCNP解析模型的結(jié)果符合較好,表明Magic的計算結(jié)果是正確的。

4 基于Magic的BNCT動態(tài)硼濃度分析

中子注量率及硼濃度共同決定了BNCT沉積在組織中的劑量,因此在計算劑量時,這二者對劑量評估的準(zhǔn)確性有很大的影響。目前在進(jìn)行理論BNCT的治療計劃評估時,對腫瘤和正常組織的硼濃度設(shè)置都是時間固定不變的,但實際給藥后組織和血液中的硼濃度會因新陳代謝隨時間而變化,因此為臨床劑量評估的準(zhǔn)確性,有必要輸運計算時建立10B的藥代動力學(xué)模型。本文基于已驗證正確性的蒙特卡羅劑量計算程序Magic定量研究動態(tài)硼濃度下BNCT劑量學(xué)規(guī)律,為臨床應(yīng)用的安全性提供理論基礎(chǔ)。

目前,在BNCT含硼藥物BSH和BPA上,所使用的血硼濃度藥代動力學(xué)連續(xù)模型為兩房室及三房室模型,而通常組織硼濃度是通過血硼濃度和靜態(tài)組織與血液硼濃度之比的乘積作為實際組織硼濃度的代替物來建模的[21]。本文模擬含腫瘤(膠質(zhì)腫瘤)Snyder模型,選取BPA-F介導(dǎo)的多形性膠質(zhì)母細(xì)胞瘤BNCT情況進(jìn)行分析研究,圖8為硼藥(BPA-F)靜脈注射后顱內(nèi)血液中10B的開放式二房室模型。

其中：中央室(central compartment)代表血液和注射分布良好的組織;周邊室(peripheral compartment)代表注射較差的組織,如大腦;I為藥物BPA-F注射的速率,mg·kg-1·min-1;C1,C2分別為各室的10B濃度;V1,V2分別為各室的分布容積;K12,K21,K10為速率常數(shù),min-1。

結(jié)合二室模型,并以哈佛-麻省理工學(xué)院I期臨床試驗患者的藥代動力學(xué)研究數(shù)據(jù)[22]做相應(yīng)推導(dǎo),得到中央室血硼濃度C1隨時間的變化關(guān)系,如圖9所示。根據(jù)圖9,在BPA-F介導(dǎo)的多形性膠質(zhì)母細(xì)胞瘤BNCT中,腫瘤與血液的硼濃度之比通常假定為3.5∶1,腦部與血液的硼濃度之比假定為1∶1,頭皮與血液的硼濃度之比假定為1.5∶1,從而利用血硼濃度變化曲線來預(yù)測其他組織中的硼濃度變化[22]。本文選取上述曲線硼濃度上升階段5個時間點,下降階段10個時間點,共15個時間點,結(jié)合各組織倍數(shù)關(guān)系,利用Magic計算10B藥代動力學(xué)模型下腫瘤的總相對生物效應(yīng)劑量率隨時間的變化關(guān)系,如圖10所示。

圖9 中央室血硼濃度隨時間的變化關(guān)系

圖10 10B藥代動力學(xué)模型下腫瘤的總相對生物效應(yīng)劑量率隨時間的變化關(guān)系

圖10中陰影部分的時間段,可視為BNCT治療時持續(xù)照射的時間,紅色陰影部分為建立了10B的藥代動力學(xué)模型,考慮動態(tài)硼濃度變化后實際所受的劑量;藍(lán)色陰影處為目前理論BNCT穩(wěn)態(tài)下硼濃度所受的劑量。由圖10可直觀比較出二者的差異,對陰影面積(劑量值)進(jìn)行偏差分析,二者陰影區(qū)域相對偏差為11.78%。因此,在BNCT臨床治療過程中,若不考慮患者體內(nèi)硼濃度的動態(tài)變化過程,則可能導(dǎo)致實際受照劑量與處方劑量存在較大偏差,影響治療效果。

4 結(jié)論

本文基于自主開發(fā)的BNCT蒙卡劑量計算程序Magic對國際Snyder修正模型進(jìn)行劑量計算分析,通過MCNP得到深度-劑量率結(jié)果為標(biāo)準(zhǔn)進(jìn)行對比,得出Magic程序計算的結(jié)果與參考文獻(xiàn)結(jié)果相符,硼劑量率、熱中子劑量率、超熱中子及快中子劑量率、次級光子劑量率和總的相對生物效應(yīng)劑量率的相對偏差均小于5%,滿足臨床治療對計算精度的要求,驗證了Magic程序正確性,同時證明Magic程序具有構(gòu)建毫米量級千萬網(wǎng)格數(shù)模型的承載力。基于Magic建立10B的開放式二房室藥代動力學(xué)模型,對動態(tài)硼濃度下BNCT劑量學(xué)進(jìn)行定量研究,根據(jù)血硼濃度變化曲線及各組織倍數(shù)關(guān)系,計算得出不同時間點下各硼濃度的腫瘤隨時間變化的總相對生物效應(yīng)劑量率變化擬合曲線。選取BNCT持續(xù)照射時間下,考慮動態(tài)硼濃度變化后,實際的所受的真劑量值及穩(wěn)態(tài)下的硼濃度所受的假劑量值,計算得出二者的相對偏差為11.78%,可為臨床更加準(zhǔn)確的劑量計算及下一步完善含時蒙特卡羅劑量計算程序Magic提供參考。