 Research
 Open Access
 Published:
On the relationship between inhibition and receptor occupancy by nondepolarizing neuromuscular blocking drugs
Theoretical Biology and Medical Modelling volumeÂ 18, ArticleÂ number:Â 15 (2021)
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 concentrationeffect relationships of NDNBs.
Background
Nondepolarizing neuromuscular blocking drugs (NDNBs) inhibit neuromuscular transmission by competing with acetylcholine (ACh) for binding sites at the postjunctional 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., [2â€“10]). In particular, in [7, 10], in vitro experiments have been conducted through macroscopic current recordings from outsideout patches of BOSC23 cells or Xenopus oocytes where human adult (rather than mouse adult or embryonic) muscletype 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 IC_{50}, the drug concentration needed to produce a 50% inhibition of the AChinduced 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 IC_{50} estimated by regression analysis using the Hill equation for the relative current I_{antag}/I_{0} given by
where I_{0} and I_{antag} stand for the currents measured in the absence and presence of the drugs, respectively. The parameter n_{H} 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 [4â€“8], the relative current versus concentration curves were analyzed based on the twosite receptor binding model given by the following equation:
where K_{D1} and K_{D2} stand for the dissociation equilibrium constants for NDNBs binding to the first and second sites of an AChR, respectively. Since the righthand 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 [4â€“8], 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, 13â€“15] 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 3letter 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 k_{assocAi} and k_{assocDi}, for site #i(i=1, 2), respectively. Similarly, k_{dissAi} and k_{dissDi} 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 k_{open} and k_{close}, respectively. Then, the concentrations of these complexes can be calculated as a function of time by solving the following ordinary differential equations:
where [ D] stands for the concentration of NDNB at the effect site, and k_{decay} 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
where [ R]_{total} stands for the concentration of the postjunctional 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Ã—10^{7}) at the end plates of human deltoid muscle reported in [16] and the volume of the synaptic cleft (4.5Ã—10^{13}L) 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Ã—10^{6} 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 k_{decay} was determined such that the halflife 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.
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^{âˆ’8}M, and 1.0Ã—10^{âˆ’7}M. 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].
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 I_{antag}/I_{0} 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:
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 k_{dissD1} and k_{dissD2} of an NDNB are much smaller than other rate constants such as k_{dissA1},k_{dissA2} and k_{decay}. 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, x_{A}, can be defined as
where K_{A1} stands for the dissociation equilibrium constant for ACh with the site #1 given by k_{dissA1}/k_{assocA1}. Dimensionless variables \(x^{\ast }_{\text {ARA}}\) and x_{DRD} for the concentrations of the complex ARA^{âˆ—} and DRD are given by
Similarly, the dimensionless variables x_{ARA},x_{ARD},x_{ARD},x_{DRA},x_{ARO},x_{ORA},x_{DRO}, and x_{ORD} 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):
where Îº_{A1},Îº_{A2}, Îº_{open} and Îº_{close} are the dimensionless parameters representing normalized rate constants given by
The parameters Î»_{A1},Î»_{A2}, and Î¼_{A}are the dimensionless parameters representing affinities of ACh to the binding sites of the AChR given by
with K_{A2}:=k_{dissA2}/k_{assocA2}, and finally, Î¼_{D} and Î´ are the dimensionless parameters representing the siteselectivity and concentration of NDNB, respectively, given by
with K_{D1}:=k_{dissD1}/k_{assocD1} and K_{D2}:=k_{dissD2}/k_{assocD2}. Furthermore, O_{DRD},O_{ORD} and O_{DRO} are functions of Î´ standing for the fractions of the complex DRD, ORD, and DRO before the release of ACh given by
The total occupancy O_{total} of the AChRs by NDNB is given by
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 siteselectivity 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}:=k_{dissD1}/k_{decay} and Îº_{D2}:=k_{dissD2}/k_{decay} 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 O_{total} is the same, the fraction of the activated AChRs will decrease with the increase in the occupancies O_{ORD} and O_{DRO} because free ACh will decrease due to association with ORD and DRO. Note that the ratio of ORD and DRO to O_{total} is determined by the parameter Î¼_{D}. Finally, the difference between the receptor binding model (2) and the competitionbased models (3) and (8) can be clarified. By using the notation introduced above, the model (2) can be rewritten as follows:
Thus, the model (2) is clearly based on the assumption that the inhibition is proportional to the total occupancy O_{total}, while the models (3) and (8) describe the effect of the partial occupancies O_{ORD} and O_{DRO} 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^{âˆ’14}M and 1.0Ã—10^{âˆ’5}M that are spaced equidistantly on a logarithmic scale, and [ D]=0) and plotted versus the total occupancy O_{total}. 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 k_{dissA1},k_{dissA2},k_{assocA1}, and k_{assocA2} presented in Table 1 were determined by experiments using mouse adult AChRs. However, it is known that the EC_{50}, the concentration of ACh where the half of the AChRs are activated, for human adult AChRs is smaller than that for mouse adult AChRs: EC_{50}=8.48Ã—10^{âˆ’6}M or 7.91Ã—10^{âˆ’6}M for human adult AChRs [10] and EC_{50}=1.6Ã—10^{âˆ’5}M 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 K_{A1} and K_{A2} by changing k_{assocA1} and k_{assocA2} while k_{dissA1} and k_{dissA2} 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 K_{A1} and K_{A2} are equal. Thus, we investigate the effect of varying the values of K_{A1}=K_{A2}=K_{A}, 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 k_{dissD1} and k_{dissD2} are small enough. The value of parameter Î¼_{D} was assigned by changing the value of k_{assocD2}, while k_{dissD1},k_{dissD2} and k_{assocD1} kept constant at the values shown in Table 1. Also, the values of the dissociation constants k_{dissD1} and k_{dissD2} are varied to explore the difference between the original model (3) and the reducedorder 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^{âˆ’5}M and [ A]_{init}=7.75Ã—10^{âˆ’6}M, 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 [4â€“8] such that nearly 50% (nearly EC_{50}) 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 O_{total} evaluated under various settings of the dissociation equilibrium constants for ACh. The values of K_{A1}=K_{A2}=K_{A} is one of 1.0Ã—10^{âˆ’4}M,3.16Ã—10^{âˆ’5}M and 1.0Ã—10^{âˆ’5}M, 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 O_{total} in all the cases. Furthermore, it can be confirmed that the extent of nonlinearity increases with the decrease in the equilibrium constant K_{A}. In the subsequent analyses, to highlight the effects of varying other parameters, we assigned the value of 1.0Ã—10^{âˆ’5}M to K_{A}, 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 siteselectivity 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 O_{total} are nonlinear in all the cases, and the extent of nonlinearity decreases with the increase in Î¼_{D}.
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 O_{total} 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 O_{total} 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 concentrationresponse 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 O_{total}. 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 O_{total} becomes linear.
Finally, the effect of varying the dissociation rate constants k_{dissD1} and k_{dissD2} of NDNB was examined under the setting of identical rate constants, i.e. k_{dissD1}=k_{dissD2}=k_{dissD}. The upper panel of Fig. 6 shows the relationship between the fraction of activated AChRs, [ARA^{âˆ—}]_{peak}/[ARA^{âˆ—}]_{max}, and the receptor occupancy O_{total} 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 k_{dissD} is as large as k_{dissA1},k_{dissA2}, and k_{decay}. This implies that the model simplification performed in this paper can be validated for a wide range of settings of k_{dissD}. 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 k_{dissD}. Furthermore, the relationships between [ARA^{âˆ—}]_{peak}/[ARA^{âˆ—}]_{max} and O_{total} are no longer linear when the rate constant k_{dissD} 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 k_{dissD}.
Discussion
We theoretically and numerically investigated the relationship between inhibition and receptor occupancy by NDNBs. While the twosite receptor binding model (2), which assumes a linear relationship between the inhibition and the receptor occupancy by an NDNB, has been statistically tested in [4â€“8] 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, 13â€“15] 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 siteselectivity 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 postsynaptic 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 siteselectivity 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 IC_{50} for several clinically used NDNBs, cisatracrium, 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 IC_{50} 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 O_{total} 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 IC_{50} increases with the decrease in the concentration of ACh. This shift in the IC_{50} 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 twosite 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 [4â€“8], 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 concentrationeffect 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 siteselectivity 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 modelorder 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, x_{A}, can be defined as
where K_{A1} stands for the dissociation equilibrium constant for ACh with the site #1 given by k_{dissA1}/k_{assocA1}. Dimensionless variables \(x^{\ast }_{\text {ARA}}\) and x_{DRD} for the concentrations of the complex ARA^{âˆ—} and DRD are given by
Similarly, the dimensionless variables x_{ARA}, x_{ARD},x_{ARD},x_{DRA},x_{ARO},x_{ORA},x_{DRO}, and x_{ORD} 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:
where Îº_{A1},Îº_{A2},Îº_{D1} and Îº_{D2} are the dimensionless parameters representing the rates of the both dissociation and association of the ligands given by
and Î»_{A1},Î»_{A2}, and Î¼_{A}are the dimensionless parameters representing affinities of ACh to the binding sites of the AChR given by
with K_{A2}:=k_{dissA2}/k_{assocA2}, and finally, Î¼_{D} and Î´ are the dimensionless parameters representing the siteselectivity and concentration of NDNB, respectively, given by
with K_{D1}:=k_{dissD1}/k_{assocD1} and K_{D2}:=k_{dissD2}/k_{assocD2}. 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 (k_{dissD1},k_{asoocD1},k_{dissD2},k_{assocD2}) 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 multipletimescale property of the model. Such a multiscale property arises when the dissociation rate constants k_{dissD1} and k_{dissD2} of an NDNB are much smaller than other rate constants such as k_{dissA1},k_{dissA2} and k_{decay}. By considering the limit Îº_{D1}, Îº_{D2}â†’0, the following equations hold from the Eq. 17:
Thus, the values of x_{DRD},x_{ORD}+x_{ARD} and x_{DRO}+x_{DRA} do not change in the reducedorder 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 x_{DRD},x_{ORD},x_{DRO},x_{ARD} and x_{DRA} are given by
where O_{DRD},O_{ORD} and O_{DRO} stand for the fractions of the complex DRD, ORD, and DRO, respectively, at the steady state. Also, the total occupancy O_{total} of the AChRs by the NDNB is given by
By using the equations from (21) to (23), the model (3) can be reduced to the following form:
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
Pardo MC, Miller RD. Basics of Anesthesia, 7th edn. Philadelphia: Elsevier; 2018.
Yost CS, Winegar BD. Potency of agonists and competitive antagonists on adult and fetaltype nicotinic acetylcholine receptors. Cell Mol Neurobiol. 1997; 17(1):35â€“50. https://doi.org/10.1023/A:1026325020191.
Paul M, Kindler CH, Fokt RM, Dresser MJ, Dipp NCJ, Yost CS. The potency of new muscle relaxants on recombinant muscletype acetylcholine receptors. Anesth Analg. 2002; 94(3):597â€“603. https://doi.org/10.1097/0000053920020300000022.
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.
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.
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.
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.
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.
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 Trainoffour Fade. Anesthesiology. 2006; 105(3):521â€“33. https://doi.org/10.1097/0000054220060900000016.
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.
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.
Gabrielsson J, Weiner D. Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications, 5th edn. Stockholm: Swedish Pharmaceutical Press; 2016.
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/S00063495(96)79834X.
Auerbach A, Akk G. Desensitization of mouse nicotinic acetylcholine receptor channels : A twogate mechanism. J Gen Physiol. 1998; 112(2):181â€“97.
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.
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.
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/S00063495(79)852492.
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.
Hobbiger F. Neuromuscular Junction In: Zaimis E, Maclagan J, editors. Berlin: Springer: 1976. p. 487â€“581.
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.
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.
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/articlepdf/103/1/108/17424609/aep150.pdf.
Langtangen HP, Pedersen GK. Scaling of Differential Equations. Cham: Springer; 2016.
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.
Kuehn C. Multiple Time Scale Dynamics. Cham: Springer; 2015.
Acknowledgements
We thank the anonymous reviewers whose comments helped improve and clarify this manuscript.
Funding
This work was supported in part by GrantinAid for Scientific Research (KAKENHI) from the Japan Society for Promotion of Science (#20K04553).
Author information
Authors and Affiliations
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
Ethics declarations
Ethics approval and consent to participate
The present study was primarily a computerbased 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.
About this article
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/s1297602100147w
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1297602100147w
Keywords
 Anesthesia
 Neuromuscular blockade
 Receptor binding models
 Nonlinear relationship