odeintでMAPK cascadeのシミュレーション

odeintを試す題材としてMAPK cascadeのモデルを利用してみる。

成果物はここ → GitHub - mythosil/odeint-mapk

モデルの入手

モデルはBioModelsから取得。
BIOMD0000000010 - Kholodenko2000 - Ultrasensitivity and negative feedback bring oscillations in MAPK cascade

結果確認に利用するため、Online Simulatorでシミュレーション結果を入手。
Actions → BioModels Online Simulationと選択していけば良い。
Select all speciesにチェックを入れ、シミュレーションを実行する。
結果をgnuplotで表示するとこのようになる。

シミュレータ実装

まず必要なものをinclude。

#include <iostream>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>

次にsystemを作る。BioModelsから数式を写し書きする作業。

/**
 * @see https://www.ebi.ac.uk/biomodels-main/BIOMD0000000010
 */
struct MapkSystem {

  using StateType = boost::array<double, 8>;

  static constexpr StateType InitialState = {
      90.0,  // Mos
      10.0,  // Mos-P
      280.0, // Mek1
      10.0,  // Mek1-P
      10.0,  // Mek1-PP
      280.0, // Erk2
      10.0,  // Erk2-P
      10.0   // Erk2-PP
  };

  /* 係数定義は略 */

  void operator()(const StateType &x, StateType &dxdt, double t) {
    // d[Mos]/dt =
    //     - uVol*V1*Mos/((1+(Erk2-PP/Ki)^n)*(K1+Mos))
    //     + uVol*V2*Mos-P/(KK2+Mos-P)
    dxdt[0] =
        - V1 * x[0] / ((1.0 + pow(x[7] / Ki, n)) * (K1 + x[0]))
        + V2 * x[1] / (KK2 + x[1]);

    // d[Mos-P]/dt =
    //     uVol*V1*Mos/((1+(Erk2-PP/Ki)^n)*(K1+Mos))
    //     - uVol*V2*Mos-P/(KK2+Mos-P)
    dxdt[1] =
        V1 * x[0] / ((1.0 + pow(x[7] / Ki, n)) * (K1 + x[0]))
        - V2 * x[1] / (KK2 + x[1]);

    // d[Mek1]/dt =
    //     - uVol*k3*Mos-P*Mek1/(KK3+Mek1)
    //     + uVol*V6*Mek1-P/(KK6+Mek1-P)
    dxdt[2] =
        - k3 * x[1] * x[2] / (KK3 + x[2])
        + V6 * x[3] / (KK6 + x[3]);

    // d[Mek1-P]/dt =
    //     uVol*k3*Mos-P*Mek1/(KK3+Mek1)
    //     - uVol*k4*Mos-P*Mek1-P/(KK4+Mek1-P)
    //     + uVol*V5*Mek1-PP/(KK5+Mek1-PP)
    //     - uVol*V6*Mek1-P/(KK6+Mek1-P)
    dxdt[3] =
        k3 * x[1] * x[2] / (KK3 + x[2])
        - k4 * x[1] * x[3] / (KK4 + x[3])
        + V5 * x[4] / (KK5 + x[4])
        - V6 * x[3] / (KK6 + x[3]);

    // d[Mek1-PP]/dt =
    //     uVol*k4*Mos-P*Mek1-P/(KK4+Mek1-P)
    //     - uVol*V5*Mek1-PP/(KK5+Mek1-PP)
    dxdt[4] =
        k4 * x[1] * x[3] / (KK4 + x[3])
        - V5 * x[4] / (KK5 + x[4]);

    // d[Erk2]/dt =
    //     - uVol*k7*Mek1-PP*Erk2/(KK7+Erk2)
    //     + uVol*V10*Erk2-P/(KK10+Erk2-P)
    dxdt[5] =
        - k7 * x[4] * x[5] / (KK7 + x[5])
        + V10 * x[6] / (KK10 + x[6]);

    // d[Erk2-P]/dt =
    //     uVol*k7*Mek1-PP*Erk2/(KK7+Erk2)
    //     - uVol*k8*Mek1-PP*Erk2-P/(KK8+Erk2-P)
    //     + uVol*V9*Erk2-PP/(KK9+Erk2-PP)
    //     - uVol*V10*Erk2-P/(KK10+Erk2-P)
    dxdt[6] =
        k7 * x[4] * x[5] / (KK7 + x[5])
        - k8 * x[4] * x[6] / (KK8 + x[6])
        + V9 * x[7] / (KK9 + x[7])
        - V10 * x[6] / (KK10 + x[6]);

    // d[Erk2-PP]/dt =
    //     uVol*k8*Mek1-PP*Erk2-P/(KK8+Erk2-P)
    //     - uVol*V9*Erk2-PP/(KK9+Erk2-PP)
    dxdt[7] =
        k8 * x[4] * x[6] / (KK8 + x[6])
        - V9 * x[7] / (KK9 + x[7]);
  }

};

途中経過を表示するためobserverを用意する。標準出力にCSV形式で出力させる。

struct StdoutObserver {

  void operator()(const MapkSystem::StateType &x , const double t) {
    std::cout << t;

    auto size = x.size();
    for (auto i = 0; i < size; i++) {
      std::cout << " " << x[i];
    }

    std::cout << std::endl;
  }

};

RK4でシミュレーション実行。

int main(int argc, char **argv) {
  const double start = 0.0;
  const double duration = 1000.0;
  const double interval = 0.1;

  MapkSystem system;
  StdoutObserver observer;
  odeint::runge_kutta4<MapkSystem::StateType> stepper;

  auto state = MapkSystem::InitialState;

  odeint::integrate_const(stepper, system, state, start, duration, interval, std::ref(observer));

  return 0;
}

RK-Dopriならこう書く。ステップサイズをadaptiveにすることで、細かく刻むよりも実行時間を大幅に短縮できる。

int main(int argc, char **argv) {
  const double start = 0.0;
  const double duration = 1000.0;
  const double interval = 0.1;

  MapkSystem system;
  StdoutObserver observer;
  odeint::runge_kutta_dopri5<MapkSystem::StateType> stepper;
  auto controlledStepper = odeint::make_controlled(0.0001, 0.0001, stepper);

  auto state = MapkSystem::InitialState;

  odeint::integrate_adaptive(controlledStepper, system, state, start, duration, interval, std::ref(observer));

  return 0;
}

結果

Online Simulatorとほぼ同じ結果が得られている。

Next Step

陰解法も比較的容易に実装できるので試す。ヤコビ行列を用意する必要あり。
また、途中でイベントが挟まるようなモデルをodeintでどう扱うかも調べる。公式にはサポートされていないようなので、odeint自体に手を入れる必要がありそう。