A multiscale cell‐based model of tumor growth for chemotherapy assessment and tumor‐targeted therapy through a 3D computational approach

Abstract Objectives Computational modeling of biological systems is a powerful tool to clarify diverse processes contributing to cancer. The aim is to clarify the complex biochemical and mechanical interactions between cells, the relevance of intracellular signaling pathways in tumor progression and related events to the cancer treatments, which are largely ignored in previous studies. Materials and Methods A three‐dimensional multiscale cell‐based model is developed, covering multiple time and spatial scales, including intracellular, cellular, and extracellular processes. The model generates a realistic representation of the processes involved from an implementation of the signaling transduction network. Results Considering a benign tumor development, results are in good agreement with the experimental ones, which identify three different phases in tumor growth. Simulating tumor vascular growth, results predict a highly vascularized tumor morphology in a lobulated form, a consequence of cells' motile behavior. A novel systematic study of chemotherapy intervention, in combination with targeted therapy, is presented to address the capability of the model to evaluate typical clinical protocols. The model also performs a dose comparison study in order to optimize treatment efficacy and surveys the effect of chemotherapy initiation delays and different regimens. Conclusions Results not only provide detailed insights into tumor progression, but also support suggestions for clinical implementation. This is a major step toward the goal of predicting the effects of not only traditional chemotherapy but also tumor‐targeted therapies.

simulations of the main physical mechanisms underlining tumor development leads to a real possibility for evaluating pre-clinical drug design opportunities and helping the optimization of drug delivery. 1,2 The development of benign tumors is caused by excessive cell proliferation, which is commonly limited by space, or more usually, by the nutrient availability in the tissue. This early phase of tumor development, labeled avascular growth, has tumor cells in the hypoxic cell state, where they are able to survive with a lower nutrient concentration. Hypoxic cells within a tumor express the hypoxiainducible factor-1 (HIF-1), upregulating pro-angiogenic factors and triggering tumor vascularization. A denser vasculature gives access to an additional supply of nutrients, driving the tumor growth to the vascular phase. Tumor vascular growth is a feature of malignancy, enabling tumor cells to invade other tissues by entering the circulation via nearby blood vessels (metastasis). 3,4 Various mathematical techniques have been used to simulate tumor growth and associated processes, being applied to tumors both in the avascular [5][6][7] and at the vascular stages. [8][9][10][11] Indeed, most of the early mathematical models of tumor growth address avascular tumor morphology. Developing mathematical models of tumorinduced angiogenesis permits a more realistic description of nutrient availability in tumors. [12][13][14][15] See [16][17][18][19] for complete reviews of angiogenesis models. Recently, multiscale approaches have been introduced to reproduce the biological and physical mechanisms in tumor growth and angiogenesis. These models consider both subcellular and tissue scales. [19][20][21][22][23] Alarcón et al. [24][25][26] introduced a hybrid structured lattice-based model to simulate vascular tumor growth, which considers blood flow and oxygen transport in a tissue scale and accounts for cellular interactions and progress in cellular and intracellular scales. The model investigates the effects of nutrient spatial heterogeneity on the evolution and invasion of cancerous tissue, and the emergent growth laws. The diffusive transport of oxygen and VEGF within the tissue is described through reaction-diffusion equations. In a subsequent study, the authors considered the diffusion of standard cytotoxic drugs as a treatment and investigated the effects of vessel normalization on chemotherapy. 26 The authors concluded that vessel normalization improves the efficiency of chemotherapeutic drugs, observing also the decrease in the prevalence of hypoxia.
Subsequently, Owen et al. 27  The authors reported that intraperitoneal infusion is the preferable route in the initial growth phase, when the tumor is still small and avascular. Jafari Nivlouei and co-workers 36 proposed a 2D multiscale agent-based model, addressing two distinct phases in tumor growth.
In each stage, tumor progression is considered with and without normal healthy cells. The authors reported the formation of a dense intra-tumoral vascular network throughout the entire tumor mass as a sign of a high malignancy grade.
In what concerns cancer therapy, computational studies focusing on tumor response to therapy are fundamental tools that facilitate the understanding of drug's mechanism of action, helping to determine the most effective treatment protocols. Different approaches to chemotherapy modeling have been proposed, including continuous 37 and hybrid discrete-continuum models, in which the model describes the effect of interstitial fluid pressure and lymphatic drainage on drug delivery to tumors. [38][39][40] This class of studies permits the evaluation of the parameters that limit the delivery of nutrients and therapy. More recently, multicellular and multiscale techniques, which incorporate drug therapy at the extracellular level, are becoming increasingly important in tumor treatment simulation. Wang et al. 41 developed a multiscale agent-based model to simulate the melanoma tumor vascular growth and to study the response of tumor to combined drugs treatments. The authors reported that the interruption in the communications between melanoma cells and the vasculature might increase the drugs' effectiveness.
Targeted therapy is a novel type of treatment which reduces systemic drug toxicity by inducing modifications in the tumor microenvironment and not in normal cells. 42,43 Targeted drugs are characterized by the binding of their therapeutic molecules to specifically expressed receptors on the tumor cells' membranes. Kim and co-workers 44 modeled targeted therapy by focusing on specific intracellular signaling pathways that prevent cancer cells' abnormal behavior and finally induce cell apoptosis or suppress cell growth.
They developed a hybrid model where they targeted the MAPK and PI3K-AKT signaling pathways, which are activated in lung cancer, and used it to investigate the effects of this pathway inhibition under different microenvironmental conditions. The model uses a CA approach to describe the cellular process and a set of ODEs to address tumor response to the targeted therapy. The authors suggested a new treatment combination strategy based on the predicted cell signaling responses. More recently, a new class of targeted therapies has been developed, targeting cells in the hypoxic regions. In hypoxia-activated pro-drugs (HAPs), the cytotoxic agents are released under low oxygen pressure. 45 Hong et al. 46  This provides a significant and novel contribution to the field of simulating tumor growth and different methods of cancer treatment in a simplified way. The model integrates all the information received from each spatial and temporal scale to predict the system response.
It enables us to survey cell phenotypic alterations by considering the interaction of signaling molecules and the signaling pathways. This helps to explore the mechanism of anti-tumor and ECM-targeted strategies by inhibiting the activity of specific receptors. Results will not only provide detailed insights into tumor progression, but the model is also a step toward clinical implementation. This represents an opportunity to analyze tumor response to both treatment strategies (i.e., chemotherapy and combination therapy) and to evaluate typical clinical protocols.

| MATERIAL AND ME THODS
The model simulates the tumor development processes at intracellular, cellular, and extracellular scales. Each scale of the model is presented in the following sections to detail the implemented mechanisms of tumor progression. Oncogenic mutations not only cause the overexpression of genes, but also can produce mutated proteins. Growth factor receptor tyrosine kinases (RTKs) are often de-regulated in neoplasms.

| Intracellular scale
These are proteins involved in commonly activated survival signaling pathways, whose activity leads to stimulation of serine/threonine kinases (e.g., Raf and Akt) and lipid kinases (e.g., PI3Ks) through the activation of small GTPases (e.g., Ras). Moreover, the activation of Ras can originate from the loss of neurofibromin (NF1) protein, encoded by the NF1 gene. NF1 functions as a tumor suppressor that negatively regulates the activity of Ras. [48][49][50][51] Loss of tumor suppressors' function results in cancer initiation and progression because of their role in cell division inhibition, induction of apoptosis, and metastasis suppression. While the hyperactivation of Ras-ERK and PI3K-Akt signaling pathways can lead to excessive proliferation in tumor cells, mutations can promote the cancer phenotype by disabling cell death signaling. 52,53 For instance, p53 is known as a tumor suppressor protein whose loss through mutation can contribute to tumor development by the interruption of cell death signaling, as it regulates cell apoptosis by binding directly to Bax, a pro-apoptotic protein. 54 65,66 Furthermore, the endothelial cell (EC)specific cadherin transmembrane receptor, VE-cadherin, whose association with the protein ß-catenin facilitates its binding to the actin cytoskeleton, is responsible for the tight but dynamic connection between neighboring cells. 67 This protein plays an important role in providing a cohesive structure for the new blood vessel.
Strikingly, the cadherin-catenin adhesion system regulates cell proliferation and migration through downstream signaling effects during cancer development. E-cadherin loss of expression promotes the release of ß-catenin into the cytosol, which results in the activation of Wnt signaling. ß-catenin is the main effector of the Wnt signaling pathway, 68,69 and E-cadherin negatively regulates the Wnt/ ß-catenin signaling. Nevertheless, the loss of expression of cadherin by itself is not sufficient for the activation of ß-catenin signaling. 70 Recent research reveals that the loss of APC function is associated with increased levels of ß-catenin. 71,72 APC is a tumor suppressor localized inside the cells' nucleus, which regulates cell proliferation by inhibiting Wnt/ß-catenin signaling, and facilitates cell apoptosis to suppress tumor progression and metastatic cell spread. 73 The reintroduction of APC into mutant cells, in order to restore its function in Wnt/ß-catenin signaling, has been investigated in several therapeutic treatments. [74][75][76] Experimental observations demonstrate that the early loss of APC function and the activation of ß-catenin can be followed by the later loss of E-cadherin, leading to the cell's invasive behavior. 70 Considering the described key events in tumor growth and angiogenesis, a signaling cascade is modeled based on the cross talk between the main regulators of growth factors (RTKs), integrin, cadherin, and Wnt ( Figure 1A). Moreover, different experimental studies are used to integrate the information of the most important effectors that play a key role in cell cycle regulation, as presented in Table 1. The dependences between the network nodes are specified by the arrows, which indicate the activation of the corresponding effector. On the contrary, an inhibitory effect is pictured as barheaded lines. The aim of modeling the intracellular scale is to determine the cell phenotype in response to the active signals. In the current study, the Hamiltonian is considered as a sum of four effective energy terms describing cellular adhesion, cell growth, chemotaxis, and guaranteeing cell continuity. The sum is run over all the neighboring pixels, and ′ are the cells' ID, and δ is the Kronecker symbol.
• Cell growth: describes the energy involved in cell growth and proliferation through mitosis. During the cell cycle progression, cells grow until they reach twice their initial volume, and then they divide. After mitosis, the parent cell target volume is restored to its initial volume, and the offspring cell will inherit the type and target volume of the parent cell, being assigned a new unique ID. The energy responsible for setting the target cell size in the Hamiltonian is: where v denotes the cell volume, while V T is the target volume and e is the cell elasticity. σ is a negative parameter that determines the chemotaxis intensity and is denoted as chemotactic potential. • Cell continuity: Cells are a continuous medium. To keep the continuity of lattice sites that are occupied by a single cell, a constraint term is added to the Hamiltonian. This term introduces a severe increase in the system total energy when a cell is about to rupture: where is a large penalty factor that increases the system energy when there is a difference between the current contiguous cell size (v ) and the number of lattice sites occupied by the cell with unique Therefore, the contribution of the energy terms related to cell adhesion, volume, continuity, and chemotaxis, which is referred collectively as Hamiltonian, will be as follows: Boolean dependence relations of the nodes presented in the signaling network of Figure 1A, based on the experimental data of given references. The colors correspond to the color coding of the nodes in Figure 1A Node Dependence relation Reference where C is the chemoattractant agent, replaced by the nutrient concentration n in what concerns tumor cells, or the VEGF concentration (V) for the activated endothelial cells.

| Extracellular scale
At this level, concentrations of nutrients and growth factor, and the drug distribution are described by a set of reaction-diffusion PDEs.
• where n denotes nutrient concentration, D n is its diffusion coefficient, S n refers to the process of nutrients release from vessels, and B is a function expressing the uptake rate of nutrients by tumor cells, as described below: where β is the maximum consumption rate of nutrients per cell voxel for a tumor cell. The release rate of nutrients from the endothelial cells is given by: •  where c denotes the chemotherapy drug concentration, D c its diffusion coefficient, and k c represents the drug decay rate. R is the function of drugs' uptake rate by tumor cells, as described below: ρ is the maximum consumption rate of drug per proliferative cell voxel. The drug release rate from the ECs will be: • Initial and boundary conditions: The simulation starts with a tumor, with a diameter of approximately 65 µm, at the center of the computational domain. Initially, it is assumed that the signaling from all receptors, including the RTK receptors, is active, and the concentration of nutrients is sufficient to irrigate the cells. So, the initial and the boundary conditions are: n(x, y, z, t) (x,y,z)⊂ECs = s n , n (x, y, z, 0) = S 0 = 4.6pg∕voxel.
Since the secretion of VEGF is induced in the hypoxic area of the tumor, there is no VEGF concentration in the domain until a hypoxic core is formed, which means V (x, y, z, 0) = 0. In response to hypoxia, VEGF is secreted at a rate s V inside the tumor core.

| Simulation algorithm
To All parameter values used in the model are listed in Table 2.

| Computational setup
Using the open-source CompuCell3D simulation environment Initially, the simulation starts with a tumor size of ∼ 65 μm , containing proliferating cells, surrounded by pre-existing vessels.
Nutrients are constantly diffused from the vascular network, while the diffusion of the chemotherapeutic drug is carried out according to a specific treatment protocol. Moreover, the concentration of VEGF secreted from hypoxic cells is calculated, being a driver for ECs' activation, which leads to neo-vessel growth.
A sensitivity analysis is performed to tune the model and to identify and adjust key parameters. Hence, by varying a particular parameter at a time (and keeping fixed all the other Table 2 Table 2 show a balance between the contact guidance and the cell-matrix adhesion energy.
Assessing the compressibility properties, cell size is sensitive to γ e , which when large made the cells small, resistant to deformation and requiring more energy to grow. The results are also insensitive to the value of T m until it is increased by more than two orders of magnitude. The larger values of T m , the larger the cell membrane fluctuations.
To investigate the sensitivity of the obtained results to changes on the signaling thresholds, comparisons between numerical simulations and experimental data were performed, and the main results are the following: • A higher activation threshold of a receptor means that the corresponding receptor is unlikely to be activated.
• Increasing the threshold for RTK receptor activation, a regulator of cell survival suppresses tumor progression.
• The E-cadherin threshold controls contact inhibition of growth.

| Cell phenotype assessment
The cells' dynamics are controlled by the signals they constantly receive from their microenvironment. Hence, a Boolean network model is employed to predict the cell phenotype from the various receptor activation states, dependent on the implemented signaling cascades ( Figure 1A). The input-output map extracted from this Boolean network is presented in Figure 1B. The states' activation is represented • Consistent with experimental observations reported in, 137-140 cell apoptosis is the dominant phenotype when a disruption in the activity of either RTK or integrin receptors occurs. This is independent from E-cadherin activity.
• In the presence of RTK and integrin, signaling from the cadherin regulates cell motility, confirming its role in cell-cell contacts. Cell-cell adhesion matrix • In the absence of Wnt signaling, cell migration can be blocked by cadherin (case 110).

TA B L E 2 Parameters used in the model and corresponding references
• Although the reintroduction of APC into mutant cells is explored as a therapeutic strategy to drive cell apoptosis, by inhibiting pathways activated by the loss of APC, including Wnt/β-catenin, 74,75 according to the Boolean network implemented there is no clear evidence for an APC role in apoptosis and in control of Wnt signaling. [141][142][143] • In the case of main receptors' activity, the presence of NF1 impairs further progression to malignancy by inducing apoptosis.
This shows a potential treatment strategy by restoring NF1, although there is no systemic therapy until now.
To incorporate the resulting mapping, it should be noted that the activation of signals from integrin, RTK, and Wnt is related to E-cadherin activation. [144][145][146][147] The implementation of Wnt activity depends on the E-cadherin loss, which stimulates canonical Wnt signal- ing. E-cadherin is related to cell-cell contacts; the connection of each cell with its adjacent cells is a criterion that determines the activity of E-cadherin receptor signaling. Integrin activity is associated with its role in cells binding to the ECM and, therefore, the connection between cells and ECM defines the activation of the integrin receptor.
The activation of the RTK receptor signaling is controlled by nutrients availability, due to the role of the PI3K-Akt pathway in the promotion of glycolysis, necessary for cell growth. 110 Figure 1A, since it is not detectable through the signaling network. Tumor cells are able to be in a quiescent slow-growing state in regions of hypoxia and nutrient deprivation in areas far from vessels.

| Model verification
To measure the network robustness against the fluctuations and, at the same time, investigate whether the signal transduction is independent of the initial internal states, a systematic simulation is performed for all possible 2 29 initial combinations of states of all 29 internal components of the network described in Figure 1A. Results So, there is no significant increase in tumor size in the third phase of growth. See Video S1 and Video S2 for tumor avascular growth.  Figure 4A shows the first activated EC, and quiescent and necrotic cells.  Figure 5B, the tumor center is a nutrient-depleted area with a large quiescent and necrotic region. Hence, in order to access the nutrients, the tumor core releases angiogenic factors,

| 3D tumor growth in a vascular network
such as VEGF at a constant rate to stimulate the growth of new capillaries. It leads to the accumulation of VEGF throughout the tumor, being the highest density inside it ( Figure 5B). Figure 5C shows the highly vascularized tumor with a lobulated form on day 110 from different perspectives (See also Video S4).

| Chemotherapy
Chemotherapy is the application of drugs that target, in general, rap- Referring to Equation 7, drug concentration in the tumor microenvironment is calculated to assess the response to chemotherapy.

| Chemotherapy Initiation
Several investigations reported that the delay in treatment initiation may cause recurrence of cancer and have a negative impact on overall survival. [155][156][157][158] By contrast, several other studies proved that there is not a clear relation between treatment initiation delays and outcomes. [159][160][161][162] Hence, the model analyzes the effect of delay in the initiation of treatment on outcomes. Figure [163][164][165][166] However, more investigations are needed to optimize metronomic chemotherapy for each tumor type. 167 Here, drug delivery is modeled on a 10-day periodic MTD schedule that yields a well-responding therapy and is compared to when the drug is applied on a 6-day periodic metronomic schedule ( Figure 7C).   (Table 2). In what concerns the combination therapy, the cytotoxic chemotherapy is applied to the tumor at the same time, through Equation 9.
The results on tumor response to the different therapies are presented in Figure 7D, in which chemotherapy is applied accord- In order to investigate the new vasculature behavior, a comparison of tumor response to two different chemotherapy protocols is presented ( Figure 7E-F). The first case considers applying the low dose of 5 µg/m 2 at an early stage of tumor growth but it leads to treatment failure. Under these conditions, endothelial cells are rapidly expanding. In the second case, a dose of 7.5 µg/m 2 is applied on the 90th day of tumor growth, yet it was unable to avoid tumor reoccurrence after 8 cycles. In a similar way to a previous case, the successful treatment, the growth of new arteries also stopped, to some extent, and no significant change in the number of endothelial cells was reported; however, with tumor recurrence, the angiogenic rate increased ( Figure 7F).

| DISCUSS ION
In this study, mathematical and computational modeling methods for simulation of tumor growth and angiogenesis are used to explore opportunities in the development and testing of novel treatment strategies, including targeted therapies. The control of the signals involved in cell proliferation and, even more importantly, in cell apoptosis, is still a challenge in tumor reduction and/or elimination.
Here, a multiscale model is presented to test this explicit aim in order to clarify how the signaling transduction operates and affects im- and it depends on the adhesion energy, including the E-cadherin and integrin-mediated cellular adhesion and signaling from growth factors that control the growth rate. As the tumor grows and the vascular network expands, the initial spherical tumor shape changes into a lobulated form. Limitations on resources and space cause cell migration toward less dense locales, avoiding over-crowd regions. As a result, the geometry of the tumor is more irregular and disordered than the avascular tumor, and it reflects a more complex phenomenon that cannot be reproduced in 2D modeling.
Our research reveals the capability of a multiscale and threedimensional numerical simulation of tumor progression to explore the outcomes of drug treatment. 2D models yield valuable insights into the growth and dissemination of tumors, tumor-induced angiogenesis, and vascular remodeling. As a matter of fact, the 2D assumptions are indeed acceptable when tumors are either approximately flat or when they have important symmetries (e.g., when they are spherically symmetric). Although modeling in 2D is an attractive alternative to 3D calculations, as it requires notably less computational resources, in the context of solving mechanical forces that mediate cell shape and orientation, the 2D hypothesis might not be sufficiently accurate in predicting tumor behavior as a growing mass. In what concerns the modeling of cytotoxic chemotherapy, tumor responses to multiple cycles of chemotherapy are simulated, including treatment failure, relapse at a distance from the primary tumor, and effective therapy. Different protocols aiming at treatment efficacy optimization have been investigated. The model performs a dose comparison that helps to evaluate the dosing efficacy.
Moreover, analyzing the effect on outcome of a delay in chemotherapy initiation indicates that applying therapy at the earliest stage of tumor development leads to a remarkable reduction in its size over time. In highly vascularized tumors, representing high-grade cancers, the initiation delay does not guarantee the treatment success since tumor eradication has not been observed. However, a significant effect on the long-term prospects is obtained, which reflects the longer period of tumor dormancy (76 days), and postpones recurrence of tumor for more than a month (40 days).
Since drug administration frequency is important to attain the desired pharmacologic effects, and to reduce adverse reactions, the

CO N FLI C T O F I NTE R E S T
The authors have declared that no competing interests exist.

ACK N OWLED G M ENTS
The authors gratefully acknowledge the support provided by Isfahan

DATA AVA I L A B I L I T Y S TAT E M E N T
Any data generated in this study are available from the corresponding author upon reasonable request.