A Bayesian K‐PD model for synergy: A case study

Coadministration of 2 or more compounds can alter both the pharmacokinetics and pharmacodynamics of individual compounds. While experiments on pharmacodynamic drug‐drug interactions are usually performed in an in vitro setting, this experiment focuses on an in vivo setting. The change over time of a safety biomarker is modeled using an indirect response model, in which the virtual pharmacokinetic profile of one compound drives the effect of the other. Several experiments at different dose level combinations were performed sequentially. While a traditional frequentist analysis consists of estimating the model parameters based on all the data simultaneously, in this work, we consider a Bayesian inference framework allowing to incorporate the results from a historical dose‐response experiment.


INTRODUCTION
For several years, the coadministration of multiple compounds has become a standard practice in the treatment of patients. Nevertheless, the coadministration of 2 or more compounds can potentially-intentionally or unintentionally-alter the systemic exposure and/or the pharmacological effects of the individual compounds. The effect on the systemic exposure is commonly called the pharmacokinetic (or PK) interaction, whereas the pharmacological effect of the drug on the body is called the pharmacodynamic (or PD) interaction. Experiments on PD drug-drug interactions are usually performed in an in vitro setting, but are rarely undertaken in an in vivo framework, unless at late stage drug development (Phase 3). This manuscript will focus on an exceptional preclinical in vivo example, where in vivo rat data were used to investigate the synergy of 2 compounds (a marketed compound and a novel experimental compound) at a PD level.
The first PD models were introduced by Levy, 1,2 who showed that the intensity of the effect attributed to a compound can be linearly related to the logarithm of the concentration of the compound itself. In the late 1960s, recognizing that a linear relationship was only valid within a limited concentration range, Wagner 3 proposed a nonlinear concentration-response relationship, also known as the Emax model. The sigmoid Emax model is an extension of the Emax model that allows for different slopes of the concentration-effect relationship. The above models assume that effects are proportionally mediated by drug concentrations. However, sometimes a time lag between concentration and observed effect may occur. Such a delay can be caused by different factors (eg, kinetics of receptor activation, active metabolites, or change in baseline over time). More complex models were proposed to explain this phenomenon, also called hysteresis: biophase distribution models (Sheiner et al 4 ), indirect response models (Dayneka et al 5 ), and signal transduction models (Sun and Jusko 6 ).
The goal of this paper is to integrate knowledge from sequential in vivo experiments using a Bayesian hierarchical semimechanistic model. An indirect response model with a virtual PK profile (also called K-PD model 7 ) will be used, with a PD interaction at the IC 50 (the drug amount that produces 50% of maximum inhibition). Random effects are incorporated to allow for differences between animals.
The structure of the paper is as follows. Section 2 describes the motivational case study. In Section 3, a K-PD model for drug-drug interaction is proposed, and the synergy between compounds is discussed in terms of the model. Details on the Bayesian specification of the model are given as well. Section 4 then describes the results of this method applied to the data of Section 2 and a final discussion is given in Section 5.

CASE STUDY
The case study consists of 11 experiments as well as a historical experiment to be integrated in the analysis. It was part of the preclinical safety evaluation of a new experimental compound meant to be coadministered with a compound already available on the market. 8 For the marketed compound, historical experimental data were used to assess its safety. 9 Ten doses, ranging between 0 and 10 mg/kg, were investigated. A total of 55 rats were randomly allocated to receive a single dose of the marketed compound (Table 1). For each of those animals' body temperature, the safety biomarker of interest was assessed up to 24 hours after oral administration. Dose related changes were observed for extremely high doses, indicating potential side effects ( Figure 1).  Figure 1 shows the observed body temperature profiles for 9 representative individuals. They show, for high dose levels, a rapid decline of body temperature, followed by a return to the baseline level within a few hours. Further, it is noted that both the maximal observed effect and the time to reach the maximal effect increase with increasing dose levels.
In light of the dose response of the marketed compound and with the intent to develop a new compound for coadministration, a new experiment was designed to evaluate whether an interaction between both compounds occurs. A total of 20 animals were used in the experiment, with 5 animals randomized to each of 4 treatment groups. Group 1 received a single dose of the vehicle (0 mg/kg); group 2 received the maximum dose from the historical experiment (ie, 10 mg/kg); group 3 obtained a high dose of the novel compound (40 mg/kg); and group 4 received a combination of both compounds (10 mg/kg of the marketed compound and 40 mg/kg of the new compound). Body temperature was assessed every hour after administration for 4 hours. High doses were chosen to maximize the chance of detecting a PD interaction between both compounds. A PK interaction was not formally statistically tested. However, no major PK interaction was observed in an exploratory experiment. As expected, animals in the vehicle (group 1) and novel (group 3) treatment groups did not show any effect on body temperature, whereas animals attributed to the marketed compound (group 2) showed an effect in line with the historical dose response experiment of the marketed compound. Nevertheless, it was unexpected that animals assigned to the combination group (group 4) showed a more pronounced change in body temperature compared with the group treated with the marketed compound only (Figure 2). To further understand this synergetic behavior and its potential impact on future drug development, the experiment was repeated 10 times, each time with a different combination of dose levels (Table 2, Figure 3). Because plasma samples could potentially impact body temperatures, plasma concentration-time profiles were not obtained.

METHODS
In this section, we first describe a K-PD model for the analysis of the historical experimental data of the marketed compound (Section 3.1). This model is then extended to describe the synergy of compounds in Section 3.2. Details on the Bayesian estimation are given in Section 3.3.

Analysis of historical data using a K-PD model
Let us first focus on the analysis of the historical experiment of the marketed compound. Turnover models, also referred to as indirect PD response models, cope with profiles such as the ones observed in Figure 1. These models take into account the mechanism of action of the drug and other aspects of biological reality (Dayneka 5 and Sharma and Jusko 10 ) and are therefore called semimechanistic models. An enzyme, protein, or biomarker is produced and proportionally eliminated continuously over time. In the absence of the compound, the production and elimination processes are in equilibrium. The administration of a compound results in an increase of the drug concentration in plasma over time. This drug will bind to a receptor in a concentration-dependent manner, which in turn can block or boost either of these processes pending the physiological pathway. Turnover models are appropriate when there is a time lag between exposure and the observed effect of the biomarker. For an extensive review of turnover models, we refer to Sharma and Jusko 11 and Gabrielson and Weiner. 12 Let R it denote the observed body temperature for animal i at time t. We propose the use of a turnover model with an inhibition process operating on the production of body heat, namely, Equation 2 describes the body temperature change over time.R it is the expected body temperature for animal i at time t. The parameters k in and k out are the rate for production and elimination of body heat, respectively. Inhibition of the production depends on the amount of the marketed compound in the central compartment C it and is characterized by the maximal inhibition of the production rate, I max (0 ≤ I max ≤ 1) and the drug amount to inhibit the production to 50% (IC 50 ). In what follows, we will refer to this IC 50 value as the potency of the compound. The model takes into account residual error via Equation 1, with 2 R the measurement error variance. As the amount C it for animal i at time t could not be measured, a virtual PK profile is assumed to drive the inhibition of body heat 7 : This represents a first-order one-compartment PK model with oral absorption. The amount of the marketed compound in the absorption depot is represented as A e,it . The parameter k a is the absorption constant, whereas the parameter k e is the elimination constant.
The full K-PD model is given by Equations 1 to 4. At time t = 0, A e,i0 corresponds to the dose of the marketed compound, C i0 = 0 as no drug amount is present in the central compartment before the actual administration of the compound and R i0 = k in ∕k out corresponds to the fact that body temperature is in a steady state condition prior to administration of the compound.
It should be noted that the lack of information on the exposure-time profile causes a form of model overspecification, as it does not allow the identifiability of some parameters. In particular, IC 50 could not be estimated and was set to 1, as discussed in Section 3.2. In Jacqmin et al, 7 an intravenous (IV) infusion model was used, so that the underlying exposure is fully corresponding to the dose level. In our context, an oral absorption model is implemented, where the bioavailability cannot be estimated as such. As a result, IC 50 becomes not estimable (confounded up to a factor). The absence of information on the absorption phase of the marketed compound further causes nonidentifiability of k a : The compound first needs to reach the site of action to induce a physiological downstream reaction. Hence, the unobserved absorption needs to be faster than the change in response. This implies that k out is a lower boundary of the absorption rate, so the true absorption rate is censored. To overcome this identifiability problem, k a was set to e 5 , which reflects a fast change in body temperature during the first hour. This assumption of an extremely fast absorption process makes the model comparable to the model with IV bolus administration illustrated by Jacqmin et al. 7 A sensitivity analysis confirms that the IV bolus and our proposed model give the same results.
To allow for heterogeneity amongst animals, a random effect was added for body temperature at baseline R i0 , Inclusion of random effects on other parameters (k in and k out ) was also evaluated.

Synergy
Up to this point, only the effect of the marketed compound has been accounted for. In this section, the model is extended to include the effect of the novel compound onto the marketed compound. Given the absence of a PK interaction between the 2 compounds and the inactivity of the novel compound on body temperature and given the fact that the coadministration did not affect the maximal effect (I max ) observed in the historical experiment on the marketed compound, it is assumed that the presence of the novel compound increases the potency of the marketed compound with respect to body temperature. Therefore, the following definition for the potency is proposed: where A n,i0 and A e,i0 are the doses of the novel and marketed compound for the ith individual.
Although, as mentioned above, IC 50 = exp ( 0 ) is fixed to 1 due to the lack of PK information, this parameterization allows to generalize the model to a nonvirtual PK profile, from which IC 50 can be estimated.
The main effect represents how the novel compound influences the binding of the marketed compound at the receptor. Note that in the particular scenario of A e,i0 = 0, A n,i0 > 0 and ≠ 0, one obtains IC 50new = IC 50 exp ( A n,i0 ) . Hence, in this case, IC 50new ≠ IC 50 , which makes sense since binding cannot be observed when no dose of the marketed compound is administered.
The interaction coefficient can be initially considered as controversial, as in this way the potency of a compound is made dependent on its dose. However, also the amount of marketed compound near the receptor can modify the impact of the novel compound on the potency of the marketed compound itself. This interaction is expressed by .
An extensive publication list on models for coadministration of compounds can be found in Harbron et al. 13 However, such models measure in vitro synergy. One of the most widely used methods is the Loewe additivity model 13 : where y represents the effect that 2 compounds produce when they are administered as monotherapy with doses D y,1 and D y,2 , respectively and d 1 and d 2 represent the doses of the 2 compounds that produce the same effect as D y,1 and D y,2 when both compounds are coadministered. In the specific case where one compound is inactive if administered as monotherapy, the in vitro Equation 6 simplifies to Thus, an in vitro synergistic behavior is observed if the dose of the active compound required to achieve the effect in combination is less than required in monotherapy. 14 Note that, according to the so-called Saariselkä agreement, 15 the terms synergism/antagonism (without a leading adjective) are used for an increased/decreased effect following the coadministration when only one compound is active.
The proposed parameterization given in (5) is consistent with the Loewe additivity model in the special situation of one active compound. The synergy is built into the model as follows: Based on this rearrangement of the equation and under the assumption of linear kinetics, one obtains which reflects how much one needs to adjust the dose in the combination therapy to obtain the same effect y as under monotherapy. Furthermore, one can easily see that Hence, in case A n,i0 = 0, it follows that IC 50 new = IC 50 . This is comparable with the Loewe additivity model where d 2 = 0 implies that d 1 = D y,1 or, in other words, that the compound is additive to itself.
So, when 2 compounds are coadministered, the evaluation whether the interaction is additive, synergistic, or antagonistic depends on the value of IC 50 new : • IC 50 new = IC 50 , means that the effect of the marketed compound remains independent of the novel compound and corresponds to the definition of additivity; • IC 50 new < IC 50 means that a lower exposure for the marketed compound is required in the presence of the novel compound to obtain the same effect, hence synergy; • IC 50 new > IC 50 means that a higher exposure for the marketed compound is required in the presence of the novel compound to obtain the same effect, hence antagonism.
This illustrates that the model presented in this paper can be considered as an extension of the Loewe additivity approach to an in vivo framework in the special situation when a compound is inactive if administered as a monotherapy.

Bayesian data integration
The data of all experiments were collected at the same laboratory by the same scientist. This justifies to pool the experiments and analyze all data altogether.
Initially, a likelihood-based model was used to estimate the parameters based on all pooled experiments (1-11) using NONMEM software 16 (version 7.3.0). This however resulted in convergence issues even with different starting values and different estimation methods (SAEM and Laplace). Therefore, a Bayesian model was considered as it allows to take into account prior knowledge from the historical dose-response experiment on the marketed compound, thus avoiding potential performance issues and facilitating convergence by restricting the exploration of the parameter space. Prior distributions were chosen reflecting the knowledge on the marketed compound. In particular, expected values were set equal to point estimates obtained from the historical data. The standard errors were doubled to down weight the information derived from the previous analysis (since the historical experiment was conducted 2 decades before the new experiments. 9 ) The priors for k e and k out were assumed lognormally distributed, I max was assumed to follow a beta distribution, and a uniform distribution with parameters −10 and 10 was used as a prior for and . An inverse gamma distribution was used as a prior for the measurement error variance 2 R , and lognormal and inverse gamma distributions were used as priors for the mean and variance of the random effect (R 0 and 2 R 0 ). Table 3 summarizes the set of prior distributions for the parameters of the above presented model.
The use of uninformative priors was also evaluated, but it led to convergence issues because of the high model complexity relatively to the amount of data.
The Bayesian model was estimated using the No-U-Turn Sampler algorithm in Stan 17 (rstan version 2.14.1). Specifically, 3 Markov chains were run, with a total length equal to 3000 (2000 draws for each chain, of which 1000 warmup). The correlation among Markov chain samples and the degree of convergence were also assessed through the effective sample size and the R-hat statistic, respectively.
The graphic representation was made using the statistical software package Rversion 3.3.3 (Foundation for Statistical Computing, Vienna, Austria).

RESULTS
In this section, we present the results of the above described synergy model. We first show the evaluation of model performance, followed by model parameter estimation and model prediction.
Our Bayesian model seemed to fit the data well. Individual estimates of body temperature for data of experiment 1 are shown in Figure 4. Columns 1 and 2 show the observed and estimated time profiles of rats in treatment group 1 (no treatment) and group 3 (novel compound), respectively. Clearly, when the novel compound is administered as monotherapy, there are no temperature variations over time. Column 3 displays the profiles in treatment group 2 (marketed compound), showing a substantial decrease in body temperature, which reaches its maximum 1 hour after oral administration, followed by an increase that reaches the initial steady state 4 hours after administration. Column 4 shows the profiles for rats in treatment group 4 (coadministration): It is possible to observe a later maximal effect and an increase in potency as the decrease is much more pronounced, reaching its maximum approximately 2 or 3 hours after administration. Body   temperature takes longer than 4 hours before reaching the steady state. The estimated time profiles in experiments 2 to 11 are given in the Appendix S1. Table 4 shows posterior mean, 95% credible intervals, effective sample size, and R-hat statistic for the parameters in the model. The value of R-hat is equal to unity for all parameters, meaning that the chains are well mixed and they are all exploring the same regions of parameter space. The effective sample size is equal to its maximum for all parameters except for I max , and 2 R 0 . Particularly, the lower value for indicates a higher correlation among the Markov chain samples for this parameter. This was to be expected since only the observations belonging to the combination treatment group actually contribute to the estimation of .
The maximal inhibition of the production rate is likely between 16% and 22%. The negative posterior means of and (with credible intervals not including 0) confirm the presence of a pharmacodynamic synergy between the study compounds, as IC 50 new < IC 50 .
The model estimates are used to illustrate IC 50 new ( Figure 5) and maximum predicted temperature change ( Figure 6) as a function of both the marketed and the novel compound. It can be observed that the synergistic behavior affects body temperature only for high doses of the marketed compound (higher than 2.5 mg/kg), whereas low doses produce negligible deviations in body temperature independently of the novel compound dose. This illustrates that the synergistic behavior between the 2 compounds is only observed at dose level combinations well above the clinical range (until 0.1 mg/kg) and that clinically relevant dose levels remain without effect.

DISCUSSION
In this paper, a novel K-PD model for the estimation of synergy is presented. The use of a physiologically inspired model enables translating the scientific questions of interest directly into tangible parameters, leading to clearly interpretable conclusions. Our proposed model directly quantifies the impact of the new compound on the safety of the existing compound in 1 or 2 parameters. It is therefore easier to interpret when compared with a model such as a spline or polynomial. It has been demonstrated that the model is an extension of the in vitro synergy methodology to in vivo data using the Loewe additivity model in the situation when a compound is inactive if administered as a monotherapy.
Unlike the frequentist framework, the Bayesian framework allows to borrow information from the historical experiment, thus facilitating convergence. The model has proven to work well when all data coming from different experiments conducted in different time periods were pooled. Further work is being devoted to keeping the sequential nature of the experiments, ie, fitting a Bayesian model in a sequential manner, so that the posterior distributions resulting from an experiment are used to determine the hyperparameters of the prior distributions of the experiment which follows.
As mentioned previously in Section 3, several models were performed by including random effects on different parameters. The choice of random effects highly affected the model performance: As an example, substantially biased parameter estimates were observed when the random effect was allocated to the parameter k out . A possible explanation of such behavior could be related to the highly correlated parameter space, which leads to overcompensated estimates. A highly correlated parametric space is a typical property of nonlinear hierarchical models. Nevertheless, the parameter correlation becomes particularly high in models with virtual variables. 18 Prior elicitation represented another challenge in fitting the presented model: The use of uninformative priors led to convergence issues due to the combination of a complex, but physiologically inspired model and a highly correlated parameter space. When data are scarce, there is a lack of information that leads to a higher correlated parametric space, where parameters compensate each other. Increasing the informativeness of priors, on the other hand, gives the possibility to add prior information, causing more stability and less convergence problems. In such context, the information coming from the historical experiment was particularly useful, because it helped choosing the prior distributions for analyzing the data under study. These evaluations will be more deeply explored in future work. The stiffness of prior choice did not allow the consideration of multivariate prior distributions, which incorporate the intercorrelations among parameters, as this would further increase the computational burden and convergence issues.
It should be noted that only a limited design space was assessed in this case study. This is due to its explorative nature: In a preclinical context, practical constraints such as limitation of animal facilities and workload of laboratory technicians prohibit the design of an optimal set of dose combinations. Furthermore, the studies were designed without the consultation with a statistician. The case study illustrated in this manuscript highlights the importance of involving statisticians at the experimental design phase, so that both scientists and statisticians can benefit from such a collaboration. Nevertheless, the compound doses assessed in this study were sufficient for the calculation of accurate parameter estimates and good model fitting.
Lastly, it needs to be emphasized that expressing the potency in function of treatment doses is solely used in the absence of the plasma concentration-time profiles. The latter would be a more suitable covariate to describe the activity at the receptor and therefore the synergy. Hence, the dose levels are only an approximation of the true exposure at the receptor levels. As previously described, the lack of PK information causes nonidentifiability of IC 50 . It should be noted, however, that the scientific question of interest focuses on the change of IC 50 in the presence of the novel compound, rather than quantifying the value itself. This proposed method enables to answer this question through a restricted number of parameters.