2017-08-12

最近の老化研究

2017年3月の『Cell』の論文

Targeted Apoptosis of Senescent Cells Restores Tissue Homeostasis in Response to Chemotoxicity and Aging

http://www.cell.com/fulltext/S0092-8674(17)30246-5

大雑把には、(マウスの)"老化細胞"(Senescent cells)でアポトーシスを誘導してやって除去すると、若返るよという話(例えば、Fig7(D)では、毛が明らかに増量している)。ここ数年、老化細胞の除去による"若返り"の報告はいくつかあるが、いずれも老化細胞の除去に遺伝子改変を必要とするものばかりだったのに対して、この論文の方法では、遺伝子改変を必要としない(マウスの実験であるから、ヒトなどでは、別途検討が必要ではあるが)

cf)細胞老化が生体機能に及ぼす影響とその機構に関する研究

http://www.ncgg.go.jp/ncgg-kenkyu/documents/25-18.pdf

http://www.ncgg.go.jp/research/news/20160810.html

肺組織で特異的に老化細胞を除去できるように、遺伝子改変を行ったマウスでの類似実験を行っている

cf)Naturally occurring p16Ink4a-positive cells shorten healthy lifespan

https://www.nature.com/nature/journal/v530/n7589/full/nature16932.html


背景。正常な細胞は、老化刺激(テロメア短小化、DNA損傷、酸化ストレスetc)を受けると、恒久的に細胞周期を停止した細胞老化と呼ばれる状態になることが知られている。こうしてできた老化細胞は、増殖を再開することはないと考えられているが、ただ機能を失って存在するだけでなく、炎症性サイトカイン、増殖因子、細胞外マトリックス分解酵素など、様々な分泌因子を分泌するSASPと呼ばれる現象を起こし、周辺組織に炎症状態を誘発する。SASPは、周囲の細胞で細胞老化や癌を誘導するなど生体にとって有害な作用をもたらすと現在考えられている。

cf)老化した細胞ががん化を促進する仕組みをハエで解明(プレスリリース)

http://research.kyoto-u.ac.jp/research/141011_1/


論文では、老化細胞でFOXO4の発現が亢進することを確認し(Fig1(D))、また、FOXO4を阻害すると、老化細胞の生存率が下がる(Fig(1(H))ことから、FOXO4が老化細胞のアポトーシスを抑制し、老化細胞の生存・維持に寄与していると推測し、FOXO4のp53結合部位を含む部分ペプチドのD-retro inverso (DRI)-isoform(FOXO4-DRI)を設計し、マウスに腹腔投与(i.p. injection)したという感じ(5mg/kgを一日一回x3回)。

#D-retro inverso-isoformというのは初耳だったのだけど、配列とアミノ酸のキラリティ(つまりD型アミノ酸から構成される)が逆になってるような異性体のことらしい。なんだかよく分からないけど、DRI-isoformはしばしば、元のペプチドと同じ活性を保持し、かつ、安定であるらしい(多分プロテアーゼなどで分解されにくい?)。今の場合は、p53とFOXO4の結合を競合阻害するのが目的だろうと思う(部分ペプチドなので、p53を抑制する機能はない)


投与を行った"老いたマウス"の年齢は104〜110週齢程度のよう(Fig7 captionなど)。20ヶ月齢のマウスが人間でいうと60歳相当程度らしいので、人間換算で30歳弱と若干若い感じがする。マウスを直接扱ったことがないので、何とも言えないけど、老化研究では、もう少し高齢のマウスを対象にすることが多いように思う。一応、Fig6(B)の腎臓組織のSA-β-GAL染色では、130week miceでも、結構、老化細胞が増えているらしきことが見て取れるし、最初の方のcfであげた

細胞老化が生体機能に及ぼす影響とその機構に関する研究

http://www.ncgg.go.jp/ncgg-kenkyu/documents/25-18.pdf

でも、12ヶ月齢のマウスが使われている(これは、肺組織に於ける老化細胞の除去実験)。


毒性について。FOXO4は、FOXO1とFOXO3に比べると、注目されることが少なく、非老化細胞ではあまり重要な役割を持っていないらしい。実際にFOXO4をノックアウトしたマウスは、目立った表現型を示さない("FOXO4 knockout mice do not show a striking phenotype")らしいので、FOXO4の阻害は、非老化細胞に対する影響は少ないと予想される。


[疑問]FOXO4の阻害が老化細胞のアポトーシスを引き起こすなら、FOXO4ノックアウトマウスでは、老化細胞の数が少なく(論文の結果が正しければ)見た目の老化も抑制されそうに思えるが、目立った表現型は示さないという報告とは矛盾する気もする。老化が遅いなどの表現型は注意して観察しなければ気付かない微妙なもの(致死とかに比べればstriking phenotypeでない)で単に見落とされているという可能性もあるけど


[補足]癌抑制遺伝子として有名な転写因子p53は、細胞周期を停止し細胞増殖の抑制を行うだけでなく、細胞老化を引き起こし、アポトーシスを誘導することも知られている。これらの機能は、独立しているわけでなく、細胞やDNAがダメージを受けた時に細胞周期を停止し修復を試みると同時に、ダメージが重度で(癌化しそう)であればアポトーシスを起こし、中度であれば細胞老化を引き起こすのだろうと思われる(回復可能な程度に軽度であれば、細胞周期を再開する)。というわけで細胞老化を誘導する薬剤を高濃度で与えると、アポトーシスを誘導するっぽい。培養細胞に於いて、p53の発現を亢進すれば、癌は減るが老化は進行し(寿命も縮むらしい)、一方、抑制すれば、老化は抑えられるが癌は増えるらしい(マウス個体でp53をノックアウトすると早期に死亡するので解析は難しいが、培養細胞での実験と矛盾しない結果が知られている)。このような癌化と細胞老化のトレードオフは、テロメラーゼの高発現/ノックアウトでも観察することができる(マウスはテロメアが長く、テロメラーゼノックアウトマウスで個体老化が観察できるようになるまで数世代かかる)らしい

cf)癌抑制遺伝子p53 による老化の制御

http://www.pharmacol.or.jp/fpj/topic/topic121-2-130.htm


#細胞周期の一時停止、細胞老化の誘導、アポトーシスの選択がどのようにされているのか詳細は不明。p53には、リン酸化部位が沢山あり、中でもSer46がリン酸化される(このリン酸化を行う酵素はHIPK2,ATK kinase,DYRK2など複数ある?)と、アポトーシス関連遺伝子(AIP1という遺伝子が知られる)の発現が誘導されるらしい。一方、Ser15のリン酸化(AMPKなどによる?)は細胞老化関連遺伝子の発現と関連するらしい。DNA損傷が起きた時、Ser46のリン酸化はSer15のリン酸化に比べると、少し後になってから起こるらしい

The p53 response to DNA damage.

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

Ser46 phosphorylation of p53 is not always sufficient to induce apoptosis: multiple mechanisms of regulation of p53-dependent apoptosis.

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

Palmdelphin, a novel target of p53 with Ser46 phosphorylation, controls cell death in response to DNA damage

http://www.nature.com/cddis/journal/v5/n5/full/cddis2014176a.html?foxtrotcallback=true


素朴な疑問として、何のために細胞老化などという機構があるのか。癌抑制のためだけであれば、とっととアポトーシスを起こすなりして除去されてもよさそうなものだけど、わざわざアポトーシスを抑制してまで生存を図る意味が分からない。最近では、創傷治癒に於いて細胞老化の果たす役割が報告されたりしていて、他にもまだ知られてない生物学的機能が存在する可能性はある

cf)Cellular senescence controls fibrosis in wound healing

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



それらの話とはまた別に組織幹細胞の減少と機能低下も個体老化の一因として挙げられることがある。細胞老化と幹細胞は、共に創傷治癒に寄与する、という共通点があるっぽい。7月のNatureに

Hypothalamic stem cells control ageing speed partly through exosomal miRNAs

http://www.nature.com/nature/journal/v548/n7665/full/nature23282.html

という論文が出ていた。これは、視床下部神経幹細胞(htNSC)の減少が老化のトリガーとなっており、特に、exosomal miRNAの放出量が減ってしまうことが主要因なのでないかとのこと(生物種はマウス)。

ちなみに、論文中では"mid-aged mice (11–16 months old)" "aged mice (22 months and older)"と書かれていて、やっぱCell論文の110wk miceとかは中年手前なんだな、という感じ。こっちの論文のFig1(a)などを見ると、16 months-old miceでも、視床下部幹細胞の減少が見て取れるし、高齢マウスでは殆ど消失しかかっている。Cell論文に比べると、視床下部幹細胞の注入("炎症"のせいか、単純に移植しても、すぐに死んでしまうので、ドミナント・ネガティブIκB-αを発現する細胞を作ったと書いてある)などにより、老化の進行が抑制されて寿命が延長してる(16ヶ月齢雄マウスに移植して4ヶ月後の生理機能(Fig3(c))や18ヶ月齢雄マウスに移植して寿命計測(Fig3(d))


この論文では、視床下部幹細胞に注目しているが、

Muscle-derived stem/progenitor cell dysfunction limits healthspan and lifespan in a murine progeria model

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3272577/

のような論文も以前に出ている(タイトルを直訳すると、『筋由来幹/前駆細胞の機能不全が、早老症モデルマウスの健康寿命と寿命を制限する』)。若いマウスのmuscle-derived stem/progenitor cellを移植することにより、早老症モデルマウスの寿命を大幅に延長できたようである(Fig4(a))。通常の老齢マウスに対して、同様の実験を行って寿命延長などが起こるかは書いていない。まぁ、視床下部幹細胞が特別というわけでもなさそうではある。また、移植後、比較的短期間で効果が出ていることから、移植した幹細胞そのものが分化しているわけではなく、何らかの分泌因子が重要らしい(が、分泌因子の正体は不明)というところで、終わっている。視床下部幹細胞の方も事情は似ていて、移植した細胞が直接分化するわけではなく、こちらは、何らかのmiRNAが重要らしいということになっている

物凄く単純化していうと、老化細胞除去は、"悪いもの"を減らすという考え方で、一方、若い幹細胞を移植するというのは、"良い物"を増やすという考え方なので、割と真逆ではある。もし後者が重要であるならば、単純に老化細胞を除去しても、機能低下した幹細胞が回復するとは限らず、高齢になってからの老化細胞除去は、効果が薄いか、手遅れである可能性もある。老化細胞除去が比較的若いマウスに対して行われている理由かもしれないけど、どうだろう。

[] 体積粘性係数の分子動力学計算

等方的な流体の場合、粘性係数には、shear viscosityとbulk viscosity(体積粘性係数)の2種類がある。水のような液体の場合、非圧縮性近似をするのが一般的なので、Navier-Stokes方程式で体積粘性係数がかかる項は0になる。気体の場合も、流速が音速よりずっと遅い場合には非圧縮近似はよい近似になっているとされている。固体の場合は、弾性が支配的なので、通常、粘性そのものが考慮されない。


Wikipedia情報によれば、体積粘性係数の値は多くの場合、あんまり知られてないらしい

Volume viscosity/Measurement

https://en.wikipedia.org/wiki/Volume_viscosity#Measurement

The volume viscosity of many fluids is inaccurately known, despite its fundamental role for fluid dynamics at high frequencies. The only values for the volume viscosity of simple Newtonian liquids known to us come from the old Litovitz and Davis review, see References. They report the volume viscosity of water at 15 °C is 3.09 centipoise.


そもそも、通常の流体運動の範疇では体積粘性係数項は消えてしまって観測量と関わってこない(というか、NS方程式から消えてるだけだけど)から、測定が難しいのは仕方ない。多分、体積粘性係数を考慮に入れるべき身近な現象は音の減衰だと思う。なので、体積粘性係数を測定する方法として、音を利用するのは自然に思える。

Bulk viscosity and compressibility measurement using acoustic spectroscopy

http://aip.scitation.org/doi/abs/10.1063/1.3095471?journalCode=jcp

という論文には"To date, measurement of ultrasound attenuation is the only known way of experimentally studying bulk viscosity."と書かれている(Section Bの冒頭)。この論文によれば、25度に於ける水の体積粘性係数は、2.4cPであるらしい(TABLE III)

医療分野では、流体ではない生体組織の弾性や粘性(の分布)を測定するエラストグラフィと総称される手法がいくつか存在し、知る限り、いずれの方法でも超音波を何らかの形で利用している。

Elastography

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

elastographyという言葉は、名前の通り、弾性の計測を主要な目的としている。癌などの腫瘍組織は一般に周囲の組織より硬くなる(乳がんが、しこりとして検出できる理由でもある)ことから、病変検出に有用と考えられている。一方、音の伝播の方程式に粘性項を入れるのは、固体であっても流体と変わるところはないので、同じ方法で粘性係数を決めることも理屈上はできる。実際に、そのような研究は存在するものの、弾性に比べると粘性を決定する医学的意義はよく分からないので、広く行われているとは言いがたいっぽい。しかし、何はともあれ、固体・流体を問わず、実験的に体積粘性係数を知る方法は、超音波の伝播や減衰を調べる以外にないらしい



現在は21世紀なので、簡単な分子の場合は、実験によらず、分子動力学で計算するという道がありえる。原理的には、shear viscosityも体積粘性係数も、Green-Kubo公式で計算できるはずなので。shear viscosityの計算は色々やられていて、手法もいくつか存在する

[1]Determining the Bulk Viscosity of Rigid Water Models

http://pubs.acs.org/doi/pdf/10.1021/jp211952y

は、分子動力学で、shear viscosityとbulk viscosityを決める論文。水のモデルとしては、TIP4P/2005がよいらしい。これは、他の論文でも、精度高いと書いてある


というわけで、Gromacsを使って、上の論文と同じことができないか試した。Gromacsのバージョンは5.1.2(aptで入れたら、これになった)。マニュアルには、gmx viscosityコマンドで、上の論文の手法が実装されてるよって書いてあるけど、そんなコマンドは存在しない。というわけで、手動で頑張ってデータを解析するしかない

また、水のモデルは、デフォルトでは、TIP4P/2005が入ってない

GROMACS files for the TIP4P/2005 model

http://www.sklogwiki.org/SklogWiki/index.php/GROMACS_files_for_the_TIP4P/2005_model

にあるのが、TIP4P/2005モデルらしいけど、面倒なので、普通のTIP4Pを使うことにした。



論文[1]によれば、

At the beginning, a MD simulation was performed at constant temperature and pressure conditions (NPT ensemble) for a period of 2 ns, in order to determine the density of liquid water as it is predicted by each water model.

The average energy ⟨E⟩ of the system at that density was computed from a second constant volume and temperature (NVT ensemble) simulation of 2 ns.

For the calculation of the bulk and shear viscosity at each temperature, we have performed a set of 100 independent constant volume and energy (NVE ensemble) simulations. For each trajectory, initial atomic momenta were sampled from a Maxwell distribution at temperature T.

らしいので、NPTアンサンブルで2ns→NVTアンサンブルで2ns→独立なNVE計算を100回(10psの平衡化の後、300psの計算をやると書いてある)という順番で計算するらしい。NVE計算では、長い時間やると、enerygy driftによって、エネルギー保存が成立しなくなるので短時間の計算を沢山やるということのよう。平衡化する場合、NVT→NPTの順が一般的な気がするけど、今回は逆


論文では、水256分子で計算しているので、適当に250分子を用意した。まあ、大体、以下のような感じのコマンドを叩く。エネルギー最小化は論文には書いてないので、特にやらなくてもいいかもしれない。NVE計算は、面倒なので100回やらずに30回しかやってない。

#-- We will get 250 molecules
gmx solvate -cs tip4p -o water.gro -box 2.1 2.0 2.0 -p conf/topol.top

#-- Energy minimization
gmx grompp -f conf/em1.mdp -o em1 -c water.gro -p conf/topol.top
gmx mdrun -s em1.tpr -deffnm em1 -v 
gmx grompp -f conf/em2.mdp -o em2 -pp em2 -po em2 -c em1 -t em1 -p conf/topol.top
gmx mdrun -s em2.tpr -deffnm em2

##gmx energy -f em1.edr -o min-energy.xvg
##gmx energy -f em2.edr -o min2-energy.xvg 

#-- NPT ensemble (to determine the density of liquid water)
gmx grompp -f conf/npt.mdp -o nptin -c em2.gro -p conf/topol.top -maxwarn 1
gmx mdrun -s nptin.tpr -deffnm nptout

#-- check density ,pressure, energy, temparature
echo -e "7\n8\n10\n14\n15\n" | gmx energy -f nptout.edr -o nptout.xvg

#-- NVT ensemble (to determine the average energy)
gmx grompp -f conf/nvt.mdp -o nvtin -c nptout -t nptout -p conf/topol.top
gmx mdrun -s nvtin.tpr -deffnm nvtout

#-- check energy
echo -e "7\n9\n11\n" | gmx energy -f nvtout.edr -o nvtout.xvg

#-- NVE ensemble
for i in `seq 1 30`;do
  gmx grompp -f conf/nve.mdp -o nvein -c nvtout -t nvtout -p conf/topol.top
  gmx mdrun -s nvein.tpr -deffnm nveout_$i
  echo -e "7\n8\n10\n" | gmx energy -f nveout_$i.edr -o nveout_$i.xvg    #--energy,temperature,pressure
done

使用したmdpファイルとtopol.topは、以下に置いた

https://gist.github.com/vertexoperator/aebb9676509e5eba608dfff6ff74bfc5




NPT計算の結果

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Total Energy               -9772.19        4.7    154.428    -16.766  (kJ/mol)
Temperature                 298.147      0.062    11.6618   -0.23701  (K)
Pressure                   -4.63437          2    868.369   -8.01756  (bar)
Volume                      7.51647     0.0023   0.124593 -0.00247449  (nm^3)
Density                     995.296        0.3    16.3923   0.194343  (kg/m^3)

とかなった。密度が大体正しそうなので多分OK。温度は298.15K=25度を指定しているので、そのへんで動く。論文[1]のTable Iでは、25度でTIP4Pを使用した時、0.9953(g/cm^3)という値を得たと書いてあるので、妥当な数値っぽい


NVT計算の結果は

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Total Energy               -9755.22        3.2    134.505    -16.928  (kJ/mol)
Temperature                 298.143       0.11    11.1177   -0.10815  (K)
Pressure                   -386.763        9.7    777.601    18.4497  (bar)

まぁ、よさげ?

Gromacsには、tcafコマンドというのがあって、私が理解してない謎の手法で、shear viscosityを計算できるらしい。

gmx tcaf

http://manual.gromacs.org/programs/gmx-tcaf.html


それで、NVT計算の結果とかに対して、

gmx tcaf -f nvtout.trr -s nvtin.tpr

とかやると、visc_k.xvgに(一部ヘッダ省略)

@    title "Fits"
@    xaxis  label "k (nm\S-1\N)"
@    yaxis  label "\xh\f{} (10\S-3\N kg m\S-1\N s\S-1\N)"
@TYPE xy
@    s0 symbol 2
@    s0 symbol color 1
@    s0 linestyle 0
 3.113 0.51114
 3.269 0.493288
 3.269 0.48944
 4.514 0.386852
 4.514 0.387317
 4.514 0.398109
 4.514 0.389984
 4.623 0.392635
 4.623 0.385249
 5.573 0.329574
 5.573 0.330219
 5.573 0.326714
 5.573 0.325943
 6.226 0.296062
 6.538 0.28417
 6.538 0.284974

などと出力される(デフォルトでは、visc_k.xvgにkとetaのみの出力がされる)。マニュアルによると、これを適当にfittingすると、shear viscosityが計算できるらしい

以下のようなコードで、実際にfittingしてみた

import numpy as np
from scipy import stats

k,y=[],[]
for line in open("visc_k.xvg"):
    if line[0]=="@":continue
    if line[0]=="#":continue
    ls = line.strip().split()
    k.append(float(ls[0]))
    y.append(float(ls[1]))


x=np.array(k)**2
y=np.array(y)
r=stats.linregress(x,y)
print(r.intercept)

計算すると、例えば0.540536022685という値を得たので、0.54e-3(kg/(m s))という感じであるらしい。水の粘性係数は、25度で0.000894(Pa s)らしいので、2/3くらいの値になっているが、論文[1]の方法では、TIP4Pで計算した場合、25度で、(4.83±0.09)e-4 (Pa s)という値になったらしいし、まぁTIP4Pの粘性に関する精度はこんなものであるらしい。NPT,NVEの結果から、tcafで計算しても、同じような結果を得た。TIP4P/2005なら、もっとよい結果が得られるらしい


Green-Kubo公式で、体積粘性係数を計算するために

B_{ACF}(t) = ¥langle ¥delta P(t) ¥delta P(0) ¥rangle

を計算する。

After an equilibration period of 10 ps, a configuration with total energy equal to the ⟨E⟩ decided before was selected. For that configuration and for the next 300 ps, the instantaneous stress tensor elements were stored in the disk every 20 fs for further analysis.

とか書いてある(<E>はNVT計算の結果を使うらしい)けど、<E>に到達するのに、ものすごく時間がかかるっぽかったので無視した(なんか勘違いしてる?)


よくやるように、

B_{ACF}{(t) = ¥langle ¥lim_{T ¥to ¥infty} ¥frac{1}{T} ¥int_{0}^{T} dt’ ¥delta P(t+t’) ¥delta P(t’) ¥rangle

で計算する。これが得られたら、論文の式(4)でfittingを行う。次のpythonコードが実装例

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt


def func(t , C , tau_f , tau_s , omega , beta_f , beta_s):
   fast_term = (1-C)*np.exp(-np.power((t/tau_f) , beta_f))*np.cos(omega*t)
   slow_term = C*np.exp(-np.power(t/tau_s , beta_s))
   return (fast_term + slow_term)


BACF = []
for i in range(1,31):
   pressures = []
   for line in open("nveout_{}.xvg".format(i)):
        if line[0]=="@" or line[0]=="#":continue
        ls  = line.strip().split()
        assert(len(ls)==4)
        t,E,_,P = ls
        pressures.append( float(P) )
   Pbar = sum(pressures)/len(pressures)
   deltaP = np.array([P-Pbar for P in pressures])
   ACF = []  #--auto-correlation function
   for t in range(len(deltaP)):
      if t==0:
          sv = np.sum(deltaP*deltaP)/(len(deltaP))
      else:
          sv = np.sum(deltaP[t:]*deltaP[0:-t])/(len(deltaP) - t)
      ACF.append( sv )
   if len(BACF)==0:
      BACF = np.copy(ACF)
   else:
      BACF = BACF+ACF


dt = 0.02  #-- pico seconds
xdata = np.arange(len(BACF))*dt
ydata = BACF/30

popt , pcov = curve_fit(func , xdata , ydata/ydata[0])
plt.plot(xdata[:1000], func(xdata[:1000], *popt))
plt.show()

#--積分値(total relaxation time)
print( np.sum(func(xdata , *popt)) * dt)

あと、単位に注意する(粘性係数の単位をPa sとなるように合わせる)。上のdtはpsなので、10^{-12}(s)で、Gromacsが出力する圧力の単位は、bar(1bar = 10^5 Ps)。体積はNVT計算の平均値、温度も同様にして、計算する。上のコードで得られた諸々の数値をもとに、

(7.51647e-27/(1.38e-23 * 298.15))*(ydata[0]*1.0e10)*np.sum(func(xdata , *popt)) * dt*1.0e-12

で、求めたい体積粘性係数の値が出る。私の場合は、0.0018268112792585417という値を得た。論文のTABLE2によれば、298KでTIP4Pを使った場合、0.001173という値を得ている。オーダーは合っている。ずれは、アンサンブル数の少なさ、方法の若干の違いによるものと思われる。最初のほうで見た2.4cPという実測値と比較しても、悪くない。論文を信じれば、TIP4P/2005の方が精度はよいらしい。おしまい

Boltzmannの壺とEhrenfestの壺

最近読んだ論文の読書感想文コーナー

[1]Statistical mechanics of money

https://arxiv.org/abs/cond-mat/0001432


N個の壺とE個のボールがあって、ボールはどれかの壺に入っているとする(論文には、ボールではなく、moneyと書いてあるし、壺とは書いてないけど、確率論に於いて、壺とボールはポピュラーなモデルなので壺にする)。(少なくとも一個以上のボールが入っている)壺をランダムに一つ選んで、やはりランダムに選んだ別の壺に移すということを繰り返すと、壺に入っているボールの平均数は当然E/Nだけど、Nがある程度大きい時、壺に入っているボールの数xの従う確率分布は、C exp(-x/T)で近似できる(Tは平均がE/Nになるように決める。Cは規格化定数)。ボール数をエネルギーと思い、壺と壺の間のボールの交換は、粒子と粒子の間のエネルギー交換の単純化と思うと、C exp(-x/T)はBoltzmann-Gibbs分布と見なせるというようなことが論文に書いてある。


統計力学のボルツマンの原理の類似を信じると、ボールがi個入っている壺の個数をn_iとした時

N=¥sum_{i=1}^B n_i

E=¥sum_{j=1}^B j n_j

という条件下で、

W = ¥frac{N!}{n_1! ¥cdots n_c!}

を最大化するような[n_1,...n_E]が、平衡状態では実現される。この条件から、Nが大きい時は、

 n_j/N ¥propto e^{-¥beta j}

であることが示せるというのは統計力学の教科書に書いてある通りで、dynamicsの話がなければ、まぁそういうもんかという感じであるものの、[1]の論文のように時間発展の仕方が定まれば、定常分布で巨視的状態[n_1,...n_E]が実現される確率は(原理的には)厳密に計算できる量となる(NとEが有限であれば、ただの有限Markov連鎖であるので)。その場合でも、(熱力学極限では)平衡状態が教科書通りに実現されるのかは自明でない気がする。で、[1]の論文にはこれが正しいことの説明があるが、数学の論文でなく、厳密な証明にはなってないと思うし、厳密な証明が知られているのかどうかは不明


このモデルを最初に考えたのは誰かと思ったら

非結合的代数によるランダムな衝突モデルの表現

https://ismrepo.ism.ac.jp/?action=pages_view_main&active_action=repository_view_main_item_detail&item_id=32615&item_no=1&page_id=13&block_id=21

によると、Boltzmann自身が、ほぼ同じモデルを考えていたっぽい。ただ、Boltzmannの結論だと¥sqrt{x} exp(-x/T)が出るらしい。ミクロカノニカル分布を仮定しても、この形は出ない。これは、本文中の条件(1.2)を仮定したためだと思うのだけど、ボールを一個ずつ交換する単純なモデルでは、これが成立する根拠はない。代わりに

A_{¥kappa ¥lambda}^{kl} = A_{kl}^{¥kappa ¥lambda}

を仮定して、Boltzmannと同様の議論を行えば、C exp(-x/T)は出る。これは良さそうな気もするけど、そもそも、このような計算が厳密には正しくないので微妙。例えば、N=4,E=10として、4つの壺に入っているボールの数(微視的状態)を(n0,n1,n2,n3)で表すとして、(1,2,7,0)->(0,3,7,0)と(1,2,3,4)->(0,3,3,4)の2種類の遷移を考える。1つめと2つ目の壺では、いずれも(1,2)->(0,3)という遷移が起きているけども、(1,2,7,0)から(0,3,7,0)へ遷移する確率と(1,2,3,4)から(0,3,3,4)へ遷移する確率は等しくない。従って、上のAがきちんと決まらない。


#元の論文では、ボールの数をmoneyと解釈しているが、現実の取引のモデルとしての妥当性は疑問(例えば、このモデルでは一旦金持ちとなっても、比較的容易に没落する)



一応、簡単な数値実験。

import numpy as np
from scipy.special import gammaln


N = 4  #--number of urns
urns = np.ones(N , dtype=np.int)*3
E = sum(urns)
p = {}
Q = {}          #-- N matrix
Nstep = 1000000

for i in range(Nstep):
    n0 = np.random.randint(0 , len(urns))
    if urns[n0]==0:continue
    n1 = np.random.randint(0 , len(urns))
    if n0==n1:continue
    bef = (urns[n0],urns[n1])
    urns[n0]-=1
    urns[n1]+=1
    aft = (urns[n0],urns[n1])
    #-- 本当は十分定常になってからcountすべき
    Q[(bef,aft)] = Q.get((bef,aft),0) + 1
    macro_state = tuple([sum(urns==b) for b in range(E+1)])
    p[macro_state] = p.get(macro_state,0) + 1



vtot = sum(p.values())
ctot = np.exp(np.sum(np.log(np.arange(N,E+N))) - np.sum(np.log(np.arange(1,E+1))))
w = np.zeros(E+1)
for k,v in p.items():
   t = np.array(k)
   c0 = np.exp(np.sum(np.log(np.arange(1,N+1))) - np.sum(np.where(t>0 , gammaln(t+1) , 0)))
   #--macrostate , observed probablitity , expected probability by microcanonical ensemble
   print( k , v/vtot , c0/ctot)
   w += t*(v/vtot)


for (s0,s1),c in Q.items():
    print (s0,s1,Q.get(((s1[1],s1[0]),(s0[1],s0[0])),0),c)


#AlderとWainwrightが1957年に始めた分子動力学計算(hard-sphere MD)では、衝突が起きるまで直進運動し、衝突時に保存則を用いて衝突後の速度を計算するというものであったらしい。MaxwellやBoltzmannが仮定していた分子運動論のイメージに近いのだと思う。この系は、決定論的でエネルギーも連続的であるので、Boltzmannの離散的なモデルよりは、現実寄りのモデルではあると思うけど、衝突によるエネルギー交換という単純なルールは保たれている



『Boltzmannの壺』モデルで、ある特定の壺に入っているボールの数の時間変化を追うと、時間平均はE/Nで、沢山サンプルを集めると、分布は、やはりC exp(-x/T)で近似できるようになるっぽい(数値実験の結果)。勿論、ボールの数はEが最大なので、NやEが有限の時は、C exp(-x/T)という分布は近似にしかならない。ところで、最初の分布は、集団に関するもので、後者は時間に関するものなので、この2つが同分布で近似できるというのは、"エルゴード仮説"っぽい(※)。これもやはり、E/Nを一定にして、NとEを無限大に飛ばした極限を考えると、この2つの分布のずれは0に収束していくと予想される。これも自明ではないと思うけど、厳密な証明は考えてないし、誰かが考えているかどうかも知らない


※)普通のエルゴード仮説は、任意の物理量の相空間上の空間平均と、相空間上の時間発展を追って時間平均を取ると等しいというものだけど、分子動力学屋は、よく、一粒子の物理量の時間方向の分布と、集団に対する分布が一致するという意味で、エルゴード仮説という用語を使っている気がする。MD屋のエルゴード仮説は、通常のエルゴード仮説よりも、強い主張になってる気がするけど、数値計算では成立することが多々あるらしい。Boltzmannのモデルでも、MD屋のエルゴード仮説が(近似的に)成立しているっぽい


#特定の壺のボール数の定常分布は、理屈上は、ボール数がk->k+1,k->k-1に変化する確率を求めれば厳密に計算できる。が、この計算は意外と面倒そうである(注目している壺にk個のボールが入っている時、残りの壺のいくつかは空であるかもしれないというのが面倒の原因)。別の直感的な考え方としては、Nが十分大きければ、2つの壺のボール数にほぼ相関はないと考えられる(例えば、N=2の時は、2つの壺のボール数には強い逆相関がある。Nが増えるにつれて、この相関が弱くなるのは、それほど非直感的ではないと思う)。一方、壺は、どれも特別でないので、各壺のボール数の定常分布は、どれも同じ分布に従うと考えられる。なので、集団としての分布は、独立同分布からサンプリングしたものとほぼ同じになるはず。定常分布に於いて、各壺にボールが(x_1,x_2,...,x_c)個入ってる確率をp(x_1,x_2,...,x_c)とすると、この確率分布は、exchangeable(x_1,...x_cの任意の置換によって値が不変)だから、De Finettiの定理によって、Nを無限大に飛ばした時、独立同分布に従うかのように振る舞い、この直感は正当化できそうな気がする(うさんくさい)


#等重率仮定を使った、よくあるタイプの計算「ある壺に入っているボールの数がxである確率は、残りの壺にE-x個のボールの数を分配する仕方の数W(E-x)に比例する」。例として、N=3,E=6の時を考えると、W(n)=n+1なので、p(x)=(7-x)/28が等重率原理から導かれる結果(平均はE/N=2になっている)。一方、定常分布を、厳密な遷移確率から数値的に計算すると、[p(0),p(1),p(2),p(3),p(4),p(5),p(6)] = [ 0.19047619, 0.25396825, 0.20634921, 0.15873016, 0.11111111, 0.06349206, 0.01587302]となって(検算として平均が2.0であることを確認)、当たり前だけど等重率仮定は全然成立していない。まあ、この例は、NやEが小さすぎて、そもそもBoltzmann-Gibbs分布で全然近似されてないので、問題とは言えないけど(ちなみに数値的に実験した限りでは、あるNまでは、p(0)<p(1)になっていて、Nが大きくなると、p(0)>p(1)と逆転するようである)


論文[1]には、エネルギー交換ルールの時間反転対称性を破ると、Boltzmann-Gibbsでない分布が出ると書いてある。私が最初同じモデルと誤解した、類似のモデルとして、multi-urn Ehrenfest modelがある。こっちは、毎回、ボールをランダムに一つ選んで取り出し、ランダムに選んだ壺に入れるというモデル。Boltzmannのモデルと違って、ボールをランダムに選ぶので、ボールが沢山入っている壺からはボールが取り出される確率が高くなり、少ししかボールのない壺からは、ボールが取り出されにくくなる。結果として、Boltzmannのモデルと違って、極端にボールが増えることもなければ減ることもない。十分時間が経った後は、特定の壺に入っているボールの数の時間分布は、B/Nの周りで正規分布(?)するっぽい。状態の集合は、Boltzmannのモデルと同じで、状態遷移の確率がBoltzmannのモデルと異なる。


multi-urn Ehrenfest modelは、元々は、N=2の時に、粒子の拡散の単純化したモデルとして提案されたらしい(N=2の時は単にEhrenfest modelという)。適当な容器内に沢山粒子があって、容器を2つの領域に分割し、各粒子がどっちの領域にいるかという問題を考えると、Ehrenfest modelとの類似性はイメージできる。


『Boltzmannの壺』とmulti-urn Ehrenfest modelの2つのモデルの違いはコードで書くと割と明確。

import numpy as np
import matplotlib.pyplot as plt

model_type = 0
N = 60  #--number of urns
urns = np.ones(N , dtype=np.int)*20
dist = []

for i in range(50000000):
    if model_type==0:
       n0 = np.random.randint(0 , len(urns))
    else:   #--multi-urn Ehrenfest model
       n0 = np.random.choice(np.arange(0,len(urns)) , p=urns/sum(urns))
    if urns[n0]==0:continue
    n1 = np.random.randint(0 , len(urns))
    if n0==n1:continue
    urns[n0]-=1
    urns[n1]+=1
    dist.append( urns[N//2] )  #--適当な壺の時間発展を追う


cnts = np.bincount(dist)
print( np.sum(np.arange(0,len(cnts))*cnts)/np.sum(cnts) )   #--print average

plt.plot(np.arange(0 , len(cnts)) , cnts)
plt.show()

Ehrenfest modelは、カットオフ現象というものが起きることが証明されてる。カットオフ現象は、多分1980年代に、Diaconisという数学者とその周辺の人々が発見した話らしい。multi-urn Ehrenfest modelにしても、最初のモデルにしても、有限Markov連鎖として定式化できる。有限Markov連鎖では、定常分布に到達する(定常分布に収束するとしてだが)のにかかる時間(混合時間)というものを考えることができるけど、定常分布に一様に近づいていくのかというと、必ずしもそうではなく、ある時間で急に定常分布に近付くということが起きるというのが、Diaconisらの発見らしく、カットオフ現象と呼ばれている(分布間の近さは全変動距離というものを使って測る)。


尤も、日常的な感覚からすると、ボールは区別することなく、また直接目にするのは、それぞれの壺に入っているボール数なので、それは、予想される通り、概ね単調に平均に収束していく。分布そのものの近さを測るということに意味があるのかどうかは、よく分からない


カットオフ現象を示すことが分かっているモデルはいくつかあって、最初に示されたのは多分、リフルシャッフルの数学的モデルと思う

Shuflling cards and stopping times

http://statweb.stanford.edu/~cgates/PERSI/year.html#86

マルコフ連鎖と混合時間 ー カード・シャッフルの数理 ー

http://www.kurims.kyoto-u.ac.jp/~kenkyubu/kokai-koza/H23-kumagai.pdf

2つ目の文献には、「一般のマルコフ連鎖について、カット・オフ現象が起こるための満足のいく必要十分条件は知られていないと言ってよい」と書いてある(数年前の話だけど、多分状況は変わってなさそう?)

Ehrenfest modelでカットオフ現象を示したのも、Diaconisらしい。ちなみに、どっちの解析でも、有限群の表現論が使われている。multi-urnだと、どうなるのか知らない

Asymptotic analysis of a random walk on a hypercube with many dimensions

http://onlinelibrary.wiley.com/doi/10.1002/rsa.3240010105/abstract

Ehrenfestモデルでは、"エントロピー"(正確には、定常分布との相対エントロピー)が単調増加することが厳密に証明されている。これ自体は、もっと広いクラスの有限Markov連鎖で成立するので、特筆すべき性質でもないけど、これはこれで全変動距離とは違う尺度で、定常分布への近さを測っていると見ることもできる。尺度を変えると、カットオフ現象とかがどう見えるのか(あるいは見えないのか)は知らない。ランダムに生成した連結なグラフ上のsimple random walkで実験した限り、KL divergence(負の相対エントロピー)の対数は、かなりキレイに直線的に減少していく(つまり、KL divergenceそのものは指数関数的に減衰する)(数値実験の詳細は下記コード参照)。こういう振る舞いをする理由が十分説明されているのかは不明



import networkx as nx
import numpy as np
import matplotlib.pyplot as plt


G=nx.fast_gnp_random_graph(950,0.02)
P=sorted(nx.connected_components(G), key = len, reverse=True)[0]

#--最大連結成分だけ残す
for c  in G.nodes():
   if not c in P:
       G.remove_node(c)


#-- transition matrix for simple random walk on G
W=np.zeros( (len(G.nodes()) , len(G.nodes())) )
for n , src in enumerate(G.nodes()):
   adjs = G.adj[src].keys()
   cnt = len(adjs)
   for tgt in adjs:
     m = G.nodes().index(tgt)
     W[n,m] = 1.0/cnt


#-- 定常分布の計算
eval,evec = np.linalg.eig(W.T)
idx = np.argmin(np.abs(eval-np.ones(len(G.nodes()))))  #--1.0に最も近い固有値のindex
pi = evec.T[idx]
pi /= np.sum(pi)

#--初期状態(任意)
p = np.zeros( len(G.nodes()) )
p[0] = np.random.random()
p[-1] += 1.0 - p[0]

#--適当なステップ数追跡
maxstep = 50
entropies=[]
for istep in range(maxstep):
    entropy = np.sum(p*np.where(p>0, np.log(p/pi), 0)) #-- KL divergence
    entropies.append(entropy)
    p = np.dot(W.T , p)



#--片対数グラフ
plt.plot(np.arange(0 , len(entropies)) , np.log(entropies) )
plt.show()