1. Introduction

2. Materials and Methods

2.1. Component Models

2.1.1. The Oncosimulator

The Oncosimulator is a lattice-based, discrete-event discrete-state approach that models tumor cell population kinetics at supercellular and tissue scale either under free growth or under treatment conditions. Within the framework of CHIC, two instances of the Oncosimulator concept have been implemented: the Lung Oncosimulator, a model of lung tumor response to external beam radiotherapy, and the Wilms tumor, WT, (Nephroblastoma) Oncosimulator, a model of Wilms tumor response to preoperative combined chemotherapy treatment of actinomycin and vincristine. The core algorithms of the Lung and WT Oncosimulators have been previously developed by the In Silico Oncology and In Silico Medicine Group and the implementation details can be found in [3,8]. Herein, Lung Oncosimulator was algorithmically extended to account for the effect of external beam radiotherapy as described in [9]. To enable communication with the CHIC platform and the exchange of data with the other component models, new code was developed and added according to guidelines [10] in both instances [see also section ‘Hypermodel coupling topology’]. Below, basic notions of the modelling approach, common for both Oncosimulators, are briefly presented.

The reconstructed tumor area is represented by a three-dimensional grid of cubic voxels, named Geometrical Cells (GCs). Each GC belonging to the tumor corresponds to a tissue volume of either 1 or 8 mm³ and is occupied by an inhomogeneous community of living and dead tumor cells. Specifically, the tumor is assumed to be organized as a hierarchy originating from a type of immature cells with unlimited mitotic capacity. These cells, termed cancer stem cells, may divide either symmetrically, with probability P_sym, to produce two stem cells or asymmetrically to produce a stem cell and a cell of limited mitotic potential (LIMP) that follows an aberrant differentiation process. LIMP cells are allowed to perform a specific number of divisions, N_LIMP, before entering an irreversible differentiated state (compartment of differentiated -DIFF- cells). Stem and LIMP cells may exist in a cycling or resting, G0, phase. Cycling cancer cells are distributed into four compartments corresponding to the four cell cycle phases (G1, S, G2, M). The withdraw probability, P_sleep, of cycling cells to a resting phase following mitosis, is regulated by the local conditions of nutrient and oxygen supply. Tumor cells are allowed to spend an average time, T_G0, in G0 phase. Afterwards they re-enter the cell cycle, with probability P_G0toG1, or die via necrosis. The necrotic loss of resting cells is assumed to be caused by nutrient or oxygen deprivation. Stem and LIMP cell categories may die with rate R_A, through spontaneous apoptosis. DIFF cells may undergo either apoptosis with rate R_ADiff or necrosis with rate R_NDiff. The time required for apoptotic and necrotic cells to be permanently removed from the tumor bulk is T_A and T_N respectively.

The Oncosimulator explicitly models chemotherapy and/or radiotherapy. During chemotherapeutic treatment, a fraction of stem and LIMP cells are assumed to undergo lethal damage by the drug(s). These cells follow a rudimentary cell cycle before apoptotic death through a cell cycle phase dictated each time by the mechanism of action of the specific chemotherapeutic agent. The effect of the drug is assumed instantaneous at the time of its administration. In the case of radiation therapy lethally damaged cells die through a radiation-induced necrotic mechanism. These cells enter a rudimentary cell cycle and die after undergoing a few mitotic divisions. The probability of cells to be hit by irradiation depends primarily on the phase they reside. Cell killing by irradiation is described by the Linear Quadratic or LQ Model [11]:

S(D)=exp[-(αD+βD²)],

(1)

where S(D) is the surviving fraction after a (uniform) dose D (Gy) of radiation to a population of cells. The parameters α (alpha) (Gy^-1) and β (beta) (Gy^-2) are called the radiosensitivity parameters of the LQ model.

The model integrates cytokinetic, metabolic, pharmacokinetic/pharmacodynamic and mechanical rules to simulate the dynamic behavior of the tumor over time. A computational grid (discretized mesh) is used to model the spatial distribution of cells within the tumor. Two mesh scans address the biological aspects (cytokinetics) and the spatial dynamics (tumor expansion/shrinkage) within the simulation. Specially-designed algorithms are employed to control the movement of cells within the mesh, defining tumor spatial evolution as described analytically in [3]. In the framework of the hypermodel presented herein, a biomechanical simulator governs the movement of biological cells within the discretized mesh, as detailed in the subsequent sections.

The model is implemented in C++ programming language.

Table 1. Model parameters of the Oncosimulator.

Table 1. Model parameters of the Oncosimulator.

Symbol	Description	Units
CELL PHASE DURATIONS
Tc	Cell cycle duration	hr
TG0	G0 (dormant phase) duration i.e., time interval before a dormant cell re-enters cell cycle or dies through necrosis	hr
TA	Time needed for both apoptosis to be completed and its products to be removed from the tumor	hr
TN	Time needed for both necrosis to be completed and its lysis products to be removed from the tumor	hr
CELL CATEGORY/PHASE TRANSITION RATES AND FRACTIONS¦
RA	Apoptosis rate of living stem and limp tumor cells, i.e., fraction of cells dying through apoptosis per unit time	hr -1
RADiff	Apoptosis rate of differentiated tumor cells	hr -1
RNDiff	Necrosis rate of differentiated tumor cells	hr -1
Psym	Fraction of stem cells that perform symmetric division	-
Psleep	Fraction of cells entering the G0 phase following mitosis	-
PG0toG1	Fraction of dormant (stem and LIMP) cells that re-enter cell cycle	-
MISCELLANEOUS PARAMETERS
NLIMP	Number of mitoses performed by LIMP cells before becoming differentiated	-
RADIOTHERAPY - CHEMOTHERAPY PARAMETERS
CKR	Cell kill rate: the numbers of biological cells lethally hit by treatment at each administration. In case of chemo, it is defined separately for each drug.	-
α/β	alpha to beta ratio	Gy
D	Dose of radiation to a population of cells	Gy

2.1.2. The Molecular Hypomodel

The molecular component explicitly models two types of signalling pathways, which are important determinants of tumor cell fate (death and proliferation) and treatment resistance, as well as the interfaces between them: the ErbB receptor mediated Ras-MAPK and PI3K-AKT pathways and the TP53 mediated DNA damage response pathways.

ErbB Receptor Mediated Ras/Raf/MAPK and PI3K/AKT Pathways: This part of the model has been adapted from Chen et al. [12] after suitable modifications. It is a continuum ordinary differential equations (ODE) based model with 504 distinct species, 827 elementary mass action type reactions, and 252 parameters. It consists of all the ErbB family of receptors and a subset of the homo and heterodimers. Though the model includes both epidermal growth factor (EGF) and heregulin (HRG) as the growth factors, for simplicity we only consider the effect of EGF here. The growth factors activate the receptors which, in turn, initialize the downstream signaling cascade. This component model incorporates the effect of receptor internalization and recycling and the cross-talks involved between the Ras-MAPK cascade and PI3K/AKT cascade. This has also been extended to take into account various mutant forms of EGFR receptor like L858R and deletion mutants which can be constitutively active.

TP53 Mediated DNA Damage Response Module: This part of the model has been adapted and modified from [13]. It consists of two main submodules corresponding to cell cycle progression and apoptosis. The first module consists of cyclin-CDK mediated cell cycle progression which affects the G1 to S phase entry of the cell cycle. The cell death module is the intrinsic apoptotic pathway which is mediated by key proteins like Bax and Bcl-2 which regulates the cell death protein Caspase 9. Damage to DNA activates ATM kinase which in turn activates TP53 through kinases Chk1 and Chk2. Activated TP53 in turn activates DNA damage repair pathways or cell death pathways depending on the extent of damage. The ultimate cell fate will depend on the combined interactions of all the various components of the network and their initial activation state. This module consists of 16 nodes with 160 negative and 218 positive feedbacks. The module is modeled using a discrete Boolean model. Depending on set thresholds, each node of the network can have two possible states – ON or OFF. The interaction between the nodes is also a discrete number which can be both positive and negative depending on whether it activates or represses the downstream node. These kinds of discrete models can give two possible outcomes: a) a point attractor which is a single steady state where the activation state of all the nodes in the network does not change over successive time steps and b) a cyclic attractor which corresponds to a sequence of repeating states (cycle). For the current model there were three possible outcomes corresponding to three different cell fates: a) cell cycle progression (point attractor with high cyclin-G and low p53 activity), b) apoptosis (point attractor high p53 and high caspase activity) and c) cell senescence (cyclic attractor with oscillations in p53 and Mdm2).

Module Interfaces and Hybrid Simulator Algorithm: The two modules are able to communicate through the states of the common nodes; these are Erk, Akt and PTEN. We run the modules sequentially where the final states of the interface nodes obtained from each module is fed to the other module at the start of each new time step. The ODE states of the common nodes are described by continuous time concentration functions which are discretized by applying appropriate thresholds before passing them to the Boolean module. The Boolean p53 model will pass the activation fraction of the common nodes to the ODE Ras-MAPK and PI3K/AKT module. We assume that the reactions in the Ras-MAPK and PI3K/AKT module are much faster than the p53 module. This enables us to partially uncouple these processes and pass pseudo-steady state information from the fast to the slow process. We assume that within the short time step of the Ras-MAPK and PI3K/AKT module, the state of p53 module is invariant. On the other hand, the p53 module will evolve with its own time scale but its behavior will be modified by the information about the interface species it receives from Ras-MAPK and PI3K/AKT module. The modules are run until the p53 module (slow process) converge to a steady state (point or cyclic attractor).

2.1.3. The Vasculature Hypomodel

The vasculature hypomodel describes the transport of nutrients, glucose in particular, in tumors at the tissue-scale. It uses the finite difference method to predict nutrient concentrations as a function of the three-dimensional spatial distribution of vessel volume fraction and tumor tissue, the latter being provided by other hypomodels. The vessel volume fraction is assumed fixed in time. The vasculature hypomodel returns a three-dimensional nutrient field, which can be used as input to other hypomodels.

In general, the vasculature plays a vital role in the transport of nutrients and therapeutics to tumors, with tumor size limited by its ability to co-opt and maintain a vessel network. The Angiogenesis hypomodel is motivated by the well-known model of vascularized tumor growth by Hahnfeldt et al. [14]. This model describes the rate of change of tumor volume $T$ as a function of carrying capacity $K$, assuming a spherical tumor:

$$T^{\prime }=-\alpha Tlog\left(\frac{T}{K}\right)$$

(2)

where $\alpha$ is a growth rate parameter. The adopted rate of change of the carrying capacity term is:

$$K^{\prime }=-{\alpha }_{2}K+bT-dK{T}^{\frac{2}{3}}$$

(3)

where ${\alpha }_2$, $b$ and $d$ are constants. Implicit in this form is a balance between a rate of increase in carrying capacity due to tumor stimulation of new vasculature, and a rate of decrease in carrying capacity due to increasing diffusion length scales as the tumor grows. Since 3D imaging data and other hypomodels can give a more general tumor shape as inputs to the hypomodel, it is useful to relax the spherical tumor assumption.

The vasculature hypomodel assumes steady-state, diffusion-limited transport of nutrient with concentration $c$, which is supplied by the vasculature at a rate dependent on vessel amount (volume fraction or density) $V$, and is consumed by tumor tissue at a rate proportional to the number or volume fraction of viable cells $P$. The tissue is assumed to comprise a tumor region and a non-tumor region, as shown in Figure 1.

In practice these regions are determined based on segmentations of clinical images, performed in a pre-processing step external to the vasculature component in the hypermodel execution. In the non-tumor region it is assumed that the tissue is well vascularized and nutrient concentrations are set to a reference value $c_n$. In the tumor region nutrient transport is described according to:

$$D{\nabla }^{2}c-\lambda Pc+\rho V\left(c_n-c\right)=0$$

(4)

where $D$ is the effective nutrient diffusion coefficient in the tumor tissue, $\lambda$ (redefined) is the rate of nutrient consumption and $\rho$is rate of nutrient delivery by vessels. In this simple model the nutrient concentration in the tumor will approach the value in the surrounding healthy tissue as the vessel volume fraction increases, or as the rate of consumption by cells decreases.

For the WT and NSCLC hypermodels the Vasculature component describes the transport of glucose and uses a vessel volume fraction that is fixed in time. This simple model was chosen to aid hypermodel integration and validation, as it has a favourably low number of input parameters.

The nutrient transport problem is solved on a regular finite difference grid in 3D. This method was chosen for computational efficiency and to avoid interpolation when used with the grid-based descriptions of cell growth used in other hypomodels.

The model is implemented in the Chaste [15] open-source C++ framework for soft tissue modelling. A custom Chaste build was developed to allow the incorporation of MUSCLE libraries for run-time coupling of hypomodels and also packaging as a standalone executable. The code can be found on GitHub [16]. Validity checks of the vasculature component are presented in Appendix B.

2.1.4. The Metabolic Hypomodel

Normal mammalian cells are exposed to a continuous supply of oxygen, glucose and other nutrients in circulating blood. In normal cells, glucose is taken up by specific transporters and is converted to pyruvate in the cytoplasm through glycolysis generating 2 moles of ATP per glucose. In the presence of sufficient oxygen, pyruvate is then completely oxidized in mitochondria generating additional 36 moles of ATP per glucose. However, when oxygen is insufficient, pyruvate is redirected away from mitochondrial oxidation and is converted to the waste product lactate. In contrast to normal cell metabolism, Warburg’s observations [17] showed that cancer cells produce a substantial amount of energy by inefficiently metabolizing glucose to lactate, independent of oxygen availability- a phenomenon termed the Warburg effect or aerobic glycolysis. The exact regulatory mechanisms of tumor metabolism are far from complete. The tumor microenvironment significantly affects the metabolic activity and rewiring, impacting metabolite transporters and glycolytic enzymes. Signaling pathways involving various oncogenes and tumor suppressor genes have been identified to play a role in the altered metabolism [18,19].

To model the metabolic adaptations of highly proliferating human cancer cells, Shlomi et al. [4] employed a genome-scale human metabolic network, comprising 1496 ORFs, 3742 reactions, and 2766 metabolites [20]. They introduced metabolic demands for biomass synthesis required for high proliferation rates, with the flux through biomass serving as the objective function. Additionally, they considered solvent capacity constraints to further restrict the fluxes of metabolic reactions. Their approach successfully replicated several experimentally observed metabolic characteristics during cancer development.

We extend the work of Shlomi et al. [4] and apply a metabolic strategy that allows for near-optimal growth solution, while maximizing lactate secretion. This approach aims to elucidate the high-flux mechanisms leading to a substantial increase in lactate production observed in tumor cells. Sub-optimal growth solutions have been noted to describe the metabolic capabilities of microorganisms under environmental stress and in the absence of sufficient evolutionary pressure [21] indicating that variability around optimal growth is not unexpected for biological systems, including cancer. The lactate maximization strategy is mathematically described as a two-step optimization problem, akin to the Flux Variability Analysis (FVA) method [22], which has been utilized to identify alternate optimal and sub-optimal metabolic states. The lactate maximization strategy employs an iterative procedure to pinpoint the minimal compromise in growth rate necessary to achieve lactate production [23].

Specifically, in the first step, the optimization problem is solved as described in [4] where cells are assumed to maximize their growth rate subject to flux balancing constraints, uptake bounds in the substrate reactions and the solvent capacity constraint. Within these specific constraints, the first problem aims to identify the maximum growth rate. The second optimization problem aims to maximize the lactate production rate. Additionally, it includes the constraint that the growth rate should not be less than a given percentage, k, of the optimal growth rate determined in the first problem. The second step is iteratively performed for decreasing values of k until a solution is found, as long as the lactate secretion rate remains below a specified tolerance (0.01 umol/mgDW/h). The model provides a solution closer to optimal growth by varying k from its maximum to lower values. The model has demonstrated its ability to capture several metabolic phenotypes observed experimentally in cancer. Slight deviations around the optimal growth rate (90-99%) were found to be sufficient for adequate lactate production, with increasing deviations observed at lower glucose uptake bounds.

2.1.5. The Biomechanics Simulator

The Biomechanics Simulator (BMS) aims to predict the mechanical impact of a growing tumor, its so-called “mass-effect”. Rapid cell division and increase of the number of cells in a tissue segment gives rise to mechanical forces that lead to volumetric growth. This tumor-induced strain results in mechanical stresses in the tumor and the surrounding healthy tissue.

Such tumor-induced biomechanical forces shape the tumor environment and are known to affect tumor growth and evolution [24], for example by reducing blood perfusion through compression of intratumoral vessels [25]. For certain tumor locations, such as the brain, mass-effect is also of direct clinical importance.

BMS is designed to model a tumor’s mass-effect and the resulting mechanical stress distribution on macroscopic length scales, i.e., the tissue-level. First, tumor-induced strains are computed as spatially varying volumetric growth factor ${ϵ}_{growth}\left(\mathit{x}\right)$, where $\mathit{x}$ indicates the location in space, based on volumetric considerations and the ratio of local tumor cell concentration to a reference concentration $c_0$. Incorporating these strains into a continuum mechanics model of the tumor growth domain then allows the resulting stresses to be computed. We have previously employed this approach for modelling the mechanical stresses resulting from tumor growth [26,27,28,29].

In the CHIC hypermodelling framework, mechanical information is used to identify the most likely growth direction of the tumor. Thus, the role of BMS is to compute tissue stresses resulting from the current tumor cell distribution in each time step. From this information, the direction for tumor growth and shrinkage are computed to inform the redistribution of tumor cells in the next time step.

BMS relies on the Finite Element Method (FEM) to solve the linear-momentum equilibrium equations with a given mechanical material model and in a patient-specific anatomy. Its usage involves a pre-processing step in which a personalized FE Model of the patient-specific anatomy is created and parametrized, followed by iteratively coupled execution with OS, as described in section 2.2.3. A custom pre-processing pipeline has been developed to automate the model configuration process, including the assignment of material properties and boundary conditions from simple configuration options. In combination with automatic segmentation tools, this pipeline permits rapid generation of patient-specific FEM models for personalized simulations (section 2.3.5). FEM model and pre-processing pipeline are implemented using Open-Source libraries and software packages (CGAL, VTK, FEBio).

2.2. Hypermodel Coupling Topology

2.2.1. Oncosimulator – Molecular Hypomodel

Cellular intrinsic sensitivity or resistance to treatment is a determinant of treatment outcome. The cell kill probability (CKP), i.e., the probability of a tumor cell to be lethally hit by a given therapeutic regimen is an input parameter of the Oncosimulator (Table 1) and it is explicitly computed by the Molecular model based on the molecular profile (e.g., EGFR mutations, miRNA expression data etc.) of the patient.

The unidirectional data flow between the Oncosimulator and the molecular model has been implemented by a serial coupling topology, i.e., the component models are called and executed sequentially. Within CHIC framework, the aforementioned coupling topology is orchestrated via TAVERNA workflow management system [30] which passes CKP as a command line argument to the Oncosimulator (Figure 2).

2.2.2. Oncosimulator – Vasculature – Metabolic Hypomodels

The hyper-modelling scenario dictates an iterative, tightly-coupled communication scheme between the Oncosimulator, Vasculature and Metabolic models (Figure 2).

As tumor grows, well-vascularized regions that provide ample nutrients to cancer cells may coexist with nutrient-limited areas within the tumor mass. In this work, we focus on glucose as the sole limiting resource, although oxygen and glutamine can also be considered. The dependence of glucose uptake on glucose concentration is modeled using Michaelis-Menten kinetics as illustrated in equation (5). Here, C represents the glucose concentration, V_max corresponds to the maximum rate of the process (dependent on factors like GLUT receptor concentration), and K_m is the saturation constant indicating the glucose concentration at which the uptake rate equals to V_max/2. This implies that glucose availability sets an upper limit on the glucose uptake rate. We further assume that V_bound reaches V_max when the glucose concentration C reaches its maximum observed value in tissues (C_max= 0.9 kg m^-3).

$$v_{bound}=v_{max}\frac{C}{K_m+C},$$

(5)

A relatively slow varying environment is assumed where cancer cells can operate at optimal or near optimal growth rates constrained by the current nutrient availability. The Oncosimulator computes and passes the initial and updated (during execution) tumor domain geometry and population of proliferating, quiescent, terminally differentiated, apoptotic and necrotic tumor cells at each voxel of the grid (GC-geometrical cell) to the vasculature model. The vasculature model solves the reaction-diffusion equation for glucose transport at each GC, based on the current tumor geometry and cell composition, and outputs the normalized glucose field ( $\frac{c}{c_n}$at each GC) for use by the metabolic model. The spatiotemporal-dependent inflow of glucose flux is constrained by the Michaelis-Menten kinetics model and an instantaneous optimization problem is solved for each cell position and time point by the metabolic model. During that time interval, the fluxes of the metabolic model are assumed constant. The metabolic model, given the available glucose concentration at every position in the computational grid, provides information regarding the uptake fluxes (e.g., glucose), intake fluxes (e.g., lactate) and the local proliferation rate of the tumor cells that reside within each GC. The latter is passed to the Oncosimulator to update its state.

Considering that Oncosimulator does not have as an intrinsic parameter the proliferation rate of the tumor cells as modulated by local metabolism (see model parameters in Table 1), a parameter transformation needs to take place. The local conditions of nutrient supply, such as glucose concentration, primarily regulate the withdrawal of tumor cells in a quiescent state, in an attempt by the tumor to sustain viability under conditions of reduced nutrient supply [31]. Hence, a reasonable first approximation, is to translate the proliferation rate, a, of the cell population to the fraction of newborn cells entering quiescent state, P_sleep, within the population. The following formula has been considered:

$$P_{sleep}=\frac{1-e^{a{T}_c}/2}{1-\left.{\left(P\right.}_{G0toG1}/T_{G0}\right)/\left(a+1/\left.T_{G0}\right)\right.},$$

(6)

where T_C the cell cycle duration, T_G0 the residence time of tumor cells in quiescent state and P_G0toG1 the fraction of quiescent cells re-entering cell cycle. Equation (6) has been derived from Eq. (7) in [32] that describes the proliferation rate of a tumor cell population with stem and progenitor cell hierarchy as regulated by the symmetric divisions of stem cells, the spontaneous apoptosis, the withdraw of cells in a quiescent state following mitosis and the cycle entry of quiescent cells:

$$e^{(\alpha +R_A)T_C}=\left(1+P_{sym}\right)\left(1-P_{sleep}+P_{sleep}\frac{P_{G0toG1}/T_{G0}}{R_A+1/T_{G0}+\alpha }\right),$$

(7)

where P_sym the fraction of stem cells that divide symmetrically and R_A the rate of spontaneous apoptosis. As a first approximation, we consider that metabolism has no effect on spontaneous apoptosis and symmetric divisions and, thereby, ignore stem hierarchy and apoptosis by setting P_sym=1 and R_A=0. The described parameter transformation has been implemented as intrinsic part of the Oncosimulator.

Furthermore, the proliferation rate, a, returned by the metabolic model consists an upper bound of the cell cycle (T_{C, max}=ln(2)/α), and corresponds to the case that no newborn cell withdraws to a quiescent state (P_sleep=0). For cell cycles longer than ln(2)/α, the P_sleep computed by (6) is negative and therefore biologically unrealistic.

The previously described cyclic coupling is implemented via MUSCLE platform [10]. The component models are called simultaneously and the data exchanges take place dynamically, i.e., during the runtime of the component models. Data transfer is triggered by the Oncosimulator at every predefined time interval. All the component models run until the Oncosimulation finishes executing. MUSCLE is triggered through TAVERNA.

2.2.3. Oncosimulator – Biomechanics Simulator

The hyper-modelling scenario dictates an iterative tightly-coupled communication scheme between the Oncosimulator (OS) and the Biomechanics simulator (BMS) (Figure 2 and Figure 3). The dynamic coupling is implemented via MUSCLE as previously described. The Oncosimulator computes cell proliferation in the case of free growth or cell loss in the case of treatment. Starting from an initial spatial map of cancer cell concentrations, BMS receives an (updated) tumor cell concentration map from OS in each time step. BMS and OS operate on distinct domains (OS: tumor only; BMS: tumor & healthy tissue) and discretization (OS: regular grid; BMS: unstructured mesh). Therefore, OS cell concentrations $c\left(\mathit{x}\right)$ are first mapped into the BMS simulation domain before spatial maps of local tumor-induced volumetric growth ${ϵ}_{growth}\left(\mathit{x}\right)$ can be computed. As the predicted mechanical stresses are intended to inform the spatial redistribution of existing tumor cells, we define the “growth” strain relative to a fixed maximum carrying capacity $c_0$, below (above) which additional tumor cells may still be accepted (will be rejected) in a region of interest:

$${\mathsf{ϵ}}_{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}(\mathbf{x})={\left(\frac{\mathrm{c}\left(\mathbf{x}\right)}{{\mathrm{c}}_0}\right)}^{\frac{1}{3}}-1,$$

(8)

Tumor-induced strain assumes positive (negative) values when the local tumor cell concentration exceeds (is less than) the maximum carrying capacity. From this strain map, tumor-induced mechanical stresses $\mathsf{\sigma }$(x) are computed using a linear-elastic material model with tissue-specific mechanical parameters, Young’s modulus $E$ and poisson ratio $\upsilon$. The resulting pressure field $\mathrm{p}\left(\mathbf{x}\right)=1/3\mathrm{t}\mathrm{r}\left(\mathsf{\sigma }\right(\mathbf{x}\left)\right)$ is then mapped back from the unstructured grid of the BMS simulation domain into the regular grid of the OS where the “direction of least-pressure” is computed as the normalized negative pressure gradient:

$$\mathrm{d}\left(\mathrm{x}\right)=-\frac{\nabla \mathrm{p}\left(\mathrm{x}\right)}{‖\nabla \mathrm{p}\left(\mathrm{x}\right)‖},$$

(8)

This information is returned to OS where it informs the movement and redistribution of tumor cells. Both models and their joint application (not as separate components) have been tested successfully for brain tumor simulation [26].

Figure 3. Data exchange and Computation in BMS-OS coupled execution.

Figure 3. Data exchange and Computation in BMS-OS coupled execution.

3. Results

4. Discussion

5. Conclusions

References

1. Goh WWB, Wong L. The Birth of Bio-data Science: Trends, Expectations, and Applications. Genomics, Proteomics & Bioinformatics. 2020;18(1):5-15. [CrossRef]
2. Anaya-Isaza A, Mera-Jiménez L, Zequera-Diaz M. An overview of deep learning in medical imaging. Informatics in Medicine Unlocked. 2021;26:100723. [CrossRef]
3. Stamatakos GS, Kolokotroni EA, Dionysiou DD, Georgiadi ECh, Desmedt C. An advanced discrete state–discrete event multiscale simulation model of the response of a solid tumor to chemotherapy: Mimicking a clinical study. Journal of Theoretical Biology. 2010;266(1):124-139. [CrossRef]
4. Shlomi T, Benyamini T, Gottlieb E, Sharan R, Ruppin E. Genome-scale metabolic modeling elucidates the role of proliferative adaptation in causing the Warburg effect. PLoS Comput Biol. 2011;7(3):e1002018. [CrossRef]
5. Borgdorff J, Bona-Casas C, Mamonski M, et al. A Distributed Multiscale Computation of a Tightly Coupled Model Using the Multiscale Modeling Language. Procedia Computer Science. 2012;9:596-605. [CrossRef]
6. Viceconti M, Bnà S, Tartarini D, et al. VPH-HF: A software framework for the execution of complex subject-specific physiology modelling workflows. Journal of Computational Science. 2018;25:101-114. [CrossRef]
7. http://www.chic-vph.eu.
8. Stamatakos GS, Georgiadi EC, Graf N, Kolokotroni EA, Dionysiou DD. Exploiting Clinical Trial Data Drastically Narrows the Window of Possible Solutions to the Problem of Clinical Adaptation of a Multiscale Cancer Model. Zhang L, ed. PLoS ONE. 2011;6(3):e17594. [CrossRef]
9. Dionysiou DD, Stamatakos GS, Uzunoglu NK, Nikita KS, Marioli A. A four-dimensional simulation model of tumour response to radiotherapy in vivo: parametric validation considering radiosensitivity, genetic profile and fractionation. J Theor Biol. 2004;230(1):1-20. [CrossRef]
10. Borgdorff J, Mamonski M, Bosak B, et al. Distributed multiscale computing with MUSCLE 2, the Multiscale Coupling Library and Environment. Journal of Computational Science. 2014;5(5):719-731. [CrossRef]
11. Steel GG, ed. Basic Clinical Radiobiology. 3. ed. Arnold; 2002.
12. Chen WW, Schoeberl B, Jasper PJ, et al. Input–output behavior of ErbB signaling pathways as revealed by a mass action model trained against dynamic data. Molecular Systems Biology. 2009;5(1):239. [CrossRef]
13. Choi M, Shi J, Jung SH, Chen X, Cho KH. Attractor landscape analysis reveals feedback loops in the p53 network that control the cellular response to DNA damage. Sci Signal. 2012;5(251):ra83. [CrossRef]
14. Hahnfeldt P, Panigrahy D, Folkman J, Hlatky L. Tumor development under angiogenic signaling: a dynamical theory of tumor growth, treatment response, and postvascular dormancy. Cancer Res. 1999;59(19):4770-4775.
15. Mirams GR, Arthurs CJ, Bernabeu MO, et al. Chaste: An Open Source C++ Library for Computational Physiology and Biology. Prlic A, ed. PLoS Comput Biol. 2013;9(3):e1002970. [CrossRef]
16. https://github.com/jmsgrogan/Chic.
17. Warburg O. On the origin of cancer cells. Science. 1956;123(3191):309-314. [CrossRef]
18. Cairns RA, Harris IS, Mak TW. Regulation of cancer cell metabolism. Nat Rev Cancer. 2011;11(2):85-95. [CrossRef]
19. Vander Heiden MG, Cantley LC, Thompson CB. Understanding the Warburg Effect: The Metabolic Requirements of Cell Proliferation. Science. 2009;324(5930):1029-1033. [CrossRef]
20. Duarte NC, Becker SA, Jamshidi N, et al. Global reconstruction of the human metabolic network based on genomic and bibliomic data. Proc Natl Acad Sci U S A. 2007;104(6):1777-1782. [CrossRef]
21. Ibarra RU, Edwards JS, Palsson BO. Escherichia coli K-12 undergoes adaptive evolution to achieve in silico predicted optimal growth. Nature. 2002;420(6912):186-189. [CrossRef]
22. Mahadevan R, Schilling CH. The effects of alternate optimal solutions in constraint-based genome-scale metabolic models. Metab Eng. 2003;5(4):264-276. [CrossRef]
23. E. Tzamali, V. Sakkalis, and K. Marias, “The effects of near optimal growth solutions in genome-scale humanME cancer metabolic model,” presented at the Bioinformatics & Bioengineering (BIBE), 2012 IEEE 12th International Conference, Larnaca, Cyprus, 2012.
24. Jain RK, Martin JD, Stylianopoulos T. The Role of Mechanical Forces in Tumor Growth and Therapy. Annu Rev Biomed Eng. 2014;16(1):321-346. [CrossRef]
25. Helmlinger G, Netti PA, Lichtenbeld HC, Melder RJ, Jain RK. Solid stress inhibits the growth of multicellular tumor spheroids. Nat Biotechnol. 1997;15(8):778-783. [CrossRef]
26. May CP, Kolokotroni E, Stamatakos GS, Büchler P. Coupling biomechanics to a cellular level model: an approach to patient-specific image driven multi-scale and multi-physics tumor simulation. Prog Biophys Mol Biol. 2011;107(1):193-199. [CrossRef]
27. Bauer S, May C, Dionysiou D, Stamatakos G, Buchler P, Reyes M. Multiscale Modeling for Image Analysis of Brain Tumor Studies. IEEE Trans Biomed Eng. 2012;59(1):25-29. [CrossRef]
28. Rikhtegar F, Kolokotroni E, Stamatakos G, Buchler P. A model of tumor growth coupling a cellular biomodel with biomechanical simulations. In: Proceedings of the 2014 6th International Advanced Research Workshop on In Silico Oncology and Cancer Investigation - The CHIC Project Workshop (IARWISOCI). IEEE; 2014:1-4. [CrossRef]
29. Abler D, Büchler P. Evaluation of a Mechanically Coupled Reaction–Diffusion Model for Macroscopic Brain Tumor Growth. In: Gefen A, Weihs D, eds. Computer Methods in Biomechanics and Biomedical Engineering. Lecture Notes in Bioengineering. Springer International Publishing; 2018:57-64. [CrossRef]
30. http://www.taverna.org.uk/.
31. Nik Nabil WN, Xi Z, Song Z, et al. Towards a Framework for Better Understanding of Quiescent Cancer Cells. Cells. 2021;10(3):562. [CrossRef]
32. Kolokotroni EA, Dionysiou DD, Uzunoglu NK, Stamatakos GS. Studying the growth kinetics of untreated clinical tumors by using an advanced discrete simulation model. Mathematical and Computer Modelling. 2011;54(9-10):1989-2006. [CrossRef]
33. Chou CH, Chang NW, Shrestha S, et al. miRTarBase 2016: updates to the experimentally validated miRNA-target interactions database. Nucleic Acids Res. 2016;44(D1):D239-247. [CrossRef]
34. Brenner DJ. The linear-quadratic model is an appropriate methodology for determining isoeffective doses at large doses per fraction. Semin Radiat Oncol. 2008;18(4):234-239. [CrossRef]
35. Kolokotroni E, Dionysiou D, Veith C, et al. In Silico Oncology: Quantification of the In Vivo Antitumor Efficacy of Cisplatin-Based Doublet Therapy in Non-Small Cell Lung Cancer (NSCLC) through a Multiscale Mechanistic Model. Komarova NL, ed. PLoS Comput Biol. 2016;12(9):e1005093. [CrossRef]
36. Kyroudis CA, Dionysiou DD, Kolokotroni EA, Stamatakos GS. Studying the regression profiles of cervical tumours during radiotherapy treatment using a patient-specific multiscale model. Sci Rep. 2019;9(1):1081. [CrossRef]
37. Li L, Wei Y, To C, et al. Integrated Omic analysis of lung cancer reveals metabolism proteome signatures with prognostic impact. Nat Commun. 2014;5(1):5469. [CrossRef]
38. Davidson SM, Papagiannakopoulos T, Olenchock BA, et al. Environment Impacts the Metabolic Dependencies of Ras-Driven Non-Small Cell Lung Cancer. Cell Metab. 2016;23(3):517-528. [CrossRef]
39. Greenspan HP. Models for the Growth of a Solid Tumor by Diffusion. Stud Appl Math. 1972;51(4):317-340. [CrossRef]
40. Faraldi M, Gomarasca M, Sansoni V, Perego S, Banfi G, Lombardi G. Normalization strategies differently affect circulating miRNA profile associated with the training status. Sci Rep. 2019;9(1):1584. [CrossRef]
41. Jeong J, Oh JH, Sonke JJ, et al. Modeling the Cellular Response of Lung Cancer to Radiation Therapy for a Broad Range of Fractionation Schedules. Clinical Cancer Research. 2017;23(18):5469-5479. [CrossRef]
42. Liu C, Tsao MS. Proto-oncogene and growth factor/receptor expression in the establishment of primary human non-small cell lung carcinoma cell lines. Am J Pathol. 1993;142(2):413-423.
43. Masuda N, Fukuoka M, Takada M, Kudoh S, Kusunoki Y. Establishment and characterization of 20 human non-small cell lung cancer cell lines in a serum-free defined medium (ACL-4). Chest. 1991;100(2):429-438. [CrossRef]
44. Durand RE, Sham E. The lifetime of hypoxic human tumor cells. Int J Radiat Oncol Biol Phys. 1998;42(4):711-715. [CrossRef]
45. Ginsberg T. Modellierung und Simulation der Poliferationsregulation und Strahlentherapie normaler und maligner Gewebe. Als Ms. gedr. VDI-Verl; 1996.
46. Chvetsov AV, Palta JJ, Nagata Y. Time-dependent cell disintegration kinetics in lung tumors after irradiation. Phys Med Biol. 2008;53(9):2413-2423. [CrossRef]
47. Kerr JF, Wyllie AH, Currie AR. Apoptosis: a basic biological phenomenon with wide-ranging implications in tissue kinetics. Br J Cancer. 1972;26(4):239-257. [CrossRef]
48. Gavrieli Y, Sherman Y, Ben-Sasson SA. Identification of programmed cell death in situ via specific labeling of nuclear DNA fragmentation. The Journal of cell biology. 1992;119(3):493-501. [CrossRef]
49. Bursch W, Paffe S, Putz B, Barthel G, Schulte-Hermann R. Determination of the length of the histological stages of apoptosis in normal liver and in altered hepatic foci of rats. Carcinogenesis. 1990;11(5):847-853. [CrossRef]
50. Rawlins EL, Hogan BLM. Epithelial stem cells of the lung: privileged few or opportunities for many? Development. 2006;133(13):2455-2465. [CrossRef]
51. Rawlins EL, Hogan BLM. Ciliated epithelial cell lifespan in the mouse trachea and lung. Am J Physiol Lung Cell Mol Physiol. 2008;295(1):L231-234. [CrossRef]
52. Lippmann M. Environmental Toxicants: Human Exposures and Their Health Effects. 3rd ed. J. Wiley; 2009.
53. Flindt R. Amazing Numbers in Biology. Springer-Verlag; 2006.
54. Pine SR, Ryan BM, Varticovski L, Robles AI, Harris CC. Microenvironmental modulation of asymmetric cell division in human lung cancer cells. Proc Natl Acad Sci U S A. 2010;107(5):2195-2200. [CrossRef]
55. Morrison BJ, Steel JC, Morris JC. Sphere culture of murine lung cancer cell lines are enriched with cancer initiating cells. PLoS One. 2012;7(11):e49752. [CrossRef]
56. Van Leeuwen CM, Oei AL, Crezee J, et al. The alfa and beta of tumours: a review of parameters of the linear-quadratic model, derived from clinical radiotherapy studies. Radiat Oncol. 2018;13(1):96. [CrossRef]
57. Ishizawa K, Rasheed ZA, Karisch R, et al. Tumor-Initiating Cells Are Rare in Many Human Tumors. Cell Stem Cell. 2010;7(3):279-282. [CrossRef]
58. Törmänen U, Eerola AK, Rainio P, et al. Enhanced apoptosis predicts shortened survival in non-small cell lung carcinoma. Cancer Res. 1995;55(23):5595-5602.
59. Moon SW, Kim JJ, Jeong SC, Kim YH, Han JW. Clinical significance of tumor necrosis and viability in non-small cell lung cancer. J Thorac Dis. 2022;14(4):892-904. [CrossRef]
60. O’Neill AJ, Staunton MJ, Gaffney EF. Apoptosis occurs independently of bcl-2 and p53 over-expression in non-small cell lung carcinoma. Histopathology. 1996;29(1):45-50. [CrossRef]
61. Zaider M, Minerbo GN. Tumour control probability: a formulation applicable to any temporal protocol of dose delivery. Phys Med Biol. 2000;45(2):279-293. [CrossRef]
62. Grannis FW. Limitations of molecular testing in combination with computerized tomographic for lung cancer screening. Nat Commun. 2022;13(1):3890. [CrossRef]
63. Groninger E, Meeuwsen-de Boar T, Koopmans P, et al. Pharmacokinetics of vincristine monotherapy in childhood acute lymphoblastic leukemia. Pediatr Res. 2002;52(1):113-118. [CrossRef]
64. Hill CR, Cole M, Errington J, Malik G, Boddy AV, Veal GJ. Characterisation of the Clinical Pharmacokinetics of Actinomycin D and the Influence of ABCB1 Pharmacogenetic Variation on Actinomycin D Disposition in Children with Cancer. Clin Pharmacokinet. 2014;53(8):741-751. [CrossRef]
65. Theerakitthanakul K, Saetang J, Kruatong J, et al. Senescence Process in Primary Wilms’ Tumor Cell Culture Induced by p53 Independent p21 Expression. J Cancer. 2016;7(13):1867-1876. [CrossRef]
66. Royer-Pokora B, Busch MA, Tenbusch S, et al. Comprehensive Biology and Genetics Compendium of Wilms Tumor Cell Lines with Different WT1 Mutations. Cancers. 2020;13(1):60. [CrossRef]
67. Berrebi D, Leclerc J, Schleiermacher G, et al. High Cyclin E Staining Index in Blastemal, Stromal or Epithelial Cells Is Correlated with Tumor Aggressiveness in Patients with Nephroblastoma. Westermark P, ed. PLoS ONE. 2008;3(5):e2216. [CrossRef]
68. Jurić I, Pogorelić Z, Kuzmić-Prusac I, et al. Expression and prognostic value of the Ki-67 in Wilms’ tumor: experience with 48 cases. Pediatr Surg Int. 2010;26(5):487-493. [CrossRef]
69. Craft A. Growth rate of Wilms’ tumour. The Lancet. 1999;354(9184):1127. [CrossRef]
70. Middleton P, Banieghbal B, Pitcher R, Schubert P. Radiological response and histological findings in nephroblastoma: Is the any correlation? Afr J Paediatr Surg. 2020;17(3):39. [CrossRef]
71. Lee CL, Blum JM, Kirsch DG. Role of p53 in regulating tissue response to radiation by mechanisms independent of apoptosis. Transl Cancer Res. 2013;2(5):412-421.
72. Roose T, Chapman SJ, Maini PK. Mathematical Models of Avascular Tumor Growth. SIAM Rev. 2007;49(2):179-208. [CrossRef]
73. https://ec.europa.eu/digital-single-market/en/news/great-achievements-eu-funded-chic-project.