home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / biology / gsrc208a.zip / KINETICS.C < prev    next >
C/C++ Source or Header  |  1993-07-18  |  28KB  |  924 lines
  1. #include "copyleft.h"
  2.  
  3. /*
  4.     GEPASI - a simulator of metabolic pathways and other dynamical systems
  5.     Copyright (C) 1989, 1992  Pedro Mendes
  6. */
  7.  
  8. /*************************************/
  9. /*                                   */
  10. /*      reaction rate equations      */
  11. /*                                   */
  12. /*        Zortech C/C++ 3.0 r4       */
  13. /*          MICROSOFT C 6.00         */
  14. /*          Visual C/C++ 1.0         */
  15. /*           QuickC/WIN 1.0          */
  16. /*             ULTRIX cc             */
  17. /*              GNU gcc              */
  18. /*                                   */
  19. /*   (include here compilers that    */
  20. /*   compiled GEPASI successfully)   */
  21. /*                                   */
  22. /*************************************/
  23.  
  24. #include <math.h>
  25. #include "globals.h"
  26. #include "globvar.h"
  27. #include "datab.h"
  28.  
  29. /*****************************************/
  30. /*                                       */
  31. /*              CONVENTIONS              */
  32. /*                                       */
  33. /* all functions take as arguments:      */
  34. /* double s[] - vector of concentrations */
  35. /* int r - the index of the reaction     */
  36. /*                                       */
  37. /*     all fuctions return a double.     */
  38. /*                                       */
  39. /*                                       */
  40. /*****************************************/
  41.  
  42.  
  43. /*                                   */
  44. /* user defined rate-equations         */
  45. /*                                   */
  46.  
  47. double value_of( double s[], int r, struct treet *t, int idx )
  48. {
  49.  switch( t->node[idx].item )
  50.  {
  51.   case 'N': return (double) t->constant[(int)t->node[idx].val];
  52.   case 'I': switch( (int) t->id[(int)t->node[idx].val][9] )
  53.             {
  54.              case 0: return params[r][ (int) t->id[(int)t->node[idx].val][0] ];
  55.              case 1:
  56.              case 2:
  57.              case 3: return s[ eff[r][ (int) t->id[(int)t->node[idx].val][0] ] ];
  58.             }
  59.   case 'F': switch( t->node[idx].val )
  60.             {
  61.              case '+': return value_of( s, r, t, (int) t->node[idx].left );
  62.              case '-': return - value_of( s, r, t, (int) t->node[idx].left );
  63.              case 'L': return log10( value_of( s, r, t, (int) t->node[idx].left ) );
  64.              case 'l': return log( value_of( s, r, t, (int) t->node[idx].left ) );
  65.              case 'e': return exp( value_of( s, r, t, (int) t->node[idx].left ) );
  66.              case 'S': return sin( value_of( s, r, t, (int) t->node[idx].left ) );
  67.              case 'C': return cos( value_of( s, r, t, (int) t->node[idx].left ) );
  68.             }
  69.   case 'O': switch( t->node[idx].val )
  70.             {
  71.              case '+': return value_of( s, r, t, (int) t->node[idx].left ) +
  72.                               value_of( s, r, t, (int) t->node[idx].right );
  73.              case '-': return value_of( s, r, t, (int) t->node[idx].left ) -
  74.                               value_of( s, r, t, (int) t->node[idx].right );
  75.              case '*': return value_of( s, r, t, (int) t->node[idx].left ) *
  76.                               value_of( s, r, t, (int) t->node[idx].right );
  77.              case '/': return value_of( s, r, t, (int) t->node[idx].left ) /
  78.                               value_of( s, r, t, (int) t->node[idx].right );
  79.              case '^': return pow( value_of( s, r, t, (int) t->node[idx].left ),
  80.                                    value_of( s, r, t, (int) t->node[idx].right ) );
  81.             }
  82.  }
  83. }
  84.  
  85. double translate(double s[], int r)
  86. {
  87.  return value_of( s, r, &tree[kinetype[r]-MAX_TYP], (int) tree[kinetype[r]-MAX_TYP].node[0].left );
  88. }
  89.  
  90.  
  91. /*                                   */
  92. /* the zero returning function       */
  93. /*                                   */
  94.  
  95. double ret_zero( double s[], int r, int e )
  96. {
  97.  return (double) 0;
  98. }
  99.  
  100.  
  101. /*                                   */
  102. /* derivatives for effectors         */
  103. /* by finite differences             */
  104. /*                                   */
  105.  
  106. double deff( double s[], int r, int e )
  107. {
  108.  double k, x, x1, f1, f2;
  109.  
  110.  x = k = s[e];
  111.  /*
  112.     if x is zero, we will calculate the derivative at a small positive
  113.     value (no point in considering negative values!). let's stick with
  114.     hrcz (the flux resolution)
  115.  */
  116.  if( x==0 ) x = options.hrcz;
  117.  x1 = x * 0.01;
  118.  s[e] = x - x1;
  119.  f1 = rateq[r]( s, r );
  120.  s[e] = x + x1;
  121.  f2 = (rateq[r]( s, r ) - f1)/(x*0.02);
  122.  s[e] = k;
  123.  return f2;
  124. }
  125.  
  126. /*                                   */
  127. /* A -->  or  --> B                  */
  128. /*                                   */
  129.  
  130. double i0(double s[], int r)
  131. {
  132.  return params[r][0];
  133. }
  134.  
  135.  
  136. /*                                   */
  137. /* A --> B                           */
  138. /*                                   */
  139.  
  140. double i11(double s[], int r)
  141. {
  142.  return s[eff[r][0]]*params[r][0];
  143. }
  144.  
  145. /*                                 */
  146. /* A <--> B                        */
  147. /*                                 */
  148.  
  149. double r11(double s[], int r)
  150. {
  151.  return s[eff[r][0]]*params[r][0]-s[eff[r][1]]*params[r][1];
  152. }
  153.  
  154. /*                                 */
  155. /* A + B --> C                     */
  156. /*                                 */
  157.  
  158. double i21(double s[], int r)
  159. {
  160.  return s[eff[r][0]]*s[eff[r][1]]*params[r][0];
  161. }
  162.  
  163. /*                                 */
  164. /* A + B <--> C                    */
  165. /*                                 */
  166.  
  167. double r21(double s[], int r)
  168. {
  169.  return s[eff[r][0]]*s[eff[r][1]]*params[r][0]-s[eff[r][2]]*params[r][1];
  170. }
  171.  
  172. /*                                   */
  173. /* A --> B + C                       */
  174. /*                                   */
  175.  
  176. double i12(double s[], int r)
  177. {
  178.  return s[eff[r][0]]*params[r][0];
  179. }
  180.  
  181. /*                                 */
  182. /* A <--> B + C                    */
  183. /*                                 */
  184.  
  185. double r12(double s[], int r)
  186. {
  187.  return s[eff[r][0]]*params[r][0]-s[eff[r][1]]*s[eff[r][2]]*params[r][1];
  188. }
  189.  
  190. /*                                 */
  191. /* A + B + C --> D                 */
  192. /*                                 */
  193.  
  194. double i31(double s[], int r)
  195. {
  196.  return s[eff[r][0]]*s[eff[r][1]]*s[eff[r][2]]*params[r][0];
  197. }
  198.  
  199. /*                                 */
  200. /* A + B + C <--> D                */
  201. /*                                 */
  202.  
  203. double r31(double s[], int r)
  204. {
  205.  return s[eff[r][0]]*s[eff[r][1]]*s[eff[r][2]]*params[r][0]-s[eff[r][3]]*params[r][1];
  206. }
  207.  
  208. /*                                   */
  209. /* A --> B + C + D                   */
  210. /*                                   */
  211.  
  212. double i13(double s[], int r)
  213. {
  214.  return s[eff[r][0]]*params[r][0];
  215. }
  216.  
  217. /*                                 */
  218. /* A <--> B + C + D                */
  219. /*                                 */
  220.  
  221. double r13(double s[], int r)
  222. {
  223.  return s[eff[r][0]]*params[r][0]-s[eff[r][1]]*s[eff[r][2]]*s[eff[r][3]]*params[r][1];
  224. }
  225.  
  226. /*                                 */
  227. /* A + B --> C + D                 */
  228. /*                                 */
  229.  
  230. double i22(double s[], int r)
  231. {
  232.  return s[eff[r][0]]*s[eff[r][1]]*params[r][0];
  233. }
  234.  
  235. /*                                 */
  236. /* A + B <--> C + D                */
  237. /*                                 */
  238.  
  239. double r22(double s[], int r)
  240. {
  241.  return s[eff[r][0]]*s[eff[r][1]]*params[r][0]-s[eff[r][2]]*s[eff[r][3]]*params[r][1];
  242. }
  243.  
  244. /*                                 */
  245. /* A + B + C --> D + E             */
  246. /*                                 */
  247.  
  248. double i32(double s[], int r)
  249. {
  250.  return s[eff[r][0]]*s[eff[r][1]]*s[eff[r][2]]*params[r][0];
  251. }
  252.  
  253. /*                                 */
  254. /* A + B + C <--> D + E            */
  255. /*                                 */
  256.  
  257. double r32(double s[], int r)
  258. {
  259.  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];
  260. }
  261.  
  262. /*                                 */
  263. /* A + B --> C + D + E             */
  264. /*                                 */
  265.  
  266. double i23(double s[], int r)
  267. {
  268.  return s[eff[r][0]]*s[eff[r][1]]*params[r][0];
  269. }
  270.  
  271. /*                                 */
  272. /* A + B <--> C + D + E            */
  273. /*                                 */
  274.  
  275. double r23(double s[], int r)
  276. {
  277.  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];
  278. }
  279.  
  280. /*                                 */
  281. /* A + B + C --> D + E + F         */
  282. /*                                 */
  283.  
  284. double i33(double s[], int r)
  285. {
  286.  return s[eff[r][0]]*s[eff[r][1]]*s[eff[r][2]]*params[r][0];
  287. }
  288.  
  289. /*                                 */
  290. /* A + B + C <--> D + E + F        */
  291. /*                                 */
  292.  
  293. double r33(double s[], int r)
  294. {
  295.  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];
  296. }
  297.  
  298.  
  299. /*                                 */
  300. /* derivatives for one substrate   */
  301. /*                                 */
  302.  
  303. double d1Xs1( double s[], int r, int e )
  304. {
  305.  return params[r][0];
  306. }
  307.  
  308.  
  309. /*                                 */
  310. /* derivatives for two substrates  */
  311. /*                                 */
  312.  
  313. double d2Xs1( double s[], int r, int e )
  314. {
  315.  return params[r][0]*s[eff[r][1]];
  316. }
  317.  
  318. double d2Xs2( double s[], int r, int e )
  319. {
  320.  return params[r][0]*s[eff[r][0]];
  321. }
  322.  
  323.  
  324. /*                                 */
  325. /* derivatives for three substrates*/
  326. /*                                 */
  327.  
  328. double d3Xs1( double s[], int r, int e )
  329. {
  330.  return params[r][0]*s[eff[r][1]]*s[eff[r][2]];
  331. }
  332.  
  333. double d3Xs2( double s[], int r, int e )
  334. {
  335.  return params[r][0]*s[eff[r][0]]*s[eff[r][2]];
  336. }
  337.  
  338. double d3Xs3( double s[], int r, int e )
  339. {
  340.  return params[r][0]*s[eff[r][0]]*s[eff[r][1]];
  341. }
  342.  
  343.  
  344. /*                                 */
  345. /* derivatives for one product     */
  346. /*                                 */
  347.  
  348. double dX1p1( double s[], int r, int e )
  349. {
  350.  return -params[r][1];
  351. }
  352.  
  353.  
  354. /*                                 */
  355. /* derivatives for two products    */
  356. /* (reactions with one substrate)  */
  357. /*                                 */
  358.  
  359. double dX2p1s1( double s[], int r, int e )
  360. {
  361.  return -params[r][1]*s[eff[r][2]];
  362. }
  363.  
  364. double dX2p2s1( double s[], int r, int e )
  365. {
  366.  return -params[r][1]*s[eff[r][1]];
  367. }
  368.  
  369. /*                                 */
  370. /* derivatives for two products    */
  371. /* (reactions with two substrates) */
  372. /*                                 */
  373.  
  374. double dX2p1s2( double s[], int r, int e )
  375. {
  376.  return -params[r][1]*s[eff[r][3]];
  377. }
  378.  
  379. double dX2p2s2( double s[], int r, int e )
  380. {
  381.  return -params[r][1]*s[eff[r][2]];
  382. }
  383.  
  384. /*                                 */
  385. /* derivatives for two products    */
  386. /* (reactions with 3 substrates)   */
  387. /*                                 */
  388.  
  389. double dX2p1s3( double s[], int r, int e )
  390. {
  391.  return -params[r][1]*s[eff[r][4]];
  392. }
  393.  
  394. double dX2p2s3( double s[], int r, int e )
  395. {
  396.  return -params[r][1]*s[eff[r][3]];
  397. }
  398.  
  399.  
  400. /*                                 */
  401. /* derivatives for three products  */
  402. /* (reactions with 1 substrate)   */
  403. /*                                 */
  404.  
  405. double dX3p1s1( double s[], int r, int e )
  406. {
  407.  return -params[r][1]*s[eff[r][2]]*s[eff[r][3]];
  408. }
  409.  
  410. double dX3p2s1( double s[], int r, int e )
  411. {
  412.  return -params[r][1]*s[eff[r][1]]*s[eff[r][3]];
  413. }
  414.  
  415. double dX3p3s1( double s[], int r, int e )
  416. {
  417.  return  -params[r][1]*s[eff[r][1]]*s[eff[r][2]];
  418. }
  419.  
  420. /*                                 */
  421. /* derivatives for three products  */
  422. /* (reactions with 2 substrates)   */
  423. /*                                 */
  424.  
  425. double dX3p1s2( double s[], int r, int e )
  426. {
  427.  return -params[r][1]*s[eff[r][3]]*s[eff[r][4]];
  428. }
  429.  
  430. double dX3p2s2( double s[], int r, int e )
  431. {
  432.  return -params[r][1]*s[eff[r][2]]*s[eff[r][4]];
  433. }
  434.  
  435. double dX3p3s2( double s[], int r, int e )
  436. {
  437.  return  -params[r][1]*s[eff[r][2]]*s[eff[r][3]];
  438. }
  439.  
  440. /*                                 */
  441. /* derivatives for three products  */
  442. /* (reactions with 3 substrates)   */
  443. /*                                 */
  444.  
  445. double dX3p1s3( double s[], int r, int e )
  446. {
  447.  return -params[r][1]*s[eff[r][4]]*s[eff[r][5]];
  448. }
  449.  
  450. double dX3p2s3( double s[], int r, int e )
  451. {
  452.  return -params[r][1]*s[eff[r][3]]*s[eff[r][5]];
  453. }
  454.  
  455. double dX3p3s3( double s[], int r, int e )
  456. {
  457.  return  -params[r][1]*s[eff[r][3]]*s[eff[r][4]];
  458. }
  459.  
  460.  
  461.  
  462. /*                                                   */
  463. /* irreversible Henri-Michaelis-Menten rate equation */
  464. /*                                                   */
  465.  
  466. double ihmm( double s[], int r)
  467.           /* eff[r][0] - substrate           */
  468.           /* params[r][0] - Km for substrate    */
  469.           /* params[r][1] - max forward vel.    */
  470. {
  471.  double sb;
  472.  
  473.  sb = s[eff[r][0]];
  474.  return (sb*params[r][1])  / ( sb + params[r][0] );
  475. }
  476.  
  477. double dihmms( double s[], int r, int e )
  478. {
  479.  double c;
  480.  
  481.   c = 1 + s[eff[r][0]]/params[r][0];
  482.  return params[r][1]/(c*c*params[r][0]);
  483. }
  484.  
  485.  
  486. /*                                                 */
  487. /* reversible Henri-Michaelis-Menten rate equation */
  488. /*                                                 */
  489.  
  490. double hmm( double s[], int r)
  491.           /* eff[r][0] - substrate           */
  492.           /* eff[r][1] - product             */
  493.           /* params[r][0] - Km for substrate    */
  494.           /* params[r][1] - Km for product      */
  495.           /* params[r][2] - max forward vel.    */
  496.           /* params[r][3] - max reverse vel.    */
  497. {
  498.  double nsb, npr;
  499.  
  500.  nsb = s[eff[r][0]]/params[r][0];
  501.  npr = s[eff[r][1]]/params[r][1];
  502.  return (nsb*params[r][2] - npr*params[r][3]) / (1 + nsb + npr);
  503. }
  504.  
  505. double dhmms( double s[], int r, int e )
  506. {
  507.   double c;
  508.  
  509.   c = s[eff[r][1]]*params[r][0] + s[eff[r][0]]*params[r][1] + params[r][0]*params[r][1];
  510.   return ( params[r][0]*params[r][1] * ( s[eff[r][1]]*(params[r][2]+params[r][3]) + params[r][2]*params[r][1] ) ) / ( c*c );
  511. }
  512.  
  513. double dhmmp( double s[], int r, int e )
  514. {
  515.  double c;
  516.  
  517.  c = s[eff[r][1]]*params[r][0] + s[eff[r][0]]*params[r][1] + params[r][0]*params[r][1];
  518.  return - ( params[r][0]*params[r][1] * ( s[eff[r][0]]*(params[r][2]+params[r][3]) + params[r][0]*params[r][3] ) ) / ( c*c );
  519. }
  520.  
  521.  
  522. /*                                                                     */
  523. /* reversible Henri-Michaelis-Menten rate equation w/ competitive inh. */
  524. /*                                                                     */
  525.  
  526. double hmmsi( double s[], int r )
  527.           /* eff[r][0] - substrate           */
  528.           /* eff[r][1] - product             */
  529.           /* eff[r][2] - inhibitor           */
  530.           /* params[r][0] - Km for substrate    */
  531.           /* params[r][1] - Km for product      */
  532.           /* params[r][2] - max forward vel.    */
  533.           /* params[r][3] - max reverse vel.    */
  534.           /* params[r][4] - inhibition const.   */
  535. {
  536.  double nsb, npr, in;
  537.  
  538.  nsb = s[eff[r][0]]/params[r][0];
  539.  npr = s[eff[r][1]]/params[r][1];
  540.  in = s[eff[r][3]];
  541.  return (nsb*params[r][2] - npr*params[r][3]) / (1 + nsb + npr + in/params[r][4]);
  542. }
  543.  
  544. double dhmmsis( double s[], int r, int e )
  545. {
  546.  double sb, pr, inn, c;
  547.  
  548.  sb = s[eff[r][0]];
  549.  pr = s[eff[r][1]];
  550.  inn = s[eff[r][2]]/params[r][4];
  551.   c = 1 + sb/params[r][0] + pr/params[r][1] + inn;
  552.  return (params[r][2]*(1+inn)+pr*(params[r][2]+params[r][3])/params[r][1]) / (c*c*params[r][0]);
  553. }
  554.  
  555. double dhmmsip( double s[], int r, int e )
  556. {
  557.  double nsb, pr, in, c;
  558.  
  559.  nsb = s[eff[r][0]]/params[r][0];
  560.  pr = s[eff[r][1]];
  561.  in = s[eff[r][2]];
  562.   c = 1 + nsb + pr/params[r][1] + in/params[r][4];
  563.  return - (params[r][3]*(1+in/params[r][4])+nsb*(params[r][2]+params[r][3])) / (c*c*params[r][1]);
  564. }
  565.  
  566. double dhmmsii( double s[], int r, int e )
  567. {
  568.  double nsb, npr, in, c;
  569.  
  570.  nsb = s[eff[r][0]]/params[r][0];
  571.  npr = s[eff[r][1]]/params[r][1];
  572.  in  = s[eff[r][2]];
  573.   c = 1 + nsb + npr + in/params[r][4];
  574.  return - (nsb*params[r][2] - npr*params[r][3]) / (c*c*params[r][4]);
  575. }
  576.  
  577. /*                                                                         */
  578. /* reversible Henri-Michaelis-Menten rate equation w/ non-competitive inh. */
  579. /*                                                                         */
  580.  
  581. double hmmci( double s[], int r )
  582.           /* eff[r][0] - substrate           */
  583.           /* eff[r][1] - product             */
  584.           /* eff[r][2] - inhibitor           */
  585.           /* params[r][0] - Km for substrate    */
  586.           /* params[r][1] - Km for product      */
  587.           /* params[r][2] - max forward vel.    */
  588.           /* params[r][3] - max reverse vel.    */
  589.           /* params[r][4] - inhibition const.   */
  590. {
  591.  double nsb, npr, in;
  592.  
  593.  nsb = s[eff[r][0]]/params[r][0];
  594.  npr = s[eff[r][1]]/params[r][1];
  595.  in = s[eff[r][2]];
  596.  return (nsb*params[r][2] - npr*params[r][3]) / ((1 + nsb + npr) * (1 + in/params[r][4]));
  597. }
  598.  
  599. double dhmmcis( double s[], int r, int e )
  600. {
  601.  double npr, c;
  602.  
  603.  npr = s[eff[r][1]]/params[r][1];
  604.    c = 1 + s[eff[r][0]]/params[r][0] + npr;
  605.  return (params[r][2] + npr*(params[r][2]+params[r][3])) / (c*c*params[r][0]*(1+s[eff[r][2]]/params[r][4]));
  606. }
  607.  
  608. double dhmmcip( double s[], int r, int e )
  609. {
  610.  double nsb, c;
  611.  
  612.  nsb = s[eff[r][0]]/params[r][0];
  613.  c = 1 + nsb + s[eff[r][1]]/params[r][1];
  614.  return -(params[r][3] + nsb*(params[r][2]+params[r][3])) / (c*c*params[r][1]*(1+s[eff[r][2]]/params[r][4]));
  615. }
  616.  
  617. double dhmmcii( double s[], int r, int e )
  618. {
  619.  double nsb, npr, c, d;
  620.  
  621.  nsb = s[eff[r][0]]/params[r][0];
  622.  npr = s[eff[r][1]]/params[r][1];
  623.   c = 1 + nsb + npr;
  624.   d = 1 + s[eff[r][2]]/params[r][4];
  625.  return (npr*params[r][3] - nsb*params[r][2]) / (d*d*c*params[r][4]);
  626. }
  627.  
  628. /*                                                                     */
  629. /* reversible Henri-Michaelis-Menten rate equation w/ mixed inhibition */
  630. /*                                                                     */
  631.  
  632. double hmmmi( double s[], int r )
  633.           /* eff[r][0] - substrate            */
  634.           /* eff[r][1] - product              */
  635.           /* eff[r][2] - inhibitor            */
  636.           /* params[r][0] - Km for substrate     */
  637.           /* params[r][1] - Km for product       */
  638.           /* params[r][2] - max forward vel.     */
  639.           /* params[r][3] - max reverse vel.     */
  640.           /* params[r][4] - spec. inh. const.    */
  641.           /* params[r][5] - catal.inh. const.    */
  642. {
  643.  double sb, pr, in;
  644.  
  645.  sb = s[eff[r][0]];
  646.  pr = s[eff[r][1]];
  647.  in = s[eff[r][2]];
  648.  return (sb*params[r][2]/params[r][0] - pr*params[r][3]/params[r][1]) /
  649.         (1 + in/params[r][4] + (sb/params[r][0] + pr/params[r][1])*(1 + in/params[r][5]) );
  650. }
  651.  
  652. double dhmmmis( double s[], int r, int e )
  653. {
  654.  double sb, pr, in, c, d;
  655.  
  656.  sb = s[eff[r][0]];
  657.  pr = s[eff[r][1]];
  658.  in = s[eff[r][2]];
  659.   c = 1 + pr/params[r][1] + in/params[r][4] + pr*in/(params[r][5]*params[r][1]);
  660.   d = c + sb/params[r][0]*(1 + in/params[r][5]);
  661.  return (params[r][2]*c + params[r][3]*pr*(1+in/params[r][5])/params[r][1]) / (params[r][0]*d*d);
  662. }
  663.  
  664. double dhmmmip( double s[], int r, int e )
  665. {
  666.  double sb, in, c, d;
  667.  
  668.  sb = s[eff[r][0]];
  669.  in = s[eff[r][2]];
  670.   c = 1 + sb/params[r][0] + in/params[r][4] + sb*in/(params[r][5]*params[r][0]);
  671.   d = c + s[eff[r][1]]/params[r][1]*(1 + in/params[r][5]);
  672.  return -(params[r][3]*c + params[r][2]*sb*(1+in/params[r][5])/params[r][0]) / (params[r][1]*d*d);
  673. }
  674.  
  675.  
  676. double dhmmmii( double s[], int r, int e )
  677. {
  678.  double sb, pr, in, c;
  679.  
  680.  sb = s[eff[r][0]]/params[r][0];
  681.  pr = s[eff[r][1]]/params[r][1];
  682.  in = s[eff[r][2]];
  683.   c = 1  + in/params[r][4] + (sb + pr) * (1 + in/params[r][5]);
  684.  return -(sb*params[r][2]-pr*params[r][3]) * (1/params[r][4]+(sb+pr)/params[r][5])
  685.         / (c*c);
  686. }
  687.  
  688.  
  689. /*                                                                  */
  690. /* reversible Henri-Michaelis-Menten rate equation w/ specific act. */
  691. /*                                                                  */
  692.  
  693. double hmmsa( double s[], int r )
  694.           /* eff[r][0] - substrate            */
  695.           /* eff[r][1] - product              */
  696.           /* eff[r][2] - activator            */
  697.           /* params[r][0] - Km for substrate     */
  698.           /* params[r][1] - Km for product       */
  699.           /* params[r][2] - max forward vel.     */
  700.           /* params[r][3] - max reverse vel.     */
  701.           /* params[r][4] - activation const.    */
  702. {
  703.  double sb, pr, ac;
  704.  
  705.  sb = s[eff[r][0]];
  706.  pr = s[eff[r][1]];
  707.  ac = s[eff[r][2]];
  708.  return  ac * (sb*params[r][2]*params[r][1] - pr*params[r][3]*params[r][0])
  709.          /
  710.          ( ac * ( params[r][0]*params[r][1] + params[r][0]*pr + params[r][1]* sb )
  711.            + params[r][0]*params[r][1]*params[r][4] );
  712. }
  713.  
  714. double dhmmsas( double s[], int r, int e )
  715. {
  716.  double pr, ac, c;
  717.  
  718.  pr = s[eff[r][1]];
  719.  ac = s[eff[r][2]];
  720.   c =  ac * ( params[r][0]*params[r][1] + params[r][0]*pr + params[r][1]* s[eff[r][0]] )
  721.        + params[r][0]*params[r][1]*params[r][4];
  722.  return ac*params[r][0]*params[r][1] * (  params[r][1]*params[r][2]*(ac+params[r][4])
  723.                                         + ac*pr*( params[r][2]+params[r][3])
  724.                                        )
  725.        / ( c * c );
  726. }
  727.  
  728. double dhmmsap( double s[], int r, int e )
  729. {
  730.  double sb, ac, c;
  731.  
  732.  sb = s[eff[r][0]];
  733.  ac = s[eff[r][2]];
  734.   c =  ac * ( params[r][0]*params[r][1] + params[r][0]*s[eff[r][1]] + params[r][1]*sb )
  735.        + params[r][0]*params[r][1]*params[r][4];
  736.  return - ac*params[r][0]*params[r][1] * (  ac*sb*(params[r][2]+params[r][3])
  737.                                          + params[r][0]*params[r][3]*(ac+params[r][4])
  738.                                         )
  739.         / ( c * c );
  740. }
  741.  
  742. double dhmmsaa( double s[], int r, int e )
  743. {
  744.  double sb, pr, ac, c;
  745.  
  746.  sb = s[eff[r][0]];
  747.  pr = s[eff[r][1]];
  748.  ac = s[eff[r][2]];
  749.   c =  ac * ( params[r][0]*params[r][1] + params[r][0]*pr + params[r][1]* s[eff[r][0]] )
  750.        + params[r][0]*params[r][1]*params[r][4];
  751.  return params[r][0]*params[r][1]*params[r][4] * (  sb*params[r][1]*params[r][2]
  752.                                                   - pr*params[r][0]*params[r][3]
  753.                                                  )
  754.        / ( c * c );
  755. }
  756.  
  757.  
  758. /*                                                                         */
  759. /* reversible Henri-Michaelis-Menten rate equation w/ catalytic activation */
  760. /*                                                                         */
  761.  
  762. double hmmca( double s[], int r )
  763.           /* eff[r][0] - substrate            */
  764.           /* eff[r][1] - product              */
  765.           /* eff[r][2] - activator            */
  766.           /* params[r][0] - Km for substrate     */
  767.           /* params[r][1] - Km for product       */
  768.           /* params[r][2] - max forward vel.     */
  769.           /* params[r][3] - max reverse vel.     */
  770.           /* params[r][4] - activation const.    */
  771. {
  772.  double sb, pr, ac;
  773.  
  774.  sb = s[eff[r][0]];
  775.  pr = s[eff[r][1]];
  776.  ac = s[eff[r][2]];
  777.  return   ac * ( sb*params[r][1]*params[r][2] - pr*params[r][0]*params[r][3] )
  778.         / (   (ac+params[r][4])
  779.             * (sb*params[r][1] + pr*params[r][0] + params[r][0]*params[r][1])
  780.           );
  781. }
  782.  
  783. double dhmmcas( double s[], int r, int e )
  784. {
  785.  double ac, pr, c;
  786.  
  787.  pr = s[eff[r][1]];
  788.  ac = s[eff[r][2]];
  789.   c = s[eff[r][0]]*params[r][1] + pr*params[r][0] + params[r][0]*params[r][1];
  790.  return   ac*params[r][0]*params[r][1] * (  pr*(params[r][2]+params[r][3])
  791.                                           + params[r][1]*params[r][2]
  792.                                          )
  793.         / ( c * c * ( ac + params[r][4] ) );
  794. }
  795.  
  796. double dhmmcap( double s[], int r, int e )
  797. {
  798.  double sb, ac, c;
  799.  
  800.  sb = s[eff[r][0]];
  801.  ac = s[eff[r][2]];
  802.   c = sb*params[r][1] + s[eff[r][1]]*params[r][0] + params[r][0]*params[r][1];
  803.  return - ac*params[r][0]*params[r][1] * (  sb*(params[r][2]+params[r][3])
  804.                                           + params[r][0]*params[r][3]
  805.                                          )
  806.         / ( c * c * ( ac + params[r][4] ) );
  807. }
  808.  
  809. double dhmmcaa( double s[], int r, int e )
  810. {
  811.  double sb, pr, c;
  812.  
  813.  sb = s[eff[r][0]];
  814.  pr = s[eff[r][1]];
  815.   c = s[eff[r][2]] + params[r][4];
  816.  return   params[r][4] * ( sb*params[r][1]*params[r][2] + pr*params[r][0]*params[r][3] )
  817.         / ( c * c * ( params[r][0]*params[r][1] + pr*params[r][0] + sb*params[r][1]) );
  818.  
  819. /* WHAT IS THIS DOING HERE???
  820.  ( (params[r][4]*(sb*params[r][2]/params[r][0] - pr*params[r][3]/params[r][1])) / (c*c*(1+sb/params[r][0]+pr/params[r][1])) );
  821. */
  822. }
  823.  
  824.  
  825. /*                                                                     */
  826. /* reversible Henri-Michaelis-Menten rate equation w/ mixed activation */
  827. /*                                                                     */
  828.  
  829. double hmmma( double s[], int r )
  830.           /* eff[r][0] - substrate           */
  831.           /* eff[r][1] - product             */
  832.           /* eff[r][2] - activator           */
  833.           /* params[r][0] - Km for substrate    */
  834.           /* params[r][1] - Km for product      */
  835.           /* params[r][2] - max forward vel.    */
  836.           /* params[r][3] - max reverse vel.    */
  837.           /* params[r][4] - spec.act. const.    */
  838.           /* params[r][5] - cat. act. const.    */
  839. {
  840.  double sb, pr, ac;
  841.  
  842.  sb = s[eff[r][0]];
  843.  pr = s[eff[r][1]];
  844.  ac = s[eff[r][2]];
  845.  return   ac * (sb*params[r][2]*params[r][1] - pr*params[r][3]*params[r][0])
  846.         / (  params[r][0]*params[r][1]*(ac+params[r][4])
  847.            + (ac+params[r][5])*(sb*params[r][1]+pr*params[r][0])
  848.           );
  849. }
  850.  
  851. double dhmmmas( double s[], int r, int e )
  852. {
  853.  double pr, ac, c;
  854.  
  855.  pr = s[eff[r][1]];
  856.  ac = s[eff[r][2]];
  857.   c =  params[r][0]*params[r][1]*(ac+params[r][4])
  858.      + (ac+params[r][5])*(s[eff[r][0]]*params[r][1]+pr*params[r][0]);
  859.  return   ac*params[r][0]*params[r][1] *
  860.           (  params[r][1]*params[r][2]*(ac+params[r][4])
  861.            + pr*(ac+params[r][5])*(params[r][2]+params[r][3])
  862.           )
  863.         / (c * c);
  864. }
  865.  
  866. double dhmmmap( double s[], int r, int e )
  867. {
  868.  double sb, ac, c;
  869.  
  870.  sb = s[eff[r][0]];
  871.  ac = s[eff[r][2]];
  872.   c =  params[r][0]*params[r][1]*(ac+params[r][4])
  873.      + (ac+params[r][5])*(sb*params[r][1]+s[eff[r][1]]*params[r][0]);
  874.  return - ac*params[r][0]*params[r][1] *
  875.           (  params[r][0]*params[r][3]*(ac+params[r][4])
  876.            + sb*(ac+params[r][5])*(params[r][2]+params[r][3])
  877.           )
  878.         / (c * c);
  879. }
  880.  
  881. double dhmmmaa( double s[], int r, int e )
  882. {
  883.  double sb, pr, ac, c;
  884.  
  885.  sb = s[eff[r][0]];
  886.  pr = s[eff[r][1]];
  887.  ac = s[eff[r][2]];
  888.   c =  params[r][0]*params[r][1]*(ac+params[r][4])
  889.      + (ac+params[r][5])*(sb*params[r][1]+pr*params[r][0]);
  890.  return    (sb*params[r][1]*params[r][2] - pr*params[r][0]*params[r][3])
  891.          * (  params[r][4]*params[r][0]*params[r][0]
  892.             + params[r][5]*(sb*params[r][1]+pr*params[r][0])
  893.            )
  894.         / (c * c);
  895. }
  896.  
  897.  
  898. /*                    */
  899. /* Hill rate equation */
  900. /*                    */
  901.  
  902. double hill( double s[], int r )
  903.           /* eff[r][0] - substrate           */
  904.           /* params[r][0] - K for substrate    */
  905.           /* params[r][1] - max forward vel.    */
  906.           /* params[r][2] - Hill coef.      */
  907. {
  908.  double sb;
  909.  
  910.  sb = pow(s[eff[r][0]],params[r][2]);
  911.  return sb*params[r][1]/(pow(params[r][0],params[r][2])+sb);
  912. }
  913.  
  914. double dhills( double s[], int r, int e )
  915. {
  916.  double sb, kh, sb1, d;
  917.  
  918.  sb = pow(s[eff[r][0]], params[r][2]);
  919.  kh = pow(params[r][0], params[r][2]);
  920.  sb1 = pow(s[eff[r][0]], params[r][2] - 1);
  921.  d = ( kh + sb );
  922.  return (kh*sb1*params[r][2]*params[r][1]) / (d * d) ;
  923. }
  924.