我剛接觸Chaste那會(huì)兒差點(diǎn)被它“勸退”。作為一個(gè)主打細(xì)胞群體動(dòng)力學(xué)仿真的開源框架它的文檔散、示例雜、編譯方式又跟普通程序不一樣光是把一個(gè)demo跑起來就折騰了我一整天。但一旦摸清楚它的操作套路你會(huì)發(fā)現(xiàn)Chaste在處理細(xì)胞增殖、遷移、群體行為建模這些事上真的是利器——尤其當(dāng)你需要把細(xì)胞尺度上的個(gè)體行為跟組織尺度上的群體動(dòng)力學(xué)耦合在一起的時(shí)候它提供的現(xiàn)成模塊能幫你省掉大量造輪子的時(shí)間。這篇是“Chaste實(shí)戰(zhàn)”系列的第4篇重點(diǎn)講基本仿真操作從環(huán)境準(zhǔn)備、工程結(jié)構(gòu)、第一個(gè)仿真跑通到參數(shù)調(diào)優(yōu)、結(jié)果輸出和常見異常排查。適合已經(jīng)了解Chaste能做什么、但還沒真正上手跑過仿真的讀者。前面幾篇聊過框架背景和核心數(shù)學(xué)模型這篇我們直接動(dòng)手。1. 環(huán)境準(zhǔn)備先把“能跑起來”這件事搞定1.1 版本選擇和依賴安裝的坑Chaste對(duì)編譯環(huán)境的要求比較“復(fù)古”官方主推Linux環(huán)境下的安裝方式。我第一次在Mac上折騰光是OpenMPI的版本兼容問題就把我卡了半死。如果你不想在環(huán)境上耗太多時(shí)間直接上Ubuntu22.04 LTS實(shí)測(cè)很穩(wěn)裝好基礎(chǔ)依賴就能少走彎路。需要裝的東西說來也不復(fù)雜核心依賴就這幾類編譯器工具鏈g、gcc、make、cmake并行計(jì)算庫OpenMPI數(shù)值計(jì)算庫PETSc、HDF5、Boost、Eigen可視化輔助ParaView看結(jié)果用在Ubuntu上一條命令能解決大部分依賴sudo apt-get install g gcc make cmake libopenmpi-dev openmpi-bin petsc-dev libhdf5-dev libboost-all-dev libeigen3-dev這里有個(gè)關(guān)鍵細(xì)節(jié)PETSc版本和Chaste版本要配套。Chaste 3.4以上版本對(duì)PETSc 3.x支持比較友好但如果你用的PETSc版本太新或太老編譯時(shí)會(huì)出現(xiàn)莫名其妙的頭文件報(bào)錯(cuò)。我踩過最典型的一個(gè)PETSc 3.16某個(gè)版本改了Vec相關(guān)的接口定義導(dǎo)致Chaste編譯到一半直接掛掉。提示如果編譯報(bào)錯(cuò)指向petscversion.h大概率是PETSc版本不匹配。建議先查Chaste官方文檔里指定的PETSc版本范圍再?zèng)Q定裝哪個(gè)版本。別圖省事直接裝最新版穩(wěn)定壓倒一切。1.2 源碼獲取和編譯的正確姿勢(shì)Chaste的源碼托管在官方倉庫獲取方式很簡(jiǎn)單git clone https://chaste.cs.ox.ac.uk/git/chaste.git注意這個(gè)倉庫是包含了Chaste所有核心模塊的“全家桶”cell_based細(xì)胞級(jí)建模的核心、heart心臟電生理、lung肺、crypt腸道隱窩等。對(duì)于做細(xì)胞群體動(dòng)力學(xué)仿真的我們來說真正核心的模塊是cell_based其他的可以暫時(shí)不編譯。編譯前需要配置hostconfig文件。這是Chaste特有的配置方式告訴構(gòu)建系統(tǒng)“你機(jī)器上有哪些依賴、裝在哪個(gè)位置”。最簡(jiǎn)單的方式是復(fù)制一份默認(rèn)配置cd chaste cp hostconfig/settings.example hostconfig/settings.sh然后編輯settings.sh在里面指定PETSc的安裝路徑export CHASTE_BUILD_TYPERelease export PETSC_DIR/usr/lib/petsc export PETSC_ARCHlinux-gnu-c-opt我建議編譯時(shí)用Release模式Debug模式雖然方便調(diào)試但仿真速度會(huì)慢很多——細(xì)胞群體動(dòng)力學(xué)仿真通常要跑數(shù)千步時(shí)間迭代性能差距不是開玩笑的。編譯整個(gè)Chaste需要的時(shí)間比較長(zhǎng)就算機(jī)器好也得半小時(shí)往上。建議編譯時(shí)指定要編譯的模塊scons build_cell_based1這樣能省掉不少時(shí)間。第一次編譯成功那一刻你會(huì)覺得前面折騰環(huán)境的時(shí)間都值了——因?yàn)楹竺嫠蟹抡娌僮鞫冀⒃谶@個(gè)編譯鏈路上。2. 工程結(jié)構(gòu)Chaste項(xiàng)目應(yīng)該怎么組織2.1 一個(gè)仿真項(xiàng)目的標(biāo)準(zhǔn)目錄Chaste不像普通C工程那樣隨意組織源文件。它有自己的“規(guī)矩”理解這套規(guī)矩你才能在各模塊之間復(fù)用代碼。標(biāo)準(zhǔn)工程目錄長(zhǎng)這樣projects/ └── my_cell_project/ ├── src/ │ ├── cell_models/ │ ├── simulations/ │ └── utilities/ ├── test/ │ ├── simulations/ │ └── utilities/ └── apps/ └── src/src目錄放你自己寫的源碼模塊test目錄放可執(zhí)行的仿真測(cè)試——在Chaste的世界里“測(cè)試”這個(gè)詞的含義很廣單元測(cè)試是測(cè)試完整跑一個(gè)仿真也算“測(cè)試”。你寫一個(gè)仿真入口文件本質(zhì)上就是寫一個(gè)繼承自TestClass的測(cè)試類。這也是很多新手困惑的地方“我就是想跑個(gè)仿真干嘛還要寫成測(cè)試”Chaste的構(gòu)建系統(tǒng)設(shè)計(jì)成“一切以測(cè)試驅(qū)動(dòng)”仿真入口天然復(fù)用測(cè)試框架好處是能白嫖一堆斷言工具、臨時(shí)目錄管理和日志系統(tǒng)。理解了這一點(diǎn)你就不會(huì)覺得別扭了。2.2 我的第一個(gè)仿真文件從模板改起來與其從零手寫不如拿官方示例改。我最常用的是test目錄下的TestNumericalMethodsTutorial.hpp它的結(jié)構(gòu)清晰適合作為起點(diǎn)。仿真文件的核心結(jié)構(gòu)無非這幾塊#include cxxtest/TestSuite.h #include CellsGenerator.hpp #include HoneycombMeshGenerator.hpp #include NodeBasedCellPopulation.hpp #include OffLatticeSimulation.hpp #include FixedDurationGenerationBasedCellCycleModel.hpp class TestMyFirstSimulation : public CxxTest::TestSuite { public: void TestBasicSimulation() throw(Exception) { // 1. 建立網(wǎng)格 // 2. 生成細(xì)胞 // 3. 配置細(xì)胞周期模型 // 4. 設(shè)置仿真參數(shù) // 5. 運(yùn)行并輸出 } };這里的命名規(guī)范很Chaste類名以Test開頭公共方法以Test開頭。因?yàn)闃?gòu)建系統(tǒng)靠約定識(shí)別哪些是需要執(zhí)行的測(cè)試不遵守這套命名規(guī)則你的仿真根本不會(huì)被編譯進(jìn)可執(zhí)行文件里。3. 第一個(gè)仿真跑起來一個(gè)細(xì)胞群體增殖小實(shí)驗(yàn)3.1 從網(wǎng)格到細(xì)胞群體我習(xí)慣用一個(gè)最經(jīng)典的例子入門在二維規(guī)則網(wǎng)格上讓一圈細(xì)胞按固定周期增殖觀察群體擴(kuò)張。這類模型對(duì)應(yīng)的是腫瘤球生長(zhǎng)、傷口愈合等場(chǎng)景的基礎(chǔ)原型。建立網(wǎng)格和細(xì)胞群體代碼非常簡(jiǎn)潔// 建立一個(gè)12x12的蜂窩網(wǎng)格 HoneycombMeshGenerator generator(12, 12, 0); MutableMesh2,2* p_mesh generator.GetMesh(); // 生成細(xì)胞 std::vectorCellPtr cells; CellsGeneratorFixedDurationGenerationBasedCellCycleModel, 2 cell_generator; cell_generator.GenerateBasicRandom(cells, p_mesh-GetNumNodes()); // 組裝細(xì)胞群體 NodeBasedCellPopulation2 cell_population(*p_mesh, cells);這段代碼里HoneycombMeshGenerator是Chaste內(nèi)置的網(wǎng)格生成器蜂窩網(wǎng)格在二維空間里能實(shí)現(xiàn)最均勻的鄰居關(guān)系分布——這對(duì)模擬細(xì)胞接觸抑制和機(jī)械力傳導(dǎo)至關(guān)重要。比較有意思的是GenerateBasicRandom這個(gè)方法它給每個(gè)細(xì)胞隨機(jī)分配一個(gè)細(xì)胞周期狀態(tài)。如果不做這一步所有細(xì)胞的周期完全同步仿真出來的群體生長(zhǎng)曲線會(huì)呈現(xiàn)詭異的“階梯狀”而不是平滑的指數(shù)增長(zhǎng)。實(shí)際生物學(xué)里細(xì)胞群體本來就處于不同周期階段的混合狀態(tài)這個(gè)隨機(jī)化操作就是模擬這種異質(zhì)性。3.2 配置細(xì)胞周期和仿真參數(shù)細(xì)胞周期模型選FixedDurationGenerationBasedCellCycleModel意思是每個(gè)細(xì)胞的周期時(shí)間固定分裂一次后進(jìn)入下一輪周期。這個(gè)模型最簡(jiǎn)單適合作為基線。一個(gè)細(xì)胞周期的默認(rèn)配置是FixedDurationGenerationBasedCellCycleModel* p_model static_castFixedDurationGenerationBasedCellCycleModel*(cell.GetCellCycleModel()); p_model-SetCellCycleDuration(12.0); // 單位小時(shí)這里的時(shí)間單位需要注意Chaste中所有時(shí)間默認(rèn)單位是小時(shí)。細(xì)胞周期設(shè)為12小時(shí)意味著在仿真時(shí)間內(nèi)每個(gè)細(xì)胞大約每12小時(shí)分裂一次。然后設(shè)置仿真總時(shí)長(zhǎng)和時(shí)間步長(zhǎng)SimulationTime* p_time SimulationTime::Instance(); p_time-SetEndTimeAndTimeStep(120.0, 0.1); // 仿真120小時(shí)步長(zhǎng)0.1小時(shí)步長(zhǎng)選0.1小時(shí)是基于穩(wěn)定性和效率的折中。Chaste默認(rèn)采用Euler-Maruyama法求解隨機(jī)微分方程步長(zhǎng)太大會(huì)導(dǎo)致數(shù)值發(fā)散步長(zhǎng)太小則計(jì)算量暴漲。0.1小時(shí)在細(xì)胞周期為12小時(shí)的場(chǎng)景下相當(dāng)于每個(gè)細(xì)胞周期內(nèi)做120次積分精度完全夠用。3.3 運(yùn)行仿真并輸出結(jié)果組裝好細(xì)胞群體和仿真配置后創(chuàng)建一個(gè)OffLatticeSimulation對(duì)象來驅(qū)動(dòng)整個(gè)仿真流程OffLatticeSimulation2 simulator(cell_population); simulator.SetOutputDirectory(TestBasicSimulation); simulator.SetEndTime(120.0); simulator.SetSamplingTimestepMultiple(10); // 每10步記錄一次結(jié)果運(yùn)行仿真只需一行simulator.Solve();編譯和運(yùn)行這個(gè)仿真在projects/my_cell_project目錄下執(zhí)行scons my_cell_project/build/test/TestMyFirstSimulation ./build/test/TestMyFirstSimulation仿真運(yùn)行過程中會(huì)在testoutput/TestBasicSimulation目錄生成結(jié)果文件。跑完后你能看到類似這樣的輸出目錄結(jié)構(gòu)testoutput/TestBasicSimulation/ ├── results.viznodes ├── results.vizcelltypes ├── results.viznodeorigins ├── simulation_time.dat ├── results_from_time_0/ │ ├── results.pvd │ ├── results_0.vtu │ ├── results_10.vtu │ └── ... └── cell_population.dat重點(diǎn)說兩個(gè)文件results.pvd這是ParaView的時(shí)間序列索引文件。用ParaView打開這個(gè)文件能按時(shí)間軸播放群體擴(kuò)張動(dòng)畫直觀看到細(xì)胞增殖過程。cell_population.dat記錄每個(gè)時(shí)間點(diǎn)的細(xì)胞坐標(biāo)和狀態(tài)是后續(xù)做數(shù)值分析如計(jì)算群體半徑、細(xì)胞密度分布的原始數(shù)據(jù)來源。我第一次跑通這個(gè)例子的時(shí)候最震撼的是看ParaView里的動(dòng)畫細(xì)胞從一個(gè)點(diǎn)開始一圈圈向外擴(kuò)就像洋蔥的橫截面在生長(zhǎng)。那種“模型跑起來”的直觀反饋比任何數(shù)值指標(biāo)都更能增強(qiáng)你對(duì)仿真框架的信心。4. 參數(shù)調(diào)優(yōu)讓你的仿真具備生物學(xué)意義4.1 時(shí)間尺度與空間尺度參數(shù)背后的真實(shí)含義很多人調(diào)參數(shù)時(shí)只關(guān)心“能不能跑”不關(guān)心“跑出來的結(jié)果有沒有物理意義”。這其實(shí)是個(gè)大坑。Chaste里的參數(shù)不是隨便填的數(shù)字它在建模之初就有明確的生物學(xué)對(duì)應(yīng)。最重要的三個(gè)尺度參數(shù)參數(shù)默認(rèn)值生物學(xué)含義調(diào)參影響細(xì)胞周期時(shí)長(zhǎng)12小時(shí)一次有絲分裂所需時(shí)間決定群體擴(kuò)張速率時(shí)間步長(zhǎng)0.1小時(shí)數(shù)值積分的分辨率過大會(huì)震蕩過小會(huì)浪費(fèi)算力網(wǎng)格密度由網(wǎng)格生成器決定初始細(xì)胞間的間距影響細(xì)胞間機(jī)械力傳遞調(diào)細(xì)胞周期參數(shù)時(shí)我一般遵循一個(gè)原則先定時(shí)間尺度再定計(jì)算精度。以腫瘤生長(zhǎng)建模為例如果你關(guān)注的是“數(shù)天到數(shù)周”的宏觀生長(zhǎng)模式細(xì)胞周期設(shè)為12-24小時(shí)都是合理的但如果你關(guān)注的是“細(xì)胞分裂事件”本身的過程時(shí)間步長(zhǎng)就得縮到分鐘級(jí)別仿真成本會(huì)急劇上升。4.2 隨機(jī)性控制的細(xì)節(jié)讓結(jié)果可復(fù)現(xiàn)細(xì)胞群體動(dòng)力學(xué)仿真本質(zhì)上是一個(gè)隨機(jī)過程——細(xì)胞周期的初始化是隨機(jī)的力模型里可能引入隨機(jī)漲落。這帶來一個(gè)問題每次跑出來的結(jié)果都不一樣那怎么做可復(fù)現(xiàn)實(shí)驗(yàn)Chaste提供了隨機(jī)數(shù)生成器的種子控制接口RandomNumberGenerator* p_rng RandomNumberGenerator::Instance(); p_rng-Reseed(42);在仿真開頭顯式調(diào)用Reseed能把隨機(jī)數(shù)序列固定下來。這樣即便換了機(jī)器只要代碼和參數(shù)相同仿真結(jié)果也能完全重演。對(duì)于后續(xù)做參數(shù)敏感性分析、模型驗(yàn)證這是不可或缺的一步。提示多組對(duì)比仿真時(shí)不同組別要用不同的種子否則隨機(jī)差異可能被誤讀為參數(shù)效應(yīng)得出完全錯(cuò)誤的結(jié)論。4.3 仿真規(guī)模的合理選擇二維仿真的計(jì)算量通常還能接受但一上三維規(guī)模會(huì)呈數(shù)量級(jí)增長(zhǎng)。我?guī)湍抽_發(fā)者調(diào)試過一個(gè)三維腫瘤模型網(wǎng)格規(guī)模從2D的144個(gè)節(jié)點(diǎn)增加到3D的近萬個(gè)節(jié)點(diǎn)同樣的仿真時(shí)長(zhǎng)耗時(shí)從10分鐘漲到幾乎一天。規(guī)模選擇建議初步驗(yàn)證階段用2D小規(guī)模10x10網(wǎng)格跑通邏輯參數(shù)調(diào)優(yōu)階段2D中規(guī)模20x20跑多組對(duì)比正式實(shí)驗(yàn)階段3D規(guī)?？紤]并行計(jì)算Chaste的并行支持比我預(yù)期的成熟用MPI跑分布式仿真沒有想象中復(fù)雜。啟動(dòng)命令差不多是這樣mpirun -np 8 ./build/test/TestMyLargeSimulation不過要提醒一句不是所有模塊都能無縫并行。OffLatticeSimulation的并行效果尚可但涉及細(xì)胞機(jī)械力全局求解隱式時(shí)間積分的部分并行效率會(huì)打折扣。初次嘗試并行時(shí)建議先用小規(guī)模數(shù)據(jù)測(cè)試加速比再?zèng)Q定要不要上大規(guī)模。5. 結(jié)果可視化看懂你的仿真到底在干什么5.1 ParaView查看細(xì)胞運(yùn)動(dòng)的正確姿勢(shì)ParaView是Chaste結(jié)果可視化的標(biāo)配工具。但直接一堆vtu文件拖進(jìn)去默認(rèn)顯示方式可能讓你一頭霧水——滿屏的網(wǎng)格點(diǎn)根本看不出細(xì)胞動(dòng)態(tài)。我的習(xí)慣是這樣做打開results.pvd加載整個(gè)時(shí)間序列將渲染方式切到Surface默認(rèn)是Points視覺效果不好給細(xì)胞團(tuán)塊加一個(gè)Threshold濾波器過濾掉背景區(qū)域用Color by按細(xì)胞類型或細(xì)胞周期階段上色這樣你能清晰看到細(xì)胞群體從初始分布到增殖、擴(kuò)張、甚至競(jìng)爭(zhēng)排斥的動(dòng)態(tài)過程。如果想把結(jié)果導(dǎo)出成動(dòng)畫可以設(shè)置ParaView的Animation面板逐幀保存圖片再合并成視頻。我用這個(gè)方法給不少組會(huì)匯報(bào)做過仿真動(dòng)畫效果比純數(shù)據(jù)圖表直觀太多。5.2 數(shù)值指標(biāo)分析不能只看動(dòng)畫好看的可視化能幫你定性判斷“模型行為是否合理”但定量分析才是體現(xiàn)仿真價(jià)值的關(guān)鍵。我最常做的分析有三個(gè)群體細(xì)胞總數(shù)隨時(shí)間變化從cell_population.dat里提取每個(gè)時(shí)間節(jié)點(diǎn)的細(xì)胞數(shù)畫出生長(zhǎng)曲線看是否符合指數(shù)增長(zhǎng)預(yù)期。群體半徑變化計(jì)算細(xì)胞群體的均方根半徑分析擴(kuò)張速率。細(xì)胞密度分布統(tǒng)計(jì)單位面積內(nèi)細(xì)胞數(shù)量判斷是否有局部擁擠或空隙。Chaste原生提供了一些數(shù)據(jù)查詢接口比如unsigned num_cells cell_population.GetNumTissueCells();但復(fù)雜的自定義統(tǒng)計(jì)我習(xí)慣直接把cell_population.dat導(dǎo)出來用Python或R做后處理。把仿真引擎和數(shù)據(jù)分析解耦能避免頻繁重編譯項(xiàng)目的痛苦。5.3 跨時(shí)間片的細(xì)胞追蹤有時(shí)你關(guān)心某個(gè)特定細(xì)胞比如某個(gè)突變細(xì)胞克隆在群體中的位置變化。Chaste提供了基于節(jié)點(diǎn)索引的追蹤能力——每個(gè)細(xì)胞在cell_population.dat里都有固定ID只要記錄它在不同時(shí)間片的坐標(biāo)就能重建它的運(yùn)動(dòng)軌跡。實(shí)操思路是跑完仿真后用腳本解析results.viznodes文件按ID篩出目標(biāo)細(xì)胞繪制其坐標(biāo)隨時(shí)間變化的曲線。這種追蹤對(duì)研究細(xì)胞運(yùn)動(dòng)性和克隆擴(kuò)張模式特別有用。6. 從2D到3D復(fù)雜度躍遷的幾個(gè)注意點(diǎn)6.1 3D網(wǎng)格生成的差異從2D升級(jí)到3D第一個(gè)感受就是網(wǎng)格生成的邏輯變了。2D用的HoneycombMeshGenerator在3D世界里沒有直接對(duì)應(yīng)物需要改用HoneycombMeshGenerator結(jié)合分層堆疊的方式或者直接用CuboidalMeshGenerator生成規(guī)整的立方體網(wǎng)格。代碼大致是這樣CuboidalMeshGenerator generator; std::vectorunsigned dimensions; dimensions.push_back(5); // x方向節(jié)點(diǎn)數(shù) dimensions.push_back(5); // y方向節(jié)點(diǎn)數(shù) dimensions.push_back(5); // z方向節(jié)點(diǎn)數(shù) boost::shared_ptrMutableMesh3,3 p_mesh generator.GetMesh(dimensions);注意維度參數(shù)的單位是節(jié)點(diǎn)數(shù)而不是長(zhǎng)度實(shí)際物理尺寸需要乘上網(wǎng)格間距。很多人第一次轉(zhuǎn)3D時(shí)在這里算錯(cuò)導(dǎo)致細(xì)胞間距不符合預(yù)期。6.2 三維場(chǎng)景下的力模型差異在2D中細(xì)胞間的機(jī)械力傳導(dǎo)相對(duì)簡(jiǎn)單通常取近鄰細(xì)胞對(duì)被限制在平面內(nèi)。到了3D細(xì)胞間的力學(xué)相互作用多了一個(gè)自由度接觸面積、體積守恒約束、粘彈性模型的選擇都會(huì)顯著影響群體動(dòng)態(tài)。我在3D仿真中遇到過一個(gè)問題默認(rèn)的簡(jiǎn)單線性彈簧力模型在3D中表現(xiàn)極不穩(wěn)定細(xì)胞會(huì)明顯相互穿透。這是因?yàn)?D空間中每個(gè)細(xì)胞的鄰居數(shù)遠(yuǎn)多于2D疊加的力效應(yīng)被放大單純線性彈簧不夠用了。后來換用了Chaste內(nèi)置的GeneralisedLinearSpringForce并調(diào)整了彈簧剛度參數(shù)問題才緩解。這給我的教訓(xùn)是從2D遷移到3D不是把坐標(biāo)改成3維就完事力學(xué)參數(shù)通常需要重新標(biāo)定。7. 常見運(yùn)行錯(cuò)誤的排查思路7.1 編譯期錯(cuò)誤頭文件找不到怎么辦編譯錯(cuò)誤里最常見的就是頭文件缺失。Chaste的頭文件依賴關(guān)系比較隱晦——很多類是放在深層目錄里的光靠#include直覺經(jīng)常碰壁。我的排查流程是先看報(bào)錯(cuò)信息里fatal error: xxx.hpp: No such file or directory在Chaste源碼目錄下全局搜索這個(gè)頭文件find . -name xxx.hpp找到后修改include路徑或者引入正確的模塊依賴另一個(gè)典型的編譯錯(cuò)誤是類名拼寫或命名空間錯(cuò)誤。Chaste的類和命名空間組織非常嚴(yán)格CellBasedSimulationArchiver這類工具類往往藏在深層命名空間里用之前最好先查一下文檔。7.2 運(yùn)行期錯(cuò)誤仿真跑到一半崩潰運(yùn)行期最讓人頭疼的問題之一是“細(xì)胞數(shù)量爆炸”。有些細(xì)胞周期模型配置不當(dāng)會(huì)導(dǎo)致細(xì)胞指數(shù)級(jí)增殖最終內(nèi)存被塞滿。典型場(chǎng)景你設(shè)置了StochasticDurationGenerationBasedCellCycleModel但沒限制細(xì)胞周期的最短和最長(zhǎng)時(shí)長(zhǎng)結(jié)果某些細(xì)胞的周期短到離譜繁殖速度快得失控。解決思路p_model-SetMinCellCycleDuration(4.0); p_model-SetMaxCellCycleDuration(20.0);給細(xì)胞周期設(shè)置合理邊界是保證仿真穩(wěn)定基礎(chǔ)的常用手段。7.3 數(shù)值發(fā)散步長(zhǎng)過大惹的禍仿真結(jié)果出現(xiàn)“細(xì)胞飛走”或者坐標(biāo)變成NaN一般都能歸因到時(shí)間步長(zhǎng)過大。癥狀是前幾十步還好好的某一步突然所有細(xì)胞坐標(biāo)全變成亂數(shù)。遇到這種情況把時(shí)間步長(zhǎng)縮小一個(gè)數(shù)量級(jí)再試。實(shí)時(shí)上固定步長(zhǎng)0.1小時(shí)在2D模型里一般夠用但一旦引入強(qiáng)機(jī)械力或更密集的網(wǎng)格就可能需要縮小到0.01小時(shí)。注意除了縮小步長(zhǎng)還要檢查力模型的參數(shù)是否合理。有些力模型在參數(shù)設(shè)置不當(dāng)?shù)那闆r下即便步長(zhǎng)很小也會(huì)發(fā)散——那種情況通常是模型本身不穩(wěn)定建議換個(gè)力模型或在參數(shù)空間里搜索穩(wěn)定區(qū)間。8. 從仿真到科研落地的經(jīng)驗(yàn)之談8.1 仿真結(jié)果如何與實(shí)驗(yàn)數(shù)據(jù)對(duì)照跑通仿真只是第一步真正的挑戰(zhàn)是讓仿真結(jié)果跟實(shí)驗(yàn)數(shù)據(jù)對(duì)得上。我常用的對(duì)照思路是不直接比“絕對(duì)數(shù)值”而是比“無量綱趨勢(shì)”。比如實(shí)驗(yàn)里觀察到細(xì)胞群體擴(kuò)張速度隨時(shí)間遞減你的仿真模型是否也在后期出現(xiàn)擴(kuò)張速率放緩這種趨勢(shì)層面的對(duì)比允許模型參數(shù)在一定范圍內(nèi)波動(dòng)更關(guān)注行為模式的匹配比強(qiáng)行讓仿真曲線精確貼合實(shí)驗(yàn)曲線現(xiàn)實(shí)得多。做趨勢(shì)對(duì)照時(shí)可以在同一坐標(biāo)下把實(shí)驗(yàn)數(shù)據(jù)和仿真數(shù)據(jù)標(biāo)準(zhǔn)化到相同尺度然后疊加畫圖。如果兩條曲線的形態(tài)走向一致說明你的模型抓住了核心機(jī)制如果形態(tài)不一致大概率是模型設(shè)置遺漏了關(guān)鍵環(huán)節(jié)——比如沒有考慮細(xì)胞凋亡、忽視了接觸抑制、沒有加入營(yíng)養(yǎng)限制等。8.2 參數(shù)敏感性分析的必要性一套模型參數(shù)往往有十多個(gè)如果不做敏感性分析很難說你得出的結(jié)論是依賴某個(gè)具體參數(shù)取值還是模型本身的魯棒行為。做參數(shù)敏感性分析最常見的方式是單參數(shù)掃掠一個(gè)參數(shù)變化其余固定觀察關(guān)鍵輸出指標(biāo)的變化。跑多組對(duì)比仿真時(shí)務(wù)必每組用不同的隨機(jī)種子前面提過并把仿真統(tǒng)一納入自動(dòng)化腳本中。手動(dòng)一個(gè)個(gè)改參數(shù)跑仿真又慢又容易出錯(cuò)?，F(xiàn)在我的工作流是寫好一個(gè)仿真模板參數(shù)全部用命令行傳參控制再用Python腳本批量生成參數(shù)組、批量跑仿真、批量收集結(jié)果。一套流程跑下來幾十組仿真參數(shù)一夜之間就能完成。一個(gè)人工改參數(shù)、人工跑仿真的實(shí)驗(yàn)相比之下效率差太多了。8.3 Chaste還可以往哪個(gè)方向擴(kuò)展如果你已經(jīng)能熟練跑通基礎(chǔ)仿真接下來可以嘗試的方向很多。我比較看好的幾個(gè)方向多尺度耦合把細(xì)胞級(jí)模型與組織級(jí)反應(yīng)-擴(kuò)散方程耦合模擬營(yíng)養(yǎng)梯度對(duì)細(xì)胞增殖的影響。Chaste有現(xiàn)成的PDE模塊可以嘗試打通細(xì)胞模塊和PDE模塊的聯(lián)合仿真。突變克隆競(jìng)爭(zhēng)在細(xì)胞群體里引入具有不同增殖優(yōu)勢(shì)的突變細(xì)胞觀察克隆擴(kuò)張和競(jìng)爭(zhēng)動(dòng)態(tài)。這對(duì)腫瘤演化建模非常有用。接觸抑制的機(jī)理建模在細(xì)胞周期模型中引入基于鄰居密度的反饋控制讓細(xì)胞在擁擠環(huán)境中停止增殖。這是組織穩(wěn)態(tài)調(diào)控的關(guān)鍵機(jī)制。這些方向在Chaste里都有一定的官方支持但需要你自己拼裝組合。老實(shí)說上手階段每天都會(huì)遇到各種小問題但堅(jiān)持跑上兩三個(gè)完整項(xiàng)目你會(huì)對(duì)細(xì)胞群體動(dòng)力學(xué)仿真有真正的“手感”——那種參數(shù)一變就能預(yù)判結(jié)果走勢(shì)的直覺比任何教程都珍貴。