Oxygen transport and stem cell aggregation in stirred-suspension bioreactor cultures.

Stirred-suspension bioreactors are a promising modality for large-scale culture of 3D aggregates of pluripotent stem cells and their progeny. Yet, cells within these clusters experience limitations in the transfer of factors and particularly O2 which is characterized by low solubility in aqueous media. Cultured stem cells under different O2 levels may exhibit significantly different proliferation, viability and differentiation potential. Here, a transient diffusion-reaction model was built encompassing the size distribution and ultrastructural characteristics of embryonic stem cell (ESC) aggregates. The model was coupled to experimental data from bioreactor and static cultures for extracting the effective diffusivity and kinetics of consumption of O2 within mouse (mESC) and human ESC (hESC) clusters. Under agitation, mESC aggregates exhibited a higher maximum consumption rate than hESC aggregates. Moreover, the reaction-diffusion model was integrated with a population balance equation (PBE) for the temporal distribution of ESC clusters changing due to aggregation and cell proliferation. Hypoxia was found to be negligible for ESCs with a smaller radius than 100 µm but became appreciable for aggregates larger than 300 µm. The integrated model not only captured the O2 profile both in the bioreactor bulk and inside ESC aggregates but also led to the calculation of the duration that fractions of cells experience a certain range of O2 concentrations. The approach described in this study can be employed for gaining a deeper understanding of the effects of O2 on the physiology of stem cells organized in 3D structures. Such frameworks can be extended to encompass the spatial and temporal availability of nutrients and differentiation factors and facilitate the design and control of relevant bioprocesses for the production of stem cell therapeutics.

Introduction

Realization of the therapeutic potential of pluripotent stem cells (PSCs) including embryonic stem cells (ESCs) and induced pluripotent stem cells (iPSCs) hinges on the development of platforms for large-scale PSC expansion and directed differentiation. The culture of PSCs has been demonstrated in laboratory scale stirred-suspension bioreactors (SSBs) as cell aggregates [1] or on various scaffold types [2], [3]. Compared to static culture, SSBs afford higher cell densities and tighter control of process variables leading to more efficient utilization of media and growth factors [1].

Despite the advantages of culturing PSCs as aggregates in stirred-suspension vessels, less attention has been devoted to limitations in the transfer of medium components to and within PSC clusters. Concentration gradients in the bioreactor environment contribute to stem cell population heterogeneity leading to variable responses to self-renewal or differentiation stimuli. Such unevenness is more pronounced for oxygen characterized by low solubility in aqueous solutions and increased rate of consumption by metabolically active stem cells. Oxygen availability affects directly the viability, proliferation and differentiation propensity of stem cells in vivo and in vitro [4]. Cultivated human ESCs (hESCs) are less prone to spontaneous differentiation and chromosomal aberrations at hypoxic (2–5%) than normoxic (21%) O2 level without significantly reduced proliferation [5], [6]. Lower O2 tension (pO2) in culture also predisposes stem cells to commit along particular lineages including endothelial cells [7] and chondrocytes [8] similar to in vivo processes [4]. More details emerge about O2 modulation of the activity of stem cell-fate controlling pathways such as the canonical Wnt/β-catenin [9] and Notch cascades [10]. These effects are largely mediated by transcriptional regulators such as the hypoxia-inducible factors (HIF) interacting with a wide gamut of genes including the pluripotency marker POU5F1 (Oct4) [11]. Hence, knowledge of O2 profile among stem or progenitor cells is important for cell fate prediction and control.

The distribution of O2 was reported for single hESC aggregates under static conditions [12], [13]. The findings from these studies however may not apply directly to aggregate cultures of ESCs in SSBs where the transfer of O2 through multiple interfaces complicates the analysis of the O2 tension (pO2) that each cell experiences. In most small-scale setups, O2 is transferred via the gas/liquid interface (headspace aeration) to the culture medium under agitation. From the medium bulk O2 fluxes to cells via a boundary layer surrounding each aggregate and a pore network within each cluster. These ultrastructural characteristics of aggregates have also not been considered explicitly to date. And although a steady state assumption simplifies the mathematical framework for O2 transfer in stem cell aggregate culture, it may not proper in the case of SSB cultures. This is because the continuous cell proliferation and dynamics of agitation-induced aggregation result in the temporal modulation of O2 levels. Furthermore, eliminating time as a variable in the analysis of mass transfer does not allow the calculation of the duration that a particular cell (or fraction) experiences a pO2 below a certain threshold. Differences in the “residence time” at a certain O2 concentration among cells may be consequential for the differentiation and/or self-renewal of the stem cell population.

Here, the distribution of O2 was calculated from experimental data linked with a mathematical model for mouse ESC (mESC) and hESC aggregates cultivated in dishes and spinner flasks. Aggregates cultured in spinner flasks under different agitation rates were analyzed at different time points and ultrastructural attributes such as porosity and tortuosity were determined for the first time. The effective diffusivity of O2 within the aggregates and parameters associated with the specific rate of O2 consumption were also computed by coupling measurements with a transient diffusion-reaction model. This model was subsequently paired to a population balance equation (PBE) depicting the time evolution of aggregate size distribution of proliferating ESCs under agitation. The multiscale PBE/diffusion-reaction model facilitated the calculation of the O2 profile in the medium bulk and inside the aggregates in stirred suspension. Moreover, the fraction of cells experiencing hypoxia was predicted. The computational model described in this study can be utilized in the development of differentiation strategies and the design and control of relevant bioprocesses for the production of stem cell therapeutics.

Material and Methods

Experimental methods

Embryonic stem cell culture

Human ESCs (H9 cells, WiCell, Madison, WI; passages 30–60) were maintained on dishes coated with Matrigel (BD Biosciences, San Jose, CA) and grown in 5% CO2/95% air at 37°C with mTESR1 medium (Stem Cell Technologies, Vancouver, BC) replaced daily. Four hours before harvesting, cells were treated with the Rho-associated kinase inhibitor (ROCKi) Y-27632 (10 µM, EMD Chemicals, Gibbstown, NJ). Colonies were dissociated with Accutase (Innovative Cell Technologies, San Diego, CA) into single hESCs. For static aggregate culture, cells were transferred to poly HEMA-treated (Sigma-Aldrich, St. Louis, MO) Petri dishes with mTESR1. Dispersed hESCs (105/ml) were seeded in 125-ml ProCulture spinner flasks (Corning, Corning, NY) with mTESR1 medium and 10 µM ROCKi. Agitation was kept constant during each experiment and medium was replaced daily thereafter without the addition of ROCKi.

Mouse E14Tg2a ESCs (Mutant Mouse Regional Resource Centers (MMRRC), University of California-Davis, CA) were maintained in defined serum-free medium (DSFM) as described [14] (see ). For stirred-suspension culture [14], mESCs were transferred in DSFM to 125-ml ProCulture spinner flasks at 5×104 cells/ml. The agitation was kept constant (60–100 rpm) throughout each run.

Oxygen consumption

Cell aggregate samples of known volumes were transferred to 15 ml-centrifuge tubes and allowed to settle. Without removing aggregates, the medium was replaced by PBS (pH 7.4) containing 0.35 g/L HEPES, 0.5 g/L bovine serum albumin, and 300 mg/dL glucose (Sigma-Aldrich, St. Louis, MO) with a Bunsen solubility coefficient of αm = 1.27×10−9 mol/(cm3 mm Hg) at 37°C [15], [16]. The cell suspension was transferred to a 30-ml cylindrical container (2.3 cm wide) with a magnetic cylindrical stirring bar (1.59 cm×0.79 cm). The container was set on a submersible magnetic stirrer plate in a water bath at 37°C bath and agitation was set to 60 rpm. A dissolved O2 (DO) polarographic probe (VWR, Bridgeport, NJ) was inserted in the container through the cap with a rubber seal ensuring airtight fitting. DO data (mg/l) were logged on a computer connected to the DO probe/meter. Then, aggregates were transferred to a 96-well plate and images were taken and processed (see below) to obtain the total number of aggregates and pertinent size distributions [17].

Stem cell aggregate micrograph acquisition and analysis

Samples of known volume from spinner flasks, dishes or the DO chamber were transferred to a 96-well plate and images from each well were acquired at 4x magnification with a camera (Moticam 2300, Motic, Richmond, BC) connected to the microscope. Raw images were imported into ImageJ (U.S. National Institutes of Health, Bethesda, Maryland, USA, http://imagej.nih.gov/ij/) and for counting total numbers of aggregates. After background subtraction (rolling ball method), the diameter of each aggregate was calculated as the mean of two perpendicular diameters and size distributions were generated [17].

Porosity ε and tortuosity τ were determined in aggregates transferred to culture medium containing 10 mg/ml fluorescein isothiocyanate (FITC)-dextran (MW: 4.4 kDa, Sigma). The aggregate suspension was then transferred to a microscopy slide with a well, sealed with a coverslip and visualized by confocal laser scanning microscopy (CLSM) with a Zeiss LSM 510 Metal NLO system (Zeiss, Thornwood, NY). Aggregates were optically sectioned (10–12 sections per aggregate) at 5-µm steps and the micrographs were processed in MATLAB (Mathworks, Natick, MA) to obtain porosity and tortuosity values ( ; see ).

Figure 1

Structural characteristics (tortuosity and porosity) of mESC and hESC aggregates under different conditions.

(A) Optical sections were acquired for each aggregate via CLSM and processed to obtained tortuosity and porosity values. For each condition, at least 10 aggregates were imaged by CLSM (30 sections/aggregate) for data analysis. Tortuosity (lines) and porosity (bars) data are shown for (B) mESC aggregates at different time points at 100 rpm, (C) mESC aggregates at different agitation rates on day 3, and (D) hESC aggregates at different agitation rates on day 4. The values in (B–D) are shown as mean±st. dev.

Statistical Analysis

Data were expressed as mean±st.dev. unless stated otherwise. ANOVA and the post hoc Tukey test were performed using Minitab (Minitab Inc, State College, PA) with p<0.05 considered as significant.

Modeling methods

Transient reaction-diffusion model

For the transient reaction-diffusion model employed in this study, the aggregates were assumed to be spherical and thus spherical coordinates were employed to describe the O2 concentration (C) profile in each aggregate:where r is the aggregate radius, t is the time and is the effective diffusion coefficient. Taking into account the porosity ε and tortuosity τ, is calculated aswhere Dt is the diffusion of O2 in the tissue (aggregate).

The consumption of O2 (in mol/(cm3.s)) was described by Michaelis-Menten (M-M)-type kinetics:where Vmax is the maximum O2 consumption rate and KM is the M-M constant. Because the system comprises aggregates of multiple sizes, cell clusters were classified into groups or bins. Here, the utilization of 7 bins provided sufficient resolution of the aggregate size distribution. Equation 1 was solved for each group [18] (see ) subject to appropriate initial and boundary conditions.

where R is the radius for the aggregates in the l th bin. Calculation of the value of the O2 transfer coefficient k from the bulk liquid to the aggregates was based on the Frössling equation [19] relating the particle Sherwood () and Schmidt () numbers shown in the section.

The total rate of O2 transfer from the bulk to all aggregates ( groups) is:where and are the fraction of aggregates with average diameter and the total number of aggregates, respectively. Integrals were calculated via the trapezoidal rule using C values from solving Equation 1 and the Euler method was implemented yield estimates of the O2 concentration in the medium bulk at the -time step in Equation 7.

Values for , and were obtained by minimizing the objective function [20] comparing the calculated values with the experimental ones at each time point:

Development of PBE model for the bioreactor culture of ESC aggregates

A size-structured PBE model was employed to simulate the temporal evolution of size of ESC aggregates cultured in spinner flasks [1] considering the contributions of ESC proliferation and collisions between particles to aggregate size growth:

The density function is defined such that is the fraction of aggregates with sizes between and at time per unit volume of the culture. The size corresponds either to the aggregate volume or mass since the buoyant density of cells does not vary significantly during the cell cycle [21], [22]. The mass of individual cells was taken as and . Aggregate breakage was not observed particularly after the first day of culture and therefore was not considered.

The rate of aggregate size change, , due to cell proliferation was modeled after the Gompertz equation [23]:where M is the aggregate mass reached as and is a constant characteristic of the cell proliferation. Both parameters were estimated from temporal data of aggregates in static cultures.

Different aggregation kernels have been reported in various systems [24], [25]. Given the non-rigid nature of cells/aggregates, we applied here an aggregation kernel (Eq. 11) for liquid-liquid dispersions assuming dynamic deformation of the colliding bodies [26]. The coalescence efficiency (exponential term) decreases with increasing average size of the colliding aggregates. The parameters , and are obtained from ESC aggregate size distributions measured in spinner flask cultures.

The PBE was expressed with respect to the aggregate radius, which was measured experimentally, assuming spherical aggregates. However, the density function was transformed in terms of the logarithm of the aggregate radius normalized to that of a single cell (∼7.5 µm) to allow for an extended range of aggregate sizes to be considered (see [17] and ).

The integro-partial differential PBE was solved numerically by finite (backward) differences over equally spaced nodes in the lnr domain. Kernel parameters for different agitation rates were obtained from aggregate size data. The total number of aggregates and biomass were calculated as the zeroth and first moments of the cell distribution.

Integration of PBE and reaction-diffusion models

The PBE and reaction-diffusion models were combined for predicting the O2 profile in stirred-suspension cultures of ESC aggregates ( ). The total number of aggregates , size distribution of aggregates , oxygen concentration in the medium bulk , working volume of bioreactor (assumed equal to that of the medium) and agitation rate were used for initialization of the integrated model (see ). Cell growth and aggregation were modeled (Eq. 10–11) and the consumption of O2 by all the aggregates in the bioreactor was determined (Eq. 7). This required knowledge of the O2 concentration distribution within individual aggregates (Eq. 1). Subsequently, the bulk O2 level was updated based on the net effect of O2 consumption and supply:

The rate of O2 supply through headspace aeration was:where is the O2 partial pressure in the gas phase (95% air/5% CO2, 37°C, 760 mm Hg), is the oxygen partial pressure in bulk medium, is the mass transfer coefficient at the gas-liquid interface. Values of were determined experimentally (see ).

Results

Porosity and tortuosity of ESC aggregates

The transport of O2, nutrients and pluripotency-maintaining or differentiation-inducing factors to and within stem cell aggregates is greatly influenced by the aggregate ultrastructure. Moreover, whether the intra-aggregate ultrastructure is strongly affected by culture conditions in SSBs -particularly agitation- remains unknown. Thus, before proceeding with the solution of the diffusion-reaction model for ESC aggregates, we used CLSM and image analysis to determine two ultrastructural descriptors –porosity (ε) and tortuosity (τ)– for mESC and hESC aggregates formed in dishes or spinner flasks at different stirring rates ( ). The porosity of mESC aggregates in spinner flasks at 100 rpm was 0.229±0.016 at 24 hours ( ) and decreased slightly after 3 (0.218±0.011) and 5 days (0.220±0.011) of culture suggesting that aggregates were undergoing compaction. For 3-day cultures the ε values for 60 and 80 rpm were 0.217±0.02 and 0.210±0.025, respectively ( ). When cultured in dishes, mESCs formed looser aggregates with ε = 0.259±0.041 (day 3). Yet, no statistically significant differences were noted in the porosity between aggregates cultured in dishes and spinner flasks at different agitation rates. Human ESCs in static culture formed aggregates with a porosity of 0.254±0.017 ( ) whereas hESC aggregates from spinner flasks at 45 and 60 rpm displayed porosities of 0.258±0.011 and 0.270±0.007, respectively.

The tortuosity of mESC aggregates at 100 rpm was 1.427±0.083 at 24 hours and did not change significantly after 3 (1.504±0.093) and 5 days (1.547±0.181) of culture ( ). Moreover, clusters formed at 60 rpm (1.460±0.035), 80 rpm (1.534±0.086) and in dishes (1.457±0.422) exhibited similar tortuosity. Human ESC aggregates were characterized by slightly higher average τ values than mESC clusters ( ). After 4 days in stirred suspension an average τ of 1.612±0.135 (45 rpm) or 1.551±0.086 (60 rpm) was determined whereas the tortuosity of dish aggregates was 1.638±0.136.

Our findings provide a first account of porosity and tortuosity values determined directly for mESC and hESC aggregates formed in dish and stirred-suspension cultures at different times and under typical agitation rates. The data did not show strong dependence of the ultrastructural descriptors considered here on the cultivation mode, i.e. dish vs. stirred suspension nor on the stirring speeds tested.

Effective diffusion and consumption of O2

A transient reaction-diffusion model was employed to determine the values of the diffusion coefficient and consumption () of O2 in the ESC aggregates from experimentally measured O2 concentration values in the medium. First, the mean size was determined of aggregates in dishes and stirred-suspension cultures at different agitation rates ( ). At the end of the culture, mESC aggregates at 100 rpm were smaller (mean diameter: 113.04±15.64 µm, p<0.05) than those at 60 rpm (173.8±23.9 µm) and 80 rpm (169.01±16.13 µm). The mean diameter of hESC aggregates at 45 and 60 rpm at day 4 was 163.8 and 192.45 µm, respectively and similar to those in static cultures (166.61±12.4 µm) ( ). For each sample, the number fraction distribution was determined ( ).

Figure 2

Experimental measurement of ESC aggregate size and O2 consumption profiles for static and spinner flask cultures.

Mean aggregate diameter data are shown for dish and bioreactor cultures at different agitation rates for (A) mESCs (*p<0.05 vs. 60 and 80 rpm) and (B) hESCs. (C) A representative ESC aggregate size distribution is shown (mESC culture at 100 rpm). (D) Experimental (solid line) and simulation (dotted line) results for the consumption of O2 over time. Parameter (Vmax, D, and KM) values were acquired by fitting the transient reaction-diffusion model to the experimental data for different agitation rates (Table 1).

Table 1

Parameter values with a critical concentration of 8 mm Hg.

Culture condition	Dt (×105 cm2/s)	Vmax (×108 mol/cm3.s)	KM (×108 mol/cm3)
mESCs
Static culture (aggregates)	1.631±0.19 [0.242±0.024]	2.615±0.997 [4.732±0.318]	2.81±1.904 [0.518±0.165]
60 rpm	1.716±0.223 [0.364±0.02]	2.667±0.378*,# [5.055±0.163]	2.721±0.385 [0.782±0.329]
80 rpm	1.283±0.864 [0.418±0.05]	3.103±1.03 [3.602±1.06]	4.575±1.915 [0.392±0.233]
100 rpm	2.375±0.199 [0.386±0.114]	4.831±0.593 [7.322±0.83]	4.063±0.693 [0.563±0.06]
hESCs
Static culture (aggregates)	1.20±0.170 [0.217±0.022]	0.995±0.464 [2.029±0.478]	3.473±0.232 [0.062±0.033]
45 rpm	1.536±0.322 [0.547±0.093]	1.16±0.355 [1.591±0.963]	3.23±0.257 [1.564±0.155]
60 rpm	1.163±0.454 [0.69±0.223]	1.36±0.427 [1.266±0.542]	4.126±1.138 [0.876±0.443]

Values in brackets are calculated with a 0 mm Hg critical concentration.

*p = 0.005 [0.019] for Vmax at 60 vs. 100 rpm,

p = 0.032 [0.001] for Vmax for mESCs vs. hESCs at 60 rpm.

Cells within aggregates consumed O2 resulting in a gradual decrease of the DO level in the medium ( ). Typically, DO measurements were performed until the O2 concentration in the medium dropped below 10 mm Hg. Solution of the model coupled with experimental data on DO, yielded values for and ( ). It should be noted that parameter values were calculated for a minimum (critical) O2 concentration at the aggregate center of either (i) 0 mm Hg, or (ii) 8 mm Hg [27]. Irrespective of the minimum O2 level considered, the values of (see Eq. 2) were lower than that of the diffusion coefficient of O2 in pure water at 37°C (2.78×10−5 cm2/s). As a measure of maximum O2 consumption rate, Vmax was higher for mESCs in suspension culture at 60 rpm than for hESCs (p<0.05). Mouse ESCs also displayed greater Vmax at 100 vs. 60 rpm suggesting an effect of the agitation on O2 metabolism (p<0.05). The mean Vmax for hESCs at 45 and 60 rpm did not differ significantly.

Oxygen profile within ESC aggregates in static and stirred-suspension cultures

Once the values of and were determined, the reaction-diffusion model was implemented to calculate the O2 profile for ESC aggregates of different size. The diffusion of O2 becomes limited as the aggregate radius increases leading to reduced O2 levels for cells residing closer to the core. To better illustrate this, the O2 distribution was calculated in hESC aggregates sampled from spinner flask cultures assuming a constant O2 concentration in the bulk of 148 mm Hg (saturation) ( ). Because hypoxia studies typically employ pO2 of ∼30 mm Hg (4–5% O2) [5], we calculated the fractions of ESCs exposed to O2 below this level for different aggregate sizes. Human ESC aggregates smaller than 200 µm were normoxic throughout their volume. In contrast, almost 23% of the hESCs in aggregates larger than 400 µm were under hypoxia. This fraction increased to 70% for aggregates with 1,000 µm diameter ( ). The O2 profiles in hESC aggregates of various sizes maintained in suspension (60 rpm) or static culture are shown in . Similarly, the O2 distributions within mESC aggregates cultured at 60 rpm or statically are also shown ( ). It should be noted that the results shown in were obtained assuming that all aggregates in culture are the same size. Simulation results for additional agitation rates are shown in .

Figure 3

Steady-state O2 profile within ESC aggregates of varying sizes.

(A) Three-dimensional distribution of O2 and hypoxic cell fraction of hESC aggregates with radius ranging from 20 µm to 500 µm. The bulk oxygen level was taken as 148 mm Hg. Fractions of ESCs exposed to O2 below a particular threshold (y-axis) for different size aggregates cultured at 60 rpm or in dishes (static) shown for (B) hESCs and (C) mESCs. Results were obtained assuming that all cultured aggregates were the same size. The minimum intra-aggregate O2 level was fixed at 8 mm Hg (critical concentration).

In general, the fractions of cells experiencing an O2 concentration below a particular threshold were lower in stirred suspension compared to dish cultures for given size clusters. For instance, the hypoxic cell fractions are larger in the plots from dishes compared to spinner flask cultures for the same size aggregates ( ). This is most likely due to the faster transfer () of O2 from the medium bulk to each aggregate in agitated vessels compared to static cultures (Eq. S2). Practically however, we rarely observed aggregates with a radius over 300 µm and only in static cultures. The higher hindrance for O2 (and nutrient) transport to large clusters possibly impacts ESC viability and reduces proliferation by stimulating differentiation. This and the agitation-induced shear in spinner flask cultures may limit further growth of bigger aggregates.

Aggregates from stirred-suspension cultures showed similar O2 profiles across all the agitation rates we tested (). As the aggregate size increased, cells positioned over 100 µm from the aggregate surface experienced less than 10 mm Hg of O2. Compared to hESCs, the fractions of mESCs under the same O2 concentration threshold were greater for equal size aggregates owing most likely to the higher Vmax of the latter. The difference was most significant for aggregates with a radius between 100–250 µm.

Nevertheless, the expansion of mESCs and hESCs in spinner flasks gave rise to aggregate diameters mostly below 200–300 µm ( ). Thus, if a normoxic concentration of O2 can be maintained in the bulk only 2.2% of the cells (located in the aggregate center) experience O2 tension below 30 mm Hg ( ). In fact, when we analyzed cells by qPCR for HIF1α expression () there was no notable differences between static and stirred-suspension cultures. Also, compared to monolayer cultures at normoxia, HIF1α expression increased only slightly for mESCs (<2.0-fold) and hESCs (<1.5-fold) within aggregates.

PBE model for ESC aggregation in stirred-suspension culture

For the O2 profiles presented above, the O2 concentration in the medium was assumed to be fixed at saturation (∼148 mm Hg). This assumption is valid if O2 can be supplied efficiently and continuously to the cultivation system or the culture is sparse so that O2 consumption is trivial. However, stirred-suspension cultures accommodate high cell densities resulting in fluctuating O2 levels in the bioreactor. Agitation and proliferation also drive the formation and growth of aggregates over time with concomitant changes in the intra- and inter-aggregate concentrations of O2. Given the multitude of effects of O2 on stem cell physiology, predicting the time-variant distribution of O2 in scalable ESC aggregate cultures is highly desirable.

To that end, a PBE model was developed for the temporal evolution of the ESC cluster size distribution in an SSB. The growth rate of ESC aggregates due to cell proliferation was described by the Gompertz equation. For this purpose, the growth of sparsely cultured mESC aggregates in dishes was monitored and the Gompertz parameters were evaluated (M = 9.71×106, αG = 5.72×10−3 hr−1) based on the recorded size changes ( ). The Gompertz model accurately described the growth of mESC aggregates (R2 = 0.99, n = 20). The maximum attainable size for mESC aggregates was calculated at 1600 µm ().

Figure 4

PBE model for the ESC aggregate culture in spinner flasks.

(A) Mouse ESC aggregates in dishes were monitored and their size was measured over several days. The experimental data for aggregate growth due to mESC proliferation were compared with the Gompertz model. Relevant parameters were determined: M = 9.71×106, αG = 5.72×10−3 hr−1. (B) Representative data of size (radius) distribution for mESC aggregates cultured at 60 rpm (day 1). Histogram and a continuous distribution curve are shown. (C) Model size distributions (lines) shown for a 96-hour culture at 60 rpm. The original distributions (dash-circle lines) for each time point are shown for comparison.

Parameters for the aggregation kernel were determined from size data of mESC aggregates cultured in spinner flasks ( ). It should be noted that the kernel parameters were calculated from data of cell cultures at 0–72 hours ( ). At longer times, aggregation was negligible and cell proliferation drove the increase in cluster size. With increasing agitation rate, the parameter value trends are such that the coalescence frequencies (i.e. collisions leading to aggregate formation) become lower. At 60 rpm, the mean aggregate radius increased to ∼105 µm on day 2 (mainly due to collision-driven aggregation) in line with the model predictions ( ). Subsequent increase in aggregate size was less pronounced as the mean radius was approximately 150 µm (day 4). The corresponding average size at higher agitation rates was lower (e.g. ∼90 µm at 100 rpm, ). The simulation results were corroborated by observations that a lower agitation rate generally promotes formation of larger aggregates. It should also be noted that the cell viability was not affected significantly by agitation (>90%).

Table 2

Kernel parameter values at different agitation rates.

Agitation rate (rpm)	k (hr−1 µm−7/3)	k1 (µm−3α)	α (unitless)
60	1.26±0.77×10−3 *	1.94±1.12×10−4	8.06±0.79×10−1 *
80	6.36±4.47×10−6	4.62±3.25×10−4	5.96±0.085×10−1 *
100	6.53±6.23×10−6	4.11±0.63×10−5	9.91±0.65×10−1 *

The parameter values were shown as mean ± st.dev. calculated by model simulation to experimental data from mESC cultures (at least 3 independent runs per agitation rate).

*p<0.05 when compared to the mean values obtained for the other agitation rates.

Taken together, the PBE model developed here allowed the prediction of the temporal evolution of the ESC aggregate size distribution (and the total aggregate number and biomass; see below). Simulation results were in good agreement with the experimental observations from stirred-suspension ESC culture at different agitation rates.

Prediction of O2 profile and hypoxia fraction in ESC aggregate stirred-suspension culture

With the PBE model of ESC aggregation in the SSB in place, the O2 distribution in the medium and within ESC aggregates was computed by coupling the diffusion reaction and PBE models. For this purpose, the coefficient was measured for the O2 transfer via headspace aeration to the bulk medium (air-liquid interface, ). The slope () was 4.87±0.13 hr−1 yielding a of 16.57±0.43 cm/hr for our cultivation system. Simulations were run for spinner flask cultures of mESCs at 60 rpm and a seeding density of 5×104 cells/ml. The temporal size distribution of mESC aggregates is shown in . During the first day, there was a significant increase in mean cluster size mainly due to aggregation. By day 3, mESC aggregates had a mean diameter of around 200 µm in agreement with our cell culture data. The total number of aggregates dropped substantially during the first day due to aggregation but remained relatively steady afterwards ( ). At day 4, the calculated total number of aggregates was around 2.7×104/50 ml of culture similar to our experiments and the corresponding biomass had increased 160-fold.

Figure 5

Model prediction of cluster size and O2 concentration distribution in a stirred suspension bioreactor of cultured mESC aggregates.

(A) Density plot of mESC aggregate size distribution during a 4-day culture. (B) Total aggregate number and biomass of mESC aggregates over 4 days of culture. (C) Oxygen concentration in the medium bulk and fraction of ESCs experiencing pO2 below 30 mm Hg. The model was run for a 50-ml working volume and 60 rpm with mESCs seeded at 5×104 cells/ml. Specific model parameters are from Tables 1–2. (D) Distribution of residence time of ESCs exposed to ≤30 mm Hg pO2 during a 4-day culture.

The fractions of cells exposed to various levels of O2 in spinner flask cultures were also predicted. Here, the fractions were computed of cells experiencing O2 tension below 30 mm Hg, which may be favorable for maintenance of ESC pluripotency ( ). We also calculated the fractions for extremely low pO2 (<10 mm Hg), which hampers ESC self-renewal and induces death [5]. The pO2 in the medium decreased to 82 mm Hg after 4 days with 58% of cells experiencing O2 less than 30 mm Hg. The corresponding fraction of ESCs at 10 mm Hg or less was 32%. If the O2 concentration in the culture bulk was fixed at 148 mm Hg throughout the cultivation period, only 33% of cells would be at pO2<30 mm Hg by day 4 (data not shown).

In addition, the residence time distribution ( ) shows that mESCs entered low O2 level (<30 mm Hg) at different time points. This variation is attributed to the distributed size of cell clusters and random aggregation under stirring. The data clearly illustrate that although a group of cells is exposed to a particular O2 range (e.g. below 30 mm Hg), not all cells within the group experience the O2 level for the same amount of time. Therefore, information about the fraction of cells at different pO2 levels should be coupled to residence times provided by the model presented here. Our findings show that stem cell aggregate size heterogeneity is linked to population heterogeneity with respect to the duration of exposure at different levels of O2 concentration despite the fairly homogeneous environment induced by mixing in SSBs.

Discussion

The transport of factors among cells within 3D structures such as aggregates is an important aspect of the performance of cultivation systems. The importance becomes even more pronounced for stem cell culture modalities given the central role of O2 in pluripotency and lineage specification. We employed a transient reaction-diffusion model to predict the O2 profile inside ESC aggregates of varying sizes. The effective diffusivity and consumption rate parameters of O2 in hESC and mESC aggregates were determined for the first time based on experiments without a steady state constraint. A PBE model was further developed for the SSB culture of ESC aggregates. The combined diffusion-reaction/PBE framework was utilized to predict the aggregation outcome and the O2 profile inside ESC clusters and the medium bulk in stirred-suspension over time. As a result, the fractions of cells exposed to particular levels of O2 and pertinent residence times were determined allowing to draw conclusions about the ensuing population heterogeneity and state of the cells (e.g. viability, differentiation potential).

The oxygen level to which ESCs are exposed to significantly affects cell fate decisions. Hypoxia (∼30 mm Hg) does not modulate adversely the proliferation rate of hESCs, which exhibit reduced spontaneous differentiation and chromosomal abnormalities [5], [6], [28]. Human ESCs and iPSCs exhibit differentiation efficiencies along particular lineages dependent on O2 availability. For instance, low pO2 enhances hESC differentiation to chondrogenic, endothelial and cardiac cells [7], [8], [29]. Mouse ESC EBs generate significantly more myeloid and erythroid progenitors when cultured at 3% O2 instead of 21% O2 [30]. However, extremely low pO2 can be detrimental for cells. Thus, delicate pO2 control is essential in culture systems to avoid hampering stem cell growth, viability and self-renewal or commitment. This is particularly relevant in scaling up stem cell cultivation where local concentration gradients and transfer of factors are more pronounced. The model presented in this study permits the prediction of the actual O2 level facilitating informed decisions about dynamic pO2 control and optimal aggregate size in systems for stem cell expansion.

For the results reported here, an implicit assumption was made that the O2 consumption rate (i.e. KM and Vmax) does not vary with time or the location of cells within each cluster. Although this is a premise that is common with previous reports, we cannot exclude the possibility of ESC metabolism adaptation to available O2 levels. For example, Vmax was estimated to vary with pO2 albeit within the same order of magnitude [31]. Human ESCs reportedly lower their O2 uptake rate when the dissolved O2 decreases from 21% to at 5% and 1% O2 [32]. Such adjustment may aid the cells to maintain the concentration of available O2 above zero at the core of aggregates with a diameter >200 µm thereby reducing the likelihood for necrosis ( in the absence of the constraint of 8 mm Hg critical pO2). However, further studies are warranted on the adaptation of O2 consumption by pluripotent stem cells in response to microenvironment O2 levels. With more information becoming available, our model can accommodate spatially (e.g. radial distance from the aggregate surface) and temporally changing O2 consumption rates.

Our study also raises the issue of ‘residence time’ (duration of exposure) of cells especially to low O2 levels. To our knowledge, this has not been investigated previously and even in reports of stem cell aggregate cultivation in stirred-suspension vessels, emphasis is placed on the actual pO2 the cells experience. It is important to note that even after the first 5 minutes of mESC culture in 2.2% O2, increases are observed in HIF-1α expression and reactive O2 species [33]. At 12–24 hours of hypoxia, cell-cycle regulatory proteins such as cyclins D1 and E, CDK2 and CDK4 along with retinoblastoma phosphorylation increase yielding a larger fraction of mESCs in the S phase and higher overall cell numbers. Accelerated cell proliferation in aggregate regions of hypoxia will lead to increased O2 consumption further depressing pO2 levels. As our model indicates, an increasing fraction of cells will reside in the hypoxic region of the aggregates which grow continuously. The ‘cut-off’ time of cell exposure to low pO2 before hypoxia-induced changes (e.g. cell differentiation, apoptosis etc.) become irreversible is unknown. Such information will be essential in formulating strategies (e.g. by modulating agitation ( )) to achieve ESC aggregate sizes so that a high cell fraction remains pluripotent for downstream differentiation to a desired phenotype.

The effects of O2 on gene expression are mediated largely by HIF transcription factors and as noted, mESCs upregulate HIF-1α within 5 minutes of exposure to 2.2% O2. Cells in our cultures did not exhibit significant differences in HIF-1α expression under different agitation rates (). This may be attributed to the average cluster diameter, which generally ranged between 200–300 µm ( ), leading to rather limited fractions of cells exposed to 2.2% O2 or less. Yet, we acknowledge that HIF-1α transcripts (probed by qPCR) may not be representative of the corresponding protein amount and activity levels. Moreover, the qPCR results are population ‘averages’ and HIF expression variation at locales with different pO2 cannot be discerned. Lastly, changes due to fluctuating O2 can be brought about by HIF-independent mechanisms such as the environmental sensing mammalian target of rapamycin (mTOR) [34].

Although this study focused on the transport of O2 in ESC aggregates in dish and stirred-suspension cultures, the concurrent availability and consumption of other components such as nutrients and factors for self-renewal or differentiation should be considered. A steady state analysis of the diffusion of glucose and (generic) cytokines (in addition to O2) in human EBs was recently presented [12]. The diffusivities of glucose and cytokine were lower than that of O2 suggesting that their transport places further limits on ESC aggregates to avoid cell starvation and death. Such limitations on transport may actually hinder the growth of aggregates and cause their size to stabilize. Interestingly enough, Cameron et al. [35] reported that over 21 days hESC aggregates displayed a maximum size at day 10 ranging between 400 and 500 µm without further increase. These observations support the choice of the Gompertz equation for modeling aggregate growth due to cell proliferation. The Gompertz equation in our study suggested a maximum radius larger than our experimental observations. This discrepancy may be explained by the fact that experimental data were obtained from static cultures with relatively constant bulk concentration of substrates (oxygen, glucose, etc.) due to daily medium changes whereas realistically these concentrations should decrease markedly over time. The maximum radius attained by ESC aggregates may be lower if substrate-dependent growth kinetics are considered. Additional mass balance equations for nutrients and growth factors can be incorporated in the model presented here to provide a more accurate representation of the bioreactor stem cell culture.

The diffusion and consumption of O2 in mESC and hESC cultures has been investigated in various studies [12], [13], [31], [36], [37] although only in a subset of these reports parameters were determined from experiments. The values reported for Vmax and KM are the same order of magnitude (10−8 mol.cm−3.s−1 and 10−8 mol.cm−3, respectively) as those presented here. The O2 consumption rate for mESCs was previously [38] measured at 4×10−17 (mol/cell.s), i.e. close to our observed rate assuming cells are spherical with a 15 µm diameter. Values for the diffusion coefficient based on experimentally determined ultrastructural ESC characteristics of aggregates such as porosity and tortuosity were not available to date in the literature for comparison. In fact, the ε and τ values did not vary significantly with the culture modality, i.e. dish or stirred suspension and the agitation speeds tested. This could suggest that cells adopt particular arrangements in the aggregate interior without significant influence from external mechanotransduction but this warrants further investigation. Such organizations may be mediated by cytoskeletal actions, which are central to the assembly of other cell types into aggregates [39].

It should be noted that no overt changes were observed in the values of Vmax, KM or Dt at different agitation rates. In part this can be explained by the fact that the majority of cells are within aggregates not experiencing the different levels of shear induced by the stirring rates tested here. Moreover, the DO measurements were performed in PBS-based solution with a known Bunsen solubility coefficient. This facilitates the comparison with findings from previous studies and eliminates potential discrepancies due to differences in culture media utilized in those reports. Even though we cannot rule out potential changes in parameter values if the measurements were performed in culture medium, these changes are expected to be minimal as shown for other systems [36], [37], [40].

A multiscale approach was described here with the coupling of the transient diffusion-reaction model for individual clusters with the PBE for the temporal evolution of size for ESC aggregate population in the bioreactor due to cell proliferation and collisions. The kernel used in this study comprises two parts: (i) The sum of cubic roots of volumes of the aggregates colliding raised to the power of 7/3 is applicable for a shear field with non-linear velocity profile [25], [26] as in the case of bioreactors with impellers and baffles. This is in contrast to the Smoluchowski kernel [41] used for cell aggregation (e.g. of platelets [42]) assuming a linear shear flow. (ii) A multiplicative factor dependent on operating conditions including the agitation rate represents the coalescence frequency written as the product of the collision frequency and the coalescence efficiency. Use of the kernel assumes that the aggregates are not rigid spheres but deform during collision. More importantly, the coalescence efficiency (collisions leading to new aggregate formation) decreases with increasing aggregate size in line with our experiments. Because of these attributes, the kernel is deemed suitable for describing the aggregation of stem cells in stirred-suspension vessels despite its original derivation for droplet coalescence in liquid-liquid dispersions. Explicit functions of the kernel parameters dependent on agitation rates can be estimated from the data ( ), and additional operating variables may be included (e.g. cell seeding density which was kept constant for all experiments here). Since transient size distributions of ESC aggregates cultivated in spinner flasks are available, a proper approach to determining the kernel for the PBE model is to solve the inverse problem [43]. Work is underway for deriving kernel functions via solving the inverse problem and comparing the solution to the functional utilized in this study.

The approach described here is certainly not without limitations but amenable to improvements. As already mentioned, the transport and consumption or secretion of factors other than O2 were not considered. Fluctuations in the consumption of nutrients and/or accumulation of waste will almost certainly affect the microenvironment within aggregates. These can be considered by expanding the model, for instance, to include reaction-diffusion equations (see Equation 1) for additional species. Moreover, O2 consumption was assumed invariant with the position of each cell within aggregates. Pertinent variations can be introduced by considering different expressions for Vmax and KM based on spatial criteria or other conditions (e.g. concentration of particular nutrients). In the context of envisioned stem cell process designs, the model in its current or expanded form can be used in conjunction with investigating relevant bounds in culture variables.

Nonetheless, the framework presented here can be applied widely to address both basic questions of stem cell physiology and issues surrounding relevant bioprocesses utilizing SSBs. The model can be employed to generate maps of concentrations of O2 (and other factors) expediting research on the role of hypoxia on stem cell self-renewal or commitment along particular lineages. Gaining a better understanding of the effects of O2 on fate selection is critical in developing optimal strategies for stem cell differentiation. The design and control of stem cell bioprocesses will also benefit from quantitative multiscale approaches such as the one described here, thereby bringing us closer to the realization of the potential of stem cell therapeutics.

Steady-state O The distance in the horizontal axis is measured from the aggregate center. The O2 level in the medium bulk was kept fixed at 148 mm Hg for all conditions. The aggregate radius ranged from 30 µm to 500 µm. Values for Vmax, KM and were taken from Table 1.

(TIF)

Click here for additional data file.

HIF1α gene expression in (A) mESCs and (B) hESCs cultured as aggregates in dishes and spinner flasks.

(TIF)

Click here for additional data file.

Comparison of experimental data (dash-circle lines) and model (solid lines) results for different agitation rates at day 4 of culture: (A) 60 rpm, (B) 80 rpm, (C) 100 rpm.

(TIF)

Click here for additional data file.

Measurement of O (A) Oxygen concentration during the experiment. The dissolved O2 in the medium bulk was first depleted by gassing with nitrogen. After the depletion phase, air was allowed to fill the space above the medium and the O2 level was measured. (B) Linear regression of the recovery phase for O2 marked in (A). The standard deviation was calculated from at least three independent experiments. The O2 transfer coefficient was determined from the slope of the curve.

(TIF)

Click here for additional data file.

Supplemental information on cell culture, image analysis, modeling methods and measurements of bioreactor O

(DOCX)

Click here for additional data file.