朱 健 白 曈 李寶生舒華忠 Antoine Simon Renaud de Crevoisier

(1東南大學(xué)影像科學(xué)與技術(shù)實(shí)驗(yàn)室, 南京 210096)(2山東省腫瘤防治研究院, 濟(jì)南 250117)(3雷恩第一大學(xué)信號與圖像實(shí)驗(yàn)室, 法國雷恩 35510)

腫瘤放療并發(fā)癥概率預(yù)測模型參數(shù)擬合方法

朱 健1,2白 曈2李寶生1,2舒華忠1Antoine Simon3Renaud de Crevoisier3

(1東南大學(xué)影像科學(xué)與技術(shù)實(shí)驗(yàn)室, 南京 210096)(2山東省腫瘤防治研究院, 濟(jì)南 250117)(3雷恩第一大學(xué)信號與圖像實(shí)驗(yàn)室, 法國雷恩 35510)

為了建立具有群體特異性的腫瘤放療NTCP預(yù)測模型,提出了一種模型參數(shù)擬合方法.首先,基于NTCP模型的特點(diǎn)構(gòu)建最大似然函數(shù);然后,分別采用確定性優(yōu)化方法和隨機(jī)性優(yōu)化方法對最大似然函數(shù)進(jìn)行優(yōu)化,分析優(yōu)化過程的時間成本及優(yōu)化結(jié)果,探討用于擬合NTCP模型參數(shù)的最優(yōu)方法.實(shí)驗(yàn)結(jié)果表明,用于擬合NTCP模型參數(shù)的最大似然函數(shù)是非凸的,存在局部最優(yōu)解;遺傳算法是一種最穩(wěn)定的最大似然函數(shù)優(yōu)化方法,其運(yùn)行時間比模擬退火算法短,而且可以在每次優(yōu)化結(jié)束后給出全局最優(yōu)解,以作為NTCP模型參數(shù).所提方法可以幫助腫瘤放療工作者在臨床隨訪數(shù)據(jù)的基礎(chǔ)上建立具有群體特異性的放療并發(fā)癥預(yù)測模型.

腫瘤;放射治療;并發(fā)癥;NTCP模型

正常組織并發(fā)癥概率(normal tissue complication probability, NTCP)模型是建立在放療中劑量-體積關(guān)系上的一種數(shù)學(xué)模型.該模型可以通過調(diào)整參數(shù)來描述不同器官在接受一定劑量照射后出現(xiàn)放療并發(fā)癥的概率,從而對放療毒性反應(yīng)進(jìn)行預(yù)測[1-3],也可根據(jù)模型預(yù)測結(jié)果對不同的治療方案進(jìn)行生物效應(yīng)量化對比[4-6].選取模型對放療并發(fā)癥進(jìn)行預(yù)測前,需要對NTCP模型進(jìn)行建模,根據(jù)所適用的患者群體出現(xiàn)并發(fā)癥的隨訪數(shù)據(jù),來擬合具有群體特異性的NTCP模型參數(shù)[7].據(jù)筆者所知,目前國內(nèi)外文獻(xiàn)未見針對NTCP建模方法的介紹.鑒于此,本文提出了一種針對NTCP模型的擬合方法.首先,構(gòu)建出最大似然函數(shù);然后,設(shè)計實(shí)驗(yàn)以尋找最優(yōu)的參數(shù)擬合方法.

1 最大似然函數(shù)的構(gòu)建

采用最大似然估計法來擬合NTCP模型參數(shù).最大似然估計法是一種統(tǒng)計方法,用于求解樣本集中相關(guān)概率密度函數(shù)的參數(shù).給定一個概率分布D,其概率密度函數(shù)(連續(xù)分布)或概率聚集函數(shù)(離散分布)為fD.從該分布中抽出k個采樣{x1,x2,…,xk},其概率為

p(x1,x2,…,xk)=fD(x1,x2,…,xkθ)

(1)

式中,θ為分布參數(shù).

盡管采樣數(shù)據(jù)來自于分布D,但θ可能未知.一旦獲得x1,x2,…,xk,便能從中找到一個關(guān)于θ的估計.最大似然估計法是指在所有可能的θ取值中,尋找使該采樣可能性最大化的解.

在數(shù)學(xué)上實(shí)現(xiàn)最大似然估計法,首先需要定義最大似然函數(shù)為

l(θ)=fD(x1,x2,…,xkθ)

(2)

使該采樣可能性最大的θ值即為其最大似然估計.

在NTCP模型參數(shù)擬合的應(yīng)用中,對于患者i,無論采用哪種模型,其NTCP值Ni均可表述為如下的函數(shù)形式[8]:

Ni=F(P,di,vi)

(3)

式中,P為模型參數(shù);di為微分后的劑量;vi為體積單元.

對于用來擬合模型參數(shù)的一組患者數(shù)據(jù),其對數(shù)似然函數(shù)可表示為

(4)

式中,Ri為參數(shù),如果患者i在臨床觀察或隨訪中出現(xiàn)放療并發(fā)癥,則Ri=1,否則Ri=0.

利用優(yōu)化方法將最大似然函數(shù)取最大值,便可擬合得到各模型參數(shù).

2 實(shí)驗(yàn)

通過對第1節(jié)中的最大似然函數(shù)進(jìn)行優(yōu)化,便可快速準(zhǔn)確地得到NTCP模型參數(shù)擬合結(jié)果.優(yōu)化方法可以分為兩大類:確定性優(yōu)化方法和隨機(jī)性優(yōu)化方法.前者的特點(diǎn)是速度快,但容易陷入局部極小;后者雖然不易陷入局部極小,從概率上可以收斂到全局最優(yōu),但是收斂速度比較慢[9].據(jù)筆者所知,目前國內(nèi)外文獻(xiàn)未見關(guān)于NTCP模型是否為凸函數(shù)的討論,即不確定其是否存在局部最小值.因此,考慮到確定性優(yōu)化方法在運(yùn)行速度方面的優(yōu)勢,選擇該方法開展實(shí)驗(yàn).

式(4)即為優(yōu)化的目標(biāo)函數(shù).考慮到該函數(shù)的復(fù)雜性,對其是否為凸函數(shù)的數(shù)學(xué)證明也會比較復(fù)雜.鑒于此,采用反證法進(jìn)行討論.即借助確定性優(yōu)化方法的特點(diǎn),從反面證明優(yōu)化目標(biāo)函數(shù)是否為凸函數(shù).若目標(biāo)函數(shù)為凸函數(shù),即解空間中不存在任何一個點(diǎn)使函數(shù)陷入局部最優(yōu),則無論取何值為解空間起始點(diǎn),均能找到一個唯一不變的最優(yōu)點(diǎn),該點(diǎn)即為全局最優(yōu)點(diǎn);若目標(biāo)函數(shù)為非凸函數(shù),則可能在某一次優(yōu)化過程中將局部最優(yōu)點(diǎn)作為返回結(jié)果.因此,采用確定性優(yōu)化方法時,如果在若干次重復(fù)的優(yōu)化過程中得到了不同的解,則表明該目標(biāo)函數(shù)存在局部最小值,即該函數(shù)為非凸函數(shù).

在文獻(xiàn)[10]的數(shù)據(jù)基礎(chǔ)上,以前列腺癌放療引起直腸晚期并發(fā)癥(大于等于二級損傷,LENT/SOMA標(biāo)準(zhǔn))為例,考察5年隨訪時間中257例患者出現(xiàn)該并發(fā)癥的情況.目標(biāo)函數(shù)中的NTCP模型以LKB模型為例[8],某正常器官受照射至出現(xiàn)50%并發(fā)癥概率所需的劑量T∈[50,100],體積效應(yīng)因子n∈(0,1),斜率因子m∈(0,1).實(shí)驗(yàn)的具體步驟如下:

① 選擇確定性優(yōu)化方法中的單純形法進(jìn)行優(yōu)化.每次優(yōu)化開始之前,對T,n,m隨機(jī)取值,作為單純形優(yōu)化的起始點(diǎn).

② 調(diào)用Matlab軟件運(yùn)行環(huán)境中的fminsearch函數(shù),實(shí)現(xiàn)對目標(biāo)函數(shù)的優(yōu)化.重復(fù)運(yùn)行10次,記錄每次優(yōu)化后得到的解、目標(biāo)函數(shù)值和運(yùn)行時間.

③ 在相同患者數(shù)據(jù)基礎(chǔ)上,利用模擬退火算法和遺傳算法對目標(biāo)函數(shù)分別優(yōu)化10次,求得全局最優(yōu)解,并與根據(jù)確定性優(yōu)化方法得到的解進(jìn)行對比.為保證隨機(jī)性優(yōu)化方法能獲得全局最優(yōu)解,將模擬退火算法的初始溫度設(shè)置為100 ℃,每次降溫的衰減系數(shù)為0.9,搜索步長為0.1,每一溫度下的迭代次數(shù)為104;為確保跳出局部極值、找到全局最優(yōu)解,將遺傳算法中的進(jìn)化代數(shù)設(shè)置為104.

3 結(jié)果與分析

將單純形算法重復(fù)運(yùn)行10次且每次都隨機(jī)選取優(yōu)化起始點(diǎn),優(yōu)化結(jié)果見表1.由表可知,在第8次和第9次優(yōu)化過程中,單純形法于開始階段便陷入了局部最優(yōu)解,得到了與其他情況下截然不同的解和最大似然函數(shù)值.因此,按照最大似然估計法構(gòu)成的優(yōu)化目標(biāo)函數(shù)是非凸函數(shù),存在局部最優(yōu)解,不適于用確定性優(yōu)化方法.

表1 單純形算法的優(yōu)化結(jié)果

當(dāng)n≈0.02,參數(shù)T和m在解空間變化時,最大似然函數(shù)取到最小值的解不是唯一的,而是一個區(qū)間.當(dāng)n≈0.02,T≈77.97時,參數(shù)m與最大似然函數(shù)值的關(guān)系見圖1.由圖1(a)可知,當(dāng)m=

(a) m∈(0,1)

(b) m→0時的局部放大圖

0.14時,最大似然函數(shù)取得全局最小值,此時的m值即為最大似然函數(shù)的全局最優(yōu)解;當(dāng)m→0時,曲線存在一個平臺,對其放大之后可見一個最大似然函數(shù)的局部最小值(見圖1(b)),此時的m值即為單純形算法第8次優(yōu)化后得到的局部最優(yōu)解.

綜上所述,從優(yōu)化結(jié)果和圖形表達(dá)2個方面證實(shí)了由NTCP模型構(gòu)成的最大似然函數(shù)作為優(yōu)化目標(biāo)函數(shù)是存在局部最小值的,該目標(biāo)函數(shù)為非凸函數(shù),不適于用確定性優(yōu)化方法.因此,需要考慮隨機(jī)性優(yōu)化方法.選取模擬退火算法和遺傳算法,重復(fù)運(yùn)行10次后的優(yōu)化結(jié)果分別見表2和表3.由表可知,模擬退火算法運(yùn)行時間較長,弱初始值選取不合適,則可能陷入局部最小值,溫度無法下降,從而導(dǎo)致優(yōu)化過程結(jié)束;遺傳算法的優(yōu)化結(jié)果則更穩(wěn)定,每次優(yōu)化后都能找到全局最優(yōu)解,且每次的運(yùn)行時間基本相同.因此,在對NTCP建模時,建議使用遺傳算法擬合模型參數(shù),確保在最短的時間內(nèi)獲得全局最優(yōu)解,使所建模型對并發(fā)癥預(yù)測更準(zhǔn)確.

表2 模擬退火算法的優(yōu)化結(jié)果

表3 遺傳算法的優(yōu)化結(jié)果

為了證實(shí)遺傳算法所得的優(yōu)化結(jié)果為全局最優(yōu)解,可通過全面搜索算法進(jìn)行驗(yàn)證.在搜索步長足夠小的情況下,全面搜索算法可以確保找到全局最優(yōu)解,但其完成一次全局優(yōu)化的時間較長.以優(yōu)化相同樣本為例,在運(yùn)行環(huán)境為Intel Core2 Duo CPU,T5670,1.80 GHz,2 GB內(nèi)存的計算平臺上,利用全面搜索算法,計算出最大似然函數(shù)值的時間約為0.023 6 s.令T∈[50,100],n∈(0,1),m∈(0,1),取值步長為0.01,則完成全局優(yōu)化需要進(jìn)行5×108次計算,大約需要13 d.因此, 全面搜索算法僅適于對其他算法進(jìn)行驗(yàn)證,不適用于常規(guī)的參數(shù)優(yōu)化.

4 結(jié)語

本文通過實(shí)驗(yàn)方法證明了擬合NTCP模型所構(gòu)建的最大似然函數(shù)是非凸函數(shù),存在局部最優(yōu)解;對比了不同優(yōu)化算法在擬合NTCP模型參數(shù)方面所表現(xiàn)出的特點(diǎn).實(shí)驗(yàn)結(jié)果表明,遺傳算法的優(yōu)化結(jié)果和運(yùn)行時間較穩(wěn)定,且每次優(yōu)化結(jié)束后都可以給出全局最優(yōu)解.因此,在對NTCP建模時,建議使用遺傳算法擬合模型參數(shù),確保在最短的時間內(nèi)獲得全局最優(yōu)解,使所建模型對并發(fā)癥預(yù)測更準(zhǔn)確.所提方法可以幫助腫瘤放療工作者在臨床隨訪數(shù)據(jù)的基礎(chǔ)上建立具有群體特異性的放療并發(fā)癥預(yù)測模型.

References)

[1]Tucker S L, Li M, Xu T, et al. Incorporating single-nucleotide polymorphisms into the Lyman model to improve prediction of radiation pneumonitis[J].IntJRadiatOncolBiolPhys, 2013, 85(1): 251-257.

[2]Strigari L, Pedicini P, D’Andrea M, et al. A new model for predicting acute mucosal toxicity in head-and-neck cancer patients undergoing radiotherapy with altered schedules[J].IntJRadiatOncolBiolPhys, 2012, 83(5): e697-e702.

[3]Gulliford S L, Partridge M, Sydes M R, et al. Parameters for the Lyman Kutcher Burman (LKB) model of normal tissue complication probability (NTCP) for specific rectal complications observed in clinical practise[J].RadiotherapyandOncology, 2012, 102(3): 347-351.

[4]Fellin F, Azzeroni R, Maggio A, et al. Helical tomotherapy and intensity modulated proton therapy in the treatment of dominant intraprostatic lesion: a treament planning comparison[J].RadiotherapyandOncology, 2013, 107(2): 207-212.

[5]Amin N P, Miften M, Thornton D, et al. Effect of induction chemotherapy on estimated risk of radiation pneumonitis in bulky non-small cell lung cancer[J].MedicalDosimetry, 2013, 38(3): 320-326.

[6]de Sanctis V, Bolzan C, D’Arienzo M, et al. Intensity modulated radiotherapy in early stage Hodgkin lymphoma patients: is it better than three dimensional conformal radiotherapy?[J].RadiationOncology, 2012, 7: 129-1-129-9.

[7]Zhu J, Zhang Z C, Li B S, et al. Analysis of acute radiation-induced esophagitis in non-small-cell lung cancer patients using the Lyman NTCP model[J].RadiotherapyandOncology, 2010, 97(3): 449-454.

[8]朱健,李寶生,舒華忠,等.正常組織并發(fā)癥概率模型綜述[J].中國生物醫(yī)學(xué)工程學(xué)報,2014,33(2):233-240. Zhu Jian, Li Baosheng, Shu Huazhong, et al. Review of normal tissues complication probability models[J].ChineseJournalofBiomedicalEngineering, 2014, 33(2): 233-240. (in Chinese)

[9]袁亞湘,孫文瑜.最優(yōu)化理論與方法[M].北京:科學(xué)出版社,2001:50-51.

[10]Zhu J, Simon A, Ospina J D, et al. Predictive models of bladder toxicity in prostate cancer radiotherapy[J].EurJCancer, 2011, 47(S1): S486.

Parameter fitting method of NTCP predictive model in radiation oncology

Zhu Jian1,2Bai Tong2Li Baosheng1,2Shu Huazhong1Antoine Simon3Renaud de Crevoisier3

(1Laboratory of Image Science and Technique, Southeast University, Nanjing 210096, China)(2Shandong Cancer Hospital and Institute, Jinan 250117, China)(3Laboratoire Traitement du Signal et de l’Image, Université de Rennes 1, Rennes 35510, France)

To establish the population specific NTCP (normal tissue complication probability) prediction model in radiation oncology, a parameter fitting method is proposed. First, the maximum likelihood function is constructed based on the characteristic of the NTCP model. Then, the deterministic optimization method and the stochastic optimization method are used to optimize the maximum likelihood function, respectively. The time cost and the optimization results are analyzed to find the better method for fitting the parameters of the NTCP model. The experimental results show that the maximum likelihood function for fitting the parameters of the NTCP model is non-convex, indicating that there exist the local optimal solutions. The genetic algorithm is the most stable optimization algorithm for fitting the NTCP model, and the running time is less than that of the simulated annealing algorithm. In this algorithm, the global optimal solutions, which are regarded as the parameters of the NTCP model, can be obtained after each optimization. The proposed method can help the researchers in radiation oncology establish the population specific NTCP predictive models based on clinical follow-ups.

tumor; radiotherapy; complication; NTCP (normal tissue complication probability) model

10.3969/j.issn.1001-0505.2015.02.011

2014-10-01. 作者簡介: 朱健(1980—)，男，博士，助理研究員；李寶生(聯(lián)系人)，男，博士，研究員，博士生導(dǎo)師，baoshli@yahoo.com.

國家自然科學(xué)基金資助項(xiàng)目(61271312，81272501，81301298)、國家重點(diǎn)基礎(chǔ)研究發(fā)展計劃(973計劃)資助項(xiàng)目(2011CB707904).

朱健,白曈,李寶生,等.腫瘤放療并發(fā)癥概率預(yù)測模型參數(shù)擬合方法[J].東南大學(xué)學(xué)報:自然科學(xué)版,2015,45(2):256-259.

10.3969/j.issn.1001-0505.2015.02.011

TP391.9

A

1001-0505(2015)02-0256-04