metconan.c

GESPI 2.0动态系统模拟工具

#include "copyleft.h"

/*
    GEPASI - a simulator of metabolic pathways and other dynamical systems
    Copyright (C) 1989, 1992  Pedro Mendes
*/

/*************************************/
/*                                   */
/*     metabolic control analysis    */
/*        using Reder's method       */
/*                                   */
/*          MICROSOFT C 6.00         */
/*           QuickC/WIN 1.0          */
/*             ULTRIX cc             */
/*              GNU gcc              */
/*                                   */
/*   (include here compilers that    */
/*   compiled GEPASI successfully)   */
/*                                   */
/*************************************/


#include <stdio.h>
#include <math.h>
#include <float.h>
#include "globals.h"                       /* global symbols                 */
#include "globvar.h"                       /* extern global variable         */
#include "matrix.h"                        /* symbols & function declarations*/
#include "datab.h"

#define MCA_OK 0
#define MCA_NOCC 1

void calc_Dxv( void )
{
 register int i, j;

 for( i=0; i<nsteps; i++ )
  for( j=0; j<nmetab; j++ )
   Dxv[i][j] = partder[i][j]( xss, i, j );
}

int calc_Gamma( void )
{
 register int k, j;
 int i;

 for(i=0;i<indmet;i++)                               /* aux1 = rstoi * Dxv   */
  for(j=0;j<nmetab;j++)
  {
   aux1[i][j] = (double) 0;
   for(k=0;k<nsteps;k++) aux1[i][j] += (double) rstoi[i][k] * Dxv[k][j];
  }
 for(i=0;i<indmet;i++)                               /* aux2 = aux1 * ml     */
  for(j=0;j<indmet;j++)
  {
   aux2[i][j] = (double) 0;
   for(k=0;k<nmetab;k++) aux2[i][j] += aux1[i][k] * (double) ml[k][j];
  }
 if( lu_inverse( (double(*)[MAX_MET][MAX_MET])aux2,
                 (double(*)[MAX_MET][MAX_MET])aux1,
                 indmet
               ) == NON_INVERT )	                 /* aux1 = 1 / aux2      */
 {
  printf("  * non-invertible Dxv\n");
  return MCA_NOCC;
 }
 for(i=0;i<nmetab;i++)								 /* aux2 = ml * aux1     */
  for(j=0;j<indmet;j++)
  {
   aux2[i][j] = (double) 0;
   for(k=0;k<indmet;k++) aux2[i][j] += (double) ml[i][k] * aux1[k][j];
  }
 for(i=0;i<nmetab;i++)                               /* Gamma = -aux2*rstoi  */
  for(j=0;j<nsteps;j++)
  {
   Gamma[i][j] = (double) 0;
   for(k=0;k<indmet;k++) Gamma[i][j] -= aux2[i][k] * (double) rstoi[k][j];
  }
 return MCA_OK;
}

void calc_C( void )
{
 register int i, j;
 int k;

 for(i=0;i<nsteps;i++)                               /* aux = Dxv * Gamma    */
  for(j=0;j<nsteps;j++)
  {
   aux1[i][j] = (double) 0;
   for(k=0;k<nmetab;k++) aux1[i][j] += Dxv[i][k] * Gamma[k][j];
  }
 for(i=0;i<nsteps;i++)                               /* C = I + aux          */
  for(j=0;j<nsteps;j++)
   C[i][j] = aux1[i][j] + ( i==j ? (double) 1 : (double) 0 );
}

void calc_tt( void )
{
 register int i,q;

 ttt = (double) 0;                                   /* reset total t.time   */
 for( i=0; i<nmetab; i++)
 {
  tt[i] = (double) 0;                                /* reset trans. time    */
  for( q=0; q<nsteps; q++ )
   if ( stoi[i][q] > 0 )                             /* sum of positive fl.  */
    tt[i] += fabs( rateq[q](xss, q) );
  if ( tt[i] == (double) 0 )                         /* if no positive fl... */
  {
   for( q=0; q<nsteps; q++ )
    if ( stoi[i][q] < 0 )                            /* sum of negative fl.  */
     tt[i] += fabs( rateq[q](xss, q) );
  }
  if ( tt[i] != (double) 0 )
  {
   tt[i] = xss[i] / tt[i];
   ttt += tt[i];
  }
  else
  {
   if (xss[i] == x0[i]) tt[i] = (double) 0;      /* disguised extern metablt */
   if (xss[i] != x0[i]) ttt = tt[i] = INFINITY;  /* if rate=0 then tt=+oo    */
  }
 }
 for( i=nmetab; i<totmet; i++ ) tt[i] = (double) 0;
}

int equilibrium( void )
{
 register int i;
 double hr;

 for( i=0, hr = (double) 0; i<nsteps; i++ )
  if ( rateq[i](xss, i) > hr )
   hr = rateq[i](xss, i);

 return ( hr <= options.hrcz );
}

void scale_coef( void )
{
 int i, j;

 for( i=0; i<nsteps; i++ )
  for(j=0;j<totmet;j++)
  {
   if( fabs( flux[i] ) >= options.hrcz )
    Dxv[i][j] *= xss[j] / flux[i];
   else Dxv[i][j] = FLT_MAX;
  }

 for( i=0; i<totmet; i++)
  for( j=0; j<nsteps; j++)
  {
   if( fabs( xss[i] ) >= options.hrcz )
    Gamma[i][j] *= flux[j] / xss[i];
   else
    Gamma[i][j] = FLT_MAX;
  }

 for( i=0; i<nsteps; i++)
  for( j=0; j<nsteps; j++)
  {
   if(
      ( fabs( flux[i] ) >= options.hrcz ) &&
      ( fabs( flux[j] ) >= options.hrcz )
     )
    C[i][j] *= flux[j] / flux[i];
   else C[i][j] = FLT_MAX;
  }
}