The Simulation of Dynamic Systems - Journal of Chemical Education

Information • Textbooks • Media • Resources

The Simulation of Dynamic Systems Sidney Toby* and Frina S. Toby Department of Chemistry, Rutgers, the State University of New Jersey, 610 Taylor Road, Piscataway, NJ 08854-8087; *[email protected]

Most introductions to chemical kinetics emphasize the steady-state assumption in the derivation of rate laws for complex mechanisms because the assumption that d[intermediates]/dt equals zero usually greatly simplifies the mechanistic algebra. Although the steady-state assumption is quite accurate in some systems, in others it is not, and the question of when it may or may not be valid in complex reacting systems is not easily answered. On the other hand, the exact integration of a series of rate laws is often difficult and various simplifying approaches have been described (1–3). Computer simulation has made the steady-state assumption generally unnecessary and Pavlis (4) has recently discussed the basis for the programs used. Packaged programs for desktop computers that make possible the simulation of chemically reacting systems have two benefits. First, they enable students who do not have an extensive knowledge of differential equations to analyze almost any mechanism of known rate constants and to obtain concentration–time plots for all reactants, intermediates, and products. Second, the students can use the techniques they learned in chemical kinetics to model dynamic systems outside of chemistry and analyze areas such as population trends, oscillating systems, and the viral infection of humans. Of several commercial programs available, one of us has used the one described here for several years in a course in Chemical Kinetics. The program, ACUCHEM/ACUPLOT, was devised by Braun, Herron, and Kahaner (5) at the former National Bureau of Standards (now the National Institute of Standards and Technology) and is free of copyright. The program, which uses the Runge–Kutta method, was designed for an IBM-type computer containing a mathematics coprocessor, and although it shows signs of age, it is easy to use, fast, and available (it may be downloaded from the Internet (6).) A comparison with a commercial program costing $495 (plus site license fees) showed that ACUCHEM/ ACUPLOT was faster in calculating concentration–time matrices by about a factor of 8. The absolute speed depends on the computer used and the complexity of the mechanism to be integrated. Oscillating systems take longer for simulation than other kinds of reacting processes and a comparison of a typical oscillation mechanism discussed in this article was made on three different computers. The times needed for computation on a 286, a 486DX, and a Pentium-based computer were 698, 34, and 6 s respectively. The built-in plotting program, while inelegant, gives clear, scaled graphs; however, for display in this article the concentration–time matrices were exported to another plotting program. Questions for students are marked with an asterisk when it is desirable to import the ACUPLOT file into another plotting program. The application of ACUCHEM/ACUPLOT to simulate oscillating systems has already been discussed by one of us (6 ). Since the programs use “E” for exponentiation, this terminology is used throughout the paper. Unless otherwise specified, integration tolerance should be set at 0.001. 1584

“Warm-Up” Problem We begin with a two-step system adapted from a problem in the kinetics text by Gardiner (7). A species of animal has an initial population of 1.0 animal/km2 and the population growth follows the first order reaction A→A+A

(1)

with a growth rate of 10.0% in 10.0 years. (a) How long will it take for the population to reach 1000 animals/km2? (Ans: 725 years) (b) The species becomes aggressive when meeting its own kind and destroys itself in a bianimalar reaction: A+A→C+A

(2)

with a rate constant of 5.0E{5 km2 animal {1 year{1. This prevents the population explosion that reaction 1 would otherwise cause. Assume a steady state in [A] and calculate this steady state concentration. (Ans: 191 animal/km2 ) (c) Simulate reactions 1 and 2 with the given initial concentration and rate constants and compare the resulting plot of [A] with the predicted steady-state value. How long will it take to reach 95.0% of the steady state concentration? (Ans: 860 years) Radioactive Decay and the Formation of Radon Table 1 gives the mechanism by which 238U (which has a half-life comparable to the age of the solar system) changes to the stable 206Pb. Within the sequence there are the two Table 1. Decay of 238U to 206Pb Reaction

Half-Life of Reactant

238

U→

234

Th →

234

Pa →

234

234

234

Th + 4He

4.5E9 y

Pa + β{

24 days

U + β{

1.2 min

234

U→

230

Th →

226

226

Ra →

222

222

Rn →

218

230

Th + 4He

2.5E5 y

Ra + 4He

8.0E4 y

Rn + 4He

1.6E3 y

Po + 4He

3.8 days

Po (99.96%) → 214Pb + 4He

3.0 min

218

Po (0.04%) → 218At + β{

3.0 min

Pb →

214

27 min

218

At →

214

214

Bi (99.96%) →

214

Bi (0.04%) → 210Tl + 4He

20 min

210

Tl →

1.0 min

218 214

210

Bi + 4He

210

Pb →

210

210 210

Bi →

Po →

210

2.0 s Po + β

214

Pb + β{

Po →

214

Bi + β{

210

{

20 min

Pb + 4He

1.6E{4 s

Bi + β{

22 y

Po + β{

5.0 days

206

Pb + 4He

138 days

Journal of Chemical Education • Vol. 76 No. 11 November 1999 • JChemEd.chem.wisc.edu


branches: one involving 218Po and the other, 214 Bi. To simulate the mechanism, the rate constants, obtained from the half-lives, should be multiplied by the appropriate branch percentages, although in this case the very small branch components may be neglected. Starting with an arbitrary amount of 238U, the buildup of radon can then be easily plotted. Good review articles on radioactivity (8) and radon formation (9) are available and a discussion with students of the risks from radon in basements and from nuclear power plants is well worthwhile.

Questions for Students 1. Simulate the system given in Table 1, setting the initial concentration of 238 U to an arbitrary value (say, 1.00). The α and β particles may be omitted from the rate equations for simplicity. How long does it take for [ 238U] to equal [206Pb]? 2. Does the concentration of radon reach a peak or reach a steady state? You may need to try differing time spans to answer this question. 3. Why is the formation of radon important even though its concentration is much less than that of radium? 4. (More difficult). Assuming the steady state assumption (ss) for all intermediates, calculate [Ra]ss and [Rn]ss. By comparing these values with those of the simulation at different time spans, comment on the validity of the steadystate assumption for this system. Enzyme Reactions: The Michaelis–Menten Mechanism The Michaelis–Menten mechanism in its simplest form may be written: E + S → E?S (k1) E?S → E + S (k{1) E?S → E + P (k2) where enzyme E reacts with substrate S to form complex E?S and product P. To obtain a rate law not involving the concentration of free enzyme it is necessary to make the steady-state assumption in [E?S] and also to assume that the initial enzyme concentration, [E0], equals the sum of the free and bound enzyme concentrations: [E0] = [E] + [E?S]. This gives the usual form of the rate law:

{ d[S] d[P] k 2[S][E0] = = [S] + K m dt dt

(3)

where the Michaelis constant Km = (k {1 + k 2)/k1. The restrictions on the rate law seem severe but are justified by the extremely small concentrations of free and bound enzymes compared to those of the substrate. The use of a spreadsheet to simulate enzyme kinetics has recently been described by Bruist (10), but ACUCHEM/ACUPLOT makes the simulation much easier. The resulting concentration–time matrix does not require the steady-state assumption but is not immediately verifiable. Concentrations are usually easier to measure than rates and we may check enzyme kinetics predictions by using the integrated form of eq 3: [S 0] – [S] + Kmln{[S0]/[S]} = k2[E0]t

(4)

Equation 4 predicts that a plot of [S] + Kmln[S] against time

should be linear with a slope of {k2[E0] and a y-intercept of [S0] + Kmln[S0].

Questions for Students 1. Simulate the Michaelis–Menten mechanism with k 1 = 1.06E+7 M{1s{1, k{1 = 14.9 s{1, and k2 = 1600 s {1. Take the initial concentrations of enzyme and substrate as 5.0E{10 and 0.010 M, respectively, and make a plot of [E], [S], [E?S], and [P] against time. 2. How long does it take for [P] to reach 99.0% of its final value? 3. At what time will [S] = [P]? 4. What is [E] at the end of the reaction? 5.* If we make the steady-state assumption in [E?S] and assume that [E0] = [E] + [E?S], we may integrate the resulting rate law to give [S0] – [S] + Kmln{[S0]/[S]} = k2[E 0]t. Use the values of [S] from your simulation to plot [S] + Kmln[S] against time. Evaluate the slope and intercept. Comment to instructor: Depending on their level, students may be asked to derive eq 3 or 4. When plotting eq 4 over a large time range (say 20,000 s) it is important to use adequate significant figures, or large rounding errors will develop. Population Growth—Human We are concerned here with a comparison of the population dynamics in developed and developing countries. This is a huge topic but the article by Swiegers (11) is an excellent introduction. A simplified model for human population growth can be written in two steps: P→P+P

birth occurs with rate constant kb birth rate is first order in [P]

P→C

death occurs with rate constant kd death rate is first order in [P]

In a typical developed country the average life expectancy is 75 years and an average couple produces 2.15 children in their lifetime. So the rate constant is (2.15 people/2 people) per 75 years or kb = 0.0143 y {1 and kd = 1/75 = 0.0133 yr{1. In a typical developing country the average life expectancy is 50 years and an average woman will have 6.2 children in her lifetime (or 3.1 children/parent in 50 years). Surprisingly, this crude model gives reasonable results, as described by Swiegers (11), who then suggests refinements. The refinements involve the steady-state assumption, which computer simulation makes unnecessary.

Questions for Students Use the above model for life and death rates in a developed (P) and developing (Q) country with initial conditions P0 = 5.0E7 people and Q0 = 1.0E6 people. 1. Calculate the time it will take for the population to double in (i) a developed country and (ii) a developing country. 2. When will the population of the developing country overtake that of the developed country? 3. Discuss briefly the deficiencies of the model and suggest ways by which the model could be made more accurate. The article by Swiegers (11) may be helpful.

JChemEd.chem.wisc.edu • Vol. 76 No. 11 November 1999 • Journal of Chemical Education

1585


Population Growth—Insect

Table 2. Model of an Insect Population

For this topic we are concerned with the topic of eradicating an insect that spreads disease. We then have the problem of which insecticide to use: one that incapacitates the juvenile insect or the adult insect. As in the previous topic, Swiegers’ paper (11) is a good introduction and we may then set up and solve the following problem. An insect population consists of juveniles (J) and adults (A). The population behavior may be represented as shown in Table 2.

Questions for Students Using the above model for describing insect life cycles: 1. Plot adult and juvenile concentrations ([A] and [J]) vs time for 400 days with k1 = 0.050 day{1, k2 = 7.0E{4. (This is a 2nd-order rate constant with the unusual units of a insect{1 day{1, where a ≡ are. Few people know what an are is, or what are are). Take k2B = k2C = k2D = k2E = k2F = 10 day{1, k3 = 0.075 day{1, k4 = 0.0020 day{1. Initial insect concentrations are (no pun intended) [A0 ] = 2.5 insects/a, [J0 ] = 0.50 insects/a. 2. The adult insects spread disease and it is proposed that one of two insecticides be used. At a fixed application rate, one insecticide (“J-DED”) kills juveniles with a rate constant k5 = 1.5 day{1 and the other insecticide (“A-DED”) kills adults with a rate constant of k6 = 1.5E{2 day{1. Decide which insecticide would be the more effective by calculating the adult population after 300 days, using first J-DED and second A-DED. Plot [A] and [J] vs time for the two cases. 3. Discuss briefly the sensitivity of the system with regard to insecticide efficiency. If you double k5 and k6, what happens to the respective values of [A] and [J] after 300 days? Formation of Acid Rain Acid rain is due to both sulfuric and nitric acid. A good starting reference is the paper by Weinstein-Lloyd and Lee (12); secondary references are contained within this article.

Reaction

Comment

J→A

First-order rate constant k1 , juvenile changes to adult via step 1st order in [ J].

(A + A → 50 J) Two adults die via a 2nd-order process after producing on average 50 juveniles. Kinetics simulation programs cannot in general cope with this level of fecundity, but we make a good approximation by replacing this step with the following: A+A→J+K

Second-order rate constant k2 ; K is a dummy product.

K→J+L

First-order constant k2 A .

L→J+M

k2 B ; the dummy rate constants ki are made large so these processes are essentially instantaneous.

M→J+N

k2 C

N→J+O

k2 D

O→J+P

k2 E

P→J+J

k2 F . We’ll stop here; 8 progeny are representative.

J→C

k3 , 1st-order rate constant. Juvenile dies from predation.

A→C

k4 , 1st-order rate constant. Adult gets swatted. (We don’t distinguish corpses.)

The formation of H2SO4 occurs in the liquid phase (reaction 19 in Table 3) because the gas-phase reaction between H2O2 and SO2 is immeasurably slow at room temperature. By contrast, HNO3 is formed in the gas phase (reaction 8). The initial concentrations (molecule/cm3) are: SUN, 1.0E20; ENGINE, 1.0E20; O 2, 1.0E19; O3, 3.0E14; NO, 1.0E13; hν, 1.0E18; H2O, 1.5E17; CO, 1.0E15; SO2, 2.0E12 (case 1) and 4.0E14 (case 2). First-order units are s {1 and secondorder units are (molecules/cm3) {1 s{1. The concentration unit molecules/cm3 is somewhat clumsy but is frequently used in atmospheric chemistry. ACUCHEM will not accept more than two reactants or products (note that it treats A → 2B as different from A → B + B), so reactions 3, 5, and 11 in the table are simplified to second order; and since the total pressure is constant at 1.0 atm, no error results. A concentration–time plot for this system is shown in Figure 1.

Table 3. Reactions in the Formation of Acid Rain Reaction

Comment

SUN → hν

1.0E{20

2

ENGINE → NO + SO2

1.0E{15

combustion produces NO and SO2

3

NO + NO (+ O2) → NO2 + NO2

2.0E{19

actually 3rd order

4

NO2 + hν → NO + O

1.0E{20

photolysis of NO2 produces ground-state O (3P)

5

O + O2 → O3

6.0E{15

actually 3rd order

6

O3 + hν → O2 + O*

1.0E{21

photolysis of O3 produces excited O (1D)

7

O* + H2O → OH + OH

1.0E{12

electronically excited O-atom reacts with H2O

8

OH + NO2 → HNO3

2.4E{11

9

HNO3 + H2O → HNO3(aq)

1.0E{15

OH + CO → CO2 + H

1.2E{13

10

1586

Rate Constant

1

11

H + O2 → HO2

6.0E{13

12

HO2 + NO → OH + NO2

8.2E{12

13

HO2 + HO2 → H2O2 + O2

1.7E{12

14

HO2 + H2O → HO2(aq)

1.0E{15

15

HO2(aq) → HO2 + H2O

1.0E{14

16

H2O2 + H2O → H2O2(aq)

1.0E{15

17

H2O2(aq) → H2O2 + H2O

1.0E{14

18

HO2(aq) + HO2(aq) → H2O2(aq) + O2(aq)

1.0E{12

19

H2O2(aq) + SO2 → H2SO4(aq)

1.0E{18

nitric acid rain formed actually 3rd order

sulfuric acid rain formed


Information • Textbooks • Media • Resources 3 10

Concn / (molecule cm-3)

CAR/1011

Concn / (molecule cm-3)

11

NO2/10

O3/(5x1013)

8

HNO3 /(2x1012) 13

H2SO4/10

6

4

2

0

O3/(5x1014) NO/1010 2

1

0 0

500

1000

1500

SMOG/(5x1011)

2000

0

5000

Time / s

10000

15000

20000

Time / s

Figure 1. Plot of the concentrations of sulfuric acid rain, nitric acid rain, and two intermediates in the simulation of acid rain.

Figure 2. Plot of the concentrations of smog, cars, and some intermediates in the simulation of smog formation.

Questions for Students

Questions for Students

Sulfuric acid rain (H2SO4[aq]) comes from SO2, which in turn comes mostly from engine exhausts (case 1) or from burning high-sulfur coal (case 2). Answer the following both for case 1 and for case 2. 1. Do nitric acid rain (HNO3[aq]) and H2SO4(aq) reach steady states? If so, how long does it take them to reach 95% of their final values? 2. Which is more abundant, H2SO4(aq) or HNO3(aq)? 3. Describe the behavior of [NO2 ] as a function of time.

1. Plot a graph showing how the concentrations of CAR, NO, HC, and SMOG vary over a period of several hours. If any of these concentrations goes through a maximum, give the value and the time it takes to reach the maximum. 2. Taking O, NO, HC, O3, and NO3 as intermediates and ignoring the nonchemical steps 1 and 2, what is the overall stoichiometry of the reaction? 3. Assuming steady states in the intermediates of question 2, obtain an expression for [NO] as a function of [CAR], [O 2], and rate constants only. 4. Compare your value of [NO]ss with the values from the simulation as the reaction proceeds. Comment briefly.

Formation of Smog The formation of photochemically formed urban smog occurs via a very complex mechanism (see for example Crutzen [13]). We give a simplified version in Table 4. The model simulates smog formation during the morning hours. It will not be valid when the incident light changes greatly. For the chemical reactions, first-order rate constants have units of s{1 and second-order units are cm3 molecule{1 s{1. Initial concentrations are GARAGE 1.0E12, O2 5.0 E18. A plot of this system is shown in Figure 2.

HIV Dynamics A good introduction to a very complex area is provided by Nowak and McMichael (14). A more advanced discussion is given by Coffin (15). An individual exposed to HIV may show symptoms of the virus for a few weeks. The symptoms then normally subside for several years, after which AIDS usually appears. We attempt to simulate these dynamics with

Table 4. Model for the Formation of Smog Reaction

Rate Constant

Comments

1

GARAGE → CAR

2.0E{4

2

CAR → PARK

3.0E{4

Reactions 1 and 2 "generate" cars, whose number goes through a maximum during the morning rush hour.

3

CAR → CAR + NO

2.0E{4

Cars (+ fuel) emit NO.

4

CAR → CAR + HC

2.0E{4

Cars (+ fuel) emit unburned hydrocarbons, HC.

5

NO + O2 → NO3

3.5E{19

6

NO + NO3 → NO2 + NO2

1.0E{10

Our simulation program will not accept termolecular reactions such as 2NO + O2 → 2NO2 , so it is necessary to write two steps for reactions 5 and 6.

7

NO2 → NO + O

10

This is actually a photochemical reaction, but assume that the incident light is constant during the morning.

8

O + O2 → O3

2.8E{12

A third-body reaction, written as bimolecular for simplicity at a total pressure of 1.0 atm.

9

NO + O3 → NO2 + O2

1.0E{12

—

HC + O3 → SMOG

1.0E{13

We write a single smog-forming reaction and omit the formation of PAN, etc.

10


1587


the reaction mechanism shown in Table 5. Units of rate constants are day{1 (first order) and day{1 µ L{1 (second order). Initial concentrations are in particles per microliter of blood: [B] = 5.0E12, [C] = 1.0E6, [BT] = 1.0E3, [HT] = 1.0E3, [KT] = 1.0E3, [M] = 1.0E3, [V] = 1.0E{1. A long-term plot is shown as Figure 3.

Questions for Students 1. Run the HIV mechanism for a short time span, say 25 days. Can you simulate the early peak in virus concentration? What else changes? 2. Carry out an HIV simulation over a long time span, say 7500 days. During this time what happens to the concentrations of the following species: helper T cells (HT), killer T cells (KT), antibodies (A), and virus (V)? 3. What effect does reducing [V 0], the initial virus concentration, by one and two orders of magnitude have on the time for the initial virus peak? 4. What effect does variation of k12, the rate at which the infected cell produces viruses, have on the time for the initial viral peak and the delayed buildup?

Table 5. Model of Dynamics of HIV Infection Reaction

Rate Constant

B →C

1.0E{14

2

C →D

5E{7

Cell dies.

3

B → BT

1.0E{4

Body produces B cells.

4

BT → D

1.0E{3

—

5

B → HT

9.0E{6

Body produces helper T cells.

6

HT → D

1.0E{7

—

7

B → KT

1.0E{4

Body produces killer T cells.

8

KT → D

5.0E{4

—

9

B →M

1.0E{4

Body produces macrophages.

10

M →D

1.0E{4

—

11

V + C → CV

5.0E{3

Free virus infects cell.

12

CV → CV + V

4.0E2

Infected cell produces virus. Macrophage displays epitopes.

Body produces nonspecific cells.

13

M + V → ME

1.0E{10

14

HT + ME → HTME + P

1.0E{7

Helper sees epitopes, releases protein.

15

P + BT → A + BT

1.0E{7

Protein releases antibodies from B cells.

16

A + V → AV

1.0E{8

Antibodies bind to virus.

17

M + AV → M + D

1.0E{10

Macros kill bound viruses.

18

CV + P → CVP

1.0E{8

Protein binds to infected cell.

19

KT + CVP → KT + D

1.0E{7

Killer cell destroys bound infected cell.

20

V + HT → V + D

1.0E{8

Virus kills helper cells.

21

KT + HT → KT + D

1.0E{10

Killer cells kill helper cells.

Computer Simulation of Oscillating Systems

Lotka 1: Damped Oscillations Lotka’s first mechanism (20) had as the overall reaction A → P, where reactant A was continuously generated and the mechanism involved an intermediate X, which was formed autocatalytically. The mechanism may be represented by: →A k1; in this zero-order reaction the reactant has been omitted A + X → 2X k2; second order X→P k3; second order This mechanism can be simulated using k1 = 1.0, k2 = 1.5, and k3 = 10, where the rate constants have designated units of time in s and concentration in mmol/L. Initial concentra-

8 HT/(2x10 3 ) KT/10 11

Concentration

6

M /(2x10 10 ) V / 1 0 11

4

2

0

0

2000

4000

6000

Time / day Figure 3. Plot of the concentrations of virus (V), helper T-cells (HT), killer T-cells (KT), and macrophages (M) in the simulation of HIV infection. 80 A X/10–2 60

Concentration

Oscillating systems are important in both chemistry and biology and have an extensive literature (16 ). A laboratory experiment in physical chemistry based on a nitrogen gas oscillator produced by the decomposition of ammonium nitrite has been described (17 ). Oscillations in carbon dioxide pressure occur in a soda bottle with a leaky seal and these oscillations have been discussed (18) and simulated using the Stella II program (19). The modeling of chemical oscillators is complex and has been discussed recently by one of us (6 ). We consider here oscillations in biological populations, which may be modeled using only three steps (for damped oscillations) or four steps (for undamped) and are thus more suitable as an introduction to the topic. The mathematical description and solutions to the systems that give rise to oscillating reactions were first described by Lotka for damped (20) and undamped (21) oscillations. The original papers by Lotka, involving second-order differential equations, modeled no chemical system known at that time or today. Instead the model was derived to explain biological population systems where oscillations in predator and prey populations were known to occur with time.

1588

Comment

1

P/2

40

20

0

0

20

40

60

Time / s Figure 4. Lotka 1 mechanism showing damped oscillations with time.



to a total of 50 points, but when doubt occurs about a particular region, the time scale may be changed to magnify that region. When one species in a chemically reacting system oscillates, at least one other species will oscillate in an outof-phase fashion. By plotting one oscillating concentration against another a time-free projection called a limit plot is obtained; this will show the characteristics of damped or undamped oscillation. When a limit plot of [A] vs [X] is made, using the data from Figure 4, the values spiral into a point with the result shown in Figure 5. If this graph is plotted manually in real time, the spiral is clear and the plot contrasts with a similar plot for undamped oscillations to be considered shortly. We may obtain information on the stoichiometry of the system and also obtain values for the steady state concentrations of [A] and [X] after the oscillations have died away. The Lotka 1 mechanism gives rise to the following rate equations:

[A ]

8

6

4

0

10

20

[X ]/0.1 Figure 5. Lotka 1 mechanism: limit plot of A vs X. 150

(5)

d[X]/dt = k2[A][X] – k3[X]

(6)

d[P]/dt = k3[X]

(7)

Summing and integrating eqs 5–7 gives the stoichiometry as

100

Concentration

d[A]/dt = k1[W] – k2[A][X]

[A] + [X] + [P] = [A]0 + [X] 0 + [P] 0 + k1t 50

where [X] 0 and [P] 0 are negligible. If we set d[A]/dt and d[X]/dt = 0, then steady-state concentrations for [A] and [X] are obtained: [A]ss = k3 /k2 and [X]ss = k1/k3. Questions for Students

A X/10–3 P/400

0

0

1000

2000

3000

4000

Time / s Figure 6. Lotka 2 mechanism showing undamped oscillations with time. 101

1. What would be the effect on the oscillations of (i) increasing and (ii) decreasing k1? 2.* Show that a plot of [A] vs [X] gives a spiral. 3.* Plot [A] + [X] + [P] vs time. What are slope and intercept of the resulting straight line? 4. Confirm that when the oscillations decay, [A] and [X] reach the expected steady-state value.

Lotka 2: Undamped Oscillations

[A ]

An ecological analogy follows where deer eat grass, reproduce and are eaten in turn by wolves. Wolves

100

Earth

Grass

Wolves + Wolves Deer + Deer

Deer Corpse 99

40

80

[X ] x 1000

120

160

Figure 7. Lotka 2 mechanism: limit plot of A vs X.

tions are [A] = 10 and [X] = 1.0 × 10{11 mmol/L with an integration tolerance of 10{4. Using the above values, Lotka 1 was simulated for a 75-s “experiment” and the results are shown in Figure 4. The concentrations of both A and X initially undergo large oscillations, which are rapidly damped, and the product concentration initially increases in a stepwise fashion. The program is limited

The population of both deer and wolves (X and Y in the mechanism) can oscillate. Lotka’s second mechanism (21) also had A → P as the overall reaction where reactant A was continuously generated. Two intermediates were involved, however, X and Y, both formed autocatalytically. The mechanism is: →A k4 A + X → 2X k5 X + Y → 2Y k6 Y→P k7 This mechanism was simulated using k4 = 10, k5 = 1.0, k6 = 1.0, and k7 = 0.1, where units are as before. The initial concentrations are [A] = 10 and both [X] and [Y] = 1.0 × 10{11


1589


mmole/L; the integration tolerance is set at 10{3 . Under these conditions, the concentrations of the intermediates X and Y increase rapidly and, after a period of initial instability, maintain uniform, undamped oscillations. This is shown in Figure 6 together with the product concentration. The basic difference between the Lotka 1 and 2 mechanisms appears in the limit plots. When the oscillations are damped, the limit plot spirals to a central point, seen in Figure 5. When the oscillations are undamped, the systems move through repeating values of X and Y and the limit plot eventually is a closed loop, as seen in Figure 7. Figure 7 is plotted with data from Figure 6, and is in contrast with Figure 5. Questions for Students 1. What is the stoichiometry relating A, X, Y, and P in Lotka 2? 2. Ignoring the initial region of instability, what is the effect on the period of the oscillations when k5 and k6 change? 3.* Using the above conditions for Lotka 2, are circular limit plots obtainable for [A] vs [X]? For [A] vs [Y]? 4. Derive expressions for the steady-state values of [A], [X], and [Y] and show that [A], [X], and [Y] approach their expected steady-state values. Literature Cited 1. Eberhart, J. G.; Levin, E. J. Chem. Educ. 1995, 72, 193.

1590

2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16.

17. 18. 19. 20. 21.

Gellene, G. I. J. Chem. Educ. 1995, 72, 196. Ramachandran, B. R.; Halpern, A. M. J. Chem. Educ. 1997, 74, 975. Pavlis, R. R. J. Chem. Educ. 1997, 74, 1139. Braun, W.; Herron, J. T.; Kahaner, D. K. Int. J. Chem. Kinet. 1988, 20, 51. Toby, S. Chem. Educator 1996, 1(4). Gardiner, W. C. Rates and Mechanisms of Chemical Reactions; Benjamin: New York, 1969. Hutchison, S. G.; Hutchison F. L. J. Chem. Educ. 1997, 74, 501. Atwood, C. H. J. Chem. Educ. 1992, 69, 351. Bruist, M. F. J. Chem. Educ. 1998, 75, 372. Swiegers, G. F. J. Chem. Educ. 1993, 70, 364. Weinstein-Lloyd, J; Lee, J. H. J. Chem. Educ. 1995, 72, 1053. Crutzen, P. J. Angew. Chem. Int. Ed. Engl. 1996, 35, 1758. Nowak, M. A.; McMichael, A. J. Sci. Am. 1995, 273, 58. Coffin, J. M. Science, 1995, 267, 483. See, for example, Oscillations and Traveling Waves in Chemical Systems; Field, R. F.; Burger, M., Eds.; Wiley: New York, 1985. Gray, P.; Scott, S. K. Chemical Oscillations and Instabilities; Oxford University Press: New York, 1990. Melka, R. F.; Olsen, G.; Beavers, L.; Draeger, J. A. J. Chem. Educ. 1992, 69, 596. Soltzberg, L. J.; Bowers, P. G.; Hofstetter, C. J. Chem. Educ. 1997, 74, 711. High Performance Systems: Hanover, NH. Lotka, A. J. J. Phys. Chem. 1910, 14, 271. Lotka, A. J. J. Am. Chem. Soc. 1920, 42, 1595.


The Simulation of Dynamic Systems - Journal of Chemical Education

Recommend Documents