test1.cpp

Go to the documentation of this file.

#include "indii/fmri/hemodynamic/FlowBalloonModel.hpp"
#include "indii/fmri/hemodynamic/BOLDCalculator.hpp"
#include "indii/ml/ode/AdaptiveRungeKutta.hpp"

#include <iostream>
#include <string.h>

using namespace std;
using namespace indii::fmri::hemodynamic;
using namespace indii::ml::ode;

/**
 * @file test1.cpp
 *
 * Basic test of FlowBalloonModel.
 *
 * This test configures a FlowBalloonModel with parameters and
 * functions set according to the first experiment in Buxton et
 * al. (1998). Results are output on stdout, tab delimited, with
 * columns representing:
 *
 * \li \f$t\f$; time
 * \li \f$f_{in}(t)\f$
 * \li \f$f_{out}(t)\f$
 * \li \f$CMRO_2(t) = \frac{E(t)}{E_0}f_{in}(t)\f$ (Buxton et al. 2004).
 * \li \f$q\f$
 * \li \f$v\f$
 * \li \f$y\f$; BOLD response (%)
 *
 * Results are as follows:
 *
 * \image html test1.png "Results, c.f. Buxton et al. (1998), Figure 1"
 * \image latex test1.eps "Results, c.f. Buxton et al. (1998), Figure 1"
 */

/**
 * \f$\tau_0\f$
 */
static const double TAU_0 = 2.0;

/**
 * \f$E_0\f$
 */
static const double E_0 = 0.4;

/**
 * \f$V_0\f$
 */
static const double V_0 = 0.01;

/**
 * \f$\alpha\f$
 */
static const double ALPHA = 0.5;

/**
 * Duration of simulation (s).
 */
static const double END = 20.0;

/**
 * \f$f_{in}(t)\f$; Trapezoidal function with a rise time of 4 s,
 * duration of 4 s and fall time of 4 s.
 */
double F_IN(double t, const double y[], void* o) {
  double f;
  if (t < 4) {
    f = (0.7/4.0)*t + 1;
  } else if (t < 8) {
    f = 1.7;
  } else if (t < 12) {
    f = (-0.7/4.0)*(t - 8) + 1.7;
  } else {
    f = 1.0;
  }
  return f;
}

/**
 * @brief \f$f_{out}(t) = \frac{7}{3}v - \frac{4}{3}\f$;
 *
 * i.e. linear relative to \f$v\f$ with a gradient of
 * \f$\frac{7}{3}\f$, \f$f_{out} = 1\f$ when \f$v = 1\f$.
 */
double F_OUT(double t, const double y[], void* o) {
  return (7.0/3.0)*y[FlowBalloonModel::V] - 4.0/3.0;
}

/**
 * Run tests.
 */
int main(int argc, const char* argv[]) {
  FlowBalloonModel balloonModel;
  BOLDCalculator boldCalculator(&balloonModel);
  AdaptiveRungeKutta ode(&balloonModel, balloonModel.suggestInitialState());

  /* set up balloon model */
  balloonModel.setFunction(FlowBalloonModel::F_IN, F_IN);
  balloonModel.setFunction(FlowBalloonModel::F_OUT, F_OUT);
  balloonModel.setParameter(FlowBalloonModel::TAU_0, TAU_0);
  balloonModel.setParameter(FlowBalloonModel::E_0, E_0);
  balloonModel.setParameter(FlowBalloonModel::V_0, V_0);
  balloonModel.setParameter(FlowBalloonModel::ALPHA, ALPHA);

  double cmro2;

  double t = 0.0;
  while (t < END) {
    cmro2 = balloonModel.getFunction(FlowBalloonModel::E) /
        balloonModel.getParameter(FlowBalloonModel::E_0) *
        balloonModel.getFunction(FlowBalloonModel::F_IN);
    cout << t << "\t";
    cout << balloonModel.getFunction(FlowBalloonModel::F_IN) << "\t";
    cout << balloonModel.getFunction(FlowBalloonModel::F_OUT) << "\t";
    cout << cmro2 << "\t";
    cout << ode.getVariable(FlowBalloonModel::Q) << "\t";
    cout << ode.getVariable(FlowBalloonModel::V) << "\t";
    cout << boldCalculator.calculate(ode.getState()) * 100.0 << "\t";
    cout << endl;

    t = ode.step(END);
  }

  return 0;
}

---