ARTICLE pubs.acs.org/jchemeduc
Nonisothermal Analysis of Solution Kinetics by Spreadsheet Simulation Robert de Levie* Chemistry Department, Bowdoin College, Brunswick, Maine 04011, United States
bS Supporting Information ABSTRACT: A fast and generally applicable alternative solution to the problem of determining the useful shelf life of medicinal solutions is described. It illustrates the power and convenience of the combination of numerical simulation and nonlinear least squares with a practical pharmaceutical application of chemical kinetics and thermodynamics, validated by comparison with literature data. KEYWORDS: Graduate Education/Research, Upper-Division Undergraduate, Physical Chemistry, Computer-Based Learning, Aqueous Solution Chemistry, Drugs/Pharmaceuticals, Kinetics, Reactions
T
extbooks are traditionally and justifiably conservative, and kinetic analysis is typically illustrated by isothermal kinetic measurements that validate kinetic laws and provide the corresponding reaction rate constants. By combining these rate constants for several such measurements on the same reacting system at different temperatures, one can then compute the corresponding activation energy, provided that the mechanism stays the same. The preceding article1 illustrates the proper analysis of isothermal reaction kinetics, which is not always quite as straightforward as it is sometimes made out to be. Industrially, however, the above, two-stage approach is often considered too time-consuming, and more efficient shortcuts are sought. One of these is to use nonisothermal kinetic analysis of, for example, medicinal solutions, for which the decomposition kinetics need be known in order to establish their useful shelf life. This is now a well-established approach, introduced to pharmacology by Rogers,2 who used an inverse-logarithmic temperaturetime profile that made the analysis of first-order irreversible kinetics mathematically tractable. It was soon followed by the use of a reciprocal heating regime.3 Further development included solving the equation for a linear temperature increase with time, and then fitting the data to resulting sets of model curves based on assumedly precise first and last data points.4,5 Alternatively, the problem was approached by numerical integration6 or by differentiation.7,8 A brief overview of the early development of these different approaches was given by Kipp.9 The combination of digital simulation with nonlinear leastsquares optimization yields yet another approach, with the advantages of conceptual simplicity, direct comparison of the model and the raw experimental data, complete flexibility in terms of temperature profile and model assumptions, and ease of implementation. Moreover, it is a good illustration of the general power of the combination of numerical simulation with nonlinear Copyright r 2011 American Chemical Society and Division of Chemical Education, Inc.
least squares, a quite general approach that can solve many mathematically well-posed problems that lack suitable closedform analytical solutions. As our software we will use the Excel spreadsheet, underscoring the inherent simplicity and general availability of this approach. We will consider a unidirectional (i.e., irreversible) first-order decomposition of the active ingredient of initial concentration c0 in a medicinal solution, so that dc ¼ kc dt
ð1Þ
where c is concentration, t is time, and k is the rate constant of the reaction, with c = c0 at t = 0. This model is not chosen for its mathematical simplicity (which is irrelevant in a digital simulation), but because all literature examples we have found happen to involve this simplest of all kinetic mechanisms. There is therefore no need to go beyond first-order kinetics in this article, but also no difficulty to extend the model to more complicated kinetics if this were required. We assume that the absolute temperature T is kept constant for a time period Δt, and is then raised in stepwise fashion to T + ΔT. For the concentration c, we will use the semi-implicit Euler method,10 which approximates the concentration c on the righthand side of eq 1 during that period Δt as c + Δc/2, that is, as the average of the initial and final reagent concentrations, c and c + Δc, respectively, in that short interval. The implied assumption here is that, over a sufficiently short interval, the concentration change is a nearly linear function of time t. The slightly simpler explicit Euler approximation,10 which takes the average concentration c over that small interval Δt as equal to its initial value, is Published: November 04, 2011 79
dx.doi.org/10.1021/ed100948n | J. Chem. Educ. 2012, 89, 79–86
Journal of Chemical Education
ARTICLE
not used here because it sometimes caused oscillatory behavior. Not only is it less stable, it is also less accurate than the semiimplicit method. In both cases, eq 1 is approximated by the corresponding difference equation Δc/Δt = kc. We therefore use Δc Δc ¼ k c þ ð2Þ Δt 2 so that identifying c with cold and c + Δc with cnew yields c þ Δc ¼
1 kΔt=2 c 1 þ kΔt=2
or 1 kΔt=2 cold 1 þ kΔt=2
ð3Þ
Figure 1. The relationship between temperature and time reported by Kipp (open circles) and a straight line (drawn) fitting these data using least-squares analysis.
where, in this case of decomposition, Δc will be a negative quantity. For the dependence of the rate constant k on the absolute temperature T, we will assume the Arrhenius equation E k ¼ k0 exp ð4Þ RT
the final temperature or a constant (in this case: zero) concentration c. The required relations only need to connect the molecular response in one time interval to that during the immediately preceding period, a convenient situation for a numerical simulation.
cnew ¼
where k0 is a constant (sometimes called the “frequency factor”), E is the activation energy of the reaction, and the gas constant R has the value 1.9858775 103 kcal mol1 K1. Now assume that the temperature, after having been kept constant for a period Δt, is raised abruptly to T + ΔT. To compute the rate constant kT+ΔT at that new temperature, we write eq 4 twice, first for T and then for T + ΔT, and find its ratio as E 0 k exp kTþΔT RðT þ ΔTÞ ¼ E kT 0 k exp RT E 1 1 ¼ exp ð5Þ R T þ ΔT T so that
E 1 1 R T T þ ΔT E 1 1 ki ¼ k0 exp R T0 Ti
kTþΔT ¼ kT exp
’ A FIRST EXAMPLE In a study of the acid-catalyzed hydrolysis of p-nitrophenyl acetate, Kipp9 took occasional samples of the reaction mixture at different times and temperatures for subsequent analysis by liquid chromatography, which yielded separate peaks for p-nitrophenyl acetate and for its hydrolysis product, p-nitrophenol. From these, he determined the concentration fractions c/C of the unreacted p-nitrophenyl acetate, where c is its concentration and C the sum of the concentrations of p-nitrophenyl acetate and p-nitrophenol. We will here use the set of experimental data reported by Kipp in his Table 3. For our simulation, it is convenient to place, at the top of an Excel spreadsheet, say in cells C3:C5, the three assumed parameters, that is, in cell C3 the initial concentration fraction c0/C of p-nitrophenyl acetate (where c0 is the value of the concentration c at the beginning of the experiment), in cell C4 the reaction activation energy E (in units of kcal mol1), and in C5 the reference rate constant k25 at 25 °C (in min1). Because a convenient layout of the spreadsheet is a significant organizing (and time-saving) part of the procedure, we will describe it here in some detail. Leaving rows 615 open for subsequent results and column headings, we copy the data of Kipp’s Table 3 to the spreadsheet, using A16:A46 for time t, with the corresponding temperatures, in °C, ranging from 35.9 to 64.1 °C, in B16:B46. In C16:C46, we then compute 1/RT, where R = 0.0019858775 kcal mol1 K1, and T is 273.15 plus the temperature in °C. Moreover, it is convenient to compute 1/RT25 =1/(0.0019858775*298.15) in cell G3, with a label in cell F3. In D16, we place the instruction for the rate constant k as =$C$5*EXP($C$4*($G$3-C16)), and again copy this down all the way to row 46. Finally, we enter Kipp’s experimental concentration ratios c/C in E16:E46. The dependence of temperature (in B16:B46) on time (in A16:A46) is illustrated in Figure 1. Linear least-squares analysis shows that it is described quite well by a linear dependence of the listed temperature (in °C) on reaction time t, that is, as a0 + a1t. The parameters found are a0 = 35.667 ( 0.063 and a1 = (0.06926 ( 0.00027) °C1, with a standard deviation of the fit, sf, of
or ð6Þ
where the reference temperature k0 is often taken at a convenient temperature. By making ΔT and Δt sufficiently small, this stepwise, discontinuous approach can also be used to approximate a gradual, continuous temperature change. Below we will therefore use a sequence of discrete periods of constant duration Δt, during each of which the temperature is kept constant at a value ΔT above that applied during the previous interval Δt. We start with the initial concentration c0 and rate constant k0 at the initial temperature T0, which is maintained for a time Δt, during which the amount decomposed is given by Δc = k0cΔt/[1 + k0Δt/2]. For the next interval we repeat this, but with c0 + Δc as our new starting concentration, and with the value of kT+ΔT appropriate for the new temperature T + ΔT. This iterative process is repeated until we have reached 80
dx.doi.org/10.1021/ed100948n |J. Chem. Educ. 2012, 89, 79–86
Journal of Chemical Education
ARTICLE
Figure 2. The experimental data reported by Kipp (open circles) and the corresponding, simulated data after (A) rough visual adjustment of the assumed fitting parameters to c0/C = 0.95, E = 30, and k25 = 1.5 104 and (B) after minimizing the sum of squares of the residuals with the nonlinear least-squares macro Solver. Figure 3. The top of the spreadsheet (A) annotated with some labels and dimensions, after use of Solver and SolverAid and (B) after the refinements described in the text.
0.182 °C. A quadratic term a2 (the usual way to represent slight nonlinearity) is unwarranted, as it would have a standard deviation s2 = 2.4 106, which is larger than the absolute value |a2| = 2.1 106 of the coefficient a2 itself. For the inverse relation of time t as a function of temperature, we find t in °C as (514.7 ( 2.8) + (14.432 ( 0.056) (T 273.15). Returning to the top of the spreadsheet, in cell F4 place the label Δt/2=, and in cell G4 deposit the corresponding instruction =(A60-A59)/2. We are now ready to start the simulation. Separated by 12 empty rows, we extend the second column, again starting with 35.9 °C in cell B59 and ending well beyond 64.1 °C, say at 70.0 °C in cell B400, but now with constant increments of 0.1 °C. In A59:A400, we then compute values for t = 514.7 + 14.432 times the corresponding temperatures in °C, in C59:C400 the associated values of 1/RT as in C16:C46, and in D59:D400 the corresponding k values as in D16:D46. To keep round-off errors to a minimum, we use the spreadsheet-generated results for the best-fitting least-squares line directly (which happens to be in cells $R$73 and $S$73 on my spreadsheet) rather than those rounded values of 514.7 + 14.432. In F59, we place =C3, in F60 the instruction =F59*(1-D59*$G$4)/ (1+D59*$G$4) for the concentration fractions c/C based on eq 3 and the assumed initial concentration fraction c0/C. Copy this down to F61:F400. We now refer in G16:G46 to the simulated c/C values in F59: F400 at the corresponding temperatures in B59:B400. We have used the 12-row offset to facilitate connecting temperatures in B16:B46 with the corresponding values in B59:B400. In G16 and G17 (for a temperature of 35.9 °C in B16 and B17), place the instruction =F59, in G18 for 37.1 °C use =F71, in G19 for
38.0 °C deposit =F80, and so on. Then, plot the experimental data in F16:F400 and the calculated ones in E16:E400 as a function of temperature in B16:B400. We can use this graph to manually adjust the assumed values for c/C, E, and k25, in order to make the simulated curve visually coincide roughly with the experimental data, as illustrated in Figure 2A. We can easily obtain a better manual fit than is shown here, by fiddling some more with the parameters, but a crude manual fit such as shown here is usually all that is required. After all this preliminary work, we are now ready for the actual data fitting. In cell G5, we use the Excel function SumXMY2 (F16:F46,G16:G46), an instruction that is easy to remember as sum of (x minus y) to the power 2, to compute SSR, the sum of the squares of the residuals, where the residuals are the differences between the experimental and corresponding simulated y values. We use Solver to minimize that sum of squares of residuals in G5, the Solver “target”, by adjusting the parameter values for c0/C, E, and k25 in C3:C5. Because the adjustable parameters differ by orders of magnitude, we must use Solver’s automatic scaling option. Starting from different, visually plausible values results in essentially identical answers. The result is shown in Figure 2B. We now use SolverAid11,12 to provide the standard deviation of the overall fit of the model to the data (in cell H5); the standard deviations of c0/C, E, and k25 (in D3:D5); and the associated covariance matrix and correlation coefficients (in B6: D8 and B9:D11, respectively). The top of the spreadsheet will now look like that shown in Figure 3A. 81
dx.doi.org/10.1021/ed100948n |J. Chem. Educ. 2012, 89, 79–86
Journal of Chemical Education
ARTICLE
D400. Note that this has only a minor effect on the resulting k values in Figure 3B, because ku,25 is about 20 times smaller than our estimate for kc,25 and, more importantly, Eu is also smaller than Ec so that, relative to the catalyzed process, the uncatalyzed reaction becomes even less important at higher temperatures. With the above changes, run Solver again, followed by SolverAid. We now find for the effective, overall rate constant k25 the value (7.154 ( 0.032) 104 min1, and for the effective activation energy E = 16.681 ( 0.035 kcal mol1, only slightly different from our earlier result.
Table 1. The Results of Changing the Number of Equidistant Time and Temperature Steps, Ns, within the Spreadsheet Intervals Used in Rows 59400, which for Ns = 1 Were Δt = 1.4432 min and ΔT = 0.1 °C c0/C
Ec /(kcal mol1)
kc,25/(104 min1)a
1
1.00549 ( 0.00043
16.681 ( 0.035
7.154 ( 0.032
10
1.00549 ( 0.00046
16.695 ( 0.038
7.109 ( 0.034
100
1.00549 ( 0.00046
16.695 ( 0.038
7.112 ( 0.034
Ns
a
In the present case, using 100 times smaller steps only has a minor effect on kc,25.
’ USING SMALLER STEPS IN TEMPERATURE AND TIME The second refinement concerns the step size Δt. We already replaced the unevenly spaced experimental data (in rows 1646) by data at regular time and temperature intervals (in rows 59400), and in principle could do the same at smaller intervals, by extending the column lengths. It is far easier, however, to leave the spreadsheet intact, and instead to use a simple custom function, as shown in the Appendix (see the Supporting Information), to perform the computation in smaller time increments. We substitute this function, which we will call NStep1, for the instruction in E60, then copy it down to row 400. The actual instruction entered in E60 should read =NStep1($G$9,A59,A60, $R$73,$S$73, D59, E59,$C$4,$C$5,$G$7,$G$8). Here cell G9 contains the value 1 for the number of steps per interval, Ns, and $R$73 and $S$73 refer to the cells containing the least-squares slope and intercept of the dependence of time (in minutes) on temperature (in °C) on my spreadsheet (these latter two addresses are most likely different on yours). The activation energy Ec and reference rate constant kc,25 of the catalyzed reaction are in cells C5 and C6, and the corresponding parameters Eu and ku,25 for the uncatalyzed process in cells G7 and G8. After copying this instruction down to row 400, try it. Nothing should happen, because with the first value in the argument of Nstep1, Ns = 1, you specified just one step, in which case the result should be identical to that obtained without the function. Then, try Ns = 10, use Solver and SolverAid, and see whether the result has changed significantly. If so, try Ns = 100, and so forth, until you get a stable reading of all the significant figures in your result. (With so many extra steps, the computation will slow down noticeably.) You can of course use different integer numbers of steps, such as Ns = 2, 4, 8, 16, .... Table 1 summarizes some results, and Figure 4 shows the final spreadsheet for Ns = 100. In principle, we could also use the spreadsheet to determine Eu and ku,25 from the experimental data, but in practice that would be nearly impossible, because Eu , Ec and ku,25 , kc,25, especially in view of the pronounced collinearity of the system, see Figure 4.
’ REFINING THE MODEL BY INCORPORATING A COMPETING REACTION The above model might be incorrect or applied incorrectly. In this case, a more complete model should include the uncatalyzed reaction and make sure that the intervals Δt and ΔT used are small enough to avoid significant systematic errors. Only after we have tied down these two possible “loose ends” can we be reasonably confident of the results. Below we will make the corresponding adjustments. Note that we need not use global weighting13 because the parameters c0/C, E, and k25 are determined directly using nonlinear least-squares, so that no transformation of variables is involved. To correct for the uncatalyzed hydrolysis of p-nitrophenyl acetate, Kipp merely subtracted the rate constant k25 of the uncatalyzed reaction, as determined by Tucker and Owen,14 from his computed result to obtain corrected rate constants. This cannot serve as a general method, although it is most likely adequate in this particular case, because Tucker and Owen reported the rate constant of the uncatalyzed reaction at 25 °C as 1.096 105 s1, and its activation energy as 14.29 kcal mol1, which are both smaller than the corresponding values for the catalyzed process. The lower activation energy of for the uncatalyzed process is of course improbable, and most likely reflects the different experimental conditions used by the different authors. Nonetheless, we here will use the same numerical values that Kipp used (even if their numerical compatibility is questionable) to indicate how, in general, an additional parallel reaction can be incorporated properly in the simulation. We therefore consider two reactions, with in this case are the (pseudo-)first-order catalyzed reaction with a rate constant kc, and the first-order uncatalyzed process with rate constant ku, together with their activation energies Ec and Eu, respectively. We replace k in eqs 13 by kc + ku, as appropriate for two parallel reaction pathways, and use Ec 1 1 ki ¼ kc, 25 exp R T25 Ti Eu 1 1 þ ku, 25 exp ð7Þ R T25 Ti
’ COMPARISON WITH LITERATURE DATA To compare our result with those in the literature, divide the rate constant kc,25 = (7.112 ( 0.034) 104 min1 by 60 s min1, and by the reported hydrogen ion concentration [H+] of (800/810) 0.25 = 0.2469 mol L1, to yield k/[H+] = (4.801 ( 0.023) 105 L mol1 s1. As can be seen in Table 2, the results obtained here agree well with the literature values. The correlation coefficients in cells C11 and D10 show significant collinearity between E and k25, as is typical of this type of problems.1 For ease of comparison, we have listed in Table 2 the activation energy and rate constant of the catalytic reaction. While these and their counterparts for the uncatalyzed reaction
to compute their values at different temperatures Ti. We use the same values as used by Kipp, viz., Eu = 14.29 kcal mol1 and ku,25 = 3.7914 105 min1 derived from measurements by Tucker and Owen at 65 °C.14 To implement this, we enter these values in cells G7 and G8, respectively, with corresponding labels in F7:F8, then modify the instruction in cell D16 to read =$C$5*EXP($C$4*($G$3-C16))+$G$8*EXP ($G$7*($G$3-C16)), and copy this to D17:D46 and to D59: 82
dx.doi.org/10.1021/ed100948n |J. Chem. Educ. 2012, 89, 79–86
Journal of Chemical Education
ARTICLE
are the chemically most meaningful quantities, their combined effect is probably more relevant for defining the solution stability, perhaps in the form of “effective” kinetic parameters Eeff and keff. A final comment: the pseudo-first-order rate constant k25/[H+] of the acid-catalyzed reaction in the first column of Table 2 has a value similar to the first-order rate constant k25 reported by Bruice and Schmir,16 and their activation energy is actually significantly higher. These measurements, obtained near room temperature in a 5.4 mM phosphate buffer at pH 8 containing 28.5% (v/v) ethanol, are clearly incompatible with the values listed in Table 2. The results of Kipp (which were the basis of our analysis) were for aqueous 0.247 M HCl containing only 1.25% (v/v) ethanol. This difference in the solvent used is possibly responsible for the discrepancy. Our analysis of Kipp’s experimental data exhibits considerable collinearity between the rate constant k and the activation energy E, see Figure 4 or cells C11 and D10 in Figure 3. Negative correlation coefficients close to 1 were also reported by Kipp.9 Figure 4 shows the corresponding error surface,1 that is, log(SSR) for various values of E and k0, where SSR, the quantity minimized by least-squares and the “target” of Solver, again denotes the sum of squares of the residuals.
fairly regular time intervals. These samples were quenched in acetic acid, and then assayed spectrophotometrically. Figure 5 shows a plot of their temperaturetime profile, which leastsquares analysis shows to be expressible quite well by the quadratic relationship T 273.15 = (20.92 ( 0.27) + (0.4456 ( 0.0098)t (0.000503 ( 0.000073)t2 where t is time in minutes. Madsen et al.6 used a seventh-order polynomial, which has an even smaller value of sf though a higher F value, as readily established with the macro LSPoly.12 For our purpose, such a high-order fit is not only unnecessary but also undesirable, because it leads to oscillatory interpolations. Both methods kept residuals to within (0.6 °C. Using reaction times t from 0 to 165 min, with increments of 0.5 min, we calculate the temperature using this quadratic relationship, and proceed as before, again initially assuming that the temperature remains essentially constant during those intervals Δt. The spreadsheet layout and data analysis procedure are therefore similar to that used previously, the only difference being that Madsen et al. used a standard reference temperature of 20 °C, which we will therefore adapt here as well. Starting from a crude first manual parameter adjustment of A0 = 0.7 absorbance units, E = 19 kcal mol1, and k20 = 0.001 min1, see Figure 6A, Solver finds A0 = 0.640385 ( 0.000051, E = 20.2574 ( 0.0079 kcal mol1, and k20 = 0.00064698 ( 0.00000098 min1, as shown in Figure 6B. The resulting values for E and k20 again exhibit a considerable collinearity, with a linear correlation coefficient of approximately 0.99. Figure 7 displays the corresponding error surface, log(SSR) as a function of E and k20, with a collinearity that is even more obvious than in Figure 4. Madsen et al. provided equivalent information in numerical form in their Table 5. Given the more than 5 times faster temperature change than in the previous section, we slightly modify the function Nstep1 to
’ A SECOND EXAMPLE As a second example, we consider the data of Madsen et al.6 on the alkaline degradation of riboflavin, aka vitamin B2. Madsen et al. used a 104 M solution of riboflavin in 0.1 M NaOH, and took temperature readings (to (0.5 °C) and reaction samples at
Figure 4. A gray scale error surface (where darker gray represents lower error) for the model used here with the experimental data reported by Kipp.
Figure 5. The relationship between temperature and time reported by Madsen et al. (open circles) and a quadratic curve (drawn) fitting these data.
Table 2. Comparison of Reported Results for the Acid-Catalyzed Hydrolysis of p-Nitrophenyl Acetate k25/[H+]/(105 L mol1 s1)
a
E/(kcal mol1)
Authors
Commentsa
Ref
5.3
—
Connors
15
—
—
18
Eriksen and Stelmach
3
i
5(1
21 ( 2
Eriksen and Stelmach
3
ni
4.75
17.1
Tucker and Owen
14
i and ni average.
4.1 ( 2.1
17.8 ( 5.7
Kipp
9
ni, derivative method
4.4 ( 0.1
17.4 ( 0.5
Kipp
9
ni, integration method
4.80 ( 0.02
16.69 ( 0.04
this analysis
ni, numerical simulation
i denotes isothermal and ni nonisothermal. 83
dx.doi.org/10.1021/ed100948n |J. Chem. Educ. 2012, 89, 79–86
Journal of Chemical Education
ARTICLE
Figure 7. A gray scale map of the error surface for the model used here with the experimental data reported by Madsen et al.
Table 3. The Results of Changing the Number of Equidistant Time and Temperature Steps, Ns, within the Spreadsheet Intervals Used in Rows 40 through 340, which for Ns = 1 were Δt = 0.5 min and ΔT = 0.3701 °C Ainit
E/(kcal mol1)
k20/(104 min1)a
1
0.6404 ( 0.0013
20.26 ( 0.23
3.23 ( 0.14
10
0.6404 ( 0.0013
20.28 ( 0.23
3.19 ( 0.13
100
0.6404 ( 0.0013
20.28 ( 0.23
3.19 ( 0.13
Ns
Figure 6. The experimental data reported by Madsen et al. (open circles) and the corresponding, simulated data after (A) rough visual adjustment of the parameters to A0 = 0.7, E = 19, and k20 = 0.001 and (B) after minimizing the sum of squares of the residuals with the nonlinear least-squares routine Solver.
a
Using more steps than with Ns = 100 has only a minor effect on k20, but increases its uncertainty (because of accumulative truncation errors) and computation time.
Nstep2 in order to see whether this result needs refinement. Note that such refinement only concerns the computation, and yields an estimate of the random errors (as a standard deviation) but does not answer questions such as whether, for example, the temperature in the reaction vessel was uniform, or whether other systematic errors might have been involved. Madsen et al. provide a comparison with literature data. For ease of comparison, we multiply our final value in Table 4 for k20 by 60 min h1 and divide by [OH] = 0.1 mol L1 as appropriate, and made similar adjustments in some of the other reported data.
(2.353.12) 1012 min1. Our single result falls well within these ranges. Moreover, Hempenstall et al. listed a number of values for three different data sets (all apparently different from the above-quoted one), for both differential and integral nonisothermal analysis, the latter using the method of Madsen et al.,6 as well as for isothermal measurements. They reported values for E ranging from 86.5 to 94.5 kJ mol1, and for k0 from (0.438.41) 1012 min1, where both ranges again widely bracket our result. As in examples 1 and 2, we observe a linear correlation coefficient between E and k0 of near unity, in this case of 0.99989, suggesting why these numbers may vary so much depending on experimental noise in individual measurements, because collinearity makes this analysis method quite sensitive in this respect. The error surface is illustrated in Figure 9, and may appear to show a different trend because the dependence of k0 on E in the “trench” is opposite from that of k20 or k25 on E. This paradox indeed disappears when we instead compute and plot k20 or k25 versus E.
’ A THIRD EXAMPLE As our final exampl,e we will use the data set reported by Hempenstall et al.8 who used the nonisothermal derivative approach to study the stability of a solution of potassium phenoxymethylpenicillin buffered at pH 9, and published one of their data sets. Use of LSPoly1 shows that its temperature can be fitted quite well by a relatively low-order polynomial in time, as illustrated in Figure 8, as long as the first data point is excluded, as was also done by Hempenstall et al. With that polynomial, we then generated a table of temperatures at 1-min intervals, and subsequently interpolated it using smaller time intervals with the slightly adjusted version Nstep3, which uses a fourth-order polynomial to compute the temperature as a function of time, and is based directly on the Arrhenius eq 4. We find c0/C = 0.99962 ( 0.00070 (which is not statistically different from 1), E = 91.50 ( 0.24) kJ mol1, and k0 = (2.93 ( 0.26) 1012 min1. For this particular data set, Hempenstall et al. reported E = 91.3 kJ mol1 within a range from 90.9 to 91.7 kJ mol1, and a k0 -value of 2.71 1012 min1 within the range
’ DISCUSSION In this article, we have emphasized the spreadsheet layout more than the details of the specific software used, because of the visual nature of the spreadsheet and the importance of a transparent data organization. The software used has been described earlier,10 is freely available for noncommercial use,12 and is self-documented and open-access, which means that it can be examined in detail by anyone so inclined, can be modified to adapt to specific needs, and can even be pilfered for useful parts. We have illustrated a fast, efficient, quite general, and readily accessible method to estimate the kinetic temperature stability of 84
dx.doi.org/10.1021/ed100948n |J. Chem. Educ. 2012, 89, 79–86
Journal of Chemical Education
ARTICLE
Table 4. Comparison of Reported Rate Parameters for the Base-Catalyzed Degradation of Riboflavin k20/[OH]/ (L mol 1 h1)
E/ (kcal mol1)
Authors
Ref
Comments
0.18
17.85
Rogers
2
ni*
0.16
19.2
Guttmann
17
i*
0.14 ( 0.02
20.39 ( 0.05
Madsen et al.
6
ni, integration method
0.22 ( 0.005
20.12 ( 0.04
Madsen et al.
6
ni, integration method*
0.192 ( 0.008
20.3 ( 0.1
this analysis
ni, numerical simulation
i denotes isothermal, ni nonisothermal. Base used 0.1 M (where indicated with *) or 0.05 M NaOH.
Figure 8. The temperaturetime profile used by Hempenstall et al.8 (open circles), and a fourth-order power series with a0 = 26.88 ( 0.11, a1 = 0.650 ( 0.011, a2 = (1.19 ( 0.34) 103, a3 = ( 4.17 ( 0.40) 105, and a4 = (1.46 ( 0.16) 107 fitted to these data (small solid points connected by thin line).
Figure 9. The error surface for our analysis of the penicillin data of Hempenstall et al.
but these do not appear to present any serious impediments in the cases shown here. The use of nonlinear least squares avoids the need for special derivations, and circumnavigates the constraints posed by linear least-squares methods. The present approach is also sufficiently transparent to be practicable in the undergraduate physical chemistry laboratory. In general, one can do far more with the combination of numerical simulation and nonlinear least-squares than one can ever hope to achieve with closed-form solutions and linear least squares, and future chemists should be exposed to these easily applied yet very powerful data analysis tools, which are unlikely to disappear any time soon, but instead are certain to evolve further, because the combination of simulation and nonlinear least squares can form an intuitive, highly flexible, readily available yet very potent general approach to data analysis, especially now that simple visual tools are available to guard against local minima.1
solutions of biological or medicinal importance. Already in 1976, Kulshreshtha18 listed well over 400 publications in this field, which has rapidly expanded since that time. The method proposed here is fast, readily adaptable to any reaction mechanism or combination thereof, of any (integral or other) order, can be used with any temperaturetime profile, is conceptually simple and easily implemented, has built-in checks for numerical accuracy, and provides estimates of the precision of its results. It is described here mainly to illustrate the general power and convenience of numerical simulation coupled with nonlinear leastsquares analysis. It can of course be implemented in any numerical software system, and is here illustrated with Excel, at present the most ubiquitous numerical software available. When used routinely, it might be worthwhile to use a more accurate fourth-order RungeKutta algorithm instead of the semiexplicit Euler method used here for reader transparency, in which case we may not even need a custom function such as NStep. In the above examples, we have assumed first-order kinetics because these were known to be applicable in these cases, and therefore sufficed to demonstrate that the proposed method works. Moreover, the method will often yield better or worse fits depending on the appropriateness of the kinetic model used, but that may not always be a reliable criterion. When the kinetics are unknown, users would therefore do well to perform preliminary isothermal experiments to establish the type of kinetic model applicable to their system before attempting a nonisothermal approach. Some problems inherent in experimental studies of chemical kinetics, such as the collinearity illustrated in the companion publication (DOI: 10.1021/ed100947d),1 cannot be avoided,
’ ASSOCIATED CONTENT
bS
Supporting Information Appendix showing how to use a simple custom function. This material is available via the Internet at http://pubs.acs.org, and in ref 12.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
’ ACKNOWLEDGMENT I thank Theresa Ziemba of Bristol-Myers Squibb, a participant in one of my PittCon short courses, for calling my attention to 85
dx.doi.org/10.1021/ed100948n |J. Chem. Educ. 2012, 89, 79–86
Journal of Chemical Education
ARTICLE
this problem, and Panos Nikitas and Carl Salter for their very helpful comments.
’ REFERENCES (1) de Levie, R., J. Chem. Educ. DOI: 10.1021/ed100947d. (2) Rogers, A. R. J. Pharm. Pharmacol. 1963, 15 (Suppl.), 101T–105T. (3) Eriksen, S. P.; Stelmach, H. J. J. Pharm. Sci. 1965, 54, 1029–1034. (4) Zoglio, M. A.; Windheuser, J. J.; Vatti, R.; Maulding, H. V.; Kornblum, S. S.; Jacobs, A.; Hamot, H. J. Pharm. Sci. 1968, 57, 2080–2085. (5) Maulding, H. V.; Zoglio, M. A. J. Pharm. Sci. 1970, 59, 333–337. (6) Madsen, B. W.; Anderson, R. A.; Herbison-Evans, D.; Sneddon, W. J. Pharm. Sci. 1974, 63, 777–781. (7) Waltersson, J.-O.; Lundgren, P. Acta Pharm. Suec. 1982, 19, 127–136. (8) Hempenstall, J. M.; Irwin, W. J.; Li Wan, Po, A.; Andrews, A. H. J. Pharm. Sci. 1983, 72, 668–673. (9) Kipp, J. E. Internat. J. Pharm. 1985, 26, 339–354. (10) de Levie, R. Advanced Excel for Scientific Data Analysis, 3rd ed.; Atlantic Academic LLC: Orr's Island, ME 2011. (11) de Levie, R. J. Chem. Educ. 1999, 76, 1594–1598. (12) These Excel macros are part of the MacroBundle, freely downloadable from http://www.bowdoin.edu/∼rdelevie/excellaneous. (accessed Oct 2011). (13) de Levie, R. J. Chem. Educ. 1986, 63, 10–15. (14) Tucker, I. G.; Owen, W. R. Int. J. Pharm. 1982, 10, 323–337. (15) Connors, K. A. Interest 1963, 2, 51–62, as quoted by Kipp.8. (16) Bruice, T. C.; Schmir, G. L. J. Am. Chem. Soc. 1957, 79, 1663–1667. (17) Guttmann, D. E. J. Pharm. Sci. 1962, 51, 1162–1166. (18) Kulshreshtha, H. K. Def. Sci. J. 1976, 26, 189–204.
86
dx.doi.org/10.1021/ed100948n |J. Chem. Educ. 2012, 89, 79–86
