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

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., [2][3][4][5][6][7][8][9][10]). 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 IC 50 , 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 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][5][6][7][8], the relative current versus concentration curves were analyzed based on the two-site 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 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 [4][5][6][7][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][14][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.

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

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   [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 * 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 site-selectivity 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 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 := 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 competition-based 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)

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 (  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 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 −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][5][6][7][8]

Results
The upper panel of Fig. 4 shows the relationship between the fraction of activated AChRs given by   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

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 [4][5][6][7][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][14][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 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 IC 50 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 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 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 [4][5][6][7][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 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.