| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |





* Department of Chemistry and Biochemistry,
Department of Mathematics,
Department of Pharmacology, and the Howard Hughes Medical Institute, University of California, San Diego, La Jolla, California 92093; and
Department of Biochemistry and Molecular Biophysics and the Center for Computational Biology, Washington University School of Medicine, St. Louis, Missouri 63110
Correspondence: Address reprint requests to Kaihsu Tai, Dept. of Chemistry and Biochemistry, University of California at San Diego, 9500 Gilman Drive, La Jolla, CA 92093-0365. Fax: 858-534-7042; E-mail: ktai{at}mccammon.ucsd.edu.
| ABSTRACT |
|---|
|
|
|---|
| INTRODUCTION |
|---|
|
|
|---|
| The neuromuscular junction |
|---|
|
|
|---|
ACh diffuses across the synaptic cleft and potentially binds to acetylcholine receptors, ion channels embedded in the postsynaptic folds. AChR ion channels are multiprotein membrane-spanning complexes found at packing densities of up to 10,000 µm-2 at the crests of the postsynaptic folds. In the absence of ACh, an AChR channel is impermeable to ion flow. However, once two ACh molecules bind to it, AChR opens and ions flow through the muscle cell membrane: sodium inward, potassium outward. This ion flow defines an endplate current (EPC) across the muscle membrane that, when strong enough, induces contraction. Roughly 1 ms after activation, ACh molecules are released back into solution and AChR ion channels close.
ACh molecules remain in the cleft until they are hydrolyzed by AChE, the biomolecular off-switch for synaptic transmission. AChE is present as clusters of three tetramers suspended by collagen stalks bound to the muscle membrane at varying density (600 µm-2 to 2500 µm-2) throughout the postsynaptic folds. As an extremely fast enzyme capable of destroying ACh molecules at rates approaching theoretical limits, AChE provides a very efficient mechanism to terminate synaptic transmission for subsequent signaling (Taylor, 1996
; Shen et al., 2002
).
Experimental data are available on the ultrastructure and activity of NMJs of different muscle types to guide the initial development of mathematical models (Land et al., 1981
, 1984
; Kandel and Siegelbaum, 1991
; Anglister et al., 1995
; Miyazawa et al., 1999
; Brejc et al., 2001
). Ultrastructural differences between the NMJs in vertebrate fast (twitch or extensor digitorum longus) and slow (tonic or soleus) muscles have been observed (Ellisman et al., 1976
; Gisiger and Stephens, 1982
; Florendo et al., 1983
; Fahim et al., 1984
). Also, geometric and reactive deviations in mouse NMJs due to muscular dystrophy have been documented (Shalton and Wareham, 1980
; Ellisman, 1981
; Tremblay et al., 1988
; Gisiger and Stephens, 1988
). Local measurements of miniature endplate current (mEPC) are used to infer the functional implications of such structural differences on an NMJ's efficiency. With our simulation infrastructure, we have begun to recreate the effects of the various parameters of different muscle types in silica. It is evident that experimental data directly coordinated with simulations are needed to bring increasingly realistic models to maturity.
| MATHEMATICAL SETTING |
|---|
|
|
|---|
, with appropriate conditions on its boundary, 
. Formally, the concentration varies from its initial state, u(x, 0), according to
![]() | (1) |
![]() | (2) |
is the outward unit normal, and 
act

denotes the cumulative reactive surface of AChE with specific reactivity kact (Smart and McCammon, 1998
act, for which a linear reaction scheme yielding a Robin boundary condition is now assumed. Estimation of the specific reactivity, kact, is described in each example. Attempts to lift this linear assumption on AChE binding and use more involved boundary conditions will be the focus of future research.
Numerical solution and visualization
Several methods of lines to solve Eqs. 1 and 2 have been built in the context of FEtk (Holst, 2001
), a finite element package written in "clean object-oriented" C language. Details on the resulting linear systems at each timestep, using a backward Euler method of lines, may be found in the Appendix. At each timestep, the conjugate gradient method is used with termination chosen such that time-truncation error is less than 10-10. Given a surface triangulation of an NMJ boundary with realistic folds, a robust and flexible mesh generation package, NETGEN (Schöberl, 1997
), develops a structured volume mesh that has been interfaced with FEtk (Figs. 1 and 8). A structured mesh fine enough to sufficiently describe the constraining surface geometry is used in the examples that follow, and refining this mesh (i.e., halving edge length) has no appreciable effect. Similarly, the timestep has been chosen to adequately resolve a postsynaptic response curve; e.g., the range 10-1 µs to 101 µs was studied to ensure that the simulation is not overly sensitive to the chosen timestep, 1 µs. Computing 1000 timesteps for a system of
33,000 vertices (such as the rectilinear synapse model) takes less than 1 h on an Intel (Santa Clara, CA) Xeon 1.80 GHz personal computer. The classical duo of backward Euler method of lines combined with conjugate gradient is nearly optimal for ordinary diffusion. Moreover, FEtk memory usage scales linearly with the number of vertices. The FEtk output formats are readable using MATLAB (MathWorks, Natick, MA), GMV (Los Alamos National Laboratory, Los Alamos, NM), and OpenDX (International Business Machines, White Plains, NY).
|
|
By defining the postsynaptic detection level, L(t), as the weighted surface integration
![]() | (3) |
R(x) is AChR density, we are able to observe general timescale trends of observed mEPC, such as rise time and signal duration. Toward more accurate estimations of mEPC, we account for ACh-AChR binding in the model to properly estimate the time course of open AChR ion channels. This entails appealing to the equilibrium mechanism
and assuming that AChR is saturated with ACh. For the equilibrium constant, we have
where
is the fraction of open AChR channels. This implies
leading to a new measure of postsynaptic response
![]() | (4) |
(t) does not consider the effect of ACh retention when bound to AChR. Such will require altering the zero-flux condition on postsynaptic folds with an integral to record accumulations of ACh bound to the ion channels. For brevity, we do not present plots of
(t) in all the examples that follow, as they cannot yet be compared with experimental data. | NMJ MODELS |
|---|
|
|
|---|
NMJ geometry and initial conditions
In Smart and McCammon (1998)
a simplified rectilinear model of an NMJ is considered with three identical secondary clefts and the following measurements (Fig. 1): primary cleft length = 2.0 µm; primary cleft width = secondary cleft width = 0.05 µm; separation between secondary clefts = 0.5 µm; secondary cleft depth = 0.8 µm.
We include a vesicle as a sphere of radius 0.024 µm fused to the middle of the presynaptic membrane. Its center is placed 0.016 µm above this membrane, leaving a circular area of radius 0.018 µm as the pore opening. The initial ACh concentration is only non-zero inside the vesicles, where it is 300 mM = 1.8 x 108 µm-3 (Schwartz, 1991
; Zimmermann, 1993
; Smart and McCammon, 1998
; Südhof and Scheller, 2000
). With the vesicle size stated above, this translates to 6.06 x 103 ACh molecules present. We recognize that this is a rather important yet simplistic aspect of our simulation, and sensitivity to the initial distribution of ACh will be included in future studies.
AChE clusters are treated as three-dimensional "holes" in the NMJ mesh with boundary condition Eq. 2. The previous finite element simulation (Smart and McCammon, 1998
) uses an AChE cluster density of
240 µm-2. By comparison, in Stiles and Bartol (2000)
, an active site density of 1800 µm-2 in the primary cleft (equivalent to a cluster density of
150 µm-2) and twice that in the secondary cleft is used. The length of the collagen stalk that connects AChE from the synaptic basal lamina is on the scale of 0.05 µm, depending on the species (Rotundo et al., 1997
). As done previously, we place AChE midway between the presynaptic and postsynaptic membranes in our rectilinear models. We have chosen to use cubic boxes with 0.02 µm sides spaced 0.1 µm apart in a square array spanning the primary and secondary clefts to afford a cluster density of 100 µm-2.
Reactivity of AChE
The reactivity of an AChE cluster in our model is relative to the surface area used to represent it. Following Smart and McCammon (1998)
, we determine
by conducting several separate but simple simulations. We represent an AChE cluster as a single cube with 0.02 µm edges, and place it at the center of a cube with 0.30 µm edges. The boundary condition for 
act in Eq. 2 is used on the surface of the interior cube, the AChE cluster. On the exterior cube boundary 
- 
act, the concentration of ACh is held at a different constant value for each "trial" simulation. The "current" of ACh consumed by AChE is defined to be
![]() | (5) |
E + P, the flux into the domain is estimated in linear proportion to the concentration of ACh maintained on the boundary (Rice, 1985
![]() | (6) |
3 x 109 M-1 · min-1 = 5 x 10-6 µm3 · µs-1 has been determined by experiment (Radi
et al., 1995
. FEtk with backward Euler timestepper is used to compute steady-state ACh diffusion between the concentric cubes, with validation from an explicit steady-state solver. Sampling the range 1.00 x 105 µm-3 to 1.00 x 107 µm-3 (that is, 0.166 mM to 16.6 mM) as the imposed ACh concentration on the outer boundary, we find
. This is slightly larger than the value used previously (Smart and McCammon, 1998
Solution and postsynaptic detection
The time course of ACh diffusion in the rectilinear NMJ model has been solved using backward Euler and FEtk (as outlined in the Appendix) with timesteps of 1 µs. Fig. 2 shows the total number of ACh molecules over 400 µs; the same decay has been shown in the solid curve in Fig. 3 of Smart and McCammon (1998)
.
|
|
R (Stiles and Bartol, 2000
R ranges from 7000 µm-2 to 10,000 µm-2. It is 2000 µm-2 to 3000 µm-2 another 0.2 µm to 0.3 µm from the crests, below which density falls off dramatically. Here we use a piecewise constant density of
R(x) = 8500 µm-2 from the crests down 0.25 µm into the secondary clefts,
R(x) = 2500 µm-2 an additional 0.25 µm into the secondary clefts, and then zero down to the troughs. The resulting postsynaptic detection level L(t) is shown in Fig. 3. The rise time (duration from the release of the vesicle to the peak) is
20 µs, much shorter than that of mEPC. However, since we do not yet consider ACh binding of AChR, this is reasonable. Moreover, it is comparable to the previous observation (Smart and McCammon, 1998
Variations in muscle type
Fast-twitch and slow-twitch NMJs
Simple rectilinear NMJ models for fast-twitch and slow-twitch muscle have been built according to the primitive dimensions listed in Table 1. Each model contains one vesicle, three identical secondary clefts, and the AChR distribution stated above. Also, as before, AChE clusters are 0.02 µm edged cubes placed 0.025 µm above the postsynaptic membrane with density 100 µm-2.
|
|
(t) is shown for the slow- and fast-twitch NMJ models. The peak-amplitudes are delayed and decay is accelerated in comparison to L(t). As expected, the more rapid decay of fast-twitch muscle signal is preserved. There are miniature endplate potential (mEPP as opposed to mEPC) data on fast- and slow-twitch muscle in mouse (Florendo et al., 1983
|
|
|
Mesh geometry and initial condition
We are grateful to have received a realistic three-dimensional surface-mesh model of an NMJ from Stiles and co-workers, formed from morphing a two-dimensional electron micrograph (Stiles and Bartol, 2000
; Stiles et al., 2001
). The resulting NETGEN volume mesh (Fig. 8) includes the vesicle, while realizing AChE clusters requires the special attention described below. In total, there are 175,609 vertices and 771,408 simplices in the volume mesh. Edges range from 2.45 x 10-3 µm to 1.05 x 10-1 µm in length, simplex volume ranges from 1.87 x 10-7 µm3 to 2.70 x 10-4 µm3, and the worst (largest) edge ratio in any simplex is 14.4.
The vesicle is a triangulated sphere of radius 0.024 µm placed roughly in the center of the presynaptic membrane. The pore opening is a single triangle, with area 0.00054 µm2, at the locus of vesicle attachment. Initially, the ACh concentration is zero everywhere but inside the vesicle, where it is 300 mM = 1.8 x 108 µm-3. This translates to
6.9 x 103 ACh molecules present at the beginning of simulation.
AChE clusters and reactivity
Placing AChE holes into the tetrahedral mesh is a nontrivial process that needs streamlining before high-throughput studies can be conducted. Currently, using the APBS (Baker et al., 2001
) and FEtk infrastructure, we repeatedly refine, down to a surface area threshold of 0.0024 µm2, the simplices that contain the positions of AChE monomers specified by Stiles and co-workers. The tetrahedral simplices containing the AChE positions are removed by marking adjacent faces of neighboring simplices as boundary. The resultant mesh has 55,589 tetrahedral AChE clusters, distributed as shown in Fig. 9. The average surface area is 0.0017 µm2, representative of 59,072 AChE molecules. This average cluster surface area must be accounted for when estimating the specific reactivity
Since Eq. 6 scales linearly with surface area of the cluster, the
estimated in the previous section can be adjusted according to the smaller average surface area in the realistic model (by a factor of
1.4). This simplification leads to
= 2.753 x 10-3 µm · µs-1.
|
|
| DISCUSSION AND CONCLUSIONS |
|---|
|
|
|---|
This research will continue to evolve on two fronts. First, with improved understanding of reaction kinetics, simulations will be more directly related to observable mEPC and mEPP. Second, increasing availability of ultrastructural data from electron tomography will offer the capacity for extensive comparative studies (Harlow et al., 2001
; Gustafsson et al., 2002
). Only in this manner can parameter space, including NMJ geometry, vesicle placement, receptor binding, and reactivity, be comprehensively explored to better understand neuromuscular diseases that affect NMJ processes.
The current study illustrates some of the effects that are associated with differences in the geometry of the NMJ, including the lengthening of the lifetime of ACh in the NMJ in the case of decreased number and depth in secondary folds. It sets the stage for more elaborate and realistic models that will be explored in the future.
| APPENDIX: DISCRETE FORMULATION |
|---|
|
|
|---|
![]() | (7) |
![]() | (8) |
u(x, tn) and define
![]() | (9) |
![]() | (10) |

, the weak form becomes
![]() | (11) |
1,...,
N}
V. The Galerkin approximation
![]() | (12) |
![]() | (13) |
To formulate Eq. 13 into a matrix equation, we write the terms
![]() | (14) |
![]() | (15) |
![]() | (16) |
![]() | (17) |
the mass matrix
and the solution vectors u = [ui], u° = [u°i]. | ACKNOWLEDGEMENTS |
|---|
|
|
|---|
K.T. is a La Jolla Interfaces in Science Fellow, supported partly by the Burroughs Wellcome Fund. This project is supported in part by the Howard Hughes Medical Institute, the W. M. Keck Foundation, the San Diego Supercomputer Center, the National Biomedical Computation Resource, the National Science Foundation, and the National Institutes of Health.
Submitted on August 19, 2002; accepted for publication December 12, 2002.
| REFERENCES |
|---|
|
|
|---|
Baker, N. A., D. Sept, S. Joseph, M. J. Holst, and J. A. McCammon. 2001. Electrostatics of nanosystems: application to microtubules and the ribosome. Proc. Natl. Acad. Sci. USA. 98:1003710041.
Braess, D. 1997. Finite Elements. Cambridge University Press, Cambridge.
Brejc, K., W. J. van Dijk, R. V. Klaassen, M. Schuurmans, J. van der Oost, A. B. Smit, and T. K. Sixma. 2001. Crystal structure of an ACh-binding protein reveals the ligand-binding domain of nicotinic receptors. Nature. 411:269276.[Medline]
Ellisman, M. H. 1981. The membrane morphology of the neuromuscular junction, sarcolemma, sarcoplasmic reticulum and transverse tubule system in murine muscular dystrophy studied by freeze-fracture electron microscopy. Brain Res. 214:261273.[Medline]
Ellisman, M. H., J. E. Rash, L. A. Staehelin, and K. R. Porter. 1976. Studies of excitable membranes. II. A comparison of specializations at neuromuscular junctions and nonjunctional sarcolemmas of mammalian fast and slow twitch muscle fibers. J. Cell Biol. 68:752774.
Fahim, M. A., J. A. Holley, and N. Robbins. 1984. Topographic comparison of neuromuscular junctions in mouse slow and fast twitch muscles. Neuroscience. 13:227235.[Medline]
Florendo, J. A., J. F. Reger, and P. K. Law. 1983. Electrophysiologic differences between mouse extensor digitorum longus and soleus. Exp. Neurol. 82:404412.[Medline]
Ghaffari-Farazi, T., J.-S. Liaw, and T. W. Berger. 1999. Consequence of morphological alterations on synaptic function. Neurocomputing. 2627:1727.
Gisiger, V., and H. Stephens. 19821983. Correlation between the acetylcholinesterase content in motor nerves and their muscles. J. Physiol. (Paris). 78:720728.
Gisiger, V., and H. R. Stephens. 1988. Localization of the pool of G4 acetylcholinesterase characterizing fast muscles and its alteration in murine muscular dystrophy. J. Neurosci. Res. 19:6278.[Medline]
Gustafsson, J. S., A. Birinyi, J. Crum, M. Ellisman, L. Brodin, and O. Shupliakov. 2002. Ultrastructural organization of lamprey reticulospinal synapses in three dimensions. J. Comp. Neurol. 450:167182.[Medline]
Harlow, M. L., D. Ress, A. Stoschek, R. M. Marshall, and U. J. McMahan. 2001. The architecture of active zone material at the frog's neuromuscular junction. Nature. 409:479484.[Medline]
Holst, M. 2001. Adaptive numerical treatment of elliptic systems on manifolds. Adv. Comput. Math. 15:139191.
Hosaka, Y., T. Yokota, Y. Miyagoe-Suzuki, K. Yuasa, M. Imamura, R. Matsuda, T. Ikemoto, S. Kameya, and S. Takeda. 2002.
1-Syntrophindeficient skeletal muscle exhibits hypertrophy and aberrant formation of neuromuscular junctions during regeneration. J. Cell Biol. 158:10971107.
Kandel, E. R., and S. A. Siegelbaum. 1991. Directly gated transmission at the nerve-muscle synapse. In Principles of Neural Science. E. R. Kandel, J. H. Schwartz, and T. M. Jessell, editors. Appleton and Lange, Norwalk, Connecticut. pp. 135152.
Land, B. R., E. E. Salpeter, and M. M. Salpeter. 1981. Kinetic parameters for acetylcholine interaction in intact neuromuscular junction. Proc. Natl. Acad. Sci. USA. 78:72007204.
Land, B. R., W. V. Harris, E. E. Salpeter, and M. M. Salpeter. 1984. Diffusion and binding constants for acetylcholine derived from the falling phase of miniature endplate currents. Proc. Natl. Acad. Sci. USA. 81:15941598.
Miyazawa, A., Y. Fujiyoshi, M. Stowell, and N. Unwin. 1999. Nicotinic acetylcholine receptor at 4.6 Å resolution: transverse tunnels in the channel wall. J. Mol. Biol. 288:765786.[Medline]
Naka, T., and N. Sakamoto. 2002. Simulation analysis of the effects of the simultaneous release of quanta of acetylcholine on the endplate current at the neuromuscular junction. Math. Comput. Sim. 59:8794.
Radi
, Z., D. M. Quinn, D. C. Vellom, S. Camp, and P. Taylor. 1995. Allosteric control of acetylcholinesterase catalysis by fasciculin. J. Biol. Chem. 270:2039120399.
Rice, S. A. 1985. Diffusion-Limited Reactions. Elsevier, Amsterdam, The Netherlands.
Rotundo, R. L., S. G. Rossi, and L. Anglister. 1997. Transplantation of quail collagen-tailed acetylcholinesterase molecules onto the frog neuromuscular synapse. J. Cell Biol. 136:367374.
Schöberl, J. 1997. NETGEN: an advancing front 2D/3D-mesh generator based on abstract rules. Comput. Visual Sci. 1:4152.
Schwartz, J. H. 1991. Synaptic vesicles. In Principles of Neural Science. E. R. Kandel, J. H. Schwartz, and T. M. Jessell, editors. Appleton and Lange, Norwalk, Connecticut. pp. 225234.
Shalton, P. M., and A. C. Wareham. 1980. Some factors affecting spontaneous transmitter release in dystrophic mice. Muscle Nerve. 3:120127.[Medline]
Shen, T., K. Tai, R. Henchman, and J. A. McCammon. 2002. Molecular dynamics of acetylcholinesterase. Acc. Chem. Res. 35:332340.[Medline]
Smart, J. L., and J. A. McCammon. 1998. Analysis of synaptic transmission in the neuromuscular junction using a continuum finite element model. Biophys. J. 75:16791688.
Stiles, J. R., and T. M. Bartol. 2000. Monte Carlo methods for simulating realistic synaptic microphysiology using MCell. In Computational Neuroscience: Realistic Modeling for Experimentalists. E. De Schutter, editor. CRC Press, New York. pp. 87127.
Stiles, J. R., T. M. Bartol, M. M. Salpeter, E. E. Salpeter, and T. J. Sejnowski. 2001. Synaptic variability: new insights from reconstructions and Monte Carlo simulations with MCell. In Synapses. W. M. Cowan, T. C. Südhof, and C. F. Stevens, editors. Johns Hopkins University Press, Baltimore, Maryland. pp. 681731.
Südhof, T. C., and R. H. Scheller. 2000. Mechanism and regulation of neurotransmitter release. In Synapses. W. M. Cowan, T. C. Südhof, and C. F. Stevens, editors. Johns Hopkins University Press, Baltimore, Maryland. pp. 177215.
Taylor, P. 1996. Anticholinesterase agents. In The Pharmacological Basis of Therapeutics. J. G. Hardman, L. E. Limbird, P. B. Molinoff, R. W. Ruddon, and A. G. Gilman, editors. McGraw-Hill, New York. pp. 161176.
Tremblay, J. P., L. Grégoire, R. Sasseville, G. Guay, and C. Belhumeur. 1988. Reduction of postjunctional fold density and depth in dystrophic mice. Synapse. 2:148156.[Medline]
Zimmermann, H. 1993. Synaptic Transmission: Cellular and Molecular Basis. Thieme, Stuttgart, Germany.
This article has been cited by other articles:
![]() |
C. P. Bolter and D. J. English The effects of tertiapin-Q on responses of the sinoatrial pacemaker of the guinea-pig heart to vagal nerve stimulation and muscarinic agonists Exp Physiol, January 1, 2008; 93(1): 53 - 63. [Abstract] [Full Text] [PDF] |
||||
![]() |
Y. Cheng, J. K. Suen, D. Zhang, S. D. Bond, Y. Zhang, Y. Song, N. A. Baker, C. L. Bajaj, M. J. Holst, and J. A. McCammon Finite Element Analysis of the Time-Dependent Smoluchowski Equation for Acetylcholinesterase Reaction Rate Calculations Biophys. J., May 15, 2007; 92(10): 3397 - 3406. [Abstract] [Full Text] [PDF] |
||||
![]() |
Y. Song, Y. Zhang, C. L. Bajaj, and N. A. Baker Continuum Diffusion Reaction Rate Calculations of Wild-Type and Mutant Mouse Acetylcholinesterase: Adaptive Finite Element Analysis Biophys. J., September 1, 2004; 87(3): 1558 - 1566. [Abstract] [Full Text] [PDF] |
||||
![]() |
Y. Song, Y. Zhang, T. Shen, C. L. Bajaj, J. A. McCammon, and N. A. Baker Finite Element Solution of the Steady-State Smoluchowski Equation for Rate Constant Calculations Biophys. J., April 1, 2004; 86(4): 2017 - 2029. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |