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

<!--
新型コロナウイルス感染症（COVID-19）を受けてMathWorks社は**まだキャンパスライセンスを契約していない大学でも** 大学の学生・教職員向けに期間付きでライセンスを開放してくれているようです。
> COVID-19 Alert: To enable remote work/online classes, MATLAB and add-on products are available to use through August 31, 2020.  Download to your computer or access MATLAB Online and MATLAB Drive from your browser.

- ライセンスは2020年8月31日までの期間限定です。
- ****-u.ac.jp などの大学ドメインのついたメールアドレスでユーザー登録する必要があるようです。

自分の大学が対応しているかは以下のサイトから確認できます。
https://jp.mathworks.com/academia/tah-support-program/eligibility.html


または各大学のポータルサイトでも確認できます。（以下は大阪大学の例です。）
https://jp.mathworks.com/academia/tah-portal/osaka-university-31487301.html
大学の英語表記 + MathWorks などで検索すれば引っかかると思います。[^covid-19]

[^covid-19]: 他にもいくつかの大学でライセンスの提供があります。（確認したところのみ）<br>大阪大学　　https://jp.mathworks.com/academia/tah-portal/osaka-university-31487301.html<br>九州大学　　https://jp.mathworks.com/academia/tah-portal/kyushu-university-31486098.html<br>京都大学　　https://jp.mathworks.com/academia/tah-portal/kyoto-university-31485310.html<br>東北大学　　https://jp.mathworks.com/academia/tah-portal/tohoku-university-31485743.html<br>名古屋大学　https://jp.mathworks.com/academia/tah-portal/nagoya-university-31486892.html

<!-- Available to use through August 31, 2020 -->


# はじめに

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

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


## [Hodgkin-Huxley方程式] とは

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

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



## [クラス] とは

クラスとは、データとそのデータに関連する処理・操作を関連づけて定義したものです。
プログラミング言語ごとにその思想は異なりますが、matlabにおけるクラスの考え方は「[MATLAB でのオブジェクト指向デザイン]」を確認してください。[^class]  
[^class]: 詳しくは[MATLAB でのオブジェクト指向デザイン]の２項目にまとめられています。<br>[オブジェクト指向の設計を使用する理由]<br>[クラスの MATLAB での役割]


# チャネル クラスの実装

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

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

```matlab:Channel.m
classdef Channel
  ...
end
```

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

<details>
<summary>:point_right: <code><b>%%%%%%%%【詳細説明を表示する】%%%%%%%%</b></code> :point_left: ◀</summary><div>

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

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

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


```math
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)
$$

</div></details>

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

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

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

## プロパティの実装

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

```matlab:Channel.m
classdef Channel
  properties
    Gx
    Vx
  end
end
```

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

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

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

```matlab
  properties
    PI =  3.141592653589793
  end
```

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

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


## メソッドの実装

クラスの機能は[Methods ブロック][メソッドと関数]の中に実装します。

### [コンストラクター]の実装

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

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

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

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

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

```matlab
>> chX = Channel();
```

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

[^constArg]: matlabではコンストラクターが入力引数がない場合でも実行できるように実装することを推奨しています。この場合、 [`nargin`] を利用して入力引数によって処理を切り替える必要があります。[入力引数なしでコンストラクターを呼び出す条件](https://jp.mathworks.com/help/matlab/matlab_oop/class-constructor-methods.html#btn2kiy)を参照してください。


```matlab: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つのプロパティに値が渡されるのがわかります。


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

### イオン電流関数の実装


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

```math
I_X = G_X \left( V - V_X \right)
```

```matlab: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()`)を用いて実装した関数を呼び出します。

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

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

```matlab
>> chX.Gx
ans =
     5

>> chX.Vx
ans =
     7

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

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

<details>
<summary>:point_right: <code><b>%%%%%%%%【コードを表示する】%%%%%%%%</b></code> :point_left: ◀</summary><div>

```matlab: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
```

</div></details>


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

一番実装が簡単なものから始めましょう。
[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`）を用意して以下のとおりに記述します。

```matlab:ChLeak.m
classdef ChLeak < Channel
end
```

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

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

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

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

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

```matlab
SubClassObject@SuperClassName(...)

% Example
obj@Channel(...)
```

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

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

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

詳細は、[サブクラス コンストラクター](https://jp.mathworks.com/help/matlab/matlab_oop/class-constructor-methods.html#brd2m9e-6)を参照してください。

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

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



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

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

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

```matlab: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)$ を用いて以下のように表記します。[^gate]

[^gate]: 論文中にはゲートの開口確率という概念は出てきません。論文の式の中では単純に変数 $n$ として表記されていますが、ここでは表記上の見通しをよくするため、$n$ を $P$ に添えて表記します。
 
$$
G_{K} = G_{K}^{max} {P_{n}(V)}^4
$$

## ゲートの2状態モデル

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


<details>
<summary>:point_right: <code><b>%%%%%%%%【詳細説明を表示する】%%%%%%%%</b></code> :point_left: ◀</summary><div>

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

```math
\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}$ の時間遷移は以下の微分方程式で表すことができます。

```math
\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)$ を考えると以下のようになります。

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

</div></details>

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

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

### ゲート開口確率の追加

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

```matlab:ChPotassium.m
classdef ChPotassium < Channel
  properties
    Pn
  end
  ...
end
```

### 状態遷移速度の追加

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

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

```math
\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 ブロック][メソッドと関数]のみです。

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

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

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


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

```matlab: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
```

また、クラスの実体を伴わずにメソッドの実行を行うことができます。
　
<details>
<summary>:point_right: <code><b>%%%%%%%%【コードを表示する】%%%%%%%%</b></code> :point_left: ◀</summary><div>

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

```matlab
>> 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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/159851/c417c835-fb73-ecd6-cfeb-b0ade189781f.png)

</div></details>


### [コンストラクター]の変更

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

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

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


```matlab: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`）では以下のように定義されています。


```matlab: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`）の **同名の** メソッドを呼び出す際は`@`を用いて、以下の様に記述します。

```matlab
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$ （の４乗）で補正します。

```matlab: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$ の時間変化は以下の微分方程式で表記できます。

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

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

```matlab: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
```

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

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

<details>
<summary>:point_right: <code><b>%%%%%%%%【コードを表示する】%%%%%%%%</b></code> :point_left: ◀</summary><div>

```matlab: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
```

</div></details>

# ナトリウムチャネル

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

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

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


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

```math
\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}
```

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

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

<details>
<summary>:point_right: <code><b>%%%%%%%%【コードを表示する】%%%%%%%%</b></code> :point_left: ◀</summary><div>

```matlab: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
```

</div></details>


# HHM クラスの実装

## Hodgkin-Huxley Model

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

<details>
<summary>:point_right: <code><b>%%%%%%%%【詳細説明を表示する】%%%%%%%%</b></code> :point_left: ◀</summary><div>

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

```math
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$ をもちいて以下のように表記できます。

```math
\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$ についての微分方程式を得ます。

```math
\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$ の大きさから 電圧の時間発展を求めることができます。

</div></details>

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

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

|| [Hodgkin and Huxley (1952d)][HH(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}$）です。 

```math
E_M = E_R - V
```

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

```math
I = - I_S
```


## プロパティの実装

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

```matlab: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()` を定義する必要があります。


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

詳しくは、[依存プロパティ][依存プロパティの set および get メソッド]を参照してください。


## メソッドの実装

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

### コンストラクターの実装

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

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

```matlab: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
```

### 微分方程式の実装

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

```math
\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 $ を加えることのできるように準備しておきます。

```matlab: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の基本的な実装はできました。
あとは数値計算によって値の時間変化を計算していくだけです。

### オイラー法の実装

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

```math
\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クラスのメソッドとしてオイラー法を実装していきます。

```matlab: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`）も更新したいことと思います。

<details>
<summary>:point_right: <code><b>%%%%%%%%【補足説明を表示する】%%%%%%%%</b></code> :point_left: ◀</summary><div>

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

```matlab: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`）の値が変更されていないことがわかります。

```matlab
>> 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` で値を確認しても値は元の初期値のままとなります。

</div></details>

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

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

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

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

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


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

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

```matlab: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`）のプロパティの両方に代入しています。

```matlab: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$ を渡せるようにしてありました。これは中点法や古典的ルンゲ=クッタ法を使用する際に使います。

<details>
<summary>:point_right: <code><b>%%%%%%%%【古典的ルンゲ=クッタ法の実装】%%%%%%%%</b></code> :point_left: ◀</summary><div>

```matlab: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
```

</div></details>


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

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

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

```matlab: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 の長さ注入した場合の膜電位の変化を計算します。

```matlab
>> 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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/159851/052b4a93-4f21-4069-0a2f-946a7642430c.png)

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

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

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


<details>
<summary>:point_right: <code><b>%%%%%%%%【コードを表示する】%%%%%%%%</b></code> :point_left: ◀</summary><div>

```matlab: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つの変数の値を調節してください

```matlab
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
```

</div></details>

![run.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/159851/1c8e0c31-d5f3-12f0-1d53-7d384098d3f6.gif)

<details>
<summary>:point_right: <code><b>%%%%%%%%【HHM.m のコード】%%%%%%%%</b></code> :point_left: ◀</summary><div>

```matlab: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
```

</div></details>

# まとめ

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

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

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

[HH(1952d)]:
https://doi.org/10.1113/jphysiol.1952.sp004764


[Hodgkin-Huxley方程式]:
https://bsd.neuroinf.jp/wiki/Hodgkin-Huxley%E6%96%B9%E7%A8%8B%E5%BC%8F

[コンダクタンス]:
https://ja.wikipedia.org/wiki/%E3%82%B3%E3%83%B3%E3%83%80%E3%82%AF%E3%82%BF%E3%83%B3%E3%82%B9


[クラス]:
https://jp.mathworks.com/help/matlab/object-oriented-programming.html

[Classdef ブロック]:
https://jp.mathworks.com/help/matlab/matlab_oop/classdef-block.html

[クラス プロパティ]:
https://jp.mathworks.com/help/matlab/matlab_oop/properties.html

[プロパティの定義]:
https://jp.mathworks.com/help/matlab/matlab_oop/specifying-properties.html

[依存プロパティの set および get メソッド]:
https://jp.mathworks.com/help/matlab/matlab_oop/access-methods-for-dependent-properties.html

[メソッドと関数]:
https://jp.mathworks.com/help/matlab/matlab_oop/specifying-methods-and-functions.html

[コンストラクター]:
https://jp.mathworks.com/help/matlab/matlab_oop/class-constructor-methods.html

[通常のメソッド]:
https://jp.mathworks.com/help/matlab/matlab_oop/ordinary-methods.html

[静的メソッド]:
https://jp.mathworks.com/help/matlab/matlab_oop/static-methods.html#responsive_offcanvas

[メソッドの属性]:
https://jp.mathworks.com/help/matlab/matlab_oop/method-attributes.html

[クラスの階層 — 概念]:
https://jp.mathworks.com/help/matlab/matlab_oop/hierarchies-of-classes-concepts.html

[サブクラス オブジェクトでのスーパークラス メソッドの呼び出し]:
https://jp.mathworks.com/help/matlab/matlab_oop/calling-superclass-methods-on-subclass-objects.html


[MATLAB でのオブジェクト指向デザイン]:
https://jp.mathworks.com/help/matlab/object-oriented-design-with-matlab.html

[オブジェクト指向の設計を使用する理由]:
https://jp.mathworks.com/help/matlab/matlab_oop/why-use-object-oriented-design.html

[クラスの MATLAB での役割]:
https://jp.mathworks.com/help/matlab/matlab_oop/handle-objects.html

[`classdef`]:
https://jp.mathworks.com/help/matlab/ref/classdef.html

[`nargin`]:
https://jp.mathworks.com/help/matlab/ref/nargin.html

[`fplot`]:
https://jp.mathworks.com/help/matlab/ref/fplot.html

[`handle`]:
https://jp.mathworks.com/help/matlab/ref/handle-class.html

[`deal`]:
https://jp.mathworks.com/help/matlab/ref/deal.html