kinetics.c

GESPI 2.0动态系统模拟工具

#include "copyleft.h"

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

/*************************************/
/*                                   */
/*      reaction rate equations      */
/*                                   */
/*        Zortech C/C++ 3.0 r4       */
/*          MICROSOFT C 6.00         */
/*          Visual C/C++ 1.0         */
/*           QuickC/WIN 1.0          */
/*             ULTRIX cc             */
/*              GNU gcc              */
/*                                   */
/*   (include here compilers that    */
/*   compiled GEPASI successfully)   */
/*                                   */
/*************************************/

#include <math.h>
#include "globals.h"
#include "globvar.h"
#include "datab.h"

/*****************************************/
/*                                       */
/*              CONVENTIONS              */
/*                                       */
/* all functions take as arguments:      */
/* double s[] - vector of concentrations */
/* int r - the index of the reaction     */
/*                                       */
/*     all fuctions return a double.     */
/*                                       */
/*                                       */
/*****************************************/


/*                                   */
/* user defined rate-equations		 */
/*                                   */

double value_of( double s[], int r, struct treet *t, int idx )
{
 switch( t->node[idx].item )
 {
  case 'N': return (double) t->constant[(int)t->node[idx].val];
  case 'I': switch( (int) t->id[(int)t->node[idx].val][9] )
            {
             case 0: return params[r][ (int) t->id[(int)t->node[idx].val][0] ];
             case 1:
             case 2:
             case 3: return s[ eff[r][ (int) t->id[(int)t->node[idx].val][0] ] ];
            }
  case 'F': switch( t->node[idx].val )
            {
             case '+': return value_of( s, r, t, (int) t->node[idx].left );
             case '-': return - value_of( s, r, t, (int) t->node[idx].left );
             case 'L': return log10( value_of( s, r, t, (int) t->node[idx].left ) );
             case 'l': return log( value_of( s, r, t, (int) t->node[idx].left ) );
             case 'e': return exp( value_of( s, r, t, (int) t->node[idx].left ) );
             case 'S': return sin( value_of( s, r, t, (int) t->node[idx].left ) );
             case 'C': return cos( value_of( s, r, t, (int) t->node[idx].left ) );
            }
  case 'O': switch( t->node[idx].val )
            {
             case '+': return value_of( s, r, t, (int) t->node[idx].left ) +
                              value_of( s, r, t, (int) t->node[idx].right );
             case '-': return value_of( s, r, t, (int) t->node[idx].left ) -
                              value_of( s, r, t, (int) t->node[idx].right );
             case '*': return value_of( s, r, t, (int) t->node[idx].left ) *
                              value_of( s, r, t, (int) t->node[idx].right );
             case '/': return value_of( s, r, t, (int) t->node[idx].left ) /
                              value_of( s, r, t, (int) t->node[idx].right );
             case '^': return pow( value_of( s, r, t, (int) t->node[idx].left ),
                                   value_of( s, r, t, (int) t->node[idx].right ) );
            }
 }
}

double translate(double s[], int r)
{
 return value_of( s, r, &tree[kinetype[r]-MAX_TYP], (int) tree[kinetype[r]-MAX_TYP].node[0].left );
}


/*                                   */
/* the zero returning function       */
/*                                   */

double ret_zero( double s[], int r, int e )
{
 return (double) 0;
}


/*                                   */
/* derivatives for effectors         */
/* by finite differences             */
/*                                   */

double deff( double s[], int r, int e )
{
 double k, x, x1, f1, f2;

 x = k = s[e];
 /*
    if x is zero, we will calculate the derivative at a small positive
    value (no point in considering negative values!). let's stick with
    hrcz (the flux resolution)
 */
 if( x==0 ) x = options.hrcz;
 x1 = x * 0.01;
 s[e] = x - x1;
 f1 = rateq[r]( s, r );
 s[e] = x + x1;
 f2 = (rateq[r]( s, r ) - f1)/(x*0.02);
 s[e] = k;
 return f2;
}

/*                                   */
/* A -->  or  --> B                  */
/*                                   */

double i0(double s[], int r)
{
 return params[r][0];
}


/*                                   */
/* A --> B                           */
/*                                   */

double i11(double s[], int r)
{
 return s[eff[r][0]]*params[r][0];
}

/*                                 */
/* A <--> B                        */
/*                                 */

double r11(double s[], int r)
{
 return s[eff[r][0]]*params[r][0]-s[eff[r][1]]*params[r][1];
}

/*                                 */
/* A + B --> C                     */
/*                                 */

double i21(double s[], int r)
{
 return s[eff[r][0]]*s[eff[r][1]]*params[r][0];
}

/*                                 */
/* A + B <--> C                    */
/*                                 */

double r21(double s[], int r)
{
 return s[eff[r][0]]*s[eff[r][1]]*params[r][0]-s[eff[r][2]]*params[r][1];
}

/*                                   */
/* A --> B + C                       */
/*                                   */

double i12(double s[], int r)
{
 return s[eff[r][0]]*params[r][0];
}

/*                                 */
/* A <--> B + C                    */
/*                                 */

double r12(double s[], int r)
{
 return s[eff[r][0]]*params[r][0]-s[eff[r][1]]*s[eff[r][2]]*params[r][1];
}

/*                                 */
/* A + B + C --> D                 */
/*                                 */

double i31(double s[], int r)
{
 return s[eff[r][0]]*s[eff[r][1]]*s[eff[r][2]]*params[r][0];
}

/*                                 */
/* A + B + C <--> D                */
/*                                 */

double r31(double s[], int r)
{
 return s[eff[r][0]]*s[eff[r][1]]*s[eff[r][2]]*params[r][0]-s[eff[r][3]]*params[r][1];
}

/*                                   */
/* A --> B + C + D                   */
/*                                   */

double i13(double s[], int r)
{
 return s[eff[r][0]]*params[r][0];
}

/*                                 */
/* A <--> B + C + D                */
/*                                 */

double r13(double s[], int r)
{
 return s[eff[r][0]]*params[r][0]-s[eff[r][1]]*s[eff[r][2]]*s[eff[r][3]]*params[r][1];
}

/*                                 */
/* A + B --> C + D                 */
/*                                 */

double i22(double s[], int r)
{
 return s[eff[r][0]]*s[eff[r][1]]*params[r][0];
}

/*                                 */
/* A + B <--> C + D                */
/*                                 */

double r22(double s[], int r)
{
 return s[eff[r][0]]*s[eff[r][1]]*params[r][0]-s[eff[r][2]]*s[eff[r][3]]*params[r][1];
}

/*                                 */
/* A + B + C --> D + E             */
/*                                 */

double i32(double s[], int r)
{
 return s[eff[r][0]]*s[eff[r][1]]*s[eff[r][2]]*params[r][0];
}

/*                                 */
/* A + B + C <--> D + E            */
/*                                 */

double r32(double s[], int r)
{
 return s[eff[r][0]]*s[eff[r][1]]*s[eff[r][2]]*params[r][0]-s[eff[r][3]]*s[eff[r][4]]*params[r][1];
}

/*                                 */
/* A + B --> C + D + E             */
/*                                 */

double i23(double s[], int r)
{
 return s[eff[r][0]]*s[eff[r][1]]*params[r][0];
}

/*                                 */
/* A + B <--> C + D + E            */
/*                                 */

double r23(double s[], int r)
{
 return s[eff[r][0]]*s[eff[r][1]]*params[r][0]-s[eff[r][2]]*s[eff[r][3]]*s[eff[r][4]]*params[r][1];
}

/*                                 */
/* A + B + C --> D + E + F         */
/*                                 */

double i33(double s[], int r)
{
 return s[eff[r][0]]*s[eff[r][1]]*s[eff[r][2]]*params[r][0];
}

/*                                 */
/* A + B + C <--> D + E + F        */
/*                                 */

double r33(double s[], int r)
{
 return s[eff[r][0]]*s[eff[r][1]]*s[eff[r][2]]*params[r][0]-s[eff[r][3]]*s[eff[r][4]]*s[eff[r][5]]*params[r][1];
}


/*                                 */
/* derivatives for one substrate   */
/*                                 */

double d1Xs1( double s[], int r, int e )
{
 return params[r][0];
}


/*                                 */
/* derivatives for two substrates  */
/*                                 */

double d2Xs1( double s[], int r, int e )
{
 return params[r][0]*s[eff[r][1]];
}

double d2Xs2( double s[], int r, int e )
{
 return params[r][0]*s[eff[r][0]];
}


/*                                 */
/* derivatives for three substrates*/
/*                                 */

double d3Xs1( double s[], int r, int e )
{
 return params[r][0]*s[eff[r][1]]*s[eff[r][2]];
}

double d3Xs2( double s[], int r, int e )
{
 return params[r][0]*s[eff[r][0]]*s[eff[r][2]];
}

double d3Xs3( double s[], int r, int e )
{
 return params[r][0]*s[eff[r][0]]*s[eff[r][1]];
}


/*                                 */
/* derivatives for one product     */
/*                                 */

double dX1p1( double s[], int r, int e )
{
 return -params[r][1];
}


/*                                 */
/* derivatives for two products    */
/* (reactions with one substrate)  */
/*                                 */

double dX2p1s1( double s[], int r, int e )
{
 return -params[r][1]*s[eff[r][2]];
}

double dX2p2s1( double s[], int r, int e )
{
 return -params[r][1]*s[eff[r][1]];
}

/*                                 */
/* derivatives for two products    */
/* (reactions with two substrates) */
/*                                 */

double dX2p1s2( double s[], int r, int e )
{
 return -params[r][1]*s[eff[r][3]];
}

double dX2p2s2( double s[], int r, int e )
{
 return -params[r][1]*s[eff[r][2]];
}

/*                                 */
/* derivatives for two products    */
/* (reactions with 3 substrates)   */
/*                                 */

double dX2p1s3( double s[], int r, int e )
{
 return -params[r][1]*s[eff[r][4]];
}

double dX2p2s3( double s[], int r, int e )
{
 return -params[r][1]*s[eff[r][3]];
}


/*                                 */
/* derivatives for three products  */
/* (reactions with 1 substrate)   */
/*                                 */

double dX3p1s1( double s[], int r, int e )
{
 return -params[r][1]*s[eff[r][2]]*s[eff[r][3]];
}

double dX3p2s1( double s[], int r, int e )
{
 return -params[r][1]*s[eff[r][1]]*s[eff[r][3]];
}

double dX3p3s1( double s[], int r, int e )
{
 return  -params[r][1]*s[eff[r][1]]*s[eff[r][2]];
}

/*                                 */
/* derivatives for three products  */
/* (reactions with 2 substrates)   */
/*                                 */

double dX3p1s2( double s[], int r, int e )
{
 return -params[r][1]*s[eff[r][3]]*s[eff[r][4]];
}

double dX3p2s2( double s[], int r, int e )
{
 return -params[r][1]*s[eff[r][2]]*s[eff[r][4]];
}

double dX3p3s2( double s[], int r, int e )
{
 return  -params[r][1]*s[eff[r][2]]*s[eff[r][3]];
}

/*                                 */
/* derivatives for three products  */
/* (reactions with 3 substrates)   */
/*                                 */

double dX3p1s3( double s[], int r, int e )
{
 return -params[r][1]*s[eff[r][4]]*s[eff[r][5]];
}

double dX3p2s3( double s[], int r, int e )
{
 return -params[r][1]*s[eff[r][3]]*s[eff[r][5]];
}

double dX3p3s3( double s[], int r, int e )
{
 return  -params[r][1]*s[eff[r][3]]*s[eff[r][4]];
}



/*                                                   */