#!/usr/bin/perl -w
use strict;

# Models glycolysis according to the method and values given
# by Eisenthal & Cornish-Bowden [(1998) J Biol Chem 273:5500-5505],
# with some help from Bakker et al [(1997) J Biol Chem
# 272:3207-3215]. The branch going off into glycerol metabolism
# has been ignored. I also fuzzed over the distinction between
# the cytoplasm and the glycosome.

################## INITIALIZATION ###################

             ###### Metabolites ######
# @values holds the current values of the variables

my $Glcx    = 5;             # glucose (external)
my $Glc     = 0;             # glucose (internal)
my $Glc6P   = 0;             # glucose-6-phosphate
my $Fru6P   = 0;             # fructose-6-phosphate
my $FruP2   = 0;             # fructose-1,6-bisphosphate
my $Gla3P   = 0;             # glyceraldehydes-3-phosphate
my $DHAP    = 0;             # dihydroxyacetone phosphate
my $GbisP   = 0;             # 1,3-diphosphoglycerate
my $ThreePG = 0;             # 3-phosphoglycerate
my $PEP     = 0;             # phosphoenolpyruvate
my $Pyr     = 0;             # pyruvate (internal)
my $Pyrx    = 0;             # pyruvate (external)
my $Pg      = 4;             # High energy P state (=2[ATP] + [ADP])
my $AXP     = 3.9;           # Total adenine nucleotides (=[AMP]+[ADP]+[ATP])
my $ATP;                     # Calculation deferred until necessary rate
my $ADP;                     # constants defined
my $AMP;
my $NADX    = 4;             # NAD + NADH pool
my $NAD     = $NADX * 0.9;   # Aribitrary choice (figured that $NAD > $NADH
my $NADH    = $NADX - $NAD;  #   under aerobic conditions

my @values = ($Glcx, $Glc, $Glc6P, $Fru6P, $FruP2, $Gla3P, $DHAP, 
              $GbisP, $ThreePG, $PEP, $Pyr, $Pyrx, $Pg);  

         ###### Collections of rate constants ######
#  @RxnN_constants holds all the constants for Reaction N
#  $RatioN holds the equilibrium constant for Reaction N

my @Rxn1_constants;
my @Rxn2_constants;
my $Ratio3;
my @Rxn4_constants;
my @Rxn5_constants;
my $Ratio6;
my @Rxn7_constants;
my @Rxn8_constants;
my $Ratio9;
my @Rxn10_constants;
my @Rxn11_constants;
my $Ratio12;
my $Ratio13;
my @slopes;

            ###### Other constants ######

my $dt = 0.0001;             # time step
my $time_steps = 3000;
my $maxnum = 10**100;        # a very large number

Set_Kinetic_Constants();
($AMP,$ADP,$ATP) = Calculate_AXP(@values);

###################### FILES ########################

my $output_path = 'glycolysis.txt';
open OUTPUT_FILE, ">$output_path" or die "Can't open $output_path: $!\n";

################### MAIN PROGRAM ####################

Output_Header();
Output(@values);
foreach my $time_step (1 .. $time_steps) {
   @slopes = model(@values);
   my @increments = multiply($dt,@slopes);
   @values = add(\@values,\@increments);
   @values = Equilibrate_Pools(@values);
   ($AMP,$ADP,$ATP) = Calculate_AXP(@values);
   Output(@values);  
}

#################### SUBROUTINES ####################

sub Set_Kinetic_Constants {
# Constants given: Vf, Vr, Ks1, Ks2, Kp1, Kp2, pwr
#    If Ks1 = -1 then use additional constants a, b, c, d, giving:
#       Ks1 = a ( 1 + ATP/b + ADP/c + AMP/d)
#    Extra constant for Rxn 1 is alpha, used per Eisenthal 1998
#    Extra constant for Rxn 5 is Ki1, used per Eisenthal 1998
# For reactions in equilibrium, ratio of P/S is given
# All nonequilibrium reactions use v_of() to calculate v, except:
#    Rxn1 and Rxn4, which use v_of_Rxn1() and v_of_Rxn4()

  #1. Glucose transport: Glcx <--> Glc
   @Rxn1_constants = (106.2, 106.2, 2, $maxnum, 2, $maxnum, 1, 0.75);
  #2. Hexokinase: ATP + Glc --> ADP + Glc6P
   @Rxn2_constants = (625, 0, 0.116, 0.1, 0.126, $maxnum, 1);
  #3. Hexose-P isomerase: Glc6P <--> Fru6P (equilibrium)
   $Ratio3 = 0.29;          
  #4. Phosphofructokinase: Fru6P + ATP <--> FruP2 + ADP
   @Rxn4_constants = (780, 0, 0.82, 0.026, $maxnum, $maxnum, 1.2); 
  #5. Aldolase: FruP2 <--> Gla3P + DHAP
   @Rxn5_constants = (184.5, 220, -1, 1, 0.067, 0.015, 1,
                         0.009, 0.68, 1.51, 3.65, 0.098);
  #6. Triose-P isomerase: DHAP <--> Gla3P (equilibrium)
   $Ratio6 = 0.045;         
  #7. Glyceraldehyde 3-P dehydrogenase: Gla3P + NAD <--> GbisP + NADH
   @Rxn7_constants = (1470, 980, 0.15, 0.45, 0.1, 0.02, 1);
  #8. 3-P-Glycerate kinase: GbisP + ADP <--> 3PG + ATP
   @Rxn8_constants = (640, 18.5, 0.05, 0.1, 1.62, 0.29, 1);
  #9. P-Glycerate mutase, enolase, PEP transport: 3PG <--> PEP (equilibrium)
   $Ratio9 = 1.25;
  #10. Pyruvate kinase: PEP + ADP <--> Pyr + ATP
   @Rxn10_constants = (2600, 0, -1, 0.114, $maxnum, $maxnum, 2.5, 
                          0.34, 0.57, 0.64, $maxnum);
  #11. Pyruvate export: Pyr <--> Pyrx
   @Rxn11_constants = (160, 0, 1.96, $maxnum, $maxnum, $maxnum, 1);
  #12. Adenylate kinase (myokinase): ADP <--> ATP + AMP (equilibrium)
   $Ratio12 = 1.25;
  #13. Growth (ATP utilization)
   $Ratio13 = 50;
}

sub Calculate_AXP {
# Taken from Bakker et al (1997), Equations 41, 21, and 1
# Based on:  2ADP <--> ATP + AMP    (leading to equation 41)
#            Pg = 2ATP + AMP        (equation 21)
#            AXP = AMP + ADP + ATP  (equation 1)

   my ($Glcx, $Glc, $Glc6P, $Fru6P, $FruP2, $Gla3P, $DHAP,  
              $GbisP, $ThreePG, $PEP, $Pyr, $Pyrx, $Pg) = @_;  
   my $a = 1 - 4 * $Ratio12;
   my $b = $AXP - $Pg * $a;
   my $c = - $Ratio12 * $Pg * $Pg;
   my $ATP = (-$b + sqrt($b*$b - 4*$a*$c)) / (2*$a);
   my $ADP = $Pg - 2 * $ATP;
   my $AMP = $AXP - $ATP - $ADP;
   return ($AMP, $ADP, $ATP);
}

sub model {

   my ($Glcx, $Glc, $Glc6P, $Fru6P, $FruP2, $Gla3P, $DHAP,  
              $GbisP, $ThreePG, $PEP, $Pyr, $Pyrx, $Pg) = @_;  

   my $v_Rxn1 = v_of_Rxn1($Glcx, -1, $Glc, -1, @Rxn1_constants);
   my $v_Rxn2 = v_of($ATP, $Glc, $ADP, $Glc6P, @Rxn2_constants);
   my $v_Rxn4 = v_of($Fru6P, $ATP, -1, -1, @Rxn4_constants);
   my $v_Rxn5 = v_of_Rxn5($FruP2, -1, $Gla3P, $DHAP, @Rxn5_constants);
   my $v_Rxn7 = v_of($Gla3P, $NAD, $GbisP, $NADH, @Rxn7_constants);
   my $v_Rxn8 = v_of($GbisP, $ADP, $ThreePG, $ATP, @Rxn8_constants);
   my $v_Rxn10 = v_of($PEP, $ADP, $Pyr, $ATP, @Rxn10_constants);
   my $v_Rxn11 = v_of($Pyr, -1, $Pyrx, -1, @Rxn11_constants);
   my $v_Rxn13 = $Ratio13 * $ATP/$ADP;

   my $dGlcx_dt = -$v_Rxn1;
   my $dGlc_dt = +$v_Rxn1 - $v_Rxn2;
   my $dGlc6P_dt = $v_Rxn2;
   my $dFru6P_dt = -$v_Rxn4;
   my $dFruP2_dt = +$v_Rxn4 - $v_Rxn5;
   my $dGla3P_dt = +$v_Rxn5 - $v_Rxn7;
   my $dDHAP_dt = +$v_Rxn5;
   my $dGbisP_dt = +$v_Rxn7 - $v_Rxn8;
   my $d3PG_dt = +$v_Rxn8;
   my $dPEP_dt = -$v_Rxn10;
   my $dPyr_dt = +$v_Rxn10 - $v_Rxn11;
   my $dPyrx_dt = +$v_Rxn11;
   my $dPg_dt = -$v_Rxn2 - $v_Rxn4 + $v_Rxn8 + $v_Rxn10 - $v_Rxn13;

   return ($dGlcx_dt, $dGlc_dt, $dGlc6P_dt, $dFru6P_dt, $dFruP2_dt, $dGla3P_dt,
      $dDHAP_dt, $dGbisP_dt, $d3PG_dt, $dPEP_dt, $dPyr_dt, $dPyrx_dt, $dPg_dt);  
}

sub v_of { 
   my ($S1, $S2, $P1, $P2, $Vf, $Vr, $Ks1, $Ks2, $Kp1, $Kp2, $pwr, 
       $a, $b, $c, $d) = @_;
   if ($Ks1 == -1) {$Ks1 = $a * (1 + $ATP/$b + $ADP/$c + $AMP/$d)}
   my $term1 = ($Vf * $S1**$pwr) / ($Ks1**$pwr); 
   if ($S2 >= 0) {$term1 =  $term1 * $S2 / $Ks2}
   my $term2 = $Vr * $P1 / $Kp1;
   if ($P2 >= 0) {$term2 = $term2 * $P2 / $Kp2}
   my $numerator =  $term1 - $term2;
   my $denominator = (1 + $S1**$pwr/$Ks1**$pwr + $P1/$Kp1)
                   * (1 + $S2/$Ks2 + $P2/$Kp2);
   return $numerator/$denominator;
}

sub v_of_Rxn1 { 
   my ($S1, $S2, $P1, $P2, $Vf, $Vr, $Ks1, $Ks2, $Kp1, $Kp2, $pwr, 
       $a) = @_;
   my $numerator =  ($Vf * $S1**$pwr) / ($Ks1**$pwr)
                  - ($Vr * $P1) / ($Kp1);
   my $denominator = 1 + $S1/$Ks1 + $P1/$Kp1 + $a*$S1*$P1/($Ks1*$Kp2);
   return $numerator/$denominator;
}

sub v_of_Rxn5 { 
   my ($S1, $S2, $P1, $P2, $Vf, $Vr, $Ks1, $Ks2, $Kp1, $Kp2, $pwr, 
       $a, $b, $c, $d, $Ki1) = @_;
   if ($Ks1 == -1) {$Ks1 = $a * (1 + $ATP/$b + $ADP/$c + $AMP/$d)}
   my $numerator =  ($Vf * $S1) / ($Ks1)
                  - ($Vr * $P1 * $P2) / ($Kp1 * $Kp2);
   my $denominator = 1 + $S1/$Ks1 + $P1/$Kp1 + $P2/$Kp2
                       + $S1*$P1/($Ks1*$Ki1) + $P1*$P2/($Kp1*$Kp2);
   return $numerator/$denominator;
}

sub Equilibrate_Pools {
   my ($Glcx, $Glc, $Glc6P, $Fru6P, $FruP2, $Gla3P, $DHAP, 
              $GbisP, $ThreePG, $PEP, $Pyr, $Pyrx, $Pg) = @_;  

   my $Sum3 = $Glc6P + $Fru6P;
   $Glc6P = $Sum3/(1 + $Ratio3);
   $Fru6P = $Sum3 - $Glc6P;
   my $Sum6 = $DHAP + $Gla3P;
   $DHAP = $Sum6/(1 + $Ratio6);
   $Gla3P = $Sum6 - $DHAP;
   my $Sum9 = $ThreePG + $PEP;
   $ThreePG = $Sum9/(1 + $Ratio9);
   $PEP = $Sum9 - $ThreePG;

   return ($Glcx, $Glc, $Glc6P, $Fru6P, $FruP2, $Gla3P, $DHAP,  
              $GbisP, $ThreePG, $PEP, $Pyr, $Pyrx, $Pg);
}   

sub Output {
  print join("\t", @_), "\n";
  print OUTPUT_FILE join("\t",@_), "\n";
}

sub Output_Header {
  print OUTPUT_FILE join("\t", ("Glcx", "Glc", "Glc6P", "Fru6P", 
        "FruP2", "Gla3P", "DHAP", "GbisP", "3PG", "PEP", "Pyr", 
        "Pyrx", "Pg")), "\n";  
}

sub multiply {
  my ($scalar,@array) = @_;
  my @product;
  foreach my $i (0 .. @array - 1) {
     $product[$i] = $scalar * $array[$i];
  }
  return @product;
}

sub add {
# Adds two arrays together. Negative numbers set to 0.
  my ($array1ref,$array2ref) = @_;
  my @array1 = @$array1ref;
  my @array2 = @$array2ref;
  my @sum;
  foreach my $i (0 .. @array1 - 1) {
     $sum[$i] = $array1[$i] + $array2[$i];
     if ($sum[$i] < 0) {$sum[$i] = 0}
  }
  return @sum;
}
