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.
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 . 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 .
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 . 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 . In recent years, microfluidic cell culture devices have been used to study cancer cell growth, cell cycle, apoptosis, and other mechanisms . 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 . 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.  hypothesized about a switchlike behavior in HIF activity under decreasing oxygen gradient which has been later verified by Qutub and Popel . Nguyen et al.  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 . 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.
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.
where is the instantaneous spatial concentration of oxygen, is the medium fluid velocity, is the diffusion coefficient for oxygen in the medium, and 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 .
where is the instantaneous spatial concentration of ascorbate, is the diffusion coefficient for ascorbate in the medium, and 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.
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 . 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.
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 . However, hypoxic tumor cells have been known to reduce this demand of oxygen by selectively shutting off the mitochondrial activity . 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 . 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 . 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.
Here, and denote the rate of ascorbate uptake and the extracellular ascorbate concentration, respectively, and 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.
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.
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. . The sequential binding of collagen prolyl hydroxylase during hydroxylation has been identified to proceed in the order of iron → 2-oxoglutarate → oxygen → peptide . The same sequence of binding has been suggested for HIF hydroxylation by PHD . 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).
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 “ ,” and the width of the cell block is “ .” The ratio of to 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.
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  and the cytosol . 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. ), in the present study we chose to use the experimentally reported oxygen consumption rates from Ref.  as the sink term for oxygen in the cell-block region after multiplying it with the volumetric cell density (). 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 , 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.
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 -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 . 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 -direction. Our results, in the dynamic equilibrium, were compared with model results from Ref. , 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. .
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 .
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.
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. . 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.
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.
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.
National Science Foundation (Grant No. DMS-1317671).