J. Phys. Chem. B 2008, 112, 4283-4289
4283
The Subdiffusive Targeting Problem Joel D. Eaves and David R. Reichman* Department of Chemistry, Columbia UniVersity, 3000 Broadway, New York, New York 10025 ReceiVed: June 22, 2007; In Final Form: October 25, 2007
Recent experiments in living cells have observed subdiffusive motion, where the mean squared displacement of a molecule grows sublinearly with time as 〈x 2〉 ∼ tR. Motivated by these experiments, we develop a theory for subdiffusion-limited bimolecular enzyme kinetics. As in normal diffusion, the statistics of the reaction times depend on the ratio between the distance to the target and the reaction rate in a scale invariant manner. Contrary to some studies, we find no critical value of R and therefore no evidence to suggest that it has been a target of evolutionary optimization.
I. Introduction Simple models of enzyme kinetics, based on phenomenological rate equations, fit bulk data from experiments in dilute solutions well. Such models have provided valuable insight into how protein structure, for example, influences reaction rates.1,2 In contrast to “well-stirred” solutions where the majority of in vitro experiments are performed, the interior of a cell appears more akin to a gel, where disorder and frustration are rules rather than exceptions. The physical environment inside a cell is a crowded ionic solution dense with proteins, organelles, and carbohydrates. Biochemists have appreciated this distinction as early as 1981, when it became apparent that DNA replication would not proceed in vitro without an inert polymer in solution to mimic the effect of cellular crowding.3 In these experiments, the crowding elements function to stabilize the native state of a protein or enzyme, but cytoplasmic crowding also has consequences for reaction kinetics. Recent observations from the dynamics of transport in living cells has put one of the fundamental assumptions underlying simple kinetic models based on the laws of mass-action into question.4 Is it reasonable to expect a crowded system to remain “well-mixed” on all time and length scales? Inside a cell, the concentration of macromolecules is large (hundreds of grams per liter) and restricts the available free volume.5 Particles navigating thorough this labyrinth can display mean squared displacements that grow sublinearly in time, as 〈x 2〉 ∼ tR, in contrast to Brownian motion where R ) 1. At first glance, this seems like a small quantitative detail. But particles moving subdiffusively display qualitatively different trajectories than diffusive ones do. A particle undergoing subdiffusive motion dwells for long periods of time in small regions of space. Long dwell times are so frequent in subdiffusion that the mean dwell time diverges.6 Conversely, divergent mean dwell times imply that the underlying trajectories are subdiffusive.23 If a signaling molecule negotiating its way through a crowded environment has to locate a distant target, it may become trapped and never encounter the target in its lifetime. Similarly, a signaling protein that begins close to the target might find it with higher probability and in less time than it would if motion was diffusive. Figure 1 A and B illustrates these two scenarios. Recent in vitro experiments performed in solutions with high random coil polymer concentrations to duplicate the effect of * Corresponding author. Electronic address:
[email protected].
macromolecular crowding in cells have measured subdiffusive motion for proteins between 2 and 20 nm in diameter.5 In another experiment, measured subdiffusive motion of colloids of ∼100 nm (about the size of cellular organelles) in diameter in a dense crosslinked network of F-actin. These experiments showed that R scaled as the diameter of the particle to the mean pore diameter of the actin network. Another study found subdiffusive motion for proteins between 10 kDa and 2 MDa in HeLa cells.8 Thus, although the mechanism and length scale for crowding in these cases can be different, the distributions of dwell times distributed according to power laws over the dynamic range of the experiment with associated subdiffusion exponents (R) ranging from ∼0.7 to 0.9.7,9 These observations cast doubt on the “fundamental hypothesis of chemical kinetics” when applied to in vivo cellular reactions.4 In a pioneering experiment, Golding and Cox9 observed subdiffusive motion and power law waiting times for the motion of an mRNA molecule in liVing bacterial cells. They made a simple argument based on scaling concepts10 to suggest that a subdiffusion exponent of R ) 2/3 implies that two molecules would always react. In gene regulation, for example, a small number of molecules (tens to hundreds) of repressors and transcription factors often modulate the state of the genetic switch. The state of the switch is sensitive to the number fluctuations of regulating molecules11 and is not always deterministic. In their measurements, R was conspicuously close to 2/3, which might imply that subdiffusion can reduce these fluctuations drastically, perhaps rationalizing the observation that regulatory sequences are often only tens of nanometers from the gene sequence of the regulating protein.12 Golding and Cox9 further showed that the subdiffusive exponent was not sensitive to severe changes in the cellular environment, including the disruption of actin-like cytoskeletal proteins. Although the subdiffusion constant, DR ∝ 〈x 2〉/tR, varied as much as an order of magnitude in different conditions, R was robust to physiological changes. A critical value of R might imply that it has been optimized under evolutionary pressure.9 On mesoscopic length and time scales (several nanometers in length and hundreds of picoseconds in time), one often describes reacting systems with spatially varying concentrations by developing partial differential equations, called reactiondiffusion equations. On these scales of measurement, the concentration of reacting molecules is large enough that it can be modeled as a continuous function in space and time. When
10.1021/jp0749017 CCC: $40.75 © 2008 American Chemical Society Published on Web 03/18/2008
4284 J. Phys. Chem. B, Vol. 112, No. 14, 2008
Eaves and Reichman
Figure 1. (A) A concentration of signaling molecules (yellow background), IPTG, for example, transfers through the cell wall of an E. coli cell (blue), navigates through a crowded environment, and modulates translation through the genetic switch (red dot) on the DNA (black squiggle). These molecules subdiffuse over distances that are hundreds of nanometers. A genetic switch (B) encodes for the regulatory sequence of repressors that control the translation of genes related to the stimulus in A. All dotted lines represent passive steps that can occur via subdiffusion, but in contrast to A, the length scales here are on the order of tens of nanometers. Either case can be modeled as a subdiffusive targeting problem (C). Each lattice site contains the number of particles that occupy a region of space of linear size . Transitions between sites occur in continuous time. The molecules in solution (white circles) wait for a time chosen from a waiting time distribution with fat fails, ψ(t), before moving to an adjacent site. The site at the origin (green circle) is the reaction zone. When molecules reach it, they can either react to form product (red circle) with an exponential waiting time dependence (red line) or re-engage in the random walk. The dashed line denotes the reactive flux of molecules that have not reacted and have moved back into solution. The waiting time distribution for this step is different than that in the subdiffusive case (solid arrows). These considerations lead to the master equations developed in this manuscript.
the transport is diffusive, the density of particles is a continuous field that evolves according to the diffusion equation. The reactions impose boundary conditions that are easy to understand and simple to derive. The same is not true of reaction-subdiffusion systems because two types of waiting time distributions appear in these problems: the one characterizing subdiffusion has fat tails and no intrinsic time scale, while the reaction time distribution has a finite mean. In contrast to reaction-diffusion models, the reaction and diffusion terms have different statistics (they are no longer both composite Markov processes). As a result, two types of boundary conditions have been analyzed in the literature: a time non-local boundary condition, which exhibits a dynamical critical point with R ) 2/3 in three dimensions,13 and a timelocal boundary condition, which has no critical value of R.14 Because phenomenology alone cannot distinguish between these two types of reactive boundaries, Section II details our derivation of the reaction-subdiffusion equation for targeting. Our approach extends Sokolov and co-workers’15 intuitive theory of sudden modulation to the bimolecular targeting problem. Once we have derived the set of master equations that describes the reaction-subdiffusion process, we develop a systematic asymptotic expansion, similar to that in ref 16, to arrive at a continuous time and space differential equation with associated reactive boundary conditions. We then show that a continuous-time Monte Carlo simulation of the discrete system agrees with the continuum approximation when two conditions are met: the particles are born several molecular diameters from one another, and the reaction rate is fast. Finally, in Section III
we discuss the implications of our results for the dynamics of chemical reactions inside the cell. II. Model Consider the binding and reaction kinetics between an enzyme and a substrate in a well-stirred solution. When the enzyme (E) and substrate (S) come together, they form an intermediate (I). This step is reversible. After binding, the intermediate can either transform the substrate into a product (P) and recycle the enzyme catalyst or unbind back to E and S. The chemical kinetic equations follow from17
k/b E + S 98 I ku I 98 E + S kf I 98 E + P
(1) (2) (3)
(E) is assumed to be massive enough that it does not diffuse on the time scale that the substrate targets. One typically assumes that the intermediate is in steady state, which implies that E is also time-independent. Equation 1 is then pseudo-first-order, with kb ) Ek/b. With these assumptions, the rate of change of the substrate concentration obeys the linear differential equation
dS ) -k S dt where k ) kbkf /(ku + kf).
(4)
Subdiffusive Targeting Problem
J. Phys. Chem. B, Vol. 112, No. 14, 2008 4285
Let us now couple this reaction scheme to the subdiffusion of S. To simplify matters, let us consider the problem in one spatial dimension and use the continuum results to generalize the results in arbitrary dimensions. By coarse graining space into cells of linear size that contain enough substrate molecules that number fluctuations can be ignored, the mass action law from eq 4 applies at the level of an individual cell (Figure 1C). Assuming that the enzyme concentration is low enough that two E molecules do not interact, let us fix E to the site at the origin, and let pn(t) represent the number of molecules of S at spatial lattice site n at time t (Figure 1C). Imagine that S executes a random walk in the cells outside the origin, n g 1. The density of S outside the reaction zone (n g1) obeys the discrete equation of continuity
dpn(t) ) J+ n (t) - Jn (t), dt
ng1
(5)
where J( n (t) are the gain (loss) fluxes per unit time. The master equation is “one step” in the sense that site n is connected only to adjacent sites in the lattice. Therefore, the gain flux at n is related to the loss flux
1 1 J+ n (t) ) Jn-1(t) + Jn+1(t) 2 2
(6)
It is simple to show that in Laplace space15 the loss flux is ˜ (s))/(1 simply related to the probability density, J˜n (s) ) s[(ψ ψ ˜ (s))]p˜ n(s), where s is the Laplace frequency conjugate to time t. ψ ˜ (s) is the characteristic function to the distribution of waiting times, ψ(t). Expanding the characteristic function, ψ ˜ (s) ∼ 1 (sτ)R, where 0 < R < 1, and τ is an intrinsic time scale. Typically τ is the time scale for molecular collisions, and for t/τ . 1 in the time domain
dpn(t) D ˆ 1-R ) (pn-1(t) + pn+1(t) - 2pn(t)), dt 2τR
ng2
(7)
obeys the differential equation
dN(t) ) -k p0(t) dt
For k * 0, the density of particles is not conserved. For these types of problems, one typically analyzes the distribution of first passage times. But because not every molecule reacts upon encounter, we have to instead examine the distribution of first reaction times. The number of particles to leave the system between time t and t + dt is related to the distribution of first reaction times, J(t)
J(t) ) -
D ˆ
1 d f(t) ) Γ(R) dt
f(t′)
∫0 dt′ (t - t′)1-R t
(8)
dp1(t) D Rˆ 1-R ˆ 1-R (p (t) 2 p (t)) + p0(t) ) 2 1 dt 2τR τR
(9)
dp0(t) D Rˆ ˆ p1(t) + k p0(t) ) R dt 2τ τR
(10)
(
1-R
)
s p˜ n(s) - pn(0) )
s (p˜ n+1(s) + p˜ n-1(s) - 2 p˜ n(s)) (14) 2(τ s)R
Let us think of p˜ n(s) as discrete samples from the continuous probability density, P ˜ (x, s), where x ) n and pn(t) ) P(x, t). In this picture, is as a “small” control parameter for a perturbation expansion
˜ (1)(x, s) + 2 P ˜ (2)(x, s) + ... P ˜ (x, s) ) P ˜ (0)(x, s) + P
1 d +k dt Γ (R)
(
(16)
s f s/λ
(17)
The continuous distribution in the coarse grained time is L (P(x, λt)) ) P ˜ (x, s/λ)/λ. Equation 14 now becomes an equation for P ˜ (x, s)
sP ˜ (x, s) - P(x, 0) )
λRs1-R (P ˜ (x + , s) + 2τR P ˜ (x - , s) - 2 P ˜ (x, s)) (18)
Expanding the right-hand side of eq 14 shows the relationship between λ and 2 (0) ˜ (x, s) 2 1-R ∂ P s R 2τ ∂2 x
(19)
Defining λ-R ≡ 2/(2τR) turns eq 19 into the fractional diffusion equation
-k(t-t′)
) ∫0t dt′ (te- t′)1-R f(t′)
(15)
t f λt
˜ (0)(x, s) - P(0)(x, 0)) ) λ-R(s P
where
Rˆ 1-R f(t) )
(13)
Next, we coarse grain time by a scale factor λ
The waiting time distribution at site n ) 0 is different than ψ(t) because at this site the molecule can undergo reaction. At this site, ψ(t) gets replaced by e-ktψ(t), the probability for not reacting and surviving until time t.15 When considering the flux back to site n ) 1 from the origin, we obtain the equations for the boundary
1-R
dN(t) dt
Equations 7, 9, and 10 can be analyzed to obtain the reactionsubdiffusion equations for targeting. A. Fokker-Planck Equation. To generalize the results of the previous section to higher dimensions and make the connection to reaction-subdiffusion equations, we now pass from the discrete space master equation to a continuum space Fokker-Planck equation. The lattice spacing, , provides a natural parameter for expansion, provided that the spatial variations of the density are slow on this length scale.16 In Laplace space, the master equation for sites n g 2 is
D ˆ 1-R is a Riemann-Liouville fractional derivative operator of order 1 - R that when acting on a function of time is 1-R
(12)
(11)
is the reactive flux operator. The fraction of particles that have remained in the system ∞ by time t is the survival probability, N(t) ) ∑n)0 pn(t). N(t)
∂ P(0)(x, t) ∂2P(0)(x, t) )D ˆ 1-R ∂t ∂2x
(20)
A similar equation holds for P(1). The evolution equation for P(2), however, contains higher-order spatial derivatives that account for more rapid spatial fluctuations. The two equations,
4286 J. Phys. Chem. B, Vol. 112, No. 14, 2008
Eaves and Reichman
eqs 9 and 10, form the boundary condition for the differential equation. For sites outside the reaction zone, we expect that spatial variations are slow compared to . The continuum solution cannot be extended to n ) 0 because at this site reaction and diffusion meet and can generate discontinuities. Taking this into account, the equations that determine the boundary condition, eqs 9 and 10 in coarse grained space and time are
(s P ˜ (, s) - P(, 0)) ) ˜ (2, s) - 2P ˜ (, s)) + 2(s + c)1-R p˜ 0(s)/ (21) s1-R(P 2
(s p˜ 0(s) + c p˜ 0(s) - p0(0)) ) ˜ (, s) - 2(s + c)1-R p˜ 0(s)/ (22) s1-R P In the above, the reaction rate constant in the coarse grained time is c ≡ kλ. The appearance of the small parameter in the denominator of eqs 21 and 22 makes the perturbation scheme trivially divergent for finite values of p˜ 0(0)(s); hence expanding to order O(-1), we find that p˜ 0(0)(s) ) 0. Expanding to O() yields the reflecting boundary condition
∂P(0)(x, s) ∂x
|
x)
)0
(23)
˜ (x, s)|x) ) [(s P ˜ (0) - P(0)(, 0)) - P(0)(, 0) + s1-R∂x P (s + c)Rs1-R P ˜ (0)(, s)] + O(2) (24) In the above, we have used eq 20 to eliminate the spatial second derivative and used the first-order result, p0(1)(0) ) limsf∞ s p˜ 0(1)(s) ) P(0)(, 0)/2. The first two terms grouped in parenthesis on the right-hand side of eq 24 are the Laplace transform of the time derivative of the probability at the boundary. The third term is the density at the boundary at time zero. The continuum solution should be independent of , so these terms must be zero. The first requirement, that the time derivative of the density is zero on the boundary, is consistent with our analysis of the chemical kinetics in steady state. The requirement that we do not begin with the system at the boundary is natural because this would impart rapid spatial and temporal fluctuations outside the scope of a simple continuum approximation. The fourth term contains the factor c and can be scaled so that it is independent of the lattice spacing. As the length over which absorption occurs tends to zero, it is natural to expect the absorption rate to become infinite so that the absorption rate per unit density is constant. The time scales that we are interested in are much larger than the absorption time (s , c), so we expand (s + c)R on the righthand side of eq 24 to obtain
|
x)
≈ cRP(, t)
(25)
or
2 ∂P(x, t) 2τR ∂x
|
x)
≈
DR
∂P(x, t) ∂x
|
x)0
) γRP(0, t)
(27)
It is the generalization of the radiation boundary condition in reaction-diffusion problems and is one principal result of this work. As we will explore later, this type of boundary condition contains symmetries that are not present the time non-local one analyzed by Sung et al.13 The first reaction distribution is the flux to the origin from molecules originating at J˜(s) ) DRs1-R∂x P ˜ (x, s)|x)0 R/2
J˜(z;κ) )
e-z 1 + (z/κ)R/2
(28)
where z ) sτD is the frequency normalized by natural units of time, τRD ) x02/DR, and κ is the ratio of the reaction rate to the diffusion rate to the reaction rate, κ ≡ τD((γR)/(xDR))2/R. Equation 28 can be expanded in a power series in terms of zR/2 ∞
At this order in the expansion, there is no feedback between reaction and diffusion. Expanding to O(2) and using the ˜ (x, z)|x) ) 0 + s1-R∂x P ˜ (1)(x, s)|x) perturbation series, s1-R∂x P + ..., we find
∂P(x, t) ∂x
holding DR ≡ 2/(2τR) and γR ≡ kR/2 constant, where DR is the subdiffusion coefficient and γR is the generalized bimolecular rate. This lowest-order non-trivial continuum approximation yields the boundary conditions Seki et al.14 have derived using a different method
kR P(, t) 2
(26)
After ordering the series in , we finally pass to the continuum limit by taking the limits that and τR both go to zero while
J˜(z;κ) )
∑ cn zRn/2 n)0
(29)
n where the coefficients cn are cn ) ∑k)0 ((-1)n)/((n - k)!κRk/2) with radius of convergence κ > z. For z near zero, J˜(z) ∼ 1 |c1|zR/2. To order zR/2, the exponent of the z dependence of this series is identical to the exponent of the first passage distribution, e-zR/2. Therefore, at long times J(t) has fat tails that go to zero as J(t) ∼ t-1-R/2. The walkers that have survived in the system to asymptotically long times have never seen the origin and the exponent of the asymptotic time dependence matches that of the first passage distribution. We emphasize that the continuum approximation is asymptotic and is most accurate when the (activation limited) reaction rates are faster than both 1/t and 1/τD. Before examining some results and implications for biological reactions in three spatial dimensions, we turn to numerical simulations of the master equation to check our asymptotic results. B. Numerical Simulations. We have simulated random walkers in one spatial dimension by using the “first reaction method” in continuous-time Monte Carlo.4 Consider a random walker on the lattice outside of the origin. Away from the origin, the walker executes a random walk, moving to the left or right with equal probability. During each step, the walker dwells for a duration chosen from the waiting time distribution, ψ(t). When the walker hits the origin, it may react with probability per unit time k or it may move back to site 1. At the origin, one draws two times, one from the distribution of reaction times, k e-kt, and one from the distribution of waiting times, ψ(t). We have chosen the one-sided Levy distribution, the limiting distribution for a sum of statistically independent processes that have tails that scale as ∼t-1-R, to generate the samples from ψ(t). If the time drawn for the reaction time is less than that drawn for the waiting time, then the walker exits the system. Otherwise, it moves back to site 1 and re-engages in the walk. In this way, one generates a “stochastically correct” chronological sequence of events for each walker.4 The density, pn(t) in the master equation, is the weight of many statistically independent walkers at time t, having initiated from the same initial position.
Subdiffusive Targeting Problem
J. Phys. Chem. B, Vol. 112, No. 14, 2008 4287
Figure 2. (A) Cartoon of our numerical procedure. Three walkers evolve by continuous-time Monte Carlo (see text). Each walker roams on the lattice until it reacts with the target and exits the system at time te. The survival probability (B) is the average of a large number of trajectories depicted in A. It is the probability that the walker remains in the system between time 0 and t. The probability per unit time of leaving the system is first reaction distribution, which is essentially the derivative of B and therefore the histogram of exit times. (C) The targeting problem in three dimensions. Because each walker is statistically independent and the motion is isotropic, the statistics of the arrival times only depend on the starting distance r0 relative to the reaction cross section, σ. In simulations it is easier to consider several independent trials, whereas in theory it is more convenient to analyze the density of walkers remaining at time t.
The fraction of walkers in the system at time t is
N(t) ) 〈θ (te - t)〉
(30)
where the 〈...〉 denotes an average over many random walks, θ(x) is the Heaviside step function, and te is the exit time. From eq 13, the first reaction probability is simply a histogram of the exit times
J(t) ) 〈δ(t - te)〉
(31)
J(t) is simple to compute. The survival probability is the cumulative time integral of the histogram of exit times and converges much more rapidly in simulations than J itself. Figure 2 is an illustration of the calculation. Figure 3 a shows the results of the calculation comparing two different starting lattice sites. In one set of trials, walkers begin at lattice site n ) 25, while in the other trial they start at n ) 50. We determined τD by running simulations on the entire real axis and determining DR from a plot of 〈x2〉t-R versus t. In both sets of trials, k > τD-1. C. Scaling. The statistics of the lattice random walkers are functions of the initial lattice site, n0, and the reaction rate, k. The continuum solution relies on limits taken to achieve an asymptotic approximation in space and time. As a result, there is no simple correspondence between the n0 and k parameters and the τ and κ parameters. Though τ and κ might be extracted by a fit, there is a better way to validate the continuum approximation. In the continuum limit, J(t) depends on the reaction time scale, k, and the initial position, x0, only through the dimensionless parameter κ. If the continuum approximation is accurate, then J(t) should depend on the starting lattice point, n0, and the reaction rate, k, through the parameter
Y ) n 0 kR
(32)
Figure 3. (A) Continuous-time Monte Carlo simulations for random walkers in one spatial dimension. The histogram of exit times, J(t), is the probability per unit time that a walker exits the system via reaction. The beginning lattice site is a factor of 2 different in these cases, but we have scaled k so that Y is constant. (B) Rescaled first reaction probabilities show excellent data collapse, in accordance with the continuum approximation, over several decades in time.
Equation 28 leads to the scaling hypothesis
J(t, k, n0) ) τD J (t/τD, Y )
(33)
where J is a scaling function. Different Y values describe different master curves for the reaction distributions and associated survival probabilities. More generally, lines in the (n0, k-R) plane that pass through the origin are families of parameters that describe statistically similar J(t). The curves in Figure 3A are simulated J(t) for starting positions that differ by a factor of 2. We have scaled k to keep Y constant for R ) 0.7. The results are in Figure 3 and exhibit excellent data collapse to a master curve over several decades in time. The first reaction distribution exhibits self-similarity, where a realization of one set of parameters can be mapped onto another one by scaling. A walker that begins close to the target but has a higher probability per encounter of reacting with the target has a statistically similar distribution of reaction times as a walker that begins far from the target but has a lower probability to react. This behavior is a direct consequence of the scale-free property of the waiting time distribution, ψ(at) ) a-1-Rψ(t). The specific values of k and τD describe quantitatively different curves with the same qualitative features: a rise from zero to a maximum and turn over into a power law with exponent J(t) ∼ t-1-R/2.
4288 J. Phys. Chem. B, Vol. 112, No. 14, 2008
Eaves and Reichman
III. Implications For Biochemistry In Vivo Golding and Cox9 have suggested that R ) 2/3 could be a critical value of the subdiffusion exponent because this value of R prevents the reactant and product from subdiffusing away from one another. In the previous section, we showed that finite reaction rates do not modify the reaction distribution qualitatively. It is therefore fruitful to analyze the κ f ∞ limit of J(t) or the first passage problem in three dimensions. If one is only interested in the distances that walkers will travel to encounter, then the problem is isotropic and is only a function of the radial distance between the particles. The walkers will either hit either the target whose reaction cross section with the walker is σ or it will wander off to infinity, never to find the target (Figure 2C). In three dimensions, the Laplacian in spherical coordinates replaces the second derivative in eq 20. The density of walkers obeys the fractional diffusion equation, eq 20, which in Laplace space is
g˜(r, s|r0) -
δ(r - r0) ) ξ˜ (s)2 ∂r2 g˜ (r, s|r0) 4π r0s
(34)
Here, g˜ (r, s|r0) ≡ rP(r, s|r0) is the radial density function, ξ˜ (s)
) xDR/sR is the frequency-dependent correlation length, and the initial condition describes a uniform distribution of particles initiating on a spherical shell a distance r0 - σ from the surface of the target. In first passage problems, g˜ (σ, s|r0) ) g˜ (R, s|r0) ) 0 where R is the escape radius, and R > r0 > σ. Solving eq 34 for the flux into σ after taking the limit that R f ∞ yields the distribution of first passage times
F ˜ (s) )
σ -(sτD)R/2 e r0
(35)
where F ˜ (s) is the characteristic function for the density of first passage times, F(t), and τR/2 ) (r0 - σ)/xDR. We are D interested in the fraction of random walkers that hit σ as time ˜ (s) ) (σ)/(r0). This is becomes infinite, which is limsf0 F independent of the subdiffusion exponent and is identical to the “diffusion to capture” result.18 At infinite times, the transient nature of the random walks in three dimensions determines the hitting probability. In dimension two and lower, all walkers eventually hit the surface and the survival probability tends to zero. But in three dimensions, there is a finite probability of never finding the target. In either the diffusive or subdiffusive case, the probability that a walker will visit the surface once is p ) (σ)/(r0). The probability that two particles meet once, return to their initial separation, and then wander away to infinity is p(1 - p). The probability of executing N round trips before wandering apart is pN(1 - p). The average number of round trips is then18 ∞
〈N〉 )
∑ npn(1 - p) n)0
(36)
Equation 36 can be summed exactly to yield 〈N〉 ) (σ)/(r0 σ), which tends to infinity as r0 tends toward σ. Following Berg’s argument,18 even if the surface of the target enzyme is patchy, a walker beginning close to the surface will visit it many times and will hit the reaction zone with high probability. Because p is independent of R, this argument is identical in the subdiffusion case. This gives physical insight into why the partially absorbing boundary is not qualitatively different from the totally absorbing boundary, even in the subdiffusive case.
Figure 4. Distribution of first passage times, or κ f ∞ limit of J(t) in one spatial dimension. In three spatial dimensions, these curves differ only by a constant amplitude factor, σ/r0 (see text). These functions are one-sided Levy distributions with index R/2 and have no critical point. Instead, the distributions of arrival times broaden as R decreases.
Figure 4 shows the first passage time distribution for various values of R for R f ∞. Note that the distribution of passage times is not qualitatively different for the diffusive case. If we consider the system to be in two states, solvated or bound, then J(t) is the waiting time distribution of the solvated state. For asymptotically long times, (t . τD), J(t) ∼ t-1-R/2, for example, in Figure 1B. The mean first passage time diverges even in the diffusive case. If the targeting molecules must travel a distance comparable to the diameter of the cell, as in Figure 1A, then the outer boundary condition rules the asymptotics and J(t) ∼ t-1-R. In this case, it can only be the fat-tailed distribution of waiting times that causes the mean first passage time to diverge.24 Whether or not a given reaction will ever achieve asymptotic behavior depends on several problem specific parameters, such as the lifetime of the proteins relative to the time required to reach the cell wall. If subdiffusion plays a significant role, then it must do so dynamically. The qualitatiVe difference between diffusion and subdiffusion is in the distribution of arrival times. Unfortunately, such dynamical effects are not apparent in experiments that look at the rested state of a genetic switch such as the lac operon.19 If a molecule passes through the cell wall and binds to a repressor molecule on the nucleoid, then it has a much broader distribution of arrival times if the transport is subdiffusive. This might increase the bandwidth of the switching times but should have little effect on experiments conducted in steady state.11,19 A recent study found that allowing transcription factors to diffuse significantly increased the intrinsic noise level of a gene circuit.20 When the motion is diffusive, the effect of diffusion can be modeled by renormalizing the binding affinity of the repressor and studying the well-stirred problem.20 Subdiffusion might change this picture qualitatively. The rules employed in our continuous-time Monte Carlo algorithm apply equally well to a multiple species reaction-diffusion system, as in the spatial Gillespie algorithm.21 In this algorithm, the system is decomposed into many small boxes with a size chosen so that on the length scale of an individual unit the system is well-stirred. Consider a molecule of RNA polymerase (RNAp) competing with a repressor for the operator binding site on the DNA in one of these boxes. In this case, binding takes place with constant probability per unit time and these times are drawn from a distribution with finite mean. The molecules will continue to compete for the site until one of the molecules diffuses to the next lattice site. Imagine that a value tr is drawn for the
Subdiffusive Targeting Problem RNAp binding time and tp is drawn for the repressor. The probability that either molecule will bind to the operator instead ∞ dt′ψ(t′), which increases draof diffusing away is ∫max(t r,tp) matically if ψ(t) has fat tails. Similarly, a multimeric repressor, such as the lac repressor, that can bind to several different sites of the DNA has a higher probability of staying attached in the subdiffusive case. These more complicated many-body effects are beyond the scope of the current manuscript but are the focus of future work. IV. Summary Motivated by enzyme kinetics and recent experiments in vivo, we have developed a model for the subdiffusive targeting problem. Our model begins with the development of a master equation with intuitive dynamical rules. The continuous time and space approximation to the lattice problem leads to a reaction-subdiffusion fractional differential equation with a mixed boundary condition, generalizing the radiation boundary condition in the diffusive targeting problem. Though scaling arguments9 and a time non-local boundary condition13 argue for a critical point at R ) 2/3, our results for the distribution of reaction times do not exhibit a critical point at any value of R.14 Instead, J(t) makes a smooth transition between diffusive and subdiffusive regimes. J(t) exhibits scaling, where the finite reaction rate and initial separation between the substrate and enzyme target can be scaled to collapse on a single master curve. The statistics of reaction times along lines that pass through the origin in the (n0, k-R) plane differ only by a change of scale. Because it appears that there is no critical value of R in isotropic three-dimensional targeting, it is not likely that evolution has driven the value of R close to 2/3. The argument pursued in ref 9 equates ((V/(DR3/2))(2/(3R)) with the aVerage time that a particle spends in a small volume of space V. Of course, the mean time that a particle spends in that region of space diverges in the subdiffusive problem. In both diffusive and subdiffusive transport, a protein locates its target with higher probability per unit time if it begins close to it. But it is likely that transcription factors and repressors lie close to their target sites to guard against spurious copying errors during the evolution, growth, and division of the bacterium. The subdiffusive physical environment in the interior of these cells is thus an indirect consequence of evolutionary pressures that incorporate a large and dense number of biochemical pathways into a small cellular volume. Acknowledgment. J.D.E. thanks Jun Shimada and Peter Mayer for stimulating discussions, as well as Michael Lomholt and Ralf Metzler for correspondence. We also thank J. Klafter
J. Phys. Chem. B, Vol. 112, No. 14, 2008 4289 for comments on the manuscript. The National Science Foundation (NIRT - 0210426) funded this work. References and Notes (1) Saunders, W. S.; Koshland, D.; Eshel, D.; Gibboms, I. R.; Hoyt, M. A. J. Cell Biol. 1995, 128, 617. (2) Lowrey, P. L.; Shimomura, K.; Antoch, M. P.; Yamazaki, S.; Zemenides, P. D.; Ralph, M. R.; Menaker, M.; Takahashi, J. S. Science 2000, 288, 483. (3) Kornberg, A. J. Bacteriol. 2000. (4) Gillespie, D. Annu. ReV. Phys. Chem. 2007, 58, 35. (5) Banks, D. S.; Fradin, C. Biophys. J. 2005, 89, 2960. (6) Metzler, R.; Klafter, J. Phys. Rep. 2000, 339, 1. (7) Wong, M.; ost Gardel, I.; Reichman, D. R.; Weeks, E. R.; Valentine, M.; Bausch, A.; Weitz, D. A. Phys. ReV. Lett. 2004, 92, 178101. (8) Weiss, M.; Elsner, M.; Kartberg, F.; Nilsson, T. Biophys. J. 2004, 87. (9) Golding, I.; Cox, E. C. Phys. ReV. Lett. 2006, 96, 098102. (10) Halford, S. E.; Marko, J. F. Nucleic Acids Res. 2004, 32, 3040. (11) Becskei, A.; Kaufmann, B. B.; van Oudenaarden, A. Nat. Genet. 2005, 37, 937. (12) Warren, P.; Rein ten Wolde, P. J. Mol. Biol. 2004, 342, 1379. (13) Sung, J.; Barkai, E.; Silbey, R. J.; Lee, S. J. Chem. Phys. 2002, 116, 2338. (14) Seki, K.; Wojcik, M.; Tachiya, M. J. Chem. Phys. 2003, 119, 7525. (15) Sokolov, I. M.; Schmidt, M. G. W.; Sagues, F. Phys. ReV. E 2006, 73, 031102. (16) van Kampen, N.; Oppenheim, I. J. Math. Phys. 1972, 13, 842. (17) Michaelis, L.; Menten, M. Biochem. Z. 1913, 49, 333. (18) Berg, H. C. Random Walks in Biology; Princeton University Press: Princeton, NJ, 1993. (19) Ozbudak, E. M.; Thattai, M.; Lim, H. N.; Shraiman, B. I.; van Oudenaarden, A. Nature 2004, 427, 737. (20) van Zon, J. S.; Morelli, M. J.; Tanase-Nicola, S.; Rein ten Wolde, P. Biophys. J. 2006, 4350-5367. (21) Andrews, S.; Arkin, A. Curr. Biol. 14, 16. (22) Lomholt, M. A.; Zaid, I. M.; Metzler, R. Phys. ReV. Lett. 2007, 98, 200603. (23) Assuming, of course, that the second cumulant of the jump length distribution is also finite. There is currently no experimental evidence to suggest otherwise. (24) In the process of completing this manuscript, we became aware of a similar work by Lomholt et al.22 Although our master equations are almost identical, our method of analysis leads to different conclusions. In particular, our eq 27 coincides with their eq 8 only if we allow P˜ (, s) to be nonanalytic, which would violate consistency. More specifically, both the boundary condition and subdiffusion equation are explicit asymptotic expressions for long time scales and large length scales and must be independent of the lattice spacing. The physical requirement that we do not allow random walkers to start on the origin arises naturally from the analysis of eq 24. However, the boundary conditions in eq 27, ref 22, and ref 14 are equivalent when walkers do not begin at the origin. Although our approach allows for the derivation of corrections to the Fokker-Planck equation from fast spatial fluctuations, as in ref 16, such corrections have little bearing on asymptotically slow dynamics. Indeed, the continuum limit is not useful if there are rapid spatial or temporal fluctuations. It is not presently clear how the differences in the boundary conditions would alter the picture of targeting in biological problems in the terms of weak ergodicity breaking, as presented in ref 22. This is a topic worthy of further investigation.
