test2.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>
#include <math.h>

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

/**
 * @file test2.cpp
 *
 * Basic test of FlowBalloonModel.
 *
 * This test configures a FlowBalloonModel with parameters and
 * functions set according to the second 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 test2.png "Results, c.f. Buxton et al. (1998), Figure 2"
 * \image latex test2.eps "Results, c.f. Buxton et al. (1998), Figure 2"
 *
 * @note The only difference between this test and that in test1.cpp
 * is the definition of F_OUT -- linear in the case of test1.cpp and
 * nonlinear in the case of test1.cpp.
 *
 * @note Look out for a double line in the graph for \f$f_{out}(v)\f$
 * which may indicate that the error bounds need to be tightened or
 * that discontinuities in \f$f_{in}(t)\f$ are not being handled
 * well. As the volume increases, peaks and then decreases with time
 * in this test, this graph has two lines drawn. They should be
 * exactly the same.
 */

/**
 * \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}{80}\left[\frac{20}{3}(v -
 * 1)\right]^{3} + 1\f$;
 *
 * i.e. nonlinear relative to \f$v\f$. The exact form of the function
 * used in Buxton et al. (1998) Figure 2 is unknown, but this seems a
 * reasonable approximation.
 */
double F_OUT(double t, const double y[], void* o) {
  return 0.7/8.0 * pow((y[FlowBalloonModel::V] - 1.0)*2.0/0.3, 3.0) + 1.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;
}

---