zundairo.net
zundairo.net > zunda framework > チュートリアル

zunda framework チュートリアル

さいしょに

まず，あたらしいプロジェクトを作成して，フレームワークのDLLかプロジェクトを参照に追加します。このチュートリアルでは，Hodgkin-Huxleyモデルを構築します。

膜電位V は，膜容量CMと入力電流IM，ナトリウム電流INa，カリウム電流IK，リーク電流IL によって次式で表されます。INa，IK，IL は，それぞれ最大コンダクタンスgと平衡電位VNa，VK，VL，また，INaの活性化パラメータm，不活性化パラメータh，IKの活性化パラメータnによって表されます。φは温度定数を表します。

数式１

数式２

チャネルを作成する

リークチャネル電流

まず，リークチャネルを作成します。チャネルはCurrentChannelクラスを継承して作ります。

リークチャネルは，コンダクタンスと平衡電位で次のように記述できますIl = -gl (Vm - Vl)。ここでは，コンダクタンスは0.3mS/cm2，平衡電位は-49.4mVとします。チャネルで流れる電流はGetCurrentメソッドをオーバーライドすることで指定します。例では，Enableを外部で指定してチャネルの有効化を設定できるようにしてあります。

/// <summary>
/// リーク電流
/// </summary>
public class LeakCurrent : CurrentChannel {
    /// <summary>
    /// LeakCurrentインスタンスを初期化します．
    /// </summary>
    /// <param name="selfNeuron"></param>
    public LeakCurrent(Neuron selfNeuron)
        : base(selfNeuron) {
        Conductance = 0.3;
        ReversalPotential = -49.4;
    }

    /// <summary>
    /// リーク電流
    /// </summary>
    /// <returns></returns>
    public override double GetCurrent() {
        if(Enable) {
            return -Conductance * (SelfNeuron.Potential - ReversalPotential);
        } else {
            return 0.0;
        }
    }
}

カリウムチャネル電流

つぎにカリウムチャネルのクラスを作ります。これも同じようにCurrentChanelクラスから継承してチャネルを作ります。

カリウムチャネルの電流は，平衡電位VKと膜電位依存性活性化変数nによって式(3)によって表されます。ここでは，nは時間で微分する必要があるのでNeuronクラスにある微分変数リストへ式と共に登録する必要があります。登録はAddODEメソッドで行い，ODEActivateKのメソッドを(デリゲートまはたラムダ式でもよい)を初期値と名前とともに渡します。

あとは，リーク電流と同じようにGetCurrentメソッドをオーバーライドしてカリウム電流を取得する式を記述します

public class KCurrent : CurrentChannel {
    private int activateIdx;

    /// <summary>
    /// インスタンスの初期化．
    /// </summary>
    /// <param name="selfNeuron"></param>
    public KCurrent(Neuron selfNeuron)
        : base(selfNeuron) {
        Conductance = 36.0;
        ReversalPotential = -72.0;

        this.activateIdx = selfNeuron.Length;
        selfNeuron.AddODE(ODEActivateK, 0.0, "ActivateK");

        Length = 1;
    }

    /// <summary>
    /// 活性化変数の取得と設定をします．
    /// </summary>
    public double ActivateK {
        get { return SelfNeuron[activateIdx]; }
        set { SelfNeuron[activateIdx] = value; }
    }

    /// <summary>
    /// 微分方程式：K活性化変数n
    /// </summary>
    /// <returns></returns>
    public double ODEActivateK() {
        double alpha, beta;
        alpha = (-0.01 * (SelfNeuron.Potential + 50)) / (-1 + Math.Exp(-(SelfNeuron.Potential + 50) / 10));
        beta = 0.125 * Math.Exp(-(SelfNeuron.Potential + 60) / 80);
        return SelfNeuron.Temperature * (alpha * (1 - ActivateK) - beta * ActivateK);
    }

    /// <summary>
    /// カリウム電流
    /// </summary>
    /// <returns></returns>
    public override double GetCurrent() {
        if(Enable) {
            return -Conductance * ActivateK * ActivateK * ActivateK * ActivateK * (SelfNeuron.Potential - ReversalPotential);
        } else {
            return 0.0;
        }
    }
}

ナトリウムチャネル電流

最後にナトリウムチャネルを作ります。カリウムチャネルと同様に膜電位依存性の活性化変数のmがあります。ただし，ここではもう一つ膜電位依存性の不活性変数も必要です。微分する変数がもう一つ増えたのでコンストラクタで追加する変数も同じように増やします。

public class NaCurrent : CurrentChannel {
    private int activateIdx, inactivateIdx;

    /// <summary>
    /// インスタンスの初期化．
    /// </summary>
    /// <param name="selfNeuron"></param>
    public NaCurrent(Neuron selfNeuron)
        : base(selfNeuron) {
        Conductance = 120.0;
        ReversalPotential = 55.0;

        this.activateIdx = selfNeuron.Length;
        selfNeuron.AddODE(ODEActivateNa, 0.0, "ActivateNa");
        this.inactivateIdx = selfNeuron.Length;
        selfNeuron.AddODE(ODEInactivateNa, 0.0, "InactivateNa");

        Length = 2;
    }

    /// <summary>
    /// 活性化変数の取得と設定．
    /// </summary>
    public double ActivateNa {
        get { return SelfNeuron[activateIdx]; }
        set { SelfNeuron[activateIdx] = value; }
    }

    /// <summary>
    /// 不活性化変数の取得と設定．
    /// </summary>
    public double InactivateNa {
        get { return SelfNeuron[inactivateIdx]; }
        set { SelfNeuron[inactivateIdx] = value; }
    }

    /// <summary>
    /// 微分方程式：Na活性化変数m
    /// </summary>
    /// <returns></returns>
    public double ODEActivateNa() {
        double alpha, beta;
        alpha = (-0.1 * (SelfNeuron.Potential + 35)) / (-1 + Math.Exp(-(SelfNeuron.Potential + 35) / 10));
        beta = 4 * Math.Exp(-(SelfNeuron.Potential + 60) / 18);
        return SelfNeuron.Temperature * (alpha * (1 - ActivateNa) - beta * ActivateNa);
    }

    /// <summary>
    /// 微分方程式：Na不活性化変数h
    /// </summary>
    /// <returns></returns>
    public double ODEInactivateNa() {
        double alpha, beta;
        alpha = 0.07 * Math.Exp(-(SelfNeuron.Potential + 60) / 20);
        beta = 1 / (1 + Math.Exp(-(SelfNeuron.Potential + 30) / 10));
        return SelfNeuron.Temperature * (alpha * (1 - InactivateNa) - beta * InactivateNa);
    }

    /// <summary>
    /// ナトリウム電流
    /// </summary>
    /// <returns></returns>
    public override double GetCurrent() {
        if(Enable) {
            return -Conductance * ActivateNa * ActivateNa * ActivateNa * InactivateNa * (SelfNeuron.Potential - ReversalPotential);
        } else {
            return 0.0;
        }
    }
}

ニューロンを作る

チャネルができて準備ができたのでいよいよニューロンを作ります。Neuronクラスから継承します。まず膜電位を微分する必要があるのでODEPotentialを微分変数に追加します。ODEPotentialは式(1)がすでに定義されていて，各電流は追加したチャネルから取得されます。

つぎにさっき作ったチャネルのインスタンスを作成します。インスタンスを作るときに現在のニューロンをコンストラクタに設定することでチャネルの関連付けは完了します。

つぎに膜容量や温度定数，各チャネルの開口率などを初期化してAddCurrentで式(1)の電流を設定します。

ほかに外部からチャネルの状態を見るためにアクセサなどを設定してもいいかもしれません。

式(1)のIMの外部からの入力電流は，ここでは定義せずメソッドのみ追加しておきます。

public class HHNeuron : Neuron {
    private LeakCurrent leak;
    private NaCurrent na;
    private KCurrent k;

    /// <summary>
    /// Hodgkin Huxley ニューロンモデルを作成します．
    /// </summary>
    public HHNeuron() {
        //膜電位微分方程式を追加します
        AddODE(ODEPotential, -60.0, "Potential"); //potential
        
        //各チャンネルを作ります
        leak = new LeakCurrent(this);
        na = new NaCurrent(this);
        k = new KCurrent(this);

        //膜容量・温度定数の初期化
        //初期化しない場合は1.0が代入されます
        Capacitance = 1.0;
        Temperature = 1.0;

        //各チャンネルの開口率などの初期化
        na.ActivateNa = 0.0529342;
        na.InactivateNa = 0.596111;
        k.ActivateK = 0.317682;

        //膜電位微分方程式に各チャンネルの電流を追加します
        AddCurrent(na.GetCurrent, k.GetCurrent, leak.GetCurrent, GetStimulusCurrent);
    }

    #region 各チャンネルのアクセサ
    /// <summary>
    /// リークチャンネルを取得します．
    /// </summary>
    public LeakCurrent Leak { get { return leak; } }

    /// <summary>
    /// Naチャンネルを取得します．
    /// </summary>
    public NaCurrent NaChannel { get { return na; } }

    /// <summary>
    /// Kチャンネルを取得します．
    /// </summary>
    public KCurrent KChannel { get { return k; } }
    #endregion
}

シミュレーションする

さあ，これで準備が全て整いました。シミュレーション呼び出すメインを作っていきます。コードのキモは，Calculatorで微分計算機を初期化してそれに，新しいニューロンのインスタンスを追加します。そして，そこに刺激入力を定義し，メインループで回して行きます。ログの取得や発火率の計算もできるようになってるのでサンプルファイルを見ながらやってみてね！

static void Main(string[] args) {
    SimulationLog simlog = new SimulationLog("simlog.txt"); //計算時間などのログを取ります
    Calculator calc = new Calculator(ODEType.RungeKutta4);//微分方程式の計算方法を指定します
    SimulationTime.Delta = 0.005;//計算のきざみ幅を設定します
    SimulationTime.EndTime = 500;//終了時間を設定します
    SimulationTime.Interval = 50;//出力ファイルのインターバルを設定します
    ProgressBar pb = new ProgressBar(ProgressBar.Type.Percents);//プログレスバーを表示するようにします

    HHNeuron n = new HHNeuron();//あたらしいHHNeuronを作ります．

    n.StimulusCurrentFunction = delegate(double t) {//刺激入力を定義します
        if(10 < t && t <= 210) {
            return 7.0;
        } else if(210 < t && t <= 410) {
            return 18.0;
        } else {
            return 0.0;
        }
    };

    MeasureFiringRate fr = new MeasureFiringRate(n);//nの発火率を計測するfrを作ります．
    calc.AddNeuron(n);//微分計算機にnを追加します

    try {//計算途中に値がトンだ場合に備えて
        //出力ファイル名の設定をします
        using(DataWriterBinder logs = new DataWriterBinder("data0.txt", "data1.txt", "rate.txt")) {
            while(SimulationTime.IsSimulation()) {//メインループ
                //ログ書き込みのメソッドは下記のWriteとかWriteWithTime以外にもいろいろあるので
                //試してみてね
                logs[0].WriteWithTime(n.Potential, n.NaChannel.ActivateNa);//ログの書き込み
                logs[1].Write(n);//ログの書き込み
                logs[2].WriteAsyncTimeData(fr.GetFiringRate());//発火率の書き込み
                calc.Calc();//計算を1ステップ進めます
            }
        }
    } catch(InvalidOperationException ex) {//値が飛んだ場合は強制終了
        Console.Error.Write(ex.ToString());
    }
    pb.Stop();//プログレスバーを停止します．
    simlog.Finish(calc);//シミュレーションのログを取ります
}