Aritalab:Lecture/NetworkBiology/Barabasi-Albert Model

スケールフリー性

多くの研究論文で「スケールフリー」という言葉が使われます。もともとは「桁 (scale) に依存しない (free) 性質」を意味し、昔から物理学の世界で知られていたものです。日本語ではべき則と呼ばれます。ただし、ネットワークについては1999年から言われはじめました。ネットワークの頂点に接続する辺の数を次数といい、全頂点のなかで次数 k の頂点の占める割合を p(k) と書きます。ネットワークがスケールフリーであるとは次数の分布がべき則に従うことを意味し、式で書くと p(k) ∝ k^-λ になります。

この定義は1999年の論文 Barabási A-L, Albert R "Emergence of Scaling in Random Networks" Science 286(5439):509-512 をきっかけに広まりました。

自然界にはべき則が普通にみられ、以下の例がよく知られています。（この表におけるL, Cの意味はリンク解析の項を参照。）

Network (size)	λ	L	C
インターネット (1.5 x 10⁵)	λ_in=2.1, λ_out=2.7	3.1	0.11	Barabasi and Albert (1999)
部分的インターネット (3000-6000)	2.16	3.7	0.18-0.3	Faloutous et al. (1999)
映画俳優の共演関係 (2.2 x 10⁵)	2.3	3.7	0.79	Watts and Strogatz (1998)
PubMedデータベースの共著関係 (1.5 x 10⁶)	?	4.6	0.066	Newman (2000, 2001)

しかし、べき則やそれを生み出すメカニズムは最近わかった事実ではありません。科学界ではおよそ50年周期で話題になってきました。べき則を示す有名な例として、ジップの法則、パレートの法則などがあります。

20世紀初頭
- Pareto, V. (1896, 1897) Cours d’economie politique. Reprinted as a volume of Oeuvres Complètes (Droz, Geneva, 1896, 1965). Pareto, V. Cours d’Economique Politique (Macmillan, Paris,　1897) 所得の分布がべき則に従うことを示す
- Yule, G. (1924) A mathematical theory of evolution, based on the conclusion of Dr. J.C. Willis. F.R.S. Phil. Trans. R. Soc. Lond Ser. B 213, 21–87 優先的選択のメカニズムを示す
20世紀半ば
- Zipf, G.K. (1949) Human Behavior and the Principle of Least Effort, Addison-Wesley 単語の頻度分布がべき則に従うことを示す
- Simon, H.A. (1955) On a class of skew distribution functions. Biometrika 42, 425–440 乗算過程が漸近的にべき則に従うことを示す
21世紀初頭
- Barabási A-L, Albert R (1999) これまでの研究をネットワークの言葉に置き換える

またべき則に近い形である対数正規分布とべき則の関係についても多くの研究がなされました。

乗算過程による対数正規分布の一般性（20世紀初頭）
- Kapteyn, J.C. (1903) Skew frequency curves in biology and statistics in Astronomical Laboratory, Noordhoff, Groningen

様々な事象が対数正規分布をなすことを示した論文

- Gibrat, R. (1931) Les Inegalites Economiques, Libraire du Recueil Sirey, Paris

企業サイズが対数正規分布をなすことを示した論文。その後、都市のサイズ等が対数正規分布に従うことが Gibrat's law と呼ばれる。ただ、都市サイズが power law (zipf's law) に従うとする論文もある。(英語版ウィキペディア Gibrat's law）

乗算過程がべき則を生むメカニズム（20世紀半ば）
- Kesten, H. (1973) Random difference equations and renewal theory for products of random matrices. Acta Mathematica CXXXI, 207–248

対数正規分布をつくるランダムウォークでサイズが小さい部分に反射壁があるとべき分布になることを数理的に示した論文

最近 (21世紀初頭）
- Reed, W.J. and Hughes, B.D. (2002) From gene families and genera to incomes and internet file sizes: why power-laws are so common in nature. Phys. Rev. E 66, 067103

バラバシの論文やレビューにはネットワークのスケールフリー性を自分たちがあたかも初めて見つけたかのように書いていますが、その姿勢が不誠実であることは筆者を含む複数の論文が指摘しています。

Arita, M. (2005 Jul.) Scale-freeness and biological networks. J Biochem. 138(1):1-4.
Keller, EF. (2005 Oct.) Revisiting "scale-free" networks. Bioessays 27(10):1060-8.

Barabási-Albert Model

Barabási, Albertがネットワークの成長モデルとして提唱したのは以下の過程です。(正確には初期条件が異なるが本質は同じ。）

時間ステップ 0 における 1 頂点からスタートする。
単位時間毎に頂点を1つずつ追加し、既に存在する頂点と $m$ 本の辺でつなぐ。
新しい辺はそれぞれ確率 $\textstyle p_i = k_i / \sum_j k_j$ (ここで $k_i$ は頂点iの次数) で接続先を決定する。

この過程は優先的選択(preferential attachment)と呼ばれ、一般にrich gets richerメカニズムなどともいわれます。

べき則のパラメータ $\lambda$ (ラムダ) について

BAモデルから $\lambda = 3$ を簡単に導出できます。初期条件として時刻 0 に 1 頂点の場合、 $t$ 時間後には頂点数 $t+1$ 、辺数 $mt$ です。頂点 i の次数 $k_i$ は時間が経つと単調増加します。ここで、次数を $k_i(t)$ という連続関数として考えてみます。頂点 i の次数は単位時間あたり次数 $k_i$ に比例して増加するので

$\frac{\partial k_i}{\partial t} = m \frac{k_i}{\sum_j k_j} = \frac{k_i}{2t}$

これを解くと $\textstyle k_i(t) = m(t/t_i)^{1/2}$ 、ただし $t_i$ は頂点 $i$ が追加された時間で $k_i(t_i)=m$ となる値です（初期条件）。これを式変形すれば、頂点の次数がある値 $k$ になる時間は $t=t_i(k/m)^2\,$ であることがわかります。頂点 $i$ の次数（期待値）が $k$ のとき、 $t_i$ 以前に追加された頂点の次数は全て $k$ 以上になっています。

頂点の次数分布 p(k) を求めるために、まず頂点の次数が $k$ より大きくなる割合（累積分布関数）を求めます。次数が k より大きくなるような頂点全体の割合は、時刻 $t(m/k)^2\,$ より前に追加された頂点全体の割合に等しくなります。

$\mbox{Pr}[k_i(t) > k] = \mbox{Pr}[t_i < t(\frac{m}{k})^2] = (\frac{m}{k})^2$

左の等号は、頂点 $i$ が追加された時間 $t_i$ が、時刻 $t$ に対して $(m/k)^2\,$ よりも早ければ（小さければ）、次数が $k$ を超えることを示しています。頂点が常に一定量ずつ追加されることを考えると、全体サイズを1とするなら、その割合は単純に $(m/k)^2\,$ としてよいことがわかります。

次数分布 $\mbox{Pr}[k_i(t) = k]\,$ は $\mbox{Pr}[k_i(t) > k] - \mbox{Pr}[k_i(t) > k+1]\,$ と書けるので、その差はkを変数と考えて $\textstyle - \frac{d}{d k}\mbox{Pr}[k_i(t) > k]$ になります。実際に $k$ で微分すると $\mbox{Pr}[k_i(t) = k] = \frac{2m^2}{k^3}$ 。すなわち辺の次数は $k^{-3}$ に比例します。

優先的選択でない場合

新しい辺が結びつく先が次数に比例する値ではなく、一様分布に従うと仮定してみます。頂点 i は頂点数に正比例して辺を得る確率が減っていきます。

$\frac{\partial k_i}{\partial t} = m/(n-1) = m/t$

これを解くと $k_i(t) = m \log t + const.$ です。初期条件として頂点 i　が追加された時間 $t_i$ に辺の数が m であること $k_i(t_i) = m$ をいれると $k_i(t) = m\Big(\log(\frac{t}{t_i})+1\Big)$ 。頂点 i の次数が k になる時間は $k_i \exp(k/m -1)$ であることがわかります。こちらも累積分布関数を計算してみましょう。

$\mbox{Pr}[k_i(t) > k] = \mbox{Pr}\Big[t_i < t\Big(\frac{1}{\exp(k/m - 1)}\Big)\Big] = ee^{-k/m}$

これを上のセクションと同様に微分することで $\mbox{Pr}[k_i(t)=k] = \frac{e}{m}e^{-k/m}$ 、すなわち辺の次数は指数的に減少することになります。

何がべき分布を作るのか

バラバシ－アルバートモデルについて言われる批判のひとつに、自然界においてはべき則のパラメータλの多くが2から3の間をとるのにλ=3のときしか説明できないというものがあります。これはどちらかというと的を外した意見です。

べき則のパラメータが定数になる理由のエッセンスは、優先的選択において次数の時間変化を規定する微分方程式が $\textstyle \frac{dy}{dx} = \frac{y}{2x}$ つまり $y/x$ の形をとることです。この方程式を解くと $y=cx^{1/2}$ という答えが得られ（ $c$ は適当な定数）、この係数1/2が $\lambda=-3$ を作り出します。だから異なる $\lambda$ の値を作り出すには $y/2x$ における比の2を他にずらせばよいのです。