Narrow Escape - Biologie ENS

3 downloads 408 Views 411KB Size Report
Sep 11, 2006 - The narrow escape problem consists in computing the mean time it takes a ... the time it activates a given protein on the cell membrane.
Narrow Escape: Theory and Applications to Cellular Microdomains D. Holcman ∗†, Z. Schuss‡, A. Singer§ September 11, 2006

Abstract We consider a Brownian particle (molecule, protein) confined to a bounded domain (a cell or a compartment) by a reflecting boundary, except for a small window through which it can escape. The narrow escape problem is to calculate the mean sojourn time before it escapes. This time diverges as the window shrinks, thus rendering the calculation of the mean escape time a singular perturbation problem. We summarize the various asymptotic formulas we obtained for several cases, including regular domains in two and three dimensions and some singular domains in two dimensions. The escape time is linked to many applications, because it corresponds to the mean time it takes for a molecule to hit a target binding site. We review applications in cellular biology.

1

Introduction

The narrow escape problem consists in computing the mean time it takes a diffusing particle to reach a small absorbing portion of an otherwise reflecting boundary of a given domain. As the absorbing boundary shrinks to zero, the ∗

Department of Mathematics, Weizmann Institute of Science, Rehovot 76100, Israel ([email protected]) † D´epartement de Math´ematiques et de Biologie, Ecole Normale Sup´erieure, 46 rue d’Ulm 75005 Paris, France. D. H. is supported by the program “Chaire d’Excellence”. This research is supported by the HFSP Grant 0007/2006-C. ‡ Department of Mathematics, Tel-Aviv University, Ramat-Aviv, Tel-Aviv 69978, Israel. This research was partially supported by a research grant from the US-Israel Binational Science Foundation ([email protected]) § Department of Mathematics, Program in Applied Mathematics, Yale University, 10 Hillhouse Ave. PO Box 208283, New Haven, CT 06520-8283 ([email protected])

1

mean time to absorption diverges to infinity, rendering the narrow escape a singular perturbation problem. Once formulated in terms of boundary value problems for partial differential equations, their singular perturbation analysis yields explicit asymptotic expressions for the mean escape time, depending on the diffusion coefficient, the ambient potential, dimensions, local and global geometrical properties of the domain and its boundary. The narrow escape problem is ubiquitous in cellular biology, because it concerns the random time between the release of a given particle in a cell and the time it activates a given protein on the cell membrane. In this context, the reciprocal of the escape time is the forward binding rate, used in homogeneous chemical reactions. It can also be used as a controlling rate in a Markov model of a chemical reaction. Recent experiments, using the patch clamp or calcium dyes, reveal some of the complexity of molecular signaling in cellular compartments. These are submicron domains that may contain only a small number of molecules, of the order between just a few and up to hundreds. This is the case in microdomains such as synapses, sensory compartments of cells such as outer segment of photoreceptors, and many others. At these low molecular concentrations, adding external chemical compounds, necessary for experimental purposes, alters the signaling pathway and thus modifies the physiological phenomenology. This calls for physical and mathematical modeling as a fundamental tool for both the quantitative and qualitative study of chemical reactions in microdomains. Indeed, the function of biological microdomains and specifically neurobiological microstructures is largely unknown and much effort has been expended for the last 20 years to unravel the molecular pathways responsible for the maintenance or modulation of cellular functions, and ultimately, to extract fundamental principles [1, 2, 3, 4, 5, 6]. Thus, it is still unclear what are the precise rules underlying memory at a synaptic level; we still don’t know how the synaptic weight is modulated by the neuronal activity [2, 3, 7]. It is also not clear how neurotransmitter release is modulated by synapses and how action potentials inside a dendrite are generated from local voltage summation. Cellular biology is rife with such examples. To obtain quantitative information about chemical processes occurring in microdomains, modeling and simulations seem to be inevitable. However, most of them are based on continuum concepts, which means that the medium is assumed homogeneous and the number of molecules involved is assumed sufficiently large, which is not the case here. In order to derive principles and to estimate the role of a microdomain, we have developed various computations to quantify precisely the role of the geometry in particle diffusion. These lead to the development of Markov models of chemical reactions and accounts for the small number of molecules involved. For example, the forward binding rate of a chemical reaction does not make 2

much sense when only a few molecules are involved. The widely used Smoluchowski formula k = 2πD[X] [8] for the binding rate k of Brownian particles with diffusion coefficient D and concentration [X], improved in [9, 10] assumes an infinite medium and spherical absorbing or partially absorbing surface. This computation does not account for the restricted geometry of microdomains. Our explicit computations give k for microdomains with various geometrical properties. specifically, we give several explicit analytical formulas for the mean time it takes a diffusing molecule to reach a small absorbing (binding) portion of the boundary of a microdomain. We discuss several applications to the forward binding rate to chemical reactions, receptor trafficking on the surface membrane of cells, and an application to vesicle trafficking in neuronal growth.

2

Formulation and mathematical results

We assume that a particle (e.g., receptor, molecule) in a biological medium is a Brownian particle in a field of force, so their motion can be described by the overdamped Langevin equation (known as the Smoluchowski limit), √ 1 ˙ x˙ + ∇φ(x) = 2D w, γ

(2.1)

where D=

kB T , mγ

(2.2)

γ is the dynamics friction coefficient, φ(x) is potential per unit mass, T is absolute temperature, kB is Boltzmann’s constant, m is the mass of the par˙ is a vector of independent δ-correlated Gaussian white noises, ticle, and w which represent the effect of the thermal motion and random collisions of the molecule with the molecules of the surrounding medium. The derivation of the Smoluchowski equation (2.1) is given in [12, 13, 14] for the three-dimensional motion of a molecule in solution, where Einstein’s equation (2.2) is valid. In two dimensions, such as for diffusion of a cylinder in the surface of a membrane, the diffusion coefficient is given by the Saffman-Delbrock formula [15] ¶ µ kB T µh D= (2.3) log 0 − γE 4πµh µa where µ is the viscosity of the solution, µ0 is the viscosity of the membrane, h the length of the cylinder, a is its radius, and γE is Euler’s constant 0.5772... . A generic problem in cellular biochemistry is to estimate the mean sojourn time of a Brownian particle in a bounded domain Ω, before it escapes through 3

a small absorbing window ∂Ωa in its boundary ∂Ω. The remaining part of the boundary ∂Ωr = ∂Ω−∂Ωa is assumed reflecting for the particle. The reflection may represent a high potential barrier on the boundary, or an actual physical obstacle. The opening may represent a narrow corridor in the barrier, or a defect in the physical obstacle (see figure 1). The biological interpretation of the mean sojourn time is discussed below. The escape time can be estimated asymptotically in the limit ε=

|∂Ωa | 0

(2.7)

∂pε (x, t) pε (x, t) ∂φ(x) + = 0, ∂n γ ∂n

for x ∈ ∂Ωr , t > 0,

(2.8)

where ρ0 (x) is the initial pdf (e.g, ρ0 (x) = 1/|Ω| for a uniform distribution or ρ0 (x) = δ(x − y), when the molecule is initially located at position y). In the later case, the function Z Z ∞ uε (y) = dx pε (x, t | y) dt, (2.9) Ω

0

where pε (x, t | y) is the pdf conditioned on the initial position, represents the mean conditional sojourn time. It is the solution of the boundary value problem (see [12]) D∆uε (y) − ∇φ(y) · ∇uε (y) = −1,

for y ∈ Ω

(2.10)

uε (y) = 0,

for y ∈ ∂Ωa

(2.11)

∂uε (y) = 0, ∂n

for y ∈ ∂Ωr .

(2.12)

4

Equation (2.12) is the adjoint boundary condition to (2.8). The survival probability is Z Sε (t) = pε (x, t) dx, (2.13) Ω

where

Z pε (x, t) =

pε (x, t | y)ρ0 (y) dy.

(2.14)



The density pε (x, t) can be computed by expanding in eigenfunctions pε (x, t) =

∞ X

ai (ε)ψi,ε (x)e−λi (ε)t ,

(2.15)

i=0

where λi (ε) (resp. ψi,ε ) are the eigenvalues (resp. eigenfunctions) of the Laplacian, and the coefficients ai (ε) depend on the initial function ρ0 (y). When ε > 0, the eigenvalues are strictly positive so the steady state is lim pε (x, t) = 0

t→∞

(the Brownian particle escapes in finite time with probability 1) and the mean sojourn time is asymptotically 1/λ0 (ε), when ε >

1 , λ1 (ε)

(2.16)

which shows that under the small hole approximation, the survival probability is exponential, Sε (t) ≈ e−λ0 (ε)t because

for t >>

Z

1 , λ1 (ε)

(2.17)

Z

a0 (ε) =

ρ0 (x) dx = 1, Ω

ψ0 (x) dx = 1,

(2.18)



due to normalization. Thus the stream of particles absorbed in the narrow opening is Poissonian. If the medium contains initially N0 independent particles, the population N (t) decays exponentially with rate λ0 (ε), N (t) = N0 Sε (t) ∼ N0 e−λ0 (ε)t

for ε ¿ 1,

t>

1 . λ1 (ε)

(2.19)

The asymptotic solution of equations (2.10)-(2.12) depends on the dimension (2 or 3) and the local geometry near the small opening [17, 18, 19, 20]. 5

When the boundary of the for ε ¿ 1 by      uε (y) =    

domain is regular, the escape time uε (y) is given |Ω| 1 ln + O(1) for n = 2 πD ε

(2.20)

|Ω| [1 + o(1)] for n = 3, 4aD where a represents the small radius of a geodesic disk located on the surface of the domain Ω. The function uε (y) does not depend on the initial position y, except for a small boundary layer near ∂Ωa , due to the asymptotic form found in [17, 18, 19, 20].

In dimension 2, the first order term matters, because, for example, if ε ≈ 1 10−1 , then ln ≈ 2.3 so the second term in the expansion (2.20) is comparable ε to the leading term. The second term can be found when Ω is a circular disk of radius R and the particle starts at the center [19, 20], as · ¸ R2 1 1 E[τ | x(0) = 0] = log + log 2 + + O(ε) . (2.21) D ε 4 The escape time, averaged with respect to a uniform initial distribution in the disk, is given by · ¸ R2 1 1 Eτ = log + log 2 + + O(ε) . (2.22) D ε 8 The geometry of the small opening can affect the escape time: if the absorbing window is located at a corner of angle α, then · ¸ |Ω| 1 Eτ = log + O(1) . (2.23) αD ε More surprising, near a cusp in a two dimensional domain, the escape time Eτ grows algebraically, rather than logarithmically: in the domain bounded between two tangent circles, the escape time is µ ¶ |Ω| 1 Eτ = + O(1) , (2.24) (d − 1)D ε where d > 1 is the ratio of the radii. Finally, when the domain is an annulus, the escape time to a small opening located on the inner circle involve a second R1 parameter which is β = < 1, the ratio of the inner to the outer radii, the R2 escape time, averaged with respect to a uniform initial distribution, is Eτ =

(2.25)

· ¸ (R22 − R12 ) 1 1 1 1 R22 2 log − R22 + O(ε, β 4 )R22 . log + log 2 + 2β + 2 D ε 21−β β 4 6

Equation (2.25) contains two terms of the asymptotic expansion of Eτ . The case β ≈ 1 remains open, and for general domains, the asymptotic expansion of the escape time remains an open problem. So does the problem of computing the escape time near a cusp point in three-dimensional domains.

3

Calcium diffusion in dendritic spines

Dendritic spines are small protrusions located on the surface of a neuronal dendrite; they receive most of the excitatory inputs. The physiological role of the spines is still unclear. The number and the shapes of the spines are highly correlated with cortical and synaptic plasticity [21, 22]. Calcium dynamics in dendritic spines is a fundamental signal, which can trigger physiological changes involved in remodeling the synaptic weight [4, 5]. Calcium diffusion involves many pathways (see figure 2): it can bind to calcium buffers or calcium binding molecules, enter the endoplasmic reticulum, or be pumped out through exchangers or pumps. The escape time determines much of the time course of calcium dynamics and its dependence on the geometry of a dendritic spine. In [23], we approximated the shape of a dendritic spine by a spherical head Ω connected to the rest of the dendritic shaft by a cylindrical neck of length L and radius a. In this approximation the radius of the neck is small relative to that of the spine head, so the mean time for a diffusing ions to escape a dendritic spine by diffusion alone can be decomposed into the mean time to find the small absorbing window (the neck) and the mean time to escape into the dendritic shaft. The former is approximated by formula (2.20), while the second time is L2 /D, which gives the escape time as τ≈

|Ω| L2 + . 4aD D

(3.1)

According to equation (2.17), the time course of calcium dynamics can be approximated by a single exponential with rate constant λ = 1/τ . The influence of pumps on calcium dynamics has been investigated in [24]. It was shown that they affect the dynamics by shifting the extrusion rate. To compute the shift, consider N pumps with an identical extrusion rate ξ, uniformly distributed along the neck. Then the total extrusion rate is λ≈

1 + ξN. τ

(3.2)

It is known that the geometry of the dendritic spines changes, depending on variables such as exposure to calcium concentration [25, 21, 26]. That is, dendritic spines can regulate dynamically their geometry, and possibly the 7

distributions and the number of pumps. Thus they regulate dynamically the fraction of calcium reaching the dendrite. The ratio of pumped to arriving ions can be evaluated according to several parameters. We demonstrated in [24] that the distribution of pumps along the spine neck can modulate such a ratio. On the other hand, given a fixed number of pumps, there exists a critical length, such that above it the spine head is effectively isolated, and below it, it conducts calcium [25]. Finally, the radius of the spine neck does not play a significant role in calcium diffusion inside a thin spine neck, when it is approximated as a small cylinder, because this parameter enters only in the second exponential in the sum (2.15), while the neck length enters directly the first one. Using the narrow escape computation, it is possible to compute the effect of crowding ions in a dendrite, which modifies the calcium time course by changing directly the expression of the first exponential decay rate.

4

Markov model of chemical reactions with a small number of molecules

In a closed microdomain Ω, a chemical reaction that involves only a few species, such as k1 M +S * ) SM, k−1

(4.1)

can be described in terms of a Markov process. We assume that M molecules are diffusing inside a domain Ω, and the substrate S consists of S0 binding molecules of small size a, located on the boundary and are well separated. 1 The the mean time for a diffusing molecule to bind is and depends only k−1 on the local potential [11], [12]. It is given by the Arrhenius law k−1 = Ce−∆E/kB Te ,

(4.2)

where C is a constant that depends on the temperature Te , the electrostatic potential barrier ∆E generated by the binding molecule, and the friction coefficient, while k1 represents the forward rate, which is given by the small window approximation as 4aS0 D . k1 = |Ω| The probabilities pk (t) = Pr{SM (t) = k} that there are exactly k SM 8

molecules produced at time t satisfy the master equations p˙k (t) = − (kk−1 + (S0 − k)k1 ) pk (t) + [k−1 (k + 1)] pk+1 (t) + [k1 (S0 − k + 1)] pk−1 (t) for k ≥ 1,

(4.3)

while for k = 0, p˙0 (t) = −S0 k1 p0 (t) + k−1 p1 (t) and for k = S0 , p˙S0 (t) = −S0 k−1 pS0 (t) + k1 pS0 −1 (t). These equations are derived by writing that the probability of formation of a new bond between times t and t + ∆t is k1 (S0 − k + 1) ∆t + o(∆t), while the probability of breaking a bond is k−1 (k + 1) ∆t + o(∆t). We get Pr{SM (t + ∆t) = k} = Pr{SM (t) = k + 1} [k−1 ∆t (k + 1)] + Pr{SM (t) = k − 1} [k1 ∆t (S0 − k + 1)] + Pr{SM (t) = k} (1 − kk−1 ∆t − (S0 − k)k1 ∆t) + o(∆t) for S0 > k > 1. The mean and the variance of pk are defined, respectively, as M (t) =

S0 X

kpk (t),

2

σ (t) =

S0 X

k 2 pk (t) − M 2 (t).

k=1

k=1

For example, the steady state mean and variance are computed by solving S0 X pk = 1. The directly the recurrence (4.3) with the normalization condition k=1

steady state probabilities are (CSk0 are the binomial coefficients) µ ¶k CSk0 k−1 pk = pk (∞) = µ , ¶S k−1 0 k1 1+ k1 and the moments are M (∞) =

S0 X

kpk (∞) = S0

k=1

k−1 = S0 k−1 + k1

k−1 4aS0 D k−1 + |Ω|

4aS0 D k−1 k1 |Ω| σ 2 (∞) = S0 = S0 µ ¶2 . 2 (k−1 + k1 ) 4aS0 D k−1 + |Ω| k−1

9

This framework is sufficiently general to include the Michaelis-Menten reaction theory and modulation of the reaction due to push-pull reactions. When the microdomain Ω is open and diffusing molecules arrive at random times, the time course of any chemical reaction inside Ω is affected, yet the properties of the moments can be estimated analytically [27].

5

Receptor trafficking on a neuronal membrane near a synapse: contribution to the synaptic weight

The synaptic weight between a pre- and a postsynaptic neurons depends in part on the number of postsynaptic receptors. In general, starting from a presynaptic neuron, the postsynaptic current depends on the type k, the number Nk , the conductance γk of the receptors, and the probability pk (t|f ) that channel k opens, conditioned on the firing rate f . The conditional probability pk (t|f ) depends on such diverse biophysical parameters as the presynaptic terminal dynamics and up to the gating properties of the channel, which depend on its molecular structure. The conductivity γk controls the flux of ions through a channel of type k and it depends mainly on the molecular properties of the primary subunit structure. The postsynaptic current is given by X γk Nk pk (t|f ). (5.1) I(t) = k

There is experimental evidence that channels are not static, so the number of channels Nk in a given domain is not a constant, but rather fluctuates in time. Glutamatergic and GABA-ergic receptors on the surface of neurons, as well as other receptors, traffic in and out of a fundamental microstructure called the Postsynaptic Density (PSD) [28, 29, 2]. The receptor movement has been approximated so far as mostly Brownian, with measurable sojourn times (confinements) in small subdomains (see figure 3). When the confining restriction is due to a corral zone around the receptor, the mean escape time, given by formula (2.20), is a good approximation of the measured confinement time [27]. However, inside the PSD region, receptors can be anchored to the membrane when they bind to scaffolding proteins. Figure 4 represents a coarsegrained model of the PSD into three compartments, where the free binding sites have been combined into the green annulus, the bound sites - into the orange disk, the grey annulus represents the domain of free diffusion of the receptors, and the scaffolding molecules, as well as many other structures, such as transmembrane molecules, sub-membranous cytoskeleton constituting obstacles and fences are grouped into the brown boundary. The total number 10

of scaffolding molecules (orange plus green) is constant, but the proportion of bound and unbound depends on the number of receptors in the PSD. The peripheric brown fence is connected to the extra-synaptic region through a small hole that restricts the dynamics of the receptors. A random trajectory of a receptor (blue) has been drawn in both pictures (black broken line). K1 is the forward binding rate of a receptor to the scaffolding molecules (which depends on the total number of scaffolding molecules and the mean time it takes to enter this domain), K−1 is the backward binding rate. The Markov model, combined with the narrow escape approximation, can account for the interactions of receptors with the scaffolding molecules and give an estimate of the mean and variance of the number of bound receptors at synapses, as well as of the chemical reactions described in the previous paragraph [30]. The computation is based on four assumptions: 1) the stream of receptors entering the synaptic region is Poissonian with rate equal to the inward flux J. 2) The escape rate of receptors from the PSD is the recip|Ω| 1 rocal of the escape time τ ∼ ln , given by eq.(2.20). This assumption πD ε is motivated by the fact that transmembrane molecules, such as the one involved in adhesion, or such that bind receptors, may act as pickets, and that sub-membranous molecules can create a fence (see figure 4) [31]. Finally, 3) a receptor can bind a free scaffolding molecule at rate k2 , according to the standard law of chemical reactions, and 4) they can dissociate from a scaffolding molecule with a rate k−2 . Using assumptions 1-4, the steady state Markov model balances the flux J with the number of escaping receptors. The number of bounded receptors and their mean can be estimated from the model as a function of the flux (see figure 5) [30]. The number of bound receptors Rb is given by Rb =

Jin k2 τ1 S0 . |Ω| 1 Jin k2 ln + k−2 πD ε

(5.2)

Recent experimental findings (Choquet, Triller personal communications) have demonstrated that the receptor movement is modulated by neuronal activity and depends on the type of receptors, which results in a modulation of the flux Jin .

6

Residence time of a receptor in the synapse

The calculation of the mean time a Brownian molecule spends inside a microdomain Ω before it escapes through one of the small holes on the boundary, when it can be caught and released by a scaffolding molecule inside Ω, is a generalization of the narrow escape problem described above. This mean 11

time depends on several parameters such as the backward binding rate (with the agonist molecules), the mean escape time from the microdomain and the mean time it takes for the molecule to reach the binding sites (forward binding rate). This mean time is usually called the dwell time and experimental measurements, such as FRAP (Fluorecent Recovery After Photobleaching) or FCS (Fluorescence correlation spectroscopy), are used for its estimate. We assume for simplicity that the microdomain Ω can be decomposed into a part where the receptor is diffusing freely and a part where it is bound to the scaffolding molecule (see figure 3). Under these assumptions the dwell time can be computed by counting all the possibilities of a Brownian molecule to exit after 0,1,... bindings. Summing the probabilities of all these events, we obtain a geometric series and the dwell time E(τD ) is calculated to be µ ¶ 1 − mδ 1 E(τD ) = hτ i + hT i + , (6.1) mδ k−1 where hτ i is the mean time to exit when no binding occurs, hT i is the mean time to enter the binding site area, mδ is the probability to bind before exit, given a uniform initial distribution, and k−1 is the backward binding rate. An explicit expression for the dwell time formula (6.1) can be obtained in term of the geometry, by estimating asymptotically the probability pδ (x) that a Brownian molecule starting at x exits before entering a binding disk of radius δ. The estimate was obtained so far when the binding domain is a disk of radius δ and the free diffusing space is an annulus Aδ of outer radius R. This estimate was obtained under the assumption that the ratio β = δ/R ¿ ε, where ε is the ratio of the small opening to the total length of the boundary. The other cases remain open [32]. The results are Z 1 pδ (r, θ)r dr dθ log β A(δ) Z mδ = hpδ i = = + o(1). 1 1 log + 2 log + 2 ln 2 r dr dθ β ε A(δ)

12

Finally the dwell time formula is µ ¶ 1 − mδ 1 E(τD ) = hτ i + hT i + mδ k−1 (6.2) µ µ ¶ ¶ R 1 = ln + ln 2 + D ε µ ¶   1 ¾ ½ 2 ln + 2 ln 2  1  1 R2 ε  + o(1).  µ ¶ µ ¶ + log  1 1 β k−1  ln 16D log + ln 2 β ε 2

In addition, the mean Mb and the variance Vb of the number of bounds made by a single molecule before it exits Ω are given by 1 1 − mδ ε + O(1) (6.3) = = 1 mδ log β µ ¶  µ ¶ 1 1 µ ¶ µ ¶ ln + 2 ln 2 + 2 ln 1 − mδ 1  ε β    (6.4) = ≈ 2 ln 2   1 mδ ε log2 β 2 log

Mb

Vb

These expressions are valid for ε