Population Pharmacodynamic Modeling and Covariate Detection for Central Neural Blockade

The data we used for this analysis came from a study investigating the difference in the effect between two different preparations of plain bupivacaine. One of the two preparations (Astra, Sodertalje, Sweden) contained 5 mg bupivacaine HCl, and 8 mg NaCl per milliliter and the other (Sintetica, Mendrisio, Switzerland) 5 mg bupivacaine HCl and 9 mg NaCl per milliliter. The additional mg/ml NaCl increases the baricity, so that it becomes isobaric at 37 degrees Celsius. [7] 

The study was approved by the Institutional review board and informed consent was obtained from each patient. Consecutive patients scheduled for elective lower limb surgery or inguinal hernia repair were assigned to one of two different preparations of plain bupivacaine. Patients were excluded from the study in case of a contraindication to spinal anesthesia. The spinal puncture was performed using a 25-G needle with the patient in the lateral decubitus position. The amount of injected local anesthetic was chosen by the anesthesiologist to be between 3 and 4.5 ml (15 and 22.5 mg). The anesthesiologist also was free to choose the site of injection, which varied between L1-L2 and L5-S1. The local anesthetic was injected over 30 s. The patient was then turned supine. The level of the central neural blockade was assessed by testing analgesia to pinprick in the mid clavicular line on both sides. During the first hour, block level was recorded at 2, 5, 10, 15, 30, and 60 min in most subjects. Subsequently, the block level was assessed at 90, 120, 150, 180, 210, 240, 270, 300, 330, 360, 390, and 420 min, provided it did not interfere with surgery. All assessments of level of analgesia in a single subject were performed by the same anesthesiologist. The dermatome level measurements were expressed as increasing numbers between 1 (sacral 5) and 30 (cervical 1). For the further analysis the average of the measured levels of both sides were used. The observation was generally discontinued when the patient left the operating room. In some patients, the level of analgesia was tested at a later time on the ward.

This data set had many characteristics of observational data, [4]which made it especially suitable for applying a mixed-effects pharmacodynamic model. Ninety-six subjects with a total of 714 measurements (between 4 and 10 measurements per subject) were available for this analysis.

The development of the structural model was based on the understanding of the physiology of the spread of local anesthetics in the cerebrospinal fluid and on the evaluation of level of analgesia over time. Both the observations of the time course of local anesthetic concentration in the cerebrospinal fluid [8]and the level of analgesia, [3]describe a curve with a relatively fast onset followed by a plateau and a more gradual decline of the response. It is assumed that spinal anesthesia is uniform and ascends from the lowest level, i.e., that whenever a particular level is blocked, all levels beneath will be blocked. Two models, which differed in structure and the number of parameters, were derived. The first model is a combination of a biexponential pharmacokinetic model, describing the onset and offset of effect, and a linear pharmacodynamic model. The second model combines a biexponential pharmacokinetic model with an E^maxtype pharmacodynamic model. With the second, more complicated model it was possible to describe a plateau of the effect. A detailed description of the structural, variance, and residual error model is given in the Appendix.

From the continuous description of the time course of the effect, different clinically important endpoints were estimated for each individual in the study. We estimated the time to reach T10 (onset time), time to reach the maximal level, duration of analgesia (time of level being above T10), and the maximally reached level. The endpoints were predicted using the empirical Bayes estimates of each individual's parameters with a dose of 20 mg bupivacaine. Thus, an estimate was obtained from the dose-normalized effect profiles.

The influence of the covariates age, height, weight, injection site, and preparation of drug on the clinical endpoints were explored. A stepwise search with a generalized additive model [9]allowed for detection of linear and nonlinear relationships between parameters and covariates. A potential nonlinear relationship was modeled with a natural cubic spline with knots at the 33rd and the 66th percentile of the data. In the stepwise search the best model was selected according to the Akaike Information Criterion (AIC). [10]We normalized the covariates with respect to the median to obtain intercepts of the regression line which are the values for the average patient. That is, the clinical endpoints (i.e., duration of analgesia) were modeled to covariates minus the median covariate (i.e., ageⁱ. . .n --median(age), where n is number of individuals).

We tested whether data gathered in a specific patient during the onset of spinal analgesia could predict the offset of spinal analgesia for that patient. Bayes parameters for each individual were calculated using the population parameters and the measurements of level of analgesia up to 30 min. Then the time course of the effect for each subject was predicted up to 420 min. The same analysis was performed using the data up to 60 min. From the subjects with measured data points after 60 min, the association between the measurement and the prediction was assessed by the coefficient of determination (R²). The performance of the population typical parameter set and the two Bayesian parameter sets (30 and 60 min) were compared by examining the prediction of the duration of anesthesia. This was done by calculating the duration of anesthesia (defined as a level of analgesia above T10) with each of these parameter sets and comparing these predictions with the "true" duration. The latter was calculated using the Bayes estimates of the individual parameters based on the complete data set. The statistical inferences were based on an F test.

(Figure 1) shows the raw data. More data points were available from the first 60 min than from later time points. For the initial model-building procedure all of the data were included. For testing correlation of covariates with clinical endpoints the data of one subject was excluded because the level of analgesia did not reach T10. The demographic data are summarized in Table 1.

The data were best described by the "E^max" type model 2. This model, which allows the spinal analgesia to reach a plateau, resulted in a significantly better fit (P < 0.001) compared to the linear pharmacodynamic model 1 based on the likelihood ratio test (see Appendix). Figure 2(A) shows a fit through the data of a representative subject using model 1. The same subject's data fitted with the E^maxtype model 2 are depicted in Figure 2(B). Figure 3shows the relationship between measured and predicted level of analgesia for all subjects. The close correlation between measured and predicted response shows that the population model adequately characterizes the response in each patient. The model parameters are listed in Table 2. The interindividual variability on the parameter F was estimated to be very small and adding the term to the model did not significantly improve the fit. Figure 4shows the time course of all subjects calculated from the empirical Bayes estimates with dose normalized to 20 mg. The maximum level in these 96 subjects varied between T10 and T1. The standard deviation of the residual, unexplained error was 1.4.

The good description of each patient's data by the model and that 4-10 measurements were obtained for each individual allowed accurate estimates of the clinical endpoints for each patient. The distribution of the different measures across the patient population is shown in Figure 5. The mean time to reach T10 was 4.2 min (median 3.0 min). The persons with the long onset time between 12 and 18 min were statistically not different from the others with respect to investigated covariates. An average level of analgesia of T6 was reached (median T6). The mean for the duration was 172 min (median 174 min) and the time to reach the maximal level was 35.5 min (median 33.6 min). Among the eight subjects with the long onset time, only three had a short (< 60 min) duration. The large variance estimates indicate the wide distribution of the clinical response. Covariate analysis aims to explain this variability in the response with patient characteristics and other possible influences on the drug effect.

The covariate analysis is summarized in Table 3. No significant nonlinear relationships were found. The onset time linearly increased with height by about 2 min per 10 cm but decreased with increasing weight by about 45 s per 10 kg. Thus, a short heavy person is expected to have a more rapid onset of spinal blockade than a tall thin one. An 80-yr-old person will reach the maximal level 6 min earlier than a 20-yr-old one. An important part of the variability in the duration also is explained by variability in height. The duration between the extremes of a 150-cm short and a 200-cm tall person would differ by 65 min. Thus, a major part of the interindividual variability in the effect (Figure 5) is explained with the covariates.

Neither the injection site nor the preparation of bupivacaine 0.5% revealed a correlation with the clinical endpoints. The finding that injection site did not prove to be a major covariate is probably owing to the small number of patients (10 of 96) that received an injection at a site other than at L3-L4 or L4-L5.

Model 2 was used to evaluate the success of Bayesian forecasting. The top panel in Figure 6shows a plot of observed levels of analgesia (all levels > 60 min) versus the predicted levels using Bayesian parameter estimates for each subject obtained with the observations up to 30 min. The coefficient of determination (R²) was 0.71. The lower panel in Figure 6with observed versus predicted levels from the data up to 60 min shows a similar strong association. The R²for the after 60-min predictions was 0.72. With the prediction using the population mean parameters an R²of 0.62 was calculated. With the Bayesian parameters obtained from all the data the R²value was 0.9 (Figure 3).

The variance of the difference between the calculated duration and the true duration improved with the parameters based on the onset data (P = 0.0191 for Bayes parameters from data up to 30 min and P = 0.0048 for Bayes parameters from data up to 60 min). Best, median, and worst predicted effect profiles are given in Figure 7for the 30-min Bayes prediction and in Figure 8for the 60-min prediction.

By looking at the cloud of dots representing the measured levels of analgesia over time of all subjects (Figure 1), an underlying pattern appears. Modeling techniques attempt to express this pattern in mathematical formulas. Thus, models are objects that imitate the properties of other real objects in a simpler and more convenient way. Holford [11]compares pharmacokinetic and pharmacodynamic models with road maps, which are useful but limited representations of the actual geography. The population model derived in this study characterized both the analgesic response to intrathecal injection of plain bupivacaine in each patient and the inter-individual variability in response in the patient population. The mixed-effects modeling approach used to derive the population model is especially useful when the data is unbalanced, that is, when a different number of data points for each subject enter the analysis, when the measurements for each subject are taken at different time points and when not enough data points are available in all patients to fit each patient's analgesic response profile independently. [12] 

We treated the data as representing a continuous measurement. Strictly speaking, the levels of analgesia to pinprick represent a 30-level categorical scale. However, with so many categories, treating the measurements as continuous is a statistically reasonable approximation. The statistical model assumes that an increment of the blockade from T2 to T1, for example, is the same as an increment from S1 to S2.

The first model represents a biexponential pharmacokinetic model combined with a linear pharmacodynamic model. One exponential term describes the onset and the other describes the offset of the effect. The additional parameter CO is a coefficient that scales the response. In the second model, the scaling coefficient is called F, because its meaning is mathematically different. The incorporation of an additional parameter "EM" (maximal effect or level) in model 2 improved the fit statistically and graphically as well. This more complicated model resembles a biexponential pharmacokinetic model combined with an E^maxpharmacodynamic model. In our implementation the parameter EM forms the plateau in the individual response. Because this factor limits the maximal effect that the model can predict, it can only be used for the dose range of bupivacaine that was used in the study; that is, between 15 and 22.5 mg. Theoretically, the maximal response to a high dose of bupivacaine will be complete blockade in all patients.

In principle, the models can be built using any mathematical equation. However, the purpose of the structural model development in this study was to obtain an accurate description of each patient's time course of analgesia using a pharmacologically realistic model. The success of the modeling is clearly shown in Figure 2for a representative patient and in Figure 3for all patients. Each curve in Figure 4represents the continuous prediction for one individual, showing the large interindividual variability in the response. The level of analgesia at each point in time can be calculated and interesting clinical values such as maximally attained level or duration of effect can be obtained from the model. The interindividual variability in Figure 4with its large distribution is impressive but expected for bupivacaine. [13]The residual, unexplained error is graphically represented in Figure 3The error includes measurement error and the part of the data that is not explained by the model. For instance, brief positional changes can influence the spread of the local anesthetic. [14]Because this event is not part of the model and was not recorded, it would cause a deviation of the observation from the prediction and become part of the residual error. Nevertheless, it can be assumed that the bulk part is measurement error. Thus, the obtained value indicates that the measured level of analgesia was within an accuracy of about 1.4 dermatomes.

Observational data can be gathered during normal clinical practice. Their weakness is also their strength because these data represent "real-life" data without controlling the independent variables by selecting a subpopulation of patients or restricting the clinical methods. Hence, the conclusions, when drawn carefully, are practically applicable and relevant. Although the data for this study were obtained in a prospective clinical trial, designed to investigate the influence of two different preparations of plain bupivacaine, they have many characteristics of observational data.

Because the patients received different doses of bupivacaine and a correlation between dose and height seemed to exist, the clinical endpoints were calculated for each patient as if they had received the same dose of 4 mg. The possibility of dose normalization and the possibility of predicting the whole time course of effect in subjects with only a minimal number of data points are unique to modeling approaches.

With the generalized additive model nonlinear relationships between the effect and the covariate are detectable. It is an extension to linear regression models, such as multiple regression models, which also have been used for the analysis of data from central neural blockade. [15]Height and weight were the most important covariates. The height of the patients in this study ranged from 148 cm to 193 cm. The predicted difference of all the clinical parameters (Table 3) with this difference in height of 45 cm is clinically significant. This agrees with the results of controlled studies. [13,16,17]The linear combination of height and weight was preferred over the simple regression of body mass index by the Akaike information criterion we used for testing the models. [18]In our analysis, age was correlated with the time required to reach the maximal level. Other studies found increasing duration of blockade with age [19]and increased maximum spread. [20]Although other studies [13,21]have shown that a higher site of injection is associated with a more rapid onset of spinal anesthesia, in our investigation the injection site was not a significant covariate. We performed most of the injections at the L3-L4 and L4-L5 vertebral levels (86 of 96). Thus, although a puncture site of L2-L3 might yield a faster onset than a puncture site of L4-L5, this was not detectable as a statistically significant covariate from our data.

The use of Bayesian forecasting for the prediction of individual pharmacokinetics was proposed by Sheiner et al. [22]based on digoxin data. Maitre and Stanski [23]have shown that with this method the prediction of intraoperative plasma concentrations of alfentanil improves, by using one measured plasma concentration after 80 min of infusion. Obviously, an observation of the response in an individual contains significant information about that person's pharmacodynamic parameters. This information is especially important if the variability of the predictions based on the population parameters is large. Bayesian regression balances the uncertainty in the individual's parameters against the uncertainty in the measured response. With only a few data points per person it is possible to obtain unique estimates of the parameters.

We used information from the onset of effect to predict the offset of effect. The calculated coefficients of determination of 0.71 and 0.72 confirm that the offset of effect is well predicted from the initial measurements. This is a significant improvement over the value obtained from the predictions with the mean population parameters, which is the "best" prediction one would use if no measurements are available for a specific subject.

That the same data set was used to derive the population parameters and for the Bayesian prediction of offset of analgesia may have biased this result to some extent. However, it was the goal to demonstrate the application of Bayesian forecasting. Dividing our data into a learning and test set or using prospectively gathered data under the same experimental conditions would most likely not change the result qualitatively.

We have derived a mixed-effects model for the description of the time course of central neural blockade. Based on the population model a continuous effect profile over time was obtained for each individual, even for those in whom only a small number of data points were available. Clinically important endpoints such as onset time, maximal level of analgesia, time to reach maximal level, and duration were correlated to patient covariates such as age, height, weight, puncture site, and kind of preparation of bupivacaine used. We have found that from the observation of the onset of spinal anesthesia the offset can be predicted.

Two models, differing in structure and the number of parameters, were derived. Model 1 [1.0] described an onset and an offset of effect linearly related to dose by means of two exponential terms: Equation 1

L^pis the level of neural blockade predicted by the model, alpha¹and alpha²are the first-order rate constants of offset and onset, respectively, and CO is a scaling coefficient. This model is equivalent to a combination of a biexponential pharmacokinetic model (such as used to describe plasma concentration time profiles after extravascular administration or effect site concentration time profiles after intravenous drug administration) and a linear pharmacodynamic model. The second model is a nonlinear transformation of model 1, which allows the spinal analgesia to reach a plateau: Equation 2

This model is equivalent to a combination of a biexponential pharmacokinetic model and an E^maxtype pharmacodynamic model. In this model, alpha¹and alpha²describe the rate of onset and offset similar to the first model. alpha¹could be interpreted as a parameter that describes longitudinal migration of the drug within the liquor space and alpha²describes the first-order elimination of drug from the spinal fluid by absorption. F is a coefficient that scales the response, but nonlinearly. The interindividual variability of the model parameters (EM, F, alpha¹, alpha²) was modeled by an exponential variance model, assuming the parameters are log normally distributed: Equation 3where phiⁱis the vector of model parameters for the ith subject, phi the vector of the typical (mean) parameter values for the population, and e^etaiexpresses the random difference between the parameters values of the ith subject and the population mean parameters. Values of etsⁱare assumed to be identically independently normally distributed with mean zero and diagonal variance-covariance matrix Omega with diagonal elements (omega sub 1², . . ., omega^m2). The log normal distribution of the parameters is a common assumption in pharmacokinetic and pharmacodynamic modeling because it prevents the parameters from having nonphysiologic negative values and because the distribution of the parameters tends to be skewed, it is a better approximation than the normal distribution. This diagonal variance covariance matrix assumes that there is no correlation between the parameters.

The residual error, which includes the measurement error, was characterized with an additive model as follows: Equation 4where L sub o,ij is the jth observed level of analgesia from the ith individual, L^p,ij the predicted level for this individual according to model 1 or 2 and the values of epsilon^ijare assumed to be identically, independently and normally distributed with mean zero and variance sigma sup 2. This so called homoscedastic error model assumes that the residual error is of the same magnitude over the whole range of measurements.

The model was fitted to the data using the mixed effect modeling program NONMEM version IV. [24],* The first-order conditional estimation method was used to estimate the values of the population parameters, phi, Omega, sigma². Based on the population parameters and the individual's data, empirical Bayes estimates were obtained of the values phiⁱfor each individual. [22,25] 

The choice between models 1 and 2 was done by visual inspection of the population fit, the individual fits based on the empirical Bayes parameter estimates of each individual and by a likelihood ratio test. Minus twice the log of the likelihood (-2LL) ratio of a reduced and a full model is approximately distributed as a chi²random variable with degrees of freedom equal to the number of parameters being added. A decrease of more than 6.6 points in the likelihood ratio is significant at the P < 0.01 Level.

*Beal SL, Sheiner LB: NONMEM User's Guide. San Francisco, University of California San Francisco, 1979.