Agent-Based Models of the Tumor Microenvironment: Predicting Flow-Mediated Invasion and Cancer Viability in Response to Interstitial Flow, Chemotherapy, and Stromal Cell Density

Therapeutic delivery within the tumor is attenuated by
static pressure and the permeability of the tissue, which is
mediated by the extent of ECM remodeling. Additionally, cellular
signaling from stromal fibroblasts limits the efficacy of
treatment. This phenomenon is of particular relevance at the tumor
border, where there is a transitional region from predominantly
cancer cells to predominantly stromal cells. These factors promote
cancer growth and select for a potentially metastatic
subpopulation. Furthermore, interstitial fluid flow at the tumor
border affects the migration characteristics of cancer cells.
Interstitial fluid flow has been implicated in an increase in
cancer cell invasion, which is the major effector of cancer spread
and decreased patient viability. In brain cancer, two different
mechanisms have been implicated in this increased invasion, and
this invasion is mediated by both CXCR4 and CD44. Therefore, TME
complexity necessitates utilization of robust models to elucidate
the effects of the TME on cancer development and progression. The
current work defines two novel agent-based models that describe and
predict cancer specific outcomes within in vitro TME mimetic
systems. Our models indicate that brain cancer cell invasion is
increased in the presence of interstitial flow, and this increased
invasion is driven by two specific mechanisms, CXCR4-CXCL12
autologous chemotaxis and CD44 mechanotransduction. Additionally,
in vitro and in silico models of the tumor border transition zone
predict that regional viability within the breast tumor border is
affected by both fibroblast signaling and transport mechanisms, and
this affect could be selecting for important resistant
subpopulations.


CHAPTER 1: Introduction
The tumor microenvironment (TME) is composed of varied factors. Cancer cells display unique characteristics that lead to the development of disordered vasculature, hypoxia, elevated pH, and necrosis.
Additionally, the developing tumor is affected by several local environment factors including cancer associated cell signaling, extracellular matrix remodeling, and interstitial fluid flow. Together, these TME factors lead to pro-metastatic phenomena.
Therapeutic delivery within the tumor is attenuated by static pressure and the permeability of the tissue, which is mediated by the extent of ECM remodeling. Additionally, cellular signaling from stromal fibroblasts limits the efficacy of treatment. This phenomenon is of particular relevance at the tumor border, where there is a transitional region from predominantly cancer cells to predominantly stromal cells. These factors promote cancer growth and select for a potentially metastatic subpopulation.
Furthermore, interstitial fluid flow at the tumor border affects the migration characteristics of cancer cells.
Interstitial fluid flow has been implicated in an increase in cancer cell invasion, which is the major effector of cancer spread and decreased patient viability. In brain cancer, two different mechanisms have been implicated in this increased invasion, and this invasion is mediated by both CXCR4 and CD44.
Therefore, TME complexity necessitates utilization of robust models to elucidate the effects of the TME on cancer development and progression. The current work defines two novel agent-based models that describe and predict cancer specific outcomes within in vitro TME mimetic systems. Our models indicate that brain cancer cell invasion is increased in the presence of interstitial flow, and this increased invasion is driven by two specific mechanisms, CXCR4-CXCL12 autologous chemotaxis and CD44 mechanotransduction. Additionally, in vitro and in silico models of the tumor border transition zone predict that regional viability within the breast tumor border is affected by both fibroblast signaling and transport mechanisms, and this affect could be selecting for important resistant subpopulations.
Together, our data describe varied TME factors that mediate cancer cell response, and the utilization of agent-based models has facilitated the prediction, characterization, and explanation of in vitro findings.

The Tumor Microenvironment
The tumor microenvironment is generally defined as the local factors that influence the growth and progression of cancer. While the components of the tumor microenvironment are cancer type and context specific, there are several common constituents inherent to tumors. Cancer associated cells ( Figure 1F/E), the physical adherent matrix or extracellular matrix ( Figure 1D), and biophysical forces ( Figure 1H) are several key components of the tumor microenvironment 1,2 , while vasculature is additionally important for drug delivery considerations ( Figure 1A). Characterization of the tumor microenvironment provides insight into the mechanisms of cancer migration, proliferation, metastasis, and chemotherapy resistance, thus facilitating better standard of care and increased treatment efficacy.

Tumor Characteristics
The most abundant component of the tumor microenvironment is cancer cells ( Figure 1B). These cells are defined by several universal hallmarks, classically defined by Hanahan and Weinberg, including growth suppression evasion, mutation, replicative immortality, inflammation, invasion mechanisms, and stimulating angiogenesis 3,4 . These ubiquitous cancer hallmarks facilitate the development of unique tumor characteristics. Unchecked growth coupled with replicative immortality yields increased metabolic activity and associated decrease in pH 5 within the tumor microenvironment. Increased metabolic activity additionally causes oxygen depletion 6 and regions of hypoxia. Therefore, unchecked growth creates a detrimental environment for the cancer cells and subsequent regions of necrosis ( Figure 1C).
Additionally, this environment creates a selective pressure that yields regional subpopulations that are more likely to proliferate and invade 7 .
The lack of oxygen within tumors causes cancer cells to secrete pro-angiogenic factors to stimulate the growth of new blood vessels 8 ( Figure 1A), however, this tumor vasculature, or neovasculature, is significantly altered from normal tissue blood vessels. Neovasculature is defined by abnormal blood vessel shapes, and large fenestrations caused by misalignment of endothelial cells 9 . These gaps allow cells to intravasate into the vasculature, which has been proposed as a mechanism of metastasis 10 . Interestingly, these neovasculature abnormalities also lead to a distinct phenomenon known as the enhanced permeability retention (EPR) effect. Treatment of tumors with large molecule drugs is enhanced over small molecule treatment due to accumulation into the tumor by extravasation from the leaky vasculature into the tumor bulk, while simultaneously being blocked from extravasation from the vasculature into other organs due to normal epithelial junction size exclusion 10 .

Tumor Associated Cells
While cancer cells are the most abundant cell type in the TME, they do not act alone. Propinquity to surrounding stroma yields infiltration of several cell types into the tumor bulk, which then interact with the cancer cells. These cell types include immune cells, such as T lymphocytes and B lymphocytes, myeloid cells, such as neutrophils and macrophages, adipocytes, and fibroblasts 11 . T lymphocytes and B lymphocytes either help promote or hinder the growth of cancer depending on the particular secreted factors and the cancer context 12,13 . Tumor associated macrophages are generally considered to be protumorgenic 13 , and these macrophages have been shown to enhance cancer cell migration 14 . Adipocytes recruit cancer cells in specific via adipocyte specific secreted factors 15 .
However, perhaps the most important cell type in the development and progression of cancer is the fibroblast ( Figure 1G). Normally, fibroblasts are active in the native wound healing response. After an injury, platelets in the blood activate fibroblasts by releasing transforming growth factor beta (TGF-β); the fibroblasts then facilitate the remodeling and repair of local tissue. Activated fibroblasts secrete proinflammatory and pro-angiogenic cytokines that further stimulate remodeling and repair, and these fibroblasts deposit the structural components of the tissue, or the extracellular matrix (ECM) 16 . After the wound has been healed, the fibroblasts are inactivated or removed from the injury site.
Contrastingly, in chronic wounds such as ulcerative colitis, where fibroblasts remain activated for extended periods of time, the risk of cancer development is much greater than in normal tissue due to the pro-tumorgenic properties of these cells 17 . Interestingly, developing tumors displays many of the characteristics of a chronic wound as inflammation, increased angiogenesis, and matrix remodeling by fibroblasts ( Figure 1D) are all present in tumor growth 18 . Cancer associated fibroblasts (CAFs) are the most abundant cell type in the tumor stroma, and cancer cells have been shown to secrete fibroblast activating factors including TGF-β. Activated fibroblasts in the TME ( Figure 1F), in turn, secrete mitogenic growth factors ( Figure E), and secrete their own TGF-β, which suppresses cells in the immune system 19 . Finally, fibroblast-cancer interactions have been directly implicated in decreased chemotherapy efficacy 20 , leading to sustained tumor growth and the potential for subsequent metastasis.

Extracellular Matrix
The extracellular matrix (ECM) is defined by the various macromolecules that the give tissue structural integrity and rigidity. These molecules include proteins, glycoproteins, proteoglycans, and polysaccharides, with unique physical properties that allow for chemical crosslinking and tissue cellular development 21 . The ECM provides a unique biological scaffold for cell adherence, growth, and development. Cells utilize several signaling mechanisms to respond to biomechanical forces and sense the rigidity of the matrix. The extent of ECM elasticity allows cells to sense external stresses and modulate their behavior, and the elasticity is also driver for cell differentiation and stem cell maintenance 22 25 .
In tumors, the ECM is significantly altered from normal tissue ECM. Normal regulatory function is disrupted in CAFs, and these fibroblasts actively deposit collagen in the TME, leading to fibrosis and increased stromal density 26 . In addition to upregulation of collagen deposition, matrix degradation enzymes are concurrently downregulated, which leads to further accumulation of ECM in the TME. For example, breast cancer tissue is 10 times stiffer than normal tissue 27 . The orientation of fibers is also altered in tumors. While normal tissue is characterized by random non-oriented fibers, the ECM in tumors is linearized parallel to epithelium, which has been hypothesized to support cellular migration and invasion. In fact, metastatic cancer progression is defined by upregulation of ECM sensing mechanisms such as CD44, and cancer cell populations with these receptors are more likely to survive and spread 28 .

Biophysical Forces
Increased ECM deposition, unchecked cellular growth, and compromised vasculature all lead to elevated pressure in the tumor bulk. Interestingly, as classically described by Jain and colleagues 29,30 , this pressure is consistently high throughout the tumor bulk, as compared to the tumor border. At the tumor border, the pressure is normal compared to the elevated pressure in the tumor bulk, which drives interstitial fluid from the tumor bulk into the surrounding stroma ( Figure 1H). The extent of pressure gradient formation is modulated by the leakiness of the vasculature, and vascular normalization has been shown to decrease IFP 29 . This interstitial flow has implications for disease spread and metastasis as well as therapeutic efficacy. Both of these topics are discussed below.

Barriers to Therapeutic Efficacy
Tumors are notoriously difficult to treat. The combination of leaky vasculature, elevated interstitial fluid pressure, abundant extracellular matrix deposition, adaptive cellular resistance, and CAFtumor cell interactions all lead to attenuated chemotherapeutic efficacy. In general, the obstacles to therapy can be described by physical barriers and cell-mediated barriers. In order to be effective, therapeutic agents must reach the tumor cells; however, several factors mitigate the transport of solute in the TME. Therapy that does reach tumor cells must then overcome barriers of adaptive tumor cell immunity and tumor-associated cell signaling. Together, these factors lead to detrimental patient outcomes and the possibility for disease recurrence.

Physical Barriers to Therapeutic Efficacy
As discussed previously, matrix deposition, rapid cellular growth, and disrupted vasculature lead to elevated pressure within tumors. This pressure leads to two distinct physical barriers. Firstly, while large molecule drugs, such as nanoparticles or antibodies, accumulate in the tumor due to the EPR effect, this accumulation is mostly perivascular. Drug penetration into the tumor from the vasculature is limited by a relatively constant pressure throughout the tumor bulk, meaning that diffusion dominates convection in this context. Since large molecules do not readily diffuse, the distribution is localized, and particles are more likely to get trapped in the matrix 31,32 . However, it should be noted that tumor heterogeneity could yield a moderate amount of drug penetration depending on the distribution of vasculature and the shape of the tumor 33 .
Secondly, at the tumor border the differential in pressure between the tumor bulk and the surrounding stroma yields fluid efflux from the tumor into the surrounding stroma 30,34,35 . This means that any drug that diffuses to the tumor border, or any drug released by vasculature near the border, will leave the tumor via convective transport. As large molecule delivery is dominated by convection, this is especially problematic for many chemotherapeutics delivered via nanoparticles or targeted antibodies.
However, the ECM composition also affects the extent of interstitial flow independently of the interstitial pressure, as permeability and porosity of the tissue of interest affect the likelihood of fluid transport. Neoplastic tissue ECM displays remarkable heterogeneity 36 , and the resistance to convective flow in the interstitium is dependent on the abundance and type of matrix molecule in the particular TME context. For example, transport in chondrosarcoma, a cartilage tumor, vs hepatoma, a liver tumor, would be significantly different due to greater relative abundance of ECM molecules in cartilage as compared to the liver 34 . Mass transfer of a molecule of molecular weight of 10 3 g/mol would be dominated by diffusion in chondrosarcoma and convection in hepatoma due solely to the abundance of ECM components. Increased ECM secretion by fibroblasts in the TME, therefore, significantly hinders solute transport.

Cellular Resistance
In addition to physical barriers, therapeutic efficacy is attenuated by cellular dynamics. Cancer cells can be inherently resistant to drugs due to mutations and adaptations. For example, the above mentioned interstitial fluid flow at the tumor border causes gradients of drugs to form near the tumor edge. These drug gradients provide a selective pressure for adaptive immunity 37 . Several mechanisms have been implicated for multi-drug resistance in cancer cells including upregulation of DNA repair mechanisms, drug metabolism, and drug alteration 38 , however, two of the major effectors of classical drug resistance in cancer is increased drug efflux or decreased drug uptake 39 . ATP binding cassette transporters or ABC transporters have been characterized as major mediators of drug removal from the cell, thus limiting cytotoxic effects. ABC transporter expression is upregulated in multi-drug resistance, and this mechanism non-specifically yields broad resistance to multiple and varied drugs 40 . Additionally, some cancer cells have observable decrease in drug uptake due to specific downregulation of plasma membrane binding proteins 41 .
While, cancer cells can develop resistance independently, cancer associated cells also play a role in tumor resistance to chemotherapy, and fibroblasts, specifically, have been implicated in cancer cell chemotherapy immunity. Marusyk and colleagues describe cancer cell insensitivity to lapatanib due to coculture and proximity to fibroblasts. This resistance is mediated by decreased drug accumulation and drug interaction with fibroblast secreted factors 42 . Cancer associated fibroblasts also secrete pro-survival factors that limit cancer cell apoptosis and could act to protect cancer cells from chemotherapeutic effects 20 .
In the context of neoadjuvant breast cancer therapy, the presence of cancer associated fibroblasts signatures has been directly implicated in poorer survival outcomes 43 ,and targeted therapy against these stromal cells has been shown to improve intratumoral uptake of doxorubicin within in vivo murine models 44 .

Mechanisms of Invasion and Migration
Cancer cells are not stationary. These cells utilize several mechanisms to interact with the ECM and migrate. More severe forms of cancer are characterized by the invasive nature of the cancer cells, as cancer progression by metastasis is mediated by cellular invasion and motility. Gliomas, brain tumors, are composed of particularly invasive cell populations, which are highly metastatic, leading to extremely poor patient outcomes. However, the mechanisms of this migration are relatively unknown, necessitating further research and experimentation.

Cancer Cell Motility
In general, the motility of cancer is mediated by a distinct phenotypic switch from epithelial cell types to mesenchymal cell types. This switch is defined by lack of E-cadherin expression. Cadherins are cell-cell adhesion proteins that allow cells to adhere to surrounding cells, and disruption of this particular cadherin leads to disruption of tight junctions and cell polarity, which frees the cancer cell from close contact with surrounding cells 45 . Disruption of E-cadherin function has also been implicated with the loss of cell-matrix adhesion mechanisms, thus freeing the cancer cells and promoting unhampered migration and associated metastasis 46 .
In addition to a decrease in adhesion, migratory cancer cells also display an alteration of gene expression related to cancer cell invasion. CXCR4, a receptor associated with directed cellular motility toward the complementary ligand CXCL12, is over-expressed in many cancers 47

Interstitial Flow and Increased Invasion
While many mechanisms have been implicated with migration in general, research indicates that interstitial flow enhances the invasion of cancer cells 53,54,55 as compared to normal motility levels. Since interstitial flow is dominant at the tumor border, this increased invasion could be a major factor mediating cancer metastasis 53 . Two mechanisms have been proposed for this increased invasion with interstitial flow. Firstly, it is known that the migratory receptor CXCR4 expression is upregulated in invasive cancer cells 55,47 , however, certain cell populations have both the receptor (CXCR4) and secrete the ligand for that receptor (CXCL12). In flow conditions, it is hypothesized that secreted CXCL12 is convectively driven away from the cell secreting the ligand, thus forming a gradient of CXCL12 in front of the cell in the flow direction ( Figure 2A) 54,55 . If a cell has the receptor for CXCR4, it will migrate in the direction of flow in response to the stimulatory gradient. This directed invasion is termed autologous chemotaxis 54 .
It has also been shown that CD44 expression is upregulated in migratory cancer cells 51 . CD44 mediates migration via signal transduction due to interaction with the ECM, and it has been shown that blocking of this receptor can alleviated increased invasion in the presence of interstitial flow 56 . Rousseau further describes the importance of CD44 as a mediator of dimerization and activation of other receptors, notably CXCR4, through interaction with the ECM and mechanotransduction 57 . It is therefore possible that increased invasion in the presence of interstitial flow is also mediated by a mechanotransduction mechanism and concurrent phosphorylation of CXCR4 ( Figure 2B).

3D In Vitro Models of the Tumor Microenvironment
As the tumor microenvironment is very complex, and there are many factors that influence the progression and spread of cancer, it is necessary to develop relevant experimental models to study cancer in an appropriate environment. 2D culture models are the easiest method to study cancer, however, research has shown that 2D culture models are not always relevant for the unique 3D environment observed in vivo. 58 Murine models are the current standard for experimental and preclinical cancer research, however, these models are complex and are less controllable than in vitro models. Therefore, bridging this experimental gap with simple and relevant 3D culture models provide medium throughput results comparable to in vivo conditions. Several in vitro systems have been developed to study cancer within the 3D context.

Simple 3D Models
Perhaps the simplest 3D models utilized to study cancer are tumor spheroids. Spheroids are 3D masses or cellular clusters within a cell culture dish that are grown from single cell suspensions. These masses replicate many of the properties inherent to in vivo tumors including cell to cell interactions, regions of hypoxia, elevated pH, and natural ECM deposition 42 . Spheroids can also be utilized for coculture studies as other cell types can be incorporated into the tumor spheroids 59 . Spheroids are especially optimal for drug screening studies, as the 3D spheroids display regional heterogeneity replicative of the actual tumor 60 .
Despite several advantages, tumor spheroids are also limited. Firstly, several types of cancer do not form spheroids 60 ; therefore, these models cannot be utilized. Secondly, while ECM components can be incorporated into spheroids, the abundance and type of components are not always replicative of the real tumor environment 61 . Additionally, the morphology and structure of the spheroids is highly variable.
While this is a benefit in terms of developing varied tumor-like structures, it is not advantageous if the composition of the tumor is to be tuned experimentally 60 . Finally, the method of spheroid formation does not account for interstitial flow dynamics.

ECM Hydrogel Models
In order to address the limitations of simple 3D models, hydrogel scaffolds have been developed to replicate the 3D tumor microenvironment. Natural hydrogel scaffolds are characterized by the utilization of physiologically relevant ECM components to create a structurally stable 3D growth medium for cancer cell development and study 62 . These models are robust, as the ECM composition can be tuned to relevant concentrations and the dominant types of ECM components can incorporated to create hydrogels with mechanical properties similar to the tissue of interest 63 . Additionally, cancer associated cell types can be co-cultured in these hydrogels for cell to cell interaction characterization, and the ratios of cells can be tuned to values observed in vivo 64 . Finally, incorporation of fluid on top of these hydrogels introduces a flow rate through the gel comparable to the interstitial flow rate observed within tumors, thus providing an experimental model for the study of interstitial flow effects on cancer cells 54,55,64 .

Agent-Based Models
The same complexities that necessitate robust in vitro models also demand appropriate in silico modeling platforms. In silico research methodologies provide a supplementary data analysis context that can enhance and predict biological phenomena. Agent-based models (ABMs) are particularly suited for biological applications, as they can predict and describe complex spatial and temporal biological interactions 65,66 . In particular, the benefit of agent-based models is the emergent behavior that is possible through incorporation of specific model components. In essence, the sum of many components yields results that would not be predicted by the individual components separately.

An Introduction to Agent-Based Models
Agent-based models are composed of several key components. The most important components of these models are the agents, and there are two types of agents utilized in these models, patches and turtles. Patches are equivalently sized squares that make up the coordinate grid of the agent-based model.
The user defines a spatial framework for the model that has dimensions N x N where N is the number of patches in either a row or a column. For example, a spatial framework of 4 x 4 would have a total of 16 patches in a symmetric square grid. Each patch is unique, and can be utilized within the model. For example, a patch can be programmed to store a specific variable and the user can reference that variable later.
The second component of these models are turtles. Turtles are the objects that interact on top of the patch grid. Turtles are dynamic objects that can move, grow, shrink, die, etc based on the desired outcome and the programmed rules of the model. Turtles can also store variables, like the patches, and the turtles can also receive variables from patches beneath them or transfer them to other turtles. The user can specify as many types of turtles as desired, and each individual turtle or each type of turtle can be programmed with different governing behavior.
The behavior of the agents is updated at successive intervals determined by the user and outputted to a display screen. Much like individual frames in a movie or animation, the location, behavior, size, etc of the turtles and/or patches is updated dynamically according to the specific programming methods. This yields a viewable dynamic model with temporal aspects controlled by the degree of change between frames. For example, if the user desired to create an agent-based model of a car moving rightward at a speed of 20mph, the user could specify a patch grid of 20 x 20 with each patch being 5 miles long. At each individual frame, the user specified program could tell the car turtle to move rightward by a one patch distance on the coordinate grid. If the time between two frames was specified as 15 minutes, iterating this method for four frames would move the car rightward by 20 miles in 60 minutes or 20mph.
This method could then be iterated as necessary.

Agent-Based Models of Cellular Motility
As agent-based models are inherently spatial, agent-based models characterizing cellular motility have been previously utilized. Emmonet et al. describe an agent-based model of bacterial chemotaxis with and without chemo-attractant using the predictive equations developed from classical microbiology. Here, motility in the absence of attractant is defined by random motility synonymous to chemical diffusion 67 .
Enderling and colleagues developed an agent-based model of tumor development that incorporates migration programming methods in order to predict tumor morphology 68 . Interestingly, Heiko and Enderling also explore tumor development as a result of random migration and chemoattractant-mediated migration. Here the morphology is dependent on the directionality of migration 69 .

Agent-Based Models in Pharmacology
While agent-based models are not utilized as readily for the dynamics of chemotherapy, as other modeling methodologies are more ubiquitous, model-based descriptive techniques are being utilized in the area of pharmacology. These models help to facilitate hypothesis testing and characterize complex dynamics in a simple and replicable platform 70 . For example, an agent-based model has been used to predict the effect of novel drug targets on the fungal infections based on either inhibition of fungal growth or phenotypic transition 71 . Agent-based models have also been utilized concurrently with descriptive pharmacokinetics and pharmacodynamics to predict the effect of drugs at the cellular level 72 . In general, however, it appears that many models are context and hypothesis specific, and the utilization of agentbased models to describe treatment dynamics is an emerging field 70 .

Motivation and Conclusion
The tumor microenvironment is very diverse; the combination of cellular factors, ECM, biophysical forces, and therapeutic barriers all contribute to complexity, which is difficult to characterize and difficult to study. Therefore, robust research methods must be employed to elucidate the effects of this unique environment on the development and spread of cancer. In vitro hydrogel systems utilized concurrently with robust in silico agent-based models facilitate the discovery of TME specific tumor outcomes in the context of interstitial flow and chemotherapy treatment. Here we describe two agentbased models that predict the effects of interstitial flow on brain cancer invasion and the effects of fibroblasts co-culture on cancer cell viability in the presences of chemotherapy.
While agent-based models have been utilized to describe cellular motility 67, 68, 69 , these models have not incorporated cellular invasion as a specific outcome measure. Additionally, interstitial flow was not a factor in any of the above models. Here we describe a novel agent-based model that replicates our in vitro invasion assay and predicts mechanisms of interstitial flow-mediated increased invasion of glioma cells. models provide insight into complex systems, and can be utilized in this context 70 . Here we describe a novel agent-based model that predicts the region specific response of cancer cells to chemotherapy within the tumor border to stroma transition region based on the kinetics of therapeutic transport and the protective effect of fibroblasts. We have concurrently developed a layered hydrogel system of the tumor border to stroma transition zone. Regional variability in this in vitro model, in response to chemotherapy, further elucidate the importance of the TME in mediating cancer cell resistance to chemotherapy.

Introduction
Cancer cell motility is linked to decreased patient survival due to metastasis and evasion of treatment. Glioblastoma, brain cancer, is particularly deadly as the cancer cells are characteristically invasive 73,74 . This effect is enhanced by interstitial flow at the tumor border. Interstitial flow develops due to a pressure differential between the tumor bulk and the surrounding tissue, which leads to fold increases in flow velocities 75 . This interstitial flow has been shown to increase the invasion of the cancer cells 55 .
Since the standard of care involves surgical resection of the primary bulk, which often misses cells at the invasive edge of the tumor, this population is optimally primed for further invasion 76 .
Increased invasion due to interstitial flow has been predicted by two distinct mechanisms. Firstly, autologous chemotaxis, mediated by cell populations expressing both the migratory receptor CXCR4 and secreting the ligand to that receptor CXCL12, has been proposed by Munson et al. 55 as a dominating mechanism for increased invasion. These cells secrete CXCL12 which is convectively drive by interstitial flow. This mediates the creation of autologous gradients of CXCL12, which the cell can respond to with its CXCR4 receptor.
Tarbell and colleagues have also described an alternative mechanism in several cancer cell lines 56 . This mechanism involves CD44 mechanotransduction, where cancer cells respond to interstitial flow indirectly by mechanical sensing by the CD44 receptor. It has further been shown that CD44 can phosphorylate other membrane proteins, including CXCR4, which indicates that the increase in flowmediated invasion could also be mediated by CXCR4 phosphorylation by CD44. Here, both in vitro and in silico experimental models have been utilized to describe flow-mediated invasion, and they have been utilized to predict the dominating mechanism of this motility. Additionally, specific subpopulations have been implicated in mediation both mechanotransduction and autologous chemotaxis.

Glioma Stem Cell Culture
G34 cells were primary derived as described and were provided as a kind gift from Dr Jakub Godlewski (Cleveland Clinic) and have been previously published and characterized. 77 Cells were cultured following the previously described protocol. 77 Briefly, GSCs were cultured in Neurobasal media (ThermoFisher) with N2 and B27 without vitamin A supplements (ThermoFisher), human recombinant bFGF and EGF (50 ng/mL, ThermoFisher), Glutamax (ThermoFisher) and Penicillin-Streptomycin

Three-Dimensonal (3D) In Vitro Interstitial Flow Model
Hyaluronan is the most abundant extracellular matrix element in the brain and was therefore chosen to be a primary matrix component of our 3D in vitro model. 24 This matrix gel replicates both the stiffness and composition of the brain, but also allows us to achieve comparable flows that have been reported in vivo. 55,78 Glioblastoma stem cells were resuspended in hyaluronan-collagen matrix at a density of 1 million cells/ml as previously described. 55 Briefly, cells were resuspended in 50 µl of 0.2% hyaluronan (ESIBIO) /0.12% rat tail collagen I (Corning) and cross-linked with PEGDA (ESIBIO). This solution was applied to an 8-micron pore 96-well tissue culture insert system (Corning), and allowed to gel for 2 hours. 55  [Methods modified from Kingsmore et al.]

Agent-Based Cell Migration Model Construction
The agent-based model (ABM) was constructed using the Repast Java framework and the Relogo Java package. Cell types were determined by the permutation of all receptor and/or ligand conditions (Table   1).  Semi-directional invasion was implemented by introducing one movement in the flow direction for every 13 random movements ( Figure 5). Each tic of the scheduling method was defined as one minute, and invasion percentages were produced after an 18 hour timespan. The velocity of the cells was defined as 7 µm/min as calculated by Munson et al. previously, and the 2D coordinate grid was initiated to model a vertical cross-section of a tissue culture insert with dimensions of 4.5 mm x 3.5 mm (7 by 9 patch box grid with side length equal to 0.5mm). 55 The number of total cells (2800) was established by the number of cells in a finite slice of an experimental gel. The percentage of each type of receptor positive/negative cell was determined from flow cytometry data, and the baseline brain cancer model was created from the range of the percentages for all glioma stem cells.

Figure 5| Specific ABM methods flow chart.
Programming methods for each type of motility utilized in the agent-based model.

Random Motility Validation
As random motility was an essential component of the computational model, the random motility programming method, described by facing cells in random directions at each scheduling tic and moving them forward a distance related to their average velocity, was validated using protocols established to describe bacterial motility in the absence of attractant.
In 1967 Adler and Dahl pioneered a method for determining bacterial motility by relating displacement to a motility coefficient synonymous to a diffusion coefficient. 80 The authors described a general equation for the diffusion of glucose in a gel based on experimentally observed behaviors where 0 is the initial concentration, c is the concentration at distance x measured from the origin at time t, and D is the diffusion coefficient of the substance. The authors noted that, in the case of their experimental setup, the first term of the equation was negligible, and the equation could be written as where the variables are described above. The authors additionally noted that the migration of bacteria within gels matched the trends observed with glucose diffusion. This lead to an equation of squared displacement of bacteria where M is motility coefficient synonymous to a diffusion coefficient, c is the minimum number of detectable bacteria, and 0 is the number of bacteria added to the assay. where C is the concentration of cells at radial distance r from the origin at time t, 0 is the number of cells present in the model, and D is the random motility coefficient. This equation can be rearranged into the form ln ( which is almost identical to the equation determined by Adler and Dahl except for an additional term (3/2) to account for the 3 rd dimensionality of the agent-based model. This equation simplifies to assuming that the second term in Equation (5) is negligible, as stated by Adler and Dahl. This expression is equivalent to where <r 2 > is ensemble mean squared displacement of all cells and d is the dimensionality of the system.
where 2 (∆ ) is the mean squared displacement at a particular time step 2 (∆ ) is the squared displacement of a cell from the origin at a particular time step. It should be noted that the squared displacement 2 (∆ ) is calculated as follows in order to normalize displacement of the cells at the zero coordinate in an x-y grid, where xt is the x coordinate at time t, yt is the y coordinate, x0 is the x coordinate at time zero, and y0 is the y coordinate at

Statistical Analysis
Statistical analyses were run using Graphpad Prism. Paired t-tests and two-way ANOVA were used for analysis of same subject groups. Unpaired t-tests and two-way ANOVA were used for analysis of independent experimental groups and computational data. . Multivariate ANOVA was used to compare generated curves for sensitivity analyses of the agent-based model. For all data, *p<0.05, **p<0.01, and ***p<0.001, and graphs are given as mean +/-standard error of the mean.

Simulated random motility is representative of experimental cellular motility in the absence of stimulant
Comparison of the simple random bacterial ABM to the AgentCell model developed by Emonet et al. showed that at a simulated velocity of 0.058 mm/s, the random motility coefficient was determined to be 3.99 *10 -6 cm^2/s ( Figure 6A). This was compared to the value of 4.21*10 -6 cm^2/s reported by to size limitations of output excel files, and the difference in slope is accounted for by the dimensionality difference (2D vs 3D). The histogram of displacement r was comparable to the normal trend observed in the AgentCell model ( Figure 6B). The trend followed a normal distribution, with mean of 1.39 mm, as confirmed by a normal probability plot. After validation of the random motility with bacterial cells, the same method was used to determine the motility coefficient for glioma cells in the absence of fluid flow for both simulated displacements ( Figure 7A) and experimental cell displacement ( Figure 7C). The agent model motility coefficient obtained from calculating the mean squared displacement for all cells (Equations 7-9) at a cell velocity of 1.077 um/s was determined to be 2.39*10 -9 cm 2 /s. This was compared to the experimentally determined motility coefficient of 2.38*10 -9 cm 2 /s. The percentage of cells above 1.2 um in the positive y direction was also determined for the agent model ( Figure 7B) and experimental conditions ( Figure 7D).

Mechanotransduction and autologous chemotaxis mediate flow responsive invasion
After validation of random motility, blocking studies were performed both in vitro and in silico to simulate disruption of specific receptors (CD44/CXCR4) and the specific ligand (CXCR4) hypothesized to be important in flow stimulated invasion. Migration simulations were performed using the percentages of expected cell types based on flow cytometry data ( Figure 8A). Simulated invasion and experimental invasion results were comparable ( Figure 8B). facilitated mechanotransduction and CXCR4-CXCL12 migratory signaling, as blocking of either of these mechanisms removed the flow response. Neither mechanism was sufficient to elicit flow responsive invasion alone, but taken together, the effect is observe. It was further hypothesized that autologous chemotaxis was the major effector of this increased invasion in relation to the CXCR4/CXCL12 receptorligand paired signaling.

Two subpopulations account for increased motility in the presence of interstitial flow
Since it was hypothesized that autologous chemotaxis and mechanotransduction were both necessary for flow responsive invasion, a sensitivity analysis of the ABM was performed to test the importance of specific cell types on the overall invasion response. Specifically, the cell types of relevance to autologous chemotaxis (CXCR4+/CXCL12+, Figure 9A

Conclusion
An agent-based model of glioma cell invasion was developed using the methods described by others, and its utilization confirms the propriety of agent-based models utilizing motility algorithms.
Random motility was validated using established equations of cellular motility. Incorporation of both directed motility and random motility rules into the agent-based model elucidated the importance of both mechanisms to increased glioma invasion in the presence of interstitial flow. Results indicate that two important subpopulations could be driving this response, and selection of either population could increase overall likelihood of metastatic spread.

Introduction
Chemotherapy is still a near-ubiquitous treatment approach for multiple forms of solid tumor cancers. Though it works in many cases, chemotherapy often fails against cancer, resulting in poor prognosis across multiple forms of the disease. Doxorubicin (DOX) is a commonly used chemotherapy against multiple cancers including breast, bladder, ovarian, and lung, alone or in concert with other treatments. [82][83][84] Doxorubicin is a common chemotherapeutic delivered prior to, or after surgery in cases of the deadliest form of breast cancer, triple negative. In addition to its clinical relevance, as a drug with wellunderstood pharmacokinetics, doxorubicin provides an optimal model drug for probing dynamics of chemotherapeutic treatment.
Systemic chemotherapeutic delivery has been limited, first and foremost, by transport restrictions to and within tumors 32 . Specifically, the tumor microenvironment (TME), defined as the tumor cells, tumor-associated cells, extracellular matrix, and biomechanical forces 1,2 , that interact in and around the tumor, provides two distinct barriers to drug delivery. Firstly, delivery of chemotherapeutic agents to the tumor bulk from the circulation is attenuated by interstitial fluid pressure, which is elevated in tumors. 85 This pressure can limit the transvascular movement of small molecule drugs into the interstitial spaces leading to retention of therapy, particularly larger molecules, near the tumor-associated blood vessels. Drug that does move into the interstitial tumor space is subject to a number of further constraints that limit the transport through the tissue. Limited transport of therapeutic agents results from matrix deposition and crosslinking by cancer associated fibroblasts, uptake by therapeutics by stromal cells, and reduced overall void space due to unrestricted cell growth 1 . This limited drug distribution is thought to particularly reduce therapeutic access to invading cancer cells at the tumor border. These cells are especially deadly as they are thought to be responsible for subsequent metastasis of cancers, which is a leading cause of death.
Interestingly, at these border regions of the tumors, where cells are invading, there is increased interstitial fluid flow. The higher interstitial fluid pressure in the tumor bulk relative to the normal pressure in the surrounding stroma yields an efflux of fluid from the tumor to the surrounding tissue 76 . These forces are known to alter cellular invasion, promoting movement towards draining lymphatics. However, the interaction of these particular forces with chemotherapeutic transport and delivery is not explored in the context of this complex transitional microenvironment.
Secondly, intercellular interactions in the tumor microenvironment can contribute to reduced therapeutic response. Several TME-mediated factors have been implicated in the development of chemotherapy resistance within solid tumors 86 . These factors include hypoxia 87 , reduced pH 88 , nutrient deprivation 89 , adhesion-mediated resistance 90,91 , and drug gradient formation 37 . However, cancer-stromal interactions, specifically between cancer cells and fibroblasts in breast cancer, appear to be a dominating factor in acquired chemotherapy resistance 92 . Cancer associated fibroblasts secrete pro-survival factors that limit cancer cell apoptosis and could act to protect cancer cells from chemotherapeutic effects 20 . In the context of neoadjuvant breast cancer therapy, the presence of cancer associated fibroblasts signatures has been directly implicated in poorer survival outcomes 43 ,and targeted therapy against these stromal cells has been shown to improve intratumoral uptake of doxorubicin within in vivo murine models 44 .
In relation to these barriers to chemotherapy treatment, the tumor border is a unique environment, as stromal interactions, chemotherapy gradients, and interstitial flow are all present in this region.
Interstitial fluid velocity and pressure are the greatest at the tumor border 29 , thus mediating convectiondriven chemotherapy transport through this region 35 . Additionally, gradients of doxorubicin have been observed at the breast tumor border in vivo 37 , providing a potential for adaptive chemotherapy resistance.
Finally, cancer cell propinquity to stromal fibroblasts, coupled with zones of cellular transition 93 from the tumor to the surrounding stroma, yields stromal heterogeneity. At the tumor border, the tissue transitions from regions of very few fibroblasts in the tumor bulk to regions of very few cancer cells relative to stromal cells; therefore, this transition region is defined by cellular gradients, which could mediate further cancer cell insensitivity to chemotherapy due to fibroblast interactions. Collectively, this distinct environment couples several TME-specific factors that could lead to chemotherapy resistance and decreased patient survival.
Due to the complexity of this region, robust in vitro and in silico models are necessary to probe the effects of therapy coupled with interstitial flow. 3D collagen hydrogels have been utilized in several cancer contexts as in vitro models of the tumor microenvironment 76 . These models can incorporate extracellular matrix proteins and pressure driven flow, thus replicating biophysical parameters inherent to in vivo tumors and providing a culture medium that is more replicative of real tissues than 2D culture systems 94 . In silico research methodologies provide a supplementary data analysis context that can enhance and predict biological phenomena. Agent-based models (ABMs) are particularly suited for biological applications, as they can predict and describe complex spatial and temporal biological interactions 65,66 .
Here, we establish in vitro and in silico models that elucidate the effects of fluid and solute transport, cellular heterogeneity, and fibroblast interactions on breast cancer viability following doxorubicin treatment. Our study utilizes a novel 3D in vitro tumor bulk to stroma transition model (TSTM), as well as concurrent 2D culture systems, to predict the regional variations in viability that occur within the microenvironment at the tumor border. Additionally, in silico methodologies predict the dominant fluid dynamic properties that influence doxorubicin treatment efficacy within the tumor border transitional region. Together our data provide evidence of a unique fibroblast protective effect that yields varied resistance to chemotherapy in a cancer to fibroblast ratiometric dependent manner, which is sustained in several experimental models. These findings illuminate the importance of regional TME heterogeneity in selecting for viable populations that could be impacting breast cancer progression.

MDA-MB-231 human breast adenocarcinoma (Tumor cells: TCs) and human dermal fibroblasts
(Fibroblasts: Fbs) were obtained from the American Type Culture Collection (ATCC). Both cell types were cultured in Dulbecco's modified Eagle's medium (DMEM, Gibco), supplemented with 10% fetal bovine serum (FBS, Seradigm). Cells were passaged weekly and grown at 37°C in a sterile incubator (5% CO2 and 95% Oxygen) in gamma irradiated tissue culture treated flasks.

Two-Dimensional (2D) Hanging Well Coculture Assay
MDA-MB-231 cells were seeded into a 24 well tissue culture companion plate for hanging cell culture inserts (VWR International) at a density of 40,000 cells per well. HDFs were seeded into hanging culture inserts (VWR international), with 1.0 micron membrane pore size, at a density of 10,000 cells per insert. For the control condition, HDF seeding was neglected. Cells were grown in serum free DMEM for 24 hours. Prior to introduction of chemotherapy, the cell culture insert was removed for the conditioned experimental group. Doxorubicin was introduced, to bring the final concentration to 10uM. TC percent live and doxorubicin accumulation were determined as described.

Two-Dimensional (2D) Ratiometric Coculture Assay Constant Total Cell Number Ratio Experiment
For the first experiment, HDFs and Cell Tracker (ThermoFisher) deep red labelled MDA-MB-231 cells were introduced into the same well of a 48 well tissue culture treated cell culture plate at varied ratios of cancer cells to fibroblasts (4:1, 2:1, 1:1, 1:2, 1:4). Here, the total cell density was held constant at 30,000 cells, and cell numbers were varied internally. For example, a ratio of 4:1 had 24,000 cancer cells and 6,000 fibroblasts, while a ratio of 2:1 had 20,000 cancer cells and 10,000 fibroblasts. Control conditions were seeded at cancer cell densities comparable to experimental conditions, however, fibroblasts were not introduced. The cells were incubated for 24 hours at 37°C, followed by introduction of doxorubicin HCL diluted in serum free DMEM (10uM).

Doxorubicin Accumulation Analysis
For doxorubicin accumulation experiments, the cells were trypsinized (0.25% trypsin, Gibco) at two hour time points (at 6 hours post doxorubicin introduction for the 2D hanging well experiments), centrifuged at 2000 rpm for two minutes, and subsequently lysed with RIPA Buffer (Thermo Scientific) for 15 minutes on ice. Cellular lysate fluorescence intensity was determined using a fluorescent plate reader (Omega FLUOstar) at 495nm excitation and 590nm emission. Fluorescence intensity was compared to a doxorubicin standard curve for concentration determination. The gels were removed from the inserts and imaged using fluorescent microscopy (EVOS FL).

Three-Dimensional (3D) Homogenous In vitro Interstitial Flow Model
Five images were taken per experimental technical replicate and averaged for statistical analysis. Dead cancer cells were determined by colocalization of blue, deep red, and green fluorescent markers.

Generation of IC50 Curves
HDFs and MDA-MB-231 cells were seeded into separate 96 well tissue culture plates at a density of 10,000 cells per well. Cells were incubated for 24 hours at 37°C. Afterward, doxorubicin solutions at varying concentrations diluted from 100uM were introduced into appropriate wells. Dead cells were determined as previously described. IC50 values were determined using Matlab curve fitting algorithms.

Development of Concentration Gradients using Comsol
Spatial concentration profiles were developed using a finite element based, partial differential boundary solver in Comsol Multiphysics. 2D Geometry of the hydrogel spatial grid (4.5 mm by 3.  Table 2). The velocity in the hydrogel is further described by where ( 2 ) is the matrix permeability coefficent, ( * ) is the dynamic viscosity of the fluid, and ( ) is the pressure. The permeability coefficient ( Table 2) (Table 1).

Agent-based Model Construction
Comsol time-dependent convection diffusion concentration gradients (Comsol Multiphysics 5.1) were developed using Darcy's law of fluid flow coupled with simple diffusion dynamics, representative of doxorubicin transport within in vitro breast mimetic collagen hydrogel/BME hydrogels 95 (Supplemental Methods). Physical model parameters were determined from experimental methods and literature values (Table 3). Concentration profiles were outputted as text files for use in agent-based model (ABM) simulations (Figure 10a).

Table 3| Agent-based model equations and parameters
The agent-based model (ABM) was constructed using the Repast Java framework and the Relogo Java package. A 2D coordinate grid of 64 x 64 patches was initiated to represent a 4.0mm x 3.5mm finite slice of experimental gel 54 . For homogenous gel simulations, the total number of cells was 2800, and the fibroblasts and/or cancer cells were spawned randomly within the gel at the user specified ratios ( Figure   10a). For the layered tumor transition model, the grid space was split into 5 equivalent segments ( Figure   10c). Fibroblasts and cancer cells were spawned randomly within the sections at user specified ratios for a total cell number of 2800. For the single culture tumor transition model, cancer cells were spawned within segments at densities comparable to experimental conditions. The total cell number was 1400 (Figure 9b).   Table 2).
The average concentration encountered by the cells was updated based on the stored value. At specified time intervals, the expected live percentage, based on IC50 curves and the average concentration encountered by the cells, was determined. This value was scaled to account for cellular density (Table 2).
If the modeling condition accounted for both fibroblast and cancer cells, the expected live percentage was scaled to account for the fibroblast protective effect based on the ratio of cancer cells to fibroblasts. The cells were removed or kept in the simulation after comparing a randomly generated number to the expected live percentage. Percent live of remaining cells was then calculated.

Sample Selection
Patient samples were accessed through the University of Virginia Biorepository and Tissue Research Facility. These samples were selected from patients with a definitive diagnosis of node-negative breast cancer and who received no treatment prior to tumor resection. Samples were de-identified before use. All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional review board of the University of Virginia and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards.

Immunohistochemistry
Formalin-fixed, paraffin-embedded sections were deparaffinized with xylene and rehydrated in graded ethanols and citrate-based antigen retrieval was performed (Vector Labs, Burlingame, CA). Samples were permeabilized (0.01% Triton) and blocked in goat serum. Based on markers previously established in the literature, breast cancer cells were identified by anti-pan-cytokeratin staining and cancer-associated fibroblasts were identified by anti-alpha-smooth muscle actin staining. Samples were incubated with pancytokeratin antibody (Thermoscientific) followed by secondary Cy5-goat anti-mouse (Thermoscientific).
These steps were then repeated for the TRITC-conjugated alpha-smooth muscle actin antibody (ebioscience). The samples were incubated with DAPI (Sigma-Aldrich, St. Louis, MO) and mounted with Fluoromount-G (SouthernBiotech). All antibodies were used at dilutions recommended by the manufacturer for paraffin-embedded tissues. Stained slides were imaged with an EVOS fluorescent microscope (Thermoscientific). Random non-overlapping 856 × 476 µm (407, 465 µm2) regions at the tumor-stroma border were selected for imaging. Images were processed using ImageJ and Photoshop.

Statistical Analysis
Statistical analyses were run using Graphpad Prism. Paired t-tests and two-way ANOVA were used for analysis of same subject groups. Unpaired t-tests and two-way ANOVA were used for analysis of independent experimental groups and computational data. MANOVA analysis using the SPSS software package was utilized for normalized distance comparisons within experimental gels for both computational and experimental conditions. All assays were performed with a minimum of three biological replicates.
P<0.05 was considered statistically significant for all statistical tests.

Fibroblasts reduce cancer cell death in response to doxorubicin and lowered levels of doxorubicin accumulation
The effects of contact independent fibroblast signaling on cancer cell response to doxorubicin were first tested by treating MDA-MB-231 cells with either HDF conditioned media or comparable control media ( Figure 12A) in 2D for 24 hours. Following doxorubicin treatment, tumor cells grown in conditioned media had a greater percentage of live cells compared to unconditioned controls ( Figure 12B). Furthermore, tumor cells treated with either Fb-conditioned or unconditioned media showed no difference in percent live tumor cells in the absence of doxorubicin. Notably, at 2 and 6 hours, the mass of doxorubicin present in cell lysate was significantly decreased following pre-incubation with fb-conditioned media prior to treatment compared to unconditioned media controls ( Figure 12C).
To further test contact independent fibroblast signaling, fibroblasts were seeded on a hanging porous insert with tumor cells seeded onto the bottom of the well plate at a ratio of 4:1 tumor cells to fibroblasts ( Figure 12D). While preconditioning with fibroblasts showed no effect, co-culturing with fibroblasts during treatment yielded a significant increase in the percentage of live tumor cells compared to However, there was no difference in doxorubicin accumulation in the conditioned group as compared to the tumor cells alone control. Examination of the fibroblasts also revealed that the co-culture condition not only led to changes in the tumor cell population, but also led to opposite effects in fibroblasts. Co-culture led to increased death of fibroblasts and increased uptake of doxorubicin by fibroblasts ( Figure S1). This decrease in doxorubicin accumulation correlated with increased percentage of live cells in comparable conditions after doxorubicin treatment, indicating the potential for fibroblast-mediated alterations in tumor cell drug accumulation and reduced treatment response.

Cancer cell response to doxorubicin is dependent on the ratio of tumor cells:fibroblasts
The tumor border includes a transition from the tumor bulk, with relatively few fibroblasts, to the surrounding stroma, with very few tumor cells. Therefore, it was desirable to determine the extent of the fibroblast protective effect within the context of varied cancer cell to fibroblast ratios and varied total cell densities, of relevance to the tumor-stroma transition (TST) zone.
We  To further probe this phenomenon, the number of tumor cells was held constant, and the cell ratio was altered by adjusting only the total number of fibroblasts in coculture ( Figure 13C), prior to doxorubicin treatment. A control group was also made by seeding tumor cells at the same number as all other experimental conditions without the addition of fibroblasts. A comparable viability effect was observed as seen in the hanging well cultures: A ratio of 4:1 tumor cells to fibroblasts yielded the greatest viability of tumor cells and the greatest difference in viability as compared to the control condition ( Figure 13D). These data indicate that the fibroblast protective effect is not dependent on the total number of tumor cells; rather, it depends on the ratio of tumor cells to fibroblasts. Additionally, while fibroblast-derived factors need to be present to observe the effect, the ratio of tumor cells to fibroblasts and not the total number of fibroblasts is the determining factor.

Ratiometric effects are conserved when cells are cultured in a breast-mimetic 3D microenvironment
Since the in vivo microenvironment is a 3-Dimensional space with which the cells can readily interact, and it has been shown that dimensionality can affect therapeutic outcomes in cancer, comparable in vitro hydrogel doxorubicin treatment experiments were developed to extend the coculture experimental findings. Collagen-basement membrane hydrogels were developed with varying ratios of tumor cells and fibroblasts found to be relevant (4:1, 1:1, 1:4) in 2D ( Figure 14A,B). Control gels were seeded with tumor cells, alone, at the same total cell number as the coculture hydrogels ( Figure 4B). Media with and without doxorubicin was applied to the gels for 18 hours. Total tumor cell viability in the absence of doxorubicin was approximately 90% for all conditions ( Figure S3, 13C). Analysis of the total tumor cell viability within the gels after treatment indicated that a ratio of 4:1 tumor cells to fibroblasts yielded a significantly higher viability than the comparable single culture condition.

Gradients of chemotherapy form across invading edges of tumors
In addition to the cellular heterogeneity at the tumor edge, increased fluid pressure at the tumor border drives interstitial fluid flow into the adjacent stroma generating a gradient of chemotherapeutic 37 .
Collagen hydrogels have been utilized in the past to mimic this region. Here, we developed in silico theoretical concentration profiles within collagen hydrogel mimetics based on the geometry of the hydrogels, Darcy's law of fluid transport, and the general diffusion equation (Table 1). Because the fluid velocity is very low (~0.5 um/s), concentration gradients develop between the top and the bottom of the gel ( Figure 15A). Notably, the concentration of the drug was greatest at the top of the gel relative to the bottom.
However, it should be noted that the concentration became homogenous in the hydrogel after 3 hours of simulated flow. Our ABM predicted more significant results than our 3D hydrogel experiments, however, the viability trend was comparable between in vitro and in silico models. Viability of tumor cells in each cocultured condition was approximately half of the ABM predicted viability percentage. This is likely due to a reduction in the fibroblast protective effect within 3D systems, as compared to 2D culture system, as a result of the extracellular matrix components and greater dispersal of cells.

Development of an in silico tumor to stromal transition (TST) model ABM shows response across the invasive edge of tumors
Following confirmation of the fibroblast protective effect in 3D, an in silico ABM tumor to stroma transition (TST) model (Figure 11), representing the in vivo tumor bulk to stroma transition area ( Figure   16A), was utilized to predict the region specific effects of doxorubicin treatment within a heterogeneous 3D-hydrogel

A corresponding in vitro model of the invasive edge indicates similar region specific cancer cell response
A penta-layered hydrogel in vitro TST model was created by successively depositing hydrogel solutions with increasing cancer cell to fibroblast ratios from bottom to top of a cell culture insert ( Figure   16C). Tumor cells in alternating layers were labelled with different cell tracker dyes ( Figure 16D). Distinct regions were confirmed, via confocal imaging, prior to introduction of flow. These regions were confirmed following flow application ( Figure 16D).
The region-specific viability within the in vitro TST model confirmed the trends predicted by the comparable ABM ( Figure 16E), with the addition of fibroblasts increasing viability overall. However, fibroblasts did not confer significantly increased viability in the lowest layer of the TST model as was  Altering the diffusion coefficient, analogous to changing the size of the drug, while maintaining convective parameters, yielded no significant difference in the general region specific viability of tumor cells ( Figure 17C), except for a 100-fold increase in the diffusion coefficient. Similarly, overall percent live tumor cells within the hydrogels was only significantly altered from the control condition in the case of 100-fold increase in the diffusion coefficient ( Figure 17D). Analysis of the concentration profiles indicated that convection was dominating diffusion in this context ( Figure 18). Alterations of the permeability coefficient, analogous to altering the matrix density, and associated alterations in velocity ( Figure 18D), while maintaining the diffusion parameters, yielded varied alterations in region specific and overall cancer cell viability ( Figure 17E). While the viability at the top of the gels was always greater than the bottom of the gel, for comparable conditions, increasing the permeability decreased the overall viability ( Figure 17F). At low values of permeability (0.01K and 0.1K) there was no difference in region specific or overall viability, however, the viability was not zero. This was due to the fact that the gel becomes saturated with drug almost instantly at these values; therefore, the viability effects become dominated by the drug dosage and the fibroblast protective effect. At very low values of permeability, the region specific and overall viability was comparable to diffusion only conditions, indicating that diffusion was dominating convection in these conditions.

Mechanisms behind fibroblast induced chemotherapy resistance
Multidrug resistance in cancers is a documented phenomenon, and can be meditated by several cellular mechanisms; however, anthracycline resistance is dominated primarily by classical multidrug resistance 40 . Classical multidrug resistance is described by acquired resistance to chemotherapy by lowered intracellular drug concentration; this can be mediated by increased drug efflux or decreased drug uptake.
For anthracyclines, such as doxorubicin, intracellular drug accumulation is decreased by increased drug efflux via ATP-dependent transporters 96 97,98 , and this saturation occurs even at a dosed doxorubicin concentration as low as 5μM 97 . Therefore, it is not likely that the reduced viability in coculture conditions is simply due to uptake by stromal fibroblasts. That being said, localized effects of this accumulation could have further reaching effects on total tumor treatment response based on observations of gradient formation of doxorubicin in patients and resistance to treatment 99 .
One of our more perplexing findings is that there is a protective effect at greater ratios of tumor cells to fibroblasts, yet mechanisms underlying this response are more challenging to parse. From our data, we interpret that a concentration-dependent, fibroblast-derived signaling molecule is necessary for the protective effect to occur, and this effect is observed at higher ratios of tumor cells to fibroblasts, regardless of the overall cell numbers. In other words, increasing the number of fibroblasts relative to tumor cells decreases the protective effect. Fibroblast activation, a major hallmark of cancer-associated fibroblasts, has been attributed to TGFβ, reactive oxygen species, and modifications to the extracellular matrix, among other mechanisms. Activation of fibroblasts leads to further secretion of TGFβ, as well as increased secretion of other cytokines such as IL-6 and extracellular matrix molecules such as Tenascin C 100 . TGFβ and IL-6 both increase resistance to anthracycline chemotherapy in breast cancer 101,102 , and Tenascin C is associated with chemoresistance in lung cancer 103 . It's possible that in our system, the ratiometric effect is due to a unique balance between signaling molecules from cancer cells to fibroblasts on a per cell basis needing to be higher than signaling molecules from fibroblasts to cancer cells, but the effect is amplified regardless when cells are touching or sharing an extracellular matrix as in our co-culture systems.
Fibroblasts have also been associated with creating niches in which cancer stem cells, a particularly drug resistant population of cancer cells, reside and proliferate 104 . The ratio effect may be due to finding a particular sweet spot for creation of such a niche for these chemoresistant cells.
Regardless, the fibroblast protective effect is a dominating factor for increased cancer cell viability in our 2D and 3D in vitro systems. In silico homogenous gel models predict that region specific cancer cell viability is affected by chemotherapy gradients; however, our in vitro and in silico models predict that the overall viability is dominated by the dependent fibroblast protective effect. This indicates that viability trends within homogenous tumor regions are predicted by concentration gradients, yet the magnitude of the viability percentage may be dominated by the relative abundance of fibroblasts.
However, within regions of cellular heterogeneity, the region-specific viability response to doxorubicin treatment is not solely predicted by the average doxorubicin concentration that the tumor cells experience ( Figure 5). While tumor cells at the top of the hydrogels, representative of the tumor bulk region, experience a greater average concentration of drug as compared to the bottom of the gel, representative of normal stroma, the viability in this region is greater. This indicates a potentially important subpopulation for cancer progression and recurrence.

Drug transport in the tumor microenvironment
Fluid dynamics can alter drug concentration gradients and affect region specific and overall viability trends in the tumor border region 95,34 . Reversed flow enhances the cancer cell viability in the bulk region and simultaneously decreases the viability in the stromal region, while additionally increasing the overall viability. This condition is of import, as the direction of fluid velocity in vivo is not known, and regional tumor heterogeneity could yield solute influx into the tumor from the surrounding stroma 34 .
Diffusion, in the absence of convection, does not alter viability at 18 hours of simulation due to lack of solute penetration into the hydrogel. These results indicate that convection is a dominating factor for concentration profile development and corresponding cancer viability response 37 . This is further confirmed by alterations of the diffusion coefficient, while simultaneously keeping convection constant. In these simulations, the diffusion coefficient only affects the viability at very large values. This confirms that convection is the driving force for doxorubicin concentration development, and diffusion dominates only when molecules are very small.
These results are comparable to those reported by Jain and colleagues 29,30,34,35 ; the authors highlight the importance of convection and diffusion for solute transport in neoplastic tissues. Specifically, convection is dominant in these tissues when the molecular weight of the solute is significantly large, while diffusion dominates with very small particles 32 . This rationale was further used to explore the transport effects of large molecules within the tumor interstitium and corresponding implications for antibody-based drug delivery 30 . The results indicate that high molecular weight drug efficacy, in whole tumors, is attenuated by extravasation from the tumor bulk. While doxorubicin would not be considered a large molecule drug, it is well within the size range for convection driven flow in neoplastic tissues 34 .
However, convection dominance is also mediated by the permeability of the tissue of interest. As breast tumor heterogeneity affects the permeability and corresponding fluid velocity 105 in vivo, characterizations of alterations in permeability were explored to determine the effect of tumor permeability heterogeneity on cancer cell viability. Physiologic breast tumor interstitial fluid velocities are predicted to be between 0-10 µm/s as compared to normal interstitial velocities between 0.1-1 µm/s 53 . The permeability sensitivity analysis accounted for velocity alterations within these ranges. Diffusion dominates convection when interstitial velocities are decreased to normal tissue values (0.1K and 0.01K), as is apparent from comparison to the diffusion alone concentration profile. Similarly, regional viability was comparable between the low permeability and diffusion only conditions. As permeability is increased, viability decreases due to rapid homogenization of concentration gradients. Therefore, decreasing the permeability would decrease the extravasation rate from the tumor bulk and subsequently increase the efficacy of overall drug penetration into the tumor, previously reported as a method to increase whole tumor therapeutic efficacy 106 . However, regional drug efficacy at the tumor border could be potentially attenuated with decreased permeability, as dispersed concentration profiles would develop in this region.

Implications of our models to therapeutic efficacy and clinical response
Together our data indicate that several TME-specific factors affect cancer cell viability within the tumor stoma transition region. Interestingly, we see that it is the cells nearest the modeled invading stroma that are more susceptible to therapy than those in the tumor bulk. There is a great deal of research focused on targeting or preventing development of invading cells within these stromal zones 55,107-109 as a means for preventing metastasis. Examination of tumor cells live in vivo indicates heterogeneity in the subsets of cells that invade through tissues 110 . Though it is intuitive to think that invaded cancer cells will not be accessed by therapies and thus will continue to invade and lead to systemic metastases, our data here may indicate that cells closer to the bulk may resist treatment better and potentially become invading cells posttherapeutic administration. Further, our results are limited by only examining one breast cancer cell line, MDA-MB-231 and one fibroblast type. Due to the ubiquitous presence of fibroblasts in carcinoma develop in multiple organs, our findings may translate in part, though broader interpretation requires further investigation to understand the interaction of intercellular tumor-stroma signaling and transport. This is particularly true for highly desmoplastic cancers, such as pancreatic carcinoma, which are known to be highly transport limited with extreme ratios of stroma to cancer compared to other solid tumors 111,112 .

Conclusion
Here we find that there is a distinctive interaction between the contributions of fibroblasts and transport to the activity of a common chemotherapeutic, doxorubicin, in a simulated breast cancer microenvironment. Our findings indicate that specific ratios of tumor cells to fibroblasts result in increased drug resistance and that these forces are dominant over transport limitations within our particularly in vitro and in silico models. These results indicate that it is important to study transport and cellular interactions simultaneously when examining cancer-related therapy as both factors play an important role in the in vivo response.

CHAPTER 5: Conclusion
The current work describes in vitro and in silico models utilized to probe the effects of the TME on cancer progression and metastasis. Together, these models characterize TME-specific effectors of cancer cell response and motivate further explorative research. Two distinct agent-based models were developed to analyze and predict the factors that correlate with cancer development and progression. The first model provides evidence for two major mechanisms of increased invasion in the brain, in the presence of interstitial flow. The second model predicts a specific cellular population that could be resistant to chemotherapy at the tumor border due to unique TME factors.
Additionally, in vitro 2D and 3D models elucidate the effects of fibroblast on cancer cell resistance to doxorubicin therapy. Novel in vitro layered transition hydrogel models were utilized to probe the region specific effects of fibroblast resistance. Due to therapeutic transport factors and the protective effect, regional variations in viability in response to flow-delivered therapy were observed, as predicted by agent-based models. Together these data indicate the importance of TME factors in altering cancer cell behavior, and predictive models further confirm these effects in a replicable manner.
Our in vitro layered model replicates the morphology of the tumor border transition and could be utilized for further probing of the population specific phenomenon observed through interactions with fibroblast. The specific model could be tuned to incorporate varied cell densities, cancer subtypes, and chemotherapy drugs. The layered hydrogel model can further be utilized to probe the effects of transport on therapeutic efficacy by modulating local ECM composition.
Most importantly, characterization of two separate in silico models with distinct cancer types indicates the transferability and utility for agent-based model use in multiple tumor microenvironments.
Components of these model can now be utilized to predict cancer response in more complex contexts.
While migration was not incorporated into the ABM of the breast cancer TME due to the minimal cancer cell velocity observed with live imaging 53 , research has shown that overall cancer cell migration is increased with fibroblast coculture 113 . This indicates a potential for probing the mechanism of increased invasion within breast cancer coculture conditions using the methodologies established for brain cancer migration. Further probing these migration characteristics in the presence of interstitial flow delivered chemotherapy will provide a broader understanding of invasive subpopulations within the breast cancer microenvironment.
Alternatively, little is known about the effects of chemotherapy on the invasive characteristics of glioma cells. Standard of care treatment for glioma patients includes adjuvant temozolomide treatment 114 , therefore, interstitial flow delivered chemotherapy is also of importance in the brain. Incorporating chemotherapy treatment into the established agent-based model of invasion mechanisms could provide insight into resistant subpopulations within the diffuse invasive edge, characteristic of brain cancers following surgical resection. Furthermore, it has been established that tumor associated cells in brain cancer alter the migration characteristics of glioma cells, and the presence of astrocytes within the brain TME limits patient survival in vivo 64 . Thus, it is possible that a similar cancer protective effect could be observed in brain tumors as in breast tumors, and the tumor border transition area could, likewise, be an important niche for protected cancer subpopulations.
Altogether, our research provides support for the pertinence of agent-based models for TME applications. Results obtained from in silico models demonstrate the usefulness of agent-based models for predicting the behavior of cancer cells in the TME, meditated by complex environmental factors. The versatility of this modeling strategy is confirmed through the predictive simulations that represent in vitro findings in two varied TME contexts. Thus, agent-based models are an optimal modeling platform for synthesizing complex interactions into tunable, robust representations of physical processes, and the techniques developed in this work could be utilized to predict similar effects in comparable settings.