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. Electronic supplementary material The online version of this article (10.1007/s10928-018-09617-0) contains supplementary material, which is available to authorized users.


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][2][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][17][18][19][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 concentrationtime 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 interindividual 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.

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 lg kg -1 , 1 ? 0.7 lg kg -1 and 10 ? 7 lg 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.
C p and C t denote drug concentration in central (plasma) and peripheral compartments. V c , V t , Cl and Cl d 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.
R is cortisol concentration, k out the fractional turnover rate and R eq stands for the expression in Eq. 5. The oscillatory behaviour of turnover rate was modelled by Eq. 3. Note that in the following x = 2p/24 h -1 .
k avg is positive and corresponds to the average turnover rate. t 0 is the phase shift between -12 and 12 h, and a, a number between 0 and 1, describes the amplitude of the oscillations as a proportion of k avg . Choosing a 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 I max is maximum inhibitory capacity, IC 50 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 p, eq , the cortisol response is given by 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 R eq 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 p (t ij ) and R(t ij ) 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 ij . c C p ðt ij Þ andRðt ij Þ are the predicted concentrations for subject i at time point t ij . e ij as well as s ð1Þ ij and s ð2Þ ij were assumed to be normally distributed with zero mean and respective standard deviations e as well as r 1 and r 2 .

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., logðhÞ $ Log-Normal ðl; s 2 Þ ; l $ Normal ðm 0 ; 1Þ ; s $ Student-t ð4; s 0 ; 0:25Þ ; where m 0 and s 0 are prior parameters and h stands for Cl, Cl d , V c and V t . Some parameters in the cortisol turnover model were modelled with correlations as where X = LDL T , D ¼ diagðs 2 1 ; s 2 2 ; s 2 3 ; s 2 4 ; s 2 5 Þ and L is a lower-triangular matrix. In this representation the matrix D contains the variances and LL T is the correlation matrix. In addition n $ Normal ðl n ; s n Þ ; Hyperparameters l in Eqs. 10, 11 and the diagonal elements of D as well as s in Eq. 11 are distributed as where m 0 and s 0 are prior parameters, v = 2.5 for hyperparameters related to a as well as t 0 , and v = 1 otherwise. The three-parameter Student-t distributions used in Eqs. 9, 12 for non-negative s are truncated distributions. The correlation matrix LL T 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-errormodel standard deviations e and r 1 and r 2 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 lg 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 lg 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 lg 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-RubinR 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.

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).
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 I max 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)/k out ) 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).
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 d (Table 1) as well as a, t 0 and IC 50 ( Table 2). Uncertainty in IIV standard deviations was comparably larger, with 95% 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 d were largest (less than 0.9). Predicted IIV for the parameters in the cortisol model was larger, particularly for parameters k avg , a and IC 50 , 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 k avg , a, t 0 , k out and IC 50 . 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 k avg and a, k out and a, and IC 50 and a, 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.

Simulation of two overnight DST protocols
We simulated time courses for healthy horses at four different dose levels (10, 20, 30 and 40 lg 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.
The variability in cortisol response around the threshold value (10 lg L -1 ) proposed in Dybdal et al. [1] was high after administration of 10 lg kg -1 but decreased steadily for increasing dose levels (Fig. 4). After administration of the 40 lg 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 lg 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 Values reported as median and 95% credible interval 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 lg 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 lg 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).

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 a 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 k out 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 k avg , a and IC 50 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 d , a, t 0 and IC 50 . 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 d . Parameters a and t 0 are mostly informed by baseline measurements, as their influence is diminished by dexamethasone administration. However, as reported in the  [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 lg L -1 in blood plasma 19 or 24 h after administration of 40 lg 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 lg 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][43][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 nonlinear mixed effects modelling allowed estimation of variation between individuals finding IIV in parameters k avg , a and IC 50 . 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 lg 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.
Acknowledgements Felix Held was supported by a research student grant from Grünenthal GmbH. Parameter inference was performed on resources at Chalmers Centre for Computational Science and Engineering (C3SE) provided by the Swedish National Infrastructure for Computing (SNIC). The experimental costs and drug analyses were funded by the Swedish-Norwegian Foundations for Equine Research. This work was also partially funded by the Swedish Foundation for Strategic Research.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creative commons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Cortisol concentration for constant dexamethasone concentration
The solution to the cortisol model equations (Eqs. 2-4) can be found analytically for constant dexamethasone concentration C p = C p, eq . Mathematically, Eq. 2 describes a first-order differential equation driven by an external oscillation (the turnover rate). Using methods described by Farkas [45], it can be shown that there exists a unique oscillating solution to Eq. 2. This solution will have the same period as the driving turnover rate and any other solution will tend towards this oscillating solution. Assume that the model equations are solved by where x = 2p/24 h -1 and t 0 is the phase shift in the turnover rate (Eq. 3). Solving for a, b and c gives To determine characteristic parameters for the equilibrium oscillation, it is more convenient to write Eq. 13 in the form R eq ðtÞ ¼ A þ B cosðx ðt À CÞÞ where A, B and C are as in Eq. 6. This representation can be found by standard arguments for the sum of trigonometric functions. It can be seen from these equations that A and B depend on the drug concentration, while C does not.

Determination of prior parameters
Prior beliefs about means and variances for hierarchically distributed parameters were quantified through a metaanalysis of Ekstrand et al. [12]. Therein, point estimates of parameter values and their relative precision (coefficient of variation, CV) in percent were reported. If h is the reported point estimate and CV its relative precision then standard errors r were calculated as r = (CVÁh)/100. To arrive at an estimate for the population mean and population variance for a parameter, e.g., Cl, weighted averages of the point estimates were formed. Each estimate was weighted by its squared standard error. Let h i and r i be the estimate and standard error determined for the ith horse. Then wherel andx 2 are the population mean and variance, respectively. Parameters were modelled hierarchically either through a log-normal distribution (e.g., Cl), a (scaled) logit-normal distribution (for a, t 0 and I max ) or a normal distribution (for n). For all parameters the prior location parameters m 0 and s 0 were determined by matching mean and variance of the respective distribution with those in Eq. 15.
For a normal distributionl andx 2 can be used directly. Analytical formulas for the mean and variance of the lognormal distribution exist. It holds that where m and s are the location and scale parameters of the log-normal distribution. l and x 2 were matched withl and x 2 in Eq. 15 to determine m and s. Values can be found in the table below. In case of the (scaled) logit-normal distribution no closed-form solutions for mean and variance exist. Therefore, moments were calculated by numerical integration using SciPy's fixed_quad function [46]. Numerical optimization was used to determine location and scale parameters. The difference between numerically approximated moments and those in Eq. 15 was used as a target function. SciPy's fsolve, which implements methods from MINPACK [47], was used for numerical optimization. The resulting optimized parameters were m and s. Values can be found in the table below.
The parameter controlling the average turnover rate k avg was not directly given by Ekstrand et al. [12]. Instead, the parameter R ij = k avg, ij /k out, i for individual i and occasion j was given. 1 Note that the parameter k avg was allowed to vary by experimental occasion in Ekstrand et al. [12]. In the model presented herein, the parameter k avg was modelled with IIV instead. Individual-and occasion-specific parameter estimates for k avg, ij could easily be recovered as k avg, ij = R ij Ák out, i . Relative precisions were used to calculate the standard errors r R ij and r k out;i of each individual-/ occasion-dependent estimate. Standard errors for k avg, ij were then approximated by propagation of uncertainty [48] r k avg;ij % k avg;ij ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Summing over all individuals and occasions, weighted averages, as in Eq. 15, were formed. These were matched with mean and variance of a log-normal distribution to get location m and scale s for k avg . Values can be found in the table below.
Location parameters m and s for prior distribution of hyper-parameters l and s.

Calculation of simulation of sensitivity and specificity
By sampling j = 1,…, M individual parameter vectors, called h i j for the ith set of hyper-parameters, sensitivity Se (i) and specificity Sp (i) can be calculated as Here, U is the cumulative distribution function of the standard normal distribution, andRðh ðiÞ j Þ is the predicted cortisol measurement obtained in the DST. In the case of sensitivity, sick subjects were simulated andRðh ðiÞ j Þ was therefore predicted at baseline without drug. For specificity, healthy subjects were simulated andRðh ðiÞ j Þ was predicted after administration of 40 lg kg -1 dexamethasone. Repeating this procedure for many sets of hyper-