Agent-based modeling in cancer biomedicine: applications and tools for calibration and validation

ABSTRACT

Computational models are not just appealing because they can simulate and predict the development of biological phenomena across multiple spatial and temporal scales, but also because they can integrate information from well-established in vitro and in vivo models and test new hypotheses in cancer biomedicine. Agent-based models and simulations are especially interesting candidates among computational modeling procedures in cancer research due to the capability to, for instance, recapitulate the dynamics of neoplasia and tumor – host interactions. Yet, the absence of methods to validate the consistency of the results across scales can hinder adoption by turning fine-tuned models into black boxes. This review compiles relevant literature that explores strategies to leverage high-fidelity simulations of multi-scale, or multi-level, cancer models with a focus on verification approached as simulation calibration. We consolidate our review with an outline of modern approaches for agent-based models’ validation and provide an ambitious outlook toward rigorous and reliable calibration.

KEYWORDS:

• Agent-based modeling
• multi-scale
• multi-level
• calibration
• validation
• optimization
• biomechanics
• biophysics
• cancer simulation
• precision oncology

1. Introduction

Cancer initiation and progression is often catalyzed by genomic instability and chronic inflammation (see Figure 1(a)). These processes set in motion a complex multistage trajectory characterized by a series of random events that confer neoplastic capabilities to normal cells. Through intricate interactions with neighboring cell populations, neoplastic cells orchestrate the formation of tissue within its own microenvironment and the remodeling of the microvasculature and the lymphatic system. Notably, this microenvironment evolves dynamically, influenced by the active contributions of unaware normal cell populations.¹

Figure 1. (a) Diseases pathogenesis. Environmental factors, ionising radiation, and an unhealthy lifestyle (e.g. smoke, obesity) are often regarded as triggers for the pathogenesis of, among others, cardiovascular and neurodegenerative diseases, cancer and diabetes. Together with chronic inflammation, where inflammatory cells (macrophages, T lymphocytes and neutrophils) secrete large amounts of reactive/nitrogen oxygen species and cytokines, these factors can impair the normal functioning of cells. As a consequence, apoptosis and necrosis of cardiomyocytes, beta cells in the islets of Langerhans and neurons can lead to cardiomyopathy, diabetes, and Alzheimer’s disease, respectively. The aforementioned risk factors can also cause genetic mutations and chromosome instability in somatic cells which can, in turn, switch to a hyperplastic, malignant phenotype and become cancerous. As these cells proliferate and consume nutrients, carcinogenesis is initiated, which stimulates angiogenesis and tissue remodelling. Finally, if circulating tumour cells extravasate the blood vessels’ epithelium, other organs can be invaded. [ROS: Reactive Oxygen Species, NOS: Nitrogen Oxygen Species] (b) Calibration process formalism of multi-scale biophysical processes: from cell-level homeostatic dynamics, to competition and cooperation at population level, and up to organ-level phenotypic interactions.

Understanding the underlying mechanisms of cancer and devising effective therapeutic strategies are formidable challenges, further exacerbated by the disease’s complexity and heterogeneity. With a myriad of cancer types exhibiting distinct molecular profiles and clinical behaviors, the quest for universally effective treatments remains elusive. Moreover, the significance of patient-specific parameters, including genetic variations, Tumor Microenvironment (TME) heterogeneity characteristics, and individual responses to treatment, further complicates the therapeutic landscape. Developing and testing new treatments in the clinical setting is a protracted endeavor, often spanning years of rigorous experimentation and clinical trials. However, mathematical and computational models – most commonly referred as in silico models – offer a promising avenue to expedite drug discovery and development process.² By leveraging computer simulations, investigators can systematically probe the mechanisms of cancer mechano-bio-physics, explore in silico new hypothesis for cancer biology while also develop novel procedures for disease stratification and prognosis, contribute toward reducing animal testing and reduce the costs associated with clinical trials, and promote personalized therapy of cancer.

In cancer biomedicine, mathematical models based on Ordinary Differential Equations (ODEs) and Partial Differential Equations (PDEs) formulations serve as indispensable tools for simulating the dynamics of neoplasia across different scales. A comprehensive overview is given in the excellent reviews by Byrne et al.,³ Yin et al.⁴ and Hadjicharalambous et al.⁵ These modeling approaches offer efficient means to capture the complex interactions within the tumor – host ecosystem. ODEs are ideal mathematically simplifying the dynamics of large populations of cells (e.g., the growth dynamics of tumors, for instance, can be realistically approximated by a Gompertz function) by averaging parameters of cell clusters. Compartmental models are a class of ODE-based models – they are commonly used in pharmacokinetic – pharmacodynamic analysis of drugs – that are capable of representing phenomena evolving at various scales and driven by different kinetics. They can capture dynamic fluctuations between different cell populations within the TME, as well as, for instance, changes in the occupancy of cellular receptors or transitions between transcription factors. With regard to the analysis of drugs, such models are being used to investigate a drug’s path through the body, and its resulting concentration in different body compartments, or with the physiological and behavioral consequences produced by a drug molecules, e.g., in the cell viability or its phenotypic response. On the other hand, PDE-based models are well-established in mathematically replicating growth phenomena in developmental biology as well as in cancer growth, the mass balance of cell populations (including transport of cells as in tumor metastasis and in cell differentiation) and in the delivery of anti-cancer drugs, e.g.^6–10 These processes are typically modeled using Turing’s reaction – diffusion model in two or multiple dimensions in space – time. By considering both spatial and temporal changes, PDEs offer pointwise information about the distribution and dynamics of substances within the TME. For instance, regions characterized by necrotic cells and hypoxia can be accurately modeled, providing insights into the spatial heterogeneity of tumor growth and response to therapy. Together, mathematical models governed by differential equations of cancer mechano-bio-physics can provide complementary approaches for understanding cancer dynamics and can offer insights into the spatial and temporal complexities inherent in tumor systems. Furthermore, these systems can often be efficiently implemented as matrices and solved rapidly, benefiting from many fast matrix solving algorithms, thus reducing the need for extensive computational resources. However, their ‘modelling resolution’ is limited by assumptions intrinsic to their continuous nature – they can mathematically describe bio-physical phenomena via averaging of the dynamical changes occurring in cell populations or the balance of biochemical cues, etc. This offers therefore room for the development of particle-based models to simulate biological cells (or other pertinent entities of particulate matter) as discrete entities.

Agent-Based Modeling (ABM) is a powerful simulation modeling technique that has seen increasing research and development in biology and in biomedical applications.^11,12 ABM represents processes and phenomena in terms of computational agents. Agents can be denoted as particles to represent cells in biology, or segments to represent neurites or vessels, which can reside in space either in a structured lattice or following an off-lattice approach. In ABM, agents (cells) are programmed with respect to their behavior and interaction with other agents, modeled as Markov processes (i.e., such that the transition to a new state depends solely on the previous one) using mathematical rules for their description. When the agents’ decisions arise from probabilistic reasoning, stochastic systems can be simulated and complex, higher-scale behaviors emerge as the simulation clock ticks. The agents’ rules can, under suitable conditions, enable multi-cellular systems to self-organize into highly nonrandom structures. Thus, ABM simulations in cell biology are characterized by the dynamics of autonomous and heterogeneous entities whose phenotypic behavior is explicitly modeled (see Figure 2(b)), yet cell behavior can adapt in time based on cell-to-cell and cell-to-microenvironment interactions (see Figure 2(c)).

Figure 2. Approaches to calibration and validation of agent-based models. (a) Prototypical representation of the calibration and validation process and its connection to the agent-based model. High-level methodological approaches for ABM calibration and validation include optimization-based, Machine Learning (ML)-based, and hybrid modeling approaches. These techniques facilitate the translation of performance measures into parameter configurations of an agent-based model. (b) Prototypical representation of an agent’s interaction with the surrounding environment in an agent-based model. The agent perceives the environment and processes the gathered information. The agent’s internal rules, defining its behaviors, can update its internal state and trigger actions on the environment, resulting in local changes. (c) A common (not exhaustive) set of cell behavior mechanisms includes: ‘Divide,’ ‘Migrate,’ ‘Change morphology’ (grow, shrink, polarize), ‘Die’ (programmable, stimulated), ‘Differentiate’ from one cellular phenotype (red) to another (purple), ‘React to a biochemical cue,’ ‘Secrete a substance,’ and ‘Interact’ with other cells or the extracellular matrix. (d) An ABM software can simulate multicellular populations to resolve the interactions amongst different cell phenotypes as well as the interactions of cells with blood vessels, the lymphatics, and the immune system. In the above schematic, biochemical interactions are indicated as a cloud of particles or protrusions to indicate lipids/proteins, while biomechanical interactions are depicted with springs.

In the context of cancer biomedicine, ABM offers a unique capability to simulate single-cell behaviors and interactions. Figure 3 illustrates several examples of agent-based models of two- or three-dimensional mono- or multicultures of cells that simulate the growth of tumor spheroids, the development of the brain cortex, the lung airways at the alveolus-level, or hybrid models that couple ABM with Finite Element Method (FEM) to simulate glioma development. Given the pivotal role of single-cell interactions in driving tumor progression and heterogeneity, agent-based models provide a powerful framework for predicting the evolution of cancer. Particularly significant is ABM’s ability to model phenotypic changes and mechanical interactions, mirroring the complex dynamics observed in vivo. Furthermore, advancements in imaging techniques enable the tracking of single cells, facilitating the conversion of imaging data into mechanistic rules for ABM simulations. This integration of experimental data with computational models not only enhances our understanding of cancer biology but also enables the design of more targeted and informative experiments. By leveraging ABM, researchers can generate hypotheses, test various scenarios, and refine experimental designs in silico before conducting costly or time-consuming in vitro or in vivo studies (see Figure 2(a)). Thus, ABM serves as an invaluable tool for bridging the gap between experimental observations and theoretical predictions in cancer research, ultimately driving the development of more effective therapeutic strategies.

Figure 3. Collage of still images from agent-based model simulation examples,^13–17 depicting either off-lattice ABM or cellular automata approaches, which demonstrate the applicability of the method to simulate mono- and multicultures of cells in two and three dimensions, across different spatial scales as well as being coupled with other numerical methods (e.g., Finite Difference Method, FEM).

In this review, we begin by surveying a comprehensive yet non-exhaustive list of agent-based models addressing various aspects of cancer biomedicine. Additionally, we offer a comparison with existing approaches and explore hybrid multi-scale adaptations. Finally, we delve into top-down and bottom-up strategies for calibrating and validating agent-based models (see Figure 1(b)).

2. Agent-based modeling in cancer biomedicine

2.1. ABM as a mechanistic approach (2010–present)

Intrinsically decentralized, inherently interactive, and multi-entity, biological systems are ideal to simulate using ABM (see Figures 2(c,d)). Among other diseases, cancer biology has witnessed an increasing number of modeling efforts via ABM over the years; ABM, alone or coupled with other in silico modeling techniques, has been employed to explore the dynamics of the TME, to probe and identify new therapeutic agents, to promote clinical translation by aiding the development of new diagnostic tools, and to bridge the gap between animal and human data.^5,18,19 The short survey of Metzcar et al.² presents the state-of-the-art ABM simulations related to cancer hypoxia and necrosis, tumor-induced angiogenesis (i.e., the sprouting of new vessels from the existing ones triggered by tumor signals), invasion, stem cell dynamics and immuno-surveillance in neoplasia during the 2000s and the first half of the previous decade. In this first section, a survey of the major agent-based approaches of the last decade to model cancer development and potential therapeutic strategies is presented (see Table 1).

Table 1. Overview of the relevant literature from the last decade in ABM procedures in cancer biomedicine.

Model focus	Main topics/findings	References
Carcinogenesis	Inflammation	^Citation20
	Genomic instability and karyotype	^Citation21,Citation22
Tumor growth	Avascular/Vascular	^Citation14,Citation23–25,Citation28,Citation37,Citation60
	Role of hypoxia	^Citation45,Citation46,Citation61
	Cell heterogeneity and phenotypic switches	^Citation23,Citation29,Citation30,Citation32,Citation34
Spatial characterization	ECM and cell migration	^Citation43,Citation55
	Proliferation vs migration	^Citation31,Citation32,Citation38,Citation39,Citation44
	Shape metrics for tumor morphology	^Citation39,Citation40
	Tumor cells’ flows	^Citation41,Citation42
Somatic cells interactions	Immune cells	^Citation44,Citation47,Citation56
	Connective tissue	^Citation55,Citation57
Metastatic processes	Main contributing factors	^Citation8,Citation36,Citation39
Treatment strategies	Interrupt cell-cell communication	^Citation25,Citation26,Citation54–56
	Reactivation of senescence	^Citation35
	Combined treatments	^Citation16,Citation27,Citation58,Citation59,Citation62
	Oxygenation level	^Citation33,Citation45,Citation46
	Morphology and packing density	^Citation40,Citation41
	Immunotherapy	^Citation44,Citation47–52,Citation55
	Oncolytic viruses	^Citation53

2.1.1. Pathogenesis

Despite being just nodes in a much broader network, chronic inflammation, unhealthy lifestyle, and environmental factors have long been regarded as triggers in the pathogenesis of cancer (see Figure 1(a)). An et al.²⁰ introduced an inflammation and cancer development ABM platform. Built on a hierarchy of relationships between cancer hallmarks, in the model higher-order processes are driven by lower-order ones. The basal DNA damage rate of healthy cells was exacerbated by reactive oxygen species secreted by inflammatory cells (e.g., neutrophils, monocytes, and macrophages). If the DNA-repair rate was exceeded, an abstract genome (i.e., a simplified genetic profile where different gene locations are associated with different cellular pathways) was impacted and cells’ functions were impaired. Damaged cells could also recruit more inflammatory cells, hence, leading to the establishment of a positive feedback loop that highlights the potential role of anti-inflammatory therapies in cancer care. Thus, in the work by An et al.²⁰ cancer was presented as an evolutionary process, and an interesting analogy between oncogenesis and evolution was highlighted, both requiring a basal level of genetic instability. Contextually, Araujo²¹ and Lynch²² and their colleagues emphasized on the role of genomic instability in carcinogenesis. In the first model, the impact of each gene was proportional to the number of copies in the genome, where abstracted genes (i.e., proxies for realistic gene networks that regulate a single, specific cellular function) were used for division, death, and segregation mechanisms. When in homeostasis, the effects of proto-oncogenes (linked to the cell growth) and tumor suppressor genes were balanced. However, cell duplication might randomly lead to aneuploidy in daughter cells causing hyperplastic growth. Different treatments were simulated, which eventually led to the conclusion that the best prognosis arises from a combination of chemotherapy and surgery. Besides, the simulations generated novel genotypes that resemble patterns found in cancer patients. Similarly, the model presented in the work by Lynch et al.²² linked chromosomal instability with the cells’ karyotype, i.e., the whole set of cellular chromosomes. A framework to measure the chromosomal instability level was built and three models for the selection of a cell after the division were compared. Following the division, the two cells underwent a selection procedure and, if not fit, they were removed from the agent-based simulation. The total fitness is then computed as a sum of the single chromosomes’ contributions. Phylogenetic trees were used in combination with approximate Bayesian computation (a computational method based on the Bayes theorem that can generate posterior probability density functions for the models’ parameters) to estimate the probability of missegregation from an observed population of cells. The results showed that sampling karyotypes in a cell population doesn’t allow direct determination of the chromosomal instability, as other factors play key roles. However, selection and missegregation shape the karyotype diversity in a population of aneuploid cells.

2.1.2. Neoplastic cell dynamics

Cells in neoplasia can undergo phenotypic switching multiple times during their lifetime. This process is typically driven by endogenous or/exogenous factors, such as the lack of nutrients or mechanical stimuli. The ability to change their behavior and adapt to the surroundings can be effectively modeled by agent-based models, as demonstrated for example in the paper of Chen et al.²³ They simulated the avascular growth of an in vitro tumor spheroid via a two-dimensional (2D) agent-based model, where nutrients are supplied solely by the existing environment. Interestingly, some modeling concepts were borrowed from the evolutionary game theory in that phenotype switching (mostly proliferation – migration and vice versa) is regulated via a payoff matrix. The phenotype-to-phenotype competition was modeled as a rewarded game where the environment (i.e., the available resources) influences the cells’ phenotype, the tumor rate of growth and the surface roughness. The latter, in particular, can be quantified. In the papers of Kareva et al.²⁴ and Phillips et al.,²⁵ vascular tumor growth was simulated using ABM driven by tumor hypoxia and the subsequent secretion of pro-angiogenic cytokines (e.g., vascular endothelial growth factor). The models showed that the rate of diffusion and consumption of the growth factors is key in the angiogenesis. The competition game is reiterated, although in this case it concerns high and low-affinity angiogenesis regulators. Reducing the tumor – endothelium communication would therefore allow control over the angiogenesis mechanisms as shown, for instance in Wang et al.²⁶ There, the Loewe combination index²⁷ (which quantifies the extent of drug synergy in combination therapies) revealed that merging cell-killing drugs and vascular endothelial growth factor inhibitors provides the best treatment for melanoma. These mechanisms were further investigated in the work of Lima and his colleagues,²⁸ where a novel coarse-grained modeling approach was introduced. Multiple cells with the same phenotype were modeled using a single agent and the total computational time dropped by 93–97% with little difference with respect to the cell agents model. A multi-scenario (i.e., with different initial conditions) calibration of the agent-based model was performed via moment-based Bayesian inference. Finally, a time-dependent global sensitivity analysis allows the authors to identify the increase in the death rate due to lack of glucose as the most influential parameter. Overall, all spheroid model simulations showed that, eventually, the core of the tumor was dominated by necrotic cells, while the ones alive were located at the rim of the tumor – typically referred as the tumor enhancing region.

Predicting cell heterogeneity and understanding the relative impact of intrinsic versus environmental factors on its emergence is of paramount importance; this is often seen as the most influential factor in multi-drug resistance, treatment failure and relapse. Although hardly achievable with in vivo models, in silico models provide frameworks to quantitatively measure these relative weights. Gallaher et al.²⁹ developed a model of intra-tumor heterogeneity in Glioblastoma (GBM) and used the in silico model to show that some level of intrinsic heterogeneity is required to capture the migration behavior observed in single-cell data. Interestingly, the environmental heterogeneity alone proved to be insufficient. In their model, GBM growth was driven by the platelet-derived growth factor that acted both via paracrine and autocrine signaling (i.e., by interacting with the surrounding cells and regulating the secreting cells themselves, respectively). Combining treatments that inhibit cancer cell proliferation and suppress cell invasiveness (i.e., migration), led to an increased efficacy with respect to monotherapies. Moreover, the model not only proved the predictive capacity of single-cell data in silico, but it also emphasized its importance by showing that cell populations with heterogeneous phenotypes displayed similar growth dynamics and final density distributions. Interestingly, Greene and his colleagues³⁰ further explored the role of cell heterogeneity on tumor growth using a continuous-time Markov chain model (i.e. a stochastic process where the transition to a new state depends only on the previous one). The model described cell states’ transitions and parameter estimation was performed using constrained trust region algorithms (i.e., optimization methods focused on a small region of the parameter space and constrained by, for example, biological observations) with nonlinear least squares. The phenotypic trade-off was similarly investigated by Gerlee et al.,³¹ which focused on GBM, and later on by Gallaher et al.³² Here, the migration (favored by selection at the early stages of tumor development) and the proliferation (favored more in the later phases) capabilities improved simultaneously up to a certain coordinate of the trait space. From there, the increase in one of the rates could only come at the expense of the other. Different shapes of the trait space were analyzed and an inverse proportionality between the rate of cell turnover and phenotypic variability among cancer cells was found. Modeling of GBM has also been the focus of,³³ where special emphasis was put on the interplay between hypoxia and cancer progression. The results of their study illustrated differences in terms of the spatial distribution of oxygen/nutrients within the TME. This, in turn, can affect not only the rate of growth of the carcinoma but also the migrational capacity of the glioma cells. Notably, the authors draw attention toward hypoxia as a catalyst for dangerous mutations at a higher rate. In an effort to model the phenotypic transitions of cancer cells, Axenie and Kurz³⁴ proposed a model that learnt the mechanistic rules governing cancer’s cells phenotypic staging from quiescent to proliferation and from proliferation to apoptosis. Using a typical cancer phenotypic state space, quiescent cancer cells (Q) can become proliferative (P) or apoptotic (A). Non-necrotic cells become hypoxic when oxygen drops below a threshold value, while hypoxic cells can recover to their previous state or become necrotic. The transitions among these states are stochastic events generated by Poisson processes. Although trained on a limited time series of raw immunohistochemistry and morphometric data from 17 patients, the lightweight machine learning system was also able to accurately predict: tumor volume evolution, the dependency between histopathological and morphological data (such as nutrient diffusion penetration length within the breast tissue), the ratio of cell apoptosis to proliferation rates and the radius of the tumor. More recently, Gazeli and his colleagues¹⁶ employed ABM to simulate in vitro experiments on melanoma (B16F10) cancer cells monolayer growth when treated with Doxorubicin alone, or in combination with a treatment based on cold atmospheric pressure plasma jet. Their model was designed to probe the mechanisms of action of each therapeutic approach (cytotoxic drug or/with plasma); this was characterized through model-derived probabilities of the melanoma cell apoptosis and division. They presented an interesting approach that combines in silico with in vitro, and they demonstrated how simulations can help to speed up laboratory work and, thus, reduce the costs for cancer drug/treatment testing.

An interesting finding is reported in Poleszczuk et al.,³⁵ where the authors observed a reduction in the tumor volume by exploiting the induction of cancer cells’ senescence. They presented an ABM procedure to probe the competition of cell subpopulations. Cancer stem cells were allowed to “fight” against cancer progenitor cells to ensure vital space and nutrient supply. Their simulations showed that the tumor growth is regulated in two distinctive phases, in which an initial increase of the cancer cell population was followed by a reduction in the proliferation rate and eventually tumor control. The latter was in fact caused by the reactivation of the senescence program in the progenitor cells that constrained the stem cells in the tumor core, thereby limiting their proliferation. Similar in silico results concerning neoplasia spatial inhibition were presented by Norton et al.³⁶; their focus was to investigate the role of cancer cell seeding in metastatic tumor progression. In particular, the authors modeled the metastatic cells dissemination in two potential scenarios: the ‘site’ seeding, where cancer-stem cells were injected from a single direction due to a breach in the vasculature, and the ‘volume’ seeding, where seeding was allowed to happen at random locations in the metastatic tumor region. While migration promoted tumor growth in every scenario, even when the cell division rate was considered high, volume seeding enhanced tumor growth. However, the impact of the seeding procedure was reported higher when cells’ quiescence inhibited spatial growth. Moreover, by extending the simulator into a three-dimensional (3D) ABM simulator, the authors were able to recapitulate the visual differences in the tumor morphology with respect to the proliferation and migration parameters of their stochastic model. Recently, de Montigny and his colleagues¹⁴ proposed a hybrid approach that uniquely combines agent-based and finite element modeling to simulate GBM progression and bridge the gap between continuum-based and discrete system dynamics. While the transport of nutrients and intra-cellular signaling were simulated using the FEM, the cells were modeled as agents, with volume averaging used to interface the two spatial scales (tissue and cell scale). Their methodology helped reduce by several orders of magnitude the number of simulated agents and, consequently, bring down the total simulation time. Notwithstanding this, the hybrid model can replicate growth patterns of both low- and high-grade tumors and assesses the role of platelet-derived growth factor on the tumor shape and size at later time points. A somewhat similar methodology was presented by Rahman et al.,³⁷ where GBM growth was modeled at multiple spatial and temporal scales, ranging from sub-cellular signaling pathways to the progression of the tumor tissue. By combining PDE solvers, an ABM simulator and ODE solvers for the tissue, cellular and intra-cellular scales, respectively, the authors provided a coupled cancer model where inter-compartmental communication ensures synchronization. The model was tested to replicate experimental findings concerning tumor growth and cell proliferation both in physiological conditions and following the administration of biochemical compounds (i.e. rapamycin).

2.1.3. Spatial characterisation

Tumor diagnosis together with the best treatment strategies can benefit from the characterization of the shape and spatial features of a solid tumor. Structural imaging, such as Magnetic Resonance Imaging (MRI) and histology images, and in silico modeling, such as ABM simulations, can provide spatial information and make predictions of the geometric features of carcinoma as they develop over time. Focusing on ABM (see Table 1), local modifications in the TME can be linked to outcomes at the tumor scale, overcoming the limitations imposed by reaction-diffusion equations that take average parameters from MRI or Computed-Tomography scans. As argued by Klank and her colleagues,³⁸ modeling GBM growth as a Brownian motion via mechanistic rules of an agent-based formulation is ample to replicate the development of a highly packed tumor core with an enhancing (proliferative) boundary at the tumor – host interface. In fact, the overall expansion speed of the lesion depends on single-cell migration rates. In Norton³⁹ and Karolak et al.,⁴⁰ shape metrics such as the mean chord length, the moment of inertia, the radius of gyration and the accessible surface area are used to characterize tumor morphology and packing density. The results show that this data can be extracted from diagnostic images and allow for tumor invasiveness and cancerous nature predictions. And while the more generic and commonly used tumor diameter doesn’t provide information about the tumor architecture, the aforementioned metrics could be adopted to supplement this limitation. In fact, at similar tumor sizes substantial morphological differences may be concealed. In this regard, computational models could be used to map anatomical compositions (in terms of shape metrics) to the corresponding effective drug penetration rates. These, in turn, could be employed to adapt the therapeutic protocols (e.g., drug doses and schedule) prior to treatment. Additionally, Norton et al.³⁹ reported that tumor morphology and invasiveness are directly linked to cancer cell phenotypic ratios, the level of hypoxia and the number of chemokine receptors. Indeed, reduced tumor growth was observed following the virtual administration of a drug that impaired cancer cell proliferation. Similarly, Bull et al.⁴¹ modeled via ABM the advective flow of microspheres from the tumor rim to the tumor core that results from the outer pressure in tumor spheroids. Their simulations illustrated distinct spherical shells in which the cells’ movement was either dominated by Brownian motion or advection. More specifically, the diffusion of the micro-beads located in the tumor rim was Brownian dominated, as the parent proliferating cells, located in this outer shell, placed their daughter cells randomly upon division. On the contrary, dying cells in the necrotic core left empty spaces. This, in turn, led to a depression that was counterbalanced by the cells in the outer shell, resulting in an advective motion. The authors introduced novel metric parameters (such as the waiting time of cells in the proliferating rim or the radial infiltration velocity in the shell between the rim and the necrotic core) whose values were mapped to the composition of the simulated spheroids. This, in turn, provided a new way to infer the underlying morphology (e.g. the quiescent area) from measurements of the microbead’s trajectories. Another approach to modeling tumor cells’ flows was presented by Jamous et al.,⁴² which simulated oncostreams (i.e., cells migration in opposite directions) and flocks (i.e., cells migration in the same direction) in 2D and 3D. They reported that the presence of oncostreams correlates with tumor progression, while they also interrogated in silico the parameter space impacting the mode of cancer cell migration. The authors showed that the formation of flocks in 2D simulations augmented as the cells’ shape was shifted from round to ellipsoid. Using the total polarization of the configuration as a proxy, the higher steering capability of the cells (which is correlated with the eccentricity of their shape) was found to be the reason behind the increased flock formation. As the simulation domain was extended from 2D into 3D, cells were provided with an additional degree of freedom. While both streams and flocks emerged at low cell numbers, only the streams were observed at higher densities. Moreover, the authors found that the tumor dynamics is strongly affected by the cell density and that both flocks and streams emerge when the ability of the cells to steer drops; thus, dismantling of oncostreams was proposed by Jamous et al.⁴² as a new therapeutic approach. Another agent-based model of tumor cell movement was presented in the work of Suveges and colleagues,⁴³ which emphasized the role of the Extracellular Matrix (ECM). More specifically, the authors developed a hybrid multi-scale model to investigate if and how the ECM could impact the cell invasion patterns of cancer cells. Cells were modeled using an agent-based model which was linked to, affected and was affected by a continuous model of the ECM. To simulate the adhesive interactions between the cells and the ECM, which don’t occur at a single point of contact, the authors employed non-local adhesion integrals. These allowed long-distance interactions to be taken into account by defining a sensing region over which the adhesion strengths were integrated. Their model demonstrated that aligned ECM fibers are necessary for tumor aggregations to move, while tumor invasion is impaired when the matrix fibers are aligned in parallel to the tumor margin. Importantly, their in silico findings were confirmed with experimental results, and they argued that the tumor expansion speed could be predicted from the orientation of the ECM fibers.

The impact of neoplastic cell heterogeneity on the tumor shape was further investigated by Gong et al.,⁴⁴ where a spatially resolved agent-based model was combined with a Quantitative Systems Pharmacology model (QSP). The QSP, an immuno-oncology mathematical model of ODEs, was used to simulate interactions among multiple compartments at the tissue scale. Interestingly, the propagating front of the tumor-enhancing region and the tumor necrotic core were simulated by two different agent-based models. The effect of immune checkpoint inhibitors (such as antiPD-1) on the tumor growth was also simulated. Virtual patients and clinical trials were probed using sets of parameters generated via Latin Hypercube Sampling (LHS), and their agent-based model aided in identifying predictive biomarkers for the tumor diameter, anti-PD-1 responsiveness, and time to cancer progression. Spatial heterogeneities within the tumor volume are not limited to cell morphologies, but rather encompass substance concentrations (e.g., glucose and oxygen). Modeling of these transients is well suited to agent-based models and allows, among other uses, inspecting the effect of local concentration changes on the cell cycle. Representative works that focused on the latter aspect were those of Hong et al.⁴⁵ and Kempf et al.⁴⁶ These models provided insights into the role of hypoxic conditions in cancer treatments, assessing the capability of hypoxia-activated pro-drugs in killing heterogeneous bystander cells otherwise unreachable. Moreover, they highlighted the importance of timing and hypoxic sensitizers to maximize the efficacy of radiotherapy.

2.1.4. Somatic cells’ role

Inherent components of the environment that surrounds the tumor tissue, somatic cells can contribute to neoplasia and promote its development. Amongst them, immune cells and, in particular, lymphocytes, play key roles in the process and thus provide an interface for testing new treatment modalities. Contextually, this section summarizes some of the latest attempts at modeling the interplay between somatic and cancer cells, while potential therapeutic approaches are outlined below.

In the work of Gong et al.,⁴⁷ they emphasized the spatial patterns of ligand PDL1 that is secreted by immune cells (as a way to inhibit excessive activity), but also by cancer cells after protracted exposure. The high spatial resolution of the agent-based model allowed the in silico framework to correlate pre-treatment immune architecture, patients’ features and immune checkpoint inhibitor outcomes. Moreover, they attempted to predict treatment responders using a threshold on the distance between the PDL1-positive cells and the tumor surface. The model was further extended by the authors in subsequent studies^44,48 (see last paragraph in subsection 2.1.3) with the introduction of a new module. Notably, the module was used to simulate the human body in a 4-compartment model, where the bloodstream served as a source for T-cells and myeloid-derived suppressor cells. While they employed ODEs to simulate the dynamics of the whole tumor at the tissue scale, the agent-based model replicated local changes at the cell scale in a small representative region. However, an ODE version of the agent-based model was built to keep consistency between the agent-based and the other models. The two were sequentially solved and used to update each other, with the input values scaled properly. The effects of immunotherapy (i.e., anti-PD1) and different values of the migration and proliferation rates of the cells (encoded by adimensional parameters) on the tumor morphology were investigated. Additionally, an innovative use of a 2D Gaussian kernel density to smooth the discrete spatial distribution of the cells allowed the authors to introduce a new way to locate the boundaries of the tumor-invasive front from digital pathology images. A different therapeutic approach, namely Chimeric Antigen Receptor (CAR) T-cell therapy, was modeled and investigated by Prybutok et al.⁴⁹ Simulations of both a dish and a tissue (where nutrients are thus provided by the vasculature) resulted in the identification of the best treatment strategy. This maximizes cancer cell death by CAR T-cells while minimizing the elimination of low-level antigen-expressing healthy cells. The studies of Beck et al. and Khazen et al.,^50,51 provide tools to characterize the immune response mediated by Cytotoxic T Lymphocytes (CTLs). In the first one, a simple model with space competition was used to prove that both CTLs contact and cytokine secretion are needed for tumor cell killing. In the second, where CTLs could adhere to tumor cells, the optimal effector/target ratios for tumor control were found and the human CTLs killing per capita was quantified. Under certain circumstances, immunotherapy alone might not be enough to eradicate solid tumors. The work of Kather et al.⁵² reports that for patients affected by microsatellite-stable colorectal cancer, effective immunotherapy strategies don’t exist. However, the agent-based model implemented by the authors shows that a combined therapy aimed at boosting the immune system while targeting the stroma can eradicate the simulated tumors in 75% of the runs. A permeable stroma allows the lymphocytes, whose number is elevated following an external injection, to effectively counteract the immune evasion of the cancer cells and avoid the inhibition of the cell-killing mechanisms. Finally, viral-infected cells are used as a proxy to activate the CD4+ T-cells in the model outlined by Jenner et al.,⁵³ where the authors examined the impact of the relative density of the stromal cells on the efficacy of oncolytic viruses for GBM treatment. The oncolytic viruses, simulated as a diffusing field, are uptaken by both cancer and stromal cells and the intracellular dynamics is modeled with ODEs. However, while replication and subsequent lysis occur in cancer cells, the stromal ones act as sinks and reduce the overall viral infiltration. As a consequence, cytokines are not released and the response of CD8+ T-cells is hampered. The simulations, whose outcomes were validated against heterogeneous patient samples, showed that high viral biding rates could be ineffective if the relative density of GBM cells is low. Only an increase in the number of CD8+ T-cells led to a significant reduction of the tumor size.

The interaction between myeloma cells (i.e., malignant plasma cells that undergo uncontrolled growth and proliferation) and bone marrow stromal cells was investigated by Su et al. and Ji et al.^54,55 Myeloma cells are thought to closely collaborate with bone marrow stromal cells in a positive-feedback loop that leads to niche stiffening and mechanical protection from drugs. As multiple myeloma has proven to be able to develop multi-drug resistance and evade the host immune response and relapse, combining multiple therapies could lead to improved outcomes. Both the aforementioned models are robust in silico procedures to test the joint efficacy of different anti-cancer drugs. To quantitatively measure the synergistic effects of the drugs, the authors employed the Loewe combination index.²⁷ Moreover, both the models are hybrid and multi-scale in that ODEs are used to simulate intracellular dynamics, while agents simulate the cells. Besides, Ji et al.⁵⁵ built upon the work of Su et al.⁵⁴ by implementing an immune system within the model. The authors simulated drugs that could target: the myeloma cells, their immune tolerance, the biomechanical phenotype of the bone marrow stromal cells and the communication between the latter and the myeloma cells. The models, which successfully replicated the tumor growth and interactions with the host cells, provide valuable resources to determine the efficacy of multi-drug treatments and the most promising dose combinations.

Other noteworthy models of interactions between cancer and blood cells, as well as other host cells, were presented in the papers of Uppal et al.⁵⁶ and Heidary et al.⁵⁷ There, the role of platelets, key players in metastasis, and fibroblasts, turned into cancer-associated cells, was explored. In,²⁴ an ABM was employed to compare physiological wound healing and tumor-induced angiogenesis to interrogate the interplay between cancer cells and platelets. The model suggests that, by disrupting the physiological setting, tumor edges interfere with the well-orchestrated release of angiogenesis inhibitors. The resulting condition resembles a ‘wound that never heals'. The role of the cross-talk between the endothelium and cancer was further explored by Yan et al.,⁵⁸ who presented a hybrid model of GBM progression. Of note, the model features both normal endothelial cells and trans-differentiated vascular endothelial cells, together with neoplastic stem and differentiated cells. The in silico results illustrated that the combination of therapies traditionally used in isolation can lead to enhanced results in GBM treatment. Modeling of the endothelium and the perivascular niche of the GBM was also explored by Randles et al.,⁵⁹ where the authors employed ABM to optimize an existing therapeutic regimen. Scalable simulations were combined with simulated annealing (i.e., an optimization algorithm used to efficiently explore the parameters space) to infer the best timing for both chemo- and radiotherapy administration. The parameters obtained were then used to implement and test the schedule in vivo. This results in an improved outcome and thus provides experimental evidence for the initial assumptions regarding the stem-like cancer cell differentiation and translocation mechanisms.

2.2. Machine learning approaches and ABM

Typically formulated as an optimization problem, simulation calibration has become a very interesting candidate in the ML community. This is because, on one side, in silico models can generate large quantities of data and, on the other side, biological or/and medical data is sometimes hard to collect or very expensive to acquire. Learning-based approaches offer an attractive alternative to optimization-based calibration approaches. They are especially interesting, as they also need to update behavior rules embedded in ABM, as shown in Figure 2(c). Such approaches tackle the realistic reproduction of mechanistic dynamics in biological systems by learning the mapping from clinical data and model parameters to a performance metric (see Figure 2(a)). In this space, the approaches are rather diverse. For instance, the work of Barde et al.⁶³ presented the first empirical application example of a novel probabilistic model calibration methodology. This is a systematic study on methodologies to calibrate and benchmark agent-based models. Their focus is on the comparison of ABM simulations carried out with standardized criteria based on accepted information-theoretical measures for forecasting performance and explanatory power. The purpose was to show that there are efforts in using ML to improve ABM calibration. Taking agent-based models closer to the data remains an open challenge, especially when considering biological processes. This involves calibrating the agent-based models to minimize the gap between simulation and reality. Computational intelligence, as a sub-field of artificial intelligence, is an exceptional candidate for achieving this due to its focus on computational principles and knowledge extraction based on biological principles, e.g., evolutionary programming, fuzzy logic, artificial immune systems etc. ML algorithms have the power to extract insights from experimental data, hence capturing those peculiarities that a model doesn’t typically capture. Data distribution changes, bias, drift, and other phenomena can be learnt and used to judiciously parameterize agent-based models and ABM simulations. This aspect was the direct focus of the research carried out in the work of Lamperti et al.⁶⁴ The authors explicitly tackled parameter space exploration and calibration of agent-based models combining supervised ML and intelligent sampling. This way the researchers could build a surrogate meta-model that provided a fast and accurate approximation of simulated model behaviors, dramatically reducing computation time. More precisely, the ML surrogate (i.e., an adaptive twin) facilitated large-scale explorations of the parameter-space, while providing a powerful filter to gain insights into the complex functioning of agent-based models capturing complex dynamics across scales. Using computational intelligence for learning ABM simulation parameters, Singh et al.⁶⁵ employed artificial life optimization. Their in silico framework implemented a hybrid model using micro-simulation and ABM techniques to generate an artificial society. The agents in this model derive their decisions and behaviors from real data (i.e., a micro-simulation feature) and interact among themselves (i.e., an ABM feature) to proceed in the simulation realization. Such approaches have been reported to map very well on the problem structure, as it is typically found in cancer cell biology, where local cell behaviors propagate in upper tissues or organ properties change.

The work of Niida et al.⁶⁶ proposed a very computationally powerful and parallelized approach to deal with model uncertainty (and its impact on calibration) in ABM. They highlighted the role of interactive visualization to help identify suitable model parameters. This is crucial in handling the highly nonlinear dynamics of cellular interactions and cancer evolution. In this context, the calibration process of ABM simulations is structured around the concept of adaptability. By adaptability, we refer to the fact that model parameters can impact one another, either via direct relationships (e.g., diffusion constant and physical properties of the extracellular matrix) or implications on system dynamics (e.g., vascularization and cellular metabolism can both affect the growth dynamics of cancer). More precisely, as the model complexity increases, so do the constraints among the model parameters and their inter-dependencies with respect to given summary statistics.

Using empirical priors and statistical learning, Lima and his colleagues²⁸ proposed a moment-based Bayesian inference to account for the stochasticity of the coarse-grained ABM in a tumor growth multi-scale model. The approach presented very clever methods for quantifying uncertainties due to limited temporal observational data of cancer growth and staging at different spatial and temporal levels. Overall, the approach reduced the computational time of ABM simulations while reliably/realistically capturing tumor dynamics and its inherent nonlinearity. Using a hierarchical optimization simulation for calibrating the agent-based model, Amaran et al.⁶⁷ integrated optimization techniques into simulation analysis. The primary goal of simulation-based optimization is to improve the performance of the models through Monte Carlo processes. More specifically, the Monte Carlo simulation allowed the system to find the optimal set of parameters for a given criterion based on a modular thresholding method. Here, by performing sensitivity analysis in conjunction with parameter optimization, the authors gained a comprehensive understanding of the process dynamics. Ultimately, this empowered the researchers to drive continuous simulation improvement and achieve superior quality. The work of Akasiadis et al.⁶⁸ stands out through the ingenious use of a typical calibration methodology. Their in silico method only considered numerical optimization and the goodness-of-fit only captured the quantitative aspects of the calibration. Finally, neural learning techniques were employed for comparing spatial simulations to tumor imaging, going beyond basic metrics retrieved from tumor images and ABM simulations. In this manner, such algorithms may evaluate the model fit quantitatively. More recently, the work of Cess and Finley⁶⁹ employed representation learning and a neural network to project an input into low-dimensional space, which is a representative example. The authors utilize a neural network to represent the ABM simulations and tumor images as low-dimensional points, with the distance between them, serving as a quantitative indicator of their differences. The main contribution of the study is that it draws on neural networks to create low-dimensional single-valued representations of both tumor image volume calculation and ABM simulations. The distance between these single point values serves as a quantitative indicator of the differences between the two (i.e. metric of performance of the ABM simulation at the core of a calibration process). In order to guide the reader, we tabulate the relevant machine learning approaches and their successful application for ABM in Table 2.

Table 2. Machine learning tools for ABM and overview of the relevant state-of-the-art literature.

Approaches	References
Computational Intelligence	^Citation65,Citation66,Citation68,Citation69
Probabilistic techniques	^Citation63,Citation64,Citation67
Bayesian and Statistical Learning	^Citation28,Citation70,Citation71

2.3. Multi-scale or multi-level numerical methodologies in ABM

Multi-scale or, as some investigators prefer to use as a term, multi-level in silico modeling systems in biology and pharmacology are robust methodologies that can aid researchers in understanding and probing the fundamental mechanisms of biological phenomena and in clinical applications. The excellent review of Morvan⁷² lays out clearly the definition of multi-scale methodologies using agent-based models, although in his survey he adopts the term ‘multi-level agent-based modelling’ instead of ‘multiscale’ since the latter, according to him and to Gil-Quijano et al. (see reference in⁷²), has a restrictive meaning as it focuses on the spatiotemporal extent of levels and not on the interactions and organization within the biological system under investigation. We will thus use the term of a ‘multi-level agent-based model’ in this survey whereas, for semantic reasons, we will preserve the term ‘multi-scale model’ if used in the cited works below. Some examples of early attempts at multi-level approaches using ABM that considered coupling cell-scale to molecular-scale dynamics include the paper of Athale et al.⁷³ They presented a model of gene – protein interactions integrated in an agent-based model system to probe the ability of brain cancer cells to ‘switch’ between migrating and proliferating phenotypes, while test how molecular species interact with other molecules within and across sub-cellular compartments. However, in order to keep the review relatively short, we limit the depth of our survey to the relevant multi-scale/multi-level ABM papers in biology and biomedicine of cancer presented during the last decade. Older review articles, such as those of Deisboeck,⁷⁴ Stamatakos,⁷⁵ Walpole⁷⁶ and their colleagues, present the highlights of published multi-scale in silico models of biological systems in cancer, cardiovascular and biomedicine, and demonstrate those early successes using agent-based models in the respective context.

2.3.1. Multi-scale/multi-level methods that couple ABM with differential equations-based solvers

There is a fair list of published papers that demonstrate the coupling between discrete systems – as in the case of agent-based models – and continuum-based models (i.e., these usually formulated using numerical methods that solve PDEs) or network models (i.e., such models are typically described using a set of coupled or non-coupled ODEs) to describe the biological cross-talk amongst different spatial scales by considering pertinent numerical algorithms and techniques suitable for the modeling task. Carcinogenesis dynamics unfolds across both spatial and temporal scales. ABM approaches are a suitable candidate to capture not only the molecular dynamics but also the aggregated dynamics of communication and regulatory networks, as shown by the very interesting multi-level ABM tools developed by Montagna and her colleagues.⁷⁷ Here, drosophila embryo cells are modeled as agents that divide, move, secrete/uptake substances, while Montagna’s method encompasses also the balance of the molecules that mediate cell-to-cell communication and a gene regulatory network to simulate the molecular biology of the cells (i.e., the reactions taking place inside the cells). Their multi-level ABM methodology was implemented via a multi-threaded discrete event scheduler using software Repast Simphony to simulate the expression patterns of the embryo cells against experimental evidence from the FlyEx database. Zhang et al.⁷⁸ proposed an agent-based brain tumor model that encompasses intercellular level to describe cell – cell interactions intracellular-scale dynamics by employing a system of ordinary differential equations to describe selected molecular pathways relevant to glioblastoma multiforme pathophysiology (i.e., phenotypic switches in cells from migration to proliferation), and the ‘tissue-scale’ to model the balance of chemoattractant concentration (through isotropic diffusion, secretion, and consumption by the cells). The main focus of this work was in the design and development of their in-house C++/CUDA implementation of the multi-scale agent-based model, which was parallelized with respect to both the chemo-attractant diffusion and the intracellular signaling processes using graphics processing units (GPUs) computing. They reported a considerable computational speed-up of the GPU-based design of the multiscale ABM simulator compared to the one of the sequential design – this was amongst the few early works that demonstrated the potential of multi-scale ABM to simulate real-time cancer progression. Cai et al.⁷⁹ developed a three-dimensional hybrid cellular automata model, which is part of the family of agent-based models, to study the dynamics of tumor spheroids and probe the effect of hypoxia, cell phenotypic behavior due to microenvironment biochemical factors. They solved three coupled reaction-diffusion equations to simulate the dynamics of the ECM, oxygen and ECM-degrading enzymes and communicated in a partitioned fashion the solution to an ABM simulator of the cell dynamics similar to the methodological approach of Zhang. Alfonso and his colleagues⁸⁰ presented a comprehensive cancer model to study in silico immune cell infiltration and interactions in the breast ductal lobular epithelium. Following an ABM formulation, their multi-scale model accounts for myoepithelial, luminal, and immune cells (each type allowed to reside in a separate lattice in 2D space) whose behavior included immune cell trafficking, cell migration, immunosuppression, epithelial cell proliferation, damage, programmed, and induced cell death, and cell lysis. Also, they modeled the transport and secretion of chemokines that control the induction of an immune response in the terminal ductal lobular units of the breast epithelium. They calibrated their model from imaging data of immuno-histochemical epithelial, vascular, and immune cell markers from healthy women, and they investigated recurrent inflammation during physiological menstrual cycles and normal hormone levels, while they analyzed in silico parameter perturbations that can lead to carcinogenesis. Interestingly in the latter part, as suggested by the agent-based model simulation results, they observed that epithelial damage induced higher variations in immune cell infiltration. Later, Gong et al.⁴⁷ presented a multi-compartmental, multi-scale model of tumor development and anti-tumor immune response, which included interleukin-2 (a cytokine attributed to immunological homeostasis and classification), cytotoxic T lymphocytes and neoplastic cells. Cells were set to interact in an off-lattice 3D space and follow a set of rules including division, migration, cytotoxic killing, and immune evasion. They developed an in-house C++ code for their ABM simulator to explore spatio-temporal tumor immune response to PD1 and PDL1 inhibition. To reduce the number of simulations used to quantify simulation uncertainty, they employed a particular statistical methodology for Monte Carlo simulation based on the LHS. This approach permitted to efficiently explore the parameter space and conduct a comprehensive sensitivity analysis with a small number of samples, where each potential parameter combination was ranked based on the accuracy of the model predictions and performance. The in silico results generated interesting findings on the spatial patterns of different cell types without treatment that resembled patterns reported in cancer patient biopsies, and that ABM-simulated response to anti-PDL1 treatment is affected by the neoantigen characteristics of a patient. Letort and her colleagues⁸¹ presented an open-source simulator, PhysiBoSS, which combined intracellular signaling using Boolean modeling and multi-cellular dynamics and behavior using ABM. As a use-case to demonstrate their modeling tool was modeling cell-fate decisions in response to treatment of cytokine Tumor Necrosis Factor (TNF) in order to illustrate the cell – cell communications. They also explored in silico the effect of different treatments and the behavior of several resistant mutants, while also testing the dynamics of cancer cell population with respect to the spatial heterogeneity of biochemical cues and resources (i.e., oxygen). Pally et al.⁸² presented in a paper both their experimental and computational work to interrogate cancer cell migration into cellular interactions with the basement membrane and its remodeling, the transition from basement membrane to type-I collagen, and the subsequent remodeling of, and migration within, type-I collagen in the context of early breast carcinomatosis. They built a multi-scale 3D organo- and pathotypic experimental assay, with the ABM implementation based on a cellular Potts model using the open-source software CompuCell3D. The model encompassed cancer cell proliferation and apoptosis, cell adhesion with the ECM, and the TME with respect to matrix remodeling through reaction – diffusion – based morphogen dynamics of Matrix Metalloproteinases (MMPs) and tissue inhibitors of MMPs. Pally and colleagues designed a culture model of MDA-MB-231 cells to form reconstituted basement membrane-coated suspended clusters to mimic the invasion patterns of breast cancer cells in vivo and probe how ECM density, MMPs and N-linked glycosylation concentration impacts cancer cell invasiveness. Their approach considers a hybrid system of partial differential equations, for the spatio-temporal evolution of the densities of the epithelial-like cancer cells, ECM and the MMPs, and stochastic differential equations, for the time evolution of the mesenchymal-like cancer cells including their migration along ECM gradients with these cells described as particles. Sfakianakis and his colleagues⁸³ presented a multi-scale modeling framework for cancer invasion of the ECM by the combined action of epithelial-like cancer cells and mesenchymal-like cancer cells. The implicit/explicit Runge-Kutta finite volume numerical method was employed to solve the continuum-based part of the model, and an explicit Euler-Maruyama scheme was employed for the solution of the stochastic particle-based part of the model. In the same year, Macnamara and her colleagues⁸⁴ proposed a multi-level in silico method to simulate the growth of a solid tumor, migration of cancer cells within heterogeneous tissue, and the effects of fiber and vascular structure in cancer development. In contrast to previous works, they modeled each cell as a viscoelastic deformable sphere with cell dynamics governed by ordinary differential equations to describe cell – cell forces, cell random migration and cell to ECM friction, and cell behavior (i.e., cell cycle, growth, and birth) as in agent-based models. Cancer cells were allowed to occupy voxels within a Cartesian grid, which was also populated with additional agents that were contiguously structured such that they formed capillary vessels. They also accounted for the balance of nutrients and oxygen at the tissue level in the form of reaction-diffusion equations that were discretized with finite elements. Macnamara built an in-house C++ solver for cell dynamics and used the FreeFem++ platform for as an FEM solver, while they ran (non-cancer specific) simulations of tumor growth around an arbitrarily defined blood-vessel network.

2.3.2. Hybrid formulations of multi-scale/multi-level procedures

In addition to the multi-scale/multi-level numerical methodologies that couple ABM with continuum-based techniques (e.g., finite elements, finite differences), several hybrid in silico cancer modeling procedures have demonstrated the capacity to simulate neoplasia development from micro to meso to macro-scales, as well as hybrid techniques that incorporate machine learning and optimization algorithms for model parameter exploration and verification. The following paragraphs outline the most notable examples of such novel hybrid formulations that relate to ABM.

De Montigny et al.¹⁴ proposed a hybrid multi-level in silico cancer model that was tailored to simulate glioma growth. Their off-lattice agent-based model encompasses (host and neoplastic) cell growth, division, migration and adhesion, the dynamics of the extracellular matrix, the effects of oxygen and nutrient availability in cell survival, or the switch of cancer cells into a hypoxic or necrotic state, and the signaling triggered by chemical cues and growth factors. In contrast to all the above-mentioned papers, the multi-level formulation in¹⁴ uniquely couples a continuum-based finite element model for the solution of reaction-diffusion equations (i.e., to predict the balance of cytokines, growth factors, and oxygen) at the tissue scale with an agent-based model via the volume averaging method. The hybrid simulator was developed by coupling an in-house FEM solver FEB3 and the open-source platform BioDynaMo for the ABM simulations. Using the in silico framework, they examined the impact of cell – cell and cell – ECM interactions in (macroscopic) tumor growth, brain tissue perfusion and tumor necrosis, as well as they assessed computationally the differences between low- and high-grade glioma growth, vascularization and necrosis and compared to experimental data from the literature.

Lima et al.²⁸ presented a reduced-order ABM methodology coupled with Bayesian inference modeling for parameter calibration to manage the stochasticity of the agent-based model. In terms of the multi-scale modeling, ABM was used for the phenotypic behavior and the geometric properties description of the cells, while the dynamics of nutrients was modeled at the tissue scale as a reaction – diffusion process. Their hybrid model was demonstrated to simulate the development of a BT-474 human breast cancer cells in vitro, using time-resolved microscopy data, and employed a moment-based Bayesian inference method to quantify the uncertainties owing to limited temporal observational data of carcinoma’s growth.

Ponce-de-Leon and colleagues⁸⁵ presented a multi-scale model of cancer cell dynamics with signaling for TNF-receptor dynamics, as in their previous intracellular signaling work in Boolean modeling.⁸¹ Cancer cells were modeled as agents residing in 2D or 3D lattice, which accounts for the presence of oxygen and the cytokine TNF. Their in silico model was built by combining the open-source ABM simulator PhysiCell and the software PhysiBoSS that was developed by the authors. Subsequently, they integrated their simulator with an Extreme-scale Model Exploration with the Swift platform to carry out exploration tests of the agent-based model parameter space which was ultimately deployed to optimize dosage-specific treatments for tumor regression. They probed for the effect of the spatial distribution of cancer cells on the treatment parameters optimizing the supply strategies in cell monolayers and three-dimensional tumor spheroids; similarly, they interrogated the robustness of the effective treatments with respect to the cell population heterogeneity of the cancer cells. Following the modeling work in,⁸⁵ Ruscone et al.⁸⁶ proposed an enhanced multi-scale model to interrogate possible targets that can help block or suppress the invasive phenotypes of cancer cells. More specifically, the improvements are focused at the intracellular scale where they incorporated mechanisms of epithelialto-mesenchymal transition and cell metastasis. They used the in silico model to test the role of tumor protein 63 and metalloproteinase MT1-MMP in tumor invasion, as well as that of the tyrosine kinase protein SRC in an epithelial monolayer, while they also tested possible drug candidates to block migration in the ECM of cells that have undergone epithelial-to-mesenchymal transition. Tsingos and her colleagues⁸⁷ presented a spatially inhomogeneous cellular Potts model to simulate cell migration in a fibrous matrix. To overcome the substrate homogeneity of the Potts model, they coupled it with a ‘background’ bead-spring biomechanical model of the ECM where fiber networks were modeled using molecular dynamics. A unique feature of their multi-level approach was the incorporation of contractile pulling by the cells through discrete focal adhesion-like sites on the fiber network. Despite their simulation experiments coming from the angle of morphogenesis and tissue healing, their in silico method can be adapted to study cancer cell infiltration and invasiveness. Also very recently, Miller et al.⁸⁸ presented a multi-scale modeling approach to evaluate the effect of chemotherapeutics on patient tumors based on metabolomic analysis results of lung cancer biopsy data. Despite their multi-scale model being based on previous work from the same group (see reference³² therein), they integrated metabolomic analysis evidence (from patient tumors) and modeling. Interestingly, they built a synthetic dataset using Monte-Carlo by resampling selected parameter values of the multi-scale model to simulate chemotherapy, while they considered seven evaluation metrics to quantify the tumor response; however, as the authors suggest, their results need further validation with metabolic evidence from different and larger patient datasets.

2.3.3. Towards robustness and high fidelity of multi-scale/multi-level methodologies

The literature survey laid out in the previous subsections illustrates a remarkable development track record of multi-scale and multi-level methodologies using ABM in oncology over the past decade. The computational cost and complexity to test the numerical stability of a multi-scale/multi-level in silico procedure, to probe the sensitivity of the numerical schemes involved, and to calibrate the models across multiple spatial scales remains a challenge however. In addition to the computational cost emerging when, it comes to simulating thousands or even millions of agents and the inherent stochasticity of an agent-based model, this stipulates realizing a great number of simulations to accurately represent the statistical features of an in silico cancer model. Thus, the effort toward achieving a high level of robustness and fidelity in a multi-scale/multi-level ABM formulation elevates tremendously in proportion to increasing the size of the biological system under investigation, as well as with the quantity and the modalities of data coming from the laboratory or the clinic.

Interestingly, as it is presented in subsection 2.2, investigators have attempted to amalgamate sophisticated machine learning and optimization algorithms for learning the simulation parameters, to quantify the model uncertainty and its impact on calibration on agent-based models (e.g., as in⁶⁵ and²⁸). However, the majority of multi-scale and multi-level methodologies have used data to constraint and validate the cancer model predictions on a single scale, e.g., usually at the tissue level through tumor size measurements, or in the order of the largest temporal scale, typically at the order of the time-duration of a preclinical cancer experiment. The cited papers that follow illustrate the gradual evolution of relevant cancer in silico models to accomplish calibration across the scales. Among the early attempts to simulate multicellular tumor spheroids that mimic the TME dynamics was that of Cai et al.,⁷⁹ who modeled cells as agents while the balance of nutrients and enzymes was modeled following a continuum-based method. Despite their method encompassed multiple spatial scales, i.e., the dynamics of cells and proteins, the simulation predictions were tested using history plots of the tumor spheroids with respect to size. Later, following a similar modeling path, Mao et al.⁸⁹ presented a hybrid continuum/agent-based model for HCT116 tumor spheroids to simulate hypoxia-dependent interactions between ionizing radiation and a hypoxia-activated prodrug; this in silico tool was used by Hong et al.⁴⁵ to build a pharmaco-kinetic/pharmaco-dynamic model and probe the bystander effects of hypoxia-activated prodrugs in cancer cell killing. To inform their model at the continuum scale (i.e. the average rate of diffusion of the drug molecules in the medium) they adopted parameters from the literature, while to calibrate the agent-based model they used in vitro data coming from flow cytometry analysis, confocal microscopy imaging data of the spheroids, and fluorescent staining of the cancer cells to mark their protein expression. In,⁵² they presented an in silico model of lymphocyte – tumor – stroma interactions to interrogate the response to immunotherapy and stroma-targeting therapies on human colorectal cancers. As in,⁴⁵ Kather et al. employed data from horizontal in vitro migration experiments on lymphocytes to inform the agent-based model, as well as ex vivo measurements (based on morphological processing of single slice images) on histological human tumor tissue samples to estimate the proliferation, apoptosis and distance to necrosis parameter values, and provide a quantitative basis for the cell-scale modeling. Rahman et al.³⁷ coupled in ‘space and time’ multi-scale cancer model spanned from tissue (using FEM) to cellular (using ABM) and subcellular scale, with the latter being represented by signaling pathways. In a similar fashion, de Montigny et al.¹⁴ integrated ABM with the FEM using a volume-averaging formulation to build a multi-level brain cancer simulator. Both models, however, were calibrated and tested from observed data at a single scale (tissue level), i.e., history plots of the tumor volume and average volume fraction of cell groups. With an exception to the model in¹⁴ where tissue-scale (FEM) parameters for tumor and host cell dynamics were inferred (a process often called data upscaling) from the cell-scale (ABM) simulation predictions, in³⁷ the modeling parameters on the lower scales (cellular and subcellular) were either adapted from the literature or estimated. The multi-scale approach of Lima²⁸ employed in vitro data to separately inform (at the cell scale) the mechanistic agent-based model of human breast cancer cells’ phenotypic behavior, and (at the continuum tissue scale) the transport and balance of glucose concentrations and cytokines in general. The important contribution of their paper concerns the Bayesian inference concept applied for the time-dependent sensitivity analysis of the in silico model and to interrogate the model parameters’ space. Using lightweight neural networks, Axenie et al.³⁴ extracted the mechanistic relations governing phenotypic staging and tumor volume development.⁹⁰ Very recently, Cesaro and her colleagues⁹¹ demonstrated their multi-level TME simulator that couples mechanistic agent-based models with PDE-based solvers in two dimensions. An innovative feature of their paper was the data-driven strategy they adopted to inform the agent-based model using bulk gene expression data from The Cancer Genome Atlas database. They also used evidence from the scRNA-seq dataset of human colorectal cancers to calibrate with respect to the tumor mutational burden and the inhibitory immune checkpoint that suppresses T-cell activation and to estimate the cell (HCA, T-reg, CAF, immune) fraction and in their model. The multi-scale model of Ponce-de-Leon⁸⁵ considered the Covariance Matrix Adaptation Evolutionary Strategy for the numerical optimization of their agent-based model to analyze the treatment parameters of the tumor necrosis factor cytokine and its effect on cancer regression. Despite their model being inherently multi-scale, their in silico experiments were calibrated on macroscopic quantities, i.e., drug dose, injection time and duration, and tumor size. In the same year however, Miller et al.⁹² presented a multi-level approach that is unique in that it proposes to link tumor metabolomic measurements from patients into the mathematical model for tissue-scale behavior of a carcinoma progression or control, the development of angiogenesis, the effect of chemotherapy, etc. However, an important limitation of the cancer model concerns that its behavior depends mainly on the metabolomic data available, and how they are appropriately weighted and combined to determine the effect on the (mechanistic) model parameters. Finally, Ruiz-Martinez and colleagues⁴⁸ proposed a hybrid ODE-based/ABM simulation tool investigate therapeutic strategies related to anti-cancer immunity and immune checkpoint inhibition. The rules for the cancer (stem-like, progenitor, and senescent) and the immune (CD8+ T and Treg) cells were defined in the agent-based model provided in vivo evidence from the literature,⁹³ while for the calibration of the differential and algebraic equations (120 in total) they employed a model. Their, a mechanistic modeling method often used for drug discovery, was based on a relevant one for non-small cell lung cancer and incorporated data from single-cell RNA sequencing.

2.4. ML-assisted calibration: a two-sided sword

Agent-based numerical procedures that are supported by ML models, as has been briefly outlined in subsection 2.2, can significantly contribute to the generation of agent-based models with suitable model parameters. However, automated calibration can also entail challenges and difficulties. This section elaborates on the comparison between ABM with and without ML assistance. Notably, we provide further explanatory comments to establish the wider context of ML and its involvement in ABM.

From a modeler’s perspective, a pure ABM procedure is, as is more generally mechanistic modeling, often seen as demanding with regard to the determination of model variables. As explained above, ML techniques can be employed to efficiently search the parameter space of mechanistic models and determine optimized model parameters. However, an important criterion in biomedical models is that parameters need to be ‘biologically plausible,’ if not (ideally) directly experimentally informed model parameters. An agent-based model where model parameters were estimated from the literature, without ML-assisted optimization, is the study of Macklin et al.⁹⁴ However, due to limitations in data availability, it is usually impossible to infer all model parameters from the literature. ML can help address this problem. For instance, Demetriades and his colleagues⁶² employed ML to infer various parameters on the pharmacological impact of cancer drugs. B¨orlin and his colleagues⁹⁵ employed model parameters obtained both from the literature as well as ML-derived ones. Jalalimanesh et al.⁹⁶ make use of reinforcement learning to optimize radiation treatment. Overall, the need to account for biomedical plausibility in model parametrization highlights the importance of interdisciplinary collaboration particularly for mechanistic, multi-scale models⁹⁷ (more about multi-scale models in subsection 2.3).

A crucial criterion in computational modeling is ‘explainability.’ This aspect has recently gained much attention due to the fact that large language models are very problematic with regard to gaining insights into the human-understandable causes of outputs. Especially when it comes to biomedical applications, explainability is of utmost importance. Clinicians need to understand exactly the reasons behind their findings to adequately inform their decision-making process, given the potentially life altering impacts of these. Mechanistic models, as in for example the studies of Macklin et al.⁹⁴ and de Montigny et al.,¹⁴ produced experimentally verifiable hypotheses that can ultimately lead to deeper insights on glioma cancer growth and cell necrosis. Notably, many ML methods suffer from the same issues as large language models as they can be perceived as black boxes that may have excellent performance with regard to a given biomedical problem, but nevertheless limited clinical benefits. Therefore, modelers should be aware of the potential pitfalls when interfacing mechanistic models with ML algorithms, especially when it comes to practical impact in the clinical setting. However, it is notable here that ML methods come in different types, i.e., fundamentally black-box and white-box models. Indeed, ML approaches that enable explainability (or interpretability) exist, as described in the work of Linardatos et al.⁹⁸

In view of the trend of in silico models’ rising complexity and the number of parameters they consist of, then comes the need to rigorously assess ‘robustness’ and ‘sensitivity’ of the model, for instance with regard to the model parameters and/or different initial conditions. Given that no two biological systems are exactly identical, a theoretical model must tolerate changes to model parameters, at least within reasonable boundaries. Due to this inherent variability, biological systems in cancer usually comprise redundancy, checkpoints, and control loops that permit for changes without endangering important outcomes. In the context of computational modeling, limited changes to model parameters should not lead to implausible in silico outcomes.

Additionally, a crucial goal of computational modeling is the generation of hypotheses and experimentally verifiable predictions. To this end, ABM sensitivity analysis can help gain insights into the impact of the model parameters. A direct but basic way to accomplish such tests for ‘robustness’ and parameter ‘sensitivity’ is to execute the model with different model parameters in a grid-like manner – in this methodology, ML cannot be deemed pertinent. However, it may come to a case where some parameters of an agent-based model may have a stronger impact than others. In this direction, ML can be a helpful tool to efficiently sample, identify, and rank model parameters in terms of their importance when considering these as ‘features.’ This is called Variable Importance Analysis and, for instance, the random forest ML method is commonly used for such analysis, as demonstrated in the work of Pereda et al.⁹⁹ In a similar approach, Retzlaff and his colleagues¹⁰⁰ used decision trees for Variable Importance Analysis in their agent-based model; the authors indicated that cell cycle duration and motility in the context of solid tumor metastasis are the most important factors with regard to therapy resistance.

An aspect that requires consideration for the usage of mechanistic modeling and ML-assisted mechanistic modeling is ‘scalability.’ Given that mechanistic modeling should be, for the sake of biological plausibility, based on local information exchange only, its simulation can naturally make use of parallelized and distributed computing.^101,102 For numerous ML methods, their adaptation for large-scale applications requires often custom efforts, since every algorithm has a distinct communication pattern, as demonstrated in the work of Verbraeken et al.¹⁰³ Along those lines, synchronization requirements among nodes can vary across ML methods, as well as suitability for specific hardware (e.g., CPUs versus GPUs). Overall, the smooth and efficient interfacing with mechanistic modeling constitutes nevertheless a challenging task that remains to be addressed in the future.

Table 3. Comparison of pure MM versus ML-MM.

Criterion	Well-suited for	References
Identification of model parameters	ML-MM	^Citation62,Citation69,Citation85,Citation95,Citation96
Accordance with experimentally measured parameters	MM	^Citation94,Citation104
Explainability	MM	^Citation14,Citation85,Citation94
Robustness and sensitivity	ML-MM	^Citation99,Citation100
Computational scalability	MM	^Citation101,Citation102

Clearly, specific problems in cancer biomedicine require consideration of the associated advantages and disadvantages, and there is no one-size-fits-all approach for any given scientific quest. ML can be a highly valuable asset for researchers employing ABM. However, the combination of these distinct approaches can also entail challenges, and so their symbiotic application is not necessarily warranted. Nevertheless, specific aspects that need to be considered on a case-by-case basis can be appreciated (Table 3).

2.5. Fusing mechanistic and learning approaches (physics-informed systems)

The in silico modeling approaches discussed above focus on the calibration of Markovchain/ODE-/PDE-based models to simulate processes involved in neoplasia. Such mechanistic models account for assumptions about the dynamics of the systems in both the temporal and spatial dimensions, as outlined in subsection 2.3. Advanced numerical methods and high-performance computing enable high-fidelity simulations of such calibrated mechanistic models to run at scale. However, most approaches for calibrating biological systems’ models focus on fit quality. The common noun in the current approaches landscape demonstrates that the calibration error, which varies depending on the optimization approach, reaches an insurmountable barrier that can result in a standstill in selecting the “optimal” model. This model selection process reaches another stale point when considering capturing corner cases in the spectrum of the modeled system’s behaviors. The current approaches look at modifications of the models themselves with adaptive features or the analysis of the influence of the system’s characteristics on the system’s behavior in corner-case situations. Learning-based approaches offer an attractive alternative to these optimization-based calibration approaches. An additional advantage is the explanatory power that physics can offer when building ABM behavior rules for the simulation. Such approaches tackle the realistic reproduction of the biological system’s behaviors by learning the underlying mapping from data and model parameters to a performance metric or goodness-of-fit criteria of plausibility. But, in order to gain the best of the two worlds, mechanistic ‘biases’ can be ‘injected’ into ML models and leverage the power of learning from large amounts of data through a ‘directed’ search for the solution, in other words, the realm of physics-informed ML modeling. Fusing mechanistic biophysical models and learning algorithms amounts to introducing appropriate observational, inductive or learning biases that can direct the learning process toward reaching physically plausible solutions. This new conceptual framework of physics-informed ML framework coined by Karniadakis and his team¹⁰⁵ proposed training ML models from additional information obtained by enforcing the physical laws (for example, at random points in the continuous space-time domain). Such physics-informed learning integrates (noisy) data and mathematical models, and implements them through neural networks or other kernel-based regression networks for calibrating or optimizing ABM parameters. Practically, this can be done by introducing inductive, observational, or learning biases in the learning process, under the form of a loss function, regularization term, or event calibration metric. Multiple candidate approaches that focus on physics-informed learning for simulation calibration have been proposed, each one focusing on a different component of the overall problem. Rutter et al.⁷⁰ proposed an incremental mixture approximate Bayesian computational procedure for colorectal cancer simulation calibration. They used a simulated sample from the posterior distribution of model parameters given calibration targets to inform national cancer screening guidelines. For instance, to achieve computational gains in large-scale simulations the work of Wood et al.¹⁰⁶ developed a novel computationally efficient method for direct generalized additive model smoothness selection. Designed as a highly stable, but carefully structured, calibration system, the proposed approach achieved a computational efficiency that led, in simulations, to lower mean computation times than the schemes that are based on working model smoothness selection.

Because it may be challenging to abstract and define the rules that control an agent-based model from experimental data, at least in an objective manner, there is a particularly synergistic potential to utilize ML to help infer the most effective, system-specific ABM rules, as shown in the work of Sivakumar et al.¹⁰⁷ Once such rule sets are developed, a large volume of ABM simulations can produce a plethora of data, and ML can be used in that setting as well. For instance, statistical measures that accurately and meaningfully characterize the stochastic outputs of a system and its features are one use of ML in this context. ABM simulations can produce credible (realistic) datasets to subsequently use for training ML algorithms (e.g., for regularization, to prevent overfitting), as an example of synergy in the other direction (from ABM to ML). In an effort to develop a general-purpose computational framework, Spolaor et al.¹⁰⁸ introduced a novel approach for the analysis of hybrid models consisting of a quantitative (or mechanistic) module and a qualitative module that can reciprocally control each other’s dynamic behavior through a common interface. This qualifies as a mix of inductive and observational biases. More precisely, the system of Spolaor and colleagues took advantage of precise quantitative information about the temporal evolution of the modeled system through the definition and simulation of the mechanistic module. At the same time, it described the behavior of biophysical model components and their interactions that are not known in full detail, by exploiting fuzzy logic in the definition of the qualitative module. Such approaches are deemed to be suitable for the analysis of cancer morphogenesis, an intricate chain of biological mechanisms that enable cell populations to reproducibly self-organize into specific shapes or patterns. Through physics-informed simulations, a modeler can modulate the state through controlled signal transduction on a range of spatial and temporal dimensions that include a variety of mechanisms and systems, as demonstrated in the work of Glen et al.¹⁰⁹ The path of hybrid approaches is further strengthened by the work of Ward et al.¹¹⁰ which proposed a dynamic calibration of agent-based models using data assimilation. More precisely, investigators tackled the question of how such models can be dynamically calibrated using the ensemble Kalman filter, a standard method of data assimilation. The work developed a new type of Kalman filter-based system in a simple setting for data assimilation and fusion in ABM calibration for cancer development. Combining probabilistic machine learning in a physics-augmented framework, the work of Moon et al.¹¹¹ claimed a new technique dedicated to improved calibration and validation of agent-based models. The framework identified periods of deviation between the simulation and the observation with the Hierarchical Dirichlet process hidden Markov model. This allows the framework to automatically calibrate the temporal macro parameters by searching parameter spaces with a broader likelihood of validation for tumor growth under a compromised immune system.

When considering clinical sequencing of surgery and chemotherapy, the work of Axenie and Kurz³⁴ illustrated how a physics-informed ML system can extract the pharmacokinetics of a common breast cancer chemotoxic medication while also concurrently learning the patterns of tumor development in a variety of breast cancer cell lines. In a very recent study, Beik et al.⁷¹ introduced a Bayesian multi-model inference methodology for a dual purpose. On one end, the model quantified how mechanistic hypotheses can explain given experimental datasets, basically by attaching the probabilistic explanation to data peculiarities. On the other, the model demonstrated how each dataset informs a given model hypothesis, thus, enabling hypothesis space exploration in the context of available data. The approach was successfully used to probe standing questions about heterogeneity, lineage plasticity, and cell – cell interactions in tumor growth mechanisms of small-cell lung cancer. The methodological approach in⁷¹ complements the physics-informed ML approaches with a strong probabilistic framing of hypothesis testing and variable interactions in cancer modeling.

On the other hand, when considering the translation of in silico models in cancer progression, physics-informed ML approaches have been successful in tumor volume prediction after learning without supervision tumor phenotypic stages from breast cancer cell lines (e.g., see³⁴) as another simultaneous task. For the purpose of optimization-free calibration of ABM simulations, the work of Axenie et al.¹¹² introduced a physics-informed fuzzy logic calibration system. Using spatiotemporal models of agents’ interactions, the in silico system could regress, based on human experts, the plausible solutions of the goodness-of-fit metric (e.g., Akaike Information Criterion, Root Mean Squared Error (RMSE), etc.). Benefiting from expert knowledge, known physics models and inference capabilities, the calibration framework in¹¹² provided a very good trade-off between plausible/realistic reproduction of real dynamics, plausible choice of model parameters, and a very fast calibration procedure. From a methodology point of view, this approach is superior to those using solely optimization algorithms, as in the paper of Akasiadis and his colleagues⁶⁸ for instance. They only considered numerical optimization, where the goodness-of-fit only captured the quantitative aspects of the tumor growth agent-based model calibration. Yet, when considering the plausibility of the candidate solution, the numerical approach might offer a, sometimes, non-intuitive or plausible parametrization of the agent-based model behavior rules. This can be overcome by infusing a mathematical description of tumor growth covariates or other mechanistic dependencies in the learning function (i.e., modeling a loss function of the log-likelihood distance from the data to the mechanistic model). The different approaches for physics-informed learning calibration of ABM simulations demonstrate the potential such an approach has to leverage known models and learning algorithms and demonstrate how their combination is beneficial to achieve plausible, realistic simulations. Thus, physics-informed ML offers an attractive numerical procedure for extracting an accurate human-understandable representation of the underlying dynamics of physical interactions crucial to typical oncology problems, as demonstrated by the very encouraging results from multiple predictive tasks instantiations in oncology, as shown in the work of Kurz et al.¹¹³

This overview highlights the way ML systems may enhance clinical decision-making using effective computational techniques that benefit from embedding priors in the learning processes in order to guide their convergence toward plausible solutions. In order to do this, we think that such platforms provide a link between the modeler, the data scientist, the data, and the practising physician. In order to guide the reader, we offer a synthetic overview on the successful recipes for agent-based models’ calibration that fuse machine learning and mechanistic modeling in Table 4.

Table 4. Fusing mechanistic and machine learning approaches for ABM calibration; a synthetic overview on relevant state-of-the-art methodologies.

Approach	Reference(s)
Physics-informed Neural Networks	^Citation34,Citation105,Citation113
Bayesian Techniques	^Citation70,Citation71
Probabilistic Techniques	^Citation68,Citation106,Citation110,Citation111
Rule-based Techniques	^Citation68,Citation107,Citation108,Citation112
Hybrid Modelling Techniques	^Citation108–110

3. Strategies for agent-based models calibration and validation

3.1. Calibration as a multi-stage validation

Independent of the underlying parameter inference approach, a calibration procedure typically comprises multiple stages: data acquisition, scale choice (global or local), performance metrics (M) definition, the definition of the metrics for the goodness-offit (G), and the choice of an optimization algorithm. A first, and crucial question is what input data is available for an agent-based model calibration? In vitro cell line data are typically available, but poised by small sample size, uneven sampling, and multi-modality. On the other hand, in vivo data, although collected under strict (and thus more easily reproducible) protocols, often provides limited coverage of the study’s phenomena, focusing only on specific aspects. Although Equation 1 defines the calibration problem, we now exemplify how this is applied in practice. Consider the problem of phenotypical transitions in tumors. Due to their spatio-temporal modeling capabilities, ABM simulations can capture such dynamics. For instance, if we consider solid tumors, the typical cancer phenotypic state space is composed of three states quiescent (Q), proliferative (P) and apoptotic (A). Non-necrotic cells become hypoxic (H) when oxygen drops below a threshold value. Hypoxic cells can recover to their previous state or become necrotic (N). The transitions among these states are stochastic events generated by Poisson processes. Mapping this example onto Equation 1: β can be composed, for instance, of cells’ cycle time value, cells apoptosis time, proliferation index, and apoptosis index; F is a mapping function from β to states transition rates; x can be composed of the mean Q-P transition rates, and Q-A transition rates, f is a function of goodness-of-fit which can be, for instance, RMSE and symmetric Mean Absolute Percentual Error of F to the real timeseries of raw immuno-histochemistry and morphometric data of solid tumors. Second, the discrepancy in calibration is, typically, measured by the goodness-of-fit of the simulated parameters to the real parameters. This is supplemented by a series of calculated quantities by the metrics of performance in both the real world and simulation. Every goodness of fit evaluation can be performed either globally (e.g., using Least-squared errors) or locally (e.g., using maximum likelihood). Since at its core, model calibration is basically an optimization problem, the underlying algorithm aims to converge to a solution that is close to the global minimum of the goodness-of-fit metric, G, while obeying imposed constraints on parameters’ values. The formulation of the calibration procedure can be synthetically described as:

$\begin{array}{cc}\underset{\mathit{Msim},\mathit{Mreal}}{\mathit{min}}& \mathit{f}\left(\mathit{Mreal},\mathit{Msim}\right)\end{array}$

(1) $\mathrm{with}\mathit{f}\left(.\right)=\mathit{G}\left(.\right),M_{\mathit{sim}}=\mathit{F}(\mathit{x}|\mathit{\beta }),M_{\mathit{real}}=\mathit{F}\left(\mathit{x}\right)$(1)

$\mathrm{s}.\mathrm{t}.\mathit{\beta min}\le \mathit{\beta }\le \mathit{\beta max},$

where F is a mapping function of the biophysical system’s temporal trajectory data x provided a set of model parameters β calculated after simulating the model to extract the simulated values of the performance metrics M_sim, with M_real being the observed values of the performance metrics calculated by F from data x; f is a function representing the goodness-of-fit G (i.e., the realism of the simulation or closeness of M_real to M_sim), and f(M_real,M_sim) is the objective function to be optimized calculated from x – in principle a function describing the discrepancy between simulation and reality (i.e., M_sim ≈ M_real). A mapping of this formalism to the disease trajectory is illustrated in Figure 1(b). This formalism was adopted by many researchers in their attempts to calibrate ABM simulations. We chose to analyze two very relevant candidates that currently capture the state-of-the-art strengths and limitations of such methods. Cess and Finley⁶⁹ presented a novel approach that applied neural networks to represent both tumor images and ABM simulations as low dimensional points, with the distance between points acting as a quantitative measure of the difference between the two. This enabled the authors to extract a quantitative comparison of tumor images and ABM simulations, where the distance between simulated and experimental images can be minimized using standard parameter-fitting algorithms. We can see that each of the quantities in Equation 1(1) $\mathrm{with}\mathit{f}\left(.\right)=\mathit{G}\left(.\right),M_{\mathit{sim}}=\mathit{F}(\mathit{x}|\mathit{\beta }),M_{\mathit{real}}=\mathit{F}\left(\mathit{x}\right)$(1) can take arbitrary dimensions, but there will always be a ‘distance-based’ goodness-of-fit G to characterize the plausibility of the simulation after calibrating the models’ parameters β. But, as the literature survey shows, so far there has been no shared view on the quantification of validity in agent-based simulations. Very recently, Troost et al.¹¹⁴ conceptualized validation by systematically substantiating the premises on which conclusions from simulation analysis for a particular modeling context are built. They provided a formal extension to the classical approach in Equation 1(1) $\mathrm{with}\mathit{f}\left(.\right)=\mathit{G}\left(.\right),M_{\mathit{sim}}=\mathit{F}(\mathit{x}|\mathit{\beta }),M_{\mathit{real}}=\mathit{F}\left(\mathit{x}\right)$(1) . To this extent, a strict definition of the parameters β is the problem dependent and the choice of functions f depends on the trade-off between explanatory power, predictive accuracy, and the plausibility captured by the bounded β values. They proposed an assessment of the validity of agent-based models by incorporating valid conclusions from simulation analysis in a context-adequate method that touched model construction, model and parameter inference, uncertainty analysis, and the simulation process itself. The authors have built a comprehensive understanding of validation for agent-based models with a modular structure for versatile applications across several fields, where measures of performance are typically problem-specific. Introducing the framework multiple dimensions were defined to characterize the modeling context and purpose with a clear stage-wise buildup: premises of common validation approaches – a detailed protocol to guide context-adequate model construction and review and, finally, a consistent tracking of uncertainty propagation through the modeling process.

We now turn our attention to how optimization-based and optimization-free multistage methods complete the landscape of relevant approaches for multi-stage validation. In research targeting the estimation of the parameters of a stochastic process model for a macroparasite population within a host, the team of Drovandi et al.¹¹⁵ employed approximate Bayesian computation to model the immunity of the host as an unobserved model variable. Despite the very limited data, the authors had available, the process rate’s time constants were inferred reasonably precisely with a grounded plausibility proof. The approach involved a three-stage Markov process for which the observed data likelihood was computationally intractable. The proposed algorithm was validated on an autologistic model prior to parameters inference from experimental data. Interestingly, the model also captured the extra-binomial variation of the immune system. The results were also supported by the study of Carr et al.¹¹⁶ who presented a similar framework of Bayesian modeling and inference as in the work of Jorgensen et al.¹¹⁷ They proposed an efficient Bayesian inference method for a stochastic agent-based model. The study mitigated the use of the Bayesian setting (a) by constructing lightweight surrogate models to substitute the simulations used for inference, and (b) by circumventing the need for Bayesian sampling schemes and directly estimating the posterior distribution. This multi-staged approach demonstrated realistic results in tumor growth prediction. Considering a similar scenario of tumor growth curve extraction, the work of Wang et al.¹¹⁸ proposed a method for calibration of a Voronoi cell-based model for tumor growth using approximate Bayesian computation. Interestingly, the work involved as well estimating the distribution of parameters that govern cancer cell proliferation (i.e., the distribution of β in Equation 1(1) $\mathrm{with}\mathit{f}\left(.\right)=\mathit{G}\left(.\right),M_{\mathit{sim}}=\mathit{F}(\mathit{x}|\mathit{\beta }),M_{\mathit{real}}=\mathit{F}\left(\mathit{x}\right)$(1) ) and recovering outputs that match the experimental data. Their results showed that the proposed approach, and its multi-stage extension, provided insights into tumor growth and a good quantification of this process uncertainty. Multi-stage calibration describes a very promising avenue to explore, also because it is also supported by work that fuses mechanistic modeling and ML. A very good candidate subsuming these principles is the work of Axenie et al.,¹¹² where the optimization step was reduced to a simple feed-forward inference through a logic model of spatio-temporal interactions of the agents. This modeling stage allowed the system to ‘inject’ proper biases in the learning process of the mapping F in the generic process captured by Equation 1(1) $\mathrm{with}\mathit{f}\left(.\right)=\mathit{G}\left(.\right),M_{\mathit{sim}}=\mathit{F}(\mathit{x}|\mathit{\beta }),M_{\mathit{real}}=\mathit{F}\left(\mathit{x}\right)$(1) while keeping inference efficient. Additionally, the employed metrics of performance demonstrated a tight coupling among the micro- and macro-scopic dynamics of the agents. Interestingly, the calibration of microscopic parameters took into account also the aforementioned scale coupling for a plausible candidate parameter configuration. Table 5 summarizes the most relevant publications that categorize the most important modeling approaches for ABM calibration as a multi-stage validation process.

Table 5. Calibration as a multi-stage validation. An overview on relevant state-of-theart methodologies that follow the formal description in Equation 1(1) $\mathrm{with}\mathit{f}\left(.\right)=\mathit{G}\left(.\right),M_{\mathit{sim}}=\mathit{F}(\mathit{x}|\mathit{\beta }),M_{\mathit{real}}=\mathit{F}\left(\mathit{x}\right)$(1) .

Approach	Reference(s)
Neural Networks	^Citation69,Citation112
Generic Optimisation Techniques	^Citation114
Bayesian Inference	^Citation114–118

3.2. Models comparison — benchmarking

Given the availability of different computational techniques and approaches, the scientific community acknowledges the need to assess and compare different explanatory computational models. In particular, this abundance of in silico models gives rise to the question of what criteria investigators should consider when making decisions of appropriateness. One commonly used process to address this question is model benchmarking. Model benchmarking stands for the assessment and quantification of in silico models according to well-specified criteria. This provides the means to objectively infer metrics and uses these to compare the models, enabling a framework for improved selection.

Notably, without proper standardization the benchmarking criteria can be strongly influenced by the research aims of the individual investigator(s), institutions or national funding policies, thus, causing biases and misconceptions. Here, we elaborate on the set of criteria that we believe to be crucial for ABM simulations in biomedicine and hence well-suited for benchmarking studies in the future.

Naturally, computational models’ primary objective is to achieve a high-fidelity predictive performance. To accomplish this, specific measured data need to be reproduced by the in silico model in a repeatable fashion (‘replicability’). The capability in doing so is quantified using measures of ‘accuracy.’ In the biomedical context, such actions could relate to various anatomical, physiological, omics or other types of biological information. For instance, the study of Borlin et al.⁹⁵ presented an agent-based model that quantitatively captures several biological measurements of autophagy, which plays crucial roles in cellular and organismal homeostasis, including the response to diverse stresses. Ideally, the model performance refers to some type of agreement with such experimental data; for instance, this could be the number of cells or subcellular vesicles of a certain type. Well established measures, such as Mean Squared Error (MSE), RMSE, the F-score and the Receiver Operating Characteristic (ROC) curve, are commonly used to quantify the agreement between computational predictions and experimentally observed data. Nevertheless, oftentimes such agreement may be only of a qualitative nature, which can still exhibit significant explanatory power. For instance, a computational model that captures up- or down-regulation of certain metabolic pathways in the right context, has value, even if the magnitude of change is not quantitatively accounted for. However, in the presence of multiple plausible models, quantitative accuracy is a crucial factor and may determine which models to select over others.

As previously discussed (section 2.4), computational models and in particular mechanistic models require parametrization using biological information. These should be ideally experimentally measured, and at least based upon plausible evidence. Given that such parametrization may be difficult to conduct based on real-world data, it is best to minimize the presence of estimated model parameters while also reducing to the absolute necessary ones, based on the design of the agent-based model, the list of assumptions. Arguably, models that require fewer parameters and less guesswork can generally be considered superior to those where significant subjective inputs are necessary. When multiple models confer equivalent agreement with experimental data, models with comparably reduced ‘complexity’ are preferable to others, in accordance with Occam’s Razor.¹¹⁹ Thus, with regard to mode benchmarking, it is important to consider and compare model complexities. At the very least, model complexity should reflect the complexity and disparity of the empirical data used to inform and calibrate the in silico model. It is beyond the scope of this review, however, to give a comprehensive overview of model complexity measures, which comprise for instance the number of model parameters, number of agents, the ontology graph, the Bayesian Information Criterion, and many others. The interested reader is recommended however to read.^119–122

Explanatory power is, as mentioned earlier, a particularly important aspect of an in silico model in biomedicine. A hypothetical model A, which is considered in this example as a black-box one, may perform significantly better in reproducing experimental data than a hypothetical mechanistic model B. However, the real-world value of a model B may still be higher. No patient will agree to surgery without a proper explanation of the underlying reasoning. Indeed, model A’s applicability is a-priori limited in contexts where the stakes are high and interventions can have long-term implications. Nevertheless, to date, the explainable category of a model B has seen very limited translation and thus application in the clinical setting has been rudimentary.

However, explanatory power is also difficult to quantify and so remains rarely assessed in existing benchmarking studies. A potential framework to assess and grade a model’s explainability is as follows: in its most basic form, a model could be associated with a binary variable for being of the black-box or white-box type. Black-box models are entirely data-driven, while white-box models constitute theoretical models that describe the mechanisms of a given system. In many cases, this basic differentiation could be extended to a spectrum of explainability, associated with gray-box models. Pure ABM is a mechanistic approach, hence it is a white-box type. Given that models exist where ABM interacts with statistical approaches (e.g., Monte Carlo methods¹⁷), it is conceivable that a more fine-grained explainability measure can be constructed. Explainability is an active area of artificial intelligence research and is particularly important to consider when assessing the usage of computational models in biomedical applications. Explainability needs to be accounted for in accordance with experimental data across multiple spatial and temporal scales. Again, the extent to which different scales are validated can be associated with explainability.¹²³ For instance, the study of Bauer et al.¹³ combined ABM of single-cell behavior in terms of cell death, with developmental outcomes in terms of large-scale tissue anatomy. The larger the range of spatial and temporal scales that an agent-based model allows to capture and reproduce, the higher its explanatory power. Along those lines, Wijeratne and Vavourakis⁷ studied the impact of drug-borne nanoparticles on multiple levels including the fluid dynamics, biochemics, vasculature, and overall solid tumor, while Ruscone et al.⁸⁶ also proposed an agent-based model and probe dynamics on the scales of genetics, metabolism, cell behavior and population.

Biomedical applications are in many regards very challenging, particularly when it comes to personalized medicine. Many diseases can progress very differently depending on lifestyle, age, genetics and comorbidities. In order for a computational model to have predictive power, it requires the capability to account for variation as observed in a given system (e.g., interpersonal variation, age-dependent changes, drug-induced changes, etc.). This aspect could therefore be accounted for in the benchmarking of computational models, for instance by producing distributions of outputs for multiple runs with different initial settings and/or stochastic dynamics. At the same time, computational models also need to be consistent, so that a given result depends on very specific initial conditions. Such variation and consistency can be captured using ‘sensitivity’ analysis.^28,66,124

It is not astounding that most bench-marking studies focus on computational performance, which depends on ‘computational efficacy.’ Needless to say, the more efficiently and faster a model can be simulated, the better suited it is for parameter optimization, refinement, and adaptation. While this does not necessarily entail enhanced explanatory power, it can lead to an edge when it comes to real-world applications. This is particularly true in the biomedical context where data acquisition has reached very high throughput levels, requiring frequent analysis, fast processing, and repeated updating and adaptation of computational models. Advanced code optimization techniques can significantly reduce computational time.¹²⁵ In addition to run-time, hardware resource demands that the model implementation imposes, such as usage of memory and energy, are factors to be considered. Moreover, the volume of code (lines of code) and its complexity¹²⁶ play significant roles when assessing computational models. Moreover, ML methods often out-compete mechanistic, multi-scale models in terms of computational efficacy.¹²⁷ Therefore, the incorporation of techniques ML can be advantageous for certain applications, such as for the prognosis of patients who may suffer from time-critical issues.¹²⁸

To conclude this section, several factors that are relevant for comparing computational models are outlined in Table 6 – the table is not meant however to provide an exhaustive list of the relevant references. These factors can serve as a stepping stone to formulate an ecosystem comprising the multitude of computational models, and establish their strengths and weaknesses in a systematic benchmarking approach. Naturally, the factors need to be considered on a case-by-case basis, and their relative importance varies according to context and the stakes at hand.

Table 6. Benchmarking multi-scale computational models. Here, important criteria for benchmarking are shown, along with some example quantities that can be objectively measured and relevant references.

Criterion	Quantities	Reference(s)
Accuracy	MSE, RMSE, F-score, ROC curve	^Citation28,Citation95
Complexity	# of model parameters, # of agents, connectivity density of ontology graph, Bayesian Information Criterion	^Citation119–122
Explanatory power	black box vs. white box model, accordance with experimental data across scales	^Citation7,Citation13,Citation17,Citation86
Model robustness/sensitivity	coefficient of variation, sensitivity indices	^Citation28,Citation66,Citation124
Computational efficacy	lines of code, run-time, memory usage, energy demand	^Citation125–128

4. Discussion and outlook

4.1. Modern mechanistic modelling approaches in cancer biomedicine

Agent-based models have proven to be invaluable tools for providing new insights and testing hypotheses on complex heterogeneous systems, including biological systems. In the context of cancer biomedicine, this novel mechanistic approach offers a unique way (a) to fuse information from in vitro and in vivo data including for instance gene expression and medical imaging,²⁹ (b) to develop novel in silico procedures for patient diagnosis and neoplasia stratification,³⁹ (c) to identify and compare the efficacy of potential treatment strategies⁵⁹ as well as to narrow the cardinality of the set of possible experiments down to the most promising ones while also (d) to aid reducing animal testing and clinical trials on humans.⁴⁷

Recent efforts aimed at employing ABM to model tumor tissues have been presented in subsection 2.1, which unveils the versatility of the approach. In fact, at multiple spatial and temporal scales, the models examined have demonstrated how various aspects of cancer dynamics can be captured, spanning from the pathogenesis to the growth and interaction with somatic cells. Furthermore, ABM offers a versatile framework that can be seamlessly integrated with other modeling approaches, facilitating validation against large-scale experimental datasets. ABM’s granular representation of individual entities and their interactions allows for the incorporation of detailed biological mechanisms, capturing emergent properties and system-level behaviors. This flexibility enables ABM to complement both traditional modeling techniques such as equation-based models which excel in capturing system-level dynamics and spatial gradients, and more recent data-driven approaches (see, for example, the excellent review by Metzcar et al.¹²⁹). This integrative approach offers mutual benefits, as ABM provides detailed mechanistic insights that inform the development and refinement of equation-based models, while equation-based models provide a framework for scaling up ABM simulations to larger spatial and temporal scales. By combining the strengths of both approaches, investigators can construct hybrid models that bridge the gap between micro- and macroscopic phenomena, enabling a more comprehensive understanding of system-level dynamics and facilitating the design of targeted therapeutic interventions. Hence, agent-based models are excellently suited to produce explanatory powerful computational models that can integrate various data modalities, heterogeneous information as well as background information that may not be readily available in the form of experimental measurements.

4.2. Calibration ‘recipes’ for validating the simulations of agent-based models

Calibration is the methodological procedure to compute model parameter configurations when the in silico model can produce realistic simulation outputs. This process can be approached from multiple perspectives and there is a plethora of tools capable to optimize the parameters and, thus, to improve the fidelity of large-scale agent-based models in cancer biomedicine. Independent of their inner workings, the calibration process is a ‘search’ in a large parameter space. This search can be guided (as in the case of ML-based approaches) by embedding information about the physics of the process being modeled or rules of the evolution of the process dynamics. While most of the calibration approaches formulate the ‘search’ as an optimization method (i.e., to minimize a distance among some performance indices), there is much more to gain from embedding rules and physics in the optimization process. This has been formally described in subsection 2.5, where a comparative benchmark has been developed to offer the readers a ‘recipe’ for approaching the validation of agent-based models through calibration. However, we acknowledge that there is no ‘one size fits all’ solution, and we thus recommend formulating the calibration process as suggested in subsection 3.1 followed by a benchmarking as suggested in subsection 3.2 and in Table 6. Besides, choosing the best fitting candidate from the “models zoo” is thoroughly covered in subsection 2.2 and in the excellent review of Metzcar and his colleagues.² Starting there the reader could get a grasp on the necessary functionality and capabilities of already successful simulation tools and machine learning and physics-informed learning algorithms in cancer biomedicine. Notably, it is currently not viable to quantitatively compare different ABM studies due to the model heterogeneity and the varying experimental data used for validation. We note that increasing usage of common experimental datasets and computational frameworks will help establish ways to compare agent-based models. Such benchmarks could then be used to justify the usage of specific ABM works over others, which will ultimately help bridge the gap toward clinical application.¹³⁰ The models we described in our review are mainly preclinical, computational and biological. It has been very well emphasized by Markowetz¹³⁰ that there is a wealth of models excelling in preclinical setups, but out of these only a few have ever been independently validated, let alone been adopted in the clinic. As we rolled out an interdisciplinary ABM-based tools overview that applies to cancer biomedicine, we still propose and encourage: (a) model developers to have a roadmap to clinical practice, (b) provide a model implementation and regulation plan to prepare the model for the clinic, and (c) regularly interact and exchange ideas with clinicians to stay abreast of their views and new developments in clinical practice.

4.3. Outlook to ‘non-invasive’ calibration of agent-based models

While a computer simulation is a cost-effective and safe way to evaluate hypotheses on new cancer therapies, immune – tumor interactions, and bio-physics of tumor invasion, existing simulators are often impractical due to inefficient and/or impractical control interfaces. This hinders the data exchange of the simulation evolution with the ‘external’ world where signals, data, or decisions to update the temporal trajectory of the simulation might be applied. A mature simulator, under the control of external systems, can be used to run the models being tested within a closed-loop environment. However, state-of-the-art simulators are often impractical for the task at hand due to their inefficient control interfaces. The challenges can be divided into two primary areas: the requirement for efficient synchronization, and the need for flexible and scalable data processing during runtime. In terms of the synchronization model, state-of-the-art ABM simulators primarily use a conservative step-based approach, which needs synchronization at every interaction point. However, the high frequency of interactions with the simulation process caused by data collection (i.e. probing) can result in an overwhelming amount of message flow, leading to significant overhead in the simulation execution. Alternatively, an optimistic approach that allows for full asynchronicity and causality violations exists, but implementing a rollback mechanism to address missing decision-making or data collection times is complex, error-prone, and not supported by most ABM simulators. Processing simulation data during runtime poses other challenges, e.g., storage and large volumes of simulation outputs, secure and fast transfer of sensitive clinical data, extracting valuable insights from a limited data amount, making predictions based on an unevenly distributed sample etc. The objective is to dynamically instruct the running simulation, allowing it to store and output the data of interest while avoiding any redundancy. An essential yet challenging aspect of this problem is performing temporal operations and on-the-fly data processing within queries, for example, filtering or computing data to retrieve results over a specific time interval.

List of Acronyms

ABM	=	Agent-Based Modeling
CAR	=	Chimeric Antigen Receptor
CTLs	=	Cytotoxic T Lymphocytes
ECM	=	Extracellular Matrix
FEM	=	Finite Element Method
GBM	=	Glioblastoma
LHS	=	Latin Hypercube Sampling
ML	=	Machine Learning
ML-MM	=	Machine Learning-assisted Mechanistic Modeling
MM	=	Mechanistic Modeling
MMPs	=	Matrix Metalloproteineases
MSE	=	Mean Squared Error
ODE	=	Ordinary Differential Equation
PDE	=	Partial Differential Equation
QSP	=	Quantitative Systems Pharmacology model
RMSE	=	Root Mean Squared Error
ROC	=	Receiver Operating Characteristic
TME	=	Tumor Microenvironment
TNF	=	Tumor Necrosis Factor

Acknowledgments

V.V. acknowledges the financial support from the Cyprus Cancer Research Institute for the project “PROTEAS” (Call: CCRI Bridges in Research Excellence; Funding Agreement No. CCRI 2021 FA LE 105), and the European Innovation Council for the project “IgnitePLASMA” (Call: HORIZON-EIC-2023-PATHFINDEROPEN-01; Grant agreement ID: 101129853).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

R.B. was supported by the Engineering and Physical Sciences Research Council of the UK (grant no. EP/S001433/2). V.V. was supported by the Cyprus Cancer Research Institute [CCRI_2021_FA_LE_105] and the European Innovation Council [101129853].