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%) O 2 level without significantly reduced proliferation [5,6]. Lower O 2 tension (pO 2 ) 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 O 2 modulation of the activity of stem cell-fate controlling pathways such as the canonical Wnt/b-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 O 2 profile among stem or progenitor cells is important for cell fate prediction and control.
The distribution of O 2 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 O 2 through multiple interfaces complicates the analysis of the O 2 tension (pO 2 ) that each cell experiences. In most small-scale setups, O 2 is transferred via the gas/liquid interface (headspace aeration) to the culture medium under agitation. From the medium bulk O 2 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 O 2 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 O 2 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 pO 2 below a certain threshold. Differences in the ''residence time'' at a certain O 2 concentration among cells may be consequential for the differentiation and/or self-renewal of the stem cell population.
Here, the distribution of O 2 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 O 2 within the aggregates and parameters associated with the specific rate of O 2 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 O 2 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.

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% CO 2 /95% air at 37uC 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 mM, 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 (10 5 /ml) were seeded in 125-ml ProCulture spinner flasks (Corning, Corning, NY) with mTESR1 medium and 10 mM ROCKi. Agitation was kept constant during each experiment and medium was replaced daily thereafter without the addition of ROCKi.
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 a m = 1.27610 29 mol/ (cm 3 mm Hg) at 37uC [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 cm60.79 cm). The container was set on a submersible magnetic stirrer plate in a water bath at 37uC bath and agitation was set to 60 rpm. A dissolved O 2 (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 96well 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 e and tortuosity t 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-mm steps and the micrographs were processed in MATLAB (Mathworks, Natick, MA) to obtain porosity and tortuosity values ( Fig. 1A; see Materials S1).
Statistical Analysis. Data were expressed as mean6st.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 O 2 concentration (C) profile in each aggregate: where r is the aggregate radius, t is the time and D Ã t is the effective diffusion coefficient. Taking into account the porosity e and tortuosity t, D Ã t is calculated as where D t is the diffusion of O 2 in the tissue (aggregate). The consumption of O 2 (V O2 in mol/(cm 3 . s)) was described by Michaelis-Menten (M-M)-type kinetics: where V max is the maximum O 2 consumption rate and K M is the M-M constant. Because the system comprises aggregates of multiple sizes, cell clusters were classified into L 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 Materials S1) subject to appropriate initial and boundary conditions.
LC Lr r~0~0 symmetry about the aggregate center ð Þ ð 5Þ where R l is the radius for the aggregates in the l th bin. Calculation of the value of the O 2 transfer coefficient k l from the bulk liquid to the aggregates was based on the Frössling equation [19] relating the particle Sherwood (Sh) and Schmidt (Sc) numbers shown in the Materials S1 section. The total rate of O 2 transfer from the bulk to all aggregates (L groups) is: where f l and z total are the fraction of aggregates with average diameter 2R l 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 O 2 concentration C tz1 M D calc in the medium bulk at the (tz1)-time step in Equation 7.
Values for D Ã t , V max and K M were obtained by minimizing the objective function [20] comparing the calculated values C tz1 M D calc with the experimental ones C tz1 M D exp 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 n x,t ð Þ is defined such that n x,t ð Þdx is the fraction of aggregates with sizes between x and xzdx at time t per unit volume of the culture. The size x 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 x o and x C~x {x 0 . Aggregate breakage was not observed particularly after the first day of culture and therefore was not considered.
The rate of aggregate size change, Lx Lt , due to cell proliferation was modeled after the Gompertz equation [23]: where M is the aggregate mass reached as t?? and a G is a constant characteristic of the cell proliferation. Both parameters were estimated from temporal data of aggregates in static cultures. Different aggregation kernels K x x 0 j ð Þ 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 k, k 1 and a 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 mm) to allow for an extended range of aggregate sizes to be considered (see [17] and Materials S1).
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 O 2 profile in stirred-suspension cultures of ESC aggregates (Fig. 1). The total number of aggregates z total , size distribution of aggregates f l , oxygen concentration in the medium bulk C M , working volume of bioreactor V M (assumed equal to that of the medium) and agitation rate N RPM were used for initialization of the integrated model (see Materials S1). Cell growth and aggregation were modeled (Eq. 10-11) and the consumption of O 2 by all the aggregates in the bioreactor was determined (Eq. 7). This required knowledge of the O 2 concentration distribution within individual aggregates (Eq. 1). Subsequently, the bulk O 2 level C M was updated based on the net effect of O 2 consumption and supply: The rate of O 2 supply through headspace aeration was: where C Ã is the O 2 partial pressure in the gas phase (95% air/5% CO 2 , 37uC, 760 mm Hg), C M is the oxygen partial pressure in bulk medium, k l Ã a is the mass transfer coefficient at the gas-liquid interface. Values of k l Ã a were determined experimentally (see Materials S1).

Porosity and tortuosity of ESC aggregates
The transport of O 2 , 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 descriptorsporosity (e) and tortuosity (t)-for mESC and hESC aggregates formed in dishes or spinner flasks at different stirring rates (Fig. 1A). The porosity of mESC aggregates in spinner flasks at 100 rpm was 0.22960.016 at 24 hours (Fig. 1B) and decreased slightly after 3 (0.21860.011) and 5 days (0.22060.011) of culture suggesting that aggregates were undergoing compaction. For 3day cultures the e values for 60 and 80 rpm were 0.21760.02 and 0.21060.025, respectively (Fig. 1C). When cultured in dishes, mESCs formed looser aggregates with e = 0.25960.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.25460.017 (Fig. 1D) whereas hESC aggregates from spinner flasks at 45 and 60 rpm displayed porosities of 0.25860.011 and 0.27060.007, respectively.
The tortuosity of mESC aggregates at 100 rpm was 1.42760.083 at 24 hours and did not change significantly after 3 (1.50460.093) and 5 days (1.54760.181) of culture (Fig. 1B). Moreover, clusters formed at 60 rpm (1.46060.035), 80 rpm (1.53460.086) and in dishes (1.45760.422) exhibited similar tortuosity. Human ESC aggregates were characterized by slightly higher average t values than mESC clusters (Fig. 1D). After 4 days in stirred suspension an average t of 1.61260.135 (45 rpm) or 1.55160.086 (60 rpm) was determined whereas the tortuosity of dish aggregates was 1.63860.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 O 2
A transient reaction-diffusion model was employed to determine the values of the diffusion coefficient D Ã t and consumption (V max , K M ) of O 2 in the ESC aggregates from experimentally measured O 2 concentration values in the medium. First, the mean size was determined of aggregates in dishes and stirred-suspension cultures at different agitation rates ( Fig. 2A). At the end of the culture, mESC aggregates at 100 rpm were smaller (mean diameter: 113.04615.64 mm, p,0.05) than those at 60 rpm (173.8623.9 mm) and 80 rpm (169.01616.13 mm). The mean diameter of hESC aggregates at 45 and 60 rpm at day 4 was 163. 8 and 192.45 mm, respectively and similar to those in static cultures (166.61612.4 mm) (Fig. 2B). For each sample, the number fraction f l distribution was determined (Fig. 2C).
Cells within aggregates consumed O 2 resulting in a gradual decrease of the DO level in the medium (Fig. 2D). Typically, DO measurements were performed until the O 2 concentration in the medium dropped below 10 mm Hg. Solution of the model coupled with experimental data on DO, yielded values for D Ã t ,V max and K M ( Table 1). It should be noted that parameter values were calculated for a minimum (critical) O 2 concentration at the aggregate center of either (i) 0 mm Hg, or (ii) 8 mm Hg [27]. Irrespective of the minimum O 2 level considered, the values of D t (see Eq. 2) were lower than that of the diffusion coefficient D O2 of O 2 in pure water at 37uC (2.78610 25 cm 2 /s). As a measure of maximum O 2 consumption rate, V max was higher for mESCs in suspension culture at 60 rpm than for hESCs (p,0.05). Mouse ESCs also displayed greater V max at 100 vs. 60 rpm suggesting an effect of the agitation on O 2 metabolism (p,0.05). The mean V max 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 D Ã t , V max and K M were determined, the reaction-diffusion model was implemented to calculate the O 2 profile for ESC aggregates of different size. The diffusion of O 2 becomes limited as the aggregate radius increases leading to reduced O 2 levels for cells residing closer to the core. To better illustrate this, the O 2 distribution was calculated in hESC aggregates sampled from spinner flask cultures assuming a constant O 2 concentration in the bulk of 148 mm Hg (saturation) (Fig. 3A). Because hypoxia studies typically employ pO 2 of ,30 mm Hg (4-5% O 2 ) [5], we calculated the fractions of ESCs exposed to O 2 below this level for different aggregate sizes. Human ESC aggregates smaller than 200 mm were normoxic throughout their volume. In contrast, almost 23% of the hESCs in aggregates larger than 400 mm were under hypoxia. This fraction increased to 70% for aggregates with 1,000 mm diameter (Fig. 3A).
The O 2 profiles in hESC aggregates of various sizes maintained in suspension (60 rpm) or static culture are shown in Fig. 3B. Similarly, the O 2 distributions within mESC aggregates cultured at 60 rpm or statically are also shown (Fig. 3C). It should be noted that the results shown in Figs. 3B-C were obtained assuming that  all aggregates in culture are the same size. Simulation results for additional agitation rates are shown in Figure S1.
In general, the fractions of cells experiencing an O 2 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 (Figs. 3B-C). This is most likely due to the faster transfer (k l ) of O 2 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 mm and only in static cultures. The higher hindrance for O 2 (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 O 2 profiles across all the agitation rates we tested (Fig. S1). As the aggregate size increased, cells positioned over 100 mm from the aggregate surface experienced less than 10 mm Hg of O 2 . Compared to hESCs, the fractions of mESCs under the same O 2 concentration threshold were greater for equal size aggregates owing most likely to the higher V max of the latter. The difference was most significant for aggregates with a radius between 100-250 mm.
Nevertheless, the expansion of mESCs and hESCs in spinner flasks gave rise to aggregate diameters mostly below 200-300 mm (Fig. 2). Thus, if a normoxic concentration of O 2 can be maintained in the bulk only 2.2% of the cells (located in the aggregate center) experience O 2 tension below 30 mm Hg (Fig. 3A). In fact, when we analyzed cells by qPCR for HIF1a expression (Fig. S2) there was no notable differences between static and stirred-suspension cultures. Also, compared to monolayer cultures at normoxia, HIF1a 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 O 2 profiles presented above, the O 2 concentration in the medium was assumed to be fixed at saturation (,148 mm Hg). This assumption is valid if O 2 can be supplied efficiently and continuously to the cultivation system or the culture is sparse so that O 2 consumption is trivial. However, stirred-suspension cultures accommodate high cell densities resulting in fluctuating O 2 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 O 2 . Given the multitude of effects of O 2 on stem cell physiology, predicting the time-variant distribution of O 2 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.71610 6 , a G = 5.72610 23 hr 21 ) based on the recorded size changes (Fig. 4A). The Gompertz model accurately described the growth of mESC aggregates (R 2 = 0.99, n = 20). The maximum attainable size for mESC aggregates was calculated at 1600 mm (lnr max~5 :4).
Parameters for the aggregation kernel were determined from size data of mESC aggregates cultured in spinner flasks (Fig. 4B). It should be noted that the kernel parameters were calculated from data of cell cultures at 0-72 hours ( Table 2). 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 mm on day 2 (mainly due to collision-driven aggregation) in line with the model predictions (Fig. 4C). Subsequent increase in aggregate size was less pronounced as the mean radius was approximately 150 mm (day 4). The corresponding average size at higher agitation rates was lower (e.g. ,90 mm at 100 rpm, Fig. S3). 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%).
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

Prediction of O 2 profile and hypoxia fraction in ESC aggregate stirred-suspension culture
With the PBE model of ESC aggregation in the SSB in place, the O 2 distribution in the medium and within ESC aggregates was computed by coupling the diffusion reaction and PBE models. For this purpose, the coefficient k l Ã was measured for the O 2 transfer via headspace aeration to the bulk medium (air-liquid interface, Fig. S4). The slope (k l Ã : a=V M ) was 4.8760.13 hr 21 yielding a k l Ã of 16.5760.43 cm/hr for our cultivation system. Simulations were run for spinner flask cultures of mESCs at 60 rpm and a seeding density of 5610 4 cells/ml. The temporal size distribution of mESC aggregates is shown in Fig. 5A. 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 mm 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 (Fig. 5B). At day 4, the calculated total number of aggregates was around 2.7610 4 /50 ml of culture similar to our experiments and the corresponding biomass had increased 160-fold. The fractions of cells exposed to various levels of O 2 in spinner flask cultures were also predicted. Here, the fractions were computed of cells experiencing O 2 tension below 30 mm Hg, which may be favorable for maintenance of ESC pluripotency (Fig. 5C). We also calculated the fractions for extremely low pO 2 (,10 mm Hg), which hampers ESC self-renewal and induces death [5]. The pO 2 in the medium decreased to 82 mm Hg after 4  In addition, the residence time distribution (Fig. 5D) shows that mESCs entered low O 2 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 O 2 range (e.g. below 30 mm Hg), not all cells within the group experience the O 2 level for the same amount of time. Therefore, information about the fraction of cells at different pO 2 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 O 2 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 O 2 in pluripotency and lineage specification. We employed a transient reaction-diffusion model to predict the O 2 profile inside ESC aggregates of varying sizes. The effective diffusivity and consumption rate parameters of O 2 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 O 2 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 O 2 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 O 2 availability. For instance, low pO 2 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% O 2 instead of 21% O 2 [30]. However, extremely low pO 2 can be detrimental for cells. Thus, delicate pO 2 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 O 2 level facilitating informed decisions about dynamic pO 2 control and optimal aggregate size in systems for stem cell expansion. For the results reported here, an implicit assumption was made that the O 2 consumption rate (i.e. K M and V max ) 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 O 2 levels. For example, V max was estimated to vary with pO 2 albeit within the same order of magnitude [31]. Human ESCs reportedly lower their O 2 uptake rate when the dissolved O 2 decreases from 21% to at 5% and 1% O 2 [32]. Such adjustment may aid the cells to maintain the concentration of available O 2 above zero at the core of aggregates with a diameter .200 mm thereby reducing the likelihood for necrosis (Fig. S1 in the absence of the constraint of 8 mm Hg critical pO 2 ). However, further studies are warranted on the adaptation of O 2 consumption by pluripotent stem cells in response to microenvironment O 2 levels. With more information becoming available, our model can accommodate spatially (e.g. radial distance from the aggregate surface) and temporally changing O 2 consumption rates.
Our study also raises the issue of 'residence time' (duration of exposure) of cells especially to low O 2 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 pO 2 the cells experience. It is important to note that even after the first 5 minutes of mESC culture in 2.2% O 2 , increases are observed in HIF-1a expression and reactive O 2 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 O 2 consumption further depressing pO 2 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 pO 2 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 ( Table 2)) to achieve ESC aggregate sizes so that a high cell fraction remains pluripotent for downstream differentiation to a desired phenotype.
The effects of O 2 on gene expression are mediated largely by HIF transcription factors and as noted, mESCs upregulate HIF-1a within 5 minutes of exposure to 2.2% O 2 . Cells in our cultures did not exhibit significant differences in HIF-1a expression under different agitation rates (Fig. S2). This may be attributed to the average cluster diameter, which generally ranged between 200-300 mm (Fig. 3), leading to rather limited fractions of cells exposed to 2.2% O 2 or less. Yet, we acknowledge that HIF-1a 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 pO 2 cannot be discerned. Lastly, changes due to fluctuating O 2 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 O 2 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 O 2 ) in human EBs was recently presented [12]. The diffusivities of glucose and cytokine were lower than that of O 2 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 mm 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 O 2 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 V max and K M are the same order of magnitude (10 28 mol . cm 23 . s 21 and 10 28 mol . cm 23 , respectively) as those presented here. The O 2 consumption rate for mESCs was previously [38] measured at 4610 217 (mol/cell . s), i.e. close to our observed rate assuming cells are spherical with a 15 mm 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 e and t 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 V max , K M or D t 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 ( Table 2), 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 O 2 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, O 2 consumption was assumed invariant with the position of each cell within aggregates. Pertinent variations can be introduced by considering different expressions for V max and K M 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 O 2 (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 O 2 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. Figure S1 Steady-state O 2 profile for ESC aggregates of different sizes cultured under different conditions. Oxygen concentration profile for (A) hESC aggregates cultured in static or stirred suspension (45 rpm and 60 rpm) culture, and (B) mESC aggregates cultured in static or stirred suspension (60 rpm, 80 rpm and 100 rpm) culture. The distance in the horizontal axis is measured from the aggregate center. The O 2 level in the medium bulk was kept fixed at 148 mm Hg for all conditions. The aggregate radius ranged from 30 mm to 500 mm. Values for V max , K M and D Ã t were taken from Table 1