A Cellular Automata Model of Diffusion in Aqueous Systems

A Cellular Automata Model of Diffusion in Aqueous Systems LEMONT B. KIER*X, C.-K. CHENG†, BERNARD TESTA‡,

AND

PIERRE-ALAIN CARRUPT‡

Received February 4, 1997, from the The Departments of *Medicinal Chemistry and †Mathematical Sciences, Virginia Commonwealth University, Richmond, Virginia 23298, and the ‡Institute of Medicinal Chemistry, University of Lausanne, CH-1015 Lausanne, Switzerland. Final revised manuscript received April 1, 1997. Accepted for publication April 9, 1997X. Abstract 0 A cellular automata model of a solute diffusing in water has been created and studied for the influential attributes. The results with this model are in agreement with experimental results; that is, that lipophilic solutes diffuse faster than do polar solutes. The model reveals that a solution composed of a relatively lipophilic solute permits a greater extent of diffusion of another solute. This Åobservation is in agreement with the model showing a diffusion preference of a solute between two solutions made up of differing polarities. The solute diffuses farther into the lipophilic solution. A temperature−lipophilicity phase diagram shows the influence of these two attributes on the rate of diffusion. A model of diffusion through solutions containing stationary ingredients reveals a faster rate when the ingredient is lipophilic. We are led to a conclusion that the relative lipophilicity of solutes or stationary ingredients in a solution has a direct influence on the rates of diffusion of other solutes in their midst.

Introduction The diffusion of small, organic nonelectrolytes in water is an important process linked to biochemical and pharmacological effects. Diffusion is the process whereby ligands pass through bulk water to the vicinity of an effector. There is considerable information on aqueous diffusion scattered through the literature, but it is usually in the form of studies on a few molecules with different protocols for measurement. Thus, useful information for systematic studies is scarce. Two studies, however, have addressed this issue and have presented data suitable for analysis.1,2 We will refer to these studies later to compare their results with our models. As an alternative to information derived from experiments, there is increasing attention being paid to the use of models to simulate some dynamic process and to generate information that may be difficult to acquire experimentally. These models, sometimes referred to as research in silico, are computer simulations using molecular dynamics, Monte Carlo, or cellular automata dynamics. We have been interested in the potential of this latter method to study water and solution phenomena.3-8 This method appears to offer possibilities to study the influence of variable lipophilicity changes on the diffusion of small organic molecules in water. The rationale for the use of this approach to model aqueous diffusion is based on its demonstrated potential to simulate the physical properties emerging from the dynamic interactions of subsystem ingredients. Complex systems arise from bottom-up encounters based on forces that may be represented by rules. Because the complex system is more than the sum of its parts, a dynamic synthesis, as with cellular automata, creates a process that may produce a useful model of such a system. Clearly, water and aqueous systems are complex and so we seek to explore the potential for this model.

Theoretical Section Cellular AutomatasCellular automata are computergenerated dynamical systems that are discrete in space, time, X

Abstract published in Advance ACS Abstracts, June 1, 1997.

774 / Journal of Pharmaceutical Sciences Vol. 86, No. 7, July 1997

state and whose behavior is specified completely by rules governing local relationships. It is an effort to simplify the often numerically intractable dynamic simulations into a set of rules that mirror intuition and are easy to compute. As an approach to the modeling of emergent properties of complex systems, it is beneficial to be visually informative about the progress of dynamic events at the molecular system level.7 Since the development by von Neumann,9 a variety of applications ranging from gas phenomena to biological applications have been reviewed by Ermentrout and EdelsteinKeshet.10 We have used cellular automata to advance our understanding of water and solution phenomena and have embarked on a series of studies with this goal in mind.3-8 Our model is composed of a grid of squares called cells (Figure 1). Each cell, i, has four adjoining neighbors, j, and four extended neighbors, k, beyond j. This configuration is called an extended von Neumann neighborhood. Each cell is assigned a state governing whether it is empty or is occupied by a water molecule or other ingredients. The contents of a cell may break away from an occupied neighboring cell or move to join a neighboring cell that is occupied (see equations in the Appendix). These trajectories are assigned as probabilistic rules at the beginning of the dynamics to reflect an anticipated relationship among the ingredients in the system. The rules are then applied randomly to each cell in turn until all cells have computed their state and trajectory. This procedure is one iteration of time. The rules are applied uniformly to each cell of the same state. The initial state of the system is random, hence it does not determine subsequent configurations at any iteration. The same set of rules do not yield the same configurations except in some average sense. The configurations after many iterations reach a collective organization that possesses a relative constancy in appearance and in reportable counts of cells and are called attributes. These attributes are the emergent characteristics of a complex system. Traditional studies of solution phenomena have been conducted at the bulk system level using kinetic and thermodynamic approaches. In these cases, there are an Avagodro's number of molecules in the system. Any mechanism purported to be at the molecular level is the result of an averaging process assigning the attributes found to molecules as a probability. All events at the molecular level can only be viewed as being uniform from such a process. More recently, the use of quantum mechanical methods, such as molecular orbital theory, has come into practice to study single molecules and very small numbers of molecules. These calculations, treating the ensemble as a united molecule, have attempted to model a phenomena from this sparse collection. Extrapolation of results to much larger numbers in the system has constituted an effort to model the large system. Cellular automata studies of dynamic systems make it possible to model the intermediate level in which the properties are manifested as measurable phenomena but at which the individual roles of the ingredients are discernible in the model. This approach makes it possible to model the progress of encounters in the system that give rise to the measured properties, thereby making it possible to approach the design

S0022-3549(97)00051-8 CCC: $14.00

© 1997, American Chemical Society and American Pharmaceutical Association

Figure 1sThe extended von Neumann neighborhood used in these studies.

of molecular level mechanisms without invoking an averaged behavior. Alternatives to cellular automata to study complex systems at this level include molecular dynamics and Monte Carlo simulations. The problems with these alternatives are the vast amounts of computer time needed to generate very brief intervals of dynamic events and the limited variety of ingredients and environments that can be employed in the simulations. We believe from our experience that cellular automata has advantages of speed and versatility in creating simulations of dynamic events at the molecular level. The objective of this study and of the others in this series is twofold. First, we want to explore the potential for cellular automata models to simulate enough of the physical properties of the molecular systems so that some degree of confidence can be achieved. In each of our studies, we began by running the dynamics and recording emerging attributes. These data have been interpreted and compared with experimental evidence for agreement, and correspondence was achieved in each study. From this preliminary base of confidence, we proceeded to look for new and unexpected attributes from the dynamics. These observations form the basis of predictions or possibilities of phenomena as yet unstudied or unforeseen. These new attributes and their physical interpretations form the basis of ideas about molecular systems that might be explored experimentally. Thus, our studies are designed to offer possible mechanisms for molecular system events and to create a pool of idea leads. We can say that these models, if valid in preliminary simulations, may offer answers to unasked questions. The Modeling RelationshipsUsing cellular automata dynamics we have recently studied models of water,3 solution phenomena,4 the hydrophobic effect,5 solute dissolution,6 immiscibility and solute partitioning,7 and micelle formation.8 Within the context of these studies it should be made clear what the cells, the configurations generated, and the cellular automata model represent. This definition is important to derive any understanding of the results of a simulation and to dispel misunderstanding based on direct comparisons with molecular models. A cell with a state value encrypting occupation by a particular ingredient is not a model of a molecule with specific electronic and steric features, but is a statement of the existence in space-time of an object that has certain rules governing its trajectory. These rules govern the transaction of the ingredients in a cell with all other ingredients in neighboring cells, the von Neumann neighborhood. Each occupied cell is interpreted to be a molecule, specified or not, that is defined only by its state and transition functions. The electronic and topological characterizations are considered to be subsumed into these rules. The transactions in the dynamics result in configurations after several iterations that involve more than a single molecule. There is simulated after some number of iterations, the emergent behavior of a molecular system. We define here a molecular system as the minimum number of molecules (or their

Figure 2sThe phase diagram depicting the influence on diffusion of solute lipophilicity and water temperature. The diffusion extent is in terms of the average count of cells separating the diffusing solute from a central point on the grid, at some iteration time. The temperature of the water is expressed by the breaking probability of water, PB(W). The lipophilicity of the solute, S, is expressed as the breaking probability of a solute-water pair, PB(WS).

modeled surrogates) necessary to model a phenomenon which is recognizable as emergent behavior. This emergent behavior is epitomized by visual patterns or calculable attributes, such as the average counts or size of configuration features. The molecular system is intermediate between the isolated molecule level and the bulk phase models of systems. Molecular systems may be modeled with molecular dynamics, Monte Carlo, or cellular automata dynamics. Molecular level phenomena are modeled with quantum mechanical-based methods, whereas bulk phase phenomena are modeled with statistical and thermodynamic methods.

Results First Study: Solute Polarity and DiffusionsThe relationship of a solute to the solvent, water, in terms of relative polarity or lipophilicity is certainly an influential factor in the diffusion through the solvent. In this study we evaluated this relationship with a series of rules while monitoring the dynamics in a cellular automata simulation. The polarity of the solute, S1, was encoded into the water-solute breaking probability, PB(WS1), so that a high value encodes a relatively nonpolar solute. At the same time, a low value of the joining parameter, J(WS1), encodes a low affinity between these two ingredients. Twenty solute molecules were introduced into the central 200 square area of the grid. The net average distance traveled by the 20 solute molecules from a central point after 1000 iterations was recorded. This distance was averaged over 10 runs and was expressed as the diffusion experienced by this particular solute type. The water breaking probability, PB(W), was varied from 0.20 to 0.70, modeling a system temperature of 20 to 70 °C.3 Several values of the PB(WS1) and J(WS1) rules were used in the dynamics to generate a profile of the polarity influence. The diffusion was expressed as the average difference between the starting, random average distance from the center at iteration zero, and the average distance traversed after 1000 iterations. These results are shown in a phase diagram in Figure 2. The results indicate that nonpolar solutes, in our model, diffuse farther than polar solutes. This result was a conclusion that we reached in an earlier cellular automata model of the dynamics of solute dissolution from a solid mass.6 The conclusion is consistent with experimental studies in which

Journal of Pharmaceutical Sciences / 775 Vol. 86, No. 7, July 1997

Table 1sDiffusion Coefficients of Molecules of Comparable Volume Molecule

MW

D24 × (106) cm2/sa

D24 × (106) cm2/sb

Ethanol Formamide 2-Propanol Acetamide Butanol t-Butanol Propionamide Pentanol (15 °C) Ethyl acetate (15 °C)

46 46 60 60 74 74 74 88 88

8.25 11.30 7.14 8.30 6.68 s 7.31 4.44 6.24

1.26 1.72 1.07 1.32 s 0.98 1.20 s s

a

b

Data from ref 1. Data from ref 2.

the relative polarity of solutes was evaluated as an influence on diffusion1,2 as shown in Table 1. In these two studies, conducted under different conditions, the less hydrogenbonding molecule (less polar) of the same molecular weight diffused the farthest in each study. In each case, this molecule was the amide or ester rather than the alcohol. With increasing volume, as exemplified by the molecular weight, the diffusion coefficient was reduced. In our simulation, all cells were modeled to be of the same size, and our results (Figure 2) are in agreement with this experimental evidence. We found that less polar molecules, with higher PB(WS1) values, diffuse farther. Second Study: Water Temperature and DiffusionsIn this study we varied the temperature of the water and modeled the diffusion of a solute as a function of its relative polarity. The concentration of the S1 solute was 20 cells randomly scattered in the center of the grid. The polarity of the solute relative to water was modeled with the parameter, PB(WS1). The water temperatures were modeled by changing PB(W). The average distances of diffusion, expressed as cell counts, were recorded at 1000 iterations with multiple runs of 10. The effect of temperature change coupled with the relative polarity of the solute on the diffusion is reflected in the phase diagram of Figure 2. Increasing the temperature resulted in an increase in the rate of diffusion as described by the Einstein equation.11 This finding is in agreement with much experimental evidence. Third Study: Diffusion Through a SolutionsA Single Solution PhasesAn important aspect of diffusion that might lend itself to study using cellular automata is the diffusion of a solute, S1, through a solution containing a cosolvent, S2. Our interest is in the effect of the polarity of S2 on water properties that influence the diffusion of S1. To model this situation, we introduced 280 cells at random in the water, representing a cosolvent, S2. There were 20 S1 solutes introduced into the center of the grid as in the previous study. The relationships between the solute, S1, and any other S1 were encoded into the rules, PB(S1) and J(S1). The interaction rules for the two different solutes were encoded into PB(S1S2) and J(S1S2). Of particular interest in this study are the relative polarities of the two solutes. These were encoded into the water-solute parameters, PB(WS1) and PB(WS2) and the corresponding J parameters. The dynamics were run 10 times for 5000 iterations each. Several variations of the polarity of the two solutes were included in the studies. The results are shown in Table 2. In a separate study we observed that the concentration of the S2 solute is inversely related to the diffusion rate, a finding coinciding with the solvent viscosity influence. The dominant influence on the diffusion of solute S1 is its polarity. The less polar (more lipophilic), the greater the extent of diffusion, as is shown in Table 2. The polarity of the cosolvent, S2, has some influence on the extent of diffusion of S1. A less polar S2 results in a greater diffusion of S1. In all cases, the diffusion of S1 was more extensive in pure water. 776 / Journal of Pharmaceutical Sciences Vol. 86, No. 7, July 1997

Table 2sEffect of Relative Polarities on the Diffusion of a Solute, S1, Through a Solution

PB(WS1) )a

PB(WS2

0.25

0.50

0.75

0.25 0.90

3.69 3.95

8.71 10.11

9.62 11.35

a When no solute S is present and P (WS ) ) 0.75, the average distance of 2 B 1 diffusion of S1 is 12.67 cells.

Table 3sAverage Extent of Diffusion into a Polar and a Nonpolar Solution Solution

Average Distance of Diffusion of S1 (cells)

Average Number of Cells Diffusing

Nonpolar Polar

16.81 14.75

9.52 7.96

Two Bulk-Phase SolutionssIn a second study to evaluate the influence of polarity on the extent of diffusion, we created a model of two solutions competing for the diffusion of a solute. A thin layer of solute S1 was randomly placed in a horizontal layer between two solution phases composed of solutes S2 and S3, respectively. The S2 solute was made relatively nonpolar with the parameters PB(WS2) ) 0.90 and J(WS2) ) 0.25. The S3 solute was made relatively polar with the parameters PB(WS3) ) 0.25 and J(WS3) ) 1.0. The diffusing solute S1 was given an intermediate polarity with the parameters PB(WS1) ) 0.50 and J(WS1) ) 0.50. The system consisted of a grid of 3025 cells, with 1880 being designated as water, 20 as S1, and 100 each as S2 and S3. The dynamics were run 10 times for 2000 iterations each to produce average values of the distance S1 diffused into each solution. The number of S1 molecules diffused into each solution was also recorded. The results are shown in Table 3. More S1 molecules diffused farther into the nonpolar solution S2 than into the more polar solution of S3 solutes. This result is in agreement with the previous study with a single solute with varying polarity. Fourth Study: Diffusion into Realms of Stationary IngredientssIn this study we used cellular automata to model the diffusion of a solute, S1, through several realms of water containing stationary ingredients of varying polarity. The purpose of this study was to simulate a situation in which a solute is exposed to stationary ingredients in water, such as the array of amino acid side chains on the surface of a protein molecule. A solution of 20 S1 solutes randomly scattered within 100 cells containing water were positioned in the center of a 3025-cell grid. In four regions (A, B, C, and D) surrounding this central area were placed stationary ingredients of varying polarities. The layout of the different regions are shown in Figure 3 and are described as follows: (A) 20 stationary cells, S2, randomly scattered in an area of 150 cells containing water [the polarity of S2 is encoded by PB(WS2) ) 0.20 and J(WS2) ) 2]; (B) 10 stationary cells, S2, and 10 stationary cells, S3, randomly scattered in an area of 150 cells containing water [the polarity of S3 is encoded by PB(WS3) ) 0.90 and J(WS3) ) 0.25]; (C) 20 stationary cells, S3, randomly scattered in an area of 150 cells containing water; (D) water cells with no ingredients. The mobile S1 cells representing a mobile solute were encoded with intermediate polarity values [PB(WS1) ) 0.50, J(WS1) ) 0.50]. These cells were allowed to respond to the dynamics for a period of 2000 iterations, starting from a random distribution. The dynamics were repeated 50 times to obtain average values of solute movement into the various regions. These values are recorded in Table 4. The solute,

Figure 3sA model exposing different regions (A, B, C, and D) containing solutions with different lipophilicity characteristics to a mobile solute, S. Table 4sDiffusion of a Solute Through Stationary Ingredients in Water

Region

Average Distance Traveled (Cells)

Maximum Distance Traveled (Cells)

Number of Solutes in Region

Boundary Crossings

Center A B C D

3.85 7.93 8.31 9.73 11.58

5.29 9.08 9.78 13.31 14.94

5.00 1.50 1.98 2.94 2.98

s 57 92 115 118

S1, diffuses farther into the pure water region as we have observed before. Among the other regions in Figure 3, the solute diffused farther and in greatest numbers into the most lipophilic region, C, and diffused the least into the most hydrophilic region, A.

solutes with varying degrees of lipophilicity. This simulation creates an interesting and possibly important model simulating possible influences on this process. Another model is created in which a solute is exposed to two possible solution realms into which it may diffuse. The solute preferentially diffused into the solution realm composed of a more lipophilic solute. This model could be likened to the influence of a concentration gradient on the diffusion of a solute. In this case, it is the general polarity of a second solute influencing a preferred path of diffusion of another solute. This scenerio raises the possibility of selected diffusion within cells governed by factors other than concentration. The final study was an attempt to model an environment simulating stationary ingredients with varying lipophilicities. Such a situation would arise on the surface of a protein molecule where a field of side chains intrudes into bulk water. Our model shows a result that is interpreted as a favored path of a mobile solute through a field of mostly lipophilic stationary ingredients. This scenerio could be of interest in assessing the diffusion of a drug through bulk water to a receptor. These are models in silico, and they cannot replace experiments in situ. However these models are sources of new ideas and possibilities that may emerge from a simulation that reproduces some elements of reality. A complex system defies full understanding even from the most sophisticated experimental design. Its attributes emerge from many complicated interactions of ingredients.12 Experimental results can only focus on a few of these events; we use these events to make a provisional model of the whole. A cellular automata model may be looked on the way we regard a quantum mechanical model of a molecule. The model simulates systems and their events, yielding attributes that we interpret to explain observations, which is all that we can ever do to understand nature. Cellular automata is thus put to the test in this study to see how far it can go in the quest for understanding.

Discussion This study was undertaken for two purposes. The first, a general one, was to extend our evaluation of the method of dynamic simulation using cellular automata. In a series of papers we have shown that some water and aqueous solution phenomena can be reproduced using stochastic rules relating cells designed to represent water and solutes. The initial challenge in any of these studies has been to establish some credibility for the models in regard to experimental events. We have achieved some measure of success in these studies in this regard. In the current study we began by reinforcing our belief that we can mirror reality by again showing that diffusion of a solute is influenced by the water temperature and its lipophilicity. A more specific objective in each study, and in this one as well, is to discover if there are some attributes emerging from the dynamics that present new sets of conditions. These attributes may be interpretable as unexplored or unforeseen physical phenomena that are a source of new insight into the behavior of the system and a source of new ideas for further experimental design. This potential for discovery is one of the values of in silico models. The present study reveals that a model simulating a solute, a, exhibits a reduced diffusion in the presence of another solute, b. The extent of that decrease is dependent on the lipophilicity of the second solute, b; the more hydrophilic is b, the greater the decrease in diffusion of a. Whether this relationship is generally observed is not clear, but the idea is of interest because no biological milieu is pristine, with only one solute. It is of further interest when one considers that the passage of a neurotransmitter across a synapse is by diffusion. Certainly within the synapse is a collection of

Appendix Two parameters are selected for a simulation to govern the probabilities for movement of the contents of the cells in the grid. The breaking probability, PB(XY) is the probability for molecule, X, in cell i (Figure 1) to detach from the molecule of type Y, at some j cells when there is exactly one occupied j cell. The value for PB(XY) lies in the closed unit interval, and when X and Y are the same molecular types, it is denoted by PB(X) or PB(Y). The second parameter, J(X), describes the movement of the molecule, X, in cell i toward or away from the molecule Y, in cell k in the extended von Neumann neighborhood when the intermediate j cell is vacant. This second parameter represents the ratio of the probability that a molecule of type X at cell i will move toward a k cell occupied by a molecule of type Y while the intermediate j cell is vacant, and the probability that a molecule of type X at cell i will move toward a vacant k cell while the intermediate j cell is vacant. J is a positive real number. When J(XY) ) 1, the molecule in cell i has the same probability of movement toward an occupied cell k as when cell k is empty. When J(XY) > 1, the molecule in cell i has a greater probability of movement toward an occupied cell k than when cell k is empty. When J(XY) < 1, that molecule i has less probability of movement toward an occupied cell k than when cell k is empty. When X and Y are the same molecular types, the second parameter J is denoted by J(X) or J(Y).

References and Notes 1. Horowitz, S. B.; Fenichel, I. R. J. Phys. Chem. 1964, 68, 33783385.

Journal of Pharmaceutical Sciences / 777 Vol. 86, No. 7, July 1997

2. Gary-Bobo, C. M.; Weber, H. W. J. Phys. Chem. 1969, 73, 11551156. 3. Kier, L. B.; Cheng, C.-K. J. Chem. Inf. Comput. Sci. 1994, 34, 647-652. 4. Kier, L. B.; Cheng, C.-K. J. Chem. Inf. Comput. Sci. 1994, 34, 1334-1337. 5. Kier, L. B.; Cheng, C.-K.; Testa, B.; Carrupt, P.-A. Pharm. Res. 1995, 12 615-620. 6. Kier, L. B.; Cheng, C.-K. Pharm. Res. 1995, 12, 1521-1525. 7. Cheng, C.-K.; Kier, L. B. J. Chem. Inf. Comput. Sci. 1995, 35, 1054-1059. 8. Kier, L. B.; Cheng, C.-K. Pharm. Res. 1996, 13, 1418-1421. 9. Von Neumann, J. In Theory of Self-Reproducing Automata; Burks, A., Ed.; University of Illinois: Urbana, IL, 1966.

778 / Journal of Pharmaceutical Sciences Vol. 86, No. 7, July 1997

10. Ermentrout, F. B.; Edelstein-Keshet, L. J. J. Theor. Biol. 1993, 160, 97-133. 11. Einstein, A. Elecktrochemistry 1908, 14, 235-239. 12. Kier, L. B.; Testa, B. Adv. Drug Res. 1995, 26, 1 - 38.

Acknowledgments The senior author thanks the Fondation Herbette (University of Lausanne) for support and the University of Lausanne for the position of Visiting Professor. The authors thank Matthew Dowd for technical assistance. The calculations were carried out with the program DING-HAO.

JS9700513