Secondary Logo

Journal Logo

Anesthetic Pharmacology: Research Report

Simulation of the Kinetics of Neuromuscular Block

Implications for Speed of Onset

Dilger, James P. PhD

Author Information
doi: 10.1213/ANE.0b013e31827ee17f
  • Free

Although the vertebrate neuromuscular junction is the most studied biological synapse, its behavior in the presence of neuromuscular blocking drugs is not fully understood. One outstanding issue is the relationship between the speed of onset of paralysis and the potency of nondepolarizing muscle relaxants. Kopman1 quantified the relationship between onset speed and potency in humans. Gallamine, a low potency relaxant (ED95 = 2.4 mg/kg), has a 2-fold faster onset than pancuronium, a high potency relaxant (ED95 = 0.07 mg/kg). This pattern is seen with most currently used nondepolarizing muscle relaxants.2–4 A similar trend is seen for onset of paralysis by 20 relaxants in cats.5

Two explanations have been proffered to account for this: pharmacokinetics6–11 and buffered diffusion.5,12,13 Pharmacokinetic explanations consider that the onset time is primarily determined by factors such as drug elimination rate and blood perfusion of muscle tissue. Indeed, the efficiency of blood perfusion of different muscles may account for differences in onset time.6,11 The buffered diffusion explanation considers that the time needed for drug diffusion into the neuromuscular junction (the effect compartment) is rate limiting. This is termed buffered diffusion because of the high density of drug binding sites (the muscle nicotinic acetylcholine receptor, nAChR) within the narrow gap between nerve and muscle. Drug molecules diffusing into the gap have a high probability for binding to receptors; the receptors buffer the free drug concentration within the neuromuscular junction. A quick calculation reveals the magnitude of this effect. The density of nAChRs within the junction is approximately 10,000 µm−2. Consider a 1 µm2 muscle membrane surface area within the 0.05 µm junctional gap. For a free concentration of 1 µM drug molecules (corresponding to an intermediate potency drug), there will be

within this volume. However, 9000 additional drug molecules must enter the volume to bind to and inhibit 90% of the AChRs. A more potent relaxant, administered at a correspondingly lower dose, will provide a smaller flux of drug molecules into the synapse and require a longer equilibration time.

The time needed for equilibration by buffered diffusion was examined by experiments in which relaxant was perfused directly onto isolated frog neuromuscular junctions.14,15 This approach circumvents pharmacokinetic steps. Onset times were assessed by measuring the time course of muscle depolarization in response to direct application of ACh. The 50% onset times were in the range of 0.3 to 3.6 seconds over a 42-fold range of relaxant potencies.15

A theoretical analysis of buffered diffusion into the neuromuscular junction has been performed but this assumed rapid binding of antagonists to receptors and did not consider synapse geometry.14 Diffusion of molecules within a geometrically realistic synapse but without receptor binding was explored using Monte Carlo simulation.16 Attempts have been made to incorporate buffered diffusion into pharmacokinetic models.7,10,12,17 Measurement of the kinetics of antagonist (muscle relaxant) binding and dissociation for AChRs18–20 have facilitated Monte Carlo simulations to examine factors that affect buffered diffusion into a simple rectangular slab model of a synapse.21 The current paper uses Monte Carlo simulation using a realistic computer model of the mammalian neuromuscular junction22 to provide a rigorous analysis of the time course of buffered diffusion. The aim is to make quantitative predictions about the role of antagonist affinity in onset of paralysis and to compare these with clinical observations.


Monte Carlo simulations with an instantaneous concentration change were performed using MCell 3.1 (Pittsburgh Supercomputing Center, Pittsburgh, PA; Monte Carlo Simulator of Cellular Microphysiology,,23 running on an Apple MacBook Pro with 2.4 GHz Intel Core 2 Duo (Mac OSX 10.5.6, Apple, Cupertino, CA ). A typical simulation required approximately 24 hours. Simulations using a time-dependent concentration change were performed on the Stony Brook University Seawulf Supercomputer Cluster. These simulations required up to 3 weeks to complete. Molecule state occupancy values were stored every 0.01 to 1.0 milliseconds.

A realistic model of the neuromuscular junction based on electron micrographs of rat diaphragm muscle was used23–25 (Fig. 1A). This represents a 3 µm long (x-dimension) by 2 µm wide (y-dimension) by 1 µm high (z-dimension) slice from a junction that extends further in both x-directions (dashed red lines). The muscle surface (blue) has multiple folds giving a total surface area of approximately 32 µm2. AChRs are distributed on the muscle membrane at a density of 10,000 µm−2 on the tops and the upper third of the folds and a density of 2970 µm−2 on the middle third.26 The bottom third contains only the voltage-gated sodium channels that generate the muscle action potential (not modeled). The nerve surface (yellow) has an area of 5.7 µm2 and includes a 6 × 5 array of ACh-containing vesicles above the muscle folds (orange spheres). The ACh release time course is modeled by passive diffusion of 6200 ACh molecules through an expanding pore.22 Acetylcholinesterase (AChE) molecules (density 1800 µm−2) are distributed over the basal lamina between the nerve and muscle membranes27 (not shown): molecules diffuse through the basal lamina without hindrance.

Figure 1:
A, Realistic model of the neuromuscular junction. This model23 , 24 represents a slice of a neuromuscular junction and is based on electron micrographs The junction would actually continue along the x-axis as indicated by the dashed red lines. The overall dimensions are 3 × 2 × 1 µm and the gap between the nerve and muscle is approximately 0.05 µm. The upper (green) regions of the front and rear x-z surfaces are clamped to the desired antagonist concentration; antagonists enter the synapse through these regions as indicated by the arrows. These regions absorb molecules of acetylcholine (ACh). All other surfaces reflect ACh and antagonists. The nerve surface (yellow) includes a 6 × 5 array of vesicular release sites for ACh above the muscle clefts (orange spheres). Receptor molecules are located on the tops and upper parts of the junctional folds (blue). Acetylcholine esterase molecules are located on a basolaminar matrix between the nerve and muscle membranes (not shown). B, Acetylcholine receptor (AChR) kinetic scheme. The acetylcholine receptor (R) (R represents the unliganded state) has 2 binding sites for acetylcholine (A) and competitive antagonists (C). The sole conducting state in this model ARA* occurs when 2 ACh molecules bind and the receptor undergoes a conformational change (rate constant β from ARA to ARA*). Open channels usually close by returning to state ARA, but can also undergo another conformational change to a desensitized state ADA. All of the kinetic rate constants have been determined experimentally and are listed in Table 1. Inset: AChE kinetic scheme. ACh binds to a high affinity site on the esterase (E → AE) but can also bind with low affinity to the acetate-bound state of the enzyme (AcE → AAcE). Choline (Ch) is generated at the hydrolysis step (AE → AcE). Acetate (Ac) is released during the process AAcE → AE. The simulations do not track choline or acetate. The rate constants are listed in Table 1.

The gap between nerve and muscle is approximately 0.05 µm. The front and rear x-z surfaces of this gap (green) are considered to be in contact with an external reservoir of an antagonist at concentration [C]. Antagonists diffuse into this gap (green arrows) but cannot directly enter the folds on the x-z surfaces because these are within the muscle. The y-z surfaces of the gap are considered to be in contact with adjacent slices of the neuromuscular junction (dashed red lines): those surfaces reflect diffusing molecules and external antagonists cannot enter through them.

Figure 1B shows the model describing the interactions of ACh (A) and antagonist (C) with the AChR (R). The rate constants (Table 1) are well established experimentally. The AChR has 2, nonidentical sites for ligand binding.34 Occupancy of either site (or both sites) by antagonist prevents channel opening. The range of antagonist dissociation rates used is indicated in Table 1. These kinetic parameters are based on experimental values for embryonic18,19 and adult35 mouse receptors. Many antagonists are selective between the 2 binding sites.30,31,34,36 A range of values for the high L1eq and low L2eq affinity sites were considered (Table 1). “Standard” conditions refer to




is the antagonist dissociation rate for the high affinity site, L1eq and L2eq are the equilibrium constants for antagonists at the high and low affinity sites respectively and [C] is the antagonist concentration.

Table 1:
Parameters Used in the Simulations

The kinetics of activation of adult mouse AChR were used.28 ACh does not discriminate between the 2 binding sites. After 2 molecules of ACh bind, the receptor can undergo a conformational change to the open state, ARA*. Desensitization proceeds primarily from the open state.29,37Figure 1B (inset) describes hydrolysis of ACh to choline and acetate by AChE.32 The fates of choline and acetate are not simulated.

A simulation time step of 0.2 microseconds was used; the average lifetime of every molecule was >50 time steps. Some simulations were repeated with a 0.1 microseconds step time; these did not produce significantly different results. Simulations were performed using 2 different concentration profiles. In both cases, the initial conditions were that all 128,000 receptors were unliganded and there were no free ligands inside the synapse. (1) The drug concentration immediately outside the neuromuscular junction, [Csynapse], was clamped instantaneously to some nonzero value. (2) The results of pharmacokinetic modeling11,38 were used to provide the external concentration as a function of time. The plasma concentration as a function of time (Fig. 2A) was based on published measurements after a bolus injection11 and a linear interpolation between measured points. The Sheiner model38 was used to calculate the effect-site concentration assuming a drug elimination rate of 0.01 min−1 and keo = 0.15, 0.30, or 0.60 min−1 (Fig. 2B). In the Sheiner model, keo is the rate of movement of drug from plasma to the neuromuscular junction, the effect site. In this paper, I use Figure 2B to represent the time dependence of drug concentration immediately outside the synapse. Diffusion into the synapse is then determined by Monte Carlo simulation.

Figure 2:
A, The plasma concentration of vecuronium (open circles) is based on Figures 1 and 5 of reference.11 The solid line is the plasma concentration profile assumed for the simulations. B, The relative concentration of drug at the effect site predicted by the profile in panel A using the Sheiner model38 with keo = 0.15, 0.3, and 0.6 min−1. These profiles are used as the input for simulations in which the drug concentration immediately outside the synapse.

Simulations were divided into approximately 100 segments of varying duration with a constant drug concentration. The fractional drug concentration for each segment was changed from the previous value by ±0.01 so as to model the synaptic concentration profile of the Sheiner model (Fig. 2B). In all simulations, the fraction of unliganded receptors R(t) was converted to twitch amplitude, T1(t), using:

where γ is the Hill factor describing the steepness of the twitch-concentration curve and R50 is the fraction of unliganded receptors at 50% twitch suppression. Literature values6 of γ = 6.16 and R50 = 0.09 were used, and this produced a minimal twitch amplitude of 0.1 (used in all simulations). This value of R50 when combined with the standard simulation parameter values (Table 1) gives an effective concentration at 50% twitch height, EC50, of 178 nM.

Simulations proceeded until the fractional equilibrium receptor occupancy, R, was achieved.

The antagonist concentration was chosen so that, at equilibrium, approximately 94% of the receptors have at least 1 antagonist molecule bound, that is, R = 0.0625. This corresponds to a steady-state twitch amplitude of 0.1. Table 2 lists the published clinical parameters, the time to 50% twitch height, and EC50 that are used for comparison with the simulations in this paper. Analysis calculations were performed using Igor Pro 6.04 (WaveMetrics Inc., Lake Oswego, OR).

Table 2:
Clinical Onset Times and Concentrations Used for Comparison with Simulations


Model Validation

The ability of the morphological model (Fig. 1A) and MCell, to describe the time course and variability of miniature endplate currents (mepc, the current produced by the release of ACh from 1 vesicle), has been demonstrated.22,25 To validate the model with respect to competitive antagonism, mepcs were simulated under conditions of different receptor densities. The results were compared with published experimental data from mouse diaphragm muscle where (+)-tubocurarine was used to reduce the effective receptor density.51Figure 3A shows simulated mepcs using the standard receptor density of 10,000 µm−2 and a density of 2000 µm−2. This illustrates one aspect of the safety margin at the neuromuscular junction: when only 20% of the receptors remained, the mepc amplitude was still 40% of control. There was excellent agreement between simulations and experimental results both when AChE was active (Fig. 3B) and inhibited (not shown). The simulations also recapitulated the decrease in decay time constant seen with increasing concentrations of (+)-tubocurarine when AChE is inhibited51 (not shown). The high safety margin is not diminished by altering the parameters in the model. Figure 3C shows the effect of 2-fold increases or decreases in receptor gating kinetics, ACh binding, the amount of ACh contained in a vesicle, and the receptor density. The figure shows the ratio of mepc amplitude when 50% of the receptors are inhibited to the mepc amplitude without inhibition. All of these changes preserve a margin of safety because the ratio is >0.5.

Figure 3:
A, Examples of miniature endplate currents (mepcs) simulated by the release of a single vesicle containing 6200 acetylcholine (ACh) molecules. Two receptor densities are shown. The black lines are the averages of 9 separate mepcs from 9 of the 30 release sites, the gray lines are the 9 individual mepcs. B, The dependence of mepc amplitude on receptor density from published experimental data51 (solid circles) and from simulations (open circles and solid line). The error bars correspond to the standard deviations of 9 simulations. The downward curvature of the mepc amplitude versus receptor density graph is one manifestation of the margin of safety at the neuromuscular junction. In particular, when 50% of the receptors were inhibited, the mepc amplitude was still 70% of control (dashed lines). C, The effect of changing model parameters on the margin of safety. For each condition, the margin of safety is assessed as the ratio of mepc amplitude when 50% of the receptors are inhibited to that when no receptors are inhibited. Using the standard parameters (Table 1), this is 0.70 ± 0.08 (dashed line, also refer to panel B). Separate simulations were done in which a single parameter was divided by 2 (light gray bars) and multiplied by 2 (dark gray bars). The parameters were the channel opening rate (β), the channel closing rate (α), the ACh dissociation rate (k−1), the ACh content of the vesicle and receptor density. A high margin of safety (>0.5) is maintained for all of these changes.

Instantaneous Change of Drug Concentration Near the Synapse

Figure 4 illustrates a simulation using the standard antagonist (Table 1). This EC50 of 178 nM closely corresponds to the properties of pancuronium at human AChRs.31,35 Equilibration of the synapse with an antagonist concentration of 237 nM results in 94% occupancy of the receptors and an endplate current (epc, the current produced by release of ACh from all 30 vesicles) whose amplitude is 10% of control. Figure 4A shows the time course of state occupancy for the unliganded receptor, R, and the antagonist-bound states, RC, RC and CRC. Figure 4B shows the time course of the free antagonist concentration within the synapse. The figures illustrate 2 general features of buffered diffusion. (1) Equilibration is limited by the time required for sufficient antagonist molecules to enter the synapse. Here, 168,500 antagonist molecules are required (64,000 for CR, 2 × 49,000 for CRC, and 6,500 for RC). The initial flux of antagonist into the synapse is 45,000 s−1. If this flux persisted (it does not because the concentration gradient dissipates with time), approximately 4 seconds would be needed to satisfy the receptors. In the absence of receptors, the synapse would equilibrate with 220 antagonist molecules in approximately 5 milliseconds. Thus, the presence of the receptors slows the effective diffusion of the antagonist by a factor of approximately 1000. (2) With buffered diffusion, the time required for 50% receptor occupancy is shorter than for 50% free antagonist equilibration (3.6 seconds vs 9.3 seconds in Fig. 4, A and B respectively). This illustrates the high buffering capacity of the receptors. Figure 4C shows the time dependence of twitch amplitude calculated from Equation (1). Due to the steep dependence of twitch on receptor occupancy (γ = 6.16) and the high margin of safety (R50 = 0.09), twitch depression occurred over a narrow range of unliganded receptors (R). It required 4.2 seconds to go from 80% to 20% twitch; during this time, the number of unliganded receptors decreased from 12% to 7.6%. Twitch amplitude decreased to 50% at 12.6 seconds.

Figure 4:
Simulation for an instantaneous change of drug concentration immediately outside the synapse. A, Simulation of the receptor state occupancy (number of receptors) versus time for an antagonist with effective concentration at 50% twitch height (EC50) = 178 nM (see Table 1 for other parameter values). The labels are: R = unliganded receptors, CR = receptors with an antagonist bound to the low affinity site, RC = receptors with an antagonist bound to the high affinity site, CRC = receptors with antagonists bound to both sites (see Fig. 1B). At t = 0 s, the concentration of antagonist outside the synapse was clamped to 237 nM. Half of the receptors have at least 1 molecule of antagonist bound by t = 3.6 s (dashed line) and equilibrium is reached by approximately 22 s when 94% of the receptors were occupied. At this point, the endplate current was 10% of control. B, Antagonist in the synapse as a function of time. The free concentration of antagonist is 50% maximum at t = 9.3 s (dashed line) and saturates at 221 molecules corresponding to 237 nM. The gray line shows all data points at 1-ms intervals; the black line is a running average of 100 points. C, Twitch amplitude (black line) calculated from the simulated fraction of unliganded receptors (gray line) using Equation 1. The times to 80%, 50%, and 20% twitch amplitude are indicated by the dashed lines. D, The time dependence of twitch amplitude for drugs of different affinity (L1eq) values ranging from 10 to 3000 nM. The EC50 values are indicated on each trace. The drug concentration for each simulation is chosen to produce the same equilibrium conditions (94% of receptors blocked, 10% twitch). Dashed lines indicate the times to 50% twitch.

Dependence on Antagonist Affinity, Selectivity, and Concentration

Because the rate of buffered diffusion is determined by the flux of antagonist molecules into the synapse, a higher affinity antagonist is expected to exhibit a longer equilibration time. Such an antagonist would be administered at a lower concentration, giving rise to a smaller concentration gradient and a lower flux. Figure 4D shows the results of simulations with different values of antagonist dissociation rate,

, (1–300 s−1). The antagonist concentration, [C], was changed in proportional to

so that that the equilibrium occupancy of R and equilibrium twitch were the same for each of the simulations. The corresponding EC50 values are 59 to 17800 nM; encompassing the range of values for clinically used muscle relaxants (Table 2). The logarithmic time scale of Figure 4D shows that a 3-fold decrease in dissociation rate led to a 3-fold increase in time to 50% twitch.

So far, the simulations considered situations for which the antagonist was selective for one binding site over the other by a factor of 10. The site selectivity of muscle relaxants ranges from nonselective to >100-fold selective.31 For a nonselective drug, L1eq = L2eq = 79 nM and [C] = 3 × L1eq = 237 nM produces the same equilibrium level of receptor occupancy as the “standard” antagonist and concentration (Table 1). Buffered diffusion of a nonselective drug is slower than for a highly selective drug. Time to 50% receptor occupancy was 4.2 seconds (compared with 3.6 seconds), and the time to 50% twitch was 15 seconds (compared with 12.6 seconds). The slower equilibration time is due to 1.4-fold more bound drug molecules being required, mostly due to a larger number of doubly bound receptors (CRC), and these must be supplied by the same external concentration. Increasing the selectivity factor from 10- to 100-fold did not produce a significant change in twitch onset times.

For a given affinity antagonist, equilibration time depended on the concentration of the antagonist used; higher concentrations led to faster equilibration. For a simple rectangular box model of a synapse, it was shown with high concentrations of antagonist (>L1eq) that the reciprocal of the time to 50% equilibration is directly proportional to antagonist concentration.21 This also applies to the model used in this paper (not shown).

Time-Dependent Change in Plasma Concentration of Drug

After a bolus injection of a drug, the plasma concentration increases to a peak and fluctuates as the drug circulates through the body.52 A 2-compartment model has been used to describe the onset of paralysis by vecuronium.11,38,53 In this model, the plasma concentration of the drug varies as shown in Figure 2A. The drug concentration at the synapse is determined by the movement of drug from plasma into the effect compartment (Fig. 2B). I assumed that this relationship describes the drug concentration immediately outside the synapse and used this to simulate the diffusion of drug into the neuromuscular junction per se.

Figure 5 shows the results of a simulation using a concentration change profile corresponding to a slowly perfused muscle, the adductor pollicis, (keo = 0.15 min−1)11 and the “standard” antagonist (EC50 = 178 nM). Figure 5B shows how the drug concentration in the interior of the synapse (dark line) did not increase as quickly as the external drug concentration (thin line) for the first 70 seconds. A similar lag is seen (Fig. 5A) in the simulation of unliganded receptors (R, dark line) compared with what would be expected if buffered diffusion did not occur (dashed line). However, by the time the number of unliganded receptors decreased to levels that produce twitch inhibition (Fig. 5C), the lag was nearly gone. Buffered diffusion caused an increase of only 12 seconds in 50% twitch onset time. A larger shift was observed for more potent drugs, 18 seconds for EC50 = 59 nM and 80 seconds for EC50 = 18 nM (Fig. 5D).

Figure 5:
Simulation for a time-dependent change of drug concentration in the plasma. A, The number of unliganded receptors as a function of time determined from the simulations (heavy line) with keo = 0.15 min−1 and the situation for no buffered diffusion (thin, jagged line). The antagonist effective concentration at 50% twitch height (EC50) = 178 nM and other parameters were the same as those used in Figure 4. B, Free antagonist molecules in the synapse as a function of time. The simulation results are displayed as a 100-point sliding average (heavy line). The thin line represents the antagonist concentration immediately outside the synapse and also the effect site concentration in the absence of buffered diffusion. C, Twitch amplitude calculated from the number of unliganded receptors. The results (100-point sliding average) are shown as a heavy line. The time dependence of twitch in the absence of buffered diffusion is shown as a thin jagged line. The dashed lines indicate the times to 80%, 50%, and 20% twitch amplitude. In this simulation, buffered diffusion delayed the time to 50% twitch by 12 s (6%). D, Twitch amplitude resulting from simulations with 2 higher potency antagonists. The effect of buffered diffusion increases with drug potency. For EC50 = 59 nM, the delay was 18 s (9%). For EC50 = 18 nM, the delay was 80 s (42%). The drug concentration for each simulation is chosen to produce the same equilibrium conditions (94% of receptors blocked, 10% twitch). The thin jagged represents the time dependence of twitch in the absence of buffered diffusion.

Figure 6 summarizes the results of simulations with 4 different values of keo (keo = ∞ corresponds to an instantaneous concentration change at the synapse, Fig. 4) and compares them with clinical onset times. The contribution of buffered diffusion can be seen as the increase in onset time from the high concentration (low potency) asymptote of each curve. Although an instantaneous concentration change at the synapse produced the steepest concentration dependence of onset times (keo = ∞), these times were much faster than clinical onset. For the slowest concentration versus time profile (keo = 0.15 min−1), onset times were essentially determined by the rate of drug perfusion of the muscle; buffered diffusion played a very small role except for an antagonist with EC50 = 28 nM. Higher equilibration rates (keo = 0.3 and 0.6 min−1) resulted in a greater fractional increase in onset time due to buffered diffusion for the most potent drugs. However, no single value of keo can account for all of the clinical data.

Figure 6:
A comparison of clinical onset times with onset times derived from simulations. In the simulations, the time course of drug concentration immediately outside the synapse followed 4 different profiles as indicated by the values of keo, the rate constant for movement of drug from plasma to immediately outside the synapse. keo = ∞ indicates an instantaneous concentration change outside the synapse. The drugs shown are C = cisatracurium, M = mivacurium, V = vecuronium, P = pancuronium, Ro = rocuronium, T = (+) tubocurarine, Ra = rapacuronium, and G = gallamine. The symbols for mivacurium and vecuronium have been separated for clarity; their onset times are nearly identical (Table 2).

Effects of Synapse Geometry

The model of Figure 1 was constructed from measurements of rat neuromuscular junction. Equally detailed models for human neuromuscular junction are not available, so it is essential to determine the extent to which differences in anatomy may alter the results of buffered diffusion simulations. Simulations were performed with the width of the nerve-muscle contact (nominally about 2 µm, y-dimension of Fig. 1) being doubled or tripled while the receptor density was held constant. A value of keo = 0.6 min−1 was chosen for these simulations because this approximates the onset speed for low potency drugs (Fig. 6). Increasing the nerve-muscle contact width slowed the rate of antagonist binding to receptors (Fig. 7A), the rate at which the free antagonist concentration equilibrated (Fig. 7B) and the rate of twitch suppression (Fig. 7C). Figure 7D summarizes the results of simulations using a range of EC50 values and compares these with the clinical data. Onset times for nerve-muscle contact widths of 4 or 6 µm (2× or 3×) encompass most of the clinical onset times. It would be expected that a similar result would be seen if the distance between the nerve and muscle (z-dimension in Fig. 1) were decreased by a factor of 2 or 3. Conversely, a decrease in nerve-muscle contact area or an increase in the width of the synaptic cleft would be expected to decrease the importance of buffered diffusion.

Figure 7:
The effect of the width of the nerve-muscle contact on buffered diffusion into the synapse. A, The relative number of unliganded receptors as a function of time for simulations with effective concentration at 50% twitch height (EC50) = 178 nM (the standard conditions of Table 2) and keo = 0.6 min−1. The thin, jagged line indicates how this varies in the absence of buffered diffusion. The nerve-muscle contact width was 1, 2, or 3 times the standard width (2 µm, Fig. 1A). Because the receptor density is constant, the number of receptors in the 2× and 3× simulations was larger by a factor of 2 and 3, respectively. The narrow gap between the nerve and muscle cells limits the flux of antagonists into the synapse, and this translates into a longer receptor equilibration times. B, The free agonist concentration within the synapse as a function of time for the 3 different nerve-muscle contact widths. As the width was increased, the drug required a longer time to equilibrate. The thin line represents the antagonist concentration immediately outside the synapse and also the effect site concentration in the absence of buffered diffusion. C, Twitch amplitude calculated from the number of unliganded receptors. The wider the synapse, the longer it takes to depress twitch amplitude. The time to 50% twitch amplitude and the percent of this time caused by buffered diffusion was: no buffered diffusion (68 s, 0%), 1× width (77 s, 13%), 2× width (90 s, 32%), and 3× width (118 s, 74%). D, A comparison of clinical onset times with onset times derived from simulations. The drug abbreviations are given in the legend to Figure 6. The simulations suggest that if the nerve-muscle contact width were between 2 and 3 times the standard width (Fig. 1A), buffered diffusion could account for the differences in clinical onset times assuming a single value of keo.


The goal of this investigation was to determine whether buffered diffusion might be expected to affect the onset of paralysis by muscle relaxants. I considered a range of muscle relaxant potencies. I considered 4 scenarios for the rate at which drug arrives immediately outside the synapse (neuromuscular junction): instantaneously (Fig. 4), and with the rate of transfer variable (keo) being 0.15, 0.3, and 0.6 min−1 (Figs. 5 and 6). I considered a standard model of the rat neuromuscular junction (Fig. 1A) and some modifications of this morphology (Fig. 7). I conclude that for keo = 0.6 min−1 and a 2- to 3-fold wider (than the standard rat model) nerve-muscle contact width, buffered diffusion could account for the differences in paralysis onset times among clinically used muscle relaxants (Fig 7D). Other morphological differences, particularly a smaller gap between nerve and muscle, may act in conjunction with the contact width to increase the role of buffered diffusion.

Critique of the Model

A complex computer model was used in this study. Modeling is crucial for understanding phenomena that occur under nonequilibrium conditions. One must, however, assess the sensitivity of the results to the assumptions used. The kinetic model of muscle nAChR activation, desensitization, and competitive inhibition (Fig. 1B) is well supported by experiments. The kinetic parameters are appropriate for room temperature rather than 37°C; very little kinetic data are available for 37°C. Increasing temperature decreases antagonist affinity20 and increases diffusion coefficients. Thus, buffered diffusion would be faster at higher temperatures.

The morphological model of the neuromuscular junction (Fig. 1A) is based on rat diaphragm muscle. Vertebrate species differ in the nerve-muscle contact area, degree of muscle membrane folding, quantal content, distribution of contact areas, and transverse width of individual nerve-muscle contacts.54,55 These studies indicate that the human neuromuscular junction consists of a 20 to 30 µm wide cluster of strands. However, each strand within the cluster may be only several microns in width. The extent to which Schwann cells extend over the cluster is unknown but even if they formed a narrow barrier around the strands, this would not significantly impede the movement of antagonist molecules because the path would not be lined with receptors. I performed simulations assuming a range of widths for the nerve-muscle contact region (Fig. 7). The human neuromuscular junction has 7 to 8 times more muscle area than nerve area.54 Thus, the model, with an area ratio of 5.7:1, is reasonable. Species differences in nerve-muscle contact area and in quantal content may compensate for each other.54 The length of the contact region (i.e., the extent along the x-direction of Fig. 1A) does differ among species, but this does not affect the time required for antagonists to diffuse into the interior as they enter only through the narrow upper regions of the front and rear x-z surfaces (green surface and arrows in Fig. 1A).

The simulations show that the onset time for receptor equilibration was inversely proportional to L1eq when the concentration change was instantaneous (Fig. 4D). This indicates that antagonist diffusion into the synapse is strongly buffered by the receptors. Deviations from this rapid binding scenario occur when (a) the diffusion coefficient is 100-fold larger, (b) the antagonist association constant is 100-fold smaller, (c) the receptor density is 100-fold smaller, or (d) the synaptic gap is 10-fold larger.21 None of these alternate scenarios appears very likely.

The models used in some pharmacokinetic/pharmacodynamic studies appropriately consider buffering effect of receptors on the concentration of drug in the neuromuscular junction. However, they did not consider either the way in which the narrow junction limits the flux of drug molecules or the kinetics of drug association and dissociation. In these models, the receptor concentration is a fitting parameter and the resulting values grossly underestimate the actual receptor concentration (280 nM7,12 or 1.5 µM10 compared with approximately 320 µM, calculated from 104 receptors/µm2/0.05 µm). This paper demonstrated that Monte Carlo simulation of diffusion and pharmacodynamics is a robust method for determining the time course of receptor binding in a realistic synapse.

Buffered Diffusion with an Instantaneous Concentration Change Near the Synapse

The first set of simulations was performed using the assumption that the volume surrounding the synapse was clamped instantaneously to a given concentration of antagonist. The slope of the log–log plot of onset time against L1eq is −1.0 (not shown). Experimentally, concentration clamp conditions were approximated with iontophoretic application of antagonists and ACh to isolated frog neuromuscular junctions.14,15 (+)-Tubocurarine applied at 4 × Leq (Leq = 500 nM for frog muscle AChR) inhibited endplate potentials with a time constant of 0.5 seconds.14 This is within a factor of 2 of our simulated value. Similar experiments using 7 antagonists with a 50-fold range in potency found the slope of a log(onset time) versus log (Leq) plot between −0.7 and −0.8.15 This is close to the theoretical relationship especially considering the experimental difficulties and the potential confounding effect of different binding site selectivities. In contrast, the dependence of onset time for human paralysis on L1eq is weaker, giving a slope of −0.22 on a log–log plot (not shown).

Buffered Diffusion with a Time-Dependent Change in Plasma Concentration

The antagonist concentration immediately outside the synapse was calculated from published values of plasma vecuronium concentrations11 and a 2-compartment pharmacokinetic model38 (Fig. 2). The simulations indicate that buffered diffusion slows the initial rate of equilibration of antagonist molecules with the receptor (Fig. 5B). Due to the high safety margin, twitch depression does not occur until t >80 seconds. By this time, however, buffered diffusion had only a small effect on receptor occupancy, antagonist concentration and twitch (Fig. 5, A–C). The effect was significantly increased for simulations using a 3- to 10-fold higher potency drug (Fig. 5D). Figure 6 shows that there is little agreement between simulated and experimentally observed twitch onset times. In particular, simulations of buffered diffusion using any single value of keo cannot account for the differences among the drugs.

The morphological model of Figure 1, based on rat neuromuscular junction, may not faithfully represent the geometry of the human junction. Modification of the model to allow for a 2- or 3-fold wider extent of nerve-muscle contact provided a better congruence between simulated and observed onset times (Fig. 7D). A similar result would be obtained if it were assumed that the width of the junction was unchanged, but the width of the synaptic cleft was decreased by a factor of 2 or 3. Results from simulations assuming a 2- to 3-fold wider synapse (4–6 µm), and keo = 0.6 min−1 encompass all of the clinical data (Fig. 7D). Until more specific anatomical information about the human neuromuscular junction is available, it will not be possible to know for certain whether buffered diffusion can account for the dependence of onset time on drug potency.


Name: James P. Dilger, PhD.

Contribution: This author designed and conducted the study, analyzed the data, and wrote the manuscript.

Attestation: James Dilger has seen the original study data, reviewed the analysis of the data, approved the final manuscript, and is the author responsible for archiving the study files.

This manuscript was handled by: Tony Gin, FANZCA, FRCA, MD.


The author thanks Dr. Aaron Hartman, graduate of the Medical School Program, Stony Brook University, for assistance with writing and testing MCell3 model files; Dr. Jeevendra Martyn, Department of Anesthesia, Critical Care and Pain Medicine, Massachusetts General Hospital for his helpful advice on an early version of the manuscript; and Dr. Mihai Sadean, Department of Anesthesiology, Stony Brook University for use of his spreadsheet to make pharmacokinetic/pharmacodynamics calculations.


1. Kopman AF. Pancuronium, gallamine, and d-tubocurarine compared: is speed of onset inversely related to drug potency? Anesthesiology. 1989;70:915–20
2. Kopman AF, Klewicka MM, Kopman DJ, Neuman GG. Molar potency is predictive of the speed of onset of neuromuscular block for agents of intermediate, short, and ultrashort duration. Anesthesiology. 1999;90:425–31
3. Kopman AF, Klewicka MM, Ghori K, Flores F, Neuman GG. Dose-response and onset/offset characteristics of rapacuronium. Anesthesiology. 2000;93:1017–21
4. Kopman AF, Klewicka MM, Neuman GG. Molar potency is not predictive of the speed of onset of atracurium. Anesth Analg. 1999;89:1046–9
5. Bowman WC, Rodger IW, Houston J, Marshall RJ, McIndewar I. Structure:action relationships among some desacetoxy analogues of pancuronium and vecuronium in the anesthetized cat. Anesthesiology. 1988;69:57–62
6. Bragg P, Fisher DM, Shi J, Donati F, Meistelman C, Lau M, Sheiner LB. Comparison of twitch depression of the adductor pollicis and the respiratory muscles. Pharmacodynamic modeling without plasma concentrations. Anesthesiology. 1994;80:310–9
7. Bhatt SB, Amann A, Nigrovic V. Onset-potency relationship of nondepolarizing muscle relaxants: a reexamination using simulations. Can J Physiol Pharmacol. 2007;85:774–82
8. Beaufort TM, Nigrovic V, Proost JH, Houwertjes MC, Wierda JM. Inhibition of the enzymic degradation of suxamethonium and mivacurium increases the onset time of submaximal neuromuscular block. Anesthesiology. 1998;89:707–14
9. Proost JH, Wierda JM, Meijer DK. An extended pharmacokinetic/pharmacodynamic model describing quantitatively the influence of plasma protein binding, tissue binding, and receptor binding on the potency and time course of action of drugs. J Pharmacokinet Biopharm. 1996;24:45–77
10. De Haes A, Houwertjes MC, Proost JH, Wierda JM. An isolated, antegrade, perfused, peroneal nerve anterior tibialis muscle model in the rat: a novel model developed to study the factors governing the time course of action of neuromuscular blocking agents. Anesthesiology. 2002;96:963–70
11. Fisher DM, Szenohradszky J, Wright PM, Lau M, Brown R, Sharma M. Pharmacodynamic modeling of vecuronium-induced twitch depression. Rapid plasma-effect site equilibration explains faster onset at resistant laryngeal muscles than at the adductor pollicis. Anesthesiology. 1997;86:558–66
12. Donati F, Meistelman C. A kinetic-dynamic model to explain the relationship between high potency and slow onset time for neuromuscular blocking drugs. J Pharmacokinet Biopharm. 1991;19:537–52
13. Donati F. Onset of action of relaxants. Can J Anaesth. 1988;35:S52–8
14. Armstrong DL, Lester HA. The kinetics of tubocurarine action and restricted diffusion within the synaptic cleft. J Physiol (Lond). 1979;294:365–86
15. Min JC, Bekavac I, Glavinovic MI, Donati F, Bevan DR. Iontophoretic study of speed of action of various muscle relaxants. Anesthesiology. 1992;77:351–6
16. Lacks DJ. Tortuosity and anomalous diffusion in the neuromuscular junction. Phys Rev E Stat Nonlin Soft Matter Phys. 2008;77:041912
17. Proost JH, Houwertjes MC, Wierda JM. Is time to peak effect of neuromuscular blocking agents dependent on dose? Testing the concept of buffered diffusion. Eur J Anaesthesiol. 2008;25:572–80
18. Wenningmann I, Dilger JP. The kinetics of inhibition of nicotinic acetylcholine receptors by (+)-tubocurarine and pancuronium. Mol Pharmacol. 2001;60:790–6
19. Demazumder D, Dilger JP. The kinetics of competitive antagonism by cisatracurium of embryonic and adult nicotinic acetylcholine receptors. Mol Pharmacol. 2001;60:797–807
20. Demazumder D, Dilger JP. The kinetics of competitive antagonism of nicotinic acetylcholine receptors at physiological temperature. J Physiol (Lond). 2008;586:951–63
21. Dilger JP. Monte Carlo simulation of buffered diffusion into and out of a model synapse. Biophys J. 2010;98:959–67
22. Stiles JR, Van Helden D, Bartol TM Jr, 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–52
23. Stiles JR, Bartol TM JrDe Schutter E. Monte Carlo methods for simulating realistic synaptic microphysiology using MCell. Computational Neuroscience: Realistic Modeling for Experimentalists. 2001 Boca Raton, FL CRC Press:87–127
24. Bartol TM Jr, Land BR, Salpeter EE, Salpeter MM. Monte Carlo simulation of miniature endplate current generation in the vertebrate neuromuscular junction. Biophys J. 1991;59:1290–307
25. Stiles JR, Bartol TM, Salpeter MM, Salpeter EE, Sejnowski TJCowan WM, Südhof TC, Stevens CF. Synaptic variability. New insights from reconstruction and Monte Calro simlations with MCell. Synapses. 2001 Baltimore, MD The Johns Hopkins University Press:681–731
26. Salpeter MMSalpeter MM. Vertebrate neuromuscular junctions: general morphology, molecular organization, and functional consequences. The Vertebrate Neuromuscular Junction. 1987 New York, NY Alan R. Liss:1–54
27. Salpeter MM. Electron microscope radioautography as a quantitative tool in enzyme cytochemistry. I. The distribution of acetylcholinesterase at motor end plates of a vertebrate twitch muscle. J Cell Biol. 1967;32:379–89
28. Akk G, Auerbach A. Inorganic, monovalent cations compete with agonists for the transmitter binding site of nicotinic acetylcholine receptors. Biophys J. 1996;70:2652–8
29. Dilger JP, Liu Y. Desensitization of acetylcholine receptors in BC3H-1 cells. Pflugers Arch. 1992;420:479–85
30. Fletcher GH, Steinbach JH. Ability of nondepolarizing neuromuscular blocking drugs to act as partial agonists at fetal and adult mouse muscle nicotinic receptors. Mol Pharmacol. 1996;49:938–47
31. Liu M, Dilger JP. Site selectivity of competitive antagonists for the mouse adult muscle nicotinic acetylcholine receptor. Mol Pharmacol. 2009;75:166–73
32. Rosenberry TL. Quantitative simulation of endplate currents at neuromuscular junctions based on the reaction of acetylcholine with acetylcholine receptor and acetylcholinesterase. Biophys J. 1979;26:263–89
33. Krnjevic K, Mitchell JF. Diffusion of acetylcholine in agar gels and in the isolated rat diaphragm. J Physiol (Lond). 1960;153:562–72
    34. Sine SM, Taylor P. Relationship between reversible antagonist occupancy and the functional capacity of the acetylcholine receptor. J Biol Chem. 1981;256:6692–9
    35. Dilger JP, Vidal AM, Liu M, Mettewie C, Suzuki T, Pham A, Demazumder D. Roles of amino acids and subunits in determining the inhibition of nicotinic acetylcholine receptors by competitive antagonists. Anesthesiology. 2007;106:1186–95
    36. Liu M, Dilger JP. Synergy between pairs of competitive antagonists at adult human muscle acetylcholine receptors. Anesth Analg. 2008;107:525–33
    37. Auerbach A, Akk G. Desensitization of mouse nicotinic acetylcholine receptor channels. A two-gate mechanism. J Gen Physiol. 1998;112:181–97
    38. Sheiner LB, Stanski DR, Vozeh S, Miller RD, Ham J. Simultaneous modeling of pharmacokinetics and pharmacodynamics: application to d-tubocurarine. Clin Pharmacol Ther. 1979;25:358–71
    39. Tran TV, Fiset P, Varin F. Pharmacokinetics and pharmacodynamics of cisatracurium after a short infusion in patients under propofol anesthesia. Anesth Analg. 1998;87:1158–63
    40. Roy JJ, Varin F. Physicochemical properties of neuromuscular blocking agents and their impact on the pharmacokinetic-pharmacodynamic relationship. Br J Anaesth. 2004;93:241–8
    41. Laurin J, Donati F, Nekka F, Varin F. Peripheral link model as an alternative for pharmacokinetic-pharmacodynamic modeling of drugs having a very short elimination half-life. J Pharmacokinet Pharmacodyn. 2001;28:7–25
    42. Cameron M, Donati F, Varin F. In vitro plasma protein binding of neuromuscular blocking agents in different subpopulations of patients. Anesth Analg. 1995;81:1019–25
    43. Alloul K, Whalley DG, Shutway F, Ebrahim Z, Varin F. Pharmacokinetic origin of carbamazepine-induced resistance to vecuronium neuromuscular blockade in anesthetized patients. Anesthesiology. 1996;84:330–9
    44. Cronnelly R, Fisher DM, Miller RD, Gencarelli P, Nguyen-Gruenke L, Castagnoli N Jr. Pharmacokinetics and pharmacodynamics of vecuronium (ORG NC45) and pancuronium in anesthetized humans. Anesthesiology. 1983;58:405–8
    45. Wood M, Stone WJ, Wood AJ. Plasma binding of pancuronium: effects of age, sex, and disease. Anesth Analg. 1983;62:29–32
    46. Stanski DR, Ham J, Miller RD, Sheiner LB. Pharmacokinetics and pharmacodynamics of d-tubocurarine during nitrous oxide-narcotic and halothane anesthesia in man. Anesthesiology. 1979;51:235–41
    47. Ghoneim MM, Kramer E, Bannow R, Pandya H, Routh JI. Binding of d-tubocurarine to plasma proteins in normal man and in patients with hepatic or renal disease. Anesthesiology. 1973;39:410–5
    48. Dragne A, Varin F, Plaud B, Donati F. Rocuronium pharmacokinetic-pharmacodynamic relationship under stable propofol or isoflurane anesthesia. Can J Anaesth. 2002;49:353–60
    49. Schiere S, Proost JH, Schuringa M, Wierda JM. Pharmacokinetics and pharmacokinetic-dynamic relationship between rapacuronium (Org 9487) and its 3-desacetyl metabolite (Org 9488). Anesth Analg. 1999;88:640–7
    50. Ramzan IM. Molecular weight of cation as a determinant of speed of onset of neuromuscular blockade. Anesthesiology. 1982;57:247–8
    51. Pennefather P, Quastel DM. Relation between subsynaptic receptor blockade and response to quantal transmitter at the mouse neuromuscular junction. J Gen Physiol. 1981;78:313–44
    52. Fisher DM. (Almost) everything you learned about pharmacokinetics was (somewhat) wrong! Anesth Analg. 1996;83:901–3
    53. Ducharme J, Varin F, Bevan DR, Donati F. Importance of early blood sampling on vecuronium pharmacokinetic and pharmacodynamic parameters. Clin Pharmacokinet. 1993;24:507–18
    54. Wood SJ, Slater CR. Safety factor at the neuromuscular junction. Prog Neurobiol. 2001;64:393–429
    55. Slater CR, Fawcett PR, Walls TJ, Lyons PR, Bailey SJ, Beeson D, Young C, Gardner-Medwin D. Pre- and post-synaptic abnormalities associated with impaired neuromuscular transmission in a group of patients with ‘limb-girdle myasthenia’. Brain. 2006;129:2061–76
    © 2013 International Anesthesia Research Society