Stem cell biomanufacturing under uncertainty: A case study in optimizing red blood cell production

As breakthrough cellular therapy discoveries are translated into reliable, commercializable applications, effective stem cell biomanufacturing requires systematically developing and optimizing bioprocess design and operation. This article proposes a rigorous computational framework for stem cell biomanufacturing under uncertainty. Our mathematical tool kit incorporates: high‐fidelity modeling, single variate and multivariate sensitivity analysis, global topological superstructure optimization, and robust optimization. The advantages of the proposed bioprocess optimization framework using, as a case study, a dual hollow fiber bioreactor producing red blood cells from progenitor cells were quantitatively demonstrated. The optimization phase reduces the cost by a factor of 4, and the price of insuring process performance against uncertainty is approximately 15% over the nominal optimal solution. Mathematical modeling and optimization can guide decision making; the possible commercial impact of this cellular therapy using the disruptive technology paradigm was quantitatively evaluated. © 2017 American Institute of Chemical Engineers AIChE J, 64: 3011–3022, 2018


Figure 11: Example of Computational Fluid Dynamics Visualization
We designed a 3D bioreactor representation in the computational fluid dynamics (CFD) software COMSOL; Figure 11 diagrams an example system. We tested many possible geometries of hollow fiber including the examples illustrated in Figure 12 with 6, 8 , and 10 hollow fibers; we also tried varying configurations of 3 PAN / 3 CRM, 4 PAN / 4 CRM, and 5 PAN / 5 CRM hollow fibers ( Figure 13). We designed our mesh with 9.37 ⇥ 10 3 elements (average quality: 0.238) and used a direct solver within COMSOL. The results in Figure 14 represent the bioreactor fraction which is glucose limited (LF GLU 2 [0, 1]) for the 30 days of the culture. Observe that more variability is induced via our inability to predict hollow fiber configuration a priori (Figure 14) than by the Krogh approximation ( Figure 9). : Simulations representing the glucose-limited bioreactor fraction over a 30-day culture. These simulations were run using the geometry specifications in Figure 12 and the hollow fiber placements in Figure 13.

Baseline Optimization Problem
Misener et al. 24 introduce the full mathematical model describing this system; the following model summarizes that work and highlights in red bold the uncertain parameter values whose value is fixed at a baseline here. Recall that this model approximates the dynamic system by discretizing the ODEs into 5 weeks, w 2 {1, ..., 5}.
Our optimal design is based on the Colijn and Mackey 31 hematopoiesis model 24 ; a discretized version is labeled Proliferation & Differentiation in the preceding model. But there are alternate models of hematopoiesis; we have no guarantee that the Colijn and Mackey 31 model is correct. The Lobato da Silva et al. 30 model for cellular growth, proliferation, and differentiation, illustrated in Figure 8, has the following form: where H h, w is the concentration of each cell type; [P h ] is the concentration of the parent cell type; g is a factor limiting expansion; k e h , k d h , and k k h represent expansion, differentiation, and death, respectively. Lobato da Silva et al. 30 further define a maturation map (illustrated in Figure 8) between the cell types: pluripotent stem cells (PSC) become stem cells (SC); SC become myeloid stem cells (MySC) or lymphoid stem cells (LySC); MySC become colonyforming unit -granulocyte, macrophage (CFU-GM), colony-forming unit -megakaryocyte (CFU-MEG), burst-forming unit -erythroid (BFU-E), or colony-forming unit -eosinophils (CFU-EO); BFU-E become colony-forming unit -erythroid (CFU-E). Figure 8 shows the path for driving stem cells towards becoming RBC precursor CFU-E 30 .
To test the alternate model of hematopoiesis, we substitute the equations for maximum cell density (H MAX ) and the reaction rate (V O 2 ) of O 2 with equivalent equations which are functions of the cell types in the Lobato da Silva et al. 30 model (PSC, SC, MySC, LySC, CFU-GM, CFU-MEG, BFU-E, CFU-EO); the new coefficients for the reaction rates (V O 2 , h ) and initialization count with respect to each different cell type are taken from Table 2. Additionally, the following equations for the Proliferation & Differentiation equations in the baseline optimization problem change:      This section considers the 30 uncertain parameters in Table 4; we define the single-and multi-variate sensitivity analysis. Parameters with known error bars vary within their expected uncertainty levels; the remaining parameters were allowed to take values 50% (L1), 90% (L2), 110% (U1), and 150% (U2) of their nominal values. Define the optimal point of the optimization algorithm to be f ⇤ (x, u) where variables x are the design decisions and u ⇢ R p are uncertain parameters. Eq. (2) introduces D f ⇤ i (x, u, d i,`) as a quantity similar to an elementary effect and normalizes it as the percent difference from the nominal optimal value; d i,`p erturbs uncertain parameter from its nominal u i value with respect to levels`2 {L1, L2, U1, U2}. We also define the sensitivity index S i in Eq. (3): Table 5 lists the uncertain parameters with sensitivity indices S i > 10%. We also define 2D parameter variation D
Calculating S i, j , we find that 41 combinations of 2 variables induce > 10% change beyond what we expect from the linear analysis and 12 combinations inducing > 15% change. For a linear optimization model we would expect the sensitivity S i, j to be approximately 0; significantly large values of S i, j indicate that nonlinearity in the optimization model is affecting the final outcome.

Robust Optimization Counterpart
This section first develops a robust counterpart model for oxygen consumption in the BR; we explain the implications in some level of detail. We also briefly develop the remaining robust counterparts for the uncertain parameters in Table 4; Figure 7 diagrams the ideas.
Oxygen Consumption: Oxygen consumption V O 2 is the sum of V O 2 , h for each type of hematopoietic cell h weighted by cellular concentration H h : Robust optimization only applies to inequalities, so the equality in Eq. (6) becomes ; the BR is O 2 limited and therefore a global optimization algorithm will drive the total oxygen consumption V O 2 as low as possible and just including the side is effectively the same as the equality because any BR not consuming the maximum amount of oxygen available would automatically improve itself by consuming more oxygen. We subsequently define the (e, d )-Interval Robust Counterpart 65 : Setting (e = 0, d = 0) in Eq. (7) corresponds to the nominal optimization model which does not incorporate uncertainty; a worst-case error analysis corresponds to setting e to its maximum value and d = 0. There are robust optimization methods which select levels of e and d to achieve a desired trade-off between a nominal and worst-case model, but we do not have enough data for the RBC BR to accurately parameterize e and d ; we therefore consider the nominal and worst-case models separately and in comparison to each other. Glucose Consumption V Glc, w / Lactate Production V Lac, w : Original equation in model formulation is based on experimental data: To transform into the robust counterpart, we begin by observing that the equalities can equivalently be the inequalities that make the equation as bad as possible: the baseline noticeably (by more than ⇡ 10%) when varying parameters, so we will not consider robustifying the design with respect to these parameters. Cell Flux Across Membrane: This parameter strongly affects final price, so we take the original equation: g P, w =Ĵ Cellŝ DH Cells · p · N HF, CRM · R 2, CRM · L Vol R 8 w 2 {1, ..., 5} and replace it with an equivalent inequality: g P, w Ĵ Cellŝ DH Cells · p · N HF, CRM · R 2, CRM · L Vol R 8 w 2 {1, ..., 5} This works because the optimization model is going to want to make g P, w as large as possible. Rearranging: g P, w Ĵ Cells J U Cellŝ DH Cells · p · N HF, CRM · R 2, CRM · L Vol R  d 8 w 2 {1, ..., 5}