FormalPara Key Points

Using previously collected data and a sequential modelling approach, we were able to describe dexmedetomidine absorption kinetics and further disposition together with its sympatholytic and haemodynamic effects.

Accurate models for haemodynamic effects were obtained by combining norepinephrine-dependent sympatholysis with either dexmedetomidine-evoked hypertension or changes in central neural activity.

Our final models precisely describe dexmedetomidine pharmacokinetics and accurately predicts dexmedetomidine-induced sympatholysis and other pharmacodynamic effects.

1 Introduction

Patients in palliative care frequently present with pain and anxiety. These symptoms are commonly managed with opioids, benzodiazepines, antidepressants and antipsychotic agents [1, 2], all of which have clinically significant adverse effects. Alternatives to the current regimes are therefore needed to achieve clinically adequate analgesia and anxiolysis, while avoiding undesired drug effects.

Dexmedetomidine is a selective and potent sedative used in anaesthesia and pain medicine for its anxiolytic, sedative and analgesic properties [3]. Dexmedetomidine mediates its effects via α2-adrenoceptors to inhibit the release of norepinephrine (NE) from presynaptic nerve terminals and to hyperpolarise postsynaptic neurons [4]. Its unique mechanism of action causes dose-dependent sedation without a significant risk of respiratory depression, making dexmedetomidine a promising agent for palliative sedation.

Although dexmedetomidine is only registered for intravenous (IV) administration, many alternative routes of administration have been documented for off-label use [3,4,5]. Our previous work indicates that subcutaneous (SC) infusion of dexmedetomidine is a feasible alternative to the IV route, especially when venous access is unattainable [6]. In our article, we have documented that the plasma concentrations of dexmedetomidine increase to clinically adequate sedative concentrations after SC administration, while the sympatholytic effects are significantly attenuated as compared with IV dosing [6]. The sedative profile of dexmedetomidine has been very well characterised [7, 8], but only a few studies have evaluated the pharmacokinetic–pharmacodynamic (PK–PD) relationship and quantitatively described the various clinically observed effects of dexmedetomidine [9,10,11,12,13], especially after extravascular dosing [9].

Population-based PK analysis of dexmedetomidine has been reported in the literature, with a range of study subjects from healthy adult volunteers [9, 13] to hospitalised children [14,15,16], and critically ill adult patients [17, 18]. Most of these studies were conducted using data collected after IV dosing, and PD data were included only in a few studies [9,10,11,12,13]. To our knowledge, population-based PK-PD analysis of subcutaneously administered dexmedetomidine has not been previously described.

Our aim was to develop semi-mechanistic population models using healthy volunteers’ data to characterise the pharmacokinetics of dexmedetomidine and its sympatholytic and haemodynamic effects after IV and SC dosing. A suitable dataset was available from our previously published clinical study [6]. We hypothesised that our population-based PK-PD models would be able to: (1) accurately describe absorption of subcutaneously administered dexmedetomidine, (2) account for the disposition kinetics of dexmedetomidine, and (3) predict the effects of dexmedetomidine on plasma NE levels, heart rate (HR), and systolic and diastolic arterial blood pressure (SAP and DAP, respectively) after IV and SC dosing.

2 Materials and Methods

2.1 Study Protocol

We have reanalysed our data from a previously published clinical trial [6] (registered as ClinicalTrials.org identifier NCT02724098). The study protocol (EudraCT 2015-004698-34) was approved by the Ethics Committee of the Hospital District of Southwest Finland (184/1800/2015) and by the Finnish National Agency for Medicines (3/2016). Detailed information about the study protocol has been provided in the Electronic Supplementary Material (ESM). Briefly, a two-period randomised crossover study was conducted to investigate the PK and pharmacological effects of subcutaneously and intravenously administered dexmedetomidine in ten (two dropouts for the IV arm of the study) healthy human volunteers. We administered either a single IV or SC dose of 1 µg/kg of dexmedetomidine during a 10-min infusion in a randomised order. Arterial blood samples were collected pre-dose and then at 5, 10, 15, 20, 30 and 45 min and 1, 1.5, 2 and 3 h after the start of the infusions into EDTA tubes for determination of concentrations of dexmedetomidine, NE and epinephrine in plasma. Additional blood samples were collected at 4, 5, 6, 8 and 10 h for determination of dexmedetomidine. Heart rate, SAP, DAP and peripheral oxygen saturation were recorded at times of blood sampling. Plasma concentrations of dexmedetomidine were determined with high-performance liquid chromatography-tandem mass spectrometry and concentrations of NE and epinephrine in plasma were determined with high-performance liquid chromatography and coulometric electrochemical detection as described previously [6, 19].

2.2 Data Analysis

Non-linear mixed effects modelling was conducted using the NONMEM® (version 7.3.0) software package (ICON Development Solutions, Ellicott City, MD, USA) [20] and supported by the Perl-Speaks-NONMEM toolkit (PsN version 4.6.0) [21], for model development, execution and evaluation. RStudio (version 1.0.153) using R (version 3.6.0) [22] was used for diagnostic model evaluation and graphical visualisations.

2.3 Model Development

Population models were specified with differential equations that were solved with ADVAN13 subroutine and parameters estimated using first-order conditional estimation with interaction. A sequential PK-PD approach was used as described previously [23, 24]. In the first step, a PK model was developed using IV administration data only, followed by the implementation of a semi-mechanistic absorption model for SC administration. Next, a PD model was developed for NE/epinephrine concentration–time data. Finally, effect models for the PD parameters, i.e. HR, SAP and DAP were developed.

2.3.1 Structural Models

For dexmedetomidine pharmacokinetics, we started with a simple one-compartment empirical model, and nested models were built in a stepwise manner using two- and three-compartmental models. Zero-order input of IV dexmedetomidine, and linear pharmacokinetics with first-order elimination clearance were assumed. A semi-mechanistic model describing dexmedetomidine absorption from the SC site of administration during a SC infusion were added to the disposition model. The inhibitory effect of dexmedetomidine on the spill-over of NE and epinephrine into plasma was modelled using indirect response models. Decreases in catecholamine concentrations and their effects were then used as surrogate measures of the changes in subsequent pharmacodynamic endpoints.

2.3.2 Stochastic Model

The between-subject variability on the fixed-effect parameters in the PK and PD models was specified using exponential model,

$${\Phi}_{i}=\theta \bullet {\text{e}}^{{\eta }_{i}},$$

where \({\Phi}_{i}\) is the fixed effect parameter for an individual, \(\theta\) is the population estimate of the parameter and \({\eta }_{i}\) describes the between-subject variability of the ith individual from the population estimate.

Between-subject variability was added to the fixed-effect parameters in the models gradually in nested models and was retained if it led to a significant reduction in the objective function value (OFV), did not cause numerical instability in the models, and if the shrinkage of the added between-subject variability was less than or equal to 25%.

For the PK models of intravenously and subcutaneously administered dexmedetomidine, the residual variability between observations was coded into the models with a constant coefficient of variation error model,

$${C}_{ij,x}={\widehat{C}}_{ij,x}\bullet \left(1+{\varepsilon }_{ij,x}\right).$$

An additive error model was used for NE and epinephrine concentrations, and for the remaining PD variables,

$${C}_{ij,x}={\widehat{C}}_{ij,x}+{\varepsilon }_{ij,x,}$$

where x is one of the dexmedetomidine/NE/epinephrine concentrations or one of the remaining PD endpoints, \({C}_{ij,x}\) is the jth observation for the ith individual at the xth time point, \({\widehat{C}}_{ij,x}\) is the corresponding model-predicted value and \({\varepsilon }_{ij,x}\) is a normally distributed random variable to quantify the residual variability in the values of x. \(\eta\) and \(\varepsilon\) were assumed to be normally distributed random-effect parameters.

2.3.3 Model Evaluation

The primary criterion for discerning amongst nested models was the difference in the OFV (ΔOFV) at a significance level of 0.05, which equates to a reduction of 3.84 points in the OFV for every added degree of freedom. Akaike information criteria (AIC) was used for non-nested models. Relative standard errors and retrieval of covariance matrices were used to check model precision and parameter identifiability. Standard goodness-of-fit plots and visual predictive checks [25] were routinely used to test model appropriateness and predictive performance. The precision in the final model parameters was estimated using the sampling importance resampling (SIR) procedure [26]. To achieve robustness in the final SIR results, the procedure was run with 40,000 samples and 2000 resamples. Covariance matrices from the final models were used as a sample proposal distribution. Output from these SIR runs were evaluated using the 95% confidence intervals in the parameter estimates, OFV distribution and graphical output from the SIR package in PsN.

2.3.4 Simulations

We also assessed model performance with simulations using the final parameter estimates for the fixed and random effects in 1000 virtual subjects with NONMEM. The first set of simulations was conducted with the final PK model only, to evaluate model output for dexmedetomidine concentration–time profiles for short- and long-term infusions using the IV and SC routes. Two clinically relevant infusion time scales of 8 h and 48 h were used at dexmedetomidine doses of 0.25, 0.5, 1, 2, 3 and 4 µg/kg/h. These simulations were conducted to evaluate the time required for achieving steady-state plasma concentrations and to compare relative areas under the curve achieved with SC administration as compared to the IV route. Additionally, we evaluated the deposition of dexmedetomidine in the SC fat tissue after SC dosing to investigate possible accumulation after SC drug administration for possible clinical implications.

Next, the same dosing schemes were used to simulate the effect of dexmedetomidine on the NE concentration–time course in plasma. These simulations were used to evaluate decreases in NE levels during continuous infusions in a dose-dependent manner. The aim was to check a model’s performance in dynamically predicting the inhibition of NE spill-over over increasing doses and in a temporal manner, and to identify the dose required for a near-complete cessation of NE spill-over into the plasma.

Third, different dexmedetomidine dosing schemes and the resulting NE profiles were used to evaluate drug effects on haemodynamic parameters. Because dexmedetomidine influences blood pressure by two distinct mechanisms, model simulations were used to provide information on the development of NE-associated hypotension (sympatholysis) and dexmedetomidine-evoked hypertension (direct vasoconstriction) over an increasing dose scale. Simultaneously, HR model simulations were conducted to evaluate the dose dependency of HR inhibition by dexmedetomidine, as well as the interplay between the inhibition of NE release and other central nervous system (CNS) effects contributing to dexmedetomidine-evoked bradycardia. Last, the ncappc package of PsN [27] was used as an internal evaluation method to perform a non-compartmental analysis (NCA)-based comparison between the observed values in the dataset to model-derived predictions for the PK parameters.

3 Results

A schematic of the final PK–PD model is shown in Fig. 1 and Table 1 shows the comparison between different models tested during model development. A more detailed description of the model development including the set of mass transfer differential equations is provided in the ESM. The model estimates in all the final PK and PD models were computed without bias, and numerical and graphical diagnostics demonstrated adequate precision and robustness of the PK (Table 2, Fig. 2) and PD models (Table 3, Figs. 2, 3, 4).

Fig. 1
figure 1

Schematic of the final semi-mechanistic model for dexmedetomidine pharmacokinetics (a), its effect on norepinephrine (NE) release (b), systolic and diastolic blood pressures (c), and heart rate (d). A amount, C1 dexmedetomidine concentration on central compartment, CE amount of (heart rate/systolic or diastolic blood pressure), CNS central nervous system, EC effect compartment, EC50 concentration causing 50% effect from EMAX, EMAX maximum effect, Fspill fraction on NE spilled out from the release compartment, HPN hypotension, HR heart rate, HTN hypertension, I(t) dosing event, IV intravenous, k1e, k5e and k7e rate constants for distribution to the effect compartment, kIN,R rate constant for NE release, ka absorption rate constant, ke0 elimination rate constant for the effect, kOUT,P rate constant for NE elimination from the plasma, P plasma, R release, SC subcutaneous, \(\gamma\) hill parameter

Table 1 Comparison of different models tested during model development
Table 2 Parameter estimates from the final pharmacokinetic model for dexmedetomidine (DEX) and norepinephrine (NE) with median and 95% confidence intervals (CIs) from the sampling importance resampling (SIR) procedure with 40,000 final samples and 2000 resamples
Fig. 2
figure 2

Visual predictive checks based on 1000 simulations showing dexmedetomidine and norepinephrine after a 1-µg/kg intravenous (a, c) or subcutaneous (b, d) dexmedetomidine infusion. The black solid and dashed lines are the observed percentiles (10th, 50th and 90th percentiles) respectively, and the blue ribbon is the corresponding median predictive interval. Light blue ribbons are predicted percentiles. Black circles are individual observations

Table 3 Pharmacodynamic parameters for blood pressure and heart rate (HR) models
Fig. 3
figure 3

Visual predictive checks based on 1000 simulations showing heart rate (upper row), systolic blood pressure (middle row) and diastolic blood pressure (lower row) after a 1-µg/kg intravenous (a, c, e) or subcutaneous (b, d, f) dexmedetomidine infusion. The black solid and dashed lines are the observed percentiles (10th, 50th, and 90th percentiles), respectively, and the blue ribbon is the corresponding median predictive interval. Light blue ribbons are predicted percentiles. Black circles are individual observations

Fig. 4
figure 4

Simulated dexmedetomidine concentration profiles in plasma using the final mechanism-based model and 8-h (upper row) and 48-h (lower row) intravenous (a, c) and subcutaneous (b, d) continuous infusions of varying dexmedetomidine doses

3.1 Pharmacokinetic Model

Pharmacokinetic modelling of dexmedetomidine concentration–time data began with an empirical one compartment model (OFV = \(-306\)), but the model fit was poor and therefore two- and three-compartment mammillary models were considered. Although the OFV reduced significantly for a three-compartment model (ΔOFV = \(-348\)) as compared with a two-compartment model (ΔOFV = \(-325\)), bootstrap analysis demonstrated wide confidence intervals for the added model parameters, and therefore the two compartmental system was retained as the final structural PK model.

A first-order absorption process was initially tested to describe the absorption of dexmedetomidine after a short SC infusion (ΔOFV = \(-659\), AIC = \(-1276\)), as previously reported for nasal mucosal bioavailability [9], but visual predictive checks demonstrated inconsistent output. Multiple different solutions were tested (see the ESM, p. 3) without success. Finally, a semi-mechanistic biphasic absorption model was constructed (AIC = \(-1436\)), such that a fast absorption process \(({k}_{a,\text{F}\text{A}\text{S}\text{T}})\) governed the influx of drug from a SC depot into the systemic circulation and a fat layer permeation constant \(({k}_{\text{F}\text{A}\text{T}})\) was added to describe drug distribution into the local SC fat tissue. The addition of bioavailability parameters in the final PK model for the SC dosing phase improved the model fit (\(\Delta \text{O}\text{F}\text{V}= -201, \text{A}\text{I}\text{C}= -1469\)).

Fractional permeation of subcutaneously administered dexmedetomidine was estimated in the PK model (\({F}_{\text{F}\text{A}\text{T}}=0.55\)), while the fast absorption process from the depot (\({K}_{a, \text{F}\text{A}\text{S}\text{T}}\)) was linked with a bioavailability parameter (\({F}_{\text{D}\text{E}\text{P}\text{O}\text{T}}=0.79\)) that governed the systemic procurement of dexmedetomidine only via the fast process. A slow absorption process \(({k}_{a,\text{S}\text{L}\text{O}\text{W}})\) was included to describe drug release from the SC fat tissue. Although this scenario represented a biphasic absorption process, no strict criterion was coded into the model to describe the temporal shift between the two processes, and both were modelled simultaneously to represent a dynamic interplay between fast and slow absorption processes from the injection site. An attempt to model a bioavailability parameter associated with slow absorption from a fat compartment failed with numerical instability and over-parametrisation. Net bioavailability of dexmedetomidine (0.89) through SC administration was estimated assuming complete recovery of the drug from the fat layer into the blood circulation.

The final PK model produced biologically plausible parameter estimates (Table 2). Our results demonstrate that inclusion of a fat compartment into the model allowed adequate description of the data. The rate constants governing SC absorption were estimated with high precision (Table 2) and the visual predictive checks for the final PK model (Fig. 2) show good predictive performance. In addition, the results indicate that the ratio of fast to slow absorption of dexmedetomidine from the local injection depot and from the SC fat tissue compartment is approximately 1:9, while the capacity of permeation into the SC fat tissue compartment is about three-fold compared to the capacity of the simultaneous fast absorption process from the depot compartment into the systemic circulation (Table 2).

The final PK model was used to calculate the area under the curve \(({\text{A}\text{U}\text{C}}_{0-\infty }=\left(\text{D}\text{O}\text{S}\text{E}\times F\right)/\text{C}\text{L})\) levels for the study population, and the levels of \({\text{A}\text{U}\text{C}}_{0-\infty }\) were also simulated using the NCA package in PsN [21]. The model calculated levels of the \({\text{A}\text{U}\text{C}}_{0-\infty }\) were \(2.01\) and \(1.82\) \(\text{n}\text{g} \cdot \text{h}/\text{m}\text{L}\) for the IV and SC phase, respectively, while the simulated values of these variables were synchronous with calculated results, at \(1.94\) and \(1.78 \; \text{n}\text{g} \text{h}/\text{m}\text{L}\) respectively using the NCA package. Although our results demonstrate a lower relative bioavailability of dexmedetomidine in the SC dosing phase, comparable values of \({\text{A}\text{U}\text{C}}_{0-\infty }\) establish the strong clinical utility of the SC route of dexmedetomidine administration.

3.2 Norepinephrine Model

An indirect response model of NE was used, and baseline levels of NE were estimated for each individual in the final model with an inter-individual variability parameter. Plasma levels of NE represent an approximately 15% spill-over fraction of the amount of NE released from sympathetic nerve endings [28], and our NE model was specified as a composite of two separate compartments, representing sympathetic nerves (release compartment) and plasma (circulatory compartment), respectively. Thereafter, the initial conditions for the NE release compartment were specified with the rate constants for the resting-state release and spill-over of NE from the release compartment. Because dexmedetomidine is known to reduce the rate of the resting-state release of NE from sympathetic nerve endings in vivo, an inhibitory effect on the release rate constant (\({k}_{\text{I}\text{N},\text{R}}\)) at \(\text{t}\text{i}\text{m}\text{e}=t\) was coded into the model using a indirect response model,

$${k}_{\text{I}\text{N},\text{R}}={{k}_{\text{I}\text{N},\text{R}}}^{t=0}\bullet \left(1-\frac{{C}_{1}}{{\text{E}\text{C}}_{50,\text{D}\text{E}\text{X}}+{C}_{1}}\right),$$

where \({{k}_{\text{I}\text{N},\text{R}}}^{t=0}\) defines the resting-state rate of NE release at \(t=0\), \({C}_{1}\) is the concentration of dexmedetomidine in the central compartment and \(\text{E}{\text{C}}_{50, \text{D}\text{E}\text{X}}\) is the concentration of dexmedetomidine that causes 50% of the maximal inhibitory effect on NE release in vivo.

Dexmedetomidine \({\text{E}\text{C}}_{50}\) for NE release was estimated at approximately \(0.33 \text{n}\text{g}/\text{m}\text{L}\) (Table 2), which is numerically close to the maximum concentrations achieved in the study subjects with SC dosing of 1 µg/kg, and much lower than the drug concentrations achieved after IV dosing. The values of the rate constant for NE spill-over from the release compartment \(({k}_{\text{O}\text{U}\text{T},\text{R}})\) and the rate constant for NE elimination from the plasma \(({k}_{\text{O}\text{U}\text{T},\text{P}})\) were \(9.03 \; {\text{h}}^{-1}\) and \(11.4\; {\text{h}}^{-1}\), respectively, which are in accordance with a recent report modelling the inhibition of NE release following dexmedetomidine challenge in humans [9].

3.3 Epinephrine Model

Our data for epinephrine concentrations included multiple below limit of quantification values. All attempts to model epinephrine failed, indicating that these models were over-parametrised to provide a numerically stable predictive model output.

3.4 Systolic/Diastolic Arterial Blood Pressure Models

Decreases in SAP and DAP after dexmedetomidine administration were assumed to result from reduced NE levels in the release compartment, and a sigmoidal \({E}_{\text{M}\text{A}\text{X}}\) model with a biophase was utilised to estimate individual SAP and DAP values. The hypertensive effect of dexmedetomidine through α2-adrenoceptors has been documented before [13], and therefore we assumed that a binomial response, i.e. combined, and dose-dependent hypertensive (α2-adrenoceptor mediated) and hypotensive (NE mediated) effect was applicable to both SAP and DAP. The hypotensive effects of NE \(({E}_{\text{B}\text{P}, \text{H}\text{P}\text{N}})\) concentration fluctuations and hypertensive blood pressure effects due to plasma dexmedetomidine concentrations \(({E}_{\text{B}\text{P}, \text{H}\text{T}\text{N}})\) on SAP and DAP were defined in the model as,

$${E}_{\text{B}\text{P},\text{H}\text{P}\text{N}}=1+\left(\frac{{E}_{\text{M}\text{A}\text{X}}\bullet {{C}_{\text{E},\text{N}\text{E}}}^{{\gamma }_{\text{H}\text{P}\text{N}}}}{{{\text{E}\text{C}}_{50, \text{N}\text{E}}}^{{\gamma }_{\text{H}\text{P}\text{N}}}+{{C}_{\text{E},\text{N}\text{E}}}^{{\gamma }_{\text{H}\text{P}\text{N}}}}\right),$$
$${E}_{\text{B}\text{P},\text{H}\text{T}\text{N}}=1+\left(\frac{{E}_{\text{M}\text{A}\text{X}}\bullet {{C}_{\text{E},\text{D}\text{E}\text{X}}}^{{\gamma }_{\text{H}\text{T}\text{N}}}}{{{\text{E}\text{C}}_{50,\text{D}\text{E}\text{X}}}^{{\gamma }_{\text{H}\text{T}\text{N}}}+{{C}_{\text{E},\text{D}\text{E}\text{X}}}^{{\gamma }_{\text{H}\text{T}\text{N}}}}\right),$$

where \({E}_{\text{M}\text{A}\text{X}}\) is the maximum effect on SAP or DAP and was fixed to 1 for both of these circumstances, \({C}_{\text{E}}\) is the apparent concentration of either NE or dexmedetomidine in the effect compartments, and \({\text{E}\text{C}}_{50}\) is the concentration of NE or dexmedetomidine resulting in half-maximal effects for hypotension and hypertension, respectively. \(\gamma\) is the hill coefficient that signifies the slope of the concentration–response curve for the development of hypotensive or hypertensive effects.

Then, the combined hypotensive and hypertensive effect of NE and dexmedetomidine on SAP or DAP, respectively, was modelled as,

$${E}_{\text{B}\text{P}}= {E}_{\text{M}\text{I}\text{N},\text{B}\text{P}}\bullet \left({E}_{\text{B}\text{P},\text{H}\text{P}\text{N}}+ {E}_{\text{B}\text{P},\text{H}\text{T}\text{N}}\right),$$

where \({E}_{\text{M}\text{I}\text{N}, \text{B}\text{P}}\) is the minimum blood pressure estimated as a model parameter.

A binomial response model for SAP and DAP data was clearly preferable over the models that used NE as the only surrogate measure of these PD endpoints. Models using NE as the only measure of blood pressure development produced final OFV values of \(-567\) and \(-828,\) respectively. A binomial response model for SAP and DAP data with hypertensive action of dexmedetomidine resulted in a significant OFV reduction (− 622 and − 920 for SAP and DAP, respectively).

The effect-site rate constant for hypotensive NE effect on SBP levels \(({k}_{\text{e}0,\text{H}\text{P}\text{N}})\) was estimated at \(22.1\; {\text{h}}^{-1}\) and for DBP at \(24.3\; {\text{h}}^{-1}\) (Table 3), showing a fast equilibration between the release compartment and the effect site, which is biologically plausible considering the action of NE on blood pressure in humans. The model-predicted \(\text{E}{\text{C}}_{50,\text{N}\text{E}}\) for SAP and DAP was \(4.56 \; \text{n}\text{g}/\text{m}\text{L}\) and \(4.32\; \text{n}\text{g}/\text{m}\text{L}\), respectively, and the values of \({\upgamma }_{\text{H}\text{P}\text{N}}\) were \(6.01\) and \(4.\)99 for SAP and DAP, respectively, showing that small changes in NE levels can have marked changes in blood pressure.

In comparison, the values of effect-site rate constants for a dexmedetomidine-induced hypertensive effect \(\left({k}_{\text{e}0,\text{H}\text{T}\text{N}}\right)\) on SAP and DAP were estimated at \(22.1\; {\text{h}}^{-1}\) and \(30 \; {\text{h}}^{-1}\), respectively, which indicates rapid equilibration between the plasma and effect compartment but also signifies very short half-lives for these processes in vivo. The values of these rate constants were fixed to their respective population estimates to ease up the estimation routines of the final models, considering that dexmedetomidine-induced vasoconstriction is a documented phenomenon [13]. Interestingly, the \(\text{E}{\text{C}}_{50,\text{D}\text{E}\text{X}}\) of dexmedetomidine for SAP and DAP was noticed to be considerably different, at \(2.37 \; \text{n}\text{g}/\text{m}\text{L}\) and \(0.81 \; \text{n}\text{g}/\text{m}\text{L}\), respectively, while \({\upgamma }_{\text{H}\text{T}\text{N}}\) for SAP and DAP was \(1.52\) and \(1.99,\) respectively, indicating that the hypertensive effect of dexmedetomidine on blood pressure is less pronounced compared with the hypotensive effect mediated by NE inhibition.

3.5 Heart Rate Model

Bradycardia after dexmedetomidine dosing was initially modelled with the sympatholytic effect of NE release inhibition only, as previously documented [9]. Inhibition of NE release was able to account for the general trend in the HR data, but a clear model misspecification was observable in the model predictions, especially after SC dexmedetomidine administration. We added an inhibitory effect of SAP on HR, mediated by a baroreceptor reflex in the CNS, to account for the misfit in the HR model, which resolved the model misspecification and resulted in an unbiased output.

The effect of SAP on CNS neural activity was specified using an indirect response model as follows:

$${k}_{\text{i}\text{n},\text{C}\text{N}\text{S}}= {{k}_{\text{i}\text{n}, \text{C}\text{N}\text{S}}}^{t=0}\bullet \left(1-\left(\frac{{\text{I}\text{P}\text{R}\text{E}\text{D}}_{\text{S}\text{A}\text{P}}}{{\text{E}\text{C}}_{50,\text{S}\text{A}\text{P}}+{\text{I}\text{P}\text{R}\text{E}\text{D}}_{\text{S}\text{A}\text{P}}}\right)\right)$$

where \({k}_{\text{i}\text{n},\text{C}\text{N}\text{S}}\) is the physiological rate of development for CNS neural activity, \({{k}_{\text{i}\text{n}, \text{C}\text{N}\text{S}}}^{t=0}\) are the initial conditions for the neural reflex at time = 0, \({\text{I}\text{P}\text{R}\text{E}\text{D}}_{\text{S}\text{A}\text{P}}\) is the individual predicted SAP levels from the final model, and \({\text{E}\text{C}}_{50,\text{S}\text{A}\text{P}}\) is the level of SAP required to cause a 50% change in central neural activity.

We assumed that the overall HR effect caused by dexmedetomidine is the sum of the HR-modulating reflex (\({E}_{\text{H}\text{R},\text{C}\text{N}\text{S}}\)) and sympatholytic (\({E}_{\text{H}\text{R},\text{N}\text{E}}\)) mechanisms as pointed out recently [13]. The former was considered to result from an initial vasoconstrictive hypertension after dexmedetomidine administration, which thereby modulates central neural activity through the baroreceptor reflex, causing a diminished HR effect. The latter driver of the HR effect was coded using NE as a surrogate measure, while reduced NE release and resulting sympatholysis were anticipated to cause bradycardia.

Hence, the two components of the overall effect were estimated simultaneously in the final HR model as follows:

$${E}_{\text{H}\text{R},\text{C}\text{N}\text{S}}=1-\left(\frac{{E}_{\text{M}\text{A}\text{X},\text{H}\text{R}}\bullet {{C}_{\text{E},\text{C}\text{N}\text{S}}}^{{\gamma }_{\text{N}\text{R}}}}{{{\text{E}\text{C}}_{50,\text{C}\text{N}\text{S}}}^{{\gamma }_{\text{N}\text{R}}}+{{C}_{\text{E},\text{C}\text{N}\text{S}}}^{{\gamma }_{\text{N}\text{R}}}}\right),$$
$${E}_{\text{H}\text{R},\text{N}\text{E}}=1+\left(\frac{{E}_{\text{M}\text{A}\text{X}, \text{H}\text{R}}\bullet {{C}_{\text{E},\text{N}\text{E}}}^{{\gamma }_{\text{N}\text{E}}}}{{{\text{E}\text{C}}_{50,\text{N}\text{E}}}^{{\gamma }_{\text{N}\text{E}}}+{{C}_{\text{E},\text{N}\text{E}}}^{{\gamma }_{\text{N}\text{E}}}}\right),$$
$${E}_{\text{H}\text{R}}={E}_{\text{M}\text{I}\text{N},\text{H}\text{R}}\bullet \left({E}_{\text{H}\text{R},\text{C}\text{N}\text{S}}+{E}_{\text{H}\text{R},\text{N}\text{E}}\right).$$

This approach led to a numerically stable model and resulted in physiologically plausible parameter estimates and corrected the previous misspecification in the HR model (Table 3). Norepinephrine \(\text{E}{\text{C}}_{50,\text{N}\text{E}}\) for HR effect was estimated at \(3.12 \; \text{n}\text{M}\), while \({\gamma }_{\text{N}\text{E}}\) was \(2.85\), showing a steep concentration–response relationship for NE on HR development. \({k}_{\text{e}0,{\text{H}\text{R}}_{\text{N}\text{E}}}\) for NE at the HR effect compartment was \(6.15 \;{\text{h}}^{-1}\), suggesting that the effect is slower as compared with the NE on SAP/DAP but has a longer half-life.

Additionally, \({k}_{\text{e}0,\text{N}\text{R}}\) for the effect of neural reflex on HR was \(0.438\; {\text{h}}^{-1}\), which demonstrates a long half-life for this process, which is biologically plausible. The final parameter estimate for \({\gamma }_{\text{N}\text{R}}\) is \(12.4,\) indicating that minute changes in neural responses can have an almost all-or-none effect on HR.

3.6 Simulations

Our simulations with the final PK model indicate that both IV and SC administration of dexmedetomidine doses of 0.5–1 µg/kg resulted in previously reported therapeutic plasma concentrations [7]. The NCA indicated that the model-predicted peak concentrations and AUCs were consistent with previously reported results (ESM) [6]. The time to achieve steady-state plasma concentrations with a constant-rate SC infusion was 15–20 h, which is considerably longer compared with 3 h after IV dosing (Fig. 4).

Fig. 5
figure 5

Simulated norepinephrine concentration profiles in plasma using the final mechanism-based model and 8-h (upper row) and 48-h (lower row) intravenous (a, c) and subcutaneous (b, d) continuous infusions of varying dexmedetomidine doses

Our model results show that SC administration of dexmedetomidine results in a milder inhibition of NE release from the release compartment at the same dose concentrations, which translates into an attenuation of an otherwise pronounced haemodynamic effect with the IV route. A dose concentrations of 4 µg/kg/h or above can result in lowering NE levels in the release compartment below 0.1 nM, leading to very low circulating NE levels, both in IV as well as SC administration (Fig. 5). The simulations with our final SAP and DAP models show a dynamic interplay between the hypotensive and hypertensive effects, depending upon the concentrations of NE and dexmedetomidine. Subcutaneous dosing of 0.25–1 µg/kg/h does not evoke notable hypertensive effects, and a unanimously downward trend can be seen in blood pressure simulations at these dose concentrations. After IV dosing, initial high dexmedetomidine concentrations result in a short hypertensive effect (\({E}_{\text{B}\text{P},\text{H}\text{T}\text{N}}\)), which is swiftly followed by a hypotensive effect (\({E}_{\text{B}\text{P},\text{H}\text{P}\text{N}}\)) due to reduced NE release (ESM). Furthermore, our simulations with IV doses suggest that the hypertensive effect is minimal at lower doses but increases considerably with increasing doses and surpasses the hypotensive effect at dose concentrations above 2 µg/kg. However, because of the rapid disappearance of the hypertensive effect of dexmedetomidine, the hypotensive response due to reduced NE release causes blood pressure to fall (ESM).

The addition of an effect from changes in central neural activity following SAP perturbations helped to simulate the dynamic interplay between these effect components on HR over a range of doses and dosing schemes and our simulations suggest a dose-dependent decline. The effects are clearly seen to be more pronounced for the IV route over the SC route and the slope of decline after IV administration is much steeper as compared with SC administration, owing to high initial concentrations of dexmedetomidine (ESM).

4 Discussion

Our primary aim was to develop semi-mechanistic models to characterise dexmedetomidine concentration– and response–time courses after SC and IV dosing. Using previously collected data and a sequential PK–PD modelling approach, we were able to describe dexmedetomidine absorption kinetics and further disposition together with its sympatholytic and haemodynamic effects. Our models were able to capture the data trends with relatively high precision, and the model predictions show good agreement with the actual observations. Internal evaluation confirmed the predictive performance of the models and the accuracy of the final model parameter estimates, with narrow confidence intervals.

A first-order absorption model has previously been used for nasal dexmedetomidine administration [9], but it did not accurately explain the time course and trends in dexmedetomidine absorption after SC dosing. The original data indicated fast absorption after SC dosing, followed by a long equilibration phase. A wide range of absorption models was tested, but a more mechanistic solution was needed. Dexmedetomidine is highly lipid soluble [29], which pointed towards developing an absorption model with an additional SC fat compartment describing a slower absorption process. As the SC fat compartment is perfused with its own blood supply, we assumed direct absorption from the fat layer into peripheral circulation with a slow absorption rate constant. A redistribution model from the SC fat layer into the injection depot before absorption into plasma was also tested. The former approach was chosen for being physiologically more plausible, and it was also supported by the data. The final absorption model showed an adequate fit for the data at hand. Furthermore, results from our NCA simulations indicate that our model predicted similar AUC as found previously [6].

We postulated that SC fat could act as a slow-release reservoir of dexmedetomidine during prolonged exposure. The results from model simulations indicate that at higher SC doses, the peak concentration in plasma is postponed, while the total systemic exposure remains nearly unchanged. Simulations also suggest that plasma concentrations remain higher for longer time after SC dosing compared with IV dosing. Our results propose that constant-rate infusions of 8 and 48 h and at dexmedetomidine dose rates between 0.25 and 0.5 µg/kg/h would be adequate to reach previously suggested therapeutically effective plasma concentrations of 0.3–0.8 ng/mL [7] after both IV and SC administration. This implies that adequate palliative sedation might be achieved using low dosing rates and the convenience of an extravascular SC route of administration. However, it needs to be highlighted that according to our results, the time needed to achieve steady-state plasma concentrations was higher in SC (15–20 h) vs IV (5–10 h) administration, and this may be worth a consideration in clinical settings.

Our final PK model shows consistency with previous reports. A two-compartmental PK model was best suited to explain the disposition kinetics of the drug, which is in line with previous population-based modelling reports [9]. Some previous studies have used three-compartmental models, but our initial effort showed a two-compartment model more stable with more precise parameter estimates and smaller standard errors. Furthermore, parameter estimates from our final PK model agree with previous studies reporting elimination clearance estimates between 41 and 57 L/min [9, 30]. Interestingly, the estimate of clearance from our study with SC dosing is completely in line with a previous report using an intranasal route of administration (i.e. 52.6, 44.0 and 44.5 L/min after IV, intranasal and SC dosing, respectively) [9, 13].

The estimation of bioavailability parameters in our final model demonstrated that the fat compartment receives approximately 55% of the drug from the depot, and of the remaining 45%, a further 79% reaches the systemic circulation via fast absorption. The resulting net bioavailability of dexmedetomidine via the SC route value is 89% (\(F= 0.89\)), which is slightly higher than our previous estimation using NCA. Bioavailability of dexmedetomidine after intranasal, intramuscular, oral and buccal administration has been investigated previously [6, 31, 32]. Our current findings suggest approximately similar bioavailability of dexmedetomidine after SC and intranasal administration and confirm that the bioavailability is sufficient for clinical usefulness also after SC administration. Compared to intranasal dosing, SC dosing seems more appropriate for palliative care, as the elimination is much slower.

We measured NE with dense sampling to characterise the effect on sympathetic nervous system after dexmedetomidine dosing. We assumed that for modelling purposes, NE levels at the sites of release could be used as a surrogate biomarker of sympatholysis. An indirect response model [33] could capture the sympatholytic effect of dexmedetomidine on NE release with more precise model estimates, relative standard errors and narrower confidence intervals compared with previous results [9]. Our simulations suggest doses above 4 µg/kg nearly abolish plasma NE levels, which may have clinical implications to be considered when using dexmedetomidine at high SC infusion rates. However, it should be emphasised that much lower doses are clinically effective, indicating that the SC dosing route is a safe alternative for palliative sedation.

Dexmedetomidine exerts biphasic effects on blood pressure [34,35,36]. At high initial concentrations and especially after IV dosing, a fast hypertensive effect predominates due to an action on vascular α2-adrenoceptors [37, 38]. Sympatholysis develops somewhat later when the drug distributes into the CNS to regulate sympathetic nervous activity [7, 36, 39]. Our model estimate for \({\text{E}\text{C}}_{50}\) of dexmedetomidine for causing a rise in blood pressure was 2.45 ng/mL. It is therefore understandable that the initial hypertensive response is much more prevalent after IV than after SC dosing as the former leads to sudden high plasma concentrations of dexmedetomidine. Considering fragile patients requiring palliative care may have concomitant cardiovascular problems, SC administration should thus attenuate the haemodynamic challenges associated with dexmedetomidine administration.

The role of NE in the modulation of HR after dexmedetomidine administration has been previously documented and also utilised in PD modelling of the HR response after IV dosing [9]. Yoo et al. studied extravascular dosing via intranasal administration and used an indirect response model to estimate NE effect on HR by employing plasma NE concentrations and the NE biophase [9], but high between-subject variability was seen. A recent study of Colin et al. used dexmedetomidine plasma concentrations as a surrogate marker for HR and showed that HR model estimates followed predicted dexmedetomidine plasma concentration linearly [11]. Their results indicate that the lower boundary for HR variability after IV dexmedetomidine is 22 beats/min, but again the inter-individual variability was quite high (48%). We also first tested an indirect response model to account for the NE effect on HR, but a clear misspecification was observed in the model predictions. Changes in blood pressure are known to result in modulation of the baroreceptor reflex and changes in its activity help to maintain haemodynamic homeostasis. The effect of the baroreflex attributable to dexmedetomidine pharmacokinetics was recently shown to account for initial vasoconstriction after dexmedetomidine dosing [13].

We hypothesised, based on this finding [13] and our results, that changes in central neural reflexes that drive the HR effect and depend on baroreflex activity are essential for accurate prediction of the HR time course. The HR model was significantly improved after including an effect component dependent on blood pressure levels, which mediates the effects of the central neural activity. Simulations with our final HR model indicate that the dexmedetomidine has a more pronounced effect on HR after intravascular dosing, and SC dosing seems to be also safer in this regard. Our results agree with a previous finding demonstrating that extravascular dosing may avoid strong acute haemodynamic changes seen after intravascular use.

Our study has several limitations. Although we aimed to implement mechanistic modelling solutions, many assumptions regarding underlying physiological processes had to be made during the development of our models for simplification. These assumptions override the deeply complex and multi-layered nature of the effects of dexmedetomidine on haemodynamic parameters and subjective effects. Additionally, we could not achieve a stable epinephrine model because of over-parameterisation issues because of the large amount of concentrations below the limit of quantification. However, our current models show reasonable performance with low standard errors and narrow confidence intervals. This indicates that although data from only ten subjects were studied, our dense sampling could provide enough information for model development. Finally, the data used in the modelling were collected from healthy volunteers, thus our findings should be translated to clinical practice only with appropriate caution.

5 Conclusions

Our semi-mechanistic models were able to capture the PK and PD profiles of both intravenously and subcutaneously administered dexmedetomidine, describe the inhibition of NE release, and evaluate the haemodynamic and subjective drug effects. The PK model provides a mechanistic exploration of the absorption process and disposition after SC dexmedetomidine with excellent precision. Endogenous mediators were utilised in an indirect response model for inhibition of NE release to control the diverse effects of the drug. Norepinephrine effects on HR, SAP and DAP were linked via a biophase compartment, which reflects the effects of endogenous mediators and their access to the heart and blood vessels. Our current model provides accurate and clinically plausible predictions, which can provide additional information for drug therapy used in palliative care.