help button home button Biophys. J.
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS

This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Tai, K.
Right arrow Articles by McCammon, J. A.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Tai, K.
Right arrow Articles by McCammon, J. A.
Biophysical Journal 84:2234-2241 (2003)
© 2003 The Biophysical Society

Finite Element Simulations of Acetylcholine Diffusion in Neuromuscular Junctions

Kaihsu Tai*, Stephen D. Bond*,{dagger}, Hugh R. MacMillan*,{dagger}, Nathan Andrew Baker{ddagger}, Michael Jay Holst{dagger} and J. Andrew McCammon*,§

* Department of Chemistry and Biochemistry, {dagger} Department of Mathematics, § Department of Pharmacology, and the Howard Hughes Medical Institute, University of California, San Diego, La Jolla, California 92093; and {ddagger} 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
 TOP
 ABSTRACT
 INTRODUCTION
 The neuromuscular junction
 MATHEMATICAL SETTING
 NMJ MODELS
 DISCUSSION AND CONCLUSIONS
 APPENDIX: DISCRETE FORMULATION
 ACKNOWLEDGEMENTS
 REFERENCES
 
A robust infrastructure for solving time-dependent diffusion using the finite element package FEtk has been developed to simulate synaptic transmission in a neuromuscular junction with realistic postsynaptic folds. Simplified rectilinear synapse models serve as benchmarks in initial numerical studies of how variations in geometry and kinetics relate to endplate currents associated with fast-twitch, slow-twitch, and dystrophic muscles. The flexibility and scalability of FEtk affords increasingly realistic and complex models that can be formed in concert with expanding experimental understanding from electron microscopy. Ultimately, such models may provide useful insight on the functional implications of controlled changes in processes, suggesting therapies for neuromuscular diseases.


    INTRODUCTION
 TOP
 ABSTRACT
 INTRODUCTION
 The neuromuscular junction
 MATHEMATICAL SETTING
 NMJ MODELS
 DISCUSSION AND CONCLUSIONS
 APPENDIX: DISCRETE FORMULATION
 ACKNOWLEDGEMENTS
 REFERENCES
 
The neuromuscular junction (NMJ) is the point of communication between neurons and muscle fiber in the orchestration of muscle contraction. This study encompasses the release of neurotransmitter acetylcholine (ACh), its hydrolysis with acetylcholinesterase (AChE) clusters, and its reactive presence near acetylcholine receptor (AChR) molecules. Past computational modeling of synaptic transmission on this scale has involved either differential equations governing continuum reaction-diffusion (Smart and McCammon, 1998Go; Ghaffari-Farazi et al., 1999Go) or particle methods such as Brownian dynamics and Monte Carlo (Stiles and Bartol, 2000Go; Stiles et al., 2001Go). Here, we present an improved finite element framework to solve continuum reaction-diffusion using FEtk (Holst, 2001Go), an efficient platform for adaptive multiscale modeling. The importance of this framework is its flexibility to evolve alongside improved understandings of reaction kinetics. Now, with the capacity to represent realistic NMJs, comparative studies can be conducted with the guidance of coordinated experimental data. This should provide insight on the functional implications of NMJ variations associated with neuromuscular diseases, such as muscular dystrophy and myasthenia gravis, ultimately suggesting potential therapeutic intervention.


    The neuromuscular junction
 TOP
 ABSTRACT
 INTRODUCTION
 The neuromuscular junction
 MATHEMATICAL SETTING
 NMJ MODELS
 DISCUSSION AND CONCLUSIONS
 APPENDIX: DISCRETE FORMULATION
 ACKNOWLEDGEMENTS
 REFERENCES
 
A typical NMJ features a smooth presynaptic neuron membrane and a folded postsynaptic muscle surface of crests and troughs. When an action potential reaches the end of the nerve, it results in a localized influx of calcium ions that, in turn, causes vesicles to fuse to the neuron membrane and release ACh. In our model, vesicles are treated as having just opened and we do not yet account for subsequent relaxation of the neuron membrane. Eventually, this effect may be incorporated into FEtk, along with spatial control of vesicle placement according to the varying presence of calcium ions.

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, 1996Go; Shen et al., 2002Go).

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., 1981Go, 1984Go; Kandel and Siegelbaum, 1991Go; Anglister et al., 1995Go; Miyazawa et al., 1999Go; Brejc et al., 2001Go). 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., 1976Go; Gisiger and Stephens, 1982Go; Florendo et al., 1983Go; Fahim et al., 1984Go). Also, geometric and reactive deviations in mouse NMJs due to muscular dystrophy have been documented (Shalton and Wareham, 1980Go; Ellisman, 1981Go; Tremblay et al., 1988Go; Gisiger and Stephens, 1988Go). 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
 TOP
 ABSTRACT
 INTRODUCTION
 The neuromuscular junction
 MATHEMATICAL SETTING
 NMJ MODELS
 DISCUSSION AND CONCLUSIONS
 APPENDIX: DISCRETE FORMULATION
 ACKNOWLEDGEMENTS
 REFERENCES
 
Continuum formulation
In a continuum model, the concentration of ACh, u(x, t), is assumed to satisfy the time-dependent diffusion equation in the synaptic cleft, {Omega}, with appropriate conditions on its boundary, {partial}{Omega}. Formally, the concentration varies from its initial state, u(x, 0), according to

(1)

(2)
where D = 4.0 x 10-4 µm2 · µs-1 is the diffusion coefficient, is the outward unit normal, and {partial}{Omega}act {partial}{Omega} denotes the cumulative reactive surface of AChE with specific reactivity kact (Smart and McCammon, 1998Go). Note that Eq. 2 states a zero-flux condition on all boundaries except the surfaces representing AChE clusters, {partial}{Omega}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, 2001Go), 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, 1997Go), 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).



View larger version (61K):
[in this window]
[in a new window]
 
FIGURE 1  Three views (edges, outside, inside) of the finite element mesh for the rectilinear model of the neuromuscular junction, with three secondary clefts and a spherical vesicle fused to the presynaptic membrane. The cubes represent acetylcholinesterase. 1°, primary cleft; 2°, secondary cleft.

 


View larger version (77K):
[in this window]
[in a new window]
 
FIGURE 8  Three views of the realistic mesh: overview; enlarged box containing the vesicle on the presynaptic membrane; and a cross section of the secondary cleft (shaded). The AChE "holes" are the patches of white within. The height of the clefts from crest to trough is ~1 µm.

 
Postsynaptic detection and AChR binding
Eqs. 1 and 2 quantify the presence of ACh in the NMJ, but relating this to experimentally observable mEPC involves estimating not only the density of AChR, but also the stochastic nature of its binding and retention. We now present two approximations for detection, and stress that improved realism will be the subject of future work.

By defining the postsynaptic detection level, L(t), as the weighted surface integration

(3)
where {gamma}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 {theta} is the fraction of open AChR channels. This implies leading to a new measure of postsynaptic response

(4)
According to Naka and Sakamoto (2002)Go, K = 3.6 x 102 mM-2 = 1.0 x 10-9 µm6. Note that {Lambda}(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 {Lambda}(t) in all the examples that follow, as they cannot yet be compared with experimental data.


    NMJ MODELS
 TOP
 ABSTRACT
 INTRODUCTION
 The neuromuscular junction
 MATHEMATICAL SETTING
 NMJ MODELS
 DISCUSSION AND CONCLUSIONS
 APPENDIX: DISCRETE FORMULATION
 ACKNOWLEDGEMENTS
 REFERENCES
 
Rectilinear benchmarking
For our first model, we have created a simplified rectilinear NMJ similar to that in Smart and McCammon (1998)Go. We discuss the details of implementing this model before presenting preliminary studies on the functional implications of variations in geometry and kinetics. Rectilinear models have also been studied with MCell (Stiles and Bartol, 2000Go; Stiles et al., 2001Go). All simulations are in micrometers and microseconds, with concentrations in particles per volume and densities in particles per area.

NMJ geometry and initial conditions
In Smart and McCammon (1998)Go 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, 1991Go; Zimmermann, 1993Go; Smart and McCammon, 1998Go; Südhof and Scheller, 2000Go). 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, 1998Go) uses an AChE cluster density of ~240 µm-2. By comparison, in Stiles and Bartol (2000)Go, 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., 1997Go). 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)Go, 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 {partial}{Omega}act in Eq. 2 is used on the surface of the interior cube, the AChE cluster. On the exterior cube boundary {partial}{Omega} - {partial}{Omega}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)
Naturally, J must also equal the flux of ACh into the domain due to imposing constant ACh concentration on the outer boundary. Assuming a nonsaturated steady state and the catalytic scheme E + S -> E + P, the flux into the domain is estimated in linear proportion to the concentration of ACh maintained on the boundary (Rice, 1985Go; Smart and McCammon, 1998Go). Writing this proportionality constant as k, we set

(6)
and consider simulations with different sustained steady-state levels of ACh on the outer boundary. Since k = kcat/KM {approx} 3 x 109 M-1 · min-1 = 5 x 10-6 µm3 · µs-1 has been determined by experiment (Radic et al., 1995Go), enforcing Eq. 6 for each steady state leads to an estimation of . 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, 1998Go), due to the fact that we are using here a different area for the active site of AChE. With the enzymatic activity of AChE approximated, Eqs. 1 and 2 can now be solved for the rectilinear NMJ.

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)Go.



View larger version (7K):
[in this window]
[in a new window]
 
FIGURE 2  Total number of ACh molecules in the rectilinear NMJ model decreases over time.

 


View larger version (7K):
[in this window]
[in a new window]
 
FIGURE 3  The postsynaptic detection level L(t) in the rectilinear NMJ model.

 
Determining the postsynaptic detection level in Eq. 3 entails defining AChR density. MCell models typically use a graduated receptor density function {gamma}R (Stiles and Bartol, 2000Go; Stiles et al., 2001Go). From the crest to 0.2 µm to 0.3 µm into the fold, {gamma}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 {gamma}R(x) = 8500 µm-2 from the crests down 0.25 µm into the secondary clefts, {gamma}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, 1998Go).

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.


View this table:
[in this window]
[in a new window]
 
TABLE 1 
 
In Fig. 4, the rise time of L(t) for the fast-twitch NMJ model is 45 µs, and a peak value of 4.25 x 108 µm3 is obtained. In contrast, the slow-twitch geometry leads to a more gradual response, reaching the peak value 3.12 x 108 µm3 in 43 µs. The decay time (duration from vesicle release until the decay reaches half of the peak L(t)) for the fast-twitch muscle is 213 µs, significantly shorter than that displayed in the slow-twitch geometry, 270 µs.



View larger version (10K):
[in this window]
[in a new window]
 
FIGURE 4  L(t) for NMJ models of fast-twitch and slow-twitch muscles.

 
In Fig. 5, {Lambda}(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., 1983Go). Signal decay according to this data is longer, possibly due to oversimplification in the rectilinear model.



View larger version (11K):
[in this window]
[in a new window]
 
FIGURE 5  {Lambda}(t) for NMJ models of fast-twitch and slow-twitch muscles.

 
Dystrophic NMJs
Some dystrophic muscle displays improperly developed NMJs in terms of ultrastructure (Hosaka et al., 2002Go), but again specific data are limited. For the rectilinear dystrophic NMJ model described in Table 1, the rise time, decay time, and amplitude of L(t) are all slightly lower than normal (Fig. 6). Of course, this simply states the obvious: such muscle does not function as efficiently! It is likely that NMJ malformation is less a cause than an effect in the pathology of dystrophy.



View larger version (11K):
[in this window]
[in a new window]
 
FIGURE 6  Normal versus dystrophic fast-twitch muscle.

 
Reduced AChE density
As could be said for the distribution of AChR, aberrations in the presence of AChE can have functional implications. For the slow-twitch muscle NMJ, lowering AChE density to 50 µm-2 leads to a slightly longer rise time, 46 µs, with a lower peak level of 2.81 µm3 (Fig. 7). The decay time is about twice as long, at 505 µs. Specific experimental information on variations in AChE density among different muscle types is limited. As such, this example merely serves to illustrate that an abnormal density can have an appreciable effect on the efficiency of synaptic transmission across an NMJ.



View larger version (10K):
[in this window]
[in a new window]
 
FIGURE 7  Comparing L(t) for different AChE densities. High density, 100 µm-2; low density, 50 µm-2.

 
Realistic NMJ model
Finally, we demonstrate the capacity of FEtk to model ACh transmission within a realistic NMJ geometry drawn from electron microscopy. Given this realism, our infrastructure is poised to accompany experiments.

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, 2000Go; Stiles et al., 2001Go). 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., 2001Go) 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.



View larger version (8K):
[in this window]
[in a new window]
 
FIGURE 9  Surface area distribution of the AChE holes in the realistic mesh.

 
AChR density and postsynaptic detection level
We use approximately the same AChR density criterion as the above rectilinear models. Fig. 10 shows the postsynaptic detection level L(t) for the more realistic NMJ model. The rise time is 10 µs, about half as long as that in the previous rectilinear NMJ model, but comparison at this stage is somewhat unfounded.



View larger version (8K):
[in this window]
[in a new window]
 
FIGURE 10  The postsynaptic detection level L(t) over 150 µs in the realistic NMJ model.

 

    DISCUSSION AND CONCLUSIONS
 TOP
 ABSTRACT
 INTRODUCTION
 The neuromuscular junction
 MATHEMATICAL SETTING
 NMJ MODELS
 DISCUSSION AND CONCLUSIONS
 APPENDIX: DISCRETE FORMULATION
 ACKNOWLEDGEMENTS
 REFERENCES
 
We have shown that simple rectilinear models are capable of capturing electrophysiological trends observed in NMJs from different muscle types. Difference in decay time between slow- and fast-twitch muscles is observed and dystrophic muscle exhibits a muted postsynaptic response. Our effort to make physiologically relevant connections is an advancement from previous continuum models (Smart and McCammon, 1998Go). As our primary contribution, we have created a finite element infrastructure capable of dealing with a realistic mesh based on electron microscopy data. This motivates future partnerships with coordinated experimental investigations on the relationship between muscle function and NMJ architecture.

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., 2001Go; Gustafsson et al., 2002Go). 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
 TOP
 ABSTRACT
 INTRODUCTION
 The neuromuscular junction
 MATHEMATICAL SETTING
 NMJ MODELS
 DISCUSSION AND CONCLUSIONS
 APPENDIX: DISCRETE FORMULATION
 ACKNOWLEDGEMENTS
 REFERENCES
 
We use the backward Euler algorithm to define a method of lines and reduce Eqs. 1 and 2 to the following elliptic problem at time tn, knowing the previous concentration, u(x, tn-1):

(7)

(8)
It is convenient to write un colone u(x, tn) and define

(9)
Integrating by parts, we obtain the weak form of Eq. 8,

(10)
where V is the test space (Braess, 1997Go; Holst, 2001Go). Enforcing boundary condition on {partial}{Omega}, the weak form becomes

(11)
For a discrete solution to Eq. 11, we employ a finite element space Vh = span{{phi}1,...,{phi}N} V. The Galerkin approximation

(12)
satisfies the discrete weak form

(13)
knowing the discrete solution from the previous timestep, To formulate Eq. 13 into a matrix equation, we write the terms

(14)

(15)

(16)
It follows that

(17)
where the stiffness matrix the mass matrix and the solution vectors u = [ui], u° = [u°i].


    ACKNOWLEDGEMENTS
 TOP
 ABSTRACT
 INTRODUCTION
 The neuromuscular junction
 MATHEMATICAL SETTING
 NMJ MODELS
 DISCUSSION AND CONCLUSIONS
 APPENDIX: DISCRETE FORMULATION
 ACKNOWLEDGEMENTS
 REFERENCES
 
Our thanks go to Dr. Burak Aksoylu, Prof. Mark H. Ellisman, Dr. Jason Smart, Dr. Jason Ka-Chun Suen, and Prof. Palmer W. Taylor for useful discussions. We are grateful to Prof. Tom Bartol and Prof. Joel Stiles for providing the realistic NMJ surface mesh.

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
 TOP
 ABSTRACT
 INTRODUCTION
 The neuromuscular junction
 MATHEMATICAL SETTING
 NMJ MODELS
 DISCUSSION AND CONCLUSIONS
 APPENDIX: DISCRETE FORMULATION
 ACKNOWLEDGEMENTS
 REFERENCES
 
Anglister, L., J. R. Stiles, B. Haesaert, J. Eichler, and M. M. Salpeter. 1995. Acetylcholinesterase at neuromuscular junctions. In Enzymes of the Cholinesterase Family. D. M. Quinn, A. S. Balasubramanian, B. P. Doctor, and P. Taylor, editors. Plenum, New York. pp. 277–285.

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:10037–10041.[Abstract/Free Full Text]

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:269–276.[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:261–273.[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:752–774.[Abstract/Free Full Text]

Fahim, M. A., J. A. Holley, and N. Robbins. 1984. Topographic comparison of neuromuscular junctions in mouse slow and fast twitch muscles. Neuroscience. 13:227–235.[Medline]

Florendo, J. A., J. F. Reger, and P. K. Law. 1983. Electrophysiologic differences between mouse extensor digitorum longus and soleus. Exp. Neurol. 82:404–412.[Medline]

Ghaffari-Farazi, T., J.-S. Liaw, and T. W. Berger. 1999. Consequence of morphological alterations on synaptic function. Neurocomputing. 26–27:17–27.

Gisiger, V., and H. Stephens. 1982–1983. Correlation between the acetylcholinesterase content in motor nerves and their muscles. J. Physiol. (Paris). 78:720–728.

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:62–78.[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:167–182.[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:479–484.[Medline]

Holst, M. 2001. Adaptive numerical treatment of elliptic systems on manifolds. Adv. Comput. Math. 15:139–191.

Hosaka, Y., T. Yokota, Y. Miyagoe-Suzuki, K. Yuasa, M. Imamura, R. Matsuda, T. Ikemoto, S. Kameya, and S. Takeda. 2002. {alpha}1-Syntrophin–deficient skeletal muscle exhibits hypertrophy and aberrant formation of neuromuscular junctions during regeneration. J. Cell Biol. 158:1097–1107.[Abstract/Free Full Text]

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. 135–152.

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:7200–7204.[Abstract/Free Full Text]

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:1594–1598.[Abstract/Free Full Text]

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:765–786.[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:87–94.

Radic, Z., D. M. Quinn, D. C. Vellom, S. Camp, and P. Taylor. 1995. Allosteric control of acetylcholinesterase catalysis by fasciculin. J. Biol. Chem. 270:20391–20399.[Abstract/Free Full Text]

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:367–374.[Abstract/Free Full Text]

Schöberl, J. 1997. NETGEN: an advancing front 2D/3D-mesh generator based on abstract rules. Comput. Visual Sci. 1:41–52.

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. 225–234.

Shalton, P. M., and A. C. Wareham. 1980. Some factors affecting spontaneous transmitter release in dystrophic mice. Muscle Nerve. 3:120–127.[Medline]

Shen, T., K. Tai, R. Henchman, and J. A. McCammon. 2002. Molecular dynamics of acetylcholinesterase. Acc. Chem. Res. 35:332–340.[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:1679–1688.[Abstract/Free Full Text]

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. 87–127.

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. 681–731.

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. 177–215.

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. 161–176.

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:148–156.[Medline]

Zimmermann, H. 1993. Synaptic Transmission: Cellular and Molecular Basis. Thieme, Stuttgart, Germany.




This article has been cited by other articles:


Home page
Exp PhysiolHome page
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]


Home page
Biophys. JHome page
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]


Home page
Biophys. JHome page
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]


Home page
Biophys. JHome page
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]


This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Tai, K.
Right arrow Articles by McCammon, J. A.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Tai, K.
Right arrow Articles by McCammon, J. A.


HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
Copyright © 2003 by the Biophysical Society.