In previous studies, the authors reported on the absorption and disposition kinetics of levobupivacaine and ropivacaine. The current study was designed to develop a population pharmacokinetic-pharmacodynamic model capable of linking the kinetic data to the analgesic effects of these local anesthetics (i.e., sensory neural blockade).
A disposition compartmental model was fitted to concentration data of the intravenously administered deuterium-labeled anesthetics, and a model consisting of two parallel absorption compartments and the identified disposition compartments was fitted to concentration data of the concomitantly epidurally administered unlabeled anesthetics. The epidural segments were modeled by individual central and peripheral absorption compartments and effect sites, which were fitted to the simultaneously acquired pinprick data. A covariate model incorporated the effects of age.
The threshold for epidural anesthesia increased from the lower to the higher segments. The central effect compartment equilibration half-lives were approximately 15 min for levobupivacaine and 25 min for ropivacaine. For levobupivacaine, age reduced the equilibration half-lives at all segments; for ropivacaine, age increased the anesthetic sensitivity at segments T12 and higher.
A population pharmacokinetic-pharmacodynamic model was developed that quantitatively described sensory blockade during epidural anesthesia, including the effects of age. The model may be useful to individualize dose requirements, to predict the time course of sensory blockade, and to study new local anesthetics.
EPIDURAL anesthesia is obtained by the injection of a local anesthetic drug into the epidural space. The clinical characteristics of the ensuing sensory neural blockade, such as onset time, intensity, and duration, depend directly on the changes in the concentration of the local anesthetic at the axonal membrane, which are dependent on pharmacokinetic factors. The rate of systemic absorption of local anesthetics gives some indication of the relation between neural blockade and the amount of drug remaining at or near the site of injection. In the past two decades, we indirectly assessed epidural drug concentrations by estimation of the time course of systemic absorption from the epidural space in humans.1–5
In our previous studies, the pharmacokinetic data were analyzed for each individual separately. However, population pharmacokinetic–pharmacodynamic (PK-PD) analysis of epidural anesthesia is important because it enables the description of both within- and between-subject variability, enables the development of predictive models, and, consequently, improves therapeutic outcome.6Schnider et al. 7developed a population pharmacodynamic model to describe the time course and blockade level of spinal anesthesia. The objective of the current study was to develop a population PK-PD model of epidural anesthesia based on pharmacokinetic absorption profiles. In a first attempt, we applied the model to levobupivacaine and ropivacaine data and determined the impact of age on the model parameters.
Materials and Methods
Study Design
For the development of the epidural PK-PD model, data were used from two previous studies3,5on the epidural injection of levobupivacaine and ropivacaine (see table 1for patient characteristics). In these studies, arterial blood samples were obtained before, during, and after the epidural injection of the local anesthetic (at the L3–L4 interspace; table 2). The maximum duration of sampling was 24 h. When after epidural injection satisfactory anesthetic conditions (i.e. , the presence of a bilateral sensory blockade, assessed by pinprick) were obtained (usually 15–25 min after the injection), the same but now deuterium (2H3)–labeled local anesthetic was administered intravenously (table 2). Blockade assessments were made every 5 min during the first 30 min, every 15 min for the next 3.5 h, and subsequently every 30 min for the next 6 h. Refer to our previous publications for further details. For the current study, there were no additional protocols that needed approval of an institutional review board.
Development of the Pharmacokinetic Model
The reanalysis of the pharmacokinetic data from the previous studies had two aims: (1) to obtain population pharmacokinetic models for the local anesthetics levobupivacaine and ropivacaine together with a covariate analysis (to enable individualization at the pharmacokinetic level) and (2) to obtain optimal individual absorption parameters for a pharmacodynamic model that describes the effects of those anesthetics on sensory blockade (and to enable individualization on the pharmacodynamic level). The pharmacokinetic analysis was performed in four steps: (1) A three-compartment structural model was fitted to the concentration data of the intravenously administered labeled local anesthetics (fig. 1). (2) A model including covariate age was developed. (3) A model, consisting of two parallel absorption compartments and the model of disposition developed in step 2, was fitted to concentration data of the epidurally administered unlabeled local anesthetic (fig. 1; note that the disposition parameters were fixed to their empirical bayesian [see Statistical Analysis section] values obtained in step 2). (4) A model including covariate age was developed.
Fig. 1. Schematic diagram of the pharmacokinetic–pharmacodynamic model. By intravenous administration of a labeled local anesthetic (LA), its disposition is identified and described by pharmacokinetic parameters V 1(central volume) and k ij(distribution rate constants between the compartments; V 2and V 3are derived parameters). By epidural administration of the unlabeled LA, its absorption is identified and described by parameters F 1, F 2, k a1, and k a2. Assumption 1: At time T = 0, the unlabeled LA administered at the L3–L4 interspace spreads instantaneously across all segments (for segments T9 and T10, model details are given). Consequently, there is only drug distribution between the central absorption compartments for T = 0, not for T > 0 ( dashed arrows ), and there is no distribution between the effect sites via the cerebrospinal fluid (CSF) ( dotted arrows ). Assumption 2: For T > 0, the local absorption profiles are identical to the identified global absorption profile. Consequently, the two parallel absorption compartments are equivalent with the central and peripheral absorption compartments as described by parameters V c, k abs, k cp, and k pc( V pis a derived parameter). When the effect site concentration C e( t ) is above a threshold C thr, a stimulus S is not perceived because it is blocked ( e.g. , at T10); when C e( t ) is below C thr, it is not blocked and can be perceived ( e.g. , at T9).
Fig. 1. Schematic diagram of the pharmacokinetic–pharmacodynamic model. By intravenous administration of a labeled local anesthetic (LA), its disposition is identified and described by pharmacokinetic parameters V 1(central volume) and k ij(distribution rate constants between the compartments; V 2and V 3are derived parameters). By epidural administration of the unlabeled LA, its absorption is identified and described by parameters F 1, F 2, k a1, and k a2. Assumption 1: At time T = 0, the unlabeled LA administered at the L3–L4 interspace spreads instantaneously across all segments (for segments T9 and T10, model details are given). Consequently, there is only drug distribution between the central absorption compartments for T = 0, not for T > 0 ( dashed arrows ), and there is no distribution between the effect sites via the cerebrospinal fluid (CSF) ( dotted arrows ). Assumption 2: For T > 0, the local absorption profiles are identical to the identified global absorption profile. Consequently, the two parallel absorption compartments are equivalent with the central and peripheral absorption compartments as described by parameters V c, k abs, k cp, and k pc( V pis a derived parameter). When the effect site concentration C e( t ) is above a threshold C thr, a stimulus S is not perceived because it is blocked ( e.g. , at T10); when C e( t ) is below C thr, it is not blocked and can be perceived ( e.g. , at T9).
The absorption compartments were characterized by the following parameters: F 1, F 2, t ½,a1, and t ½,a2, denoting the fractions absorbed and absorption half-lives of the two compartments. We used a parameterization with F and F′ 1where F1=F ·F′ 1and F 2=F · (1 −F′ 1), because F , the total absorption, should be around 1, and interindividual variability in F 1is likely to be counteracted by approximately the same variability in F 2.
Development of the Pharmacodynamic Model
The Epidural Segments.
Each segment was modeled by its own central and peripheral absorption compartments (fig. 1). After a unit dose, the amount of anesthetic in a central absorption compartment A c(t ), is given by
with k a1= log(2)/t ½,a1and k a2= log(2)/t ½,a2for convenience. The concentration Cc (t ) in a central absorption compartment after a dose consisting of amount A (fraction of total dose; see Main Assumptions section) of anesthetic is then described by
with f 1=F 1·A /V cand f 2=F 2·A /V c, where V cis the volume of that central absorption compartment. Effect sites were postulated to exist at each segment to account for lags between the absorption compartment concentration and effect, quantified by equilibration rate constant k e0. The effect site concentration C e(t ; k e0) is then given by
Epidural blockade (analgesia with respect to pinprick) was assumed to occur when C e(t ; k e0) exceeds a concentration threshold C thr(fig. 1).
Types of Observations.
The acquired pharmacodynamic data contained three types of observations:
an interval (t on, t off) with C e(t on; k e0) ≥C thrwhere a blockade occurred, and C e(t off; k e0) < C thrwhere the blockade disappeared;
as type 1, but at the last assessment time a blockade was still present; and
there was no blockade at all assessment times.
Although this might seem counterintuitive, a type 3 observation does provide information, because no blockade occurs because C thris relatively high, V cis large, and/or A is small. With only type 3 observations in the data set, only a lower limit of C thrcan be determined for a fixed k e0. Together with type 1 and type 2 observations, type 3 observations do allow for more accurate estimates of C thrand k e0. Evidently, also the reverse is true: When type 3 observations would be discarded, estimates of C thrwould be biased downward.
Main Assumptions.
According to the pharmacodynamic model, each segment is described by the following six parameters: F 1, F 2, k a1, k a2, k e0, and C thr. In addition, rate constants should be defined that describe the transport of the anesthetic between segments. However, there were at most only two informative observations per dermatome, and the pharmacokinetic analysis provides global rather than local absorption process parameters. So there are more parameters than observations. The following assumptions allowed us to proceed (fig. 1):
The longitudinal spread of the anesthetic solution across the segments of the epidural space occurs instantaneously.
The parameters F 1, F 2, k a1, and k a2are equal for each segment, albeit that F 1and F 2are interpreted with respect to the fraction of the dose that is present in each segment after the initial spread.
Under these assumptions, both the local (segmental) and global (total epidural; obtained by adding all local absorption profiles) processes are described by the same parameters as those obtained from the pharmacokinetic analysis. If the local absorption or transport processes would be subject to large intrasegmental variability, it is unlikely that the global absorption profile would display such a clear biphasic pattern.1,8Furthermore, distribution via the cerebrospinal fluid can be assumed to be minimal.9We therefore believe that the postulated assumptions are reasonable.
After the initial spread, the amounts (A ) of anesthetic present in each of the segments, and the central absorption volumes (V c) at those segments, are not simultaneously identifiable. Only a combination parameter A thr=C thr·V c/A can be estimated. For parameter estimation, we therefore used (cf. equations 1 and 3):
and assumed that a blockade occurs when A e(t ; k e0) ≥A thr. Parameter A thris a measure of the anesthetic sensitivity. With two observations per dermatome, the remaining unknown parameter, k e0, can be estimated. Parameters A thrand k e0were assumed to be lognormally distributed across the population. Because of the assumptions postulated, they can be estimated for each of the segments independently. Further details are given in the appendix.
Covariate Analysis
The pharmacokinetic and pharmacodynamic parameters (θ, indexed by i ) were multiplied by factors containing log-transformed and centered (at 50) covariate age:
where θpop,i denote population (typical) values; c = log(age/50), αAGE,i denote covariate coefficients; and ηi denote residual interindividual variabilities.
Statistical Analysis
Statistical analysis of the pharmacokinetic data were performed with the NONMEM version VI software package (a data analysis program for nonlinear mixed effects modeling)10using its routine ADVAN11 for three-compartmental models of the intravenous data and analytical expressions (for processor time reduction) in $PRED for the five-compartmental model of the epidural data. Nonzero variance terms of the ηs were determined for the population model without covariate age; the structure was assumed to be the same for the model with the covariate included (mainly a decrease in variances is to be expected). Constant relative errors were assumed for the concentration measurements. All numerical results are given with three significant digits (conform NONMEM output). The empirical bayesian estimates (provided by NONMEM) of the disposition parameters were used for the estimation of the absorption parameters; the empirical bayesian estimates of the absorption parameters were used for the estimation of the pharmacodynamic parameters.
Statistical analysis of the pharmacodynamic data was performed using software written in the computer language C by one of the authors (E.O.) using the free GNU Scientific Library (GSL).∥Additional information regarding this is available on the Anesthesiology Web site at http://www.anesthesiology.org. With distributions of parameters A thrand k e0postulated, the probability of observing the data per dermatome can be computed by integration across the (A thr, k e0) plane where it is possible that the data are as observed. The likelihood of observing the data are the product of the probabilities of all observations per dermatome. By maximizing this likelihood, the parameters of the lognormal distributions of A thrand k e0can be obtained. (Note that NONMEM cannot be used when the integrand used for integration across the (A thr, k e0) plane is not continuous.) The boundaries of the 95% confidence intervals of the parameter estimates were determined by those parameter values that increase the objective function by 3.84 points (i.e. , the likelihood profile method). Further details are given in the appendix.
The improvement of the model fit by inclusion of separate parameters for the different anesthetics and covariate age via forward selection was based on Akaike’s Information-Theoretic Criterion (AIC).11The best model was defined as the one with the lowest value of the criterion, where AIC = MVOF + 2p , where MVOF is the −2 log-likelihood and p is the number of parameters. The model with the lowest value of the criterion is assumed to have the best predictive properties.12,13However, to present a classic measure of the significance of difference from zero of each of the final covariate coefficients of the pharmacokinetic models, likelihood ratio tests were performed (Δobjective function value > 3.84 corresponds to P < 0.05; Δobjective function value > 6.63 corresponds to P < 0.01). To present a classic measure of significance of the selected covariate coefficients of the pharmacodynamic models, their 95% confidence intervals were graphically displayed together with zero ordinate lines.
Except for F′ 1, lognormal distributions were postulated for all parameters. The value of F′ 1was constrained between 0 and 1, which was accomplished by using the inverse logit transformation on its η.
With distributions of A thrand k e0available, the probability of blockade at each segmental level can be determined as a function of time and age. Blockade probabilities were computed with respect to a standardized dose of 100 mg for ropivacaine and 125 mg for levobupivacaine.
Results
Pharmacokinetic Analysis
Best, median, and worst model fits of the labeled local anesthetic concentration data (panels 1, 2, and 3, respectively) according to the coefficient of determination, together with the corresponding model fits of the unlabeled local anesthetic concentration data, are shown in figure 2. The worst model fits of the unlabeled local anesthetic data still provide adequate inputs for the pharmacodynamic models.
Fig. 2. Best, median, and worst model fits ( filled circles ; panels 1, 2, and 3, respectively) according to the coefficient of determination ( R 2) for unlabeled levobupivacaine and ropivacaine concentration data (panels A and B , respectively), together with the corresponding model fits of the labeled local anesthetic concentration data ( open circles ). The numbers in the subplots denote patient identification number/ R 2unlabeled/ R 2labeled.
Fig. 2. Best, median, and worst model fits ( filled circles ; panels 1, 2, and 3, respectively) according to the coefficient of determination ( R 2) for unlabeled levobupivacaine and ropivacaine concentration data (panels A and B , respectively), together with the corresponding model fits of the labeled local anesthetic concentration data ( open circles ). The numbers in the subplots denote patient identification number/ R 2unlabeled/ R 2labeled.
Levobupivacaine and ropivacaine disposition parameter estimates are presented in table 3. The only parameter that differed significantly between the two anesthetics was k 31. A significant effect of age on V 1, and therefore on elimination clearance, was detected. In addition, a significant effect of age on intercompartmental distribution rate constant k 12was detected.
Levobupivacaine and ropivacaine absorption parameter estimates are presented in table 4. Parameters t ½,a1and t ½,a2were significantly different between the two anesthetics. Significant effects of age on t ½,a1and F′ 1were detected. Parameter F was approximately 10% higher than 1 for both anesthetics (and had small interindividual variability around it).
Pharmacodynamic Analysis
Figure 3presents the population estimates of t1/2,ke0 (converted from the k e0values for convenience of interpretation) and A thrwith their 95% confidence intervals for each segment separately. Both parameters were significantly different between levobupivacaine and ropivacaine. At the L2 segment level, parameter t1/2,ke0 was approximately 5 min for levobupivacaine and 9 min for ropivacaine, which approximately doubled at higher and tripled at lower dermatome levels. Parameter A thrincreased from S5 to T5, with the highest overall level for levobupivacaine. The variabilities in k e0and A thrwere approximately 45% and 25%, respectively (data not shown).
Fig. 3. Equilibration half-lives and absorption thresholds, with their 95% confidence intervals, at segments S5–T5 of the two local anesthetics under study. To guide the eye , dashed lines were drawn at the parameter values at level L2 (epidural injection was at the L3–L4 interspace).
Fig. 3. Equilibration half-lives and absorption thresholds, with their 95% confidence intervals, at segments S5–T5 of the two local anesthetics under study. To guide the eye , dashed lines were drawn at the parameter values at level L2 (epidural injection was at the L3–L4 interspace).
Figure 4presents the effects of age on the population values. For levobupivacaine, age increase the equilibration rate constants k e0at all segments. For ropivacaine, age increased the anesthetic sensitivity A thrat segments T12 and higher.
Fig. 4. Effect of age on k e0and A thrat segments S5–T5 as quantified by the respective covariate coefficients ( equation 5) together with their 95% confidence intervals, for the two local anesthetics under study.
Fig. 4. Effect of age on k e0and A thrat segments S5–T5 as quantified by the respective covariate coefficients ( equation 5) together with their 95% confidence intervals, for the two local anesthetics under study.
The nature of the pharmacodynamic data are such that the model fits exactly for each individual and dermatome (with two measurements (onset and offset of effect), and two parameters (t1/2,ke0 and A thr)); only their probability distributions across the population can be estimated. The probabilities of blockade (after a dose of 100 mg ropivacaine and 125 mg levobupivacaine) as a function of time and age of the patient are presented in figure 5. The effects of increasing age are clearly visible at the higher segments (T9 and higher) where the probability of block dramatically increases, and by the increase of duration of block. To facilitate the interpretation of figure 5, a three-dimensional view of the blockade probabilities for levobupivacaine in a patient aged 50 yr is given in figure 6. Additional information regarding this is available on the Anesthesiology Web site at http://www.anesthesiology.org.
Fig. 5. Blockade probabilities for the two local anesthetics (ropivacaine dose 100 mg, levobupivacaine dose 125 mg) and the effect of age. The size of the dots is proportional to the probability of blockade; in addition, 90%, 75%, and 50% isoeffect lines were drawn ( solid , heavily , and lightly dashed lines , respectively).
Fig. 5. Blockade probabilities for the two local anesthetics (ropivacaine dose 100 mg, levobupivacaine dose 125 mg) and the effect of age. The size of the dots is proportional to the probability of blockade; in addition, 90%, 75%, and 50% isoeffect lines were drawn ( solid , heavily , and lightly dashed lines , respectively).
Fig. 6. Blockade probabilities for levobupivacaine for age 50 yr as a function of time and segment level.
Fig. 6. Blockade probabilities for levobupivacaine for age 50 yr as a function of time and segment level.
Discussion
We developed a predictive population PK-PD model of epidural anesthesia and determined its parameters, taking into account the effects of age, for the local anesthetics levobupivacaine and ropivacaine. The pharmacokinetic analysis yielded similar population parameter estimates with respect to our previous studies. The PK-PD model was able to predict the probability of block and duration of anesthesia per segment. The covariate analysis showed that age had a significant and clinically important effect on the spread and duration of analgesia.
We previously determined absorption kinetics after epidural administration of bupivacaine, levobupivacaine, and ropivacaine.1–5Stable isotope–labeled analogs were intravenously infused to obtain disposition kinetics. Subsequently, with the use of the disposition kinetics and the plasma concentration of the epidurally administered unlabeled local anesthetic, the absorption kinetics could be obtained by deconvolution. The plasma concentration–time curves of individual patients were adequately described by fitting directly the aggregated model of two parallel first-order absorption compartments and the disposition profile (a two- or three-compartmental model). These analyses allowed the determination of the profiles of the different local anesthetics and the effects of age.
However, these individual analyses may be hampered by interindividual variability, caused by sex differences and genetic, environmental, and pathophysiologic factors.14–16A population approach of these data may be attractive because it can explain a part of the wide variability by incorporating covariates, such as the type of local anesthetic and age. In addition, it enables the development of models that are able to make predictions about certain clinically important endpoints, such as maximal level and duration of blockade. This may, consequently, improve therapeutic outcome in future patients.6
Pharmacokinetic Analysis
The current analysis yielded similar results with respect to our previous studies; small differences, however, were apparent. For example, we now observed that the clearances of the two anesthetics decrease significantly with increasing age.
The observed differences between the absorption parameters of the two anesthetics can be explained from their vasoconstrictory and vasodilatory properties.5The effects of those differences on the concentration in the epidural space after a standardized bolus dose are shown in figure 7. Levobupivacaine concentration levels are the highest, and this anesthetic could be called the most efficient in terms of concentration. This observation has to be linked to the pharmacodynamic potency to determine whether it is also effect efficient (see Discussion, Pharmacodynamic Analysis section).
Fig. 7. Comparison of absorption profiles of levobupivacaine and ropivacaine in the epidural space (absorption parameters from table 4).
Fig. 7. Comparison of absorption profiles of levobupivacaine and ropivacaine in the epidural space (absorption parameters from table 4).
Parameter F , which represents the fraction of the administered dose that is present in the patient’s blood, should be 1. It was greater than 1 for levobupivacaine and ropivacaine. Possible explanations are given in our previous articles3–5; in brief, bioavailability may be overestimated because of the finite experiment time from which the disposition kinetics (in particular elimination clearance) are determined.
Pharmacodynamic Analysis
We observed that the sensitivity to the two anesthetics, quantified by an increase of parameter A thr, decreased with segment height. Assuming that the sensitivity of the nerves to the presence of an anesthetic, and the volume of the central absorption compartments do no differ between segments, this decrease is due to a smaller amount of anesthetic reaching the higher segments. The overall A thris highest for levobupivacaine. So although the absorption characteristics cause this anesthetic to be concentration efficient, this is counteracted by its lower sensitivity (figs. 3 and 7). Ropivacaine was found to have a lower speed of onset and offset (as expressed by t1/2,ke0 ) than levobupivacaine. This may be explained by the fact that the physiochemical properties of ropivacaine are similar to those of levobupivacaine, apart from its lower lipid solubility. Because lipid solubility is related to potency, ropivacaine may be less potent.17The smaller t1/2,ke0 in the neighborhood of the L2 segment is probably related to differences in the distribution of the anesthetics in the epidural space close to the site of injection.
Anatomical and physiologic changes associated with advancing age may affect the pharmacokinetics and the nerve block characteristics after epidural administration of local anesthetics. A declining number of neurons, deterioration in myelin sheaths in the dorsal and ventral roots, changes in the anatomy of the spine, and intervertebral foramina may contribute to altered nerve block characteristics after epidural anesthesia.18,19Furthermore, the number of axons in peripheral nerves decreases with advancing age, and the conduction velocity diminishes, particularly in motor nerves.19–21With increasing age, changes in the connective tissue ground substances and epidural fat may result in changes in local distribution, i.e. , in the distribution rate of the local anesthetic from the site of injection (the epidural space) to the sites of action.18We observed important age effects for ropivacaine and levobupivacaine, respectively: (1) A thrdecreased at the higher segment levels (T12 and higher); and (2) onset and offset of the sensory blockade, as expressed by parameter t1/2,ke0 was faster with increased age. The most probable mechanism of an increased anesthetic sensitivity with age is an age-dependent change in the longitudinal spread of the anesthetics in the epidural space. This is evident because we observed an increased sensitivity with age mainly at the higher segments. The reduced loss of anesthetic via sclerotic intervertebral foramina during its rostral spread with increasing age may explain this observed age effect. The higher level of analgesia in older patients may as well be attributed to a greater sensitivity (a lower A thr) such that with the same local anesthetic concentrations at higher thoracic segments, blockade occurs in older but not in younger patients. Consequently, to obtain comparable epidural blocks, smaller doses of local anesthetic solutions should be administered to older as compared with younger patients. The effect of age on the speed of onset and offset of effect is probably related to changes in epidural fat. Uptake into extraneural tissues, such as epidural fat, limits the rate and extent of drug distribution to the nerves and thereby changes the time profile of clinical potency.22
PK-PD Models of Spinal Anesthesia
To the best of our knowledge, this is the first PK-PD model developed for lumbar epidural administration of local anesthetics. There are, however, PK-PD models for spinal anesthesia in pigs and in humans.7,23,24Shafer et al. 23performed PK-PD modeling of intrathecal neostigmine using cerebrospinal fluid samples and analgesia (visual analog scale) scores. The pharmacokinetic part of their model included a component to account for the spatial dimension in the intrathecal administration–analgesia relation. The pharmacokinetic model developed by Ummenhofer et al. 24consists of four spatially interconnected subunits of four serially connected compartments to describe drug distribution in each subunit. From microdialysis data, it was not possible to estimate all model parameters for all subunits, and some additional assumptions had to be made. For example, exchange between spinal levels was assumed to occur only via the cerebrospinal fluid and not via the epidural space. For our model of epidural anesthesia, we assumed that no exchange occurs via the cerebrospinal fluid9and epidural space, except in the latter during epidural administration. Schnider et al. 7developed a population pharmacodynamic model of spinal anesthesia from pinprick assessments. Although their model is applicable to epidural anesthesia, it is unable to predict sensory blockade at levels lower than the maximum. Furthermore, because they did not take into account pharmacokinetic data (on disposition and absorption) and hence the local anesthetic concentration, their pharmacodynamic analysis is affected by pharmacokinetic variability. Although this seems of limited importance when treating and predicting epidural anesthesia in a new patient, the inclusion of covariates in the pharmacokinetic model will improve the prediction of sensory anesthesia in this particular patient.
Critique on Methods
A potential drawback of our model is the underlying essential assumption that rostral and caudal spread of the anesthetic in the epidural space is instantaneous and subsequently remains unchanged (apart from absorption). In reality, the local anesthetic spreads with a certain delay. Incorporation of segment-dependent delay is not feasible, because we have only two measurements per dermatome (onset and offset times of blockade). Note that in reality t1/2,ke0 is not the delay to the segments but the delay from the segment to the effect site (i.e. , the spinal nerve roots). Consequently, we may have slightly overestimated the value of t1/2,ke0 .
Furthermore, the results that we present here are valid for the specific volume of the local anesthetic given as well as the location of the epidural puncture (L3–L4 interspace). It has been shown that in contrast to dose, volume per se has little or no effect on sensory block height and quality of anesthesia.25,26We expect, however, that the site of injection affects the parameter values of our model.
The way covariate age was included in the models allows for a covariate effect in one way only (either a decrease or an increase). In our previous articles, patients were assigned to age groups; however, this allowed for physiologically unrealistic covariate effects.
There was no intraindividual variability in the pharmacodynamic data because of the nature of the binary assessments. There was, however, intraindividual uncertainty in the actual times of onset and offset of effect, depending on the sampling times according to the protocol; this uncertainty will bias the estimates of interindividual variability (ωAthr2and ωk e02) upward.
Conclusions
In conclusion, we developed a predictive PK-PD model of epidural anesthesia in humans. We were able to demonstrate the importance of age on the spread of the sensory blockade as well as on the speed of onset/offset of blockade. The model allows the study of various important factors in epidural anesthesia, such as the site of injection, the volume of the injected fluid, combined spinal–epidural injections, anesthetic-opioid interaction, and pregnancy. Furthermore, the model may be used in the development and study of new local anesthetics for epidural use.
References
Appendix
Likelihood Function for the Set of Observations
For a type 1 observation, pinprick assessment times yield intervals (t on,1, t on,2) and (t off,1, t off,2) that constitute the closest intervals around the actual onset and offset times of effect (fig. 8). The remaining assessments do not contain any additional information, because they were always 0 for t ≤t on,1and t ≥t off,2, and 1 for t on,2≤t ≤t off,1(a result of the choice of assessment times with respect to the rate of change in the effect site absorption profile). Therefore, the uncertainty about whether there is a blockade cannot be estimated and has to be neglected, but there is uncertainty with respect to the actual onset and offset times. This uncertainty is taken into account by the postulated normal probability distributions for the logarithms of A thrand k e0with variances ωAthr2and ωk e02, respectively. Thus, they are a measure of the combination of intraindividual and interindividual variabilities (which cannot be separated). Figure 9depicts the (A thr, k e0) area for which it is possible that the observed intervals (t on,1, t on,2) and (t off,1, t off,2) are observed. The likelihood of an observation in individual i is given by the (double) integral of the probability distributions of A thrand k e0across this area:
Fig. 8. Central and effect site absorption profiles in the epidural space ( solid and hatched curves , respectively) given by equations 1 and 4with parameter estimates from patient 4210 substituted. The dots denote the effect site absorption A e( t ; k e0) at onset and offset times t on,1, t on,2 , t off,1, and t off,2, respectively. The anesthetic sensitivity A thr( horizontal line ) and central effect site equilibration rate constant k e0were approximated from the average ( t on , t off) (see text).
Fig. 8. Central and effect site absorption profiles in the epidural space ( solid and hatched curves , respectively) given by equations 1 and 4with parameter estimates from patient 4210 substituted. The dots denote the effect site absorption A e( t ; k e0) at onset and offset times t on,1, t on,2 , t off,1, and t off,2, respectively. The anesthetic sensitivity A thr( horizontal line ) and central effect site equilibration rate constant k e0were approximated from the average ( t on , t off) (see text).
Fig. 9. The area of integration is the area bounded by the effect site absorption profiles for t on,1 , t on,2 , t off,1, and t off,2(when observed) as functions of k e0(given by equation 4with fixed absorption parameters). The dot in the center denotes the ( A thr , k e0) estimate from the average ( t on , t off).
Fig. 9. The area of integration is the area bounded by the effect site absorption profiles for t on,1 , t on,2 , t off,1, and t off,2(when observed) as functions of k e0(given by equation 4with fixed absorption parameters). The dot in the center denotes the ( A thr , k e0) estimate from the average ( t on , t off).
where NL denotes the fact that A thrand k e0are normally distributed in the log-domain and μ and ω2denote, in NONMEM parlance, typical population values and variabilities. The integration interval for A thrdepends on k e0, and the integration interval for k e0is determined by nonempty A thrintegration intervals (fig. 9). The integrals are evaluated in the log-domain; the inner integral can be numerically approximated by one of the GSL gsl_cdf_normal routines, and the outer integral can be numerically approximated by one of the GSL gsl_integration_qag routines, depending on whether the k e0range is finite or infinite. The likelihood of the set of observations is given by the product of equation 7for each individual.
For a type 2 observation, a blockade is still present at the last assessment time, so A e(t end) =A e(t off,1) ≥A thrand t off,2is not observed, which means that the information about A thrand k e0is not bounded by the line A e(t off,2, k e0) in figure 9, and the area of integration therefore also includes the horizontally shaded part (with k e0↓ 0).
For a type 3 observation, A e(t ) < A thrfor all assessment times; only the maximum value of A e(t max) gives information, because if A e(t max) < A thr, A e(t ) < A thrfor all t ≠t max(the observations are not independent, and function A e(t ) has only one maximum).
The implementation of the integrals in the C program were checked using Monte Carlo simulation, i.e. , by counting the number of samples from the lognormal distributions of A thrand k e0that fall within the region of integration.
The model parameters were estimated by the method of maximum likelihood, i.e. , they are given by those values that maximize L =∏iL i(for all individuals i , equation 6); the maximum was numerically determined by the GSL routine gsl_multimin_fminimizer_nmsimplex.
Initial Parameter Values
Suppose that t on,1≈t on,2and t off,1≈t off,2, so that they can be approximated by their averages t onand t off. Then there exists a unique relation between (t on, t off) and (A thr, k e0); the latter are given by those values that satisfy A e(t on; k e0) =A e(t off; k e0) ≡A thr(using equation 4). They are unique because A e(t on; k e0)↑ when k e0↑ while A e(t off; k e0)↓ when k e0↑. The equality is numerically found by the GSL root-finding routine gsl_root_fsolver_brent. Initial parameter values were calculated as the means and variances of the logarithms of the A thrand k e0estimates at each dermatome for all individuals from the type 1 observations.
Blockade Probabilities
Once the distributions of A thrand k e0are known (their parameters estimated), the probability of blockade at time t is given by:
where the population values μAthrand μke0are functions of age as given by equation 5, and the parameters of A e(t ; k e0) (equation 4) are given by population estimates of F 1, F 2, k a1, and k a2, which are also functions of age. When the population distributions of the absorption parameters are postulated correctly and/or their variances are small, the bias caused by taking their population estimates instead of in addition integrating across their respective distributions is assumed to be negligible.
