A Model-Based Pharmacokinetic/Pharmacodynamic Analysis of the Combination of Amoxicillin and Monophosphoryl Lipid A Against S. pneumoniae in Mice

Combining amoxicillin with the immunostimulatory toll-like receptor 4 agonist monophosphoryl lipid A (MPLA) represents an innovative approach for enhancing antibacterial treatment success. Exploiting pharmacokinetic and pharmacodynamic data from an infection model of Streptococcus pneumoniae infected mice, we aimed to evaluate the preclinical exposure-response relationship of amoxicillin with MPLA coadministration and establish a link to survival. Antibiotic serum concentrations, bacterial numbers in lung and spleen and survival data of mice being untreated or treated with amoxicillin (four dose levels), MPLA, or their combination were analyzed by nonlinear mixed-effects modelling and time-to-event analysis using NONMEM® to characterize these treatment regimens. On top of a pharmacokinetic interaction, regarding the pharmacodynamic effects the combined treatment was superior to both monotherapies: The amoxicillin efficacy at highest dose was increased by a bacterial reduction of 1.74 log10 CFU/lung after 36 h and survival was increased 1.35-fold to 90.3% after 14 days both compared to amoxicillin alone. The developed pharmacometric pharmacokinetic/pharmacodynamic disease-treatment-survival models provided quantitative insights into a novel treatment option against pneumonia revealing a pharmacokinetic interaction and enhanced activity of amoxicillin and the immune system stimulator MPLA in combination. Further development of this drug combination flanked with pharmacometrics towards the clinical setting seems promising.


Pharmacometric Analysis
All modelling and simulation tasks were performed using NONMEM ® 7.4.1 (ICON Clinical Research LLC, Gaithersburg, MD, USA) which was executed by PsN (version 4.7.0 [6]). As compiler, GFortran (version 6.3.0, GNU Compiler Collection) for macOS High Sierra (Apple Inc., Cupertino, CA, USA) was used. Additional software tools to record, visualize, and evaluate the modelling results were Pirana (version 2.9.7 [7]) including XQuartz (version 2.7.11), RStudio (version 1.2.1335), and R (version 3.6.0) with packages xpose4 (version 4.6.1), ggplot2 (version 3.2.1), and vpc (version 1.1.0). For parameter estimation, Laplacian estimation with interaction and due to stability reasons first-order estimation were used for the pharmacometric PK/PD and TTE model, respectively. To handle samples with concentrations below the LLOQ, the so-called M3 method, which estimates the likelihood that the concentration is below the LLOQ, was used [8].
For comparison between nested or non-nested models, analysis of the objective function value (OFV) or Akaike Information Criterion (AIC) [9], respectively, was performed. A decrease of more than 3.84 (p ≤ 0.05 at α = 0.05, df = 1) of the OFV was defined as statistically significant. In addition, model evaluation was performed by considering goodnessof-fit (GOF) plots, simulation-based diagnostics (e.g., VPC), log-likelihood profiling (LLP) and bootstrap analysis [10]. For the VPC, 1000 simulations including interindividual variability (IIV (PK data)) or unexplained variability (PK/PD data) based on the original dataset were performed resulting in a 90% prediction interval and a 90% confidence interval (CI) around the predicted percentiles (5tt, 50th, and 95th percentile). For final models, a nonparametric bootstrap analysis with 1000 bootstrap samples with replacement stratified to the respective study groups was performed to determine the precision and accuracy (95% CI) of the bootstrap parameter estimates [11]. These were compared to the final model-predicted parameter estimates of the original dataset.

Pharmacokinetic Submodel
For AMX concentrations in serum, one-, two-, and three-compartment models after oral administration with linear and nonlinear absorption and elimination processes and with or without lag time were investigated. Potential covariates such as mouse type (RjOrl:Swiss, Balb/cJRj), MPLA coadministration (yes/no) and the dose of AMX (0.4 mg/kg, 14 mg/kg) were investigated: Their influence on the model parameters clearance, volume of distribution, intercompartmental clearance or absorption rate constant as categorial (mouse type, MPLA coadministration) or continuous (AMX dose) covariate was evaluated. The influence of MPLA coadministration on certain model parameters was evaluated by allowing MPLA as binary covariate to change respective PK model parameters proportionally with the AMX dose: where θ1 is the typical model predicted value of a PK parameter P in mice that were not treated with MPLA and θ2 is a value for the influence altering θ1 in case of MPLA coadministration and also depending on the dose of AMX ("DOSE"). Sampling twice per individual mouse due to the developed highly sensitive and little sample volume requiring bioanalytical method, enabled estimation of IIV for the PK data. This source of variability was included for various parameters in an exponential stochastical submodel, assuming PK parameters to be lognormally distributed. The RUV characterising, e.g., experimental imprecision, e.g., in the bioanalytical method and blood sampling was investigated as additive, proportional and combined RUV model.

Bacterial Disease Submodel
To describe the growth of bacteria in untreated mice, a bacterial disease submodel was developed. Based on the visual time course of the bacterial data, the number of bacteria in the lung (Nbacteria, lung) was characterized by a time-delayed (klag) first-order growth rate constant (kg) describing the natural growth of the bacterial strain [12]. To also capture treatment-unrelated killing and natural death processes of bacteria, an additional firstorder rate constant (kkill,lung) term was introduced [13]: Additional implementations such as limiting bacterial growth to a maximum value or different models, e.g., bacterial killing as Emax model [12,13], introducing prebacterial compartments (i.e., bacteria which do not replicate yet) [13] or a simplification without klag, were investigated and evaluated for improvement in OFV or AIC.
In the spleen, implementation of an analogous growth model as in lung was investigated in detail by applying growth, kill, and death kinetics. The transit process of bacteria from lung to spleen was investigated by transit models being able to delay the transition from one to another comportment more physiologically than simple lag times. The firstorder transit rate constant (ktr) was computed by dividing the actual number of transit compartments (n) by the estimated mean transition time (MTT): To more specifically estimate the number of transit compartments, an advanced approach was utilized as proposed by Savic et al. [14] using the Stirling approximation.

Effect Compartment Submodel
To link the PK submodel to the time-delayed PD effect of AMX at target sites in lung and spleen, effect compartments were introduced [13,15,16]. Thus, by allowing the transfer of AMX from serum to a lung or spleen compartment, a concentration-dependent time delay was introduced between the respective compartments and different in-and outflow kinetics (first-order, Michaelis Menten) were investigated. By, e.g., using first-order inand outflow, the concentration in the effect compartment Ce was determined as follows: where ke0 is a first-order rate constant for the antibiotic transfer and Cserum and Ce the concentrations of AMX in serum and the effect compartment, respectively.

Disease and Treatment Submodel
The effect of AMX on certain bacterial disease parameters (kg, klag, kkill,lung) or as bacterial killing process was characterized by investigating different effect models in the lung and spleen compartment such as a simple linear, a power, an ordinary Emax or a sigmoidal Emax model: In comparison, the sigmoidal Emax model served as reference and the alternative models were evaluated for worsening in the model performance: where E(Ce) is the drug effect of a defined effect compartment antibiotic concentration, Emax is the maximum effect, EC50 is the AMX concentration to achieve half-maximum effect, Ce is the apparent AMX concentration in either lung or spleen and H is the Hill factor [13]. Since no PK information for MPLA was collected and only one dose of MPLA was investigated, the influence of MPLA treatment was analyzed as binary covariate (yes/no) and studied for its ability to increase the killing effect characterized by, e.g., kkill, lung (Equation (2)), mathematically implemented as proportionality factor.
To characterize if MPLA had a significant effect on the efficacy of AMX, the potential of MPLA to modify the effect of AMX was investigated, e.g., by implementing an influence on the maximum effect E max or the concentration of AMX to achieve half of the maximum effect EC 50 according to the General Pharmacodynamic interaction model proposed by Wicha et al. [17] as fractional change of respective effect parameters. Determined 95% CI of the estimated change parameter were evaluated, where the inclusion of 1 in the 95% CI indicated no significant change.
RUV was estimated and implemented additively on a logarithmic scale. For the pooled studies, implementation of IIV was impossible due to the fact that only one sample was collected for quantification of bacteria per organ from each individual mouse. This setting did not allow to reliably estimate more than one hierarchical level of variability simultaneously.

Time-to-Event Analysis
In order to investigate survival, a TTE analysis was performed based on reported survival of mice. Assuming a specific distribution that is able to describe the risk of observing an event, in this case death, at a specific time point a hazard function can be generated quantifying the risk of a mouse changing from alive to dead [18]. Of several existing parametrisations being used as hazard functions h(t), e.g., constant hazard distribution [19], Gompertz [19,20] and Weibull hazard distribution functions [19,21,22], log-logistic function [23], log-normal function [23,24], or surge function [25], the surge function was chosen as the most suited dynamic risk descriptor for the survival data: where SA is the surge amplitude, PT is the peak time, SW is the surge width at half-maximum intensity, and γ is a shape parameter of the peak. Variability in the hazard can be explained by potential covariates allowing differentiation between the study groups. Given the hazard function h(t), the survival function S(t) indicating the probability of being alive until a certain study time t can be derived. The same baseline hazard was assumed for all individuals, as only one observation per individual mouse existed due to the nature of the study and no IIV could thus be estimated. Potential covariates were included exponentially as change in hazard on overall survival to generate a link between the PK/PD model and the TTE model [18,26,27]: where βi (i=1, …, n) are coefficients characterising the extent of potential covariate effects of different covariate values xi (i=1, …, n). As PK predictors, mono-or combination therapy, the dose of AMX, Cmax of AMX, the area under the concentration-time curve (AUCAMX) of AMX, or the actual AMX concentration at a certain time t (C(t)) and respective combinations of single covariates were investigated. PD predictors included bacterial numbers at 12 h or 36 h in lung or in spleen (Clung, time, Cspleen, time), as well as the area under the bacterial number-time curve in lung or spleen (AUClung, AUCspleen) or the actual bacterial number at time t (CFU(t)). In addition, PK/PD parameters as T>MIC, T>EC50, the AUCAMX divided by the MIC (AUCAMX/MIC) or Cmax divided by the MIC (Cmax/MIC) were evaluated.

Pharmacokinetic Submodel
A two-compartment model for serum concentrations of AMX ( Figure S1), parameterized in terms of clearance (CLAMX = 124 mL/h), intercompartmental flow (Q = 71.9 mL/h), volume of distribution of the central (Vc = 15.4 mL), and peripheral (Vp = 50.7 mL) compartment with first-order absorption kinetics (ka = 5.04 h -1 ) including a lag time (tlag = 0.125 h) as well as first-order elimination kinetics best described the data. Parameter estimates for mice were plausible and are shown in Table S1. Investigation of other compartmental models or nonlinear absorption or elimination kinetics did not improve the model performance. Parameter estimates were estimated with acceptable RSEs (<25.5%) and a proportional RUV model (28.2% coefficient of variation). IIV was included on CLAMX (23.4% CV) and Q (25.7% CV). Prediction-corrected VPCs displayed a good predictive performance of the PK model over the entire investigated time period and also captured fractions of samples that were below the LLOQ ( Figure S2) as well as model estimates were within the 95% CI of the bootstrap analysis (Table S1). Still, higher variability was observed at higher AMX concentrations between 0 and 0.5 h after treatment ( Figure S2A).

Bacterial Disease Submodel
As described in Figure S1 (right), the final bacterial disease submodel consisted of a time-delay (klag = 0.0595 h -1 ), bacterial growth (kg = 0.477 h -1 ) and treatment-unrelated killing and natural death effects (kkill,lung = 0.274 h -1 ) in lung. The initial bacterial number of 6.12 log10 CFU/lung in lung was precisely estimated and in the range with the intranasally administered number of S. pneumoniae (≈6 log10 CFU). Here, a simultaneous and at the same time reliable estimation of k lag along with k g was challenging during the entire development process of the PK/PD model with a higher uncertainty especially for k lag (46.2%). However, given the underlying data of the bacterial growth submodel, simplifications or other implementations were not possible without neglecting important trends of initial killing or delayed growth in the data. Here, all parameters were essential to characterize the bacterial growth kinetics sufficiently, and hence a higher uncertainty was accepted.
The transit of bacteria from lung to spleen was best described by exactly estimating the number of transit compartments leading to a MTT of 42.0 h. Due to the marked delay in transition of bacteria determined in spleen about 3 h after treatment, a high number of compartments (n = 23) as well as the concomitant high MTT were plausible, especially considering the different barriers that bacteria had to overcome to transit from lung to spleen. In the spleen, the increase in bacterial numbers was only driven by the inflow from lung leading to accumulation of bacteria. The reliable estimation of bacterial growth as well as natural death and treatment-unrelated killing effects as a first-order killing rate constant kkill,spleen was not supported by the underlying data. Due to the setting of the study (e.g., study period of 48 h), bacteria reducing effects only occurred in presence of at least one drug.

Effect Compartment Submodel
Two separate effect compartments with first-order in-and outflow kinetics allowed a virtual transfer of AMX into lung (ke0, lung = 0.125 h -1 ) and spleen (ke0, spleen = 0.0435 h -1 ), ending up in a delayed increase and slower elimination compared to serum (t1/2,spleen = 15.9 h; t1/2,lung = 5.54 h; t1/2,serum = 0.0861 h) with lower maximum concentrations ( Figure S5). Hence, the antibiotic effect lasted over a longer time than AMX was actually above the MIC in serum.

Disease and Treatment Submodel
The effect of AMX on S. pneumoniae in the lung was best characterized by a sigmoidal Emax model as drug-dependent bacterial killing process. The concentration-effect relationship was steep indicated by the estimated high Hill factor, which was not quantifiable with high precision. Nevertheless, H was significantly higher than 1 as indicated by a LLP (95% confidence interval: 1.96-105) and was fixed to 20 to increase model stability [28]. Contrarily, a power model characterized by a first-order killing rate constant (kAMX) with a separate Hill factor representing the slope of the AMX effect compartment concentration and effect relationship in spleen was most stable and best captured killing effects in spleen. More complex implementations of the AMX effect in spleen, such as a sigmoidal Emax model, did not improve the model performance. Comparison of the AMX effect in lung and spleen did not only reveal that different mathematical implementations were better able to capture the effect of AMX in the respective compartments, but also the killing effect was substantially increased in spleen. In addition, no other killing effects were present within the structure of the PK/PD model in case of AMX monotherapy. Whereas in the lung, the effect of MPLA was implemented as its ability to enhance the efficacy of kkill,lung, in spleen the MPLA-related killing was estimated as a separate first-order killing rate constant and independent of natural death and treatment-unrelated killing effects. These were not able to be included in the spleen given manifold processes affecting the bacterial numbers in spleen. In case of combination therapy with MPLA neither the implementation of MPLA as influence on the maximum effect of AMX (∆E max , p = 0.944, df = 1, α = 0.05) nor on the potency of AMX (∆EC 50 , p = 0.272, df = 1, α = 0.05) revealed a significant change of the OFV in a LRT. In addition, the 95% CI of a LLP included parameter estimates of 1 for ∆E max (∆E max = 1.01; 95%CI 0.875-1.16) and ∆EC 50 (∆EC 50 = 0.772; 95% CI 0.415-1.23), respectively, and no improvement in GOF plots was observed indicating only additive effects.