以原子級幾何模型模擬脫氧核糖核酸（DNA）分子損傷的遍歷算法目前存在著計算時間長、收斂速度慢、計算設備要求高等缺點，因此本文提出一種基于能量沉積點和羥自由基（·OH）位置信息的含噪聲密度聚類（DBSCAN）算法，該算法結合概率統計方法，能快速地得到 DNA 鏈斷裂產額，幫助進一步研究 DNA 集簇性損傷的變化規律。首先，課題組依托放射生物學蒙特卡羅模擬軟件 Geant4-DNA，模擬質子和次級粒子在細胞核中的輸運以及它們對水分子電離激發的過程，獲得粒子能量沉積點和羥自由基的位置信息。然后在簡化的細胞核幾何模型中，利用損傷發生概率函數計算得到 DNA 損傷點的空間分布數據集，確定損傷點密度，由 DSBCAN 算法計算 DNA 單鏈斷裂（SSB）和雙鏈斷裂（DSB）產額。最后分析 DNA 鏈斷裂產額隨傳能線密度（LET）值變化的趨勢，總結出不同損傷簇的變化規律。實驗結果表明，與遍歷算法相比，本文新算法模擬速度快且模擬精度高，模擬結果同其他實驗和模擬方法相比較，具備良好的一致性。實現了在分子水平上快速獲得質子輻射致 DNA 集簇性損傷更為精細的信息，為輻射生物損傷機制的研究提供了一種重要和有力的研究手段。
引用本文： 唐菁, 張鵬程, 肖欽峰, 李杰, 桂志國. 基于含噪聲密度聚類算法的質子輻射誘導的 DNA 集簇性損傷研究. 生物醫學工程學雜志, 2019, 36(4): 633-642. doi: 10.7507/1001-5515.201812008 復制
引言
20 世紀 60 年代，Kuzin[1]的研究證明，細胞核中的脫氧核糖核酸（deoxyribonucleic acid，DNA）損傷會影響遺傳信息的復制和轉錄，如果不能及時修復損傷，則會引起基因突變、染色體畸變、細胞失活等嚴重的生物學后果。因此，DNA 分子損傷是輻射生物效應研究的一個重要課題，而蒙特卡羅徑跡結構法是研究 DNA 損傷最主要的模擬方法。目前國際上已有納米劑量級別的蒙特卡羅模擬程序，包括：蒙特卡羅徑跡軟件 KURBUC（Karolinska Institutet，瑞典）和粒子徑跡軟件 PARTRAC（Institute of Radiation Protection，德國）[2-3]，但都未公開發布。一套由歐洲核子研究組織（European organization for nuclear research，CERN）開發的開源幾何徑跡模擬軟件 Geant4（CERN，歐洲），能夠模擬包含電子、質子，重離子等各種粒子在物質中的輸運過程。而用于 DNA 的放射生物學蒙特卡羅模擬軟件 Geant4-DNA（National Institute for Nuclear and Particle Physics：IN2P3，法國）則是幾何徑跡模擬軟件 Geant4（CERN，歐洲）的擴展包，可以研究電離輻射在 DNA 上引起的早期生物損傷。在模擬中，DNA 損傷的空間分布與能量沉積點及羥自由基（·OH）有著密切關系。目前，大多數研究均基于原子級的 DNA 幾何模型，每一個損傷簇的獲取都會遍歷細胞核中所有原子，并計算能量沉積點或·OH 與周圍原子的距離[2-3]，因此在解決實際問題中，從徑跡結構計算到化學反應，再到 DNA 鏈斷裂計算的整個過程模擬非常耗時，需要在復雜的計算機集群上運行。于是有一些研究小組從概率統計和數據挖掘的角度探索新的方法，如 Nikjoo 等[4]利用 K 均值算法（K-means）研究粒子的能量沉積簇。Francis 等[5]利用 Ester 等[6]提出的含噪聲密度聚類（density-based spatial clustering of applications with noise，DBSCAN）算法研究 DNA 損傷簇，其在計算時間和收斂方面非常高效。
研究對比顯示，DBSCAN 算法與 K-means 算法不同，DBSCAN 算法不需要提前設置簇的數量，且不受簇形狀的影響。Francis 等[5]將 DBSCAN 算法應用于放射生物學中，通過電離粒子引起的能量沉積分布和損傷發生概率函數，確定損傷聚集的簇位置，得到 DNA 雙鏈斷裂（double strand break，DSB）與單鏈斷裂（single strand break，SSB）的比率。Dos Santos 等[7]在不同細胞核的幾何模型上，利用 DBSCAN 算法加速獲得 DSB 和 SSB 產額（yield），評估 DNA 密度與放射性損傷的相關性。然而上述研究都只考慮了粒子放射后的物理過程，得到的 DNA 損傷只與粒子徑跡中的能量沉積點相關，忽略了水分子的電離激發產生的化學產物對 DNA 的影響。實驗證明，大部分 DNA 鏈斷裂是由化學過程中·OH 攻擊 DNA 分子的間接損傷誘發的[8]。因此在蒙特卡羅模擬中，化學過程是不能忽略的，它對損傷簇的產額有著極其重要的影響。本文基于放射生物學蒙特卡羅模擬軟件 Geant4-DNA（IN2P3，法國），在已有研究基礎上，對 DBSCAN 算法的應用進行擴展，考慮損傷不僅與粒子的能量沉積點相關，還與·OH 的位置有關。通過模擬獲得水輻解產物的分布，由·OH 誘導的損傷發生概率計算 DNA 的間接損傷產額，研究 DSB 與 SSB 的比率以及 DSB 產額在加入化學模塊后的變化規律。因為目前的技術難以在真實的實驗環境中獲取更多的信息，因此本文通過與其他研究小組的實驗和模擬方法相對比，驗證本文方法的準確性。
1 方法
1.1 蒙特卡羅徑跡結構方法
de Nardo 等[9]和 Grosswendt 等[10]的實驗得到粒子與納米級媒質分子相互作用后能量沉積點的空間分布，而此能量沉積的簇分布成為研究粒子放射最基本且最關鍵的信息。用蒙特卡羅徑跡結構方法模擬質子及產生的次級粒子在輸運過程中發生電離和激發的物理過程，將能量轉移到 DNA 或水分子上，獲得能量沉積點的分布，其中沉積在 DNA 上的能量會誘發 DNA 的直接損傷。
本文假設直接損傷發生概率隨能量沉積值的增加而增加[5]，如式（1）所示，采用線性損傷概率函數的方法：當能量 E 低于下限值 Emin = 5 eV 時，誘發損傷的概率 P 為零；概率隨能量線性增加，達到或超過上限值 Emax = 37.5 eV 時，概率 P 為 1。
 |
另一部分轉移到水分子的能量使得水分子電離或激發，在能量沉積處生成·OH、H·、H2、H2O2、H3O+、
共 6 種輻解產物。參照粒子徑跡軟件 PARTRAC（Institute of Radiation Protection，德國）中的參數值[11]，設置不同水輻解過程的分支比、各產物相互反應比例系數。基于粒子—連續體表示法（particle-continuum）[12]，將每一個化學反應物看作包含在水連續體中的一個粒子，按照單位時間步長 1 ps，計算每一個粒子的位置。然后根據 Ballarini 等[13]研究的方法確定反應半徑，當反應物之間的距離小于反應半徑時發生化學反應，反應物被反應產物所替代。最后這些粒子隨時間逐漸擴散，呈現出自由基清除效應，化學反應基本完成。整個過程時長 1 ns。
化學過程中 DNA 間接損傷由·OH 攻擊 DNA 分子產生。Milligan 等[14]在實驗中測得，造成 DNA 鏈斷裂的·OH 數量占·OH 總產額的 13%，即總產額 20% 的·OH 與糖磷酸組反應，其中 65% 的反應會使得 DNA 發生斷裂[15]，因此本文將間接損傷的損傷概率 Rind 設為固定值 0.65。由此根據能量沉積點和·OH 的空間分布，采用不同的損傷發生概率函數計算得到 DNA 損傷點的空間分布數據集。
1.2 聚類算法與概率統計
針對模擬中生成的位置信息數據量大的特點，本文采用 DBSCAN 的加速算法來提高計算的效率。DBSCAN 算法是一種基于密度的聚類算法，它基于一組鄰域（neighborhood）參數（鄰域 ε，最小樣本個數 m）來刻畫樣本分布的緊密程度。給定數據集 D，算法原理如圖 1 所示。
圖1
DBSCAN 算法原理
Figure1.
Principle of DBSCAN algorithm
（1）首先找出數據集 D 中各樣本點的 ε 鄰域，如果樣本點 O 的 ε 鄰域至少包含 m 個樣本，則 O 是一個核心對象（core object），確定 D 中所有的核心對象以組成一個核心對象集合 Ω；
（2）接著從 Ω 隨機選取一個核心對象作為種子（seed），找出由其密度可達（density-reachable）的樣本生成聚類簇 C1；
（3）然后將 C1 中包含的核心對象從 Ω 中去除，再從更新后的 Ω 中隨機選取另一個核心對象作為種子生成下一個聚類簇 C2。上述過程不斷重復，直至 Ω 為空，迭代完成，生成聚類簇集合 C = {C1，C2，
，Ck}。沒有被任何聚類簇包含的點，即為噪聲點（noise）。
本文的鄰域參數（ε，m）分別代表了可以形成損傷簇的相鄰點間的最大距離和形成 DSB 損傷簇的最小 SSB 數量，它們的取值分別是 3.2 nm 和 2[5]。由此所獲得的簇損傷分為兩種：SSB 和在 DNA 兩個不同螺旋鏈上斷裂形成的 DSB。具體算法流程如圖 2 所示。
圖2
DBSCAN 算法流程
Figure2.
Flow chart of DBSCAN algorithm
本文采用的是近似幾何模型，且模擬是假設在均勻的水中進行，因此需要設置 Pp 和 Pc 兩個參數來表示粒子或化學產物擊中敏感區域的均勻概率。Francis 等[16]曾得出粒子電離轉變為鏈斷裂的概率在 9%～12% 之間，其他影響占 4% 的結論。因此 Francis 等[16]的 Pp 是 16%。而這個方法不能準確地體現出化學過程中水輻解產物造成鏈斷裂的影響。本文計算粒子電離轉變為鏈斷裂的概率 Pp 如式（2）所示。
 |
參考 Meylan 等[17]堿基對的幾何數據，細胞核內 6 × 109個堿基對（base pair，bp）的總體積約為 3.57 μm3，含水合層 DNA 體積是 28.57 μm3，細胞核的總體積為 732 μm3。因此細胞核中含水合層 DNA 的體積與細胞核體積的比值為 0.039。因為考慮化學過程的影響，所以本文粒子擊中敏感區域造成鏈斷裂的均勻概率要低于文獻[16]中 9% 的下限值。經過調試，發現 Pp = 5.5% 時模擬效果最佳，即將 0.039 的體積比乘以 1.4 倍的調節系數，并且考慮到能量沉積點同時作用在水分子上，即生成化學產物并擊中敏感區域的概率與 Pp 密切相關，因此 Pc 也設為 5.5%。
2 模擬
在放射生物學蒙特卡羅模擬軟件 Geant4-DNA（IN2P3，法國）平臺上構建軸半徑依次為 9.85、7.1、2.5 μm 的橢圓球體狀的成纖維細胞核模型。模擬質子初始能量在 0.5～50 MeV 之間，初始粒子的入射位置為橢球體細胞核表面上隨機的一個點，入射方向面向球心，如圖 3 所示。圖右上角的局部放大圖中，不同的點代表不同粒子和化學產物的位置信息。
圖3
成纖維細胞核的幾何模型
Figure3.
Geometric model of a human fibroblast nucleus
計算 DNA 損傷簇的流程圖如圖 4 所示。采用多線程同時并發模擬多個事件，每個事件單獨地模擬物理過程和化學過程，記錄 DNA 損傷點的數量，判斷并匯總 SSB 和 DSB 總個數，計算 SSB 產額（以符號 YSSB 表示）和 DSB 產額（以符號 YDSB 表示）如式（3）、（4）所示。
圖4
DNA 損傷的流程圖
Figure4.
Flow chart of DNA damage simulation
 |
 |
其中 NSSB 和 NDSB 分別是 SSB 和 DSB 的總個數；Dose 是質子在細胞核內的吸收劑量，單位是戈瑞（gray，Gy）；Nbp 為細胞核內所含堿基對，為 6 × 109 bp；YSSB 和 YDSB 分別是 SSB 和 DSB 產額，單位是 Gy?1·bp?1·10?9。
簡化的幾何模型沒有考慮細胞核內的組蛋白以及 DNA 分子的纏繞、扭曲等信息，因此生成的 SSB 和 DSB 產額會增大相應的倍數[18]。本文假設這些信息對 SSB 和 DSB 的影響相同，即產額都增大了 10 倍，將 DBSCAN 算法計算的 SSB 和 DSB 產額都乘以 0.1 的調節因子來抵消沒有考慮細胞核內各種組分的影響，最后進行模擬結果的對比與分析。
3 結果與分析
模擬生成的 DNA 損傷產額結果如表 1 所示，其中直接損傷總產額 Ydir 是由能量沉積誘導的 SSB 和 DSB 產額之和，間接損傷總產額 Yind 是由·OH 攻擊 DNA 生成的 SSB 和 DSB 產額之和，損傷總產額 Ytot 是 Ydir 與 Yind 之和。
表1
不同能量的 DNA 損傷產額
Table1.
DNA damage yields of different energies
3.1 SSB 和 DSB 的產額對比
如圖 5 所示為 SSB 和 DSB 產額依傳能線密度（linear energy transfer，LET）變化的對比圖。LET 表示每單位徑跡長度介質吸收的能量。分別選用 Friedland 團隊開發的粒子徑跡軟件 PARTRAC（Institute of Radiation Protection，德國）[19-20]以及其他研究人員如 Francis 等[5]、Meylan 等[17]、Lampe 等[21-22]和 Rosales 等[23]的模擬結果進行對比。本文的結果與這些研究基本吻合。存在的差異主要是因為模擬的細胞核幾何體模型、間接損傷發生概率、化學反應時間等設置的不同。
圖5
SSB 和 DSB 產額對比圖
Figure5.
Comparison chart of SSB and DSB yields
粒子徑跡軟件 PARTRAC（Institute of Radiation Protection，德國）是 DNA 集簇性損傷模擬研究的主要參考對象。本文的 SSB 產額在 10 keV·μm?1后，略高于 Friedland 等[19]于 2003 年所得的研究結果。因為后者詳細的細胞核幾何模型不僅包含了 DNA 聚集度信息，還包含了組蛋白。Nygren 等[18]的實驗結果總結，DNA 聚集度對 SSB 和 DSB 的影響相同，而組蛋白的去除使得 SSB 的增長比 DSB 要高 2 倍，而本文的假設并沒有考慮組蛋白的影響差異，因此 SSB 會略高。Friedland 等[20]于 2017 年指出，粒子徑跡軟件 PARTRAC（Institute of Radiation Protection，德國）最新的細胞核模型因其更復雜的 DNA 結構和媒質密度的不同，導致 SSB 產額總體降低。
相對于 Francis 等[5]的 DBSCAN 算法，本文增加了粒子對水分子電離激發后發生的化學過程，·OH 與糖磷酸組反應使得 DNA 鏈斷裂的總數上升，SSB 和 DSB 產額有不同程度的增加。但當 LET 值在 26～40 keV·μm?1之間時，本文的 SSB 產額比 Francis 等[5]的下降得更快，自由基清除效應產生，即對于高 LET 值的輻射粒子，水輻解生成的自由基密度增加，彼此間相互反應的概率也相應增加，DNA 損傷減少。
對比而言，Meylan 等[17]在其構建的 DNA 幾何模型中，6 個球形幾何體組成一個核苷酸組合，包括 2 個脫氧核糖，2 個磷酸和一對堿基。根據文獻[24]的研究，脫氧核糖的 5 個氫原子反應點，·OH 誘導損傷的主要占據 2 個，因此該課題組設定的間接損傷發生概率為 0.4，小于本文的 0.65。另外 Meylan 等[17]誘導 DNA 鏈斷裂的能量值是一個固定值 17.5 eV，化學階段時間范圍截止到 2.5 ns，這些參數設置都是導致 DSB 產額較低的原因。
Lampe 等[21-22]借鑒了 Meylan 等[17]的 DNA 幾何模型及 0.4 的間接損傷概率，但為了加速化學過程的模擬，設置了一個距離參數 rkill，即在反應物發生化學反應前，距離 DNA 分子 rkill 以外的自由基都將提前被清除，因此 SSB 和 DSB 產額都較低。
3.2 直接損傷和間接損傷對鏈斷裂的影響
本文進一步研究了直接損傷和間接損傷對 DNA 鏈斷裂的影響。如圖 6 所示，直接損傷的產額不受 LET 值變化的影響，基本保持在 60 Gy?1·bp?1·10?9上下。當 LET 值高于 10 keV·μm?1的質子照射時，DNA 斷裂總產額隨 LET 增加而降低，這種減少是由間接作用誘導的鏈斷裂產額明顯降低引起的。因此，隨著 LET 的增加，直接損傷產額對 DNA 鏈斷裂總產額的占比從低 LET 輻射的 32% 上升到了高 LET 輻射的 37%。這種趨勢在 Friedland 等[19]和 Nikjoo 等[25]的研究中也出現過。
圖6
直接損傷和間接損傷對 DNA 鏈斷裂的影響
Figure6.
Influence of direct and indirect damages on DNA strand breaks
3.3 DSB 產額隨 LET 變化的規律
DSB 損傷是輻射所致生物效應中最重要的原初損傷，質子等高 LET 輻射粒子由于具有布拉格（Bragg）峰這樣獨特的物理性質，因此誘導生成的 DSB 變化規律與 X 射線等低 LET 輻射有著本質的不同。如圖 7 所示，可以看出 DSB 產額總體上隨著 LET 的增加而增加。但在 1～10 keV·μm?1的區間內，本文的 DSB 產額出現了明顯的先下降再上升的趨勢。這說明了受 LET 影響的 DSB 產額值不僅與未參加活性產物重組過程的自由基數量有關，還與自由基聚集和電離作用的總體水平相關[8]。隨著 LET 的增加，較高電離密度的局部區域變得越來越密集，增加了電離簇和·OH 簇的密度，于是逐漸密集的空間位置提高了·OH 與其他活性產物的重組機會，減少了·OH 數量，降低了 DSB 的產額。但是 LET 的繼續增大，更大尺寸的電離簇促使粒子能量沉積更加密集，電離產生更多的·OH 攻擊 DNA 分子，·OH 與其他活性產物的重組不再是影響 DSB 的重要原因，DSB 產額值升高，因此產生了 DSB 產額先降低后升高的凹曲線。
圖7
DSB 產額隨 LET 變化
Figure7.
The yields of DSB versus LET
3.4 DSB/SSB 變化曲線
在脈沖凝膠電泳分析實驗中，研究人員通過比較直接損傷引起的鏈斷裂產率來估計間接損傷的相對貢獻，DSB/SSB 比值就成為一個有用的研究參數。Leloup 等[8]分別測量了在含有 2 mmol 和 200 mmol 甘油溶液中被照射質粒的 DNA 損傷產額。其中，在含有 2 mmol 甘油的低清除效應環境中，·OH 間接作用貢獻了大約 99% 的鏈斷裂損傷；而含有 200 mmol 甘油的環境類似于細胞環境，直接損傷約占總損傷的 35%。
如圖 8 所示，Leloup 等[8]的 2 mmol 甘油實驗和本文模擬得到的間接作用 DSB/SSB 曲線最陡，斜率逐漸增長。由前面提到的間接損傷總額隨著 LET 值的增加而呈下降的趨勢，本文得出結論：在輻射的間接作用中，隨著 LET 值的增加，DSB 產額逐漸增加，SSB 產額則逐漸減少。圖 8 中，直接作用的 DSB/SSB 曲線要平緩一些，這與直接損傷不受 LET 值的影響而產額總量基本保持不變有關。此外，本文同時考慮能量沉積誘導的直接損傷和·OH 誘導的間接損傷后的 DSB/SSB 曲線變化趨勢與 Leloup 等[8]的 200 mmol 甘油實驗結果相似。
圖8
DSB/SSB 比率曲線圖
Figure8.
DSB/SSB ratio
4 討論
本文通過能量沉積點和·OH 位置信息研究 DNA 集簇性損傷，由于大多數生物物理化學模型都是從半經驗基礎開始的，簡單的幾何更有助于在納米尺度上突出相互作用模式的重要性以及 DNA 中的集簇性損傷[26-28]。因此本文基于放射生物學蒙特卡羅模擬軟件 Geant4-DNA（IN2P3，法國），構造簡化的細胞核幾何模型，并利用軟件中的 G4EmDNAPhysics 和 G4EmDNAChemistry 函數構造出模擬的物理和化學模型[29-33]。為了獲得良好的統計取樣，首先保證初始質子發射位置隨機均勻，且根據 Nikjoo 等[34]基于平均頻率分布的特定能量均值比較法所得出的結論決定模擬事件數，本文模擬了 1 000 次質子輸運徑跡，使得每個質子在單位微米之內的損傷點數量誤差控制在 2% 以內。
模擬結果顯示，在細胞核環境中，直接作用的能量沉積對直接損傷產額和間接損傷位置的決定有著重要作用；·OH 間接作用則對 SSB 和 DSB 產額以及鏈斷裂復雜性的增加有明顯影響。由于當前放射生物學蒙特卡羅模擬軟件 Geant4-DNA（IN2P3，法國）的化學過程版本還存在一定的局限性，例如：忽略了高 LET 值的自由基快速重組過程，不能完整體現高濃度產物之間的反應過程[3, 30]。因此，如圖 5 所示的當 LET 高于 40 keV·μm?1時，SSB 產額出現了與其他結果不一樣的向上趨勢，以及如圖 6 所示的間接損傷對 DNA 鏈斷裂影響的研究中，本文的模擬并沒有如粒子徑跡軟件 PARTRAC（Institute of Radiation Protection，德國）模型[19]那樣下降至 50%。
相對于高精度的原子級 DNA 幾何模型，本文的簡化幾何模型還不能如實體現出 DNA 聚集程度對鏈斷裂產額的影響，在實驗中，DNA 的折疊狀態和組蛋白是·OH 清除能力急劇下降的主要因素之一[35]。本文假設上述信息對 SSB 和 DSB 的影響相同，同時選用了 0.1 的產額調節因子，然而實際上 DNA 的幾何形態對 SSB 和 DSB 的影響是有區別的。本課題組下一步的研究會區分這種不同。
本文 DSB 產額比實驗值高的原因，需要進一步擴展聚類算法，展開對 DSB 復雜度及其碎片分布的研究，特別是在實驗中觀察到的非隨機 DSB 分布現象[36]，有助于為質子、重離子等高 LET 值致 DNA 和生物體損傷機制提供支撐，為將來 DNA 修復模擬的研究提供可能性。
5 結論
本文在前人的基礎上，補充了化學過程中·OH 對 DNA 的間接損傷過程，提出了一種基于能量沉積點和·OH 位置信息的 DBSCAN 算法，為了快速得到 SSB、DSB 的產額，沒有考慮細胞核內具體的 DNA 組成信息，但是利用 DBSCAN 算法以及在模擬流程中選擇的參數值保證了本文與其他實驗結果的一致性。并且，這些參數是在可以解釋或證明的實際范圍內選擇的。本文由于補充了相應的化學過程，能夠研究不同能量質子的徑跡結構對 DNA 直接和間接損傷的影響，準確地再現了隨著 LET 的增加，自由基清除效應對 DNA 損傷簇的影響，特別是對 DSB 產額實現了更加準確的評估，為將來進一步使用 DBSCAN 算法研究相對生物學效應、DNA 修復模擬等提供基礎。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
20 世紀 60 年代，Kuzin[1]的研究證明，細胞核中的脫氧核糖核酸（deoxyribonucleic acid，DNA）損傷會影響遺傳信息的復制和轉錄，如果不能及時修復損傷，則會引起基因突變、染色體畸變、細胞失活等嚴重的生物學后果。因此，DNA 分子損傷是輻射生物效應研究的一個重要課題，而蒙特卡羅徑跡結構法是研究 DNA 損傷最主要的模擬方法。目前國際上已有納米劑量級別的蒙特卡羅模擬程序，包括：蒙特卡羅徑跡軟件 KURBUC（Karolinska Institutet，瑞典）和粒子徑跡軟件 PARTRAC（Institute of Radiation Protection，德國）[2-3]，但都未公開發布。一套由歐洲核子研究組織（European organization for nuclear research，CERN）開發的開源幾何徑跡模擬軟件 Geant4（CERN，歐洲），能夠模擬包含電子、質子，重離子等各種粒子在物質中的輸運過程。而用于 DNA 的放射生物學蒙特卡羅模擬軟件 Geant4-DNA（National Institute for Nuclear and Particle Physics：IN2P3，法國）則是幾何徑跡模擬軟件 Geant4（CERN，歐洲）的擴展包，可以研究電離輻射在 DNA 上引起的早期生物損傷。在模擬中，DNA 損傷的空間分布與能量沉積點及羥自由基（·OH）有著密切關系。目前，大多數研究均基于原子級的 DNA 幾何模型，每一個損傷簇的獲取都會遍歷細胞核中所有原子，并計算能量沉積點或·OH 與周圍原子的距離[2-3]，因此在解決實際問題中，從徑跡結構計算到化學反應，再到 DNA 鏈斷裂計算的整個過程模擬非常耗時，需要在復雜的計算機集群上運行。于是有一些研究小組從概率統計和數據挖掘的角度探索新的方法，如 Nikjoo 等[4]利用 K 均值算法（K-means）研究粒子的能量沉積簇。Francis 等[5]利用 Ester 等[6]提出的含噪聲密度聚類（density-based spatial clustering of applications with noise，DBSCAN）算法研究 DNA 損傷簇，其在計算時間和收斂方面非常高效。
研究對比顯示，DBSCAN 算法與 K-means 算法不同，DBSCAN 算法不需要提前設置簇的數量，且不受簇形狀的影響。Francis 等[5]將 DBSCAN 算法應用于放射生物學中，通過電離粒子引起的能量沉積分布和損傷發生概率函數，確定損傷聚集的簇位置，得到 DNA 雙鏈斷裂（double strand break，DSB）與單鏈斷裂（single strand break，SSB）的比率。Dos Santos 等[7]在不同細胞核的幾何模型上，利用 DBSCAN 算法加速獲得 DSB 和 SSB 產額（yield），評估 DNA 密度與放射性損傷的相關性。然而上述研究都只考慮了粒子放射后的物理過程，得到的 DNA 損傷只與粒子徑跡中的能量沉積點相關，忽略了水分子的電離激發產生的化學產物對 DNA 的影響。實驗證明，大部分 DNA 鏈斷裂是由化學過程中·OH 攻擊 DNA 分子的間接損傷誘發的[8]。因此在蒙特卡羅模擬中，化學過程是不能忽略的，它對損傷簇的產額有著極其重要的影響。本文基于放射生物學蒙特卡羅模擬軟件 Geant4-DNA（IN2P3，法國），在已有研究基礎上，對 DBSCAN 算法的應用進行擴展，考慮損傷不僅與粒子的能量沉積點相關，還與·OH 的位置有關。通過模擬獲得水輻解產物的分布，由·OH 誘導的損傷發生概率計算 DNA 的間接損傷產額，研究 DSB 與 SSB 的比率以及 DSB 產額在加入化學模塊后的變化規律。因為目前的技術難以在真實的實驗環境中獲取更多的信息，因此本文通過與其他研究小組的實驗和模擬方法相對比，驗證本文方法的準確性。
1 方法
1.1 蒙特卡羅徑跡結構方法
de Nardo 等[9]和 Grosswendt 等[10]的實驗得到粒子與納米級媒質分子相互作用后能量沉積點的空間分布，而此能量沉積的簇分布成為研究粒子放射最基本且最關鍵的信息。用蒙特卡羅徑跡結構方法模擬質子及產生的次級粒子在輸運過程中發生電離和激發的物理過程，將能量轉移到 DNA 或水分子上，獲得能量沉積點的分布，其中沉積在 DNA 上的能量會誘發 DNA 的直接損傷。
本文假設直接損傷發生概率隨能量沉積值的增加而增加[5]，如式（1）所示，采用線性損傷概率函數的方法：當能量 E 低于下限值 Emin = 5 eV 時，誘發損傷的概率 P 為零；概率隨能量線性增加，達到或超過上限值 Emax = 37.5 eV 時，概率 P 為 1。
|
另一部分轉移到水分子的能量使得水分子電離或激發，在能量沉積處生成·OH、H·、H2、H2O2、H3O+、 共 6 種輻解產物。參照粒子徑跡軟件 PARTRAC（Institute of Radiation Protection，德國）中的參數值[11]，設置不同水輻解過程的分支比、各產物相互反應比例系數。基于粒子—連續體表示法（particle-continuum）[12]，將每一個化學反應物看作包含在水連續體中的一個粒子，按照單位時間步長 1 ps，計算每一個粒子的位置。然后根據 Ballarini 等[13]研究的方法確定反應半徑，當反應物之間的距離小于反應半徑時發生化學反應，反應物被反應產物所替代。最后這些粒子隨時間逐漸擴散，呈現出自由基清除效應，化學反應基本完成。整個過程時長 1 ns。
化學過程中 DNA 間接損傷由·OH 攻擊 DNA 分子產生。Milligan 等[14]在實驗中測得，造成 DNA 鏈斷裂的·OH 數量占·OH 總產額的 13%，即總產額 20% 的·OH 與糖磷酸組反應，其中 65% 的反應會使得 DNA 發生斷裂[15]，因此本文將間接損傷的損傷概率 Rind 設為固定值 0.65。由此根據能量沉積點和·OH 的空間分布，采用不同的損傷發生概率函數計算得到 DNA 損傷點的空間分布數據集。
1.2 聚類算法與概率統計
針對模擬中生成的位置信息數據量大的特點，本文采用 DBSCAN 的加速算法來提高計算的效率。DBSCAN 算法是一種基于密度的聚類算法，它基于一組鄰域（neighborhood）參數（鄰域 ε，最小樣本個數 m）來刻畫樣本分布的緊密程度。給定數據集 D，算法原理如圖 1 所示。
圖1
DBSCAN 算法原理
Figure1.
Principle of DBSCAN algorithm
（1）首先找出數據集 D 中各樣本點的 ε 鄰域，如果樣本點 O 的 ε 鄰域至少包含 m 個樣本，則 O 是一個核心對象（core object），確定 D 中所有的核心對象以組成一個核心對象集合 Ω；
（2）接著從 Ω 隨機選取一個核心對象作為種子（seed），找出由其密度可達（density-reachable）的樣本生成聚類簇 C1；
（3）然后將 C1 中包含的核心對象從 Ω 中去除，再從更新后的 Ω 中隨機選取另一個核心對象作為種子生成下一個聚類簇 C2。上述過程不斷重復，直至 Ω 為空，迭代完成，生成聚類簇集合 C = {C1，C2，，Ck}。沒有被任何聚類簇包含的點，即為噪聲點（noise）。
本文的鄰域參數（ε，m）分別代表了可以形成損傷簇的相鄰點間的最大距離和形成 DSB 損傷簇的最小 SSB 數量，它們的取值分別是 3.2 nm 和 2[5]。由此所獲得的簇損傷分為兩種：SSB 和在 DNA 兩個不同螺旋鏈上斷裂形成的 DSB。具體算法流程如圖 2 所示。
圖2
DBSCAN 算法流程
Figure2.
Flow chart of DBSCAN algorithm
本文采用的是近似幾何模型，且模擬是假設在均勻的水中進行，因此需要設置 Pp 和 Pc 兩個參數來表示粒子或化學產物擊中敏感區域的均勻概率。Francis 等[16]曾得出粒子電離轉變為鏈斷裂的概率在 9%～12% 之間，其他影響占 4% 的結論。因此 Francis 等[16]的 Pp 是 16%。而這個方法不能準確地體現出化學過程中水輻解產物造成鏈斷裂的影響。本文計算粒子電離轉變為鏈斷裂的概率 Pp 如式（2）所示。
|
參考 Meylan 等[17]堿基對的幾何數據，細胞核內 6 × 109個堿基對（base pair，bp）的總體積約為 3.57 μm3，含水合層 DNA 體積是 28.57 μm3，細胞核的總體積為 732 μm3。因此細胞核中含水合層 DNA 的體積與細胞核體積的比值為 0.039。因為考慮化學過程的影響，所以本文粒子擊中敏感區域造成鏈斷裂的均勻概率要低于文獻[16]中 9% 的下限值。經過調試，發現 Pp = 5.5% 時模擬效果最佳，即將 0.039 的體積比乘以 1.4 倍的調節系數，并且考慮到能量沉積點同時作用在水分子上，即生成化學產物并擊中敏感區域的概率與 Pp 密切相關，因此 Pc 也設為 5.5%。
2 模擬
在放射生物學蒙特卡羅模擬軟件 Geant4-DNA（IN2P3，法國）平臺上構建軸半徑依次為 9.85、7.1、2.5 μm 的橢圓球體狀的成纖維細胞核模型。模擬質子初始能量在 0.5～50 MeV 之間，初始粒子的入射位置為橢球體細胞核表面上隨機的一個點，入射方向面向球心，如圖 3 所示。圖右上角的局部放大圖中，不同的點代表不同粒子和化學產物的位置信息。
圖3
成纖維細胞核的幾何模型
Figure3.
Geometric model of a human fibroblast nucleus
計算 DNA 損傷簇的流程圖如圖 4 所示。采用多線程同時并發模擬多個事件，每個事件單獨地模擬物理過程和化學過程，記錄 DNA 損傷點的數量，判斷并匯總 SSB 和 DSB 總個數，計算 SSB 產額（以符號 YSSB 表示）和 DSB 產額（以符號 YDSB 表示）如式（3）、（4）所示。
圖4
DNA 損傷的流程圖
Figure4.
Flow chart of DNA damage simulation
|
|
其中 NSSB 和 NDSB 分別是 SSB 和 DSB 的總個數；Dose 是質子在細胞核內的吸收劑量，單位是戈瑞（gray，Gy）；Nbp 為細胞核內所含堿基對，為 6 × 109 bp；YSSB 和 YDSB 分別是 SSB 和 DSB 產額，單位是 Gy?1·bp?1·10?9。
簡化的幾何模型沒有考慮細胞核內的組蛋白以及 DNA 分子的纏繞、扭曲等信息，因此生成的 SSB 和 DSB 產額會增大相應的倍數[18]。本文假設這些信息對 SSB 和 DSB 的影響相同，即產額都增大了 10 倍，將 DBSCAN 算法計算的 SSB 和 DSB 產額都乘以 0.1 的調節因子來抵消沒有考慮細胞核內各種組分的影響，最后進行模擬結果的對比與分析。
3 結果與分析
模擬生成的 DNA 損傷產額結果如表 1 所示，其中直接損傷總產額 Ydir 是由能量沉積誘導的 SSB 和 DSB 產額之和，間接損傷總產額 Yind 是由·OH 攻擊 DNA 生成的 SSB 和 DSB 產額之和，損傷總產額 Ytot 是 Ydir 與 Yind 之和。
表1
不同能量的 DNA 損傷產額
Table1.
DNA damage yields of different energies
3.1 SSB 和 DSB 的產額對比
如圖 5 所示為 SSB 和 DSB 產額依傳能線密度（linear energy transfer，LET）變化的對比圖。LET 表示每單位徑跡長度介質吸收的能量。分別選用 Friedland 團隊開發的粒子徑跡軟件 PARTRAC（Institute of Radiation Protection，德國）[19-20]以及其他研究人員如 Francis 等[5]、Meylan 等[17]、Lampe 等[21-22]和 Rosales 等[23]的模擬結果進行對比。本文的結果與這些研究基本吻合。存在的差異主要是因為模擬的細胞核幾何體模型、間接損傷發生概率、化學反應時間等設置的不同。
圖5
SSB 和 DSB 產額對比圖
Figure5.
Comparison chart of SSB and DSB yields
粒子徑跡軟件 PARTRAC（Institute of Radiation Protection，德國）是 DNA 集簇性損傷模擬研究的主要參考對象。本文的 SSB 產額在 10 keV·μm?1后，略高于 Friedland 等[19]于 2003 年所得的研究結果。因為后者詳細的細胞核幾何模型不僅包含了 DNA 聚集度信息，還包含了組蛋白。Nygren 等[18]的實驗結果總結，DNA 聚集度對 SSB 和 DSB 的影響相同，而組蛋白的去除使得 SSB 的增長比 DSB 要高 2 倍，而本文的假設并沒有考慮組蛋白的影響差異，因此 SSB 會略高。Friedland 等[20]于 2017 年指出，粒子徑跡軟件 PARTRAC（Institute of Radiation Protection，德國）最新的細胞核模型因其更復雜的 DNA 結構和媒質密度的不同，導致 SSB 產額總體降低。
相對于 Francis 等[5]的 DBSCAN 算法，本文增加了粒子對水分子電離激發后發生的化學過程，·OH 與糖磷酸組反應使得 DNA 鏈斷裂的總數上升，SSB 和 DSB 產額有不同程度的增加。但當 LET 值在 26～40 keV·μm?1之間時，本文的 SSB 產額比 Francis 等[5]的下降得更快，自由基清除效應產生，即對于高 LET 值的輻射粒子，水輻解生成的自由基密度增加，彼此間相互反應的概率也相應增加，DNA 損傷減少。
對比而言，Meylan 等[17]在其構建的 DNA 幾何模型中，6 個球形幾何體組成一個核苷酸組合，包括 2 個脫氧核糖，2 個磷酸和一對堿基。根據文獻[24]的研究，脫氧核糖的 5 個氫原子反應點，·OH 誘導損傷的主要占據 2 個，因此該課題組設定的間接損傷發生概率為 0.4，小于本文的 0.65。另外 Meylan 等[17]誘導 DNA 鏈斷裂的能量值是一個固定值 17.5 eV，化學階段時間范圍截止到 2.5 ns，這些參數設置都是導致 DSB 產額較低的原因。
Lampe 等[21-22]借鑒了 Meylan 等[17]的 DNA 幾何模型及 0.4 的間接損傷概率，但為了加速化學過程的模擬，設置了一個距離參數 rkill，即在反應物發生化學反應前，距離 DNA 分子 rkill 以外的自由基都將提前被清除，因此 SSB 和 DSB 產額都較低。
3.2 直接損傷和間接損傷對鏈斷裂的影響
本文進一步研究了直接損傷和間接損傷對 DNA 鏈斷裂的影響。如圖 6 所示，直接損傷的產額不受 LET 值變化的影響，基本保持在 60 Gy?1·bp?1·10?9上下。當 LET 值高于 10 keV·μm?1的質子照射時，DNA 斷裂總產額隨 LET 增加而降低，這種減少是由間接作用誘導的鏈斷裂產額明顯降低引起的。因此，隨著 LET 的增加，直接損傷產額對 DNA 鏈斷裂總產額的占比從低 LET 輻射的 32% 上升到了高 LET 輻射的 37%。這種趨勢在 Friedland 等[19]和 Nikjoo 等[25]的研究中也出現過。
圖6
直接損傷和間接損傷對 DNA 鏈斷裂的影響
Figure6.
Influence of direct and indirect damages on DNA strand breaks
3.3 DSB 產額隨 LET 變化的規律
DSB 損傷是輻射所致生物效應中最重要的原初損傷，質子等高 LET 輻射粒子由于具有布拉格（Bragg）峰這樣獨特的物理性質，因此誘導生成的 DSB 變化規律與 X 射線等低 LET 輻射有著本質的不同。如圖 7 所示，可以看出 DSB 產額總體上隨著 LET 的增加而增加。但在 1～10 keV·μm?1的區間內，本文的 DSB 產額出現了明顯的先下降再上升的趨勢。這說明了受 LET 影響的 DSB 產額值不僅與未參加活性產物重組過程的自由基數量有關，還與自由基聚集和電離作用的總體水平相關[8]。隨著 LET 的增加，較高電離密度的局部區域變得越來越密集，增加了電離簇和·OH 簇的密度，于是逐漸密集的空間位置提高了·OH 與其他活性產物的重組機會，減少了·OH 數量，降低了 DSB 的產額。但是 LET 的繼續增大，更大尺寸的電離簇促使粒子能量沉積更加密集，電離產生更多的·OH 攻擊 DNA 分子，·OH 與其他活性產物的重組不再是影響 DSB 的重要原因，DSB 產額值升高，因此產生了 DSB 產額先降低后升高的凹曲線。
圖7
DSB 產額隨 LET 變化
Figure7.
The yields of DSB versus LET
3.4 DSB/SSB 變化曲線
在脈沖凝膠電泳分析實驗中，研究人員通過比較直接損傷引起的鏈斷裂產率來估計間接損傷的相對貢獻，DSB/SSB 比值就成為一個有用的研究參數。Leloup 等[8]分別測量了在含有 2 mmol 和 200 mmol 甘油溶液中被照射質粒的 DNA 損傷產額。其中，在含有 2 mmol 甘油的低清除效應環境中，·OH 間接作用貢獻了大約 99% 的鏈斷裂損傷；而含有 200 mmol 甘油的環境類似于細胞環境，直接損傷約占總損傷的 35%。
如圖 8 所示，Leloup 等[8]的 2 mmol 甘油實驗和本文模擬得到的間接作用 DSB/SSB 曲線最陡，斜率逐漸增長。由前面提到的間接損傷總額隨著 LET 值的增加而呈下降的趨勢，本文得出結論：在輻射的間接作用中，隨著 LET 值的增加，DSB 產額逐漸增加，SSB 產額則逐漸減少。圖 8 中，直接作用的 DSB/SSB 曲線要平緩一些，這與直接損傷不受 LET 值的影響而產額總量基本保持不變有關。此外，本文同時考慮能量沉積誘導的直接損傷和·OH 誘導的間接損傷后的 DSB/SSB 曲線變化趨勢與 Leloup 等[8]的 200 mmol 甘油實驗結果相似。
圖8
DSB/SSB 比率曲線圖
Figure8.
DSB/SSB ratio
4 討論
本文通過能量沉積點和·OH 位置信息研究 DNA 集簇性損傷，由于大多數生物物理化學模型都是從半經驗基礎開始的，簡單的幾何更有助于在納米尺度上突出相互作用模式的重要性以及 DNA 中的集簇性損傷[26-28]。因此本文基于放射生物學蒙特卡羅模擬軟件 Geant4-DNA（IN2P3，法國），構造簡化的細胞核幾何模型，并利用軟件中的 G4EmDNAPhysics 和 G4EmDNAChemistry 函數構造出模擬的物理和化學模型[29-33]。為了獲得良好的統計取樣，首先保證初始質子發射位置隨機均勻，且根據 Nikjoo 等[34]基于平均頻率分布的特定能量均值比較法所得出的結論決定模擬事件數，本文模擬了 1 000 次質子輸運徑跡，使得每個質子在單位微米之內的損傷點數量誤差控制在 2% 以內。
模擬結果顯示，在細胞核環境中，直接作用的能量沉積對直接損傷產額和間接損傷位置的決定有著重要作用；·OH 間接作用則對 SSB 和 DSB 產額以及鏈斷裂復雜性的增加有明顯影響。由于當前放射生物學蒙特卡羅模擬軟件 Geant4-DNA（IN2P3，法國）的化學過程版本還存在一定的局限性，例如：忽略了高 LET 值的自由基快速重組過程，不能完整體現高濃度產物之間的反應過程[3, 30]。因此，如圖 5 所示的當 LET 高于 40 keV·μm?1時，SSB 產額出現了與其他結果不一樣的向上趨勢，以及如圖 6 所示的間接損傷對 DNA 鏈斷裂影響的研究中，本文的模擬并沒有如粒子徑跡軟件 PARTRAC（Institute of Radiation Protection，德國）模型[19]那樣下降至 50%。
相對于高精度的原子級 DNA 幾何模型，本文的簡化幾何模型還不能如實體現出 DNA 聚集程度對鏈斷裂產額的影響，在實驗中，DNA 的折疊狀態和組蛋白是·OH 清除能力急劇下降的主要因素之一[35]。本文假設上述信息對 SSB 和 DSB 的影響相同，同時選用了 0.1 的產額調節因子，然而實際上 DNA 的幾何形態對 SSB 和 DSB 的影響是有區別的。本課題組下一步的研究會區分這種不同。
本文 DSB 產額比實驗值高的原因，需要進一步擴展聚類算法，展開對 DSB 復雜度及其碎片分布的研究，特別是在實驗中觀察到的非隨機 DSB 分布現象[36]，有助于為質子、重離子等高 LET 值致 DNA 和生物體損傷機制提供支撐，為將來 DNA 修復模擬的研究提供可能性。
5 結論
本文在前人的基礎上，補充了化學過程中·OH 對 DNA 的間接損傷過程，提出了一種基于能量沉積點和·OH 位置信息的 DBSCAN 算法，為了快速得到 SSB、DSB 的產額，沒有考慮細胞核內具體的 DNA 組成信息，但是利用 DBSCAN 算法以及在模擬流程中選擇的參數值保證了本文與其他實驗結果的一致性。并且，這些參數是在可以解釋或證明的實際范圍內選擇的。本文由于補充了相應的化學過程，能夠研究不同能量質子的徑跡結構對 DNA 直接和間接損傷的影響，準確地再現了隨著 LET 的增加，自由基清除效應對 DNA 損傷簇的影響，特別是對 DSB 產額實現了更加準確的評估，為將來進一步使用 DBSCAN 算法研究相對生物學效應、DNA 修復模擬等提供基礎。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
