Mechanistic characterization of oscillatory patterns in unperturbed tumor growth dynamics: the interplay between cancer cells and components of the tumor microenvironment.
Aymara Sancho-Araiz (1,2), Zinnia P Parra-Guillen (1,2), Jean Bragard (3,6), Sergio Ardanza (3,6), Victor Mangas-Sanjuan (4,5), Iñaki F. Troconiz (1,2,6)
(1) Pharmacometrics & Systems Pharmacology Group, Department of Pharmaceutical Technology and Chemistry, School of Pharmacy and Nutrition, University of Navarra, Pamplona, Spain. (2) IdiSNA, Navarra Institute for Health Research, Pamplona, Spain. (3) Department of Physics and Applied Math. University of Navarra, Pamplona, Spain (4) Department of Pharmacy and Pharmaceutical Technology and Parasitology, Faculty of Pharmacy, University of Valencia, 46100 Burjassot, Valencia, Spain. (5) Interuniversity Research Institute for Molecular Recognition and Technological Development, 46100 Burjassot, Valencia, Spain (6) Institute of Data Science and Artificial Intelligence, DATAI, University of Navarra, Pamplona, 31080, Spain
Objectives:
The adequate characterization of the tumor dynamics is crucial to understand cancer progression and response to treatment. In this sense, different empirical models are being used to describe unperturbed tumor growth (TG) in xenograft experimental settings [1]. However, in most cases, the tumor size/volume profiles show oscillatory trajectories which are ignored by the standard TG models and treated as random noise.
In this work, we investigate unperturbed non-monotonic growth profiles and hypothesized that they result from the interaction between the tumor and its microenvironment. We strive to simplify and as a main objective, develop a computational structure resembling the complexity of longitudinal tumor profiles and allowing the analyst to obtain information at the individual level.
Methods:
Tumor volume (TV) over time from control groups of different xenografts preclinical studies (85 mice including 10 cell lines and 6 tumor types) [2] were pooled together and analyzed with Monolix 2021R1.
First, an exhaustive statistical analysis was performed to discharge the possibility that the oscillatory patterns observed in the raw data were due to random variations (Kolmogorov-Smirnov (KS) test). Then, different models, ranging from empirical to more mechanistic, accounting for tumor dynamics were fitted to the data. However, all of the structures predicted a monotonic (smoothly continuous) increase in tumor volume over time.
Therefore, since none of the models provided a fair description of the data, we develop a novel TG model based on the prey-predator paradigm largely used in ecology [3] where the main assumptions are: (i) available resources (RES) such as nutrients or oxygen increase the tumor growth capacity (GC), promoting proliferation, and (ii) the increase in TV is accompanied with higher demand and consumption of RES, which eventually leads to a decreased GC that could translate into tumor shrinkage.
Model performance of the aforementioned model was compared with standard unperturbed TG models commonly used in population pharmacokinetic/pharmacodynamics analysis using the Akaike information criterion (AIC), and goodness of fit (GOF) plots including autocorrelation plots. Additionally, a local sensitivity analysis in which one parameter at a time was increased or decreased a 50% was performed to evaluate the impact on the predicted RES, GC and tumor size on day 18 after tumor inoculation.
Furthermore, a simulation study was undertaken to investigate the impact of virtual dosing schedules and combination therapies. Anticancer treatment was mimicked through a constant change in prey-predator parameters, one at a time and in multiple combinations, during 21 days (starting a day 7, 15, or 21 after tumor inoculation). Subsequently, the area under the tumor size vs time curve (AUTC) for each scenario was computed and compared to untreated AUTC.
Results:
The possibility that the oscillatory patterns were a consequence of random variations was statistically discarded (KS test p< 10-9). In addition, the statistical analysis also revealed that the characteristics of the oscillations were evenly distributed across the different cell lines without any pattern suggesting that a particular tumor type was more prone to show non-constant growth rates.
In the final structure, the RES-GC interaction is governed by: k1 and k3 (1st order rate constants) representing the growth rate of RES and the decrease in the GC, respectively; and k2 (2nd order) controlling the consumption rate of RES and the rise of GC. Additionally, this interplay impacts the rate at which TV grows controlled by kge (2nd order proliferation constant). A spontaneous cell death (kdeath, 1st order rate constant) is also assumed. Moreover, at the time of cell inoculation, TV0 and GC0 were estimated, whereas RES0 was fixed to 1. Note that the GC-RES system is uncoupled with tumor dynamics.
All model parameters were estimated with adequate precision (<30%), being the typical values of k1, k2 and k3, 0.061 day-1, 0.055 au-1×day-1 and 0.32 day-1, respectively. Moreover, the proposed model showed increased adequacy in all the numerical (ΔAIC > -87.19 and Δerror > -0.03 compared to the best standard model) and graphical diagnostics performed (p>0.05 in the lag plot and no tendency in the weighted residuals versus time plot). In this respect, the lag plots, scarcely used as GOF plots to compare among competitor models, emerge as a powerful diagnostic tool to support the selected model. In addition, the system showed to be more sensitive to the parameters controlling the oscillation behavior (i.e. k1 and k2). In this sense, higher values of k1 result in a higher number of tumor volume oscillations per unit of time (higher oscillatory periods), and the parameter k2, on the other hand, will mainly govern the amplitude of those oscillations.
The findings of the simulation studies showed that, besides chemotherapeutic agents increasing the killing rate of cancer cells, anti-angiogenic drugs able to decrease resources proliferation (k1) (able to reduce AUTC more than 50%), combined with immunotherapies that increase the degradation of the growing capacity (k3 or k2), have a significant impact on tumor shrinkage. Furthermore, the exploration allowed us to identify the optimal dosing schedule (simultaneous versus sequential administration) that could greatly increase efficacy in future combination approaches.
Conclusions:
In summary, the work presents a new semi-mechanistic model able to describe the oscillatory tumor size profiles over time in the absence of treatment administration, as a result of the interaction between cancer cells, the growing capacity of the tumor and the available resources in the tumor microenvironment. Although this work is not the first attempt of studying this complex pattern, to the best of our knowledge, it is the first data-driven modeling exercise that describes non-monotonic unperturbed TG through a simple system of ordinary differential equations, easily implemented in any modeling platform. The relevance of this approach lies in its potential translational impact at the time to design drug combinations strategies, emerging the possibility of locating anticancer drug effects on different parameter targets.
References:
[1] Benzekry S, Lamont C, Beheshti A, Tracz A, Ebos JMLL, Hlatky L, et al. Classical Mathematical Models for Description and Prediction of Experimental Tumor Growth. PLoS Comput Biol. Public Library of Science; 2014;10:e1003800.
[2] Parra-Guillen ZP, Mangas-Sanjuan V, Garcia-Cremades M, Troconiz IF, Mo G, Pitou C, et al. Systematic modeling and design evaluation of unperturbed tumor dynamics in xenografts. J Pharmacol Exp Ther. 2018;366:96–104
[3] Kareva I, Luddy KA, O’Farrelly C, Gatenby RA, Brown JS. Predator-Prey in Tumor-Immune Interactions: A Wrong Model or Just an Incomplete One? Front Immunol. Frontiers Media S.A.; 2021;12:3391