T. Bartol; T. Sejnowski; B. Land; E. Salpeter ; M. Salpeter

(based on Neurosciences poster 422.17 presented at the Nov 2000 meeting)


The simplest chemical kinetics scheme for the development of miniature endplate currents (mepcs) at the neuromuscular junction (NMJ) is:

To explain mepc characteristics (e.g. rise-time, amplitude, fall-time) with a simulation one has to assign numerical values to the four kinetics parameters, plus a diffusion constant (D) and quantal packet size. Previous measurements on lizard neuromuscular junctions (NMJ's) under various experimental conditions have led to one set of MCell simulation parameters (http://www.mcell.cnl.salk.edu) that match observed mepc characteristics. A major question of interest now is uniqueness of these simulation parameters and their relevance to other synapses. We carried out 47000 MCell simulations of a simplified NMJ to derive a table of 1568 sets of simulation parameters and their predicted mepc characteristics. With a slight loss of accuracy, only four parameters: b/a, b/k_, k+/k_ and D/k_ needed to be varied, allowing a wide range to be covered. The simulations were performed at the San Diego Supercomputer Center on BlueHorizon, an IBM SP with 1152 processors. The table is summarized in graphical form. An algorithm is presented for searching and interpolating the table for determining the range of simulation parameters which will equally well match mepc characteristics and for comparisons to other synapses. Supported by: NIH NS-09315 NSF IBN-9603611, and by NSF to the San Diego Supercomputer Center.


Glossary nmj = neuromuscular junction
mepc = miniature endplate current

Physical input parameters:

Kp = forward binding constant (in 1/mol sec)
Km = unbinding constant (in 1/sec)
b = channel opening rate (in 1/sec)
a = channel closing rate (in 1/sec)
D = diffusion constant (in cm2/sec)
N = number of ACh molecules in a quantal packet (0.5N is maximum number of channel openings)
h = height of the synaptic cleft (in microns) (bounded by two infinite planes)
s = density of ACh receptors (micron-2)
so = normal value of s

Three Timescales:

Binding time: tb = 1.39*h / (s *Kp)
Diffusion time for a 'saturated disk' of area (N/s): td = 0.8*N / (4*p*s*D)
Mean burst duration for open channels: mbd = (1+b/2Km) / a

Three simulated conditions for mepcs (see also Fig 1 below):
(1) 'normal' : s = so and acetylcholinesterase (AChE) intact.
(AChE parameters from Stiles, et.al. Biophys J. 77, 1177, 1999)
(2) 'noAChE': esterases omitted.
(3) 'reduced': s = 0.3*so and esterases omitted.

Figure 1. Simulated mepcs under three conditions: (1) at normal AChR density, (2) with AChE removed (3) with AChE removed and AChR density of 1/3 normal. All three were computed using the same input parameters from the middle of the ranges described in the Glossary

Output quantities from the mepc simulations:

tr = rise time (20% to 80% of peak during rising phase)
tf = fall time (90% to 33% of peak during falling phase) plateau time (80% during rising to 90% of peak during falling phase)
P = plateau time / (tr tf ),0.5
hch = peak mepc amplitude, expressed as number of open channels efficiency = (hch /0.5N)x100 , expressed in percent (Measurement of rise time, fall time and plateau are also shown in Fig 2).

Figure 2. The measurement scheme for simulated mepcs. Averaged traces were used in all cases. The peak value was determined from a quadratic fit. Rise time, fall time, and plateau time were measured as shown between the 20-80% rise, 90-30% fall and 80% rise-90% fall respectively.

Experimental Target Values: Three dimensionless ratios predicted by Monte Carlo simulation results are:

(tr normal) / (tf normal)
(tf noAChE) / (tf normal)
(tr reduced) / (tr noAChE)

Experimental values of these ratios depend on the particular muscle and temperature, but typical values are: tf normal ~= 1200 msec
tr normal ~= 90 msec
(tr normal) / (tf normal) ~= 13.3
(tf noAChE) / (tf normal) ~= 3.63
(tr reduced) / (tr noAChE) ~= 1.90.

MCell Simulations

Monte Carlo simulations of the reaction scheme outlined in the introduction were performed using MCell, a general simulation system for coupled reaction-diffusion systems with arbitrary bounding geometry. For the MCell simulations we fixed four quantities: h=0.05 microns, so=7500 AChR/micron2, N=5000 and mbd=1000 msec. These values imply tb=279/ Kp and td=212/D msec. We carried out 1568 sets of simulations separately for the three conditions 'normal', 'noAChE' and 'reduced'. These simulations covered all combinations of:

b/a = 100, 50, 20, 10;
b/2Km = 8, 4, 2, 1, 0.5, 0.25, 0.125;
tb = 279, 139, 70, 35, 17, 9, 4;
td = 200, 141, 100, 71, 50, 35, 25, 18;

To minimize stochastic noise, each set of simulation parameters was run 20 times and the results averaged. A quadratic curve fit was done near the peak to obtain an estimate of the amplitude. The estimated rise and fall times (Fig 2) were then obtained by line fits to appropriate portions of the mepc waveform. Each simulation set yields the output values tr, tf, hch, and P defined in the Glossary for each of the three conditions. The first version of the catalog (see below) consists of the input parameters and the corresponding outputs. The values of hch vary greatly and the values of tf vary slightly over the whole simulation catalog.

Some of the results are shown in Fig 3-5 below.

Figure 3. Model predictions for (a) rise time, (b) efficiency, (c) fall time, plotted against binding time tb. The simulations were for normal conditions with diffusion time td of 35 (green), 71 (blue) and 140 (red) msecs. In all cases b/a=50 and b/2Km=1. The thin,black diagonal lines represent the effect of changing tb and td by the same ratio. (a) Note that rise time varies remarkably little when tb is varied considerably with td constant. Thus a pathological condition which increases cleft-height h only, such as myasthenia gravis, would have little effect on the rising phase, if s stayed the same. However, if td is increased by dropping s (with BTX or a pathological condition) then the rise time increases. (b) The efficiency decreases rapidly with increasing tb, when tb is large. A disease which increases h/s (such as myasthenia gravis) and thus decreases tb would have a less serious effect on mepc amplitude if Kp were large. (c) For medium-to-large tb, the fall time is close to the mean burst duration. Only for very small tb is the binding to receptors so efficient that (even in the presence of esterase) there is significant rebinding during the falling phase.


Figure 4. Model predictions for (a) rise time and (b) efficiency plotted against b/2Km. The plots are for values of the isomerization ratio b/a of 10 (blue), 20 (green) and 100 (red). Note again that the rise time varies very little with varying b/2Km, especially for the most likely intermediate values of b/a. The efficiency is somewhat small when b/2Km is very small, i.e. when unbinding inhibits channel openings.


Figure 5. Two model mepcs (normal condition) with different input parameters, but giving the same values for tr and tf. The blue curve comes from a model with large Kp and small D which gives a high efficiency of 54% while the red curve comes from a model with small Kp and large D with quite low efficiency of 6% (requiring a large N). The two mepcs are almost indistinguishable. The shape of a mepc under normal conditions alone cannot give information on whether the kinetics parameters lead to high or low efficiency. The ratio (tf noAChE) / (tf normal) is appreciably different for the two models, 51% for the blue mepc and 5.6% for the red mepc. The relationship between this ratio and the efficiency is shown in Figure 6.


Figure 6. (a) A scatter diagram of efficiency versus tb/td. Efficiency generally decreases as tb/td increases, but as the color indicates, three ranges of tb+td are mixed together. (b) A scatter diagram of (tf noAChE) / (tf normal) versus tb/td. (c) A scatter diagram of efficiency versus the target ratio(tf noAChE) / (tf normal). The correlation between the two quantities is fairly strong, since both decrease with increasing tb/td. The curvature is due to the fact that (tf noAChE) / (tf normal) is always greater than 1 and efficiency is always less than 100%. For all panels, the colors indicate values for tb+td of >150 (red), 75 to 150 (green), and <75 (black).


The Output Catalogs

Scaling Prescriptions

The six rate parameters plus the quantal packet size and diffusion constant of ACh are needed for a Monte Carlo simulation (see Glossary). There are two scaling arguments which greatly reduce the number of simulations required:

  1. Given one simulation with mbd=1000 msec, if a, b, Kp, Km, and D are multiplied by the same factor x, then for all three conditions, tr and tf are multiplied by x, and hch and P are unchanged.
  2. If N and D are multiplied by the same factor , then amplitude hch is multiplied by that factor but efficiency and all rise and fall times are unchanged.

The model output must be converted from normalized, model units back to physical units for comparison with real data. If we start with an MCell simulation which produces some values for hch and tf (normal condition), we can convert to real units by defining two scaling factors: The amplitude scaling factor is defined as Sa=(actual amplitude)/(model amplitude) and the time scaling factor is defined as St=(actual falltime)/(model falltime). Multiplying the modeled waveform amplitude by Sa and the modeled waveform time by St scales the model waveform to the actual data. Removing the amplitude scaling for N and D consists of multiplying the model values by Sa. Removing the time scaling for a, b, Kp, Km, and D consists of multiplying by 1/St.

Range of Models which Fit the Three Target Values

Another version of the simulation catalog consists of the 1568 models which have been scaled so that all have amplitude hch=872 and fall time tf =1200 msec. With this scaling, the full range of parameter values explored by the models is; a =880 to 61000 1/sec, b=9100 to 2.7x106 1/sec, Kp=0.62 to 430x107 1/mol*sec, and Km=4900 to 4.9 x105 1/sec. Of the 1568 models, 73 predict all three target ratios (see Glossary) within 20%. Table 1a gives the range of kinetic parameters for this subset. With four input parameters to be chosen, and only three output quantities to fit to target values, we could not obtain input values uniquely. One therefore has to specify at least one input parameter (from experiment). For instance, many people have measured b although there seems to be a fairly wide range in reported values. In principle, by choosing some exact value for b and requiring an exact fit to three target values, there should be a unique set of reaction parameters. In practice, one has to accept an uncertainty range on each target value. The catalog user can specify the four uncertainty ranges and then obtain the set of models which satisfy the ranges. One example is in Table 1b and 1c. Two different ranges of b were chosen from the literature. Out of the 73 best fitting models, 9 have b within the range (60,000 to 100,000) shown in Table 1b, and 9 have b within the range (40,000 to 60,000) shown in Table 1c. Table 1b and 1c give the range of a, Kp, and Km for each range of b. If one specifies the value of a second input parameter to within some uncertainty range, the set of models is further narrowed. For instance, some measurements of a have also been done. If one specifies a to be in the range 2900 to 5000 (with b in the range in Table 1c), then only one of the 1568 models remaining is b/a=20, b/2Km=2, tb=70, td=100. On the other hand if one chooses a in the range 1500 to 2500, the only remaining model is b/a=20, b/2Km=1, tb=35, td=71.

Another way of looking at the wide range of parameters which fit the data is to plot surfaces representing the locus of parameters which yield each of the three target ratios in input-parameter space. The following figure has axes proportional to the logs of D, Kp, and b/2Km. Since the green and blue surfaces (optimal
(tf noAChE) / (tf normal) and (tr reduced) / (tr noAChE) respectively) are almost conicident throughout the parameter range, we lose power to uniquely fit data.


  1. To specify a model uniquely eight parameters have to be specified. In our catalog only four dimensionless parameters are varied; b/a, b/2Km, tb and td, while four others are kept fixed. However scaling prescriptions are given to obtain results for other values, without requiring further simulations.
  2. Although one would expect the rise time to be directly correlated with the binding time, tr varies remarkably little over a wide range of tb, if td is held constant. This near-constancy is due to some accidental cancellations: As tb increases, the efficiency decreases, which results in the peak occurring earlier in the binding phase.
  3. With b and three target values specified to within 20%, the resultant range of uncertainty in a, Kp, and Km is enormously larger than 20%. To obtain a moderately unique model, one has to specify the value of a second parameter (out of a, Kp, Km, N and D). The large range of uncertainty is connected to the cancellation noted in (2) above. A range of 20% in tr results in virtually the whole range of tb we simulated (See Figure 3a).
  4. The curvature of the rising phase depends on whether the ratio of b/a is moderate or large (in the range of 10 to 100). If very accurate mepc waveforms should become available, the rising phase could give information on b/a.