Modelling Oxygen Transport, Microcarrier Aggregation, and Hydrodynamic Constraints in Stirred Bioreactors for Scalable Developmental Engineering

1. Introduction

Developmental engineering (DE) has emerged as a promising bottom-up strategy for generating functional tissues through the progressive assembly of modular tissues (MTs) [1,2]. In this approach, MTs are formed by culturing multiple cell types on microscale three-dimensional (3D) modular scaffolds (MSs), which are subsequently assembled in a controlled, layer-by-layer manner [3]. Unlike conventional top-down tissue engineering strategies, which rely on pre-fabricated macroscopic scaffolds, DE offers an inherently scalable route to tissue fabrication by alleviating mass transport limitations—particularly oxygen diffusion—that constrain tissue size, viability, and functional maturation [4].

Despite its conceptual advantages, the translation of DE to large-scale applications remains constrained by the lack of robust bioprocesses capable of producing functional MTs at scale. While suitable microscale MSs can be manufactured using established biomaterials fabrication techniques [5,6], the scalable culture of cells on these scaffolds to form viable, well-functioning MTs remains a significant challenge. In particular, achieving consistent MT size, morphology, and cellular performance under bioreactor conditions represents a critical bottleneck. Without reliable and scalable strategies for MT production, the subsequent assembly of larger, functional tissue constructs is fundamentally limited, restricting the practical implementation of DE for regenerative medicine and advanced cell-based therapies.

Microcarriers—typically spherical beads with diameters below 500 μm—have long been employed as microscale supporting matrices for adherent mammalian cell cultures [7]. By providing extensive surface area for cell attachment, microcarriers enable high-density cultures in wave and stirred-tank bioreactors and underpin many industrial bioprocesses [8,9]. Consequently, microcarrier-based technologies play a central role in the global biopharmaceutical industry, driven by advances in disease biology, cell culture technologies, and increasing therapeutic demand [10,11].

Given the strong parallels between microcarriers and microscale MSs used in DE—particularly in terms of size, material requirements, and suitability for cell attachment—microcarrier-based bioprocesses provide a relevant and well-established framework for MT production. However, a key challenge in adapting these systems lies in microcarrier aggregation during culture. In conventional mammalian cell culture, microcarrier aggregation is typically regarded as undesirable due to its adverse effects on mass transport and culture homogeneity [12,13]. Microcarrier aggregation increases diffusion distances for oxygen and nutrients, leading to spatial gradients that impair cellular proliferation, differentiation, and viability. Oxygen transport is especially limiting, as diffusion distances exceeding approximately 200 μm result in hypoxic conditions that markedly alter cellular behaviour [14,15]. Given that individual microcarriers commonly range from 80 to 500 μm in diameter, microcarrier aggregation readily produces structures exceeding this critical length scale, often leading to hypoxic cores and cell necrosis [16].

In the context of DE, however, controlled aggregation of MS-supported MTs may represent an opportunity rather than a limitation. Aggregated MTs with defined sizes and morphologies could serve as higher-order building blocks, enabling hierarchical tissue assembly while preserving the advantages of modularity. The central challenge is therefore not the formation of aggregates per se, but the lack of quantitative understanding of how aggregate geometry, packing configuration, and hydrodynamic conditions influence oxygen transport and, consequently, MT viability and function.

Although several strategies have been explored to mitigate excessive aggregation and oxygen limitation in microcarrier based cell cultures—such as optimising inoculation protocols, microcarrier concentration, and impeller configuration [12]—their application to MT production remains limited. While impeller design and agitation speed have been shown to significantly influence aggregate size and cell density [17], quantitative links between aggregate morphology and oxygen availability within MT aggregates remain poorly characterised [18].

This study addresses these gaps by systematically investigating the effects of microcarrier aggregation on oxygen transport and cell culture performance. A computational simulation framework was employed to model microcarrier cluster formation and quantify the influence of aggregate geometry on the minimum oxygen concentration (C_min). The simulations examined the effects of microcarrier diameter, packing configuration, and aggregate morphology on oxygen diffusion, while hydrodynamic factors—including impeller design, agitation speed, microcarrier material properties, and medium composition—were evaluated for their roles in aggregate formation, suspension behaviour, and cell viability.

By integrating oxygen transport modelling with hydrodynamic analysis, this study provides a mechanistic understanding of the interplay between oxygen limitation, shear-induced stress, and cell damage in microcarrier-based tissue cultures. Importantly, this work reframes microcarrier aggregation—traditionally viewed as a processing drawback—as a controllable design parameter for the scalable manufacture of MTs with tunable sizes and shapes. These findings establish a rational basis for bioprocess optimisation aimed at enabling large-scale MT production and their subsequent assembly into functional tissues using DE principles.

2. Methodology

2.1. Modeling Oxygen Diffusion in 3D Tissue Constructs

To simulate oxygen diffusion through cells or tissues cultured on individual or aggregated microcarriers, the parameters used for different cell types, microcarriers, and media are listed in Appendix A. The thickness of each cell was assumed to be 10 micrometres (t = 10 µm). The cell culture media were treated as Newtonian fluids, and the microcarriers were assumed to be non-porous, perfectly spherical particles. Microcarrier oxygen permeability was not considered in this research, as each simulation evaluated the most demanding diffusion pathway, defined as transport occurring entirely through oxygen-consuming cells or tissues. McMurtrey developed a diffusion model to describe nutrient concentration within tissue constructs [19]. Based on a modified form of Fick’s second law (Equation 1), analytical solutions were derived for both steady-state and transient diffusion in 3D constructs. Oxygen concentration (C) was modelled as a function of tissue thickness (t), diffusion coefficient (D), the cellular oxygen consumption rate (φ), and the overall tissue thickness (T). The term φ (Equation 2) is defined as the product of the single-cell oxygen consumption rate (sOCR) and the cell density within the tissue constructs (ρ_cells).

$${\displaystyle \frac{\partial C}{\partial t}}=\phi -D\ast {\displaystyle \frac{{\partial }^{2}C}{\partial T^2}}$$

(1)

$$\phi =sOCR\ast {\rho }_{cells}$$

(2)

Equation 3 represents the solution reported by McMurtrey for oxygen diffusion from the surface to the centre of a spherical tissue construct under quasi-steady-state conditions, where C(r) is the Oxygen concentration, r is the depth of the tissue, R is the maximum depth or thickness of the tissue, C₀ the oxygen concentration at the deepest point of the tissue (R). A one-dimensional (1D) diffusion model was considered sufficient, as multiple studies have demonstrated its accuracy for describing oxygen transport in bioprocesses [20,21,22].

$$C\left(r\right)={\displaystyle \frac{\phi }{6D}}\ast \left(r^2-R^2\right)+C_o$$

(3)

The model simulated the stationary phase of cell growth, as this phase contains the highest cell numbers and therefore the greatest likelihood of microcarrier aggregation [23]. During this phase, both the number of viable cells and the cell density were assumed to remain constant, with ρ_cell = 1.00×10¹² cells/m³ [24,25]. By assuming that all oxygen is consumed at the deepest point of the tissue (C₀ = 0), C₀ represents the minimum oxygen concentration (C_min) required to sustain all cells within a tissue of maximum depth R. The potential effects of hyperoxia [26,27] were not considered in this model. A full derivation of Equation 3 is provided in Appendix B.

2.2. Geometric Modelling of Microcarrier Aggregation and Diffusion Distance Estimation

Equation 3 shows that oxygen concentration has an exponential dependence on R, which was hypothesised to be determined by the morphology and size of aggregated microcarriers during cell cultures. To estimate R, the model represented aggregated microcarriers in a three-dimensional (3D) space. Assuming identical spherical microcarriers, the centre-to-centre distance between two tangential microcarriers is equal to the diameter of a single microcarrier (d_mc). This distance is also referred to as the Euclidean distance (E) and is calculated from the coordinates of the centres of any two spheres, i and j, as defined in Equation 4.

$$E=d_{mc}=\sqrt{{\left(x_i-x_j\right)}^2+{\left(y_i-y_j\right)}^2+{\left(z_i-z_j\right)}^2}$$

(4)

In this framework, R is defined as the distance from the centroid of an aggregate, composed of n microcarriers, to the surrounding culture medium. The centroid ($\overline{x}, \overline{y}, \overline{z}$), representing the geometric centre of the aggregate, corresponds to the location furthest from the oxygen source in the surrounding medium in all directions. Accordingly, for an aggregate consisting of n microcarriers with radius ($r_{mc}$), $R$ was calculated using Equation 5, assuming that cells form a confluent monolayer on the surface of each microcarrier.

$$\begin{gathered} If E\ne 0 and\underset{1\le i\le n}{\min }\left\{E\right\}>r_{mc} where j=centroid coordinates and \\ i=all microcarrier coordinates, \\ R\approx \underset{1\le i\le n}{\min }\left\{\left|\underset{1\le i\le n}{\max }{x}_i-\overline{x}\right|,\left|\underset{1\le i\le n}{\min }{x}_i-\overline{x}\right|,\left|\underset{1\le i\le n}{\max }{y}_i-\overline{y}\right|,\left|\underset{1\le i\le n}{\min }{y}_i-\overline{y}\right|,\left|\underset{1\le i\le n}{\max }{z}_i-\overline{z}\right|,\left|\underset{1\le i\le n}{\min }{z}_i-\overline{z}\right|\right\}+r_{mc}+t_{cell} \end{gathered}$$

(5)

This equation is only applicable when the centroid lies outside all microcarriers $\left(\underset{1\le i\le n}{\min }\left\{E\right\}>r_{mc}\right)$, as cells cannot grow within the microcarriers under the non-porous assumption. If the centroid lies within a microcarrier $\left(\underset{1\le i\le n}{\min }\left\{E\right\}<r_{mc}\right)$, the model relocates it to the nearest point on the microcarrier surface using a normalised vector (Equation 6). The resulting surface coordinates are then substituted for the centroid coordinates in Equation 5 to calculate the updated value of R.

$$\begin{gathered} If E\ne 0 and\underset{1\le i\le n}{\min }\left\{E\right\}<r_{mc} where j=centroid coordinates and \\ i=all microcarrier coordinates, \\ {x}_{surf}=\left({\displaystyle \frac{\overline{x}}{\underset{1\le i\le n}{\min }\left\{E\right\}}}\ast r_{mc}\right)+r_{mc}, y_{surf}=\left({\displaystyle \frac{\overline{y}}{\underset{1\le i\le n}{\min }\left\{E\right\}}}\ast r_{mc}\right)+r_{mc}, z_{surf}=\left({\displaystyle \frac{\overline{z}}{\underset{1\le i\le n}{\min }\left\{E\right\}}}\ast r_{mc}\right)+r_{mc} \end{gathered}$$

(6)

Equations 5 and 6 are valid when the centroid does not coincide with a microcarrier centre ($E\ne 0$), in which case the nearest centroid position can be extrapolated. When $E=0$, however, the model repositions the centroid along the $\pm x, \pm y$ and $\pm z$ axes and evaluates each case to determine the maximum value of $R$ (Figure 1.A). This procedure is implemented using Equation 7, where L denotes the number of layers within an aggregate, represented by the number of distinct z-coordinates.

$$\begin{gathered} If E=0 where j=centroid coordinates and i=all microcarrier coordinates, \\ R\approx \underset{1\le i\le n}{\max }\{\left|\underset{1\le i\le n}{\max }{x}_i-(x_{centre}+r_{mc})\right|,\left|\underset{1\le i\le n}{\min }{x}_i-(x_{centre}-r_{mc})\right|,\left|\underset{1\le i\le n}{\max }{y}_i-(y_{centre}+r_{mc})\right|, \\ \left|\underset{1\le i\le n}{\min }{y}_i-(y_{centre}-r_{mc})\right|,\left|\underset{1\le i\le n}{\max }{z}_i-(z_{centre}+r_{mc})\right|, \\ \left|\underset{1\le i\le n}{\min }{z}_i-(z_{centre}-r_{mc})\right|\}+t_{cell}+\left\{\begin{array}{c}{r}_{mc} if L\ge 2\\ 0 if L=1\end{array}\right. \end{gathered}$$

(7)

Excluding experimental observations, which appear largely random, there is limited precedent regarding the orientation and sequence of microcarrier aggregation beyond overall size distributions [17,18]. Consequently, two key assumptions were adopted: (1) denser microcarrier packing results in a greater diffusion distance to the centroid, and therefore (2) the worst-case aggregation sequence is the one that produces the densest packing. The first assumption restricted the geometric variations of microcarrier aggregates to the three densest known packing arrangements for identical, perfect spheres: hexagonal close packing (HCP), face-centred cubic (FCC), and body-centred cubic (BCC) packing [28,29]. The second assumption substantially reduced the number of possible aggregation sequences, enabling the simulations to construct aggregates either conventionally (layer-by-layer along the +z axis) or alternately (adding microcarriers on either side of a central layer along the ±z axis) to allow cell growth in each packing configuration.

Equation 8 was used to calculate the volume of individual aggregates (V_aggregate) based on microcarrier morphology and the packing fraction ($\eta$). The packing fraction $\eta$ is defined as the ratio of the volume occupied by solid particles (microcarriers) to the volume of the corresponding unit cell (microcarriers plus viable cells). BCC packing has a packing fraction of 0.68 (${\eta }_{\mathrm{BCC}}=0.68$), whereas both HCP and FCC packing are denser, each achieving a packing fraction of 0.74 (${\eta }_{\mathrm{HCP}}={\eta }_{\mathrm{FCC}}=0.74$) [28,29]. This calculation slightly overestimates aggregate volume because it assumes microcarriers are separated by twice the cell thickness ($E=2r_{\mathrm{mc}}+2t_{\mathrm{cells}}$). The overall aggregate density (${\rho }_{\mathrm{aggregate}}$) was then determined from the volume fractions of the microcarriers and attached cells, together with their respective material densities.

$${\displaystyle \frac{n\ast {\left(r_{mc}+t_{cells}\right)}^3}{\eta }}=V_{aggregate}$$

(8)

2.3. Modeling Hydrodynamics and Cell–Microcarrier Interactions in Stirred Bioreactors

It is well-established that increasing the revolutions per minute (RPM) of the impellers or using more aggressive impeller designs in stirred reactors leads to smaller, more favourable aggregations [17]. Similarly, microcarrier type and the hydrodynamic properties of the medium can significantly affect aggregation geometry and overall culture performance [30,31]. However, these factors also influence other culture conditions, including microcarrier suspension, cell damage, and mixing efficiency. To quantify these effects, the model employs experimental correlations and equations. The Zwietering correlation (Equation 9) is a widely accepted method for predicting solids suspension and provides the minimum stirring speed, $N_{min} \left(RPS\right)$, required to prevent any particle from remaining stationary at the bottom of the vessel for more than 2 seconds [32,33]. All symbols retain their standard meanings, except for the empirically derived constant S, which depends on impeller and reactor design, and the mass fraction of microcarriers in the medium, X. The diameter of the aggregate ($d_{aggregate}$) was obtained from each simulation as the largest aggregation for a given n, estimated using the maximum Feret’s diameter (m). Values for S were based on a flat-bottomed, cylindrical vessel with four baffles of width 1/10 the tank diameter and a liquid height equal to the tank diameter [34,35]. All relevant data are provided in Appendix C.

$$N_{min}=\left({\displaystyle \frac{S\ast v_{media}^{0.1}\ast d_{aggregate}^{0.2}\ast {\left[\frac{g\ast \left({\rho }_{aggregate}-{\rho }_{media}\right)}{{\rho }_{media}}\right]}^{0.45}\ast X^{0.13}}{d_{impeller}^{0.85}}}\right)$$

(9)

Based on the agitation speed, the power requirement in Watts (P) was calculated using Equation 10, where P₀ is the power number—an empirically derived constant that varies with impeller type and flow regime. For turbulent flow, P₀ is constant, with values listed in Appendix C. In the transient region, however, P₀ depends on both impeller type and geometry, with the full calculation detailed in Appendix D. An agitator efficiency ${\eta }_{effiency}$ of 80% was assumed for all simulations.

$$P={\displaystyle \frac{P_0\ast {\rho }_{media}\ast {\left(N\right)}^3\ast d_{impeller}^5}{{\eta }_{effiency}}}$$

(10)

Promoting turbulent flow enhances mass transfer and ensures homogeneity within the culture, improving nutrient distribution and facilitating waste removal. Maintaining homogeneity is critical for producing large scale cultures with higher specific growth rates over extended periods [36]. According to Nienow [37], the Reynolds number for Newtonian fluids in a stirred bioreactor is defined as:

$$Re={\displaystyle \frac{{\rho }_{media}\ast N\ast d_{impeller}^2}{{\mu }_{media}}}$$

(11)

where all symbols have their usual meanings. For turbulent flow in a bioreactor, Re > 2 × 10⁴ [37]. Equation 11 determines the flow regime around the impeller; in larger reactors, greater energy dissipation can create stagnant regions. However, the model assumes turbulence throughout the reactor, as Re > 2 × 10⁴ exceeds the threshold, and impeller geometry is scaled accordingly.

While turbulence improves mass transfer and homogeneity, it also generates kinetic energy that can damage cells. Kolmogorov’s universal equilibrium theory of local isotropic turbulence was used to calculate the smallest characteristic eddy size (${\lambda }_K$) as a function of the medium’s kinematic viscosity (${\upsilon }_{media}$) and the total energy dissipated per unit mass of fluid (${\epsilon }_T)$ [34]. Equation 12 provides this calculation, comparing ${\lambda }_K$ to the diameter of the selected microcarrier to assess cell damage risk. Values of ${\lambda }_K$below ${\displaystyle \frac{2}{3}}{d}_{mc}$ moderately reduce cell growth, while values below ${\displaystyle \frac{1}{2}}{d}_{mc}$ are detrimental to cell health and can disrupt the cell membrane [38,39].

$${\lambda }_K={\left({\displaystyle \frac{{\upsilon }_{media}^3}{{\epsilon }_T}}\right)}^{{\displaystyle \frac{1}{4}}} where, {\epsilon }_T={\displaystyle \frac{P}{{\rho }_{media}\ast V}}$$

(12)

3. Results and Discussions

3.1. Influence of Aggregate Geometry and Microcarrier Size on Oxygen Limitation

Equations 1–7 were applied to simulate the number of microcarriers (n) within an aggregate and evaluate the resulting minimum oxygen concentration (C_min) across different microcarrier types. As shown in Figure 1.B, the simulations demonstrate that C_min increases with increasing individual microcarrier diameter (d_mc) within aggregates. For example, Cultisphere G, S, and GL (d_mc = 255 μm) reached peak C_min values of ~0.24 mM at n = 17, whereas Cytodex 1 (d_mc = 190 μm) reached only ~0.14 mM at the same aggregation number. This outcome is consistent with Equation 3, as larger d_mc increases the aggregate radius (R), which is exponentially related to C_min, and aligns with experimental trends reported by Preissmann et al. [40] and Rafiq et al. [41].

Beyond simple size effects, the simulations demonstrate that $C_{\mathrm{m}\mathrm{i}\mathrm{n}}$is strongly governed by aggregate geometry and packing configuration. While larger microcarriers are generally associated with higher $C_{\mathrm{m}\mathrm{i}\mathrm{n}}$, the increase in $C_{\mathrm{m}\mathrm{i}\mathrm{n}}$does not scale proportionally with $n$, indicating a significant dependence on how the microcarriers are arranged within the aggregate. This effect was first observed in Figure 1.B for 2≤ n ≤ 3. When n = 1, the diffusion radius R was equal to t_cell due to the assumption of a confluent cell monolayer. In contrast, when n = 2, cell bridging occurred between adjacent microcarriers, displacing the centroid into the interstitial “valley” between the particles. Similarly, at n = 3, the densest packing configuration was achieved, with tangential microcarriers arranged within a single layer or plane. As a result, for any configuration forming a single layer, or for any orientation where 2≤ n ≤ 3$,$ oxygen must diffuse through a cell mass with an effective diffusion length of t_cell + r_mc. The addition of a fourth microcarrier enabled the formation of a multilayer tetrahedral structure in a hexagonal close-packed (HCP) configuration, resulting in an increase in $C_{\mathrm{m}\mathrm{i}\mathrm{n}}$. Figure 2.A presents the corresponding $R$ values for all packing configurations, plotted against their respective number of the microcarriers in each aggregate or the unit cell numbers. Similarly, the incorporation of a fifth microcarrier permitted the formation of a third layer, yielding a triangular bipyramidal structure (also consistent with an HCP configuration) and a further increase in $C_{\mathrm{m}\mathrm{i}\mathrm{n}}$. However, the model predicts a plateau in $C_{\mathrm{m}\mathrm{i}\mathrm{n}}$for $6\le n\le 10$(Figure 1.B), arising from the formation of a square bipyramidal structure (Figure 2.B). This plateau indicates that increasing $n$does not necessarily lead to higher $C_{\mathrm{m}\mathrm{i}\mathrm{n}}$when aggregates adopt the densest packing configurations. Importantly, this observation is directly relevant to the manufacture of modular tissues for developmental engineering, as it suggests that larger modular tissues can be formed through the controlled aggregation of modular scaffolds. In particular, increasing the number of modular scaffolds ($n$) within specific ranges (e.g. $6\le n\le 10$) does not necessarily reduce oxygen availability, since $C_{\mathrm{m}\mathrm{i}\mathrm{n}}$can remain effectively constant under these packing configurations.

Figure 1.B further illustrates variations in $C_{\mathrm{m}\mathrm{i}\mathrm{n}}$gradients for aggregates with $11\le n\le 17$across different microcarrier types. For Ti-doped phosphate glass (Ti-7) microcarriers, a plateau in $C_{\mathrm{m}\mathrm{i}\mathrm{n}}$is observed, indicating the emergence of a dominant aggregation geometry at $n=13$. This configuration corresponds to the kissing number, defined as the minimum number (12) of identical spheres that can simultaneously contact and surround a central sphere. In contrast, Cytodex 2 microcarriers—whose diameter is approximately 98.2% larger—exhibit an increase of approximately 28.9% in $C_{\mathrm{m}\mathrm{i}\mathrm{n}}$over the range $13\le n\le 17$. These results demonstrate that the microcarrier diameter ($d_{mc}$) influences $C_{\mathrm{m}\mathrm{i}\mathrm{n}}$both directly, as described by Equation 3, and indirectly by governing packing configurations that determine aggregate geometry. These findings establish a novel mechanistic relationship between microcarrier size, aggregate geometry, and oxygen limitation, demonstrating the feasibility of culturing larger modular tissues through the controlled aggregation of suitable modular scaffolds without an associated increase in hypoxic risk.

3.2. Effects of Microcarrier Properties and Aggregation on Suspension and Cell Damage

Figure 3.A illustrates the influence of microcarrier density on the minimum stirring speed required for suspension ($N_{\mathrm{m}\mathrm{i}\mathrm{n}}$), as predicted by the Zwietering correlation (Equation 9). Aggregates with lower packing fractions contained a higher proportion of microcarriers relative to cells, resulting in a greater effective aggregate density and explaining the observed reduction in $N_{\mathrm{m}\mathrm{i}\mathrm{n}}$for $6\le n\le 10$. In addition, increasing microcarrier density led to higher suspension requirements overall. For instance, Ti-doped phosphate glass microcarriers (${\rho }_{mc}=2.75\times {10}^3{\mathrm{kg} \mathrm{m}}^{-3}$) required stirring speeds of approximately 70–100 RPM to remain suspended, whereas other microcarriers typically required only 10–30 RPM. Among the latter, Hillex II-170 exhibited relatively higher $N_{\mathrm{m}\mathrm{i}\mathrm{n}}$values, attributable to its denser polystyrene core compared with biopolymer-based microcarriers such as dextran, cellulose, and gelatin. Consequently, biopolymer microcarriers are often preferred, as lower $N_{\mathrm{m}\mathrm{i}\mathrm{n}}$values are associated with reduced energy consumption and a lower risk of shear-induced cell damage. In addition, their biodegradability and inherent biocompatibility promote efficient cell attachment and facilitate enzymatic cell harvesting [42]. However, their natural origin can introduce batch-to-batch variability and susceptibility to premature degradation, which may result in inconsistent culture performance [18,43].

Figure 3.B illustrates the combined effects of microcarrier diameter ($d_{mc}$) and the number of microcarriers per aggregate ($n$) on the cell population per aggregate, represented by the calculated total cell volume (Equation 8). In contrast to Figure 1.B, which exhibits a plateau in $C_{\mathrm{m}\mathrm{i}\mathrm{n}}$between $6\le n\le 10$, Figure 3.B shows a continued increase in cell population with increasing $n$. This divergence indicates that larger cellular volumes can be accommodated within aggregates without a corresponding increase in the effective diffusion path length ($R$). These results therefore support the earlier conclusion that maximal cell populations can be achieved through controlled aggregation while maintaining oxygen transport constraints within acceptable limits.

Although cell population generally increases with increasing $n$, a pronounced reduction is observed at $n=11$. This occurs because the densest packing configuration transitions from body-centred cubic (BCC) to face-centred cubic (FCC) packing (Figure 2.A), which increases the packing fraction and consequently reduces the accessible microcarrier surface area available for cell attachment. It should be noted, however, that the model may overestimate cell growth, as it assumes that cells readily bridge gaps between tangential microcarriers or between microcarriers aligned along the same axis, as demonstrated in related experimental studies [44]. Moreover, although aggregates formed from smaller microcarriers ($d_{mc}$) display lower maximum cell populations as a result of their reduced absolute surface area for cell attachment (Figure 3.B), they possess higher surface area-to-volume ratios, which can facilitate more efficient cell culture on a mass basis. This trend is consistent with observations reported by Ng et al [45]. Nevertheless, the same study also demonstrated that larger microcarriers may support more than threefold higher cell-per-bead loading, an effect attributed in part to microcarrier porosity. Such internal structural features are not incorporated into the present model and may therefore contribute to discrepancies between predicted and experimentally observed cell densities.

Table 1 summarises the influence of different stirring speeds or energy inputs required to suspend various types of microcarriers in DMEM cell culture medium supplemented with 0% or 20% fetal bovine serum (FBS) in stirred tank reactors with an impeller-to-tank diameter ratio of 0.33 ($0.33d_{tank}$) on the potential damage to cells cultured on the suspended microcarriers, as evaluated using Equation 12. The results clearly indicate that increasing microcarrier diameter ($d_{mc}$) reduces the rotational speed (RPM) required to induce a risk of cell damage. For instance, at 94 RPM, Cytopore 1 and 2 ($d_{mc}\approx 240\mu \mathrm{m}$) fall within the high Kolmogorov risk regime, whereas Cytodex 2 ($d_{mc}\approx 167.5\mu \mathrm{m}$) remains within the low-risk category. This size-dependent susceptibility to shear-induced damage is consistent with previous experimental and theoretical studies [34,39,46]. Additionally, increasing the medium density (${\rho }_{\mathrm{media}}$) reduces the minimum stirring speed ($N_{\min }$) required for suspension, as smaller density differences between microcarriers and the medium lower the hydrodynamic forces. This also decreases the total energy dissipation (${\epsilon }_T$), thereby reducing the risk of cell damage. Table 1 illustrates this effect: for all microcarriers, the RPM at which cell damage occurs is higher in DMEM supplemented with 20% FBS (${\rho }_{\mathrm{media}\text{@}20\%}=1.02\times {10}^3{\mathrm{kg}/\mathrm{m}}^3$) compared to DMEM without FBS (${\rho }_{\mathrm{media}\text{@}0\%}=1.00\times {10}^3{\mathrm{kg}/\mathrm{m}}^3$). Liste-Calleja et al [31] similarly reported improved cell densities and viabilities for HEK-293 cultures in higher FBS concentrations, attributing this primarily to enhanced nutrient availability. However, the present model suggests an additional hydrodynamic contribution, where increased medium density mitigates shear-induced damage, further supporting improved culture outcomes.

3.3. Effects of Agitation and Impeller Design on Cell Damage Risk

The effect of energy input in stirred tank bioreactors on physical damage to the cells cultured on microcarriers was evaluated by converting agitation speed (RPM) into power dissipation and eddy characteristics, which were subsequently related to the Kolmogorov length scale (λₖ) (Equations 10–12). Figure 4.A illustrates the variation in λₖ as a function of impeller type, impeller diameter (d_impeller), and RPM. Increasing RPM leads to a reduction in λₖ, reflecting higher energy dissipation and the breakdown of larger eddies into smaller turbulent structures. At lower RPM values, a local decrease in λₖ is observed for each impeller configuration. This behavior arises because the Reynolds number surpasses the turbulent transition threshold, prompting a change in the impeller power number, as described by Nienow [37]. Smaller eddies generate steep velocity gradients at the microcarrier surfaces, thereby increasing the likelihood of cell damage [47]. While comparisons between eddy size and cell diameter are commonly used to assess shear-induced damage [46], in microcarrier-based cell culture systems λₖ should instead be evaluated relative to the microcarrier diameter (d_mc), as this represents the most conservative or “worst-case” hydrodynamic condition [34].

Figure 4.A further demonstrates the reduction in λₖ with increasing impeller diameter (d_impeller). For instance, a Rushton impeller with d_impeller = 0.25 d_tank generates eddies of approximately 86.8 μm at maximum RPM. In contrast, a Rushton impeller with d_impeller = 0.5 d_tank produces a comparable eddy size at only 47 RPM (${\lambda }_{K@150RPM}$ ≈ 87.14 μm), and at its maximum RPM generates eddies less than half that size. Although this trend is supported by Cherry and Papoutsakis [48] and Grein et al [46], both studies also report that cell damage may result from microcarrier–microcarrier and microcarrier–impeller collisions, the frequency of which increases with agitation speed. Consequently, relying solely on the Kolmogorov length scale as an indicator of cell damage oversimplifies the combined effects of RPM and impeller geometry. For example, increasing d_impeller enlarges the impeller surface area, thereby increasing the probability of microcarrier–impeller interactions and associated mechanical damage.

As summarised in Figure 4.A, simulations based on Equations 10–12 and Appendix D indicate that λₖ varies with impeller type. The Rushton impeller generated the smallest eddies (λₖ at 150 RPM = 36.50 μm), owing to its flat, vertical blades, which promote high turbulence intensity, elevated local shear stresses, and predominantly radial flow [49]. For this reason, Rushton impellers are widely applied in fermentation processes, where bacterial cells—characterised by small size and relatively robust cell walls—can better tolerate low λₖ conditions [34]. In comparison, the six-blade, 45° pitched-blade turbine operating in upward pumping mode produced the second smallest eddies (${\lambda }_{K@150RPM}=36.50\mu m$). The inclined blade configuration results in comparatively lower turbulence intensity and energy dissipation than the Rushton design [50]. Although the λₖ profiles for upward and downward pumping turbines of identical diameter are effectively the same, Ibrahim et al. [51] reported that downward pumping configurations are particularly advantageous in shear-sensitive microcarrier-based cell culture systems.

More aggressive impeller geometries, characterised by larger diameters and flatter blades, generally exhibit lower minimum agitation speeds (N_min) because the turbulent energy dissipation rate (ε_T) is higher. However, these designs also demand greater power input as a consequence of increased drag and hydrodynamic resistance [52,53]. Simulation results derived from Equations 9–11, summarised in Table 2, indicate that increasing d_impeller reduces both N_min and the agitation speed required to reach fully turbulent conditions (N_turbulence), owing to improved bulk circulation and higher power transmission to the fluid. Consistent with this observation, Ramírez et al. [53] reported that larger impellers generate elevated shear stresses, which in turn decrease the extent of cell aggregation, as also noted by Zhang et al [17].

3.4. Model Assumptions, Limitations, and Applicability

Although the model incorporates both theoretical formulations and empirically derived correlations, it demonstrates strong agreement with experimental observations. For instance, the Zwietering correlation has been reported to exhibit a strong confidence interval of 95% ±13% [32]. Similarly, deviations in the impeller power number are relatively minor, with Rushton impellers showing a reported relative error of approximately 3.4% [54]. Nevertheless, the model assumes a uniform microcarrier diameter, which may not reflect practical conditions. Manufacturing techniques such as emulsion solvent evaporation are known to produce broad particle size distributions [55]. This variability is particularly pronounced in bio-based microcarriers, including Cultispher G, S, and GL (d_mc ± 125 μm), compared with more uniform synthetic alternatives such as Hillex II-170 (d_mc ± 10 μm).

Because mammalian cells are highly sensitive to hydrodynamic stress, bioprocesses for their cultivation commonly employ surface aeration as a low-shear oxygen transfer strategy [56]. In this approach, oxygen diffuses from the headspace into the culture medium across the free liquid surface rather than being introduced through sparging or dispersed bubbles. The present model assumes that surface aeration alone is adequate to sustain cell growth; however, this assumption is often not valid at industrial scale [57]. Although the model does not provide quantitative predictions for gassed systems, their effects are extensively reported in the literature. Sparged or bubbled systems increase the gas–liquid interfacial area, thereby enhancing the volumetric mass transfer coefficient (k_La) [58]. Improved oxygen transfer supports higher cell densities and enables larger microcarrier aggregates to receive sufficient oxygen to their core. Gas sparging can also promote turbulence and assist in microcarrier suspension in low-viscosity media—an effect not accounted for in the Zwietering correlation [34]. However, excessive gas flow rates may generate cavities behind impeller blades, decreasing effective fluid–impeller contact and consequently requiring higher agitation speeds to maintain suspension and mixing. Furthermore, bubbled systems pose a risk of cell damage due to steep velocity gradients generated during bubble rise and rupture [59,60]. The use of more aggressive impeller geometries can mitigate some of these effects by reducing bubble size, lowering associated shear damage, and improving the surface-area-to-volume ratio, thereby achieving higher k_La values [61].

In each simulation, the calculated R value was used to represent the most conservative diffusion pathway—namely, a route assumed to remain in continuous contact with cells rather than passing through a microcarrier. Accordingly, the minimum oxygen concentration (C_min) was determined under the assumption that oxygen consumption occurred along the entire diffusion distance without interruption by non-consuming regions. This conservative framework enabled accurate estimation of R for microcarrier aggregates with spherical and cubic geometries (Figure 1.A). However, minor deviations arose for microcarrier aggregates with triangular or tetrahedral morphologies. The limitation becomes more evident in Figure 4.B, where a large diamond-shaped aggregate is shown. In this configuration, the true shortest diffusion pathway follows an angular trajectory rather than a straight line, meaning the model slightly overestimates R relative to the physical system. Despite this, a modest overestimation is advantageous in practice, as it yields a C_min value marginally above the actual cellular oxygen requirement, thereby providing a conservative safety margin.

4. Conclusions

This study utilised theoretical simulations to characterise microcarrier aggregation morphologies in three-dimensional space, establishing a novel framework for analysing the impact of aggregation on oxygen transport. The findings indicate that the minimum oxygen concentration (C_min) is strongly governed by packing configuration. Within a defined range, this implies that larger, cell-populated tissue constructs may be generated through aggregates comprising a greater number of microcarriers (n), without proportionally increasing the effective diffusion path length (R). Through combined literature evaluation and supplementary simulations, the study also investigated the influence of mixing and agitation in stirred bioreactors on microcarrier aggregation dynamics. Elevated agitation speeds and more aggressive impeller designs—characterised by increased diameter and flatter blade geometry—were shown to improve mixing efficiency and lower the minimum suspension speed (N_min). However, these advantages are associated with a heightened risk of shear-induced cellular damage. Furthermore, media with higher density and viscosity were found to attenuate mechanical stress on cells, whereas increased microcarrier density necessitates higher agitation rates to maintain suspension, thereby intensifying hydrodynamic forces. Collectively, these results underscore the complex coupling between hydrodynamics, microcarrier aggregation behaviour, and oxygen transport. They further highlight the importance of predictive modelling in evaluating whether microcarrier aggregation—traditionally regarded as a limitation in cell culture—could instead be strategically harnessed as a modular tissue manufacturing approach within developmental engineering.

Although the study provides novel theoretical insight into microcarrier aggregation, the model has not yet been experimentally validated. Future work should therefore focus on verifying the predicted aggregate geometries and diffusion path lengths (R) using controlled stirred bioreactor experiments and quantitative imaging. Direct comparison between simulated and measured aggregate size distributions and oxygen profiles would strengthen model reliability. While aggregation was shown to decrease with increasing RPM, the quantitative relationship between agitation and aggregate stability remains undefined. Future research should incorporate force-balance analysis and, where possible, couple population balance modelling with hydrodynamic characterisation to establish numerical relationships. These advances would improve bioprocess design and enable controlled production of modular tissues by regulating aggregate size through modular tissue and operating conditions. These modular tissues could subsequently be assembled into larger, spatially organised constructs following developmental engineering principles, facilitating the fabrication of functional tissues for applications in regenerative medicine, disease modelling, drug screening, and potentially organ-scale tissue reconstruction.