A solvable model for the diffusion and reaction of neurotransmitters in a synaptic junction

Jorge L Barreda; Huan-Xiang Zhou

doi:10.1186/2046-1682-4-5

Research article
Open Access
Published: 02 March 2011

A solvable model for the diffusion and reaction of neurotransmitters in a synaptic junction

Jorge L Barreda¹ &
Huan-Xiang Zhou¹

BMC Biophysics volume 4, Article number: 5 (2011) Cite this article

5674 Accesses
4 Citations
0 Altmetric
Metrics details

Abstract

Background

The diffusion and reaction of the transmitter acetylcholine in neuromuscular junctions and the diffusion and binding of Ca²⁺ in the dyadic clefts of ventricular myocytes have been extensively modeled by Monte Carlo simulations and by finite-difference and finite-element solutions. However, an analytical solution that can serve as a benchmark for testing these numerical methods has been lacking.

Result

Here we present an analytical solution to a model for the diffusion and reaction of acetylcholine in a neuromuscular junction and for the diffusion and binding of Ca²⁺ in a dyadic cleft. Our model is similar to those previously solved numerically and our results are also qualitatively similar.

Conclusion

The analytical solution provides a unique benchmark for testing numerical methods and potentially provides a new avenue for modeling biochemical transport.

1. Background

In intercellular and intracellular spaces, passive transport of biomolecules is a common phenomenon. Because such processes are difficult to probe directly by experiments, numerical modeling is increasingly used to gain insight. Two processes that have been extensively modeled are the diffusion and reaction of the transmitter acetylcholine in a neuromuscular junction [1–6] and the diffusion and binding of Ca²⁺ in the dyadic cleft of a ventricular myocyte [7, 8]. In contrast to previous numerical approaches, here we present an analytical solution of a model for the diffusion and reaction of acetylcholine in a synaptic cleft (or Ca²⁺ in a dyadic cleft). Our model is similar to those previously solved numerically; hence our analytical solution potentially provides a new avenue for modeling biochemical transport. More importantly, an analytical solution provides a unique benchmark for testing numerical methods. Such a solution has been lacking up to now; the present work fills this gap.

Neuromuscular junction refers to the cleft between a motor neuron and a muscle fiber. As illustrated in Figure 1, the neuronal signal for muscle contraction is mediated by acetylcholine. These neurotransmitter molecules are initially inside vesicles located in the pre-synaptic axon terminal. When an action potential reaches the axon terminal, the vesicles release acetylcholine molecules into the synaptic cleft. These molecules then diffuse toward the post-synaptic membrane and bind to acetylcholine receptors in the membrane. Acetylcholine binding activates these ligand-gated ion channels, allowing Na⁺ to flow in and generating an action potential along the muscle fiber. Finally excess acetylcholine molecules around the post-synaptic membrane are broken down by acetylcholinesterase to prevent continued activation of acetylcholine receptors.

A related system is a dyadic cleft, which spans the gap between the cell membrane in a transverse tubule and the membrane of a sarcoplasmic reticulum. As Figure 2 illustrates, Ca²⁺ can enter the cell through L-type Ca²⁺ channels on the cell membrane in response to the arrival of an action potential. The ions then diffuse to reach and activate ryaonodine receptors in the membrane of the sarcoplasmic reticulum. The activated ryaonodine receptors release Ca²⁺ from the sarcoplasmic reticulum, which ultimately lead to muscle contraction.

Here we propose a simple but not unrealistic model for the diffusion and reaction of acetylcholine in a neuromuscular junction. The model also applies to the diffusion and binding of Ca²⁺ in a dyadic cleft. We are able to find an analytical solution for this model. For convenience we describe our model in the language of neuromuscular junction. As shown in Figure 3, we model the synaptic cleft as the space between two infinite, parallel planar membranes. In the pre-synaptic membrane, there is a periodic array of circular openings, via which acetylcholine molecules enter the cleft. In the post-synaptic membrane, there is a periodic array of circular disks, where the acetylcholine molecules are absorbed. The quantity of interest is the total flux, J(t), at time t of acetylcholine molecules across the post-synaptic membrane. This model allows us to use periodic boundary conditions in the transverse direction. Our results are qualitatively similar to those obtained previously by numerical solutions [1–7]. More realistic ingredients can be added to our model and still permit analytical solutions.

2. Methods

We set up a coordinate system such that the x and y axes are parallel to the pre- and post-synaptic membranes. The synaptic junction has depth L_z . The junction is periodic in the x and y directions, with periodicities of L_x and L_y , respectively. In each "unit cell", a synaptic vesicle bursts at time t = 0, releasing the neurotransmitters into the cleft. We model the release of the neurotransmitters as a transient flux, u(t), that is confined to a circular opening with radius R. We place the synaptic vesicle at the center of the pre-synaptic face of the unit cell. After diffusing to the post-synaptic membrane, neurotransmitters are absorbed by a circular disk in each unit cell, with radius a. We place this "sink" also at the center of the post-synaptic face of the unit cell. The exact shapes and locations of the pre-synaptic opening and the post-synaptic sink are not essential for the analytical solution of our model. The quantity of interest is the total flux, J(t), through the post-synaptic face of each unit cell.

We place the origin of the Cartesian coordinate system at the center of the pre-synaptic face, with the z axis pointing toward the post-synaptic face. The concentration of neurotransmitters at position r and at time t is C(r, t). Within the synaptic junction, it is governed by the diffusion equation,

∂ C (r, t) ∂ t = D ∇ 2 C (r, t)

(1)

where D is the diffusion constant. The initial condition is

C (r, 0) = 0

(2)

The boundary condition at z = 0 is

D ∂ C (r, t) ∂ z = u (t) when ρ < R

(3a)

D ∂ C (r, t) ∂ z = 0 when ρ > R

(3b)

where ρ = (x² + y²)^1/2 is the distance to the z axis. The boundary condition at z = L_z is

C (r, t) = 0 when ρ < a

(4a)

D ∂ C (r, t) ∂ z = 0 when ρ > a

(4b)

We solve the problem in Lapace space. For a function f(t) of t, we denote the Laplace transform as $f^(s) ≡ ∫ 0 ∞ d t e − s t f (t)$ . The Laplace transform of Eq. (1), using the initial condition of Eq. (2), is

s C^(r, s) = D ∇ 2 C^(r, s)

(5)

The solution appropriate for the periodic boundary conditions in the x and y directions has the form

C^(r, s) = ∑ l, m = 0 ∞ [α l m (s) e γ l m (s) z + β l m (s) e − γ l m (s) z] cos (2 l π x / L x) cos (2 m π y / L y)

(6)

where

γ l m (s) = [s / D + (2 l π / L x) 2 + (2 m π y / L y) 2] 1 / 2

(7)

The coefficients α_lm (s) and β_lm (s) are to be determined from boundary conditions.

To make use of the boundary condition of Eq. (3), we need the 2-dimensional cosine transform of a function in x and y that has value 1 when ρ < R and value 0 when ρ > R. The coefficient of the cos(2lπ/L_x )cos(2mπ/L_y ) term is

q l m = 4 ε l ε m L x L y ∬ ρ < R d x d y cos (2 l π x / L x) cos (2 m π y / L y)

(8)

where ε₀ = 1/2 and ε_l = 1 for l > 0. For later reference, we note that the function having value 1 when ρ < a and value 0 when ρ > a has the following coefficients in its 2-dimensional cosine transform:

p l m = 4 ε l ε m L x L y ∬ ρ < a d x d y cos (2 l π x / L x) cos (2 m π y / L y)

(9)

Applying the boundary condition of Eq. (3), we find

D γ l m (s) [α l m (s) − β l m (s)] = u ˆ (s) q l m

(10)

For the boundary condition of Eq. (4) at z = L_z , we use the constant-flux approximation [9]:

D ∂ C^(r, s) ∂ z = Q (s) when ρ < a

(11a)

D ∂ C^(r, s) ∂ z = 0 when ρ > a

(11b)

where the quantity Q(s) is to be determined by the condition

< C^(r, s) > ≡ (π a 2) − 1 ∬ ρ < a d x d y C^(r, s) = 0

(12)

Equation (11) leads to

D γ l m (s) [α l m (s) e γ l m (s) L z − β l m (s) e − γ l m (s) L z] = Q (s) p l m

(13)

With Eqs. (10) and (13), we solve for coefficients α_lm (s) and β_lm (s). The results are

α l m e γ l m (s) L z = 1 D γ l m (s) − q l m u ˆ (s) + p l m Q (s) e γ l m (s) L z e γ l m (s) L z − e − γ l m (s) L z

(14a)

β l m e − γ l m (s) L z = 1 D γ l m (s) − q l m u ˆ (s) + p l m Q (s) e − γ l m (s) L z e γ l m (s) L z − e − γ l m (s) L z

(14b)

Inserting Eq. (6) in Eq. (12) and using Eqs. (14) and (9), we find Q(s) as

Q (s) = u ˆ (s) ∑ l, m = 0 ∞ p l m q l m ε l ε m γ l m (s) sinh [γ l m (s) L z] ∑ l, m = 0 ∞ p l m 2 cosh [γ l m (s) L z] ε l ε m γ l m (s) sinh [γ l m (s) L z]

(15)

Finally the total flux through the sink on the post-synaptic face is

J^(s) = π a 2 Q (s)

(16)

The flux accumulated over all times,

I = ∫ 0 ∞ d t J (t)

(17a)

is of interest. In our model,

I = J^(0) = π R 2 u^(0)

(17b)

which is independent of a. Equation (17b) is simply a consequence of ligand conservation, i.e., the total number of ligands released from the synaptic vesicle is the same as the total number of ligands absorbed by the receptors.

The analytical solution derived above can be implemented on any function u(t) modeling neurotransmitter release from the synaptic vesicle. We focused on an exponentially decaying function:

u (t) = e − t / t 0

(18a)

Its Laplace transform is

u^(s) = 1 s + 1 / t 0

(18b)

Our model then contains two parameters related to time: D and t₀. For two sets of these parameters, e.g., D₁ and t₀₁ in one and D₂ and t₀₂ in the other, it can be shown that the corresponding response functions satisfy

J 1 (D 2 t) = J 2 (D 1 t) when D 1 t 01 = D 2 t 02

(19)

We now briefly describe the details of our implementation of the analytical solution. To calculate the q_lm and p_lm coefficients of Eqs. (8) and (9), we first carried out the integration over y analytically. The remaining integration over x was done numerically using the Gauss-Legendre quadrature with two points. The summations over l and m in Eq. (15) were truncated at l = m = 40. The Laplace transform of $J^(s)$ was inverted by the Stehfest algorithm [10]. A Fortran90 code for the implementation is available upon request.

3. Results

We now present some illustrative results. The parameters of our model are as follows: L_x = L_y = 500 nm; L_z = 50 nm; R = 20 nm; a varied from 2.5 to 40 nm; t₀ varied from 1 to 10 ms; and D varied in (0.4-4) × 10⁵ nm²/ms.

Our single absorbing disk is used to model ligand binding to multiple receptors. Increasing a thus mimics an increasing number, N_rec, of receptors per unit cell. We expect a to increase linearly with increasing N_rec when N_rec is small; the increase in a then slows down at higher N_rec (see Discussion). Figure 4 displays the dependence of the response function J(t) on a (and hence N_rec). As a (i.e., N_rec) increases, J(t) more and more quickly reaches a higher and higher maximum and decays faster and faster. This is the expected behavior. Moreover, according to Eqs. (17b) and (18b), the areas under the curves for different a values are all equal to πR²t₀. Note also that all our J(t) curves have the familiar shape seen in previously numerical studies [1–7].

Figure 5 shows the change in the response function when the speed at which the neurotransmitters are released into the synaptic cleft is varied. To make a fair comparison, J(t) is scaled by t₀ so that, effectively, the total number of neurotransmitters entering the synaptic cleft is fixed. It is clear that, as t₀ increases (i.e., as the speed of neurotransmitter release decreases), the response function rises and decays more slowly. This behavior has previously been specifically modeled by Stiles et al. [1].

There is some uncertainty on the diffusion constant of acetylcholine in neuromuscular junctions. In previous models [1–6], the value of D ranged from (0.25-6) × 10⁵ nm²/ms. In Figure 6 we display the change in the response function when D is varied from (0.4-4) × 10⁵ nm²/ms. As expected, when neurotransmitter diffusion slows down, the response function is also delayed. In our model, decreasing D has the same effect on the response function as increasing t₀ [see Eq. (19)].

4. Discussion and Conclusion

We have presented an analytical solution to a model for the diffusion and reaction of acetylcholine in a neuromuscular junction. The model also applies to the diffusion and binding of Ca²⁺ in a dyadic cleft. Our results are qualitatively similar to those obtained previously from models solved numerically [1–7].

Perhaps the greatest value of our analytical solution is that it provides a benchmark for testing numerical methods. Diffusion and reaction of ligands in intercellular and intracellular spaces have been modeled either on a particle description or a concentration description. The former type of models have been solved by Monte Carlo simulations [1, 7], while the latter type of models have been solved by either finite-difference [2, 6] or finite-element [3–5] methods. The two types of models have been shown to give equivalent results [8]. The level of realism of our model approaches those of the models solved numerically; hence our analytical solution will be able to serve as a good benchmark for the numerical methods.

Our model has room for increasing the level of realism and still allows for analytical solution. For example, we modeled ligand binding to receptors as absorbing. The binding can be modeled as partially absorbing if binding does not occur at every ligand-receptor encounter. The partial absorption condition takes the form

D ∂ C (r, t) ∂ z = κ C (r, t) when ρ < a

(20)

where the reactivity κ controls the degree of partial absorption. When κ = 0, the region ρ < a becomes reflecting and binding cannot occur. When κ → ∞, the region ρ < a becomes fully absorbing and Eq. (20) reduces to the absorbing condition of Eq. (4a). Note that the bimolecular rate constant k for binding to the partially absorbing disk is given by [11]

1 k = 1 4 D a + 1 π κ a 2

(21)

One can thus parameterize a and κ by matching k with experimental data for the ligand-receptor binding rate constant.

Another simplification of our model is that a single absorbing disk is used to represent the receptors. We accounted for the presence of multiple receptors per unit cell by increasing the radius a of this disk. One can identify a by requiring that the bimolecular rate constant k calculated for the single absorbing disk is the same as that for the N_rec receptors. For the single absorbing disk we have k = 4Da [see Eq. (21)]. If each of the N_rec receptors is modeled as an absorbing disk with a small radius a₀, then k ≈ 4N_recDa₀ when N_rec is small [12]. Therefore a ≈ N_reca₀ for small N_rec. As N_rec increases, the rate constant k for the multiple receptors and correspondingly the radius a of the equivalent disk also increase, but the increase slows down as N_rec becomes large [12, 13]. Instead of using a single equivalent absorbing disk, analytical solution is actually still permitted when multiple receptors, each represented by a (partially) absorbing small disk, are present at arbitrary positions on the post-synaptic face. An alternative way to model the presence of multiple receptors is to assume that the whole post-synaptic face is partially absorbing, with an reactivity given by [12, 13]

K = N rec k 0 S (1 − f rec)

(22)

where S = L_xL_y is the area of the post-synaptic face, f_rec = N_recπa₀²/S is the fraction of the post-synaptic face that is covered by the receptors, and k₀ is given by

1 k 0 = 1 4 D a 0 + 1 π κ a 0 2 (1 − f rec)

(23)

We have modeled ligand-receptor binding as irreversible. This is somewhat justified for modeling the neuromuscular junction, in which acetylcholinesterase can break down acetylcholine molecules newly released from the receptors. No such mechanism is present for Ca²⁺ in the dyadic cleft. Reversible binding can be treated by appropriate boundary conditions [14] on the post-synaptic face. Another important detail is that both acetylcholine and ryanodine receptors have multiple binding sites for their ligands so that there are multiple ligand-occupation states for the receptors. Again, these can be treated by appropriate boundary conditions.

The geometries of some of the models previously solved numerically are more sophisticated than that of our model. In particular, secondary folds of the neuromuscular junction has been included in some of the previous models [1, 3–5]. A formalism for treating ligand binding to a site buried in a narrow tunnel has been developed [15] and may be adopted for treating the narrow secondary folds in the neuromuscular junction. However, analytical solution requires idealized geometries; the kind of realistic geometries drawn from electron microscopy that can be handled by a finite-element method [4, 5] is beyond the reach of analytical solution. Nevertheless, with all the new ingredients outlined above, analytical solution will potentially provide a new avenue for modeling biochemical transport.

References

1.
Stiles JR, Van Helden D, Bartol TM, Salpeter EE, Salpeter MM: Miniature endplate current rise times less than 100 microseconds from improved dual recordings can be modeled with passive acetylcholine diffusion from a synaptic vesicle. Proc Natl Acad Sci USA. 1996, 93: 5747-5752. 10.1073/pnas.93.12.5747.
Article Google Scholar
2.
Naka T, Shiba K, Sakamoto N: A two-dimensional compartment model for the reaction-diffusion system of acetylcholine in the synaptic cleft at the neuromuscular junction. Biosystems. 1997, 41: 17-27. 10.1016/S0303-2647(96)01659-0.
Article Google Scholar
3.
Smart JL, McCammon JA: Analysis of Synaptic Transmission in the Neuromuscular Junction Using a Continuum Finite Element Model. Biophysical journal. 1998, 75: 1679-1688. 10.1016/S0006-3495(98)77610-6.
Article Google Scholar
4.
Tai K, Bond SD, MacMillan HR, Baker NA, Holst MJ, McCammon JA: Finite element simulations of acetylcholine diffusion in neuromuscular junctions. Biophys J. 2003, 84: 2234-2241. 10.1016/S0006-3495(03)75029-2.
Article Google Scholar
5.
Cheng Y, Suen JK, Radic Z, Bond SD, Holst MJ, McCammon JA: Continuum simulations of acetylcholine diffusion with reaction-determined boundaries in neuromuscular junction models. Biophys Chem. 2007, 127: 129-139. 10.1016/j.bpc.2007.01.003.
Article Google Scholar
6.
Khaliq A, Jenkins F, DeCoster M, Dai W: A new 3 D mass diffusion-reaction model in the neuromuscular junction. Journal of Computational Neuroscience. 2010, 1-17.
Google Scholar
7.
Koh X, Srinivasan B, Ching HS, Levchenko A: A 3 D Monte Carlo analysis of the role of dyadic space geometry in spark generation. Biophys J. 2006, 90: 1999-2014. 10.1529/biophysj.105.065466.
Article Google Scholar
8.
Hake J, Lines GT: Stochastic binding of Ca2+ ions in the dyadic cleft; continuous versus random walk description of diffusion. Biophys J. 2008, 94: 4184-4201. 10.1529/biophysj.106.103523.
Article Google Scholar
9.
Shoup D, Lipari G, Szabo A: Diffusion-controlled bimolecular reaction rates. The effect of rotational diffusion and orientation constraints. Biophysical Journal. 1981, 36: 697-714. 10.1016/S0006-3495(81)84759-5.
Article Google Scholar
10.
Stehfest H: Algorithm 368: Numerical inversion of Laplace transform. Communication of the ACM. 1970, 13: 47-49. 10.1145/361953.361969.
Article Google Scholar
11.
Zhou H-X, Szabo A: Theory and simulation of the time-dependent rate coefficients of diffusion-influenced reactions. Biophysical Journal. 1996, 71: 2440-2457. 10.1016/S0006-3495(96)79437-7.
Article Google Scholar
12.
Shoup D, Szabo A: Role of diffusion in ligand binding to macromolecules and cell-bound receptors. Biophysical Journal. 1982, 40: 33-39. 10.1016/S0006-3495(82)84455-X.
Article Google Scholar
13.
Zwanzig R: Diffusion-controlled ligand binding to spheres partially covered by receptors: an effective medium treatment. Proc Natl Acad Sci USA. 1990, 87: 5856-5857. 10.1073/pnas.87.15.5856.
Article MathSciNet Google Scholar
14.
Gopich IV, Szabo A: Kinetics of reversible diffusion influenced reactions: the self-consistent relaxation time approximation. Journal of Chemical Physics. 2002, 117: 507-517. 10.1063/1.1482701.
Article Google Scholar
15.
Zhou H-X: Theory of the diffusion-influenced substrate binding rate to a buried and gated active site. Journal of Chemical Physics. 1998, 108: 8146-8154. 10.1063/1.476255.
Article Google Scholar

Download references

Acknowledgements

This work was supported in part by Grant GM58187 from the National Institutes of Health.

Author information

Affiliations

Department of Physics, Institute of Molecular Biophysics, Tallahassee, Florida, 32306, USA
Jorge L Barreda & Huan-Xiang Zhou

Authors

Jorge L Barreda
View author publications
You can also search for this author in PubMed Google Scholar
Huan-Xiang Zhou
View author publications
You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Huan-Xiang Zhou.

Additional information

Authors' contributions

JLB implemented the analytical solution, did the calculations, and prepared the figures. HXZ derived the analytical solution and wrote the paper. Both authors read and approved the final manuscript.

Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

Authors’ original file for figure 1

Authors’ original file for figure 2

Authors’ original file for figure 3

Authors’ original file for figure 4

Authors’ original file for figure 5

Authors’ original file for figure 6

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and Permissions

About this article

Cite this article

Barreda, J.L., Zhou, H. A solvable model for the diffusion and reaction of neurotransmitters in a synaptic junction. BMC Biophys 4, 5 (2011). https://doi.org/10.1186/2046-1682-4-5

Download citation

Received: 09 February 2011
Accepted: 02 March 2011
Published: 02 March 2011
DOI: https://doi.org/10.1186/2046-1682-4-5

Keywords

Acetylcholine
Sarcoplasmic Reticulum
Synaptic Vesicle
Neuromuscular Junction
Synaptic Cleft