基于機器學(xué)習(xí)軟骨細(xì)胞的時間依賴性力學(xué)行為及本構(gòu)參數(shù)反演1)

魏新宇 桑建兵 張睿琳 王靜遠(yuǎn) 劉寶友

(河北工業(yè)大學(xué)機械工程學(xué)院,天津 300401)

引言

軟骨細(xì)胞稀疏分布于構(gòu)成軟骨的細(xì)胞外基質(zhì)ECM 內(nèi),負(fù)責(zé)ECM 的維護(hù)和修復(fù)[1-2],其感知到的機械信號能夠調(diào)節(jié)關(guān)節(jié)軟骨的生理機能[3-4].此外,已有研究表明,軟骨細(xì)胞在不同的病理或生理過程中表現(xiàn)出一定的特征反應(yīng),如骨關(guān)節(jié)炎軟骨細(xì)胞較正常軟骨細(xì)胞的黏彈性特性存在明顯差異[5],以及在自然生長過程中的軟骨細(xì)胞整體黏彈性和恢復(fù)變形的能力受年齡因素影響發(fā)生改變[6].因此,準(zhǔn)確量化單個軟骨細(xì)胞的力學(xué)特性及其對機械負(fù)荷的響應(yīng)對于揭示軟骨生物力學(xué)特性與骨性關(guān)節(jié)炎病理十分重要.

近年來,力學(xué)表征實驗技術(shù)的發(fā)展為推導(dǎo)描述軟骨細(xì)胞材料特性的本構(gòu)模型提供了數(shù)值基礎(chǔ).研究表明,軟骨細(xì)胞表現(xiàn)出高度非線性彈性[7-8].Florea 等[9]基于納米壓痕實驗和計算建模,揭示了軟骨細(xì)胞隨時間變化的力學(xué)行為的關(guān)鍵本構(gòu)參數(shù).文獻(xiàn)[10-11]基于微管抽吸實驗的數(shù)值模擬及理論分析探討了軟骨細(xì)胞半空間模型、不可壓縮球體模型和可壓縮球體模型的彈性參數(shù)的差異.文獻(xiàn)[12-13]基于無側(cè)限壓縮試驗探究了軟骨細(xì)胞孔隙彈性本構(gòu)模型的滲透率數(shù)值差異性及參數(shù)可靠性問題.關(guān)于軟骨細(xì)胞本構(gòu)模型的研究取得了一定程度的進(jìn)展,但其力學(xué)響應(yīng)與計算模型本構(gòu)參數(shù)之間的高度復(fù)雜非線性使得軟骨細(xì)胞的本構(gòu)參數(shù)反求仍是一個具有挑戰(zhàn)性的問題.

人工智能的飛速發(fā)展為其他學(xué)科領(lǐng)域提供了解決問題的新思路,能夠高效模擬輸入和輸出變量之間的非線性相關(guān)性的機器學(xué)習(xí)方法逐步應(yīng)用于生物力學(xué)領(lǐng)域[14-15].研究表明,利用可靠的正演模型來進(jìn)行反演分析是系統(tǒng)識別材料性質(zhì)的有效手段[16-17].李洋等[18]集成有限元與智能算法實現(xiàn)了骨骼肌本構(gòu)參數(shù)的反演方法研究.Arbabi 等[19]基于人工神經(jīng)網(wǎng)絡(luò)和軟骨孔隙彈性壓痕有限元模型,確定了軟骨的力學(xué)和物理性能.Wang 等[20]基于紅細(xì)胞的光鑷?yán)爝^程的有限元模型利用BP 神經(jīng)網(wǎng)絡(luò)實現(xiàn)了紅細(xì)胞膜本構(gòu)參數(shù)識別,探索了原子力顯微鏡壓痕下紅細(xì)胞的局部變形機理.然而,取決于參數(shù)組的維度以及生物代理模型的復(fù)雜度,如何在有限的實驗及仿真成本下提高生物材料參數(shù)反演的效率和精度仍是值得探討的問題.Liu[21]提出了一種用于正反向力學(xué)問題實時計算的雙向深度網(wǎng)絡(luò)TW-Deepnets (twoway deepnets)模型,該模型被證明可以在基于物理定律的小樣本問題上取得更高的預(yù)測精度[22].同時,在有關(guān)中小型結(jié)構(gòu)表格數(shù)據(jù)方面,決策樹算法具有良好的預(yù)測效果,隨機森林(random forest,RF)是一種基于決策樹的回歸技術(shù).在計算生物學(xué)工作中,RF 已被證明在處理小樣本量,高維特征空間和復(fù)雜數(shù)據(jù)結(jié)構(gòu)方面具有獨特的優(yōu)勢[23].

因此,本論文基于軟骨細(xì)胞無側(cè)限壓縮實驗有限元模型以及TW-Deepnets 模型、RF 模型,分別提出了兩種反演方法,用于識別軟骨細(xì)胞黏彈性本構(gòu)參數(shù).同時,結(jié)合數(shù)據(jù)采樣對比和智能算法超參數(shù)優(yōu)化過程,擬合Nguyen 等[24]的實驗曲線來驗證和對比所提出的反演方法的有效性,并通過使用統(tǒng)計指標(biāo)R2對比了這兩種方法對骨細(xì)胞本構(gòu)參數(shù)的預(yù)測能力.

1 有限元建模

基于Nguyen 等[24]的實驗,使用有限元軟件ABAQUS(6.14)模擬了單個軟骨細(xì)胞的無側(cè)限壓縮松弛實驗,通過微操作技術(shù)在懸浮液中對分離出的單軟骨細(xì)胞的機械響應(yīng)進(jìn)行了實驗研究.因此,本論文中的骨細(xì)胞有限元模型沒有考慮細(xì)胞與基質(zhì)間的黏附效應(yīng).軟骨細(xì)胞建模為一個直徑為9 μm 的球形,該值是基于Nguyen 等[24]測定的牛軟骨細(xì)胞的平均高度.如圖1(a)所示,建立的有限元模型由剛性壓縮板和球形細(xì)胞組成,考慮到幾何形狀和邊界條件的對稱性,采用了八分之一模型進(jìn)行計算,由2420個線性六面體單元和11 139 個節(jié)點組成,單元類型采用C3D8R 單元.沿對稱面施加對稱邊界條件,細(xì)胞和板之間定義無摩擦的面面接觸.在有限元分析過程中建立了兩個分析步: 第一個分析步為0.75 s,在這一分析步中,如圖1(b)剛性壓縮板勻速施加豎直向下的2.25 μm 的位移;第二個分析步為2.25 s,在這一分析步中,剛性壓縮板位移保持不變,模擬軟骨細(xì)胞應(yīng)力松弛過程.

圖1 軟骨細(xì)胞壓縮有限元模型Fig.1 Finite element model of chondrocyte compression

2 MSnHS 本構(gòu)模型

黏彈性模型是細(xì)胞力學(xué)研究中最常用的力學(xué)模型之一[25].Nguyen 等[26]提出了一種修正的非線性黏彈性本構(gòu)模型(MSnHS),解釋了單細(xì)胞的應(yīng)變率依賴行為,但研究只限于軟骨細(xì)胞的壓縮過程.本研究中,基于MSnHS 本構(gòu)結(jié)合無側(cè)限壓縮實驗的壓縮及松弛過程來探究單個軟骨細(xì)胞的時間依賴性全局力學(xué)性能,并分析細(xì)胞應(yīng)力松弛機制對各本構(gòu)參數(shù)的依賴性.軟骨細(xì)胞的非線性、彈性力學(xué)響應(yīng)采用可壓縮Neo Hookean 超彈性模型,該模型的應(yīng)變能密度表達(dá)式為

式中,C10為與剪切模量相關(guān)的本構(gòu)模型參數(shù);W為單位體積應(yīng)變能;D1為材料的可壓縮性;J為材料的體積比;為左柯西?格應(yīng)變張量的第一應(yīng)變不變量.

黏性模型的剪切松弛模量以及體積松弛模量的一項Prony 級數(shù)展開表達(dá)式為

其中,G(t)和K(t)為剪切松弛模量和體積松弛模量;G0和K0為瞬時剪切模量體積模量;g1,k1,τ1為Prony 級數(shù)展開式的材料常數(shù).將式(3)的黏性材料特性與式(1)描述的超彈性材料特性合并即為描述軟骨細(xì)胞的MSnHS 本構(gòu)模型.

3 雙向深度神經(jīng)網(wǎng)絡(luò)與隨機森林

雙向深度神經(jīng)網(wǎng)絡(luò)模型,基于物理定律的模型,如有限元方法FEM,光滑F(xiàn)EM (S-FEM)和無網(wǎng)格模型進(jìn)行訓(xùn)練,用于材料和結(jié)構(gòu)的正、反力學(xué)問題的實時計算[21].有限元系統(tǒng)矩陣的良好特性,使得有限元?人工智能深度網(wǎng)絡(luò)能夠有效地訓(xùn)練正向力學(xué)問題,為高成本的反向傳播提供了堅實的基礎(chǔ).在本論文中,基于有限元收集的數(shù)據(jù),正問題神經(jīng)網(wǎng)絡(luò)建立了一個由材料參數(shù)預(yù)測軟骨細(xì)胞壓縮接觸力隨時間歷程變化的代理模型,反問題神經(jīng)網(wǎng)絡(luò)則結(jié)合實驗數(shù)據(jù)及正問題代理模型進(jìn)行MSnHS 本構(gòu)參數(shù)反求,利用雙向深度神經(jīng)網(wǎng)絡(luò)進(jìn)行軟骨細(xì)胞參數(shù)反求問題流程圖如圖2 所示.

圖2 本研究中提出的本構(gòu)參數(shù)識別方法流程圖Fig.2 The flow chart of the parameter identification method proposed in this study

隨機森林是一種經(jīng)典的Bagging 模型,其弱學(xué)習(xí)器為決策樹模型,在回歸問題上具有強大的非線性關(guān)系學(xué)習(xí)能力[27-28].RF 回歸模型基于原始數(shù)據(jù)的隨機抽樣創(chuàng)建大量不相關(guān)的決策樹,并根據(jù)每個決策樹的平均值獲取最終結(jié)果.如圖3 所示,與其他算法不同,RF 回歸模型將原始數(shù)據(jù)通過隨機抽樣分為兩部分,即袋內(nèi)訓(xùn)練數(shù)據(jù)和袋外驗證數(shù)據(jù)OOB.使用OOB 數(shù)據(jù)而非外部數(shù)據(jù)用于評估學(xué)習(xí)性能,使得RF 模型具有更強的回歸預(yù)測能力.利用隨機森林進(jìn)行軟骨細(xì)胞參數(shù)反求問題流程圖如圖2(b)所示.

圖3 隨機森林模型預(yù)測結(jié)構(gòu)示意圖Fig.3 The structure of the RF model

4 預(yù)測結(jié)果與討論

4.1 數(shù)據(jù)準(zhǔn)備

MSnHS 本構(gòu)參數(shù)空間的必要參數(shù)為五個,分別為C10(取值范圍700～ 2630),D1(取值范圍50～2260),g1(取值范圍0.25～ 0.7),k1(取值范圍0～ 1),τ1(取值范圍0～ 1),參數(shù)范圍的確定基于相關(guān)文獻(xiàn)對軟骨細(xì)胞本構(gòu)模型的研究.研究表明,均勻采樣和拉丁超立方采樣(LHS)在采樣空間的均勻性及空間填充性上具有各自的優(yōu)勢[29].本研究中,為了提高基于采樣的分析效率、收斂性和魯棒性,使用將均勻采樣與LHS 結(jié)合的混合抽樣方法,即在參數(shù)空間內(nèi)分別采用LHS 及均勻抽樣法各采集200 組材料參數(shù)繼而合并得到采樣量為400 的材料參數(shù)組.接下來,基于軟骨細(xì)胞的有限元模型,獲取材料參數(shù)與細(xì)胞和壓板間的接觸響應(yīng)力的對應(yīng)數(shù)據(jù)集,用于搭建能夠反求本構(gòu)參數(shù)的TW-Deepnets 模型及RF 模型.為了消除特征間的量級差異性,使得回歸問題能夠更加準(zhǔn)確快速的收斂,本研究對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理

式中,X*為標(biāo)準(zhǔn)化后的值,X為原始值;Xmax,Xmin分別為特征數(shù)據(jù)范圍內(nèi)的最大值和最小值.

為了測試基于混合抽樣法所得數(shù)據(jù)集的分析效率,基于均勻抽樣法,拉丁超立方抽樣法以及混合抽樣法的數(shù)據(jù)集,對所搭建的軟骨細(xì)胞神經(jīng)網(wǎng)絡(luò)代理模型的預(yù)測效果進(jìn)行了對比.代理模型的輸入為本構(gòu)參數(shù),輸出接觸力,以均方誤差函數(shù)MSE為損失函數(shù),表達(dá)式如式(5)所示

式中,n為樣本量,ypred為預(yù)測值,yi為真實值.

如圖4,基于混合采樣數(shù)據(jù)的軟骨細(xì)胞壓縮實驗神經(jīng)網(wǎng)絡(luò)代理模型的預(yù)測誤差MSE明顯小于基于其他兩種單一采樣方法數(shù)據(jù)的預(yù)測誤差,說明本研究采用的混合采樣法在數(shù)據(jù)驅(qū)動預(yù)測問題中的高效性.

圖4 軟骨細(xì)胞代理模型的誤差對比Fig.4 Comparing the MSE of neural network proxy model under different sampling methods

4.2 貝葉斯超參數(shù)優(yōu)化

神經(jīng)網(wǎng)絡(luò)模型的主要超參數(shù)包括學(xué)習(xí)率、隱藏層數(shù)、隱藏單元數(shù)、激活函數(shù),通常使用經(jīng)驗公式或交叉驗證方法[30]對超參數(shù)進(jìn)行估計.本文基于貝葉斯(Bayesian)優(yōu)化構(gòu)造了神經(jīng)網(wǎng)絡(luò)超參數(shù)搜索空間模型,用高斯過程填充樣本點之間的區(qū)域,尋找驗證集上均方誤差函數(shù)MSE最小值的參數(shù)組合.神經(jīng)網(wǎng)絡(luò)最優(yōu)超參數(shù)組搜索空間為: 學(xué)習(xí)率(10?6～10?1)、隱藏層數(shù)(1～ 10)、隱藏單元數(shù)(5～ 50)、激活函數(shù)(Relu/Sigomid).圖5 為優(yōu)化后雙向神經(jīng)網(wǎng)絡(luò)的網(wǎng)絡(luò)架構(gòu).同樣使用貝葉斯優(yōu)化對隨機森林的重要超參數(shù)進(jìn)行優(yōu)化,隨機森林主要模型參數(shù)包括決策樹的數(shù)量、RF 最大深度、內(nèi)部節(jié)點分裂所需的最小樣本.表1 為Bayesian 優(yōu)化后的TW-Deepnets模型及RF 模型超參數(shù)值.

圖5 TW-Deepnets 模型架構(gòu)Fig.5 TW-Deepnets model architecture

表1 Bayesian 優(yōu)化算法確定的模型超參數(shù)Table 1 Hyperparameter values obtained by Bayesian optimization

4.3 參數(shù)預(yù)測結(jié)果分析

有限元收集到的初始樣本量以7:3 的比例對訓(xùn)練集和測試集進(jìn)行分配.TW-Deepnets 模型使用均方誤差(MSE)作為損失函數(shù),通過真實值與預(yù)測值的距離來指導(dǎo)模型的收斂方向.其中正問題神經(jīng)網(wǎng)絡(luò)在測試集上的MSE為2.076 2 × 10?6,反問題神經(jīng)網(wǎng)絡(luò)在測試集上的MSE為1.4 × 10?4,正反問題模型在測試集上的良好表現(xiàn)保證了模型的泛化能力及預(yù)測的準(zhǔn)確性.對于隨機森林模型的評估,采用均方根誤差函數(shù)(RMSE)測量預(yù)測值與真實值之間的偏差,RF 模型在測試集上的RMSE為0.041,偏差程度較小.均方根誤差函數(shù)公式為

式中,N為樣本量,ypred為預(yù)測值,yi為真實值.

為了更好地對比兩種模型對軟骨細(xì)胞本構(gòu)參數(shù)的預(yù)測能力,通過決定系數(shù)R-square(R2)來評估參數(shù)C10,D1,g1,k1,τ1真實值與預(yù)測值的擬合程度.R2利用數(shù)據(jù)的平均值作為誤差基準(zhǔn),觀察預(yù)測誤差的偏離程度.R2定義為

式中,yi為真實值,ypred為預(yù)測值,y為真實值的平均值.

圖6 為TW-Deepnets 模型和RF 模型對10 組參數(shù)預(yù)測的R2結(jié)果對比,兩種模型預(yù)測所得各個參數(shù)均能達(dá)到較高的擬合優(yōu)度(R2均值 >0.80),而TW-Deepnets 模型總體情況好于RF 模型.可以觀察到C10,g1的預(yù)測效果擬合優(yōu)度總是很高(R2>0.9),而D1,k1,τ1的擬合優(yōu)度會低一些,這說明在利用軟骨細(xì)胞壓縮實驗挖掘其全局力學(xué)性能與MSnHS 材料參數(shù)復(fù)雜的非線性關(guān)系中,C10,g1起重要影響,而τ1,k1的影響占比較小,D1最小,這與文獻(xiàn)[9]的研究結(jié)論一致.圖7 為RF 由參數(shù)預(yù)測壓縮力時的Scikitlearn 特征重要性,代表各個特征對于預(yù)測力的貢獻(xiàn)大小,即更加直觀地驗證了MSnHS 本構(gòu)各參數(shù)對描述軟骨細(xì)胞時間依賴性力學(xué)行為的影響占比.

圖6 本構(gòu)參數(shù)預(yù)測性能檢驗Fig.6 Test of constitutive parameters prediction

圖7 MSnHS 本構(gòu)參數(shù)特征重要性占比Fig.7 Importance ratio of MSnHS constitutive parameters

利用文獻(xiàn)[24]中0～3 s 內(nèi)單個軟骨細(xì)胞受到50%壓縮程度下與壓縮板間的實驗接觸力數(shù)據(jù),TWDeepnets 模型和RF 模型各預(yù)測得到一組MSnHS本構(gòu)參數(shù)值.表2 為兩個參數(shù)反求模型的本構(gòu)參數(shù)預(yù)測結(jié)果.將預(yù)測所得的本構(gòu)參數(shù)輸入有限元模型,對比有限元響應(yīng)曲線與實驗曲線的擬合程度來驗證預(yù)測所得參數(shù)的可靠性.圖8(a)為分別基于兩個參數(shù)反求模型預(yù)測所得有限元模型響應(yīng)曲線與實驗曲線對比.其中,由TW-Deepnets 模型預(yù)測所得的MSnHS 有限元響應(yīng)可以很好地和實驗曲線擬合(R2=0.987),RF 模型的預(yù)測響應(yīng)在壓縮階段可以與實驗曲線很好地擬合,而應(yīng)力松弛階段存在輕微的滯后現(xiàn)象,需要達(dá)到最終松弛狀態(tài)的時間更長,壓縮反作用力更小,與實驗曲線的擬合程度為R2=0.912.因此,本文提出的集成有限元模型和TW-Deepnets模型以及集成有限元模型和RF 模型進(jìn)行參數(shù)反求的過程都是預(yù)測軟骨細(xì)胞的MSnHS 本構(gòu)參數(shù)的有效方法,其中TW-Deepnets 模型表現(xiàn)出更高的準(zhǔn)確度.

表2 MSnHS 本構(gòu)模型參數(shù)預(yù)測結(jié)果Table 2 Parameters prediction results of MSnHS

將本文中利用TW-Deepnets 模型預(yù)測所得的軟骨細(xì)胞壓縮有限元響應(yīng)與文獻(xiàn)[13]利用傳統(tǒng)反有限元方法進(jìn)行參數(shù)預(yù)測所得的軟骨細(xì)胞壓縮力學(xué)響應(yīng)曲線進(jìn)行對比.從圖8(b)中可以直觀看出,無論是基于黏彈性本構(gòu)模型還是基于孔隙黏彈性本構(gòu)模型,傳統(tǒng)反有限元方法預(yù)測所得的響應(yīng)力最大值點都存在較大偏差,并且基于黏彈性本構(gòu)的預(yù)測模型對于細(xì)胞松弛過程中響應(yīng)力松弛速率與實驗存在一定差異.這說明,相較于傳統(tǒng)的反有限元方法,基于TWDeepnets 模型的參數(shù)反演方法能夠預(yù)測更加可靠的軟骨細(xì)胞本構(gòu)參數(shù),從而能夠建立更加準(zhǔn)確的計算模型進(jìn)行力學(xué)性能研究.同時,MSnHS 本構(gòu)模型也被證明能夠很好地對軟骨細(xì)胞的時間依賴性全局力學(xué)性能進(jìn)行描述.

圖8 模型預(yù)測與實驗的時間?力曲線對比Fig.8 Model predictions are compared with experimental time-force curves

軟骨細(xì)胞通常表現(xiàn)出應(yīng)變率依賴性行為,其力學(xué)響應(yīng)也與加載速率有關(guān),如圖9 為軟骨細(xì)胞有限元模型在不同壓縮加載速率下的力響應(yīng)曲線,壓縮速度越快,最終變形時細(xì)胞與壓板間的接觸力越大.然而在細(xì)胞松弛后,即在細(xì)胞保持?jǐn)?shù)秒后,所有力都達(dá)到相同的平衡力,這與Nguyen 等[24]的實驗結(jié)論一致,表明本研究中的軟骨細(xì)胞本構(gòu)關(guān)系及有限元模型能夠捕捉單軟骨細(xì)胞的應(yīng)變率依賴性行為.圖10 為軟骨細(xì)胞受到50%壓縮后3 s 應(yīng)力松弛前后的Mises 應(yīng)力云圖(單位為MPa).正如松弛前的應(yīng)力云圖所示,軟骨細(xì)胞受到壓縮后的最大應(yīng)力值點集中在細(xì)胞的中心,應(yīng)力值由中心向外圍不斷減小.隨著松弛時間的增加,細(xì)胞整體應(yīng)力值逐漸減小,在3 s 后達(dá)到最終狀態(tài),此時應(yīng)力最大處的應(yīng)力值減少了4.969 kPa.

圖9 軟骨細(xì)胞不同壓縮速度下的時間?力曲線Fig.9 Time-force curves of chondrocytes at different compression rates

圖10 軟骨細(xì)胞3 s 應(yīng)力松弛前后Mise 應(yīng)力云圖Fig.10 Mises stress cloud diagram of chondrocytes before and after 3 s stress relaxation

5 結(jié)論

本論文提出了分別利用TW-Deepnets 模型和RF 模型結(jié)合有限元系統(tǒng)來識別軟骨細(xì)胞的MSnHS本構(gòu)參數(shù)的反演方法.建立了軟骨細(xì)胞的無側(cè)限壓縮有限元模型用于數(shù)據(jù)收集,利用決定系數(shù)R2對兩種模型對各參數(shù)的預(yù)測性能進(jìn)行了對比評估,分析了參數(shù)特征重要性占比.最后,利用文獻(xiàn)[24]的實驗曲線驗證了反演方法的有效性,得出如下結(jié)論.

(1) 基于建立的軟骨細(xì)胞無側(cè)限壓縮實驗有限元模型,結(jié)合數(shù)據(jù)采樣對比和超參數(shù)優(yōu)化,利用TWDeepnets 模型和RF 模型確定了能夠描述單個軟骨細(xì)胞時間依賴性力學(xué)特性的本構(gòu)參數(shù)取值,預(yù)測所得有限元響應(yīng)與實驗曲線吻合良好.

(2) MSnHS 本構(gòu)模型能夠很好地描述包括壓縮過程和后續(xù)的應(yīng)力松弛過程中軟骨細(xì)胞的時間依賴性力學(xué)性能.其中,利用軟骨細(xì)胞無側(cè)限壓縮實驗挖掘其全局力學(xué)性能與MSnHS 材料參數(shù)復(fù)雜的非線性關(guān)系中,C10,g1起重要影響,而τ1,k1的影響占比較小,D1最小.

(3) TW-Deepnets 模型和RF 模型對于預(yù)測軟骨細(xì)胞黏彈性本構(gòu)參數(shù)都具有很大的潛力,就本研究來說,TW-Deepnets 模型表現(xiàn)出更高的準(zhǔn)確度和穩(wěn)定性.該方法也可進(jìn)一步推廣到其他維度更復(fù)雜的本構(gòu)參數(shù)反求及其他類型細(xì)胞在負(fù)載下更復(fù)雜的細(xì)胞微觀力學(xué)現(xiàn)象.