Reviews in Computational Chemistry

234 Cellular Automata Models of Aqueous Solution Systems

Relative Lipophilicity

Figure 20 The flux of water through the membrane as a function of changing solute lipophilicity.

observation led to the speculation that the increasing concentration of solute molecules within the membrane should retard somewhat the flux of water through the membrane. Retardation was indeed found to result from increasing the PB(WS) value of the solute molecules beyond the critical point (Figure 20). At the critical PB(WS) value, the flux of water through the membrane went from a relatively constant amount to a lower level, decreasing with increasing PB(WS) values. This behavior has apparently not been noted or studied, but it provokes interest to seek experimental confirmation.

Acid Dissociation

The success with the variegated cell in creating a model of an amphiphile,62 described above, opened the possibility of using a similar type of cell to model an organic acid molecule. The dissociation of an acid and the influence of the environment on this process was investigated67 using the general reaction scheme shown in Figure 21. A cell modeling an organic acid molecule was divided into two parts, one face representing the carboxyl group Y, and the other three faces, X, representing the nondissociating, nonpolar parts of the acid molecule. The rules assigned to each face corresponded

Figure 21 The scheme employed in the modeling of the dissociation of an organic acid. Scheme (a) denotes the neutral organic acid with the Y sector representing the carboxyl group and X the remainder of the molecule. Scheme (b) denotes the hydrated organic acid. Scheme (c) denotes the water–acid ion pair. Scheme (d) denotes the dissociated acid and a hydronium ion.

236 Cellular Automata Models of Aqueous Solution Systems

from one ingredient to another across or through the system without interruption. Some objects are distributed over a space in a random fashion. Because of the scarcity of these objects, little or no physical contact is encountered. No information is exchanged within the system. If enough additional objects are randomly added to the system, there arises a finite probability that some of these objects may be associated to form clusters. There is some exchange of information within the clusters, but they are isolated and so the information exchange is confined within each cluster. If enough objects are randomly added to the system, the possibility arises that some clusters may appear as a single cluster that spans the entire length or width of the system. This spanning cluster produces a conduit over which an uninterrupted flow of information is possible across the system. This flow of information takes place via a process called percolation. The minimum number of objects in the system necessary to have a finite probability of percolation occurring is called the percolation threshold or percolation point.

Percolation is widely observed in chemical systems. It is a process that can describe how small, branched molecules react to form polymers, ultimately leading to an extensive network connected by chemical bonds. Other applications of percolation theory include conductivity, diffusivity, and the critical behavior of sols and gels. In biological systems, the role of the connectivity of different elements is of great importance. Examples include selfassembly of tobacco mosaic virus, actin filaments, and flagella, lymphocyte patch and cap formation, precipitation and agglutination phenomena, and immune system function.

The onset of percolation and the conditions that produce this phase-tran- sition phenomenon have been of considerable interest to many disciplines of science. The classical studies have focused on immobile ingredients in a system that increase in concentration by randomly adding to the collection of particles. It is possible to estimate the conditions leading to the onset of percolation under these circumstances. When the ingredients are in motion, this estimation is far more difficult. It is an obvious challenge that was tackled using cellular automata.

The studies69 were designed to reveal several attributes of the dynamic system. The first of these is the configuration of the system as additional cells were added to the grid. Of particular interest were the relative proportions of the five possible configurations of the cells, described by the concentrations of each configuration, fx, defined earlier in this review. The second attribute studied was the growth of the size of the largest cluster as the number of cells in the system increases. Of prime interest is the concentration of occupied cells at which the percolation phenomenon occurs. Also of interest were the occupied cell concentration when the onset of percolation begins and the occupied cell concentration where percolation is occurring, on average, 50% of the time. The study also looked at the information content of the system as a function of concentration change and at possible relationships among attributes of the system.

A series of runs were made starting with a concentration of 100 occupied cells in a grid of 3025 cells. Successive increases in the occupied cell concentration were introduced in increments of 100. The concentrations of the five configurations, f0 through f4, were recorded for each occupied cell concentration. The size of the largest cluster was also recorded at each occupied cell concentration. Finally, the occupied cell concentrations at which the onset of percolation takes place and at which there is a 50% probability of percolation were recorded.

The percolation point in these studies was calculated from the concentration producing a 50% probability of percolation. A comparison between these percolation onset concentrations and the concentrations corresponding to the maximum Shannon information content, Imax, revealed a close correspondence.69 The Imax concentration appears to be very close to the percolation onset concentration for each of the three parameter sets. The largest cluster size at a concentration producing 50% percolation was about the same for all parameter sets. This size is a cluster size occupying close to 20% of the grid area in these studies. Stated another way, as the concentration increases for any parameter set study, there is a relatively common cluster size corresponding to a 50% probability that percolation will occur. This interpretation is an alternative definition of a percolation point when the system is dynamic. It appears from this modeling that a maximum of diversity among the fx types is necessary to create conditions in which the percolation process may begin.

SOLUTION KINETIC MODELS

First-Order Kinetics

Many important natural processes ranging from nuclear decay to unimolecular chemical reactions are first order, or can be approximated as first order, which means that these processes depend only on the concentration to the first power of the transforming species itself. A cellular automaton model for such a system takes on an especially simple form, since rules for the movements of the ingredients are unnecessary and only transition rules for the interconverting species need to be specified. We have recently described such a general cellular automaton model for first-order kinetics and tested its ability to simulate a number of classic first-order phenomena.70

The prototype first-order transition is radioactive decay, A ! B, in which the concentration [A] of a species A decreases according to the rule that each A ingredient has a probability Pt(A,B) per unit time (here, per iteration) of converting to some other form B. For small numbers of A ingredients, the actual decay curve observed for [A] is rather jagged and only roughly exponential, as a result of the irregular decays expected in this very finite, stochastic system, as illustrated in Figure 22. However, as the number of decaying ingredients is increased, the decay curve approaches the smooth exponential fall-off

238 Cellular Automata Models of Aqueous Solution Systems

Number of Iterations

Figure 22 Exponential decays (offset) for first-order models with (a) 100 cells, (b) 400 cells, and (c) 2500 cells. The transition probability is 0.04 per/iteration.

expected for a deterministic system obeying the rate equation

This smooth decay is also illustrated in Figure 22.

When a reverse transition probability Pt(B,A) for the transition B ! A is included, the model simulates the first-order equilibrium:

A ! B

Here too, the finite size of the system causes notable fluctuations, in this case in the value of the equilibrium constant K, which fluctuates with time about the deterministic value

240 Cellular Automata Models of Aqueous Solution Systems

The series model can be extended to longer series and to the inclusion of reversibility to illustrate a variety of fundamental kinetic phenomena in an especially simple and straightforward manner. Depending on the relative rates employed, one can demonstrate the classic kinetic phenomena of a rate-limit- ing step and preequilibrium,72 and one can examine the conditions needed for the validity of the steady-state approximation commonly used in chemical kinetics.70

Kinetic and Thermodynamic Reaction Control

Parallel competing reactions

can be simulated.70 An especially interesting example is present when the reactions are reversible:

With properly chosen transformation probabilities, cellular automata can be used to examine the conditions governing thermodynamic and kinetic control of reactions.73 Knowledge of such conditions is an important element in organic synthesis74,75 and is an important consideration for many industrial reactions.76,77 Figure 24 shows a reaction being modeled with the parameters Pt(A,B) ¼ 0.01, Pt(A,C) ¼ 0.001, Pt(B,A) ¼ 0.02, and Pt(C,A) ¼ 0.0005, and a set of 10,000 ingredients. Starting with reactants A, the kinetically favored product B is produced in excess in the initial stages of the reaction, whereas at later times the thermodynamically favored product C gains dominance. The number of interacting species (10,000) is rather large in this illustration, and the results of the cellular automata model are in good agreement with those found in a deterministic, numerical solution for the same conditions.78 For example, the cellular automata model yields final, equilibrium concentrations for species B and C of [B] ¼ 0.1439 0.0038 and [C] ¼ 0.5695 0.0048 compared to the reported deterministic values of 0.14 and 0.571, respectively.

Excited-State Kinetics

Another important application of the first-order model is the examination of the groundand excited-state kinetics of atoms and molecules.79 These systems are characterized by competing first-order transitions representing both radiative and nonradiative processes. The radiative processes normally

Figure 24 Illustration of kinetic and thermodynamic reaction control: B is the kinetically favored product (higher probability of formation from A), and C is the thermodynamically favored product (greater equilibrium constant with A).

include light absorption from the ground state as well as fluorescence and phosphorescence originating from the excited states. The nonradiative processes include spin-allowed internal conversions and spin-forbidden intersystem crossings. The ‘‘species’’ in this case are the different atomic or molecular electronic states involved, represented in the visualizations by different colors. Cellular automata models can be constructed to represent both pulse conditions, which are normally used, for example, to determine excited-state lifetimes, and steady-state conditions, which are customarily employed to determine spectra and the quantum yields of the processes taking place. Pulse conditions are simulated by starting all the ingredients in a particular excited state or in a distribution of excited states characteristic of a starting distribution created by some process of interest; the states are allowed to decay according to the transition probabilities appropriate to that system, and

242 Cellular Automata Models of Aqueous Solution Systems

the time evolutions of their populations are followed. Under pulse conditions, the decays eventually lead to a condition in which all of the ingredients are in the ground state (in the absence of additional trapping states). For steadystate conditions, additional transitions from the ground to the excited states are introduced with appropriate transition probabilities.

One illustration of the excited-state cellular automata model is the dynamics of the excited-state transitions of oxygen atoms.80 The oxygen atom has a 3P ground state and 1S and 1D excited states. Emissions from the latter two excited states play an important role in the dramatic light dis- plays—the aurora borealis or ‘‘northern lights’’—seen under certain conditions in the northern polar skies and similar emissions have been detected in the atmospheres of Mars and Venus. The excited states are believed to be produced mainly by dissociative recombinations of ionized oxygen molecules and

electrons generated in the atmosphere by ultraviolet (UV) bombardment during the daylight hours.81,82

In this equation, the species O* and O** are unspecified atomic oxygen states, 3P, 1S, or 1D. The most prominent feature in the atmospheric displays is nor-

1 1 ˚

mally the green, spin-allowed S ! D transition appearing at 5577 A. Using transition probabilities taken from the compilation of Okabe,83

we simulated the dynamics associated with these atomic transitions under both pulse and steady-state conditions. For the pulse simulations, two starting conditions were examined: the first in which all ingredients started in the upper 1S excited state, and the second in which the ingredients started in a distribution believed characteristic of that produced by the dissociative recombination process in Eq. [12].82 The simulations yield excited-state lifetimes and luminescence quantum yields consistent with the experimental observations for these properties. An interesting feature arising from these studies is the possibility that a 1D/3P population inversion could occur under certain atmospheric conditions, as a result of the exceedingly long lifetime of the 1D state, estimated to be about 200 s.82

Second-Order Kinetics

Several groups have developed cellular automata models for particular reaction–diffusion systems. In particular, the Belousov–Zhabotinsky oscillating reaction has been examined in a number of studies.84–86 Attention has also been directed at the A þ B ! C reaction, using both lattice-gas models 87–90 and a generalized Margolus diffusion approach.91 We developed a simple, direct cellular automaton model92 for hard-sphere bimolecular chemical reactions of the form

Figure 25 Irreversible second-order reaction A þ B ! C þ D, with Pr(A,B) ¼ 0.1 and initial conditions [A]0 ¼ 100 cells, [B]0 ¼ 200 cells.

As before, the different species are assigned different colors in the visualization. In this model, the reactant and product species diffuse about the grid in random walks. When the species A and B encounter each other (come to adjacent cells) on the grid, the probability that these species transform to C and D is determined by an assigned reaction probability Pr(AB). The simulations take place on a toroidal space such that ingredients leaving the grid on one side appear at the opposite edge. Initially, the ingredients are placed randomly on the grid.

The production of species C over time for starting counts of 100 A and 200 B ingredients on a 100 100 ¼ 10,000 cell grid is shown in Figure 25. The expected second-order rate law