Pharmacodynamic Modeling of Muscle Relaxants: Effect of Design Issues on Results

Our simulations were based on an investigation by Bergeron et al. , 4who gave bolus doses of cisatracurium ranging from 75 to 300 μg/kg (approximately 1.5–6.0 times the ED⁹⁵reported in adults during barbiturate, N²O–opioid anesthesia. 8,9Bergeron et al.  4sampled arterial blood at 1, 2, 3, 4, and 6 min after bolus administration of cisatracurium, followed by increasing intervals until 480 min. Their pharmacokinetic analysis was based on a two-compartment model that assumed that cisatracurium Cp peaked immediately after bolus administration and then decreased monotonically.

Because Bergeron et al.  4did not report individual Cp and twitch tension data, we reconstructed their Cp and effect data for the 75- and 300-μg/kg doses based on their publication. We used their mean pharmacokinetic parameters to simulate a single Cp-versus -time curve (at intervals of 3 s) for each of the 75- and 300-μg/kg doses, assuming monotonic decay. These Cp-versus -time curves and the values of Bergeron et al.  4for each dose for the equilibration rate constant (k^e0) between Cp and those at the effect site for each dose were then used to simulate concentrations of cisatracurium (Ce) at the effect site (the neuromuscular junction). The resulting values of Ce and the values of Bergeron et al.  4for C⁵⁰and the Hill factor (γ) governing the sigmoidicity of the concentration–effect relation for each dose were then used to simulate the resulting twitch depression (effect) data.

To evaluate whether these simulated effect data were consistent with those obtained clinically by Bergeron et al. , 4we determined the time to 98% twitch depression and compared it to the value for onset of Bergeron et al.  4. Our values for time to 25 and 75% twitch recovery for each of the two doses were compared to those of Bergeron et al.  4 

The time at which cisatracurium Cp peaks after bolus administration is not known. Previous studies show that Cp of vecuronium 10and atracurium 11peaks approximately 0.6 min (range, 0.4–0.9 min) after their bolus administration. Therefore, we tested seven additional models in which cisatracurium's Cp peaked at 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, and 1.5 min after its bolus administration. To create these Cp-versus -time curves, we assumed that the model with monotonic decay of Cp yielded accurate values for Cp at all times after Cp peaked (fig. 1); we also assumed that Cp increased in a linear manner from time 0.0 min to the time of its peak.

We then used a semiparametric convolution approach ( Appendix) to relate each of the eight Cp-versus -time curves (i.e. , Cp peaking at 0.0 min plus the seven models in which Cp peaked at 0.2–1.5 min) to the corresponding effect curve for each of the two doses. Each analysis yielded values for C⁵⁰and k^e0. For each cisatracurium dose, these values were plotted against the time to peak Cp.

In muscle relaxant studies, twitch depression is typically quantified by recording the evoked mechanical (or electrical) response of the adductor pollicis muscle on a strip chart. 12To synchronize twitch recording with administration of the muscle relaxant, the investigator might mark the strip chart recording when the drug is administered. These markings might be erroneous; we tested systematic errors of −6 to +6 s.

These analyses were based on the simulated twitch tension data generated for the 75- and 300-μg/kg doses of cisatracurium. Four new twitch tension-versus -time datasets were created by shifting the time base −6, −3, +3, and +6 s. Plasma concentration of cisatracurium was assumed to peak at 0.2, 0.8, or 1.5 min after injection, as above. We then used the semiparametric convolution approach described earlier to relate each of the three Cp-versus -time curves (i.e. , Cp peaking at 0.2, 0.8, or 1.5 min) to each of the five effect curves (four time-shifted effect curves plus one with no shift in time base) for each of the two doses. Each analysis yielded values for C⁵⁰and k^e0. For each cisatracurium dose, these values were plotted against the shift in time base.

A typical strip chart recorder (TA240; Gould Electronics, Valley View, OH) has a 40-mm width for full-scale recording. The width of a recorded signal is approximately 1 mm (2.5% of full scale). The baseline signal may fluctuate several percent of full scale, a result of venous pulsations in the hand and movement artifact. In that investigators typically determine twitch depression by measuring distance from the twitch baseline to the peak of each evoked twitch, small inaccuracies in these measurements are inevitable. We were concerned that administering too small or too large a dose might impact on accuracy of the pharmacodynamic estimates. If the dose is small, the resulting “signal” is small (i.e. , peak effect is < 25% twitch depression); this error (“noise”) represents a larger fraction of the signal than if peak effect were larger. If the administered dose is sufficiently large, few measurements are available during onset. Also, if twitch is ablated for a long period, changes in twitch tension that occur during the drug's distribution phase (e.g. , a decrease from 99.9% twitch depression to 99.1%) cannot be measured accurately. Hence, error in effect measurements might affect estimates of k^e0and, therefore, C⁵⁰. Therefore, we examined whether the magnitude of peak effect (ranging from 20% twitch depression to ablated twitch for a prolonged period) affects the reliability of the pharmacodynamic estimates. Based on the pharmacokinetic–pharmacodynamic data of Bergeron et al.  4for the 75-μg/kg dose, we estimated that the doses producing 20% (ED²⁰), 50% (ED⁵⁰), 80% (ED⁸⁰), and 99% (ED⁹⁹) effect were approximately 30, 37.5, 45, and 75 μg/kg, respectively. Using the pharmacokinetic and pharmacodynamic parameters for the 75-μg/kg dose, we simulated the time course of Cp, Ce, and effect for each of these doses as well as 2 × ED⁹⁹(150 μg/kg) and 4 × ED⁹⁹(300 μg/kg). To simulate different degrees of “noise,” 10 random sets of values of homoscedastic error ‡with an SD of either 2.5 or 10% of full scale (representing 1 or 4 mm, respectively, on a 40-mm twitch recording scale) were generated (Excel; Microsoft, Redmond, WA), i.e. , a different error was simulated for each time interval from 0.0 min to complete recovery of twitch tension. These errors were added to the “true” effect measurements. If the resulting twitch was less than 0 (i.e. , > 100% twitch depression), it was set to 0 (100% twitch depression).

Pharmacodynamic parameters were estimated for each simulated dataset, assuming that plasma concentration of cisatracurium peaked 0.0 min after injection. For each level of error, mean and SD were determined. Then, bias and variability (both expressed as a percentage of the “true” value) in the estimates of C⁵⁰and k^e0were plotted against dose.

Pharmacokinetic–pharmacodynamic simulations and analyses were performed for each individual dataset using NONMEM. 13Values are reported as mean ± SD.

The pharmacokinetic and pharmacodynamic parameters used to simulate the onset and recovery data yielded values consistent with those reported by Bergeron et al.  4(table 1).

As the simulated time to peak Cp increased from 0.0 to 1.5 min, estimates of C⁵⁰for the 300-μg/kg dose decreased; values of C⁵⁰for the 75-μg/kg dose varied minimally (fig. 2). For the 300-μg/kg dose, values for k^e0increased more than twofold as the time of peak Cp increased from 0.0 to 1.5 min. The magnitude of increase was smaller with the 75-μg/kg dose.

Estimates for C⁵⁰resulting from a shift of pharmacodynamic data by −6, −3, +3, and +6 s vary more with the 300-μg/kg dose than with the 75-μg/kg dose (fig. 3). The effect of the timing error is similar for Cp peaking at 0.2, 0.8, and 1.5 min after injection. Similarly, shifting the time base of the pharmacodynamic data by ± 6 s affected k^e0more with the large dose than with the small dose (fig. 4).

The ED²⁰dose was associated with more bias (fig. 5) and variability (fig. 6) in both C⁵⁰and k^e0compared with larger doses. There was minimal bias and variability with doses ranging from ED⁵⁰to 2 × ED⁹⁹. The largest dose, 4 × ED⁹⁹, yielded larger bias and variability than doses ranging from ED⁵⁰to 2 × ED⁹⁹. Both variability and bias were larger with the larger magnitude of error.

Studies of muscle relaxants have provided important insights into various issues of pharmacokinetic–pharmacodynamic modeling, including the need for an effect compartment to accommodate the delay between plasma concentrations and effect and the need to consider that a polyexponential equation describes the initial plasma concentration-versus -time course poorly. We were concerned that certain methodologic issues in pharmacokinetic–pharmacodynamic modeling could impact on the results of these analyses. We used a published dataset to examine several issues, including the impact of different assumptions about the Cp-versus -time profile during the period immediately after drug administration, the impact of systematic error in the timing of effect data, and the impact of dose magnitude in estimating potency.

Our simulations permitted the time of peak Cp to vary from 0.0 to 1.5 min. The latter of these values is larger than that reported for any of the muscle relaxants that have been sampled intensively during the initial minute after bolus dosing. For vecuronium, 10atracurium, 11mivacurium, 14and doxacurium, 3individual time to peak Cp ranges from 0.42 to 0.92 min, and mean values for each drug range from 0.47 to 0.67 min. Thus, it is likely that time to peak Cp for cisatracurium does not exceed 1.0 min. Regardless, as we permitted the time at which Cp peaked to increase from 0.0 to 1.0 min, estimates for C⁵⁰changed markedly with the 300-μg/kg dose of cisatracurium. In contrast, the same change in time at which Cp peaked after a 75-μg/kg dose of cisatracurium had a much smaller effect on C⁵⁰. As a result, for those analyses in which Cp was assumed to peak 1.0 min after bolus administration, the difference in C⁵⁰for the two doses was small; as Cp peaked earlier, the differences became larger (fig. 2). Therefore, if Cp actually peaks as late as 1.0 min after bolus administration of cisatracurium, it is likely that there is, at most, a minimal effect of dose size on C⁵⁰, a finding that contradicts that of Bergeron et al.  4If Cp peaks at 0.6–0.8 min (as is the case for muscle relaxants for which early intense sampling has been performed 3,10,11,14), then C⁵⁰may differ between doses. However, we speculate that the magnitude of difference predicted in our analyses (130 ng/ml for the 75-μg/kg dose and 178 ng/ml for the 300-μg/kg dose, assuming that Cp peaks at 0.6 min) is too small to detect with an unpaired study, as was performed by Bergeron et al.  4When we used the flawed assumption that Cp peaks at 0.0 min (as is typical in most pharmacokinetic–pharmacodynamic studies), the difference in C⁵⁰between the two doses (136 ng/ml for the 75-μg/kg dose and 209 ng/ml for the 300-μg/kg dose) was sufficiently large for Bergeron et al.  4to detect differences between groups. Our simulations suggest that if cisatracurium's C⁵⁰varies with dose and if the analysis of Bergeron et al.  4were performed using a more appropriate method of analysis (i.e. , assuming that Cp peaked at a time later than 0.0 min), their study may not include sufficient subjects to support their conclusion. These results confirm the claim of Ducharme et al.  10that estimation of pharmacodynamic parameters depends on an accurate description of the early time course of Cp. For example, to demonstrate that vecuronium's C⁵⁰varied with dose (as was suggested by Bragg et al. , 5who modeled pharmacodynamics without plasma concentration data), Fisher et al.  6sampled arterial plasma at 0.5 min (in addition to a sampling regimen similar to that of Bergeron et al.  4) and analyzed the data using the “reasonable” assumption that vecuronium's Cp peaks at 0.5 min.

We observed that estimates of C⁵⁰for the smaller cisatracurium dose varied minimally as a function of the time at which Cp peaked. We presume that this occurs because the time at which twitch is ablated with this dose is sufficiently late (4.8 ± 2.3 min) that there is adequate information regarding the plasma concentration time course and effect time course before effect peaks so as to accurately describe the relation between Cp, Ce, and effect. In contrast, with the larger cisatracurium dose, twitch is ablated at 1.8 ± 0.5 min, so that patients studied by Bergeron et al.  4typically had only one plasma sample obtained before twitch was ablated. In turn, the incorrect (but typical) assumption regarding the time at which Cp peaked markedly influenced the input function for the effect compartment, leading to flawed estimates of the pharmacodynamic parameters.

Our simulations indicate the importance of early samples when effect peaks early. If early samples cannot be obtained, pharmacodynamic modeling may be flawed. Another design issue that could lead to incorrect modeling of the early plasma concentration-versus -time course is the use of venous samples. For example, Donati et al.  15demonstrated that atracurium's arterial Cp is markedly larger than venous Cp during the initial 2 min. In that arterial Cp accurately describes the input to the neuromuscular junction, use of venous samples may lead to inaccurate estimates of pharmacodynamic parameters. The inaccuracy of pharmacodynamic parameters is likely to be largest for those drugs with the largest difference between arterial and venous Cp values. If arterial blood cannot be sampled (e.g. , for ethical reasons), then the dosing regimen should be designed so as to minimize the difference between arterial and venous Cp during times critical for the pharmacodynamic analysis. This can presumably be accomplished by administering the muscle relaxant as a brief infusion, as was suggested originally by Sheiner et al.  2 

We assumed that Cp increased in a linear manner from time 0 to the time at which it peaked and then decreased in a monotonic manner. This assumption is slightly flawed. First, intensive sampling during the first minute after bolus administration of several muscle relaxants 3,10,11,14reveals that concentrations of these drugs are not detectable in arterial blood for approximately the first 10 s. Second, the increase in Cp from the first detectable concentration to the peak is not exactly linear. Third, after the peak occurs, there may be several oscillations in Cp, presumably the result of recirculation. Regardless, our approach is markedly closer to the actual Cp-versus -time course than that used by most investigators. First, the commonly-used monotonic decay approach assumes that Cp is maximal at time 0. Second, we speculate that although oscillations exist, their magnitude, coupled with the damping effect of drug transfer to the effect compartment, is probably insufficient to impact on the time course of effect.

To examine the influence of systematic misspecification of the timing of effect data on pharmacodynamic parameter estimation, we simulated a shift of the effect measurements by ± 3 or ± 6 s. Such a shift would occur if the investigator systematically misrecorded the time of drug administration. This might occur under several circumstances. For example, the investigator might administer the drug, then make the notation on the recorder. Also, certain recorders are designed in a manner that prevents access to the paper that is presently being recorded; thus the investigator must wait until the strip chart advances and then make a mark at the estimated time of drug administration.

Although this timing error is trivial compared with the several-hour duration of a neuromuscular study (and it is even small compared with the 1.8 min to twitch ablation with the larger dose of cisatracurium), it influenced the estimates of the pharmacodynamic parameters for the larger dose of cisatracurium. In contrast, misspecification of the timing of effect data with cisatracurium administration influenced the pharmacodynamic estimates minimally with the small dose of cisatracurium, for which twitch was ablated much later, i.e. , at 4.8 min. The findings with the large dose of cisatracurium indicate the importance of accurate timing of dosing events in pharmacodynamic studies, particularly when the drug's onset is rapid.

In that all physiologic measurements involve error, one responsibility of an investigator is to maximize the signal-to-noise ratio, thereby minimizing the impact of noise on parameter estimates. For neuromuscular studies, this noise can be minimized by accurate application of stimulation electrodes, establishing a control response that varies less than 2% for at least 3 min, and maintaining core temperature of 35°C or more and peripheral temperature of 32°C or more. 12We were concerned that the impact of noise would be largest when the maximal effect was small, e.g. , if a small dose of cisatracurium produced only 20% twitch depression. Our simulations support this speculation. However, we also note that as the cisatracurium dose increased beyond the ED⁹⁹, thereby producing a prolonged period during which twitch was ablated, both bias and variability increased (figs. 5 and 6). We offer two possible explanations. First, as the dose increases, the period during which twitch is ablated increases; therefore, there is minimal “information” regarding the relation between changing Cp and effect (i.e. , Cp and Ce may be decreasing, but changes in effect cannot be measured). Second, as the dose increases, onset time shortens so the quantity of effect data during onset decreases (e.g. , during a 1.8-min onset period, there are only 9 measurements of twitch at 12-s intervals, §whereas the 4.8-min onset period permits 24 measurements). Our analysis did not permit us to evaluate the impact of two additional factors that might influence the accuracy of pharmacodynamic estimates with larger doses. First, several investigators 16,17demonstrated that an inadequate stabilization period before muscle relaxant administration results in twitch tension recovery exceeding the baseline value. Whether or not, and the manner in which, an investigator adjusts the data to correct for this overrecovery influences the pharmacodynamic parameters. Second, the longer the period from drug administration to complete recovery, the greater the likelihood is that the twitch tension signal will be unstable, e.g. , by movement of the arm or by changes in body temperature. In that a larger dose yields a longer time to complete recovery, this suggests a disadvantage to administration of doses larger than those necessary. Our simulations indicate that doses ranging from the ED⁸⁰to the ED⁹⁹are probably optimal. This is in contrast to recommendations by other authors, e.g. , Bergeron et al.  4recommend that “doses relevant to the anesthetic practice [presumably in the range of 2 × ED⁹⁹] be used for the estimation of EC50 [termed C⁵⁰in the current article] values.”

One issue of our study design warrants comment. For all pharmacodynamic analyses, we used abundant and “perfect” Cp data. This contrasts to the real-world situation in which many factors contribute to variability in Cp data. These include errors in the timing of blood samples, blood samples being obtained over lengthy intervals (so that their assigned time is not truly representative), contamination of blood samples by intravenous fluids, assay errors, and other factors. It is likely that including this variability in our analyses would have affected our results. In particular, estimates of variability in the estimation of the optimal dose for potency determination are likely to be larger than those reported here.

In summary, we examined several issues that potentially influence pharmacodynamic estimates in neuromuscular studies. We demonstrate that the typical assumption that Cp peaks at 0.0 min after bolus administration results in flawed pharmacodynamic estimates and may lead to misleading conclusions regarding the effect of dose magnitude on these parameters. We demonstrate that timing of effect measurements must be precise, particularly if the onset of effect is rapid. Finally, we demonstrate that doses outside of the range of ED⁸⁰to ED⁹⁹are more likely to yield flawed estimates of the pharmacodynamic parameters than doses within that range.

The authors thank Lewis Sheiner, M.D. (Professor of Laboratory Medicine at the University of California San Francisco, San Francisco, California), for his seminal contributions to pharmacodynamic modeling.

Semiparametric modeling of effect data. A typical approach to pharmacokinetic–pharmacodynamic modeling of muscle relaxants starts by using a compartmental model to describe the plasma concentration (Cp)-versus -time curve. Parameters from this compartmental model and effect data are then used to estimate the pharmacodynamic parameters. In this second step, the pharmacokinetic parameters are used to generate an idealized Cp-versus -time curve. A rate constant k^e0“acts” on these Cp values to generate a effect compartment concentration (Ce)-versus -time curve; the Ce values are then manipulated mathematically using the Hill equation (which involves C⁵⁰and γ) to estimate an effect-versus -time curve. The process starts with initial parameter estimates that are revised in successive iterations until the fit of the model to the data are optimized. The first of these steps, describing the Cp-versus -time curve, cannot always be optimized using a simple compartmental (polyexponential) model in which Cp decreases monotonically. Thus, the compartmental model in the first step can be replaced by one or more functions that describe the Cp-versus -time course. We chose a linear representation for the initial increase in Cp; other nonlinear approaches would also be acceptable. The remainder of the steps in fitting the effect data are identical with the traditional parametric approach and the semiparametric approach.