indii/fmri/hemodynamic/FlowBalloonModel.hpp

#ifndef INDII_FMRI_HEMODYNAMIC_FLOWBALLOONMODEL_HPP
#define INDII_FMRI_HEMODYNAMIC_FLOWBALLOONMODEL_HPP

#include "BalloonModel.hpp"

#include "indii/ml/aux/vector.hpp"

namespace indii {
  namespace fmri {
    namespace hemodynamic {

/**
 * Balloon %model driven by blood flow.
 *
 * @author Lawrence Murray <lawrence@indii.org>
 * @version $Rev: 285 $
 * @date $Date: 2007-07-20 17:25:40 +0100 (Fri, 20 Jul 2007) $
 *
 * This class implements the Balloon Model described by @ref
 * Buxton1998 "Buxton et al. (1998)". This system models a venous
 * compartment as a balloon, described by a volume and deoxyhemoglobin
 * (dHb) content. These combine to produce the Blood Oxygenation Level
 * Dependent (BOLD) signal.
 *
 * The %model is driven by blood flow into the venous
 * compartment. This causes the compartment to expand like a balloon,
 * increasing the pressure exerted by it and consequently the blood
 * flow out of the compartment. Concurrently, oxygen extraction, and
 * therefore the rate of dHb production, is more efficient at slower
 * flow rates. These simple interactions simulate a biologically
 * plausible BOLD signal.
 *
 * @section FlowBalloonModel_details Details
 *
 * The state variables \f$q\f$, representing normalised dHb content,
 * and \f$v\f$, representing normalised venous volume, are updated
 * according to the system of differential equations given by @ref
 * Buxton1998 "Buxton et al. (1998)":
 *
 * @anchor diffsystem
 * \f{eqnarray*}
 *   \frac{dq}{dt} &=& \frac{1}{\tau_0}\left[
 *     f_{in}(t)\frac{E(t)}{E_0} -
 *     f_{out}(v)\frac{q}{v}
 *     \right] \\
 *   \frac{dv}{dt} &=& \frac{1}{\tau_0}\left[
 *     f_{in}(t) - f_{out}(v)
 *     \right]
 * \f}
 *
 * Other symbols represent biophysical parameters and functions of the
 * system, handled using indii::ml::ode::ParameterCollection and
 * indii::ml::ode::FunctionCollection, with default values given in
 * FlowBalloonModelDefaults.
 *
 * Singularities are handled using limits. See
 * FlowBalloonModelDefaults::DQDT.
 *
 * @section FlowBalloonModel_usage Usage
 *
 * The basic usage idiom is as follows. Firstly create an instance of
 * FlowBalloonModel and set up the initial state:
 *
 * @code
 * FlowBalloonModel bm;
 * indii::ml::aux::vector y(bm.suggestInitialState());
 * @endcode
 *
 * The %model may then be simulated by solving the differential system
 * using indii::ml::ode::AdaptiveRungeKutta:
 *
 * @code
 * indii::ml::ode::AdaptiveRungeKutta solver(&bm, y);
 * double end = 60.0;
 * double t = 0.0;
 *
 * while (t < end) {
 *   t = solver.step(end);
 * }
 * @endcode
 *
 * In order for the %model to do something interesting, however, you
 * will want to at least override the default FlowBalloonModel::F_IN
 * function with your own custom function that controls flow. The
 * following function, for example, generates a flow of 2.0 for 1 s,
 * followed by a constant flow of 1.0:
 *
 * @code
 * double F_IN(double t, const double y[], void* o) {
 *   if (t < 1.0) {
 *     return 2.0;
 *   } else {
 *     return 1.0;
 *   }
 * }
 * @endcode
 *
 * Note the signature of this function, described by
 * indii::ml::ode::f_t. The %model can be made to use this function by
 * calling setFunction() before the first call to step():
 *
 * @code
 * bm.setFunction(FlowBalloonModel::F_IN, F_IN);
 * @endcode
 *
 * Other biophysical parameters and functions may be customised using
 * the setParameter() and setFunction() methods inherited from
 * indii::ml::ode::ParameterCollection and
 * indii::ml::ode::FunctionCollection.
 *
 * Finally, you will want to actually retrieve the state of the %model
 * after each call to step() using getVariable() or getState(). For
 * example, modifying the loop above:
 *
 * @code
 * while (t < end) {
 *   t = solver.step(end);
 *
 *   std::cout << t << ',';
 *   std::cout << solver.getVariable(FlowBalloonModel::Q) << ',';
 *   std::cout << solver.getVariable(FlowBalloonModel::V) << std::endl;
 * }
 * @endcode
 *
 * or:
 *
 * @code
 * while (t < end) {
 *   t = solver.step(end);
 *   y = solver.getState();
 *
 *   std::cout << t << ',';
 *   std::cout << y(FlowBalloonModel::Q) << ',';
 *   std::cout << y(FlowBalloonModel::V) << std::endl;
 * }
 * @endcode
 *
 * The values of biophysical functions, which are time dependent, may
 * be retrieved using getFunction() if they are of interest also.
 */
class FlowBalloonModel : public BalloonModel {
public:
  /**
   * State variable indices.
   */
  enum StateVariable {
    /**
     * \f$q = \frac{Q(t)}{Q_0}\f$; normalised dHb content of venous
     * compartment at the current time.
     */
    Q,

    /**
     * \f$v = \frac{V(t)}{V_0}\f$; normalised volume of the venous
     * compartment at the current time.
     */
    V
  };

  /**
   * Parameter indices.
   */
  enum BiophysicalParameter {
    /**
     * \f$V_0\f$; volume of the venous compartment at rest.
     */
    V_0,
  
    /**
     * \f$E_0\f$; oxygen extraction rate of the venous compartment at
     * rest.
     */
    E_0,

    /**
     * \f$\tau_0 = \frac{V_0}{F_0}\f$; mean transit time through the
     * venous compartment at rest.
     */
    TAU_0,

    /**
     * \f$\alpha\f$; Grubb's exponent.
     */
    ALPHA
  };

  /**
   * Function indices.
   */
  enum BiophysicalFunction {
    /**
     * \f$f_{in}(t) = \frac{F_{in}(t)}{F_0}\f$; normalised flow into
     * the venous compartment.
     */
    F_IN,

    /**
     * \f$f_{out}(t) = \frac{F_{out}(t)}{F_0}\f$; normalised flow out
     * of the venous compartment.
     */
    F_OUT,

    /**
     * \f$E(t)\f$; oxygen extraction rate of the venous compartment.
     */
    E
  };

  /**
   * Create new %model with the default biophysical parameters and
   * functions defined in FlowBalloonModelDefaults.
   */
  FlowBalloonModel();
  
  /**
   * Destructor.
   */
  virtual ~FlowBalloonModel();

  virtual indii::ml::aux::vector suggestInitialState();

  /**
   * @see indii::ml::ode::DifferentialModel
   */
  virtual void calculateDerivatives(double t, const double y[], double dydt[]);

private:
  /**
   * Number of state variables.
   */
  static const unsigned int DIMENSIONS = 2;

  /**
   * Number of parameters.
   */
  static const unsigned int NUM_PARAMETERS = 4;

  /**
   * Number of functions.
   */
  static const unsigned int NUM_FUNCTIONS = 3;

};
  
    }
  }
}

#endif

---