A two-step deconvolution-analysis-informed population pharmacodynamic modeling approach for drugs targeting pulsatile endogenous compounds

Pharmacodynamic modeling of pulsatile endogenous compounds (e.g. growth hormone [GH]) is currently limited to the identification of a low number of pulses. Commonly used pharmacodynamic models are not able to capture the complexity of pulsatile secretion and therefore non-compartmental analyses are performed to extract summary statistics (mean, AUC, Cmax). The aim of this study was to develop a new quantification method that deals with highly variable pulsatile data by using a deconvolution-analysis-informed population pharmacodynamic modeling approach. Pulse frequency and pulse times were obtained by deconvolution analysis of 24 h GH profiles. The estimated pulse times then informed a non-linear mixed effects population pharmacodynamic model in NONMEM V7.3. The population parameter estimates were used to perform simulations that show agonistic and antagonistic drug effects on the secretion of GH. Additionally, a clinical trial simulation shows the application of this method in the quantification of a hypothetical drug effect that inhibits GH secretion. The GH profiles were modeled using a turnover compartment in which the baseline secretion, kout, pulse secretion width, amount at time point 0 and pulse amplitude were estimated as population parameters. Population parameters were estimated with low relative standard errors (ranging from 2 to 5%). Total body water (%) was identified as a covariate for pulse amplitude, baseline secretion and the pulse secretion width following a power relationship. Simulations visualized multiple gradients of a hypothetical drug that influenced the endogenous secretion of GH. The established model was able to fit and quantify the highly variable individual 24 h GH profiles over time. This pharmacodynamic model can be used to quantify drug effects that target other endogenous pulsatile compounds. Electronic supplementary material The online version of this article (doi:10.1007/s10928-017-9526-0) contains supplementary material, which is available to authorized users.


Introduction
Many drugs exhibit wanted and/or unwanted effects on the secretion of endogenous pituitary hormones such as growth hormone (GH), prolactin or luteinizing hormone (LH). These hormones are secreted in 'bursts' or pulses and, particularly in the case of GH, vary in amplitude, time of secretion and may be influenced by circadian rhythmicity and sleep [1][2][3][4]. When studying drug effects on the secretion of GH, a non-compartmental analysis on the plasma concentration-time profile is commonly performed. This results in mean and maximum plasma GH concentrations over a specified time interval or the area under the plasma GH concentration-time curve (AUC) to be used for comparison between groups [5][6][7]. Importantly, by reducing the complex concentration-time profile of GH to summary statistics, the use of time as a continuous variable is lost. Therefore, the analysis results are highly dependent on the timeframe of observations and sampling interval, e.g. the endogenous GH AUC of 0-12 h can yield different results than the AUC of 12-24 h within the same individual due to the secretion of variable pulses. This contributes to high variability in these summary statistics. Furthermore, the commonly used single GH measurement or multiple-point mean [8,9] does not capture the total secretion profile of an individual and thereby limits the correct quantification of a possible drug effect over time.
More advanced analysis methods, such as deconvolution analysis, have been developed to extract more information from a pulsatile profile [10,11]. With deconvolution analysis, the observed concentrations are treated as the product of secretion and elimination processes. The underlying pulsatile secretion processes are estimated as Gaussian shaped events. The time points of these events are optimized to fit the data via multiple in-and exclusion steps [3,10,12]. Deconvolution analysis provides information on the regularity, the frequency, the amplitude, baseline secretion and the secretion width of pulses on an individual level [12,13]. Even though this increases the amount of information that is retrieved from pulsatile GH concentration-time profiles, time cannot be used as a continuous variable in the reported tables. Thus, deconvolution results are still dependent on the study design, and it therefore has limited utility in comparing results between studies where different dosing regimens or study designs are used.
In drug development, population approach non-linear mixed effects (NLME) models are often used to study the pharmacokinetics and pharmacodynamics (PK/PD) of drugs over time. For example, direct effect, turnover or pool models are commonly used to describe endogenous pituitary hormone secretion over time [14][15][16]. However, such pharmacodynamic models cannot account for a highly variable pulsatile secretion. In the literature, the implementation of pulsatile functions in NLME models has been limited to compounds with a low number of pulses (melatonin [17], LH [2], ACTH [18]). When a higher number of pulses is observed, the numerical complexity of the model increases and the model stability decreases when the pulse location, duration and amplitude of multiple pulses need to be estimated.
The aim of this study is to develop a new method to quantify and model pulsatile data for the development of a pharmacodynamic model that is able to quantify highly variable pulsatile secretion patterns over time by combining deconvolution analysis techniques and NLME modeling. Deconvolution analysis was previously performed on data from a clinical study where GH concentrations were sampled every 10 min over a 24 h period in normal weight, upper body obese (UBO) women with large visceral fat areas and lower body obese (LBO) women with small visceral fat areas, before and after weight loss [19]. This resulted in the identification of reduced GH secretion in UBO subjects and no significant difference in GH half-life or volume of distribution. This densely sampled and highly variable dataset was used for the model development of this study. The deconvolution method used in the study of Pijl et al. [19] has been improved and applied in a new software package, AutoDecon, which was implemented in this study [3]. Furthermore, GH concentration-time profiles after the administration of hypothetical drugs that have antagonistic or agonistic properties were simulated. As a proof of concept of the application of this method on the quantification of a drug effect, a simulated drug effect was re-estimated using data from a clinical trial simulation.

Study design
Data were obtained from a clinical study which has been reported previously [19]. In short, 16 women (8 UBO and 8 LBO subjects) followed a weight loss diet to study GH kinetics, before and after weight loss, compared to normal weight control subjects (N = 8). Blood samples for GH analysis were taken at 10 min intervals for 24 h, resulting in a maximum of 144 samples per individual per occasion. The lower limit of quantification (LLOQ) for the GH immunofluorometric assay was 0.03 mU/L.

Deconvolution analysis
Individual deconvolution analysis of the 24 h GH concentration time profiles was performed in AutoDecon [3]. This software provides a ''nonsubjective, standardized, and completely automated algorithm'' [3] for the analysis of pulsatile profiles. During data assembly, the expected measurement times with 10 min intervals were used, as this program requires regularly spaced time input. Separate data files were created for each individual in which the mean GH observation of the duplicate measurement, the standard error of the mean and the number of analysis replicates for each time point were included. The previously reported coefficient of variation of the observations over the concentration range was used to calculate the standard error of the mean for each observation [10,19]. For AutoDecon to function, initial estimates of the secretion width and an initial half-life of GH (5 and 15 min, respectively) were used as input. These parameters were optimized for each individual during the deconvolution analysis. The convolution integral that is implemented in AutoDecon is depicted in Eq. 1 [3].
where C(t) is the GH concentration over time, consisting of the integral of all secretion events (S n (s)), plus the concentration at time point zero (C(0)). The elimination function (E) can follow a 1-or 2-compartment disposition for GH, determined in AutoDecon. Deconvolution analysis resulted in the estimation of an individual's pulse frequency, pulse times, secretion width, baseline secretion, pulse amplitude and half-life. After deconvolution analysis, statistical comparison of the results between the LBO and UBO groups compared to normal weight control subjects was done using an independent 2-group t test.

Model development
For each individual, the pulse frequency and the individual pulse times, retrieved from the deconvolution analysis, were included in the dataset for NLME modeling (NON-MEM V7.3 [20]). The random effects structure was incorporated in the model by a ln-normal transformation of the random effects (g) on the population parameters [21]. Significant inter-individual variability (IIV) on population parameters was included in the model following a forward inclusion method (p \ 0.05). The residual error distribution (e) was drawn using an additive, proportional or combined (additive and proportional) residual error structure using parameters from a normal distribution. Various types of variance-covariance matrices were tested for the correlation between the random effects when identified by Pearson correlation plots. Due to high IIV and intra-individual variability in the height of the pulses (pulse amplitude) within a 24 h concentration-time profile, each pulse was estimated as a different occasion (BOV) which enabled the estimation of GH pulses of different heights within one individual. The amplitude of a pulse was modeled according to Eq. 2.
where Amplitude n is the amplitude of pulse n, h population is the population amplitude parameter, g is the random effects distribution of the IIV and j n is the BOV, the intra-individual variability, for pulse n. Two equations (Eqs. 3, 4) were tested to fit a pulsatile event in the NLME model. Equation 3 is adapted from the documentation of AutoDecon [3], Eq. 4 is adapted from a previous publication that models a single Gaussian shaped pulse in NONMEM [2].
where S n (t) is the secretion over time for pulse n, Amplitude n is the amplitude of pulse n, PulseTime n is the time which corresponds with the maximum secretion time for pulse n (retrieved from deconvolution analysis), and Se-cretionSD corresponds with the width of the pulses. In Eq. 3, the exponential transformation limits the S n (t) to positive values only. In Eq. 4, the exponents 2 and 4 were evaluated during model development.

Covariates
The following covariates were explored: weight, height, age, lean body mass (LBM), total body water (TBW), total body fat, percentage fat mass, percentage LBM and percentage TBW. All weight related covariates were calculated using Bioelectrical Impedance Analysis (Bodystat 1500, Bodystat Ltd., Isle of Man, UK). Covariate relationships were explored using visual exploration of the Pearson correlation plots and their correlation coefficients. When a correlation coefficient was 0.5 or higher, covariate relationships were formally tested for significance in the structural model using linear, power and exponential relationships. Circadian rhythmicity was explored as a covariate on the Amplitude parameter using a cosine function with a 24 h acrophase or as a day/night effect.
Covariates were included using a forward inclusion method (p \ 0.05) combined with backward elimination (p \ 0.01) and centered around their mean values.

Model evaluation
Models were evaluated on basis of objective function value (OFV, which approximates -2*Log Likelihood), goodness of fit (GOF) plots and numerical evaluation [22]. Model hypothesis testing was done under the assumption that the difference in OFV is v 2 -distributed with the degrees of freedom determined by the number of additional parameters in the more complex model. A drop in OFV of more than 3.84 points (p = 0.05) resulted in accepting the model with one additional degree of freedom. For backward elimination of a covariate, an increase of less than 6.6 points in OFV (p = 0.01) was needed for exclusion. Model comparison implementing Eqs. 3 and 4 was done using the Bayesian information criterion (BIC) due to changes in the structural model. GOF plots consisted of population-and individual model predictions versus observations, conditional weighted residuals with interaction (CWRESI) versus clock time and population predictions. Numerical evaluation was based on the relative standard error (RSE) for population parameters, normalized prediction distribution error (NPDE) analysis, coefficient of variation for random effects (CV%) and the condition number [22,23].

Simulation
Multiple simulations with the developed model were performed to visualize the effect of hypothetical drugs on the GH concentration-time profile targeting GH secretion. For these simulations, the parameter estimates of the developed model were used to simulate a typical individual with a pulse interval of 1.57 h (estimated mean pulse interval).
The half-life of the hypothetical drug was fixed to 6 h to simulate a short-term effect which can still be observed in a 24 h period. The drug effect was implemented using an E max relationship driven by the amount of the hypothetical drug where the maximum effect was reached immediately after bolus dose administration (10 mg). The effect of the drug on the pulsatile secretion was modeled using Eqs. 5 and 6.
where E max varied between -1 and 10 to simulate inhibitory and stimulatory effects and was 0 for simulations of the typical individual. A(t) is the amount of the hypothetical drug remaining in the body. For simulations of the inhibitory drug effect, the E max was equal to -1, -0.75 and 0. For simulations of the stimulatory drug effect, the E max was equal to 0, 2, 5 and 10. EA 50 was fixed to 2 mg with c equal to 5 to simulate a relatively fast offset of the drug effect within 18 h after dosing. Covariate relationships were simulated at their mean values. A new clinical trial was simulated to investigate whether the applied drug effects could be correctly re-estimated using the developed method proposed in this study. Simulations were performed including the identified IIV and residual error structure of the developed model. Pulses were simulated at regular time intervals. Additionally, a CV% of 10% was added on the EA 50 . A total of 5 cohorts of 8 subjects were simulated (placebo, 2.5, 5, 7.5 and 10 mg) with an E max of -0.9 (inhibiting the GH secretion by 90%) and an EA 50 of 3 mg with c equal to 5. The parameters estimated in the deconvolution-analysis-informed population model were fixed to their corresponding values. The re-estimated drug effect parameters were then compared to the simulated 'true' model values to judge the ability to correctly recover the drug effect given this true model.

Results
A total of 2 subjects (1 from the LBO and 1 from the UBO group) did not complete the occasion after weight loss. For these individuals, only data from the visit before weight loss were included in model development. Data below the LLOQ (n = 52, \1%) were excluded from this analysis. A total of 5377 GH observations were used for model development. Figure 1 visualizes the GH concentrationtime profiles during the 24 h observation period for three individuals. High intra-and inter-individual variability can be observed in the pulse amplitude and the time of GH secretion. Table 1 shows the summary of the deconvolution analysis results for each group, before and after weight loss. A 1-compartment elimination function was identified as the best fit for GH disposition. The secretion width of the pulses was found to be significantly lower in LBO subjects before weight loss and in all UBO subjects compared to normal weight subjects. The UBO subjects before weight loss showed a reduction in the baseline secretion compared to normal weight subjects. No significant differences in the pulse frequency, half-life, amplitude or pulse interval were identified between normal weight and obese subjects.

Model development
The GH observations were best described using a turnover compartment, as depicted in Fig. 2. The baseline secretion was modeled as a steady-state condition, using a continuous zero-order input (k in ) in the central compartment and first-order elimination (k out ) that describes the elimination of GH from the body. The k in was estimated as Baseline Á k out so that the Baseline parameter is estimated as mU/L. The pulsatile secretion is the sum of the secretion of all pulses at a certain time point (Eq. 7) where npulses is the pulse frequency of an individual. The differential equation of the central GH compartment is presented in Eq. 8.
A large proportion of the subjects showed GH concentrations above baseline at the start of the measurement period (e.g. Figure 1, black line) due to the secretion of an endogenous pulse of GH prior to the start of the observation period. To account for these initial concentrations, the central compartment was initialized at an initial concentration (A_0). The estimation of BOV between the GH profiles before and after weight loss within one individual resulted in numerical difficulties due to the large number of random effects (54?) on the amplitude parameter to be estimated. Therefore, the two occasions of one individual were stratified using unique subject identifiers. Thereby reducing the number of random effects within one individual but losing the ability to identify intra-individual variability between occasions before and after weight loss.

Covariate analysis
The estimation of a 24 h cosine function or a day/night effect on Amplitude did not result in the improvement of the OFV, indicating that no circadian rhythmicity could be identified on this data. Inspection of the correlation plots of the x 2 distribution identified covariate relationships between the distribution of Baseline, Amplitude and Se-cretionSD with the TBW (%), as shown in Fig. 3. The TBW (%) also showed a high degree of correlation between the weight of a subject in which heavier subjects had a lower percentage total body water. A power covariate relationship showed to be superior over other tested relationships. The inclusion of this covariate relationship on the Baseline, Amplitude and SecretionSD parameters resulted in a significant decrease (DOFV = -48) in OFV. Backward elimination did not result in the removal of a covariate in this model. The previously observed covariate correlations between the random effects were reduced to a random scatter around 0 after inclusion of the covariate, indicating that part of the identified variability could be explained by the included covariate. The correlation plots of TBW(%) with the x 2 distribution of Baseline, Amplitude and SecretionSD after inclusion of the covariate relationship in the structural model, can be found in online resource I.

Model evaluation
For visualization purposes, individual model fits over time, for one individual per group, are depicted in Fig. 4.   Table 2. The population parameters show low RSEs (ranging from 2 to 5%), indicating high accuracy in the estimation of the population parameters. The CV% was moderate ranging from 26.9 to 70.8% for most population parameters. High CV% was identified for the A_0 and the j on Amplitude with a CV% of 521 and 302% respectively. The high CV% for these parameters originates from the differences in initial concentration (the possible occurrence of a pulse before the start of the observation period) and the high variability in pulse amplitudes between the separate pulses within one individual. The condition number of the final model was 11.3, indicating no ill conditioning [23]. NPDE analysis results can be found in online resource II. The NONMEM code of the final deconvolution-analysis-informed population model (FOCE?I, ADVAN = 13, TOL = 6, SIGL = 6, NSIG = 3) can be found in online resource III.

Simulation
The simulated GH concentration-time profiles are depicted in Fig. 6 and stratified on stimulatory (Fig. 6a) and inhibitory (Fig. 6b)     the GH secretion returned back to a pulsatile profile. The inhibition of GH secretion by 75% showed reduced GH concentrations but did not completely counter the pulsatile secretion of GH in this scenario. The drug effect over time is depicted in online resource IV. The simulated clinical trial resulted in a total of 5800 observations which were then fitted to the true simulated model structure (including the drug effect) and to a reduced model structure (excluding the drug effect). Modeling the GH secretion while excluding a drug effect resulted in an OFV of -1871 points. The inclusion of an E max relationship for the drug effect resulted in a significant drop in OFV (DOFV = -234). The re-estimated model parameters are presented in Table 3 and are near identical to the parameter values used for simulation of the clinical trial data. The low RSE indicates high accuracy in the re-estimation of the drug effect parameters given the true simulated model. High shrinkage (51%) was observed on the random effects of the EA 50 . The NONMEM model codes used for simulation and re-estimation, including the GOF plots for the re-estimated model, are provided in online resource V.

Discussion
The developed method, using NLME modeling with NONMEM, has proven to adequately describe variable endogenous pulsatile GH concentration-time profiles in women. This resulted in the identification of three  covariate relationships and the quantification of BOV and IIV on the amplitude of pulses. In addition, simulations have been performed to show hypothetical drug effects on GH concentration-time profiles. The correct re-estimation of the drug effect in the simulated clinical trial data shows the applicability of this method in the analysis of variable pulsatile data. The RSEs of the population parameters were low (\10%), indicating stable predictions of population parameters. The CV% on A_0 and Amplitude were high, which is in agreement with the observed high variability in the observed concentration at time point 0 and between the amplitude of pulses within an individual, respectively. The high GH concentrations at the first measurement indicate that sampling should be performed at several pre-dose time points to be able to quantify GH pulses that may have occurred before dose administration. The population halflife of GH was estimated as 15 min (ln(2)/k out ) in this population and can be compared with the half-life reported in literature [27][28][29]. No circadian rhythmicity was identified in this population. This could be due to potential differences between men and women in regards to nocturnal GH secretion [30].
The study of Pijl et al. [19] resulted in the identification of a stratification between the GH secretion in UBO compared to LBO and normal weight subjects. The deconvolution analysis performed in this study did also result in the identification of a lower baseline secretion in UBO subjects before weight loss. No changes in the pulse amplitude were identified between groups. In this study, significant continuous covariate relationships have been identified in which the total body water was the best covariate explaining the IIV for Baseline, Amplitude and SecretionSD. Furthermore, the possibility to follow the GH concentration over time in a NLME model increased the information that is retrieved from these dense sampled observations.
During model development, it was required to stratify the two occasions, before and after weight loss, of one individual with unique identifiers to prevent numerical instability of the model. This resulted in the loss of estimation of the intra-individual variability between these separate occasions. Due to the high variability in the GH profiles that were observed between the occasions before and after weight loss, the long duration between the occasions (up to 6 months) and the nature of this paper to establish a new quantification method using NLME modeling, this stratification was deemed appropriate.
The Gaussian shaped events that are fitted in deconvolution analysis have the advantage to be applicable on a wide range of pulsatile compounds [31,32]. However, this generalization inherently brings the disadvantage that no mechanistic information of the biological regulation of the compound is incorporated in the structural model. Especially when multiple components are involved in the regulation of the secretion, e.g. stimulation by growth hormone releasing hormone and inhibition by somatostatin in the case of GH, more mechanistic insights can be beneficial in explaining the variability between individuals and possible extrapolation to disease states.
Compared to previously used non-compartmental analysis, the current analysis enables the quantification of a drug effect targeting GH secretion over time in which an E max or other PK/PD relationships can be determined using NLME in NONMEM. The correct estimation of the PK/PD relationship parameters on a pulsatile compound will be mainly dependent on 3 components: (1) the frequency of sampling, which is critical for compounds with a short halflife [11], (2) the duration of the observation period should preferably include the on-and offset of the drug effect and (3) the expected variability in response should be taken into account when choosing the study sample size.
The method proposed in this study can be applied to analyze or simulate different pulsatile profiles. From the analysis step, information can be obtained on the pulsatile behavior (pulse frequency, width and amplitude) in a population and on the PK/PD relationship. From a clinical perspective, this provides more insight in the relationship between the pulsatile secretion of a compound and the effect that a drug has over time compared to differences in the AUC or C max when comparing placebo with treated individuals. The provided simulations in this study visualize what effect hypothetical drugs can have on a pulsatile profile. This was simulated for drugs having agonistic or antagonistic properties. After administration of an agonistic compound with an E max of 10, the maximum GH concentration reached was *50 mU/L. These simulations are physiologically plausible since the simulated concentration are in agreement with data from acromegalic patients in which GH concentrations are in the same range [33]. The provided antagonistic simulations suggest a reduction to baseline after fully blocking the pulsatile secretion (E max = -1). These antagonistic drug simulations (Fig. 6b) correspond with the clinical response as shown after treatment with Ocreotide, a somatostatin analogue, which can completely block the endogenous pulsatile secretion of GH for 5 h after dosing [7]. Blocking the pulse amplitude by 75% lowers the pulsatile GH secretion but did not fully block the pulsatile profile. The pulse interval can be changed in the model code to simulate a different pulse frequency, including variability on the pulse interval. Simulations incorporating the established E max and/or EC 50 values of a known drug on GH secretion can then be used to support dose selection or to investigate novel dosing regimens for the studied drug. The observation period of a study can be expanded if simulations show that the full effect is reached at a later stage than expected (e.g. indirect response or delayed release formulation of the drug) or that a fast offset of the drug effect suggests the need for more selective monitoring for multiple hours at an early stage compared to a full 24 h period.
Instead of using 'simple' summary statistics (mean, AUC, C max ) for complex and variable pulsatile profiles, it is now possible to quantify and model a pulsatile profile over time and to identify what effect a potential drug has on the secretion. This two-step deconvolution-analysis-informed population pharmacodynamic model enables the possibility to analyze highly variable pulsatile profiles with a NLME PK/PD model in NONMEM.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Ethical approval All procedures performed in this study involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. For the original clinical study, approval by the ethics committee of the Leiden University Medical Center was obtained. For this retrospective study, formal consent was not required.
Informed consent Informed consent was obtained from all individual participants included in the study.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.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.