数理生物学に入門してみる（１）力学系

　東京大学教養学部後期課程には「数理生物学」という授業があります。自分はその授業を受講しているのですが、内容についてまとめることで学びになるだろうと思い、ブログの記事にしていくことにしました。なお、この授業の教科書として『細胞の理論生物学ダイナミクスの視点から』（東京大学出版会, 2020）が出版されたので、その内容も（怒られない程度に？）踏まえていくことにします。

細胞の理論生物学: ダイナミクスの視点から

作者:金子邦彦,澤井哲,高木拓明,古澤力
発売日: 2020/04/06
メディア: 単行本

　初回の今回は、教科書の第1章、授業では1, 2回目で扱った力学系についてまとめます。主に常微分方程式の解析に関係した話になっているようです。

力学系
- 力学系とは
- 例
力学系の安定性
力学系を描く
- matplotlib
補足（末尾）
- 断熱消去
- 線形微分方程式の解法について

　生体内で起こる時間変化をシミュレーションするためには、それらを時間変化する複数の変数ととらえ、その動きを考えることが必要になります。即ち、状態に対応した変数 $x_1(t), x_2(t), \cdots , x_n(t)$ の形にモデル化し、それらの時間変化を解析することになります。このとき、各変数の値が時間に応じてどのように変化していくのか（時間発展といいます）を知ることが出来れば、先の変数の状態も計算で予測することが出来ます。つまり、変数の初期値と時間発展の規則が求まれば、その後の全ての時間の状態を求めることが出来ます。

　適度に選んできた変数＝状態の組（状態空間と呼びます）に対し、初期条件と時間発展の規則を与えると、数値的に解析が可能になるということです。例えば、変数が3つ存在したとして、それぞれの数値を表すための軸を用意すれば、3次元の空間で状態を考えることが出来ます。つまり、

$(x_1, x_2, x_3) = (1, 1, 0)$

のように状態を3次元空間上の1点と同じものとみなせます。こう考えたとき、時間発展によってこの点は様々な方向へと動き、軌跡を描くことが想像できます。こうした「状態空間における点の動きを記述するシステム」を力学系と呼びます。

　一方、時間発展の規則をどのように記述するのかも重要です。たとえば、ある位置に存在する点が次の瞬間にどこへ向かうのかが規則にあたります。一番簡単に考えられるのが、原点なら (1, 0, 0) の方向へ向かう、などと「その位置に来たら次の瞬間に向かうべき方向が決まっている」規則です。すなわち、各点での速度ベクトルが固定されており、変数の値の変数で記述できるということです。数式で表すと

$\displaystyle \frac{d \boldsymbol{x}}{d t} = \boldsymbol{f}(\boldsymbol{x})$

と書けます。ここで、 $\boldsymbol{x} = (x_1, x_2, \cdots , x_n)$ は列ベクトルとして扱い、 $\boldsymbol{f}$ は点 $\boldsymbol{x}$ の速度は依存しない関数（のベクトル）です。この形の微分方程式で書ける力学系は、自励系とも呼ばれています。各点で速度のベクトルの規則が決められているため、ベクトルが各点に住んでいるように捉えることが出来ます。この状態はベクトル場と呼ばれており、住み着いているベクトルの規則が $\boldsymbol{f}(\boldsymbol{x})$ として表されています。

f:id:cakkby6:20201010152621p:plain — 自励系の例：各点に灰色の速度ベクトルが住み着いており、赤点で表された初期状態から各点での速度ベクトルに沿って緑の軌跡を描く

例

　自励系の例として、”トグル・スイッチ”、或いは Bistability と呼ばれる系を考えてみます。
　生体の中には、2つの遺伝子がお互いを抑制し合う関係にあるものがあります。つまり、遺伝子1, 2において、各々の発現するタンパク質 $p_1, p_2$ が互いの発現を抑制し合う関係にあるということです。

　タンパク質の発現は、各遺伝子のDNAから mRNA が転写され、それが翻訳されることで行われます。つまり、この系の中には mRNA が含まれており、本来はそれも状態空間の中に組み込むのがよいと言えます。しかし、ｍRNAは即座に合成・分解される「速い変数」にあたります。このような速い変数を無視して考えても誤差はあまりなく、一方で解析は楽に行えるようになります。こうして速い変数を消去することは、断熱消去と呼ばれます（記事の末尾で断熱消去の実際の様子を説明します）。

　mRNAの断熱消去により考えるのはタンパク質の発現量 $p_1, p_2$ のみで良くなりました。次に、これらの関係を化学反応のシミュレーションと同様にして考えていきます。

　 $p_1$ は、何も抑制されていないときは、単位時間あたり $\alpha_1$ の量だけ発現され、 $\beta_1$ の比で自然に分解されていくものとします。これらのパラメータはどんどんと多くなってしまうため、分解される量を 1 とするように $\displaystyle \alpha = \frac{\alpha_1}{\beta_1}$ と正規化します。また、 $p_2$ も同様だとします。この時点で

$\displaystyle \frac{d p_1}{dt} = \alpha - p_1$
$\displaystyle \frac{d p_2}{dt} = \alpha - p_2$

と表されることが分かります。分解される量は比×発現量となっています。

　そしてこのモデルにお互いを抑制する状況を追加します。様々な追加方法があるのですが、ここでは「タンパク質が2量体を形成して、相手の遺伝子の発現を抑制する」状況にします。このとき、タンパク質 $p_1$ の発現が抑制されるためには、 $p_2$ が 2量体を形成する必要があります。 2量体の形成は化学反応として考えられるため、質量作用の法則（化学平衡の法則）により $p_2^2$ に比例した分だけ起こります。そしてこの分だけ $p_1$ の発現を抑制します。

　発現の抑制を式に組み込むのには、極限的な状況を考えてみると良さそうです。つまり、 $p_2^2$ が 0 かと無限大の状況でどのようになるのかを考えます。 $p_2^2 = 0$ の時は、抑制がないため元のモデルの通りになり、 $p_2^2 = \inf$ のときは、全く発現できない状態になると考えます。それを満たすように式にすると、次のようになります。

$\displaystyle \frac{d p_1}{dt} = \frac{\alpha}{1 + p_2^2} - p_1$
$\displaystyle \frac{d p_2}{dt} = \frac{\alpha}{1 + p_1^2} - p_2$

抑制する比率や2量体の形成確率は、ここでは1にしました。こうしてモデルを表す自励系の数式が完成しました。

これを実際に（ベクトル場と共に）グラフにすると以下のようになります。

f:id:cakkby6:20201010160200p:plain — トグル・スイッチのモデルの図示： $\alpha = 4$ としている。3つの初期状態からの軌跡を表示している。

トグル・スイッチのモデルの図示： $\alpha = 4$ としている。3つの初期状態からの軌跡を表示している。

力学系の安定性

安定性

　上の図にて水色とピンクで表された線は、ヌルクラインと呼ばれているものです。ヌルクラインとは、住み着いている速度ベクトルの各要素のうち、いずれかが 0 となるような点の集合で、

$\displaystyle \frac{d x_i}{dt} = f_i(\boldsymbol{x}) = 0$

で表される $\boldsymbol{x}$ の集合のことです。上の図では、

$\displaystyle \frac{d p_1}{dt} = \frac{\alpha}{1 + p_2^2} - p_1 = 0 \Rightarrow p_1 = \frac{\alpha}{1 + p_2^2}$

$\displaystyle \frac{d p_2}{dt} = \frac{\alpha}{1 + p_1^2} - p_1 = 0 \Rightarrow p_2 = \frac{\alpha}{1 + p_1^2}$

の2つの方程式からなります（水色が前者、ピンクが後者）。ヌルクラインの上では、速度ベクトルの要素が0の方向には動かないようになることがわかります。

　ヌルクラインのうち、さらにすべての要素 $f_i$ において0となる点は、固定点（平衡点・不動点など）と呼ばれています。固定点は速度ベクトル自体が0なので、その点から動くことはなくなります。"トグル・スイッチ"の場合は、3点あることが確認できます。

　上の例だと、全ての軌道がいずれかの固定点に向けて動いていく様子が分かります。このように、固定点から少し離れた場所を初期状態としたときに、固定点の方へ向かうのか、それとも固定点から離れていく動きをするのか、といったことは固定点の安定性と呼ばれています。

　安定性は、例えば地図上の山と窪みの地形を考えてみると分かりやすいです。頂点も窪みも両方その場ではボールがとどまることが出来るため、固定点と言えます。窪みから少し離れた場所にボールを置いてみると、窪みの方へ落ちていき、窪みの中に居座り続けるようになります。一方で、山の頂点から少し離れたところにボールを置くと、頂点からは離れていく方向に動いていきます。

　窪みのように安定な点は沈点、頂点のように離れていく不安定な点は湧点（上から見ると湧きだしているように見える）と呼ばれています。他にも、固定点の周囲を回転しながら近づいていく・離れていく点や、回転しながら一定の距離を保つ点、ある方向からは近づけるが、他の方向では離れていく点（鞍点）などがあります。

　このような固定点の安定性を数学的に考えるためには、一種の近似を用います。

　地形の例で考えてみます。各点で速度ベクトルは複雑に変化している場合もありますが、固定点の安定性のみを考える時には、固定点のごく周辺の地形の凹凸だけを考えればよいことが分かります。丸みを帯びた山だったり、かなり急峻な谷かもしれないですが、そのような地形の細かい違いを除いた凹凸のみを見ればよく、それは「地形を平面で近似」したものでも可能です。

　このように、地形を平面で近似するというのは、考えている点での接平面と勾配ベクトルを考えることになります。線形近似と呼ばれる手法です。リミットサイクルの例でいえば、固定点が $p_1, p_2$ の2変数の時間微分に相当する $f_1, f_2$ の関数において、極大か極小かが分かれば安定性が分かります。より具体的には、固定点を $\boldsymbol{x'}$ とおいて

$\displaystyle J = \begin{pmatrix} \frac{\partial f_1(\boldsymbol{x'})}{\partial p_1} & \frac{\partial f_2(\boldsymbol{x'})}{\partial p_1} \\ \frac{\partial f_1(\boldsymbol{x'})}{\partial p_2} & \frac{\partial f_2(\boldsymbol{x'})}{\partial p_2} \end{pmatrix}$

と書かれる、ヤコビ行列を用います。ヤコビ行列の中身は、速度ベクトルを表す関数の各変数による偏微分であり、接平面に対応しています。

　線形近似を行ってから、どのように安定性を考えるのかは次の節で扱います。

安定性と固有値

　今考えている微分方程式は線形近似しているため、

$\displaystyle \begin{pmatrix} \frac{d p_1}{d t} \\ \frac{d p_2}{d t} \end{pmatrix} = \begin{pmatrix} \frac{\partial f_1}{\partial p_1} & \frac{\partial f_2}{\partial p_1} \\ \frac{\partial f_1}{\partial p_2} & \frac{\partial f_2}{\partial p_2} \end{pmatrix} \begin{pmatrix} p_1 \\ p_2 \end{pmatrix}$

と書けます。これは行列Aによって書ける

$\displaystyle \frac{d \boldsymbol{x}}{dt} = A \boldsymbol{x}$

という形の線形微分方程式です。この方程式に対して、固定点（ヤコビ行列にするために原点へ平行移動している）の安定性を解析するためには、実際に解いてしまえばよいです。線形近似すると嬉しい理由というのは、この線形微分方程式の解法が良く知られており、解析が可能だからです。

　解法全体を紹介するのは末尾で行うとして、実際に重要なのは固有値の正負と虚部だと言えます。すなわち、行列Aを対角化して得られる各固有値 $\lambda_i$ の正負・虚部の存在によって固定点の安定性が場合分けできます。ここでは2変数の場合のまとめだけ出すことにします。

両方実数で負：沈点（吸い込まれる）安定

両方実数で片方が正、片方が負：鞍点（片方は吸い込まれ、片方は離れていく） 不安定

両方実数で正：湧点（離れていく）不安定

両方複素数で実部が負：焦点（回転しながら吸い込まれる） （漸近）安定

両方複素数で実部が正：不安定焦点（開店しながら離れていく）不安定

両方純虚数：渦心点（一定の距離で点の周りを回転する）安定

　例えば"トグル・スイッチ"の場合は、 $p_1 = p_2$ の位置の固有点が不安定（鞍点）、残りが安定（沈点）だと言えます（ヤコビ行列を考えると分かります）。

　なお、固有値は大雑把には固有点からの距離の増減に対応しており、負ならば距離が減る=吸い込まれる方向、正なら距離が増える=離れていく方向だとわかります。また、固有値が複素数になる場合は、（波の複素表示などと同様に） $\cos, \sin$ の形に変形されるため、回転の運動が絡むようになると直感的に理解できます。

アトラクター

　固定点のほかにも、周期的な軌道となる場合も存在し、（安定）リミットサイクルと呼ばれています。こうして、いかなる初期条件にあっても、十分に時間がたった後には固定点や安定な軌道だけに存在することが分かると思います。この十分時間がたった後に通る軌道は、アトラクターと呼ばれています。

　Poincare-Bendixson の定理によると、軌道が有限の範囲にあるならば、2変数の自励系でのアトラクターは固有点かリミットサイクルしかないと言えます。よって、力学系においては、十分な時間がたって普段考えるべき状態（アトラクター）として主に固定点とリミットサイクルなどを考えればよいと言えます。

　また、各初期条件から複数のアトラクターへと移っていくとき、どのアトラクターになるのかという範囲が存在することが想像できます。この範囲をベイスンと呼びます。

分岐

　リミットサイクルの時の $\alpha$ のように、微分方程式の中にパラメータが含まれていることが多いです。このパラメータを動かすと、固定点やアトラクターが変化していきます。パラメータの値によってはアトラクターの安定性が変化することもあります。こうした定性的な変化を分岐(bifurcation)と呼びます。

(1) transcritical bifurcation
　パラメータによって安定性が切り替わる場合です。例えば

$\displaystyle \frac{d u}{d t} = a \frac{u}{u+1} - u$

と書かれる場合、固定点は $0, a$ の2点であり、 $a > 0$ では $x = 0$ が不安定、 $x = a$ が安定となりますが、 $a < 0$ では $x = 0$ が安定、 $x = a$ が不安定に切り替わります。このような1次元の場合の考察は、 $\dot{u}-u$ 平面を描いてみると分かりやすくなります。

f:id:cakkby6:20201011204119p:plain — 安定性の考察： $\dot{u}$ が正ならば $u$ を増加させる動き、負ならば減少させる動きをするため、点Aでは少しずれたときにAへと戻る動きから安定だと分かり、点BはBから離れる動きをするため不安定だと分かる

安定性の考察： $\dot{u}$ が正ならば $u$ を増加させる動き、負ならば減少させる動きをするため、点Aでは少しずれたときにAへと戻る動きから安定だと分かり、点BはBから離れる動きをするため不安定だと分かる

(2) saddle-node
　不安定な固定点と安定な固定点が存在し、パラメータを変えることによってそれらが衝突して対消滅する場合です。例えば次の自励系で生じます。

$\displaystyle \frac{d u}{d t} = a - u^2$

(3) pitchfork
　パラメータを変えると、安定だった1つの固定点が不安定になり、さらに2つ安定な固定点が生じるものです。"トグル・スイッチ"はこのパターンで $\alpha < 2$ では固定点が1つ、 $\alpha > 2$ では固定点が3つで最初に示した図のようになります。

(4) Hopf
　複素数を用いたモデルなどで、pitchforkのと同様に安定だった固定点がパラメータの変化によってリミットサイクルを生じる場合です。

　このように、多くの種類の分岐が存在します。こうした分岐は一見して異なるシステムでも同様の者に分類することができ、メカニズムを等価なものとみなすことが出来るようになります。

以上で教科書の第1章に相当する範囲が終わりました。第2章ではこれを細胞の実際の系について具体的に考えて行くことになります。

cake-by-the-river.hatenablog.jp

力学系を描く

matplotlib

途中で現れたグラフの図示には、python の matplotlib (とnumpy) と呼ばれるモジュールを用いました。2次元版と3次元版とのコードを載せるので、環境を持っている人は試してみるといいと思います。

　微分方程式を数値的に解くために、オイラー法と呼ばれる方法を用いました。微分を行う代わりに、微小な範囲の差分で代替する方法です。簡単に実装できますが、誤差がある程度生じます。

　 2次元版のコード（"トグル・スイッチ"の実装）

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
import numpy as np

plt.figure()

LX, LY = 5.01, 5.01  # メッシュの大きさ
gridwidth = 0.1
# メッシュの各点を作成
vX, vY = np.meshgrid(np.arange(-LX, LX, gridwidth), np.arange(-LY, LY, gridwidth))

# 軌跡をしまうリスト
X = []
Y = []

# 初期条件
X1_s, Y1_s = 0, 0.5
X2_s, Y2_s = 3, 2.6
X3_s, Y3_s = 4, 4

# 計算に用いる変数
X1, Y1 = X1_s, Y1_s
X2, Y2 = X2_s, Y2_s
X3, Y3 = X3_s, Y3_s

t = 0
# 時間の刻み幅
h = 0.0005

# パラメータと方程式
alpha = 4
f_x = lambda x, y: alpha/(1 + y**2) - x
f_y = lambda x, y: alpha/(1 + x**2) - y

# オイラー法
for i in range(32000):
    X.append(X1)
    Y.append(Y1)
    X.append(X2)
    Y.append(Y2)
    X.append(X3)
    Y.append(Y3)

    X1 += h * f_x(X1, Y1)
    Y1 += h * f_y(X1, Y1)
    X2 += h * f_x(X2, Y2)
    Y2 += h * f_y(X2, Y2)
    X3 += h * f_x(X3, Y3)
    Y3 += h * f_y(X3, Y3)
    t += h

# ベクトル場の作成
U = f_x(vX, vY) / 2
V = f_y(vX, vY) / 2

plt.quiver(vX, vY, U, V, color="#33333333", scale_units='xy', scale=5.0)

# ヌルクラインを描く
nX = np.linspace(0, 6, 2000)
plt.scatter(nX, alpha/(1+nX**2), s=1, color='#ffaaaa22')
nY = np.linspace(0, 6, 2000)
plt.scatter(alpha/(1+nY**2), nY, s=1, color='#00ccff22')

plt.scatter(X, Y, s=1, color="#00aa33")
plt.scatter(X1_s, Y1_s, s=5, color='red')
plt.scatter(X2_s, Y2_s, s=5, color='red')
plt.scatter(X3_s, Y3_s, s=5, color='red')

plt.xlim([-0.3, LX])
plt.ylim([-0.3, LY])

plt.show()

　3次元版のコード（とあるベクトル場で）

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
import numpy as np

fig = plt.figure()
ax = Axes3D(fig)

LX, LY, LZ = 4.01, 4.01, 4.01  # メッシュの大きさ
gridwidth = 0.1
vX, vY, vZ = np.meshgrid(np.arange(-LX, LX, gridwidth), np.arange(-LY, LY, gridwidth),np.arange(-LZ, LZ, gridwidth))

# 軌跡をしまうリスト
X = []
Y = []
Z = []

# 初期条件
X1_s, Y1_s, Z1_s = -1, 1, 0
X1, Y1, Z1 = X1_s, Y1_s, Z1_s

t = 0
# 時間の刻み幅
h = 0.0005

# パラメータと方程式
a = 0.1
b = 0.1
r = 0.1
r_d = 0.1
f_x = lambda x, y, z: a - r * x * y
f_y = lambda x, y, z: b - r * x * y
f_z = lambda x, y, z: r * x * y - r_d * z


# オイラー法
for i in range(32000):
    X.append(X1)
    Y.append(Y1)
    Z.append(Z1)

    X1 += h * f_x(X1, Y1, Z1)
    Y1 += h * f_y(X1, Y1, Z1)
    Z1 += h * f_z(X1, Y1, Z1)
    t += h

# ベクトル場の作成
U = f_x(vX, vY, vZ)
V = f_y(vX, vY, vZ)
W = f_z(vX, vY, vZ)

ax.quiver(vX, vY, vZ, U, V, W, color="#33333333", length=1, normalize=False)
ax.scatter3D(X, Y, Z, color="#00aa33")
ax.scatter3D(X1_s, Y1_s, Z1_s, "o", s=40, color='red')

ax.set_xlim([-LX, LX])
ax.set_ylim([-LY, LY])
ax.set_zlim([-LZ, LZ])

plt.show()

補足（末尾）

断熱消去

　 $p_1, p_2$ に対応するmRNAの量を $m_1, m_2$ とおきます。互いの遺伝子の発現抑制は、実際にはタンパク質がお互いのmRNAの発現を抑制しているため、mRNAの発現速度の項に抑制がかかります。（両方の遺伝子で同じ速度として、）mRNAの合成速度を $S$ 、分解速度を $\gamma$ とおき、タンパク質の合成・分解速度は 1 とすると

$\displaystyle \frac{d m_1}{dt} = \frac{S}{1 + p_2^2} - \gamma m_1$
$\displaystyle \frac{d m_2}{dt} = \frac{S}{1 + p_1^2} - \gamma m_2$

$\displaystyle \frac{d p_1}{dt} = m_1 - p_1$
$\displaystyle \frac{d p_2}{dt} = m_2 - p_2$

とかけます。抑制の項はタンパク質自体で考察した時と同じです。ここで、翻訳よりも転写の方が早く行われる、つまり $S, \gamma \gg 1$ とすると、タンパク質の量の変化に対してmRNAの方は十分素早く変化するため、考えたいタイミングではおよそ一定の値を取った中間の役割をしているだろうと考えられます。

　そこで、 $m_1, m_2$ が一定、つまり $\displaystyle \frac{d m_1}{dt} = \frac{d m_2}{dt} = 0$ が成り立つとして変数を消去します。

$\displaystyle m_1 = \frac{S}{\gamma (1 + p_2^2)}$
$\displaystyle m_2 = \frac{S}{\gamma (1 + p_1^2)}$

この $m_1, m_2$ を先ほどの4式の後半に代入して

$\displaystyle \frac{d p_1}{dt} = \frac{S}{\gamma (1 + p_2^2)} - p_1$
$\displaystyle \frac{d p_2}{dt} = \frac{S}{\gamma (1 + p_1^2)} - p_2$

この式の $\displaystyle \frac{S}{\gamma}$ を $\alpha$ と置きなおせば、上述した式の形になります。

線形微分方程式の解法について

　今回は対角化できる行列のみを扱うことにします。そうでない例に関しては、広義固有空間やジョルダン標準形によって解けますが、なかなか長くなってしまうので今回は扱いません。

$\displaystyle \frac{d \boldsymbol{x}}{dt} = A \boldsymbol{x}$

　（天下り的ですが）行列で書かれたこの微分方程式は、もし $\boldsymbol{x}, A$ がベクトル・行列じゃなかったとしたらすぐに指数関数 $x=e^{tA}$ が解だと分かります。そこで、もはやぶっちゃけて $A$ が行列の場合も、 $x = e^{tA}$ が解だと考えてしまいましょう。

　このように、指数関数の肩に行列が乗ったものが行列の指数関数というものです。実際の定義はテイラー展開を用いて次のように書ける行列です。

$\displaystyle e^{tA} = I + tA + \frac{t^2}{2}A^2 + \frac{t^3}{6}A^3 + \cdots$

べき級数の形で書け、それぞれを微分することが可能なため、実際にこの行列の指数関数が先の微分方程式の（一種の）解になることが分かります。また、初期条件 $\boldsymbol{x'}$ に対し、

$\boldsymbol{x} = e^{(t - t_0)A}\boldsymbol{x'}$

が一般解を成すことが言えます（1次元の例でも直感的に分かる）。問題は、この行列の指数関数をそのまま考えるのが面倒だと言うことです。ところが、行列 $A$ を対角化できるならば、これをすっきりした形で考えることが可能になります。

　行列の指数関数の特徴として、次の式が成り立つことが言えます。

$D = P^{-1} A P$

と対角化できるとき、

$e^{tA} = e^{tPDP^{-1}} = P e^{tD} P^{-1}$

と書ける。

この式は、 $D^n = (PAP^{-1})^n = P A^n P^{-1}$ により簡単に示せます。そしてこの式によって $e^{tD}$ の部分のみを考えればよくなるのですが、 $D$ が対角化された行列であることから、

$\displaystyle D = \begin{pmatrix} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda_n \end{pmatrix} \Rightarrow e^{tD} = \begin{pmatrix} e^{t \lambda_1} & 0 & \cdots & 0 \\ 0 & e^{t \lambda_2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & e^{t \lambda_n} \end{pmatrix}$

と書けます。また、 $P$ などは解空間（解となるベクトル全体がなす空間）を変えるのではなく、その基底となるベクトルの向きを変更させるに過ぎず、それらを考慮して新しく $\boldsymbol{v_i$ ] という基底の下で考えると、微分方程式の一般解を

$\boldsymbol{x} = c_1 e^{\lambda_1 t} \boldsymbol{v_1} + c_2 e^{\lambda_2 t} \boldsymbol{v_2} + \cdots + c_n e^{\lambda_n t} \boldsymbol{v_n}$

と書けます。時間変化に関係する部分は指数関数で表されています。そしてその大きさは各固有値に、その方向（基底）は $\boldsymbol{v_i}$ という固定されたベクトルで表されます。

　このことから、固有値が解の安定性を説明することが分かります。もし固有値が正ならば、指数関数の肩が常に正となり、指数関数の値は1以上＝その方向のベクトルが伸びていくため、原点（固定点）から離れていきます。逆に負ならば指数関数の値が減っていくために近づいていきます。虚数ならば減衰振動のようになることも分かります。

以上が対角化可能な場合の線形微分方程式の解き方です。対角化不可能な場合はジョルダン標準形という形で計算することで似たようなことが可能になります。