Unravelling the impact of obstacles in diffusion and kinetics of an enzyme catalysed reaction
Received
10th June 2013
, Accepted 10th October 2013
First published on 11th October 2013
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.1 Although extensively applied to model biochemical reactions inside cells, the law of mass action assumes that reactions occur in homogeneous, well stirred volumes.2 This is inconsistent with the heterogeneous environment that characterizes intracellular reactions.3 Intracellular environments possess a high macromolecular content, also known as macromolecular crowding, in the range of 5% to 40% of the total cellular volume.4 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.11 These three phases have been shown to exist in the diffusion of Cajal bodies in the nucleus of HeLa cells18 and in the diffusion of gold-labeled dioleoyIPE in the plasma membrane of fetal rat skin keratinocyte cells.19 Tracer-obstacle simulations have also shown that the anomalous diffusion slope and occurring period increase with increasing obstacle density.7 More recently, Vilaseca16 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, Kopelman22,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, Grima26 showed that deviations from diffusional isotropy significantly reduce the diffusion coefficient of tracer particles in crowded environments. Also in 2010, McGuffee and Elcock27 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.28 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: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:29where Γ 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 D0 as the diffusion coefficient of a particle in an unobstructed system (θ = 0), we simplify eqn (1) by making r = lr*, t = ωt*, D = D0D*(θ) and l2 = 2dD0ω. We then obtain the following dimensionless equation:7
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 (tCR). The slope α − 1 of the initial time steps is used to calculate α. D* is the value of diffusion as time approaches infinity. tCR 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.7
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): | ![image file: c3cp52417e-t1.tif](/image/article/2014/CP/c3cp52417e/c3cp52417e-t1.gif) | (4) |
where k1, k2 and k3 are rate constants. When few obstacles are present in the reaction media, k1 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 k1 decays over time.22,23 As a consequence, log
(k1) decays linearly at long times in a logarithmic time scale:where k0 is the rate at t = 0 and h is the fractal-like kinetic parameter. Notice that when h is 0, k1 is no longer time-dependent. Therefore, h can be used to measure deviations from the law of mass action.24 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.24
Unless otherwise stated, the initial fraction density of E, S, C and P in the lattice are [E0] = 0.01, [S0] = 0.1, [C0] = 0 and [P0] = 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 k1, k2 and k3 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):
![image file: c3cp52417e-f1.tif](/image/article/2014/CP/c3cp52417e/c3cp52417e-f1.gif) |
| 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 tCR
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 tCR 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 tCR, 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 k1, k2 and k3
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 k1 by the following integral expression: | ![image file: c3cp52417e-t2.tif](/image/article/2014/CP/c3cp52417e/c3cp52417e-t2.gif) | (6) |
We deduce k1 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 k2 and k3:3,25 | ![image file: c3cp52417e-t3.tif](/image/article/2014/CP/c3cp52417e/c3cp52417e-t3.gif) | (7) |
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.30 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: | ![image file: c3cp52417e-t4.tif](/image/article/2014/CP/c3cp52417e/c3cp52417e-t4.gif) | (8) |
In eqn (8), w1 and w2 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 w1 to w2. This approach is somehow similar to the approach of Saxton,13 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 tCR 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.18 and Saxton.11D* 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.
![image file: c3cp52417e-f2.tif](/image/article/2014/CP/c3cp52417e/c3cp52417e-f2.gif) |
| Fig. 2 Log–log plots of simulated substrate diffusion and forward rate k1 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 [E0] = 0.01, [S0] = 0.1, f = 1, r = 0.02 and g = 0.04. | |
The log–log plot of k1, 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 tCR 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,16 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.
![image file: c3cp52417e-f3.tif](/image/article/2014/CP/c3cp52417e/c3cp52417e-f3.gif) |
| Fig. 3 Anomalous diffusion exponent (α) (A), diffusion values at long times (D*) (B), crossover time (tCR) (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 [E0] = 0.01, [S0] = 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.16 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 Saxton7 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).
![image file: c3cp52417e-f4.tif](/image/article/2014/CP/c3cp52417e/c3cp52417e-f4.gif) |
| 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 [E0] = 0.01, [S0] = 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 (k2 and k3) and the second-order (k1) rate coefficients of the MM mechanism on θ and δ. Because we are interested in k1 at long times, we obtain the average and standard deviations of the last 100 time steps of k1. As expected, k1 decreases with increasing θ for particles of size 1 and 16 (Fig. 5A). When obstacles of size 1 (δ = 1) are present in the system, k1 is ∼0.9. At the maximum density of obstacles (θ = 0.4), k1 is reduced to an average of ∼0.6. When larger obstacles (δ > 16) are present, k1 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, k2 and k3 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 k2 and k3 become less dependent on θ as δ increases. This is the expected outcome when the diffusion of particles, other than O, increases with the increasing δ.
![image file: c3cp52417e-f5.tif](/image/article/2014/CP/c3cp52417e/c3cp52417e-f5.gif) |
| Fig. 5
k
1, k2 and k3 mean and standard deviations as a function of obstacle density (θ) and obstacle size (δ) in the Michaelis–Menten mechanism. (A) k1 is inversely proportional to obstacle density for δ = 1 and δ = 16. However, when δ = 49 and δ = 100, k1 increases with increasing obstacle density. The latter behaviour shifts at high densities. (B and C) k2 and k3 are inversely proportional to obstacle density, but proportional to the obstacle size. We obtain the mean and standard deviations of k1 by averaging the last 100 iterations of each simulation. We obtain the mean and standard deviations of k2 and k3 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 [E0] = 0.01, [S0] = 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 k2 and k3 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 k2 and k3 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:
| ![image file: c3cp52417e-t5.tif](/image/article/2014/CP/c3cp52417e/c3cp52417e-t5.gif) | (9) |
In the above system,
k1 is represented by the traditional fractal-like rate coefficient
k0t−h, where
h is the fractal-like kinetic parameter.
k2 and
k3 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
![[thin space (1/6-em)]](https://www.rsc.org/images/entities/char_2009.gif)
(
k1) on a logarithmic timescale note that log
![[thin space (1/6-em)]](https://www.rsc.org/images/entities/char_2009.gif)
(
k1) 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
0], [S
0] and reaction probabilities applied by Berry
25 and obtained the same outcome. In addition, we tested several other different conditions with no observed changes (results not shown).
![image file: c3cp52417e-f6.tif](/image/article/2014/CP/c3cp52417e/c3cp52417e-f6.gif) |
| 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 [E0] = 0.01, [S0] = 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.3 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 k1 and h (see Fig. 7). At the initial transient, log
(k1) 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).
![image file: c3cp52417e-f7.tif](/image/article/2014/CP/c3cp52417e/c3cp52417e-f7.gif) |
| Fig. 7 Log–log plot of rate k1 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 [E0] = 0.01, [S0] = 0.1, [O] = 0.4, f = 1, r = 0.02 and g = 0.04. | |
![image file: c3cp52417e-f8.tif](/image/article/2014/CP/c3cp52417e/c3cp52417e-f8.gif) |
| 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 [E0] = 0.01, [S0] = 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.
![image file: c3cp52417e-f9.tif](/image/article/2014/CP/c3cp52417e/c3cp52417e-f9.gif) |
| 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 [E0] = 0.01, [S0] = 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 w1 to w2. 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 w1 = 10w2, h is ∼0. This can be seen in the long time behaviour of the log–log plot of k1 over time (Fig. 10A).
![image file: c3cp52417e-f10.tif](/image/article/2014/CP/c3cp52417e/c3cp52417e-f10.gif) |
| Fig. 10 Log–log plot of rate k1 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 k1 when η = 0, and k1 when η = 1 with w1 = 10w2, where w1 and w2 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 w2 = 10 w1, no differences are seen as η is increased. When w1 ≥ w2, 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 [E0] = 0.01, [S0] = 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 w1 to w2. When w2 = 10w1, h is approximately constant (∼0.28) as η is increased (Fig. 10B). No significant differences appear, since there is a very low trapping probability. When w1= w2, h appears to decrease as η increases to 0.4 and 0.6, so to maintain or slightly increase afterwards. When w1 = 10w2, 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 w1 to w2 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, Kopelman23 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.24 In 2010, Grima28 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.31 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.4 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.7 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,27 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 Vmax and KM 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.34 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, Saxton13 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
- P. Waage and C. M. Gulberg, J. Chem. Educ., 1986, 63, 1044–1047 CrossRef.
-
B. G. Cox, Modern liquid phase kinetics, Oxford University Press, Oxford, New York, 1994 Search PubMed.
- S. Schnell and T. E. Turner, Prog. Biophys. Mol. Biol., 2004, 85, 235–260 CrossRef CAS PubMed.
- D. Hall and A. P. Minton, Biochim. Biophys. Acta, 2003, 1649, 127–139 CrossRef CAS.
- S. B. Zimmerman and A. P. Minton, Annu. Rev. Biophys. Biomol. Struct., 1993, 22, 27–65 CrossRef CAS PubMed.
- A. P. Minton, J. Biol. Chem., 2001, 276, 10577–10580 CrossRef CAS PubMed.
- M. J. Saxton, Biophys. J., 1994, 66, 394–401 CrossRef CAS.
- M. J. Saxton, Biophys. J., 1997, 72, 1744–1753 CrossRef CAS.
- M. Weiss, M. Elsner, F. Kartberg and T. Nilsson, Biophys. J., 2004, 87, 3518–3524 CrossRef CAS PubMed.
- D. S. Banks and C. Fradin, Biophys. J., 2005, 89, 2960–2971 CrossRef CAS PubMed.
- M. J. Saxton, Biophys. J., 2007, 92, 1178–1191 CrossRef CAS PubMed.
- J. A. Dix and A. S. Verkman, Annu. Rev. Biophys., 2008, 37, 247–263 CrossRef CAS PubMed.
- M. J. Saxton, Biophys. J., 2008, 94, 760–771 CrossRef CAS PubMed.
- I. Pastor, E. Vilaseca, S. Madurga, J. L. Garces, M. Cascante and F. Mas, J. Phys. Chem. B, 2010, 114, 4028–4034 CrossRef CAS PubMed.
- E. Vilaseca, I. Pastor, A. Isvoran, S. Madurga, J.-L. Garcés and F. Mas, Theor. Chem. Acc., 2011, 128, 795–805 CrossRef CAS.
- E. Vilaseca, A. Isvoran, S. Madurga, I. Pastor, J. L. Garces and F. Mas, Phys. Chem. Chem. Phys., 2011, 13, 7396–7407 RSC.
- M. J. Saxton, Biophys. J., 2012, 103, 2411–2422 CrossRef CAS PubMed.
- M. Platani, I. Goldberg, J. R. Swedlow and A. I. Lamond, J. Cell Biol., 2000, 151, 1561–1574 CAS.
- 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 CrossRef CAS PubMed.
- A. Lin, R. Kopelman and P. Argyrakis, Phys. Rev. E, 1996, 53, 1502–1509 CrossRef CAS.
- A. L. Lin, M. S. Feldman and R. Kopelman, J. Phys. Chem. B, 1997, 101, 7881–7884 CrossRef CAS.
- R. Kopelman, J. Stat. Phys., 1986, 42, 185–200 CrossRef.
- R. Kopelman, Science, 1988, 241, 1620–1626 CAS.
- R. Grima and S. Schnell, Biophys. Chem., 2006, 124, 1–10 CrossRef CAS PubMed.
- H. Berry, Biophys. J., 2002, 83, 1891–1901 CrossRef CAS.
- R. Grima, S. N. Yaliraki and M. Barahona, J. Phys. Chem. B, 2010, 114, 5380–5385 CrossRef CAS PubMed.
- S. R. McGuffee and A. H. Elcock, PLoS Comput. Biol., 2010, 6, e1000694 Search PubMed.
- R. Grima, J. Chem. Phys., 2010, 132, 185102 CrossRef.
-
D. Ben-Avraham and S. Havlin, Diffusion and reactions in fractals and disordered systems, Cambridge University Press, Cambridge, New York, 2000 Search PubMed.
- R. Grima, J. Chem. Phys., 2010, 133, 035101 CrossRef CAS PubMed.
- R. Grima and S. Schnell, ChemPhysChem, 2006, 7, 1422–1424 CrossRef CAS PubMed.
- M. Jiang and Z. Guo, J. Am. Chem. Soc., 2007, 129, 730–731 CrossRef CAS PubMed.
- I. Pastor, E. Vilaseca, S. Madurga, J. L. Garces, M. Cascante and F. Mas, J. Phys. Chem. B, 2011, 115, 1115–1121 CrossRef CAS PubMed.
- M. J. Saxton, Methods Mol. Biol., 2007, 400, 295–321 CAS.
|
This journal is © the Owner Societies 2014 |