# **Bioprocess Systems Engineering Applications in Pharmaceutical Manufacturing**

Edited by Ralf Pörtner and Johannes Möller

Printed Edition of the Special Issue Published in *Processes*

www.mdpi.com/journal/processes

# **Bioprocess Systems Engineering Applications in Pharmaceutical Manufacturing**

# **Bioprocess Systems Engineering Applications in Pharmaceutical Manufacturing**

Editors

**Ralf P ¨ortner Johannes M ¨oller**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Ralf Portner ¨ Bioprocess and Biosystems Engineering Hamburg University of Technology Hamburg Germany

Johannes Moller ¨ Bioprocess and Biosystems Engineering Hamburg University of Technology Hamburg Germany

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Processes* (ISSN 2227-9717) (available at: www.mdpi.com/journal/processes/specialissues/ PharmaceuticalManufacturing).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-5210-1 (Hbk) ISBN 9978-3-0365-5209-5 (PDF)**

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

# **Contents**

# *Editorial* **Special Issue "Bioprocess Systems Engineering Applications in Pharmaceutical Manufacturing"**

**Ralf Pörtner \* and Johannes Möller \***

Institute of Bioprocess and Biosystems Engineering, Hamburg University of Technology, 21073 Hamburg, Germany

**\*** Correspondence<EMAIL_ADDRESS>(R.P<EMAIL_ADDRESS>(J.M.)

Biopharmaceutical and pharmaceutical manufacturing are strongly influenced by the process analytical technology initiative (PAT) and quality by design (QbD) methodologies, which are designed to enhance the understanding of more integrated processes. The major aim of this effort can be summarized as developing a mechanistic understanding of a wide range of process steps, including the development of technologies to perform online measurements and real-time control and optimization. Furthermore, minimization of the number of empirical experiments and the model-assisted exploration of the process design space are targeted. Even if tremendous progress has been achieved so far, there is still work to be carried out in order to realize the full potential of the process systems engineering toolbox.

Within this Special Issue of *Processes*, an overview of cutting-edge developments of process systems engineering for biopharmaceutical and pharmaceutical manufacturing processes is given, including model-based process design, Digital Twins, computer-aided process understanding, process development and optimization, and monitoring and control of bioprocesses. The biopharmaceutical processes addressed focus on the manufacturing of biopharmaceuticals, mainly by Chinese hamster ovary (CHO) cells, as well as adenoassociated virus production and generation of cell spheroids for cell therapies.

*Both model-based process designs and the Digital Twin concept* have gained increasing interest for the development and optimization of biopharmaceutical production processes. Such methods are still not state-of-the-art for cell culture processes during development or manufacturing, although first approaches have been proposed. This highlights a need for improved methods and tools for optimal experimental design, optimal and robust process design, and process optimization for the purposes of monitoring and control during manufacturing. Three contributions within this Special Issue address this topic, which are highlighted as follows:

Bayer et al. [1] present a digital bioprocess twin used for a model-based design of experiment (DoE) to accelerate the design space exploration and thereby decrease the time needed to identify the optimum combination of critical process parameters (CPP) for the variables of interest. This Digital Twin simultaneously delivers additional process understanding while accelerating bioprocess development and optimization by applying in silico simulations and only perform the recommended experiments. A structured workflow is presented using different initial data sets to reduce experimental efforts, evaluate the results, and additionally to investigate the applicability of an intensified DoE (iDoE) for such a model-based DoE. By this, the best CPP combinations in a design space are identified with the highest space-time yield.

Strategies for multi-objective optimization of industrial biopharmaceutical processes are addressed by Hernández Rodríguez et al. [2]. In industrial applications, it is typically desirable to optimize several conflicting objectives at a time, leading to suitable tradeoffs and compromises. However, multi-objective optimization is more complex and its application is still not state-of-the-art in the context of cell culture processes (probably

**Citation:** Pörtner, R.; Möller, J. Special Issue "Bioprocess Systems Engineering Applications in Pharmaceutical Manufacturing". *Processes* **2022**, *10*, 1634. https:// doi.org/10.3390/pr10081634

Received: 8 August 2022 Accepted: 12 August 2022 Published: 18 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

due to a lack of related studies and instructions). The contribution presents a conceptual workflow which couples uncertainty-based upstream simulation and Bayes optimization using Gaussian processes. Its application is demonstrated in a simulation case for the design of a robust cell culture expansion process (seed train), meaning that despite uncertainties and variabilities concerning cell growth, low variations of viable cell density during the seed train are obtained. This approach provides the potential to be used in the form of a decision tool (e.g., for the choice of an optimal and robust seed train design or for further optimization tasks within process development).

The selection of appropriate clones of production cells is essential for the optimized manufacturing of complex biopharmaceuticals using bioreactors. A simulation study on model-assisted clone selection for CHO-cells is presented by Hernández Rodríguez et al. The authors of [3] address the question, if clonal cell populations showing high cell-specific growth rates are more favourable than cell lines with higher cell-specific productivities (or vice versa). A mechanistic cell culture model was adapted to the experimental data of such clonally-derived cell population. Uncertainties and prior knowledge concerning model parameters were considered using Bayesian parameter estimations. This model was then used to define an inoculum train protocol. It could be shown that growth rates have a higher impact on overall process productivity and for product output per year, whereas cells with higher productivity can potentially generate higher product concentrations in the production vessel.

Besides the use of mechanistic process models, computational fluid dynamic studies can support bioprocess design and optimization. One paper in particular focuses on *computer-aided process understanding*. Bioreactor design and scale-up in today's biopharma industry relies mostly on empirical correlations, experience, and engineering heuristics, which can hardly provide the link between hydrodynamics within the bioreactor and biological process behaviour. Freiberger et al. [4] investigated cellular effects and locally resolved hydrodynamics in stirred bioreactors for impellers with different spatial hydrodynamics. Therefore, the hydrodynamics, mainly flow velocity, shear rate, and power input in a single- and a three-impeller bioreactor setup were analyzed by means of CFD simulations, and cultivation experiments with antibody-producing CHO cells were performed at various agitation rates in both reactor setups. It could be shown that behaviour of the cells in the different reactor set-ups cannot be linked to parameters commonly used to describe shear effects on cells such as the mean energy dissipation rate or the Kolmogorov length scale, even if this concept is extended by locally resolved hydrodynamic parameters. Alternatively, the hydrodynamic heterogeneity was statistically quantified by means of variance coefficients of the hydrodynamic parameters fluid velocity, shear rate, and energy dissipation rate. The calculated variance coefficients of all hydrodynamic parameters were higher in the setup with three impellers than in the single impeller setup, which might explain the rather stable process behaviour in multiple impeller systems due to the reduced hydrodynamic heterogeneity.

Three contributions are dedicated to aspects of *process development and optimization*. Müller et al. [5] studied seed train intensification using an ultra-high cell density cell banking process. A frequently used approach for seed train intensification uses N − 1 perfusion, in which perfusion cultivation is carried out as the final step of inoculum production to generate ultra-high cell banks exceeding 100 × 106 cells mL−1. These cells can subsequently be used to inoculate a production bioreactor. On the one hand, the inoculum production steps can be reduced, and on the other hand a continuous process or a high-seed fed-batch process can be directly implemented with these cells instead of the otherwise usual low-seed fed-batch process. Within the study, an ultra-high cell density working cell bank was established for an immunoglobulin G-producing CHO cell line. A comparison with the standard approach shows that cell growth and antibody production are comparable, but time savings of greater than 35% are possible for inoculum production.

Ladd et al. [6] developed a process for continuous transfection for adeno-associated virus production in microcarrier-based culture. Adeno-associated virus vectors have great potential for gene therapy. However, a major bottleneck for this kind of therapy is the limitation of production capacity. Higher specific AAV vector yield is often reported for adherent cell systems compared to cells in suspension, and a microcarrier-based culture is well established for the culture of anchored cells on a larger scale. The purpose of the present study was to explore how microcarrier cultures could provide a solution for the production of adeno-associated virus vectors based on the triple plasmid transfection of HEK293T cells in a continuously operated stirred tank bioreactor. The present investigation provided a proof-of-concept of a continuous process based on microcarriers in a stirred-tank bioreactor.

Petry and Salzig [7] developed a large-scale production process of size-adjusted β-cell spheroids. The large and growing number of patients living with diabetes has generated interest in the promise of β-cell therapy to restore lost β-cell mass. For β-cell replacement therapies, one challenge is the manufacturing of a sufficient number of functional β-cells manufactured as 3D constructs, known as spheroids with a controlled size. For this, a process in a fully controlled stirred bioreactor systems was established using the INS-1 β-cell line as a model for process development. Specifically, the dynamic agglomeration of β-cells to determine minimal seeding densities, spheroid strength, and the influence of turbulent shear stress was investigated in order to generate spheroids of a defined size. The process developed in shaking flasks was successfully transferred to a stirred bioreactor, and it could be shown that functional β-cell spheroids sufficient for β-cell therapy applications can be generated.

Two contributions are dedicated to the aspects of *monitoring and control of bioprocesses*. Reyes et al. [8] provide a review on modern sensor tools and techniques for monitoring and control, addressing especial technological innovation directed towards online in situ continuous monitoring of quality attributes that could previously only be estimated offline. These new sensing technologies when coupled with software models have shown promise for unique fingerprinting, smart process control, outcome improvement, and prediction. State-of-the-art sensing technologies and their applications in the context of cell culture monitoring are reviewed with an emphasis on the coming push towards industry 4.0 and smart manufacturing within the biopharmaceutical sector. Additionally, perspectives concerning how this can be leveraged to improve both understanding and outcomes of cell culture processes are discussed.

A new and promising biosensor technology is introduced by Gaudreault et al. [9], a surface plasmon resonance (SPR)-based biosensor. SPR-based biosensors can play a role in enabling the development of improved bioprocess monitoring and control strategies. In this review, the applications of SPR that are or could be related to bioprocess monitoring in three spheres are examined such as biotherapeutics production monitoring, vaccine monitoring, and bacteria and contaminant detection. These applications mainly exploit SPR's ability to measure solution species concentrations, but performing kinetic analyses is also possible and could prove useful for product quality assessments. SPR-based biosensors exhibit potential as product monitoring tools from early production to the end of downstream processing, paving the way for more efficient production control. However, more work needs to be carried out to facilitate or eliminate the need for sample pre-processing and to optimize the experimental protocols.

We thank the authors, reviewers, and editors who have contributed to the success of this Special Issue.

**Author Contributions:** Conceptualization, J.M. and R.P.; Writing—original draft preparation, J.M. and R.P.; All authors have read and agreed to the published version of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

# *Article* **Digital Twin Application for Model-Based DoE to Rapidly Identify Ideal Process Conditions for Space-Time Yield Optimization**

**Benjamin Bayer 1,2,†, Roger Dalmau Diaz 1,2,†, Michael Melcher 1, Gerald Striedner 1,2 and Mark Duerkop 1,2,\***

**Abstract:** The fast exploration of a design space and identification of the best process conditions facilitating the highest space-time yield are of great interest for manufacturers. To obtain this information, depending on the design space, a large number of practical experiments must be performed, analyzed, and evaluated. To reduce this experimental effort and increase the process understanding, we evaluated a model-based design of experiments to rapidly identify the optimum process conditions in a design space maximizing space-time yield. From a small initial dataset, hybrid models were implemented and used as digital bioprocess twins, thus obtaining the recommended optimal experiment. In cases where these optimum conditions were not covered by existing data, the experiment was carried out and added to the initial data set, re-training the hybrid model. The procedure was repeated until the model gained certainty about the best process conditions, i.e., no new recommendations. To evaluate this workflow, we utilized different initial data sets and assessed their respective performances. The fastest approach for optimizing the space-time yield in a threedimensional design space was found with five initial experiments. The digital twin gained certainty after four recommendations, leading to a significantly reduced experimental effort compared to other state-of-the-art approaches. This highlights the benefits of in silico design space exploration for accelerating knowledge-based bioprocess development, and reducing the number of hands-on experiments, time, energy, and raw materials.

**Keywords:** *Escherichia coli*; hybrid modeling; machine learning; model-assisted DoE; quality by design; upstream bioprocessing

**1. Introduction**

For the production of biopharmaceuticals, it is of high importance to guarantee a specified product quality for patient safety. Raw materials, process deviations, and unrecognized faults may result in altered quality, and finally in batch rejection [1]. Process characterization in the biopharmaceutical industry has long been known and emphasized by the authorities, thus, processes must be closely monitored and well understood to ensure robust and uniform product quality. The most prominent guidance is the process analytical technology (PAT) guide by the US federal drug administration (FDA). Additionally, the quality by design (QbD) initiative [2] greatly emphasizes process understanding during the development of a bioprocess to guarantee a stable and uniform product quality output and fewer rejected batches [3]. To achieve these objectives, the statistical design of experiments (DoE) and advanced online monitoring are highlighted. The herein experimentally investigated design space is built by different combinations of critical process parameters (CPP) and critical material attributes (CMA), which affect the target parameters and the

**Citation:** Bayer, B.; Dalmau Diaz, R.; Melcher, M.; Striedner, G.; Duerkop, M. Digital Twin Application for Model-Based DoE to Rapidly Identify Ideal Process Conditions for Space-Time Yield Optimization. *Processes* **2021**, *9*, 1109. https:// doi.org/10.3390/pr9071109

Academic Editors: Ralf Pörtner and Johannes Möller

Received: 2 June 2021 Accepted: 23 June 2021 Published: 25 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

critical quality attributes (CQA) [4]. For such a design space exploration, different DoEs can be applied, e.g., full factorial, fractional factorial, Box–Behnken, Doehlert, and hypercubes, differing in the number of required experiments and the amount of information generated [5]. Besides the increased process understanding, for process optimization of the target molecule, it is still important to quickly find the best CPP combination in the design space, at which the production process will be performed, e.g., biomass, product titer, or space-time yield [6]. Such DoE studies are combined with process modeling to generate added value and further accelerate these tasks [7].

The most common techniques for bioprocess modeling are data-driven (black box) and mechanistic (white box) approaches, each with their own characteristics, advantages, and disadvantages [8]. Since the parameters in data-driven models do not have a physical meaning, no further process knowledge is needed, enabling a fast and easy implementation of these model types. Currently, various regression algorithms are available and commonly used, e.g., partial least squares, random forests, support vector machines, artificial neural networks (ANN) [9], and many more. However, such models are based on correlation and do not mandatorily imply causality, which can lead to inaccurate or even incorrect model predictions and conclusions. Contrarily, mechanistic models are based on theoretical considerations, i.e., the parameters have a physical meaning, and therefore ensure causality. Since these model predictions follow a purely mechanistic trend, temporary process deviations and unknown CPP impacts are not considered, which also interferes with the model performance and accuracy. To exploit the advantages of each individual model structure, a combined approach can be considered, called hybrid modeling (grey box) [10]. Since both models can complement each other in this combined structure, more precise predictions are anticipated. Such a hybrid model can be built in a parallel or serial structure, e.g., first, the data-driven part estimates parameters used in the mechanistic part, which otherwise would have to be assumed. Thereby, it is possible to also incorporate the CPP's impact into the hybrid model, which significantly strengthens the explanatory power of the model [11]. Additionally, to have assurance about the model performance and the risk of model mispredictions, typically cross-validation is performed to reduce variance, avoid overfitting, and investigate how the model performs when applied to new data [12]. A similar approach with a higher degree of freedom for creating the final model is model averaging from a leave-one-batch-out cross-validation, i.e., several developed models are averaged to improve the model stability and accuracy [13]. Even though this hybrid modeling approach has been the state-of-the-art in other industries for many years, due to the higher complexity of biological processes, it has only gained interest during the last few years [14]. Even though hybrid modeling is increasingly adopted for downstream applications [15–17], the response surface modeling of process endpoints is still more commonly applied [18], and the full potential of hybrid process modeling applications in bioprocessing has not yet been realized.

The high added value and the benefits of hybrid modeling for upstream bioprocessing become tangible when considering three major aspects of progressing towards digital biomanufacturing, i.e., delivering an increased process understanding, accelerating bioprocess development, and enabling advanced process control [19]. For all these components, various tools with different levels of complexity can be considered. Herein, soft sensors are frequently used, i.e., advanced online sensor systems such as spectrometry [20] or spectroscopy [21] in combination with a software algorithm to estimate the variables of interest in real-time, without any sampling and analytical time delay [22]. Depending on the used model structure, such soft sensors can be descriptive or predictive. While the descriptive model type can only be used to get estimated values up to the current time point, predictive models can also predict future values with a degree of uncertainty and therefore can additionally be used for process control [23]. Along with process models for the variables of interest, model-based methods for the optimization of process parameters such as the gained process information, the maximum amount of cells, or productivity were also introduced [24]. A highly interesting concept for accelerating bioprocess development and optimization in combination with model-related DoE approaches is a digital bioprocess twin [25,26]. Based on a minimal number of experiments, a hybrid model can be developed and subsequently be applied as a digital bioprocess. This digital twin then enables the simulation of further experiments, i.e., in silico exploration of the design space to shed light on the process behavior, without any additional laboratory experiments. This can be used to investigate the impact of the CPPs on the desired output, and thereby recommend the best CPP combination that maximizes it. A validation experiment at the recommended CPPs can be performed and compared to the simulation [27]. Subsequently, this digital twin model can be re-trained with the new experimental data, improving its performance by gaining a higher understanding of the process, and allowing it to explore a potential new optimum [28]. Once the recommendation of the digital twin converges at the process optimum, no new CPP combination will be proposed. Such model-based DoE and process modeling to find the best CPP combination in a design space saves raw materials and additionally operates more quickly and is cheaper compared to approaches in which experiments are only performed in the laboratory [29].

To accelerate the design space exploration and thereby greatly decrease the time needed to identify the optimum CPP combination for the variables of interest, we present a digital bioprocess twin used for model-based DoE [30]. This digital twin simultaneously delivers additional process understanding, while accelerating bioprocess development and optimization by applying in silico simulations that only perform the recommended experiments. We were particularly interested in determining the minimal number of required experiments for developing an initial digital twin, recommending further experiments to rapidly identify the best CPP combination in the design space. Such an iterative approach towards digitalization leads to a reduced experimental effort and saves various propositions of economic value while tackling current shortcomings for the implementation of such novel and promising tools [31]. Therefore, we present our structured workflow using different initial data sets to reduce experimental effort, evaluate the results, and additionally to investigate the applicability of an intensified DoE (iDoE) [32] for such a model-based DoE, to rapidly find the best CPP combinations in a design space and obtain the highest space-time yield.

#### **2. Materials and Methods**

#### *2.1. Experimental Design*

The experimental data set was derived from *E. coli* (HMS174 (DE3)) (Novagen, Germany) fed-batch cultivations at 20 L scale. For the workflow and the evaluation, a design space with three CPPs, each at three levels, was considered: the feed controlled specific growth rate *μ* (0.10, 0.15, and 0.20 h−1), the cultivation temperature *T* (30, 34, and 37 ◦C), and the induction strength *I* (0.2, 0.5, and 0.9 μmol IPTG g−1 cell dry mass), respectively. The variables of interest to be modeled were the biomass concentration (g L−1) and the space-time yield (g L−1 h−1) of the soluble fraction of the expressed protein, recombinant human superoxide dismutase. The biomass was analytically measured by thermogravimetric analysis [33] once before induction and then hourly, and the soluble product titer was measured every 2 h from the time point of induction to the last sampling at the end of the process by ELISA [34]. The fed-batch phase was carried out for four doubling times, and induction of the cells took place after the first doubling time, i.e., product formation took place for the remaining three doubling times. The values for the online measurements were available every minute and included the pH (controlled by the addition of 12.5% NaOH), off-gas (%), cultivation temperature (◦C), inlet air (slpm), dissolved oxygen (%), stirrer speed (rpm), base consumption (L), accumulated feed (L), inducer (kg), and head pressure (bar). More details about the applied exponential feeding strategy for the fed-batch phase, the utilized *E. coli* strain, the expression vector system, the online monitoring, and the offline measurements have already been presented elsewhere [35–37].

To receive meaningful information about the performance of the different digital twins and model-based DoE approaches, the design space was completely characterized. Once by common static cultivations (one CPP combination per experiment, i.e., 27 cultivations to cover all CPP combinations) and by iDoE cultivations (three CPP combinations per experiment, i.e., nine cultivations covering all 27 CPP combinations).

The intra-experimental CPP shifts in the intensified fed-batch fermentations were performed after each theoretical cell doubling post-induction of the cells, with a temporarily increased sampling interval, and executed by adjusting the setpoint value of the feed controlled specific growth rate and cultivation temperature in the process control system. Additionally, the feasibility of these shifts and the exclusion of a potential memory effect on the cells is presented in detail elsewhere [38]. A list of all the performed experiments used for comprehensive comparison is given in Appendix A.1 (Tables A1 and A2). Moreover, for the static cultivations, the maximum experimental values of the variables to be modeled are indicated. For the intensified cultivations, the maximum values were not conclusive, due to the intra-experimental shifts and the resulting multiple characterized CPP combinations, and therefore are not displayed. The two complete DoE and iDoE data sets are presented extensively and available for download as supporting information for an earlier publication [38].

#### *2.2. Data Sets*

For the initial hybrid model building and the model-based DoE, different initial data sets were used, and the respective performances for identifying the best CPP combination, obtaining the highest space-time yield were compared. These data sets were assembled out of the presented static and intensified fed-batch fermentations:

#### *2.3. Hybrid Model Development*

#### 2.3.1. Model Building

For initial model training, the different data sets were considered. To deal with the small initial data sets, avoid loss of information, and provide a more robust basis for the digital twin simulations, for each practically performed experiment, two additional in silico experiments were generated, i.e., each performed experiment was available in triplicate. For these in silico experiments, an appropriate level of analytical error was considered as random noise for the biomass (up to 5%) and the soluble product titer (up to 10%). As model inputs, the cultivation temperature (◦C), the accumulated feed (L), and the accumulated inducer (kg) were chosen to estimate the two response variables: the biomass (g L−1) and the space-time yield (g L−1 h−1). Prior to model building, the input variables were standardized using the z-score. To predict the response variables, a serial hybrid model structure was implemented. The data-driven model, an ANN, embedded in the hybrid model, and applying a Levenberg–Marquardt regularization algorithm, was chosen to estimate the specific growth rate *μ* and the soluble product formation rate *vp*/*x* as propagated predictions for the mechanistic part. The ANN consisted of three layers. The nodes of the hidden layer used hyperbolic tangent transfer functions, while the output layer used linear transfer functions. The values derived from the ANN were subsequently used in the mechanistic model, as shown in Equations (1) and (2), where *X* is the biomass concentration (g L−1), P is the soluble product titer (g L−1), I*y*/*n* is the inducer switch (zero for no induction or one for induction), and D is the dilution rate (h−1). Herein, D is used as the comprehensive term to describe the ratio between the flow of all volume additions into the reactor (L h−1), i.e., substrate feed, inductor feed and base, and the

overall reactor volume (L), which comprises the initial volume and all the added volumes. Consequently, in Equation (3), the space-time yield (*STY*) was calculated with the soluble product titer (g L−1) divided by the current utilization time of the bioreactor (h). This Bioreactor Utilization Time comprised the duration of the sterilization in place, inoculum, batch, harvest, cleaning, and the respective feed time.

$$\frac{d\mathbf{X}}{dt} = \boldsymbol{\mu} \cdot \mathbf{X} - \mathbf{D} \cdot \mathbf{X} \tag{1}$$

$$\frac{\text{dP}}{\text{dt}} = v\_{p/\text{x}} \cdot X \cdot \text{I}\_{y/n} - \text{D} \cdot \text{P} \tag{2}$$

$$STY = \frac{\text{P}}{\text{Bioreactor Ultilization Time}} \tag{3}$$

#### 2.3.2. Model Validation

For validation of the model performance, leave-one-batch-out cross-validation was performed, i.e., the initial model was built on all but one experiment, and the parameters were optimized by applying them to the experiment left out. Once no further improvement was observed, the model training stopped. To find the optimal setting to fit the experimental data, the number of neurons and hidden layers were varied. While the number of neurons was individually adapted for each data set, a single hidden layer delivered the best performance in all cases with respect to the normalized root mean square error (NRMSE) in Equation (4), where *y* is the analytical value, *yˆ* is the estimated counterpart for each sampling point (t), *y¯* is the mean of the analytical values, and N the total number of observations.

$$\text{NRMSE}\left[\%\right] = \frac{\sqrt{\frac{1}{N} \cdot \sum \left(y\_{(t)} - \hat{y}\_{(t)}\right)^2}}{\overline{y}} \cdot 100\tag{4}$$

#### 2.3.3. Model Averaging

To assess the risk of model misprediction, averaging of the individual models was performed. This averaging of the estimations from multiple models represents a robust way to deal with model uncertainties. This approach allows selecting a single model from each of the cross-validations. Depending on the initial data set, the averaged hybrid models consisted of three to five individual models. To validate this averaged model performance and its uncertainty, the NRMSE was taken into account, along with its standard deviation (SD) (Equation (5)) and the prediction interval (PI) (Equation (6)), where *yˆ*average is the estimation of the averaged model, *yˆ*model is the estimation of the respective model, *i* the index of these models, and *n* is the number of observations for each time point.

$$\text{SD}\_{\text{(t)}} = \sqrt{\frac{1}{n-1} \cdot \sum \left( \mathcal{G}\_{\text{average}(\text{t})} - \mathcal{G}\_{\text{model}(\text{i})\_{\text{(t)}}} \right)^{2}} \tag{5}$$

$$\text{PI}\_{\text{(t)}} = \text{\mathcal{Y}}\_{\text{average}} \pm \text{SD}\_{\text{(t)}} \tag{6}$$

Subsequently, the final averaged hybrid models were transferred to a digital twin environment.

#### *2.4. Digital Twin Application*

The developed hybrid models were implemented as digital twins to simulate all experiments in the given design space. Therefore, the accumulated feed, the inducer, and the inducer switch were simulated according to the feeding strategy and process time of the individual constant CPP levels, according to the desired design space boundaries. Once the simulations were performed by the digital twin, a lookup table could be used to individually evaluate the digital twin simulations. This lookup table provides the options for investigating the simulations, i.e., find the minimum or maximum values for

the response variables and their respective associated CPP combination along with the process time duration. For this case study, the lookup table was used to find the optima (maximum value) for the space-time yield in all simulated experiments, i.e., recommending the CPPs to obtain this simulated value. To validate the derived recommendation of the digital twin, a laboratory experiment with the respective settings was performed. The new experiment was then added to the previous data set and the hybrid model was re-trained including the new setup and its findings. This model-based DoE for optimizing the spacetime yield was repeated until the digital twin identified the best CPP combination and no new CPP combination was recommended. The entire workflow of the model-based DoE is presented in Figure 1. This workflow was carried out for all of the different initial data sets presented before, to evaluate the possible minimum number of required experiments for each case.

The hybrid model development, digital twin simulation, and model-based DoE were accomplished in the Novasign GmbH (Vienna, Austria) hybrid modeling toolbox.

**Figure 1.** Schematic workflow of optimizing the space-time yield using model-based DoE. Starting with an initial set of experiments from a given design space (I), a hybrid model is developed (II) and transferred to a digital twin environment. Based on the hybrid model, the digital twin simulates all experiments of the design space and recommends the best CPP combination in the design space to obtain the maximum value of the variable of interest (space-time yield) (III). In the case of a new CPP recommendation, the experiment is performed, added to the training data, and utilized to re-train the hybrid model with the new process information (IV). Once no new CPP recommendation is obtained, the digital twin identifies the best CPP combination to maximize the space-time yield and the optimization stops (V).

#### **3. Results**

#### *3.1. Analytical Space-Time Yield Maxima in the Design Space*

To confirm the simulated values and correctness of the CPP recommendation by the digital twin, the space-time yield of each CPP combination was investigated. The analytical maximum space-time yield of each cultivation is presented as a response surface in Figure 2. For simpler visualization, the results are separated into the three levels of induction strength.

**Figure 2.** Response surfaces of the analytical space-time yield maxima in the design space. The maximum value of each CPP combination is displayed as a function of the specific growth rate and the cultivation temperature for each induction level: *I* = 0.2 (**A**), *I* = 0.5 (**B**), and *I* = 0.9 (**C**). The color indicates the values of the space-time yield from dark blue (lowest value) to red (highest value).

The graphical investigation of the analytical space-time yield of each CPP combination in Figure 2 reveals the local and global optima in the design space. While at induction level *I* = 0.2, the local maximum was found at 0.0726 g L−1 h−1 (*μ* = 0.10 h−1 and *T* = 34 ◦C), and the induction level *I* = 0.5 contained the global maximum at the center point (*μ* = 0.15 h−1, *T* = 34 ◦C, and *I* = 0.5) with 0.0997 g L−1 h−1. The local maximum at induction level *I* = 0.9 resulted in 0.0915 g L−1 h−1 (*μ* = 0.10 h−1 and *T* = 34 ◦C). This visualization demonstrates that a cultivation temperature of 34 ◦C seems to be highly favorable for product formation, along with a trend towards slower specific growth rates.

#### *3.2. Initial Training Data for the Model-Based DoE*

The objective for this model-based DoE for parameter optimization was to quickly identify the best CPP combination for the highest space-time yield in the design space. To determine the minimum number of required experiments to develop meaningful hybrid models, and applied as digital twins recommending the next experiments, different initial data sets were utilized (Section 2.2 Data sets). These comprised either static or intensified cultivations, as presented in Figure 3.

As presented in Figure 2 and Table A1, the best CPP combination in the design space to maximize the space-time yield was obtained at the center point. However, there was also a local maximum with a high space-time yield at the highest induction level, which is assumed to be challenging not to become trapped in. For the design space investigation and determination of this CPP combination, different approaches can be consulted, as presented in Figure 3. First, experiments at each CPP combination were performed, characterizing the entire space without comprehensive process modeling (Figure 3A). Using this approach, the optimum in the design space was found, but this was paired with a high experimental effort and therefore time and costs. This experimental effort can be reduced by selecting a fractional factorial design and process modeling, i.e., only certain CPP combinations are performed. For this comparison, three fractional factorial designs were performed with the center point and the corners of the design space, either using nine (Figure 3B), five (Figure 3C), or only three initial experiments to build the hybrid model (Figure 3D). Since the iDoE concept proved to be suitable for accelerating the process characterization, this approach was additionally considered. Therefore, a complete set of iDoE experiments (Figure 3E) and three fractional iDoE approaches (Figure 3F–H) were used. The initial experiments of these last seven approaches were used in combination with process modeling to find the optimal CPP combination for obtaining the highest space-time yield as fast as possible, and using the workflow presented in Figure 1.

**Figure 3.** Approaches with varying initial data sets to set up the digital twin for model-based DoE. To find the best CPP combination in the given design space, different approaches with varying initial numbers of experiments were used (blue circles and lines). The full factorial DoE without the need for comprehensive modeling (**A**) was consulted in addition to fractional factorial DoEs with nine (**B**) and five (**C**), as well as a minimal approach using three (**D**), initial static cultivations. Model-based DoE approaches using iDoE cultivations were performed with the complete iDoE data set (**E**) and three fractional iDoEs (**F**–**H**).

#### *3.3. Digital Twin Simulations of the Model-Based DoE*

Out of all the presented initial data sets for the model-based DoE parameter optimization, the fractional factorial DoE with five initial static cultivations performed best, i.e., the fewest total experiments were needed by the digital twin to identify the CPP optimum for the space-time yield. A graphical presentation of this model-based DoE is presented in Figure 4. The step-by-step progression of the recommended experiments in the design space along with the simulated values compared to the analytical values for each re-trained digital twin are shown.

The model-based DoE quickly recommended the best CPP combination to obtain the highest space-time yield (Figure 4A). The correct induction level was already found after implementing the gained process knowledge from the first recommended experiment and the correct cultivation temperature after the second re-training of the digital twin. Even though the specific growth rate was the most difficult to properly assert, after two additional cultivations the optimum in the design space was found, identifying the center point CPPs as the optimum process conditions, which were already present in the initial training data. This resulted in nine performed experiments instead of twenty-seven, highlighting the advantages of knowledge-based bioprocess development. However, with this small initial data set, the simulated biomass of the first recommended experiment (Figure 4B) almost matched the analytical results, and the space-time yield was highly overestimated. Likewise, high overestimations were observed for the second (Figure 4C) and the third recommendation (Figure 4D). By adding these new recommended experiments to the initial data set, the resulting retrained hybrid model iteratively gained knowledge about the process for the next recommendation. Already, after only these three re-trainings, the fourth simulation almost converged on the analytical values (Figure 4E). The digital twin gained precision and certainty at the fifth and final recommendation (Figure 4F). Since this recommended experiment had already been performed, the model-based DoE stopped, i.e., the best CPP combination was identified, and the biomass and space-time yield of the process were accurately simulated.

**Figure 4.** Step-by-step progression of the model-based DoE and the performance compared to experimental data. The initial data set (blue circles), the recommended experiments (orange circles), and the temporal order (orange arrows) are given (**A**). For each recommended experiment (**B**–**F**), the simulated biomass (green lines) and the simulated space-time yield (blue lines) are presented along with the PI (shaded area). The time point of induction (dashed grey line) and the mean analytical values for the biomass (green diamonds) and the space-time yield (blue triangles) are indicated along with the SD (error bars).

With five initial static experiments, the digital twin simulated the biomass concentration with an appropriate accuracy from the beginning, but highly overestimated the experimental values of the space-time yield. By consecutively adding the four recommended experiments, and extending the initial data set, precise simulations were obtained. This fast convergence of the simulated space-time yield on the analytical values, along with the SD, is displayed in Table 1.

**Table 1.** Progression of the model-based DoE until the optimum was found using five initial experiments.

As seen in Table 1, the obtained recommendations of the digital twin, at which CPP combination the next experiment should be performed, converged at the best CPP combination in the design space after five recommended experiments, i.e., no new recommendation was derived. Moreover, a steep learning curve of the hybrid model was observed when the new experiments were added for re-training the digital twin. While the simulated spacetime yield of the first recommended experiment, derived from the information gained from the initial five experiments, resulted in an 8.68-fold deviation compared to the analytical value, this factor quickly decreased after including the respective validation experiments in the training data and subsequent re-training of the hybrid model. For example, the simulation of the second recommendation already displayed a decreased deviation of only 1.75-fold compared to the analytical value, while the third simulation was down to a

1.59-fold deviation. The fourth simulation only displayed a deviation from the analytical value by 1.12-fold, and the final simulation of the fifth recommendation was highly precise, displaying a simulated maximum of 0.98-fold the analytical value. This demonstrates that with only five initial experiments to start the model-based DoE, the hybrid model promptly gained process knowledge and its digital twin was able to provide the best CPP combination to obtain the highest space-time yield.

A complete quantitative and qualitative performance comparison of all the presented approaches (Figure 3) is given in Table 2. Herein, the three different fractional iDoE approaches are summarized.

**Table 2.** Performance summary of the model-based DoE approaches. Capital letters in brackets represent the DoE conditions from Figure 3.

Table 2 presents the quantitative effort and qualitative performance of each initial data set. With respect to the total required time for each presented approach, only the duration of the practical experiments (including pre- and post-processing) was taken into account for the evaluation, since using our setup, an entire experiment takes approximately one working week. However, the computational time for the hybrid model training and subsequent re-training can be neglected, since it ranges between half an hour and three hours, and depending highly on the performance of the utilized computer. While the number of required experiments remains unchanged, the needed experimental time can further be reduced by the utilization of multiple bioreactors or parallel bioreactor systems.

Since in the full factorial DoE all experiments are performed, comprehensive process modeling is redundant to find the best CPP combination for the highest space-time yield in the design space. By using this approach, the optimum was found, but paired with the highest experimental effort. For the other initial data sets, model-based DoE was applied to reduce the required number of experiments. For the fractional factorial DoEs, the number of recommended experiments increased until the optimum was found when decreasing the number of initial experiments. Herein, the fastest approach was the fractional factorial DoE with five initial experiments and four validation experiments required, i.e., only 9/27 experiments had to be performed. Moreover, in all cases, the optimum was identified. However, in this case study, the utilization of initial iDoE cultivations for model-based DoE did not lead to the identification of the best CPP combination in the design space. Regardless of selecting the entire iDoE data set or varying fractional iDoEs, the model-based DoE ended up at different locations in the design space than the optimum CPP combination. Herein, the final recommendations by the digital twin were all located at *μ* = 0.10, *I* = 0.9 and either 30 ◦C or 34 ◦C, indicating a model bias towards slow specific growth rates and temperatures, apart from 37 ◦C, where a high value or local maximum of the space-time yield is located. A more detailed progression of the recommended experiments in the design space for each of the other six model-based DoEs is shown in Appendix A.2, Figure A1 (excluding the full factorial DoE).

#### **4. Discussion**

The prominent emerging concept of model-based DoE for parameter optimization is an interesting, and yet not completely explored, topic. To accelerate this identification of optimum process conditions is of great interest for manufacturers, to reduce bioprocess development timelines. Typically, by performing all experiments in a design space, these optimum process conditions can be found, but with high experimental effort. Herein, we challenged this approach by investigating the minimum requirements for such a modelbased DoE workflow (Figure 1) to rapidly and properly discover the best CPP combinations in a design space (Figure 2), utilizing varying numbers of initial experiments (Figure 3). We demonstrated with our case study that the fastest approach to identifying the best process conditions for the highest space-time yield was an initial fractional factorial DoE with five static cultivations and four consecutively performed recommendations from the digital twin (Figure 4 and Table 1). In case scientists are limited to certain time slots for further experiments, the best x-recommendations from the digital twin can be used in the next campaign to obtain the maximum learning, according to the experimental possibilities. Interestingly, all model-based DoEs using initial iDoE cultivations failed to find the global maximum in the design space (Table 2), and recommending an incorrect optimal CPP combination after a few iterations (Figure A1). It has already been demonstrated that iDoE is favorable for accelerating process characterization. Here, a trade-off between decreased experimental effort and reduced process information can be accepted. This consideration must be handled with care when iDoE is used for process optimization, i.e., an increased model uncertainty due to decreased process information may result in divergent optima, as was the case herein. To the best of our knowledge, this iDoE concept has not been well investigated and little literature is available as a reference for microbial, and even less for mammalian, systems. Additionally, several degrees of freedom are introduced by iDoE, e.g., the number and duration of the intra-experimental CPP shifts, as well as how these should be performed. Therefore, before reliably applying iDoE for such model-based DoE approaches, more research should be performed on this subject.

Furthermore, the identification of optimum process conditions for the response to be optimized in design spaces with a higher dimensionality, as in our case study (>3 CPPs), could lead to new challenges, e.g., the occurrence of various local optima, which complicate the accurate identification of the global optimum. The robustness and applicability of digital twins to also perform reliably when confronted with this higher complexity must be further investigated. Moreover, our findings demonstrate that bioprocess modeling is not an all-in-one solution, eliminating all current limitations and obstacles; showing that it is important to consider many potentially influencing factors [39].

For instance, it is advisable for the initially used data set to introduce every CPP level to the hybrid model training, i.e., the minimal fractional factorial DoE with three initial cultivations in our case study. Otherwise, the hybrid model will be biased towards the included CPP levels in the training data and potentially would not recommend the missing setting, since the ability to correctly determine these causal relationships is lost. This bias towards CPP levels should be considered when initially investigating a design space, for which no prior process knowledge about process behavior and the responses is available, i.e., the CPPs and the appropriate levels should be well-considered and not too far apart. Hereby, the accidental generation of independent data sets, becoming missing and getting trapped in local optima, can be avoided at the start. Since this case study mainly focused on the practical application of digital twins, more detailed theoretical analysis should be performed in future studies. However, it might be desirable to re-define the CPP levels and look for new, more beneficial settings in the design space, e.g., with smaller intervals of the cultivation temperatures simulated by the digital twin. However, if a digital twin recommends an experiment next to the identified optimum CPP combination, but with a 0.5 ◦C decreased cultivation temperature and an increased space-time yield by 0.3%, the execution of this cultivation should be critically questioned. Additionally, for some CPPs, such simulated intervals are not always practically feasible, e.g., steps of 0.5 ◦C for the cultivation temperature, which might be adjustable but difficult to precisely control. This exemplary scenario demonstrates that such approaches must still be guided by human knowledge, rather than completely trusting an algorithm.

Herein, it has been demonstrated that such digital solutions enable a new knowledgebased perspective on bioprocess development and optimization, and to get more out of the

available data. Even though several of these advantages have already been recognized and discussed, much more research will be required to fully implement and exploit the potential of digitalization in the biopharmaceutical industry [40]. For instance, an up-and-coming area for future application of model-based DoE, hybrid modeling, and digital twins is found in simulating new CPP combinations out of the design space, i.e., extrapolation where appropriate. However, this again poses new challenges, such as how to validate this new setting outside the design space, e.g., an additional smaller design space with the new CPP combination as the center point could potentially be performed. Besides the validation issue, the stability of the digital twin and the underlying hybrid model structure must also be ensured. Additionally, if the mechanistic relationships are known and understood, such digital twins could be used as a basis to initially simulate new bioprocesses with similar product properties without prior experiments, e.g., product size and cytotoxicity supporting platform approaches.

#### **5. Conclusions**

In silico design space exploration using a digital bioprocess twin increases the process understanding for QbD; the impact of the CPPs on the variables of interest can rapidly be investigated. The presented workflow enabled us to quickly find process optima in a design space despite using only a small initial experimental setup. Moreover, this approach to decreasing the number of required practical experiments for process optimization becomes even more advantageous for larger design spaces. Even though, herein the dimensionality and complexity increase, which will lead to new challenges, model-based DoE has the potential to significantly lower the experimental effort; saving money, time, raw materials, and other propositions of economic value for later stages.

**Author Contributions:** Conceptualization, B.B., M.D., R.D.D.; methodology, B.B., R.D.D.; software, R.D.D.; validation, B.B., R.D.D., M.M.; formal analysis, B.B., R.D.D.; investigation, B.B.; resources, G.S.; data curation, B.B., R.D.D.; writing—original draft preparation, B.B.; writing—review and editing, B.B., R.D.D., M.M., G.S., M.D.; visualization, B.B.; supervision, M.D.; project administration, M.D.; funding acquisition, G.S., M.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Austrian Research Promotion Agency (FFG), grant number 859219.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available in the supporting information of a previous publication (https://doi.org/10.1002/biot.202000121) (accessed on 24 June 2021).

**Conflicts of Interest:** Benjamin Bayer, Roger Dalmau Diaz, Mark Duerkop, and Gerald Striedner hold shares of Novasign GmbH.

#### **Abbreviations**

#### **Appendix A**

#### *Appendix A.1 CPP Settings of All Experiments Used for Model-Based DoE*

The design space, the herein investigated CPPs (and respective levels), and cultivation approaches are introduced in the Materials and Methods section of the main manuscript. A detailed list of all performed experiments of the comprehensive comparison for the applicability of the model-based DoE workflow (Figure 1, main manuscript) is given below. Table A1 provides information about the experiments performed with one CPP combination, and Table A2 contains the intensified experiments (three CPP combinations per cultivation) and the herein performed CPP shifts. For all static experiments, the maximum experimental values of the variables modeled (biomass and space-time yield) are provided for easier comparison. For the intensified experiments, these maximum experimental values are not indicated, because these quantities are not meaningful due to multiple characterized CPP combinations per experiment. The highest space-time yield in the entire design space was obtained at CPP combination #14 (μ = 0.15 h−1, T = 34 ◦C, and I = 0.5), reaching 0.0997 g L−1 h−1 in the performed cultivation. Subsequently, the different initial data sets were evaluated in the model-based DoE (Figure 3, main manuscript), considering the number of required recommendations by the digital twin until certainty about the best CPP combination is gained.

**Table A2.** CPP combinations of the intensified experiments for the model-based DoE approach.

#### *Appendix A.2 Progression of the Recommended Experiments by Each Model-Based DoE Approach*

Out of all presented initial data sets for the model-based DoE in Figure 3 (Results section of the main manuscript), the fractional factorial DoE with five initial static experiments proved to be the fastest for identifying the best CPP combination for the highest space-time yield. This detailed progression until the optimum was found is presented in Figure 4 and Table 1 (Results section of the main manuscript). For the other six data sets used for the model-based DoE (excluding the full factorial DoE), Figure A1 presents an overview of the respective progressions, including the initially performed experiments, as well as the recommended experiments.

Besides the best performing model-based DoE with five initial static cultivations, the two other initial fractional factorial DoEs also performed well. The approach with nine initial static cultivations (Figure A1A) needed two recommendations, i.e., two further experiments to gain certainty about the optimum CPP combination, resulting in a total of 11/27 cultivations. Herein, the model quickly gained certainty about the correct induction level from the beginning, and after the second experiment also about the other two CPP levels. The model-based DoE using three initial static cultivations performed seven recommendations until the optimum was identified, i.e., 10/27 cultivations (Figure A1B). Interestingly, here the induction level was also the first CPP to be correctly recommended after two additional experiments, followed by the cultivation temperature and then the specific growth rate. However, the complete iDoE as the basis for model-based DoE (Figure A1C) did not identify the optimum, and after two recommendations by the digital twin ended up recommending CPP combination #7 (*μ* = 0.10 h−1, *T* = 30 ◦C, and *I* = 0.9). Moreover, the model-based DoE based on three different fractional iDoEs was also not able to find the optimum CPP combination. Depending on the initially selected three iDoE cultivations, it took one to four recommendations by the digital twin until these modelbased DoEs also recommended CPP combination #7 (Figure A1D,F) or CPP combination #8 (*μ* = 0.10 h−1, *T* = 34 ◦C, and *I* = 0.9) (Figure A1E) as the best CPP combination for the highest space-time yield.

**Figure A1.** Step-by-step progressions of the recommended experiments by the model-based DoE, using varying initial data sets. The initial experiments (blue circles and lines) and the respective recommendations for the next experiment (orange dots), along with the temporal order (orange arrows) are given. The fractional factorial DoE with nine (**A**) and three (**B**) initial static cultivations, the complete iDoE (**C**), and the three fractional factorial iDoEs (**D**–**F**) are presented.

In conclusion, while every approach using static cultivations as a basis for modelbased DoE could identify the optimum CPP combination in the design space, all the iDoE approaches failed to do so. However, it was already shown that the concept of iDoE is advantageous for reducing the experimental effort for process characterization but, in this particular case, it was not possible for model-based DoE to accurately identify the static process conditions optimizing a certain process output. The recommended CPP combinations by the model-based DoE with the initial iDoE cultivations were becoming trapped at high values or local maxima and were highly biased towards the highest induction level, slowest growth rate, and lower cultivation temperatures.

#### **References**

# *Article* **Designing Robust Biotechnological Processes Regarding Variabilities Using Multi-Objective Optimization Applied to a Biopharmaceutical Seed Train Desig**

**Tanja Hernández Rodríguez 1, Anton Sekulic 2, Markus Lange-Hegermann 2 and Björn Frahm 1,\***

**Abstract:** Development and optimization of biopharmaceutical production processes with cell cultures is cost- and time-consuming and often performed rather empirically. Efficient optimization of multiple objectives such as process time, viable cell density, number of operating steps & cultivation scales, required medium, amount of product as well as product quality depicts a promising approach. This contribution presents a workflow which couples uncertainty-based upstream simulation and Bayes optimization using Gaussian processes. Its application is demonstrated in a simulation case study for a relevant industrial task in process development, the design of a robust cell culture expansion process (seed train), meaning that despite uncertainties and variabilities concerning cell growth, low variations of viable cell density during the seed train are obtained. Compared to a non-optimized reference seed train, the optimized process showed much lower deviation rates regarding viable cell densities (<10% instead of 41.7%) using five or four shake flask scales and seed train duration could be reduced by 56 h from 576 h to 520 h. Overall, it is shown that applying Bayes optimization allows for optimization of a multi-objective optimization function with several optimizable input variables and under a considerable amount of constraints with a low computational effort. This approach provides the potential to be used in the form of a decision tool, e.g., for the choice of an optimal and robust seed train design or for further optimization tasks within process development.

**Keywords:** Gaussian processes; Bayes optimization; Pareto optimization; multi-objective; cell culture; seed train

#### **1. Introduction**

The development and optimization of biopharmaceutical production processes with cell cultures is cost- and time-consuming, requiring substantial lab work. This necessitates thorough planning of experiments and processes, taking into account existing process knowledge. The need for model-based decision support in biopharmaceutical manufacturing has been emphasized by the US Food and Drug Administration (FDA) [1,2], including taking into account available prior know-how and experience within the decision process and uncertainties [3]. Such methods are still not state-of-the-art for cell culture processes during development or manufacturing [3,4], although first approaches have been proposed, for example, in order to optimize the titer of a mammalian cell culture process [5]. This highlights a need for improved methods and tools for optimal experimental design, optimal and robust process design and process optimization for the purposes of monitoring and controlling during manufacturing.

But also in other engineering fields such as chemical engineering or mechanical engineering, process optimization plays an important role and is the subject of current research. Some application examples rely on dynamic models, an example is the optimization of

**Citation:** Hernández Rodríguez, T.; Sekulic, A.; Lange-Hegermann, M.; Frahm, B. Designing Robust Biotechnological Processes Regarding Variabilities Using Multi-Objective Optimization Applied to a Biopharmaceutical Seed Train Design. *Processes* **2022**, *10*, 883. https:// doi.org/10.3390/pr10050883

Academic Editor: Alina Pyka-Paj ˛ak

Received: 6 April 2022 Accepted: 21 April 2022 Published: 29 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

sustainable algal production processes [6] or the improvement of the vibration performance of cold orbital forging machines [7]. Other approaches rely on machine-learning algorithms such as those reported in [8–10].

The optimization of one objective criterion (e.g., final titer) is relatively straight forward, i.e., building an objective function with a unique response variable and applying an appropriate optimization algorithm to maximize this function. However, in industry, it is typically desired to optimize several conflicting objectives at a time, leading to suitable trade-offs and compromises. For example, when trying to maximize final titer via viable cell density while minimizing cultivation time. Multi-objective optimization provides a decision-making tool for optimal decisions in the presence of trade-offs between two or more conflicting criteria.

However, multi-objective optimization is more challenging. Its application is still not state-of-the-art in the context of cell culture processes, probably due to a lack of related studies and instructions. Moreover, within the manufacturing life cycle of biopharmaceuticals, some phases are better investigated than others. Still very few investigations are reported concerning the cell expansion process (seed train). It consists of several consecutive cultivation and passaging (transfer) steps, starting with a small amount of cell suspension because cells are frozen in small vials until they are used for a production process. The goal is to expand the number of viable cells in order to reach the required amount to inoculate (start) the production bioreactor (e.g., 10,000 L at industrial scale) while keeping them in a healthy and growing state. A high amount of operational requirements and constraints have to be fulfilled and, as reported in literature [11,12], the cell expansion process critically effects product quality and the amount of product at production scale. In [12] for example, the passage duration, as well as the initial viable cell density for each passage are reported as important parameters with high impact on process time and productivity at production scale. A careful and optimal planning of a seed train is therefore essential. However, this is not a trivial task due to the inherent variability concerning cell growth (cell growth differs from cell line to cell line and also from cultivation run to cultivation run) and uncertainty about the real state of the process due to considerable measurement uncertainties. This requires the design of a reproducible process which is robust regarding viable cell density, meaning that despite (initial) variabilities concerning cell growth, low variations of viable cell density at the end of the seed train are obtained. The goal of this paper is to close the gap between state-of-the-art optimization techniques and modern techniques from machine learning to improve the biopharmaceutical production by allowing easy to use yet powerful multi-objective optimization.

In most multi-objective optimization problems, no single best (unique optimal) solution exists, instead there is a set of optimal solutions (also called Pareto optimal solutions or non-dominated solutions), meaning for each solution that one criterion cannot be improved without degrading at least one of the other criteria. So, the decision maker has to choose from the set of non-dominated solutions according to the most preferred or important objective criterion. A promising approach to optimize objective functions, which are expensive to evaluate, is Bayes optimization. The methodology of Bayes optimization dates back to the work of Harold Kushner in 1964 [13] and gained impact through the work of Jones et al. in 1998 [14]. It is a probabilistic global optimization method for finding the maximum of objective functions that are expensive to evaluate or unknown (black-box) objective functions that are approximated using simulations [15].

In practice, the objective function could be the outcome of interest of a process, for example, process productivity or control metrics to describe the quality of a product. Input parameters can be process parameters needed to be optimized. Bayesian optimization [16] creates a quick to evaluate model, the so-called surrogate model of the objective function. In order to reduce the objective function evaluations, the surrogate model is iteratively trained and updated on new data. The positions of this new data are chosen by finding a trade-off between exploration (improving the surrogate model) and exploitation (finding optimal points). Typical surrogate models are Gaussian processes.

Gaussian processes (GP's) are popular machine learning models [17] because, due to their Bayesian nature, they work well with few data points [18]. Furthermore, they allow the inclusion of expert knowledge [19,20] and can be used in dynamic systems [6,21]. GP's are very flexible non-parametric models, hence, they can approximate any function and do not assume a predefined set of modeling functions.

Bayes optimization is successfully applied in many fields of research and economics [22]. Moreover, applications of Bayes optimization in the field of bioprocess engineering were published during the last decade [6,9,23,24]. Furthermore, this methodology was shown to be efficient in solving multi-objective optimization problems [25] and has also been applied for parameter estimation of kinetic parameters [26]. However, no applications are reported so far applying model-based multi-objective Bayes optimization within biopharmaceutical process development.

This contribution aims to present the concept of a workflow which couples uncertaintybased upstream simulation and Bayes optimization using Gaussian processes and its application in the form of a simulation case study to illustrate its applicability to a relevant industrial task in process development.

This simulation case study addresses the question if a reference seed train setup comprising five shake flask scales can be optimized through varying shake flask volumes and how many shake flask scales, three, four or five, are recommendable in terms of two objective criteria, seed train duration and deviation rate. Moreover it is investigated how the results change if cells grow with 5% lower or 5% higher maximum cell-specific growth rate.

Afterwards, two more objective criteria, titer (product concentration) and viability after 8 days in the production bioreactor, are added and seed train optimization is performed regarding four objective criteria simultaneously.

Furthermore, the suitability of the proposed method and the required number of iterations is evaluated with respect to the obtained information gain.

#### **2. Methods**

The main components of the applied methodology and the corresponding tools are described.

#### *2.1. Upstream Simulation*

Upstream simulation comprises a simulation of the cell expansion process (seed train) and simulation of the production scale. The reference upstream process taken as an application example for the here presented simulation case study comprises five consecutive shake flask scales followed by three bioreactor scales and one production scale, similar to the upstream process investigated in [27]. Further specifications are listed in Table 1.

A mathematical model is required, describing cell growth and interactions with the main limiting substrates and eventually inhibiting metabolites over time. A cell growth model, a system of ordinary differential equations (ode), already adapted to an industrial cell culture upstream process using a CHO cell line [27] has been used, which describes the dynamic behavior of viable and total cell density, *X*v and *X*t, concentrations of glucose *c*Glc, glutamine *c*Gln, lactate *c*Lac, ammonia *c*Amm and product (volumetric titer) *c*titer (see Table A1 in the Appendix A).

Moreover, such an upstream process includes several constraints, operation steps and process parameters (e.g., concerning passaging intervals, substrate/nutrient concentrations, initial viable cell densities and viable cell densities before transferring cells into the next cultivation vessel, as well as the amount of cell suspension and fresh medium), which have to be considered in the simulation workflow. A detailed description of the required components and calculation routines are described in [28,29].

Besides these requirements, several passaging strategies can be applied, helping to decide at which point in time cells should be transferred from one cultivation vessel into the next larger one and how to perform these passaging steps (e.g., which amount of cell suspension should be mixed with how much fresh cell culture medium).

For the here presented simulation study, the passaging strategy for robust seed train design was chosen, where robustness refers to the reproducibility of the seed train regarding viable cell density, meaning that despite initial uncertainties and variabilities concerning cell growth, low variations of viable cell density at the end of the seed train are obtained. This strategy grounds on the objective of reaching the previously determined threshold of viable cell density and corresponding probability distributions of viable cell density at different points in time. These distributions are used in combination with a utility function following the mean-variance principle, which grounds on the Markowitz mean-variance portfolio optimization theory [30,31]: The utility function *U*(*t*) is defined as a function of viable cell density *X*v including the expected value E(*X*v) and the variance Var(*X*v) of viable cell density, as well as a risk aversion parameter *α* which controls the amount of risk (amount of uncertainty) the user is willing to bear. In the here presented example, the risk refers to the probability that viable cell density differs from the expected value (predicted mean). A risk aversion value of *α* = 1 would mean that the expected time profile minus one time the standard deviation of *X*v is considered.

The utility function is defined through:

$$\mathcal{U}(t) = \mathbb{E}(\mathbf{X}\_{\mathbf{V}}(t)) - \mathfrak{a}\sqrt{\text{Var}(\mathbf{X}\_{\mathbf{V}}(t))}\tag{1}$$

Based on the simulated time profiles of the current cultivation scale (by solving the corresponding ode system), Equation (1) is used to calculate the utility function value *U*(*t*) per hour and to check if this value reaches or exceeds the required transfer viable cell density *X*v,transfer which is necessary to inoculate (start) the next cultivation scale fulfilling the required seeding (initial) viable cell density and the filling volume.

In the next step, it is evaluated whether the calculated point in time lies within the range of practically feasible points in time for cell passaging, *T*p. Thus, the objective is to find the minimum point in time out of the set of practically feasible points in time for passaging, *T*p, which fulfills:

$$
\mathcal{U}(t) \ge \mathcal{X}\_{\text{v,transfer}} \tag{2}
$$

$$\text{subject to:} \quad t \in T\_{\text{p}}.\tag{3}$$

Based on the obtained point in time and the corresponding concentrations of viable cells, total cells, substrates and metabolites at this point in time, starting concentrations (=initial values of the system of ordinary differential equations) of the next cultivation scale are calculated based on the defined configurations and constraints (e.g., working volumes, acceptable range of seeding viable cell density and medium concentrations). This calculation has to be performed for every cultivation scale and passaging step. For more details refer to [27–29].

#### *2.2. Bayes Optimization*

A typical mathematical optimization problem is the following: Given an objective function *f* : X → R over input space X ⊆ R*d*, the aim is to find an argument *x*∗ ∈ X , which optimizes (minimizes or maximizes) *f* .

The idea behind Bayes Optimization consists of creating a simple, probabilistic and cheap to evaluate model, a so-called surrogate model (substitute model), of the objective function *f* [15,17,32]. Bayesian optimization reduces the number of evaluations of the objective *f* via the following iterative approach: Before sampling *f* at another point, we take into account a trade-off between exploration (i.e., sampling of areas of high uncertainties) and exploitation (sampling from areas which are likely to move towards the optimum), which is encoded in a so-called acquisition function. We can find such points quickly from evaluation of the surrogate model.

Within Bayes optimization the following steps are performed:

The key of Bayesian optimization is not to rely on local approximations as many other optimization algorithms and instead to have a global viewpoint of also evaluating the function at unknown positions.

The acquisition function is used to propose the next candidate point to be evaluated based on specific criteria, for example the expected improvement of the optimization criteria, and on the reduction in predictive uncertainty. As in the case of the kernels, there is also a wide variety of possible acquisition functions to choose from. In this study, the Expected Improvement (EI) acquisition function is used [33,34].

Gaussian processes (GPs) are well suited surrogate models when making few assumptions [15]. Just like a Gaussian distribution (a normal probability distribution) is fully described by its mean *m* and variance *σ*2, a GP is fully described by a mean function *m*(*x*) and a covariance function *k*(*x*; *x* ) [17]. A GP is an extension of a multivariate Gaussian (or normal) distribution to distributions of functions in the sense that if a function *y* follows a GP distribution, i.e., *y* ∼ GP(*m*, *k*), then every evaluation of the function follows a Gaussian distribution *y*(*x*) ∼ N (*m*(*x*), *k*(*x*, *x*)). In particular, a GP returns mean and variance of the possible function values (instead of just returning a scalar), and hence also provides information about the uncertainty of a prediction. Moreover, GPs can take into account uncertainty in the form of noise, the class of Gaussian processes is closed under Bayesian updates, and such updates are computationally tractable [35].

The covariance function describes the assumed characteristics such as smoothness or periodicity of the objective function *f* [16]. They are so-called positive-definite functions, often also called kernels [17,36]. It specifies the relationship between two 'points' (vector of the input space) *x* and *x* and the corresponding changes in *f* at these points. A covariance function is described by a set of parameters, also called hyperparameters, describing a specific behavior. This is how prior information is embedded in the Bayes optimization procedure. Also in this work, the most commonly used covariance function, the Squared Exponential (SE) kernel (often also referred to as Gaussian kernel) is used [32].

#### *2.3. Problem Definition and Computational Procedure*

The goal of the presented application example is to propose a concept and a numerical procedure for optimal robust seed train design, where robustness refers to the reproducibility of the seed train regarding viable cell density, meaning that despite initial uncertainties and variabilities concerning cell growth, low variations of viable cell density at the end of the seed train are obtained.

First, seed train constraints are defined based on a chosen cell line and its characteristics concerning optimal cultivation conditions and based on the operative possibilities (e.g., feasible points in time for cell passaging). Second, the optimizable input parameters and objective criteria (objective response variables) applied in this study are defined (as also illustrated in Figure 1), followed by the formulation of the mathematical optimization problem. Thereafter, the optimization problem is solved using a workflow which connects seed train simulation and Bayes optimization.

**Figure 1.** Goal of the study is to propose a concept and a numerical framework for optimal robust seed train design (blue box in the middle), including optimizable inputs (first gray box) as well as objectives (objective criteria) used in this study (right gray box).

The following objective criteria were chosen to represent an optimal seed train: (I) a minimum duration (*d*) (=required cultivation time) of the seed train and (II) a minimum deviation rate (*D*) regarding viable cell density, i.e., the probability that the seed train will run outside predefined ranges of viable cell density (for both, seeding viable cell density and transfer viable cell density) (This is important to consider because in the case that specific constraints are not fulfilled, the performance of the cells could decrease. The growth rate could decrease and, furthermore, it has been observed that the violation of constraints could also cause less viability of the cells in the production phase [12]) (compare to Figure 1 right gray box). These two attributes shall enable an optimal start of the production scale. Note that, in addition to these criteria, the growth rate is another important parameter affecting an optimal start of the production scale and the growth rate should be high until the end of the seed train. However, in this first optimization study it is not set as optimization criterion because the here defined seed train setup (in terms of medium concentrations and possible cultivation volumes per scale) together with the aim to reduce cultivation time already supports good growth during the entire cultivation. However, for other seed train setups, it might be advisable to include growth rate at the end of the seed train into the optimization problem.

After consideration of the two mentioned objective criteria, a third and fourth objective criterion, the product concentration (titer) and the viability at the end of the cultivation in the production bioreactor (in this simulation study: after 8 days in batch mode, i.e., without addition of nutrient feeds) are added to the optimization problem (see Figure 1 right gray box (III)). Note: The authors are aware of the fact that cultivation in the production vessel itself, which is often performed in fed-batch mode, is also influenced by several process parameters having an impact on product quantity and quality. Moreover, data of further attributes would be necessary to describe product quality (e.g., of a recombinant therapeutic protein or antibody) but these are not provided and therewith not considered in this study.

The input variables that can be varied to optimize the recently mentioned objective criteria, and thus the optimizable input variables, are the filling volumes in the first five shake flask scales, *V*1, ... , *V*5 (compare to Figure 1, the part of the seed train between thawing cells from a small vial and inoculation of the first biorector). These target values are important inputs of the seed train simulation process because they are used to calculate points in time for cell passaging. Volumes in the finally proposed seed train protocol (output of the seed train simulation) may vary within allowed working volume ranges and these are also presented in this work.

Formulation of the Mathematical Optimization Problem

The optimizable variables and therewith inputs of the optimization problem are the filling volumes of the *n* shake flask scales, *V*1, ... , *Vn* which are included in the input vector:

$$\mathfrak{X} = (V\_1, \dots, V\_n)^T. \tag{4}$$

Outputs of the optimization problem are the defined objective criteria. These are seed train duration *d* and deviation rate *D* for the first optimization example. Thus, the unknown objective function (which should be minimized) can be written as follows:

$$f(\mathbf{x}) = (f\_1(\mathbf{x}), f\_2(\mathbf{x}))^T \tag{5}$$

with *f*1(*x*) =ˆ *d* and *f*2(*x*) =ˆ *D*.

The second optimization example includes a third and fourth optimization criterion, product concentration and viability at the end of the production scale (here after 8 days in the production vessel). Thus *f*(*x*) expands to:

$$f(\mathbf{x}) = \left(f\_1(\mathbf{x}), f\_2(\mathbf{x}), f\_3(\mathbf{x}), f\_4(\mathbf{x})\right)^T \tag{6}$$

with *f*1(*x*) =ˆ *d*, *f*2(*x*) =ˆ *D*, *f*3(*x*) =ˆ *c*titer,end and *f*4(*x*) =ˆ Viabilityend.

#### *2.4. Connecting Seed Train Simulation and Bayes Optimization*

Uncertainty-based seed train simulation as described in Section 2.1 was coupled with algorithms for Bayes optimization as described in Section 2.2. The workflow integrating both components is illustrated in Figure 2. The inputs of the combined framework are the input variables: Boundaries for the optimizable variables (here filling volumes) and objective criteria (here seed train duration, deviation rate and in the second example also product concentration at the end of production scale) given all required seed train configuration settings and constraints (e.g., initial concentrations, practically feasible points in time for cell passaging, acceptable ranges for viable cell density, . . . ).

**Figure 2.** Scheme showing the applied computational workflow comprising: (**A**) a Bayes optimization algorithm which is coupled with (**B**) a seed train simulation routine. Input and output values are shown in the blue boxes above and below.

First points (=combinations of optimizable variables) are determined using a Latin Hypercube design distributing these points within the design space (see Figure 2, Box A). Seed train simulations are performed at these points in order to obtain the corresponding objective criteria values. Input values together with output values form a data set. An unknown model describing the relationship between inputs and outputs is approximated through a Gaussian process (GP) which has to be trained (see Figure 2, Box A) based on the given data set. Therefore, the Gaussian process proposes a point that has to be evaluated next (see Figure 2, Box B).

A robust seed train is simulated, using a mechanistic process model, and the objective criteria are calculated. This output is then returned to the Bayes optimization (Box A) to update the GP. Usually, experiments are performed to return the experimental output. The present approach instead exploits the advantages of the model-based upstream simulation in order to reduce the experimental effort to a minimum.

These steps are repeated various times, e.g., until a previously defined number of maximum iteration steps is reached. The latter depends on the resources (human and financial resources in case of laboratory experiments or computational resources in case of in silico experiments). In every iteration the Gaussian process chooses a new point aiming to move to the optimum and at the same time to reduce model uncertainty.

Results of this optimization framework are the set of Pareto optimal setups (also called Pareto front) and their corresponding response values.

#### *2.5. Numerical Solvers and Tools*

The programming language and numeric computing environment MATLAB [37] was used for the seed train simulations. The code for the optimization workflow was written in Python [38] using the MATLAB Engine API for Python to call MATLAB as a computational engine from Python code. To perform Bayes optimization within this workflow, the library GPflow [39] was used.

#### **3. Results and Discussion**

#### *3.1. Optimization of Cultivation Vessels Regarding Number of Shake Flask Scales and Filling Volumes for Five, Four and Three Shake Flask Scales*

In this section, it is investigated which cultivation filling volumes should be used for the flask scales in order to obtain optimal results in terms of seed train duration and deviation rate, here defined as the probability that the seed train will run outside the predefined acceptable ranges for initial viable cell density (VCD) and transfer VCD (final VCD before transfer into the next cultivation vessel) per scale. The latter is a measure for the robustness of the seed train regarding viable cell density.

For assessment of the optimization results, a conventional reference seed train comprising five shake flask scales was simulated based on a non-optimized design. Therefore, a common passaging interval of 3 days per cultivation scale was fixed and filling volumes were determined following a conservative layout (i.e., choosing not too huge differences between one cultivation scale and the next to ensure that enough viable cells are generated even if they grow a little bit slower than expected).

In the first step, the optimal combination of filling volumes for five shake flask scales is investigated and the results are compared to the reference seed train. Afterwards, it is investigated if a reduction in shake flask scales from five to four or three shake flask scales leads to similar or even better results in terms of seed train duration and deviation rate. The number of bioreactor scales was kept fixed. Three bioreactors with filling volumes of 40 L, 320 L and 2100 L were used as pre-stages before inoculation of the production bioreactor with 9600 L. The assumed seed train setup is given in Table 1.

**Table 1.** Specification of the exemplary seed train setup providing information concerning cultivation vessels, required viable cell densities and the transfer of cells from one cultivation vessel into the next larger one, assumed in this work.

To find the optimal solution, multi-objective Bayesian optimization coupled with uncertainty-based seed train simulation, as described in Section 2.1, was applied. First, a Latin hypercube design for *n*lhs design 'points' (combinations of filling volumes, here *n*lhs = 10) was initiated and seed train simulation was applied to calculate the objective criteria values, here, deviation rate *D* and seed train duration *d* (replacing the normally required experimental cultivation runs) at each point. Within the Bayes optimization procedure, Gaussian processes were trained based on the simulation outcomes and an acquisition function was calculated in each iteration step in order to propose which point should be evaluated next. The input space for shake flask filling volumes (here the optimizable variables) was defined as described in Table 2, assuming the possibility of using several shake flasks in parallel for one shake flask scale and also considering their working volumes ranges.

**Table 2.** Input space for the shake flask filling volumes, containing the possible filling volumes per scale, given for optimization runs with 5, 4 or 3 shake flask scales.

#### 3.1.1. Optimization of Five Shake Flask Scales

The first optimization was performed for a seed train comprising five shake flask scales. Figure 3 shows the objective criteria values for each evaluated point, whereby the outcomes based on the initial Latin hypercube space are illustrated by blue dots and the outcomes for the proposed points based on the trained Gaussian processes are illustrated through yellow crosses. The optimal solutions are those near to the lower left corner aiming to minimize seed train duration and the deviation rate. The Pareto optimal solutions, also called non-dominated solutions, are illustrated through green circles. A solution (seed

train setup/combination of filling volumes) is called non-dominated if no solution exists leading to better (here lower) objective criteria values. As described previously, several Pareto optimal solutions can be obtained because when considering two or more objective criteria then for two different solutions one criterion might have better (here lower) value then the other solution for the same objective, while the other criterion has worse (here higher) values. The set of all Pareto optimal solutions is called Pareto front.

**Figure 3.** Algorthmically determined solutions and Pareto front of the two objective criteria seed train duration and deviation rate. (optimizable variables, here combinations of 5 shake flask filling volumes). Blue dots show an initial Latin hypercube design (LHC); yellow crosses are the points proposed by the algorithm; green circles are Pareto optimal solutions (=Pareto front).

For the investigated scenario (five shake flask scales and the seed train configuration according to Table 1) five Pareto optimal solutions were obtained (see green circles in Figure 3). It can be seen that comparing two of these solutions (green circles) each, one solution has a lower (here better) seed train duration value than the other solution and the opposite holds for the deviation rate.

The corresponding values for the optimizable variables, here shake flask filling volumes (*V*1, *V*2, *V*3, *V*4 and *V*5), and the corresponding objective criteria values, here deviation rate *D* and seed train duration *d*, are listed in Table 3.

**Table 3.** Pareto optimal solutions concerning the choice of filling volumes in shake flask scales, for three scenarios for 5 flask scales. The following bioreactor filling volumes are 40 L, 320 L and 2210 L. The averaged filling volumes in L and the resulting deviation rate (*D*) in % and seed train duration (*d*) in h are listed for each Pareto optimal solution.

The filling volume of the first scale was limited to a very narrow range (14–15 mL) (A higher variation after cell thawing was not expected). Most obtained solutions start with

the maximum value of this range (see Table 3, first column). The filling volume of flask scale 2 varies between 0.065 and 0.115 L, the filling volume of flask scale 3 between 0.340 and 0.904 L, the filling volume of flask scale 4 between 1.582 and 2.355 L and of flask scale 5 between 6.85 and 7.97 L. All five combinations lead to a deviation rate *D* of less than 7% and to a seed train duration between 520 to 537 h.

A more detailed illustration of the obtained results is presented in Figures 4 and 5. For two optimizable variables and one objective criterion each (deviation rate in Figure 4 and seed train duration in Figure 5), a contour plot is shown which illustrates the objective value for each calculated point (combination of the two variables), using the trained Gaussian processes, through colored isolines.

For example, the diagram in the top left of Figure 4 shows the deviation rate for each combination of *V*1 (filling volume in flask scale 1) and *V*2 (filling volume in flask scale 2) through colors representing the corresponding values in %, as indicated in the color bar. The results obtained through seed train simulations are shown by dots. The red dots represent the non-dominated (optimal solutions), optimal with respect to the defined multi-objective optimization problem. The dark blue area indicates combinations of *V*1 and *V*2 leading to a lower deviation rate. It can be seen that values above 0.1 for *V*2 combined with any value of *V*1 (within the given range) lead to the lowest deviation rates (below 6.2%, see dark blue area). Moreover, the optimal solutions (red dots) are mostly located in the area with higher filling volumes for shake flask 2, *V*2, except one (red dot at *V*2 ≈ 0.065).

For some combinations, a closer delimitation is possible. For example, the middle diagram in the second row (*V*3 over *V*2) shows a limited region (dark blue area) and therewith a specific combination of *V*3 and *V*2 that leads to the lowest deviation rates (<5.6%). These are around 0.3 L for *V*3 and around 0.105 L for *V*2. Furthermore, two optimal solutions (red dots) out of the set of Pareto optimal solutions (considering both objective criteria, seed train duration and deviation rate) are located in this region. The remaining red dots are located outside of the dark blue regions (see turquoise regions in the same diagram), meaning that they have higher deviation rates. Analogously, Figure 5 shows the contour plots for the second objective criterion, seed train duration. The dark blue areas show the combinations with the lowest seed train durations (approximately below 528 h). It can be seen in these diagrams that most red dots are located in the dark blue regions. For some combinations the dark blue areas are wider, distributed over several possible values for one variable, e.g., the diagram in the top center, top left, center, and center right.

Other combinations show narrower regions with low seed train durations as can be seen in the diagram showing *V*4 over *V*3. The lowest seed train duration is obtained for filling volumes between 1.5 and 2.5 L for shake flask 4 in combination with filling volumes between 0.2 and 0.8 L for shake flask 3.

Overall, these diagrams give an overview of the impact of two combined optimizable variables each on a specific objective criterion.

In addition to this information, simulated time profiles (predictive mean in green, 90% prediction bands in blue) of viable and total cell density as well as concentrations of glucose, glutamine, lactate and ammonium (see Figure 6) can be obtained for each solution, as well as a seed train protocol containing information about the calculated passaging intervals, amount of medium, etc.

**Figure 4.** Contour plots showing two optimizable variables on x and y-axis and one objective (here Deviation rate *D* in %), assigned to each combination of the two variables, through colored isolines. For example, the diagram in the top left shows the deviation rate for each combination of *V*1 (filling volume in flask scale 1) and *V*2 (filling volume in flask scale 2) through colors representing the corresponding values in %, as indicated on the color bar. Moreover, the results obtained through seed train simulations are shown by dots. The red dots represent the non-dominated (optimal solutions).

**Figure 5.** Contour plots showing two optimizable variables on x and y-axis and one objective (here seed train duration (*d*) in h), assigned to each combination of the two variables, through colored isolines. For example, the diagram top left shows the deviation rate for each combination of *V*1 (filling volume in flask scale 1) and *V*2 (filling volume in flask scale 2) through colors representing the corresponding values in %, as indicated on the color bar. Moreover, the results obtained through seed train simulations are shown by dots. The red dots represent the non-dominated (optimal solutions).

**Figure 6.** Seed train showing viable cell density (VCD) and total cell density, as well as substrate (glucose and glutamine) and metabolite (lactate and ammonium) concentrations over time and over the whole seed train (5 shake flask scales and three bioreactor scales), based on the shake flask filling volumes according to solution 1. The green lines represent the mean time course and the blue lines show the corresponding 90%-prediction band (5%- and 95%-quantiles). The plot (top left) also includes the filling volumes and the acceptable ranges for seeding VCD and transfer VCD, illustrated through dashed lines.

It can be seen in the top left of Figure 6 that based on the given filling volumes in addition to the flexibility to choose individual points in time for cell passaging in each scale, it is possible to set the seeding viable cell density at the beginning of each cultivation scale on the desired value with low variability, allowing to stay within the corresponding acceptable ranges for seeding VCD (see yellow dashed lines). Moreover, transfer VCDs lie within the corresponding acceptable range with high probability (see lower boundary, gray dashed line). Moreover, it can be seen that substrate concentrations are not depleted and according to [27], values of 20 mmol/L lactate and 5 mmol/L ammonium are not yet inhibiting concentrations for this cell line.

For a better assessment, the obtained results are compared to the reference seed train which is also defined in this work for five shake flask scales and illustrated in Figure 7. It grounds on a (non-optimized) configuration setup for five shake flask scales using fixed passaging intervals of 72 h each (common practice) and filling volumes of 15 mL (flask scale 1), 80 mL (flask scale 2), 300 mL (flask scale 3), 2000 mL (flask scale 4) and 4000 mL (flask scale 5). This choice grounds on a rather conservative approach aiming to avoid the risk of reaching too low transfer cell densities at the end of a cultivation scale but without the inclusion of probabilistic simulations.

The proposed method instead includes risk calculations and a passaging strategy aiming to minimize this risk but at the same time identifying a seed train configuration which is optimal regarding further objectives such as seed train duration in the present case.

A comparison of the seed train solutions obtained after optimization and the reference seed train shows that deviation rate is much lower after optimization (4.9–6.7% instead of 41.7%) and seed train duration could be reduced by 56 h from 576 h to 520 h. Figure 7, diagram top left shows where seeding or transfer viable cell density do not lie fully within the acceptable ranges (see red circles). This is different for the optimized solutions, e.g., solution 5, as illustrated in Figure 6, where seeding VCD lies within the acceptable range and also transfer VCD lies above the lower bound of the acceptable range for transfer VCD. This significant reduction in time (≈2 days per seed train) would contribute to a meaningful acceleration of the production process.

**Figure 7.** Reference (non-optimized) seed train showing viable and total cell density, as well as substrate (glucose and glutamine) and metabolite (lactate and ammonium) concentrations over time, based on a reference configuration setup for 5 shake flask scales using passaging intervals of 72 h each. The green lines represent the mean time course and the blue lines show the corresponding 90%-prediction band (5%- and 95%-quantiles).

#### 3.1.2. Optimization of Three and Four Shake Flask Scales

In the next step, the number of shake flask scales was reduced from five to four and then to three shake flask scales and the same optimization procedure was applied. The aim was to investigate if less cultivation vessels would lead to comparable results and if so, which target and filling volumes should be chosen. This is of interest because less operations (such as transferring cells from one scale into another one) signify less risk of failure and deviations.

Figure 8 shows the obtained values for the objective criteria deviation rate and seed train duration for different combinations of filling volumes for three (left) and for four shake flask scales (right). Furthermore, here, the solutions based on the initial Latin hypercube design are shown by blue dots and Pareto optimal solutions are highlighted through green circles.

It can be seen that for both scenarios, combinations of filling volumes could be found leading to an overall seed train cultivation time between 519 and 530 h. However, the scenario of using four shake flask scales, leads to lower deviation rates (*D* < 10%) compared to the scenario of using three shake flask scales (23% < *D* < 26%).

The corresponding filling volumes and the obtained filling volumes (based on the underlying passaging strategy) of the Pareto optimal solutions are listed in Table 4 together with the results for five shake flask scales from Table 3. The results are sorted as discovered by the optimization algorithm. The obtained filling volumes for four shake flask scales are very similar, except for shake flask scale 4 ( *V*1 = 15 mL, *V*2 = 158–200 mL, *V*3 = 1.51–1.60 L and *V*4 = 4.81–7.58 L). Some of the obtained solutions would be seen or treated as equal in practice, because the differences are rather small. For example it would not be distinguished between 0.190 and 0.195 L. Probably 200 mL would be used instead. However, the applied optimization algorithm works on a continuous input space and differentiates between the solutions listed in the Table 4, even though the differences are very low. The obtained optimal filling volumes for three shake flask scales also look similar, but with a bit more variation for shake flask 3 (4.45–5.58 L).

**Figure 8.** Algorithmically determined solutions and Pareto front regarding seed train duration and deviation rate for 3 resp. 4 shake flask scales on top resp. bottom. Blue dots show an initial Latin hypercube design (LHC); yellow crosses are the points proposed by the algorithm; green circles are Pareto optimal solutions (=Pareto front).

Comparing the results for the three scenarios (three, four and five shake flask scales) endorses a decision against the three flask scales-scenario due to the higher deviation rates (>20%), which stands for less process robustness. Between the other two scenarios (four or five shake flask scales) only little differences with respect to deviation rates are observed for the determined optimal solutions (4.9–6.7% for five shake flasks, 5.9–9.2% for four shake flasks). Using five shake flask scales would lead to more or less similar cultivations times (520–537 h) but one operational step more would be required.

This information, together with the corresponding seed train protocol, provides a solid basis to take a decision for one of the proposed optimal seed trains designs, taking into account seed train duration, robustness (expressed through deviation rates) and operational steps.

**Table 4.** Pareto optimal solutions concerning the choice of filling volumes in shake flask scales, for 3, 4 and 5 shake flask scales. The following bioreactor filling volumes are 40 L, 320 L and 2210 L. The averaged filling volumes in L, the resulting deviation rate (*D*) in % and seed train duration *d* in h are listed for each solution.

#### *3.2. Application to Further Cell Lines with Potentially Different Growth Rates*

The optimization examples presented in the previous subsection were applied to a specific CHO cell line with growth characteristics described by a set of model parameters derived from an industrial cell culture process which was investigated in [27]. If a different cell line or a clonal cell population with potentially differing growth behavior is used, then the optimization has to be performed for this specific cell line. In the following simulation study, a cell line having a 5% lower and a cell line having a 5% higher maximum cell-specific growth rate compared to the reference maximum growth rate (*μ*max = 0.028 h−1 for the first bioreactor scale and *μ*max = 0.029 h−1 for the remaining seed train scales) are assumed and the optimization is applied for both scenarios.

The results for the obtained/proposed filling volumes, as well as the corresponding seed train duration and deviation rate are listed in Table 5.

As expected, cells which grow faster (higher maximum growth rate *μ*max) would require less time until reaching a specific target cell density. This can be seen in the right column of Table 5. Using five flask scales, the optimal required seed train duration would lie between 494 and 503 h for a cell line with a 5% higher growth rate compared to the reference cell line which would need 520–537 h (see Table 3). Correspondingly, cells with a 5% lower growth rate would need more time (550–568 h). The same is observed when using four or three shake flasks.

With respect to the deviation rates which represent the robustness of the seed train design regarding variability of viable cells, it can be seen that low deviation rates of between 4.1% and 11.6% can be reached when using five or four flask scales, even if the maximum growth rate varies ±5%. A critical limit was identified for the combination of using three shake flask scales for a slower growing cell line. The corresponding optimal solution shows a comparatively higher deviation rate (19.2–29.1%) together with a high seed train duration (548–552 h).

**Table 5.** Pareto optimal solutions concerning the choice of filling volumes in shake flask scales, for 3, 4 and 5 shake flask scales, for two different scenarios. Scenario 1 assumes a 5% lower and scenario 2 a 5% higher cell-specific maximum growth rate compared to the reference maximum growth rate. The bioreactor filling volumes which follow after the shake flask scales are 40 L, 320 L and 2210 L. The filling volumes in L, the resulting deviation rate (*D*) in % and seed train duration (*d*) in h are listed for each solution (several Pareto optimal solutions can be obtained per setup).

#### *3.3. Optimization Regarding Four Objectives Including Product Concentration*

To show the applicability of the proposed method to more than two objectives, a third and a fourth objective criterion, titer concentration and viability at the end of the production vessel (after 8 days) was added. Whereas the first two objective criteria (seed train duration and deviation rate) are related to the seed train itself, the third and fourth criterion refer to the generated product in the production vessel and to the viability of the cells in the production vessel. Product concentration, as well as product quality can be influenced by many factors (seeding cell density, substrate concentrations and nutrient feeds, metabolite production, temperature, pH, dissolved oxygen and carbon dioxide concentration, osmolality and more) and also by the amount and the state of the cells at the end of the seed train. Since no data describing product quality are available, product concentration and viability are considered in this study. A further simplification that was made is the assumption that the production vessel is performed in batch-mode (meaning

without any addition of nutrient feeds or medium renewals). The reason for this simplification is to avoid confounding effects. The authors are aware of the fact that many factors affect product concentration and product quality and when data of other critical process parameters or quality attributes are available, these could also be considered in the same manner. The main purpose of the present simulation example is to demonstrate how the proposed method can be applied to more than two objectives and how the corresponding results can be illustrated and interpreted.

To obtain a visual overview for multiple objective criteria in one figure, a so-called spider plot (or net plot) can be used, which is shown in Figure 9.

**Figure 9.** Spider plot showing the objective criteria values (seed train duration, deviation rate, titer and viability after 8 days in the production vessel) for the Pareto optimal solutions for 5 shake flask scales.

The horizontal axis shows the values of the deviation rate (on the right) and of the viability (on the left). The vertical axis shows the values of the seed train duration (above) and of the titer (below). The aim of the optimization was to minimize seed train duration and deviation rate and to maximize viability and titer. Each color (hyperplane) represents one of the Pareto optimal seed train configurations (based on the optimal combinations of filling volumes in shake flask scales). Since seed train duration and deviation rate should be minimal and titer and viability should be maximal, hyperplanes covering the lower left area would be desired. However, no such solution (hyperplane) was obtained. The reason is that the optimization problem contains conflicting objective criteria, meaning that an improvement of one criterion leads to a degradation of another criterion. The here presented solutions are all non-dominated (see the green circles in the figures for two objective criteria). For all shown solutions, the deviation rate is rather low (4.9–7.3%), the seed train duration lies between 521 and 562 h and a titer of approximately 430–433 mg/L (assuming here a cell-specific production rate of *q*titer,max = 3.9 × 10−10 mg cell−1 h−1, as reported in [40]) and a viability of 52–53% is reached after 8 days in the production vessel (here via batch-mode). Of course, the obtained values depend a lot on the real process conditions (production bioreactor probably performed in fed-batch model) and the model parameter values obtained after model validation. However, the presented simulation example shall illustrate how the proposed approach can be applied for risk-based decision making under consideration of several criteria that should be optimal.

#### *3.4. Impact of Performed Iterations during Bayes Optimization*

For the example of three shake flask scales, (followed by three bioreactor scales) and optimizing filling volumes for all shake flask scales with respect to the two objective criteria: seed train duration and deviation rate, the number of performed iterations during the optimization procedure was varied. First, 10 initial points (combinations of filling volumes) distributed based on a Latin hypercube design were evaluated, followed by 10 Bayes iterations, which means that 10 times the algorithm updates the black box model (the Gaussian process), calculates the acquisition function and proposes the next point based on the outcome of this calculation. Then, the optimization was performed again for the same seed train setup but using 20 and then 30 Bayes iterations. The obtained solutions are illustrated in Figure 10.

**Figure 10.** Algorthmically determined solutions and Pareto front of the two objective criteria seed train duration and deviation rate (optimizable variables, here combinations of 3 shake flask filling volumes) for 10 (**top left**), 20 (**top right**) and 30 (**bottom left**) Bayes iterations; blue dots showthe initial Latin hypercube design (LHC); yellow crosses show the points proposed by the algorithm; green circles are the Pareto optimal solutions (=Pareto front).

Increasing the number of iterations from 10 to 20 helped to identify one solution that has not been discovered when running only 10 iterations. This can be seen when comparing the green circles in the diagram top left and the green circles in the diagram top right. The solution with *D* ≈ 22 and *d* = 523 cannot been found in the diagram top left.

Increasing the number of iterations from 20 to 30 did not lead to an improved optimum as can be seen when comparing the green circles in the Figure 10 top right diagram and bottom left diagram. This underlines the efficiency of the Bayes optimization. In the present example, only 10 initial points (distributed randomly according to a Latin hypercube design) and 20 Bayes optimization iteration steps were required to obtain the results which were confirmed when applying 30 iteration steps.

#### *3.5. Summary*

The objective of the first optimization problem was to design a robust seed train (cell expansion process), which means a seed train layout (including the number of cultivation scales, filling volumes and passaging intervals) leading to a reproducible seed train with

low variability regarding viable cell density and with a minimum seed train duration. The obtained solutions were compared to a non-optimized reference seed train and a comparison showed that the deviation rate is much lower after optimization (<10% instead of 41.7%) and seed train duration could be reduced by 56 h from 576 h to 520 h, which means a significant reduction of more than 2 days.

Addressing the question of if variation of the number of shake flask scales (and therewith the number of passaging steps) would lead to similar results in terms of deviation rates and seed train duration, it turned out that a reduction to three shake flask scales, would mean an increase in deviation rate and is therefore not recommended, at least under the assumed working volume ranges.

In industrial practice, typically more than one cell line is in use (different cell lines may be used to produce different molecules/products). Since growth rates of different cell lines differ, it was investigated how optimal seed train designs would differ for cell lines with 5% higher or lower growth rates. It turned out that the same optimization procedure could be easily adapted (by modification of the model parameter maximum growth rate) and applied to the modified setup revealing critical limits, e.g., for the combination of using three shake flask scales for a slower growing cell line. The latter shows comparatively high deviation rates (19.2–29.1%) together with high seed train durations (548–552 h instead of 519–523 h for the reference growth rate).

To show the applicability of the proposed method to more than two objective criteria, a third and fourth objective criterion, product concentration (titer) and viability after 8 days in the production phase, were added and the optimization was performed regarding four objective criteria in total. These are seed train duration, deviation rate (i.e., the probability that the seed train will run outside the predefined criteria), titer and viability at the end of the production phase.

Moreover, it was investigated for one seed train configuration (three shake flask scales and two objectives using the reference cell growth rate) if increasing the number of Bayes iterations would identify different optima. A number of 20 Bayes iterations turned out to be sufficient, because running 20 or 30 Bayes iterations showed similar results, which underlines the efficiency of the Bayes optimization approach.

In the present case study, the volumes are considered as fixed after optimization. If the production process allows for more flexibility in terms of adapting the volume within a specific range in the case that cells grow slower than the expected mean, then a reduction in the deviation rate can be achieved because varying the volumes allows for regulation of the inoculum viable cell density. However, this flexibility is not always given due to regulatory requirements and therefore not considered in the present study.

#### **4. Conclusions**

A concept has been developed to use process models in combination with algorithms for Bayes optimization using Gaussian processes to solve multi-objective optimization problems in the context of biopharmaceutical production processes. To illustrate this approach, a relevant exemplary optimization problem was chosen and solved using the proposed method.

The goal was to find optimal combinations of filling volumes for the shake flask scales of a seed train leading to a minimum deviation rate regarding viable cell densities and a minimum process duration. Compared to a non-optimized reference seed train, the optimized process showed much lower deviation rates regarding viable cell densities (<10% instead of 41.7%) using five or four shake flask scales and seed train duration could be reduced by 56 h from 576 h to 520 h.

Overall, it is shown that applying Bayes optimization to a multi-objective optimization function with several optimizable input variables and under a considerable amount of constraints, lead to revealing results with a low computational effort. This approach provides the potential to be used in form of a decision tool, e.g., for the choice of an optimal and robust seed train design but also to further optimization tasks within process development.

It should be noted that Bayes optimization and the corresponding computational modules could also be applied, even if no mechanistic process model is available, following a slightly different workflow. Instead of performing model-based in silico experiments (process simulations), real lab experiments would be performed and fed back to update the black box model (here the Gaussian process). This adaptive procedure (also called Bayesian experimental design or experimental design with Bayesian optimization [41]) or further related optimization methods might be promising tools to support experimental planning, process characterization, process transfer or optimization of cell culture processes but they still require further research and being embedded in software solutions that are easy to use for operators.

**Author Contributions:** Conceptualization, T.H.R., A.S., B.F. and M.L.-H.; methodology, T.H.R., A.S. and M.L.-H.; software, A.S. and T.H.R.; validation, A.S. and T.H.R.; formal analysis, T.H.R.; investigation, T.H, A.S., M.L.-H. and B.F.; resources, M.L.-H. and B.F.; data curation, T.H.R.; writing original draft preparation, T.H.R.; writing—review and editing, A.S., B.F. and M.L.-H.; visualization, T.H.R. and A.S.; supervision, M.L.-H. and B.F.; project administration, T.H.R., B.F. and M.L.-H.; funding acquisition, B.F. and T.H.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** The article processing charge (APC) was funded partially by Ostwestfalen-Lippe University of Applied Sciences and Arts (TH OWL).

**Acknowledgments:** The authors would like to express special thanks to Christoph Posch (Novartis Technical Research and Development) for the fruitful scientific exchange regarding the here presented case study. Moreover, we acknowledge support for the Open Access fees by Ostwestfalen-Lippe University of Applied Sciences and Arts (TH OWL) in the funding program Open Access Publishing.

**Conflicts of Interest:** All authors T.H.R., A.S., M.L.-H. and B.F. do not have any conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **Appendix A**

**Table A1.** Mechanistic model [27,42–45] for description of cell growth, cell death, substrate uptake, metabolite production and antibody production applicable to batch and fed-batch mode.

#### **Appendix B**

*Application to Other Cell Lines with Potentially Higher and Lower Maximum Growth Rates*

**Figure A1.** *Cont*.

**Figure A1.** *Cont*.

**Figure A1.** Pareto solutions for 3, 4 and 5 shake flask scales and for three different growth rates, reference maximum growth rate (left column), a 5% lower (middle column) and a 5% higher growth rates (right column) showing the objective criterion seed train duration over objective criterion deviation rate, using 20 optimization iterations; Blue dots: based on the initial Latin hypercube (LHC) design; Yellow crosses: based on the proposed points (by the algorithm); Green circles: Pareto optimal solutions. (**a**) 5 sf, *μ*,max,ref. (**b**) 5 sf, *μ*,max,95%. (**c**) 5 sf, *μ*,max,105%. (**d**) 4 sf, *μ*,max,ref. (**e**) 4 sf, *μ*,max,95%. (**f**) 4 sf, *μ*,max,105%. (**g**) 3 sf, *μ*,max,ref. (**h**) 3 sf, *μ*,max,95%. (**i**) 3 sf, *μ*,max,105%.

#### **References**

# *Article* **Considerations of the Impacts of Cell-Specific Growth and Production Rate on Clone Selection—A Simulation Study**

**Tanja Hernández Rodríguez 1, Sophie Morerod 2, Ralf Pörtner 3, Florian M. Wurm 2,4 and Björn Frahm 1,\***

**Abstract:** For the manufacturing of complex biopharmaceuticals using bioreactors with cultivated mammalian cells, high product concentration is an important objective. The phenotype of the cells in a reactor plays an important role. Are clonal cell populations showing high cell-specific growth rates more favorable than cell lines with higher cell-specific productivities or vice versa? Five clonal Chinese hamster ovary cell populations were analyzed based on the data of a 3-month-stability study. We adapted a mechanistic cell culture model to the experimental data of one such clonally derived cell population. Uncertainties and prior knowledge concerning model parameters were considered using Bayesian parameter estimations. This model was used then to define an inoculum train protocol. Based on this, we subsequently simulated the impacts of differences in growth rates (±10%) and production rates (±10% and ±50%) on the overall cultivation time, including making the inoculum train cultures; the final production phase, the volumetric titer in that bioreactor and the ratio of both, defined as overall process productivity. We showed thus unequivocally that growth rates have a higher impact (up to three times) on overall process productivity and for product output per year, whereas cells with higher productivity can potentially generate higher product concentrations in the production vessel.

**Keywords:** clonal cell population; phenotypic diversity; inoculum train; uncertainty-based; cell culture model; biopharmaceutical manufacturing

#### **1. Introduction**