近年研究表明，某些藥物能改變基因調控網絡（GRNs）中的調控反應參數，促使細胞從病變狀態轉移到正常狀態。在此基礎上，有學者構建了基于吸引子和分岔理論的狀態調控框架。現有的控制策略中，控制信號規律都是自行擬定的，并且所使用的參數擾動法僅能求出單控制變量情況下狀態轉移所需的時間。因此，針對以非線性常微分方程組（ODEs）形式描述的復雜生物網絡，本文提出了一種基于動態優化算法改變網絡狀態的最優調控方法，在掌握生物網絡基本特性的基礎上建立了動態優化問題，并以簡單低維的三節點 GRN 和復雜高維的癌癥 GRN 為例，利用 MATLAB 計算得到單控制變量及多控制變量條件下的最優控制策略。本文所提方法旨在實現精確且快速的狀態調控，為試驗研究和實際治療提供參考依據。
引用本文： 揭豪, 袁美晨, 朱國柱, 洪偉榮. 基于動態優化算法的復雜生物網絡的狀態調控. 生物醫學工程學雜志, 2020, 37(1): 19-26. doi: 10.7507/1001-5515.201810008 復制
引言
基因表達過程涉及到 DNA、RNA 和蛋白質等多種物質，因而從分子生物學角度研究細胞行為，往往會描繪出一個由 DNA、RNA、蛋白質和小分子之間相互調控作用構成的復雜的基因調控網絡（gene regulatory networks，GRNs）。隨著基因微陣列芯片和轉錄組測序技術的發展，研究人員可以獲得大量高通量生物數據，因此 GRNs 的構建愈加普遍，其應用也愈發廣泛[1]。GRNs 的應用主要表現為：① 提供了生物分子相互作用的“藍圖”，能夠模擬系統的動態特性；② 根據特定表現型相對應的基因組表達情況對復雜疾病進行預測和診斷；③ 幫助識別致病基因和疾病發展路徑，為臨床給藥方案提供依據。
GRNs 的本質是在分子水平上描述系統動力學的抽象網絡模型，可以從控制理論出發研究其系統狀態的調控方法。對 GRNs 動力學的分析和干預是現代治療方式發展的新方向，其最終目標是開發干預活生物體功能的方法，以便將細胞從惡性狀態帶入良性狀態[2]。學界提出了各種控制策略使得 GRNs 動力學以期望的方式發生改變，生物學研究表明，通過引入改變細胞行為的因子或藥物，這些策略是可以實現的，有望成為未來基因治療的干預措施[2]。2006 年，Pal 等[3]針對概率布爾模型網絡，利用最優無限時域控制實現黑色素瘤 GRNs 的狀態轉變。2011 年，Liu 等[4]通過卡爾曼可控性秩條件來判斷線性常微分方程組（ordinary differential equations，ODEs）所述系統的可控性，并以 GRNs 為例進行說明。2016 年，Wang 等[5]提出了非線性 ODEs 系統的控制方法，通過參數擾動驅動系統實現狀態轉變。2017 年，Klickstein 等[6]從平衡控制能量和控制精度的角度，以新陳代謝等復雜網絡為研究對象，提出了線性 ODEs 系統的最優控制方法。2017 年，Jin 等[7]提出了一種非線性動態網絡域控制方法，實現了不同吸引域之間的轉變，并設計了控制方程的一般形式。
目前，常規用于描述 GRNs 的布爾網絡模型難以捕捉基因的連續動態特性，而線性 ODEs 模型簡化了非線性模型所特有的動態行為，因此非線性 ODEs 模型在模擬基因調控過程中的連續動態行為方面具有明顯的優勢。由于非線性 ODEs 的復雜性，目前對于該形式所表述的動態網絡的控制方法研究還不夠深入，大部分現有方法是基于物理學中的能量景觀地貌或者間接地從吸引域的角度出發。另一方面，控制變量隨時間變化的信號一般是自行制定的，而控制時長、控制效果等顯然是由控制信號規律決定的。為尋求實現精準、快速、高效控制的信號規律，本文以簡單低維的三節點 GRN 和復雜高維的癌癥 GRN 為例，針對非線性 ODEs 模型構建最優控制問題，提出了一種基于動態優化算法的復雜生物網絡狀態調控方法，能在已知吸引子網絡的基礎上直接獲得多控制變量同時作用的最優控制信號。
1 理論與方法
1.1 GRNs 建模
1.1.1 數學模型
目前，用于描述 GRNs 的模型主要有以下幾種：有向圖、布爾網絡、貝葉斯網絡、隨機方程以及 ODEs 等[8]。其中，ODEs 模型應用最為廣泛，非線性 ODEs 模型能很好地捕捉基因的連續動態特性，其一般形式如式（1）所示：
 |
其中，xi 表示 RNA、蛋白質或小分子的濃度，是隨時間 t 變化的非負變量；fi（?）表示非線性函數，除受系統狀態 x 自身影響外，還可以包含一些外部影響因素的輸入[9]；n 為涉及到的生物分子的個數。
實驗表明，GRNs 中各個分子間的相互作用表現為開關類方式，即調控函數是 S 形的，通常可用希爾函數（hill function）來描述[10]：以 
表示激勵作用，以 
表示抑制作用。對特定的 GRNs 來說，s 和 n 一般是實驗測定的常數系數[5]。有向圖由節點（頂點）和連接節點的有向邊組成，三節點 GRN 的有向圖形式如圖 1 所示，包含 3 個基因和 9 條作用邊，對應的非線性 ODE 模型如式（2）所示。該三節點 GRN 從實際生物網絡中抽象而來，具備復雜生物網絡的基本元素和特性，常用于系統動力學分析和控制方法的通用性研究，因此各變量未賦予具體的物理單位。
圖1
三節點 GRN
Figure1.
Three-node GRN
 |
其中，X1、X2、X3 代表三種基因或基因產物，x1、x2、x3 分別表示 X1、X2、X3 的濃度，
、
、
分別表示 X1、X2、X3 濃度的變化率，a1、b2、c3 分別表示 X1、X2、X3 對自身基因表達的激勵作用強度，a2、a3 分別表示 X1 對 X2、X3 基因表達的抑制作用強度，b1、b3 分別表示 X2 對 X1、X3 基因表達的抑制作用強度，c1、c2 分別表示 X3 對 X1、X2 基因表達的抑制作用強度，n 和 s 均為希爾函數中的參數，k 表示自降解速率。根據文獻[4]，參數取值為 n = 4、s = 0.5、k = 1.0，激勵作用強度 a1、b2、c3 的初始標準值為 1.0，抑制作用強度 a2、a3、b1、b3、c1、c2 的初始標準值為 0.1。
1.1.2 吸引子
基于對 GRNs 的分析，一些生物網絡擁有多個代表正常、病變和凋亡狀態的穩定穩態，或稱吸引子[11]。吸引子的本質是微分動力系統的平衡點，平衡點在微分方程定性理論中也稱為“奇點”。在平衡點上，系統狀態可以維持穩定，即平衡點 x*滿足 
。吸引子是具有一定特性的平衡點，每個吸引子都對應一個鄰域，稱為吸引域，只要系統狀態位于吸引子的吸引域內，系統就會向該吸引子收斂。三節點 GRN 的吸引子在相空間中的分布如圖 2 所示，可以看出該 GRN 存在 8 個吸引子，依次記為 A（0.105 0，1.167 0，0.105 0）、B（0.240 8，0.240 8，0.240 8）、C（1.167 0，0.105 0，0.105 0）、D（1.057 1，1.057 1，0.009 5）、E（0.009 5，1.057 1，1.057 1）、F（0.105 0，0.105 0，1.167 0）、G（1.057 1，0.009 5，1.057 1）、H（0.940 2，0.940 2，0.940 2）。
圖2
三節點 GRN 的吸引子分布
Figure2.
Attractor distribution of three-node GRN
1.1.3 狀態調控
文獻[7]研究表明，改變 ODEs 中的調控反應參數，可以改變相應吸引子的位置及吸引域，從而使原狀態落入過渡吸引子的吸引域，然后向過渡吸引子轉移；一定時間后系統狀態落入期望吸引子的吸引域，撤除參數的改變，系統將繼續向期望吸引子轉移，最終實現由原狀態向期望狀態的轉變。這些參數的生物學意義為，基因產物對自身或其他基因表達水平的激勵或抑制作用強度，可通過一組誘導或抑制的藥物進行調整，從而具有治療功能[12]。
實現狀態轉移的基本控制原理為分岔理論。對于含參數的系統，當參數變動并通過某些臨界值時，系統的定性狀態（例如：平衡狀態、周期運動的數目和穩定性等）會發生突然變化，這種變化稱為分岔，是一類常見的重要非線性現象[13]。通過對 GRNs 的模擬計算與參數擾動，可以識別出控制狀態轉移的參數及其臨界值（即分岔點）。三節點 GRN 在初始吸引子為 A、參數 b1 變化的情況下，系統穩定后狀態變量 x 的變化規律如圖 3 所示，形象地描繪了分岔現象，臨界值為 b1 = 2.756 1。在臨界值上，x 的值發生了突變，而在臨界值兩側，x 變化平穩。
圖3
狀態變量穩定值隨控制變量變化的情況
Figure3.
Variation of steady state variables when control variables varying
1.1.4 吸引子網絡
將三節點 GRN 系統的每一個吸引子（A～H）分別定義為一個節點，以有向線段表示能控制實現的狀態轉移路徑，可以繪制出吸引子網絡圖[5]，如圖 4 所示。該圖是一張連通圖，說明可以實現從任意初始狀態到任意期望狀態的轉移，該 GRN 是完全可控的。
圖4
三節點 GRN 的吸引子網絡
Figure4.
Attractor network of three-node GRN
1.2 動態優化算法
1.2.1 動態優化問題
最優化問題的數學表述一般包含目標函數和以等式或不等式表達的約束條件。根據變量是否連續、各函數是否線性、約束模型是否與時間相關等條件，可對優化問題進行分類[14]。動態優化問題是針對連續時間動態模型的最優控制問題，其數學模型如式（3）～式（8）所示[15]。
 |
 |
 |
 |
 |
 |
其中 x 為狀態變量，u 為控制變量，
為狀態變量對時間 t 的一階導數，
為終值時刻，?（?）為目標函數，F（?）為動態系統的非線性 ODEs 模型，xmax、xmin 分別為狀態變量上、下邊界，umax、umin 分別為控制變量上、下邊界，x0 為狀態變量初始時刻值。
動態優化問題的求解方法可以分為間接法、迭代動態規劃方法和直接法三種[15]。只有直接法適合于求解大規模的約束動態優化問題，近些年動態優化算法的研究也主要集中于此。直接法主要有序貫算法、聯立算法和擬序貫算法，其基本思想都是通過一定程度的離散把原來的連續時間問題轉換為離散問題，即將有限時間域分割為有限個時間單元，在時間單元端點上施加約束條件，使動態優化問題轉換為非線性規劃（nonlinear programming，NLP）問題，通過 NLP 算法迭代求解。
1.2.2 擬序貫算法
擬序貫算法結合了序貫算法和聯立算法的優點，適合求解帶路徑約束的大規模強非線性動態優化問題，分為模擬計算層和優化計算層[15-16]。將有限時間域分割為有限個時間單元，控制變量表示為分段常數函數。在模擬計算層中，狀態變量用正交配置法[17]離散，從而將非線性 ODEs 離散為非線性代數方程組，可以通過牛頓迭代法求解，得到離散點上的狀態變量值。與此同時，根據離散點上狀態變量關于控制變量的偏導數公式可計算出靈敏度矩陣供優化計算層獲取目標函數對于控制變量的梯度信息。這樣控制變量分段函數值成為了優化問題的決策變量，約束條件中不再含有 ODEs 等式約束，狀態變量僅以路徑約束的形式在不等式約束中出現，原問題簡化為 NLP 問題，可以直接應用標準的序列二次規劃（sequential quadratic programming，SQP）算法求解[18]。
1.2.3 目標函數
目標函數表征期望達到的最優性能指標，而本文需要達到兩個目標——實現狀態轉移和轉移時間最小化，并且第一個目標的優先級更高，體現在目標函數的權重系數中。動態過程被分割成了有限個時間單元，要實現轉移時間最小化，需將離散時間單元長度作為決策變量一并寫入目標函數中，故目標函數的形式選取如式（9）所示[19]：
 |
其中，?（?）為目標函數，
、
、
為終值時刻 
時的狀態變量值，
、
、
為期望吸引子的狀態變量值，w1、w2、w3 為權重系數，NL 為離散時間單元數量，
為第 i 個離散時間單元的長度。
1.2.4 程序實現
從吸引子網絡中選取一條狀態轉移路徑，可以確定初始吸引子和目標吸引子 
、
、
。將該路徑對應的調控反應參數作為控制變量，根據臨界值選取合適的邊界條件，其余參數取其標準值，得到 ODEs 約束，形成完整的動態優化問題。通過數值計算軟件 MATLAB R2018b（MathWorks Inc.，美國）實現擬序貫算法，計算得到最優控制策略。
2 結果與討論
2.1 三節點 GRN 優化調控
選取狀態轉移路徑 A→D，則 b1 臨界值為 2.756 1、c1 臨界值為 0.188 9，選擇合適的變量上下限，實際問題中可由療效和藥物用量限制得出。b1 和 c1 均作為控制變量的動態優化問題如式（10）～式（15）所示。
 |
 |
 |
 |
 |
 |
實現三節點 GRN 狀態調控 A→D 的最優控制策略如圖 5 所示。b1 單獨作用時，狀態轉移所需最小時間 T = 17.97，但最后一段時間控制變量 b1 回到了初始值 0.1，說明這段時間不需要在外部施加控制作用，故實際最優控制時間為 T = 16.79，且最優控制策略形同矩形波。c1 單獨作用時，實際最優控制時間為 T = 4.28，表明該參數的控制效率優于 b1。b1、c1 同時作用時，狀態轉移所需最小時間縮短為 T = 4.91，而實際最優控制時間縮短為 T = 2.82，相較于單控制變量情形，控制效率有顯著的提升。
圖5
狀態轉移 A→D 的最優控制策略
Figure5.
Optimal control strategies of state transition A→D
2.2 癌癥 GRN 模型
2014 年，Li 等[20]從一些癌變相關的關鍵標志基因出發尋找有相互作用的基因，基于以往的實驗結果構建了典型癌癥 GRN，如圖 6 所示。該 GRN 有 32 個節點、66 條激勵作用邊和 45 條抑制作用邊，其中節點表示基因，藍色箭頭線表示激勵作用，紅色圓頭線表示抑制作用，八邊形節點和實線線段分別表示控制癌變狀態向正常狀態轉移的關鍵基因和關鍵作用邊。根據各基因的生物學作用將節點基因分為四類：綠色代表抑癌標志基因，黃色代表凋亡標志基因，橙紅色代表癌變標志基因，紫色代表其他基因。
圖6
癌癥 GRN
Figure6.
Cancer GRN
圖 6 中各基因分別為：磷酸酯酶與張力蛋白同源物（phosphatase and tensin homolog，PTEN）、B 淋巴細胞瘤-2 基因（B-cell lymphoma 2，BCL2）、BCL2 相關死亡蛋白（BCL2 associated agonist of cell death，BAD）、血管內皮生長因子（vascular endothelial growth factor，VEGF）、BCL2 相關 X 蛋白（BCL2 associated X protein，BAX）、周期蛋白依賴性激酶 1/2/4（cyclin-dependent kinase 1/2/4，CDK1/CDK2/CDK4）、肝細胞生長因子（hepatocyte growth factor，HGF）、表皮生長因子受體（epidermal growth factor receptor，EGFR）、轉化生長因子 α/β（transforming growth factor α/β，TGFα/TGFβ）、腫瘤壞死因子（tumor necrosis factor α，TNFα）、E2F1、視網膜母細胞瘤蛋白（retinoblastoma protein，RB）、共濟失調毛細血管擴張突變基因（ataxia telangiectasia-mutated gene，ATM）、野生型 p53 誘導的磷酸酶（wild-type p53-induced phosphatase 1，Wip1）、生長素反應因子（auxin response factor，ARF）、人端粒酶逆轉錄酶（telomerase reverse transcriptase in humans，hTERT）、雄性激素受體（androgen receptor，AR）、活化 B 細胞的 κ 輕鏈增強核因子（nuclear factor kappa-light-chain-enhancer of activated B cells，NFκB）、缺氧誘導因子（hypoxia-inducible factors 1，HIF1）；鼠雙微體基因 2/X（murine double minute 2/X，MDM2/MDMX）。
如圖 6 所示的癌癥 GRN 有三個吸引子，對應于細胞的三種生物學狀態：正常狀態 N、癌變狀態 C 和凋亡狀態 A。其吸引子網絡如圖 7 所示，三種狀態可以相互轉移。將該網絡所包含的 ATM、P53、P21、PTEN、CDH1、RB、ARF、AR、MYC、AKT、EGFR、VEGF、HGF、H1F1、hTERT、MDM2、CDK2、CDK4、CDK1、E2F1、Caspase、BAX、BAD、BCL2、NFκB、RAS、TGFα、TNFα、TGFβ、Wee1、MDMX、Wip1 依次編為 1～32 號基因，具體的非線性 ODEs 模型、參數取值及各吸引子的數值詳見文獻[7, 20]。通過對模型方程的模擬計算與參數擾動，可以找到關鍵作用邊中能使細胞從癌變狀態轉移到正常狀態的可行路徑共有 4 條，以┤表示抑制、→表示激勵，則分別為 R1（CDK2 ┤RB）、R10（P53 ┤ VEGF）、A9（HGF→AKT）、A11（AKT→VEGF）。4 條可行路徑上的調控反應參數符號及對應的臨界值如表 1 所示，以路徑 R1 為例，R1 表示 CDK2 對 RB 的抑制作用，CDK2 為 17 號基因，RB 為 6 號基因，故參數符號記為 
，標準值為無外部影響時的數值。
表1
可行路徑臨界值
Table1.
Critical values of feasible paths
圖7
癌癥 GRN 的吸引子網絡
Figure7.
The attractor network of cancer GRN
2.3 癌癥 GRN 優化調控
生物學研究表明，6 號基因 RB（
）、10 號基因 AKT（
）和 21 號基因 Caspase（
）是表征細胞狀態的三個重要基因，在不同的細胞狀態下表達水平存在很大的差異，可作為癌癥 GRN 的狀態觀測項寫入目標函數中。將路徑 R1 和 A9 對應的調控反應參數同時作為控制變量，目標為控制 GRN 從癌變狀態 C 轉移至正常狀態 N 并且滿足控制時間最小化，據此建立如式（16）～式（21）所示的動態優化問題。另外兩條路徑 R10 與 A11 實現狀態轉移的效果較差，體現在其所需控制時間很長，同時考慮到程序計算效率，不將其納入最優控制策略計算中。
 |
 |
 |
 |
 |
 |
實現癌癥 GRN 狀態調控 C→N 的最優控制策略如圖 8 表示。路徑 R1 對應參數 
單獨作用時，最優控制時間 T = 16.93，因為分岔后的過渡吸引子與表示正常狀態的吸引子很接近，故沒有返回標準值的變化趨勢。路徑 A9 對應參數 
單獨作用時，實際最優控制時間 T = 44.28，相對于路徑 R1 來說所需控制時間較長。
和 
同時作為控制變量作用于網絡時，最優控制時間 T = 11.07，其中對路徑 A9 的實際最優控制時間 T = 4.20。與單控制變量相比，多控制變量共同作用時，控制效率有顯著的提升。
圖8
狀態轉移 C→N 的最優控制策略
Figure8.
Optimal control strategies of state transition C→N
3 結論
對復雜生物網絡建立非線性 ODEs 模型，通過模擬計算和參數擾動得到系統的吸引子網絡，在此基礎上建立動態優化問題，利用動態優化算法求解，最終得到實現網絡狀態調控的最優控制策略。該方法具有普適性，既適用于簡單低維的生物網絡，也適用于與復雜高維的生物網絡；既適用于求解單控制變量問題，也適用于求解多控制變量問題。本研究所列舉的計算結果多為矩形脈沖形式，這可能與網絡本身特性有關，對于其他網絡不能得出矩形脈沖信號為最優策略的一般結論，但具有參考價值，必要時可以建立相應動態優化問題求解得到最優控制策略。另外，由計算結果可知，多控制變量共同作用時的控制效率往往高于單控制變量。
本文基于動態優化算法，提出了尋找實現復雜生物網絡狀態轉移最優控制策略的一般化方法，存在一定的應用前景。非線性 ODEs 模型均為實際生物網絡經過量化后建立的，狀態變量 x 與物質的實際濃度存在確定的函數關系。藥物的促進和抑制作用強度也經過量化，故控制變量 u 與藥物使用的劑量和周期存在一一對應的關系，即建立起了數學模型與實際情況的聯系。
大部分醫學實驗和實際治療中需要制定藥物的用量和使用周期，往往依照經驗而定，亟需定量化的參考。將本文的方法結論應用于實踐，可以指導以上參數的定量選擇，實驗中減小試錯次數、提高效率，實踐中提升療效。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
基因表達過程涉及到 DNA、RNA 和蛋白質等多種物質，因而從分子生物學角度研究細胞行為，往往會描繪出一個由 DNA、RNA、蛋白質和小分子之間相互調控作用構成的復雜的基因調控網絡（gene regulatory networks，GRNs）。隨著基因微陣列芯片和轉錄組測序技術的發展，研究人員可以獲得大量高通量生物數據，因此 GRNs 的構建愈加普遍，其應用也愈發廣泛[1]。GRNs 的應用主要表現為：① 提供了生物分子相互作用的“藍圖”，能夠模擬系統的動態特性；② 根據特定表現型相對應的基因組表達情況對復雜疾病進行預測和診斷；③ 幫助識別致病基因和疾病發展路徑，為臨床給藥方案提供依據。
GRNs 的本質是在分子水平上描述系統動力學的抽象網絡模型，可以從控制理論出發研究其系統狀態的調控方法。對 GRNs 動力學的分析和干預是現代治療方式發展的新方向，其最終目標是開發干預活生物體功能的方法，以便將細胞從惡性狀態帶入良性狀態[2]。學界提出了各種控制策略使得 GRNs 動力學以期望的方式發生改變，生物學研究表明，通過引入改變細胞行為的因子或藥物，這些策略是可以實現的，有望成為未來基因治療的干預措施[2]。2006 年，Pal 等[3]針對概率布爾模型網絡，利用最優無限時域控制實現黑色素瘤 GRNs 的狀態轉變。2011 年，Liu 等[4]通過卡爾曼可控性秩條件來判斷線性常微分方程組（ordinary differential equations，ODEs）所述系統的可控性，并以 GRNs 為例進行說明。2016 年，Wang 等[5]提出了非線性 ODEs 系統的控制方法，通過參數擾動驅動系統實現狀態轉變。2017 年，Klickstein 等[6]從平衡控制能量和控制精度的角度，以新陳代謝等復雜網絡為研究對象，提出了線性 ODEs 系統的最優控制方法。2017 年，Jin 等[7]提出了一種非線性動態網絡域控制方法，實現了不同吸引域之間的轉變，并設計了控制方程的一般形式。
目前，常規用于描述 GRNs 的布爾網絡模型難以捕捉基因的連續動態特性，而線性 ODEs 模型簡化了非線性模型所特有的動態行為，因此非線性 ODEs 模型在模擬基因調控過程中的連續動態行為方面具有明顯的優勢。由于非線性 ODEs 的復雜性，目前對于該形式所表述的動態網絡的控制方法研究還不夠深入，大部分現有方法是基于物理學中的能量景觀地貌或者間接地從吸引域的角度出發。另一方面，控制變量隨時間變化的信號一般是自行制定的，而控制時長、控制效果等顯然是由控制信號規律決定的。為尋求實現精準、快速、高效控制的信號規律，本文以簡單低維的三節點 GRN 和復雜高維的癌癥 GRN 為例，針對非線性 ODEs 模型構建最優控制問題，提出了一種基于動態優化算法的復雜生物網絡狀態調控方法，能在已知吸引子網絡的基礎上直接獲得多控制變量同時作用的最優控制信號。
1 理論與方法
1.1 GRNs 建模
1.1.1 數學模型
目前，用于描述 GRNs 的模型主要有以下幾種：有向圖、布爾網絡、貝葉斯網絡、隨機方程以及 ODEs 等[8]。其中，ODEs 模型應用最為廣泛，非線性 ODEs 模型能很好地捕捉基因的連續動態特性，其一般形式如式（1）所示：
|
其中，xi 表示 RNA、蛋白質或小分子的濃度，是隨時間 t 變化的非負變量；fi（?）表示非線性函數，除受系統狀態 x 自身影響外，還可以包含一些外部影響因素的輸入[9]；n 為涉及到的生物分子的個數。
實驗表明，GRNs 中各個分子間的相互作用表現為開關類方式，即調控函數是 S 形的，通常可用希爾函數（hill function）來描述[10]：以 表示激勵作用，以 表示抑制作用。對特定的 GRNs 來說，s 和 n 一般是實驗測定的常數系數[5]。有向圖由節點（頂點）和連接節點的有向邊組成，三節點 GRN 的有向圖形式如圖 1 所示，包含 3 個基因和 9 條作用邊，對應的非線性 ODE 模型如式（2）所示。該三節點 GRN 從實際生物網絡中抽象而來，具備復雜生物網絡的基本元素和特性，常用于系統動力學分析和控制方法的通用性研究，因此各變量未賦予具體的物理單位。
圖1
三節點 GRN
Figure1.
Three-node GRN
|
其中，X1、X2、X3 代表三種基因或基因產物，x1、x2、x3 分別表示 X1、X2、X3 的濃度，、、 分別表示 X1、X2、X3 濃度的變化率，a1、b2、c3 分別表示 X1、X2、X3 對自身基因表達的激勵作用強度，a2、a3 分別表示 X1 對 X2、X3 基因表達的抑制作用強度，b1、b3 分別表示 X2 對 X1、X3 基因表達的抑制作用強度，c1、c2 分別表示 X3 對 X1、X2 基因表達的抑制作用強度，n 和 s 均為希爾函數中的參數，k 表示自降解速率。根據文獻[4]，參數取值為 n = 4、s = 0.5、k = 1.0，激勵作用強度 a1、b2、c3 的初始標準值為 1.0，抑制作用強度 a2、a3、b1、b3、c1、c2 的初始標準值為 0.1。
1.1.2 吸引子
基于對 GRNs 的分析，一些生物網絡擁有多個代表正常、病變和凋亡狀態的穩定穩態，或稱吸引子[11]。吸引子的本質是微分動力系統的平衡點，平衡點在微分方程定性理論中也稱為“奇點”。在平衡點上，系統狀態可以維持穩定，即平衡點 x*滿足 。吸引子是具有一定特性的平衡點，每個吸引子都對應一個鄰域，稱為吸引域，只要系統狀態位于吸引子的吸引域內，系統就會向該吸引子收斂。三節點 GRN 的吸引子在相空間中的分布如圖 2 所示，可以看出該 GRN 存在 8 個吸引子，依次記為 A（0.105 0，1.167 0，0.105 0）、B（0.240 8，0.240 8，0.240 8）、C（1.167 0，0.105 0，0.105 0）、D（1.057 1，1.057 1，0.009 5）、E（0.009 5，1.057 1，1.057 1）、F（0.105 0，0.105 0，1.167 0）、G（1.057 1，0.009 5，1.057 1）、H（0.940 2，0.940 2，0.940 2）。
圖2
三節點 GRN 的吸引子分布
Figure2.
Attractor distribution of three-node GRN
1.1.3 狀態調控
文獻[7]研究表明，改變 ODEs 中的調控反應參數，可以改變相應吸引子的位置及吸引域，從而使原狀態落入過渡吸引子的吸引域，然后向過渡吸引子轉移；一定時間后系統狀態落入期望吸引子的吸引域，撤除參數的改變，系統將繼續向期望吸引子轉移，最終實現由原狀態向期望狀態的轉變。這些參數的生物學意義為，基因產物對自身或其他基因表達水平的激勵或抑制作用強度，可通過一組誘導或抑制的藥物進行調整，從而具有治療功能[12]。
實現狀態轉移的基本控制原理為分岔理論。對于含參數的系統，當參數變動并通過某些臨界值時，系統的定性狀態（例如：平衡狀態、周期運動的數目和穩定性等）會發生突然變化，這種變化稱為分岔，是一類常見的重要非線性現象[13]。通過對 GRNs 的模擬計算與參數擾動，可以識別出控制狀態轉移的參數及其臨界值（即分岔點）。三節點 GRN 在初始吸引子為 A、參數 b1 變化的情況下，系統穩定后狀態變量 x 的變化規律如圖 3 所示，形象地描繪了分岔現象，臨界值為 b1 = 2.756 1。在臨界值上，x 的值發生了突變，而在臨界值兩側，x 變化平穩。
圖3
狀態變量穩定值隨控制變量變化的情況
Figure3.
Variation of steady state variables when control variables varying
1.1.4 吸引子網絡
將三節點 GRN 系統的每一個吸引子（A～H）分別定義為一個節點，以有向線段表示能控制實現的狀態轉移路徑，可以繪制出吸引子網絡圖[5]，如圖 4 所示。該圖是一張連通圖，說明可以實現從任意初始狀態到任意期望狀態的轉移，該 GRN 是完全可控的。
圖4
三節點 GRN 的吸引子網絡
Figure4.
Attractor network of three-node GRN
1.2 動態優化算法
1.2.1 動態優化問題
最優化問題的數學表述一般包含目標函數和以等式或不等式表達的約束條件。根據變量是否連續、各函數是否線性、約束模型是否與時間相關等條件，可對優化問題進行分類[14]。動態優化問題是針對連續時間動態模型的最優控制問題，其數學模型如式（3）～式（8）所示[15]。
|
|
|
|
|
|
其中 x 為狀態變量，u 為控制變量， 為狀態變量對時間 t 的一階導數， 為終值時刻，?（?）為目標函數，F（?）為動態系統的非線性 ODEs 模型，xmax、xmin 分別為狀態變量上、下邊界，umax、umin 分別為控制變量上、下邊界，x0 為狀態變量初始時刻值。
動態優化問題的求解方法可以分為間接法、迭代動態規劃方法和直接法三種[15]。只有直接法適合于求解大規模的約束動態優化問題，近些年動態優化算法的研究也主要集中于此。直接法主要有序貫算法、聯立算法和擬序貫算法，其基本思想都是通過一定程度的離散把原來的連續時間問題轉換為離散問題，即將有限時間域分割為有限個時間單元，在時間單元端點上施加約束條件，使動態優化問題轉換為非線性規劃（nonlinear programming，NLP）問題，通過 NLP 算法迭代求解。
1.2.2 擬序貫算法
擬序貫算法結合了序貫算法和聯立算法的優點，適合求解帶路徑約束的大規模強非線性動態優化問題，分為模擬計算層和優化計算層[15-16]。將有限時間域分割為有限個時間單元，控制變量表示為分段常數函數。在模擬計算層中，狀態變量用正交配置法[17]離散，從而將非線性 ODEs 離散為非線性代數方程組，可以通過牛頓迭代法求解，得到離散點上的狀態變量值。與此同時，根據離散點上狀態變量關于控制變量的偏導數公式可計算出靈敏度矩陣供優化計算層獲取目標函數對于控制變量的梯度信息。這樣控制變量分段函數值成為了優化問題的決策變量，約束條件中不再含有 ODEs 等式約束，狀態變量僅以路徑約束的形式在不等式約束中出現，原問題簡化為 NLP 問題，可以直接應用標準的序列二次規劃（sequential quadratic programming，SQP）算法求解[18]。
1.2.3 目標函數
目標函數表征期望達到的最優性能指標，而本文需要達到兩個目標——實現狀態轉移和轉移時間最小化，并且第一個目標的優先級更高，體現在目標函數的權重系數中。動態過程被分割成了有限個時間單元，要實現轉移時間最小化，需將離散時間單元長度作為決策變量一并寫入目標函數中，故目標函數的形式選取如式（9）所示[19]：
|
其中，?（?）為目標函數，、、 為終值時刻 時的狀態變量值，、、 為期望吸引子的狀態變量值，w1、w2、w3 為權重系數，NL 為離散時間單元數量， 為第 i 個離散時間單元的長度。
1.2.4 程序實現
從吸引子網絡中選取一條狀態轉移路徑，可以確定初始吸引子和目標吸引子 、、。將該路徑對應的調控反應參數作為控制變量，根據臨界值選取合適的邊界條件，其余參數取其標準值，得到 ODEs 約束，形成完整的動態優化問題。通過數值計算軟件 MATLAB R2018b（MathWorks Inc.，美國）實現擬序貫算法，計算得到最優控制策略。
2 結果與討論
2.1 三節點 GRN 優化調控
選取狀態轉移路徑 A→D，則 b1 臨界值為 2.756 1、c1 臨界值為 0.188 9，選擇合適的變量上下限，實際問題中可由療效和藥物用量限制得出。b1 和 c1 均作為控制變量的動態優化問題如式（10）～式（15）所示。
|
|
|
|
|
|
實現三節點 GRN 狀態調控 A→D 的最優控制策略如圖 5 所示。b1 單獨作用時，狀態轉移所需最小時間 T = 17.97，但最后一段時間控制變量 b1 回到了初始值 0.1，說明這段時間不需要在外部施加控制作用，故實際最優控制時間為 T = 16.79，且最優控制策略形同矩形波。c1 單獨作用時，實際最優控制時間為 T = 4.28，表明該參數的控制效率優于 b1。b1、c1 同時作用時，狀態轉移所需最小時間縮短為 T = 4.91，而實際最優控制時間縮短為 T = 2.82，相較于單控制變量情形，控制效率有顯著的提升。
圖5
狀態轉移 A→D 的最優控制策略
Figure5.
Optimal control strategies of state transition A→D
2.2 癌癥 GRN 模型
2014 年，Li 等[20]從一些癌變相關的關鍵標志基因出發尋找有相互作用的基因，基于以往的實驗結果構建了典型癌癥 GRN，如圖 6 所示。該 GRN 有 32 個節點、66 條激勵作用邊和 45 條抑制作用邊，其中節點表示基因，藍色箭頭線表示激勵作用，紅色圓頭線表示抑制作用，八邊形節點和實線線段分別表示控制癌變狀態向正常狀態轉移的關鍵基因和關鍵作用邊。根據各基因的生物學作用將節點基因分為四類：綠色代表抑癌標志基因，黃色代表凋亡標志基因，橙紅色代表癌變標志基因，紫色代表其他基因。
圖6
癌癥 GRN
Figure6.
Cancer GRN
圖 6 中各基因分別為：磷酸酯酶與張力蛋白同源物（phosphatase and tensin homolog，PTEN）、B 淋巴細胞瘤-2 基因（B-cell lymphoma 2，BCL2）、BCL2 相關死亡蛋白（BCL2 associated agonist of cell death，BAD）、血管內皮生長因子（vascular endothelial growth factor，VEGF）、BCL2 相關 X 蛋白（BCL2 associated X protein，BAX）、周期蛋白依賴性激酶 1/2/4（cyclin-dependent kinase 1/2/4，CDK1/CDK2/CDK4）、肝細胞生長因子（hepatocyte growth factor，HGF）、表皮生長因子受體（epidermal growth factor receptor，EGFR）、轉化生長因子 α/β（transforming growth factor α/β，TGFα/TGFβ）、腫瘤壞死因子（tumor necrosis factor α，TNFα）、E2F1、視網膜母細胞瘤蛋白（retinoblastoma protein，RB）、共濟失調毛細血管擴張突變基因（ataxia telangiectasia-mutated gene，ATM）、野生型 p53 誘導的磷酸酶（wild-type p53-induced phosphatase 1，Wip1）、生長素反應因子（auxin response factor，ARF）、人端粒酶逆轉錄酶（telomerase reverse transcriptase in humans，hTERT）、雄性激素受體（androgen receptor，AR）、活化 B 細胞的 κ 輕鏈增強核因子（nuclear factor kappa-light-chain-enhancer of activated B cells，NFκB）、缺氧誘導因子（hypoxia-inducible factors 1，HIF1）；鼠雙微體基因 2/X（murine double minute 2/X，MDM2/MDMX）。
如圖 6 所示的癌癥 GRN 有三個吸引子，對應于細胞的三種生物學狀態：正常狀態 N、癌變狀態 C 和凋亡狀態 A。其吸引子網絡如圖 7 所示，三種狀態可以相互轉移。將該網絡所包含的 ATM、P53、P21、PTEN、CDH1、RB、ARF、AR、MYC、AKT、EGFR、VEGF、HGF、H1F1、hTERT、MDM2、CDK2、CDK4、CDK1、E2F1、Caspase、BAX、BAD、BCL2、NFκB、RAS、TGFα、TNFα、TGFβ、Wee1、MDMX、Wip1 依次編為 1～32 號基因，具體的非線性 ODEs 模型、參數取值及各吸引子的數值詳見文獻[7, 20]。通過對模型方程的模擬計算與參數擾動，可以找到關鍵作用邊中能使細胞從癌變狀態轉移到正常狀態的可行路徑共有 4 條，以┤表示抑制、→表示激勵，則分別為 R1（CDK2 ┤RB）、R10（P53 ┤ VEGF）、A9（HGF→AKT）、A11（AKT→VEGF）。4 條可行路徑上的調控反應參數符號及對應的臨界值如表 1 所示，以路徑 R1 為例，R1 表示 CDK2 對 RB 的抑制作用，CDK2 為 17 號基因，RB 為 6 號基因，故參數符號記為 ，標準值為無外部影響時的數值。
表1
可行路徑臨界值
Table1.
Critical values of feasible paths
圖7
癌癥 GRN 的吸引子網絡
Figure7.
The attractor network of cancer GRN
2.3 癌癥 GRN 優化調控
生物學研究表明，6 號基因 RB（）、10 號基因 AKT（）和 21 號基因 Caspase（）是表征細胞狀態的三個重要基因，在不同的細胞狀態下表達水平存在很大的差異，可作為癌癥 GRN 的狀態觀測項寫入目標函數中。將路徑 R1 和 A9 對應的調控反應參數同時作為控制變量，目標為控制 GRN 從癌變狀態 C 轉移至正常狀態 N 并且滿足控制時間最小化，據此建立如式（16）～式（21）所示的動態優化問題。另外兩條路徑 R10 與 A11 實現狀態轉移的效果較差，體現在其所需控制時間很長，同時考慮到程序計算效率，不將其納入最優控制策略計算中。
|
|
|
|
|
|
實現癌癥 GRN 狀態調控 C→N 的最優控制策略如圖 8 表示。路徑 R1 對應參數 單獨作用時，最優控制時間 T = 16.93，因為分岔后的過渡吸引子與表示正常狀態的吸引子很接近，故沒有返回標準值的變化趨勢。路徑 A9 對應參數 單獨作用時，實際最優控制時間 T = 44.28，相對于路徑 R1 來說所需控制時間較長。 和 同時作為控制變量作用于網絡時，最優控制時間 T = 11.07，其中對路徑 A9 的實際最優控制時間 T = 4.20。與單控制變量相比，多控制變量共同作用時，控制效率有顯著的提升。
圖8
狀態轉移 C→N 的最優控制策略
Figure8.
Optimal control strategies of state transition C→N
3 結論
對復雜生物網絡建立非線性 ODEs 模型，通過模擬計算和參數擾動得到系統的吸引子網絡，在此基礎上建立動態優化問題，利用動態優化算法求解，最終得到實現網絡狀態調控的最優控制策略。該方法具有普適性，既適用于簡單低維的生物網絡，也適用于與復雜高維的生物網絡；既適用于求解單控制變量問題，也適用于求解多控制變量問題。本研究所列舉的計算結果多為矩形脈沖形式，這可能與網絡本身特性有關，對于其他網絡不能得出矩形脈沖信號為最優策略的一般結論，但具有參考價值，必要時可以建立相應動態優化問題求解得到最優控制策略。另外，由計算結果可知，多控制變量共同作用時的控制效率往往高于單控制變量。
本文基于動態優化算法，提出了尋找實現復雜生物網絡狀態轉移最優控制策略的一般化方法，存在一定的應用前景。非線性 ODEs 模型均為實際生物網絡經過量化后建立的，狀態變量 x 與物質的實際濃度存在確定的函數關系。藥物的促進和抑制作用強度也經過量化，故控制變量 u 與藥物使用的劑量和周期存在一一對應的關系，即建立起了數學模型與實際情況的聯系。
大部分醫學實驗和實際治療中需要制定藥物的用量和使用周期，往往依照經驗而定，亟需定量化的參考。將本文的方法結論應用于實踐，可以指導以上參數的定量選擇，實驗中減小試錯次數、提高效率，實踐中提升療效。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
