Availability of essential species like oxygen is critical in shaping the dynamics of tumor growth. When the intracellular oxygen level falls below normal, it initiates major cascades in cellular dynamics leading to tumor cell survival. In a cellular block with cells growing away from the blood vessel, the scenario can be aggravated for the cells further inside the block. In this study, the dynamics of intracellular species inside a colony of tumor cells are investigated by varying the cell-block thickness and cell types in a microfluidic cell culture device. The oxygen transport across the cell block is modeled through diffusion, while ascorbate (AS) transport from the extracellular medium is addressed by a concentration-dependent uptake model. The extracellular and intracellular descriptions were coupled through the consumption and traffic of species from the microchannel to the cell block. Our model shows that the onset of hypoxia is possible in HeLa cell within minutes depending on the cell location, although the nutrient supply inside the channel is maintained in normoxic levels. This eventually leads to total oxygen deprivation inside the cell block in the extreme case, representing the development of a necrotic core that maintains a dynamic balance with growing cells and scarce supply. The numerical model reveals that species concentration and hypoxic response are different for HeLa and HelaS3 cells. Results also indicate that the long-term hypoxic response from a microfluidic cellular block stays within 5% of the values of a tissue with the basal layer. The hybrid model can be very useful in designing microfluidic experiments to satisfactorily predict the tissue-level response in cancer research.

Introduction

Cells in every living organism depend crucially on the supply of oxygen for survival. In a typical human tissue, consisting of thousands of contiguous cells, an intricate network of blood vessels continuously supplies oxygen and other nutrients. However, tissues under cancerous growth often have chaotic vessel architecture, altered blood rheology, and defective vessel autoregulation. This leads to transient interruptions in tissue blood flow and consequently regions of tissue with very low levels of oxygen, a condition often termed hypoxia. Surprisingly, tumor cells exhibit significant resistance to cellular apoptosis that should follow hypoxia. This response has been largely attributed to the transcription factor hypoxia-inducible factor (HIF) in mammalian response to oxygen levels [1]. This factor is principally a heterodimer of bHLH-PAS proteins which has been identified as a key element in one of the most stable evolutionarily conserved pathways [2].

The first characterized member of the HIF family, HIF1, are expressed both in α and β isoforms in mammalian cells. Their activation in hypoxic conditions has been associated with a large number of transcriptional activity that influences the formation of blood vessels (angiogenesis), nutrient transport, switch in energy metabolism, cell proliferation, and migration [3]. In normal and well-oxygenated conditions (normoxic), HIF is hydroxylated by intracellular prolyl and asparaginyl hydroxylases and eventually degraded. A number of intracellular pathways investigating their dynamics have been proposed. In Vitro studies using different cell-lines and in vivo mouse xenograft models have been used extensively for the study of HIF pathways [4,5]. A number of cell culture experiments have reported the impairment of tumor growth under the ablation of HIF [6]. In recent years, microfluidic cell culture devices have been used to study cancer cell growth, cell cycle, apoptosis, and other mechanisms [7]. These devices can be used for anticancer drug testing and to study the cellular response over desired time span under easily reproducible conditions, which are not possible in traditional in vitro experiments [8]. Although in vitro experiments have characterized the primary dynamics of the HIF–prolyl hydroxylase domain (PHD) pathways, detailed mapping of intracellular interactions between different cellular compartments (e.g., mitochondria, nucleus, etc.) and their impact on tumor progression have not been explored. More importantly, the numerous extracellular influences that are common in vivo are still poorly understood.

System level biophysical models have opened up a new way of exploring the mechanistic behavior of these intracellular interactions. The mathematical models allow for the testing of biologically plausible hypotheses and suggest explanations for conflicting observations from experiments. Over the past decade, a number of models have been developed to investigate HIF interactions. The theoretical model developed by Kohn et al. [9] hypothesized about a switchlike behavior in HIF activity under decreasing oxygen gradient which has been later verified by Qutub and Popel [10]. Nguyen et al. [11] included the effect of both prolyl and asparaginyl hydroxylases in the study of intracellular interactions with compartmentalized reactions. However, the early models are limited in detailed intracellular networks and rarely considered the influence of extracellular effects. More recently, we have proposed a model considering the coupling between extracellular nutrient supply in a microfluidic setting [12]. Still, all these models focus on the dynamics of a single cell and its environment. For a growing tumor, tissue-level effects need to be included in the mathematical model for a complete description of the hypoxic environment.

In this study, we have developed a tissue-level model with contiguous cells inside a microfluidic environment. The microfluidic platform has immense potential in exploring tissue-level cellular activity in cancer research. Transport and evolution of different species in both extracellular and intracellular spaces are addressed with an appropriate mathematical description along with interaction and exchange between these two domains. Change in tissue dimensions that affect the nutrient supply has been studied systematically. Also under varying levels of hypoxia, we investigated the transcriptional response in the tissue. Parameters involving intracellular interactions, such as hydroxylation, hypoxia feedback and response, have been quantified. The detailed mathematical modeling of the microfluidic system can be useful in designing efficient systems and identify target key species.

The rest of the paper is arranged as follows: First, we present the mathematical formulations for the extracellular transport inside a microfluidic cell culture device. Then, the relations for consumption and uptake by the tissue regions are presented. Next, the intracellular reaction network is discussed. We have compartmentalized the cytosolic and nuclear interactions along with shuttling of species between the compartments. In our model, we have kept the supply of species at the normoxic level in the flow media, but cell-block thickness is increased to see hypoxic effect further inside the cellular block representing a growing tissue. We also studied the variation in hypoxic response with changing cell types. Finally, the effect of basal layer is probed in hypoxic response.

Mathematical Model

The physical system in our model consists of a microfluidic cell culture environment where the cell growth and nutrient supply are controlled externally. For the study of hypoxic response, the major nutrient of concern is oxygen dissolved in the extracellular media. Additionally, other key cofactors of the cellular hypoxia pathway, such as ascorbate (AS), iron, and glucose, can also be supplied in the extracellular environment. The control of extracellular nutrient concentration is critical in maintaining a viable environment for the cells to thrive. The spatiotemporal dynamics of the extracellular and intracellular domains and their crosstalk are presented in this section.

Extracellular Region.

The general transport of nutrients in the microfluidic chip is represented by the species mass conservation equation [13]. For oxygen transport, the governing equation can be given as
(1)

where CO is the instantaneous spatial concentration of oxygen, u¯ is the medium fluid velocity, DO is the diffusion coefficient for oxygen in the medium, and RO is the generalized source/sink contribution involving reactions. The concentration at the inlet of the flow channel maintained at a constant level, while the outlet was set as fully developed with the channel side walls as impermeable. A constant medium flow velocity of 100 μm/s was used at the center of the channel following the low-shear mass transport experiments [14].

Similar to oxygen, the transport of ascorbate through the microchannel is modeled as
(2)

where CAS is the instantaneous spatial concentration of ascorbate, DAs is the diffusion coefficient for ascorbate in the medium, and RAS is the sink term for ascorbate uptake by the cells. The ascorbate inlet concentration in the channel was held constant as specified by a Dirichlet condition. All the relevant initial and boundary conditions are listed in Table 1.

Table 1

Initial and boundary conditions for transport equations in the microfluidic extracellular domain

Initial conditionsaInletbSide wallscOutletd
OxygenCO(x,y)=207μMC0=CgivenC0y=0C0x=0
AscorbateCAS(x,y)=10μMCAS=CgivenCASy=0CASx=0
Initial conditionsaInletbSide wallscOutletd
OxygenCO(x,y)=207μMC0=CgivenC0y=0C0x=0
AscorbateCAS(x,y)=10μMCAS=CgivenCASy=0CASx=0
a

Unless otherwise mentioned.

b

Dirichlet conditions chosen according to species supply concentration. Default concentrations, reported in supplementary Table S3 which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection, are used unless otherwise mentioned.

c

The side walls are assumed impermeable to the nutrients.

d

The concentrations are lengthwise invariant, and the flow (subregion a) is fully developed at the outlet. The cell-block part of the outlet (subregion b) is maintained impermeable to the species.

Communication Between Intracellular and Extracellular Regions

Transmembrane Oxygen Transport.

Cells continuously consume oxygen in order to maintain the metabolic reactions. Inside the microfluidic chip, this consumption naturally results in oxygen deficient regions where cells have been cultured. Oxygen is among a small number of molecules that can directly diffuse through cell membrane whose diffusion parameters and consumption rates have been extensively studied in the literature [15]. In our model system, oxygen from the flow channel can diffuse through the cell membrane/basal membrane of the tissue block. Since most of the cellular mass is aqueous, the diffusivity of oxygen in water, blood, and tissue is in the same order. However, often basal membrane impedes free transport resulting in a lower diffusivity. All the relevant transport properties used in this study have been presented in Table 2.

Table 2

Transport parameters

ParameterValueReference
Oxygen diffusivity in aqueous media2.66×109 m2/s[16]
Oxygen diffusivity in basal layer(2.0±0.5)×1010 m2/s[17]
Oxygen diffusivity in tissue1.75×109 m2/s[18]
Ascorbate diffusivity in aqueous media6.63×1010 m2/s[19]
Ascorbate diffusivity in tissue2.13×1011 m2/s[20]
HeLa OCR26.86×1018 moles/(cell s)[21]
HeLa mitochondria knockout OCR12.5×1018 moles/(cell s)
HeLaS3 OCR8.42×1018 moles/(cell s)
HeLaS3 mitochondria knockout OCR5.64×1018 moles/(cell s)
ParameterValueReference
Oxygen diffusivity in aqueous media2.66×109 m2/s[16]
Oxygen diffusivity in basal layer(2.0±0.5)×1010 m2/s[17]
Oxygen diffusivity in tissue1.75×109 m2/s[18]
Ascorbate diffusivity in aqueous media6.63×1010 m2/s[19]
Ascorbate diffusivity in tissue2.13×1011 m2/s[20]
HeLa OCR26.86×1018 moles/(cell s)[21]
HeLa mitochondria knockout OCR12.5×1018 moles/(cell s)
HeLaS3 OCR8.42×1018 moles/(cell s)
HeLaS3 mitochondria knockout OCR5.64×1018 moles/(cell s)

The oxygen consumption rate (OCR) depends on the type of cells and activities. For instance, in healthy cells, mitochondria consume most of the available oxygen in meeting the energy requirements [21]. However, hypoxic tumor cells have been known to reduce this demand of oxygen by selectively shutting off the mitochondrial activity [22]. On the other hand, experimental studies have related selective damage to the mitochondrial genes (or any member of the mitochondrial electron transfer chain) directly to reduction in oxygen consumption, tumor growth, and increased glycolysis [23]. In our model, we combined these two well established lines of research to use the OCR of mitochondrial knockout cells as a baseline value to represent mitochondrial activity suppression in tumor cells. It should be remembered, however, that the oxygen consumption rate is constantly varying in cells in vivo, which is very difficult to quantify continuously in experiments.

In our model, we represent the consumption as localized sink terms for the tissue regions. In in vitro cell culture experiment, oxygen levels are commonly expressed as the volume fraction of oxygen in the air supplied to the cell medium. Thus, for an oxygen solubility of 1.30 μmol/l/mm Hg in water (at 37 °C and standard pressure), 21% ambient oxygen corresponds to 207 μM oxygen in the solution [5]. We have considered concentrations less than 5% to be hypoxic and changed oxygen consumption rates to corresponding mitochondria knock-out cells. The consumption rates used in the present model for different cellular states are also reported in Table 2.

Ascorbate Uptake.

Ascorbates are consumed by human cells in different oxidized states, since we do not produce it ourselves. These different states include ascorbic acid, ascorbate, ascorbate free radical, and dehydroascorbic acid. The primary active transport of ascorbate is due to the sodium-dependent vitamin C transporter proteins, but the glucose transporters can also play a role in ascorbate uptake from plasma [24]. Availability, selectivity, and concentration of these transporters vary with cell types. Although the uptake of intestinal ascorbate is very tightly regulated [25], recent pharmacokinetic data suggest that intravenous administration of ascorbate can result in much higher ascorbate level in plasma [26]. We have previously proposed and validated a simplified ascorbate uptake model against cell-well based experimental data [12]. The model captures the ascorbate accumulation behavior over a period of 5 h. The extracellular concentration of ascorbate is taken as substrate concentration in the Michaelis–Menten enzyme-substrate kinetic model as
(3)

Here, A˙Sin and ASext denote the rate of ascorbate uptake and the extracellular ascorbate concentration, respectively, and ρcell is the volumetric cell density. The Michaelis–Menten parameters were estimated based on experimental data and are reported in supplementary Table S1 which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection. It is important to note that the parameters will most certainly be different for different cell types, and therefore with appropriate data, they can be adjusted for different cell lines.

Intracellular Region.

The intracellular region in our model consists of arrays of contiguous cells. Each cell is represented by a 25 μm × 25 μm rectangular zone. The complex intracellular biochemistry of hypoxic response involves both intracellular and extracellular interactions. The shift in cellular equilibrium due to hypoxia manifests transcriptional activity in the nucleus. We have therefore compartmentalized the intracellular interactions to cytosolic and nuclear reactions, with different species moving back and forth between these two regions. In this model, we have primarily focused on the HIF hydroxylation pathways, as shown in Fig. 1. All the species and relevant abbreviations are listed in Table 3.

Fig. 1
Schematic of the normoxic and hypoxic reaction pathways. Hypoxia-inducible factors (HIF1α) gets hydroxylated by PHDs and the FIH in parallel pathways. PHD gets activated by a number of other cofactors such as iron (Fe), DG, and AS. Hydroxylated HIF later gets degraded by ubiquitination reaction after being tagged by vHL proteins. When there is a lack of oxygen supply, both hydroxylation pathways lose their effectiveness. HIF, in this case, accumulates in cytosol and shuttles inside the nucleus where it also faces hydroxylation pathways but it can also form a dimer (HIFd) with HIF1β. This dimer can activate the hypoxia response elements and initiate a large number of transcriptional activities. Feedback mechanisms such as elevated PHD production are also activated at this stage. Uptake and consumption of extracellular nutrients (O2 and AS) have also been considered in the model.
Fig. 1
Schematic of the normoxic and hypoxic reaction pathways. Hypoxia-inducible factors (HIF1α) gets hydroxylated by PHDs and the FIH in parallel pathways. PHD gets activated by a number of other cofactors such as iron (Fe), DG, and AS. Hydroxylated HIF later gets degraded by ubiquitination reaction after being tagged by vHL proteins. When there is a lack of oxygen supply, both hydroxylation pathways lose their effectiveness. HIF, in this case, accumulates in cytosol and shuttles inside the nucleus where it also faces hydroxylation pathways but it can also form a dimer (HIFd) with HIF1β. This dimer can activate the hypoxia response elements and initiate a large number of transcriptional activities. Feedback mechanisms such as elevated PHD production are also activated at this stage. Uptake and consumption of extracellular nutrients (O2 and AS) have also been considered in the model.
Close modal
Table 3

List of species and relevant abbreviations

NotationDescription
OOxygen
ASAscorbate
PHDProlyl hydroxylases
Fe2+/Fe3+Iron (ferrous/ferric)
PFBinding of iron and prolyl hydroxylases
DG2-oxoglutarate
PFDBinding of PF and 2-oxoglutarate
PFDOBinding of PFD and oxygen
PFDOABinding of PFDO and ascorbate
HIF1αHypoxia-inducible factor 1 α
HIF1αPhProlyl hydroxylated HIF
HIF1αAhAsparaginyl hydroxylated HIF
HIF1αPAhProlyl hydroxylated HIF1αAh
FIHFactors inhibiting HIF
HIF1βHypoxia-inducible factor 1 β
HIFdDimer of HIF1α and HIF1β
HREHypoxia response elements
HIFdHREBinding of HIFd and HRE
vHLvon Hippel–Lindau
H2O2Hydrogen peroxide
Subscript “nQuantity inside nucleus
NotationDescription
OOxygen
ASAscorbate
PHDProlyl hydroxylases
Fe2+/Fe3+Iron (ferrous/ferric)
PFBinding of iron and prolyl hydroxylases
DG2-oxoglutarate
PFDBinding of PF and 2-oxoglutarate
PFDOBinding of PFD and oxygen
PFDOABinding of PFDO and ascorbate
HIF1αHypoxia-inducible factor 1 α
HIF1αPhProlyl hydroxylated HIF
HIF1αAhAsparaginyl hydroxylated HIF
HIF1αPAhProlyl hydroxylated HIF1αAh
FIHFactors inhibiting HIF
HIF1βHypoxia-inducible factor 1 β
HIFdDimer of HIF1α and HIF1β
HREHypoxia response elements
HIFdHREBinding of HIFd and HRE
vHLvon Hippel–Lindau
H2O2Hydrogen peroxide
Subscript “nQuantity inside nucleus

Normoxic levels of oxygen allow the PHDs to hydroxylate HIF at their N-terminal transactivation domains. However, a number of experiments have suggested that this reaction is facilitated by a number of intracellular cofactors like iron (Fe), 2-oxoglutarate (DG), AS, etc. [4]. The sequential binding of collagen prolyl hydroxylase during hydroxylation has been identified to proceed in the order of iron → 2-oxoglutarate → oxygen → peptide [27]. The same sequence of binding has been suggested for HIF hydroxylation by PHD [28]. It should be mentioned that different PHD isoforms have specific nuclear or cytosolic localizations, but we have treated them as one component for simplicity.

The intracellular ascorbate performs a number of interactions. Other than being a cofactor for PHD (Eq. (S2.4)) in the PHD–HIF hydroxylation, it also reduces ferric ions (Fe3+) that were oxidized by cellular peroxides (H2O2). The oxidation of ferrous ions also produces reactive oxygen species, which contributes to a number cellular activity switches like apoptosis.

Activity of the PHD proteins is highly dependent on the oxygen concentration. But under well-oxygenated conditions (aka normoxia) the primarily unidirectional prolyl hydroxylation (HIF1αPh) ensures a steady drop in intracellular HIF concentrations. Factors inhibiting HIF represent another major HIF hydroxylation component known as asparaginyl hydroxylation (HIF1αAh). Although this pathway remains active at a lower concentration of oxygen than the PHDs, asparaginyl hydroxylation is also known to be reversible. Activated PHD isoforms also participate in the prolyl hydroxylation of the asparaginyl hydroxylases. Hydroxylation of HIF results in a binding site for the von Hippel–Lindau (vHL) tumor suppressor protein. This vHL protein tags the HIF for ubiquitylation and eventual proteasomal degradation. Hypoxia-inducible factor can also translocate into the nucleus, but the nuclear PHD hydroxylates the translocated HIF and subsequently degraded by vHL similar to cytosolic reactions.

In low-oxygen (hypoxic) environment, HIF escapes hydroxylation in both compartments and accumulates over time in the nucleus. The HIFβ isoform which is primarily located inside the nucleus binds with nuclear HIF1α and produces the HIF dimer (HIFd). This dimer can bind to hypoxia response elements (HREs) in hypoxia responsive genes and contribute to a host of transcriptional activity. A negative feedback of HIF-dependent activation involves the production of PHDs which can contribute to the reduction of HIF levels. However, cytosolic translation and degradation of HIF also play their roles in influencing this equilibrium.

A list of all the reactions described above is given in supplementary Table S2 which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection. The corresponding mathematical representation of these chemical reactions involves formulation of a system of ordinary differential equations. Mass-law kinetics and Michaelis–Menten forms have been used for these formulations, and the corresponding ordinary differential equations are reported in supplementary Table S3 which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection. For example, the temporal evolution of intracellular ascorbate (Eq. (S3.8)) involves binding with PFDO complex resulting in a drop proportional to their instantaneous concentrations (first term), a regain proportional to the amount available through the breakdown of the PFDOA complex (second term), loss in ferric to ferrous reduction reaction (third term), and overall uptake from extracellular region constituting a rise (fourth term).

Numerical Methods

We used a hybrid mathematical description in order to represent the intracellular dynamics as well as the microfluidic operations. As presented in Fig. 2, the channel width is divided into two sections: the width of the unobstructed flow area is “ a,” and the width of the cell block is “ b.” The ratio of a to b is defined as a nondimensional parameter δ. Cell boundaries were considered to be aligned with domain grids, and the length and width of each cell are 25 μm. To see the oxygen penetration across the tissue, we chose to take a tissue section at least 60 cells wide (∼1500 μm) for the standard case. Coupling that width with equally wide flow channel resulted in a 3 mm wide channel. The width of the flow channel or tissue block can be further changed by adding more computational domain. From a biological perspective, this feature of the model could be used to mimic blood vessel of different widths found in vivo. The length of the channel could be designed differently per the needs of the experiment. For our chosen flow velocity of 100 μm/s, the 1.5 cm flow length provided a satisfactory fully developed condition at the outlet.

Fig. 2
Schematic representation of the computational domain inside a microfluidic channel. Total channel width consists of the unobstructed flow region (a) and cell-block region (b). The nondimensional parameter δ represents the ratio of a to b, which varies with the thickness of cell block. The basal layer, encountered in in vivo experiments, poses an additional impediment for free entry of different species into tissue. Extracellular species like oxygen and AS are transported in the fluid channel. Small species like oxygen can diffuse through the cell block. However, AS transport inside the block requires intervention of active transporter proteins. Depending on theavailability of oxygen and other cofactors, each cell undergoes through a cascade reactions presented in Fig. 1 and supplementary Table S2 which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection. In this study, the length of the channel is 1.5 cm and a+b=const=3000μm.
Fig. 2
Schematic representation of the computational domain inside a microfluidic channel. Total channel width consists of the unobstructed flow region (a) and cell-block region (b). The nondimensional parameter δ represents the ratio of a to b, which varies with the thickness of cell block. The basal layer, encountered in in vivo experiments, poses an additional impediment for free entry of different species into tissue. Extracellular species like oxygen and AS are transported in the fluid channel. Small species like oxygen can diffuse through the cell block. However, AS transport inside the block requires intervention of active transporter proteins. Depending on theavailability of oxygen and other cofactors, each cell undergoes through a cascade reactions presented in Fig. 1 and supplementary Table S2 which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection. In this study, the length of the channel is 1.5 cm and a+b=const=3000μm.
Close modal

For each cell in the block, solving a set of ordinary differential equations (see supplementary Table S3 which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection) resulted in the instantaneous intracellular concentrations. The system of ODE is solved using a fourth-order Runge–Kutta algorithm, for which the parameter values such as rate constants and initial concentrations are reported in supplementary Tables S1 and S4, respectively, which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection. The intracellular system is compartmentalized between the cytosol and the nucleus with appropriate shuttling between compartments. However, each compartment in a cell is assumed to be well-mixed.

Due to its small size, oxygen can diffuse through the cell membrane [17] and the cytosol [18]. In case of oxygen, the whole domain (subregions a and b in Fig. 2) is considered in Eq. (1) as a continuous region. However, the diffusivity is different in each subregion and at their interface. Moreover, there is no advection in the tissue-block. Although the uptake of oxygen can be represented by an experimentally tuned function to represent the sink term (Eq. (5) in Ref. [12]), in the present study we chose to use the experimentally reported oxygen consumption rates from Ref. [21] as the sink term for oxygen in the cell-block region after multiplying it with the volumetric cell density (RO=ρcell×OCR). The consumption rates as well as the diffusivity of oxygen in different subregions are reported in Table 2. In case of ascorbate however, the primary source term was concentrated in cells (control volumes) adjacent to the blood vessel. It arose due to the active transport of ascorbate by transporter protein given by Eq. (3).

The governing partial differential equations (Eqs. (1) and (2)) are solved using the finite volume method [12,13]. These transport equations were discretized spatially following exponential (exact) scheme [29], and for the temporal evolution, a first-order implicit formulation was used. Discretized system of algebraic equations was solved iteratively in an in-house code written in FORTRAN95 with tridiagonal matrix algorithm. Since the two-dimensional matrix is solved line-by-line tridiagonal matrix algorithm approach, the solution must be iterated until a satisfactory convergence is achieved. A convergence tolerance of 10−5 was used with a time step size of 0.05 s. The numerical calculations were carried out on a Linux platform using a 2.4 GHz Intel Xeon CPU. A flowchart of the hybrid algorithm is given in Fig. 3.

Fig. 3
Flowchart for the hybrid algorithm. Each cell in the domain must be evaluated individually at each time for temporal evolution of species, which is also connected to the spatial distribution of the species presented in the extracellular domain that freely diffuses through the cellular block or, transported in actively. The spatial solution must go through iterations in order to converge at every time step.
Fig. 3
Flowchart for the hybrid algorithm. Each cell in the domain must be evaluated individually at each time for temporal evolution of species, which is also connected to the spatial distribution of the species presented in the extracellular domain that freely diffuses through the cellular block or, transported in actively. The spatial solution must go through iterations in order to converge at every time step.
Close modal

Results and Discussion

The response of cancer cells at the tissue level is studied by integrating a space and time-dependent mathematical representation for oxygen and ascorbate inside the whole computational domain along with a system of ordinary differential equations to take care of the intracellular dynamics for each cell. The thickness of tissue inside the microchannel is controlled by changing the number of cells in y -direction, while the width of the microfluidic channel (computational domain) was kept constant. We specifically studied the transport of oxygen inside the cell block by considering the appropriate oxygen consumption rate during hypoxia and normoxia by maintaining a constant concentration of oxygen (207 μM) at the microchannel entry.

Onset of Hypoxia.

While the dissolved oxygen gets transported through the microfluidic channel by both advection and diffusion, it can only move by diffusion once inside the cell block. Owing to the smaller size of the oxygen molecule, the diffused oxygen can get inside the cell without any active transport mechanism. The intracellular oxygen is continuously consumed primarily in mitochondria for energy production and for a number of other functions [21]. In a growing tumor tissue, this generally creates an oxygen lean core which might become necrotic. The presented system (Fig. 2) explores this idea by keeping the extracellular media oxygen rich, while allowing the cells to consume oxygen continuously. The oxygen concentration evolution, as shown in Fig. 4, is quick inside the microfluidic device. Within minutes, two clear zones are formed inside the device (Fig. 4(a)). While the flow channel is maintained normoxic, the cell block rather quickly becomes hypoxic and reaches a steady-state within the first 5 min. Figures 4(b) and 4(c) show two representative cases with 60 and 30 HeLaS3 cells in the y -direction. Our results, in the dynamic equilibrium, were compared with model results from Ref. [31], as a reference for both cases. They have used a hybrid stochastic model for growing the tissue away from the blood vessel and provided normalized results for the first 30 cells. Further validations for the model and detailed discussion on its suitability for the HIF pathway can be found in Ref. [12].

Fig. 4
Spatiotemporal oxygen concentration distribution with HeLaS3 cells. (a) The contour of oxygen levels in the whole computational domain after 30 min showing a jump in concentration across the basal layer of the cell block. Concentration distribution inside the block approaches dynamic equilibrium within minutes. (b) Oxygen levels across the domain at different times for δ = 1.0 and (c) δ = 3.0, where the equilibrium level inside the cell block is compared with results from Ref. [30] for model validation. At the inlet boundary, supply for oxygen was maintained at 207 μM, while the initial intracellular concentrations are reported in supplementary Table S4 which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection. The side walls of the channel are taken as impervious and the outlet is fully developed.
Fig. 4
Spatiotemporal oxygen concentration distribution with HeLaS3 cells. (a) The contour of oxygen levels in the whole computational domain after 30 min showing a jump in concentration across the basal layer of the cell block. Concentration distribution inside the block approaches dynamic equilibrium within minutes. (b) Oxygen levels across the domain at different times for δ = 1.0 and (c) δ = 3.0, where the equilibrium level inside the cell block is compared with results from Ref. [30] for model validation. At the inlet boundary, supply for oxygen was maintained at 207 μM, while the initial intracellular concentrations are reported in supplementary Table S4 which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection. The side walls of the channel are taken as impervious and the outlet is fully developed.
Close modal

Effect of Tissue Thickness on Oxygen Level.

Next, we study the oxygen concentration distribution in HeLa cell by changing the thickness of the tumor tissue. The oxygen concentration evolution in the domain is shown in Fig. 5 for four such cases by increasing the number of cells in the y-direction from 10 to 90. Figures 5(a) and 5(b) are very representative of a normoxic scenario in the tissue, where we have considered 10 and 30 cells, respectively, in the y-directions. In these two cases, the normoxic supply of oxygen is sufficient enough to maintain the oxygen balance across the whole tissue even with constant consumption. When the tissue width is increased further to 60 and 90 cells (Figs. 5(c) and 5(d)), the oxygen demand far exceeds the supply through diffusion from the blood vessel. As a result, within minutes, the distal cells become highly hypoxic. Our predictions with 90 cells resulted in severely low oxygen after 50 cells. However, that does not restrict the spread of necrotic core until the 50 cell limit. In reality, chronic hypoxia can result in a necrotic region starting even after 10–15 cells from the blood vessel depending on the stage of the tumor [32].

Fig. 5
Effect of cell-block thickness on oxygen concentration. Results are obtained at the midlength of the channel. Four δ values (a) 11.0, (b) 3.0, (c) 1.0, and (d) 0.33 were considered which correspond to 10, 30, 60, and 90 HeLa cells in the y-direction, respectively. For δ value as low as 1.0, although the cells are in hypoxia, they maintain a steady dynamic equilibrium. However, for δ less than 1.0, we see almost negligible oxygen further inside the block. For all cases, entry of extracellular medium was maintained at 207 μM for oxygen, and the initial intracellular concentrations are given in supplementary Table S3 which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection.
Fig. 5
Effect of cell-block thickness on oxygen concentration. Results are obtained at the midlength of the channel. Four δ values (a) 11.0, (b) 3.0, (c) 1.0, and (d) 0.33 were considered which correspond to 10, 30, 60, and 90 HeLa cells in the y-direction, respectively. For δ value as low as 1.0, although the cells are in hypoxia, they maintain a steady dynamic equilibrium. However, for δ less than 1.0, we see almost negligible oxygen further inside the block. For all cases, entry of extracellular medium was maintained at 207 μM for oxygen, and the initial intracellular concentrations are given in supplementary Table S3 which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection.
Close modal

Hypoxic Response Inside the Cell Block.

Figure 6 shows the hypoxic response from the intracellular reaction network across the block for the four cases presented in Fig. 5. It should be noted that this response includes a reduction in oxygen consumption due to a glycolytic switch in energy metabolism, a positive feedback producing higher transcriptional activity especially in the chronic case. Moreover, negative feedbacks resulting in the synthesis of PHD to hydroxylate free HIF are also considered in the model. Since the transcriptional activity results in activation of hundreds of genes in the nucleus, we present the generalized mRNA response from the model as a transcriptional activity. As expected, the highest transcriptional activity was observed in the case of 90 cell block and occurred at the furthest cell from the flow channel. This value has been used to normalize hypoxia response for all other cases. The accumulation and chronic activity are apparent in all four cases when compared to Fig. 5, where the oxygen concentration over the whole block reaches local equilibrium, but the HIF activity is accumulative and keeps increasing with time.

Fig. 6
Effect of cell-block thickness on the transcriptional activity. The variation of tissue thickness and initial conditions are kept exactly same as those in Fig. 5. Although oxygen concentration reaches equilibrium fairly quickly, the extent of hypoxic activity increases steadily once the cells are kept in hypoxic conditions. The highest level of activity was observed for the cells furthest away from the channel. The level found at 1 h for δ = 0.33 was used to normalize transcriptional activity for all four cases considered.
Fig. 6
Effect of cell-block thickness on the transcriptional activity. The variation of tissue thickness and initial conditions are kept exactly same as those in Fig. 5. Although oxygen concentration reaches equilibrium fairly quickly, the extent of hypoxic activity increases steadily once the cells are kept in hypoxic conditions. The highest level of activity was observed for the cells furthest away from the channel. The level found at 1 h for δ = 0.33 was used to normalize transcriptional activity for all four cases considered.
Close modal

Hypoxic Response in Different Types of Cancer Cells.

The OCR is also a very important parameter in hypoxia study. This parameter is responsible for dictating the equilibrium level reached in the tissue and will certainly result in variations depending on cell types. The effect of oxygen consumption rate is studied by considering two different types of cancer cells: HeLa and HeLaS3 cells. The OCR values in healthy and mitochondria knockout variants of both cells were reported in Ref. [21]. Owing to the difference in OCR, changes in oxygen concentration over the whole block are significantly high in early stages. However, they reach a steady difference especially inside the block as the cell blocks reach oxygen equilibrium (Fig. 7(a)). On the contrary, the difference in hypoxic response (Fig. 7(b)) remains higher in the front portion of the block over the whole time and is found lower further inside due to the similarity in hypoxic activity at chronic stages.

Fig. 7
Effect of cell type variation through a change in the OCR. Different OCR levels for HeLa and HeLaS3 were used from experimental measurements [21]. In all cases, absolute values of the difference were normalized by corresponding levels observed for HeLa cells. (a) Numerically predicted difference in oxygen concentration at different time and location. Different OCR values for HeLa and heLaS3 result in a significant variation in oxygen levels inside the block initially. However, it reaches steady levels quickly and maintains a fixed level after a distance of about 15 cell diameters. (b) Difference in transcriptional activities at different locations (along the y-direction) inside the tissue calculated at midlength. Initial difference is high until a dynamic equilibria is reached over the domain.
Fig. 7
Effect of cell type variation through a change in the OCR. Different OCR levels for HeLa and HeLaS3 were used from experimental measurements [21]. In all cases, absolute values of the difference were normalized by corresponding levels observed for HeLa cells. (a) Numerically predicted difference in oxygen concentration at different time and location. Different OCR values for HeLa and heLaS3 result in a significant variation in oxygen levels inside the block initially. However, it reaches steady levels quickly and maintains a fixed level after a distance of about 15 cell diameters. (b) Difference in transcriptional activities at different locations (along the y-direction) inside the tissue calculated at midlength. Initial difference is high until a dynamic equilibria is reached over the domain.
Close modal

Effect of Basal Layer in Microfluidic Setup.

Although a layer of epithelial cells forms the blood vessel boundary, it is not possible to exactly replicate the in vivo tissue structure in microfluidic cell culture environment. In most microfluidic studies, cellular blocks are formed by growing cells inside a stagnant or flowing media. Therefore, we studied the difference in oxygen levels and corresponding changes in hypoxic response when we consider the cellular block with and without the basal layer (Figs. 8(a) and 8(b)). It is important to note that from a mass transport point of view, the basal layer poses an impediment, since it generally has one order lower mass diffusivity for oxygen when compared to the cell block. As a result, we see a sharp jump in the difference in oxygen levels just at the beginning of the block. Over time, this difference propagates to the whole block and reaches a steady level of 20% in the distal end. More importantly, the difference in hypoxic response (Fig. 8(b)), although is very high initially, reaches a paltry level of 5% in less than 40 min over most of the block. This result signifies the potential of studying the transcriptional activity inside a block of cells in microfluidic cell culture device. Our numerical prediction indicates that the hypoxic cellular responses in the microfluidic chamber will be very close to those in vivo.

Fig. 8
Difference in (a) oxygen levels and (b) transcriptional activity with and without the artificial basal layer with HeLa cells at midlength of the channel. In both cases, absolute values of the difference were normalized by corresponding levels observed for the cell block with the base layer. The sharp difference in oxygen at the channel-block interface is due to the absence of the basal layer which has an insulator like effect on gradient driven mass transport. Initially, the cells deeper inside the block have similar levels of oxygen; hence, the percentage difference is almost zero. But within minutes the diffusion front reaches cells at the distal end in the case without the basal layer resulting in the rise in percentage difference in oxygen concentration. Again, the change becomes same after 10 min and constant after ∼15 cells from the start of the block. Due to the comparatively large difference in concentration just at the interface for the two cases, transcriptional activity is most prominent near this area. However, as we move inside the block, it reaches a rather small (∼5%) and consistent level in most of the cell block. All initial concentrations are presented in supplementary Table S4 which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection.
Fig. 8
Difference in (a) oxygen levels and (b) transcriptional activity with and without the artificial basal layer with HeLa cells at midlength of the channel. In both cases, absolute values of the difference were normalized by corresponding levels observed for the cell block with the base layer. The sharp difference in oxygen at the channel-block interface is due to the absence of the basal layer which has an insulator like effect on gradient driven mass transport. Initially, the cells deeper inside the block have similar levels of oxygen; hence, the percentage difference is almost zero. But within minutes the diffusion front reaches cells at the distal end in the case without the basal layer resulting in the rise in percentage difference in oxygen concentration. Again, the change becomes same after 10 min and constant after ∼15 cells from the start of the block. Due to the comparatively large difference in concentration just at the interface for the two cases, transcriptional activity is most prominent near this area. However, as we move inside the block, it reaches a rather small (∼5%) and consistent level in most of the cell block. All initial concentrations are presented in supplementary Table S4 which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection.
Close modal

Conclusions

A highly typical scenario that arises in growing tumor cells and in ductal carcinoma has been modeled in the current study. While tumor cells tend to proliferate indefinitely, their growth is cut short due to the lack of nutrients. The hypoxic environment forces the cells to change their collective behavior in order to survive in this kind of environments through a shift in metabolism and activation of a large number of transcriptional activities. The hypoxia-inducible factors serve as the prime proponent in activating the mRNA signaling cascades under both acute and chronic hypoxia. Although the PHD–HIF pathways have been widely studied and characterized in vitro, the collective cellular behavior in tissues is still poorly understood. Our model of the cellular block inside a microfluidic cell culture environment shows a number of interesting features that are helpful in designing new experiments for future studies.

We used general species transport description for oxygen and other nutrients inside the microchannel that also houses the cell block mimicking tissue in vivo next to a blood vessel. Oxygen can freely diffuse inside the cell block which is consumed by cell continuously, albeit at different rate. However, the consumption rate and intracellular evolution are dependent on the oxygen level available. The intracellular reaction network is used to predict cellular behavior in the temporal domain. The equilibrium concentration of oxygen across the block was validated in this study for different tissue thickness.

Variation in equilibrium oxygen levels and consequent changes in hypoxic behavior were observed for different types of cancer cells. We included important autocrine feedback mechanism which resulted in sustained transcriptional activity. Although we have presented a short timespan (∼1 h) in this study, longer periods can be studied for understanding chronic responses in future.

Our study suggests that the difference in cellular hypoxic response is rather small from the cellular block with and without a separate basal membrane, especially at cells far away from the basal membrane. However, the cellular blocks will most certainly have paracrine response influenced by the transcriptional activity of neighboring cells. This kind of effect is not yet considered in this model which is an important avenue to explore in future. Moreover, cell morphology and tissue structure in vivo are generally irregular. Tumor tissues have further chaotic distribution of blood vessels. Extending the model to three-dimensional description can make the scenario more realistic and help study cell morphology, mechanical stress, and tissue growth related issues in future.

Funding Data

  • National Science Foundation (Grant No. DMS-1317671).

References

1.
Semenza
,
G. L.
,
2000
, “
Hypoxia, Clonal Selection, and the Role of HIF-1 in Tumor Progression
,”
Crit. Rev. Biochem. Mol. Biol.
,
35
(
2
), pp.
71
103
.
2.
Kaelin
,
W. G.
, and
Ratcliffe
,
P. J.
,
2008
, “
Oxygen Sensing by Metazoans: The Central Role of the HIF Hydroxylase Pathway
,”
Mol. Cell
,
30
(
4
), pp.
393
402
.
3.
Gilkes
,
D. M.
,
Semenza
,
G. L.
, and
Wirtz
,
D.
,
2014
, “
Hypoxia and the Extracellular Matrix: Drivers of Tumour Metastasis
,”
Nat. Rev. Cancer
,
14
(
6
), pp.
430
439
.
4.
Hirsila
,
M.
,
Koivunen
,
P.
,
Gunzler
,
V.
,
Kivirikko
,
K. I.
, and
Myllyharju
,
J.
,
2003
, “
Characterization of the Human Prolyl 4-Hydroxylases That Modify the Hypoxia-Inducible Factor
,”
J. Biol. Chem.
,
278
(
33
), pp.
30772
30780
.
5.
Tuckerman
,
J. R.
,
Zhao
,
Y. G.
,
Hewitson
,
K. S.
,
Tian
,
Y. M.
,
Pugh
,
C. W.
,
Ratcliffe
,
P. J.
, and
Mole
,
D. R.
,
2004
, “
Determination and Comparison of Specific Activity of the HIF-Prolyl Hydroxylases
,”
FEBS Lett.
,
576
(
1–2
), pp.
145
150
.
6.
LaGory
,
E. L.
, and
Giaccia
,
A. J.
,
2016
, “
The Ever-Expanding Role of HIF in Tumour and Stromal Biology
,”
Nat. Cell Biol.
,
18
(
4
), pp.
356
365
.
7.
Pappas
,
D.
,
2016
, “
Microfluidics and Cancer Analysis: Cell Separation, Cell/Tissue Culture, Cell Mechanics, and Integrated Analysis Systems
,”
Analyst
,
141
(
2
), pp.
525
535
.
8.
Khanal
,
G.
,
Chung
,
K.
,
Solis-Wever
,
X.
,
Johnson
,
B.
, and
Pappas
,
D.
,
2011
, “
Ischemia/Reperfusion Injury of Primary Porcine Cardiomyocytes in a Low-Shear Microfluidic Culture and Analysis Device
,”
Analyst
,
136
(
17
), pp.
3519
3526
.
9.
Kohn
,
K. W.
,
Riss
,
J.
,
Aprelikova
,
O.
,
Weinstein
,
J. N.
,
Pommier
,
Y.
, and
Barrett
,
J. C.
,
2004
, “
Properties of Switch-Like Bioregulatory Networks Studied by Simulation of the Hypoxia Response Control System
,”
Mol. Biol. Cell
,
15
(
7
), pp.
3042
3052
.
10.
Qutub
,
A. A.
, and
Popel
,
A. S.
,
2006
, “
A Computational Model of Intracellular Oxygen Sensing by Hypoxia-Inducible Factor HIF1 Alpha
,”
J. Cell Sci.
,
119
(
16
), pp.
3467
3480
.
11.
Nguyen
,
L. K.
,
Cavadas
,
M. A. S.
,
Scholz
,
C. C.
,
Fitzpatrick
,
S. F.
,
Bruning
,
U.
,
Cummins
,
E. P.
,
Tambuwala
,
M. M.
,
Manresa
,
M. C.
,
Kholodenko
,
B. N.
,
Taylor
,
C. T.
, and
Cheong
,
A.
,
2013
, “
A Dynamic Model of the Hypoxia-Inducible Factor 1 Alpha (HIF-1 Alpha) Network
,”
J. Cell Sci.
,
126
(
6
), pp.
1454
1463
.
12.
Morshed
,
A.
, and
Dutta
,
P.
,
2017
, “
Hypoxic Behavior in Cells Under Controlled Microfluidic Environment
,”
Biochim. Biophys. Acta (BBA)—Gen. Subj.
,
1861
(
4
), pp.
759
771
.
13.
Sze
,
T. K. J.
,
Liu
,
J.
, and
Dutta
,
P.
,
2014
, “
Numerical Modeling of Flow Through Phloem Considering Active Loading
,”
ASME J. Fluids Eng.
,
136
(
2
), p. 021206.
14.
Khanal
,
G.
,
Hiemstra
,
S.
, and
Pappas
,
D.
,
2014
, “
Probing Hypoxia-Induced Staurosporine Resistance in Prostate Cancer Cells With a Microfluidic Culture System
,”
Analyst
,
139
(
13
), pp.
3274
3280
.
15.
Vaupel
,
P.
,
Mayer
,
A.
, and
Hockel
,
M.
,
2004
, “
Tumor Hypoxia and Malignant Progression
,”
Oxygen Sens.
,
381
, pp.
335
354
.
16.
Spaeth
,
E. E.
, and
Friedlander
,
S. K.
,
1967
, “
The Diffusion of Oxygen, Carbon Dioxide, and Inert Gas in Flowing Blood
,”
Biophys. J.
,
7
(
6
), pp.
827
851
.
17.
Sasaki
,
N.
,
Horinouchi
,
H.
,
Ushiyama
,
A.
, and
Minamitani
,
H.
,
2012
, “
A New Method for Measuring the Oxygen Diffusion Constant and Oxygen Consumption Rate of Arteriolar Walls
,”
Keio J. Med.
,
61
(
2
), pp.
57
65
.
18.
Grote
,
J.
,
Süsskind
,
R.
, and
Vaupel
,
P.
,
1977
, “
Oxygen Diffusivity in Tumor Tissue (DS-Carcinosarcoma) Under Temperature Conditions Within the Range of 20–40 C
,”
Pflügers Arch. Eur. J. Physiol.
,
372
(
1
), pp.
37
42
.
19.
Bassingthwaighte
,
J. B.
,
Kuikka
,
J. T.
,
Chan
,
I. S.
,
Arts
,
T.
, and
Reneman
,
R. S.
,
1985
, “
A Comparison of Ascorbate and Glucose-Transport in the Heart
,”
Am. J. Physiol.
,
249
(
1
), pp.
H141
H149
.http://ajpheart.physiology.org/content/249/1/H141.long
20.
Kuiper
,
C.
,
Vissers
,
M. C. M.
, and
Hicks
,
K. O.
,
2014
, “
Pharmacokinetic Modeling of Ascorbate Diffusion Through Normal and Tumor Tissue
,”
Free Radical Biol. Med.
,
77
, pp.
340
352
.
21.
Herst
,
P. M.
, and
Berridge
,
M. V.
,
2007
, “
Cell Surface Oxygen Consumption: A Major Contributor to Cellular Oxygen Consumption in Glycolytic Cancer Cell Lines
,”
Biochim. Biophys. Acta, Bioenerg.
,
1767
(
2
), pp.
170
177
.
22.
Gatenby
,
R. A.
, and
Gillies
,
R. J.
,
2004
, “
Why Do Cancers Have High Aerobic Glycolysis?
,”
Nat. Rev. Cancer
,
4
(
11
), pp.
891
899
.
23.
Chen
,
Y. J.
,
Cairns
,
R.
,
Papandreou
,
I.
,
Koong
,
A.
, and
Denko
,
N. C.
,
2009
, “
Oxygen Consumption Can Regulate the Growth of Tumors, a New Perspective on the Warburg Effect
,”
PLoS One
,
4
(
9
), p. e7033.
24.
Figueroa-Mendez
,
R.
, and
Rivas-Arancibia
,
S.
,
2015
, “
Vitamin C in Health and Disease: Its Role in the Metabolism of Cells and Redox State in the Brain
,”
Front. Physiol.
,
6
, p. 397.https://www.ncbi.nlm.nih.gov/pubmed/26779027
25.
Graumlich
,
J. F.
,
Ludden
,
T. M.
,
ConryCantilena
,
C.
,
Cantilena
,
L. R.
,
Wang
,
Y. H.
, and
Levine
,
M.
,
1997
, “
Pharmacokinetic Model of Ascorbic Acid in Healthy Male Volunteers During Depletion and Repletion
,”
Pharm. Res.
,
14
(
9
), pp.
1133
1139
.
26.
Padayatty
,
S. J.
,
Sun
,
H.
,
Wang
,
Y. H.
,
Riordan
,
H. D.
,
Hewitt
,
S. M.
,
Katz
,
A.
,
Wesley
,
R. A.
, and
Levine
,
M.
,
2004
, “
Vitamin C Pharmacokinetics: Implications for Oral and Intravenous Use
,”
Ann. Intern. Med.
,
140
(
7
), pp.
533
537
.
27.
Myllylä
,
R.
,
Tuderman
,
L.
, and
Kivirikko
,
K. I.
,
1977
, “
Mechanism of the Prolyl Hydroxylase Reaction
,”
Eur. J. Biochem.
,
80
(
2
), pp.
349
357
.
28.
Hirsilä
,
M.
,
2004
, “
Characterization of the Novel Human Prolyl 4-Hydroxylases and Asparaginyl Hydroxylase That Modify the Hypoxia-Inducible Factor
,”
Doctoral dissertation
, University of Oulu, Oulu, Finland.http://jultika.oulu.fi/Record/isbn951-42-7575-6
29.
Patankar
,
S.
,
1980
,
Numerical Heat Transfer and Fluid Flow
,
Hemisphere Publishing Corporation
,
New York
.
30.
Smallbone
,
K.
,
Gatenby
,
R. A.
,
Gillies
,
R. J.
,
Maini
,
P. K.
, and
Gavaghan
,
D. J.
,
2007
, “
Metabolic Changes During Carcinogenesis: Potential Impact on Invasiveness
,”
J. Theor. Biol.
,
244
(
4
), pp.
703
713
.
31.
Milotti
,
E.
, and
Chignola
,
R.
,
2010
, “
Emergent Properties of Tumor Microenvironment in a Real-Life Model of Multicell Tumor Spheroids
,”
PLoS One
,
5
(
11
), p. e13942.
32.
Gatenby
,
R. A.
,
Smallbone
,
K.
,
Maini
,
P. K.
,
Rose
,
F.
,
Averill
,
J.
,
Nagle
,
R. B.
,
Worrall
,
L.
, and
Gillies
,
R. J.
,
2007
, “
Cellular Adaptations to Hypoxia and Acidosis During Somatic Evolution of Breast Cancer
,”
Br. J. Cancer
,
97
(
5
), pp.
646
653
.

Supplementary data