#################################################################
#    Programs Related to Nonlinear Least Squares in FORTRAN77   #
#                        Version 3.5                            #
#################################################################
        
                        Kiyoshi Yamaoka
        Faculty of Graduate School of Pharmaceutical Sciences
              Kyoto University, Kyoto 606, Japan
        
[1] Introduction
We have developed several programs for curve fitting by non-linear least squares.  ALL programs
 are written in Microsoft BASIC(1-7). Unfortunately, the syntax of the BASIC on personal computers
 is different from those on UNIX machines and mainframe computers.  Therefore, these programs were
 translated from BASIC to FORTRAN.  It was confirmed all programs run on IBM compatible machines
(i.e. on Windows). 
        
[2] Comments on Program Files
This diskette (3.5') for the standard IBM-PC  (720K) contains one document file and six programs
 written in FORTRAN.
        
MULTI.TXT: this file
SULTIF.FOR:simulation program using direct definition of model equations
MULTIF.FOR:curve fitting program using direct definition of model equations
SRUNGEF.FOR:simulation program using ordinary differential equations
MRUNGEF.FOR:curve fitting program using ordinary differential equations
SFILTF.FOR:simulation program using fast inverse Laplace transform (FILT)
MFILTF.FOR:curve fitting program fast inverse Laplace transform

The first letters S and M in file name mean simulation and curve fitting, respectively.  The files,
 MULTIF.FOR, MRUNGEF.FOR and MFILTF.FOR correspond to the programs, MULTI, MULTI(RUNGE), and 
MULTI(FILT), respectively.
        
[3] Definitions of Pharmacokinetic Models
    Let's consider the following pharmacokinetic models as a example.
        

Cp=D/Vd exp(-ket)                                  (1)
CP=FDka/Vd/(ka-ke) (exp(-ket)-exp(-kat))/(ka-ke)   (2)


where Cp:plasma concentration, D:dose, Vd:volume of distribution, ke:elimination rate constant, 
F:absorption ratio,  ka:absorption, rate constant, t:time.
Eq.(1) is one-compartment model following rapid intravenous administration, and Eq.(2) is 
one-compartment model with first-order absorption.  The parameters Vd and ke are assumed to be 
common in Eqs.(1) and (2). Eqs.(1) and (2) mean that two time courses of the same drug are measured
 following rapid intravenous and oral (or intra-muscular) administration.  These models can be
 defined in the subroutine FFUNC in SULTIF.FOR or MULTIF.FOR as follows.

        
       SUBROUTINE FFUNC(J,P,T,CP)
       DIMENSION P(20)
       GOTO(100,200,300,400,500)J
   100 CP=P(1)/R(1)*EXP(-P(2)*T)
       RETURN
   200 CP=P(1)*P(3)*P(4)/(P(3)-P(2))*(EXP(-P(2)*T)-EXP(-P(3)*T))
       RETURN
   300 CONTINUE
       RETURN
   400 CONTINUE
       RETURN
   500 CONTINUE
       RETURN
       END
        

where P(1):Vd, P(2):ke, P(3):ka, P(4):F and R(1):D.
It should be noted that P(1), P(2).. are parameters to estimate, R(1),R(2).. are constans.  
Lines 100 and 200 are the equations defined by the user.
We construct the differential equations corresponding to the models represented by Eqs.(1) and (2).

        
dC1/dt = -ke C1
dC2/dt = -ke C2 + DFka/Vd exp(-kat)


and the initial conditions are


C1(0)=D/Vd
C2(0)=0
        

These differential equation can commonly be defined in the subroutine DIFF of SRUNGEF.FOR and MRUNGEF.FOR.

        
      SUBROUTINE DIFF(DC,C,T,P)
      DIMENSION C0(15),DC(15),C(15),P(20)
      COMMON NN,ND,DT,TT(2000),CC(2000,15),R(20)
      DC(1)=-P(2)*C(1)
      DC(2)=-P(2)*C(2)+R(1)P(3)*P(4)/P(1)*EXP(-P(4)*T)
      END


Initial conditions can be defined in the subroutine INI of SRUNGE.FOR and MRUNGE.FOR.
 

      SUBROUTINE INI(C0,P)
      DIMENSION C0(15),P(20)
      COMMON NN,ND,DT,TT(2000),CC(2000,15),R(20)
      C0(1)=R(1)/P(1)
      C0(2)=0.
      END
        

where DC(i) means the differential of C(i) with respect to t and C0(i) means the initial value. 
P(1):Vd, P(2):ke, P(3):F,  P(4):ka, R(1):D.
Transforming  Eqs.(1) and (2) into the Laplace domain, we  get the following equations.

        
CC1(s)=D/Vd/(s+ke)               (3)
CC2(s)=FDka/Vd/(s+ke)/(s+ka)     (4)
 

where CC1(s) and CC2(s) are the Laplace transforms of C1(t) and C2(t), respectively, s is the 
Laplace variable with respect to t.
It should be noted that the parameters to estimate are real, whereas CC1(s), CC2(s) and s are 
complex.  These equations can be defined in the subroutine FFUNC of  SFILTF.FOR and MFILTF.FOR as follows.

        
      SUBROUTINE FFUNC(J,T,P,S,CCP)
      DIMENSION P(20)
      COMMON R(20)
      COMPLEX S,CCP
      GOTO(100,200,300,400,500)J
  100 CCP=P(1)/(S+P(2))
      RETURN
  200 CCP= P(3)*P(4)*R(1)/P(1)/(S+P(2))/(S+P(4))
      RETURN
  300 CONTINUE
      RETURN
  400 CONTINUE
      RETURN
  500 CONTINUE
      RETURN
      END
        

The  new version of MFILTF.FOR becomes more compact  and  several times faster than the old 
version.
        
[4] Example Runs (SULTI and MULTI)
    The followings are example runs of SULTIF(or SRUNGEF and SFILTF).

 
SULTIF
###########################
# SIMULATION BY EQUATIONS #
###########################
NUMBER OF EQUATIONS ?
2
NUMBER OF PARAMETERS ?
4
P( 1) ?
100
P( 2) ?
0.5
P( 3) ?
2.0
P( 4) ?
1
INITIAL, AND FINAL TIME ?
0,8
TITLE OF T-AXIS ?
TIME(HR)
TITLE OF CP-AXIS ?
CP(UG/ML)
   CP(UG/ML) 
    100.000+1                                                  
           I                                                   
           I  1                                                
           I                                                   
     75.000+                                                   
           I    1                                              
           I    2 2                                            
           I        2 2                                        
     50.000+        1                                          
           I  2       1  2                                     
           I             1 2 2                                 
           I               1   2                               
     25.000+                 1 1 2                             
           I                     1 2 2                         
           I                         1 2 2 2                   
           I                               1 2 2  2 2 2 2 2    
       .000+2                                           1 1 2 2
           ++---------+---------+---------+---------+---------+
            .00      1.60      3.20      4.80      6.40      8.00
                                  TIME(HR)  
DO YOU CONTINUE(Y/N)?
N
Stop - Program terminated.


The followings are example runs of MULTIF(or MRUNGEF and MFILTF).
Table 1 presents two time courses of cefotiam, an antibiotics, following intravenous dose(50mg/kg)
 and intra-muscular dose(50mg/kg)

Table 1
----------------------------------------- 
time(min)    Cp^iv(ug/ml)   Cp^im(ug/ml) 
----------------------------------------- 
    5          102            6.03
   10           56.2          9.37
   15           24.8          13.4
   20           22.0          11.8
   40            8.04         10.8
   60            3.10          8.97
   90            --            6.31
---------------------------------------


MULTIF
#######################################
#          MULTI LINES FITTING        #
#######################################
E)XIT, I)NPUT, R)EAD, W)RITE, S)TOP
I
NUMBER OF LINES(1-5)?
2
NUMBER OF POINTS(1)?
6
T 1( 1)  , CP 1( 1) ?
5,102
T 1( 2)  , CP 1( 2) ?
10,56
T 1( 3)  , CP 1( 3) ?
15,24.8
T 1( 4)  , CP 1( 4) ?
20,22
T 1( 5)  , CP 1( 5) ?
40,8.04
T 1( 6)  , CP 1( 6) ?
60,3.1
NUMBER OF POINTS ( 2) ?
7
T 2( 1)  , CP 2( 1) ?
5,6.03
T 2( 2)  , CP 2( 2) ?
10,9.73
T 2( 3)  , CP 2( 3) ?
15,13.4
T 2( 4)  , CP 2( 4) ?
20,11.8
T 2( 5)  , CP 2( 5) ?
40,10.8
T 2( 6)  , CP 2( 6) ?
60,8.97
T 2( 7)  , CP 2( 7) ?
90,6.31
E)XIT, I)NPUT, R)EAD, W)RITE, S)TOP
W
FILE NAME TO WRITE ?
CEFO.IN
E)XIT, I)NPUT, R)EAD, W)RITE, S)TOP
E
(0)GAUSS-NEWTON         METHOD
(1)DAMPING GAUSS NEWTON METHOD
(2)MODIFIED MARQUARDT   METHOD
(3)SIMPLEX              METHOD
WHICH ALGORITHM DO YOU SELECT ?
1
SUBJECT NAME ?
TEST
NUMBER OF LINES (1-5) ?
2
NUMBER OF PARAMETERS ?
4
NUMBER OF POINTS ( 1) ?
6
WEIGHT OF DATA(0,1,2) ?
0
--- CONSTRAINT ON P(I) ---
(1) NO CONSTRAINT                        
(2) P(I)=Q(I)*Q(I)                       
(3) P(I)=B+(A-B)*SIN^2(Q(I))             
(4) P(I)=B+(A-B)*(EXP(Q(I))/(1+EXP(Q(I)))
WHICH CONSTRAINT DO YOU SELECT(1,2,3,4) ?
1
NUMBER OF CONSTANTS ?
1
R(1) ?
50
NUMBER OF PARAMETERS ?
4
INITIAL P( 1) ?
200.0
INITIAL P( 2) ?
0.1
INITIAL P( 3) ?
0.01
INITIAL P( 4) ?
1.0
INITIAL SS=  .11770E+04
LOOP=   1
DAMP= 1
P( 1)=  179.761400
P( 2)=     .113317
P( 3)=     .010926
P( 4)=     .929711
SS=  .12938E+03
LOOP=   2
DAMP= 1
P( 1)=  183.394800
P( 2)=     .118916
P( 3)=     .011166
P( 4)=     .915666
SS=  .11615E+03
LOOP=   3
DAMP= 1
P( 1)=  184.405700
P( 2)=     .119696
P( 3)=     .011126
P( 4)=     .917513
SS=  .11603E+03
**** CEFO.IN                 ****
BY DAMPING GAUSS NEWTON METHOD
WEIGHT =1/CP^(0)
CONSTRAINT ON P(I): NO CONSTRAINT                        
LOOP =    4
AIC =    69.800030
FILAL P( 1)=  184.466400   S.D.=   12.125880
FILAL P( 2)=     .119741   S.D.=     .008441
FILAL P( 3)=     .011121   S.D.=     .006667
FILAL P( 4)=     .917819   S.D.=     .433933
SS=  .11603E+03
DO YOU CONTINUE(Y/N) ?
Y
## TIME COURSE 1 ##
T 1=    5.000000 ,CP 1 =   101.368300(  102.000000)
T 1=   10.000000 ,CP 1 =    55.704090(   56.000000)
T 1=   15.000000 ,CP 1 =    30.610620(   24.800000)
T 1=   20.000000 ,CP 1 =    16.821200(   22.000000)
T 1=   40.000000 ,CP 1 =     1.533900(    8.040000)
T 1=   60.000000 ,CP 1 =      .139874(    3.100000)
MINIMUM,MAXIMUM OF T-AXIS ?
0,60
NAME OF T-AXIS ?
TIME
NAME OF CP-AXIS ?
CP(IV)
      CP(IV)  
    184.466+*                                                  
           I                                                   
           I                                                   
           I                                                   
    138.350+  *                                                
           I                                                   
           I   *                                               
           I    O                                              
     92.233+     *                                             
           I      *                                            
           I       *                                           
           I        O                                          
     46.117+         **                                        
           I           **                                      
           I             O***O                                 
           I                 *******         O                 
       .000+                        * ************************O
           ++---------+---------+---------+---------+---------+
            .00     12.00     24.00     36.00     48.00     60.00
                                  TIME      
SAVE(Y/N) ?
Y
FILE NAME ?
CEFO1.OUT
## TIME COURSE 2 ##
T 2=    5.000000 ,CP 2 =     6.871439(    6.030000)
T 2=   10.000000 ,CP 2 =    10.275770(    9.730000)
T 2=   15.000000 ,CP 2 =    11.794960(   13.400000)
T 2=   20.000000 ,CP 2 =    12.297230(   11.800000)
T 2=   40.000000 ,CP 2 =    10.966180(   10.800000)
T 2=   60.000000 ,CP 2 =     8.881471(    8.970000)
T 2=   90.000000 ,CP 2 =     6.370930(    6.310000)
MINIMUM,MAXIMUM OF T-AXIS ?
0,90
NAME OF T-AXIS ?
TIME
NAME OF CP-AXIS ?
CP(IM)
    CP(IM)    
      13.400+        O                                          
            I          *****                                    
            I       *** O   ******                              
            I      *              *O**                          
      10.050+     *O                   ****                     
            I                              ***O*                
            I    *                              *****           
            I   *                                    *****      
       6.700+                                             *****O
            I  *O                                               
            I                                                   
            I                                                   
       3.350+ *                                                 
            I                                                   
            I                                                   
            I                                                   
        .000+*                                                  
            ++---------+---------+---------+---------+---------+
             .00     18.00     36.00     54.00     72.00     90.00
                                  TIME      
SAVE(Y/N) ?
Y
FILE NAME ?
CEFO2.OUT
(0) GAUSS-NEWTON                METHOD
(1) DAMPING GAUSS NEWTON METHOD METHOD
(2) MODIFIED MARQUART           METHOD
(3) SIMPLEX                     METHOD
(-1)END
WHICH ALGORITHM DO YOU SELECT (0-3, OR -1) ?
-1
Stop - Program terminated.


The content of CEFO.IN is
-----------------------------
             2
             6
     5.00000,   102.00000
    10.00000,    56.00000
    15.00000,    24.80000
    20.00000,    22.00000
    40.00000,     8.04000
    60.00000,     3.10000 
             7
     5.00000,     6.03000
    10.00000,     9.73000
    15.00000,    13.40000
    20.00000,    11.80000
    40.00000,    10.80000
    60.00000,     8.97000
    90.00000,     6.31000
    90.00000,     6.31000
-----------------------------

The content of cefo1.out is
----------------------------
             6
     5.00000,   102.00000
    10.00000,    56.00000
    15.00000,    24.80000
    20.00000,    22.00000
    40.00000,     8.04000
    60.00000,     3.10000
            50
      .00000,   219.97080
     1.22449,   184.18580
     2.44898,   154.22230
     3.67347,   129.13330
     4.89796,   108.12580
     6.12245,    90.53586
     7.34694,    75.80742
     8.57143,    63.47502
     9.79592,    53.14886
    11.02041,    44.50257
    12.24490,    37.26286
    13.46939,    31.20092
    14.69388,    26.12513
    15.91837,    21.87507
    17.14286,    18.31643
    18.36735,    15.33670
    19.59184,    12.84171
    20.81633,    10.75261
    22.04082,     9.00337
    23.26531,     7.53870
    24.48980,     6.31230
    25.71429,     5.28541
    26.93878,     4.42557
    28.16327,     3.70562
    29.38776,     3.10279
    30.61225,     2.59802
    31.83673,     2.17537
    33.06123,     1.82148
    34.28571,     1.52516
    35.51020,     1.27705
    36.73470,     1.06930
    37.95918,      .89534
    39.18367,      .74969
    40.40816,      .62773
    41.63265,      .52561
    42.85714,      .44010
    44.08163,      .36851
    45.30612,      .30856
    46.53061,      .25836
    47.75510,      .21633
    48.97959,      .18114
    50.20408,      .15167
    51.42857,      .12700
    52.65306,      .10634
    53.87755,      .08904
    55.10204,      .07455
    56.32653,      .06242
    57.55102,      .05227
    58.77551,      .04377
    60.00000,      .03665
P( 1)=,      .22730
P( 2)=,      .14500
P( 3)=,      .48726
P( 4)=,      .00470
R( 1)=,    50.00000
AIC=    90.56375
---------------------------------
    

The content of cefo2.out
----------------------------
             7
     5.00000,     6.03000
    10.00000,     9.73000
    15.00000,    13.40000
    20.00000,    11.80000
    40.00000,    10.80000
    60.00000,     8.97000
    90.00000,     6.31000
            50
      .00000,      .00000
     1.83673,     2.96389
     3.67347,     5.29116
     5.51020,     7.10872
     7.34694,     8.51832
     9.18367,     9.60155
    11.02041,    10.42385
    12.85714,    11.03766
    14.69388,    11.48503
    16.53061,    11.79968
    18.36735,    12.00863
    20.20408,    12.13353
    22.04082,    12.19175
    23.87755,    12.19720
    25.71428,    12.16101
    27.55102,    12.09211
    29.38775,    11.99766
    31.22449,    11.88335
    33.06122,    11.75378
    34.89796,    11.61260
    36.73469,    11.46274
    38.57143,    11.30652
    40.40816,    11.14582
    42.24490,    10.98210
    44.08163,    10.81657
    45.91837,    10.65013
    47.75510,    10.48355
    49.59184,    10.31740
    51.42857,    10.15214
    53.26530,     9.98814
    55.10204,     9.82567
    56.93877,     9.66495
    58.77551,     9.50614
    60.61224,     9.34936
    62.44898,     9.19471
    64.28571,     9.04224
    66.12244,     8.89201
    67.95918,     8.74404
    69.79591,     8.59834
    71.63265,     8.45491
    73.46938,     8.31375
    75.30612,     8.17485
    77.14285,     8.03819
    78.97959,     7.90376
    80.81632,     7.77152
    82.65306,     7.64146
    84.48979,     7.51353
    86.32653,     7.38773
    88.16326,     7.26400
    90.00000,     7.14234
P( 1)=,      .27026
P( 2)=,      .12015
P( 3)=,     1.06559
P( 4)=,      .00920
R( 1)=,    50.00000
AIC=    69.92988
----------------------------------

------------------------------------------------------------------     
[5] Fast Inverse Laplace Transform (FILT)
    In the pharmacokinetic field, the Laplace transform  has been 
frequently  used for the linear models.  The  compartment  models 
described  by  simultaneous ordinary differential  equations  are 
usually solved by the Laplace transform method excepting the most 
simple  models.   The  dispersion  model  described  by   partial 
differential equations are also solved by the Laplace  transform.  
Though  the  transformed equation is ordinarily  obtained  in  the 
routine  work,  the inverse Laplace transform  from  the  Laplace 
domain  to  the  time domain is usually  difficult.    Even  if  the 
inverse transform is possible, the solution in the time domain is 
often given in the form of a complicated integral or an  infinite 
series which is difficult to handle. 
    The  fast  inverse Laplace transform (FILT) is  an  algorithm 
which  numerically  generates  the time  course  curve  from  the 
Laplace-transformed  model equation.  FILT was first proposed  by 
Hosono in the field of the wave optics [T.Hosono, Radio  Science, 
16,1015 (1981)].  MULTI(FILT) is a nonlinear regression  analysis 
program  where  FILT is combined with the least  squares  program 
MULTI [K.Yamaoka et al., J.Pharm.Dyn., 4, 879(1981)]. 
[6] INSTALLING AND RUNNING MULTI(FILT) 
    MULTI(FILT)  is  written  in  FORTRAN77  that  supports   the 
arithmetic of complex numbers.  Therefore, FORTRAN80 on CP/M  or 
the old version of MS-FORTRAN which is limited to the arithmetic 
of  real  numbers is improper to use MULTI(FILT).    The  minimum 
system requirements for MULTI(FILT) are 
1. 320 K bytes random access memories at least. 
2. Preferably two disk drives. 
3. MS-DOS(PC-DOS) version 2.0 or later. 
SFILT is a program for computer simulation using FILT  algorithm.  
The  details  about  SFILT is given  in  the  following  section. 
MULTI(FILT) and SFILT include a subroutine for the definition of 
model equations.  When using these programs, the user is supposed 
to  use  FORTRAN  system  diskette in A  drive  and  the  program 
diskette in B drive. 
[7] EXAMPLE RUN OF MULTI(FILT) 
    Initially,  the  user  must  define  the  Laplace-transformed 
pharmacokinetic  model in the subroutine FFUNC by means of an editor such
as WordStar  or WordMaster etc.  For example, the one-compartment  model 
is  given by the following equation. 
Cp(t) = D/Vd exp(-ket)             (5) 
where  Cp(t) is plasma concentration, D is dose, Vd is volume  of 
distribution, ke is elimination rate constant and t is time.  The 
Laplace-transformed equation corresponding to Eq.(1) is given by 
CCp(s) = D/Vd 1/(s+ke)              (6) 
where CCp(s) is transformed plasma concentration and s is  Laplace 
variable corresponding to t.  The defined equation in MEQ1.FOR is 
as follows. 

C  =====  FUNCTION FOR MULTI(FILT)  ====== 
      SUBROUTINE FFUNC(J,T,P,S,CCP) 
      DIMENSION P(20)
      COMPLEX S,CCP 
      GO TO (100,200,300,400,500) J 
  100 CCP=P(1)/(S+P(2)) 
      RETURN 
  200 CONTINUE 
      RETURN 
  300 CONTINUE 
      RETURN 
  400 CONTINUE 
      RETURN 
  500 CONTINUE 
      RETURN
      END 

where  CCP  is CCp(s), P(1) is D/Vd, P(2) is ke and  S  is 
Laplace variable.  It should be noted that CCP and S are  complex 
variables   and P(1) and P(2) are real variables.   When  another 
variable is adopted in the definition, the declaration of COMPLEX 
or REAL at the beginning of the subroutine is necessary.
    The following is an example to transform MEQ1.FOR to MEQ1.OBJ.
FOR1 B:MFILT,B:MFILT;
PAS2 
The next link of MULTI(FILT) and the execution are 
done  by  the  following commands,  assuming  that  MFILT.FOR  is 
 
previously compiled to MFILT.OBJ in the same way as mentioned above.
LINK B:MFILT,B:MFILT;
B:MFILT
The example RUN is as follows.

=======================================
           MULTI LINES FITTING         
   BY FAST INVERSE LAPLACE TRANSFORM   
=======================================
  (0)GAUSS NEWTON         METHOD
  (1)DAMPING GAUSS NEWTON METHOD
  (2)MODIFIED MARQUARDT   METHOD
  (3)SIMPLEX              METHOD
 
   WHICH ALGORITHM DO YOU SELECT ? 
1                          -- selection of algorithm 
   SUBJECT NAME ? 
ONE                         -- data name        
   NUMBER OF LINES (1-5)  ? 
1                           -- number of time courses
   WEIGHT OF DATA (0,1,2) ? 
0                           -- WT=CP^0 which consequently means 1
   NUMBER OF PARAMETERS   ? 
2                           -- P(1) and P(2)
   NUMBER OF POINTS ( 1)  ? 
4                           -- in this case, four points
   T 1( 1)  , CP 1( 1) ? 
1,10.1
   T 1( 2)  , CP 1( 2) ?    -- time course data 
2,5.2
   T 1( 3)  , CP 1( 3) ? 
3,2.3
   T 1( 4)  , CP 1( 4) ? 
4,1.2
 
 --- CONSTRAINT ON P(I) ---
  (1) NO CONSTRAINT                        
  (2) P(I)=Q(I)*Q(I)                       
  (3) P(I)=B+(A-B)*SIN^2(Q(I))             
  (4) P(I)=B+(A-B)*(EXP(Q(I))/(1+EXP(Q(I)))
  WHICH CONSTRAINT DO YOU SELECT (1,2,3,4)  ? 
1                          -- no constraint is selected   
  INITIAL P( 1) = ? 
20
  INITIAL P( 2) = ? 
0.7
  INITIAL SS=  .12271E+00
 
LOOP=   1
 DAMP= 1
 P( 1)=   20.747670
 P( 2)=     .710940
 SS=  .71331E-01
 
LOOP=   2
  DAMP= 1
 P( 1)=   20.377990
 P( 2)=     .701354
 SS=  .70923E-01
 
LOOP=   3
  DAMP= 1
 P( 1)=   20.572140
 P( 2)=     .706740
 SS=  .69009E-01
 
LOOP=   4
  DAMP= 1
 P( 1)=   20.532370
 P( 2)=     .705109
 SS=  .68991E-01
 
LOOP=   5
  DAMP= 1
 P( 1)=   20.510900
 P( 2)=     .704429
 SS=  .68951E-01
 
LOOP=   6
  DAMP= 3
 P( 1)=   20.551640
 P( 2)=     .705669
 SS=  .68909E-01
 
 **** ONE                  ****
 BY DAMPING GAUSS NEWTON METHOD
  WEIGHT =1/CP (0)
   - CONSTRAINT ON P(I) -
     NO CONSTRAINT                        
  LOOP =    7
  AIC =    -6.700156            -- Akaike's information criterion
  DP FOR JACOBIAN=     .001000
 FINAL P( 1)=   20.553290   S.D.=     .667971
 FINAL P( 2)=     .704723   S.D.=     .022700
 SS=  .68904E-01
 
 
 DO YOU CONTINUE (Y/N) ? 
Y
 
T 1=    1.000000 ,CP 1 =    10.148260(   10.100000)
T 1=    2.000000 ,CP 1 =     5.010694(    5.200000)
T 1=    3.000000 ,CP 1 =     2.473990(    2.300000)
T 1=    4.000000 ,CP 1 =     1.221586(    1.200000)
 
MINIMUM, AND MAXIMUM OF T-AXIS  ? 
0,4
MINIMUM, AND MAXIMUM OF CP-AXIS ? 
0,20
  NAME OF  T-AXIS ? 
TIME
  NAME OF CP-AXIS ? 
CONC
     CONC      
      20.000+**                                                
            I  *                                               
            I   *                                              
            I    *                                             
      15.000+     **                                           
            I       *                                          
            I        **                                        
            I          **                                      
      10.000+            O*                                    
            I              ***                                 
            I                 ***                              
            I                    ***                           
       5.000+                       **O*                       
            I                           ******                 
            I                                 ****O****        
            I                                          *******O
        .000+                                                  
            ++---------+---------+---------+---------+---------+
             .00       .80      1.60      2.40      3.20      4.00
                                       TIME      
 
 DO YOU CONTINUE (Y/N) ? 
N
Stop - Program terminated.

The O and * in the plots are the observed and theoretical values, 
respectively.   In  this  example,  the  single  time  course  is 
assumed.   When  several time courses are  available,  the  model 
equations may be defined at lines 100, 110, 120, ... in the subroutine
FFUNC. 
[8] EXAMPLE RUN OF SFILT
    SFILT  is a simulation program using FILT algorithm.  To  RUN 
SFILT,  the Laplce-transformed equations must be defined  in  the 
subroutine FFUNC.   Let's pick up the following model.

Cp = Co exp(-ket)                              (7)
Cp = CokaF/(ka-ke) (exp(-ket) - exp(-kat))     (8)

where  Co(=D/Vd) is  initial  concentration,  ke  is  elimination 
constant,  F  is extent of availability, ka  is  absorption  rate 
constant.   Eqs(7)  and  (8)  are  one-compartment  models   with 
intravenous(IV) and oral(PO) administrations, respectively.   The 
definition for this model in subroutine FFUNC is as follows. 

      SUBROUTINE FFUNC(J,T,P,S,CCP) 
      DIMENSION P(20)
      COMPLEX S,CCP 
      GO TO (100,200,300,400,500) J 
  100 CCP=P(1)/(S+P(2)) 
      RETURN 
  200 CONTINUE 
      RETURN 
  300 CONTINUE 
      RETURN 
  400 CONTINUE 
      RETURN 
  500 CONTINUE 
      RETURN
      END 

where  P(1)  is Co, P(2) is F, P(3) is ke and P(4)  is  ka.   The 
coefficient  10  at  line  200 is the factor  to  adjust  the  PO 
concentration to IV concentration.  COMPILE, LINK and RUN of SFILT
are done in the same way as MULTI(FILT) 
FOR1 SFILT,SFILT;
PAS2
LINK SFILT,SFILT;
SFILT
The following is the example execution of SFILT.

 =====  FUNCTION SIMULATION PROGRAM  =====
 =====                 BY FILT       =====
 NUMBER OF EQUATION   ? 
2
 NUMBER OF PARAMETERS ? 
4
P( 1)=
180
P( 2)=
0.918
P( 3)=
7.18
P( 4)=
0.667
 MINMUM, AND MAXIMUM OF T-AXIS ? 
0,1.5
 MINIMUM, AND MAXIMUM  OF Y-AXIS ? 
0,179
 NAME OF X-AXIS ? 
TIME
 NAME OF Y-AXIS ? 
CP
           CP             
     179.987+1                                                 
            I                                                  
            I                                                  
            I 1                                                
     134.990+                                                  
            I          222222                                  
            I  1    222      222222                            
            I     22               222222                      
      89.993+   12                       2222222               
            I    1                              2222222        
            I   2                                      22222222
            I  2  1                                            
      44.997+      1                                           
            I 2     11                                         
            I         11                                       
            I           11111                                  
        .000+                1111111111111111111111111111111111
            ++---------+---------+---------+---------+---------+
             .00       .30       .60       .90      1.20      1.50
                                  TIME           
 
 DO YOU CONTINUE (Y/N)  ? 
N
Stop - Program terminated.

The numbers 1 and 2 in plots correspond to the IV and PO concentrations, 
respectively.

[REFERENCES]
1. Y.Yano, K.Yamaoka and H.Tanaka: A Nonlinear Least Squares Program,
   MULTI(FILT), Based on Fast Inverse Laplace Transform for Micro-
  computers, Chem.Pharm.Bull. 37, 1035-1038 (1989).
2. K.Yamaoka, T.Nakagawa and T.Uno: Application of Akaike's 
   Information Criterion (AIC) in the Evaluation of Linar Pharmacokinetic
   Equations, J.Pharmacok.Biopharm., 6, 165-175 (1978).
3. K.Yamaoka et al.: An Analysis Program MULTI(ELS) Based on Extended
   Nonlinear Least Squares Method for Microcomputers, J.Pharmacobio-
   Dyn., 9, 161-173(1986).
4. K.Yamaoka and H.Tanaka: A New Version of MULTI(ELS) for Extended
   Nonlinear Least Squares, J.Pharmacobio-Dyn., 10, 26-34 (1987).
5. K.Yamaoka and T.Nakagawa: A Nonlinear Least Squares Program Based 
   on  Differential Equations, MULTI(RUNGE), for Microcomputers,
   J.Pharm.Dyn., 6, 595-606 (1983).
6. K.Yamaoka, Y.Tanigawara, T.Nakagawa and T.Uno: A Pharmacokinetic
   Analysis Program (MULTI) for Microcomputer, J.Pharm.Dyn., 4,
   879-885(1981).
7. K.Yamaoka, et al.: A Nonlinear Multiple Regression Program,
   MULTI2(BAYES), Based on Bayesian Algorithm for Microcomputers,
   J.Pharmacobio-Dyn., 8, 246-256(1985).