RNfb2T annot.1d1

Program:  RNfb   Version 2T   by Kurt W. Kohn,  May 8, 1997
Language: FutureBASIC II  2.1.1 for Macintosh (STAZ Software)

This is an annotated copy of the program RNfb.2T which was used to generate data presented in 
Kohn KW: Functional capabilities of molecular network components controlling the mammalian G1/S cell cycle phase transition. Oncogene  16:1065-1075, 1998. 
 
*********************************************************************
  This program was written by 
   Kurt W. Kohn,  Laboratory of Molecular Pharmacology,
   Division of Basic Sciences, National Cancer Institute,
   Building 37, Room 5C25, NIH, Bethesda, Maryland, 20892
   Telephone:  301-496-2769    FAX:  301-402-0752
   E-mail  kohnk@dc37a.nci.nih.gov
   The author hereby grants permission to use, copy, modify, or distribute
   this software and its documentation, provided that this notice is included 
   verbatim in any distributions, and that the author receives written notification 
   of any new distributions and distributed modifications.
   IN NO EVENT SHALL THE AUTHOR BE LIABLE TO ANY PARTY FOR 
   DIRECT,  INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIOAL DAMAGES
   ARISING OUT OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, 
   OR ANY DERIVATIVES THEREOF.  THIS SOFTWARE IS PROVIDED ON AN
   "AS IS" BASIS, AND THE AUTHOR HAS NO OBLIGATION TO PROVIDE
    MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
***********************************************************************

This program simulates reaction networks as defined in a reaction file that is read by the program.  The output of this program is written to a file named "RNtout" in comma-separated format.  In addition, graphs are created in which each plot is scaled to the same height. To get properly scaled graphs of multiple plots, another program (called "Plotout") is run which gets input from the RNtout file.  Instructions for use of the programs are given in the RNhelp file included on the disk. 

' -------------------- Globals --------------------------

RxnFi$=""           'Name of the reaction file to be used in the simulation.  
                           'The program will allow the user to select the desired file.

'The following variables are initialized merely to declare them as GLOBAL;  
'  the actual values used will be derived on the basis of input from the reaction file:

nrxns%=12           'Number of reactions.  
nspecies%=12        'Number of molecular species.  
dt1! =.1                'Time interval between plotted points
nsteps%=100           'Number of time steps plotted
nsubsteps%=10        'Number of calculated steps between successive plotted steps. 
                              '     (total steps = nsteps%*nsubsteps%)
                              '     (total time = nsteps%*nsubsteps%*dt = nsteps%*dt1)
dt!=dt1!/nsubsteps%    'Time step for each calculated reaction step

itemtype1%=0       'See LOCAL FN Scan2 for explanation of these variables
scanitem1%=0
mode1%=0  
itemtype2%=0     'itemtype2% and scanitem2% provide nesting capability
scanitem2%=0     '  that is not utilized in the current version of the program.
                                                
nplots%=1            'Number of plots to be recorded

'The following graphics constants will be recalculated during execution of the program:

hscale!=1
gheight1%=0

DEF FN showTopics USING gShowTopics&

v0%= 200                     'Distance of horizontal axis from top
hmax%= 1100               'Width of horizontal axis

black%=1                

tout%=1                                           'flag to record kinetic curves in RNtout.

_gtop=10                                          'Alternative values:  10  20
_gwidth=950                                                                '640  950
_gheight=850                                                             '850   575  394
glines!=(_gheight-_gtop)/12.5

xx!=0                                             'x scan parameter
yy!=0                                             'y scan parameter
iplot%=0                                         'index of plots outputed to RNtout

_dim1=26
_dim2=600
_dim3=800
DIM k!(_dim3)                       'Rate constants
DIM x!(_dim3)                       'Concentrations
DIM x0!(_dim3)                      'Initial concentrations
DIM plotsp%(_dim1)
DIM xt!(_dim1,_dim2)            'Concentration (selected species , time)
DIM xs!(_dim1)                      'List of selected species
DIM xmax!(_dim1)                  'Max conc of selected species

DIM i1%(_dim3),i2%(_dim3),i3%(_dim3)           'Input species
DIM o1%(_dim3),o2%(_dim3),o3%(_dim3)         'Output species
DIM e%(_dim3)                                                 ' "e" species
DIM vscale!(_dim3)
DIM d!(_dim3)                'List of concentration changes for the species in a calculated reaction step
DIM x$(_dim3)                'Names of molecular species
DIM tDelayk!(_dim3)       'Delay in the start of a reaction --  t when k!() stops being 0

'The following variable are not utilized by the current version of the program, but are best left as is:
tmax! =100
cells%=2     

END GLOBALS

' ----------- Functions -------------

'Input file is #1 (reaction file)
'Ouptut file is #3 (RNtout file)
'File #2 can be ignored for the current simulations

LOCAL FN ReadRxns

  PRINT #3,"Reactions:"
  nrxns%=0

  DO
    nrxns%=nrxns%+1
    INPUT#1,j%

    INPUT#1,i1%(j%),i2%(j%),i3%(j%),o1%(j%),o2%(j%),o3%(j%),e%(j%),a!,q$

                   'i1%(j%) is the id number of the first molecular input species for reaction j%
                   'Similarly for 3 input species, 3 output species, and an e% species.
                   'a! is a rate constant, later to be assigned to a k!()
                   'If a! is negative, the rate constant will be set equal to a! times the k of reaction e%;
                   '  this is a convenient way to link together the rate constants of several reactions.

    LONG IF e%(j%)<0                              '2p
      k!(j%) = a! * k!(ABS(e%(j% )))
    XELSE
      k!(j%)=a!
    END IF 
                                       
    PRINT j%;i1%(j%);i2%(j%);i3%(j%);o1%(j%);o2%(j%);o3%(j%);e%(j%);k!(j%),q$
    PRINT #3,j%;",";i1%(j%);",";i2%(j%);",";i3%(j%);",";o1%(j%);",";o2%(j%);",";o3%(j%);",";e%(j%);",";k!(j%);",";q$
  UNTIL q$="end"

  nrxns%=nrxns%-1
  PRINT "nrxns%=";nrxns%
  PRINT #3,"nrxns%=";nrxns%
END FN

LOCAL FN ReadSpecies
  PRINT #3, "Species:"
  k%=0                                      'counter for species to be plotted; k% -> nplots% & index for plotsp%()
  FOR j%=1 TO _dim1
    plotsp%(j%)=0
  NEXT j%
  q$=""
  nspecies%=0

  DO
    nspecies%=nspecies%+1
    INPUT#1,i%
    j%=ABS(i%)
    INPUT#1,x0!(j%),q!,q$               '(q! is not utilized by the current version of the program)
    x$(j%)=q$
    PRINT j%;x0!(j%);q!,x$(j%)
    PRINT #3,j%;",";x0!(j%);",";q!;",";x$(j%)
    LONG IF i%<0
      k%=k%+1
      plotsp%(k%)=j%                              'j% is species number in plot k%
    END IF
  UNTIL q$="end"

  nspecies%=nspecies%-1
  nplots%=k%
  PRINT "nspecies%=";nspecies%;"  nplots%=";nplots%
  PRINT #3,"nspecies%=";nspecies%
END FN

LOCAL FN InitializeSpecies
  FOR i%=0 TO _DIM1
    FOR j%=0 TO _DIM2
      xt!(i%,j%)=0
    NEXT j%
  NEXT i%
  FOR j%=1 TO nspecies%
    x!(j%) = x0!(j%)                              'initial concs
    IF x!(j%)=0 THEN x!(j%)=.000001
  NEXT j%
  FOR j%=1 TO nrxns%
    IF tDelayk!(j%) > 0 THEN k!(j%)=-ABS(k!(j%))
  NEXT j%
END FN

LOCAL FN Axes
  WINDOW 2
  CLS
  TEXT _monaco,10,,_srcCopy
  PLOT 0,_gtop TO _gwidth,_gtop
  FOR i%=1 TO nplots%
    k%=plotsp%(i%)                                'k% is species number in plot i%
    h%=0
    gheight1%=(_gheight-_gtop)/nplots%
    v0%= _gtop + gheight1% * i%
    PLOT 0,v0% TO _gwidth,v0%
    iline% = 1+INT(0.5+gheight1%*(i%-1)/_gheight * glines!)
    LOCATE 50,iline%
    PRINT k%;
    PRINT USING "#.###^^^^";xmax!(i%);
    PRINT USING " #.###^^^^ ";x!(k%);
    PRINT x$(k%);
    IF tout%=1 THEN PRINT #3,k%;",";
    'IF k%=endptSpecies% THEN PRINT USING"##.###"; endptTime!
  NEXT i%
  PRINT #3, " species"
END FN

LOCAL FN Splot
  WINDOW 2
  FOR i%=1 TO nplots%
    PLOT 0,v0%
    gwid!=_gwidth
    PRINT #3,"//"
    iplot%=iplot%+1
    PRINT #3,iplot%
    PRINT #3,"species";",";plotsp%(i%)
    PRINT #3,"xx! yy!= ,";xx!;",";yy!
    FOR j%=1 TO nsteps%
      nst!=nsteps%
      h%=INT(j%*gwid!/nst!)
      v0%= _gtop + gheight1% * i%
      temp1!=(xt!(i%,j%) / xmax!(i%))
      v%=v0%- temp1! * gheight1%
      
      PLOT h%,v%                      
      IF tout%=1 THEN PRINT #3,xt!(i%,j%);",";    'v%;v0%;i%;j%;xt!(i%,j%);xmax!(i%);temp1!;gheight1%;
    NEXT j%
    
    k%=plotsp%(i%)
    iline% = 1+INT(0.5+gheight1%*(i%-1)/_gheight * glines!)
    LOCATE 50,iline%
    PRINT SPACE$(30);
    LOCATE 50,iline%
    PRINT k%;
    PRINT USING "#.###^^^^";xmax!(i%);
    PRINT USING " #.###^^^^ ";x!(k%);
    PRINT x$(k%);
    
    IF tout%=1 THEN PRINT #3,-99.9
  NEXT i%
END FN

'**********************************************
'NOW FOR THE HEART OF THE COMPUTATION:

LOCAL FN DoCell
                                                  
  totalsteps%=nsteps%*nsubsteps%
  t! = 0                                          'time
  tstep%=0
  rstep%=0                                        'sequential number of recorded step
  's1%,s2%,s3%,s4% -- molecular species
  'd(rxn) -- change in conc (molec/cell)
  FOR j%=1 TO nplots%
    xmax!(j%)=0
  NEXT j%

  WHILE tstep% < totalsteps%
    t!=t!+dt!                                     'Next time step
    tstep%=tstep%+1

    'DETERMINE THE EXTENTS OF REACTIONS   d(r%):

    FOR r%=1 TO nrxns%
      LONG IF k!(r%) < 0
        IF t! > tDelayk!(r%) THEN k!(r%)=ABS(k!(r%))
      END IF
      LONG IF k!(r%) > 0
        IF r%=kscan% THEN  k!(r%) = EXP( logkmin! + (t!/tmax!) * (logkmax! - logkmin!))
        
        s1% = ABS(i1%(r%))                        '1st input species for rxn r%
        d!(r%) = k!(r%)*x!(s1%)*dt!               'every rxn is at least 1st order
        s2% = ABS(i2%(r%))
        IF s2%>0 THEN  d!(r%)=d!(r%)*x!(s2%)      'if rxn is 2nd order
        s3% = ABS(i3%(r%))
        IF s3%>0 THEN d!(r%)=d!(r%)*x!(s3%)       'if rxn has 3 input species
        s4% = e%(r%)                              'enzyme species
        IF s4%>0 THEN  d!(r%)=d!(r%)*x!(s4%)      'If rxn is enzyme mediated
      END IF
    NEXT r%

    FOR r%=1 TO nrxns%

          'DECREMENT THE INPUT SPECIES (unless flagged negative):

      LONG IF k!(r%) > 0
        ss1%=i1%(r%)                              'Input species of rnx r%
        s1%=ABS(ss1%)
        ss2%=i2%(r%)
        s2%=ABS(ss2%)
        ss3%=i3%(r%)
        s3%=ABS(ss3%)
        LONG IF ss1% >0
          IF d!(r%) > x!(s1%) THEN d!(r%) = x!(s1%) - 1e-12
        END IF
        LONG IF ss2% >0
          IF d!(r%) > x!(s2%) THEN d!(r%) = x!(s2%) - 1e-12
        END IF
        LONG IF ss3% >0
          IF d!(r%) > x!(s3%) THEN d!(r%) = x!(s3%) - 1e-12
        END IF
        
        IF ss1%>0 THEN x!(s1%)=x!(s1%)-d!(r%)     'Don't deplete if flagged negative
        IF ss2%>0 THEN x!(s2%)=x!(s2%)-d!(r%)
        IF ss3%>0 THEN x!(s3%)=x!(s3%)-d!(r%)
        
        'INCREMENT THE OUTPUT SPECIES (unless flagged negative):
        
        s1%=o1%(r%)                               'Now output species of the rxn
        s2%=o2%(r%)
        s3%=o3%(r%)
        IF s1%>0 THEN x!(s1%)=x!(s1%)+d!(r%)
        IF s2%>0 THEN x!(s2%)=x!(s2%)+d!(r%)
        IF s3%>0 THEN x!(s3%)=x!(s3%)+d!(r%)
      END IF
    NEXT r%
    
    LONG IF tstep% MOD nsubsteps% = 0
      rstep%=rstep%+1                             'sequential number of step recorded
      LOCATE 60,4
      CLS LINE 
      PRINT rstep%
      FOR j%=1 TO nplots%
        k%=plotsp%(j%)                            'k% is species number in plot j%
        xt!(j%, rstep%) = x!(k%)
        IF xmax!(j%)<x!(k%) THEN xmax!(j%)=x!(k%)
      NEXT j%
    END IF

  WEND

END FN

LOCAL FN UpdateLinkedConstants                    
  FOR i%=1 TO nspecies%
    IF e%(i%)<0 THEN k!(i%) = k!(ABS(e%(i%)))
  NEXT i%
END FN

'****************************************************************************
'THE FOLLOWING ROUTINE RUNS A SERIES OF SIMULATIONS FOR A LINEAR OR
'  EXPONENTIAL SEQUENCE OF VALUES OF A RATE CONSTANT OR INITIAL CONCENTRATION
'  OR DELAY TIME:
'  (provision is made for 2 embedded sequences, but this capability was not utilized in the current simulations) 

LOCAL FN Scan2
  WINDOW 2
  graph%=1
  FN Axes
  'WINDOW 1
  OPEN "O",2,"Scout"            'This file is irrelevant to the current simulations
  INPUT #1,sp%                    'The value of sp% is irrelevant to the current simulation

  INPUT #1,itemtype1%,scanitem1%,mode1%,x1!,x2!,xd!

     'If itemtype1%=0 then scanitem1% refers to an initial concentration
     '                      1                                       a rate constant
     '                      2                                       a delay time
     'x1! and x2! are the initial and final values of the item to be scanned
     'xd! is the step interval between successive values of the scanned item
     'If mode1%=0 then xd! is added to the old value to get the new value
     'If mode1%=1 then the old value is multiplied by xd! 

  INPUT #1,itemtype2%,scanitem2%,mode2%,y1!,y2!,yd!  'This line is not utilized further in the current simulation
    
  PRINT #3,sp%
  PRINT #3,itemtype1%;",";scanitem1%;",";mode1%;",";x1!;",";x2!;",";xd!
  PRINT #3,itemtype2%;",";scanitem2%;",";mode2%;",";y1!;",";y2!;",";yd!
  PRINT #3,"**"
  PRINT #3, RxnFi$
  PRINT #3, "itemtype,scanitem"
  PRINT #3, itemtype1%;",";scanitem1%
  PRINT #3, itemtype2%;",";scanitem2%
  
  yy!=y1!
  PRINT #2, y1!;",";
  DO
    LONG IF mode2%=1
      yy!=yy!*yd!
    XELSE
      yy!=yy!+yd!
    END IF
    PRINT #2, yy!;",";
  UNTIL yy!>=y2!
  PRINT #2,
  xx!=x1!
  IF xx!=0 THEN xx!=.000001
  DO
    yy!=y1!
    IF yy!=0 THEN yy!=.000001
    PRINT #2, xx!;",";
    DO
      IF itemtype1%=0 THEN x0!(scanitem1%)=xx!
      IF itemtype1%=1 THEN k!(scanitem1%)=xx!
      IF itemtype1%=2 THEN tDelayk!(scanitem1%)=xx!
      
      IF itemtype2%=0 THEN x0!(scanitem2%)=yy!
      IF itemtype2%=1 THEN k!(scanitem2%)=yy!
      IF itemtype2%=2 THEN tDelayk!(scanitem2%)=yy!
      
      FN InitializeSpecies
      FN DoCell
      
      IF graph%=1 THEN FN Splot
      LOCATE 60,3
      CLS LINE
      PRINT xx!;yy!,
      'PRINT #3, "end ,";xx!;",";yy!
      PRINT #2, x!(sp%);",";
      
      key$=INKEY$
      LONG IF key$="." 
        INPUT q$
        IF q$="." THEN STOP
      XELSE
        IF key$="g" THEN graph%=0
      XELSE
        IF key$="G" THEN graph%=1
      END IF
      
      LONG IF mode2%=1
        yy!=yy!*yd!
      XELSE
        yy!=yy!+yd!
      END IF
      FN UpdateLinkedConstants                    
    UNTIL yy!>y2!
    PRINT #2,
    LONG IF mode1%=1
      xx!=xx!*xd!
    XELSE
      xx!=xx!+xd!
    END IF
  UNTIL xx!>x2!
  CLOSE #2
END FN

'****************************************
LOCAL FN Start
  q$="fudge"
  WHILE q$<>""
    FN DoCell
    FN Axes
    FN Splot
    q$=""
    INPUT ""; q$
    IF q$="." THEN STOP
    IF q$="c" THEN q$="fudge"
  WEND
END FN

LOCAL FN Triage          'For the current simulations, one could skip directly to FN Scan2
  WHILE 1
    q$=""
    INPUT#1,mode1$                  'In the current simulations, mode1$ was always "scan2"
    PRINT mode1$
    SELECT
      CASE mode1$="set"
        INPUT #1,j%
        INPUT #1,tDelayk!(j%)
      CASE mode1$="rxn"
        INPUT#1,j%
        INPUT#1,k!(j%)
        PRINT j%;i1%(j%);i2%(j%);i3%(j%);o1%(j%);o2%(j%);o3%(j%);e%(j%);k!(j%);
        INPUT;fudge$
        FN Start
      CASE mode1$="species"
        INPUT#1,i%
        j%=ABS(i%)
        INPUT#1,x!(j%)
        FN Start
      CASE mode1$="continue"
        FN Start
      CASE mode1$="scan species" OR mode1$="logscan species"
        OPEN "O",2,"RNout"
        INPUT #1, sp%,x1!,x2!,xd!
        PRINT q$;sp%;x1!;x2!;xd!
        INPUT; fudge$
        FOR j%=0 TO nspecies%
          PRINT #2, USING "###     &"; j%,",",
        NEXT j%
        PRINT #2, fudge%
        xx!=x1!
        IF xx!=0 THEN xx!=.000001
        DO
          LONG IF mode1$="logscan species"
            xx!=xx!*xd!
          XELSE
            xx!=xx!+xd!
          END IF
          LOCATE 1,1,1
          PRINT xx!
          FN InitializeSpecies
          x0!(sp%)=xx!
          IF x0!(sp%)=0 THEN x0!(sp%)=.000001
          x!(sp%)=x0!(sp%)
          FN DoCell         
          PRINT #2, USING "##.### &"; x0!(sp%),",",
          FOR j%=1 TO nspecies%
            PRINT #2, USING "##.#### &"; x!(j%),",",
          NEXT j%
          PRINT #2, fudge%
        UNTIL xx!>=x2!
        CLOSE #2

      CASE mode1$="scan2"      'THIS IS THE ONLY CASE RELEVANT TO THE CURRENT SIMULATIONS
        FN Scan2
        BEEP
        INPUT; q$

      CASE ELSE                               
        CLOSE #1
        INPUT q$
        STOP
    END SELECT
  WEND
END FN

'*************************************************************
'-------Main--------
WINDOW -2,,(1,_gtop)-(_gwidth,_gheight),_dialogPlain,_updateVisRgn
WINDOW  1,,(1,_gtop)-(_gwidth,_gheight),_dialogPlain,_updateVisRgn
IF tout%=1 THEN OPEN "O",3,"RNtout"
WHILE 1
  RxnFi$=FILES$(_fOpen,,"Open a rxnFile",volRefNum%)
  OPEN "I",1,RxnFi$,,volRefNum%
  lin$=""
  LOCATE 0,2
  PRINT "OPENed  ";RxnFi$;"   volRefNum%=";volRefNum%
  
  WHILE lin$<>"*"          'Skip comment lines until a line consisting of * is detected:
    INPUT#1,lin$
    PRINT lin$
    PRINT #3, lin$
  WEND
  
  INPUT#1,dt1!,nsteps%,nsubsteps%
  
              'dt1! = time interval between points to be plotted
              'nsteps% = number of time steps of size dt1! to be iterated in the simulation 
              'Hence the time covered in each simulation is dt1!*nsteps%         
              'nsubsteps% = number of calculated time steps within each interval dt1!
              'Thus the time interval for each calculated step is dt1!/nsubsteps%

  PRINT #3,dt1!;",";nsteps%;",";nsubsteps%
  dt!=dt1!/nsubsteps%
  PRINT "dt1!,dt!=";dt1!;dt!;"  nsteps%=";nsteps%,"  nsubsteps%=",nsubsteps%
  q$=""
  
  FN ReadRxns
  PRINT "dt1!,dt!=";dt1!;dt!;"  nsteps%=";nsteps%,"  nsubsteps%=",nsubsteps%
  INPUT; q$
  IF q$="." THEN STOP
  FN ReadSpecies
  PRINT "dt1!,dt!=";dt1!;dt!;"  nsteps%=";nsteps%,"  nsubsteps%=",nsubsteps%
  INPUT; q$
  IF q$="." THEN STOP
  FN InitializeSpecies
  FN Triage
  WINDOW 1
  INPUT "END RxnNet program?";q$
  IF q$="y" THEN STOP
  CLOSE #1
  IF tout%=1 THEN CLOSE #3
  CLS
  WINDOW -2
  CLS
WEND
END




1