Using Limited Concentration Data for the Determination of Rate

annealing and tabu search, that has the ability to locate a minimum of a function near enough to the global minimum to yield meaningful results withou...
0 downloads 0 Views 98KB Size
Download PDF
Environ. Sci. Technol. 1998, 32, 3207-3212

Using Limited Concentration Data for the Determination of Rate Constants with the Genetic Algorithm KRISTA B. VON ARX, JOHN J. MANOCK, SCOTT W. HUFFMAN, AND MICHAEL MESSINA* Department of Chemistry, University of North Carolina at Wilmington, Wilmington, North Carolina 28403

A computational strategy is presented that uses limited concentration vs time data to iteratively determine rate constants for a postulated reaction mechanism. These rate constants, once found, can be used to compute complete concentration profiles as a function of time and can be transferred to other reaction mechanisms that contain common reaction steps. The computational algorithm presented here maps this iterative rate constant search onto a problem of functional optimization, where the values of the rate constants themselves are optimized in order to minimize a well-defined error function. This error function contains multiple minima, which obviates the use of conventional computational methods of functional optimization such as conjugate gradient and Simplex methods. Alternatively, the genetic algorithm is found to be extremely robust in its ability to find solutions near enough to the global minimum to yield meaningful results without getting “stuck” in local minima. The algorithm presented here uses the Genetic Algorithm to get near a minimum that is close to a useful local minimum and then uses the Simplex method to polish the solutions, i.e., to find this minimum. The usefulness of the algorithm is demonstrated on two model problems each based on an environmentally relevant chemical reaction system.

I. Introduction Chemical reactions that occur in the environment usually have complex reaction mechanisms involving multiple reaction steps that span long periods of time outside the laboratory where quantitative measurements of complete reaction profiles are prohibitive. Further, these reactions do not easily lend themselves to the controllable kinetic measurements often available in the laboratory, thus one usually has incomplete information about the reaction kinetics. This incomplete information can either be manifested as a small subset of the absolute rate constants for the aggregate reaction mechanism and/or sparsely measured concentration vs time data for the reactants and products. Thus important pieces of information are missing, these being the complete set of absolute rate constants for each step in the reaction mechanism, and the complete reaction profiles. A computational strategy for determining these rate constants is to define an error function that is the difference * Corresponding author e-mail: [email protected]; fax: (910)962-3013; phone: (910)962-3450. S0013-936X(97)00947-4 CCC: $15.00 Published on Web 09/03/1998

 1998 American Chemical Society

between the experimentally measured concentration data and computed concentrations that arise from a set of guessed rate constants and then find the best set of rate constants by minimizing the error function. In this work we present such a computational strategy by coupling a heuristic search algorithm that gets near a useful local minimum in parameter space and then uses the Simplex method (1) to find this useful local minimum using as initial guesses for the parameters the output from the heuristic algorithm. The Genetic Algorithm (GA) (2) distinguishes itself as one of several heuristic optimization methods, such as simulated annealing and tabu search, that has the ability to locate a minimum of a function near enough to the global minimum to yield meaningful results without getting stuck in unphysical local minima. The GA has been used in systems as varied as the conformational analysis of large biomolecules (3), prediction of protein folding (4), predicting structures of binary metallic alloys (5), and the fitting of molecular potential energy surfaces (6). There are many sources in the literature that describe the inner workings of the GA, and we refer the interested reader to these sources (7-12). The rate constant search algorithm is tested on two model systems. Each of these systems is based on the reaction mechanism for the removal of chlorpyrifos from the aqueous phase. These model systems provide hard tests of the search algorithm since the rate constants and the concentrations of the species involved in these two systems each vary over 4 orders of magnitude. We will show that the rate constants involved in these two model systems are accurately predicted using limited concentration data and assuming no initial knowledge of the rate constants.

II. Theory Given a general reaction mechanism containing a series of steps, we define two vectors: c˜(t), the concentration of all the species at any time and k˜ , the set of rate constants governing the elementary reaction steps. The ith vector component of c˜(t), ci(t), represents the concentration of species i at time, t. The kinetic equation that governs the time evolution of the concentration of species i is defined as

c3 i(t) ) fi[c˜(t),k˜ ,t]

(1)

where fi is a function dependent upon the instantaneous concentrations, c˜(t), the rate constants, k˜ , and time, t. The solution of the set c3 (t) yields the complete reaction profile of c˜(t) for a given set of rate constants, k˜ . Specifically, the concentration of species i at time tj is ci(tj). Let us assume that we have an incomplete set of concentration vs time data which we term the target concentrations. For species i, there may exist target concentration data dijat several discrete times, tij. The error function for species i is defined as

i(k˜ ) )

∑ j

({ } 1

dij

)

[ci(tj) - dij]

2

(2)

The error function for species i is the squared difference between the computed concentration, ci(tj), for a guessed set of rate constants, k˜ , and the target concentrations, dij at the same times, tj. We normalize each i(k˜ ) by the target concentration dij, since the concentrations of species involved in the reaction mechanism can vary over several orders of magnitude. Without the normalization factor i(k˜ )’s for species with large concentrations will unnaturally be weighted more heavily. The total error function for the complete VOL. 32, NO. 20, 1998 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

3207

volume of the water; F is the mass of the fish; and S is the mass of the soil. The values of d, V, F, and S are 88.8 cm, 3.59 × 108 mL, 8580 g, and 1.49 × 107 g, respectively (13). The second column in Table 1 shows the postulated values of these rate constants according to Neely and Blau (13) for an initial aqueous chlorpyrifos concentration of 5.75 ppb. We now consider another model, model II, somewhat based on the Neely and Blau reaction mechanism described above. The reason for introducing a second model is to have a test of the rate constant search algorithm where the elementary reactions are not all first-order in all the species concentrations. This will demonstrate that the rate constant search algorithm presented here is flexible and can be used to find rate constants in a variety of mechanisms. Model II is just a variation of model I where we redefine the rate of change of chlorpyrifos concentration in the water phase as FIGURE 1. The mechanism for the removal of chlorpyrifos from the aqueous phase. Here, Cw, Cs, and Cf are the concentrations of chlorpyrifos in water, soil, and fish, respectively. C, B, and I are the chlorpyrifos vapor, the hydrolysis product of chlorpyrifos, and the chlorpyrifos metabolite from fish, respectively. The rate constants are represented by k1 to k6. reaction then becomes

(k˜ ) )

∑ (k˜ ) i

(3)

i

The iterative rate constant search algorithm now proceeds as follows: An initial guess of the set of rates, k˜ , is made. Solving the kinetic eq 1 with this set of rate constants gives c˜(t), the set of concentrations as a function of time. The error function is then formed from eq 3. Minimizing this error function produces a new set of rate constants, k˜ , with which the process is repeated until the value of the error function approaches zero within some predefined tolerance. The minimization of the error function is the critical step here and a mechanism is needed to do this. This function is first minimized via the GA since the GA is quite efficient at getting near to useful local minima but not so efficient at finding the precise location of these minima (2). We then use the GA results as initial guesses in the Simplex algorithm which is efficient at finding the precise location of the useful local minimum from the set of guesses that are supplied by the GA.

III. Results and Discussion A. The Model Systems. In model I we consider a previously worked out reaction mechanism for the removal of chlorpyrifos from the aqueous phase. This is a complex mechanism that involves six rate constants and we show this mechanism in Figure 1 (13, 14). Since the only concentrations we are concerned with are those of chlorpyrifos in the aqueous phase, in soil, and in fish, the system can be completely described by the following three differential equations (13):

dCw k1Cw k3FCw k5SCw k6SCs )- k2Cw + (4) dt d V V V dCs ) k5Cw - k6Cs dt

(5)

dCf ) k3Cw - k4Cf dt

(6)

where Cw, Cs, and Cf are the concentrations of chlorpyrifos in water, soil, and fish, respectively. The set k˜ ) {k1, k2, ..., k6} is the set of rate constants for chlorpyrifos’ transfer between the phases; d is the depth of the water; V is the 3208

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 32, NO. 20, 1998

k1C2w k3FCw k5SC2w k6SCs dCw )- k2Cw + (7) dt d V V V The rate of change of chlorpyrifos concentration in the soil and fish compartments are unchanged from model I, while the rate of change of chlorpyrifos concentration in the water compartment is second order in some of the species. We will keep the same values of Neely and Blau for the rate constant set k˜ ) {k1, k2, ..., k6}, and the initial aqueous chlorpyrifos concentration. We have found that model II, although seemingly only modified slightly from model I, gives quite different concentration profiles for all the species. B. Computational Details. The differential equations in eqs 4-7 are solved numerically via a fourth-order RungeKutta adaptive step-size integrator (1). The differential equation solver is used to generate exact concentration vs time profiles for all species, and as a subroutine in the GA (12) and Simplex algorithm (1) to evaluate the error function. The error function is computed in the GA and Simplex algorithm by calling the differential equation solver as a subroutine. The GA algorithm was run with the following tuning parameters: a population size of 10, probability for mutation ) 0.1, probability for crossover ) 0.5. The GA was also run with uniform crossover and with creep mutations enabled. The search boundaries for each of the parameters were initially set at a minimum of zero, and a maximum of unity. As the program ran, everytime a search boundary was exceeded for a particular parameter, the search boundary for this parameter was increased by 10, until the GA yielded values of the parameter within the search boundary. The tolerance for optimization for the GA was 1 × 10-3, and whenever the error function was less than or equal to this tolerance the GA algorithm was stopped. It was found that it typically took on the order of 5000 GA iterations to converge to user tolerance. This took approximately 5 min of CPU time on an SGI Indigo2. When the tolerance for optimization was reached by the GA the best values for each of the six rate constants as determined by the GA were then used as initial guesses in a downhill Simplex algorithm. We defined the error function in the Simplex algorithm as in the GA. The tolerance for optimization for the Simplex algorithm was 1 × 10-4 or 10 000 iterations, whichever came first. C. Results. i. Model I. The reaction profiles of Cw, Cs, and Cf for model I with rate constants from column 2 of Table 1 are shown as solid lines in panels a, b, and c of Figure 2, respectively. We now iteratively solve for the rate constants, treating them as unknowns. A total of nine target concentrations, shown in Table 2, are picked as three concentrations at three distinct times for each of the three phases of concern. Three times represent a reasonable experimental strategy to gather concentration data in order to fit the rate constants,

TABLE 1. Rate Constants for the Removal of Chlorpyrifos from the Aqueous Phase GA Model I k˜

k1 k2 k3 k4 k5 k6

proposed

1, 3, and 7 days

after simplex

1 day and 3 days

1 day and 7 days

3 days and 7 days

no water information

no soil information

no fish GA model II: information 1, 3, and 7 days

0.12 cm h-1 0.33 .28 0.38 0.029 0.53 0.28 0.13 0.91 9.6 × 10-3 h-1 7.3 × 10-3 8.0 × 10-3 7.9 × 10-3 9.7 × 10-3 4.3 × 10-3 9.0 × 10-3 9.8 × 10-3 6.7 × 10-3 55.5 mL gfish-1 h-1 47.1 54.3 59.7 58.1 45.9 59.9 58.9 55.5 0.078 h-1 0.063 .076 0.082 0.084 0.067 0.010 0.082 0.21 -1 -1 0.67 mL gsoil h 0.74 .67 0.85 0.88 0.88 0.75 0.74 0.69 -1 0.052 h 0.061 .055 0.078 0.070 0.063 0.051 0.055 0.053

.15 9.4 × 10-3 49.3 6.9 × 10-2 0.65 5.1 × 10-2

TABLE 2. Computed Target Concentrations from the Exact Solution of Rate Equations for Model I with Proposed Rate Constants of Neely and Blau dij (ppb) Cw Cs Cf

1 day

3 days

7 days

3.1 35.0 2240.0

1.9 28.9 1530.0

0.9 14.1 735.0

of Neely and Blau the target rate constants. We refer to the complete reaction profiles that arise from these target rate constants as target reaction profiles. There are no initial guesses for the rate constants and the GA chooses their initial values randomly within the search boundaries, thus we are assuming no a priori knowledge of the rate constants. The algorithm is initiated by solving the set of differential equations in eqs 4-6 with the randomly chosen initial values for rate constants k1 to k6. The error function is computed according to eq 3 as the squared difference between the target concentrations and the computed concentrations; the computed concentrations at any particular time being the solutions of the set of differential equations in eqs 4-6. This function is then supplied to the GA as the external function to be minimized. The GA minimizes this function by finding a new set of rate constants which are then used in the subsequent steps and the process is repeated until the value of (k˜ ) approaches zero within the defined tolerance. The optimized parameter values from the GA are then used as initial guesses in the Simplex algorithm where they are further polished, as described above.

FIGURE 2. The target reaction profiles (solid line) compared to the reaction profiles arising from rate constants found by the search algorithm (dashed line) for model I. Panels a, b, and c show the comparison for Cw, Cs, and Cf, respectively. a data point at short, intermediate, and long times. These target concentrations are taken from the exact solution, thus we are testing the algorithm by seeing if it can reproduce the rate constants given in column 2 of Table 1 using this limited concentration data and no initial knowledge of the rate constants themselves. We call these predicted rate constants

In column 3 of Table 1, we show the rate constants obtained from the GA portion of the rate search algorithm using the limited concentration data from Table 2 with no a priori knowledge assumed about their values. These rate constants yielded a value of the error function of 1 × 10-3. These rate constants were then used as initial guesses in the Simplex algorithm where they are further polished. The polished values of these rate constants are shown in column 4 of Table 1, and they yield a value for the error function of 1 × 10-4. Comparing columns 2, 3, and 4 in Table 1 shows that the Simplex algorithm significantly improves the values of the rate constants obtained from the GA alone. As can be seen by comparing columns 2 and 4, the GA + Simplex search algorithm does extremely well in reproducing the values of all six rate constants. In fact rate constants k2 through k6 are within 90% of the target values (the values postulated by Neely and Blau in column 2 in Table 1) (13). Only k1 is not well-predicted by the search algorithm, the algorithm giving a value of 0.28 cm h-1 compared to a target value of 0.12 cm h-1. The inability of the search algorithm to predict the rate constant k1 is not too surprising based on the structure of the differential equations that govern the reaction in eqs 4-6. Since the two rate constants k1 and k2 only appear in eq 4, and both multiply the concentration of chlorpyrifos in the aqueous phase, Cw, only the sum (k1/d) + k2 is important to the reaction mechanism. This fact is borne out in the VOL. 32, NO. 20, 1998 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

3209

calculations. According to Table 1 this sum is approximately equal to 1 × 10-2 for all example problems based on model I, and further when the algorithm overestimates the value of k1 the value of k2 is underestimated. It is also interesting to point out here that for model II the rate constants k1 and k2 are distinguishable since k1 represents a second-order process in Cw while k2 represents a first-order process. Thus the algorithm is able to predict both k1 and k2 individually to good accuracy for model II, as demonstrated by the last column in Table 1. Overall the rate search algorithm does well in reproducing the target rate constant values with no a priori knowledge of their true values with a limited set of concentration vs time data. We still must determine if the values obtained via the rate constant search are sufficiently accurate enough to predict concentrations of species at times for which there is no experimental data available. To do this we use the rate constants determined by the rate constant search algorithm, i.e., column 4 in Table 1, in the rate equations, eqs 4-6, to obtain the complete reaction profile. In Figure 2 we show the comparison of the target reaction profiles determined for Cw, Cs, and Cf from the target rate constants (column 2 in Table 1) and the reaction profiles of these species when using the rate constants as found by the rate constant search algorithm (column 4 in Table 1). Figure 2 shows that the rate constants that are iteratively determined from the search algorithm accurately reproduce the whole reaction profile of chlorpyrifos in all three phases. It is also important to point out that the rate constants found from the search algorithm yield concentration profiles that are nearly indistinguishable from the exact solutions at times that are twice as long as the longest target data time. The GA component of the algorithm was key in finding a parameter set that was accurate. For example we ran the Simplex code without the GA preoptimization starting with the same initial parameter values as that given by the GA. The Simplex algorithm yielded solutions that came nowhere near the tolerance and became stuck in a local minimum, yielding the following values for the six rate constants, (7.7 × 10-3, 9.8 × 10-3, -1.1 × 10-2, 9.6 × 10-3, 0.138, -6.0 × 10-3). These rate constants are nowhere near the correct values and several are actually negative. The final Simplex polishing step of the search algorithm, on the other hand, was not as necessary for finding a good set of values of the rate constants. Rate consants determined from the GA component of the search algorithm alone were close to the true rate constants. But, as demonstrated by columns 2-4 in Table 1, the Simplex polishing step significantly improves the rate constants obtained from the GA component of the search algorithm. Further, the Simplex step is extremely inexpensive computationally. In all cases we found the Simplex algorithm was able to polish the rate constants within a few CPU seconds. The Simplex step is much more efficient than the GA step once the values of the rate constants are close to their optimal values. In fact even after 10 000 GA iterations we were not able to obtain a value for the error function less than 1 × 10-3, while the Simplex algorithm is able to obtain a value for the error function of 1 × 10-4 using as initial guesses the rate constants obtained after 5000 GA iterations. We now offer evidence of the robustness of the rate search algorithm. We must make sure that the algorithm did not serendipitously get stuck in a useful local minimum solely because of a lucky choice of initial values of parameters in the GA. To show that this was not the case we offer three diagnostic tests of the search algorithm. In the first test we restarted the GA with different initial values for the parameters. In the second test, we explore the dependence of the solutions on the values of the user defined tolerance parameters in the GA by computing the rate constants for 3210

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 32, NO. 20, 1998

FIGURE 3. The target reaction profiles to the results of the search algorithm with only the concentrations from two times available. Panel a shows the search results for Cw with data from all three times (solid line), 1 day and 3 days (solid line with squares), 1 day and 7 days (dashed line with diamonds), and 3 days and 7 days (dashed line with ×). Panels b and d show the same results for Cs and Cf, respectively. several different values of the tolerance parameter. Finally, the third test involved the stability of the search algorithm when there is an artificially induced error in the target concentration data. For each test the search algorithm was shown to be stable and yield roughly the same values for the six rate constants as given in column 4 of Table 1. The interested reader is referred to the Supporting Information for the numerical results of these diagnostic tests. It might be expected that the rate constant search algorithm degrades in accuracy when less information is available about the concentration of the species. We test the limits of the rate constant search algorithm by eliminating certain parts of the available target data. Lack of information can come in two ways: (1) there can be a lack of time

FIGURE 4. The target reaction profiles compared to the reaction profiles from the search algorithm with only two species’ concentration data available. Panels a, b, and c refer to Cw, Cs, and Cf, respectively. In each panel, the solid line refers to the target reaction profiles; the solid line with squares to, no available water data, the dashed line with diamonds, to no soil data; and the dashed line with × marks, to no fish data. information, meaning that we have some concentration data about all the species, but only at few times along the reaction profile, and (2) there can be a lack of species information, meaning that there is no concentration data available for one or more of the species involved in the reaction. We will test the effect on the accuracy of the rate search algorithm with respect to both of these limits. First we test the limitation of the search algorithm with respect to the lack of time information. To do this we eliminate one set of time information from the target concentration set in Table 2. In columns 5-7 in Table 1, we show the rate constants predicted from the search algorithm when we have information about days 1 and 3 (column 5), days 1 and 7 (column 6), and days 3 and 7 (column 7). As evidenced by Table 1, the rate constant search algorithm still

FIGURE 5. The target reaction profiles (solid line) compared to the reaction profiles arising from rate constants found by the search algorithm (dashed line) for model II. Panels a, b, and c show the comparison for Cw, Cs, and Cf, respectively. does well in predicting the rate constants when there is a lack of time information about the concentration of each of the species. To further show this we plot the complete reaction profiles as predicted by these rate constants. In Figure 3, the panels a, b, and c represent Cw, Cs, and Cf, respectively. These results compare favorably to the target reaction profiles. Next we test the limitation of the search algorithm with respect to the lack of concentration information of some of the species. In columns 8-10 of Table 1, we show the rate constants predicted from the search algorithm when we exclude the target concentration information from Table 2 of water (column 8), soil (column 9), and fish (column 10). Although these rate constants do stray some, they still give acceptable results. Figure 4 shows the reaction profiles predicted by the rate constants from columns 8-10 of Table 1. Panels a, b, and c refer to Cw, Cs, and Cf, respectively. Cw VOL. 32, NO. 20, 1998 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

3211

TABLE 3. Computed Target Concentrations from the Exact Solution of Rate Equations for Model II with Proposed Rate Constants of Neely and Blau dij (ppb)

1 day

3 days

7 days

Cw Cs Cf

1.3 18.6 1140.6

0.9 12.93 681.1

.65 8.9 489.0

and Cf are reproduced quite well by the algorithm except when there is a lack of information about these species. When there is no information about Cw the algorithm does poorly in predicting the reaction profile of Cw, likewise for Cf. On the other hand, when there is no target concentration data for soil, Cs, the algorithm does extremely well in reproducing the reaction profiles of all three compartments including Cs, this is evidenced by the dashed lines connected by open diamonds in panels a-c in Figure 4. The reason the Cs concentration profile is reproduced accurately when no Cs concentration data is available can be determined by examining the rate equations in eqs 4-6. The differential equation for Cs depletion only depend on rate constants k5 and k6. The differential equation for Cw depends on both these rate constants, so if concentration data is available for Cw the algorithm has a means to find k5 and k6 and thus information for the whole of the Cs concentration profile. ii. Model II. We now test the rate constant search algorithm for model II as defined by the rate equations given by eqs 5-7. As mentioned above the reaction in model II is no longer first-order in all species concentrations. The reaction profiles of Cw, Cs, and Cf with the Neely and Blau target rate constants from Table 1 are shown as solid lines in panels a, b, and c of Figure 5, respectively. Comparing Figures 2 and 5 show that models I and II give quite different reaction profiles for all three species. We now iteratively solve for the rate constants, treating them as unknowns. A total of nine target concentrations, shown in Table 3, are picked as three concentrations at distinct times for each of the three phases of concern. These target concentrations are taken from the exact solution thus, again, we are testing the algorithm by seeing if it can reproduce the rate constants given in column 2 of Table 1 using this limited concentration data and no initial knowledge of the rate constants themselves. In column 11 of Table 1 we show the search algorithm predictions of the rate constants. Comparing columns 2 and 11 shows that the search algorithm is able to reproduce the target rate constants quite accurately for model II. We now compare the reaction profiles that is obtained from the search algorithm rate constants to the exact profiles. To do this we use the rate constants determined by the rate constant search algorithm, i.e., column 11 in Table 1, in the rate equations, eqs 5-7, to obtain the complete reaction profile. In Figure 5 we show the comparison of the target reaction profiles determined for Cw, Cs, and Cf from the target rate constants (column 2 in Table 1) and the reaction profiles of these species when using the rate constants as found by the rate constant search algorithm (column 11 in Table 1). Again, the rate constants that are iteratively determined from the search

3212

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 32, NO. 20, 1998

algorithm produce reaction profiles that are nearly indistinguishable from the exact profiles at times that are twice as long as the longest target data time. The rate constant search algorithm presented here is shown to be robust and generally applicable to problems of environmental interest. We believe that such an algorithm can be useful for finding unknown rate constants in complex reaction mechanisms where there is only a sparse set of experimental concentration data. It should be pointed out here that the algorithm must be supplied with the reaction mechanism for the reaction of interest. Further work is being done to generalize the algorithm to allow for determination of both unknown reaction mechanisms and rate constants by treating the exponents in a postulated reaction mechanism, i.e., the reaction rate orders, as unknowns and having the search algorithm find their optimal values given a sparse set of experimental data.

Acknowledgments M. Messina would like to acknowledge summer support from the UNCW Center for Marine Science Research. This is contribution no. 175 from the UNCW Center for Marine Science Research.

Supporting Information Available The details and numerical results of the three diagnostic tests performed to test the robustness of the search algorithm presented here (5 pages). See any current masthead page for ordering information.

Literature Cited (1) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran, 2nd ed.; Cambridge: New York, 1992. (2) Goldberg, D. Genetic Algorithms in Search, Optimization, and Machine Learning; Addison-Wesley: New York, 1989. (3) Blommers, J.; Lucasius, C. B.; Kateman, G.; Kaptein, R. Biopolymers 1992, 32, 45. (4) Dandekar, T.; Argos, P. Protein Eng. 1992, 5, 637. (5) Smith, R. W. Comput. Phys. Commun. 1992, 71, 134. (6) Makarov, D. E.; Metiu, H. J. Chem. Phys. 1998, 108, 590. (7) Hasdorff, L. Gradient Optimization and Nonlinear Control; John Wiley & Sons: New York, 1976. (8) Debock, G. J. Trading on the Edge; John Wiley & Sons: New York, 1994. (9) Back, T. Evolutionary Algorithms in Theory and Practice; Oxford University Press: Oxford, 1996. (10) Karr, C. L.; Stanley, D. A.; Scheiner, B. J. Genetic Algorithm Applied to Least Squares Curve Fitting; Int. Bu. of Mines: Pennsylvania, 1991; Report No. 9339. (11) Cvijovic, D.; Klinowski, J. Science 1995, 267, 664. (12) Carroll, D. A Fortran Genetic Algorithm Driver. www.staff.uiuc. edu/∼carroll/ga.html, 1997. (13) Neely, W. B.; Blau, G. E. In Pesticides in Aquatic Environments; Khan, M. A., Ed.; Plenum Press: New York, 1977; pp 145-163. (14) Macek, K. J.; Waber, D. F.; Hogan, J. W.; Holz, D. D. Trans. Am. Fish. Soc. 1972, 3, 420.

Received for review October 28, 1997. Revised manuscript received July 7, 1998. Accepted July 14, 1998. ES970947+