LoginSignup
9
5

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@gkhtnk(Gaku Hatanaka)

（神経科学で学ぶ）matlab 初心者向け講座 HHM編

Last updated at Posted at 2020-05-06

はじめに

この講座では matlab初心者向けに、matlabの機能を神経科学に関係する内容に触れながら紹介します。
matlabの基本的な扱いを学んで、もう少し実践的な内容に触れてみたい方を対象にしています。

今回はHodgkin-Huxley方程式に従う神経細胞のモデルを作成することを題材にして、matlabのクラスの使い方を紹介します。

Hodgkin-Huxley方程式 とは

Hodgkin-Huxley方程式はAlan Lloyd HodgkinとAndrew Fielding Huxleyによって考案された神経軸索の電気的な活動を記述する方程式です。細胞膜のイオンチャネルでのイオンの流入出と、細胞膜に保持される電荷の変化を用いて入力電流に対する細胞膜の電位変化を模倣することを試みます。
原著論文はHodgkin and Huxley (1952d) J Physiol 1です。 またその詳細については日本語では脳科学辞典で丁寧に解説されています。Hodgkin-Huxley方程式は膜電位の時間変化を記述する常微分方程式と、膜電位の空間伝播を記述する偏微分方程式の２種類がありますが、本稿では常微分方程式のみを扱います。

クラス とは

クラスとは、データとそのデータに関連する処理・操作を関連づけて定義したものです。
プログラミング言語ごとにその思想は異なりますが、matlabにおけるクラスの考え方は「MATLAB でのオブジェクト指向デザイン」を確認してください。2

チャネル クラスの実装

初めに、イオンの通り道であるチャネルをmatlabのクラスとして実装します。

関数ファイルを作成するときと同様に、クラス名を付けた.mファイルを作成します。
クラスは Classdef ブロック に定義します。
関数ファイルを作成するときと同様に、クラス名とファイル名は一致する必要があります。

Channel.m
classdef Channel
  ...
end

電位依存性チャネルは電圧に対してイオンの流入出量を変化させます。

:point_right: %%%%%%%%【詳細説明を表示する】%%%%%%%% :point_left: ◀

Hodgkin-Huxley方程式では、イオンの流入出によって生じるイオン電流の大きさをオームの法則に従って表記しました。

オームの法則は電流 $I$ と電圧 $V$ 、抵抗 $R$ の関係式を定義します。
Hodgkin-Huxley方程式では抵抗 $R$ の代わりに、コンダクタンス $G$ (抵抗の逆数、電流の流れやすさを表す)を使用します。3

I = GV, \qquad \left(I = \frac{V}{R} : \text{オームの法則}, \quad G = \frac{1}{R} : \text{コンダクタンスと抵抗の関係} \right)

チャネルにかかる電圧は、電位 $V$ とイオンの濃度勾配によって生じる起電力 $V_X$ の差で表記することができます。したがってチャネルを流れるイオン電流 $I_X$ は以下のように表記できます。

$$
I_X = G_X \left( V - V_X \right), \qquad \left( V = \frac{I_{X}}{G_{X}} + V_{X} \right)
$$

電位依存性チャネルの基本的な性質として、以下の2点がまとめられます。

  • チャネルはコンダクタンス $G_X$ と 起電力 $V_X$ という要素を持つこと
  • 電位 $V$ が与えられることで大きさ $I_X$ のイオン電流を流すこと

そこでこれらの性質を Channel クラスに実装していきます。

プロパティの実装

クラスが保持するデータはクラス プロパティとして、Properties ブロックの中に実装します。
ここでは、チャネルの持つ2つのデータ $G_X$ と $V_X$ を用意します。

Channel.m
classdef Channel
  properties
    Gx
    Vx
  end
end

クラスのオブジェクト（以下では chX）はクラス名（Channel）を使用して作成されます。
クラスのプロパティへはドット（.）とプロパティ名（Gx）を用いてアクセスできます。

>> chX = Channel();
>> chX.Gx = 0;
>> chX.Gx
ans =
     0

プロパティブロックでは各要素に対して、既定値を割り当てることもできます。
例えば、以下ではプロパティ PI に16桁の精度で円周率を代入しています。

  properties
    PI =  3.141592653589793
  end

今回は次に述べるコンストラクターを用いてプロパティを初期化するので、
プロパティブロックで既定値を割り当てずにおいておきます。

詳細は、プロパティの定義を参照してください。

メソッドの実装

クラスの機能はMethods ブロックの中に実装します。

コンストラクターの実装

コンストラクターはクラスの実体（インスタンス）が作成される際に実行される特殊な関数です。
Methods ブロックの中にクラス名と同じ名前の関数として実装します。

methods
  function obj = Channel()
    ...
  end
  ...

コンストラクターの最初の出力は作成されるクラスオブジェクト（上の例では obj）になります。4

以下のように、クラスを実体化する際にクラスオブジェクト（以下では chX）を生成します。

>> chX = Channel();

メソッド内でクラスプロパティ（例えば Gx）にアクセスする際はクラスオブジェクト（obj）とドット（.）に続けてプロパティ名を指定します。
今回は、クラスを初期化する際に、2つの値（src_Gx, src_Vx）を受け取り、2つのプロパティ（Gx, Vx）に代入して初期化するようにしておきます5。

Channel.m
classdef Channel
  ...
  methods
    function obj = Channel(src_Gx, src_Vx)
      obj.Gx = src_Gx;
      obj.Vx = src_Vx;
    end
  end
end

ここまでで、Channel クラスを初期化する際に、プロパティ Gx, Vx に規定値を代入することができるようになりました。
例えば、以下のように初期値として 5, 7 を与えることで、2つのプロパティに値が渡されるのがわかります。

>> chX = Channel(5, 7)
chX = 
  Channel のプロパティ:
    Gx: 5
    Vx: 7

イオン電流関数の実装

次に、クラスの持つプロパティの値（Vx, Gx）と、入力引数で与えられる電位（$V$）の値から、チャネルのイオン電流（$I_X$）を計算する処理をメソッドとして実装します。

I_X = G_X \left( V - V_X \right)
Channel.m
classdef Channel
  ...
  methods
    ...    
    function dst_Ix = Ix(obj, src_Vm)
      dst_Ix = obj.Gx.*(src_Vm - obj.Vx);
    end
  end
end

生成されたクラスの実体（ここではchX）にドット(.)とメソッド名(Ix())を用いて実装した関数を呼び出します。

>> chX = Channel(5, 7);
>> chX.Ix(11)
ans =
    20

関数内でプロパティにアクセスしたのと同様に、クラス実体（chX）からもプロパティにアクセスできます。
関数の計算が正しいか、同じ計算をして検証します。

>> chX.Gx
ans =
     5

>> chX.Vx
ans =
     7

>> chX.Gx.*(11 - chX.Vx) % Verification
ans =
    20

ここまでで、チャネルの基本的な部分は用意できました。

:point_right: %%%%%%%%【コードを表示する】%%%%%%%% :point_left: ◀
Channel.m
classdef Channel
  properties
    Gx
    Vx
  end
  methods
    function obj = Channel(src_Gx, src_Vx)
      obj.Gx = src_Gx;
      obj.Vx = src_Vx;
    end
    
    function dst_Ix = Ix(obj, src_Vm)
      dst_Ix = obj.Gx.*(src_Vm - obj.Vx);
    end
  end
end

リークチャネル クラスの実装

一番実装が簡単なものから始めましょう。
Hodgkin-Huxley方程式ではリークチャネルは純粋に抵抗 $G_L$ と起電力 $V_L$ からなる回路を仮定します。

つまり、電圧 $V$ に対して、チャネルを流れる電流は以下のように表記できます。

$$
I_L = G_L \left( V - V_L \right)
$$

このため先ほど実装した Channel.m をそのまま使用してもよいのですが、
本稿では、クラスの練習も兼ねて、クラスの継承を実践してみます。

クラスの継承

クラスを継承することで、継承元のクラス(基底クラス)の性質、つまりプロパティとメソッドを受け継いだ新たなクラス（派生クラス）を作成することができます。

スーパークラス/基底クラス サブクラス/派生クラス
継承されたクラス 継承したクラス
Channel ChLeak

詳しくは、クラスの階層 — 概念を参照してください。

classdef 行において、< を用いてスーパークラスを指定します。
ここではスーパークラスを Channel を継承して、サブクラス chLeak を定義します。

新しいファイル（ChLeak.m）を用意して以下のとおりに記述します。

ChLeak.m
classdef ChLeak < Channel
end

スーパークラスメソッドの呼び出し

サブクラスでは継承したスーパークラスのメソッドを使用することができます。以下のように、ChLeak クラスは Channel クラスと同様にコンストラクターとイオン電流関数（Ix()）を使用することができます。

>> chL = ChLeak(5, 7);
>> chL.Ix(11);
ans =
    20

Hodgkin and Huxley (1952d)では、リークチャネルのコンダクタンスと起電力を $G_L = 0.3$ 、 $V_L = -10.613$ と求めています。これらの定数を実体を生成するごとに毎回代入することは面倒です。

ここでは、継承したコンストラクターを常に2つの値を代入するように変更します。スーパークラス（Channel）のコンストラクターはサブクラスで生成するサブクラスオブジェクト(obj)を利用して以下のように呼び出します。

SubClassObject@SuperClassName(...)

% Example
obj@Channel(...)

スーパークラス（Channel）のコンストラクターは src_Gx 、 src_Vx の2つの引数を取ります。
今回定義するサブクラス（ChLeak）のコンストラクターは引数を取らず、内部で Channel クラスのコンストラクターに2つの値を渡すように書き換えます。

クラスの基本的な性質は基底クラス（Channel）で定義されているので、派生クラス（ChLeak）での実装は追加で必要なところに留まります。

chLeak.m
classdef ChLeak < Channel
  methods
    function obj = ChLeak()
      obj@Channel(0.3, -10.613);
    end
  end
end

詳細は、サブクラス コンストラクターを参照してください。

コンストラクターの変更によって、クラス ChL を初期化する際には引数を渡さなくてもよくなりました。

>> chL = ChLeak()
chL = 
  ChL のプロパティ:
    Gx: 0.3000
    Vx: -10.6130

カリウムチャネル クラスの実装

カリウムチャネルもリークチャネルと同様に Channel クラスを継承して作成していきます。カリウムチャネルのコンダクタンスと起電力は $G_K = 36$ 、 $V_K = 12$ です。6

ChPotassium.m
% Potassium(Kalium) Ion Channel
classdef ChPotassium < Channel
  methods
    function obj = ChPotassium()
      obj = obj@Channel(36, 12);
    end
  end
end

カリウムチャネルのコンダクタンスの大きさは一定ではなく、膜電位の値に従って変化します。
Hodgkin-Huxley方程式では、電位に依存するカリウムチャネルの $n$ ゲートの開口確率 $P_n(V)$ を用いて以下のように表記します。7

$$
G_{K} = G_{K}^{max} {P_{n}(V)}^4
$$

ゲートの2状態モデル

Hodgkin-Huxley方程式のチャネルはゲートの開、閉の2状態によってコンダクタンスが変化します。

:point_right: %%%%%%%%【詳細説明を表示する】%%%%%%%% :point_left: ◀

チャネルのゲート（ここでは一般的にゲート $\xi$ とします）は開状態 $S_{\xi}^{Open}$ と閉状態 $S_{\xi}^{Close}$ からなる2状態モデルとして表現します。

\boxed{S_{\xi}^{Close}}
\begin{array}{}
\xrightarrow{\alpha_{\xi} (V)} \\
\xleftarrow[\beta_{\xi} (V)]{}
\end{array}
\boxed{S_{\xi}^{Open}}
\qquad
\begin{array}{}
P_{\xi}(V) = P \left(S_{\xi} = S_{\xi}^{Open} \middle| V \right)
\end{array}

このとき、ゲート開口確率 $P_{\xi}$ の時間遷移は以下の微分方程式で表すことができます。

\frac{dP_{\xi}(V)}{dt} = \alpha_{\xi} (V) (1 - P_{\xi}) - \beta_{\xi} (V) P_{\xi}

また、チャネルの定常状態 $\left( t \rightarrow \infty; \frac{dP_{\xi}}{dt} = 0 \right)$ を考えると以下のようになります。

P_{\xi} (t=\infty, V) = \frac{\alpha_{\xi} (V)}{\alpha_{\xi} (V) + \beta_{\xi} (V)}

ゲートの2状態モデルの実装

それでは、ゲートの2状態モデルを組み込んだカリウムチャネルクラスの実装をおこなっていきます。

ゲート開口確率の追加

さて、カリウムチャネル クラスでは、ゲート開口確率を利用する必要が出てきました。
サブクラス内でプロパティを定義すると、そのプロパティはサブクラス専用のものになります。

ChPotassium.m
classdef ChPotassium < Channel
  properties
    Pn
  end
  ...
end

状態遷移速度の追加

次に、状態遷移速度 $\alpha_n (V)$ と $ \beta_n (V)$ を実装します。

Hodgkin and Huxley (1952d)では、電位に依存する状態遷移速度 $\alpha_n (V)$ と $ \beta_n (V)$ を実験的に以下のように求めました。

\begin{aligned}
\alpha_n (V) &= \frac{0.01 (V + 10)} {\exp{\left(\frac{V + 10}{10}\right)} - 1} &
\beta_n (V)  &= 0.125 \exp{\left(\frac{V}{80}\right)}
\end{aligned}

状態遷移速度は膜電位 $V$ のみに依存し、ほかのプロパティ($G_K$, $V_K$, $P_n$)の値に依存しません。
つまり、クラスオブジェクトを使うことなく関数を実装できます。

このような関数は静的メソッドとして実装します。
静的メソッドは今回のようにクラスの操作には関係するものの、その実体（インスタンス）の状態（データ）とは関連しない関数です。

matlabではこのような性質を以下のようにmethods 行において属性として指定することができます。属性の指定を受けるのは、属性指定を伴う Methods ブロックのみです。

methods
  % 普通のメソッド
  ...
end

methods (Static)
  % 静的メソッド
  ...
end

属性には他にもアクセスを制限を指定するものなど、いくつか種類があります。
詳細は、メソッドの属性を参照してください。

静的メソッドでは、関数の最初の引数としてクラスオブジェクト（obj）を渡す必要がなくなります。

ChPotassium.m
classdef ChPotassium < Channel
  methods
    ...
  end
  ...
  methods (Static)
    function dst_an = an(src_Vm)
      dst_an = 0.01.*(src_Vm + 10)./(exp((-src_Vm + 10)/10) - 1);
    end
    
    function dst_bn = bn(src_Vm)
      dst_bn = 0.125.*exp(src_Vm./80);
    end
  end
end

また、クラスの実体を伴わずにメソッドの実行を行うことができます。
　

:point_right: %%%%%%%%【コードを表示する】%%%%%%%% :point_left: ◀

試しに論文中の図４のグラフを描画してみます。
$\alpha_n$ も $\beta_n$ も $V$ の関数として定義しており、fplotを用いて描画できます。

>> ChPotassium.an(0)
ans =
    0.0582

% Hodgkin and Huxley (1952d) Fig. 4
>> hold on;
>> fplot(@ChPotassium.an, [-110, 50]);
>> fplot(@ChPotassium.bn, [-110, 50]);

image.png

コンストラクターの変更

追加したゲート開口確率 $P_n$ の初期化処理を実装します。
ある電圧 $V$ におけるゲート開口確率 $P_n$ の定常状態の値は以下のように表記できます。

P_n (t=\infty, V) = \frac{\alpha_n (V)}{\alpha_n (V) + \beta_n (V)}

この値は初期状態における電圧に依存するため、再びコンストラクターに引数を追加します。

ChPotassium.m
classdef ChPotassium < Channel
  ...
  methods
    function obj = ChPotassium(src_Vm)
      obj = obj@Channel(36, -12);
      obj.Pn = obj.an(src_Vm)./(obj.an(src_Vm) + obj.bn(src_Vm));
    end
  end
  ...
end

イオン電流関数の変更

イオン電流関数を変更します。

スーパークラス（Channel）では以下のように定義されています。

Channel.m
    function dst_Ix = Ix(obj, src_Vm)
      dst_Ix = obj.Gx.*(src_Vm - obj.Vx);
    end

カリウムチャネルのコンダクタンスは以下の様に与えられます。

$$
G_{K} = G_{K}^{\text{max}} {P_{n}}^4
$$

したがってイオン電流の大きさはスーパークラス（Channel）のものに $P_n^4$ をかければ良いことがわかります。

サブクラス（ChPottasium）のメソッドでスーパークラス（Channel）の 同名の メソッドを呼び出す際は@を用いて、以下の様に記述します。

SuparClassMehod@SuperClassName(SubClassObject, ...);

% Example; Both OK
Ix@Channel(obj, src_Vm);
obj.Ix@Channel(src_Vm);

ここでは、obj.Ix@Channel(src_Vm) で　$G_x$ と $V_x$ と膜電位から求めたイオン電流を求め、ゲート開口確率 $P_n$ （の４乗）で補正します。

ChPotassium.m
classdef ChPotassium < Channel
  ...
  methods
    ...
    function dst_Ix = Ix(obj, src_Vm)
      dst_Ix = obj.Ix@Channel(src_Vm).*obj.Pn.^4;
    end
  end
  ...
end

詳細は、サブクラス オブジェクトでのスーパークラス メソッドの呼び出しを参照してください。

微分方程式の実装

$n$ ゲートの開口確率 $P_n$ の時間変化は以下の微分方程式で表記できます。

\frac{dP_{n}}{dt} = \alpha_{n} (V) (1 - P_{n}) - \beta_{n} (V) P_{n}

先で数値計算することを考え、ここでは引数に $P_n$ の微小変動 $\delta P_n $ を加えられるようにしておきます。

ChPotassium.m
classdef ChPotassium < Channel
  ...
  methods
    ...
    function dst_dn = dn(obj, src_Vm, del_Pn)
      tmp_Pn = obj.Pn + del_Pn;
      dst_dn = obj.an(src_Vm).*(1 - tmp_Pn) - obj.bn(src_Vm).*tmp_Pn;
    end
  end
  ...
end

カリウムチャネルのソースコード

以上をまとめると、カリウムチャネルの実装が完成します。

:point_right: %%%%%%%%【コードを表示する】%%%%%%%% :point_left: ◀
ChPotassium.m
% Potassium(Kalium) Ion Channel
classdef ChPotassium < Channel
  properties
    Pn
  end
  
  methods
    function obj = ChPotassium(src_Vm)
      obj = obj@Channel(36, 12);
      obj.Pn = obj.an(src_Vm)./(obj.an(src_Vm) + obj.bn(src_Vm));
    end
    j.
    function dst_Ix = Ix(obj, src_Vm)
      dst_Ix = obj.Ix@Channel(src_Vm).*obj.Pn.^4;
    end
       
    function dst_dn = dn(obj, src_Vm, del_Pn)
      tmp_Pn = obj.Pn + del_Pn;
      dst_dn = obj.an(src_Vm).*(1 - tmp_Pn) - obj.bn(src_Vm).*tmp_Pn;
    end
  end
  
  methods (Static)
    function dst_an = an(src_Vm)
      dst_an = 0.01.*(src_Vm + 10)./(exp((src_Vm + 10)/10) - 1);
    end
    
    function dst_bn = bn(src_Vm)
      dst_bn = 0.125.*exp(src_Vm./80);
    end
  end
end

ナトリウムチャネル

最後にナトリウムチャネルの実装を行います。
ナトリウムチャネルは $m$, $h$ の2種類のゲートを持ちます。

カリウムチャネルのコンダクタンスは以下のように与えられます。

$$
G_{Na} = G_{Na}^{\text{max}} {P_{m}}^3 {P_{h}}
$$

膜電位に依存する状態遷移速度はカリウムイオンチャネルと同様に、それぞれ以下のように求められました。

\begin{aligned}
\alpha_m (V) &= \frac{0.1 (V + 25)} {\exp{\left(\frac{V + 25}{10}\right)} - 1} &
\beta_m  (V) &= 4 \exp{\left(\frac{V}{18}\right)} \\
\alpha_h (V) &= 0.07 \exp{\left(\frac{V}{20}\right)} &
\beta_h  (V) &= \frac{1}{\exp{\left(\frac{V+30}{10}\right)} + 1} \\
\end{aligned}

ナトリウムチャネルのソースコード

２種類分のゲートを実装しないといけない点を除けば、カリウムチャネルと同じです。

:point_right: %%%%%%%%【コードを表示する】%%%%%%%% :point_left: ◀
ChSodium.m
% Sodium(Natrium) Ion Channel
classdef ChSodium < Channel
  properties
    Pm
    Ph
  end
  
  methods
    function obj = ChSodium(src_Vm)
      obj = obj@Channel(120, -115);
      obj.Pm = obj.am(src_Vm)./(obj.am(src_Vm) + obj.bm(src_Vm));
      obj.Ph = obj.ah(src_Vm)./(obj.ah(src_Vm) + obj.bh(src_Vm));
    end
    
    function dst_Ix = Ix(obj, src_Vm)
      dst_Ix = obj.Ix@Channel(src_Vm).*obj.Pm.^3.*obj.Ph;
    end
    
    function dst_dm = dm(obj, src_Vm, del_Pm)
      tmp_Pm = obj.Pm + del_Pm;
      dst_dm = obj.am(src_Vm).*(1 - tmp_Pm) - obj.bm(src_Vm).*tmp_Pm;
    end
    
    function dst_dh = dh(obj, src_Vm, del_Ph)
      tmp_Ph = obj.Ph + del_Ph;
      dst_dh = obj.ah(src_Vm).*(1 - tmp_Ph) - obj.bh(src_Vm).*tmp_Ph;
    end
  end
  
  methods(Static)
    function dst_am = am(src_Vm)
      dst_am = 0.1.*(src_Vm + 25)./(exp((src_Vm + 25)/10) - 1);
    end
    
    function dst_bm = bm(src_Vm)
      dst_bm = 4.*exp(src_Vm./18);
    end
    
    function dst_ah = ah(src_Vm)
      dst_ah = 0.07.*exp(src_Vm./20);
    end
    
    function dst_bh = bh(src_Vm)
      dst_bh = 1./(exp((src_Vm + 30)./10) + 1);
    end
  end
end

HHM クラスの実装

Hodgkin-Huxley Model

さて、以上でチャネルの用意はできました。
ここではこれらのチャネルをまとめて、細胞膜の電気活動の変化を模倣するモデルを構築します。

:point_right: %%%%%%%%【詳細説明を表示する】%%%%%%%% :point_left: ◀

HHM では、細胞膜に流れる全電流 $I$ を脂質2重膜の電荷 $Q_M$ の変化（つまり$I_c = \frac{dQ_M}{dt}$）とイオン電流 $I_i$ に分割します。細胞膜の静電容量 $C_M$ と電圧 $V$ を用いると電荷 $Q_M$ は $Q_M = C_M V$ と表すことができ、これらをまとめると以下のように表記できます。

I = C_M \frac{dV}{dt} + I_i

イオン電流 $I_i$ はナトリウムチャネル、カリウムチャネル、リークチャネル３つの異なるチャネルにおけるイオン電流（$I_{Na}, I_K, I_l$）に分割することができます。

$$
I_i = I_{Na} + I_K + I_l
$$

また、それぞれのチャネルのイオン電流は、コンダクタンス $G_X$ と起電力 $V_X$ をもちいて以下のように表記できます。

\begin{aligned}
I_{Na} &= G_{Na} (V - V_{Na}) \\
I_{K}  &= G_{K}  (V - V_{K})  \\
I_{l}  &= G_{l}  (V - V_{l})  \\
\end{aligned}

これらを合わせると、

$$
I = C_M \frac{dV}{dt} + G_{Na} (V - V_{Na}) + G_{K} (V - V_{K}) + G_{L} (V - V_{L})
$$

さらに式を変形することで $V$ についての微分方程式を得ます。

\frac{dV}{dt} = - \frac{1}{C_M} \left( G_{Na} (V - V_{Na}) + G_{K}  (V - V_{K}) + G_{L}  (V - V_{L}) - I \right)

この微分方程式を用いることで、全膜電流 $I$ の大きさから 電圧の時間発展を求めることができます。

さて、ここまで説明を避けていましたが、Hodgkin and Huxley (1952d) における電圧 $V$ と電流 $I$ の定義は現在、一般的に使用されている定義とは異なります。

Hodgkin and Huxley (1952d)では電圧は、静止膜電位を基準（0V）にしており、脱分極へ向かう方向を負としています。この膜電位の定義は現在一般的な細胞外を基準（0V）にし、脱分極の方向を正とするものとするものとはことなります。
また、電流に関しても、正電荷が膜を外から内へと流れる方向を正としています。これも正電荷が膜を外から内へと流れる方向を負とする現代のものとは正反対です。

Hodgkin and Huxley (1952d) 現代の定義 変換
電圧 $V$ "depolarization negative" $E_M$ "depolarization positive" $E_M = E_R - V$
電流 $I$ "inward current positive" $I_S$ "outward current positive" $I_S = -I$

これまでの実装は、論文のオリジナルの定義に沿ってきましたが、数値計算した結果は、普段慣れ親しんだ定義で見たいものです。

電圧に関しては以下の式を用いることで $V$ から膜電位 $E_M$ を求めます。ただし、$E_R$ は細胞外を基準(0V)とした時の細胞内の静止電位（静止膜電位：$E_R \approx -65\space\text{mV}$）です。

E_M = E_R - V

また、刺激電流にも注意が必要です。刺激としての細胞内への正電荷の注入（正の注入電流）は、膜において内側から外側への正電荷の流れ（膜電流）を作ります。論文の膜電流の定義ではこれは負の電流となりますが、現代の定義では正の膜電流となります。刺激電流（$I_S$）からモデルの内部で扱う膜電流（$I$）への変換は以下の式で与えられます。

I = - I_S

プロパティの実装

以下の様に必要なプロパティを定義します。
HHM クラスはプロパティとして３つのチャネルを定義します。

HHM.m
classdef HHM
  properties
    Cm = 1.0
    Vm
    Er
    
    chN
    chK
    chL
  end

依存プロパティの実装

さて、細胞外を基準とした膜電位 $E_M$ は、モデルの電圧 $V$ （プロパティ Vm） と静止膜電位 $E_R$ （プロパティ Er）に依存して決定されます。ここではEmを依存プロパティとして設定することにより、必要な時にその値を計算して使用できるようにします。

依存プロパティはそれ自体は値を持たず、アクセスされた際に対になるアクセスメソッドによって値を返します。今回 プロパティ Em を依存プロパティとして定義する場合、getメソッド get.Em() を定義する必要があります。

HHM.m
classdef HHM
  ...
  properties (Dependent)
    Em
  end
  
  ...
  
  methods
    ...
    function dst_Em = get.Em(obj)
      dst_Em = obj.Er - obj.Vm;
    end
    ...
  end
end

詳しくは、依存プロパティを参照してください。

メソッドの実装

まずはコンストラクターと微分方程式をメソッドに定義します。

コンストラクターの実装

コンストラクターでは、静止膜電位の大きさを指定できるようにします。
今回は真面目に nargin を用いて引数の数をチェックして、引数が与えられているときのみ静止膜電位の値を設定する様にします。（引数がない場合は-65を与えます。）

また、３つのチャネルに関してもここで初期化しておきます。

HHM.m
classdef HHM
  ...
  methods
    function obj = HHM(src_Er)      
      if 0 < nargin 
        obj.Er = src_Er;
      else
        obj.Er = -65;
      end
      
      obj.Vm = 0;
      obj.chN = ChSodium(obj.Vm);
      obj.chK = ChPotassium(obj.Vm);
      obj.chL = ChLeak();
    end
    ...
  end
end

微分方程式の実装

電圧の時間変化は、以下の式で与えられます。

\frac{dV}{dt} = - \frac{1}{C_M} \left( G_{Na} (V - V_{Na}) + G_{K}  (V - V_{K}) + G_{L}  (V - V_{L}) - I \right)

膜電流 $I$ の大きさは外部から与えることのできるよう、引数（src_Im）として定義します。
(ただし、ここでのsrc_Imの電流の方向は論文中で定義された方向です。)

また、今回も数値計算に備えて $V$ の微小変動 $\delta V $ を加えることのできるように準備しておきます。

HHM.m
classdef HHM
  ...
  methods
    ...
    function dst_dV = dV(obj, src_Im, del_Vm)
      tmp_Vm = obj.Vm + del_Vm;
      dst_dV = -(obj.chN.Ix(tmp_Vm) + obj.chK.Ix(tmp_Vm) + obj.chL.Ix(tmp_Vm) - src_Im)./obj.Cm;
    end
    ...
  end
end

微分方程式の数値計算

ここまでで、Hodgkin-Huxley Modelの基本的な実装はできました。
あとは数値計算によって値の時間変化を計算していくだけです。

オイラー法の実装

初めに、オイラー法を確認します。

\begin{aligned}
y(t + h) &= y(t) + h \times f(t, y) + \mathcal{O} (h^2) \\
f(t, y) &= \frac{d}{dt} y(t)\\
\end{aligned}

電圧 $V(t)$ に関する微分方程式の実装（HHM.dV()）を用いると、ステップ幅 $h$ = dt 進めた後の電圧 $V(t+dt)$ の数値解は以下のように求めることができます。ここでは、HHMクラスのメソッドとしてオイラー法を実装していきます。

HHM.m
classdef HHM
  ...
  methods
    ...
    function dst_Vm = Euler(obj, src_Im, dt)
      dst_Vm = obj.Vm + dt.*obj.dV(src_Im, 0);
    end
  end
end

さて、ほとんどの場合は１ステップ後の数値解を求めるとともに、クラスの電圧の値を保持するプロパティ（HHM.Vm）も更新したいことと思います。

:point_right: %%%%%%%%【補足説明を表示する】%%%%%%%% :point_left: ◀

試しのプロパティの更新に関して、以下のように実装してもうまく働きません。

HHM.m.wrong
classdef HHM
  ...
  methods
    ...
    function dst_Vm = Euler(obj, src_Im, dt)
      dst_Vm = obj.Vm + dt.*obj.dV(src_Im, 0);
      obj.Vm = dst_Vm; % This code does not work as expected
    end
  end
end

以下のように、クラスを初期化し、1ステップ分実行してもプロパティ（Vm）の値が変更されていないことがわかります。

>> hhm = HHM(-65);
>> hhm.Vm
ans =
     0
>> hhm.Euler(0, 0.01)
ans =
  -4.2237e-05
>> hhm.Vm
ans =
     0  % Expected to be -4.2237e-05

これはmatlabのデフォルトのクラス定義が値クラスと呼ばれるタイプであり、このクラスでは関数内でのオブジェクトの変更内容が関数の外へ自動的に影響を及ぼしません。つまり上の実行例でhhm.Euler(0, 0.01)の実行で関数内の obj に加えられたは関数の外部にある hhm へ影響を及ぼすことができません。この結果、関数の実行後に hhm.Vm で値を確認しても値は元の初期値のままとなります。

メソッド内でのオブジェクトの変更を反映させるためには２つの方法があります。

  1. メソッドの戻り値として変更されたオブジェクトを返す
  2. handleクラスを継承する

1.の方法はコンストラクターの実装のように関数の1つ目の戻り値として obj を返し、更新したいオブジェクトに代入します。

    ...
    function [obj, dst_Vm] = Euler(obj, src_Im, dt)
    ...

>> [hhm, dst_Vm] = hhm.Euler(0, 0.01)

今回は2. のhandleクラスを継承する方法によってメソッド内でのオブジェクトへの変更が元のオブジェクトへも反映するようにします。

以下のように classdef 行において handle クラスを継承します。

HHM.m
classdef HHM < handle
  ...
  methods
    ...
    function dst_Vm = Euler(obj, src_Im, dt)
      dst_Vm = obj.Vm + dt.*obj.dV(src_Im, 0);
      obj.Vm = dst_Vm;
    end
  end
end

詳しくは、[オブジェクトの変更]
(https://jp.mathworks.com/help/matlab/matlab_oop/matlab-vs-other-oo-languages.html#bslvcv1)を参照してください。

さて、HHMは電圧 $V$ の時間変化に関する微分方程式だけでなく、チャネルのゲートの開口確率の時間変化に関する微分方程式も含みます。

関数 deal は、入力をコピーをして、出力に代入します。
ここでは、微分方程式を使って求めた解析解を、HHM.Eular() 関数の出力引数と、オブジェクト（obj）のプロパティの両方に代入しています。

HHM.m
classdef HHM < handle
  ...
  methods
    ...
    function dst_Vm = Euler(obj, src_Im, dt)
      [dst_Vm, obj.Vm] = deal(obj.Vm + dt.*obj.dV(src_Im, 0));
      obj.chN.Pm = obj.chN.Pm + dt.*obj.chN.dm(obj.Vm, 0);
      obj.chN.Ph = obj.chN.Ph + dt.*obj.chN.dh(obj.Vm, 0);
      obj.chK.Pn = obj.chK.Pn + dt.*obj.chK.dn(obj.Vm, 0);
    end
  end
end

古典的ルンゲ=クッタ法の実装

微分方程式の引数として微小変動 $\delta$ を渡せるようにしてありました。これは中点法や古典的ルンゲ=クッタ法を使用する際に使います。

:point_right: %%%%%%%%【古典的ルンゲ=クッタ法の実装】%%%%%%%% :point_left: ◀
HHM.m
classdef HHM < handle
  ...
  methods
    function dst_Vm = RungeKutta4(obj, src_Im, dt)
      k1 = dt.*obj.dV(src_Im, 0);
      k2 = dt.*obj.dV(src_Im, k1/2);
      k3 = dt.*obj.dV(src_Im, k2/2);
      k4 = dt.*obj.dV(src_Im, k3);
      [dst_Vm, obj.Vm] = deal(obj.Vm + (k1 + 2*k2 + 2*k3 + k4)/6);
      
      k1 = dt.*obj.chN.dm(obj.Vm, 0);
      k2 = dt.*obj.chN.dm(obj.Vm, k1/2);
      k3 = dt.*obj.chN.dm(obj.Vm, k2/2);
      k4 = dt.*obj.chN.dm(obj.Vm, k3);
      obj.chN.Pm = obj.chN.Pm + (k1 + 2*k2 + 2*k3 + k4)/6;
      
      k1 = dt.*obj.chN.dh(obj.Vm, 0);
      k2 = dt.*obj.chN.dh(obj.Vm, k1/2);
      k3 = dt.*obj.chN.dh(obj.Vm, k2/2);
      k4 = dt.*obj.chN.dh(obj.Vm, k3);
      obj.chN.Ph = obj.chN.Ph + (k1 + 2*k2 + 2*k3 + k4)/6;
      
      k1 = dt.*obj.chK.dn(obj.Vm, 0);
      k2 = dt.*obj.chK.dn(obj.Vm, k1/2);
      k3 = dt.*obj.chK.dn(obj.Vm, k2/2);
      k4 = dt.*obj.chK.dn(obj.Vm, k3);
      obj.chK.Pn = obj.chK.Pn + (k1 + 2*k2 + 2*k3 + k4)/6;
    end
  end
end

数値計算によるシミュレーション

シミュレーション計算関数

上で実装したオイラー法を用いて、刺激電流の時系列データ（src_Is）と、時間ステップ幅（dt）を入力として、数値計算結果を出力する関数を作成してみます。obj.Euler(-src_Is(t), dt) と刺激電流の向きを逆転させることと、obj.Em を使用して細胞外電位を基準とする膜電位の値を取得する点に注意します。

HHM.m
classdef HHM < handle
  ...
  methods
    ...
    function [dst_Em, dst_Pm, dst_Ph, dst_Pn] = calc(obj, src_Is, dt)
      N = numel(src_Is);
      [dst_Em, dst_Pm, dst_Ph, dst_Pn] = deal(zeros(size(src_Is)));
      
      for n = 1:N
        obj.Euler(-src_Is(n), dt);
        dst_Em(n) = obj.Em;
        dst_Pm(n) = obj.chN.Pm;
        dst_Ph(n) = obj.chN.Ph;
        dst_Pn(n) = obj.chK.Pn;
      end
    end
  end
end

適当な電流刺激を設定して、数値計算を実行します。
ここでは大きさの異なる電流をそれぞれ 100 ms の長さ注入した場合の膜電位の変化を計算します。

>> hhm = HHM(-65);
>> dt = 0.05;
>> src_Tm = colon(0, dt, 900-dt)';
>> src_Im = zeros(size(src_Tm));
>> src_Im(100<=src_Tm & src_Tm<200) = 2;
>> src_Im(300<=src_Tm & src_Tm<400) = 4;
>> src_Im(500<=src_Tm & src_Tm<600) = 6;
>> src_Im(700<=src_Tm & src_Tm<800) = 8;

>> [dst_Em, dst_Pm, dst_Ph, dst_Pn] = hhm.calc(src_Im, dt);

>> subplot(311);
>> area(src_Tm, src_Im, 'FaceAlpha', 0.1);
>> set(gca, {'Box', 'TickDir'}, {'off', 'out'});

>> subplot(312);
>> plot(src_Tm, dst_Em, 'LineWidth', 1);
>> set(gca, {'Box', 'TickDir'}, {'off', 'out'});

>> subplot(313);
>> hold on
>> plot(src_Tm, dst_Pm, 'LineWidth', 1);
>> plot(src_Tm, dst_Ph, 'LineWidth', 1);
>> plot(src_Tm, dst_Pn, 'LineWidth', 1);
>> set(gca, {'Box', 'TickDir'}, {'off', 'out'});

電流の大きさ（上）に従ってが膜電位（下）が変化する様子をみることができます。

image.png

描画更新と同期したシミュレーション

matlabの描画機能と組み合わせることによって、さらに凝った実装もできます。

以下の例では標準正規分布から生成された疑似乱数と正弦波を、異なる割合で足した入力に対する出力をシミュレーションしています。
（正弦波の振幅は分散が１になるように $\sqrt{2}$ 倍に補正します。）

:point_right: %%%%%%%%【コードを表示する】%%%%%%%% :point_left: ◀
HHM.m
classdef HHM < handle
  ...
  methods
    ...
    function run(obj)
      T = 1200; % Trial Time    Set Smaller to Get Faster
      S = 5;    % Sample Step   Set Larger  to Get Faster
      R = 4;    % Refresh Rate  Set Larger  to Get Faster
      
      dt = 0.05;
      src_Tm = colon(0, dt, T-dt)';
      N = length(src_Tm);

      src_Is = zeros(N, 1);
      dst_Pm = zeros(N, 1) + obj.chN.Pm;
      dst_Ph = zeros(N, 1) + obj.chN.Ph;
      dst_Pn = zeros(N, 1) + obj.chK.Pn;      
      dst_Em = zeros(N, 1) + obj.Em;

      fig = clf;

      subplot(311);
      plt_Im = plot(src_Tm(1:S:N), src_Is(1:S:N));
      set(gca, {'Box', 'TickDir', 'YLim'}, {'off', 'out', [-40, 40]});

      subplot(312);
      plt_Em = plot(src_Tm(1:S:N), dst_Em(1:S:N));
      set(gca, {'Box', 'TickDir', 'YLim'}, {'off', 'out', [-100, 50]});
      

      subplot(313);
      hold on
      plt_Pm = plot(src_Tm(1:S:N), dst_Pm(1:S:N));
      plt_Ph = plot(src_Tm(1:S:N), dst_Ph(1:S:N));
      plt_Pn = plot(src_Tm(1:S:N), dst_Pn(1:S:N));
      set(gca, {'Box', 'TickDir', 'YLim'}, {'off', 'out', [0, 1]});     

      n = 0;
      m = 1;
      M = 4;
      while isgraphics(fig)
        n = mod(n, N) + 1;

        if n == N
          m = mod(m, M) + 1;
        end
        switch (m)
          case {1}
            src_Is(n) = 8.*randn;
          case {2}
            src_Is(n) = 6.*randn + 2.*sqrt(2).*sin(2.*pi.*src_Tm(n)./100);
          case {3}
            src_Is(n) = 4.*randn + 4.*sqrt(2).*sin(2.*pi.*src_Tm(n)./100);
          case {4}
            src_Is(n) = 2.*randn + 6.*sqrt(2).*sin(2.*pi.*src_Tm(n)./100);
          otherwise
            src_Is(n) =  0;
        end
        [dst_Em(n), dst_Pm(n), dst_Ph(n), dst_Pn(n)] = obj.calc(src_Is(n), dt);        
        
        if mod(n, S*R) == 1
          plt_Im.YData = src_Is([n:S:N, 1:S:n-1]);
          plt_Em.YData = dst_Em([n:S:N, 1:S:n-1]);
          plt_Pm.YData = dst_Pm([n:S:N, 1:S:n-1]);
          plt_Ph.YData = dst_Ph([n:S:N, 1:S:n-1]);
          plt_Pn.YData = dst_Pn([n:S:N, 1:S:n-1]);
          drawnow
        end
      end
    end
  end
end

描画の速度は環境に依存するので、以下の3つの変数の値を調節してください

T = 1200; % Trial Time    Set Smaller to Get Faster
S = 5;    % Sample Step   Set Larger  to Get Faster
R = 4;    % Refresh Rate  Set Larger  to Get Faster

run.gif

:point_right: %%%%%%%%【HHM.m のコード】%%%%%%%% :point_left: ◀
HHM.m
classdef HHM < handle
  properties
    Cm = 1.0
    Vm
    Er
    
    chN
    chK
    chL
  end
  
  properties (Dependent)
    Em
  end
    
  methods
    function obj = HHM(src_Er)      
      if 0 < nargin 
        obj.Er = src_Er;
      else
        obj.Er = 0;
      end
      
      obj.Vm = 0;
      obj.chN = ChSodium(obj.Vm);
      obj.chK = ChPotassium(obj.Vm);
      obj.chL = ChLeak();
    end
    
    function dst_dV = dV(obj, src_Im, del_Vm)
      tmp_Vm = obj.Vm + del_Vm;
      dst_dV = -(obj.chN.Ix(tmp_Vm) + obj.chK.Ix(tmp_Vm) + obj.chL.Ix(tmp_Vm) - src_Im)./obj.Cm;
    end
    
    function dst_Em = get.Em(obj)
      dst_Em = obj.Er - obj.Vm;
    end
    
    function dst_Vm = Euler(obj, src_Im, dt)
      [dst_Vm, obj.Vm] = deal(obj.Vm + dt.*obj.dV(src_Im, 0));
      obj.chN.Pm = obj.chN.Pm + dt.*obj.chN.dm(obj.Vm, 0);
      obj.chN.Ph = obj.chN.Ph + dt.*obj.chN.dh(obj.Vm, 0);
      obj.chK.Pn = obj.chK.Pn + dt.*obj.chK.dn(obj.Vm, 0);
    end
    
    function dst_Vm = RungeKutta4(obj, src_Im, dt)
      k1 = dt.*obj.dV(src_Im, 0);
      k2 = dt.*obj.dV(src_Im, k1/2);
      k3 = dt.*obj.dV(src_Im, k2/2);
      k4 = dt.*obj.dV(src_Im, k3);
      [dst_Vm, obj.Vm] = deal(obj.Vm + (k1 + 2*k2 + 2*k3 + k4)/6);
      
      k1 = dt.*obj.chN.dm(obj.Vm, 0);
      k2 = dt.*obj.chN.dm(obj.Vm, k1/2);
      k3 = dt.*obj.chN.dm(obj.Vm, k2/2);
      k4 = dt.*obj.chN.dm(obj.Vm, k3);
      obj.chN.Pm = obj.chN.Pm + (k1 + 2*k2 + 2*k3 + k4)/6;
      
      k1 = dt.*obj.chN.dh(obj.Vm, 0);
      k2 = dt.*obj.chN.dh(obj.Vm, k1/2);
      k3 = dt.*obj.chN.dh(obj.Vm, k2/2);
      k4 = dt.*obj.chN.dh(obj.Vm, k3);
      obj.chN.Ph = obj.chN.Ph + (k1 + 2*k2 + 2*k3 + k4)/6;
      
      k1 = dt.*obj.chK.dn(obj.Vm, 0);
      k2 = dt.*obj.chK.dn(obj.Vm, k1/2);
      k3 = dt.*obj.chK.dn(obj.Vm, k2/2);
      k4 = dt.*obj.chK.dn(obj.Vm, k3);
      obj.chK.Pn = obj.chK.Pn + (k1 + 2*k2 + 2*k3 + k4)/6;
    end
    
    function [dst_Em, dst_Pm, dst_Ph, dst_Pn] = calc(obj, src_Is, dt)
      N = numel(src_Is);
      [dst_Em, dst_Pm, dst_Ph, dst_Pn] = deal(zeros(N, 1));
      
      for n = 1:N
        obj.Euler(-src_Is(n), dt);
        dst_Em(n) = obj.Em;
        dst_Pm(n) = obj.chN.Pm;
        dst_Ph(n) = obj.chN.Ph;
        dst_Pn(n) = obj.chK.Pn;
      end
    end
    
    function run(obj)
      T = 1200; % Trial Time    Set Smaller to Get Faster
      S = 5;    % Sample Step   Set Larger  to Get Faster
      R = 4;    % Refresh Rate  Set Larger  to Get Faster
      
      dt = 0.05;
      src_Tm = colon(0, dt, T-dt)';
      N = length(src_Tm);

      src_Is = zeros(N, 1);
      dst_Pm = zeros(N, 1) + obj.chN.Pm;
      dst_Ph = zeros(N, 1) + obj.chN.Ph;
      dst_Pn = zeros(N, 1) + obj.chK.Pn;      
      dst_Em = zeros(N, 1) + obj.Em;

      fig = clf;

      subplot(311);
      plt_Im = plot(src_Tm(1:S:N), src_Is(1:S:N));
      set(gca, {'Box', 'TickDir', 'YLim'}, {'off', 'out', [-40, 40]});

      subplot(312);
      plt_Em = plot(src_Tm(1:S:N), dst_Em(1:S:N));
      set(gca, {'Box', 'TickDir', 'YLim'}, {'off', 'out', [-100, 50]});
      
      subplot(313);
      hold on
      plt_Pm = plot(src_Tm(1:S:N), dst_Pm(1:S:N));
      plt_Ph = plot(src_Tm(1:S:N), dst_Ph(1:S:N));
      plt_Pn = plot(src_Tm(1:S:N), dst_Pn(1:S:N));
      set(gca, {'Box', 'TickDir', 'YLim'}, {'off', 'out', [0, 1]});
      
      drawnow

      n = 0;
      m = 1;
      M = 4;
      while isgraphics(fig)
        n = mod(n, N) + 1;

        if n == N
          m = mod(m, M) + 1;
        end
        switch (m)
          case {1}
            src_Is(n) = 8.*randn;
          case {2}
            src_Is(n) = 6.*randn + 2.*sqrt(2).*sin(2.*pi.*src_Tm(n)./100);
          case {3}
            src_Is(n) = 4.*randn + 4.*sqrt(2).*sin(2.*pi.*src_Tm(n)./100);
          case {4}
            src_Is(n) = 2.*randn + 6.*sqrt(2).*sin(2.*pi.*src_Tm(n)./100);
          otherwise
            src_Is(n) =  0;
        end
        [dst_Em(n), dst_Pm(n), dst_Ph(n), dst_Pn(n)] = obj.calc(src_Is(n), dt);        
        
        if  mod(n, S*R) == 1
          plt_Im.YData = src_Is([n:S:N, 1:S:n-1]);
          plt_Em.YData = dst_Em([n:S:N, 1:S:n-1]);
          plt_Pm.YData = dst_Pm([n:S:N, 1:S:n-1]);
          plt_Ph.YData = dst_Ph([n:S:N, 1:S:n-1]);
          plt_Pn.YData = dst_Pn([n:S:N, 1:S:n-1]);
          drawnow
        end
      end
    end
  end
end

まとめ

今回はHodgkin-Huxley Modelを作成しながら、matlabのクラスの機能・実装方法を確認してみました。

クラスを用いることで、パラメータ（各チャネルのパラメータや、膜電位などの物理量）とそれに対する操作を関連付けて実装することができます。

正直なところ、Hodgkin-Huxley Model程度の実装量であればクラスのありがたみを実感することは難しいかもしれませんが、複雑な解析を実装する枠組みとして、クラスが役に立つ場面があるかと思います。

  1. Hodgkin AL and Huxley AF. (1952)
    A quantitative description of membrane current and its application to conduction and excitation in nerve.
    J Physiol 117:500-44
    https://doi.org/10.1113/jphysiol.1952.sp004764
    PMID: 12991237
    出版当時、チャネルやゲートの存在は明らかになっていませんでしたが（これがこの研究のすごいところではあるのですが）、本記事では説明の便宜上チャネル、ゲートなどの語を使用します。 ↩

  2. 詳しくはMATLAB でのオブジェクト指向デザインの２項目にまとめられています。
    オブジェクト指向の設計を使用する理由
    クラスの MATLAB での役割 ↩

  3. コンダクタンスは抵抗の逆数で、電流の流れやすさを表します。単位はジーメンス[S] またはモー[℧]です。モーと言っても記号の形が牛のように見えるからではなく、Ω【ohm】 を逆から読んだだけです。(:cow: !mho) ↩

  4. 特に決まりはないのですが、慣例的にクラスを定義するコートの中ではクラスオブジェクトとして obj を使用することが多いです。もちろんほかの名前を用いることもできます。 ↩

  5. matlabではコンストラクターが入力引数がない場合でも実行できるように実装することを推奨しています。この場合、 nargin を利用して入力引数によって処理を切り替える必要があります。入力引数なしでコンストラクターを呼び出す条件を参照してください。 ↩

  6. 念のため、日本で広く使われているカリウムはドイツ語/ラテン語（Kalium）から来ており、英名はPotassiumです。ただし、元素記号はKalium に由来する $K$ となっています。大変ややこしいです。 ↩

  7. 論文中にはゲートの開口確率という概念は出てきません。論文の式の中では単純に変数 $n$ として表記されていますが、ここでは表記上の見通しをよくするため、$n$ を $P$ に添えて表記します。 ↩

9
5
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
9
5

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?