Introduction

Cancer is considered a systemic disease with not only tumor lesions but also metabolic disorders, inflammation and other morbidities, which may lead to unexpected loss of body weight, aggravated toxicity and complications of anti-cancer therapies [1]. However, a lot of preclinical cancer studies absolutely choose tumor growth inhibition (TGI) derived from tumor size as a primary indicator to evaluate the in-vivo anti-cancer efficacy in animal xenografts, which might neglect other potential information related to safety issues and lead to a bias evaluation of anti-cancer therapies [2, 3]. Survival has been one of the most important endpoints in clinical cancer therapy, which may comprehensively reflect the overall result of disease progression, drug efficacy, toxicity, and their interactions [4]. Survival analysis in animal xenograft model that takes both drug effect and safety issues into consideration makes better use of experimental data and could provide more insights for specific therapy [5].

In addition to the commonly used techniques in survival analysis including nonparametric and semiparametric methods, such as Kaplan-Meier plot, log-rank test and Cox proportional hazard model [6, 7], parametric time-to-event (TTE) model is raising more and more interest, for its potential in associating time-course of important biomarkers with survival-related endpoints as well as its better ability in simulation [8, 9]. Combination of longitudinal tumor size model and TTE model is increasingly applied to facilitate decision making in oncology clinical trials and predicting the survival of individuals [10]. Nevertheless, except for tumor size, time-courses of other factors that may also relate to the survival are less considered. For example, nearly 80% of patients with advanced cancer experience progressive loss of body weight due to cachexia [11]. In preclinical studies, body weight of animals is also generally taken for a biomarker for systemic safety, individuals with significant change in body weight after being treated are considered to be not tolerated and not benefit from therapy [12]. However, there is little research investigating on the longitudinal body weight and its quantitative relationship with anti-tumor efficacy and survival.

Pancreatic cancer is a malignant tumor with 80%–90% risk of weight loss due to cachexia and an extremely high mortality [13]. In our previous study, dexamethasone (DEX) was found to be efficacious in inhibiting the growth and migration of pancreatic cancer in vitro and in vivo, where the blockade of NF-κB phosphorylation, reversion of EMT and down-regulation of IL-6 and VEGF were involved [14]. However, long-term and high dose of DEX may lead to muscle atrophy and aggravate the weight loss related with cancer cachexia [15,16,17].

In the current study, the efficacy and toxicity of DEX combined with gemcitabine (GEM), one of the first-line treatments for pancreatic cancer, was investigated in pancreatic cancer xenograft. Taking body weight as a systemic safety indicator, this study aimed to develop an integrated pharmacokinetic/pharmacodynamic (PK/PD) model to quantitatively demonstrating the impact of tumor size and systemic safety on animal survival.

Materials and methods

Chemicals and reagents

DEX and GEM were purchased from Meilunbio (Dalian, China). The DMEM medium was supplied by Macgene Biotech Co., Ltd (Beijing, China), fetal bovine serum (FBS) was provided by PAN-Biotech (Bavaria, Germany) and penicillin/streptomycin was provided by Gibco (CA, USA).

Cell culture

Human pancreatic cancer cell line SW1990 was supplied by American Type Culture Collection. Cells were cultured in DMEM medium with 10% FBS and 1% penicillin/streptomycin at 37 °C in a humidified 5% CO2 atmosphere.

Animals

The animal experiments were approved by the Institutional Animal Care and Use Committee of Peking University Health Science Center. Female BALB/c nude mice (5 weeks old) were provided by Department of Laboratory Animal Science, Peking University Health Science Center. The mice were maintained in individual ventilated cages with standard pathogen-free condition at 22–24 °C, 50%–60% humidity, and a 12 h/12 h light/dark cycles. Every effort was made to minimize the stress and the number of the animals used in the experiments.

PD study

In vivo study was carried out to investigate on the efficacy, toxicity as well as the impact on survival of DEX combined with GEM. SW1990 cells (3 × 106) were injected subcutaneously into BALB/c nude mice to establish the SW1990 xenograft model. The mice were randomly assigned into 8 groups and the treatment started when the tumor volumes reached 50–100 mm3. The dosage for DEX and GEM was selected according to our previous researches [14, 18, 19]. Considering the side-effect caused by long-term and high dose of DEX [16], DEX with longer dosing interval and lower total dose were also investigated in this study. Group A was set as vehicle control, Group B-D were treated with DEX 4 mg/kg by gavage every day/3 days/7 days respectively. Group E was given GEM 15 mg/kg every 3 days by tail vein injection, while Group F-H were given GEM 15 mg/kg every 3 days combined with DEX 4 mg/kg every day/3 days/7 days, respectively. DEX was dissolved in 10% hydroxypropyl-β-cyclodextrin (HP-β-CD) and GEM was in saline.

The two agents were given repeatedly until the death of animal or the end of study. In order to exclude the potential influence of solvent and administration on body weight or tumor size, 10% HP-β-CD was administered orally every day and saline was injected every 3 days for Group A. The length and width of the tumors were measured with caliper, and the body weight of the animals were recorded every other day. The tumor volumes were calculated as: tumor volume (mm3) = length × width2 × 0.5 [20]. Net body weight (NBW) was derived by subtracting tumor weight from animal’s body weight assuming that tumor has a density of 1 g/cm3, in order to exclude the numeric influence from the rapid-progressing tumor on the weight of other non-tumor parts of the animal. Finally, the survival of every mouse was recorded in a 60 days period, individuals that survived more than 60 days were considered as right-censored, which were euthanatized at the end of the study.

PK/PD model

The PK/PD models describing the time-course of both tumor size and NBW were developed in a sequential manner, respectively. Figure 1 showed the global scheme of the integrated model. Firstly, the PK model of DEX was described by a two-compartment model with first-order absorption rate, which was developed in our previous study [21]. Differential equations for the PK model are shown as follows:

$$\frac{{dX_{a,DEX}}}{{dt}} = - k_{a,DEX} \cdot X_{a,DEX},X_{a,DEX}(0) = X_{0,DEX}$$
(1)
$$\frac{{dX_{C,DEX}}}{{dt}} = \, k_{a,DEX} \cdot X_{a,DEX} \cdot F_{DEX} - Q_{DEX} \cdot C_{C,DEX} + Q_{DEX} \cdot C_{P,DEX}\\ - \, CL_{DEX} \cdot C_{C,DEX},X_{C,DEX}(0) = 0$$
(2)
$$\frac{{dX_{P,DEX}}}{{dt}} = Q_{DEX} \cdot C_{C,DEX} - Q_{DEX} \cdot C_{P,DEX},X_{P,DEX}(0) = 0$$
(3)
$$C_{C,DEX} = \frac{{X_{C,DEX}}}{{V_{C,DEX}}},C_{C,DEX}\left( 0 \right) = 0$$
(4)
$$C_{P,DEX} = \frac{{X_{P,DEX}}}{{V_{P,DEX}}},C_{P,DEX}\left( 0 \right) = 0$$
(5)

where Xa,DEX, XC,DEX, XP,DEX represent the amount of DEX in absorption, central and peripheral compartment, respectively. VC,DEX, and VP,DEX stand for apparent volume of distribution in central and peripheral compartment, while CC,DEX and CP,DEX are the concentrations of DEX in the relative compartment. ka,DEX is the first order absorption rate constant, CLDEX is the clearance of DEX and QDEX is the intercompartmental clearance. FDEX represents the bioavailability of DEX.

Fig. 1: The schematic structure of the integrated PK/PD model for the combination use of DEX and GEM.
figure 1

The schematic structure of the integrated PK/PD model for the combination use of DEX and GEM. Xa,DEX, XC,DEX, XP,DEX: amount of DEX in absorption, central and peripheral compartment; VC,DEX, VP,DEX: apparent volume of distribution in central and peripheral compartment; CC,DEX, CP,DEX: concentrations of DEX in central and peripheral compartment; ka,DEX: first order absorption rate constant of DEX; CLDEX: clearance of DEX; QDEX: intercompartmental clearance of DEX; FDEX: bioavailability of DEX; XC,GEM, XP,GEM: amount of GEM in central and peripheral compartment; VC,GEM, VP,GEM: apparent volume of distribution in central and peripheral compartment; CLGEM: clearance of GEM; QGEM: intercompartmental clearance of GEM; CC,GEM, CP,GEM: concentrations of GEM in in central and peripheral compartment; ET,max: maximum effect of DEX in tumor growth inhibiting; ECT,50: DEX’s concentration to achieve half of the maximum effect; KT,GEM: linear potency factor of GEM on damaging tumor cells; kt: transition rate between the tumor compartments; ψT: the combination index of DEX and GEM in tumor inhibition; k: first order rate of NBW decreasing; EW,max: maximum effect of DEX on NBW; ECW,50: concentration of DEX to achieve half of the maximum effect on NBW; KW,GEM: linear potency factor of GEM on NBW; ψW: combination index of DEX and GEM in NBW; kin: turnover rate of the Effect compartment (kin = kout).

Full size image

In addition, the two-compartment model of GEM’s PK was also obtained from publication [19]:

$$\frac{{dX_{C,GEM}}}{{dt}} = \, Q_{GEM} \cdot C_{P,GEM} - Q_{GEM} \cdot C_{C,GEM} - CL_{GEM} \cdot C_{C,GEM},\\ X_{C,GEM}(0) = \, X_{0,GEM}$$
(6)
$$\frac{{dX_{P,GEM}}}{{dt}} = Q_{GEM} \cdot C_{C,GEM} - Q_{GEM} \cdot C_{P,GEM},X_{P,GEM}(0) = 0$$
(7)
$$C_{C,GEM} = \frac{{X_{C,GEM}}}{{V_{C,GEM}}},C_{C,GEM}\left( 0 \right) = 0$$
(8)
$$C_{P,GEM} = \frac{{X_{P,GEM}}}{{V_{P,GEM}}},C_{P,GEM}\left( 0 \right) = 0$$
(9)

where XC,GEM, XP,GEM represent the amount of GEM in central and peripheral compartment, VC,GEM and VP,GEM stand for apparent volume of distribution in the relative compartment. CLGEM is the clearance of GEM and QGEM is the intercompartmental clearance. CC,GEM and CP,GEM are the concentrations of GEM in central and peripheral compartment, respectively.

The tumor natural growth model was characterized by Koch et al. [22], where the tumor growth was described by an exponential rate λ0 followed by a linear growth rate λ1 with no threshold point between two phases. The equation is shown as follows:

$$\frac{{dT}}{{dt}} = \frac{{2 \cdot \lambda _0 \cdot \lambda _1 \cdot T}}{{\lambda _1 + 2 \cdot \lambda _0 \cdot T}},T\left( 0 \right) = T_0$$
(10)

T indicates the tumor volume and T0 is the initial tumor volume at the beginning of treatment.

The PK/PD models describing the anti-tumor effect of DEX and GEM were developed separately and the typical values of the drug efficacy parameters were fixed when the model for their combination was established. DEX showed its effect by inhibiting the proliferation of tumor cells, which was similar with our previous study [23], and a Hill’s function was employed to describe it. GEM was assumed to convert the proliferating cells into non-proliferating cells, and the tumor cells were gradually brought to death through a series of transit compartments, describing the time delay between GEM exposure and tumor shrinkage [19]. Equations showing effect of DEX and GEM on the tumor size are shown as follows:

$$ \frac{{dT1}}{{dt}} = \frac{{2\lambda _0\lambda _1 \cdot T1^2}}{{(\lambda _1 + 2\lambda _0T1) \cdot T}} \cdot \left( {1 - \frac{{E_{T,max} \cdot C_{C,DEX}}}{{EC_{T,50} + C_{C,DEX}}}} \right) - K_{T,GEM} \cdot \psi _T \cdot C_{C,GEM} \cdot T1,\\ T1\left({t = 0}\right) = T_0$$
(11)
$$\frac{{dT2}}{{dt}} = K_{T,GEM} \cdot C_{C,GEM} \cdot T1 - k_t \cdot T2,T2\left( {t = 0} \right) = 0$$
(12)
$$\frac{{dT3}}{{dt}} = k_t \cdot (T2 - T3),T3\left( {t = 0} \right) = 0$$
(13)
$$\frac{{dT4}}{{dt}} = k_t \cdot (T3 - T4),T4\left( {t = 0} \right) = 0$$
(14)
$$T = T1 + T2 + T3 + T4$$
(15)

where T1 stands for the proliferating tumor cells and T2, T3, T4 stand for the non-proliferating cells, and T is the whole tumor including all these cells. ET,max is the maximum effect of DEX in tumor growth inhibiting with a theoretical upper limitation of 1, and ECT,50 is DEX’s concentration to achieve half of the maximum effect. KT,GEM is the linear potency factor of GEM on damaging tumor cells, and kt denotes the transition rate between the tumor compartments. ψT is the combination index of DEX and GEM in tumor inhibition, it is assumed that the two drugs showed synergistic effect if ψT was greater than 1 while they antagonize with each other if ψT was less than 1, and no interaction exists in PD level if ψT equal to 1.

The NBW of tumor-bearing mice was progressively decreasing since the tumor’s malignant proliferation required increasing nutrients from the host [24]. The unperturbed NBW decreasing was described with a first order rate as follows:

$$\frac{{dNBW}}{{dt}} = - k \cdot NBW,NBW\left( {t = 0} \right) = NBW_0$$
(16)

where k is the first order rate of NBW decreasing, and NBW0 is the initial NBW at the start of the treatment.

The models of DEX and GEM’s effect on the time-course of NBW were developed separately, too. Emax function was used to describe DEX’s effect on NBW, and GEM influenced the NBW in a linear way. The typical values of the drug-specific parameters were fixed in their combination model. The equation is shown as follows:

$$ \frac{{dNBW}}{{dt}} = k \cdot \left( {1 + \frac{{E_{W,max} \cdot C_{C,DEX}}}{{EC_{W,50} \,+\, C_{C,DEX}}} + K_{W,GEM} \cdot \psi _W \cdot C_{C,GEM}} \right) \cdot NBW,\\ NBW\left( {t = 0} \right) = NBW_0$$
(17)

where EW,max is the maximum effect of DEX on NBW, while ECW,50 is the concentration of DEX to achieve half of the maximum effect. KW,GEM is the linear potency factor of GEM on NBW, and ψW is the combination index of DEX and GEM in NBW.

Survival model

A parametric TTE model was developed to further explore the relationship between tumor size, NBW, and survival in the tumor-bearing mice. The hazard, cumulative hazard, and survival function are shown as follows:

$$h(t) = h_0$$
(18)
$$H(t) = {\int}_{\!\!\!0}^t {h(t)}$$
(19)
$$S(t) = e^{ - H(t)}$$
(20)

where h0 is the baseline hazard, and H(t) is the integration of h(t). S(t) stands for the probability of survival. In addition, probability density of observing death event was defined as the production of h(t) and S(t). Multiple probability density functions were tried for the baseline hazard including exponential, Weibull, log-logistic, and Gompertz distribution.

To investigate on the relationship between tumor size, NBW and survival, empirical Bayesian estimates were derived from the developed PK/PD models for tumor size and NBW, respectively, and longitudinal tumor size and NBW were considered as time-varing predictors in the TTE model. An indirect-response (IDR) model was used to fit the time-lag of drug effect better, which is assumed as a “Effect” compartment:

$$\frac{{{dE}}{{{{{\rm{ffect}}}}}}}{{dt}} = k_{in} - k_{out} \cdot Effect,Effect\left( {t = 0} \right) = 1$$
(21)

where kin and kout were equal, for the homeostasis was maintained when the initial Effect = 1 at time = 0 and NBW and tumor size remained stable. This homeostasis was interrupted by the implantation of the xenograft tumor and the treatment of DEX and GEM. DEX and GEM were assumed to influence the tumor size and NBW simultaneously and further affect the Effect compartment to work on the cumulative hazard, which derived survival:

$$\begin{array}{l}\frac{{{dE}}{{{{{\rm{ffect}}}}}}}{{dt}} = k_{in} \cdot \frac{{NBW\left( t \right)}}{{NBW_0}} - k_{out} \cdot \left( {1 + POW \cdot \log _{10}\left( {T\left( t \right)/T_0} \right)} \right) \cdot Effect,\\ t \, > \, 0,Effect\left( {t = 0} \right) = 1\end{array}$$
(22)

where NBW(t)/NBW0 represent the relative change of NBW at time = t verses the initial NBW. Log10(T(t)/T0) stands for the change of tumor size at time = 0, and POW is its potency factor on kout. Since the tumor size increased and NBW decreased with time, the value of Effect would gradually decrease from the initial value 1, and therefore h(t) increased until the death event occurred or censored.

Model development, evaluation, and external validation

The modeling and simulation in this study were conducted using NONMEM software (version 7.4.0, ICON Development Solutions, Ellicott City, MD, USA) and PsN (version 4.9.0, Uppsala University, Uppsala, Sweden). For the PK/PD models, first-order conditional estimation with interaction (FOCEI) method was used, and the inter-individual variability (IIV) was estimated in an exponential model. The residual variability was tested with additive, proportional, and mixture model. As for the parametric TTE model, Laplacian estimation method was used for parameter estimation. Since the variability of baseline hazard could not be estimated based on the survival data, all the individuals were assumed to have the same baseline hazard and therefore no IIV was estimated [5].

Models were evaluated and selected based on the rationality of the estimates, change in objective function value (OFV), and diagnostic plots. Relative standard error (RSE) was shown to evaluate the precision of the estimates and bootstrap based on 1000 runs was performed to access the parameter uncertainty of the tumor size model as well as the NBW model. On the other hand, sampling importance resampling (SIR) was applied to evaluate the robustness and parameter uncertainty of the TTE model. Besides, visual predictive check (VPC) based on 1000 times of simulation was executed for the PK/PD models and 100 times for the TTE model.

In order to avoid misleading exploration due to inaccurate estimation for hazard distribution and further improve the reliability of the parametric TTE model, additional two groups of tumor-bearing mice treated with different dosing schedules were employed as external validation set for the final model. Group I was treated with Schedule A (DEX 4 mg/kg at day 1, GEM 15 mg/kg at day 4, 7 days/cycle), and Group J was treated with Schedule B (DEX 1.33 mg/kg at day 1–3, GEM 15 mg/kg at day 4, 7 days/cycle). Tumor volume, NBW, and survival were recorded as mentioned above.

Simulations

The survival of the tumor-bearing mice under different NBW changes and tumor sizes changes were simulated to illustrate their impact on the survival of animals.

Results

DEX and GEM’s effect on tumor growth, NBW, and survival in SW1990 xenograft

The effects of DEX, GEM as well as their combination on tumor growth, NBW and the survival of SW1990 tumor-bearing mice were shown in Fig. 2. The number of mice in Group A (control) is 6 while other groups were 5, except a mouse died by accident in Group D. DEX inhibited tumor growth in a dose-dependent manner, while GEM significantly inhibited tumor growth both alone and in combination with DEX (Fig. 2a). NBW decreased with time for all the animals, and the NBW of mice treated with GEM decreased rapidly in spite of a rather slow tumor progression, except for Group H treated with DEX 4 mg/kg every 7 days combined with GEM 15 mg/kg every 3 days that showed a relatively stable NBW time-course after an obvious weight loss at the beginning of the treatment, and only Group H had significantly improved animals’ survival compared with control group (Fig. 2b, c), suggesting that the tumor sizes can not necessarily reflect survival results of tumor-bearing mice, and it would be more comprehensive to investigate on the survival as well as systemic safety biomarkers such as NBW of the animal xenografts.

Fig. 2: The effect of DEX, GEM and their combination on SW1990 xenograft.
figure 2

a Tumor growth, b NBW, and c the survival of SW1990 tumor-bearing mice. The statistical results are shown as *P < 0.05.

Full size image

PK/PD models for tumor size and NBW

The PK/PD models were developed sequentially and the parameter estimates, as well as the bootstrap 90% confidence intervals, are shown in Table 1. Some of the IIV in the final models were fixed to zero due to the limitation of data, resulting in an overestimation in the IIVs of the other parameters such as λ0 and KT,GEM, which were over 100%. For the tumor size model, the combination index ψT is 1.14 (90% CI is 0.78–1.42), indicating there was no obvious synergistic between DEX and GEM. On the other hand, the combination index in the NBW model ψW is 0.247, revealing the strong antagonism in the body weight loss between them. As the GOF plots shown (Supplementary Fig. S1), most of the population predictions and the individual predictions symmetrically distributed around the diagonal line, and the conditional weighted residuals (CWRES) were mostly within −4 and 4.

Table 1 Parameter estimates of the PK/PD models for tumor size and NBW.
Full size table

Individual fitting plots (Supplementary Fig. S2) indicated that the individual predictions basically agreed with the observations. The VPC results based on 1000 simulations were also performed (Supplementary Fig. S3). The medians of observations and simulations generally matched well, while the 10% and 90% percentiles of the simulations were wider than that of the observations, indicating the uncertainty of the model possibly due to the limited sample size. Overall, the PK/PD models were evaluated to be qualified to be integrated with the survival model. The main process for establishing the PK/PD models is summarized in Supplementary Table S1.

Parametric TTE model for survival

Empirical Bayesian estimates were derived from the aforementioned PK/PD models and the parametric TTE model for the survival of SW1990 tumor-bearing mice was further developed, the parameter estimates are shown in Table 2 and the main process for establishing the TTE model was summarized in Supplementary Table S1. The RSE of all the parameter estimates were under 30%, and the SIR results suggested that the parameters were estimated with reasonable uncertainty. The potency factor POW, representing the effect of tumor size on the assumed Effect compartment, was estimated to be rather small (9.4 × 10−7), suggesting that tumor size contributed less to the final survival results compared with NBW changes. The VPC results in Fig. 3 showed that most of the observed events were within the predicted 95% confidential interval, indicating that the model was with good ability to predict survival results in SW1990 tumor-bearing mice with the treatment of DEX and GEM, even though some misfits could still be found in Group D-H, which may be attributed to the limited sample size for model developing.

Table 2 Parameter estimates of the survival model.
Full size table
Fig. 3: Visual predictive check (VPC) for survival model for DEX combined with GEM in SW1990 xenograft.
figure 3

a–h Stand for related treatment groups. Group a was set as vehicle control; Group b–d were treated with DEX 4 mg/kg by gavage every day/3 days/7 days respectively. Group e was given GEM 15 mg/kg every 3 days by tail vein injection, while Group f–h were given GEM 15 mg/kg every 3 days combined with DEX 4 mg/kg every day/3 days/7 days, respectively. The blue solid line is the observations, the green shaded area stands for its 95% prediction interval.

Full size image

External validation

To further validate the predicting ability of the models, two groups of SW1990 xenografts with different dosing schedules were set as the external validations, where Group I contained 5 mice and Group J contained 4. For the PK/PD models, all the observed data were within the 95% prediction intervals (Fig. 4a, b). Post hoc Bayesian estimation was also performed for the external validation data set. The mean relative error (MRE) as well as root mean square error (RMSE) were calculated as 9.0% and 247.96 mm3 for the tumor size model, and 4.9% and 0.70 g for the NBW model, respectively. For the TTE model, the death events also fell into the prediction intervals (Fig. 4c). Basically, the external validation data sets met the model predictions well.

Fig. 4: External validation of the PK/PD models for the relative changes.
figure 4

Tumor size (a), NBW (b) and the TTE model (c) by Group I (left panel) and Group J (right panel). The shaded area represented the prediction interval, while the dots in a and b as well as lines in c represented the observations.

Full size image

Simulations

The survival of animals with different NBW reduction and tumor size progression was simulated based on the developed model. NBW reduction played an important role in animals’ survival in the current model, where individuals with less NBW reduction might have a much longer survival under the assumption that tumor size remains stable (Fig. 5). On the other hand, changes in tumor size made little difference in survival when the NBW were set to be stable, where the individuals with the same NBW reduction overlapped in spite of the different tumor ratio up to 50-fold change (Supplementary Fig. S4). These results suggested that there was a relatively weak connection between tumor size and survival in SW1990 xenografts, and NBW could possibly be a good predictor for the survival of tumor-bearing mice in preclinical studies of SW1990 xenograft model.

Fig. 5
figure 5

Simulated survival of SW1990 tumor bearing mice in different net body weight reduction.

Full size image

Discussion

Survival is accepted as a gold standard in oncology clinical trials for it could show the most clinically meaningful improvements from therapies [25]. Tumor inhibition has also been an important therapeutic target in both clinical and preclinical oncology studies. At preclinical stage, mice xenografts bearing cells or patient-derived tumors were commonly developed for researches in anti-tumor efficacy and TGI based on tumor size is often calculated to evaluate the efficacy of specific cancer therapy in mice xenografts [26]. However, it could be biased to only consider tumor size in oncology studies. It is reported that nutritional and psychological supportive care could significantly improve the prognosis for cancer patients [27], indicating that tumor size is not the only reason for cancer-related death. Similarly, it could be observed in preclinical studies that therapies significantly inhibiting tumor growth may also possibly came up with major safety issues and even lethal toxicity [28]. In this study, we developed an integrated PK/PD model considering both tumor size and systemic safety, in order to quantitatively investigate their impact on animal survival in pancreatic cancer xenograft.

Progressive weight loss is an important symptom of cancer cachexia which is often accompanied with advanced carcinoma, it has been reported that over 10% losses in body weight could be a significant risk factor for the survival of pancreatic cancer patients [29]. Besides, body weight of animals in preclinical studies is also considered as an indicator of cachexia as well as a systemic safety biomarker for certain therapy [30]. With the tumor size increasing, it would require more and more nutrition supplies from the host animals and also release more and more TGF-β, lactic acid and myostain, further leading to systemic inflammation, liver glycogen metabolism promotion and inhibiting skeletal muscle formation [24]. Although tumor weights are neglectable compared with human body weights, sometimes it could make a difference in animal xenografts where the increasing tumor weights may cover the body weight loss of the tumor-bearing mice, and therefore resulted in an underestimation of the systemic safety issues. In this study, NBW that derived from the body weight and tumor weight was selected as the systemic safety indicator to separate the tumor from the non-tumor part of the body weight, and the time course of NBW exhibited the weight loss caused by tumor progression together with anti-tumor therapies in a more direct manner.

In the PD study, DEX showed a dose-dependent tumor growth inhibition in SW1990 xenograft, which was consistent with our previous results, and therefore a Hill’s function was used to describing DEX’s inhibition on the proliferation of tumor cells based on the GR-related mechanism [14]. Linear function was used for GEM’s effect because there was only one dose level in the experimental data set. No PK interaction was assumed in the combination PK/PD model for the reason that DEX is mainly metabolized by CYP3A while GEM and its metabolites are excreted via renal, and it has been reported that no PK interaction was found between them in human [31]. In contrast to the synergistic effect in breast cancer xenograft [19], the combination index for tumor size model in SW1990 xenograft is slightly greater than 1, and the 90% confidence interval based on bootstrap indicated no significant synergistic effect was observed, which might be attributed to the relative lower expression of GR in SW1990 cells, and further resulting in a worse response to DEX therapy either alone or in combination with GEM [14].

Body weight was usually taken as a result of blood glucose changes in previous PK/PD modelling studies. For example, Iwasaki et al. established a semi-mechanism PK/PD model where exenatide influenced the body weight of mice via impacting their food intake [32]. However, the time course of body weight in tumor-bearing mice was seldomly quantitatively investigated. In this study, an empirical model for the longitudinal NBW data was developed by using a first order rate constant k to describe the NBW reduction caused by the increasing tumor size, which is simplified and fit-for-purpose. GEM sharply increased the NBW reduction rate with a potency factor KW,GEM of 11.6 mL/μg, while DEX slightly decreased the NBW by a EW,max of 0.369. The combination index ψW was 0.247, suggesting that the safety issue of weight loss could be improved when DEX and GEM were used in combination.

During the development of TTE mode, exponential, Weibull, log-logistic and Gompertz distributions for the baseline hazard were tested, and exponential distribution was selected in the final model for better fitness and lowest AIC (Supplementary Table S1). To integrate tumor size and NBW into the survival model, we firstly tried a proportional hazard style as follow:

$$h\left( t \right) = h_0 \times e^{\beta _1 \times X_1 + \beta _2 \times X_2 \cdots + \beta _n \times X_n}$$
(23)

where h0 is the baseline hazard, Xi is risk factor such as tumor size, NBW and their derived metrics, and βi is the associated potency factor. However, it did not provide a satisfying fitness. Considering the time lag between drug administration, tumor/NBW response and the cumulative hazard, the IDR model with an “Effect” compartment was utilized [5]. A difference with the model developed by Xie et al. was that the “Effect” compartment of IDR model was assumed to be a function to drug concentration in peritoneum and blood in that publication, which was considered to be related with the antitumor-efficacy and toxicity respectively [5]. In the current model, we assumed that both DEX and GEM showed their effect on tumor size and NBW, which were not only observed in the treatment groups but also in control group, so we tried this IDR model firstly in the control group and then in all the groups, rather than followed the sequence of “control group – single agent group – combination group” that we used in the development of tumor size model and NBW model.

Many different attempts were made to add tumor size into the survival model and the model was successfully minimized when the tumor size was added in a log10(T(t)/T0) manner with an extremely small potency factor POW near to zero and little improvement on fitness. The extremely small potency factor digitalized this neglectable influence of tumor size on the survival of tumor-bearing mice in this case, and therefore it was reserved in the final model to fit the purpose of the study better. Finally, the longitudinal relative changes of NBW and tumor size were integrated into the survival model as covariates. It should be emphasized that this structural model was limited to tumor-bearing mice, because it assumed that the hazard function was zero at hemostasis where the NBW and tumor size remained stable, which may not be extrapolated to other disease models. Besides, the hazard function h(t) should be limited to non-negative so as to avoid Lazarus distribution [33]. Since all the NBW predictions after treatment were less than the initial NBW and the tumor sizes were also greater than the initial volumes, no additional assumption was made to restrict h(t) to non-negative.

The simulation results showed that NBW played much more important role in the survival of the animals compared with tumor size. The reduction of NBW not only came from drug effect but also the increasing demand for nutrition of rapid-progressing tumor with not any treatment [34], the relationship between NBW and the increasing tumor size was actually neglected when an empirical first order rate was used to capture the data. However, a more complicated model could be unstable due to the limitation of our data, as a result it was not taken into consideration in the current study. On the other hand, NBW could be more directly associated with the survival of tumor-bearing mice than tumor sizes, suggesting that tumor-growth inhibition is not the only thing should be cared when evaluating the therapeutic efficacy in preclinical studies.

One of the limitations of this study is that the results were obtained from immunodeficient mice bearing SW1990 cells derived tumors, while other covariates except NBW might be more important to the survival results in other animal models and tumors, such as the immune status of the host, which required further investigation in more different xenografts. In addition, the PK models for DEX as well as GEM were adopted from previous publications [19, 21], and the individual PK for the mice was not considered in the PK/PD model, which may lead to over-estimation for the variability of the parameters in the PK/PD models. Furthermore, it would be a more complicated situation when it comes to clinical cancer therapy, and even more risk factors that may associate with the clinical survival endpoints need to be taken into consideration. Last but not least, the small sample size for model development and external validation was a significant limitation in the current study, and there were still some misfits in the model possibly due to the limited data, which may bring some uncertainties and keep us modest with the conclusions.

Nevertheless, the current study has revealed the importance of body weight monitoring and survival analysis in preclinical oncology studies, which was often left in the basket in the past but actually worth considering in future preclinical studies, and the proposed model structure with DEX and GEM could provide some reference to other modeling and simulation studies in preclinical stages. For other anti-tumor drugs with different mechanisms, tumor size and NBW might not possibly be the only factors connecting drug efficacy, toxicity and survival outcomes in preclinical studies. In future studies, the applicability of the model would be tested in various situations including different anti-cancer drugs and different cancer xenograft models, and the workflow of parametric survival analysis in preclinical studies could be further explored and optimized.

Conclusion

In summary, a parametric TTE model of pancreatic cancer SW1990 xenografts’ survival integrated with longitudinal PK/PD models based on tumor sizes and NBW was developed and validated. NBW showed a more direct and important relation with the survival results of SW1990 tumor-bearing mice compared with tumor sizes, suggesting that it could be biased to evaluate the anti-cancer therapeutic efficacy with tumor growth inhibition alone, and survival analysis together with systemic safety biomarkers such as body weight might deserve more attention in preclinical oncology studies.