Unravelling the impact of obstacles in diffusion and kinetics of an enzyme catalysed reaction

Abstract

Lattice gas automata simulations of diffusion-limited reactions in heterogeneous media exhibit fractal-like kinetics, which is a generalised mass action kinetics with time-dependent rate constants. We develop a two dimensional lattice gas automata simulation of the Michaelis–Menten mechanism in diffusion-limited conditions to investigate the effect of density and size of obstacles on reactant diffusion and rate coefficients. In order to simulate more physicochemical realistic conditions, reactants rotate and interact according to their specific orientation. We also model weak interaction forces between reactants and obstacles. Our results show that obstacle density and size affect diffusion, first- and second-order rates. We also find that particle rotations and weak force interactions among particles lead to a significant decay in the fractal-like kinetic exponent h. These results suggest that the effects of fractal-like kinetics disappear under less restricted conditions than previously believed in lattice based simulations.

---

1. Introduction

The law of mass action states that the velocity of a reaction is proportional to the product of the concentrations of the participating molecules.¹ Although extensively applied to model biochemical reactions inside cells, the law of mass action assumes that reactions occur in homogeneous, well stirred volumes.² This is inconsistent with the heterogeneous environment that characterizes intracellular reactions.³ Intracellular environments possess a high macromolecular content, also known as macromolecular crowding, in the range of 5% to 40% of the total cellular volume.⁴ Macromolecular crowding qualitatively and quantitatively affects molecule diffusion and reaction kinetics.^3–6

The role of crowding in diffusion has been extensively investigated.^7–17 Theoretical and computational studies show that in two-dimensional (2D) media, diffusion is characterized by three main distinct time periods: a very fast initial period, where diffusion is regular i.e., constant; a fast transient period, where diffusion is anomalous, i.e., not constant; and a final period where diffusion is regular again.¹¹ These three phases have been shown to exist in the diffusion of Cajal bodies in the nucleus of HeLa cells¹⁸ and in the diffusion of gold-labeled dioleoyIPE in the plasma membrane of fetal rat skin keratinocyte cells.¹⁹ Tracer-obstacle simulations have also shown that the anomalous diffusion slope and occurring period increase with increasing obstacle density.⁷ More recently, Vilaseca¹⁶ show that at long-times the diffusion coefficient increases with the size of the obstacles.

Since diffusion directly relates to the probability of encounters between two molecules in a finite space, diffusion is important in understanding reaction kinetics, particularly in crowded media. Theoretical and experimental work show that diffusion-limited elementary reactions depend on the space in which they occur.^20,21 Notably, Kopelman^22,23 shows that the law of mass action breaks down in heterogeneous crowded environments, where in the presence of obstacles or obstructions, the second-order rate decays linearly on a logarithmic time scale as time tends to infinity. Kopelman interprets the rate of this decay as a measure of the dimensionality of the system, where the higher the decay, the lower the dimensionality. This is the hallmark of fractal-like kinetics.^3,24

Numerous computer simulations, predominantly using lattice gas automata, have been implemented to investigate the role of obstacles in reaction kinetics.^{3,7,15,16,24,25} While these studies appear to support Kopelman's theory, lattice approaches restrict a particle's degree of movement and hence, may affect diffusion and reaction kinetics. In 2010, using an off-lattice model, Grima²⁶ showed that deviations from diffusional isotropy significantly reduce the diffusion coefficient of tracer particles in crowded environments. Also in 2010, McGuffee and Elcock²⁷ developed an off-lattice Brownian dynamic simulation that includes both rotational and translational diffusion of the fifty most abundant cytoplasmic E. coli proteins. Their work provides some insights on diffusion, association and disassociation rates, but does not explicitly consider reaction kinetics of proteins. The effects of obstacle size and density in reaction kinetics have also been studied using a modified chemical master equation and an off-lattice Brownian model.²⁸ However, many of these studies are under reaction-limited conditions, where reactants interact with low probability.

Here, we investigate the role of obstacle density and size on the kinetics of an enzyme-catalysed reaction in diffusion-limited conditions, where second-order rates are expected to be time-dependent and follow fractal-like kinetics. For this purpose, we develop an enhanced 2D lattice simulation with cyclic boundary conditions of the Michaelis–Menten (MM) mechanism of enzyme action. Our enhanced model of the lattice simulation explicitly models a particle's rotation and weak force interactions among particles in the lattice.

2. Theory

2.1. Diffusion

The well-known Einstein–Smoluchowski equation describes the relationship between the diffusion D and the distance r taken by a particle at time t:

〈r²(t)〉 = (2d)Dt (1)

where d is the topological dimension of the medium surrounding the particle. If diffusion is constant, the mean square displacement of a particle is proportional to time. However, in a crowded medium, the diffusion of a particle may be hindered by large macromolecules. In this case, the diffusion process may have a non-linear relationship with respect to time, and eqn (1) is generalized to:²⁹

〈r²(t)〉 = (2d)Γt^α (2)

where Γ is the anomalous diffusion coefficient, α is the anomalous diffusion exponent (0 < α < 1) and D(t) = Γt^α. Using l as the unit of length in the lattice, ω as the unit of time, θ as the density of obstacles, and D₀ as the diffusion coefficient of a particle in an unobstructed system (θ = 0), we simplify eqn (1) by making r = lr*, t = ωt*, D = D₀D*(θ) and l² = 2dD₀ω. We then obtain the following dimensionless equation:⁷

〈r*²〉 = D*(θ)t* (3)

Note that D*(0) = 1. In addition, notice that the term ‘2d’ disappears in the simplified equation and that the lattice model does not describe dynamics on timescales faster than ω. To simplify our notation, we drop the stars of eqn (3) throughout the paper.

Using a log–log plot of diffusion over time, diffusion is characterised by three distinct parameters: the anomalous diffusion exponent (α); the limit diffusion (D*); and the crossover time (t_CR). The slope α − 1 of the initial time steps is used to calculate α. D* is the value of diffusion as time approaches infinity. t_CR is the time required to crossover from anomalous to regular diffusion by intersecting the slope α − 1 of anomalous diffusion and the constant slope of regular diffusion at long times. If α = 1, diffusion is constant at all times and no crossover exists.⁷

2.2. The Michaelis–Menten mechanism and fractal-like kinetics

In the MM mechanism of enzyme action, an enzyme (E) and a substrate (S) bind reversibly to form a complex (C), which irreversibly yields a product (P):where k₁, k₂ and k₃ are rate constants. When few obstacles are present in the reaction media, k₁ remains constant over time in agreement with the law of mass action. However, in diffusion-limited conditions, the time required for any two reactants to interact increases and k₁ decays over time.^22,23 As a consequence, log (k₁) decays linearly at long times in a logarithmic time scale:

log k₁ = log k₀ − h log t (5)

where k₀ is the rate at t = 0 and h is the fractal-like kinetic parameter. Notice that when h is 0, k₁ is no longer time-dependent. Therefore, h can be used to measure deviations from the law of mass action.²⁴ Note that h is only greater than 0 in diffusion-limited reactions, i.e., when the binding probability of S and E is approximately equal to 1 and diffusion plays a significant role.^22,23

3. Methods

3.1. Lattice gas automata algorithm

Each simulation begins with E, S, and obstacle molecules (if present) randomly placed on a 2D lattice of size 100 × 100 with cyclic boundary conditions. At each time step, the number of potential moves for a particle is given by the coordinate number of the type of lattice used. In our model particles can move into six directions (coordinate number = 6). We select a coordinate number equal to 6 rather than a coordinate number equal to 4 because the former provides a better approximation to physically realistic reaction conditions.²⁴

Unless otherwise stated, the initial fraction density of E, S, C and P in the lattice are [E₀] = 0.01, [S₀] = 0.1, [C₀] = 0 and [P₀] = 0. θ varies between 0 and 0.40. Obstacles are mobile, square and of size 1 (1 × 1), 16 (4 × 4), 49 (7 × 7) or 100 (10 × 10). If the obstacle size is (l1 × l2), then its width is l1 and height is l2, which represent the number of occupied lattice sites in each direction. The rate coefficients k₁, k₂ and k₃ are modelled by the reaction probabilities f, r and g, respectively (Table 1). Here, f is the probability that E and S will react to form C (E + S → C) given that they have collided, and r and g are the probabilities that a given C will disassociate into E and S (C → E + S) or E and P (C → E + P), respectively. With only one exception (Section 4.3), the values of f, r and g are 1, 0.02 and 0.04, respectively. These values are identical to those used in previous simulations,^3,25 which makes these results suitable for comparison with published literature in the subject.

Table 1 Notation used in the model and simulations

Symbol	Description
t	Time
τ	Lifetime of a particle
E	Enzyme
S	Substrate
C	Complex
P	Product
O	Obstacle
θ	Density of obstacles
δ	Size of obstacles
α	Anomalous diffusion exponent
D*	Limiting diffusion coefficient
t CR	Crossover time from anomalous to regular diffusion
h	Fractal-like kinetic parameter
Spin	Probability of spinning per time step
η	Fraction of obstacles that interact with enzymes to form EO
k 1	Simulated rate of E + S → C
k 2	Simulated rate of C → E + S
k 3	Simulated rate of C → E + P
w 1	Probability of E + O → EO
w 2	Probability for EO → E + O
f	Probability of E + S → C
r	Probability of C → E + S
g	Probability of C → E + P

---

The movement of a particle follows Brownian motion, so at each time-step of the Monte Carlo sequence, a “subject” molecule moves and interacts according to the following rules (Fig. 1):

Fig. 1 Decision process for the movement and rotation of a particle in the simulations.

Randomly choose a “destination” site from six possible positions neighbouring the “subject” site.

(a) If the “destination” site is not vacant,

(i) If the molecule occupying the “subject” site is E or S and the molecule occupying the “destination” site is S or E respectively, generate a random number between 0 and 1. If this random number is lower than reaction probability f, replace the “destination” site molecule with C and remove the molecule occupying the “subject” site, making that a vacant site. Set y = y + 1, where y counts the number of E + S → C reactions that have occurred in the interval [0,t].

(ii) Otherwise, no movement or reaction occurs.

(b) If the “destination” site is vacant and the “subject” molecule is C,

(i) Generate a random number between 0 and 1.

(ii) If the random number is less than the reverse reaction probability r, place E on the “subject” site and S on the “destination” site. This corresponds to a reverse reaction C → E + S.

(iii) If the random number is greater than r but less than r + g, where g is the probability of a catalytic reaction, place E on the “subject” site and P on the “destination” site. This corresponds to the catalytic reaction C → E + P.

(iv) Otherwise, move C to the “destination” site.

(c) If the “destination” site is vacant and the “subject” molecule is not C,

(i) Move the molecule to the “destination” site.

At each time step, the total number of possible moves is equal to the total number of mobile molecules present in the lattice, so that on average, every molecule is moved once. Notice that in the MM mechanism, as particles bind together and break apart, the total number of mobile particles changes throughout the simulation. One simulation ends when either the density of substrates is lower than 0.01, or the number of time steps reaches 1000. For each time step of the simulation, we obtain the average mean square displacement, lifetime and fraction of occupied sites (densities) of every reacting molecule and obstacle. We run 500 replicates of each set of parameter values to obtain averages at each time step.

3.2. Estimation of α, D* and t_CR

To obtain the average diffusion of a particular type of particle, we interpret t* in eqn (3) as the mean lifetime (τ) of the particle. We obtain α, D* and t_CR by analysing the logarithm of diffusion on a logarithmic timescale. We calculate α from the slope (α − 1) of diffusion in the initial five time steps of the simulation, as seen in the log–log plot. D* is obtained by collecting the average diffusion of the last 250 simulation time steps, when diffusion is already at steady state. We obtain t_CR, the time required to crossover from anomalous to regular diffusion, by intersecting the slope line of anomalous diffusion with the constant line obtained from regular diffusion at long times.

3.3. Estimation of k₁, k₂ and k₃

As explained in Section 3.1, we use y as a counter of the forward reactions (E + S → C) that occur in the interval [0,t]. We assume that this counter is related to k₁ by the following integral expression:We deduce k₁ with a simple derivative of eqn (6). Using this derivative and making a simple algebraic rearrangement of the law of mass action equations for the substrate and product, we also obtain k₂ and k₃:^3,25

For the estimation of these parameters, we assume that fluctuations are small enough such that the mean densities are accurately given by the rate equations.³⁰ Parameter values are chosen arbitrarily in the simulations to facilitate comparison between our results and those published in the literature.

3.4. Particle initialization and orientation

Every particle, except obstacles of size greater than one, is randomly initialized with a given orientation and direction of rotation. There are six possible orientations, corresponding to the coordinate number equal to six. The direction of rotation is always clockwise (CW) or counter-clockwise (CCW). Every time a particle moves, the particle is reoriented one unit in the predetermined direction, so as to maintain an angular momentum. This reorientation occurs according to a spin probability (spin). Every time a particle ‘collides’ with another particle and no reaction occurs, both particles are reoriented in a new randomly chosen direction. E only binds to S when both particles' orientations are complementary for an interaction. The particle's movement is also always independent of its orientation, except in the reversible or catalytic step of the MM mechanism. Here, the movement of S (reversible step) or P (catalytic step) occurs in the direction that is complementary to its position of molecules in C.

3.5. Weak force interactions

In order to mimic non-specific interactions, we allow the binding of E to obstacles of similar size. E reversibly binds to an obstacle (O) the same way it binds to S, except there is no catalytic step:In eqn (8), w₁ and w₂ are the rates at which E associates to and dissociates from O, respectively. In this particular implementation, O acts as a trap to E and the time E remains trapped increases with an increasing ratio of w₁ to w₂. This approach is somehow similar to the approach of Saxton,¹³ in which tracer particles (E) have to search through a large number of nonreactive binding sites to find its target site (S). At the beginning of each simulation, we randomly assign which O may interact with E. This fraction of O (η) is predetermined and given as an input to the simulation. The rules of interaction follow the algorithm described in Section 3.1.

4. Results

4.1. As obstacle density increases, obstacle size plays a larger role in fractal-like kinetics

In order to test how diffusion and h are affected by either the size or the density of obstacles, we run simulations varying θ = [0, 0.1, 0.2, 0.3, 0.4] and the obstacle size, δ = [1, 16, 49, 100]. Using an average of 500 replicate simulations for each parameter combination, we obtain α, D* and t_CR as characteristics of the diffusion process (see Section 3). The log–log plot of diffusion over time shows anomalous diffusion. Diffusion is irregular (non-constant) during an initial transient, but then becomes regular after a longer timescale (see Fig. 2A). This behaviour is consistent with the findings of Platani et al.¹⁸ and Saxton.¹¹D* for a density of 1% of a single particle holds D* = 1 (results not shown). In the absence of O our 2D simulations of the MM mechanism exhibit constant diffusion at long times. The diffusion coefficient is close to, but not exactly 1 due to auto-crowding effects. Auto-crowding effects become visible when a fraction of sites occupied by E and S is large enough to affect reaction kinetics, even in the absence of O.

Fig. 2 Log–log plots of simulated substrate diffusion and forward rate k₁ as a function of time in the Michaelis–Menten mechanism. (A) Diffusion is anomalous for short times and regular for long times. Diffusion in a non-crowded environment is close to, but not exactly equal to 1 because of auto-crowding effects. (B) Long time behaviour shows fractal-like kinetics characteristics. We use the slope of the curve at long times to measure the fractal-like kinetic parameter h. Plotted values represent averages of 500 replicates each running for a maximum of 1000 time steps, with [E₀] = 0.01, [S₀] = 0.1, f = 1, r = 0.02 and g = 0.04.

The log–log plot of k₁, i.e., the rate of the reaction E + S → C, yields a linear decay at long times (see Fig. 2B). This outcome is consistent with fractal-like kinetics and supported by previous work.^3,25 Note that in the absence of O, fractal-like kinetics is still observed in 2D environments, contrary to what is observed in a 3D system.^3,25 A systematic analysis of the diffusion parameters of S and h as a function of θ and δ shows that α and t_CR do not significantly change with θ (Fig. 3A and C). These results contradict previous work that suggests a decrease of the anomalous diffusion exponent (higher slope) and an increase of h with increasing θ.^20–23 While it seems clear that anomalous diffusion is related or contributes to anomalous kinetics, it is not so clear in what conditions the anomalous diffusion exponent decreases with increasing θ. Our observations may in part be due to the low dimensionality of our system. In 3D simulations, θ may have lower impact on D* at the initial time steps and the system may respond differently to different θ. This 3D behaviour is seen in the work of Vilaseca,¹⁶ where α is shown to decrease with increasing θ. However, the authors only model single particle type diffusion and do not explicitly model reactions (see Discussion). Notice that α is on average lower than one in all of our results. This implies that anomalous diffusion is always present at short times.

Fig. 3 Anomalous diffusion exponent (α) (A), diffusion values at long times (D*) (B), crossover time (t_CR) (C), and fractal-like kinetic parameter (h) (D) as a function of obstacle density (θ) and obstacle size (δ). (A) The anomalous diffusion exponent shows no significant differences as a function of either obstacle density or size. (B) Diffusion is inversely proportional to obstacle density, but proportional to obstacle size. (C) Like α, the crossover time from anomalous to regular diffusion shows no significant differences as a function of either obstacle density or size. (D) As expected, h increases with increasing obstacle density. Plotted values represent averages of 500 replicates, each running for a maximum of 1000 time steps, with [E₀] = 0.01, [S₀] = 0.1, f = 1, r = 0.02 and g = 0.04.

Our results show that the limiting D* of S is dependent on θ and δ (Fig. 3B). As θ increases, D* decreases. In fact, D* can be reduced from ∼0.9 in an obstacle free environment to ∼0.5 in a highly crowded environment. The diffusion of S also increases as δ increase. These results are consistent with the recent work of Vilaseca.¹⁶ For a fixed occupied fraction of sites, increasing δ creates more space for free diffusion.

As expected, h increases with increasing θ (Fig. 3D). As θ increases, δ plays a larger role in fractal-like kinetics. While differences between obstacles of size greater than one are not very clear, there are obvious differences between the value of h obtained with obstacles of size one and greater. When no obstacles are present, h is ∼0.2. However, when θ is 0.4, h is ∼0.3 for obstacles of size 1, but ∼0.45 for obstacles of size 100. Together, these results indicate that increasing diffusion of S and lowering diffusion of O exacerbate fractal-like kinetics behaviour. As previously shown by Saxton⁷ in tracer-obstacle studies, the diffusion coefficient of O is similar to the diffusion of tracers at long times when the two particle types have similar size. However, as θ increases, the ratio of diffusion of O to the diffusion of S diminishes approximately linearly for obstacles of size greater than one (Fig. 4).

Fig. 4 Ratio of the diffusion of obstacles to the diffusion of substrates as a function of obstacle density (θ) and obstacle size (δ) in the Michaelis–Menten mechanism. When δ = 1, all particles have the same size, and hence, similar diffusion. When δ > 1, diffusion of obstacles is inversely proportional to obstacle density and diffusion of substrates proportional to obstacle density (see also Fig. 2B). Plotted values represent averages of 500 replicates, each running for a maximum of 1000 time steps, with [E₀] = 0.01, [S₀] = 0.1, f = 1, r = 0.02 and g = 0.04.

4.2. Both first- and second-order rate coefficients exhibit dependence on the density of obstacles

We evaluate the dependence of the first-order (k₂ and k₃) and the second-order (k₁) rate coefficients of the MM mechanism on θ and δ. Because we are interested in k₁ at long times, we obtain the average and standard deviations of the last 100 time steps of k₁. As expected, k₁ decreases with increasing θ for particles of size 1 and 16 (Fig. 5A). When obstacles of size 1 (δ = 1) are present in the system, k₁ is ∼0.9. At the maximum density of obstacles (θ = 0.4), k₁ is reduced to an average of ∼0.6. When larger obstacles (δ > 16) are present, k₁ seems to increase at lower densities and then decrease at higher densities. The slope increases when the higher diffusion of S and E (induced by increasing δ) trumps the lower diffusion coefficient (induced by increasing θ). This explains the decreasing slope at a higher θ for particles of size 100 (between θ = 0.3 and θ = 0.4). We also evaluate the first-order rate coefficients as a function of θ and δ. As θ increases, both reversible and catalytic rates decrease. This decrease is linear for obstacles of size 1, but appears with slight curvature for obstacles of greater size (Fig. 5B and C). For obstacles of size 1, k₂ and k₃ change from ∼0.017 and ∼0.035 when no obstacles are present in the system, to ∼0.01 and ∼0.02, respectively, when a maximum density of obstacles are present. The slope diminishes for obstacles of larger size, indicating that the values of k₂ and k₃ become less dependent on θ as δ increases. This is the expected outcome when the diffusion of particles, other than O, increases with the increasing δ.

Fig. 5 k ₁, k₂ and k₃ mean and standard deviations as a function of obstacle density (θ) and obstacle size (δ) in the Michaelis–Menten mechanism. (A) k₁ is inversely proportional to obstacle density for δ = 1 and δ = 16. However, when δ = 49 and δ = 100, k₁ increases with increasing obstacle density. The latter behaviour shifts at high densities. (B and C) k₂ and k₃ are inversely proportional to obstacle density, but proportional to the obstacle size. We obtain the mean and standard deviations of k₁ by averaging the last 100 iterations of each simulation. We obtain the mean and standard deviations of k₂ and k₃ by averaging the entire time course of each simulation. Plotted values represent averages of 500 replicates, each running for a maximum of 1000 time steps, with [E₀] = 0.01, [S₀] = 0.1, f = 1, r = 0.02 and g = 0.04.

Previous simulations of MM mechanism do not reveal dependence of the first-order rate coefficients on θ.^3,25 In our simulations (Section 3), we find that k₂ and k₃ are constant (consistent with the law of mass action), but with values that are only a fraction of the first-order rates obtained when no obstacles are present in the system (see Section 5). Since diffusion of a particle only occurs when a selected destination is vacant, we hypothesise that the fraction by which k₂ and k₃ diminish in a simulation is exactly set by the dimensionless diffusion value of the complex C at long times. To test this hypothesis, we compare the existent theoretical approaches with a modified version of the fractal-like kinetics over the time course of E, S, C and P. The corresponding governing ordinary differential equations are:

In the above system, k₁ is represented by the traditional fractal-like rate coefficient k₀t^−h, where h is the fractal-like kinetic parameter. k₂ and k₃ are multiplied by the dimensionless complex diffusion value at long times. The above governing rate – a new form of anomalous reaction kinetics – describes the governing behaviour of the time course of the four MM species with higher accuracy than classical kinetics (CK) and the fractal-like kinetics (FK) (Fig. 6). The classical kinetics (law of mass action) expression is obtained making h = 0 and having no diffusion explicitly represented in the equations. In the traditional fractal-like kinetics, the parameter h is obtained by computing the slope of a linear fit to the last 100 iterations of the log (k₁) on a logarithmic timescale note that log (k₁) decays linearly at long times in a logarithmic timescale, as shown in eqn (5). As with the classical kinetics expression, no diffusion is explicitly considered. We tested the same [E₀], [S₀] and reaction probabilities applied by Berry²⁵ and obtained the same outcome. In addition, we tested several other different conditions with no observed changes (results not shown).

Fig. 6 Lattice gas automata (LGA), classical kinetics (CK), fractal-like kinetics (FK) and the new form of anomalous kinetics (AK) over time for the enzyme (A), substrate (B), complex (C) and product (P) species in the Michaelis–Menten mechanism. This new form of anomalous kinetics explains lattice gas automata simulation results with higher accuracy. We use [E₀] = 0.01, [S₀] = 0.1, [O] = 0.4, δ = 16, f = 1, r = 0.02, g = 0.04.

4.3. Particle orientation reduces forward rate coefficient and fractal-like kinetics behaviour

Previous studies consider two different probabilities for the formation of C: (a) the probability that both E and S compete for the same lattice spot at the same time; and (b) the probability that E and S react given that they compete for the same lattice spot at the same time.³ While the former is affected by obstacles in the media that may increase the average time taken by any two particles to interact, the latter is solely affected by the affinity of the two molecules. In reality, although (a) and (b) are necessary, they are not sufficient. E and S can only react and form C if the orientation of the active site of E is complementary to the orientation of the binding site of S. The introduction of an orientation to every particle is expected to reduce the rate at which complexes are formed and as a consequence, reduce h. To test this hypothesis, we modify our algorithm to introduce an orientation to every particle in the system (see Section 3 and Fig. 1). At each time step, every particle is rotated with a ‘spin’ probability in a predetermined direction (CW/CCW). The higher the ‘spin’ probability, the faster the particles will rotate. Even at the maximum obstacle density (θ = 0.4), spin reduces k₁ and h (see Fig. 7). At the initial transient, log (k₁) reduces from ∼0.4 to ∼−0.5 by adding a ‘spin’ of 100%. We systematically analyse h as a function of ‘spin’, θ and δ. With θ = 0.1, h is clearly reduced as ‘spin’ increases. For instance, when we consider δ = 1, h is reduced from ∼0.2 to ∼0.05 (Fig. 8A). As θ increases, the effects of spinning on h diminish, except for the significant decrease one obtains by just adding 20% ‘spin’ (Fig. 8B–D). From ‘spin’ = 0 to 0.20, h reduction is ∼0.1 in Fig. 8B and C and 0.15 in Fig. 8D. No major differences appear significant between obstacles of different size. Also, and as expected, notice that as θ increases, h increases. On average, from no ‘spin’ to 100% ‘spin’, we see a reduction of ∼66% for θ = 0.1, ∼40% for θ = 0.2, and ∼50% for θ = 0.3 and θ = 0.4 (Fig. 8).

Fig. 7 Log–log plot of rate k₁ over time comparing the use of particle orientation versus the traditional method. When particle orientation is taken into account, the probability of a forward reaction event is diminished and the fractal-like kinetic parameter h decreases. Plotted values represent averages of 500 replicates, each running for a maximum of 1000 time steps, with [E₀] = 0.01, [S₀] = 0.1, [O] = 0.4, f = 1, r = 0.02 and g = 0.04.

Fig. 8 Fractal-like kinetic parameter h as a function of the probability of spinning for the obstacle size (δ) and obstacle density θ = 0.1 (A), θ = 0.2 (B), θ = 0.3 (C) and θ = 0.4 (D). It is not clear that particle size plays a significant role, but h decreases with increasing spinning. Plotted values represent averages of 500 replicates, each running for a maximum of 1000 time steps, with [E₀] = 0.01, [S₀] = 0.1, f = 1, r = 0.02 and g = 0.04.

E and S that have less binding affinity should also reduce fractal-like kinetics, as the rate of C formation is decreased. We evaluate h as a function of f, θ and δ. With or without obstacles, h increases with f (Fig. 9). When no obstacles are present (θ = 0), h is reduced from ∼0.18 (f = 1) to ∼0.03 (f = 0.02), a reduction of more than 80% (Fig. 9A). At higher obstacle densities, the reduction between f = 1 and f = 0.2 is less pronounced. The reduction is ∼0.57% and ∼0.55% with θ = 0.2 and θ = 0.4, respectively. Obstacle size does not appear to be significant, although as expected, obstacles of lower size exhibit less fractal-like kinetics behaviour.

Fig. 9 Fractal-like kinetic parameter h as a function of reaction probability f and the obstacle size (δ) for obstacle densities θ = 0 (A), θ = 0.2 (B) and θ = 0.4 (C). The fractal-like kinetic parameter h increases with increasing f. h seems distinctly lower for obstacles of size 1. No significant differences are seen for obstacles of larger size. Plotted values represent averages of 500 replicates, each running for a maximum of 1000 time steps, with [E₀] = 0.01, [S₀] = 0.1, r = 0.02 and g = 0.04. Note that f is always at least five times higher than the catalytic rate g.

4.4. Weak force interactions reduce fractal-like kinetics behaviour

We analyse the effects of weak force interactions on the fractal-like kinetics. In particular, we allow a fraction of obstacles (η) to bind and trap enzymes (see Section 3 and Fig. 1). The time an enzyme remains trapped increases with an increasing ratio of w₁ to w₂. For simplification, we restrict the interaction of E and O to cases where both particles have the same size and occupy one spot in the lattice. Also we only analyse different η fractions with a fixed obstacle density (θ = 0.4). When w₁ = 10w₂, h is ∼0. This can be seen in the long time behaviour of the log–log plot of k₁ over time (Fig. 10A).

Fig. 10 Log–log plot of rate k₁ over time (A), and plot of fractal-like kinetic parameter h as a function of the fraction of obstacles (η) with which enzymes may interact (B). (A) Comparison of k₁ when η = 0, and k₁ when η = 1 with w₁ = 10w₂, where w₁ and w₂ are probabilities of the reactions E + O → EO and EO → E + O, respectively. h is ∼0 for the latter case. (B) Fractal-like kinetic parameter h as a function of η. When w₂ = 10 w₁, no differences are seen as η is increased. When w₁ ≥ w₂, h decreases with increasing η. The rate at which h decreases is proportional to the probability of an enzyme being trapped by an obstacle. Plotted values represent averages of 500 replicates, each running for a maximum of 1000 time steps, with [E₀] = 0.01, [S₀] = 0.1, [O] = 0.4, f = 1, r = 0.02 and g = 0.04.

We systematically analyse h as a function of η and three different ratios of w₁ to w₂. When w₂ = 10w₁, h is approximately constant (∼0.28) as η is increased (Fig. 10B). No significant differences appear, since there is a very low trapping probability. When w₁= w₂, h appears to decrease as η increases to 0.4 and 0.6, so to maintain or slightly increase afterwards. When w₁ = 10w₂, h appears to decrease as η increases. It is important to notice that as obstacles trap enzymes, enzymes exist at lower densities in comparison to what is expected without the weak force interactions. The exacerbated effect on h as the ratio of w₁ to w₂ increases is due to a larger number of enzymes trapped by obstacles. As a result, fewer enzymes are available to bind and react with S. P formation is therefore reduced, and slowly increases over time.

5. Discussion

Reactions inside cells can take place in a heterogeneous crowded environment, far from the homogeneous, well mixed solution that characterizes the typical test tube experiments. Macromolecular crowding affects reaction kinetics because it influences the displacement of molecules, and hence, the probability of interactions between two reacting molecules. There is ample evidence on the impact of obstacles on diffusion and reaction kinetics.^3,7,11 Importantly, Kopelman²³ shows that in diffusion-limited elementary reactions, the law of mass action breaks down in spatially heterogeneous environments with obstacles. His new description of reaction dynamics, known as fractal-like kinetics, has been shown to be accurate predominantly using lattice based models in the diffusion-limited regime.^3,25 However, scepticism remains on these results as lattice approaches spatially restrict the degree of particle movement.²⁴ In 2010, Grima²⁸ used a modified chemical-master equation to investigate the role of obstacle crowding in reaction kinetics. His theory is consistent with a simulated dimerization reaction, which suggest that fractal-like kinetics could be an artefact of lattice-based simulations. This result was also confirmed analytically using the generalised mass action kinetics theory.³¹ However, these results are only valid in the reaction-limited regime, where the probability of a reaction is low. Here, we took a new approach and developed an enhanced 2D lattice simulation of the MM mechanism in the diffusion-limit regime, where fractal-like kinetics is predicted to take place. Our enhanced version of the lattice simulation provides a consistent diffusion coefficient for all particles of similar size and explicitly models particle rotation and weak force interactions among particles in the lattice. In this new context, we investigate the role of obstacle density and size on diffusion and rate coefficients.

We made obstacles mobile, so as to have obstacle diffusion implicitly set as a function of its size. We also measured the diffusion coefficients of every particle with respect to their lifetime, so as to understand the relations between diffusion of obstacles and any of the four constituents of the MM mechanism (Fig. 2 and 4). In this work we choose obstacle densities within the intracellular macromolecular crowding density range reported in the literature and four obstacle sizes that provide significantly distinct diffusion coefficients.⁴ In all of our simulations, E, S, C and P diffusion coefficients exhibit a crossover from anomalous to regular diffusion and approach similar values at long times (Fig. 3). This is also observed in the linear behaviour of log–log plots of the mean square displacement over lifetime, although C naturally exhibits short lifetimes in comparison to E, S and P (results not shown). Our results are consistent with previous findings.⁷ Our simulation results also show that the anomalous diffusion exponent or crossover time is independent of the density and size of obstacles (Fig. 3). If such dependence exists as previous work suggests, such behaviour can become visible in 3D simulations, where the effect of obstacle densities would have lower impact on initial time step measurements of the diffusion coefficient.^7,9,16 Our results are in good agreement with the full energy simulation model of McGuffee and Elcock,²⁷ where anomalous diffusion behaviour of several interacting molecules appears independent of their molecular weight. Our simulations show that increasing obstacle size contributes to increased diffusion, which in turn contributes to higher fractal-like kinetics when that diffusion is hampered by a higher obstacle density (Fig. 3).

Previous attempts to model biochemical reactions using lattice gas automata did not explicitly consider particle diffusion.^3,25 They also considered particles moving on four directions, rather than on six directions. However, the major difference between those models and our model resides on the behaviour of the complex particle or C. In previous models, while E, S and P are allowed to react or move in a time step using one randomly selected neighbouring position, C is allowed to move or react using any of the vacant neighbouring positions. This endows C with a special behaviour, and it is not realistic as it makes C move faster than any other MM species (results not shown). Provided that C has the same size as E, S or P and therefore has similar diffusion coefficients at long times, we apply the same rule to every particle and make C react or move only if the selected destination (one of the six possible destinations) is unoccupied (see Section 3 and Fig. 1). This behaviour explains the first-order rates dependency on obstacle density, which has not been previously reported (see Fig. 5). Importantly, it may explain the dependency of V_max and K_M on obstacle concentration observed in in vitro enzyme-catalysed reaction experiments.^32,33 We compare different theoretical models of reaction kinetics over the time course of the MM species and show that the difference between fractal-like kinetics occurring with and without obstacles is exactly given by the diffusion coefficient of C at long times (Fig. 6).

We showed here that particle rotation can reduce the fractal-like kinetic parameter by ∼50% (Fig. 7 and 8). This reduction can be further augmented in 3D simulations of biochemical reactions.³⁴ We also systematically reduced the forward probability f to show the relationship between the affinity of two particles and fractal-like kinetics. As the particle's affinity decreases, so does fractal-like kinetics (Fig. 9). In sum, even if reactants have high affinity, they may only react if their orientations are complementary, which reduces the rate coefficient. In 2008, Saxton¹³ showed that the time taken for a mobile reactant to find its target site increases with the number and affinity of nonreactive binding sites in the environment. Here, we show that weak force interactions, specifically those that cause enzyme sequestration, contribute to diminishing the density of free enzymes, resulting in slower enzymatic reactions and lower fractal-like kinetics behaviour (Fig. 10). These results suggest that fractal-like kinetics behaviour is also dependent on the ratio of enzyme to substrate densities, reducing its impact as this ratio is also reduced.

6. Conclusions

Previous research has focused on the effects of density and size of obstacles in diffusion. Here, we expand this investigation to diffusion-limited enzyme-catalysed reactions. We show that both first-order and second-order rate coefficients as well as the fractal-kinetic parameter are affected by the density and size of obstacles. We also propose a new form of anomalous kinetics that includes diffusion to explain lattice gas automata simulation results.

In addition, we investigate the effects of particle rotation and weak force interactions in the reaction kinetics. We find that both particle rotation and weak force interactions shift the balance from diffusion-limited to reaction-limited kinetics. This significantly reduces fractal-like kinetics, suggesting that in lattice based models, fractal-like kinetics generally appear under more restrictive conditions than previously believed.

Acknowledgements

The authors are grateful for Dr Ramon Grima's (University of Edinburgh) and Mark Whidden's (University of Michigan) critical comments of this paper. This work has been funded by a grant of the James S. McDonnell Foundation under the 21st Century Science Initiative, Studying Complex Systems Program.

References

1. P. Waage and C. M. Gulberg, J. Chem. Educ., 1986, 63, 1044–1047.
2. B. G. Cox, Modern liquid phase kinetics, Oxford University Press, Oxford, New York, 1994.
3. S. Schnell and T. E. Turner, Prog. Biophys. Mol. Biol., 2004, 85, 235–260.
4. D. Hall and A. P. Minton, Biochim. Biophys. Acta, 2003, 1649, 127–139.
5. S. B. Zimmerman and A. P. Minton, Annu. Rev. Biophys. Biomol. Struct., 1993, 22, 27–65.
6. A. P. Minton, J. Biol. Chem., 2001, 276, 10577–10580.
7. M. J. Saxton, Biophys. J., 1994, 66, 394–401.
8. M. J. Saxton, Biophys. J., 1997, 72, 1744–1753.
9. M. Weiss, M. Elsner, F. Kartberg and T. Nilsson, Biophys. J., 2004, 87, 3518–3524.
10. D. S. Banks and C. Fradin, Biophys. J., 2005, 89, 2960–2971.
11. M. J. Saxton, Biophys. J., 2007, 92, 1178–1191.
12. J. A. Dix and A. S. Verkman, Annu. Rev. Biophys., 2008, 37, 247–263.
13. M. J. Saxton, Biophys. J., 2008, 94, 760–771.
14. I. Pastor, E. Vilaseca, S. Madurga, J. L. Garces, M. Cascante and F. Mas, J. Phys. Chem. B, 2010, 114, 4028–4034.
15. E. Vilaseca, I. Pastor, A. Isvoran, S. Madurga, J.-L. Garcés and F. Mas, Theor. Chem. Acc., 2011, 128, 795–805.
16. E. Vilaseca, A. Isvoran, S. Madurga, I. Pastor, J. L. Garces and F. Mas, Phys. Chem. Chem. Phys., 2011, 13, 7396–7407.
17. M. J. Saxton, Biophys. J., 2012, 103, 2411–2422.
18. M. Platani, I. Goldberg, J. R. Swedlow and A. I. Lamond, J. Cell Biol., 2000, 151, 1561–1574.
19. A. Kusumi, C. Nakada, K. Ritchie, K. Murase, K. Suzuki, H. Murakoshi, R. S. Kasai, J. Kondo and T. Fujiwara, Annu. Rev. Biophys. Biomol. Struct., 2005, 34, 351–378.
20. A. Lin, R. Kopelman and P. Argyrakis, Phys. Rev. E, 1996, 53, 1502–1509.
21. A. L. Lin, M. S. Feldman and R. Kopelman, J. Phys. Chem. B, 1997, 101, 7881–7884.
22. R. Kopelman, J. Stat. Phys., 1986, 42, 185–200.
23. R. Kopelman, Science, 1988, 241, 1620–1626.
24. R. Grima and S. Schnell, Biophys. Chem., 2006, 124, 1–10.
25. H. Berry, Biophys. J., 2002, 83, 1891–1901.
26. R. Grima, S. N. Yaliraki and M. Barahona, J. Phys. Chem. B, 2010, 114, 5380–5385.
27. S. R. McGuffee and A. H. Elcock, PLoS Comput. Biol., 2010, 6, e1000694.
28. R. Grima, J. Chem. Phys., 2010, 132, 185102.
29. D. Ben-Avraham and S. Havlin, Diffusion and reactions in fractals and disordered systems, Cambridge University Press, Cambridge, New York, 2000.
30. R. Grima, J. Chem. Phys., 2010, 133, 035101.
31. R. Grima and S. Schnell, ChemPhysChem, 2006, 7, 1422–1424.
32. M. Jiang and Z. Guo, J. Am. Chem. Soc., 2007, 129, 730–731.
33. I. Pastor, E. Vilaseca, S. Madurga, J. L. Garces, M. Cascante and F. Mas, J. Phys. Chem. B, 2011, 115, 1115–1121.
34. M. J. Saxton, Methods Mol. Biol., 2007, 400, 295–321.

---