Mechanistic Multiscale Pharmacokinetic Model for the Anticancer Drug 2',2'-difluorodeoxycytidine (Gemcitabine) in Pancreatic Cancer.

The aim of this work is to build a mechanistic multiscale pharmacokinetic model for the anticancer drug 2',2'-difluorodeoxycytidine (gemcitabine, dFdC), able to describe the concentrations of dFdC metabolites in the pancreatic tumor tissue in dependence of physiological and genetic patient characteristics, and, more in general, to explore the capabilities and limitations of this kind of modeling strategy. A mechanistic model characterizing dFdC metabolic pathway (metabolic network) was developed using in vitro literature data from two pancreatic cancer cell lines. The network was able to describe the time course of extracellular and intracellular dFdC metabolites concentrations. Moreover, a physiologically-based pharmacokinetic model was developed to describe clinical dFdC profiles by using enzymatic and physiological information available in the literature. This model was then coupled with the metabolic network to describe the dFdC active metabolite profile in the pancreatic tumor tissue. Finally, global sensitivity analysis was performed to identify the parameters that mainly drive the interindividual variability for the area under the curve (AUC) of dFdC in plasma and of its active metabolite (dFdCTP) in tumor tissue. From this analysis, cytidine deaminase (CDA) concentration was identified as the main driver of plasma dFdC AUC interindividual variability, whereas CDA and deoxycytidine kinase concentration mainly explained the tumor dFdCTP AUC variability. However, the lack of in vitro and in vivo information needed to characterize key model parameters hampers the development of this kind of mechanistic approach. Further studies to better characterize pancreatic cell lines and patient enzymes polymorphisms are encouraged to refine and validate the current model.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

☑ 2’,2’‐difluorodeoxycytidine (gemcitabine, dFdC) treatment is associated with high variability in responses. Several studies suggest that individual genetic factors affecting dFdC metabolic pathway would lead to different concentrations of dFdC metabolites and consequently to different responses.

WHAT QUESTION DID THIS STUDY ADDRESS?

☑ It integrates knowledge from the literature into a mechanistic multiscale pharmacokinetic (PK) model describing dFdC metabolic pathway. The model simulates different concentration profiles of dFdC metabolites in the tumor for a virtual population of patients with pancreatic cancer.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

☑ It provides a platform for translating in vitro dFdC metabolites profiles to in vivo concentrations in humans at target sites, coupling together system PK and physiologically‐based pharmacokinetic (PBPK) modeling.

HOW MIGHT THIS CHANGE CLINICAL PHARMACOLOGY OR TRANSLATIONAL SCIENCE?

☑ This study shows that integrating in vitro metabolic network in in vivo PBPK models can be useful for the prediction of the active metabolites levels in the drug site of action. Data including the target enzyme level expression, the degree of their polymorphisms, and their interindividual variability are encouraged to refine the current version of the model.

2′,2′‐difluorodeoxycytidine (gemcitabine, dFdC) is a nucleoside antimetabolite prodrug effective against several solid tumors.1, 2, 3, 4 Treatment with dFdC represents the first‐line therapy of pancreatic cancer, that constitutes one of the most aggressive and lethal oncology diseases, with an overall 5‐year survival rate of < 5%.5

As a prodrug, dFdC has to be intracellularly metabolized to its active metabolite, dFdC triphosphate (dFdCTP), to exert its cytotoxic action.6 First, dFdC is taken into the cell by active transporters (hENTs and hCNTs)7 and then it is phosphorylated by deoxycytidine kinase (dCK) to its monophosphate form, dFdCMP. Then, dFdCMP is subsequently metabolized by nucleoside kinases to dFdC diphosphate (dFdCDP) and then to dFdCTP that binds to the DNA promoting apoptosis.8 The dFdC also suffers inactivation by cytidine deaminase (CDA), leading to inactive metabolite 2′,2′‐difluorodeoxyuridine (dFdU), which is excreted in urine.9

One of the biggest complications associated to treatment with dFdC is the variability in responses, ranging from lack of efficacy to severe toxicity.10 These different rates of responses to dFdC could be, in part, explained by individual genetic factors affecting its metabolic pathway, leading to different dFdCTP intracellular tumor concentrations.11 As an example, a high activity of CDA enzyme is related with a higher depletion of dFdC and, therefore, lower dFdCTP concentrations.12 It is also stated that treatment efficacy may be explained by a nonfunctional transport of the prodrug into the cell.7 Moreover, cells with low dCK levels are associated with resistance to dFdC.13, 14 In addition, some clinical studies in patients with pancreatic cancer treated with dFdC associated different expressions of the transporters or the target enzymes activity with a high or low survival probabilities.13, 15

The dFdC effects on pancreatic cancer have been described previously by using pharmacokinetic/pharmacodynamic (PK/PD), mechanistic, and semimechanistic models in in vitro,16, 17, 18, 19 preclinical in vivo,20, 21 and clinical stages.22, 23 However, to the best of our knowledge, the models developed in clinical stages do not consider the intersubject variability in dFdC metabolism (e.g., individual concentrations of enzymes involved in dFdC metabolism).

Information on the systemic and cellular dFdC PKs could be used to develop a quantitative model that mechanistically describes processes, such as drug distribution, metabolism, and active metabolite formation. A model like this could help in improving the knowledge of the system, for example, by understanding what subject characteristics mostly determine the predicted active metabolite exposure in the site of action.

Systems pharmacology is an approach that aims to develop multiscale mechanistic models that “span the divide between cell‐level biochemical models and organism‐level PK/PD models.”24 These models integrate the knowledge from various sources (e.g., in vitro experiments and physiological data),25 and typically represent the system more mechanistically than classical PK/PD models. The physiologically‐based pharmacokinetic (PBPK) modeling approach provides a framework for integrating drug‐specific parameters and in vitro measurements with physiological system‐specific parameters.26 This type of models can integrate interindividual variabilities in the concentrations of enzymes involved in drug metabolism and allows the simulation of drug concentration in specific tissues (e.g., pancreatic tumors).27

In this context, we built a mechanistic systems pharmacology model to describe dFdC PKs and dFdCTP tumor concentrations, in a population of patients with pancreatic cancer. The developed model was built by using data from the literature, including genetic and physiological intersubject variabilities. In summary, our aims were to (i) propose a translational multiscale system PK modeling approach for dFdC able to describe different concentrations of dFdC metabolites, (ii) to show the capabilities and the limitations of this kind of modeling strategy starting from a case study, and (iii) to understand what information is needed and what can be found or not in the literature.

Methods

The work was performed in three different steps. First, we developed a model (the so‐called metabolic network) to describe dFdC metabolism in vitro in two pancreatic cancer cell lines (PK9 and RPK9).14 Then, we developed a PBPK model to describe the dFdC PKs in a population of patients with pancreatic cancer. The in vitro derived metabolic network was coupled with the PBPK model in a compartment representing the pancreatic tumor after an appropriate in vitro to in vivo rescaling of the network parameters. This was done to describe the dFdC metabolism and predict dFdCTP concentrations in the site of action. Finally, we performed a global sensitivity analysis (GSA) on the developed model to identify what are the characteristics that mostly drive the dFdC and dFdCTP exposure variability in a population of patients.

The analyses were performed in MATLAB R2019a.28 Parameters were estimated by using the covariance matrix adaptation evolution strategy.29

Metabolic network

An extensive literature review was performed looking for knowledge and in vitro data to build a mathematical model of the dFdC metabolic pathway. The structure of this pathway, which has been defined over the years,6, 8, 9 is schematized in Figure 1 c. The mathematical model was built by assuming that enzymatic reactions were described by first order rate constants, except for those catalyzed by CDA, dCK, and hENT1 enzymes, which were described by Michaelis‐Menten kinetics. Metabolic network equations are reported in the Supplementary Material, Section .

Figure 1

Schematic structure of the whole body physiologically‐based pharmacokinetic (PBPK) model coupled with the metabolic network representing gemcitabine (dFdC) metabolism in the pancreatic tumor tissue. (a) PBPK model structure. Red arrows represent arterial blood flows, whereas blue arrows represent venous blood flows. The organs and tissues represented in boxes A and B are “in parallel” with respect to the blood flow, this means that they have separate blood inflows and outflows. Box A: Adipose tissue, bone, brain, gonads, heart, kidneys, muscles, and skin. Box B: Gut, spleen, and stomach. (b) Model structure of the pancreas. Pancreatic tumor and intracellular space share the same extracellular environment. (c) Schematic representation of dFdC metabolism network, including metabolites (dFdC, dFdCMP, dFdCDP, and dFdCTP), transporters (hENT1), and target enzymes responsible of driving the metabolism reactions (dCK, NMPK, NDPK, CDA, and dCMPD). Reactions catalyzed by unknown enzymes were named as KMPC, KDPMP, and KTPDP. KDNA represents the DNA binding rate of dFdCTP. This figure was modified from Servier Medical Art, licensed under a Creative Common Attribution 3.0 Generic License. http://smart.servier.com/.

Experimental data used to identify network parameters were taken from Ohmine et al.14 In vitro concentrations of dFdC metabolites (extracellular dFdC and dFdU, intracellular dFdC, dFdCMP, dFdCDP, dFdCTP, dFdU, and dFdUMP) for two pancreatic cancer cell lines (i.e., PK9 and its resistant version to dFdC, RPK9) were available.

Parameters were jointly estimated on both cell lines data by including the ratio of the three target enzymes (CDA, dCK, and hENT1) concentrations between the two cell lines as covariates of the model. , , and , the covariates, were set equal to 1 for PK9 and equal to 1.64, 0, and 1.35 for RPK9, respectively.14 However, with the available data, it was only possible to identify the parameters of a reduced metabolic network involving dFdC, dFdCMP, dFdCDP, and dFdCTP, but not dFdU and dFdUMP. In addition, to fit data, an extracellular dFdC deamination due to the activity of an extracellularly secreted CDA was added. The hypothesis is supported by some observations reported in the literature for other cancer cell lines.30

To correctly write mass balance equations, metabolite measurements were transformed from the concentration values (pmol/mg prot) reported in Ohmine et al. 14, to amount (pmol), as explained in the Supplementary Materials, Section .

PBPK model

A PBPK model was developed to describe dFdC distribution and metabolism in the body. The model, represented in Figure 1 a, consists of 14 organs and tissues; each of them (excluding arterial and venous blood) was described by a permeability limited model.31 This choice was supported by the hydrophilic nature of dFdC hampering distributions into the cells.7, 32 Drug‐specific parameters used in the model are listed in Table 1.

Table 1

Parameters values and their distributions in the population as used in the GSA

Parameters names	Distributions parameters	Distributions type	Units
pKa35	3.6	Fixed
B:P47	1.94	Fixed
fup35	1	Fixed
Molecular weight35	299.66	Fixed	g/mol
logPow35	–1.4	Fixed
Sexc,d	0, 1	Uniforma
Aged	20, 65	Uniforma	Years
E:Pe	1, 5	Uniforma
Tumor volume	32.3, 224.3	Uniforma	mL
k_hCNT1	920.1657 (33%)	Log‐normalb	1/minute
k_hENT1,in	20.1560 (24.3%)	Log‐normalb	1/minute
k_hENT1,out	25.6623 (24.3%)	Log‐normalb	1/minute
k_CDA	0.3312 (109.6%)	Log‐normalb	1/minute
dCK_int,vivo	0.45 (20%)	Log‐normalb	pmol/mg prot
hENT1_int,vivo	3.08 (53.4%)	Log‐normalb	pmol/mg prot
CDA_int,vivo	0.67 (109.6%)	Log‐normalb	pmol/mg prot

GSA, global sensitivity analysis.

For distribution parameters, minimum, maximum of the parameter. bFor distribution parameters, mean (coefficient of variation) of the log‐normal variable. cIf the extracted value is < 0.5 the subject is female (0), otherwise male (1). dHeight, organ volumes, and blood flows were generated by using the Willmann model37 for the population generation and are function of sex and age. We refer to the original publications for their values and distributions. eB:P calculated from E:P values.47 B:P was calculated by using the following equation48: , with H hematocrit equal to 0.47.

The dFdC is transported inside the cell by concentrative nucleoside transporters (mainly hCNT1) and equilibrative nucleoside transporters (mainly hENT1).7, 33 The activity of both hCNT1 and hENT1 was modeled as a first order reaction, namely and (Eqs. 1 and 2). The hCNT1 mediates a unidirectional flux from the extracellular to the intracellular space, whereas hENT1 was considered as a bidirectional transporter. Once inside each organ intracellular space, the drug was supposed to be metabolized by CDA. This process was modeled in Eq. 3 with a first order reaction too ().

is the relative expression of enzyme or transporter x in the tissue t. They were taken from the Open Systems Pharmacology Suite version 7.134 and are numbers always between 0 and 1; their values are reported in Table . The is the rate constant relative to the protein x activity, assumed to be equal in all the organs. and are the unbound extracellular and intracellular dFdC concentrations, respectively. and are the extracellular and intracellular volumes of the tissue t, they are calculated by multiplying the tissue volume for the extracellular and intracellular water fractions. Unbound fraction of dFdC was considered equal to 1.35 To account for the different transporters and enzymes expressions on each organ, was multiplied for . 27 The main hypothesis here is that the intracellular enzymes’ and transporters’ concentrations are proportional between the different organs. In Eqs. 1 and 2 the rate constants relative to the transport from the extracellular to the intracellular compartment are multiplied for . The rate could be written as , where and are the parameters of the Michaelis‐Menten equation. Therefore, the following relationship holds:

and are the intracellular and extracellular transporter concentration and is the turnover number. Equation (4) applies to both hCNT1 and hENT1.

Each organ modeled in the PBPK was described by using two compartments, representing intracellular and extracellular spaces. The generic dFdC tissue extracellular and intracellular unbound concentration dynamics are represented in Eq. (5).

is the tissue blood flow, corresponds to the arterial dFdC concentration, is the tissue to plasma partition coefficient, calculated as in Jamei et al.,31 and is the blood to plasma partition coefficient. These equations apply to all the tissues except arterial and venous blood, lungs, and pancreas.

After appropriate parameters rescaling, the metabolic network describing dFdC metabolism was included into a compartment representing the pancreatic tumor, as explained in Inclusion of the metabolic network in the PBPK model section below. All the PBPK model equations and parameter values are reported in the Supplementary Materials, Sections S1 and S6.

The four rate constants associated with enzymes and transporters activities in the PBPK model (, , and , in Eqs. (1), (2), (3)) were identified for a male mean subject (height 175 cm, weight 73 kg, and age 30 years) using the data simulated with a PK model taken from the literature (Zhang model, Supplementary Materials, Section ), for 30‐minute infusion of 3.34 mmol/m2 of dFdC.36

Finally, by adding variability to the physiological parameters, a population model was obtained. The variabilities in organ volumes and blood flows were modeled following Willmann et al.37 Briefly, in this model: (i) sex and age of the subjects are extracted; (ii) for each subject, the height is extracted from a distribution given sex and age; (iii) mean organs weight and blood flows are generated given the mean subject characteristics; and (iv) a residual variability is added. The , , , and were considered variable too, in order to account for the different protein expression among subjects. They were supposed log‐normally distributed with mean equal to the estimated values and coefficient of variation (CV) taken from Burt et al.38 To our knowledge, no information regarding CDA variability in tissues is present in the literature; thus, we decided to fix the CDA CV equal to the one of the pancreatic cancer, obtained from Ohmine et al.39 Parameters of the distributions are reported in Table 1.

Inclusion of the metabolic network in the PBPK model

A compartment representing the pancreatic tumor cells was included into the PBPK, as shown in Figure 1 b. The main hypothesis is that the tumor and the pancreatic intracellular space share the same extracellular environment and they compete for drug uptake. The network was included into the PBPK as follows: in vitro intracellular compartment corresponds to the PBPK tumor compartment, whereas the in vitro medium corresponds to the pancreatic extracellular space.

The parameters were appropriately rescaled considering the different volumes and enzymatic concentrations between the in vitro and in vivo situations, as described in the Supplementary Material, Section . In the publication where we took the in vitro data, the in vitro concentrations of some of the enzymes and transporters involved into dFdC metabolism, like CDA, dCK, and hENT1, were also reported.14 In another publication, the concentration of these enzymes in pancreatic tumor samples from 10 different subjects was measured.39 That information was used for the in vitro to in vivo rescaling.

Tumor volume was supposedly distributed in the population uniformly between 32.3 mL and 224.3 mL.40 Enzymatic concentrations were supposedly log‐normally distributed with the mean and CV derived from the pancreatic tumor samples.39

Global sensitivity analysis

Variance‐based GSA can be used to understand how the variability of certain model outputs (e.g., area under the curve (AUC)) is apportioned to the variability of the model parameters (e.g., organ volumes and enzymes abundance), indicated in GSA as inputs.41 This technique can help to understand the most important parameters responsible for the model outputs variability.

In variance‐based GSA, for each parameter , two sensitivity indices are derived from the output variance decomposition: the main effect and the total effect . is related to the part of the output variance explained by the variation of each taken alone, while is plus the interactions between parameters. Both and are always between 0 and 1. The highest and are, the more important the parameter is, whereas if is equal to 0, the parameter does not impact on the output variability.41

Variance‐based GSA was performed on the dFdC PBPK model coupled with the metabolic network. AUC of plasma dFdC and tumor dFdCTP concentrations were considered as outputs of interest. AUC was calculated from time 0 (dose administration time) to 7 days. The parameters that were considered different among subjects of a population are: sex, age, and height, organ volumes and blood flows, dFdC blood to plasma ratio, all the rates associated with drug transport and elimination in the PBPK, tumor volume, and tumor concentrations of the enzymes involved in dFdC metabolism. Given that and are related to the activity of the same enzyme, in the GSA, they were jointly considered (grouped): they shared the same variability and so they were considered as a unique parameter (). For readability purposes, all the organ volumes and blood flows residual variabilities were grouped too.41 The distributions of all the parameters are reported in Table 1.

The number of samples, n, extracted in the GSA was set to 5,000. The confidence intervals of the sensitivity indices were calculated using 10,000 bootstrap samples.42

Results

In vitro metabolic network

The parameters of the metabolic network were identified on the in vitro data from Ohmine et al.14 The estimated parameter values are listed in Table 2, where it is also indicated for which parameter covariates were included. In order to reduce the number of parameters to identify, and were imposed to be equal and , , and were considered equal too. In Figure 2, the results of the fitting process for each metabolite profile and each pancreatic cell line are shown. R 2 values for extracellular dFdC, intracellular dFdC, dFdCMP, dFdCDP, and dFdCTP are equal to 0.98, 0.1, 0.99, 0.61, and 0.87, respectively. The metabolic network was believed sufficiently capable of describing the time course of extracellular and intracellular concentrations of dFdC and its phosphorylated metabolites (dFdCMP, dFdCDP, and dFdCTP) for the two different cell lines.

Table 2

Metabolic network estimated parameters

Parameters names	Value	Units
V_max,dCK (xCOVdCK)	1.4470 × 105	pmol/hour
K_M,dCK a	4.6	µM
V_max,CDA,int (xCOVCDA)	1.1047 × 105	pmol/hour
K_M,CDA,int	0.4340 × 105	pmol
V_max,CDA,ext (xCOVCDA)	0.8382 × 105	pmol/hour
K_M,CDA,ext	4.1976 × 105	pmol
V_max,hENT1(xCOVhENT1)	5.4574	pmol/hour
K_M,hENT1	4.6963	pmol
K_NMPK	0.1087 × 105	1/hour
K_DCMPD	0.0273 × 105	1/hour
K_DPMP	0.1393 × 105	1/hour
K_DNA b	1 × 10–7	1/hour
K_INH	0.8072 × 105	1/pmol

Value taken from Bouffard et al.49 To include it into the model the value was multiplied for the intracellular volume. bGiven that its value is very close to zero, it was removed from the model without an impact on the simulations.

Figure 2

Fitting results of the in vitro metabolic network for PK9 and RPK9 pancreatic cancer cell lines. Blue lines are the model predictions and red stars the data from Ohmine et al.14 In panel (a), results for gemcitabine (dFdC) medium amount for PK9 are reported. In panels (b–e), the results for dFdC, dFdCMP, dFdCDP, and dFdCTP intracellular PK9 amounts are reported. In panel (f), the results for dFdC medium amount for RPK9 are reported. In panels (g–j), the results for dFdC, dFdCMP, dFdCDP, and dFdCTP intracellular RPK9 amounts are reported.

PBPK

We fitted the PBPK model coupled with the reduced metabolic network against a typical plasma profile of dFdC, given a single dose of 3.34 mmol/m2 infused in 30 minutes (standard administration in the clinical setting), simulated with the Zhang model.36 The identified parameters are reported in Table 1 (mean values of the rates) and the fitting results are shown in Figure 3. It is possible to observe that the PBPK model well reproduces the typical subject profiles provided by the Zhang model. From the logarithmic scale, it can be appreciated that the elimination rate is well captured. The value of R 2 is equal to 0.94.

Figure 3

Fitting and population simulation results of the physiologically‐based pharmacokinetic (PBPK) model. In panels (a, b), the result of the fitting process is shown, in natural and semilogarithmic scale, respectively. Continuous blue line represents the gemcitabine (dFdC) plasma concentration simulated by using the PBPK model, whereas the dashed red line represent the dFdC plasma concentration simulated with the Zhang model.36 In panel (c), the dFdC plasma concentrations are reported. In panels (d–g), the tumor dFdC, dFdCMP, dFdCDP, and dFdCTP concentrations are reported, respectively. In this case, blue line represents the median of the compound concentrations in the population, while red shaded area represents the 95% confidence interval.

Once the model parameters were identified on the typical profile, PK profiles of a population of 500 individuals were simulated. In Figure 3, the plasma dFdC concentration profiles and the dFdC and its metabolites pancreatic tumor profiles are shown, whereas profiles for all the other organs and tissues are reported in the Supplementary Materials, Section .

In Figure 3, it is possible to observe that although the dFdC plasma concentration drops to zero for at least the 95% of the subjects in almost 24 hours, the metabolite concentrations in the tumor decreases much slower. These results are qualitatively in agreement with the simulations obtained by the Zhang model, that predicts a drop to zero in 70 hours for the typical value of the dFdCTP concentration in the white blood cells (WBCs), used as a surrogate of the intracellular dFdCTP concentration.

In Table 3, the metrics simulated with the PBPK model for plasma dFdC and tumor dFdCTP concentrations are compared with those simulated with the Zhang model and with in vivo patients PK data from the literature.43, 44 For the two in vivo studies, metrics related to the tumor dFdCTP concentration were not available, thus, the ones relative to WBC dFdCTP concentration are reported. Concerning the dFdC plasma AUC and peak plasma concentration (Cmax), the results of the PBPK model show good agreement with both the Zhang model and the literature data. However, it is possible to observe that the PBPK model tends to overestimate the population variability of dFdC AUC. PBPK dFdCTP tumor AUC results slightly lower with respect to dFdCTP WBC AUC simulated with the Zhang model, but it shows good agreement with the results presented in Derissen et al.43

Table 3

PBPK and Zhang model output metrics compared with literature metrics

Metrics	Study	Mean	Mina	Maxa	CV (%)b		Units
		dFdC
Plasma AUC	PBPK model	3.57	1.14	11	78.06	mmol·min/L
Plasma AUC	Zhang modelc	2.21	1.22	3.59	28.92	mmol·min/L
Plasma AUC	Derissen et al. 43	1.73	0.67	3.42		mmol·min/L
Plasma AUC	Abbruzzese et al. 44	2.13	0.76	5	78.12	mmol·min/L
Plasma Cmax	PBPK model	0.048	0.031	0.071	22.99	mmol/L
Plasma Cmax	Zhang modelc	0.041	0.024	0.061	22.14	mmol/L
Plasma Cmax	Derissen et al. 43	0.043	0.017	0.082		mmol/L
Plasma Cmax	Abbruzzese et al. 44	0.056	0.012	0.108	63.24	mmol/L
		dFdCTP
Tumor AUC	PBPK model	343.86	31.18	1,127.3	94.37	mmol·min/L
WBC AUCd	Zhang modelc	564.95	331.25	931.17	29.12	mmol·min/L
WBC AUC	Derissen et al. 43	394.8	169.8	942		mmol·min/L
Tumor Cmax	PBPK model	0.15	0.07	0.3	42.99	mmol/L
WBC Cmax d	Zhang modelc	0.44	0.26	0.68	25.75	mmol/L
WBC Cmax	Derissen et al. 43	0.5	0.19	1.06		mmol/L
WBC Cmax	Abbruzzese et al. 44	0.224				mmol/L

AUC, area under the curve; Cmax, peak plasma concentration; CV, coefficient of variation; PBPK, physiologically‐based pharmacokinetic; WBC, white blood cell.

Minimum and maximum refer to: 2.5 and 97.5 percentiles for the simulations with PBPK and Zhang models; observation ranges for Derissen et al. 43 and Abbruzzese et al. 44. bFor all the metrics relative to Derissen et al. it was not possible to calculate the CV. For the WBC Cmax of Abruzzese et al. ranges and CV were not available. cFive hundred subjects were simulated (250 men and 250 women, age 29 years) including intersubject variability reported in Zhang model. Mean subject metrics were calculated for a 29 year old man. dValues of dFdCTP AUC and Cmax were simulated with the Zhang model in and , respectively. These values were converted to intracellular WBC concentration by dividing them for the mean volume of a neutrophil, equal to 299 fL.50

As a further simulation exercise, 3.34 mmol/m2 of dFdC were infused weekly for 20 weeks; results are shown in the Supplementary Materials, Section . From these results, it is possible to see that there is no accumulation of dFdC and dFdCTP, in agreement with De Lange et al. 45

A variance‐based GSA on the PBPK model coupled with the metabolic network was performed with the aim of understanding what the most important parameters are in explaining plasma dFdC and tumor dFdCTP concentration AUC variability in the population. Results are shown in Figure 4.

Figure 4

Global sensitivity analysis (GSA) results performed for (a) plasma dFdC area under the curve (AUC) and (b) tumor dFdCTP AUC. The parameters that are considered variables among subjects in the population are: sex, age, residual variability on height, organ volumes, and blood flows (height, volumes, and fluxes), blood to plasma ratio (BP), tumor volume (tumor vol), all the estimated rate constants relative to the transport and elimination of dFdC in the physiologically‐based pharmacokinetic model (kCNT pbpk, khENT1 pbpk, and kCDA pbpk) and the enzymes tumor concentrations (CDA tum, dCK tum, and hENT1 tum). Error bars represent the 95% confidence interval of the sensitivity indices calculated with 10,000 bootstrap samples.

The parameter that mainly explains the dFdC plasma concentration AUC is the rate relative to the dFdC elimination in tissues, . A critical aspect related with this parameter is that we have considered its CV equal to the one found in pancreatic tumor tissue samples, because of the lack of specific knowledge. The variation of is responsible of about all the dFdC AUC variance; then an adequate characterization of its variability in the population is needed for a more reliable prediction of the AUC variability.

The dFdCTP tumor concentration AUC variability is mainly due to dCK and CDA tumor concentrations. These results are in agreement with the fact that the resistance against dFdC in some cell lines is obtained by reducing the dCK levels.14

It is interesting to notice that hENT1 tumor expression is not important in determining dFdCTP AUC. This could be due to the fact that the hENT1 concentration was quite homogeneous in the population that we used to estimate its variability.39 Thus, it is possible that with these data, the hENT1 population variability was underestimated.

Discussion

A multiscale systems pharmacology model describing the dFdC metabolic pathway and predicting the levels of dFdCTP in the active site for a population of patients with pancreatic cancer was developed. This model was built by integrating different resources obtained from the literature: in vitro information regarding dFdC metabolism in pancreatic cancer cell lines14; a compartmental model describing plasma dFdC concentrations in patients with pancreatic cancer36; and physiological and genetic information of a given population.37 Finally, we performed a GSA in order to understand what parameters explain the predicted population variation of plasma dFdC and tumor dFdCTP AUC.

Regarding the dFdC in vitro metabolism, data used for developing the network were obtained from in vitro experiments performed after a single dose exposure of dFdC, collecting a single profile per metabolite for each cell line. With the data available, the results of the fitting process are biased, as it can be appreciated by looking at Figure 2. Despite this, the metabolic network was believed sufficiently capable of describing the profiles of dFdC and its phosphorylated metabolites. The development of this model presented several difficulties. First, important information regarding the experimental setup were not present in the original publication on the in vitro experiments, like the milligrams of protein per number of cells and the area of the well in which the cells were cultured. Moreover, it was not possible to describe the profiles of dFdU and its metabolites. Given that the metabolic network structure depends on the data available, a kind of “structural uncertainty” of the model is present. This uncertainty could potentially impact the in vivo predictions of dFdCTP concentration once the network is included in the PBPK model.

To model dFdC distribution and metabolism in the body, a PBPK model was developed accounting for the activity of plasmatic membrane transporters.7, 32 One of the main advantages of the developed PBPK model is that the metabolic network was easily coupled with the model, leading to the possibility of describing the active metabolite concentration in the site of action. In order to do this, a compartment representing the pancreatic tumor was introduced into the model and was supposed to share the same extracellular environment with the pancreas sane tissue. By doing this, there was a direct correspondence between in vitro and in vivo intracellular and extracellular environments. In this context, information regarding the target enzymes concentration both in in vitro cancer cell lines and in vivo tumor samples, was found to be particularly useful for the parameters rescaling and, thus, the inclusion of the metabolic network in the PBPK model.

With the model presented here, there was no need of estimating the blood flow directed to the tumor. This approach presents, however, some drawbacks. In fact, by using this model, it would be difficult to describe processes such as the angiogenesis and the effect of a potential antiangiogenic compound. Moreover, the thick stroma surrounding the tumor cells that characterize the pancreatic cancer was not modeled46 and this could potentially impact the predicted drug disposition in the tumor tissue.

Another advantage of the developed systems pharmacology model is that it includes the interpatient variability of parameters, such as organ volumes, blood flows, and abundances of enzymes and transporters. Then, by performing a GSA, it is possible to understand what are the parameters that with their variation in the population mostly explain the interpatient variability of metrics of interest, such as AUC of plasma dFdC and tumor dFdCTP. The GSA results highlight that the tumor dFdCTP AUC variability is mainly explained by the variation of CDA and dCK tumor concentration. These results are in accordance with the fact that the dFdC clinical response is probably related to the patient’s genotype and to different expressions of the transporters or target enzymes.11

The results of this modeling study suggest that individual genetic factors affecting dFdC metabolism would lead to different amounts of its metabolites and, consequently, different treatment responses, as dFdCTP exposure has been previously related to tumor response, and later, to survival.22 A recent multiscale network characterizes the effect of proteomics on dFdC mechanism of action and its signaling pathways, in combination with birinapant.17 Future integration of their results with those present in this study could provide insights to better understand dFdC variability and drug effects.

Further research should be done for characterizing in vitro different pancreatic cancer cells coming from patients receiving dFdC, measuring the target enzyme level expression and the degree of their polymorphisms. This would be key to assess and refine the current model.

Funding

No funding was received for this work.

Conflict of Interest

The authors declared no competing interests for this work.

Author Contributions

M.G.C. and N.M. wrote the manuscript. All authors designed the research. M.G.C. and N.M. performed the research.

Supplementary Information 1.

Click here for additional data file.