AUTOR: Vitaly Selivanov (Universidad de Barcelona)

Universalmente aceptan la modelización matemática como una herramienta para analizar los datos biológicos. El uso de los algoritmos de NAG simplifica mucho la modelización y tratamiento de los datos. Aunque existen programas informáticos especializados y lenguajes de programación de nivel alto como "Mathematica" (Wolfram Research) o "Matlab" que contienen algoritmos interiores que facilitan formulación matemática y resolución de problemas biológicas, sin embargo, en muchos casos las rutinas de NAG son preferibles. A menudo los programas especializados no cumplen las funciones necesarias en algunos casos particulares, y los lenguajes de programación de nivel alto son muy lentos en solución de los problemas realmente grandes; en este caso los algoritmos rápidos y eficaces de las bibliotecas de NAG Fortran o NAG C++ están por encima de cualquier competición.

Como ejemplo de tales problemas sirven modelización de:

• reacción-difusión de metabolitos intracelulares, que requiere resolución de ecuaciones diferenciales parabólicas con derivadas parciales en espacio de dos dimensiones (como se muestre en [1]);
• la dinámica de distribución de los isótopos estables, usados como marcas, para valoración de los flujos metabólicos intracelulares, que se requiere resolución los sistemas de cientos ecuaciones diferenciales ordinarias (ODE) [2];
• el funcionamiento de complejos respiratorios mitocondriales que se requiere resolución los sistemas grandes de ODE [3].

La aplicación de rutinas "D03RAF"de la biblioteca NAG Fortran resuelve problemas de difusión muy eficiente; para resolución numérica de sistemas grandes de ODE la aplicación de NAG del método DASSL en "D02NEF" es muy eficaz.

Para resolución numérica de sistemas de ecuaciones diferenciales parciales o sistemas grandes de ecuaciones diferenciales ordinarias los algoritmos NAG son vitales, porque ellos funcionan cientos veces más rápido que los programas escritos en lenguajes de nivel alto. Sin embargo, el software que usa estos algoritmos a menudo requiere más rutinas auxiliares y es muy cómodo cogerlos de la misma biblioteca de NAG Fortran. Por ejemplo el problema de optimización de parámetros de sistemas grande de ODE requiere ordenamiento ascensional de masivos de números reales, que son los mínimos locales de desviaciones entre los cálculos y datos experimentales. Los algoritmos de NAG, como "M01CAF", sirven perfectamente para ordenamiento de los vectores reales.

También para tratamiento de los datos sería útil comparar una distribución de los flujos metabólicos, que minimiza la diferencia entre cálculos y datos experimentales, con alguna distribución hipotética que optimiza algunos funciones, tales como producción de ATP o biosíntesis cualquier otro metabolito. Para encontrar la solución óptima de tal función objetivo las rutinas para programación lineal como "E04MFA" son muy útiles.

Los programas informáticos para tratamiento de datos citados son software multifunción escritos en C++, sin embargo esto no limitaba la implementación de rutinas NAG de Fortran mezclándoles con C++. Para usar C++ y Fortran conjunto solamente necesito un tratamiento simple que está explicado abajo.

Ejemplo de aplicación de rutinas NAG Fortran en un programa C++ usando compilador GNU g++

Véase a continuación un ejemplo de un programa escrito en C++ que llama los integradores D02NEF y D02MWF (que inicia los parámetros para D02NEF) de la biblioteca NAG Fortran:

void d02nefsv(const double tout, double *py,int iter) {
// definition of parameters for D02NEF and D02MWF:
int i,itask,lcom=200000,licom=500,icom[500],iuser[2],ires=0, maxord=5, itol=0, ifail=0;
double t=0.,xi=0., yprime[NN], rtol[NN], atol[NN], com[200000], ruser[2], h0=0.1e-6, hmax=22.0, to=tout;
//NN=400 is defined elsewhere, number of variables
rtol[0]=0.0005; atol[0]=3.0e-7;
char *jceval="Numeric";
//reading the initial values for 400 variables (*py) from a file:
horse.read(py);
// setup for D02NEF:
d02mwf_(NN,maxord,jceval[0],hmax,h0,itol,icom,licom,com,lcom,ifail);
//setting the time derivatives (yprime) of variables (py):
horse.distr(py, yprime);
//call integrator D02NEF:
d02nef_(NN,t,to,py,yprime,rtol,atol,itask,resd02,jacd02,icom,com,lcom,iuser, ruser,ifail );
}

Al integrador NAG, D02NEF, se puede llamarlo como cualquier otra función C++, excepto que es necesita el carácter '_'. Además es importante notar que numeración de elementos de un vector en Fortran empieza de 1 y en C++ de 0, entonces a[i] en Fortran es lo mismo que a[i-1] en C++.

Las funciones de Fortran se puede declarar en C++ como extern "C", las Funciones de Fortran D02NEF y D02MWF se declaran en C++ como:

extern "C" void d02mwf_(int& NEQ,int& MAXORD,char& JCEVAL,double& HMAX,double& H0,int& ITOL,int ICOM[],int& LICOM,double COM[],int& LCOM,int& IFAIL);

extern "C" void d02nef_(int& NEQ,double& T,double& TOUT,double Y[],double YDOT[],double *RTOL,double *ATOL,int& ITASK,
void (*RES)(int& NEQ,double& T,double Y[ ],double YDOT[],double R[],int& IRES,int *IUSER,double *RUSER),
void (*JAC)(int NEQ,double T,double Y[ ],double YDOT[],double **PD,double& CJ,int *IUSER,double *RUSER),
int ICOM[],double COM[],int& LCOM,int *IUSER,double *RUSER,int& IFAIL);

Cabe señalar que en Fortran se definen como:

SUBROUTINE D02MWF ( NEQ, MAXORD, JCEVAL, HMAX, H0, ITOL, ICOM, LICOM, COM, LCOM, IFAIL)
INTEGER NEQ, MAXORD, ITOL, ICOM(LICOM), LICOM, LCOM, IFAIL
double precision HMAX, H0, COM(LCOM)
CHARACTER*1 JCEVAL
SUBROUTINE D02NEF ( NEQ, T, TOUT, Y, YDOT, RTOL, ATOL, ITASK, RES, JAC, ICOM, COM, LCOM, IUSER, RUSER, IFAIL)
INTEGER NEQ, ITASK, ICOM(50+NEQ), LCOM, IUSER(*), IFAIL
double precision T, TOUT, Y(NEQ), YDOT(NEQ), RTOL(*), ATOL(*), COM(LCOM), RUSER(*)
EXTERNAL RES, JAC

La declaración y la definición correcta de parámetros como está realizado en el driver d02nefsv, que llama a los rutinas de Fortran, es suficiente para trabajar con ella en C++. Por supuesto, si esta usando Fortran95 es necesario compilar el programa completo con la biblioteca libgfortran (usando g++ option -lgfortran) y usando las bibliotecas distribuidas por NAG.

Hay que añadir algunas palabras sobre las rutinas exteriores RES y JAC:

• Jacobian puede ser determinada como función vacía, D02NEF puede calcularlo numéricamente:
  void jacd02(int NEQ,double t,double y[],double YDOT[],double **PD,double& cj,int *iuser,double *ruser){}
• La función RES que calcula los restos, puede llamar a otra función que automáticamente calcula las derivadas de variables por tiempo, que esta definida en un otro parte del programa:
  void resd02(int& neq,double& t,double y[],double ydot[],double r[],int& ires,int *iuser,double *ruser){
  double f[NN];
  horse.distr(y, f);
  for(unsigned i=0;i }

La biblioteca de NAG Fortran usada de este modo sirve no solo para los programas escritos en Fortran sino también para los programas en C++.

Referencias

[1] Selivanov VA, Krause S, Roca J, Cascante M. Modeling of spatial metabolite distributions in the cardiac sarcomere. Biophys J. 2007, 92, 3492-3500.
[2] Selivanov VA, Marin S, Lee PW, Cascante M. Software for dynamic analysis of tracer-based metabolomic data: estimation of metabolic fluxes and their statistical analysis. Bioinformatics. 2006, 22, 2806-2812.
[3] Selivanov VA, Votyakova TV, Zeak JA, Trucco M, Roca J, Cascante M. Bistability of mitochondrial respiration underlies paradoxical reactive oxygen species generation induced by anoxia. PLoS Comput Biol. 2009, 5(12):e1000619.