大腸菌走化性運動における最適な感知と制御

Abstract

大腸菌の走化性応答はこれまで定量的実験により詳細に特徴づけられてきたが，応答が示す確率性や非線形性からその適応的意義を理論的に調べることに困難があった．部分観測最適制御理論を活用することで，大腸菌の走化性応答の最適性を検証し，その適応的意味を示した理論を紹介する．

1.  生物の感知・制御と非線形性

生物は時々刻々と変化する環境の状態を感知し，感知情報に基づいて応答を制御することで，環境の変化に適応的に応答する．多くの生き物でその適応応答の振る舞いは極めて巧妙で，進化の過程を経て何らかの意味で最適化されていることを示唆する．しかしそのことを検証するには，最適な適応的システムが持つべき性質や振る舞いを明らかにし，実際の現象と定性的・定量的に比較する必要がある．これまでに神経科学や生物物理学の分野で，最適な感知や運動制御をベイズ推定理論や最適制御理論などを用いて導き，実際の生物の振る舞いと比較することで適応性の議論がなされてきた^1),2)．

実際の生物の中で推定や制御を実現する化学反応や細胞ダイナミクスは高い非線形性を示し，それを最適理論の枠組みで扱うことは容易でない．そのためこれまでの多くの研究は，説明する現象を非線形性が無視できる線形応答領域に絞ることで理論を適用してきた．入力が弱く感知・制御系が大きく変化しない範囲では，応答が入力の強さに比例する線形応答が仮定できる．しかし応答の性質が分子機構も踏まえて詳細に整理されているにも関わらず，線形応答を大きく超えた理論がない故に適応性の議論が行き届かない系があった．本稿で扱う大腸菌の走化性はその一例である．

大腸菌の走化性を担うシグナル伝達系は，早くから工学・物理学的な発想を取り入れた定量的な研究が進められてきた．1980年代にはすでに，弱い入力に対する応答が系統的に測定され，その領域の振る舞いは線形応答で説明できることが確認されている³⁾．2000年代には，大腸菌の示す線形応答領域の性質が推定・制御の観点から見て適応的なものであるのか，線形応答の理論に基づき議論された^4)-6)．非線形応答まで含めた性質についても，遺伝子欠損株を用いた定量測定と物理化学的な理論を組み合わせた特徴づけが進められ，2010年頃までに分子機構を踏まえた生物物理モデルに整理されている⁷⁾．しかし非線形応答の領域に見られる性質がどのような適応的意義を持つのかについては定まった見解がなく，特に最適な推定・制御の観点から説明することは未解決問題として残されていた．この課題を克服したのが本稿で紹介する著者らの研究である^8),9)．

著者らは線形応答に限らない応答の最適化を扱いうる理論を開発し，最適な推定・制御則が大腸菌走化性の分子機構を踏まえた数理モデルと対応して非線形応答まで含めて説明できることを示した．また我々は推定と制御の最適性を包括的に扱う部分観測確率制御の理論¹⁰⁾を用いることで，大腸菌は自ら下した制御が走化性運動に与える影響を織り込んだ上で最適な推定をしている可能性を見出した．本稿では未だ生物学への応用に乏しい部分観測確率制御の理論が，非線形応答を示す大腸菌走化性にどう適用されて何を明らかにしたのかを解説する．2節で大腸菌走化性の概要を説明した上で，3・4節では部分観測確率制御の理論で扱えるように定式化する過程を説明する．5節では部分観測確率制御の理論から帰結する内容を説明し，6節では実際の大腸菌と比較して明らかになったことを解説する．

2.  大腸菌のrun-and-tumble型走化性における感知と制御

大腸菌はべん毛の回転方向を切り替えることで，直進するrun状態とランダムに方向転換するtumble状態をとる．両者を組み合わせた運動パターンはrun-and-tumble運動と呼ばれ，tumbleの起こる頻度（tumbleレート）を受容体で感知するリガンド濃度の時間変化に基づいて制御することで走化性を実現している．濃度が増大するのはリガンド勾配を上っているときなので，tumbleレートを下げれば上り続けることが期待される．逆に濃度が減少するのは勾配を下っているときなので，tumbleレートを上げればさらに下るのを防ぐことが期待される．濃度の時間変化の正負に応じた制御が走性を実現する肝になるが，正負の判定はリガンド感知に伴う確率性のために容易でない．濃度勾配を上っているのにリガンドとレセプタの結合回数が減少するというような，相反する観測が確率的に得られる．感知の確率性のもとで，濃度の時間変化の正負をどうすれば良く推定できるのか，また推定の不確実性を加味してtumbleレートをどう制御すれば良いのか，という問題を生物自身が解く必要がある．

3.  状態空間モデルによる感知と制御の定式化

最適な感知・制御の問題を理論的に扱う第一歩は，推定・制御問題を定式化する一般的な枠組みである状態空間モデルに落とし込むことである．状態空間モデルでは，推定・制御したいが直接は観測できない状態変数の時間発展と，状態変数に関する観測の確率的な得られ方を定式化する（図1）．

図1

Run-and-tumble型走化性における感知・制御過程．

まず大腸菌の運動をモデル化する．リガンド濃度の時間変化の正負に基づく制御を捉えた単純な設定として，濃度が単調に変化する軸に沿った1次元空間上でのrun-and-tumble運動を考える．時刻tにおける大腸菌の位置と方向をそれぞれξ_tとX_tで表し，位置ξ_tは一定の速さvで方向X_tに動くと仮定することで，ξ_tの時間発展がdξ_t/dt = vX_tと表せる．方向X_tの時間発展については，tumbleレートR_tの頻度で確率的にtumbleが起こり，その度に方向がランダムに振り直されて確率1/2でX_t = +1かX_t = –1になると仮定する．数理的には，十分小さい時間幅dtで方向の反転する確率が

  

$$\mathbb{P}\left(X_{t+\text{d}t}=-X_t∣X_t,R_t\right)={\left(1/2\right)R}_t\text{d}t,$$（1）

に従う連続時間マルコフ連鎖でモデル化できる．

次に確率的なリガンド濃度感知を観測過程としてモデル化する．観測対象であるリガンドが拡散等で指数的に分布している状況を考えると，時刻tの大腸菌の位置ξ_tにおける濃度[L]_tについて ${\left[L\right]}_t\propto \mathrm{\text{exp}}\left(c{\xi }_t\right)$ が成り立つ．c > 0は勾配の傾きを表す定数である．大腸菌は濃度の絶対値に大きく依存せず，対数の時間微分に応答することが知られている^3),7)．そこで，対数濃度の時間微分にガウスノイズが加算された観測変数dY_t = Y_t+dt – Y_tを考える．数理的には，連続時間のガウスノイズを標準ウィーナー過程の増分dW_tで表すことで，確率微分方程式$\text{d}{Y}_t=\text{d}{\log \left[L\right]}_t+\sqrt{\sigma }\text{d}{W}_t$により定式化できる．σ > 0は2節で説明した感知の確率性に由来するノイズの大きさを表す定数である．[L]_t，ξ_t，X_t間の関係式を代入することで

  

$$\text{d}{Y}_t=\mathit{\text{cv}}{X}_t\text{d}t+\sqrt{\sigma }{\text{d}W}_t,$$（2）

が得られ，これは状態X_tが与えられた元で観測変数dY_tがとる値の確率 $\mathbb{P}\left(\text{d}{Y}_t∣X_t\right)$ を規定している．

以上で大腸菌走化性を，状態変数の時間発展 $\mathbb{P}\left(X_{t+\text{d}t}∣X_t,R_t\right)$ （式（1））と観測 $\mathbb{P}\left(\text{d}{Y}_t∣X_t\right)$ （式（2））からなる状態空間モデルにまとめることができた．大腸菌はX_tを直接知ることはできないが，dY_tの値を利用することができる．また大腸菌にとって式（1），（2）やそこに現れるパラメータで表される運動と観測の統計的な性質は既知であると仮定する．運動と観測の統計的な性質は方向X_tに比べて変動が緩やかで推定しやすい可能性が高く，その場合運動と観測の性質が持つ不確実性を無視して既知と仮定しても本稿で着目するX_tの推定には大きく影響しないことが期待される．

4.  最適化問題（制御関数と目的関数）の定式化

TumbleレートR_tを因果的に制御するために，原理的には時刻tまでの観測履歴Y_0:t = {Y_s|0 ≤ s ≤ t}を使うことができる．Y_0:tが持つ情報を大腸菌が保持できるかは非自明であるが，比較の基準となる最適な制御則を考えるため，Y_0:tに依存しうる制御 $R_t={\mathcal{R}}_t\left[Y_{0:t}\right]$ を想定する．ここで制御関数 ${\mathcal{R}}_t$ は線形応答に限定しない．

最適化の基準となる目的関数として，制御コストを低く抑えつつ濃度勾配を素早く上る設定を考える：

  

$$J=\mathbb{E}\left[{\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{-\gamma t}\left[v{X}_t-\frac{1}{\beta }\{R_t\text{log}\frac{R_t}{R^0}-\left(R_t-R^0\right)\}\right]dt\right].$$（3）

e^–γtは時間割引因子と呼ばれ，定数γ > 0が大きいほど最適化にあたりより近い未来を重視する設定を表す．第1項vX_tは時間積分されることで勾配を正味上る距離を，第2項はtumbleレートR_tの制御に伴うコストを表す．1/βは制御コストに対する重みづけである．

制御コストをモデル化するために，多くの制御問題では二乗コストが使われる．しかしR_tが大きく変動する非線形応答まで考える本稿の設定では，二乗コストのもと最適化するとR_tの非負性が破れるという不自然な結果が現れてしまう．非負性の破れを解消しつつ解釈しやすい定式化として，著者らはKullback-Leibler（KL）制御を用いた．KL制御では非制御状態の存在を仮定し，制御を加えることで非制御状態から乖離する度合いを制御コストと捉えてKLダイバージェンス ${\mathcal{D}}_{\text{KL}}\left[{\mathbb{P}}^{\mathrm{\text{制御}}}\parallel {\mathbb{P}}^{\mathrm{\text{非制御}}}\right]$ で定式化する．ここで ${\mathbb{P}}^{\mathrm{\text{制御}}}$ と ${\mathbb{P}}^{\mathrm{非制御}}$ はそれぞれ制御/非制御状態におけるtumbleの時系列の現れ方を表す確率分布である．大腸菌走化性では濃度勾配が存在しないときの自由運動を非制御状態と解釈でき，tumbleがある一定値R⁰で自発的に起きていると仮定することで式（3）の制御コストが導かれる．R_tの非負性の破れが解消されることはこのあとわかる．

5.  部分観測確率制御から導出される最適制御則

状態空間モデルに基づく以上の定式化に部分観測確率制御の理論¹⁰⁾を適用することで，式（3）の目的関数を最大にする意味で効率的に勾配を上る最適制御則が導出される．最適制御則は，時々刻々と得られる観測変数dY_tに基づいて逐次的に推定を更新する部分と，推定に基づいてtumbleレートR_tを決める部分で構成される（図1）．

まず最適な制御の実現に求められる推定は，事後確率 $Z_t=\mathbb{P}\left(X_t=-1∣Y_{0:t}\right)$ で与えられる．Z_tは時刻tまでの観測Y_0:tから推定される，勾配を下っている可能性を表す量である．Z_t = 1/2は方向がX_t = ±1のどちらか不明な無情報の状態に対応し，Z_t = 0, 1に近づくほど不確実性が低い状態を表す．推定Z_tを更新する方法は，非線形フィルタリング理論においてKushner方程式として知られるZ_tの確率微分方程式で表せる：

  

$$\text{d}{Z}_t=-R_t\left(Z_t-1/2\right)\text{d}t-K{Z}_t\left(1-Z_t\right)○\text{d}{Y}_t.$$（4）

ここで$○$は確率微分方程式の表し方の一つであるストラトノビッチ積分であり，このあと行う座標変換の結果と対応がとりやすいように採用している．第1項は，運動モデル（式（1））に基づく予測を表す項である．Tumbleが確率的に起こることで方向がランダム化されるので，観測が無ければ推定は不確実性を増し続ける．そのためdY_t = 0ならば推定は無情報の状態Z_t = 1/2に収束していき，その速さはtumbleレートR_tに比例する．第2項は式（2）の観測dY_tに基づく更新を表し，係数のK = 2cv/σは観測のシグナルノイズ比である．

tumbleレートR_tの最適な制御は，式（4）に従って観測履歴Y_0:tの持つ情報を要約した事後確率Z_tに基づいて実現される．具体的には以下の最適制御関数によりR_t = R*(Z_t)と制御できる：

  

$$R^{\mathrm{*}}\left(Z_t\right)=R^0\text{exp}\left[-\beta \left(Z_t-1/2\right)\text{d}V/\text{d}Z\left(Z_t\right)\right].$$（5）

V(Z_t)は価値関数と呼ばれ，最適制御理論でHamilton-Jacobi-Bellman方程式と呼ばれる2階の微分方程式を解いて得られる．V(Z_t)の微分がexpの引数に現れる形はKL制御に基づくコストの導入から帰結しており，最適化で現れるR_tの非負性が保証されている．V(Z_t)を数値的に求めて式に代入することで，図2Aに示すR*(Z_t)が得られる．R*(Z_t)はZ_t = 1/2でR⁰を通るZ_tの増大関数であることがわかる．これは無情報の推定状態Z_t = 1/2では非制御状態R_t = R⁰をとり，Z_t = 1に近づき勾配を下っている可能性が高いほどtumbleレートR_tを上げるという，直感的にも理解できる結果である．

図2

最適制御則に現れる関数の形．（A）tumbleレートを決める最適制御関数（B）負のフィードバックを担う関数Fの最適な形と実験データとの比較．最適制御では縦軸と横軸にF*とZ_tを，実験データでは縦軸と横軸にFとa_tをそれぞれプロットした．

6.  最適な推定・制御則と大腸菌の応答の対応

理論的に導かれた最適な推定・制御則を実際の大腸菌の応答と比較することで，大腸菌が効率的に走化性をしているかを検証できる．先述の通り大腸菌の応答は詳しく測定されており，背後の分子的機構も踏まえて定量的に応答を再現する生物物理モデルがすでに整備されている．もし最適な制御則が生物物理モデルと対応すれば，生物物理モデルで説明される大腸菌の応答は最適制御を実現する性質を持っているといえる．

著者らは最適制御則を適切に等価変形することで，実際に両者が対応することを見出した．まず式（4），（5）で表される最適制御則は，Z_tと1対1対応する新しい座標θ_t: = log{(1 – Z_t)/Z_t}を導入し，予測項を新しい変数μ_tで表すことで次のように書き直すことができる．

  

$$Z_t=1/\{1+\text{exp}\left({\theta }_t\right)\}$$（6）

  

$${\theta }_t=-\kappa {\mu }_t+K\log {\left[L\right]}_t+K\sqrt{\sigma }{W}_t$$（7）

  

$$\text{d}{\mu }_t/\text{d}t=F^{\mathrm{*}}\left(Z_t\right):=-R^{\mathrm{*}}\left(Z_t\right)\cdot \left(Z_t-1/2\right)/\left\{Z_t\left(1-Z_t\right)\right\}/\kappa$$（8）

  

$$R_t=R^{\mathrm{*}}\left(Z_t\right):\text{増加関数}$$（9）

κ > 0は任意定数である．対して大腸菌の濃度[L]_tに対するtumbleレートR_tの応答は，以下の生物物理モデル[7]で再現できる．

  

$$a_t=1/\{1+\text{exp}\left(f_t\right)\}$$（10）

  

$$f_t=-N\alpha m_t+N{\text{log}\left[L\right]}_t$$（11）

  

$$\text{d}{m}_t/\text{d}t=F\left(a_t\right):a_t=\stackrel{-}{a}\text{にゼロ点を持つ減少関数}$$（12）

  

$$R_t=R\left(a_t\right):\text{増加関数}$$（13）

このモデルは多数のレセプタ分子からなる大腸菌レセプタの状態を，活性分子の割合であるレセプタ活性a_tと各分子の平均的なメチル化度合いであるメチル化レベルm_tの2変数に縮約して説明したものである．またf_tはa_tを決める活性／非活性状態間の自由エネルギー差，N，αは生化学定数である．

生物物理モデルの4変数（a_t, f_t, m_t, R_t）は，最適制御則の4変数（Z_t, θ_t, μ_t, R_t）とそれぞれ対応している．レセプタ活性a_tは自由エネルギー差f_tを通して対数濃度log[L]_tとメチル化レベルm_tから定まり，tumbleレートを増大させる．これはZ_tがθ_tを通してlog[L]_tとμ_tで定まり，最適制御関数R*が増大関数であることと対応する．関数F(a_t)で決まるメチル化レベルm_tの時間発展は，レセプタ活性を平衡点a_t = āに収束させる負のフィードバックを構成している．これはμ_tの時間発展がZ_t = 1/2をゼロ点に持つ減少関数F*で表され，事後確率をZ_t = 1/2に収束させる負のフィードバックを構成していることと対応する．

著者らは以上の対応を足がかりに，さらに最適制御から導出されたZ_t全域にわたる関数形F*をShimizuら¹¹⁾がa_t全域で計測したFと比較し，その整合性を検証した．その結果が図2Bである．両者は平衡点（Fのゼロ点）から離れた非線形応答の領域まで含めて良く整合している．

特に着目すべき点は，Z_t = 1/2を挟んで左側Z_t < 1/2に比べて，右側Z_t > 1/2ではFの絶対値がより大きいという非対称性である．この性質は，自らの制御が走化性運動に与える影響を式（1）の運動モデルに基づいて予測している結果と解釈できる．実際，仮に制御関数R*(Z_t)が定数だった場合のF*はZ_t = 1/2に対してその絶対値が対称であり，F*の非対称性はR*(Z_t)が増大関数であることで初めて現れる．TumbleレートR_tが大きいほど方向X_tのランダム化が速く進むため，制御によりR_tを増大させるZ_t > 1/2の領域では推定の信頼性がより速く失われると予測できる．そのためZ_t > 1/2ではF*の絶対値を大きくし，Z_tをより速く無情報の状態Z_t = 1/2に収束させることで，制御の影響に関する予測を推定に織り込んでいると解釈できる．

生化学的な観点からは，Fの非対称性はレセプタのメチル化を担う酵素CheRと脱メチル化を担う酵素CheBのうち，後者のCheBのみが活性状態のレセプタ分子によるリン酸化修飾を受けることで現れると考えられている．CheRにはなくCheBのみに存在するリン酸化制御の機能的役割は未解明であるが，本モデルはtumbleレートの制御を考慮した最適な推定に貢献しているという新たな仮説を提示している．

7.  まとめと今後の課題

本稿では線形応答に限らない推定・制御則の最適性を部分観測確率制御やKL制御に基づいて定式化することで，大腸菌走化性の生物物理モデルと良く対応し非線形応答を説明する最適制御則が得られたことを紹介した．部分観測確率制御の理論は，本稿の例に限らず制御下の状態遷移 $\mathbb{P}\left(X_{t+\text{d}t}∣X_t,R_t\right)$ と観測 $\mathbb{P}\left(\text{d}{Y}_t∣X_t\right)$ をモデリングすれば広く適用できて，最適制御則が事後確率 $\mathbb{P}\left(X_t∣Y_{0:t}\right)$ の時間発展とそれを引数としてとる制御関数の組合せとして導出される．しかし，より複雑な環境変動や多数の観測情報が存在する状況では，部分観測確率制御理論から最適な制御関数を導出することが極めて困難になり，加えて導かれる最適な制御は生物が細胞内などで実現できる複雑さにとどまる保証がない．また生体が機能することに伴うコストを，例えば熱力学なども考慮してより適切に扱う必要もある．さらに高次の生体情報処理に対してその最適性を調べるために，これら生物特有の問題を扱いうる新たな理論を発展させることが今後の課題である．

文献

• 1)  田中 宏和 (2019) 計算論的神経科学：脳の運動制御・感覚処理機構の理論的理解へ，森北出版，東京．
• 2)   Mattingly,  H. H. et al. (2021) Nat. Phys. 17, 1426-1431. DOI: 10.1038/s41567-021-01380-3.
• 3)  Berg, H. C. (2004) E. coli in Motion, Springer, New York.
• 4)   De Gennes,  P. G. (2004) Eur. Biophys. J. 33, 691-693. DOI: 10.1007/s00249-004-0426-z.
• 5)   Clark,  D. A.,  Grant,  L. C. (2005) Proc. Natl. Acad. Sci. U.S.A. 102, 9150-9155. DOI: 10.1073/pnas.0407659102.
• 6)   Andrews,  B. W. et al. (2006) PLoS Comput. Biol. 2, e154. DOI: 10.1371/journal.pcbi.0020154.
• 7)   Tu,  Y. (2013) Annu. Rev. Biophys. 42, 337-359. DOI: 10.1146/annurev-biophys-083012-130358.
• 8)   Nakamura,  K.,  Kobayashi,  T. J. (2022) Phys. Rev. Res. 4, 013120. DOI: 10.1103/PhysRevResearch.4.013120.
• 9)   中村 絢斗， 小林 徹也（2022）システム／制御／情報 66, 346-351. DOI: 10.11509/isciesci.66.9_346.
• 10)   Segall,  A. (1977) IEEE Trans. Autom. Contr. 22, 179-186. DOI: 10.1109/TAC.1977.1101447.
• 11)   Shimizu,  T. S. et al. (2010) Mol. Syst. Biol. 6, 382. DOI: 10.1038/msb.2010.37.

Biographies

中村絢斗（なかむら　けんと）

理化学研究所基礎科学特別研究員

小林徹也（こばやし　てつや）

東京大学生産技術研究所教授