#!/usr/bin/perl -w
# ADP-As:
#    VERSION: Version 2 (13 October 2003)
#    PURPOSE: Simulates the decay of ADP-Arsenate to ADP + Arsenate
#             Reaction is unidirectional, so the rate equation is:
#                  d[P]/dt = -d[S]/dt = k[S]        
#             where S is ADP-Arsenate and P is ADP
#    OUTPUT FILES: Tab-delimited text file with the following three field:
#         time (in seconds)       [substrate]       [product]

# Reaction to be modeled is S --> P

############## LIBRARIES AND PRAGMAS ################
  use strict;

#################### CONSTANTS ######################

  my $k = 1;                 # Rate constant, in sec/mM
  my $dt = .5;               # Time increment, in sec
  my $duration = 5;          # Duration of simulation, in sec
  my $LF = "\n";             # Linefeed
  my $tab = "\t";            # Tab

#################### VARIABLES ######################

  my $S = 1;                 # Substrate concentration (mM)
  my $P = 0;                 # Product concentration (mM)
  my $t;                     # Current time
  my @slopes;                # Holds d[ ]/dt for each metabolite
  my @metabolites = ($S,$P); # Holds [ ]  for each metabolite  
  my @increments;            # Holds amount by which each metabolite increases

###################### FILES ########################

  my $output_path = 'ADPAs.txt';
  open OUTPUT_FILE, ">$output_path" or die "Can't open $output_path: $!\n";

################### MAIN PROGRAM ####################
# Time goes from 0 to duration of simulation
# At each time point, the rates of changes are calculated from the model,
#    then the amount of change of each metabolite is calculated, and
#    finally the concentrations of the metabolites are updated.

  Output_Header();

  for ($t=0; $t <= $duration; $t = $t + $dt) {
     Output($t,@metabolites);
     @slopes = model(@metabolites);
     @increments = multiply($dt,@slopes);
     @metabolites = add(\@metabolites,\@increments);
  }

#################### SUBROUTINES ####################

#### MODEL
#    Calculates rate of change for each metabolite
#    Reaction is S --> P
  sub model {
    my ($S,$P) = @_;
    my $S_prime = -$k * $S;
    my $P_prime =  $k * $S;
    return ($S_prime, $P_prime);
  }

#### OUTPUT_HEADER
#    Prints column headers for time an metabolites
  sub Output_Header {
    print OUTPUT_FILE "Time", $tab, "Substrate", $tab, "Product", $LF;
  }

#### OUTPUT
#    Prints time and values of metabolites
  sub Output {
    print join("\t", @_), "\n";
    print OUTPUT_FILE join("\t",@_), "\n";
  }

#### MULTIPLY (scalar, array)
#    Multiplies each element of the array by the scalar value
#    Returns array holding product of multiplication
  sub multiply {
    my ($scalar,@array) = @_;
    my @product;
    foreach my $i (0 .. @array - 1) {
       $product[$i] = $scalar * $array[$i];
    }
    return @product;
  }

#### ADD (array1, array2)
#    Adds each element of an array to the corresponding element of another array
#    Sets any negative sum to zero
#    Returns array holding sum of addition
#    Arguments are references to the two arrays
  sub add {
    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;
  }
