Modelling of oscillatory cortisol response in horses using a Bayesian population approach for evaluation of dexamethasone suppression test protocols.

Cortisol is a steroid hormone relevant to immune function in horses and other species and shows a circadian rhythm. The glucocorticoid dexamethasone suppresses cortisol in horses. Pituitary pars intermedia dysfunction (PPID) is a disease in which the cortisol suppression mechanism through dexamethasone is challenged. Overnight dexamethasone suppression test (DST) protocols are used to test the functioning of this mechanism and to establish a diagnosis for PPID. However, existing DST protocols have been recognized to perform poorly in previous experimental studies, often indicating presence of PPID in healthy horses. This study uses a pharmacokinetic/pharmacodynamic (PK/PD) modelling approach to analyse the oscillatory cortisol response and its interaction with dexamethasone. Two existing DST protocols were then scrutinized using model simulations with particular focus on their ability to avoid false positive outcomes. Using a Bayesian population approach allowed for quantification of uncertainty and enabled predictions for a broader population of horses than the underlying sample. Dose selection and sampling time point were both determined to have large influence on the number of false positives. Advice on pitfalls in test protocols and directions for possible improvement of DST protocols were given. The presented methodology is also easily extended to other clinical test protocols.

Introduction

Dexamethasone and other glucocorticoids are commonly used in equine medicine for the treatment of diseases and clinical testing, e.g., the dexamethasone suppression test (DST) [1-3]. In healthy horses, dexamethasone suppresses the cortisol response [4, 5]. The mechanism is challenged in horses affected by pituitary pars intermedia dysfunction (PPID) [6]. PPID is an age-related degenerative disease leading to a loss of dopaminergic neurons affecting the pars intermedia of the pituitary gland [7]. Similarities of PPID to Parkinson’s disease in humans have been found [8]. The most prevalent clinical signs of PPID in horses are hair coat abnormalities, laminitis and muscle atrophy [9]. The prevalence of PPID in horses aged more than 15 years is 21% [10].

The idea behind the DST is to observe whether the cortisol response is suppressed below a threshold after dexamethasone administration [1]. Different test protocols have been published under the name overnight DST [1, 11]. They differ in dexamethasone administration time as well as waiting-time until cortisol measurement. Both are designed as single-point observation test protocols which have limitations, since DST results may indicate presence of PPID in healthy individuals due to a circadian rhythm in cortisol production in addition to inter-individual variability [12]. In horses, plasma cortisol concentration displays an apparently symmetric circadian variation with peak concentration in the morning and nadir concentration in the afternoon/evening [13, 14].

Pharmacokinetic (PK) and pharmacodynamic (PD) modelling is done to quantify the relationship between a physiological response and drug exposure [15]. The disposition of dexamethasone in horses has been characterized in several studies [12, 16–20]. The relationship between dexamethasone and cortisol response has also been characterised by turnover modelling with inhibition of a constant [21] as well as an oscillating turnover rate [12]. Those studies focused on modelling individual concentration–time as well as response-time courses. Modelling studies have not estimated potential variability between animals.

A non-linear mixed-effects (NLME) approach allows simultaneous regression of all individuals and time courses [22]. A merit of this technique is the estimation of inter-individual variability (IIV) directly from data [23]. The use of mixed effects models has historically not been extensively used in veterinary science but recently has been given more attention [24]. Mixed-effects models can be formulated as hierarchical Bayesian models [25]. The Bayesian approach allows for incorporation of prior knowledge as well as modelling of all sources of variability and uncertainty. This allows for explicit propagation of uncertainty in parameter estimates and residual variability to predictions made using the final adjusted model [22]. Including all sources of uncertainty and simulating predictions from the full model, as is straight-forward in the Bayesian approach, increases the credibility of the predictions. Prior information about parameter estimates and variability is available from an earlier study [12].

In this study, we sought to analyse data from a previous study [12] by means of a NLME approach to investigate the IIV. We then used the adjusted model to scrutinize DST protocol designs, define weaknesses in two proposed test protocols and give directions for test improvement.

Materials and methods

Experimental setup and analytical method

Six Standardbred horses (four mares and two geldings) 6–20 years old and weighing 430–584 kg were included in the study and assigned to a randomised crossover design including four treatments and four periods. Each treatment started with an intravenous bolus dose immediately followed by 3 h of constant rate infusion of dexamethasone 21-phosphate disodium salt (Dexadreson 2 mg mL−1, Intervet AB, Sollentuna, Sweden). The dose levels were (bolus + infusion) 0.1 + 0.07 μg kg−1, 1 + 0.7 μg kg−1 and 10 + 7 μg kg−1 dexamethasone. For the control level 0.9% saline was used. Before the bolus dose (time = 0) a pre-dose blood sample was drawn. Additional blood samples were drawn during and after infusion at hours 1, 2, 3, 4, 5, 6, 9, 12, 18, 24, 36 and 48. A minimum of a 1 week wash-out period was allowed between drug treatments. The study was approved by the Ethics Committee for Animal Experiments, Uppsala, Sweden (C333/11). Total plasma dexamethasone and cortisol concentrations were analysed and quantified using Ultra High Performance Liquid Chromatography-Tandem Mass Spectrometry (UHPLC-MS/MS). The analytical method was described before elsewhere [21].

Dexamethasone exposure

A two-compartment model (Eq. 1, Fig. 1a) was fitted to experimental dexamethasone-time course data.

Fig. 1

a Dexamethasone disposition model after the bolus + constant rate infusion regimens. b Conceptual model of the cortisol turnover with an oscillating baseline kin(t)

C and C denote drug concentration in central (plasma) and peripheral compartments. V, V, Cl and Cl denote, respectively, the central and peripheral volumes, plasma clearance and inter-compartmental distribution parameter. Inf(t) represents the constant rate infusion regimen and D is the bolus dose administered at time t  =  0.

Cortisol turnover

Cortisol was modelled by a turnover model (Eqs. 2–4, Fig. 1b).

R is cortisol concentration, kout the fractional turnover rate and Req stands for the expression in Eq. 5. The oscillatory behaviour of turnover rate was modelled by Eq. 3. Note that in the following ω = 2π/24 h−1.

kavg is positive and corresponds to the average turnover rate. t0 is the phase shift between − 12 and 12 h, and α, a number between 0 and 1, describes the amplitude of the oscillations as a proportion of kavg. Choosing α this way ensures positivity of the turnover rate for all choices of parameters. The period was fixed at 24 h. The inhibitory dexamethasone mechanism function was modelled as

where Imax is maximum inhibitory capacity, IC50 the potency of dexamethasone and n is a Hill exponent. This model is a modified version of the single cosine model presented in [26].

Cortisol concentration under constant dexamethasone exposure

An oscillating turnover rate leads to oscillating cortisol concentration. Keeping dexamethasone exposure in Eq. 4 constant at a fixed concentration C, eq, the cortisol response is given by

where

A describes the average cortisol response, B the amplitude and C the phase shift of the oscillation. The model predicts only changes in the average cortisol response and amplitude due to changes in dexamethasone exposure. A derivation of Req in Eq. 5 as well as for A, B and C in Eq. 6 can be found in the Appendix. The ideas are similar to the calculations presented in Krzyzanski et al. [27].

Residual error variance model

Kinetic data was modelled on a log scale. For the drug exposure model in Eq. 1 it was assumed that

For the cortisol response model in Eq. 2 a combined error model with proportional and additive error was assumed. This was described by

Here, C(t) and R(t) are the jth measurement of the plasma concentration of dexamethasone in the central compartment and cortisol, respectively, measured for subject i at time point t. and are the predicted concentrations for subject i at time point t. e as well as and were assumed to be normally distributed with zero mean and respective standard deviations ε as well as and .

Statistical parameter model

IIV was modelled by making the following assumptions about the distribution of the parameters in Eqs. 1–4. The process of deciding which parameters were modelled with correlation is described in the Supplementary.

All parameters involved in the description of dexamethasone exposure were modelled independently log-normally distributed, i.e.,

where m0 and s0 are prior parameters and θ stands for Cl, Cl, V and V. Some parameters in the cortisol turnover model were modelled with correlations as

where Ω = LDL, and L is a lower-triangular matrix. In this representation the matrix D contains the variances and LL is the correlation matrix. In addition

Hyperparameters μ in Eqs. 10, 11 and the diagonal elements of D as well as τ in Eq. 11 are distributed as

where and are prior parameters, v = 2.5 for hyperparameters related to α as well as t0, and v = 1 otherwise. The three-parameter Student-t distributions used in Eqs. 9, 12 for non-negative τ are truncated distributions.

The correlation matrix LL was assumed to be distributed following a LKJ distribution [28] with concentration parameter 2. This is a prior for correlation matrices where samples resemble the identity matrix more closely for concentration parameters closer to 1. Residual-error-model standard deviations ε and and were assumed to be positive and were given half-Cauchy prior distributions [29] with scale 2.5. Prior parameters were estimated from a meta-analysis of Ekstrand et al. [12] as described in the Appendix.

Analysis of the dexamethasone suppression test protocol

We simulated two different overnight DST protocols. Each consisted of a dexamethasone administration time and a sampling time on the following day. Cortisol concentration is analysed in the sampled blood plasma and the result of the DST is positive if concentration is above a prescribed threshold. The protocols analysed more closely are described in Dybdal et al. [1] (protocol A) and Frank et al. [11] (protocol B). Both protocols assume administration of 40 μg kg−1 of dexamethasone. The protocols differ in administration route. Protocol A assumes intramuscular (im) administration whereas protocol B assumes intravenous (iv) administration. Test starting times were at 9.00 a.m. (protocol B) and 5 p.m. (protocol A). Plasma sampling times for determination of cortisol concentration were after 19 h (protocol A) and 24 h (protocol B), respectively. In both protocols, the test is positive (indicating sick individuals) if measured cortisol concentration is above a threshold of 10 μg L−1.

The DST protocols were analysed in light of two different aspects. First, we performed a Monte Carlo study to visualise cortisol time courses for horses subjected to each protocol. A sample of 10,000 horses was simulated from the adjusted model. For this, residual variance parameters and hyper-parameters (N = 1000) were taken from the estimated posterior parameter distribution. Then, individual parameters (N = 10) were simulated from hyper-parameters and the distributions in Eqs. 9–11. Dexamethasone- and cortisol time courses were simulated under the two test protocols for the new subjects using Eqs. 1–4 as well as the measurement equation (Eq. 8). The investigated protocols assume administration of 40 μg kg−1 dexamethasone and the aim of this simulation was to determine whether this amount is necessary or if lower doses could be sufficient. Predicted cortisol concentration at sampling time was then used for further analysis.

These concentrations were then used in a second step to investigate both the sensitivity of the test, i.e., the probability that the test is positive for a sick subject, as well as the specificity of the test, i.e., the probability that the test is negative for a healthy subject [30]. The distributions of sensitivity and specificity were simulated through a combination of Monte Carlo and analytical steps. See the Appendix for the formulas used. In horses with PPID the mechanism for dexamethasone suppression of cortisol is challenged [6]. To quantify sensitivity, simulations from sick horses were needed. We hypothesized that dexamethasone has no suppression effect on sick individuals and therefore these horses were sampled at baseline. The studies reporting protocol A and B [1, 11] determined sensitivity and specificity experimentally and this analysis aimed to investigate if model predicted and experimentally determined values are aligned.

Numerical analysis and parameter estimation

The software Stan version 2.18.0 [31] was used for parameter inference through the interface CmdStan. Stan implements the NUTS sampler [32] that uses Hamiltonian Monte Carlo (HMC) [33] for estimation of the posterior parameter distribution and allows models with differential equations. PK and PD parameters were estimated in two stages. First, PK parameters in Eq. 1 were estimated. Each individual’s PK parameters were then summarised and fixed to the respective conditional mean. In a second stage, the PD parameters in Eqs. 2–4 were estimated. In each stage, four Markov chains were started at random initial parameters around the prior parameter means. Each chain was run for 250 iterations in warm-up and sampling, respectively. This led to a total of N = 1000 samples from the posterior.

The convergence of the HMC algorithm was checked in multiple ways. Numerical divergences during parameter estimation were observed and appropriate choices about parameter distributions were made and Stan settings were tuned to reduce and avoid divergences [34]. The Gelman–Rubin statistic [35] and trace plots were used to ensure proper mixing of the Markov chains. The effective sample size [25] was observed to be at least 10% of total samples size (N = 1000). The energy Bayesian fraction of missing information (E-BFMI) [36] was checked to ensure that the parameter space was properly and efficiently explored. No external validation data was available and therefore internal model checking was performed through posterior predictive checks (PPCs) [25]. These visualisations are similar to visual predictive checks (VPCs) [37]. However, PPCs include parameter uncertainty by simulating the response from the full estimated posterior distribution, whereas VPCs omit this. Estimated parameters were summarised by median and 95% credible intervals (CIs) [25]. Population predicted ranges were calculated as described in the Supplementary.

Results

Regression of experimental time courses

The drug exposure model captured dexamethasone exposure across three orders of magnitude (Figs. 2, S1). In general, within and between subject variability in the dexamethasone time courses were low, which suggests that exposure of dexamethasone does not confound the cortisol response. Observed and regressed dexamethasone-time courses following the three dosing regimens for a representative horse are shown in Fig. 2. The final population model parameters as well as their predicted population range are shown in Table 1. Summaries of individual parameters per horse are reported in the Supplementary (Table S1).

Fig. 2

Dexamethasone and cortisol time course data and model predictions for a representative horse. Different dosing regimens are indicated by colour and respective dosing amounts are shown in the legend. Solid lines are time courses corresponding to the subject’s average parameters. Uncertainty in predicted time-courses is shown by shaded areas representing 95% of uncertainty

Table 1

Estimated population parameters, the IIV standard deviation and predicted population range for the dexamethasone model

Parameter (unit)	Typical valuea		IIV std. dev.a		Predicted population rangea
Cl (L h−1 kg−1)	0.50	(0.43, 0.56)	0.071	(0.0036, 0.22)	0.50	(0.39, 0.63)
Cld (L h−1 kg−1)	0.16	(0.078, 0.33)	0.40	(0.032, 0.86)	0.16	(0.049, 0.55)
Vc (L kg−1)	1.2	(0.74, 1.7)	0.15	(0.013, 0.48)	1.2	(0.62, 2.1)
Vt (L kg−1)	0.76	(0.52, 1.0)	0.14	(0.013, 0.48)	0.76	(0.42, 1.3)

aValues reported as median and 95% credible interval

The observed cortisol response was well captured by the model (Figs. 2, S2). The largest variability in cortisol response was observed in baseline data. Increasing exposure to dexamethasone suppressed both the average cortisol response and the amplitude of oscillations (Fig. 2), which was furthermore captured by the model. The suppression of cortisol by dexamethasone was almost complete, which is also evident from the estimated values of Imax close to 1 (Table 2). The typical potency value was predicted to be 37 ng L−1 varying between 22 and 56 ng L−1, and the half-life of the cortisol response (ln (2)/kout) was predicted to be 2.0 h varying between 1.1 and 3.4 h. Observed and regressed cortisol response-time courses following the three dosing regimens and baseline for a representative horse are shown in Fig. 2. The final population model parameters as well as their predicted population range are shown in Table 2. Summaries of individual parameters per horse are reported in the Supplementary (Table S2).

Table 2

Estimated population parameters, the IIV standard deviation and predicted population range for the cortisol model

Parameter (unit)	Typical valuea		IIV std. dev.a		Predicted population rangea
kavg (μg L−1 h−1)	15	(11, 22)	0.36	(0.15, 0.67)	15	(6.4, 37)
α (unitless)	0.30	(0.13, 0.58)	1.1	(0.77, 1.7)	0.31	(0.035, 0.85)
t0 (h)	− 2.9	(− 4.3, − 1.3)	0.12	(0.0052, 0.37)	− 3.0	(− 5.1, − 0.095)
kout (h−1)	0.34	(0.27, 0.46)	0.20	(0.080, 0.45)	0.34	(0.20, 0.62)
Imax (unitless)	0.99	(0.97, 1.0)	0.41	(0.050, 0.96)	0.94	(0.81, 0.98)
IC50 (ng L−1)	37	(22, 56)	0.56	(0.20, 0.97)	37	(9.1, 130)
n (unitless)	2.3	(1.5, 3.4)	0.39	(0.062, 0.99)	2.4	(1.1, 3.8)

aValues reported as median and 95% credible interval

Uncertainty in parameter estimates and IIV are included in the Bayesian posterior distribution and were readily analysed. Variability in typical values and IIV standard deviation was directly estimable. Predicted population ranges contain variability stemming from both uncertainty in population parameters and variability from the distributional assumptions in the section “Statistical parameter model”. Uncertainty in population parameters was moderate with larger amounts of variability observed in Cl (Table 1) as well as α, t0 and IC50 (Table 2). Uncertainty in IIV standard deviations was comparably larger, with 95% credible intervals often spanning double the magnitude of the median. For parameters in the dexamethasone model IIV was estimated to be small, with most credible intervals reaching close to zero. For Cl the smallest amount of variability (less than 0.2) was estimated whereas estimates for Cl were largest (less than 0.9). Predicted IIV for the parameters in the cortisol model was larger, particularly for parameters kavg, α and IC50, with credible intervals being clearly bounded from below. Estimated IIV for other parameters was considered to be small with most variability in the predicted population range stemming from uncertainty in estimated typical values.

After initial test runs and analysis of possible correlations between individual parameters, correlation estimates were included for kavg, α, t0, kout and IC50. The resulting estimates can be found in the Supplementary (Table S3). Most correlation parameters showed only a slight if any tendency towards the positive or negative. Results suggest possible negative correlations between kavg and α, kout and α, and IC50 and α, but high uncertainty makes the estimates inconclusive.

Simulation of cortisol response versus dexamethasone plasma concentration

Model simulations of the equilibrium dexamethasone-cortisol response relationship (Fig. 3a), and the amplitude of cortisol response (Fig. 3b) with increasing dexamethasone concentrations show that at concentrations around the potency value (about 20 to 50 ng L−1) small changes in dexamethasone plasma concentration produce disproportionally large changes in cortisol response and amplitude. At concentrations well above the potency value (> 100 ng L−1) small changes in dexamethasone plasma concentration produce small changes in cortisol response and amplitude. Note the almost complete suppression of cortisol response and its variability, with increasing dexamethasone concentrations.

Fig. 3

a Plot of the average response (A in Eq. 6). b Plot of the oscillation amplitude (B in Eq. 6). Variability is shown as shaded areas

Simulation of two overnight DST protocols

We simulated time courses for healthy horses at four different dose levels (10, 20, 30 and 40 μg kg−1 dexamethasone administered intravenously) following protocols A and B. The resulting Fig. 4 shows predicted variability in possible time-courses resulting from IIV and parameter uncertainty. Additionally, predicted time courses under protocol A and B for the horses involved in this study are shown.

Fig. 4

Simulated time courses including model-predicted variability for healthy horses undergoing the DST protocols A and B with the modification of varying the test dose of 10, 20, 30 and 40 μg kg−1 dexamethasone. The black dashed line represents the threshold of 10 μg kg−1 proposed by Dybdal et al. [1]. Arrows indicate time for drug administration at 9.00 o’clock (protocol B) or 17.00 o’clock (protocol A) and blood sample 19 h (protocol A) or 24 h (protocol B) after administration. Parameter uncertainty and IIV are incorporated into simulations. This variability is shown by blue shaded areas representing increasing amounts of variability (see legend). Predicted time courses for the horses in this study are shown as solid red lines (Color figure online)

The variability in cortisol response around the threshold value (10 μg L−1) proposed in Dybdal et al. [1] was high after administration of 10 μg kg−1 but decreased steadily for increasing dose levels (Fig. 4). After administration of the 40 μg kg−1 dose variability in cortisol response 19 h after drug administration was lower than after the three lower doses. The model-predicted cortisol response shows that, even after the highest dose, individuals from a healthy population might have cortisol plasma concentrations greater than 10 μg L−1 at sampling time.

Systematic differences between the simulations of protocol A and B are shown in Fig. 4. Following protocol A, dexamethasone is administered at 17.00 o’clock, which coincides with cortisol decrease following the model-predicted circadian rhythm. Protocol B assumes administration at 9.00 o’clock which coincides with a peak in predicted cortisol concentration. As can also be seen from the individual predictions (Fig. 4), cortisol concentrations at administration time are generally lower following protocol A than B. At sampling time, 19 h (protocol A) or 24 h (protocol B) after administration, predicted concentrations for protocol A were consistently lower than for protocol B. However, differences grew smaller with increasing dexamethasone administration. For the highest dose (40 μg kg−1), remaining dexamethasone concentration in blood plasma was predicted to 309 ng L−1 (varying between 101 to 694 ng L−1) at sampling time for protocol A and 155 ng L−1 (varying between 34 to 391 ng L−1) at sampling time for protocol B. Suppressed cortisol concentration after administration of the largest dose of dexamethasone (40 μg kg−1) was lowest roughly between 5.00 o’clock and 12.00 o’clock for protocol A and 19.00 o’clock and 3.00 o’clock for protocol B (Fig. 4).

Prediction of sensitivity and specificity for overnight DST protocols

Sensitivity and specificity were predicted through Monte Carlo simulations. All results are reported as median (95% CI). For protocol A, a sensitivity of 99.1% (93.1% to 99.7%) and a specificity of 83.4% (63.3% to 99.7%) were predicted. For protocol B, a sensitivity of 99.2% (94.2% to 99.8%) and a specificity of 77.8% (52.8% to 96.8%) were predicted. Sensitivities for both protocols turned out equal. Specificity tended towards lower values for protocol B but there is no substantial difference, given the spread of the credible intervals.

Discussion

In this study, we used an NLME approach to investigate the IIV in dexamethasone exposure and cortisol response in horses. The aim was to use the adjusted model to scrutinize DST protocol designs. We analysed two proposed test protocols through simulation, found systematic differences in resulting cortisol time courses and predicted sensitivity and specificity for each protocol.

The PK/PD model used in this study was based on a previously published model (originally published in [26] and used in [12]). For the purpose of this work, some parts of the original model were simplified and re-parametrised. An analytical solution of the cortisol model equations (Eqs. 2–4) was used to initialise the parameter estimation (see Eq. 5). This was necessary since the cortisol response model does not have a constant steady state, which is typically used for initialisation. The calibrated model (Eqs. 1, 2) mimicked the overall tendency of the experimental data well. Individual PK parameter estimates (Table S1) and most PD parameters (Table S2) agree with previously reported values [12]. Parameter α was not directly comparable, since a different parametrisation of the model is used in this study. In comparison to Ekstrand et al. [12], the model in this study does not use a second delay compartment for the cortisol response. This affects individual estimates for parameter kout which were systematically lower in this study, about half the previously estimated values. We found no evidence in the data to motivate the necessity of an additional delay compartment. Parameters kavg, α and IC50 showed evidence of non-zero IIV.

Large improvements in estimated precision for parameter n were achieved in this study. However, considerable uncertainty in some estimated parameters remained, notably Cl, α, t0 and IC50. Dexamethasone exposure data (Figs. 2, S1) suggest a two-compartment model, but the influence of the second compartment is subtle, which explains the large uncertainty in parameter estimates of Cl. Parameters α and t0 are mostly informed by baseline measurements, as their influence is diminished by dexamethasone administration. However, as reported in the results section (Figs. 2, S2) irregular fluctuations in baseline measurements confound the underlying oscillation of the cortisol response, explaining the challenge in precise estimation of these two parameters. Irregular fluctuations in baseline measurements can be explained by handling stress [38] as well as pulsatility in cortisol secretion [39].

The multiple sources of uncertainty, incorporation of IIV as well as the aim to utilize the calibrated model for prediction motivated our use of a Bayesian approach. The resulting posterior parameter distribution includes the general tendency of the estimated parameter values, as well as variability in form of IIV and parameter uncertainty. When IIV is included through the framework of NLME, parameter estimation can get more robust, since parameter estimates for one subject are informed by estimates of other subjects [22]. This can be seen in parameter n, which was estimated with much higher precision in this study than in Ekstrand et al. [12]. This holds for the Bayesian approach as well as the maximum likelihood approach. However, explicit inclusion of parameter uncertainty, when using a Bayesian approach, allows for straight forward model prediction, despite some remaining parameter uncertainty [25].

One clinically relevant application of modelling the cortisol response is evaluation of DST protocols. The overnight DST is commonly used when testing for PPID. Test protocols tend to be straight forward to be usable in both equine-clinics as well as in the field and typically involve only one or two samples. We investigated two overnight DST protocols through simulations, including an analysis of their sensitivity and specificity. Both protocols define a negative outcome as suppressed cortisol response to less than 10 μg L−1 in blood plasma 19 or 24 h after administration of 40 μg kg−1 dexamethasone [1, 11]. In our simulations dexamethasone plasma concentration at sampling time was was predicted to be 101–694 ng L−1 for protocol A and 34–391 ng L−1 for protocol B. This result is consistent with experimental data [16, 17]. Concentration at sampling time was above potency (22–54 ng L−1) for most individuals. However, there is a larger risk for dexamethasone concentration to fall below the potency value when following protocol B. Combined with the IIV in cortisol response, not all healthy horses will fall below the proposed threshold (Fig. 4).

All simulations were made under the assumption of intravenous administration of dexamethasone. Dybdal et al. [1] used intramuscular administration of dexamethasone. As shown by Soma et al. [17], dexamethasone plasma concentrations reach peak concentration after intramuscular administration with a slight time delay (about 15 min) compared to intravenous administration. However, concentrations at times 19 + h after administration, relevant for DST protocols, are comparable after intravenous as well as intramuscular administration. We therefore decided that it is reasonable to compare our simulations with the data published by Dybdal et al. [1].

For the purposes of model estimation only data from healthy horses was available. To simulate sensitivity for DST protocols A and B, we hypothesized that sick horses do not respond to dexamethasone administration. Their responses were therefore simulated at baseline. This was motivated by the known mechanism that cortisol suppression through dexamethasone is challenged in horses affected by PPID [6]. PPID also seems to lead to a reduction in circadian rhythm and an increase in average baseline concentrations [1]. However, numerical values quantifying these latter differences for the race of horses used in this study were not available. Incorporating these changes into the hypothesis stated above would have increased the predicted sensitivity, which was close to 100% already. To the best of our knowledge, there is no finding about differences in dexamethasone kinetics between healthy horses and those affected by PPID.

Sensitivity and specificity were determined through simulation for both protocols. The predicted specificity for protocol A, Dybdal et al. [1] reporting 100% specificity, was 63.3% to 99.7% and for protocol B, Frank et al. [11] reporting 76% specificity, it was 52.8% to 96.8%. The consistency between experimentally determined results and our simulations makes us confident that our model captured the overall cortisol response in healthy horses well. Predicted sensitivity for protocol A, Dybdal et al. [1] reporting 100% sensitivity, was 93.1% to 99.7% and for protocol B, Frank et al. [11] reporting 65% sensitivity, it was 94.2% to 99.8%. An explanation for the large discrepancy between our simulation and the experimental result for protocol B might be that our hypothesis only applies to a population of horses with clinically advanced PPID [40]. Such a population was considered by Dybdal et al. [1] and simulated sensitivity for protocol A is in agreement with experimental results. Future work might consider simulations of horses in which maximum suppression is decreased, imitating a mechanism that still reacts to dexamethasone suppression but at a reduced intensity.

DST protocols A and B use different starting times (afternoon vs. morning) and different waiting times until cortisol sampling (19 h vs. 24 h). Simulations show that it can be more effective to administer dexamethasone during a descending phase in cortisol’s circadian rhythm (Fig. 4). Also, cortisol’s circadian rhythm can interfere with the effectiveness of dexamethasone suppression (Fig. 4, protocol B, upswing after suppression coincides with natural upswing in cortisol production). However, no substantial difference in test outcome between protocols A and B was found in this study. A finding also consistent with a study by Sojka et al. [41] that investigated the influence of starting time on the DST and found no statistically significant difference between starting in the morning or the afternoon.

Simulations showed that a dose of 40 μg kg−1 dexamethasone is necessary to conduct the current test protocols, with lower doses leading to an increased number of false positives (Fig. 4). Small changes in dose around the potency value result in large reduction of cortisol average baseline and amplitude (Fig. 3), but this diminishes quickly for concentrations above the potency value. This has implications on dose selection and sampling time points in the DST. A dose-increment and consequently increased plasma exposure will decrease the number of false positive test results in healthy individuals. It is important to consider that negative adverse effects from glucocorticoids, e.g., hyperglycaemia, hyperinsulinemia and possibly laminitis, are assumed to increase with higher doses [42-44]. Shorter waiting time until sampling could further decrease the number of false positives. Cortisol is suppressed quickly by dexamethasone and peak suppression is reached earlier than 19 + h after administration (Fig. 4). However, this could mean sampling in the late evening or at night (between 19.00 o’clock to 3.00 o’clock for protocol B), which might be unacceptable in the field. Protocol A has better potential for improvement since suppression is at its maximum during early morning hours (between 5.00 o’clock to 12.00 o’clock for protocol A). Existing protocols show that practical applicability is an important factor in their design.

Ideally, serial sampling would be used. An extended sampling protocol with 2–3 additional samples would provide more information about cortisol behaviour and more reliable test results. Additional collection of an unaffected cortisol baseline would allow estimation of the typical cortisol response for the particular horse and allow for an increase in the predictive power of the test. However, this is not a suitable strategy for clinical routine, especially not in the field and the cost would also be accepted in lower extent by the owners. Given the challenge of producing a protocol that clinicians can perform, 2–3 samples after drug administration is the maximum to be collected. Therefore, some uncertainty in interpretation of test results remains and the diagnosis must be based not only on the test outcome.

Conclusion

Our study presents an improved model structure and parameter estimates for cortisol concentration in horses during intervention with dexamethasone. The use of non-linear mixed effects modelling allowed estimation of variation between individuals finding IIV in parameters kavg, α and IC50. Using a Bayesian approach allowed straight-forward propagation of uncertainty to simulations. The adjusted model was successfully used to scrutinize clinical test protocols through simulation. The model output and simulations indicated the importance of dose selection with doses below 40 μg kg−1 performing unfavourably. Sampling time was also found to be of importance and simulations showed that waiting times in the window 10 to 17 h could improve test performance. In addition, it was found that administrating dexamethasone in synchronisation with the down-swing in cortisol’s circadian rhythm can allow for a slight prolongation in waiting time.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary material 1 (DOCX 7139 kb)

Parameter θ	Cl	Cl d	V c	V t	k avg	α	t 0	k out	I max	IC 50	n
m 0	− 0.71	− 1.7	− 0.21	− 0.15	3.0	− 0.82	− 0.71	− 0.56	1.6	− 4.1	0.75
s 0	0.064	0.46	0.27	0.16	0.46	0.98	0.10	0.26	0.34	0.73	0.36