基于機器學習軟骨細胞的時間依賴性力學行為及本構參數反演1)

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

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

引言

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

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

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

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

1 有限元建模

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

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

2 MSnHS 本構模型

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

式中,C10為與剪切模量相關的本構模型參數;W為單位體積應變能;D1為材料的可壓縮性;J為材料的體積比;為左柯西?格應變張量的第一應變不變量.

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

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

3 雙向深度神經網絡與隨機森林

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

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

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

圖3 隨機森林模型預測結構示意圖Fig.3 The structure of the RF model

4 預測結果與討論

4.1 數據準備

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

式中,X*為標準化后的值,X為原始值;Xmax,Xmin分別為特征數據范圍內的最大值和最小值.

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

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

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

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

4.2 貝葉斯超參數優化

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

圖5 TW-Deepnets 模型架構Fig.5 TW-Deepnets model architecture

表1 Bayesian 優化算法確定的模型超參數Table 1 Hyperparameter values obtained by Bayesian optimization

4.3 參數預測結果分析

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

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

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

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

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

圖6 本構參數預測性能檢驗Fig.6 Test of constitutive parameters prediction

圖7 MSnHS 本構參數特征重要性占比Fig.7 Importance ratio of MSnHS constitutive parameters

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

表2 MSnHS 本構模型參數預測結果Table 2 Parameters prediction results of MSnHS

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

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

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

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

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

5 結論

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

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

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

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