INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

The classic Michaelis-Menten formalism for enzyme reactions (Michaelis and Menten, 1913; Segal, 1959) is based on mass-action laws, that primarily originate from Smoluchowski's approach of diffusion-limited reactions (Smoluchowski, 1917). Mass-action laws are mean-field approximations because they evaluate local reaction rates on the basis of average values of the reactant density over a large spatial domain. They rely on strict assumptions concerning, for instance, the characteristics of the reaction medium, which must be dilute, perfectly-mixed, three-dimensional, and homogenous. Moreover, the distribution of the reaction times must be exponential, i.e., of the Poisson type (Xie and Lu, 1999). Many of these assumptions fail in the case of biological reactions. Because of the complexity of biological molecules or molecular assemblies, reaction-time distributions may be broader than the Poisson one (Frauenfelder et al., 1999). Cellular media are not homogeneous, but highly compartmented. Moreover, they are not dilute solutions, because local concentrations are high. For instance, viscosity in the mitochondrion matrix is 25-37 times higher than that of a classical experimental buffer (Scalettar et al., 1991). High molecular crowding has important consequences as far as thermodynamics are concerned (Minton, 1993, 1998), but strongly influences diffusion processes as well. In the cytoplasm of fibroblasts, endogenous obstacles hinder the diffusion of tracer objects to such a point that objects with radius greater than 260 Å are immobilized (Luby-Phelps et al., 1987). Biomacromolecule diffusion coefficients in cytoplasm are usually 5-20 times lower than their values in saline (for a review, see Verkman, 2002). Thus biological media can be considered as disordered ones. Accordingly, measurements of diffusion of molecules into living cells have evidenced nonclassical behaviors (Schwille et al., 1999; Wachsmuth et al., 2000).

Furthermore, many cellular reactions occur in two-dimensional (2D) membranes. This is an important aspect to take into account, because diffusion is highly dependant on the Euclidean dimension d of the medium in which it occurs. In dimension d > 2, only a low fraction of the accessible volume is explored by a diffusing molecule, so that the molecule always escapes its initial position (noncompact diffusion). Whereas, in dimension d < 2, the molecule mainly remains in the neighborhood of its initial position (compact diffusion) (de Gennes, 1982). In other words, diffusion is not a perfectly mixing process in low dimension (d < 2) because the diffusing molecule will eventually return to its initial position with probability 1, whereas, for d > 2, there is a significant probability that the diffusing molecule will never return to its origin (Montroll and Weiss, 1965). This has important consequences on the mean-squared displacement of the molecule, which scales with time as R²  t. Whereas the exponent  = 1 for conventional diffusion,  < 1 in many low-dimensional media. In this case, diffusion is called anomalous.

The case of diffusion on fractal media such as percolation clusters is especially relevant to this issue because the effective dimension of fractal objects is less than the Euclidean dimension of the medium in which they are embedded (de Gennes, 1983; Havlin and Ben-Avraham, 1987). In contrast, percolation clusters are useful models of disordered medium in which diffusion is hindered by immobile obstacles (Saxton, 1994). Diffusion is anomalous both in low-dimensional homogenous media and on percolation clusters (Bouchaud and Georges, 1990). As a result of this abnormality, diffusion-dependant processes are altered in low dimensions. For example, diffusion-limited elementary chemical reactions of the type A + A   or A + B   in low-dimensional media are not adequately described by mass-action laws. The rate coefficient k relating the reaction rate r to A and B concentrations or densities (r = k_AB) is no more constant, but decreases at long reaction times as a power law of the reaction time: k  t^h (Kopelman, 1986). This, in turn, yields anomalous values of the reaction order X. For the homobimolecular diffusion-limited reaction, mass-action laws predict a reaction order X = 2, whereas the expected values are X = 2.5 for diffusion on 2D percolation cluster and X = 3 for diffusion in a one-dimensional homogenous medium (Kopelman, 1988). These nonconventional kinetics have been referred to as "fractal" kinetics. They arise from a spatial self organization of the reactants induced by the compact properties of diffusion (Argyrakis and Kopelman, 1990). In the A + B   case, the self organization can even lead to spontaneous segregation of the reactants into A-only and B-only regions, a phenomenon called the Zeldovich effect (Ovchinnikov and Zeldovich, 1978; Toussaint and Wilczek, 1983).

Enzyme kinetics are poorly understood when the reaction occurs in low-dimensional and disordered media such as biological ones. However, in vivo measurements have frequently evidenced deviations from traditional mass-action laws in the form of anomalous reaction orders and power-law reaction rates (Savageau, 1992, 1995). On a theoretical point of view, the understanding of enzyme reactions is hardly reducible to elementary bimolecular reactions. Basically, the conventional Michaelis-Menten scheme includes three elementary reaction steps. The abundant literature on elementary chemical reactions in low dimensions usually deals with diffusion-limited reactions occurring on percolation clusters at the threshold and initiated with equal reactant densities. Equal enzyme and substrate concentrations are poorly relevant to biological cases. Moreover, the percolation threshold is obtained at high obstructing obstacle densities, whereas, for biological reactions, the entire range of obstacle densities (from absence of obstacle to threshold density) must be considered. This indicates that the issue of enzyme kinetics in low-dimensional disordered media must be specifically addressed.

The purpose of this paper is thus to study Michaelis-Menten enzyme kinetics on low-dimensional lattices. Monte Carlo (MC) simulations are especially adapted to this kind of study because they allow the formulation of the reaction scheme as simple evolution rules. Furthermore, they provide representations of the spatial distribution of the molecules during the reaction, allowing direct imaging of the molecule repartition on the lattice. We thus use MC simulations to model the evolution of an isolated enzymatic system consisting of enzyme, substrate, product, and enzyme-substrate complex molecules diffusing on a 2D lattice while reacting according to the Michaelis-Menten reaction scheme. We also place immobile obstacles on the lattice to account for the high molecular crowding of biological membranes that hinders diffusion. The reaction medium can thus be continuously varied from homogeneous and 2D without obstacle to a fractal percolation cluster when obstacle density reaches the percolation threshold (Saxton, 1994).

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Monte Carlo simulations and model description

The Michaelis-Menten scheme of enzyme kinetics is a paradigmatic model in biochemistry. It consists in a set of three elementary chemical reactions

(1)

where E, S, P, and C represent enzyme, substrate, product and enzyme-substrate complex, respectively, and k_i is the rate coefficient associated with the elementary step i. According to the classical mass-action laws (mean-field kinetics), the k_i are constant throughout the reaction (Smoluchowski, 1917). The classic derivation of the equations describing reaction scheme 1 assumes mean-field chemical kinetics for each elementary step. This results in the set of ordinary differential equations,

(2)

(3)

(4)

where _i is the overall density of species i, and t is time. This is the classic textbook case. The set of Eqs. 2-4 is a singular perturbation problem: one asymptotic solution can be found in the neighborhood of t = 0 (boundary layer), and another for greater times, but there is no analytical solution valid over the entire course of the reaction (Murray, 1993). Biochemists usually use a quasi-stationary state assumption, which relies on the observation that, after an initial prestationary period, _C remains essentially constant during the rest of the reaction, provided that the ratio between initial E and S densities is low: _E(0)/_S(0) 1.

We simulated Reaction 1 using a MC algorithm on 2D square lattices with cyclic boundary conditions. Each molecule type (E, S, C, and P) is mobile on the lattice through diffusion, which is modeled by independent (nearest-neighbor) random walk of the individual molecules. When present, obstacles are represented as a fifth but immobile and nonreactive molecule type. The coordinates of the position of every molecule and the occupancy status of each lattice site are stored (Lin et al., 1996) and used for analysis. At any moment of the simulation, one given lattice site cannot be occupied by more than one molecule at a time (excluded volume condition). The rate coefficients k₁, k₁, and k₂ are modeled by the reaction probabilities f, r, and g, respectively (see below). At the beginning of each simulation, the E and S molecules and the obstacles (if present) are placed on the lattice by randomly choosing the coordinates for each of them. At each MC sequence, an occupied lattice site is chosen at random (excluding obstacle sites). The rules for movement and reaction of the associated molecule depend on the nature of this molecule, in agreement with Reaction 1:

1.	If the molecule occupying this site is an S, a destination site is chosen at random between its four nearest neighbors. If this destination site is unoccupied, the molecule moves to it directly. If the destination site is occupied by an E molecule, a random number is chosen between 0 and 1. If this number is lower than the reaction probability f, the destination site is turned to a C molecule and the initial S site becomes unoccupied. In all other cases, the S molecule remains at its initial position. Note that this is also valid if the chosen destination site is an obstacle (blind ant, see, for example, Majid et al., 1984).
2.	If the molecule occupying the chosen site is an E, the process is symmetric to the former case, i.e., depends on the occupancy status of the randomly-chosen destination site: movement if unoccupied, reaction with a probability f if occupied by S, or immobility in all other cases.
3.	If the molecule occupying the chosen site is a C, a random number is chosen between 0 and 1. If this number is lower than the reaction probability r, and provided that at least one of its nearest neighbors is unoccupied, the C molecule dissociates into an E and an S. The new E molecule is placed on the initial C site, whereas the new S molecule randomly moves to one of the unoccupied sites. Note that this step is not a blind ant process, contrary to the other steps. A more physically realistic way would be to choose a site at random for the new S molecule, move it to this site if unoccupied, and abort the decomposition process if occupied. However, we checked that this had no important consequences in our simulations, at least for the E, S, and obstacle-density ranges used. The initial C molecule dissociates into an E and a P in a similar fashion, if the random number is greater than r but lower than r + g. Finally, if the random number is greater than r + g, the C molecule is allowed to move to a randomly chosen unoccupied nearest-neighbor site.
4.	if the chosen molecule is a P, it moves to a randomly-chosen unoccupied nearest-neighbor site.

After each sequence, time is incremented by 1/, where  is the current number of molecules on the lattice (excluding obstacles), and another sequence begins. One time unit thus statistically represents the time necessary for each molecule to move once. The simulation goes on until a prescribed total time (600 time units in this study, unless specified). Lattice size was 100 × 100. Initial densities were always zero for C and P, and ranged from 0.002 to 0.02 for E, and from 0.01 to 0.2 for S. In this study, the values of the reaction probabilities f, r, and g were set to 1, 0.02, and 0.04, respectively, unless specified. The results presented are averages of 100 or 200 runs. Uniformly distributed random numbers were generated with the combined generator of L'Ecuyer and Cote (1991), as implemented in the RANLIB package (B. W. Brown and J. Lovato. Department of Biomathematics, The University of Texas, Houston). The algorithm was programmed with SCILAB (http://www-rocq.inria.fr/scilab/) and run on a beowulf cluster under Linux.

Data analysis

Throughout the simulation course, the overall density of each species (calculated over the whole lattice) is recorded. We also record (t), the total number of enzyme-substrate collisions that have effectively given rise to complex formation after time t. (t) is related to k₁ by

(5)

On the basis of the sampled (t) values, the time derivative d(t)/dt is estimated numerically after interpolation by third-order spline functions. Using enzyme and substrate densities at time t, k₁ can then be estimated from Eq. 5 by

(6)

The observed time-dependence of k₁ is accounted for in Eqs. 2-4, and the set of ordinary differential equations is integrated numerically in SCILAB with a solver suitable for stiff problems (Gear method, Isode solver of the ODEPACK package).

Reactant self-segregation is evaluated with a quantitative criteria Q_SP inspired by New house and Kopelman (1988), which is based on a quotient of effective pair-correlation functions

(7)

where N_SS is the number of S-S pairs, N_PP is the number of P-P pairs, and N_SP is the number of S-P pairs. A pair is defined here by two nearest-neighbor sites. N_SP for instance, is thus the number of site pairs of which one site is occupied by S and the other by P. A random distribution of the molecules over the lattice yields Q_SP = 1, whereas under S-P segregation, most of the S-P pairs are located on the interfaces between S-rich and P-rich domains, so that Q_SP > 1.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Anomalous diffusion

In lattices with immobile obstacle densities below the percolation threshold, the accessible sites (i.e., those that are not obstructed by an obstacle) form a percolation cluster (Saxton, 1994). Percolation clusters are fractal over distances shorter than the correlation length  and homogenous for larger ones. Diffusion on percolation clusters is thus expected to show a crossing-over from an anomalous regime (short times) to a normal one (longer times) as the distance explored by the diffusing molecule becomes greater than the correlation length. We first simulated the diffusion of a single S molecule on a 2D lattice with varying obstacle densities . Figure 1 shows the mean-squared displacement R² as a function of time. The observed power-law time dependence can be expressed by R²  t, where the anomalous diffusion exponent  is the slope of curves in log-log coordinates. Note that, in the physics literature, the d_w notation, where  = 2/d_w, is often used. Without obstacle, diffusion is normal over the entire time scale, i.e.  = 1 for  = 0. With increasing obstacle densities, diffusion becomes increasingly anomalous (i.e.,  decreases with increasing obstacle densities). We note that, with nonzero obstacle densities, diffusion remains anomalous over the entire time scale (10³ time units = 10³ random-walk steps). Taken together, these results are in agreement with previously published simulations (Saxton, 1994).

View larger version (17K): [in this window] [in a new window]   FIGURE 1   Diffusion of a single molecule on 2D lattices. A single molecule is deposited in the lattice at a randomly chosen site at t = 0, and the simulation proceeds as described in Methods. The mean squared displacement R2 of the molecule is plotted as a function of time, in log-log coordinates. The slope of the curves is , the anomalous diffusion exponent. Obstacle densities  are (from top to bottom) 0, 0.10, 0.15, 0.25, 0.35, and 0.40. Indicated are the values of /c for the two extreme densities, where c is the obstacle density at the percolation threshold (c = 0.4073 for such site-percolation problems, Sahimi (1994)). For comparison, the dotted line indicates a slope of 1. One time unit corresponds to one move of the random walker between two nearest-neighbor sites.

Fractal kinetics

To probe the influence of anomalous diffusion on Michaelis-Menten enzyme kinetics, we carried out MC simulations on 2D lattices (see Methods). Figure 2 shows the evaluation of the rate coefficients k₁ and k₂ (cf. Reaction 1) during the simulations, for different /_c values. k₁ and k₂ are evaluated from Eqs. 2-4 and 6 by

(8)

and

(9)

where time derivatives are evaluated numerically as described for d(t)/dt in Methods. As seen in Fig. 2 B, these coefficients remain unchanged throughout the simulation and do not depend on obstacle density, in agreement with classic chemical kinetics. Furthermore, because these constants are first-order, they are expected to correspond exactly to the reaction probabilities r and g of the simulations (see Methods). Accordingly, in Fig. 2 B, one obtains on average k₁ = r = 0.02 and k₂ = g = 0.04. The situation is different with the rate coefficient k₁ because k₁ is a second-order rate coefficient, and thus accounts for diffusion effects: to form an enzyme-substrate complex, one enzyme molecule and one substrate molecule must first collide in the course of their respective random walks. Figure 2 A shows the value of k₁ for different /_c values. Unlike k₁ or k₂ and unlike the prediction of classic Michaelis-Menten kinetics, k₁ is not constant throughout the reaction course. It crosses over from a constant regime at short times to a power law decrease at longer ones. This behavior is the hallmark of fractal kinetics:

(10)

The exponent h, the slope of the curves at long times in Fig. 2 A is the fractal kinetics exponent (Kopelman, 1988). Its value is zero for classic kinetics (short times in Fig. 2 A) and >0 in the anomalous regime (long times in Fig. 2 A). Note that one time unit in these simulations statistically corresponds to the time necessary for each molecule to move once. This is the same time scale as in Fig. 1. Diffusion is thus anomalous throughout the enzyme kinetics.

View larger version (28K): [in this window] [in a new window]   FIGURE 2   Rate coefficient time dependence with varying obstacle densities on 2D lattices. (A) k1 (t) and (B) k1 (t) and k2(t) are evaluated as described by Eqs. 6, 8, and 9, respectively. Initial enzyme and substrate densities are E(0) = 0.01 and S(0) = 0.10. The relative obstacle densities /c are 0 (curve 1), 0.491 (curve 2), 0.737 (curve 3), 0.859 (curve 4), and 0.999 (curve 5). In (A) the slope of each curve in the long-time regime is the value of the fractal kinetics exponent h that describes k1 anomalous time dependence. In (B), the curves for different obstacle densities are all identical to within statistical scatter, whether for k1(t) or k2(t).

Because of the anomalous behavior of the complex formation step, the entire reaction kinetics are anomalous as well. Figure 3 shows the evolution of the overall enzyme and complex densities during 2D simulations (curve MC) for two obstacle densities (/_c = 0.908 in Fig. 3, A and B, and 0.368 in Fig. 3, C and D). For comparison, the dotted lines indicate the results of numerical integration of the classical Michaelis-Menten equations, Eqs. 2-4 (curve CK). Except for the very short time regime, the results of the simulations (MC) are clearly not adequately described by the classic formalism (CK), whether for C (Fig. 3, A and C) or S densities (Fig. 3, B and D). This behavior is in agreement with k₁ evolution, which is almost normal in the short time regime, and anomalous for longer times (Fig. 2 A). E and P densities can be obtained from Fig. 3 using mass-conservation principles, i.e., _E(t) = _E(0)  _C(t), and _P(t) = _S(0)  _C(t)  _S(t). For this reason, the time evolution of these densities is not presented in Fig. 3. However, it can be deduced from these conservation laws that, like _S and _C, _E and _P kinetics do not comply with classical Michaelis-Menten equations. As mentioned above, the abnormality of diffusion and k₁ time dependence increases with increasing obstacle densities. In agreement with this behavior, the abnormality of the kinetics increases with obstacle densities (compare Fig. 3, A and B with Fig. 3, C and D). Consequently, the discrepancy between simulation results (MC curves) and the kinetics predicted by the classic equations (CK curves) is higher with high obstacle densities.

View larger version (24K): [in this window] [in a new window]   FIGURE 3   Kinetics of (A and C) enzyme-substrate complex and (B and D) substrate densities during Monte Carlo simulation of a Michaelis-Menten enzyme reaction on a 2D lattice, with  = 0.37 (A and B) or  = 0.157 (C and D). The curves labeled MC are the direct results of Monte Carlo simulations. The curves labeled CK (Classic Kinetics) are the results of numerical integration of the classic equations for enzyme reaction (Eqs. 2-4). The curves labeled FK (Fractal Kinetics) are the results of numerical integration of Eqs. 11-13 with k01 = 1.035, h = 0.549 (A and B) or k01 = 1.191, h = 0.347 (C and D). In each panel, initial densities are E(0) = 0.01 and S(0) = 0.20. k01 and h were determined from log-log plots of k1 time dependence, such as those presented in Fig. 2 A. The other parameters were set to k1 = r = 0.02 and k2 = g = 0.04, in agreement with Fig. 2 B.

An important consequence of these fractal kinetics concerns the so-called quasi-stationary-state hypothesis. This hypothesis states that the enzyme-substrate complex density is quasi-constant during the reaction and provides the biochemist with a useful simplification by canceling d_S/dt in Eq. 3. This, in turn, allows the reduction of Eqs. 2-4 to a single differential equation for the initial velocity of S consumption or P production. It can be shown that this hypothesis is increasingly valid when the initial substrate-to-enzyme density ratio _S(0)/_E(0) is increasingly high (Murray, 1993). For calculation time reasons, such high _S(0)/_E(0) values are difficult to study in MC simulations. In Fig. 3 A, the initial ratio _S(0)/_E(0) (= 20) is not sufficient to obtain a real quasi-stationary state for _C, even in the classic case (curve CK). However, it is obvious from this figure that 2D kinetics (MC results) are even more different from steady-state than those predicted by the classic equations. This indicates that the validity of the quasi-stationary state hypothesis in 2D enzyme kinetics is largely questionable, because the _S(0)/_E(0) values required are much higher than those predicted by classic Michaelis-Menten kinetics.

The classic Michaelis-Menten equations (Eqs. 2-4) can be modified to account for k₁ time dependence during fractal kinetics (Eq. 10). This yields the following equation set for long reaction times:

(11)

(12)

(13)

The classic Michaelis-Menten formalism is thus a special case of Eqs. 11-13, that corresponds to h = 0. To validate the capacity of these equations to describe 2D kinetics, we integrated them numerically in the anomalous time regime. For each simulation condition, h and k were determined by nonlinear least-square fit to the values of k₁ in the long-time domain (Fig. 2 A), according to Eq. 10. The results are presented in Fig. 3 (FK curves). In each case, the densities predicted by Eqs. 11-13 are exactly superimposed to the simulation results as far as _S is concerned (Fig. 3, B and D). Predicted _C kinetics are also in very good agreement with the simulation results (Fig. 3, A and C).

Reactant segregation

We were then interested in the spatial distribution of the E, S, C, and P molecules during 2D enzyme kinetics. In our simulations, the E, S, and obstacle distributions at t = 0 are spatially homogeneous over the lattice. Figure 4 shows the molecule distribution at longer times for increasing obstacle densities. The reaction time in each panel of Fig. 4 corresponds approximately to the mid-reaction point in each simulation. From visual inspection of this figure, S-P segregation is not obvious in the unobstructed case (Fig. 4 A). The evolution of the quantitative segregation criteria Q_SP during the reaction confirms this observation (Fig. 5). Although a random distribution yields Q_SP = 1, one expects Q_SP > 1 under S-P segregation. With low obstacle densities, Q_SP rapidly tends to a steady value. This value is slightly greater than 1 without obstacles (Q_SP = 1.18 ± 0.02 for /_c = 0). This slight deviation from 1 indicates that S-P segregation is weak in the unobstructed case. With increasing obstacle densities, S and P molecule segregation into S-rich and P-rich regions is clearly visible from Fig. 4, C and D. In these cases, Q_SP accordingly increases and is greater than 1 throughout the reaction (Fig. 5). The maximal values reached at t = 600 are 2.3 and 2.8 for /_c = 0.908 and 0.998, respectively. These values confirm the possibility of occurrence of a Zeldovich segregated regime in 2D enzyme kinetics. Note however, that segregation is here less dramatic than in the case of the A + B   reaction with equal initial densities, where Q_SP  10 in the Zeldovich regime (Newhouse and Kopelman, 1988).

View larger version (134K): [in this window] [in a new window]   FIGURE 4   Spatial distribution of the E (×), S (), C (+) and P () molecules on the 2D lattice for single-simulation realizations with homogeneous initial distributions. Simulation parameters are  = 0.000 (A), 0.150 (B), 0.250 (C), or 0.406 (D) and E(0) = 0.01, S(0) = 0.20 in each case. Snapshots were realized at t = 300 (A, B and C) or 600 (D) to obtain P and S values corresponding approximately to the mid-reaction point. Note that the symbols are not to scale, so that apparent densities are greater than actual ones. N.B.: For clarity, immobile obstacles are not represented in these figures.

View larger version (20K): [in this window] [in a new window]   FIGURE 5   Quantification of S-P molecules segregation during enzyme reaction on a 2D lattice with obstacles. Indicated above each curve is the corresponding relative obstacle density. QSP is calculated as described in Methods (Eq. 7). E(0) = 0.01, S(0) = 0.20 in each case.

h Variations with obstacle and substrate densities

As already mentioned, the fractal kinetics exponent h increases with increasing obstacle densities. However, h appears to be dependent on the initial substrate density as well. Figure 6 presents h dependence on obstacle density for different initial substrate densities. Whatever the initial _S value, h increases with obstacle density up to a maximal value at the percolation threshold. The value at the threshold (h( = _c) = 0.632 ± 0.024 in Fig. 6) does not depend on the initial substrate density. In contrast, for obstacle densities below the threshold, h increases with _S(0). As seen in this figure, the kinetics are anomalous even without immobile obstacles, i.e., on homogeneous 2D lattices, although diffusion is normal in this case. We further studied the variations of h on homogeneous lattices (h( = 0)), with two enzyme initial densities (Fig. 7). The data with both enzyme densities (_E(0) = 0.01 or 0.005) collapse on a single line in log-log coordinates, yielding a power-law behavior for h( = 0)

(14)

From Fig. 7, the prefactor  and the exponent  were evaluated to 0.6 and 0.43, respectively. Note that the enzyme densities are equally low as compared to most of the initial substrate densities in Fig. 7, so that a weak dependance of h( = 0) on _E(0) is not excluded. We tried to obtain the influence of the obstacle density only on the anomalous kinetics by subtracting the corresponding h( = 0) value form the h() values presented in Fig. 6. When plotted as h()  h( = 0) as a function of , the data of Fig. 6 mostly collapse on a single curve (Fig. 8). These results thus indicate that the respective influence of immobile obstacle density and substrate densities on Michaelis-Menten fractal kinetics in two dimensions are mainly additives, i.e.,

(15)

As can be seen from Fig. 8, this approximation partly fails when obstacle density approaches the percolation threshold. This could be due to the poor precision in the values of  and  estimated from Fig. 7. Alternatively, Eq. 15 may need corrections close to the threshold.

View larger version (17K): [in this window] [in a new window]   FIGURE 6   Fractal kinetics exponent h as a function of obstacle and initial substrate densities. h values were determined from log-log plots of k1 time-dependence, such as in Fig. 2 A. Initial substrate densities are: S(0) = 0.20 (), 0.1 (), 0.05 () and 0.01 (). In each case, E(0) = 0.01.

View larger version (12K): [in this window] [in a new window]   FIGURE 7   Fractal kinetics exponent h( = 0) for Michaelis-Menten kinetics on 2D lattice in the absence of immobile obstacles as a function of initial substrate density. Initial enzyme densities are E(0) = 0.02 (), 0.01 (), 0.005 (+) or 0.002 (×). The full line is a nonlinear least-square fit on the four data series, yielding h( = 0) = 0.60(S(0))0.43.

View larger version (13K): [in this window] [in a new window]   FIGURE 8   Plot of h()  h( = 0) for Michaelis-Menten kinetics on 2D lattices as a function of the relative obstacle density. Data are replotted from Fig. 6 after h( = 0) subtraction from the overall h() values presented in Fig.  6. h( = 0) is calculated for each initial enzyme density from Eq. 14. The symbols and initial enzyme and substrate densities are as in Fig. 6

The influence of the initial substrate density could be due to the excluded volume condition included in our simulations. Actually, because of this condition, each molecule on the lattice can be considered as a mobile obstacle with respect to another molecule. The importance of the excluded volume increases with _E(0) + _S(0), which could partly account for the observed h( = 0) variations. We carried out a modified version of the MC algorithm, in which all lattice sites (except obstacle ones) can be occupied by several molecules at a time. The simulations show that the h values with or without excluded volume condition are largely comparable (not shown). Thus the excluded-volume condition is not likely the main cause of the substrate density influence on our kinetics.

Finally, the results presented above assume that C formation from E + S is diffusion limited, as reaction probability f = 1. We studied the influence of partial limitation by the reaction (i.e., f < 1) of this part of the overall Michaelis-Menten scheme. Figure 9 shows k₁ for varying f values, with obstacle density  = 0.20, _E(0) = 0.01, and _S(0) = 0.10. The crossover time to the anomalous regime increases as f decreases. More importantly, h, the slope of the curves in the anomalous regime decreases at lower f values (its value at f = 0.30 is 55% of its value at f = 1.00). This indicates that the abnormality of the kinetics decreases as the E + S  C step becomes less diffusion controlled.

View larger version (17K): [in this window] [in a new window]   FIGURE 9   Evolution of k1(t)/k1(0) as a function of time for three values of the E + S  C reaction probability f. The value of f for each simulation condition is indicated on the figure. Reaction probabilities r and g are 0.02 and 0.04, respectively. Initial densities are S(0) =  = 0.20, and E(0) = 0.01

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Mass-action laws and enzyme reactions in low dimensions

Dimensionality plays a crucial role in many diffusion-dependant mechanisms in biology. Target-search processes such as specific DNA sequence search by transcription factors in cell nucleus, or receptor search by pheromones on insect's antenna, are dramatically accelerated when diffusion is at least partly two-dimensional (Adam and Delbrück, 1968; Holyst et al., 2000). Furthermore, conventional mean-field equations or mass-action laws are continuum approximations that can fail to describe diffusion-reaction processes in low dimensions, because they do not take spatial fluctuations into account. In contrast, spatial density fluctuations are inherent to discrete lattice simulations. Several recent studies have evidenced significant discrepancies between the predictions given by these two approaches. Strikingly, discrete approaches in low dimension frequently predict self organizations that are not predicted by the corresponding mean-field equations. Two-dimensional simulations of simple population-growth models show species survival on localized lattice regions where mean-field equations predict extinction of the entire population (Shnerb et al., 2000; Young et al., 2001). Discrepancies are also observed in more elaborate prey-predator models (Lipowski and Lipowska, 2000). For instance, in the case of the classic Lotka-Volterra prey-predator model, the elliptic fixed-point predicted by the mean-field equations is changed to a limit cycle in 2D simulations (Bettelheim et al., 2001).

The issue of enzyme reaction in low-dimensional media has rarely been addressed. Savageau (1995) examined the possible implications of anomalous reaction orders on enzyme kinetics under the quasi-stationary-state assumption. Zhdanov (2000) used MC simulations on 2D lattices to study an enzyme reaction with k₁ = 0, immobile enzyme molecules, and instantaneous complex decomposition. The reaction scheme is then E + S  E + P with immobile enzyme and unequal initial densities and is kinetically analogous to a trapping reaction with the enzyme as the trap. This study shows that, in the case of attractive interactions between the product molecules, phase separation between S and P molecules occurs provided that the rate constant is fast. Lin et al. (1997) studied experimentally the glucose oxidase enzyme reaction with immobilized enzyme and mobile substrate. Again, this reduces to a trapping reaction. They indeed show that, in a capillary (pseudo-1D medium), the reaction is fractal and in good agreement with the behavior predicted for d = 1.

In this paper, we studied the native Michaelis-Menten reaction scheme using MC simulations on 2D lattices with immobile obstacle densities varying from zero to the percolation threshold. We show that, as a result of the anomalous diffusion on these low-dimensional media, the conventional mean-field description fails because the reaction kinetics are of the fractal type (Kopelman, 1988). Of importance, we show that the conditions necessary to fulfill the quasi-stationary-state assumption are much more drastic than in the conventional case, so that this assumption is highly questionable. The differential equation set describing the kinetics can thus hardly be reduced to a single equation for the initial reaction velocity. This could be a significant indication for experimental enzymology, where the determination of the kinetic parameters is usually based on the quasi-stationary-state assumption derived rate law. In the absence of quasi-steady state, the kinetic parameters can still be obtained by a fit to the results of the numerical integration of the differential equation set (Eqs. 11-13). Furthermore, we show that the fractal characteristics of the kinetics increase with obstacle density and initial substrate concentration, the two contributions being mainly additive. Finally, as in diffusion-limited elementary chemical reactions, fractal kinetics seem to originate from reactant self organization that leads at relatively high obstacle densities to S-P segregation over the lattice.

Comparison with diffusion-limited heterobimolecular reactions

The diffusion-limited A + B   reaction has been widely studied both by simulation and theory, whether on homogeneous 2D lattices (Kang and Redner, 1984, 1985; Lindenberg et al., 1988; Ovchinnikov and Zeldovich, 1978; Toussaint and Wilczek, 1983) or on percolation clusters at the threshold (Anacker and Kopelman, 1987; Argyrakis and Kopelman, 1992; Kopelman, 1986; Newhouse and Kopelman, 1988). For simplicity reasons, these studies usually use equal initial A and B densities. In the case of diffusion on a fractal support, the reaction constant k behaves as k  t^h, with h = 1  d_s/2 before the onset of the Zeldovich regime, and h = 1  d_s/4 afterward (Argyrakis et al., 1993; Lin et al., 1996). d_s is called the spectral dimension and equals  × d_f where d_f is the fractal dimension of the support, and  is the anomalous exponent for the diffusion on the fractal. The theoretically expected h value for the diffusion-limited A + B   reaction on a 2D percolation cluster at the threshold is h =  in the segregated regime (Havlin and Ben-Avraham, 1987). In our simulations with equal initial densities, h = 0.6662 at the threshold (Fig. 6), in very good agreement with the predicted value.

To compare with our results for obstacle densities below the threshold, a major problem consists in the value of the spectral dimension as a function of obstacle density. Argyrakis and Kopelman (1984) have measured an effective spectral dimension d'_s from the average number of distinct sites visited, during random walks on lattices with obstacle densities varying from zero to the threshold. A nonlinear least-square fit of their data for d'_s (see Fig. 6 in their paper) yields d'_s = (7.57x² + 6.91x + 1.80)/(3.89x² + 3.81x + 1), where x = /_c. This approximation allows correlation of the h values with d'_s for different obstacle densities. As seen in Figs. 4 and 5, S-P segregation is weak in the unobstructed case. We thus expect h to vary from 1  d'_s/2 at low obstacle densities to 1  d'_s/4 close to the threshold. Figure 10 shows a test of these two expressions, on the basis of the h values determined in Fig. 6 and the calculated values for d'_s. The h values observed in our simulations indeed vary between h = 1  d'_s/2 at zero obstacle densities to h = 1  d'_s/4 close to the threshold.

View larger version (13K): [in this window] [in a new window]   FIGURE 10   Plots of h(, S(0)) + d'S/2 () and h(, S(0)) + d'S/4 () as a function of the relative obstacle density. Initial substrate and enzyme densities are equal: E(0) = S(0) = 0.01. The d'S values were deduced from Fig. 6 of Argyrakis and Kopelman (1984) (see text). The doted line indicates the value 1.

In contrast, the Michaelis-Menten reaction scheme is not an elementary hetero-bimolecular reaction, because the two monomolecular C decomposition reactions complicate the kinetics. These steps mainly introduce an effective mixing because E (and S) molecules are regenerated after a E + S  C reaction, possibly at a certain distance of the initial E-S meeting. In A + B   reactions, reactant segregation occurs because, on account of its compact property, the random walk in low dimension is not "mixing" enough to dissipate the formation of the A and B regions. Introducing exogenous mixing thus tends to prevent the occurrence of the Zeldovich regime. This is, for example, the case when the random walk is replaced by Lévy flights (Zumofen et al., 1996) or when A and B molecules are exogenously added to the lattice during the reaction (stationary-state conditions, Anacker and Kopelman, 1987; Argyrakis and Kopelman, 1992; Lin et al., 1996). In the latter case, segregation occurs for Euclidean dimension d  2 (for homogeneous lattices), whereas it can be observed for d  4 in "batch" conditions. In dimension d = 2, segregation is marginal under stationary-state conditions (Lindenberg et al., 1988). Thus, the mixing effect of the two monomolecular C decomposition reactions in Michaelis-Menten kinetics could explain the weak segregation on 2D lattices without obstacle. Alternatively, the crossover time between nonsegregated and segregated regimes could be greater than the time necessary to complete the reaction.

The decrease of k₁ with time in the unobstructed case can be explained by two alternative mechanisms. A first possibility is that reactant segregation in the unobstructed case is low, but not zero. In the A + B   diffusion-limited case under steady-state conditions, segregation in two dimensions depends sensitively on the detailed parameter values (Lindenberg et al., 1988). Accordingly, in a model for protein G activation by mobile membrane receptors, Shea et al. (1997) evidenced inhomogeneous reactant distribution at steady state, together with reduced rate constants as compared to the homogeneous case. In this case, one would thus expect a fractal power-law behavior for k₁, as assumed in Eq. 14. In contrast, the rate coefficient for a diffusion-limited A + B   reaction in 2D unobstructed systems is known to scale asymptotically as 1/ln t (Torney and McConnell, 1983). In our simulations, the time dependence of k₁ at  = 0 can as well be fitted to a 1/ln t law (not shown). Because of the statistical error, we could not discriminate between these two behaviors for k₁ in the unobstructed case. Eq. 14 is thus to be considered as an empirical law, and has the advantage of preserving the same functional form for zero and nonzero obstacle densities.

Fractal enzyme kinetics and spatial self-organization in vivo

The appearance of stable forms or patterns from previously homogeneous spatial conditions is a central issue in biology. The physical mechanisms underlying these symmetry breakings are still largely unknown. One possible mechanism has been uncovered by Turing (1952). Within the framework of this theory, spatial self organization results from the growth of spatial fluctuations of concentration, which can be observed in some reaction-diffusion systems (Murray, 1993). Nevertheless, this approach necessitates specific reaction or diffusion conditions and cannot account for all spatial self-organization events in biology.

Our simulations evidence that spatial self organization in membranes could occur through simple Michaelis-Menten enzyme kinetics. However, the relevance of these simulation results to real in vivo situations must be carefully interpreted. First, we studied an isolated enzyme reaction, independent of external enzyme regulations or molecule influx or efflux through metabolic pathways that could as well contribute to effective mixing. Second, each molecule in our simulations moves with the same diffusion coefficient. This could be a crude approximation for some biological systems, where enzyme and substrate sizes are very different. A direct evaluation of obstacle densities in biological membranes is a difficult problem because the membrane area occupied by a transmembrane protein obstacle is usually unknown, and because there exist several mechanisms hindering diffusion (for a review, see Saxton, 1999). Estimations for several biological membranes give lower limits for obstacle-area fractions ranging from 0.2 to 0.4 (Saxton, 1989). Wachsmuth et al. (2000) measured the diffusion of protein probes in nuclei by spatially resolved fluorescence-correlation spectroscopy, and evidenced a slightly anomalous diffusion, with   0.87. For 3D percolation clusters at the threshold,  = 0.541 (Havlin and Ben-Avraham, 1987), so that obstacle density in the nuclei would rather be far from the threshold. Alternatively, Schwille et al. (1999) showed that diffusion in cell membranes is anomalous, with   0.741, which is close to its value at the percolation threshold in two dimensions ( = 0.697). Our simulations show visible reactant segregation for obstacle densities as low as 0.25 (Fig. 4). Thus we think that our simulation conditions are compatible with in vivo conditions. Finally, the value of the fractal kinetics exponent h seems to decrease when the C formation step is decreasingly diffusion controlled, or increasingly reaction controlled (Fig. 9). Thus the occurrence of fractal kinetics and spatial segregation could be restricted to enzymes with high catalytic efficiency, or at least high affinity for their substrate. Simulations with lower f values than those presented in Fig. 9 demand extensive calculations because the overall reaction is accordingly slowed down but could more precisely state this issue by permitting a thorough study of h dependance upon f values.