Skip to main content

On the relationship between inhibition and receptor occupancy by nondepolarizing neuromuscular blocking drugs

Abstract

Background

Nondepolarizing neuromuscular blocking drugs (NDNBs) are clinically used to produce muscle relaxation during general anesthesia. To better understand clinical properties of NDNBs, comparative in vitro pharmacologic studies have been performed. In these studies, a receptor binding model, which relies on the assumption that the inhibition, i.e., the effect of an NDNB, is proportional to the receptor occupancy by the drug, has been effectively used to describe obtained experimental data. However, it has not been studied in literature under which conditions the above assumption can be justified nor the assumption still holds in vivo. The purpose of this study is to explore the in vivo relationship between the inhibition and the receptor occupancy by an NDNB and to draw implications on how in vitro experimental results can be used to discuss the in vivo properties of NDNBs.

Methods

An ordinary differential equation model is employed to simulate physiologic processes of the activation of receptors by acetylcholine (ACh) as well as inhibition by an NDNB. With this model, the degree of inhibition is quantified by the fractional amount of receptors that are not activated by ACh due to the presence of an NDNB. The results are visualized by plotting the fractional amounts of the activated receptors as a function of the receptor occupancy.

Results

Numerical investigations reflecting in vivo conditions show that the degree of inhibition is not proportional to the receptor occupancy, i.e., there is a nonlinear relationship between the inhibition and the receptor occupancy. However, under a setting of high concentration of ACh reflecting a typical situation of in vitro experiments, the relationship between the inhibition and the receptor occupancy becomes linear, suggesting the validity of the receptor binding model. Also, it is found that the extent of nonlinearity depends on the selectivity of NDNBs for the two binding sites of the receptors.

Conclusions

While the receptor binding model may be effective for estimating affinity of an NDNB through in vitro experiments, these models do not directly describe in vivo properties of NDNBs, because the nonlinearity between the inhibition and the receptor occupancy causes the modulation of the resultant concentration-effect relationships of NDNBs.

Background

Nondepolarizing neuromuscular blocking drugs (NDNBs) inhibit neuromuscular transmission by competing with acetylcholine (ACh) for binding sites at the post-junctional nicotinic acetylcholine receptors (AChRs). They are widely used during general anesthesia to produce muscle relaxation for facilitating tracheal intubation and for providing optimal surgical working conditions [1]. To describe in vivo properties of NDNBs and to understand the molecular mechanisms behind clinical observations, many in vitro studies have been conducted (see e.g., [210]). In particular, in [7, 10], in vitro experiments have been conducted through macroscopic current recordings from outside-out patches of BOSC23 cells or Xenopus oocytes where human adult (rather than mouse adult or embryonic) muscle-type AChRs were expressed. In [7], synergistic effects of pairs of NDNBs have been studied. Although it successfully identified evidence for synergy between many pairs of NDNBs at the receptor level, not all the results correlated with synergism observed in vivo. In [9, 10], it was found that the IC50, the drug concentration needed to produce a 50% inhibition of the ACh-induced current, decreases with the increase in the concentration of ACh. Although it was concluded in [9, 10] that this demonstrated the existence of a noncompetitive action of NDNBs, some researchers raised concerns [11] over the insufficiency in quantitative analysis to draw the conclusion. Thus, more investigations and considerations are needed to describe the clinical observations and to clarify underlying mechanisms of inhibition based on in vitro experimental results.

While there are several measures identifying in vitro properties of NDNBs, potency [12] is one of the most widely used measures for studying NDNBs. It is usually quantified by IC50 estimated by regression analysis using the Hill equation for the relative current Iantag/I0 given by

$$\begin{array}{*{20}l} \frac{I_{\text{antag}}}{I_{0}} = \frac{{\text{IC}_{50}}^{n_{\mathrm{H}}}}{{\text{IC}_{50}}^{n_{\mathrm{H}}} + {[\!\mathrm{D}]}^{n_{\mathrm{H}}}}, \end{array} $$
(1)

where I0 and Iantag stand for the currents measured in the absence and presence of the drugs, respectively. The parameter nH stands for the Hill coefficient, and [ D] for the drug concentration. Also, affinity [12], which is defined as the extent or fraction of drug binding to receptors at any given drug concentration, is used for characterizing the strength of the binding of a ligand to its receptors. In [48], the relative current versus concentration curves were analyzed based on the two-site receptor binding model given by the following equation:

$$\begin{array}{*{20}l} \frac{I_{\text{antag}}}{ I_{0}} = \frac{K_{\mathrm{D1}} K_{\mathrm{D2}}}{ K_{\mathrm{D1}} K_{\mathrm{D2}} + K_{\mathrm{D1}} {[\!\mathrm{D}]} + K_{\mathrm{D2}} {[\!\mathrm{D}]} + {[\!\mathrm{D}]}^{2} }, \end{array} $$
(2)

where KD1 and KD2 stand for the dissociation equilibrium constants for NDNBs binding to the first and second sites of an AChR, respectively. Since the right-hand side of Eq. (2) represents the fraction of free receptors, i.e., the fraction of receptors that are not occupied by the drugs, the model (2) implies that the potency of an NDNB can be completely characterized by the affinity of the drug. However, in general, there may exist other factors that determine the potency of a drug. The underlying assumption for the use of the model (2) is that the inhibition, i.e. the effect of an NDNB, is proportional to the receptor occupancy by the drug. While the model (2) has been statistically tested in [48], the key factors affecting the validity of this assumption have not been fully discussed in literature, and it has not been examined if the assumption holds in vivo.

The purpose of this study is to explore the in vivo relationship between the inhibition and the receptor occupancy by an NDNB and to draw implications on how in vitro experimental results can be used to discuss the in vivo pharmacologic properties of NDNBs. An ordinary differential equation model based on [5, 1315] is employed to simulate physiologic processes of the activation of receptors by acetylcholine (ACh) as well as inhibition by an NDNB. Based on the model, we examine the conditions under which the relationship between inhibition and receptor occupancy becomes nonlinear and discuss its clinically relevant implications.

Methods

Model of competition between acetylcholine and NDNB

At the neuromuscular junction, electrical impulses from motor nerve cause the release of transmitter, ACh, to the synaptic cleft. Then, part of the released ACh molecules bind to AChRs on the muscle membrane and results in a change in the conductance of the membrane due to the channel opening. This change causes movement of ions into the muscle cell that produces an action potential spreading over the surfaces of skeletal muscle fibers and causing muscle contraction. NDNBs act by competing with ACh for binding to the two sites of an AChR and preventing changes in the membrane conductance. To describe the competition between ACh and NDNB molecules for binding to AChRs, this paper uses the model developed in [15] with a modification to incorporate the dynamics of channel opening as described in [5, 13, 14]. The complexes formed by binding of ACh, denoted by A, and NDNB, by D, to AChR, by R, are represented by 3-letter symbols as shown in Fig. 1. The first and last letters denote the first and second ligands occupying the sites 1 and 2, respectively, and the middle letter represents the receptor R. Unoccupied sites are denoted by O, and ORO stands for free AChR. The kinetics of ACh and NDNB are represented using association rate constants kassocAi and kassocDi, for site #i(i=1, 2), respectively. Similarly, kdissAi and kdissDi stand for the dissociation rate constants for ACh and NDNB. The symbol ARA stands for AChRs bound with two ACh molecules but in the closed state, and the symbol ARA* for AChRs in the open state and thus activated due to the conformational change of AChRs. The rate constants of opening and closing of AChRs are represented by kopen and kclose, respectively. Then, the concentrations of these complexes can be calculated as a function of time by solving the following ordinary differential equations:

$$\begin{array}{*{20}l} & \frac{\mathrm{d}}{\mathrm{d} t}\mathrm{[\!A]} = - {k_{\text{decay}}} \mathrm{[\!A]} \\ &\hspace{13mm} + k_{\mathrm{dissA1}} (\mathrm{[\!ARA]} + \mathrm{[\!ARD]} + \mathrm{[\!ARO]}) \\ &\hspace{13mm} - k_{\mathrm{assocA1}} \mathrm{[\!A]} (\mathrm{[\!ORA]} + \mathrm{[\!ORD]} + \mathrm{[\!ORO]}) \\ &\hspace{13mm} + k_{\mathrm{dissA2}} (\mathrm{[\!ARA]} + \mathrm{[\!DRA]} + \mathrm{[\!ORA]}) \\ &\hspace{13mm} - k_{\mathrm{assocA2}} \mathrm{[\!A]} (\mathrm{[\!ARO]} + \mathrm{[\!DRO]} + \mathrm{[\!ORO]}), \end{array} $$
(3a)
Fig. 1
figure1

Diagram of the interactions between acetylcholine (A), NDNB (D), and the postsynaptic receptor (R). The complexes formed by binding of ACh and NDNB to AChR are represented by 3-letter symbols, where the first and last letters denote the first and second ligands occupying the sites 1 and 2, respectively. The state ARA* represents the activated AChR

$$\begin{array}{*{20}l} & { \frac{\mathrm{d}}{\mathrm{d} t}\mathrm{[\!ARA^{\ast}]} = - k_{\text{close}}\ \mathrm{[\!ARA^{\ast}]} + k_{\text{open}} \mathrm{[\!ARA]},} \end{array} $$
(3b)
$$\begin{array}{*{20}l} & \frac{\mathrm{d}}{\mathrm{d} t}\mathrm{[\!ARA]} = k_{\mathrm{assocA1}} \mathrm{[\!ORA][A]} - k_{\mathrm{dissA1}}\mathrm{[\!ARA]} \\ &\hspace{17mm} + k_{\mathrm{assocA2}} \mathrm{[\!ARO][A]} - k_{\mathrm{dissA2}}\mathrm{[\!ARA]} \\ &\hspace{18mm} {+ k_{\text{close}} \mathrm{[\!ARA^{\ast}]} - k_{\text{open}} \mathrm{[\!ARA]},} \end{array} $$
(3c)
$$\begin{array}{*{20}l} & \frac{\mathrm{d}}{\mathrm{d} t}\mathrm{[\!DRD]} = k_{\mathrm{assocD1}} \mathrm{[\!ORD][D]} - k_{\mathrm{dissD1}}\mathrm{[\!DRD]} \\ &\hspace{17mm} + k_{\mathrm{assocD2}} \mathrm{[\!DRO][D]} - k_{\mathrm{dissD2}}\mathrm{[\!DRD]}, \end{array} $$
(3d)
$$\begin{array}{*{20}l} & \frac{\mathrm{d}}{\mathrm{d} t}\mathrm{[\!ARD]} = k_{\mathrm{assocA1}} \mathrm{[\!ORD][A]} - k_{\mathrm{dissA1}}\mathrm{[\!ARD]} \\ &\hspace{17mm} + k_{\mathrm{assocD2}} \mathrm{[\!ARO][D]} - k_{\mathrm{dissD2}}\mathrm{[\!ARD]}, \end{array} $$
(3e)
$$\begin{array}{*{20}l} & \frac{\mathrm{d}}{\mathrm{d} t}\mathrm{[\!DRA]} = k_{\mathrm{assocD1}} \mathrm{[\!ORA][D]} - k_{\mathrm{dissD1}}\mathrm{[\!DRA]} \\ & \hspace{17mm} + k_{\mathrm{assocA2}} \mathrm{[\!DRO][A]} - k_{\mathrm{dissA2}}\mathrm{[\!DRA]}, \end{array} $$
(3f)
$$\begin{array}{*{20}l} & \frac{\mathrm{d}}{\mathrm{d} t}\mathrm{[\!ARO]} = k_{\mathrm{assocA1}} \mathrm{[\!ORO][A]} - k_{\mathrm{dissA1}}\mathrm{[\!ARO]} \\ &\hspace{17mm} + k_{\mathrm{dissA2}} \mathrm{[\!ARA]} - k_{\mathrm{assocA2}} \mathrm{[\!ARO][A]} \\ &\hspace{17mm} + k_{\mathrm{dissD2}} \mathrm{[\!ARD]} - k_{\mathrm{assocD2}} \mathrm{[\!ARO] [D]}, \end{array} $$
(3g)
$$\begin{array}{*{20}l} & \frac{\mathrm{d}}{\mathrm{d} t}\mathrm{[\!ORA]} = k_{\mathrm{assocA2}} \mathrm{[\!ORO][A]} - k_{\mathrm{dissA2}}\mathrm{[\!ORA]} \\ &\hspace{17mm} + k_{\mathrm{dissA1}} \mathrm{[\!ARA]} - k_{\mathrm{assocA1}} \mathrm{[\!ORA][A]} \\ &\hspace{17mm} + k_{\mathrm{dissD1}} \mathrm{[\!DRA]} - k_{\mathrm{assocD1}} \mathrm{[\!ORA] [D]}, \end{array} $$
(3h)
$$\begin{array}{*{20}l} & \frac{\mathrm{d}}{\mathrm{d} t}\mathrm{[\!DRO]} = k_{\mathrm{assocD1}} \mathrm{[\!ORO][D]} - k_{\mathrm{dissD1}}\mathrm{[\!DRO]} \\ & \hspace{17mm} + k_{\mathrm{dissD2}} \mathrm{[\!DRD]} - k_{\mathrm{assocD2}} \mathrm{[\!DRO] [D]} \\ &\hspace{17mm} + k_{\mathrm{dissA2}} \mathrm{[\!DRA]} - k_{\mathrm{assocA2}} \mathrm{[\!DRO][A]}, \end{array} $$
(3i)
$$\begin{array}{*{20}l} & \frac{\mathrm{d}}{\mathrm{d} t}\mathrm{[\!ORD]} = k_{\mathrm{assocD2}} \mathrm{[\!ORO][D]} - k_{\mathrm{dissD2}}\mathrm{[\!ORD]} \\ &\hspace{17mm} + k_{\mathrm{dissD1}} \mathrm{[\!DRD]} - k_{\mathrm{assocD1}} \mathrm{[\!ORD] [D]} \\ &\hspace{17mm} + k_{\mathrm{dissA1}} \mathrm{[\!ARD]} - k_{\mathrm{assocA1}} \mathrm{[\!ORD][A]}, \end{array} $$
(3j)

where [ D] stands for the concentration of NDNB at the effect site, and kdecay for the rate constant of the decay of free ACh due to hydrolysis of ACh by acetylcholinesterase and diffusion of ACh away from the synaptic cleft. The concentration [ ORO] of the unoccupied AChR is given by

$$\begin{array}{*{20}l} \mathrm{[\!ORO]} =& \mathrm{[\!R]}_{\text{total}} {- {[\text{ARA}^{\ast}]}} - \mathrm{[\!ARA]} - \mathrm{[\!DRD]} \\&- \mathrm{[\!ARD]} - \mathrm{[\!DRA]} \\ & - \mathrm{[\!ARO]} - \mathrm{[\!ORA]} - \mathrm{[\!DRO]} - \mathrm{[\!ORD]}, \end{array} $$
(4)

where [ R]total stands for the concentration of the post-junctional AChRs in the synaptic cleft.

Here we provide an example of simulation results of the model (3). The setting of the parameters is shown in Table 1. Following [15], the concentration [ R]total of AChRs was calculated using the number of AChRs (2.1×107) at the end plates of human deltoid muscle reported in [16] and the volume of the synaptic cleft (4.5×1013L) of rat diaphragm reported in [17]. The initial concentration [ A]init of the free ACh was calculated by assuming that the number of ACh molecules released is one tenth of the number of AChRs as discussed in [18, 19]. This means that the number of ACh molecules released is 2.1×106 corresponding to the release of 300 vesicles [20] with 7000 ACh molecules in each vesicle. Thus, while the concentration of ACh itself is sufficiently higher than needed for neuromuscular transmission, due to the considerable excess of AChRs, only a small fractional amount of AChRs will be activated during neuromuscular transmission. The rate constant kdecay was determined such that the half-life of free ACh becomes 58μs as calculated in [15]. The dissociation and association rate constants for ACh were tentatively set to the values obtained in [13, 14] from experiments using mouse adult AChRs. Similarly, the rate constants for NDNBs were set to tentative values based on experiments in [4, 5] using mouse embryonic and adult AChRs. Since the values of these rate constants may be different between human and mouse AChRs, these constants will be varied in a systematic way.

Table 1 Setting of parameters for numerical simulation in base case

Figure 2 shows the time courses of the concentrations of free ACh, [ A], of diliganded AChRs at the closed state, [ ARA], and of the activated AChRs, [ARA], after the arrival of an electrical impulse from motor nerve and the release of ACh at t=0. The initial concentrations for the nine complexes at t=0 were defined by the steady state concentrations before the release of ACh. The concentration of the NDNB was set to one of [ D]=0,3.0×10−8M, and 1.0×10−7M. As shown in the upper panel of Fig. 2, the concentration of ACh rapidly decreases due to the hydrolysis of ACh by acetylcholinesterase and binding of ACh to AChRs. The time course of [ ARA] is shown in the middle panel of Fig. 2. The concentration [ ARA] raises rapidly and decays after reaching the peak at around t=0.04 ms. Subsequently, as shown in the lower panel, the concentration [ARA] raises and reaches its peak at around t=0.3 ms. The highest [ARA], denoted by [ARA]max, is attained in the absence of NDNB, i.e., [ D]=0, and the peak concentration of the activated AChRs [ARA]peak decreases with the increase in the concentration of [ D].

Fig. 2
figure2

An example of simulation results of the model (3). a The time course of the concentration [ A] of free ACh. b The time course of the concentration [ ARA] of the complex ARA representing the activated AChRs

With this model, the effect of an NDNB can be quantified by a fraction of activated AChRs given by [ARA]peak/[ARA]max. This definition corresponds to the relative current Iantag/I0 used in in vitro experiments under the assumption that the membrane conductance is proportional to the number of activated AChRs, i.e. the number of opened ion channels:

$$\begin{array}{*{20}l} \frac{I_{\text{antag}}}{I_{0}} = {\frac{[\text{ARA}^{\ast}]_{\text{peak}}}{[\text{ARA}^{\ast}]_{\text{max}}}}. \end{array} $$
(5)

Theoretical analysis of the competition model

To conduct the subsequent numerical analysis in a systematic way, here we theoretically analyze the model (3). Specifically, we provide a simplified representation of the model under the assumption that the dissociation rate constants kdissD1 and kdissD2 of an NDNB are much smaller than other rate constants such as kdissA1,kdissA2 and kdecay. For cisatracurium, the dissociation rate constants are reported in [5] as 13s−1 and 34s−1 for mouse adult and embryonic AChRs, respectively. Also, for (+)-tubocurarine and pancuronium, these values are reported in [4] as 2.1s−1 and 5.9s−1, respectively. Whereas, the dissociation rate constants of ACh and the rate constant of hydrolysis are estimated as 18000s−1 and 12000s−1, respectively [13, 15].

The simplified model is given in a dimensionless representation. For example, a dimensionless time τ and a concentration of ACh, xA, can be defined as

$$\begin{array}{*{20}l} \tau := t \cdot {k_{\text{decay}}}, \quad x_{\mathrm{A}} := \frac{\mathrm{[\!A]}}{ K_{\mathrm{A1}}}, \end{array} $$
(6)

where KA1 stands for the dissociation equilibrium constant for ACh with the site #1 given by kdissA1/kassocA1. Dimensionless variables \(x^{\ast }_{\text {ARA}}\) and xDRD for the concentrations of the complex ARA and DRD are given by

$$\begin{array}{*{20}l} {x^{\ast}_{\text{ARA}} := \frac{{[\text{ARA}^{\ast}]}}{\mathrm{[\!R]}_{\text{total}}}}, \quad x_{\text{DRD}} := \frac{\mathrm{[\!DRD]}}{\mathrm{[\!R]}_{\text{total}} }. \end{array} $$
(7)

Similarly, the dimensionless variables xARA,xARD,xARD,xDRA,xARO,xORA,xDRO, and xORD for the remaining seven complexes can be defined by dividing by [ R]total. By using these dimensionless variables, the simplified model is given as follows (see Appendix for its derivation):

$$\begin{array}{*{20}l} {}&\frac{\mathrm{d}}{\mathrm{d} \tau}x_{\mathrm{A}} = - x_{\mathrm{A}} +\kappa_{\mathrm{A1}}\lambda_{\mathrm{A1}} \{ x_{\text{ARA}} - x_{\mathrm{A}} (O_{\text{ORD}}(\delta) - x_{\text{ARD}}) \} \\ {}& \hspace{16mm} +\kappa_{\mathrm{A2}}\lambda_{\mathrm{A2}} \{ x_{\text{ARA}} - x_{\mathrm{A}} (O_{\text{DRO}}(\delta) - x_{\text{DRA}}) \} \\ {}& \hspace{16mm} +\kappa_{\mathrm{A1}}\lambda_{\mathrm{A1}} \{ x_{\text{ARD}} + x_{\text{ARO}} - x_{\mathrm{A}} (1-x_{\text{ARA}}-x_{\text{ARO}} - O_{\text{total}}(\delta)) \} \\ {}& \hspace{16mm} +\kappa_{\mathrm{A2}}\lambda_{\mathrm{A2}} \{ x_{\text{DRA}}+ x_{\text{ORA}} - x_{\mathrm{A}} (1-x_{\text{ARA}} -x_{\text{ORA}}- O_{\text{total}}(\delta)) \}, \end{array} $$
(8a)
$$\begin{array}{*{20}l} {}& {\frac{\mathrm{d}}{\mathrm{d} \tau}x^{\ast}_{\text{ARA}} = - \kappa_{\text{close}} x^{\ast}_{\text{ARA}} + \kappa_{\text{open}}x_{\text{ARA}}, } \end{array} $$
(8b)
$$\begin{array}{*{20}l} & \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{ARA}} = \kappa_{\mathrm{A1}} \left(x_{\text{ORA}} x_{\mathrm{A}} - x_{\text{ARA}} \right) + \kappa_{\mathrm{A2}} \left(\mu_{\mathrm{A}} x_{\text{ARO}} x_{\mathrm{A}} - x_{\text{ARA}} \right) \\ & \hspace{16mm} {+ \kappa_{\text{close}} x^{\ast}_{\text{ARA}} - \kappa_{\text{open}}x_{\text{ARA}}}, \end{array} $$
(8c)
$$\begin{array}{*{20}l} & \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{ARO}} = \kappa_{\mathrm{A1}} \left\{ x_{\mathrm{A}} (1-x_{\text{ARA}}-x_{\text{ARO}}-x_{\text{ORA}}- O_{\text{total}}(\delta))- x_{\text{ARO}} \right\} \\ & \hspace{16mm} +\kappa_{\mathrm{A2}} \left(x_{\text{ARA}} - \mu_{\mathrm{A}} x_{\text{ARO}} x_{\mathrm{A}} \right), \end{array} $$
(8d)
$$\begin{array}{*{20}l} & \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{ORA}} = \kappa_{\mathrm{A2}} \left\{ \mu_{\mathrm{A}} x_{\mathrm{A}} (1-x_{\text{ARA}}-x_{\text{ARO}}-x_{\text{ORA}}- O_{\text{total}}(\delta)) - x_{\text{ORA}} \right\} \\ &\hspace{16mm} + \kappa_{\mathrm{A1}} \left(x_{\text{ARA}} - x_{\text{ORA}} x_{\mathrm{A}} \right), \end{array} $$
(8e)
$$\begin{array}{*{20}l} & \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{ARD}} = \kappa_{\mathrm{A1}} \left\{ x_{\mathrm{A}} (O_{\text{ORD}}(\delta) - x_{\text{ARD}}) - x_{\text{ARD}} \right\}, \end{array} $$
(8f)
$$\begin{array}{*{20}l} & \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{DRA}}= \kappa_{\mathrm{A2}} \left\{ \mu_{\mathrm{A}} x_{\mathrm{A}} (O_{\text{DRO}}(\delta) - x_{\text{DRA}}) - x_{\text{DRA}} \right\}, \end{array} $$
(8g)

where κA1,κA2, κopen and κclose are the dimensionless parameters representing normalized rate constants given by

$$\begin{array}{*{20}l} & \kappa_{\mathrm{A1}} := \frac{k_{\mathrm{dissA1}} }{ {k_{\text{decay}}} }, \quad \kappa_{\mathrm{A2}} := \frac{k_{\mathrm{dissA2}} }{ {k_{\text{decay}}} }, \\ & {\kappa_{\text{open}} := \frac{k_{\text{open}} }{ {k_{\text{decay}}} }, \quad \kappa_{\text{close}} := \frac{k_{\text{close}} }{ {k_{\text{decay}}} }, } \end{array} $$
(9)

The parameters λA1,λA2, and μAare the dimensionless parameters representing affinities of ACh to the binding sites of the AChR given by

$$\begin{array}{*{20}l} \lambda_{\mathrm{A1}} := \frac{\mathrm{[\!R]}_{\text{total}}}{K_{\mathrm{A1}} }, \quad \lambda_{\mathrm{A2}} := \frac{\mathrm{[\!R]}_{\text{total}}}{K_{\mathrm{A2}} }, \quad \mu_{\mathrm{A}} := \frac{K_{\mathrm{A1}}}{K_{\mathrm{A2}}}, \quad \end{array} $$
(10)

with KA2:=kdissA2/kassocA2, and finally, μD and δ are the dimensionless parameters representing the site-selectivity and concentration of NDNB, respectively, given by

$$\begin{array}{*{20}l} \mu_{\mathrm{D}} := \frac{K_{\mathrm{D1}}}{K_{\mathrm{D2}}}, \quad \delta := \frac{{[\!\mathrm{D}]}}{ K_{\mathrm{D1}} }, \end{array} $$
(11)

with KD1:=kdissD1/kassocD1 and KD2:=kdissD2/kassocD2. Furthermore, ODRD,OORD and ODRO are functions of δ standing for the fractions of the complex DRD, ORD, and DRO before the release of ACh given by

$$\begin{array}{*{20}l} O_{\text{DRD}}(\delta) &:= \frac{ \mu_{\mathrm{D}} \delta^{2} }{ (1 + \delta)(1 + \mu_{\mathrm{D}} \delta) }, \end{array} $$
(12a)
$$\begin{array}{*{20}l} O_{\text{ORD}}(\delta) &:= \frac{ \mu_{\mathrm{D}} \delta }{ (1 + \delta)(1 + \mu_{\mathrm{D}} \delta) }, \end{array} $$
(12b)
$$\begin{array}{*{20}l} O_{\text{DRO}}(\delta) &:= \frac{ \delta }{ (1 + \delta)(1 + \mu_{\mathrm{D}} \delta)}. \end{array} $$
(12c)

The total occupancy Ototal of the AChRs by NDNB is given by

$$\begin{array}{*{20}l} O_{\text{total}}(\delta) & := O_{\text{DRD}}(\delta) + O_{\text{ORD}}(\delta) + O_{\text{DRO}}(\delta). \end{array} $$
(13a)

From the simplified model (8), several insights can be obtained on the relationship between the receptor occupancy by NDNBs and the dynamics of activation of AChRs. First, it can be seen that the model (8) has only one parameter μD characterizing NDNBs. Thus, the properties of NDNBs are completely determined by the parameter μD representing the site-selectivity of NDNBs. This fact implies that even in the simulations of the original model (3), the results will be almost unchanged if the parameter μD kept constant and the parameters κD1:=kdissD1/kdecay and κD2:=kdissD2/kdecay are small enough to validate the simplification (see Appendix for details of the validity of the simplification). This observation facilitates the numerical analysis of the original model (3) by varying the parameters of the model in a systematic way as demonstrated in the rest of the paper. Second, from the model Eq. 8, it can be seen that the entire states during competition can be viewed as divided into 4 parts as shown in Fig. 3. Note that owing to an appropriate derivation of the dimensionless representation of the model, not only the terms for dissociation of NDNBs but also association terms are eliminated from the original model. As a result, the state DRD and the pairs of states (ORD, ARD) and (DRO, DRA) can be viewed as separated from the states describing the activation of AChRs, while interaction between these states may occur through the dynamics of free ACh given by (8a). Specifically, although the fraction of the state DRD will not change during the competition, the fractions of ORD and DRO will change due to association and dissociation with ACh molecules. Therefore, even if the total occupancy Ototal is the same, the fraction of the activated AChRs will decrease with the increase in the occupancies OORD and ODRO because free ACh will decrease due to association with ORD and DRO. Note that the ratio of ORD and DRO to Ototal is determined by the parameter μD. Finally, the difference between the receptor binding model (2) and the competition-based models (3) and (8) can be clarified. By using the notation introduced above, the model (2) can be rewritten as follows:

$$\begin{array}{*{20}l} \frac{I_{\text{antag}}}{I_{0}} = \frac{1}{1 + \delta + \mu_{\mathrm{D}} \delta + \mu_{\mathrm{D}} \delta^{2}} = 1 - O_{\text{total}}(\delta) \end{array} $$
(14)
Fig. 3
figure3

Diagram of the interactions between acetylcholine (A), NDNB (D), and the postsynaptic receptor (R) based on the simplified model (8). The state DRD and the pairs of states (ORD, ARD) and (DRO, DRA) can be viewed as separated from the states describing the activation of AChRs, while interaction between these states may occur through the dynamics of free ACh given by (8a)

Thus, the model (2) is clearly based on the assumption that the inhibition is proportional to the total occupancy Ototal, while the models (3) and (8) describe the effect of the partial occupancies OORD and ODRO on the activation of AChRs during competition.

Method of numerical analysis on the relationship between inhibition and receptor occupancy

The relationship between the inhibition and the receptor occupancy by an NDNB is studied via numerical simulations of the original model (3). To visualize the results, the fraction of activated AChRs given by [ARA]peak/[ARA]max is calculated under various concentrations of the NDNB (100 values between [ D]=1.0×10−14M and 1.0×10−5M that are spaced equidistantly on a logarithmic scale, and [ D]=0) and plotted versus the total occupancy Ototal. For calculating the peak concentration of ARA, the ordinary differential equation model (3) is numerically solved using the Fortran solver LSODA provided by the python package SciPy (Version 1.5.2).

The parameters of the model (3) are varied in a systematic way based on information provided by literature and the findings of the theoretical analysis as explained in the following. First, we investigate the effect of varying the kinetic constants for ACh and NDNB. For ACh, the values of kdissA1,kdissA2,kassocA1, and kassocA2 presented in Table 1 were determined by experiments using mouse adult AChRs. However, it is known that the EC50, the concentration of ACh where the half of the AChRs are activated, for human adult AChRs is smaller than that for mouse adult AChRs: EC50=8.48×10−6M or 7.91×10−6M for human adult AChRs [10] and EC50=1.6×10−5M for mouse adult AChRs [13]. Since this implies that the affinity of ACh is different between human and mouse AChRs, we investigated the effect of varying the constants KA1 and KA2 by changing kassocA1 and kassocA2 while kdissA1 and kdissA2 kept constant. Furthermore, since it has been reported in [13] that the two binding sites of mouse adult AChR have similar affinities, we assume that it is also the case for human adult AChR and the parameters KA1 and KA2 are equal. Thus, we investigate the effect of varying the values of KA1=KA2=KA, which corresponds to κA1=κA2 and μA=1. For the kinetic constants of an NDNB, the effect of varying the parameter μD is investigated, since this parameter completely characterizes the properties of the NDNB as far as the dissociation constants kdissD1 and kdissD2 are small enough. The value of parameter μD was assigned by changing the value of kassocD2, while kdissD1,kdissD2 and kassocD1 kept constant at the values shown in Table 1. Also, the values of the dissociation constants kdissD1 and kdissD2 are varied to explore the difference between the original model (3) and the reduced-order model (8).

Furthermore, we examine the effect of varying the initial concentration [ A]init representing the concentration of ACh immediately after the release of ACh at t=0. It has been known that the number of ACh molecules released in vivo is one tenth of the number of AChRs at the synaptic cleft [15]. In this paper, following [15], the concentrations of AChRs and the initial ACh were set as [ R]total=7.75×10−5M and [ A]init=7.75×10−6M, respectively. Under this setting, only a small fraction (at least less than one tenth) of the total receptor population will be activated in vivo. However, in many in vitro experiments, a quite high concentration of ACh is used [48] such that nearly 50% (nearly EC50) of the AChRs are activated or more than 90% of AChRs are activated. Thus, we investigated the effect of varying [ A]init in a wide range of concentrations from [ A]init=0.0562×[ R]total to 5.62×[ R]total.

Results

The upper panel of Fig. 4 shows the relationship between the fraction of activated AChRs given by [ARA]peak/[ARA]max and the receptor occupancy Ototal evaluated under various settings of the dissociation equilibrium constants for ACh. The values of KA1=KA2=KA is one of 1.0×10−4M,3.16×10−5M and 1.0×10−5M, and shown by red, green and blue lines, respectively. The dotted line in the figure shows the linear relationship corresponding to the model (2). The results clearly show that there are nonlinear relationships between [ARA]peak/[ARA]max and Ototal in all the cases. Furthermore, it can be confirmed that the extent of nonlinearity increases with the decrease in the equilibrium constant KA. In the subsequent analyses, to highlight the effects of varying other parameters, we assigned the value of 1.0×10−5M to KA, where the extent of nonlinearity was most prominent. Then, the lower panel of Fig. 4 shows the results under various settings of μD representing the site-selectivity of an NDNB. The red, green and blue lines show the results for μD=1, 10, and 100, respectively. It can be seen that the relationships between [ARA]peak/[ARA]max and Ototal are nonlinear in all the cases, and the extent of nonlinearity decreases with the increase in μD.

Fig. 4
figure4

Effect of varying kinetic constants for ACh and NDNB on the relationship between the fraction of activated AChRs, [ARA]peak/[ARA]max, and the receptor occupancy Ototal. a Results of varying the equilibrium dissociation constants for ACh under KA=KA1=KA2 (κA1=κA2 and μA=1). b Results of varying the site-selectivity μD=KD1/KD2 of an NDNB

Next, the upper and the middle panels of Fig. 5 show the relationship between the fraction of activated AChRs, [ARA]peak/[ARA]max, and the receptor occupancy Ototal under various settings of the initial concentration of ACh, [A]init. In the upper panel, the initial concentration [A]init is decreased from 5.62×[R]total to 0.316×[R]total, and it can be seen that the extent of nonlinearity increases with the decrease in [A]init. Also, in the middle panel, the initial concentration is further decreased to 0.0562×[R]total. The extent of nonlinearlity slightly decreases in this range of parameter setting, and thus the nonlinearity is most prominent at [ A]init=0.316×[R]total. Interestingly, the relationship between [ARA]peak/[ARA]max and Ototal becomes almost linear when the value of [A]init is larger than 3.16×[R]total. To clarify the meaning of this result, the lower panel of Fig. 5 shows the concentration-response relationship between the concentration [A]init and the activated AChRs [ARA]max. Note that the [ARA]max is defined for each setting of [ A]init as the highest [ARA] under various settings of the concentration of an NDNB, and it is attained in the absence of NDNB. At the setting of [ A]init=0.1×[R]total, which corresponds to the in vivo situation [15], only a fraction of AChRs are activated even in the absence of NDNBs, and result in a highly nonlinear relationship between [ARA]peak/[ARA]max and Ototal. However, at the setting of [ A]init≥3.16×[R]total, more than 92% of AChRs are activated in the absence of NDNB, and in this case, the relationship between [ARA]peak/[ARA]max and Ototal becomes linear.

Fig. 5
figure5

Effect of varying the initial concentration [ A]init of ACh on the numerical results. a The relationship between the fraction of activated AChRs, [ARA]peak/[ARA]max, and the receptor occupancy Ototal under various settings of [ A]init. b The concentration-effect relationship between the activated AChRs in the absence of NDNBs and the initial concentration [ A]init

Finally, the effect of varying the dissociation rate constants kdissD1 and kdissD2 of NDNB was examined under the setting of identical rate constants, i.e. kdissD1=kdissD2=kdissD. The upper panel of Fig. 6 shows the relationship between the fraction of activated AChRs, [ARA]peak/[ARA]max, and the receptor occupancy Ototal under the setting of [ A]init=0.1×[R]total corresponding to in vivo concentration. With this setting, the results are almost unchanged even when the rate constant kdissD is as large as kdissA1,kdissA2, and kdecay. This implies that the model simplification performed in this paper can be validated for a wide range of settings of kdissD. The lower panel of Fig. 6 shows the results under the setting of [ A]init=10×[R]total corresponding to in vitro concentration. In this case, it can be seen that the results are highly affected by the setting of the dissociation rate constant kdissD. Furthermore, the relationships between [ARA]peak/[ARA]max and Ototal are no longer linear when the rate constant kdissD becomes large even if the concentration [ A]init is large. This implies that the model simplification can be validated only for small values of the parameter kdissD.

Fig. 6
figure6

Effect of varying the dissociation rate constants kdissD1 and kkdissD2 on the relationship between the fraction of activated AChRs, [ARA]peak/[ARA]max, and the receptor occupancy Ototal. a Results under the setting of [ A]init=0.1×[R]total corresponding to in vivo concentration. b Results under the setting of [ A]init=10×[R]total corresponding to in vitro concentration

Discussion

We theoretically and numerically investigated the relationship between inhibition and receptor occupancy by NDNBs. While the two-site receptor binding model (2), which assumes a linear relationship between the inhibition and the receptor occupancy by an NDNB, has been statistically tested in [48] for several in vitro experimental settings, it has not been studied in literature under which conditions the above assumption holds nor if the assumption remains valid in vivo. To consider these problems, an ordinary differential equation model was introduced based on [5, 1315] to simulate the physiologic processes of activation of receptors by ACh as well as inhibition by an NDNB. The theoretical analysis performed in this paper clarified that under the assumption that the dissociation rate constants for an NDNB are small and with an appropriate nondimensionalization, the characteristics of an NDNB can be completely determined by a single parameter μD representing the site-selectivity of the NDNB for two binding sites of AChRs. We then performed numerical analysis of the model by plotting the fractional amounts of the activated AChRs as a function of the receptor occupancy. The numerical results show that under a setting of parameters reflecting in vivo environment, there is a nonlinear relationship between the inhibition and the receptor occupancy, indicating limitation of the applicability of the receptor binding model. Furthermore, it has been shown that the extent of nonlinearity depends on the parameters representing kinetic constants for ACh or NDNBs and the concentration of ACh.

Regarding the nonlinear relationship between the effect and the receptor occupancy by an NDNB, it has been well known that the twitch strength (the degree of muscle contraction) observed in vivo is not proportional to the receptor occupancy due to the high margin of safety [21]. The origin of the safety margin is a copious density of AChRs on the post-synaptic membrane and the fact that only a small fraction of AChRs needs to be activated to trigger the occurrence of an action potential and the contraction of the associated muscle fiber. Thus, the nonlinearity due to the safety margin means that the response of muscle fiber is not proportional to the fraction of activated AChRs, and the extent of nonlinearity would not be affected by the properties of an NDNB. However, this paper focused on the nonlinearity in the relationship between the receptor occupancy and the fraction of activated AChRs, and the extent of nonlinearity is affected by the properties of an NDNB. In particular, it has been revealed in this paper that the effect of an NDNB on the extent of nonlinearity can be characterized by a single parameter μD representing the site-selectivity of the NDNB.

Although the model used and simulations performed in this paper are intended to describe in vivo situations, our finding that the extent of nonlinearity is affected by the concentration of ACh is consistent with results observed through in vitro experiments. In [9, 10], it was found that the IC50 for several clinically used NDNBs, cis-atracrium, rocuronium, and vecronium, decreased with the increase in the concentration of ACh. With these results, it was concluded in [9, 10, 22] that this phenomenon indicates a noncompetitive component of inhibition under the idea that the enhancement of the inhibition was caused by a mechanism different from competitive one occurred by NDNBs in combination with high concentration of ACh. However, such a shift in the values of IC50 can also be explained by a change in the relationship between the receptor occupancy and the fraction of activated AChRs. As demonstrated in the upper panel of Fig. 5, the receptor occupancy Ototal at which [ARA]peak/[ARA]max (which corresponds to the relative current) takes the value of 0.5 increases with the decrease in the concentration of ACh, meaning that the IC50 increases with the decrease in the concentration of ACh. This shift in the IC50 is consistent with the observation in [9, 10] and occurs under a totally competitive mechanism of inhibition by an NDNB.

Interestingly, it was found that the relationship between the fraction of activated AChRs and the receptor occupancy became linear when the concentration of ACh was sufficiently large and the dissociation rate constants of an NDNB were sufficiently small. This finding may provide a reasonable justification of the use of the two-site model (2) in the analysis of kinetic constants for NDNBs through in vitro experiments. At least, the condition of large concentration of ACh is satisfied in the experiments reported in [48], where concentrations that opens about 93% to 95% of the AChRs were used. However, further consideration is needed to identify the conditions needed for the justification of the model (2), because the present study did not take into account the effect of desensitization of AChRs, which is the main cause of the decay in a measured current observed during in vitro experiments.

Conclusion

The relationship between the inhibition and the receptor occupancy by an NDNB was theoretically and numerically investigated. While a receptor binding model, which assumes a linear relationship, may be effective for estimating affinity of an NDNB through in vitro experiments, the model do not directly describe in vivo pharmacologic properties of NDNBs, because the nonlinearity between the inhibition and the receptor occupancy causes the modulation of the resultant concentration-effect relationships of NDNBs. It was found that the effect of characteristics of an NDNB on the extent of nonlinearity can be identified by a single parameter representing the site-selectivity of an NDNB.

Appendix

This appendix provides a detailed derivation of the simplified model (8). Specifically, we derive a dimensionless representation of the original model based on the technique of scaling [23] and perform model-order reduction based on singular perturbation theory [24, 25]. The scaling of the model (3) can be done by introducing dimensionless variables. For example, a dimensionless time τ and a concentration of ACh, xA, can be defined as

$$\begin{array}{*{20}l} \tau := t \cdot {k_{\text{decay}} }, \quad x_{\mathrm{A}} := \frac{\mathrm{[\!A]}}{ K_{\mathrm{A1}} }, \end{array} $$
(15)

where KA1 stands for the dissociation equilibrium constant for ACh with the site #1 given by kdissA1/kassocA1. Dimensionless variables \(x^{\ast }_{\text {ARA}}\) and xDRD for the concentrations of the complex ARA and DRD are given by

$$\begin{array}{*{20}l} {x^{\ast}_{\text{ARA}} := \frac{{[\text{ARA}^{\ast}]} }{ \mathrm{[\!R]}_{\text{total}} }}, \quad x_{\text{DRD}} := \frac{\mathrm{[\!DRD]} }{ \mathrm{[\!R]}_{\text{total}} }. \end{array} $$
(16)

Similarly, the dimensionless variables xARA, xARD,xARD,xDRA,xARO,xORA,xDRO, and xORD for the remaining seven complexes can be defined by dividing by [ R]total. By substituting these dimensionless variables to the model (3), the following equations can be derived:

$$\begin{array}{*{20}l} {}& \frac{\mathrm{d}}{\mathrm{d} \tau} x_{\mathrm{A}} = - x_{\mathrm{A}} + \kappa_{\mathrm{A1}} \lambda_{\mathrm{A1}} \{ x_{\text{ARA}} + x_{\text{ARD}} + x_{\text{ARO}} \\ {}& \hspace{13mm} - x_{\mathrm{A}} (x_{\text{ORA}} + x_{\text{ORD}} +x_{\text{ORO}}) \} \\ {}&\hspace{13mm} + \kappa_{\mathrm{A2}} \{\lambda_{\mathrm{A1}}(x_{\text{ARA}} + x_{\text{DRA}}+ x_{\text{ORA}}) \\ {}&\hspace{13mm} - \lambda_{\mathrm{A2}} x_{\mathrm{A}} (x_{\text{ARO}} + x_{\text{DRO}}+ x_{\text{ORO}}) \}, \end{array} $$
(17a)
$$\begin{array}{*{20}l} {}& {\frac{\mathrm{d}}{\mathrm{d} \tau}x^{\ast}_{\text{ARA}} = - \kappa_{\text{close}} x^{\ast}_{\text{ARA}} + \kappa_{\text{open}}x_{\text{ARA}},} \end{array} $$
(17b)
$$\begin{array}{*{20}l} {}& \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{ARA}} = \kappa_{\mathrm{A1}} \left(x_{\text{ORA}} x_{\mathrm{A}} - x_{\text{ARA}} \right) + \kappa_{\mathrm{A2}} \left(\mu_{\mathrm{A}} x_{\text{ARO}} x_{\mathrm{A}} - x_{\text{ARA}} \right) \\ {}&\hspace{16mm} {+ \kappa_{\text{close}} x^{\ast}_{\text{ARA}} - \kappa_{\text{open}}x_{\text{ARA}},} \end{array} $$
(17c)
$$\begin{array}{*{20}l} {}& \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{DRD}} = \kappa_{\mathrm{D1}} \left(x_{\text{ORD}} \delta - x_{\text{DRD}} \right) + \kappa_{\mathrm{D2}} \left(\mu_{\mathrm{D}}x_{\text{ORD}} \delta - x_{\text{DRD}} \right), \end{array} $$
(17d)
$$\begin{array}{*{20}l} {}& \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{ARD}} = \kappa_{\mathrm{A1}} \left(x_{\text{ORD}} x_{\mathrm{A}} - x_{\text{ARD}} \right) + \kappa_{\mathrm{D2}} \left(\mu_{\mathrm{D}} x_{\text{ARO}} \delta - x_{\text{ARD}} \right), \end{array} $$
(17e)
$$\begin{array}{*{20}l} {}& \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{DRA}} = \kappa_{\mathrm{D1}} \left(x_{\text{ORA}} \delta - x_{\text{DRA}} \right) + \kappa_{\mathrm{A2}} \left(\mu_{\mathrm{A}} x_{\text{DRO}} x_{\mathrm{A}} - x_{\text{DRA}} \right), \end{array} $$
(17f)
$$\begin{array}{*{20}l} {}&\frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{ARO}} = \kappa_{\mathrm{A1}} \left(x_{\text{ORO}} x_{\mathrm{A}} - x_{\text{ARO}} \right) + \kappa_{\mathrm{A2}} \left(x_{\text{ARA}} - \mu_{\mathrm{A}} x_{\text{ARO}} x_{\mathrm{A}} \right) \\ {}& \hspace{16mm} + \kappa_{\mathrm{D2}}\left(x_{\text{ARD}} - \mu_{\mathrm{D}}x_{\text{ARO}} \delta \right), \end{array} $$
(17g)
$$\begin{array}{*{20}l} {}& \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{ORA}} = \kappa_{\mathrm{A2}} \left(\mu_{\mathrm{A}} x_{\text{ORO}} x_{\mathrm{A}} - x_{\text{ORA}} \right) + \kappa_{\mathrm{A1}} \left(x_{\text{ARA}} - x_{\text{ORA}} x_{\mathrm{A}} \right) \\ {}& \hspace{16mm} + \kappa_{\mathrm{D1}} \left(x_{\text{DRA}} - x_{\text{ORA}} \delta \right), \end{array} $$
(17h)
$$\begin{array}{*{20}l} {}& \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{DRO}} = \kappa_{\mathrm{D1}} \left(x_{\text{ORO}} \delta - x_{\text{DRO}} \right) + \kappa_{\mathrm{A2}} \left(x_{\text{DRA}} - \mu_{\mathrm{A}} x_{\text{DRO}} x_{\mathrm{A}} \right) \\ {}& \hspace{16mm} + \kappa_{\mathrm{D2}}\left(x_{\text{DRD}} - \mu_{\mathrm{D}} x_{\text{DRO}} \delta \right), \end{array} $$
(17i)
$$\begin{array}{*{20}l} {}& \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{ORD}} = \kappa_{\mathrm{D2}} \left(\mu_{\mathrm{D}} x_{\text{ORO}} \delta - x_{\text{ORD}} \right) + \kappa_{\mathrm{A1}} \left(x_{\text{ARD}} - x_{\text{ORD}} x_{\mathrm{A}} \right) \\ {}& \hspace{16mm} + \kappa_{\mathrm{D1}} \left(x_{\text{DRD}} - x_{\text{ORD}} \delta \right), \end{array} $$
(17j)

where κA1,κA2,κD1 and κD2 are the dimensionless parameters representing the rates of the both dissociation and association of the ligands given by

$$\begin{array}{*{20}l} \kappa_{\mathrm{A1}} := \frac{k_{\mathrm{dissA1}} }{ {k_{\text{decay}}} }, \quad \kappa_{\mathrm{A2}} := \frac{k_{\mathrm{dissA2}} }{ {k_{\text{decay}}} }, \\ \kappa_{\mathrm{D1}} := \frac{k_{\mathrm{dissD1}} }{ {k_{\text{decay}}} }, \quad \kappa_{\mathrm{D2}} := \frac{k_{\mathrm{dissD2}} }{ {k_{\text{decay}}} }, \end{array} $$
(18)

and λA1,λA2, and μAare the dimensionless parameters representing affinities of ACh to the binding sites of the AChR given by

$$\begin{array}{*{20}l} \lambda_{\mathrm{A1}} := \frac{\mathrm{[\!R]}_{\text{total}}}{K_{\mathrm{A1}} }, \quad \lambda_{\mathrm{A2}} := \frac{\mathrm{[\!R]}_{\text{total}}}{K_{\mathrm{A2}} }, \quad \mu_{\mathrm{A}} := \frac{K_{\mathrm{A1}}}{K_{\mathrm{A2}}}, \quad \end{array} $$
(19)

with KA2:=kdissA2/kassocA2, and finally, μD and δ are the dimensionless parameters representing the site-selectivity and concentration of NDNB, respectively, given by

$$\begin{array}{*{20}l} \mu_{\mathrm{D}} := \frac{K_{\mathrm{D1}}}{K_{\mathrm{D2}}}, \quad \delta := \frac{{[\!\mathrm{D}]}}{ K_{\mathrm{D1}} }, \end{array} $$
(20)

with KD1:=kdissD1/kassocD1 and KD2:=kdissD2/kassocD2. Note that due to the scaling performed here, the number of parameters characterizing the properties of NDNB can be reduced from four to three: from (kdissD1,kasoocD1,kdissD2,kassocD2) to (κD1,κD2,μD).

Furthermore, the order of the model can be reduced using the technique of singular perturbation [24, 25] based on an inherent multiple-timescale property of the model. Such a multi-scale property arises when the dissociation rate constants kdissD1 and kdissD2 of an NDNB are much smaller than other rate constants such as kdissA1,kdissA2 and kdecay. By considering the limit κD1, κD2→0, the following equations hold from the Eq. 17:

$$\begin{array}{*{20}l} & \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{DRD}} = \frac{\mathrm{d}}{\mathrm{d} \tau}(x_{\text{ORD}}+x_{\text{ARD}}) = \frac{\mathrm{d}}{\mathrm{d} \tau}(x_{\text{DRO}}+ x_{\text{DRA}}) = 0. \end{array} $$
(21)

Thus, the values of xDRD,xORD+xARD and xDRO+xDRA do not change in the reduced-order model, or almost unchanged in the original model (3), from their initial values at τ=0. When the initial conditions are defined by the steady state concentrations under a given value of δ, the initial values of xDRD,xORD,xDRO,xARD and xDRA are given by

$$\begin{array}{*{20}l} x_{\text{DRD}}|_{\tau=0} &= O_{\text{DRD}}(\delta) := \frac{ \mu_{\mathrm{D}} \delta^{2} }{(1 + \delta)(1 + \mu_{\mathrm{D}} \delta) }, \end{array} $$
(22a)
$$\begin{array}{*{20}l} x_{\text{ORD}}|_{\tau=0} &= O_{\text{ORD}}(\delta) := \frac{\mu_{\mathrm{D}} \delta }{ (1 + \delta)(1 + \mu_{\mathrm{D}} \delta)}, \end{array} $$
(22b)
$$\begin{array}{*{20}l} x_{\text{DRO}}|_{\tau=0} &= O_{\text{DRO}}(\delta) := \frac{ \delta }{ (1 + \delta)(1 + \mu_{\mathrm{D}} \delta)}, \end{array} $$
(22c)
$$\begin{array}{*{20}l} x_{\text{ARD}}|_{\tau=0} & = x_{\text{DRA}}|_{\tau=0} = 0 \end{array} $$
(22d)

where ODRD,OORD and ODRO stand for the fractions of the complex DRD, ORD, and DRO, respectively, at the steady state. Also, the total occupancy Ototal of the AChRs by the NDNB is given by

$$\begin{array}{*{20}l} O_{\text{total}}(\delta) := O_{\text{DRD}}(\delta) + O_{\text{ORD}}(\delta) + O_{\text{DRO}}(\delta). \end{array} $$
(23)

By using the equations from (21) to (23), the model (3) can be reduced to the following form:

$$\begin{array}{*{20}l} {}&\frac{\mathrm{d}}{\mathrm{d} \tau}x_{\mathrm{A}} = - x_{\mathrm{A}} +\kappa_{\mathrm{A1}}\lambda_{\mathrm{A1}} \{ x_{\text{ARA}} - x_{\mathrm{A}} (O_{\text{ORD}}(\delta) - x_{\text{ARD}}) \} \\ {}& \hspace{16mm} +\kappa_{\mathrm{A2}}\lambda_{\mathrm{A2}} \{ x_{\text{ARA}} - x_{\mathrm{A}} (O_{\text{DRO}}(\delta) - x_{\text{DRA}}) \} \\ {}& \hspace{16mm} +\kappa_{\mathrm{A1}}\lambda_{\mathrm{A1}} \{ x_{\text{ARD}} + x_{\text{ARO}} - x_{\mathrm{A}} (1-x_{\text{ARA}}-x_{\text{ARO}} - O_{\text{total}}(\delta)) \} \\ {}& \hspace{16mm} +\kappa_{\mathrm{A2}}\lambda_{\mathrm{A2}} \{ x_{\text{DRA}}+ x_{\text{ORA}} - x_{\mathrm{A}} (1-x_{\text{ARA}} -x_{\text{ORA}}- O_{\text{total}}(\delta)) \}, \end{array} $$
(24a)
$$\begin{array}{*{20}l} {}& {\frac{\mathrm{d}}{\mathrm{d} \tau}x^{\ast}_{\text{ARA}} = - \kappa_{\text{close}} x^{\ast}_{\text{ARA}} + \kappa_{\text{open}}x_{\text{ARA}},} \end{array} $$
(24b)
$$\begin{array}{*{20}l} {}& \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{ARA}} = \kappa_{\mathrm{A1}} \left(x_{\text{ORA}} x_{\mathrm{A}} - x_{\text{ARA}} \right) + \kappa_{\mathrm{A2}} \left(\mu_{\mathrm{A}} x_{\text{ARO}} x_{\mathrm{A}} - x_{\text{ARA}} \right) \\ {}& \hspace{16mm} {+ \kappa_{\text{close}} x^{\ast}_{\text{ARA}} - \kappa_{\text{open}}x_{\text{ARA}}}, \end{array} $$
(24c)
$$\begin{array}{*{20}l} {}& \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{ARO}} = \kappa_{\mathrm{A1}} \left\{ x_{\mathrm{A}} (1-x_{\text{ARA}}-x_{\text{ARO}}-x_{\text{ORA}}- O_{\text{total}}(\delta))- x_{\text{ARO}} \right\} \\ {}& \hspace{16mm} +\kappa_{\mathrm{A2}} \left(x_{\text{ARA}} - \mu_{\mathrm{A}} x_{\text{ARO}} x_{\mathrm{A}} \right), \end{array} $$
(24d)
$$\begin{array}{*{20}l} {}& \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{ORA}} = \kappa_{\mathrm{A2}} \left\{ \mu_{\mathrm{A}} x_{\mathrm{A}} (1-x_{\text{ARA}}-x_{\text{ARO}}-x_{\text{ORA}}- O_{\text{total}}(\delta)) - x_{\text{ORA}} \right\} \\ {}&\hspace{16mm} + \kappa_{\mathrm{A1}} \left(x_{\text{ARA}} - x_{\text{ORA}} x_{\mathrm{A}} \right), \end{array} $$
(24e)
$$\begin{array}{*{20}l} {}& \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{ARD}} = \kappa_{\mathrm{A1}} \left\{ x_{\mathrm{A}} (O_{\text{ORD}}(\delta) - x_{\text{ARD}}) - x_{\text{ARD}} \right\}, \end{array} $$
(24f)
$$\begin{array}{*{20}l} {}& \frac{\mathrm{d}}{\mathrm{d} \tau}x_{\text{DRA}}= \kappa_{\mathrm{A2}} \left\{ \mu_{\mathrm{A}} x_{\mathrm{A}} (O_{\text{DRO}}(\delta) - x_{\text{DRA}}) - x_{\text{DRA}} \right\}. \end{array} $$
(24g)

Availability of data and materials

The data generated and/or analyzed in this research can be reproduced using numerical simulations explained in Methods section.

Abbreviations

NDNB:

Nondepolarizing neuromuscular blocking drugs

ACh:

Acetylcholine

AChR:

Acetylcholine receptor

References

  1. 1

    Pardo MC, Miller RD. Basics of Anesthesia, 7th edn. Philadelphia: Elsevier; 2018.

    Google Scholar 

  2. 2

    Yost CS, Winegar BD. Potency of agonists and competitive antagonists on adult- and fetal-type nicotinic acetylcholine receptors. Cell Mol Neurobiol. 1997; 17(1):35–50. https://doi.org/10.1023/A:1026325020191.

    Article  CAS  Google Scholar 

  3. 3

    Paul M, Kindler CH, Fokt RM, Dresser MJ, Dipp NCJ, Yost CS. The potency of new muscle relaxants on recombinant muscle-type acetylcholine receptors. Anesth Analg. 2002; 94(3):597–603. https://doi.org/10.1097/00000539-200203000-00022.

    Article  CAS  Google Scholar 

  4. 4

    Wenningmann I, Dilger JP. The kinetics of inhibition of nicotinic acetylcholine receptors by (+)-tubocurarine and pancuronium. Mol Pharmacol. 2001; 60(4):790–6. http://arxiv.org/abs/https://molpharm.aspetjournals.org/content/60/4/790.full.pdf.

    PubMed  CAS  Google Scholar 

  5. 5

    Demazumder D, Dilger JP. The kinetics of competitive antagonism by cisatracurium of embryonic and adult nicotinic acetylcholine receptors. Mol Pharmacol. 2001; 60(4):797–807. http://arxiv.org/abs/https://molpharm.aspetjournals.org/content/60/4/797.full.pdf.

    PubMed  CAS  Google Scholar 

  6. 6

    Demazumder D, Dilger JP. The kinetics of competitive antagonism of nicotinic acetylcholine receptors at physiological temperature. J Physiol. 2008; 586(4):951–63. https://doi.org/10.1113/jphysiol.2007.143289.

    Article  CAS  Google Scholar 

  7. 7

    Liu M, Dilger JP. Synergy between pairs of competitive antagonists at adult human muscle acetylcholine receptors. Anesth Analg. 2008; 107(2):525–33. https://doi.org/10.1213/ane.0b013e31817b4469.

    Article  CAS  Google Scholar 

  8. 8

    Liu M, Dilger JP. Site selectivity of competitive antagonists for the mouse adult muscle nicotinic acetylcholine receptor. Mol Pharmacol. 2009; 75(1):166–73. https://doi.org/10.1124/mol.108.051060.

    Article  CAS  Google Scholar 

  9. 9

    Jonsson M, Gurley D, Dabrowski M, Larsson O, Johnson EC, Eriksson LI. Distinct Pharmacologic Properties of Neuromuscular Blocking Agents on Human Neuronal Nicotinic Acetylcholine Receptors: A Possible Explanation for the Train-of-four Fade. Anesthesiology. 2006; 105(3):521–33. https://doi.org/10.1097/00000542-200609000-00016.

    Article  CAS  Google Scholar 

  10. 10

    Fagerlund MJ, Dabrowski M, Eriksson LI. Pharmacological characteristics of the inhibition of nondepolarizing neuromuscular blocking agents at human adult muscle nicotinic acetylcholine receptor. Anesthesiology. 2009; 110:1244–52. https://doi.org/10.1097/ALN.0b013e31819fade3.

    Article  CAS  Google Scholar 

  11. 11

    Dilger JP, Steinbach JH. Inhibition of muscle acetylcholine receptors by nondepolarizing drugs: Humans are not unique. Anesthesiology. 2010; 112:247–9. https://doi.org/10.1097/ALN.0b013e3181c5dbbc.

    Article  Google Scholar 

  12. 12

    Gabrielsson J, Weiner D. Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications, 5th edn. Stockholm: Swedish Pharmaceutical Press; 2016.

    Google Scholar 

  13. 13

    Akk G, Auerbach A. Inorganic, monovalent cations compete with agonists for the transmitter binding site of nicotinic acetylcholine receptors. Biophys J. 1996; 70(6):2652–8. https://doi.org/10.1016/S0006-3495(96)79834-X.

    Article  CAS  Google Scholar 

  14. 14

    Auerbach A, Akk G. Desensitization of mouse nicotinic acetylcholine receptor channels : A two-gate mechanism. J Gen Physiol. 1998; 112(2):181–97.

    Article  CAS  Google Scholar 

  15. 15

    Nigrovic V, Amann A. Competition between acetylcholine and a nondepolarizing muscle relaxant for binding to the postsynaptic receptors at the motor end plate: Simulation of twitch strength and neuromuscular block. J Pharmacokinet Pharmacodyn. 2003; 30(1):23–51. https://doi.org/10.1023/A:1023245409315.

    Article  CAS  Google Scholar 

  16. 16

    Fambrough DM, Drachman DB, Satyamurti S. Neuromuscular junction in myasthenia gravis: Decreased acetylcholine receptors. Science. 1973; 182(4109):293–5. https://doi.org/10.1126/science.182.4109.293.

    Article  CAS  Google Scholar 

  17. 17

    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(2):263–89. https://doi.org/10.1016/S0006-3495(79)85249-2.

    Article  CAS  Google Scholar 

  18. 18

    Salpeter MM, Eldefrawi ME. Sizes of end plate compartments, densities of acetylcholine receptor and other quantitative aspects of neuromuscular transmission. J Histochem Cytochem. 1973; 21(9):769–78. https://doi.org/10.1177/21.9.769.

    Article  CAS  Google Scholar 

  19. 19

    Hobbiger F. Neuromuscular Junction In: Zaimis E, Maclagan J, editors. Berlin: Springer: 1976. p. 487–581.

  20. 20

    Potter LT. Synthesis, storage and release of [14c]acetylcholine in isolated rat diaphragm muscles. J Physiol. 1970; 206(1):145–66. https://doi.org/10.1113/jphysiol.1970.sp009003.

    Article  CAS  Google Scholar 

  21. 21

    Paton WDM, Waud DR. The margin of safety of neuromuscular transmission. J Physiol. 1967; 191(1):59–90. https://doi.org/10.1113/jphysiol.1967.sp008237.

    Article  CAS  Google Scholar 

  22. 22

    Fagerlund MJ, Eriksson LI. Current concepts in neuromuscular transmission. Br J Anaesth. 2009; 103(1):108–14. https://doi.org/10.1093/bja/aep150. http://oup.prod.sis.lan/bja/article-pdf/103/1/108/17424609/aep150.pdf.

    Article  CAS  Google Scholar 

  23. 23

    Langtangen HP, Pedersen GK. Scaling of Differential Equations. Cham: Springer; 2016.

    Book  Google Scholar 

  24. 24

    Kokotović P, Khalil H, O’Reilly J. Singular Perturbation Methods in Control: Analysis and Design. Philadelphia: Society for Industrial and Applied Mathematics; 1999. https://doi.org/10.1137/1.9781611971118.

    Book  Google Scholar 

  25. 25

    Kuehn C. Multiple Time Scale Dynamics. Cham: Springer; 2015.

    Book  Google Scholar 

Download references

Acknowledgements

We thank the anonymous reviewers whose comments helped improve and clarify this manuscript.

Funding

This work was supported in part by Grant-in-Aid for Scientific Research (KAKENHI) from the Japan Society for Promotion of Science (#20K04553).

Author information

Affiliations

Authors

Contributions

HH designed the model and the computational framework and analyzed the data. EF verified the analytical methods and supervised the project. All authors discussed the results and contributed to the final manuscript. Both authors read and approved the final manuscript.

Corresponding author

Correspondence to Hikaru Hoshino.

Ethics declarations

Ethics approval and consent to participate

The present study was primarily a computer-based simulation study. As such, the data employed in this study did not require ethical approval.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Hoshino, H., Furutani, E. On the relationship between inhibition and receptor occupancy by nondepolarizing neuromuscular blocking drugs. Theor Biol Med Model 18, 15 (2021). https://doi.org/10.1186/s12976-021-00147-w

Download citation

Keywords

  • Anesthesia
  • Neuromuscular blockade
  • Receptor binding models
  • Nonlinear relationship