魏新宇 桑建兵 張睿琳 王靜遠 劉寶友

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

引言

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

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

人工智能的飛速發(fā)展為其他學(xué)科領(lǐng)域提供了解決問題的新思路,能夠高效模擬輸入和輸出變量之間的非線性相關(guān)性的機器學(xué)習(xí)方法逐步應(yīng)用于生物力學(xué)領(lǐng)域[14-15].研究表明,利用可靠的正演模型來進行反演分析是系統(tǒng)識別材料性質(zhì)的有效手段[16-17].李洋等[18]集成有限元與智能算法實現(xiàn)了骨骼肌本構(gòu)參數(shù)的反演方法研究.Arbabi 等[19]基于人工神經(jīng)網(wǎng)絡(luò)和軟骨孔隙彈性壓痕有限元模型,確定了軟骨的力學(xué)和物理性能.Wang 等[20]基于紅細胞的光鑷拉伸過程的有限元模型利用BP 神經(jīng)網(wǎng)絡(luò)實現(xiàn)了紅細胞膜本構(gòu)參數(shù)識別,探索了原子力顯微鏡壓痕下紅細胞的局部變形機理.然而,取決于參數(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].

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

1 有限元建模

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

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

2 MSnHS 本構(gòu)模型

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

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

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

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

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

雙向深度神經(jīng)網(wǎng)絡(luò)模型,基于物理定律的模型,如有限元方法FEM,光滑F(xiàn)EM (S-FEM)和無網(wǎng)格模型進行訓(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ù)測軟骨細胞壓縮接觸力隨時間歷程變化的代理模型,反問題神經(jīng)網(wǎng)絡(luò)則結(jié)合實驗數(shù)據(jù)及正問題代理模型進行MSnHS 本構(gòu)參數(shù)反求,利用雙向深度神經(jīng)網(wǎng)絡(luò)進行軟骨細胞參數(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ù)測能力.利用隨機森林進行軟骨細胞參數(shù)反求問題流程圖如圖2(b)所示.

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

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

4.1 數(shù)據(jù)準備

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)文獻對軟骨細胞本構(gòu)模型的研究.研究表明,均勻采樣和拉丁超立方采樣(LHS)在采樣空間的均勻性及空間填充性上具有各自的優(yōu)勢[29].本研究中,為了提高基于采樣的分析效率、收斂性和魯棒性,使用將均勻采樣與LHS 結(jié)合的混合抽樣方法,即在參數(shù)空間內(nèi)分別采用LHS 及均勻抽樣法各采集200 組材料參數(shù)繼而合并得到采樣量為400 的材料參數(shù)組.接下來,基于軟骨細胞的有限元模型,獲取材料參數(shù)與細胞和壓板間的接觸響應(yīng)力的對應(yīng)數(shù)據(jù)集,用于搭建能夠反求本構(gòu)參數(shù)的TW-Deepnets 模型及RF 模型.為了消除特征間的量級差異性,使得回歸問題能夠更加準確快速的收斂,本研究對數(shù)據(jù)進行標(biāo)準化處理

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

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

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

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

圖4 軟骨細胞代理模型的誤差對比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ù)進行估計.本文基于貝葉斯(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ù)進行優(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)練集和測試集進行分配.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ù)測的準確性.對于隨機森林模型的評估,采用均方根誤差函數(shù)(RMSE)測量預(yù)測值與真實值之間的偏差,RF 模型在測試集上的RMSE為0.041,偏差程度較小.均方根誤差函數(shù)公式為

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

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

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

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

利用文獻[24]中0～3 s 內(nèi)單個軟骨細胞受到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)象,需要達到最終松弛狀態(tài)的時間更長,壓縮反作用力更小,與實驗曲線的擬合程度為R2=0.912.因此,本文提出的集成有限元模型和TW-Deepnets模型以及集成有限元模型和RF 模型進行參數(shù)反求的過程都是預(yù)測軟骨細胞的MSnHS 本構(gòu)參數(shù)的有效方法,其中TW-Deepnets 模型表現(xiàn)出更高的準確度.

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

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

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

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

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

圖10 軟骨細胞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)來識別軟骨細胞的MSnHS本構(gòu)參數(shù)的反演方法.建立了軟骨細胞的無側(cè)限壓縮有限元模型用于數(shù)據(jù)收集,利用決定系數(shù)R2對兩種模型對各參數(shù)的預(yù)測性能進行了對比評估,分析了參數(shù)特征重要性占比.最后,利用文獻[24]的實驗曲線驗證了反演方法的有效性,得出如下結(jié)論.

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

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

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