最近見た面白い論文のコーナー

Special functions and the Mellin transforms of Laguerre and Hermite functions

https://arxiv.org/abs/math-ph/0612085

Hermite多項式とLaguerre多項式のMellin変換は、ゼロ点がRe(s)=1/2上にある(そして、全て単純零点)という事実が、1980年代から知られているらしい。

このMellin変換が超幾何関数の特殊値で表せて(パラメータ(a,b,c)がsの関数として表せる)、関数等式が超幾何関数のKummerの関数等式の帰結として出てくるというのが新しい発見のよう(そして、Kummerの関数等式は、GL(4)のWeyl群作用からの解釈がある)

そしてまた、任意次元の水素原子の固有エネルギー状態のMellin変換(この場合のMellin変換の定義は論文参照)も、Re(s)=1/2上にのみゼロ点を持ち、いずれも単純零点であるとのこと。これの背景には、調和振動子の適当な波動関数($O(n)$対称性のある既約表現に属するという条件がついている)のMellin変換を考えると、やはりRe(s)=1/2上にゼロ点があるという話があるらしい。

A Local Riemann Hypothesis, I

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.23.7311


水素原子の束縛状態の空間(3次元の場合はso(4,2)の極小表現)は、振動子表現(Weil表現/Segal-Shale-Weil表現)の部分空間としての実現を持つけど、表現の異なる実現に於いて、このような性質を共有するのは意外。KS変換する前と後でMellin変換の結果が同じものになるのかどうかは確認していないけど、もし一致するのであれば、表現の実現の仕方に依存しない形で、Mellin変換が理解できるのか興味がある。わたしは数論には興味がないし、Riemann予想は解けたら何が嬉しいのか、さっぱり分からないので、今まで興味がなかったけど、表現論的に何か面白いことがあるのかもなぁという感じ


調和振動子の対称性とspectrum generating algebra

ただのSegal-Shale-Weil表現(振動子表現という名前もある)の話。無限次元既約ユニタリー表現の研究は、Lorentz群やPoincare群で始まったけど、ある意味では、それより以前から知られていた、最初の無限次元既約ユニタリー表現の例でもある。


水素原子の表現論

http://d.hatena.ne.jp/m-a-o/20140130#p1

水素原子の表現論(1.5-2)spectrum generating algebra

http://d.hatena.ne.jp/m-a-o/20161010#p1

などで、

・3次元Kepler系の相空間には、群$SO(4,2)$が推移的にsymplectic変換として作用し、symplectic等質多様体であること

・水素原子の束縛状態の空間は、Lie環$so(4,2)$のある既約ユニタリー表現(極小表現)であり、$so(4,2)$はspectrum generating algebraという名前で呼ばれること

を書いた。また、Hamiltonianと交換する保存量として、角運動量とRunge-Lenzベクトルがあり、対応して、古典的にはエネルギー等値面上に、(エネルギーの大きさに応じて)$SO(4)$や$SO(3,1)$などの群が作用している。そして、水素原子に於ける主量子数はLie環$so(4,2)$の既約ユニタリー表現を$so(4)$に制限した時の分岐則として理解できる。というわけで、Kepler系では、相空間の幾何学的対称性とHamiltonianの対称性が、古典力学と量子力学を共に統制している様子を見ることができる。



調和振動子については、Heisenberg群/Heisenberg代数の表現論による記述がよく知られている。調和振動子の相空間$\mathbf{R}^{2n}$はHeisenberg群の余随伴軌道であって、幾何学的量子化が完璧に機能する。一方、(振動数が全て同一の$D$次元)調和振動子が、$su(D)$の(無限小)対称性を持つことも知られている。$su(D)$対称性は、(知らなくてもエネルギー固有値を計算できるし、そもそも一次元しか扱わない場合も多いので)量子力学の教科書には必ずしも書いていないことが多い気がするけど、調和振動子のエネルギーの縮退は、$su(D)$対称性の存在によって説明できる。調和振動子は明らかに回転対称性を持つが、これは、この調和振動子の対称性の一部で、このあたりの事情は、水素原子のエネルギーの縮退が、回転対称性より大きな対称性代数$so(4)$の表現によって説明されるのと似ている。


ところで、調和振動子の対称性$su(D)$はHeisenberg代数の部分Lie環になっていない。これ自体は、珍しいことでもないので、別に構わないという立場もありうるけど、「水素原子の表現論(1.5-2)」で触れたJosephの論文

Minimal realizations and spectrum generating algebras

https://projecteuclid.org/euclid.cmp/1103859776

には、symplectic代数$sp(2n,\mathbf{R})$が調和振動子のspectrum generating algebraであると書いてある。シンプレクティック群$Sp(2D,\mathbf{R})$は$\mathbf{R}^{2D}$に推移的には作用してないけど、稠密な開軌道は持つ。


Heisenberg代数とsymplectic群/代数との関係は、symplectic群がWeyl代数の環同型として作用することから来る(Weyl代数が、Euclid空間の余接空間の量子化であることを思い出せば自然)(Weyl代数の自己同型群自体は、もっとずっと大きい。 arXiv:math/0512169などを参照)。調和振動子の状態空間は、Heisenberg代数の既約ユニタリー表現空間であるが、この既約ユニタリー表現は、symplectic代数とHeisenberg代数の半直積の既約ユニタリー表現に一意に拡張でき、Segal-Shale-Weil表現と呼ばれる。この既約ユニタリー表現をsymplectic代数に制限すると、可約になって2つの既約表現に直和分解する。この2つの既約表現は極小表現になっていることが知られている。


注)Segal-Shale-Weil表現は、シンプレクティック群の射影表現にしかならない。量子力学で、幾何学的対称性の群の射影表現を考える必要が最初に生じたのは、スピン表現だと思うけど、事情は似ていて、スピン表現が$SO(3)$の二重被覆$SU(2)$の表現にはなっているように、シンプレクティック群の二重被覆であるメタプレクティック群の表現にはなる。まぁ、群で考えると、計算が面倒になるだけで特にメリットもないので、以下計算する時はLie環の方で考える。



というわけで、調和振動子の状態空間への$sp(2D,\mathbf{R})$の作用や、Hamiltonianの対称性$su(D)$の$sp(2D,\mathbf{R})$の部分Lie代数としての実現を見て、分岐則を調べれば、調和振動子のエネルギー固有値ごとに縮退した状態が、対称性$su(D)$の既約表現として得られるっぽいことを見る。この分岐則は、Segal-Shale-Weil表現の説明に必ず書いてあるようなことではあるし、特に複雑な計算もない(が、よい教科書は知らない)。以下、D>1を仮定する(D=1の時は、$su(1)$は自明になり、エネルギーに縮退もないので)


bosonの生成・消滅演算子と調和振動子のHamiltonianを、通常のように

 ¥[a_i^{-} , a_j^{+}¥] = ¥delta_{ij}

 ¥[a_i^{¥pm} , a_j^{¥pm} ¥] = 0

 H = ¥omega ¥sum_{k=1}^{D} (a_{k}^{+}a_k^{-} + ¥frac{1}{2})

で定める。容易な計算で

 X_{ij} = a_i^{+} a_j^{-}

は全てHamiltonianと可換になることが分かる。また簡単な計算によって

 ¥[X_{ij} , X_{kl} ¥] = ¥delta_{jk} X_{il} - ¥delta_{il} X_{kj}

なので、これがLie環$gl(D,\mathbf{C})$と同型なLie代数の基底をなすことが分かる。


Hamiltonianの対称性だけ見ていると、これで十分なのだけど、(後でspectrum generating algebraの基底も見ると)これでは(複素化した)spectrum generating algebraの部分Lie環にならないことが分かるので、

 X_{ij} = a_i^{+} a_j^{-} + ¥frac{1}{2} ¥delta_{ij}

とする必要がある。これも、依然としてLie環$gl(D,\mathbf{C})$の基底となるけど、おまけで付いてくる1/2のせいで、この表現は、一般線形群$GL(D,\mathbf{C})$の表現には持ち上がらない(群に拘りたい場合は、一般線形群の二重被覆metalinear群を考えることもできる)。しかし、$sl(D,\mathbf{C})$のCartan部分環は、標準的な基底として

 h_{i} = X_{ii} - X_{i+1,i+1}

を取ることができ、うまいこと1/2が相殺してくれるので、特殊線形群の表現には持ち上がる。$sl(D,\mathbf{C})$に入らない残りの一次元分の対称性は

 ¥sum_{k=1}^D X_{ii} = ¥sum_{k=1}^D (a_k^{+} a_k^{-} + ¥frac{1}{2})

で、定数倍を除けばHamiltonianそのものである。


#群に表現が持ち上がるかどうかという観点では、$SL(D,\mathbf{C})$や$SU(D)$は単連結だが、$GL(1)$や$U(1)$はそうでないというあたりから、面倒さがやってくる。$D$変数多項式環への$gl(D,\mathbf{C})$の作用として、行列単位に対して

 E_{ij} = x_i ¥partial_j + s ¥delta_{ij}

を取れる。対応する一般線形群$GL(D,\mathbf{C})$の表現は形式的には

 ¥pi_{s}(A)f(¥mathbf{x}) = ¥det(A)^{s} f(¥mathbf{x} A)

であるが、$s$が非整数の時、群の表現の方は$(det A)^{s}$が一価でない。しかし、$det A=1$の時に限定すれば問題はなくなり、特殊線形群$SL(D,\mathbf{C})$の表現は定まる。


spectrum generating algebraの実現。まず、複素化した$sp(2D,\mathbf{C})$を考えたほうがわかりやすい。次のように基底を取る

A_{ij} = ¥begin{pmatrix} E_{ij} & O ¥¥ O & -E_{ji} ¥end{pmatrix}

B_{ij} = ¥begin{pmatrix} O & E_{ij}+E_{ji} ¥¥ O & O ¥end{pmatrix}

C_{ij} = ¥begin{pmatrix} O & O ¥¥ E_{ij}+E_{ji} & O ¥end{pmatrix}

対応する演算子をbosonの生成・消滅演算子を使って

X_{ij} = a_i^{+} a_j^{-} + ¥frac{1}{2} ¥delta_{ij}

Y_{ij} = a_i^{+} a_j^{+}

Z_{ij} = a_i^{-} a_j^{-}

とする。これが、symplectic代数と同型なLie環をなすことは計算するのみ。Heisenberg代数にderivationとして作用していることも簡単な計算で分かる。Cartan対合$s$を

s(A_{ij}) = -A_{ji}

s(B_{ij}) = C_{ij}

で定める。このCartan対合に関する実型は$sp(2D,\mathbf{R})$と同型(証明略)。特に、上で作ったLie環$gl(D,\mathbf{C})$は部分Lie環になってることが分かり、このLie環の元のうち、Cartan対合で不変なもの全体は実Lie環$u(D)$をなす。


以上の定義によって、調和振動子の状態空間への$sp(2D,\mathbf{C})$と$sl(D,\mathbf{C})$(あるいは$sp(2D,\mathbf{R})$と$su(D)$)の作用が定まった(ユニタリー表現になることなどは示されてないけど)。symplectic代数の作用を見ると、振動子の個数を変化させないか、2個増減させることが分かるので、少なくとも振動子の個数の偶奇によって、2つの表現に直和分解する。


調和振動子の状態空間を、$su(D)$の表現として既約分解する。試みとして、|0>を基底状態として

 v_i = a_{i}^{+} |0 ¥rangle

を基底とする、ベクトル空間(振動子が一個の状態全体)を考える。そうすると

 a_{i}^{+} a_{j}^{-} (v_{k}) = ¥delta_{jk} v_{i}

なので、$su(D)$が自然表現として作用していることが分かる。一般に振動子$N$個の状態空間は、$su(D)$の自然表現の$N$次対称テンソル積表現となっていて、これらは全て既約表現である($N$次斉次多項式の空間と同一視できるので、適当に基底を取って、最高ウェイトetc.を直接計算すれば示せる。また、この事実から、上で書いた$sp(2D,\mathbf{R})$の2つの表現が既約であることが従う)。$su(D)$を複素化すると$sl(D,\mathbf{C})$なので、振動子の状態空間は、$sl(D,\mathbf{C})$の自然表現の対称テンソル積代数(必要なら、適当に完備化したもの)と同定できる。


Liouville可積分性について。調和振動子の場合、$su(D)$のCartan部分環は、$D-1$個の独立でinvolutiveな保存量を与える。Hamiltonian自身を合わせて、相空間の次元の半分の個数のinvolutieな保存量が得られる。各振動子の振動数ωが異なっていても、Liouville可積分ではあるが、対称性は、$su(D)$の部分代数に落ちていく。最も退化した場合は、$u(1)^{D-1}$となる


Heisenberg vs symplectic。量子Kepler系/水素原子との関連で考えると、Heisenberg代数でなく、symplectic代数をspectrum generating algebraとした方が自然ではある。水素原子については、調和振動子の$U(1)$作用によるHamiltonian reductionとして得られる(散乱状態の場合は、"虚の振動数"を持つHamiltonianを考える)けど、

M = ¥sum_{i=1}^2 X_{ii} - ¥sum_{i=3}^{4} X_{ii}

と可換な$sp(8,\mathbf{R})$の部分Lie代数が、$u(2,2)$と同型であるということが、Kepler系のspectrum generating algebraが$su(2,2)$であることの、1つの説明を与える。この$M$は$U(1)$作用の量子化運動量写像を与えると解釈でき、$u(1)$部分は、Hamiltonian reductionを取る際に0に取るので、結局$su(2,2)$が残る


#古典的な場合、Kustaanheimo-Stiefel変換と、このHamiltonian reductionの関係は

On the regularization of the Kepler problem

https://projecteuclid.org/euclid.cmp/1103921049

の3章に簡潔にまとまっている。


#また、Kepler系の相空間は、$su(2,2)$の極小軌道であった。一般に(symplectic群の)dual pairがあった場合、pairの片割れによる簡約相空間は、残りのpairの余随伴軌道(の閉包)と同型になることは、漠然と予想されているっぽい。が、この対応が、いつ正しいのか、正確なところは、よく分からない。

Examples of singular reduction

http://www.math.cornell.edu/~sjamaar/papers/lms.pdf

の4節(特に定理4.3)には、かなり限定的な条件ではあるが、dual pairの簡約相空間と余随伴軌道の対応について、証明がある(1990年頃のものなので、今ではもっと一般的な事実が知られていてもよさそうではある)。Kepler系の相空間が、極小軌道となることも、こうした事実の一類型となっている。



余談。Segal-Shale-Weil表現は、symplectic代数の表現としては既約表現になってない点に不満がある。Heisenberg代数とsymplectic代数の半直積を考えるという道もありそうではある。群としては、シンプレクティック群とHeisenberg群の半直積を考えることになり、この群はJacobi群という名前が付いている。Jacobi群という名前は

Jacobi group

https://en.wikipedia.org/wiki/Jacobi_group

によると、EichlerとZagierによる命名らしい。彼らはJacobi形式の研究のために、使ったとか。物理の方では、squeezed stateをJacobi群のgeneralized coherent statesとして理解できる、と書いている人がいる。私は、量子光学を知らないので、squeezed状態が何か知らないし、これが妥当なのかも確認していない(Segal-Shale-Weil表現は、Jacobi群の表現としては当然射影表現にしかならない)

Generalized squeezed states for the Jacobi group

https://arxiv.org/abs/0812.0717

The Jacobi group and the squeezed states - some comments

https://arxiv.org/abs/0910.5563



余談その2。私は全然興味がないけど、テータ関数の変換公式をWeil表現から導出するという話があったりするらしい

リー群の表現から見たテータ関数

http://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/48222

で、主に、$SL(2,R)$の場合に絞って考え方を理解しやすく解説してくれている


調和振動子は量子力学の模型としては簡単で面白くもない部類のものだけど、振動子表現は、表現論的になかなか面白いものである。


その他の初等的な可積分系について。

まず、(massiveな)自由粒子。自由粒子は教科書の最初に出てくる最も簡単な系にも関わらず、表現論的に理解するのは、それほど自明ではない。明らかに重要な代数として、(非相対論的な場合には)ガリレイ代数、(相対論的な場合には)ポアンカレ代数がいる。こいつらは、Hamiltonianと可換でない演算子(Galilean boost,Lorentz boost)を含み、対称性代数をこれらの部分代数と考えることができる。一方、Wignerの仕事によって、粒子状態の空間は、ガリレイ代数/ポアンカレ代数の既約ユニタリー表現と理解できるはずで、これを対称性代数に制限した時の既約分解を考えると、エネルギー固有状態に対応する表現が得られると予想される。エネルギースペクトルが離散的でないことから、分岐則は、直和の形ではなく直積分の形になると思われる。


#ガリレイ代数は、粒子の質量と関連する中心拡大を持つ。相空間に作用する幾何学的変換としては、この項は見えない。波動関数の変換としては、位相の回転と関連し、ガリレイ群の射影表現を考える必要があることの表れでもある。まぁ、群で考えると、計算が面倒になるだけで特にメリットもないのでry


#masslessな自由粒子としては、光子が代表的であるけども、これについては

Maxwell方程式の表現論(1)(2)

http://d.hatena.ne.jp/m-a-o/20150621#p1

http://d.hatena.ne.jp/m-a-o/20151106#p2

で書いたものに相当する。これらは量子力学系に相当するけど、対応する古典論は"幾何光学"と呼ばれるべきものである。通常の意味の幾何光学と同一なのか知らないけど、正準理論の起源が幾何光学にあった(らしい)ことを思うと、面白い。またmasslessなスカラー粒子は、素粒子として存在するのか知らないけど、対応する既約ユニタリー表現は、水素原子の束縛状態全体の空間と同じである。masslessスカラー粒子の場合と、水素原子の場合では、Hamiltonianは異なるので、対称性と分岐則も異なる


個人的に、よく分からない例として、自由剛体がある。可積分なコマとしては、Euler top,Lagrange top,Kovalevski topがあるけど、特にEuler topを理解したい。この系はHamiltonian以外に1つの保存量を持つLiouville可積分系であり、古典的には、楕円関数を使って解を書ける。量子論の場合、明らかに、$so(3)$の既約ユニタリー表現空間に作用し、それは全て有限次元空間なので、Hamiltonianの行列要素を具体的に書くことができ、原理的には解ける。とはいえ、有限次元空間に作用しているから解けますというのでは、Heisenberg模型(XXX模型)なんかも自明ということになってしまう。

Lamé equation, quantum top and elliptic Bernoulli polynomials

https://arxiv.org/abs/math-ph/0508068

なんかは、quantum Euler topの面白さと非自明さが現れているのでないかと思う。Euler topは、調和振動子や量子Kepler系のように、Lie環の表現の分岐則を見れば殆どの情報が得られるというような、単純な話にはなっていない。

また、以下のようなRisa/Asirコード書いて計算すると、半整数スピン(Nが偶数)の時は、ハミルトニアンの固有多項式がF(x)^2の形に因数分解できるらしい。つまり、理由はわからないけど、半整数スピンの時は、全ての固有状態が二重に縮退している。不思議なことである

N=5;  /*  2s+1  */

Emat = newmat(N,N);
Fmat = newmat(N,N);
Hmat = newmat(N,N);

Umat = newmat(N,N);

for(I=0;I<N-1;I++){
   Emat[I][I+1] = I+1;
   Fmat[I+1][I] = N-I-1;
   Hmat[I][I] = (N-1-2*I);
   Umat[I][I] = 1;
}
Hmat[N-1][N-1] = -(N-1);
Umat[N-1][N-1] = 1;

S1 = Emat+Fmat;
S2 = @i*(-Emat+Fmat);
S3 = Hmat;

H = a*S1*S1+b*S2*S2+c*S3*S3;

fctr(det(H-4*x*Umat));

(量子力学の入門的な教科書には載っていないという意味で)初等的でない例として、rational Calogero-Moser模型に興味がある(trigonometricなSutherland模型も殆ど同じことになると思うけど)。この模型は、一次元空間上の相互作用する多体粒子系と説明される。表現論的には、粒子が何次元空間に住んでるとか、多体系か否かは副次的な解釈の問題に過ぎないので、どうでもいいが。1970年代はじめに、量子可積分系として発見され、1978年には、Kazhdan-Kostant-Sternbergによって、Calogero-Moser spaceという"真の相空間"が構成されている。量子系の方は、Dunkl作用素という、調和振動子に於ける生成・消滅演算子のようなものを導入して解かれる。Dunkl operatorは1980年代終わり頃発見されたものだけど、90年代以降も、非常に多くの研究があり、多くの一般化や様々な代数との関連が議論されている。

The Quantum Calogero-Moser Model: Algebraic Structures

https://www.jstage.jst.go.jp/article/jpsj1946/62/9/62_9_3035/_article

Integrability and algebraic structure of the quantum Calogero-Moser model

http://www.sciencedirect.com/science/article/pii/0960077993900497

Exact Solvability of the Calogero and Sutherland Models

https://arxiv.org/abs/hep-th/9506105

Lectures on Calogero-Moser systems

https://arxiv.org/abs/math/0606233

Dunkl operators: Theory and applications

https://arxiv.org/abs/math/0210366




SBMLによる生化学ネットワークの記述と計算

SBMLは、Systems Biology Markup Languageの略。シグナル伝達や代謝系などの生化学ネットワークのモデルを形式的に記述するXML仕様で、2000年頃から継続的に策定と改訂が行われているらしい。生物分野で古典的な数理モデルといえば、酵素反応速度論(Michaelis-Menten式etc.)、Hodgkin-Huxley方程式、Turingの反応拡散系などがあるけど、前の2つは、SBMLで記述されるモデルのプロトタイプと言っていいと思う。これらのモデルは、数個の変数しかない単純なモデルだけど、最近では細胞内の生化学ネットワークを調整する遺伝子たちが沢山同定されて、複雑な生化学ネットワークを、そのままモデル化して調べるという道ができた。


その極致として、2012年に、500個の遺伝子を持つ細菌マイコプラズマの全細胞シミュレーションが行われている

A whole-cell computational model predicts phenotype from genotype.

https://www.ncbi.nlm.nih.gov/pubmed/22817898

いずれ、線虫(細胞数1000、遺伝子数20000弱)くらいの生物で全細胞シミュレーションができれば面白いと思うけど、全ての要素をモデル化するだけでも大変そう(マイコプラズマの場合でも、シミュレーションの結果から新しい酵素が見つかったとかあったようだし)


#関連する話題として、最近、RIKENで同じくマイコプラズマの全原子分子動力学シミュレーションを行ったという報告があった

バクテリア細胞質の全原子分子動力学計算(プレスリリース)

http://www.riken.jp/pr/press/2016/20161101_3/

こっちは、20nsを数カ月かけて計算するというものらしいので、計算量的に何かへ応用するというには、難しいけれど。


More Detailed Summary of SBML

http://sbml.org/More_Detailed_Summary_of_SBML

に、もう少し詳しいSBMLの説明と、ミカエリス・メンテン型酵素反応の例が載っている。


既に知られているモデルの場合は、SBMLを

BioModels Database

http://www.ebi.ac.uk/biomodels-main/

とかから見つけてこれるかもしれない。Curatedでないモデルは色々ダメなことが多いっぽいので、使えるモデルは600ちょっとと思っておいたほうがよさげではある。まぁ、この手のシミュレータを実装すること自体は、難しいことでもないと思うけど、何気にパラメータが多くて、論文などから適切な値を探すのが面倒だったりするので、モデルをまとめて公開してくれているとありがたい。


SBMLのシミュレータはいくつも公開されている。

SBML Software Matrix

http://sbml.org/SBML_Software_Guide/SBML_Software_Matrix

にソフトウェアのリストがある。libSBMLとかあるけど、こいつにはシミューレション機能とかはなく、SBMLファイルの読み書きをするだけで、割とどうでもいい。対応しているSBMLのバージョン、実装されている機能や計算速度、継続的にメンテナンスがされているのかなどを考慮して決めたいのだけど、いまいちよく分からないので、自分で実装した方が早そうであった



最も単純な利用法は、常微分方程式系に変換してソルバーに投げるというのだと思う。細胞は非常に小さく、分子種によっては分子数が非常に少ない場合があり、ゆらぎが無視できない場合もあると考えられており、その場合は確率モデルを用いる必要があると考えられている。その他にも、目的に応じて、適用したい解析手法は多数あると思う。

自分でシミュレータを実装する場合は

SBML Test Suite

http://sbml.org/Software/SBML_Test_Suite

から、テストデータやらを引っ張ってこれるっぽいので、確認用に使える。



例として、カルシウム振動のモデル

http://www.ebi.ac.uk/compneur-srv/biomodels-main/BIOMD0000000184

を使う。アストロサイトにおけるカルシウム振動の(多分最も単純な)モデル。変数としては細胞質のカルシウム濃度とIP3(イノシトール3リン酸)濃度と、ERのカルシウムの濃度。

(1)細胞外から細胞内へのCa2+の流入及び流出

(2)ERから細胞質へのカルシウムの放出(IP3受容体を介するもの)

(3)細胞質からERへのカルシウムの浸透

(4)ERから細胞質へのカルシウムの漏れ

(5)IP3の合成(カルシウム濃度依存。ホスホリパーゼCの活性がカルシウムによる調節を受けるらしい)

(6)IP3の分解

が考慮される。3つの変数が振動を起こすことが確認できる。計算結果では、振動周期は3分程度である。これが、カルシウム振動のモデルとして適切なのかどうかは不明


とりあえず、これくらいのモデルがPython(+numpy/scipy/matplotlib)で動くようにした

・多分、2.7と3.5の両方で動く

・Test Suiteは通してないし、仕様の確認もしてないので、仕様を網羅もしてないと思う

・凄く遅い


いうて、常微分方程式やし、簡単にできると思ったが、色々面倒くさいことがわかった。とりあえずミカエリス・メンテン型反応とカルシウム振動の計算結果を貼っておく

f:id:m-a-o:20161117232853p:image:w360

f:id:m-a-o:20161117232851p:image:w360

from scipy.integrate import ode
import xml.etree.ElementTree as et
import math
from numpy import prod


def sbml_root_fn(xs):
    if len(xs)==1:
        return math.sqrt(xs[0])
    elif len(xs)==2:
        return math.pow(xs[1] , 1.0/xs[0])
    else:
        assert(False),"invalid argument in sbml_root_fn:{}".format(str(xs))

class Closure:
    def __init__(self , _vars , fnbody):
        self.vars = _vars
        self.fndef = fnbody
        self.env = {}
    def __call__(self , args):
        assert(len(args)==len(self.vars)),"invalid closure call"
        saved_env = {}
        for (key,val) in zip(self.vars , args):
            if key in self.env:
                saved_env[key] = self.env[key]
            self.env[key] = val
        ret = eval_math(self.fndef , self.env)
        for key in self.vars:
            del self.env[key]
        for (key,val) in saved_env.items():
            self.env[key] = val
        return ret


def parse_math(e):
    def parse_cn(e):
       assert(e.attrib.get('type',None)!="rational"),"rational cn@parse_cn"
       return eval(e.text)
    def parse_apply(es):
       fn = None
       args = []
       for e0 in es:
          _ , tagname  = e0.tag.split('}')
          if tagname=="plus":
              fn = (lambda x:x[0]+x[1])
          elif tagname=="minus":
              fn = (lambda x:x[0]-x[1])
          elif tagname=="divide":
              fn = (lambda x:x[0]/x[1])
          elif tagname=="times":
              fn = (lambda x:prod(x))
          elif tagname=="power":
              fn = (lambda x:math.pow(x[0],x[1]))
	  elif tagname=="exp":
	      fn = (lambda x:math.exp(x[0]))
          elif tagname=="root":
              fn = sbml_root_fn
          elif tagname=="apply":
              args.append( parse_math(e0) )
          elif tagname=="degree":
              assert(fn==sbml_root_fn)
              assert(len(e0.getchildren())==1),"error in parse_apply@degree"
              e1 = e0.getchildren()[0]
              _ , tagname1  = e1.tag.split('}')
              if tagname1=="ci":
                 args.append( e1.text.strip() )
              elif tagname1=="cn":
                 assert(e0.attrib.get('type',None)!="rational"),"rational cn@parse_math"
                 args.append( eval(e1.text) )
              else:
                 assert(False),"failed to parse degree"
          elif tagname=="ci":
              if fn!=None:
                  args.append( e0.text.strip() )
              else:
                  fn = e0.text.strip()
          elif tagname=="cn":
              assert(e0.attrib.get('type',None)!="rational"),"rational cn@parse_math"
              args.append( eval(e0.text) )
          else:
              assert(False),"@parse_math:unknown tag {}".format(tagname)
       assert(fn!=None),et.tostring(es)
       return (fn,args)
    def parse_lambda(es):
       varnames = []
       for e0 in es:
          _ , tagname  = e0.tag.split('}')
          if tagname=="bvar":
              for it in e0.getchildren():
                  varnames.append( it.text.strip() )
          elif tagname=="apply":
              fnbody = parse_apply(e0.getchildren())
              return Closure(varnames , fnbody)
          else:
              assert(False),"unknown tag name@parse_lambda:{}".format(tagname)
    _ , tagname  = e.tag.split('}')
    if tagname=="math":
         assert(len(e.getchildren())==1),"error in parse_math"
         return parse_math(e.getchildren()[0])
    elif tagname=="lambda":
         return parse_lambda(e.getchildren())
    elif tagname=="apply":
         return parse_apply(e.getchildren())
    elif tagname=="cn":
         return parse_cn(e)
    else:
         assert(False),"unknown tagname@parse_math:{}".format(tagname)


def eval_math(e , env):
   if type(e)==tuple:
      fn,_args = e
      args = [eval_math(x,env) for x in _args]
      if type(fn)!=str:
          return fn(args)
      else:
          closure = env[fn]
          closure.env = env
          return closure(args)
   elif type(e)==str:
      if e in env:
          return env[e]
      else:
          assert(False),"Unknown variable in env:{}".format(e)
   else:
      return e



class Model:
   def __init__(self):
      self.args = []
      self.species = []
      self.namemap = {}
      self.env = {}
      self.reactions = []
      self.compartment = {}
      self.rules = {}
   def values(self):
      ret = []
      for key in self.args:
          v = self.env[key]
	  if type(v)==float or type(v)==int:
             ret.append( self.env[key] )
          else:
	     ret.append( eval_math(v ,self.env) )
      return ret
   def __call__(self , t , y):
      assert(len(y)==len(self.args)),"invalid input"
      ret = [0.0 for _ in y]
      for (name,val) in zip(self.args , y):
          self.env[name] = val
      rest = self.rules.keys()
      Nmax = len(rest)
      #-- process assignment rules
      for _ in range(Nmax):
          if len(rest)==0:break
          for (key,sm) in self.rules.items():
	      if key in rest:
	         try:
	            self.env[key] = eval_math(sm , self.env)
		    rest.remove(key)
		 except:
		    pass
      #-- process reactions and rate rules
      for (reactants , products , sm , local_env) in self.reactions:
          saved_env = {}
          for key,val in local_env.items():
              if key in self.env:
                 saved_env[key] = self.env[key]
              if type(val)==str:
                 self.env[key] = self.env[val]
              else:
                 self.env[key] = val
          try:
              v = eval_math(sm , self.env)
          except:
              assert(False),("failed to evaluate math expresion",sm)
          for r in reactants:
	      try:
                 idx = self.args.index(r)
	         if r in self.compartment:
                     ret[idx] -= v/self.env[self.compartment[r]]
	         else:
	             ret[idx] -= v
	      except:
	         pass
          for r in products:
	      try:
                 idx = self.args.index(r)
	         if r in self.compartment:
                     ret[idx] += v/self.env[self.compartment[r]]
	         else:
	             ret[idx] += v
	      except:
	         pass
          for key in local_env:
              del self.env[key]
          for key,val in saved_env.items():
              self.env[key] = val
      return ret


def readsbml(root):
   ns = {'sbml' : root.tag[1:-5] , 'mml':"http://www.w3.org/1998/Math/MathML"}
   ret = Model()
   for it in root.find('.//sbml:listOfCompartments' , namespaces=ns).getchildren():
       key,value = it.attrib['id'] ,float( it.attrib['size'] )
       ret.env[key] = value
   for it in root.find('.//sbml:listOfSpecies' , namespaces=ns).getchildren():
       if it.attrib.get('boundaryCondition' , 'false')!='true':
           ret.args.append( it.attrib['id'] )
       if 'name' in it.attrib:
           ret.species.append( it.attrib['name'] )
	   ret.namemap[ it.attrib['id'] ] = it.attrib['name']
       else:
           ret.species.append( it.attrib['id'] )
	   ret.namemap[ it.attrib['id'] ] = it.attrib['id']
       ret.compartment[ it.attrib['id'] ] = it.attrib['compartment']
       if "initialAmount" in it.attrib:
            vol = ret.env[it.attrib["compartment"]]
            key,val = it.attrib['id'] , float(it.attrib['initialAmount'])/vol
            ret.env[key] = val
       elif "initialConcentration" in it.attrib:
            vol = ret.env[it.attrib["compartment"]]
            key,val = it.attrib['id'] ,  float(it.attrib['initialConcentration'])
            ret.env[key] = val
       else:
            pass
   for it in root.findall('.//sbml:listOfParameters' , namespaces=ns):
       for r in it.getchildren():
           if 'value' in r.attrib:
               ret.env[ r.attrib['id'] ] = eval(r.attrib['value'])
   for it in root.findall('.//sbml:functionDefinition' , namespaces=ns):
       fndef = it.find('.//mml:math' , namespaces=ns)
       key = it.attrib['id']
       if fndef!=None:
           fn = parse_math(fndef)
           ret.env[key] = fn
   listOfInitialAssignments = root.find('.//sbml:listOfInitialAssignments' ,  namespaces=ns)
   if listOfInitialAssignments!=None:
       for it in listOfInitialAssignments.getchildren():
           varname = it.attrib['symbol']
	   ret.env[varname] = eval_math(parse_math(it.getchildren()[0]) , ret.env)
   listOfRules = root.find('.//sbml:listOfRules' , namespaces=ns)
   if listOfRules!=None:
       for it in listOfRules.getchildren():
           _ , tagname  = it.tag.split('}')
           key = it.attrib['variable']
	   if tagname=="assignmentRule":
	       rule = parse_math(it.getchildren()[0])
	       ret.rules[key] = rule
	   else:
	       ret.args.append( key )
	       ret.reactions.append( ([] , [key] , parse_math(it.getchildren()[0]) ,{}) )
   for r in root.find('.//sbml:listOfReactions' , namespaces=ns).getchildren():
       reactants = []
       products = []
       local_env = {}
       listOfReactants = r.find('.//sbml:listOfReactants' , namespaces=ns)
       if listOfReactants!=None:
           for it in listOfReactants.getchildren():
               reactants.append( it.attrib['species'] )
               if 'id' in it.attrib:
                    local_env[ it.attrib['id'] ] = it.attrib['species']
       listOfProducts = r.find('.//sbml:listOfProducts' , namespaces=ns)
       if listOfProducts!=None:
           for it in listOfProducts.getchildren():
               products.append( it.attrib['species'] )
               if 'id' in it.attrib:
                    local_env[ it.attrib['id'] ] = it.attrib['species']
       listOfModifiers = r.find('.//sbml:listOfModifiers' , namespaces=ns)
       if listOfModifiers!=None:
           for it in listOfModifiers.getchildren():
               if 'id' in it.attrib:
                    local_env[ it.attrib['id'] ] = it.attrib['species']
       listOfLocalParameters = r.find('.//sbml:listOfLocalParameters' , namespaces=ns)
       if listOfLocalParameters!=None:
           for it in listOfLocalParameters.getchildren():
               local_env[ it.attrib['id'] ] = eval(it.attrib['value'])
       sm = r.find('.//mml:math' , namespaces=ns)
       ret.reactions.append( (reactants , products , parse_math(sm) , local_env) )
   return ret


try:  #-- python2
   from urllib2 import urlopen
except: #-- python3
   from urllib.request import urlopen


import matplotlib.pyplot as plt


def runsbml(path , t0 , t1, dt):
   try:
      root = et.parse(path).getroot()
   except:
      req = urlopen(path)
      s = req.read()
      root = et.fromstring(s)
   m = readsbml(root)
   r = ode(m, None).set_integrator('vode', method='bdf', nsteps=1500 , with_jacobian=False)
   r.set_initial_value(m.values() , t0)
   res = []
   times = []
   while r.successful() and r.t < t1:
      _ = r.integrate(r.t+dt)
      times.append( r.t )
      res.append( r.y )
   for n,name in enumerate(m.args):
      dat = [x[n] for x in res]
      plt.plot(times , dat , label=m.namemap[name])
   plt.legend()
   plt.show()


if __name__=="__main__":
 #-- calcium oscillation
  runsbml("http://www.ebi.ac.uk/compneur-srv/biomodels-main/download?mid=BIOMD0000000184" , 0.0 , 2000.0 , 0.5)

[広告]

コメント
0件
トラックバック
0件
ブックマーク
0 users

[広告]

m-a-o
m-a-o