Modeling tumor size dynamics based on real-world electronic health records and image data in advanced melanoma patients receiving immunotherapy

Abstract

The development of immune checkpoint inhibitors (ICIs) has revolutionized cancer therapy but only a fraction of patients benefits from this therapy. Model-informed drug development can be used to assess prognostic and predictive clinical factors or biomarkers associated with treatment response. Most pharmacometric models have thus far been developed using data from randomized clinical trials, and further studies are needed to translate their findings into the real-world setting. We developed a tumor growth inhibition model based on real-world clinical and imaging data in a population of 91 advanced melanoma patients receiving ICIs (i.e., ipilimumab, nivolumab, and pembrolizumab). Drug effect was modeled as an ON/OFF treatment effect, with a tumor killing rate constant identical for the three drugs. Significant and clinically relevant covariate effects of albumin, neutrophil to lymphocyte ratio, and Eastern Cooperative Oncology Group (ECOG) performance status were identified on the baseline tumor volume parameter, as well as NRAS mutation on tumor growth rate constant using standard pharmacometric approaches. In a population subgroup (n = 38), we had the opportunity to conduct an exploratory analysis of image-based covariates (i.e., radiomics features), by combining machine learning and conventional pharmacometric covariate selection approaches. Overall, we demonstrated an innovative pipeline for longitudinal analyses of clinical and imaging RWD with a high-dimensional covariate selection method that enabled the identification of factors associated with tumor dynamics. This study also provides a proof of concept for using radiomics features as model covariates.

---

Study Highlights

• WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

Immune checkpoint inhibitors (ICIs) have become a cornerstone in melanoma treatment but only a fraction of patients benefit from this therapy. Population pharmacometrics models, which can quantify the contribution of specific factors to the efficacy of these therapies, have thus far been developed primarily using data from randomized clinical trials.

• WHAT QUESTION DID THIS STUDY ADDRESS?

How can real-world data (RWD) inform anticancer treatment decisions? Which clinically relevant factors are associated with tumor dynamics in patients with advanced melanoma receiving immune checkpoint inhibitors?

• WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

This study demonstrates that pharmacometric analyses of RWD can provide evidence of response to ICI therapy in the context of metastatic melanoma. Further, we show how machine learning techniques enable high-dimensional image-derived characteristics of each tumor lesion to be integrated with established pharmacometric approaches for covariate assessment.

• HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?

This study demonstrates the potential of high-dimensional analysis of clinical and image-based RWD to advance precision oncology toward more individualized treatment solutions.

---

INTRODUCTION

The development of immune checkpoint inhibitors (ICIs) has revolutionized cancer therapy. Ipilimumab, a CTLA-4 inhibitor, was approved as the first ICI for patients with advanced melanoma in 2011.¹ This was rapidly followed by the approval of antibodies targeting PD1 (e.g., nivolumab and pembrolizumab) or PDL1 (e.g., atezolizumab, durvalumab, and avelumab) in several indications. Multiple clinical trials are currently being conducted to evaluate the efficacy of ICIs as monotherapy or in combination with chemotherapy in various cancer types.² The efficacy of ICIs is particularly noteworthy in melanoma with a subset of patients achieving complete and durable response.³ However, only a subpopulation of patients benefits from ICIs and there is an increasing need for predictive biomarkers to guide ICI treatment individualization.

Modeling of tumor size dynamics was successfully reported in the literature to identify and quantify the sources of response variability. Tumor growth inhibition (TGI) models have been developed based on data from patients with melanoma receiving ICIs in the framework of randomized controlled trials (RCTs).^4,5 Results showed large interindividual variabilities (IIVs) in tumor dynamics, partly explained by disease characteristics or abnormal laboratory values. However, given the complex covariate structure, identification of prognostic and predictive markers of tumor dynamics requires more powerful algorithms, such as the ones developed for machine learning (ML) which are increasingly applied.^6–8 A recent analysis used such techniques to identify prognostic and predictive factors for tumor dynamics and overall survival in patients with gastric cancer receiving avelumab.⁹ In addition to standard clinical covariates, there is growing interest in medical imaging data to guide treatment algorithms via noninvasive image-based biomarkers.¹⁰ Quantitative image descriptors, so-called radiomics features, were shown to be predictive markers of overall survival in several cancer types^11–13 and have the potential to predict early response to ICIs in patients with melanoma.¹⁴ Whereas ML offers a powerful framework for dealing with large datasets, such as high-dimensional image-based covariates, no study has yet, to our knowledge, assessed their relationships to tumor dynamics in a TGI model.

In parallel, most prior TGI model developments have been conducted using data from RCTs and further studies are needed to translate these findings into the real-world setting. Although RCTs are generally considered to provide the highest level of evidence, retrospective analyses of real-world data (RWD) better reflect heterogeneous populations, diverse tumor phenotypes, and treatment combinations. Developing TGI models based on RWD, however, is challenging because the patients' overall tumor burden is not systematically quantified in clinical practice. The findings from radiological evaluations are mostly described qualitatively in natural language, providing insufficient details for quantitative assessment of tumor response. To overcome this limitation, we developed a semi-automated pipeline to collect quantitative retrospective radiological data on the evolution of metastatic cancers.¹⁵

In this work, we analyze RWD from [¹⁸F]Fluorodeoxyglucose (FDG) positron emission tomography/computed tomography (PET/CT) and contrast-enhanced CT images of patients with advanced melanoma receiving ipilimumab, nivolumab, pembrolizumab, or a combination of ipilimumab and nivolumab. We developed TGI models with the aim of evaluating the time course of tumor burden, quantifying the variability of response observed in a real-world setting, and identifying its sources. Besides accounting for clinical covariates, we also demonstrate how image-derived characteristics can be integrated in TGI models, and present an exploratory analysis on a subcohort of the study population that considers radiomic features as model covariates.

PATIENTS AND METHODS Data collection Cohort selection

Data aggregation and retrospective analysis for this study have been approved by the local ethics committee for patients who did not refuse general informed consent (CER-VD Protocol No. 2019–00448).

Analyses were performed using data from patients with advanced cutaneous melanoma receiving ICI (i.e., ipilimumab, nivolumab, or pembrolizumab) as a first treatment line at the Lausanne University Hospital. Cohort selection was conducted retrospectively using a pipeline for data extraction and curation published by Abler et al.¹⁵ Clinical data and treatment information (e.g., type of ICI, dose, and date of administration), as well as clinical covariates were extracted from electronic health records (EHRs). Patients were included in the analysis if their EHRs contained sufficient and reliable information about ICI administration and images (i.e., CT or PET/CT scan data). A thorough inspection of EHRs guided patient selection based on plausibility checks (e.g., a delay between two administrations longer than the approved interdose intervals could indicate that the patient received the treatment in an external hospital). In such cases, additional sections of the EHRs were explored to assess the plausibility of these findings and to potentially complete treatment information in our database.

Imaging data Collection

Routine clinical imaging data (FDG PET/CT and contrast-enhanced CT performed for oncologic follow-up) were retrieved for the selected patients. Both imaging modalities were included in the analysis, but for consistency, a patient was retained in the analysis only if all images were collected with the same modality. We considered baseline and follow-up images for the analysis if they were acquired within the 3 months preceding ICI initiation and until 1 year after ICI discontinuation or until the initiation of the next treatment line (if applicable), respectively. For PET/CT, only the patients with at least two images (regardless of the timing relative to treatment initiation) were included in the analysis, whereas for CT, those who had at least one baseline and one follow-up image were included. To prevent confounders from obscuring the real ICI effect, we excluded baseline images collected before cancer-related surgeries, as well as follow-up images collected after cancer-related surgeries or after the first radiotherapy session. Whole-body images acquired after local treatment of brain metastases were not affected by this selection because brain images are not optimally evaluable using PET/CT and therefore had been excluded from our analysis.

Lesion segmentation and volume quantification

All metastatic lesions were identified and segmented at each imaging timepoint in order to quantify their volume. Segmentation of the PET/CT scans relied on the PET-Assisted Reporting System prototype software (Siemens).^16,17 Data were further consolidated by focusing on lesions located in the patients' trunk and exploiting longitudinal imaging information.¹⁵ Differently, the 3D segmentation on CT scans was performed semi-automatically using the US Food and Drug Administration (FDA) approved Mint Medical (version 3.8.2; Mint Medical GmbH) by expert radiologists.¹⁸ For each patient and timepoint, the total tumor burden (TTB) was computed by summing the volumes of all individually segmented malignant lesions (including all target, non-target, and new lesions for CT scans).

Radiomics features as imaging covariates on a subcohort

Radiomics features provide a quantitative descriptor of a region of interest (ROI) in an image and thus allow for the noninvasive characterization of tumor lesions. A dataset of PET-based imaging features was constructed from the sets of radiomics features extracted from all identified tumor lesions. PET images at baseline and first follow-up were converted to standardized uptake value (SUV) maps and resampled to 1 × 1 × 1 mm³ isotropic resolution by linear interpolation before feature extraction. The ROI masks were resampled to the same resolution using nearest-neighbor interpolation. For each segmented malignant lesion, a set of 107 shape, SUV intensity, and SUV texture features, was computed from each scan's SUV map and each ROI using the Pyradiomics library¹⁹ with fixed bin width of 0.2 for intensity discretization. When multiple lesions were observed at the same imaging timepoint, their feature sets were combined into a single aggregated radiomics feature set: features corresponding to lesion volume and total metabolic activity were aggregated by summation; the values of all other features were aggregated by weighted mean using the lesion's proportion of TTB as a weighting factor. To capture changes in imaging appearance between baseline and follow-up examinations, a delta-radiomics feature set was computed by subtracting the values of the aggregated feature sets at baseline from those at the patients' first follow-up. The resulting radiomics feature set contained 216 features per patient: lesion count, delta lesion count, 107 aggregated radiomics features at baseline, and 107 radiomics delta-features.

TGI model development Base structural and statistical model development

A data-driven stepwise procedure was used to find the model that adequately fitted the data (Supplementary Material S1). Briefly, each model included at least a baseline tumor size parameter before treatment initiation (TTB₀ in cm³), a growth rate constant (k_growth in days⁻¹) and a killing rate constant (k_kill in days⁻¹). When baseline observation was missing, we kept the time zero observation with missing dependent value to further allow its empirical Bayes estimate (EBE). An exponential tumor growth with a fixed k_growth retrieved from literature in a population of melanoma patients receiving pembrolizumab⁴ was assumed (Supplementary Material S1). Models with a dedicated or shared k_kill for all drugs or for the two PD-1 inhibitors were compared. Drug effects were included in an additive form.

Different hypotheses on treatment dynamics were compared: ON/OFF treatment effect (using a 1/0 regressor in the dataset), treatment effect linearly dependent on drug dose, or kinetic-pharmacodynamic model²⁰ using linear relationship with dose. Images collected after treatment discontinuation (and prior to initiation of a subsequent treatment line) were included, and the effect of ICIs discontinuation assessed.

Parameter estimation and model selection

Tumor measurements were fitted using the Stochastic Approximation Expectation–Maximization algorithm implemented in Monolix. Model fitting quality was evaluated using diagnostic plots together with relative standard error (RSE, acceptable if <50%) of parameter estimates. The difference in objective function values (ΔOFV) of the reference and nested test models allowed their comparison, whereas Bayesian information criterion (BIC) was used between non-hierarchical models (Supplementary Material S1). Potential model misspecifications were assessed using prediction-corrected visual predictive checks (pcVPCs).

Identification of factors associated with tumor dynamics Clinical covariates

Once the best structural model was selected, a model including clinical covariates was developed on the entire dataset with the aim of assessing the impact of these factors on model parameters. The understanding of the mechanistic effect of covariates added to literature data from previous TGI models guided the preselection of potentially relevant clinical covariates. The following baseline characteristics were tested to explain variability of model parameters: gender, age, body mass index, sodium, erythrocytes, monocytes, leucocytes, neutrophil to lymphocyte ratio (NLR), platelets, creatinine (CRT), LDH, albumin (ALB), Eastern Cooperative Oncology Group performance status (ECOG-PS), prior adjuvant treatment (after surgery on the cutaneous primary tumor), BRAF mutation, and NRAS mutation (see details in Supplementary Material S1). Missing values were imputed to the median for continuous covariates or to the most represented category.

Preselection of radiomics features as covariate candidates

A separate analysis, including the ML-selected radiomics covariate candidates, was performed on the subset of patients for whom such data were available (thereafter referred to as “substudy model”). This model served as a proof of concept to demonstrate the inclusion of quantitative image-derived information as model covariates by using ML. It was compared to the base model run on this subset of patients (i.e., substudy base model).

Given the large number of radiomic features per patient, their direct assessment as model covariates using standard population pharmacometrics modeling approaches is not practically feasible due to redundancy, risk of collinearities, and excessive time/resource burdens. We used a multistep ML-based selection approach to identify a small subset of radiomic covariate candidates prior to evaluation in the tumor dynamics model. This selection approach seeks to identify features that are predictive for the EBEs of each of the model parameters. It consists of a global assessment step, followed by three feature selection steps, as illustrated in Figure 1 and described in Supplementary Material S1. All four steps (from 0 to 3) were repeated for the EBEs of each model parameter, resulting in different sets of radiomics features with variable size, typically between 10 and 20 features for each TGI model parameter.

FIGURE 1. Multistep machine-learning-based selection of radiomics features as covariate candidates. Before feature selection, we first established whether the radiomic dataset carries any relevant information for predicting each of the model parameters (step 0). If successful, we proceeded with feature selection for this model parameter: step 1 was designed to suppress redundant information by removing correlated features from the radiomics datasets. Step 2 compared the performance of multiple model-based feature selection approaches to identify the optimal supervised feature selection strategy. Step 3 used this optimal selection strategy to choose the final set of radiomics covariate candidates. Please see the Supplementary Material for further details about the individual selection steps.

Pharmacometric covariates assessment

In each model, covariates were first explored visually and then screened using the SAMBA-COSSAC algorithm implemented in Monolix (Stochastic Approximation for Model Building Algorithm – COnditional Sampling use for Stepwise Approach based on Correlation tests)^21,22 (Supplementary Material S1). This method uses samples from the a posteriori conditional distribution instead of EBEs, allowing an unbiased estimation of covariates’ effect even with sparse data prone to shrinkage. Forest plots showing the effect size of covariates (including uncertainty) were used to illustrate their influence on model parameters (Supplementary Material S1).

RESULTS Data

A total of 311 tumor size measurements (240 for PET/CT and 71 for CT) from 91 patients with melanoma were available for model building, with a median of three assessments per patient (range 2–10). For 27 patients (30%), the baseline tumor size was missing. Radiomic features at baseline and follow-up could be extracted for 38 patients from the PET/CT imaging cohort. Most of the images (n = 129, 41%) were collected on-treatment, whereas 108 of them (35%) were acquired after treatment discontinuation. Most of the patients (n = 38) received a combination of ipilimumab and nivolumab, and 18 of them had a follow-up nivolumab maintenance. Median doses administered were 225 mg (range 50–357 mg), 85 mg (range 56–300 mg), and 200 mg (range 128–256 mg) for ipilimumab, nivolumab, and pembrolizumab, respectively. Table 1 summarizes patient characteristics and Figure S2 shows tumor dynamics in our population.

TABLE 1 Characteristics of the study population.

Patients' characteristics (N = 91)	Median [IQR] or n (%)
Treatment received
Ipilimumab combined with nivolumab	20 (22)
Ipilimumab combined with nivolumab followed by a nivolumab maintenance	18 (20)
Ipilimumab alone	23 (25)
Nivolumab alone	3 (3)
Pembrolizumab	27 (30)
Demographics
Male sex	59 (65)
Age (years)	68 [60–76]
BMI (kg/m2)	26 [25–30]
Laboratory values
Sodium (mmol/L)	140.0 [138.0–145.0]
Missing values	7 (8)
Erythrocytes (T/L)	4.6 [4.2–4.9]
Missing values	8 (9)
Monocytes (G/L)	0.59 [0.45–0.73]
Missing values	8 (9)
Leucocytes (G/L)	7.0 [5.7–8.6]
Missing values	8 (9)
NLR	3.4 [2.4–4.3]
Missing values	8 (9)
Platelets (G/L)	248.0 [211.5–313.5]
Missing values	8 (9)
CRT (μmol/L)	80.0 [70.0–89.0]
Missing values	6 (7)
LDH (U/L)	209.0 [179.0–262.0]
Missing values	10 (11)
Albumin (g/L)	42.5 [40.0–45.0]
Missing values	27 (30)
Disease status
ECOG-PS
0	51 (56)
1	18 (20)
2	4 (4)
4	1 (1)
Missing values	17 (19)
Initial adjuvant treatment	13 (14)
Somatic mutation
BRAF
Non-mutated	53 (58)
V600E	26 (29)
V600K	6 (7)
G464R	1 (1)
G469R	1 (1)
L584F	1 (1)
L597Q	1 (1)
L597R	1 (1)
Unknown	1 (1)
NRAS
Non-mutated	56 (62)
Q61R	15 (17)
Q61K	13 (14)
Q61L	1 (1)
G12C	1 (1)
G12D	1 (1)
Mutated but mutation unknown	3 (3)
Unknown	1 (1)

Abbreviations: BMI, body mass index; CRT, creatinine; ECOG-PS, Eastern Cooperative Oncology Group Performance Status; IQR, interquartile range; LDH, lactate dehydrogenase; NLR, neutrophil to lymphocyte ratio.

Base model

It was not possible to estimate the growth rate constant, which was therefore fixed to a literature value obtained from a population similar to our patient's cohort.⁴ A TGI model with an ON/OFF treatment effect and a shared killing rate constant for the three drugs was the most parsimonious model providing the best description of the data. Similar or better performances were achieved with this model, as compared to more sophisticated models including additional killing rate constants and/or more complex treatment effects. A significantly lower OFV was observed in the model with an ON/OFF treatment effect compared to the one with a treatment effect linearly dependent on drug dose (∆OFV = −10.3, p < 0.05). The estimation of a single shared killing rate constant for the three drugs decreased BIC (∆BIC = −16.0 compared to a model with three killing rate constants and the associated IIV) and increased precision of parameter estimates. The same model parameters were used for modeling tumor dynamics during treatment and after ICIs discontinuation. Estimating separate residual errors for each image modality (i.e., CT or PET/CT) failed to improve the model. IIVs were estimated on TTB₀ and k_kill. Despite a statistically nonsignificant decrease in OFV when adding IIV on k_growth (∆OFV = −2.97, p = 0.09), it was retained in the base model because it makes sense that not all tumors grow at the same pace and to allow further covariate assessment on this parameter. Parameter estimates of the base model are reported in Table 2.

TABLE 2 Structural and random effects parameter estimates for the base model and for the model including clinical covariates (n = 91 patients).

Parameter	Base model		Model including clinical covariates
	Estimate	Relative standard error (%)	Estimate	Relative standard error (%)
TTB0 (cm3)	8.34	15.9	7.78	15.6
ωTTB0 (%)	200	10.8	123	13.0
β albumin			−5.27	27.1
β neutrophil-lymphocyte ratio			0.76	36.5
β ECOG >0			0.65	46.1
kgrowth (day−1)	0.0081	Fixed	0.0081	Fixed
ωkgrowth (%)	59	40.6	49	66.4
β NRAS mutated			−1.22	44.9
kkill (day−1)	0.01	18.9	0.008	22.1
ωkkill (%)	153	14.9	185	14.5
Additive residual error (log scale)	1.10	5.59	1.13	5.41

Abbreviations: β, covariate effect as defined in section “Pharmacometric covariates assessment”; ECOG, Eastern Cooperative Oncology Group; k_growth, tumor growth rate constant; k_kill, killing rate constant; TTB₀, total tumor burden at baseline; ω, interindividual variability reported as coefficient of variation (%).

Model equations:

OFF treatment: $\text{dTTB}/\mathrm{d}t=k_{\text{growth}}\times \mathrm{TTB}.$

ON treatment: $\text{dTTB}/\mathrm{d}t=k_{\text{growth}}\times \mathrm{TTB}-k_{\text{kill}}\times \mathrm{TTB}.$

Model including clinical covariates:

$\begin{array}{c}{\mathrm{TTB}}_{0,\mathrm{i}}=7.78\times {\left(\frac{\mathrm{ALB},\mathrm{i}}{42}\right)}^{-5,27}\times {\left(\frac{\mathrm{NLR},\mathrm{i}}{3.14}\right)}^{0.76}\times {\mathrm{e}}^{0.65\times {\text{ECOG}}_{{\mathrm{cat}}_{\mathrm{I}}}},\hfill \end{array}$ with ECOG_cat,_i = 1 when ECOG >0.

$k_{\text{growth},i}=0.0081\times {\mathrm{e}}^{-1.22\times {\text{NRAS}}_{{\mathrm{cat}}_{\mathrm{I}}}}$, with NRAS_cat,i = 1 for tumors bearing NRAS mutation.

Identification of factors associated with tumor dynamics Clinical covariates

The influence of several factors on model parameters was identified: ALB, NLR, CRT, and ECOG on TTB₀, and NRAS mutation on k_growth. The effect of CRT was not retained into the model given the observed narrow distribution of CRT levels limiting an appropriate assessment of the impact of this factor on the parameters. Table 2 reports the parameter estimates of the model including clinical covariates along with their precision. Figure 2 shows the forest plot for the effect size of covariates on TTB₀. Results highlight the clinically meaningful effect of ALB on TTB₀, with a median increase of 527% and a decrease of 47% in TTB₀ at low and high extremes of ALB levels, respectively. A smaller but still relevant effect was observed for NLR, with a median decrease of 50% in TTB₀ for the low extreme of NLR, and median 116% increase for the high extreme. TTB₀ in patients with ECOG greater than 0 was predicted to be 92% higher than in patients with ECOG equal to 0. Additionally, an effect of NRAS mutation was identified on k_growth, with tumors bearing NRAS mutation having a 70% lower k_growth compared to tumors without that mutation.

FIGURE 2. Forest plots of covariate effects on TTB0 for the model including clinical covariates. The gray area delimits the “no-effect” range (defined as ±20% of the point estimate of the covariate). The points represent the estimate of the parameter including the covariate effect and the whiskers show the 95% confidence interval. ECOG, Eastern Cooperative Oncology Group.

Radiomic features Covariate-candidates from ML screening

Sets of radiomic covariate candidates were identified for TTB₀ and k_growth. For k_kill, no indication for the existence of predictive information was found according to selection step 0 (Figure S3). Following the removal of highly correlated features (step 1), the optimal supervised feature selection strategy was determined by comparison of the average three-fold cross-validation performance across 10 repetitions with different cross-validation initialization seeds in selection step 2. Feature selection based on ridge-regression weights yielded the best performing feature sets for TTB₀ with average cross-validation R² (±SD) of 0.72 ± 0.05; permutation feature importance was identified as best feature selection strategy for k_growth with average cross-validation R² (±SD) of 0.36 ± 0.05.

The resulting sets of radiomics covariate candidates (step 3) contained 16 features for k_growth and 10 features for TTB₀ (Table S4).

Radiomic features inclusion in pharmacometric models

Despite limited data availability in this exploratory analysis, two radiomics features were significantly associated with TTB₀ (Gray Level Dependence Matrix [GLDM] dependence entropy, Gray Level Run Length Matrix [GLRLM] gray level nonuniformity), and three with k_growth (GLDM small dependence low gray level emphasis, GLDM delta dependence entropy, delta shape elongation) in the substudy model (Table S5). Most retained radiomic features were related to the texture of the lesions, characterized by GLDM. As examples, tumors with an elevated spatial heterogeneity in PET imaging appearance (reflected by the value of gray level nonuniformity) tend to have a higher TTB₀. In addition, a negative value of delta shape elongation was found to be significantly associated to high k_growth.

Model evaluation

A 25-point decrease in BIC was observed for the model including clinical covariates compared to the base model. Parameters of the model with clinical covariates were estimated with a good precision (RSE <46%), except for the IIV on k_growth (RSE = 66%; Table 2). The inclusion of clinical covariates explained 24% and 16% of IIV on TTB₀ and k_growth, respectively. The pcVPCs for this model demonstrated good agreement between observed data and model predictions, as shown in Figure 3.

FIGURE 3. The visual predictive checks for the model including clinical covariates. Shaded areas represent the 90% prediction intervals around the median (red) and the 10th and 90th percentiles (blue) based on 500 simulations. Blue solid lines show the median, 10th, and 90th percentiles of the observed data. Points represent individual observations.

A 54-point decrease in BIC was observed in the substudy model including radiomics features as compared to the substudy base model (Table S5). The inclusion of radiomics features explained 75% and 71% of the variability on TTB₀ and k_growth, respectively, whereas the variability on k_kill increased by 35%. This is probably attributed to the correlation between parameters which was not retained in the model. As expected, the limited number of patients in the substudy limited the accuracy of the parameter estimates (Table S5).

DISCUSSION

Application of population modeling methods to RWD to quantify disease trajectory, treatment effects, and associated sources of variability in clinical use settings is an emerging area of pharmacometrics. Although recent reports have described applications to therapeutic drug monitoring data,^23–25 real-world pharmacokinetic studies,²⁶ postmarketing surveillance,²⁷ and characterization of disease progression or drug effects in cardiovascular medicine,²⁸ women's health,^29,30 and nephrology,³¹ applications in oncology remain an untapped opportunity. The analysis of tumor dynamics in cancer patient populations based on real-world imaging data from EHRs, quantified trough lesions segmentation has not been described to our best knowledge. By applying a semi-automated analysis pipeline to routinely collected clinical tumor images, we developed a TGI model that successfully described tumor dynamics in patients with advanced melanoma. Modeling RWD is challenging due to the heterogeneous population and the sparse data available. The paucity of baseline images prevented us to estimate the growth rate constant that was fixed to a literature value obtained from a population similar to our cohort.⁴ Sensitivity analyses provided confidence regarding the suitability of the value for our analysis. Our data supported an ON/OFF treatment effect, with better model performances as compared to a model with a treatment effect linearly dependent on the doses. This is in line with literature evidence indicating that, at the dose ranges used in clinical practice, the exposure-response relationship of ICIs is at the plateau of the maximal response.^4,32

Additionally, we investigated different combinations of k_kill to deal with limited data in our population of patients receiving different ICIs. We first estimated a shared killing rate constant for nivolumab and pembrolizumab (i.e., PD-1 inhibitors) together with a dedicated killing rate constant for ipilimumab based on the respective molecular mechanism of action (anti-PD1 vs. anti-CTLA4). This choice was also based on literature data indicating that nivolumab and pembrolizumab may be interchangeable considering the available molecular, preclinical, and early clinical data.³³ The estimation of a shared killing rate constant for the three drugs was supported by our data with lower BIC and higher precision of parameter estimates. We also stratified pcVPCs by treatment without detecting any trend for worse prediction in a subset of patients, supporting the estimation of a shared k_kill for the three drugs. From a clinical point of view, although there is a lack of published comparisons between combination therapy of ipilimumab and nivolumab versus monotherapy with pembrolizumab for the treatment of melanoma, it has been demonstrated that both treatments provide comparable benefit in terms of clinical endpoint for other types of cancers, such as non-small cell lung cancer.^34,35

Several clinical covariates showed large effect size on tumor dynamics. Albumin was found as the most influential covariate on TTB₀, showing higher values in patients with low ALB levels. It has been reported that low ALB levels, reflecting cachexia, are associated with poor clinical outcomes³⁶ and thus likely present as higher disease burden in the analysis population. Additional covariates related to inflammation status and advanced disease were also identified as predictive factors of TTB₀, consistently with previously published studies.^4,37 Furthermore, our results indicated that tumors bearing NRAS mutation are associated with a lower growth rate. It has been reported that NRAS-mutated tumors are highly aggressive, nevertheless with a better response to ICIs.³⁸ Further studies with larger sample size are needed to disentangle the effects of NRAS mutation on k_growth and k_kill, and current finding has to be taken with caution. The addition of covariates explained 24% and 16% of IIV on TTB₀ and k_growth, respectively. Part of the unexplained variability could also be attributable to drug pharmacokinetics,⁵ molecular characteristics of the tumor (e.g., PD-L1 expression⁴), or histological features, such as lymphocyte infiltration status,^39,40 that were not collected in our study. Additionally, the number of lesions, as well as metrics related to their location and heterogeneity could improve TGI models.^41–43

Our proof of concept for selecting radiomic features as model covariates provides an example of “mechanistic learning”,⁸ where ML approaches are incorporated in a traditional pharmacometric workflow. There is a growing interest in ML as a powerful tool for screening of high-dimensional covariates.^6,7,44 Here, we successfully applied ML techniques to select covariate candidates from a high-dimensional dataset of image-derived radiomics features.

The field of radiomics builds on the premise that medical images contain information about the underlying pathophysiology, which, even if not visible to the human eye, can be captured via quantitative image analysis.⁴⁵ Reservations toward the use of radiomic prediction models in clinical decision making^10,46 are frequently attributed to the quality and reporting of radiomic studies, and the interpretability and reproducibility of radiomic prediction models.^47,48

In our study, we identified five radiomics features associated with estimated model parameters. However, as the feature values correspond to the weighted averages across a patient's tumor lesions at a given imaging timepoint or their change between timepoints (delta features), their physiological interpretation is challenging. Elevated spatial heterogeneity in PET imaging appearance (reflected for example by the value of gray level nonuniformity) was associated to higher TTB₀. This could be attributed to the fact that larger tumors may face challenges in maintaining homogeneous oxygen and nutrient supply, affecting the spatial pattern of cellular metabolic activity and thus FDG uptake. Indeed, it has been demonstrated that intratumoral FGD PET image heterogeneity is associated with regional tumor cellularity, hypoxia, necrosis, and angiogenesis.⁴⁹ In addition, a negative value of delta shape elongation, corresponding to a spherical lesion at baseline taking on a more “elongated” shape at the first follow-up, was found to be significantly associated to high k_growth. We hypothesize that tumors may have started elongating into the vessels and thus disseminating more easily. Overall, our results demonstrate that image-based covariates can effectively be integrated in conventional TGI models. Due to limitations of our dataset, radiomic information could only be leveraged for a subcohort of 38 patients. Further large-scale studies are needed to validate our methodology and assess the clinical relevance of radiomic features, to help advance the field towards clinical implementation.

In conclusion, our analyses established for the first time a robust workflow facilitating tumor growth inhibition model development using RWD quantified trough lesions segmentation. We demonstrated proof-of-principle in the context of therapy with immune checkpoint inhibitors in melanoma but this workflow is easily applicable to other cancer types and therapies. We identified potential factors associated with tumor dynamics and we furthermore demonstrated the feasibility of integrating high-dimensional radiomics features in a pharmacometric framework by using ML methods. Our real-world study requires further confirmation but demonstrates the potential for opportunities to advance precision oncology toward more individualized treatment solutions. Importantly, our example demonstrates the feasibility of extending the application of population modeling beyond clinical trial data to enable a quantitative analysis of longitudinal data mined from EHRs and clinical imaging to quantify real-world treatment effects and associated sources of variability. Future extensions of such research include establishment of quantitative linkage between real-world tumor dynamics and survival outcomes, as this could be valuable in furthering generation of real-world evidence of drug effectiveness in the post-approval setting and is relevant for setting benchmarks to guide clinical development of new therapies.

AUTHOR CONTRIBUTIONS

P.C., D.A., M.G., P.G., K.V., M.A.C., and N.T. wrote the manuscript. P.C., D.A., P.G., K.V., O.M., M.A.C., and N.T. designed the research. P.C., D.A., M.G., P.G., F.A., N.V.V., M.D., N.G., A.W., S.L., R.D.M., M.J., C.D., C.C., J.O.P., O.M., M.A.C., and N.T. performed the research. P.C., D.A., and M.A.C. analyzed the data.

ACKNOWLEDGMENTS

The authors thank the patients who signed the general consent for usage of their data and the medical teams who followed the patients. We are grateful to Marian Caikowski for his expertise with information systems at the CHUV Center for Precision Oncology. Access to clinical data was made possible thanks to the efforts of the CHUV Clinical Research Datawarehouse Team. We thank Mr. Frédéric Pedron of the CHUV Radiology Department for facilitating image transfers, and the Mint team for enabling batch exports of CT segmentation results. We thank Vijay Shah at Siemens for his help with PET/CT segmentation software. We thank Professor Adrien Depeursinge of HESSO Sierre for valuable input on radiomics feature computation. We thank Alain Munafo for his participation to the setup of this project.

FUNDING INFORMATION

This project was funded by the Merck Serono S.A, an affiliate of Merck KGaA, Darmstadt, Germany.

CONFLICT OF INTEREST STATEMENT

N.T. and P.G. are employees of Merck Institute of Pharmacometrics, Ares Trading S.A., Lausanne, Switzerland, an affiliate of Merck KGaA, Darmstadt, Germany. K.V. is an employee of EMD Serono. P.C. was an employee of Randstad and contributed as a paid contractor for the Merck Institute of Pharmacometrics, Ares Trading S.A., Lausanne, Switzerland at the time the paper was submitted. F.A. is an employee of the Swiss Data Science Center – EPFL, Lausanne, Switzerland, formally engaged in a research collaboration with the Merck Institute of Pharmacometrics, Ares Trading S.A., Lausanne, Switzerland at the time of the publishing of this study. All other authors declared no competing interests for this work.