Semimechanistic models to relate noxious stimulation, movement, and pupillary dilation responses in the presence of opioids.

Intraoperative targeting of the analgesic effect still lacks an optimal solution. Opioids are currently the main drug used to achieve antinociception, and although underdosing can lead to an increased stress response, overdose can also lead to undesirable adverse effects. To better understand how to achieve the optimal analgesic effect of opioids, we studied the influence of remifentanil on the pupillary reflex dilation (PRD) and its relationship with the reflex movement response to a standardized noxious stimulus. The main objective was to generate population pharmacodynamic models relating remifentanil predicted concentrations to movement and to pupillary dilation during general anesthesia. A total of 78 patients undergoing gynecological surgery under general anesthesia were recruited for the study. PRD and movement response to a tetanic stimulus were measured multiple times before and after surgery. We used nonlinear mixed effects modeling to generate a population pharmacodynamic model to describe both the time profiles of PRD and movement responses to noxious stimulation. Our model demonstrated that movement and PRD are equally depressed by remifentanil. Using the developed model, we changed the intensity of stimulation and simulated remifentanil predicted concentrations maximizing the probability of absence of movement response. An estimated effect site concentration of 2 ng/ml of remifentanil was found to inhibit movement to a tetanic stimulation with a probability of 81%.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

Movement is a known response to noxious stimulus. However, its applicability to assess intraoperative nociception is limited when neuromuscular blocking agents are used. Pupillary reflex dilation (PRD) is known to predict movement after noxious stimuli and to be useful in assessing nociception response.

WHAT QUESTION DID THIS STUDY ADDRESS?

This study aimed to quantitatively relate PRD and movement with noxious stimulation and opioids to gain an understanding of the underlying mechanisms that control both variables.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

This study provides a granular description of PRD and movement reflex to tetanic stimulation, which were characterized with semimechanistic population pharmacodynamic models. PRD and movement showed similar responses to painful stimulation and remifentanil, supporting the hypothesis that both reflexes share a common physiological mechanism of control.

HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?

This study provides evidence for the use of PRD as a surrogate marker for pain, guiding opioid administration during surgery.

INTRODUCTION

Titration of opioids during surgery must be performed according to the intensity of external stimuli to attain the right nociceptive–antinociceptive balance. The estimation of this balance depends on reflex quantification of responses to stimulation and nonspecific autonomic responses as well as previous clinical experience. The most direct and clinically observed reflex response to noxious stimulation is movement. Other indirect signs derived from the autonomic nervous system stimulation such as arterial hypertension, tachycardia, sweating, or presence of tears have also been clinically used as surrogate measures of response to noxious stimulation. These indirect signs can be affected by changes that are not related to painful stimulation. Therefore, the pursuit for the individualization of the analgesic component of general anesthesia ,  lead to the development and commercialization of different nociception indexes/monitors. One such monitor quantitates pupillary reflex dilation (PRD) to a homogenous electrical painful stimulation. PRD is a supra‐spinal parasympathetic reflex present during general anesthesia , shown to be triggered by noxious stimuli in anesthetized subjects. , This reflex is diminished by opioid administration,  suggesting it could be used for the assessment of analgesia in patients under general anesthesia.

Reflex movement to a noxious stimulus is considered to be a classical response, and its generation involves the transmission of nociceptive signals from the periphery to the central nervous system and includes the modulation by the rostral ventral medulla (RVM), where “on‐cells” and “off‐cells” are present. In the close vicinity of these areas, the control mechanisms of the PRD are located. Some authors have shown that the effect of externally administered opioids can affect the normal function of “off‐cells,” preventing the noxious stimuli to progress toward the cortex and hence blocking any movement response. , Similarly, the presence of opioids prevents the dilation mechanism of the pupil from working.

To better understand the effects of opioids on PRD and its relationship with the reflex movement response to nociception, we conducted a prospective clinical study where patients under the effects of propofol received varying target concentrations of the opioid remifentanil to observe and quantitate the movement and PRD responses elicited by a homogenous painful electrical stimulation. The main objective was to generate a population model relating remifentanil concentration to movement and to pupillary dilation. Based on the developed model, we simulated possible opioid dosing schemes, maximizing the probability of the absence of movement response, that is, cortical repercussion of pain, in unconscious patients under general anesthesia.

METHODS

Patient population and characteristics

Under Institutional Review Board and Ethics Committee approval (Hospital Clinic de Barcelona no. HCB/2016/0318) and following written informed consent, 78 female patients scheduled for gynecologic surgery under general anesthesia in the Ambulatory Surgery facility at Hospital Clinic de Barcelona, Spain, were included in the current study. Exclusion criteria were prior eye surgery, any ophthalmologic diseases besides refraction errors, prescription of drugs affecting the size or reflex of the pupil, and morbid obesity (body mass index > 35).

Study protocol

Upon arrival in the operating room, routine monitors were placed. General anesthesia was achieved and maintained with propofol and remifentanil administered using a Target‐Controlled Infusion system (Base Primea, Fresenius Kabi AG, Bad Homburg, Germany).

Loss of consciousness was induced setting the effect site concentration of propofol (Ce,Prop) between 5 and 11 µg/ml (Schnider et al. , ).

PRD measurements and reflex movement assessment

Two minutes after reaching the predicted pseudo–steady state for propofol, the PRD was elicited using the videopupillometer Algiscan® (IDMed, Marseille, France) connected to two electrodes placed on the volar surface of the right arm. The stimulus consisted of a 60 mA tetanus for a period of 5 s. Pupillary diameter (PDiam) was measured and recorded at 67 Hz for 13 s: 2 s before the stimulation, during the 5 s stimulation, and until 6 s after the stimulation ended. A rubber cup covered the measured eye, and the contralateral eye was taped closed.

The investigational procedure was video recorded. Retrospectively, two clinicians examined the video recordings to evaluate intensity and duration of the movement response. Movement intensity was classified into the following four‐level categorical scale: 0 = absence of movement, 1 = movement of only one limb, 2 = movement of two or three limbs, and 3 = movement of the whole body.

Within approximately 5 min after the first stimulation, while maintaining the same propofol effect site concentration, remifentanil administration was started, achieving predicted effect site concentrations (Ce,remi) varying between 0.5 and 6 ng/ml. Two minutes after reaching the predicted pseudo–steady state equilibration, a second and a third PRDs were elicited as described previously, and PDiam and movement response were recorded.

Then airway was then secured either by placement of a laryngeal mask or by endotracheal intubation. When intubation was required, 30 mg of rocuronium bromide were administered, two minutes before laryngoscopy. The hypnotic effect was titrated to maintain a Bispectral Index (BIS Vista, Medtronic, Ireland) value between 45 and 60. The analgesic effect was titrated at the discretion of the attending anesthesiologist.

At the end of surgery Ce,prop and Ce,remi were set to 1 µg/ml and 0 ng/ml, respectively. When the propofol Ce target was reached, one or two additional stimulations separated by at least 3 min were carried out before the patient spontaneously awakened or moved.

General description of the data

Figure 1 shows the raw data profiles for pupillary size, movement response, and drug concentrations used for the model development. Of all pupil measurements, 1082 were recorded in the presence of propofol (first stimulus), and 2775 were recorded with a combination of propofol and remifentanil (stimuli 2–4). For the movement analysis, 1092 observations correspond to the first stimulus and 2816 to stimuli 2–4. In both end points, the number of observations after surgery was lower than for those performed before surgery started as patients were recovering consciousness after the surgical procedure, and measures were interrupted when patients showed signs of arousal. Predicted propofol and remifentanil concentrations were recorded every second during the experimental procedure.

FIGURE 1

Overview of the experimental set‐up and data generated during the experimental procedure. Upper panels show the median concentrations of propofol (red) and remifentanil (green) over the four cycles of stimulation. Middle panels depict the raw pupil profiles for each stimulus. Orange line highlights the period of 5 s when electric stimulation is delivered. Lower panels represent the distribution of the movement scores per stimulus (Grade 0, red; Grade 1, blue; Grade 2, green; Grade 3, yellow). The bottom figure represents the timeline of the experimental procedure with the administration of propofol (red syringe), electric stimuli (yellow lightning), administration of remifentanil (green syringe), and surgical procedure. Times represent the mean starting time for each stimulus among all the individuals. Images were modified from Servier Medical Art and Freepik

Data analysis

PRD and movement data were fit separately. All of the analyses were carried out using the population approach with nonlinear mixed effects modeling and the software NONMEM 7.4. The first‐order conditional estimation method with interaction was used for the analysis of the PDiam (continuous data). The movement response was treated as an ordered categorical variable and analyzed using the Laplacian estimation method instead.

Interindividual variability (IIV) was modeled exponentially for parameters integrating the model of the PDiam. On the other hand, an additive model was used to describe IIV associated with the Logit (L) analysis of movement scores because regardless of the magnitude of the individual random effect, probabilities where kept between 0 and 1. PDiam data were logarithmically transformed for the analysis, and the residual error was modeled initially with an additive model in the logarithmic scale. Given the fact that PDiam were measured very frequently, autocorrelation should not be rule out, and a model for autocorrelated errors was also explored in the analysis as indicated by an anonymous reviewer. Residual variability does not apply in the case of the analysis of noncontinuous response.

Pharmacodynamic models

The raw data profiles presented in Figure 1 show two features that are relevant from a modeling point of view: (i) there is a delay between the beginning of the nociceptive electric impulse and the onset of the pupil and movement responses, and (ii) remifentanil reduces both the PDiam and the probability of movement. Consequently, models based on the concept of the indirect response, and drug inhibitory effects either on the perception and/or transduction of the nociceptive signal, were fit to the data.

Models assume that the degree of activation of unobserved levels of nociceptors (N act) drive the pupil and movement response. In absence of either electric stimulation or remifentanil, N act was arbitrarily set to a baseline value of 1 (N act,0), maintained by the balance between the zero and first‐order rate constants of turnover, K S and K D, respectively, as shown by Equation (1). At baseline, the rate of change of N act (dN act/dt) is null, and K S equal K D as N act,0 = 1.

Electric stimulation increases N act as indicated in Equation (2), where θ ES is the parameter accounting for its nociceptive effects, and ES is a variable with a value of 1 during the application of the electrical stimulus and 0 otherwise.

Analgesic effects of remifentanil were explored as a modulation of θ ES, and/or increasing K D, the first‐order rate constant of nociceptor deactivation, and are expressed in Equation (2) by f(Ce ,remi), and g(Ce,remi), respectively.

During model development, different structures for f() and g() considering the predicted remifentanil concentration in plasma (C p,remi) and effect site (C e,remi) were explored as shown in the Results section. Predicted values of C e,remi were obtained using Equation (3) representing the effect compartment model. Contribution of propofol to the activation of nociceptors was also investigated. where k e0, is the first‐order rate constant governing the drug distribution equilibrium between the central and effect site compartments.

PRD response

A turnover model was used to transduce the N act dynamics to PRD as shown in Equation (4), where K P,S and K P,D represent the zero‐order and first‐order rate turnover constants, respectively, being K P,S = K P,D × PDiam0, and PDiam0, the value of PDiam at baseline.

Additional effects of remifentanil on the turnover parameters were explored, in addition to the possibility that differences in PDiam and N act dynamics were negligible. The contribution of propofol to the activation of nociceptors and pupil turnover was also investigated. Linear and nonlinear inhibitory models as well as an effect site compartment were evaluated.

Reflex movement response

Movement response was recorded as an ordered categorical variable and was modeled using the proportional odds logistic regression approach. The cumulative conditional probability of observing a score greater than or equal to a certain m category in the ith subject at the jth observation time is denoted P(Y ≥ m|η ) and is described by Equation (5), where η is the individual random effect that belongs to a distribution of mean 0 and variance ω2.

L, the logit, brings together the parameters defining the baseline probabilities of movement, and the contribution of N act represented in Equation (6) as θ and h(N act), respectively.

Finally, the probability of observing an m score P(Y = m|η) is given by the following expression (Equation 7):

The contribution of first‐order Markov elements to the movement response was explored as well as considering the nonproportional odds approach and the inclusion of an additional turnover compartment as in the case of PDiam.

Model selection

The minimum value of the objective function, approximately equal to −2 × log‐likelihood (−2LL), was used to select between candidates during model building through the log‐likelihood ratio test (LRT). A decrease of 3.84 or 6.61 in the −2LL—corresponding to a 5% or 1% level of significance, respectively—was considered statistically significant for nested models differing in one parameter. For comparison across non‐nested models, the Akaike information criteria (AIC) was used. In addition to LRT or AIC, precision of parameter estimates expressed as the percentage of relative standard error (RSE[%])—calculated as the ratio (multiplied by 100) between the standard error of the parameter and the corresponding estimate—was also considered, together with the visual inspection of the goodness‐of‐fit (GOF) plots.

Model evaluation

Once the models for PDiam and movement responses were selected, the evaluation process was performed. Parameter precision was further evaluated by analyzing 500 nonparametric bootstrap data sets using the selected models and calculating the median and the 95% confidence intervals for each of the model parameters. Visual predictive checks (VPCs) corresponding to simulation‐based diagnostics were generated to evaluate model performance. A total of 1000 data sets of the same characteristics as the original were simulated using the structure and parameters from the selected models. For each simulated data set and each measurement time, either the 2.5th, 50th, or 97.5th percentiles of the PDiam or the percentage of movement scores of 0, 1, 2, and 3 were calculated. Then the areas covering the 95% prediction intervals of each of the simulated percentiles or each of the movement scores were generated and plotted together with the corresponding raw data.

Data corresponding to the fifth stimulus (the second after surgery) was not used for model development but as an internal validation data set. The model developed for stimuli 1–4 was used to simulate the PDiam and movement grades during the fifth stimulus, and this simulation confroned with the observations. In addition, in the case of pupil, individual model parameters were used to predict the PDiam during the fifth stimulus. Prediction errors (PE)—calculated with Equation (8)—were computed to quantify the precision of the individual predictions for the fifth stimulus and compared with the PE for the previous stimuli. where Pred and Obs , refer to the ith individual model predictions and observations obtained at time jth, respectively.

Software and tools

R Version 3.6.1 with RStudio interface (Version 1.2.5001) were used for building the NONMEM data set and generating the graphical output. GOF plots were generated with the R packages Xpose and ggplot2. VPCs, covariate evaluation, and bootstrap simulations were performed with Perl‐Speaks‐NONMEM. ,

RESULTS

Demographic characteristics of the patients included in the study are shown in Table 1. Supplementary material (SM) 1‐Figure 1 shows the distribution of the different combinations of propofol and remifentanil targeted during the study.

TABLE 1

Patient characteristics

Patients, n	78
Female, n	78
Age, years, median (range)	45 (27–85)
Height, cm, median (range)	160 (140–173)
Weight, kg, median (range)	64 (38–93)
Lean body mass, kg, median (range)	44 (31–58)
Body surface area, m2, median (range)	1.65 (1.26–2.04)

Model development

Pupillary reflex dilation

A model considering that the dynamics of the PRD are influenced by the time course of N act provided better fits compared with just one turnover compartment (∆‐2LL = −890; p < 0.01). In Figure 2a, the schematic representation of two compartments required to properly describe the PDiam observations is shown.

FIGURE 2

Schematic and mathematical representation of the pharmacokinetic/pharmacodynamics models developed for pupil diameter (a) and movement response (b). All parameters are described in the text. Note that in our model, KS and KD are the same parameter

Remifentanil showed significant (p < 0.01) effects in (i) reducing of the nociception activation triggered by the application of the electric stimulus and (ii) increasing the first‐order rate constant K D, as represented in Equation (9): where ES is a variable with a value of 1 during the application of the electrical stimulus and 0 otherwise. θ ES is the parameter accounting for the nociceptive effects of the ES, and C50 is the concentration of remifentanil (C p,remi) reducing the impact of θ ES to half its maximum value. These terms account for the induction of the PRD response by the electrical stimulation and the attenuation by remifentanil. The term θ remi refers to the slope at which the concentrations of remifentanil in the effect site increases K D and is related to the decrease of PDiam over time after administration of remifentanil. To prevent that remifentanil could decrease PDiam below physiological constraints, the term (PDiam‐θ Limit/PDiam) was incorporated into this effect. Finally, the contribution of propofol was analyzed and did not significantly modified the different processes reflected in the model (p > 0.05) (SM 1‐Figure 2). Figure 2a shows both the schematic and the full mathematical representation of the structural part of the population model for PDiam.

The model parameters, shown in Table 2, were estimated with good precision. The estimated PDiam at baseline was 3.8 mm. The estimates of the first‐order rate constants K D, K P,D, and k e0 corresponds to half‐lives of 9.5 s, 1.3 s, and 3.56 min, respectively, indicating faster turnovers of N act and PDiam than distribution of remifentanil from plasma to the effect site. The steady‐state concentration of remifentanil to reduce by half the impact of θ ES and double K D are 1.15 and 0.67 ng/ml, respectively.

TABLE 2

Model parameter estimates

Parameter	Estimate	RSE (%)	2.5th–97.5th	Shrinkage, %
Pupillary reflex dilation response
K S (s−1)	0.0726	19.6	0.0203–0.152	–
PDiam0 (mm)	3.80	3.3	3.57–4.02	–
IIV PDiam0 (%)	24.6	7.7	20.9–28.4	1
K P,D (s−1)	0.531	13.9	0.294–1.07	–
IIV K P,D (%)	54.4	22.0	31.5–107	26
θ ES	1.54	11.0	1.10–4.33	–
IIV θ ES (%)	49.6	11.5	36.9–63.7	9
Cov (ω2 PDiam0, ω2 θ ES)	−0.0574	18.7	−0.0978 to −0.0236	–
C50 (ng/ml)	1.15	14.0	0.868–1.52	–
θ remi (ml/ng)	1.50	12.9	1.10–2.39	–
IIV θ remi (%)	55.6	8.7	44.1–73.9	2
k e0 (s−1)	0.00324	8.0	(2.67–3.82) × 10−3	–
θ Limit (mm)	1.41	6.1	1.18–1.61	–
Residual error (%)	9.90	6.6	8.50–11.0	3
Movement response
K S (s−1)	0.437	13.9	0.336–0.617	–
θ ES	6.63	8.2	5.67–7.98	–
θ 1	−4.90	4.0	−5.38–−4.58	–
θ 2	−0.754	11.2	−0.916–−0.584	–
θ 3	−0.978	17.1	−1.32–−0.697	–
θ 1,1	1.76	13.1	1.30–2.17	–
θ 2,2	1.93	16.2	1.31–2.63	–
θ 3,3	2.43	16.7	1.58–3.50	–
C50 (ng/ml)	0.617	23.3	0.360–1.02	–

IIV is expressed as percentage coefficient of variation calculated as , where ω2 corresponds to the variance of the random effects. The rest of the terms are defined in the main text.

Abbreviations: Cov, covariance; RSE(%), percentage of relative standard error.

Data supported the estimation of moderate IIV on PDiam0, θ ES, K P,D, and θ remi ranging from 24% to 56% coefficient of variation. A correlation of −0.51 was found between random effects associated with PDiam0 and θ ES. As no tendencies were detected between individual estimates of model parameters and patient characteristics, no further investigation of potential covariate effects was undertaken. Despite that including autocorrelation in the residual model decreased the −2LL, the model performed worse as reflected by the individual predictions and VPCs, and therefore autocorrelation was not included in the selected model.

Figure 3 shows the GOF plots (a) and the results of the simulation‐based diagnostics (b) corresponding to all data used for model development, both indicating adequate model performance and an absence of model misspecifications. Similar results were found during internal validation using the data corresponding to the fifth stimulus, which are shown in SM 1‐Figure 3. SM 2 shows the NMTRAN code corresponding to the selected model.

FIGURE 3

Pupil diameter model evaluation. (a) Goodness‐of‐fit plots. Data are colored according to the stimulus: 1, red; 2, green; 3, blue; 4, purple. Dashed line represents the loess smoothing curve. (b) Visual predictive checks for Stimulus 1 to Stimulus 4. Median (solid line) and 2.5th and 97.5th percentiles (dashed lines) of observed data compared with 95% prediction intervals (shaded area) for the median and 2.5th and 97.5th percentiles based on 1000 simulations. CWRES, conditional weighted residuals; IWRES, individual weighted residuals

Reflex movement

In the movement response, the addition of the second turnover compartment did not improve the fit significantly (p > 0.05). Remifentanil elicited a significant effect (p < 0.01), reducing the nociceptive effects triggered by electric stimulation as it was previously described for PDiam. Figure 2b shows the schematic and mathematical representation of the model for movement reflex. Regarding the structure of the logit, the function h() (see Equation 6) has as argument N act – 1, which gives the value of zero in unperturbed conditions. The addition of a random effect on the L was not significant (p > 0.05), and as in the case of the PDiam, propofol was not shown to significantly contribute to antinociception.

The possibility that the probability of each score at a certain time is influenced at least in part on the score achieved in the previous time (first‐order Markov process) was also investigated. Addition of the parameters accounting for all transitions resulted in a significant (p < 0.01) improvement in the predictions, particularly for movement Grades 1 and 3 (SM 1‐Figure 4, Panel c).

FIGURE 4

Movement response model evaluation. (a) Observed versus predicted probability of movement for each grade and number of stimulus calculated from 500 simulations. Solid line is the identity representing perfect fit. (b) Visual predictive checks stratified by grade. Solid line corresponds to the median probability calculated from raw data, and the shaded areas represented the 95% prediction intervals computed from 1000 simulations

Table 2 also lists the model parameter estimates corresponding to the movement response. All parameters were estimated precisely as confirmed by the low RSE and 95% confidence intervals. The estimate of K D corresponds to a value of half‐life of 1.58 s, and C50 is estimated as 0.617 ng/ml. In absence of stimulus, the probability of transitioning to movement is negligible (<0.8%). Once a certain grade is reached, the probability of maintaining the movement intensity is 85%, 87%, and 92% for Grades 1, 2, and 3 and returning probabilities to baseline are 15%, 13%, and 8%, respectively.

Results shown in Figure 4, corresponding to GOF (left) and VPCs pooling all data together (right), indicate that the model described the data properly. SM 3 shows the NMTRAN code corresponding to the selected model incorporating the Markov elements.

SM 1‐Figure 5 splits the VPC per stimulus and shows, in addition, good model performance for the movement responses gathered in stimulus 5, which were used for internal validation purposes. Our model was also able to capture the number of transitions between movement grades observed in the study: 261 for the raw data and a median of 272 for the model predictions (248 to 297, 95% confidence interval).

FIGURE 5

Models exploration. (a) Dynamics of nociceptor activation. (b) Pupil diameter and probability of movement in response to tetanic stimulation in absence (blue) or presence of remifentanil at 2 ng/ml steady‐state concentration (purple). (c) Percentage change in pupil diameter and probability of movement after tetanic stimulation at different steady‐state concentrations of remifentanil

Model exploration

SM 1‐Figure 4 explores the impact of some of the key components of the selected models on the measured responses. Panels A and B show the changes in pupil response after the administration of remifentanil and graphically explores the individual contribution of the model effects to these changes. Figure 5a shows the time course of the predicted degree of nociceptor activation for movement and PRD, and despite similar trends, movement appears to be associated to a faster turnover. The relationship between the amplitude in the pupil response profile and the probability of movement versus steady‐state remifentanil concentrations are shown in Figure 5b,c. Similar patterns between the two variables can be observed again, with the degree of steepness higher for movement.

Clinical applicability

Model‐based simulations were performed to identify the concentration of remifentanil at which 80% of patients would not experience movement after the delivery of a tetanic current ranging from 1.2 to 0.8 times the one used in the experimental procedure (60 mA). The simulation assumed an intensity of 1.2 for laryngoscopy and intubation while an intensity of 0.8 was assumed to represent a less painful surgical wound closure. For this purpose, we simulated 1000 individuals receiving electric stimulus of varying intensity at steady‐state concentrations of remifentanil from 0 to 4 ng/ml and computed the percentage of individuals not experiencing movement. A steady‐state concentration of 2 ng/ml was found to inhibit movement in 81% percent of the individuals with the experimental stimulus intensity and 77.8% and 83.1% of the individuals with a 20% increase and reduction in the stimulus intensity, respectively.

This concentration (2 ng/ml steady state) was subsequently used to simulate the PDiam of 1000 individuals from baseline to the end of a PRD. PDiam before electric stimulation for 80% of the individuals ranged from 1.69 to 2.61 mm (10th and 90th percentiles) and were distinctly different from the PDiam at baseline (2.71 and 5.25 mm, 10th and 90th percentiles, respectively) due to the pharmacologic effects of remifentanil. However, no correlation could be established between PRD and the prestimulus or basal (PDiam0) pupil, thus preventing recommendations on treatment individualization based on pre‐stimulus pupil size.

DISCUSSION

We have described both the PRD and reflex movement responses to noxious stimulation using a semimechanistic population pharmacokinetic/pharmacodynamics model to investigate similarities between the underlying mechanisms governing these two responses and the effects of remifentanil. Using this approach, we demonstrated that movement and PRD are equally depressed following remifentanil infusions. These data suggest that pupillary dilation and movement are initiated from a common source which is depressed by opioids.

The afferent noxious stimulus activates several sites within the pain matrix that in turn directly or indirectly inhibits the Edinger–Westphal nucleus and dilate the pupil. The RVM, gigantocellular reticular nucleus, and the locus coeruleus are possible sources of this inhibition. ,

These same nuclei have descending projections into the spinal cord that influence the withdrawal reflex. “Off‐cells” in the RVM for example are antinociceptive and are inhibited by noxious stimuli. , This inhibition of “off‐cells” is blocked by opioids. We hypothesize that the RVM and the EW nucleus are affected by nociception and opioids in a similar manner. Which after further investigation and confirmatory studies, could thereby allow the anesthesia provider to observe PRD as a measure of opioid effect on suppression of movement.

The final models selected in the current investigation have been carefully and extensively evaluated using both GOF plots and simulation‐based diagnostics while being validated internally with a set of data not used during model building. Data were accurately described, and all parameters were estimated with good precision. The models were capable of handling data gathered under different experimental conditions such as (i) the absence or presence of the analgesic drug and (ii) before and after surgery, indicating model robustness to support the conclusions of this work.

The use of changes in pupillary parameters as a surrogate marker of analgesic effects has received considerable attention in the past. PRD was shown to be useful to assess the nociception/antinociception balance and titrate opioid administration as PRD is triggered by noxious stimuli during general anesthesia , and diminished by opioid administration. Our results support the observations by Barvais et al. for the correlation between remifentanil Ce and percentage PRD as shown in Figure 5.

Pharmacokinetic/pharmacodynamics modeling has been previously done for the pupillary light reflex, demonstrating the effects of opioid drugs, , , , , and to a lesser extent, drugs modulating noradrenergic pathways. , , However, there are no previous modeling studies for PRD and how nociception and opioids influence it. Although only using isolated PDiam measurements, instead of the light reflex, Skarke et al.  modeled independently pupil and pain tolerance in response to electric pain stimulation and showed that (i) the effect site equilibration half‐life (t1/2_ke0) obtained for pupil was almost identical to the one from pain tolerance data and (ii) changes in both variables occur in parallel. Recently, Mangas‐Sanjuan et al.  modeled the PDiam effects of the active components of axomadol and found that pupillary changes correlated linearly with the area under the curve of the cold pressor test.

Pharmacodynamics of remifentanil is represented by values of C50 ranging from 0.92 to 11.2 ng/ml for other indicators of central activity such as arterial pressure, , spectral edge frequency of the electroencephalogram, and respiratory depression. , Remarkably, the corresponding estimates obtained in the current evaluation are in the same order of magnitude as those reported previously. Onset of action of remifentanil is fast with values of t1/2_ke0 of 0.53 to 2.48 min. , , , , The value of t1/2_ke0 found in this study is 3.5 min; however, for the additional effect of K D, the model considers an instantaneous equilibrium between the central and the effect site compartment. Those discrepancy might be caused by the fact that we predicted but did not measure plasma concentrations. Alternatively, it might well be that the region in the central nervous system in which remifentanil elicits analgesia is, at least from a distribution point of view, different from those contributing to cortical effects and respiratory depression.

This study has the objective to evaluate whether PDiam and movement responses ensue from a nociceptive stimulus and if the counteraction of remifentanil shares a common underlying pathway. In this context, one important aspect that deserves consideration is the granularity at which we measured the responses. The PRD was measured 67 times per second and the movement response each second during and after electrical stimulation in several occasions, resulting in an appropriate experimental framework to achieve our objective. This is the first pharmacodynamic model focusing on the pupillary effect of opioids ever to use this time resolution.

There are several possibilities to analyze the current data, for instance a simultaneous analysis where an unobserved variable (i.e., N act) is scaled differently for each of the two clinical end points. We decided to use an unbiased approach consisting on modeling the two end points independently hypothesizing a common core mechanism (the turnover model for N act) triggering a reflex movement response. The results from our analysis suggest that PRD and movement responses do share that mechanism as the model parameters are essentially similar for the turnover rates and remifentanil related effect parameters. However, the fact that both variables were not modeled simultaneously represents a limitation. The model indicates that differences between both responses appear up‐stream because the incorporation for PDiam of an additional turnover compartment, and for movement response of the Markov element contributing to maintain nociception, provided similar temporal patterns (see Figure 5b). Markov elements have been commonly identified in the analysis of noncontinuous responses (count or ordered as well as nonordered categorical) and in this evaluation contributed to improve the description of the data. ,

This is an innovative model that attempts to simultaneously describe the antinociceptive effect of remifentanil in the transduction and transmission component using the plasma concentration and the modulation phase through the effect site concentrations while describing how both influence PRD and movement in response to a noxious stimulus.

The knowledge provided by this investigation may allow a better understanding of the underlying phenomenon of nociception/antinociception during general anesthesia and allow for rational opioid dosing to prevent the central nervous system activation that would elicit a movement response and use the PRD as a surrogate reliable measure of adequate antinociception in the clinical setting.

CONFLICT OF INTEREST

II‐A is currently employed by AstraZeneca. CG and EWJ are employees of Quantium Medical. PLG has received honoraria from Quantium Medical. All other authors declared no competing interests for this work.

AUTHOR CONTRIBUTIONS

All authors wrote the manuscript. PLG and MDL designed the research. SV, MA, SJ, and PLG performed the research. NMA, II‐A, and IFT analyzed the data. AC, AP, CG, EWJ, PC‐M, and JFV contributed new reagents/analytical tools.

Supplementary material 1: Supplementary Figures 1 to 5

Click here for additional data file.

Supplementary material 2: NONMEM code for pupil model

Click here for additional data file.

Supplementary material 3: NONMEM code for movement model

Click here for additional data file.

Supplementary material 4: NONMEM dataset for one subject

Click here for additional data file.