Mechanism-based pharmacodynamic model for propofol haemodynamic effects in healthy volunteers

Background: The adverse haemodynamic effects of the intravenous anaesthetic propofol are well known, yet few empirical models have explored the doseeresponse relationship. Evidence suggests that hypotension during general anaesthesia is associated with postoperative mortality. We developed a mechanism-based model that quantitatively characterises the magnitude of propofol-induced haemodynamic effects during general anaesthesia. Methods: Mean arterial pressure (MAP), heart rate (HR) and pulse pressure (PP) measurements were available from 36 healthy volunteers who received propofol in a step-up and step-down fashion by target-controlled infusion using the Schnider pharmacokinetic model. A mechanistic pharmacodynamic model was explored based on the Snelder model. To benchmark the performance of this model, we developed empirical models for MAP, HR, and PP. Results: The mechanistic model consisted of three turnover equations representing total peripheral resistance (TPR), stroke volume (SV), and HR. Propofol-induced changes were implemented by Emax models on the zero-order production rates of the turnover equations for TPR and SV. The estimated 50% effective concentrations for propofol-induced changes in TPR and SV were 2.96 and 0.34 mg ml , respectively. The goodness-of-fit for the mechanism-based model was indistinguishable from the empirical models. Simulations showed that predictions from the mechanism-based model were similar to previously published MAP and HR observations. Conclusions: We developed a mechanism-based pharmacodynamic model for propofol-induced changes in MAP, TPR, SV, and HR as a potential approach for predicting haemodynamic alterations. Clinical trial registration: NCT02043938

volunteers that performed equally well as empirical models.
The mechanism-based model adequately predicted haemodynamic changes during propofol infusion in healthy volunteers, and the correlations between haemodynamic variables were well described by the model. This mechanism-based model will be useful in predictive algorithms for propofol-induced hypotension in healthy individuals, and potentially adaptable to patients with comorbidities.
Propofol, a positive modulator of g-aminobutyric acid receptors, is one of the most commonly used drugs for sedation and anaesthesia in clinical practice. The hypnotic and analgesic effects of propofol and the doseeexposureeresponse relationship have been extensively studied through pharmacokineticepharmacodynamic (PKePD) modelling. 1e3 The effects of propofol on changes in haemodynamic variables has also been studied after target-controlled infusion in healthy volunteers and patients. 4,5 The most pronounced propofol-induced haemodynamic effect is a decrease in sympathetic tone resulting in vasodilation and a decrease in total peripheral resistance (TPR), which leads to a decrease in MAP. 6 There is considerable concern that intraoperative hypotension is associated with postoperative mortality. 7 Although extensive descriptive knowledge on the association between exposure and the cardiovascular safety of propofol is available, only a few mathematical models have addressed changes in haemodynamic variables during propofol infusion in humans, and they all use empirical approaches 8e10 ignoring the complex interactions in the cardiovascular system.
In contrast to empirical models, mechanism-based PKePD models can distinguish between drug-specific parameters, describing the interaction between the drug and biological system, and biological system-specific parameters describing the functioning of the biological system such as feedback mechanisms. Mechanism-based models have better extrapolation properties, 11 and therefore better characterise propofolinduced haemodynamic side-effects.
Snelder and colleagues 12,13 have developed a PKePD model in rats that integrates a quantitative description of the physiology of the cardiovascular system and the effect of cardiovascular drugs on the relationship between MAP, TPR, HR, and stroke volume (SV). This model has been extended by including heart contractility data in dogs, 14 and recent models with similar structure perform reasonably well in humans. 15,16 We aimed to develop a mechanism-based PKePD model that describes the relationship between cardiovascular variables and propofol plasma concentrations using the principles outlined by Snelder and colleagues. 13 In addition, three empirical models for MAP, HR, and pulse pressure (PP) (serving as a surrogate of SV in this study) were developed to benchmark the performance of the mechanism-based model. Finally, we performed simulations for the final mechanismbased model predicting the expected haemodynamic changes according to US Food and Drug Administrationapproved drug labelling regulations and replicated published studies to evaluate the model.

Study design
The data used in this analysis were recorded as part of standard vital signs monitoring for safety purposes in a previous four-period randomised sequence crossover study as described. 17 Screening and inclusion methods were described by Kuizenga  RUGLOOPII software (Demed, Temse, Belgium) was used to control infusion rates between 0 and 1200 ml h À1 via an RS-232 interface. After 2 min of baseline measurements, a 'staircase' step-up followed by step-down infusion of propofol was initiated. Effect-site concentration targets were set at 0.5, 1, 1.5, 2.5, 3.5, 4.5, 6, and 7.5 mg ml À1 using a three-compartment pharmacokinetic model with an effect compartment developed by Schnider and colleagues. 18,19 Drug administration was stopped when burst suppression ratio by EEG reached 40%. At that point, propofol effect-site concentration targets were decreased with the inverse steps as performed during induction followed by a single bolus. Each target was maintained for 12 min to reach equilibrium, and arterial blood samples were taken at the end of the equilibration time as described. 17,20 Pharmacodynamic measurements and data handling Each volunteer was connected to a Philips IntelliVue MP50 vital signs monitor (Philips Medizin Systeme, Boeblingen, Germany) for ECG and noninvasive blood pressure measurements. Noninvasive systolic and diastolic blood pressure were measured every minute. HR was derived from the ECG and was recorded every second. All variables were collected using RUGLOOP II software. Mean arterial blood pressure and PP were calculated from the measured systolic and diastolic blood pressure.
Because of the length of the study sessions, our dataset contained more than 10 000 combined observations for HR for each individual. To reduce the computational burden during model development, we reduced the number of HR measurements per subject. Before data reduction, we applied a median filter to reduce the influence of artifacts or outlying data. The width (span) of the median filter was 30 s. Data reduction was performed by retaining only the first out of every 60 consecutive median filtered observations in the dataset.

Population PKePD modelling
The 'individual pharmacokinetic parameter approach' was applied to develop the PKePD models. 21 Plasma concentration measurements were used to derive individual post hoc pharmacokinetic parameters (CL, Q2, Q3, V1, V2, and V3) based on the propofol pharmacokinetic model published by Eleveld and colleagues, 1 and the pharmacokinetic parameters for each individual were fixed in the subsequent pharmacodynamic modelling. Goodness-of-fit plots for the Eleveld pharmacokinetic model for describing the propofol pharmacokinetics of this particular dataset is shown in Online Supplementary 3, Figs 3e1.
The model by Snelder and colleagues 13 was used as a starting point for development of the mechanism-based model. This model assumes that MAP is the product of cardiac output (CO) and TPR and that CO is the product of HR and SV, as in equation (1). 22 As SV and TPR were not measured directly, we used PP as a surrogate of SV, taking into account an SV/PP ratio of 1.5, 23e25 and TPR was implemented as a latent variable in line with contemporary pharmacometric approaches for handling longitudinally correlated data. 26 A latent variable is a variable that is not directly observed but its dynamics can be inferred (through the model) from other observable correlated variables (much like the hypothetical effect-site concentration that is used in effect-site PKePD models).
The Snelder model consists of three turnover equations for TPR, SV, and HR and a negative feedback mechanism between these variables, through MAP, representing the homeostatic feedback of the baroreflex system. A direct inverse relationship between HR and SV was also included, representing the effect of shorter left ventricular filling time because of increased HR thereby decreasing SV (expressed as HR_SV). 13 In line with the approach of Snelder and colleagues, we explored propofol effects on the turnover equations for TPR, SV, and HR. Propofol-induced increase or decrease in zeroorder production rate constants (k in ) or first-order dissipation rate constants (k out ) of the turnover equations were explored using E max and sigmoid E max models. The parameter that describes the relationship between HR and SV (HR_SV) was fixed to the value reported by Snelder and colleagues. 13 The focus of the development of the mechanism-based model was to construct a model with as few propofol drug effect parameters as possible, and with a goodness-of-fit comparable with the empirical models (as judged by goodness-of-fit plots and prediction-corrected visual predictive checks).
For the separate empirical pharmacodynamic models for MAP, HR, and SV, the hysteresis between predicted plasma concentrations and the pharmacodynamic measurements was described by an indirect response model. 27 The concentrationeeffect relationships were explored using E max and sigmoid E max models.
As MAP and HR were elevated before propofol infusion, we explored whether including a time-dependent effect in the models could improve the goodness-of-fit. The effect was estimated by an empirical function with two parameters k and LTDE; equations (2) and (3). TDE 0 represents the magnitude of the initial time-dependent effect. LTDE described the level of this effect influencing baseline MAP and baseline HR (expressed as percentage of increased baselines) with a firstorder dissipation rate constant k. Base HR represents baseline HR.
Once the basic model was established, the correlations between post hoc predicted parameters for which betweensubject variability (BSV) was included and covariates were analysed. Subsequently, the influence of these covariates was evaluated by inclusion in the models. The inclusion of covariates was accepted only if the objective function value decreased more than 3.85 points, which is the 5% significance level critical quintile of the corresponding c 2 distribution. The covariates considered were age, weight, height, sex, and BMI.

Parameter estimation and model evaluation
The MAP, PP, and HR observations were simultaneously fitted using the first-order conditional estimation method with interaction using NONMEM (version 7.4; Icon Development Solutions, Hannover, MD, USA). BSV was modelled using exponential or additive models. Residual variability was described using additive error models, proportional error models, or both. At each step of model building, the change in objective function value and the median absolute population prediction error were compared between candidate models. Goodness-of-fit was graphically evaluated by plotting individual or population predictions vs the observations and the conditionally weighted residuals (CWRES) vs population predictions and time. Parameter uncertainty of the models was estimated using the covariance step in NONMEM or sampling importance resampling. 35 Finally, prediction-corrected visual predictive checks were used to internally evaluate the models. 28 All models and simulations were run using PsN 29 and Pirana 30 as back end, front end, or both to NONMEM. Graphical assessment of the goodness-of-fit and simulations and the construction of the prediction-corrected visual predictive checks were conducted in R (R Foundation for Statistical Computing, Vienna, Austria). 31 Simulations for the mechanism-based model The final mechanism-based PKePD model was used to predict expected changes in haemodynamic variables during general anaesthesia for the dosing recommendations described in the propofol US Food and Drug Administration-approved drug label. 32 According to the label, adults (<55 yr old) received an induction dose of 2.5 mg kg À1 and a maintenance dose of 0.2 mg kg À1 min À1 for the first 15 min followed by a reduction of 40% (0.12 mg kg À1 min À1 ) up to 30 min followed by 0.1 mg kg À1 min À1 afterwards. For older individuals (>55 yr old), the induction dose was 1.5 mg kg À1 and maintenance 0.1 mg kg À1 min À1 . For the simulations, the distribution of covariates (age, weight, and sex) was according to the covariate distribution (including correlations between covariates) in our study population. We also simulated propofol concentrationeeffect relations for MAP, HR, and PP in a 35-yr-old individual.
We performed 1000 simulations to obtain 95% prediction intervals for MAP and HR to benchmark our model against the available literature. First, we replicated the study by Fairfield and colleagues, 33 which consisted of nine patients with a mean age of 23 yr (95% confidence interval [CI], 20e26 yr) and mean weight of 67 kg (95% CI, 59e75 kg) who received an infusion of propofol 2.5 mg kg À1 i.v. for 20 s. We also simulated the study by Maneglia and Cousin, 34 which focused on eight females with a mean age of 86 (standard error of the mean [SEM], 5.4) yr and mean weight of 51.8 kg, who received an induction dose of propofol 1 mg kg À1 and a maintenance infu-

Statistical analysis
Model parameters are reported as typical values with associated relative standard errors (RSEs) derived from covariance matrix or sampling importance resampling. 35

Results
All volunteers received propofol according to a step-up and stepdown target-controlled infusion. Propofol plasma concentrations, MAP, HR, and PP observations used for modelling are shown in Online Supplementary 4, Figs 4e1, which shows that MAP decreases in the step-up phase, then increases in the stepdown phase and decreases again around the time of the bolus. The change in HR was similar, yet inverse to the change in MAP. PP appeared to decrease only mildly during propofol infusion.

Pharmacodynamic empirical models
The E max functions characterise the relationship between the predicted propofol plasma concentrations and MAP, HR, and PP. Estimated E max values for MAP and HR were e40.0% and 26.0%, respectively. For PP, the estimated maximum effect of PP (E max_PP ) was e8.9% with individual estimates ranging from e116.8% to 54.9%. We found that E max_PP negatively correlated with baseline PP (Base PP , r¼e0.72), meaning that individuals with higher than typical baseline PP were more likely to have a lower than typical E max_PP (r is the correlation coefficient). The estimated half-lives for the MAP, HR, and PP models were 3.15, 7.79, and 5.78 min, respectively, indicating that MAP changed first followed by PP and HR. The time-dependent effect in equations (2) and (3) was significant for MAP (D objective function value¼e91.0) and HR (D objective function val-ue¼e16.9), and resulted in a temporary increase in MAP by 6 mm Hg and HR by 3 beats min À1 with a half-life for attenuation of 14.74 and 5.46 min, respectively. Age was a significant covariate on Base MAP (D objective function value¼e13.3) and E max_PP (D objective function value¼e15.2), with Base MAP and E max_PP increasing with age according to equations (4) and (5) Parameter estimates for the final empirical models are shown in Online Supplementary 2, Tables 2e2. The goodnessof-fit plots for the separate MAP, HR, and PP models are shown in Online Supplementary 5, Figs 5e1,5-2, and 5-3, respectively. The prediction-corrected visual predictive check plots of the Table 1 Parameter estimates of final mechanism-based model. *Calculated as: ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e u 2 À 1 p Â 100%. k out , first-order dissipation rate constant of SV, HR, and TPR; Base SV , Base HR , and Base TPR , baseline SV, HR, and TPR, respectively; EC 50 TPR and EC 50 SV , concentrations that produce half of the maximal propofol effect on TPR and SV, respectively; E max_SV and E max_TPR , maximal effect of propofol on SV and TPR, respectively; FB, magnitude of the feedback; HR_SV, magnitude of the effect on HR on SV; k, first-order dissipation rate constant of the time-dependent effect; LTDE, percentage of increased baseline HR caused by the time-dependent effect; g, Hill coefficient for TPR sigmoid; AGE_E max_SV , covariate effect of age on the maximal propofol effect on SV; BSV, between-subject variability; r, correlation coefficient; sRUV, residual unexplained variability; RSE, relative standard error.

Mechanism-based model
Propofol plasma concentrations negatively affected the production rates of TPR (k in_TPR ) and SV (k in_SV ), which was best described using a sigmoid E max model and an E max model, respectively. Inclusion of an additional propofol effect on the turnover equation for HR did not improve the model. Initially, we estimated a separate k out for each of the turnover equations for TPR, SV, and HR. However, these parameter estimates were on the same order of magnitude (0.078, 0.072, and 0.073 min À1 for TPR, SV, and HR, respectively). Therefore, one k out was estimated instead (D objective function value¼1.8).
A time-dependent effect, as in equations (2) and (3), was significant for HR and accounted for the elevated HR before the start of the study session (D objective function value¼e537.2). The BSV for the estimated baseline parameters for SV, HR, and TPR (BSV_Base SV , BSV_Base HR , and BSV_Base TPR , respectively) were highly correlated. Inclusion of these correlations (denoted by rBase HR , rBase TPR , and rBase SV in Table 1) in the stochastic model resulted in improvement in the model's goodness-of-fit (D objective function value¼e59.5) and prediction-corrected visual predictive checks. Residual errors on SV, HR, and TPR were described using additive residual error models.
The final mechanism-based model is shown in equations (6)e (11), and the structure of the model is graphically depicted in Fig 1. In these equations, C p represents the plasma propofol concentration; k in_TPR , k in_SV , and k in_HR represent the zeroorder production rate constants and k out is the first-order dissipation rate constant of TPR, SV, and HR. E max_TPR and E max_SV are the maximum effects of propofol on TPR and SV, respectively. EC 50TPR and EC 50SV are the concentrations that produce 50% of the maximal drug effect on TPR and SV, respectively. SV* represents the SV affected by the negative change in MAP and drug effect. HR* represents the HR influenced by the feedback of MAP. RMAP is MAP normalised to baseline values. FB is power of the negative feedback of MAP on TPR, SV, and HR, respectively. HR_SV is a constant that represents the magnitude of the direct inverse effect of HR on SV. Structure of the mechanism-based model that describes the propofol effect and the interdependency between MAP, HR, SV, and TPR. MAP equals the product of TPR, HR, and SV. The influences of negative feedback through MAP on TPR, HR, and SV are represented by three linked turnover equations. Propofol decreases the zero-order production rates (k in ) of the turnover equations for TPR and SV. SV is influenced by the drug effect and negative feedback through MAP (SV*) and by direct inverse relationship through HR (represented by HR_SV). HR is described by combination effects of time-dependent effect (which will cause a temporary increase in HR, expressed as TDE) and feedback from MAP (HR*). PP is a surrogate of SV (taking into account an SV/PP ratio equals to 1.5). In this model, EFF_TPR and EFF_SV represent the drug effect on TPR and SV. k in_HR , k in_SV , and k in_TPR are the zero-order production rate constants of HR, SV, and TPR, respectively. k out is the first-order dissipation rate constant. FB is the power of the feedback on HR, SV and TPR. SV, stroke volume; TPR, total peripheral resistance; PP, pulse pressure. RMAP ¼ HR$SV$TPR Base HR $Base SV $Base TPR (11) Inclusion of age in the E max_SV resulted in a significant decrease in the model's objective function value (e1381.8). E max_SV increased with age according to equation (12). E max_SV(typ) is the population typical value of the maximal effect of propofol on SV. No other covariates were included in the final model. E max sv ¼E max svðtypÞ , e ð0:026,ðAGEÀ35ÞÞ 12 (12) Final mechanism-based model The estimated parameters and associated uncertainty in the estimates are shown in Table 1 Figure 2 shows the steady-state concentrationeeffect relationship for MAP, HR, and SV for a 35-yr-old individual according to the final mechanism-based model. MAP decreases from 86 mm Hg in the absence of propofol to 54 mm Hg at 10 mg ml À1 (relative decrease¼e37.2%). HR changes from 56 to 87 beats min À1 as plasma concentration increases from 0 to 10 mg ml À1 (relative increase¼35.6%). SV decreases slightly (from 82.9 to 75.6 ml) at low plasma propofol concentration (0e1.1 mg ml À1 ) and then increases to 91.2 ml as plasma concentration increases to 10 mg ml À1 .  36 in 44.7% of adults and 56.3% of older adults. Simulation showed that individuals with higher baseline MAP tended to have more chance of hypotension (correlation coefficients calculated from simulation between initial MAP and getting hypotension were e0.59 and e0.57 in adults and older individuals, respectively). Maximum reductions in MAP from baseline were 26.3% and 30.4% in adults and older patients, respectively. Compared with MAP, the change in HR after propofol infusion was moderate with an expected increase of 12.6% and 22.6% from baseline in adults and older individuals, respectively. At the same time, SV decreased slightly from baseline for adults (4.3%), whereas the decrease of SV was more pronounced for older individuals (31.7%).

Effects of propofol infusion in young and older adults
To benchmark the predictive performance of the mechanismbased model, we performed simulations according to published studies and then compared our model predictions with the observations reported in the studies. The comparison between the model prediction and the Fairfield and colleagues study 33 is shown in Table 2. Fairfield and colleagues found that 5 min after propofol induction with 2.5 mg kg À1 , MAP declined from 88 to 67 mm Hg (a relative decrease of e23.9%). In line with these observations, the 95% prediction interval for the mean change from baseline according to our model was e5.5% to e20.9%. In contrast, Fairfield and colleagues found that HR increased after 1 min and then decreased, whereas our simulations did not show a change in HR. This apparent stable HR in the first 12 min of our simulations is attributable to the counteracting time-dependent effect (decreased HR) and the drug effect of propofol (increased HR owing to the influence of feedback from MAP).
Our simulations differ from the observations described by Maneglia and Cousin. 34 Maneglia and Cousin found that in older people (mean age, 85.5 yr), after induction of 1 mg kg À1 and maintenance with 0.1 mg kg min À1 propofol, MAP declined from 96 to 81 mm Hg (e15.6%) at 15 min (Table 3). Owing to the strong age relationship on the maximum effect of SV (E max_SV ) (e.g. E max_SV was e86.2% for an 85-yr-old individual), our model predicts a mean change from baseline MAP of 38.5% (95% prediction interval from e31.9% to e45.1%) and a 20% increase in HR.

Discussion
We developed a mechanism-based pharmacodynamic model that describes propofol-induced changes in MAP, TPR  Table 1 for a 35-yr-old individual. SV, stroke volume; C p , plasma propofol concentration.
(implemented as a latent variable), SV (approximated by PP), and HR, and the interrelationship between these variables in healthy volunteers. According to our model, propofol decreases TPR and SV with a relatively greater sensitivity of TPR (the estimated maximum effects of propofol on TPR and SV were e85.1% and e23.5%, respectively). The reduced TPR and SV results in a decrease in MAP, which further influences SV, HR, and TPR through three negative feedback loops. Age was a significant covariate on E max_SV in the mechanism-based model, where the propofol-induced decrease in SV was greater in older volunteers. Both our mechanism-based model and our empirical models adequately describe the changes in MAP, HR, and PP after propofol infusion. We also showed the potential utility of the model by showing good agreement between simulations and the data on propofol-induced changes in MAP and HR from Fairfield and colleagues. 33 Propofol is known to reduce sympathetic tone and influence smooth muscle Ca 2þ flux, leading to vasodilation, reduced TPR, and ultimately lower blood pressure. 6,37,38 Propofol also lowers left ventricular preload, which reduces SV. 38,39 In line with these findings, we found that propofol reduces the zero-order production rates in the turnover equations for TPR and SV. The steady-state concentrationeeffect relationship of SV shows a biphasic effect resulting from a propofol-driven reduction in SV at low concentrations (<1.1 mg ml À1 ) and an increase in SV at higher concentrations because of the feedback effect from MAP.
Several studies have demonstrated that propofol reduces the tachycardia response to hypotension because of inhibition of the baroreceptor reflex. 37,38,40 Despite these observations, we found that HR significantly increased in healthy volunteers receiving propofol. A possible explanation is that the baroreflex-induced increase in HR outweighs the inhibitory effect of propofol on the baroreceptor reflex, and our model describes the net result of these conflicting mechanisms. Our simulations of changes in MAP and HR differed from the observations of Maneglia and Cousin. 34 The discrepancy between the model predictions and observations in older individuals is because our model suggests a strong relationship between age and E max_SV (e.g. E max_SV is e86.2% for an 85yr-old individual, whereas E max_SV is e23.5% for a 35-yr-old individual). According to equation (1), 22 the decrease in MAP from baseline is expected to be greater for an individual with a stronger effect on SV, and the increase in HR will be more pronounced owing to the negative feedback mechanism. 13 Nevertheless, the age relationship on E max_SV needs to be confirmed in our future studies. We identified a time-dependent effect on haemodynamic variables and hypothesise that this reflects a time-dependent adaptation of the feedback mechanisms in the cardiovascular system in response to propofol exposure. Unfortunately, our data did not allow us to implement this time-dependency in a more mechanistic manner. More data, especially more diverse data in terms of the propofol infusion regimen, are likely needed to improve our understanding of the origin of this time dependency.
Wada and colleagues 41 and Upton and colleagues 42 have shown for thiopental and propofol that drug-induced alterations in liver blood flow can impact the pharmacokinetics. When this is the case, a sequential PKePD modelling approach, as used in this analysis, is not appropriate and a simultaneous PKePD modelling approach is more useful.
We used the Eleveld propofol PK model 1 to derive the post hoc PK parameters for all subjects before PD modelling. This PK model assumes linear kinetics and does not take into account potential cardiac output or liver blood flow dependent pharmacokinetics. Nevertheless, goodness-of-fit plots (Online Supplementary Figs 3e1) showed that the Eleveld PK model accurately described the current dataset, suggesting that a propofol-induced decrease in liver blood flow is not present in our subjects and is not confounding our results. This is the first quantitative analysis of the cardiovascular adverse effects of propofol in a mechanism-based PKePD model in humans. Several studies have characterised the haemodynamic effects of propofol using more empirical models. 8e10 Jeleazcov and colleagues 8 developed a PKePD model to characterise the effects of propofol on MAP in which the concentrationeeffect relationships was described by a direct response model with two effect-site compartments. These two effect sites indicate that the effect of propofol on MAP is mediated by two pathways (e.g. changes in CO and systemic vascular resistance). When we use the model by Jeleazcov and colleagues (ignoring the differences in pharmacokinetic models) to predict the maximum decrease in MAP according to our study design, results are very similar (e38.8% for a 35-yr-old individual in our model vs e47.9% for the model by Jeleazcov and colleagues 8 ).
Wu and colleagues 10 developed a PKePD model describing the changes in MAP after propofol administration in normal weight and morbidly obese patients (BMI >35 kg m À2 ). Their model contains an effect-compartment model linked to the central compartment of a two-compartment pharmacokinetic model. We cannot directly compare parameters between the model by Wu and colleagues and our final model as there is no drug effect parameter on MAP in our mechanism-based model. However, when we simulate the maximum effect on MAP after propofol infusion according to our study (for a normal weight patient) using only the pharmacodynamic model presented by Wu and colleagues, results are in good agreement. In the model by Wu and colleagues, MAP decreased from 96 to 62 mm Hg (relative decrease of e35.1%), whereas our model predicts a decline in MAP from 90 to 55 mm Hg (relative decrease of e38.8%).
Our results are in general agreement with the work of Kazama and colleagues, 43 who developed an effect-site compartmental model for propofol on systolic blood pressure (SBP). The half-lives for dissipation of propofol-induced drug effects on SBP in patients aged 20e39, 40e59, 60e69, or 70e85 yr were 4.68, 5.92, 8.87, and 10.22 min, respectively, which is similar to the estimated value from our mechanistic model (9.19 min). In addition, Kazama and colleagues found that the effect on SBP increased with increasing age. The estimated EC 50 in these four groups was 4.61, 4.13, 3.96, and 2.09 mg ml À1 , respectively. In our model, a similar effect was found with E max_SV increasing with increasing age, resulting in a decrease in MAP during propofol infusion with increasing age.
The model by Snelder and colleagues 12,13 was used as a starting point in our mechanism-based model development.
There are several differences between the models. First, the Snelder model was developed based on data from rats whereas our study used data from humans. Although there are differences between rats and humans in cardiac electrophysiology, 44 the structure we used for our model is  physiologically reasonable. 22 Second, the Snelder model was developed using experimental data where SV was measured directly, whereas for our study SV measurements were unavailable and we used PP data as a surrogate of SV. 23e25 Third, feedback mechanisms were implemented as linear relationships between MAP and the rate constant product of HR, SV, and TPR in Snelder's model, whereas a power function describes our data better. Finally, circadian rhythm was not included in our model because of the short time course of the study (~5 h). There are also similarities between the two models. The structure of our mechanism-based model is the same as that of the Snelder model and only differs in implementation of the drug effect. Moreover, both models account for a potential time-varying effect that causes a temporary increase in HR and MAP. 12,13 The Snelder model was further extended by Venkatasubramanian and colleagues 14 based on data in dogs to include cardiac contractility (dPdt max ). In this refined model, a more detailed mechanistic description of SV was used that included cardiac contractility and HR. Despite this refinement, translation from dPdt max to any clinical measure of contractility is still under debate and warrants further investigation. 14 Chae and colleagues 15 used a similar mechanistic model to describe the haemodynamic changes induced by telmisartan in healthy male volunteers. Their model also includes the interrelationships between haemodynamic variables and baroreflex feedback control. There are still several differences between our mechanism-based model and the model developed by Chae and colleagues. Baroreflex control is described in the Chae model by a turnover function taking into account the fractional change in MAP, whereas we estimated the parameter FB representing the power of feedback similar to the Snelder model. 13 In our dataset, the approach by Chae and colleagues 15 performed worse than using the parameter FB as in our final model. In addition, the inverse relationship between HR and SV is not included in the Chae model.
Bahnasawy and colleagues 16 modified the Chae model by including the inverse relationship between HR and SV. To describe the short-term changes in haemodynamics, k out was fixed to 35 h À1 (half-life of 1.20 min), which is in the same order magnitude as estimates from our model (9.19 min). Bahnasawy and colleagues also showed that baseline SV needs to be fixed if there are only MAP and HR measurements to keep the model identifiable, whereas in our model PP was included, which served as a surrogate of SV and therefore baseline SV could be estimated.
Application of a mechanism-based PKePD model allows a more granular description to complex physiological systems, leading to better understanding of the system through separation of drug-specific parameters (e.g. E max and EC 50 ) and biological system-specific parameters (e.g. k out and FB). 11,45 This property enables better understanding of changes amongst the haemodynamic effects by simultaneously considering both drug effect and correlations between cardiovascular effects.
There are some limitations to our study. First, we used data from healthy volunteers. Although noxious stimuli were applied, we acknowledge that the level of noxious stimulation was less than that during surgery. However, use of healthy volunteer data is a strength in that it allows us to calibrate the model to data that lack confounding factors that are known to be present in patient data (medications, cardiovascular disease, etc.). It is reassuring that predictions from our final model are in reasonable agreement with results from studies in patients. 33 In a follow-up project, this model will serve as a basis to describe the haemodynamic effects of propofoleremifentanil combinations both in healthy volunteers and in patients, which will explore additional factors that can help explain the variability that we observe in the druginduced haemodynamic responses.
Second, we did not have SV measurements and so used PP as a surrogate of SV based on the assumption that total arterial compliance remained constant. 46 Nevertheless, propofol can cause an increase in compliance. 47 At the same time, the ratio of SV to PP has been shown to be related to body size, age, and heart rate. 23 A sensitivity analysis showed that our data are not informative for estimating SV/PP ratio or for applying a more mechanistic equation such as the equation presented by de Simone and colleagues. 23 As a consequence, predicted changes in SV or CO from our model remain uncertain and should be evaluated in future research.
In conclusion, we present a mechanism-based pharmacodynamic model for propofol effects on MAP, HR, SV, and TPR in healthy volunteers. The mechanism-based model performance is equally good to the performance of three empirical models for MAP, HR, and PP. The similarities between model simulations and published observations indicate that the model has reasonable predictive properties.