Research Article Wei Xie, Department of Mechanical and Industrial Engineering,
Northeastern University, Boston, MA 02115, USA<EMAIL_ADDRESS>
# Interpretable Biomanufacturing Process Risk and Sensitivity Analyses for
Quality-by-Design and Stability Control

Wei Xie  Department of Mechanical and Industrial Engineering, Northeastern
University, Boston, MA 02115, USA  Bo Wang  Department of Mechanical and
Industrial Engineering, Northeastern University, Boston, MA 02115, USA  Cheng
Li  Department of Statistics and Applied Probability, National University of
Singapore, Singapore  Dongming Xie  Department of Chemical Engineering,
University of Massachusetts Lowell, Lowell, MA 01854, USA  Jared Auclair
Department of Chemistry and Chemical Biology, Northeastern University, Boston,
MA 02115, USA

###### Abstract

While biomanufacturing plays a significant role in supporting the economy and
ensuring public health, it faces critical challenges, including complexity,
high variability, lengthy lead time, and very limited process data, especially
for personalized new cell and gene biotherapeutics. Driven by these
challenges, we propose an interpretable semantic bioprocess probabilistic
knowledge graph and develop a game theory based risk and sensitivity analyses
for production process to facilitate quality-by-design and stability control.
Specifically, by exploring the causal relationships and interactions of
critical process parameters and quality attributes (CPPs/CQAs), we create a
Bayesian network based probabilistic knowledge graph characterizing the
complex causal interdependencies of all factors. Then, we introduce a Shapley
value based sensitivity analysis, which can correctly quantify the variation
contribution from each input factor on the outputs (i.e., productivity,
product quality). Since the bioprocess model coefficients are learned from
limited process observations, we derive the Bayesian posterior distribution to
quantify model uncertainty and further develop the Shapley value based
sensitivity analysis to evaluate the impact of estimation uncertainty from
each set of model coefficients. Therefore, the proposed bioprocess risk and
sensitivity analyses can identify the bottlenecks, guide the reliable process
specifications and the most informative data collection, and improve
production stability.

###### keywords:

Bioprocess risk analysis, sensitivity analysis, manufacturing process
stability control, Bayesian network, process causal interdependence

##  1 Introduction

Biomanufacturing is growing rapidly and playing an increasingly significant
role in supporting the economy and ensuring public health. For example, the
biopharmaceutical industry generated more than $300 billion in revenue in 2019
and more than 40% of the drug products in the development pipeline were
biopharmaceuticals [39]. However, drug shortages have occurred at
unprecedented rates over the past decade. The current systems are unable to
rapidly produce new drugs to meet urgent needs in the presence of a major
public health emergency. The COVID-19 pandemic is having a profound impact
globally and caused over 111 millions confirmed cases by February, 2021. Even
COVID-19 vaccines are discovered, developing the production process and
manufacturing the billions of doses needed to immunize the world’s population
will be extremely time-consuming using existing technologies, thus lengthening
the time period of human and economic distress. It is critically important to
speed up the bioprocess development and ensure product quality consistency.

However, biomanufacturing faces several critical challenges, including high
complexity and variability, and lengthy lead time [17]. Biomanufacturing is
based on living cells whose biological processes are very complex and have
highly variable outputs. The productivity and product critical quality
attributes (CQAs) are determined by the interactions of hundreds of critical
process parameters (CPPs), including raw materials, media compositions,
feeding strategy, and process operational conditions, such as pH and dissolved
oxygen in the bioreactor. As new biotherapeutics (e.g., cell and gene
therapies) become more and more “personalized", the production, regulation
procedure, and analytical testing time required by biopharmaceuticals of
complex molecular structure is lengthy, and the historical observations are
relatively limited in particular for drugs in early stages of production
process development.

Therefore, it is crucial to integrate all sources of data and mechanism
information, provide the risk- and science-based understanding of the complex
bioprocess CPPs/CQAs causal interdependencies, and identify and control the
key factors contributing the most to the output variation. This study can
accelerate the development of productive and reliable biomanufacturing,
facilitate building the quality into the production process or quality-by-
design (QbD), support real-time monitoring and release, and reduce the time to
market.

Various Process Analytical Technologies (PAT) and methodologies have been
proposed to improve the bioprocess understanding and guide the process
development, decision making, and risk control; see the review in [46]. Most
PATs are based on multivariate data analysis; see Section 2. Ordinary or
partial differential equations (ODEs/PDEs) based mechanistic models are
developed for simulating individual biomanufacturing unit operations; see for
example [20]. On the other hand, various operations research/management
(OR/OM) methods are also proposed for biomanufacturing system analytics and
decision-making; see the review [17]. Overall, existing methodologies have the
key limitations: (1) the multivariate statistics based PAT and OR/OM
approaches focus on developing general methodologies without incorporating the
bioprocess causal relationship and structural mechanism information, which
limits their performance, interpretability, and adoption, especially with
limited data; and (2) the mechanistic models are usually deterministic and
focus on individual unit operations without providing an reliable integrated
bioprocess learning and risk management framework.

Driven by the critical challenges in the biomanufacturing industry, in this
paper, we propose a bioprocess semantic probabilistic knowledge graph,
characterizing the risk- and science-based understanding of integrated
production process, which can integrate all sources of heterogeneous data and
leverage the information from existing mechanism models and historical data.
Then, we introduce comprehensive and rigorous bioprocess risk and sensitivity
analyses, accounting for model risk, which can guide the process
specifications and most informative data collection to facilitate the learning
and improve the production reliability and stability (e.g., product quality
consistency).

The key contributions of this paper are three fold. First, by exploring the
causal relationships and interactions of many factors within and between
operation units (i.e., CPPs/CQAs), such as raw materials, production process
parameters, and product quality, we consider a bioprocess ontology based data
integration and develop a Bayesian network (BN) based bioprocess probabilistic
knowledge graph, characterizing the process inherent stochastic uncertainty
and causal interdependencies of all input and output factors. Second, building
on the process knowledge graph, we introduce a game theory – Shapley value
(SV) – based sensitivity analysis (SA), considering the complex bioprocess
interdependencies, which can correctly quantify the contribution and
criticality of each random input factor on the variance of outputs (i.e.,
productivity and product CQAs), identify the bottlenecks, and accelerate the
reliable bioprocess specifications. Third, since the coefficients of
interpretable bioprocess model or probabilistic knowledge graph are estimated
from limited real-world process data, which induces model uncertainty (MU) or
model risk (MR), we further propose Bayesian uncertainty quantification and
Shapley value based model uncertainty sensitivity analysis to support process
learning and faithfully assess the impact of estimation uncertainty from each
set of model coefficients. Thus, our study can: (1) identify the bottlenecks
of bioprocess; (2) accelerate the reliable process specifications and
development to improve the production process stability and facilitate QbD;
and (3) support the most “informative" data collection to reduce the model
risk of process probabilistic knowledge graph and improve the bioprocess
understanding.

This paper is organized as follows. In Section 2, we review the related
literature on biomanufacturing process modeling and PATs, Bayesian network,
and process sensitivity analysis. In Section 3, we present the problem
description and summarize the proposed framework. In Section 4, we develop the
Bayesian network (BN) based bioprocess probabilistic knowledge graph to
characterize the risk- and science-based process understanding. We derive the
Shapley value (SV) based bioprocess sensitivity analysis in Section 5 to
support the process specifications, improve the production stability, and
ensure product quality consistency. We further introduce the process model
coefficient uncertainty quantification and Shapley value based sensitivity
analysis studying the impact of each model coefficient estimation uncertainty
on process risk analysis and CPPs/CQAs criticality assessment in Section 6. We
conduct the empirical study on the performance of our proposed framework in
both simulation and real data analysis in Section 7, and then conclude with
some discussion in Section 8.

##  2 Background

The Process Analytical Technologies (PAT) are defined as “a system for
designing, analyzing and controlling manufacturing through timely measurements
of critical quality and performance attributes of raw and in-process materials
and processes, with the goal of ensuring final product quality"; see [36].
With the established process sensors and analyzers, such as near infrared
spectroscopy, Raman spectrocopy, and multiwavelength fluorescence, various
multivariate data analysis approaches have been used for bioprocess PATs,
including principal component analysis (PCA) [1], partial least squares (PLS)
[7], clustering [37], multilinear regression [53], artificial neural network
(ANN) [24], genetic algorithm [44], elastic net [41], support vector machines
[23], and root cause analysis [2]; see an overview in [40]. However, existing
PAT approaches are usually based on generalized multivariate “black-box"
approaches quantifying the input-output relationship without incorporating the
bioprocess mechanism information.

On the other hand, OR/OM methodology development for biomanufacturing analysis
and decision making is still in its infancy [17]. Mixed integer linear
programming [21, 22], dynamic lot size model [11], and queueing network and
simulation models [25, 19] have been developed to study resource planning,
scheduling and material consumption in biomanufacturing. Those approaches
focus on developing general methodologies without fully exploring the bio-
technology domain knowledge (e.g., causal relationship, structural information
of the bioprocess). Some recent works, e.g., [30, 29, 28], account for
physical-chemical characteristics and biology-induced randomness in either
fermentation or chromatography stage, and develop Markov decision models to
optimize the corresponding operational policies.

For complex systems, Bayesian network (BN) can be used to combine the expert
knowledge with data and facilitate data integration and process analysis in
various applications. For example, [52] proposed a BN based knowledge
management system for additive manufacturing. [48] introduced a BN that
combines evidence from gene co-expression and experimental data to predict
whether two genes are functionally related. [32] provided a BN approach for
system architecture generation and evaluation, and [47] applied probabilistic
graphical model to study how the usage context factors, including human
factors, situational factors, and product design factors, impact on the energy
consumption of the lightweight vehicle to guide usage scenarios and vehicle
designs. Furthermore, Bayesian posterior and belief propagation based risk
assessment has been studied in information system security [10], water mains
failure [16], and supply chain [33]. Motivated by these studies, we propose a
Bayesian network for modeling the complex interdependence of production
process parameters and bio-drug properties, which can fully utilize the
structural knowledge and causal relationship, and integrate the data from end-
to-end bioprocess.

Finally, we briefly discuss the related literature on sensitivity analysis;
see the review [3]. The existing sensitivity analysis studies associated with
Bayesian network tend to systematically vary one of network’s parameter at a
time while fixing the other parameters and then obtain analytic expressions
for the sensitivity functions [49, 4]. In our case, we are interested in
stochastic uncertainty contributed by each factor, which is closely related to
global probabilistic sensitivity analysis. Existing literature on global
sensitivity analysis can be divided into several categories, including: (1)
regression based methods, e.g., [15], which use the standardized regression
coefficients as sensitivity measure; (2) variance based methods [50, 43] which
assess the contribution of each random input based on expected reduction in
model output variance; (3) functional ANOVA decomposition [38] which provides
variance decomposition under independence through high dimensional model
representation theory; (4) density-based methods [54] that directly quantify
the output density without reference to a particular moment. Since the
commonly used variance-based sensitivity measures (i.e., first-order effects
and total effects) fail to adequately account for probabilistic dependence of
inputs and process structural interactions or interdependencies, [35]
introduced a new sensitivity measure based on the game theory, called the
Shapley Value (SV). [45] further analyzed this measure and proposed a Monte
Carlo algorithm for the estimation of Shapley values. [27] proposed Shapley
value based unified framework for interpreting predictions. Inspired by these
studies, building on the proposed BN-based bioprocess knowledge graph
characterizing the process causal interdependencies, we introduce SV-based
probabilistic sensitivity analysis to assess the contribution or criticality
of each random input (e.g., CPP and CQA ) on the output variance, while
accounting for the impact of model estimation uncertainty associated with each
set of model coefficients.

##  3 Problem Description and Proposed Framework

We create a probabilistic graph model characterizing the risk- and science-
based understanding of causal interdependencies between bioprocess CPPs/CQAs,
and then propose risk and sensitivity analyses for integrated biomanufacturing
process, accounting for model uncertainty. This study can: (1) provide a
reliable guidance on process specification, CPPs/CQAs monitoring, and most
informative data collection; (2) facilitate production stability control and
quality-by-design (QbD); and (3) accelerate real-time release, speed up the
time to market, and reduce the drug shortage.

An illustration of biomanufacturing process is provided in Fig. 1 with a fish
bone representation of bioprocess input factors introduced in each unit
operation impacting on the outputs. The biomanufacturing process typically has
several main unit operations, including: (1) media preparation, (2) inoculum
fermentation, (3) main fermentation, (4) centrifugation(s), (5)
chromatography/purification, (6) filtration, (7) fill and finish, and (8)
quality control. Steps (1)–(3) belong to upstream cell culture, Steps (4)–(6)
belong to downstream target protein purification, and Steps (7)–(8) are for
finished drug filling/formulation and final product quality control testing.

The interactions of many factors impact the variability of outputs (e.g., drug
quality, productivity). They can be divided into CPPs and CQAs in general; see
the definitions of CPPs/CQAs in ICH Q8(R2) [13].

  * CPP:

At each process unit operation, CPPs are defined as critical process
parameters whose variability impacts on product CQAs, and therefore should be
monitored and controlled to ensure the process produces the desired quality.

  * CQA:

A physical, chemical, biological, or microbiological property that should be
within an appropriate limit, range, or distribution to ensure the desired
product quality.

Since the raw material attributes are outputs of release materials, they
should be considered along with CPPs as impacting process variability.

Figure 1: An illustration of general biomanufacturing process and fish-bone
representation [51].

We represent the system output (e.g., product CQAs, productivity) with a
random variable, denoted by $Y$, which depends on CPPs/CQAs inputs, denoted by
$\mathbf{X}$, and other uncontrolled/uncontrollable input variables (e.g.,
contamination), modeled by residuals $\mathbf{e}$. We represent the impact of
complex interactions of input factors $(\mathbf{X},\mathbf{e})$ throughout the
production process on the response by
${Y}=g(\mathbf{X},\mathbf{e}|\boldsymbol{\theta})$, where the unknown function
$g(\mathbf{X},\mathbf{e}|\boldsymbol{\theta})$, specified by model
coefficients $\boldsymbol{\theta}$, models the complex interactions of
integrated bioprocess and characterizes the impact of random inputs
$(\mathbf{X},\mathbf{e})$ on the output $Y$. For notation simplification, we
consider the unit-variate response/output in the paper, and the proposed
framework can be naturally extended to a vector of responses.

To provide the risk- and science-based production process understanding and
guide the reliable process development, we need to correctly quantify all
sources of uncertainties. There are two types of uncertainty: (1) bioprocess
inherent stochastic uncertainty from CPPs/CQAs and other uncontrolled
variables (i.e., randomness of $\mathbf{X}$ and $\mathbf{e}$), which can be
reduced by the identification of missed CPPs and tighter specification of
selected CPPs; and (2) model uncertainty (MU) (i.e., the estimation
uncertainty of bioprocess model coefficients $\boldsymbol{\theta}$), which can
be reduced by collecting “most informative" process observations. Correctly
quantifying all sources of uncertainty can facilitate learning, guide risk
elimination/control, and improve robust, automatic, and reliable bioprocess
decision making.

###  3.1 Review of Game Theory based Sensitivity Measure - Shapley Value

In game theory, the Shapley value (SV) was originally introduced to evaluate
the contribution of a player in a cooperative game [42]. A cooperative game is
defined as a set of players $\mathcal{K}=\\{1,2,\ldots,K\\}$, with a function
$c(\cdot)$ that maps a subset of players to its corresponding payoff,
$c:2^{\mathcal{K}}\to\mathbb{R}$ with $c(\emptyset)=0$, where
$2^{\mathcal{K}}$ denotes the power set of $\mathcal{K}$ (i.e., the set of all
subsets of $\mathcal{K}$). Thus, $c(\mathcal{J})$ characterizes the total gain
that the players in subset $\mathcal{J}\subset\mathcal{K}$ can obtain by
cooperation. The SV of player $k\in\mathcal{K}$ with respect to $c(\cdot)$ is
defined by

$\mbox{Sh}_{k}=\sum_{\mathcal{J}\subset\mathcal{K}/\\{k\\}}\dfrac{(K-|\mathcal{J}|-1)!|\mathcal{J}|!}{K!}\left[c(\mathcal{J}\cup\\{k\\})-c(\mathcal{J})\right],$
(1)

where $K=|\mathcal{K}|$ is the total number of players and $|\mathcal{J}|$ is
the size of subset $\mathcal{J}$ from $\mathcal{K}/\\{k\\}$. This SV can be
interpreted as the average incremental payoff by including player $k$ over all
possible cooperation group formations, i.e.,
$\mathcal{J}\subset\mathcal{K}/\\{k\\}$, and $\mbox{Sh}_{k}$ can be used to
measure the contribution of the player $k$. This assessment approach satisfies
the “efficiency property" that the sum of the SVs of all players equals the
gain of the grand coalition, i.e.,
$c(\mathcal{K})=\sum_{k=1}^{K}\mbox{Sh}_{k}$.

The Shapley value was recently introduced for global sensitivity analysis to
measure the variance of output contributed by each random input [35]. Denote
the set of inputs as
$\mathbf{U}_{\mathcal{K}}=\\{U_{1},U_{2},\ldots,U_{K}\\}$, and model the
output $V=\eta(\mathbf{U}_{\mathcal{K}})$ as a function $\eta(\cdot)$ of the
inputs, accounting for their interactions. Two most commonly used variance-
based sensitivity measures are: (1) the first-order effect
$O_{k}\equiv\mbox{Var}(V)-\mbox{E}[\mbox{Var}(V|U_{k})]$ that considers the
variance reduction when we fix $U_{k}$; and (2) the total effect
$T_{k}\equiv\mbox{E}[\mbox{Var}(V|\mathbf{U}_{-k})]$ that considers the
expected remaining variance when all other factors, denoted by
$\mathbf{U}_{-k}$, are fixed. However, both measures fail to appropriately
quantify the sensitivity or variance contribution when there exist
probabilistic interdependence among inputs and process structural interaction
[45].

Built on the SV from game theory, given a cooperative game with inputs
$\mathbf{U}_{\mathcal{K}}$ as the players and the payoff as the incremental
variance in output $V$ induced by any index subset
$\mathcal{J}\subset\mathcal{K}$, one can define the payoff function as

${c(\mathcal{J})=\mbox{Var}(V)-\mbox{E}[\mbox{Var}[V|\boldsymbol{U}_{\mathcal{J}}]]~{}~{}\mbox{or}~{}~{}}c(\mathcal{J})=\mbox{E}[\mbox{Var}[V|\boldsymbol{U}_{-\mathcal{J}}]].$
(2)

Thus, [35] introduced a new SV-based sensitivity measure, with
$\mbox{Sh}_{U_{k},V}$ computed by Equations (1) and (2). In this paper, we use
$c(\mathcal{J})=\mbox{E}[\mbox{Var}[V|\boldsymbol{U}_{-\mathcal{J}}]]$ in
Equation (1), which can simplify the computation of the contribution from any
random input $U_{k}$ on the output variance $\mbox{Var}(V)$,
$\mbox{Sh}_{U_{k},V}=\mbox{Sh}_{k}$. The SV-based sensitivity analysis
overcomes the limitations of first-order effect and total effect measures by
accounting for the interdependence of inputs and process interactions. The
variance of output $V$ can be decomposed into the contribution from each
random input $U_{k}$ and we can define the criticality as the proportion of
$\mbox{Var}(V)$ contributed from $U_{k}$, denoted by $p_{U_{k},V}$,

$\mbox{Var}(V)=\sum_{k=1}^{K}\mbox{Sh}_{U_{k},V}~{}~{}\mbox{ and
}~{}~{}p_{U_{k},V}=\frac{\mbox{Sh}_{U_{k},V}}{\mbox{Var}(V)}.$

The main benefits of SV over first-order and total effect sensitivity measures
include: (1) the uncertainty contributions sum up to total variance of output;
and (2) SV can automatically account for probabilistic dependence and
structural interactions occurring in the complex production process.

###  3.2 Summary of Proposed Interpretable Bioprocess Model, Risk and
Sensitivity Analyses for Integrated Bioprocess Stability Control

Fig. 2 provides the flowchart of proposed risk and sensitivity analyses
framework, which can accelerate learning of the end-to-end production process
and guide the development of stable biomanufacturing. Parts I and II focus on
modeling and reducing of process stochastic uncertainty. Part III focuses on
analyzing and controlling the model risk. By exploring the causal
relationships and interactions of CPPs/CQAs of raw materials/in-process
materials/product within and between different process modules, in Section 4,
we develop ontology based data integration and create an interpretable
Bayesian network (BN) based bioprocess semantic probabilistic knowledge graph,
specified by the model coefficients $\boldsymbol{\theta}$. This knowledge
graph can characterize the risk- and science-based understanding of integrated
bioprocess and quantify the causal interdependencies of inputs
$(\mathbf{X},\mathbf{e})$ and output $Y$. It is interpretable and extendable,
which can support flexible process modular design, incorporate the existing
mechanisms from different modules and operation units, quantify the bioprocess
causal interdependencies, and greatly reduce the dimensionality of bioprocess
design space to guide the decision making.

Figure 2: The flowchart of proposed biomanufacturing process risk and
sensitivity analyses framework.

Building on this interpretable probabilistic knowledge graph, in Section 5, we
develop the SV-based process risk and sensitivity analyses studying stochastic
uncertainty, and derive variance decomposition to quantify the contribution
from each random input,

$\mbox{Var}(Y|\boldsymbol{\theta})=\sum_{X_{k}}\mbox{Sh}_{X_{k},Y}(\boldsymbol{\theta})+\sum_{e_{k}}\mbox{Sh}_{e_{k},Y}(\boldsymbol{\theta}),$

where Shapley values, $\mbox{Sh}_{X_{k},Y}(\boldsymbol{\theta})$ and
$\mbox{Sh}_{e_{k},Y}(\boldsymbol{\theta})$, measure the contributions from any
CPP/CQA, $X_{k}\in\mathbf{X}$, and residual factor, $e_{k}\in\mathbf{e}$
(representing the impact of remaining uncontrolled factors on the CQA
$X_{k}$), to the output variance $\mbox{Var}(Y|\boldsymbol{\theta})$. For any
input factor $W_{k}$ (i.e., either $X_{k}$ or $e_{k}$), the criticality,
$p_{W_{k},Y}(\boldsymbol{\theta})\equiv\mbox{Sh}_{W_{k},Y}(\boldsymbol{\theta})/\mbox{Var}(Y|\boldsymbol{\theta})$,
can be used to identify the bottlenecks that contribute the most to
$\mbox{Var}(Y)$, and guide the process specifications to efficiently improve
production process stability. The CPPs/CQAs $X_{k}$ with high criticality
requires more restrict stability control, while the residual $e_{k}$, with
high impact on the output variance, can guide us to identify the ignored CPPs.
Since the Shapley value (SV) based sensitivity analysis is developed based on
game theory, its combination with bioprocess probabilistic knowledge graph,
accounting for the complex causal interdependencies, can correctly assess the
risk effect from each set of random input factors on the output variation.

The “correct" process model coefficients, denoted by
$\boldsymbol{\theta}^{c}$, characterizing the bioprocess underlying
probabilistic interdependence, is unknown and estimated by using the real-
world process data, denoted by $\mathcal{X}$. Given limited historical process
data, the model risk or estimation uncertainty can have a large impact on the
bioprocess risk and sensitivity analyses. Since the estimation uncertainty of
model coefficients at different parts of bioprocess can be different and
interdependent, the model uncertainty (MU) is quantified with the joint
posterior distribution $p(\boldsymbol{\theta}|\mathcal{X})$. We further
develop the SV-based sensitivity analysis to study the impact of model
uncertainty from each part of process in Section 6, which can guide the “most
informative" data collection to reduce the impact from model risk and
efficiently improve the accuracy of bioprocess risk analysis and critiality
assessment, especially for those factors contributing the most to the output
variance. For any random input $W_{k}$, the Shapley value
$\mbox{Sh}_{W_{k},Y}$ is estimated with error, which can be contributed by the
model coefficients located along the paths propagating the uncertainty of
$W_{k}$ to the output $Y$, denoted by $\boldsymbol{\theta}(W_{k},Y)$. We
introduce the BN-SV-MU sensitivity analysis to provide the comprehensive study
over the impact of model uncertainty,

$\mbox{Var}^{*}[\mbox{Sh}_{W_{k},Y}|\mathcal{X}]=\sum_{\theta_{\ell}\in\boldsymbol{\theta}(W_{k},Y)}\mbox{Sh}^{*}_{\theta_{\ell}}\left[\left.\mbox{Sh}_{W_{k},Y}\left(\widetilde{\boldsymbol{\theta}}(W_{k},Y)\right)\right|\mathcal{X}\right]=\sum_{\theta_{\ell}\in\boldsymbol{\theta}(W_{k},Y)}\mbox{Sh}^{*}_{\theta_{\ell}}\left[\left.\mbox{Sh}_{W_{k},Y}\right|\mathcal{X}\right],$
(3)

where the subscript “$*$" represents any measure calculated based on the
posterior $p(\boldsymbol{\theta}|\mathcal{X})$ and
$\mbox{Sh}^{*}_{\theta_{\ell}}\left[\left.\cdot\right|\mathcal{X}\right]$
measures the contribution from coefficient estimation uncertainty of
$\theta_{\ell}\in\boldsymbol{\theta}(W_{k},Y)$. In the proposed interpretable
bioprocess model, $\theta_{\ell}$ can be interpreted as certain mechanistic
coefficients (e.g., cell growth rate in the cell culture). Thus, the
decomposition in (3) provides the detailed information on how the model
uncertainty of each part of integrated production process influences the
estimation uncertainty of $\mbox{Sh}_{W_{k},Y}$.

To illustrate the key ideas of the proposed bioprocess risk and sensitivity
analyses, we use a simplified monoclonal antibody (mAbs) drug substance
production example, including main fermentation, centrifuge, chromatograph,
and filtration; see the interactions in Fig. 3. We consider the dominant
CPPs/CQAs in each step, while the impacts of remaining factors are included in
$\mathbf{e}$. We are interested in the variance contribution (or criticality)
from each CPP to drug substance protein content $Y=X_{20}$, and also account
for the impact of model uncertainty on criticality assessment.

Figure 3: An example illustration visualization of the integrated bioprocess
sensitivity framework for criticality assessment and model uncertainty.

The results of risk and sensitivity analyses can be visualized along the graph
model for this bioprocess example; see Fig. 3. The process knowledge graph
model is specified by coefficients,
$\boldsymbol{\theta}=(\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta})$,
where $\boldsymbol{\mu}$ and $\boldsymbol{v}^{2}$ are the mean and variance
vectors of all factors (listed in the table), and $\boldsymbol{\beta}$
quantifies the effects from parents’ nodes on each child node. The darkness of
nodes indicates the criticality level of process factors, which can guide the
better bioprocess specification. The results show that $X_{4}$ contributes to
$55\%$ variance of $X_{20}$. In addition, the darkness of directed edges and
circle boundaries indicates the distribution of model uncertainty in the
process. For instance, the variance of dissolved oxygen in bioreactor,
$v_{4}^{2}$, has dominant model uncertainty impact on the estimation
uncertainty of criticality $p_{X_{4},Y}$, which suggests additional data
should be collected to improve the estimation of $v_{4}^{2}$. We will revisit
this example and use it to study the performance of proposed framework in
Section 7.1.

##  4 Interpretable Bioprocess Probabilistic Model Development

We develop an interpretable bioprocess probabilistic model, which can be
extendable to end-to-end biomanfacturing supply chain and support evidence-
based biopharmaceutical production process development. The model can
incorporate the existing mechanism models from each module and unit operation,
and facilitate the learning from distributed and heterogeneous process data.
In this section, we develop a Bayesian network (BN) based interpretable
bioprocess probabistic knowledge graph, which can characterize the complex
CPPs/CQAs causal interdependencies and support flexible modular
biomanfuacturing development.

By exploring the causal relationships and interactions, we consider bioprocess
ontology-based data integration, which can connect all distributed and
heterogeneous data collected from bioprocess; see more description in Online
Appendix A. This relational graph can enable the connectivity of end-to-end
process. Nodes represent factors (i.e., CPPs/CQAs, media feed, bioreactor
operating conditions, other uncontrolled factors) impacting the process
outputs, and the directed edges model the causal relationships. Within each
module, which can be each phase of cell culture process (such as cell growth
and production phases) or each unit operation, we can model complex
interactions, e.g., biological/physical/chemical interactions. In the
relational graph, the shaded nodes represent the variables with real-world
observations, including the testing and sensor monitoring data of CPPs/CQAs
for raw materials, operation conditions, and intermediate/final drug products.
The unshaded nodes represent variables without observations and residuals,
including quality status of intermediate and final drug products, and other
uncontrollable factors (e.g., contamination) introduced during the process
unit operations. Since bio-products have very complex structures, we cannot
observe the underlying complete quality status and the monitoring of CQAs can
carry partial information. Building on the bioprocess relational graph, we
develop a BN based probabilistic graphical model composing of random
CPPs/CQAs/residuals factors and their conditional dependencies via directed
edges. It can characterize the probabilistic causal interdependencies among
all factors of integrated bioprocess.

To make it easy to follow, we first provide a simple illustration example of
cell culture process, including two phases, to present the key ideas of
modular bioprocess modeling, and then develop the complete probabilistic
knowledge graph model for general integrated bioprocess. Specifically, we use
a simple bioreactor fermentation example with two phases (i.e., cell growth
and production phases; see Fig. 4) to illustrate the probabilistic graphical
model development. It is based on the causal relationships and interactions
between CPPs and CQAs. Each node represents a CPP/CQA with a random variable
$X$ modeling its variability. Each directed edge represents the causal impact
of parent node $X_{i}$ on child node $X_{j}$. The pattern-fill nodes
($X_{1},X_{2},X_{3}$) represent the CPPs. The solid fill nodes ($X_{6},X_{7}$)
represent the monitored CQAs of intermediate materials and drug products. The
nodes $X_{4}$ and $X_{5}$ represent the underlying status of working cells
after cell growth phase and the protein/impurity structure after cell
production phase. The CQAs $X_{6}$ and $X_{7}$ represent the partial
information of quality variables $X_{4}$ and $X_{5}$. Except the CPPs
$X_{1},X_{2},X_{3}$, the impacts from other uncontrolled factors introduced
during two phases of cell culture are modeled through $e_{4}^{\prime}$ and
$e_{5}^{\prime}$.

Figure 4: Left: knowledge relational graph; Right: simplified knowledge graph.

Since it is hard to uniquely specify the underlying cells/proteins with very
complex structures, $X_{4}$ and $X_{5}$ are hidden, which can lead to an
identification issue. Typically, the non-identifiable BN with hidden nodes is
transformed to the equivalent BN by structural simplification to avoid
analytical issues (see Chapter 19 in [18]). Thus, we simplify and transform
the relational graphical model to a graph without hidden nodes, depicted in
the right panel of Fig. 4. The new residual $e_{6}$ in updated graph accounts
for both original residual $e^{\prime}_{4}$ and also the uncertainty of
underlying cell health status, $X_{4}$, impacting on CQA $X_{6}$, similar for
new residual $e_{7}$. According to the right plot in Fig. 4, the sources of
bioprocess stochastic uncertainty impacting on the variability of $X_{7}$
include CPPs, $(X_{1},X_{2},X_{3})$, and other factors with the impact
represented by residuals $(e_{6},e_{7})$. Thus, we have CPPs $X_{1},X_{2}$ as
inputs and CQA $X_{6}$ as output for the first cell growth phase, and have CQA
$X_{6}$ and CPP $X_{3}$ as inputs and $X_{7}$ as output for the second protein
production phase. To study the impact of each CPP on the CQA of interest
(i.e., $X_{6}$ and $X_{7}$), we can decompose the variance of $X_{6}$ and
$X_{7}$ into the contributions from $X_{1}$, $X_{2}$ and $X_{3}$, and
remaining parts coming from $e_{6}$ and $e_{7}$; see the process risk and
sensitivity analyses in Section 5. In this way, we can identify the main
sources of uncertainty and quantify their impacts, which can guide the
CPPs/CQAs specifications and the quality control to improve the product
quality stability.

Now we describe the BN-based bioprocess model for general situations. Suppose
that the integrated bioprocess can be represented by a probabilistic graphical
model with $m+1$ nodes: $m$ process factors (denoted by $\mathbf{X}$) and a
single response, denoted by $Y$, such as the impurity concentration or protein
content. Let the first $m^{p}$ nodes representing CPPs
$\mathbf{X}^{p}=\\{X_{1},X_{2},\ldots,X_{m^{p}}\\}$, the next $m^{a}$ nodes
representing CQAs $\mathbf{X}^{a}=\\{X_{m^{p}+1},X_{m^{p}+2},\ldots,X_{m}\\}$,
and the last node representing the response $Y\triangleq X_{m+1}$ with
$m=m^{p}+m^{a}$. The modular bioprocess probabilistic knowledge graph can be
modeled by marginal and conditional distributions of each node as follows:

$\displaystyle\small X_{k}$ $\displaystyle\sim$
$\displaystyle\mathcal{N}(\mu_{k},v_{k}^{2})\mbox{ for CPP $X_{k}$ with
$k=1,2,\ldots,m^{p}$,}$ (4) $\displaystyle X_{k}$ $\displaystyle=$
$\displaystyle f(Pa({X}_{k});\boldsymbol{\theta}_{k})+e_{k}\mbox{ for CQA
$X_{k}$ with $k=m^{p}+1,\ldots,m+1$}$ (5)

where $\mathcal{N}(u,v^{2})$ denotes the normal distribution with mean $u$ and
variance $v^{2}$, and $Pa(X_{k})$ denotes the parent nodes of $X_{k}$. By
applying central limit theory (CLT), we assume that the residual
$e_{k}\sim\mathcal{N}(0,v_{k}^{2})$ with the conditional variance
$v_{k}^{2}\equiv\mbox{Var}[X_{k}|Pa(X_{k})]$. Since the amount of real-world
bioprocess batch data is often very limited, Gaussian distribution is used to
model the variability of each variable or node, which is often used in the
existing biopharmaceutical studies (see for example [6]). It also makes the
process risk and sensitivity analyses tractable.

The proposed probabilistic knowledge graph is a hybrid model of integrated
bioprocess, which can leverage the existing mechanisms and learn from real-
world process data. Basically, the prior of the function $f(\cdot)$ in a
generalized regression model (5) can be specified based on the existing
knowledge on underlying bioprocess mechanisms (e.g., biophysicochemical
kinetics) within each module of bioprocess; see for example [20, 26, 8]. The
unknown model coefficients $\boldsymbol{\theta}_{k}$ (e.g., cell growth rate,
media consumption rate) need be estimated from process data. In this paper, we
consider linear function accounting for the main effects, i.e.,

$X_{k}=\mu_{k}+\sum_{X_{j}\in Pa(X_{k})}\beta_{jk}(X_{j}-\mu_{j})+e_{k}\mbox{
for CQA $X_{k}$ with $k=m^{p}+1,\ldots,m+1$}$ (6)

where the coefficient $\beta_{jk}$ can be used to measure the effect from the
parent node $X_{j}$ to child node $X_{k}$.

Here, we use some illustrative examples to briefly show how the proposed
Bayesian network based process probabilistic model allows us to incorporate
the existing bioprocess mechanisms. We first consider the cell exponential
growth mechanism for the fermentation step, $x=x_{0}e^{\mu t}$, where $x_{0}$
and $x$ denote the starting and ending cell densities, and $\mu$ is the
unknown growth rate. This is a commonly used mechanism model in
biomanufacturing industry; see more information in [9]. Suppose that there is
a fixed cell culture duration $t$. By doing the log transformation and setting
$X_{k}=\log(x_{0})$, $X_{k+1}=\log(x)$ and $\beta_{0}=\mu t$, we can take the
exponential growth mechanism as prior and get the hybrid probabilistic model
for the exponential growth phase in fermentation or cell culture process,
$X_{k+1}=\beta_{0}+X_{k}+e_{k}$, where $e_{k}$ represents the residual term
characterizing the integrated effect from many other factors and it follows a
Gaussian distribution by following CLT. Notice that it is a special case of
BN-based process model (6). The similar idea can be applied to the situations
where we have PDE/ODE-based bioprocess kinetics mechanism models,

$\frac{d}{dt}x(t)=f(x(t);\theta_{t})\approx\frac{x(t_{k+1})-x(t_{k})}{t_{k+1}-t_{k}}=f(x(t_{k});\theta_{t_{k}}),$

where $x(t)$ can represent the concentrations of protein and metabolite waste
at time $t$ and $\theta_{t}$ can denote the nonstationary growth rate. We can
take the existing mechanism model as the prior knowledge of production
process. By applying the finite difference on the gradient $dx(t)/dt$ and
first-order Taylor approximation on function $f(\cdot)$, we can construct a
probabilistic hybrid model matching with the formula in Equation (6), which
can leverage the information from existing PDE/ODE-based bioprocess kinetics
mechanism models. This approximation can be very accurate if the data are
collected from real-time production process sensor monitoring with high
sampling frequency.

The complex bioprocess CPPs/CQAs causal interdependencies are characterized by
the BN-based probabilistic knowledge graph. Given the model parameters
$\boldsymbol{\theta}=(\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta})$
with mean $\boldsymbol{\mu}=(\mu_{1},\ldots,\mu_{m+1})^{\top}$, conditional
variance $\boldsymbol{v}^{2}=(v_{1}^{2},\ldots,v_{m+1}^{2})^{\top}$, and
linear coefficients
$\boldsymbol{\beta}=\\{\beta_{jk};k=m^{p}+1,\ldots,m+1\mbox{ and }X_{j}\in
Pa(X_{k})\\}$, the conditional distribution for each CQA node $X_{k}$ becomes,

$p(X_{k}|{Pa(X_{k})})=\mathcal{N}\left(\mu_{k}+\sum_{X_{j}\in
Pa(X_{k})}\beta_{jk}(X_{j}-\mu_{j}),v_{k}^{2}\right)\mbox{ for
}k=m^{p}+1\ldots,m+1.$

For any CPP node $X_{k}$ without parent nodes, $Pa(X_{k})$ is an empty set and
$P(X_{k}|Pa(X_{k}))$ is just the marginal distribution $P(X_{k})$ in (4).
Therefore, the joint distribution characterizing the interdependencies of CPPs
and CQAs involved in the production process can be written as
$p(X_{1},X_{2},\ldots,X_{m+1})=\prod_{k=1}^{m+1}p(X_{k}|{Pa(X_{k})})$.

##  5 Process Risk and Sensitivity Analyses

Given the bioprocess probabilistic knowledge graph specified by the
coefficients
$\boldsymbol{\theta}=(\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta})$,
we develop the BN-SV based sensitivity analysis for integrated production
process and quantify the criticality of each random input factor measuring its
contribution to the output variance $\mbox{Var}(X_{m+1})$. This study can
guide the CPPs/CQAs specification and improve the production process
stability. To make it easy to follow, we start with a simple illustration
example, provide the general process risk and sensitivity analyses, and then
present the algorithm at the end of this section.

We again use the simple example in Fig. 4 to illustrate the results of
proposed BN-SV based bioprocess risk and sensitivity analyses, which can
decompose the output variance of protein/impurity concentration $X_{7}$ after
fermentation process to each random input – including $X_{1}$, $X_{2}$,
$X_{3}$, $e_{6}$, $e_{7}$ – as

$Var(X_{7}|\boldsymbol{\theta})=\underbrace{(\beta_{16}\beta_{67})^{2}v_{1}^{2}}_{\text{contribution
from
$X_{1}$}}+\underbrace{(\beta_{26}\beta_{67})^{2}v_{2}^{2}}_{\text{contribution
from $X_{2}$}}+\underbrace{\beta_{37}^{2}v_{3}^{2}}_{\text{contribution from
$X_{3}$}}+\underbrace{\beta_{67}^{2}v_{6}^{2}}_{\text{contribution from
$e_{6}$}}+\underbrace{v_{7}^{2}}_{\text{contribution from $e_{7}$}}.$

(7)

The variance contribution from each random input, denoted by $W_{k}$ (i.e.,
$X_{1},X_{2},X_{3},e_{6},e_{7}$), depends on its variance $v^{2}_{k}$ and the
product of coefficients $\boldsymbol{\beta}$ located along the paths
propagating the uncertainty from $W_{k}$ to the output $X_{7}$; see Fig. 5.
The darker blue filled node (i.e., cell growth phase CPP $X_{2}$, feed rate)
contributes more to the output variance and has higher criticality. Thus, to
efficiently reduce the output variance, it requires more restrictive stability
control. The high impact of $e_{7}$ (with darker color) can guide us to
identify unrecognized or missed CPPs.

Figure 5: A simple example to illustrate BN-SV based process risk and
sensitivity analysis.

This simple example illustrates that the proposed production process BN-SV
risk and sensitivity analyses and the CPPs/CQAs criticality assessment are
based on the bioprocess probabilistic knowledge graph, characterizing the
complex CPPs/CQAs causal interdependencies and accounting for all sources of
process inherent uncertainty, which can (1) guide the process specifications;
(2) improve the product quality consistency and bioprocess stability; and (3)
advance the risk- and science-based understanding on bioprocess.

Now we present the general process risk and sensitivity analyses. We first
derive the Shapley value (SV) quantifying the contribution of each random
input factor from CPPs $\mathbf{X}^{p}$ and other factors $\mathbf{e}$ to
$\mbox{Var}(X_{m+1})$, which accounts for cases with dependent input factors.
According to the Gaussian BN model presented in (4) and (6), we can write

$X_{m+1}=\mu_{m+1}+\sum_{k=1}^{m^{p}}\gamma_{k,m+1}(X_{k}-\mu_{k})+\sum_{k=m^{p}+1}^{m+1}\gamma_{k,m+1}e_{k},$
(8)

where the weight coefficient of any CPP $X_{k}$ to CQA $X_{n}$ with $k\leq
m^{p}<n\leq m+1$,

$\gamma_{kn}=\beta_{kn}+\sum_{m^{p}<\ell<n}\beta_{k\ell}\beta_{\ell
n}+\sum_{m^{p}<\ell_{1}<\ell_{2}<n}\beta_{k\ell_{1}}\beta_{\ell_{1}\ell_{2}}\beta_{\ell_{2}n}+\ldots+\beta_{k,m^{p}+1}\beta_{m^{p}+1,m^{p}+2}\ldots\beta_{n-1,n},$
(9)

the weight coefficient of any $e_{k}$ to a CQA node $X_{n}$ with
$m^{p}<k<n\leq m+1$,

$\gamma_{kn}=\beta_{kn}+\sum_{k<\ell<n}\beta_{k\ell}\beta_{\ell
n}+\sum_{k<\ell_{1}<\ell_{2}<n}\beta_{k\ell_{1}}\beta_{\ell_{1}\ell_{2}}\beta_{\ell_{2}n}+\ldots+\beta_{k,k+1}\beta_{k+1,k+2}\ldots\beta_{n-1,n};$
(10)

and $\gamma_{nn}=1$ for any $n$; see the derivation for (8) in Appendix B. The
weight coefficient $\gamma_{kn}$ is the product sum of $\boldsymbol{\beta}$
located along the paths from node $X_{k}$ to node $X_{n}$ in the graph model.
Let
$\boldsymbol{W}=\\{X_{1},\ldots,X_{m^{p}},e_{m^{p}+1},\ldots,e_{m+1}\\}\triangleq\\{W_{1},W_{2},\ldots,W_{m+1}\\}$
represent all random input factors, with the index set
$\mathcal{K}=\\{1,2,\ldots,m+1\\}$. Then, the SV for the $k$-th factor $W_{k}$
is,

$\mbox{Sh}_{W_{k},X_{m+1}}=\sum_{\mathcal{J}\subset\mathcal{K}/\\{k\\}}\dfrac{(m-|\mathcal{J}|)!|\mathcal{J}|!}{(m+1)!}\left[c(\mathcal{J}\cup\\{k\\})-c(\mathcal{J})\right].$

Based on (8), we compute the cost function,

$c(\mathcal{J})=\mbox{E}[\mbox{Var}[X_{i}|\boldsymbol{W}_{-\mathcal{J}}]]=\sum_{k\in\mathcal{J}}\gamma_{k,m+1}^{2}\mbox{Var}(W_{k})+2\sum_{k_{1}<k_{2}\in\mathcal{J}}\gamma_{k_{1},m+1}\gamma_{k_{2},m+1}\mbox{Cov}(W_{k_{1}},W_{k_{2}}).$

The random input factors,
$\boldsymbol{W}=(X_{1},\ldots,X_{m^{p}},e_{m^{p}+1},\ldots,e_{m+1})$,
including CPPs and residual terms introduced at each CQA nodes, are often
independent as the real biomanufacturing process specification is often based
on each CPP or CQA. To make the proposed framework general, we consider the
potential interdependence between some inputs $W_{k_{1}}$ and $W_{k_{2}}$ with
$k_{1}\neq k_{2}$, and the covariance $\mbox{Cov}(W_{k_{1}},W_{k_{2}})$ can be
estimated by using the process data.

Then, for each $W_{k}$ and $\mathcal{J}\subset\mathcal{K}/\\{k\\}$, we can
obtain

$c(\mathcal{J}\cup\\{k\\})-c(\mathcal{J})=\gamma_{k,m+1}^{2}\mbox{Var}(W_{k})+2\sum_{\ell\in\mathcal{J}}\gamma_{k,m+1}\gamma_{\ell,m+1}\mbox{Cov}(W_{k},W_{\ell}).$

Given the BN-based bioprocess knowledge graph model parameters
$\boldsymbol{\theta}$, by applying (1), we can derive the Shapley value,
$\mbox{Sh}_{W_{k},X_{m+1}}(\boldsymbol{\theta})$, characterizing the
contribution from any input factor $W_{k}$ to the output variance,

$\displaystyle\mbox{Sh}_{W_{k},X_{m+1}}(\boldsymbol{\theta})=\gamma_{k,m+1}^{2}\mbox{Var}(W_{k})+\sum_{\ell\neq
k}\gamma_{k,m+1}\gamma_{\ell,m+1}\mbox{Cov}(W_{k},W_{\ell}).$ (11)

The derivation of (11) is provided in Appendix C. Therefore, we can decompose
the variance of output $X_{m+1}$ and estimate the contribution from each
random input from $\mathbf{X}^{p}$ and $\mathbf{e}$,

$\mbox{Var}(X_{m+1}|\boldsymbol{\theta})=\sum\mbox{Sh}_{W_{k},X_{m+1}}(\boldsymbol{\theta})=\sum_{k=1}^{m^{p}}\mbox{Sh}_{X_{k},X_{m+1}}(\boldsymbol{\theta})+\sum_{k=m^{p}+1}^{m+1}\mbox{Sh}_{e_{k},X_{m+1}}(\boldsymbol{\theta}).$
(12)

Equation (12) can be used to identify the dominant factors in $\mathbf{X}^{p}$
and $\mathbf{e}$ contributing the most to the output variance, which can guide
the CPPs identification and process specification to improve the process
stability and quality consistency. As a result, the criticality of any input
factor $W_{k}$ can be calculated as
$p_{W_{k},X_{m+1}}(\boldsymbol{\theta})\equiv{\mbox{Sh}_{W_{k},X_{m+1}}(\boldsymbol{\theta})}/{\mbox{Var}(X_{m+1}|\boldsymbol{\theta})}$.
Notice that for any independent input factor $W_{k}$, the SV in Equation (11)
is reduced to
$\mbox{Sh}_{W_{k},X_{m+1}}(\boldsymbol{\theta})=\gamma_{k,m+1}^{2}v_{k}^{2}$.
Under the case that all input factors $\boldsymbol{W}$ are mutually
independent, the variance decomposition Equation (12) can be written as
$\mbox{Var}(X_{m+1}|\boldsymbol{\theta})=\sum\gamma_{k,m+1}^{2}v_{k}^{2}$,
which gives the example results in Equation (7).

This risk and sensitivity analyses can be applied to any part of production
process including one or multiple modules. Under this situation, the input
factors $\boldsymbol{W}$ include those nodes without parent node within the
considered range of production (i.e., CPPs, CQAs or uncontrolled factors), and
output of interest $X_{i}$ is certain CQA at the end of the procedure. For
example, in Fig. 5, we consider the subgraph, including
$\\{X_{3},X_{6},X_{7}\\}$, for cell production phase with the starting CQA
$X_{6}$ carrying the information from previous cell growth phase. We can study
the impacts of $X_{6}$ and CPP $X_{3}$ on the variability of CQA $X_{7}$. The
SV of any input $W_{k}$ and the variance decomposition of $X_{i}$, still
follow Equations (11) and (12) by replacing the output $X_{m+1}$ with $X_{i}$.
The criticality of $W_{k}$ on $X_{i}$ can be measured by proportion
$p_{W_{k},X_{i}}(\boldsymbol{\theta})=\mbox{Sh}_{W_{k},X_{i}}(\boldsymbol{\theta})/\mbox{Var}(X_{i}|\boldsymbol{\theta})$.

Given the BN parameters $\boldsymbol{\theta}$, we summarize the procedure for
production process BN-SV based sensitivity analysis in Algorithm 1, in which
we consider several consecutive operation steps, and our objective is to
quantify the contribution of each random factor in $\boldsymbol{W}$ to
$\mbox{Var}(X_{i}|\boldsymbol{\theta})$.

Input: BN parameters $\boldsymbol{\theta}$, group of input factors
$\boldsymbol{W}$, response node $X_{m+1}$.

Output: Variance decomposition of $X_{i}$ in terms of all random inputs within
$\boldsymbol{W}$.

(1) Calculate the Shapley value
$\mbox{Sh}_{W_{k},X_{m+1}}(\boldsymbol{\theta})$ by using Equation (11), which
measures the contribution from $W_{k}$ to the variance of response CQA
$X_{m+1}$;

(2) Provide the variance decomposition of
$\mbox{Var}(X_{m+1}|\boldsymbol{\theta})$ by using Equation (12), and obtain
the criticality of $W_{k}$ on the variance of $X_{m+1}$:
$p_{W_{k},X_{m+1}}(\boldsymbol{\theta})=\mbox{Sh}_{W_{k},X_{m+1}}(\boldsymbol{\theta})/\mbox{Var}(X_{m+1}|\boldsymbol{\theta})$.

Algorithm 1 Procedure for Production Process BN-SV based Sensitivity Analysis

##  6 Sensitivity Analysis for Model Risk Reduction

Since the underlying true process model coefficients $\boldsymbol{\theta}^{c}$
are unknown, given finite real-world data $\mathcal{X}$, there exists the
model uncertainty (MU) characterizing our limited knowledge on the
probabilistic interdependence of integrated bioprocess. To study the impact of
MU on the production process risk and sensitivity analyses for stochastic
uncertainty and further assess CPPs/CQAs criticality, we propose the BN-SV-MU
based uncertainty quantification and sensitivity analysis, which can guide the
process monitoring and “most informative" data collection. In Section 6.1, we
develop the posterior $p(\boldsymbol{\theta}|\mathcal{X})$ and a Gibbs sampler
to generate posterior samples, $\widetilde{\boldsymbol{\theta}}^{(b)}\sim
p(\boldsymbol{\theta}|\mathcal{X})$ with $b=1,2,\ldots,B$, quantifying the
model uncertainty, and then we quantify the overall impact of model
uncertainty on the process risk analysis and CPPs/CQAs criticality assessment.
In Section 6.2, we propose the BN-SV-MU based sensitivity analysis, which can
study the impact of each model coefficient(s) estimation uncertainty on the
process risk analysis and criticality assessment; see the result visualization
in Fig. 3.

###  6.1 Bayesian Learning and Model Uncertainty Quantification

We consider the case with $R$ batches of complete production process data,
denoted as
$\mathcal{X}=\\{(x_{1}^{(r)},x_{2}^{(r)},\ldots,x_{m+1}^{(r)}),r=1,2,\ldots,R\\}$.
Without strong prior information, we consider the following conjugate (vague)
prior (with initial hyperparameters giving relatively flat density),

$p(\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta})=\prod_{i=1}^{m+1}p(\mu_{i})p(v_{i}^{2})\cdot\prod_{i\neq
j}p(\beta_{ij}),$ (13)

with
$p(\mu_{i})=\mathcal{N}(\mu_{i}^{(0)},\sigma_{i}^{(0)2}),p(v_{i}^{2})=\mbox{Inv-}\Gamma\left(\dfrac{\kappa_{i}^{(0)}}{2},\dfrac{\lambda_{i}^{(0)}}{2}\right)$
and $p(\beta_{ij})=\mathcal{N}(\theta_{ij}^{(0)},\tau_{ij}^{(0)2})$, where
$\mbox{Inv-}\Gamma$ denotes the inverse-gamma distribution. Given the data
$\mathcal{X}$, by applying the Bayes’ rule, we can obtain the posterior
distribution

$p(\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta}|\mathcal{X})\propto\prod_{r=1}^{R}\left[\prod_{i=1}^{m+1}p(x_{i}^{(r)}|x_{Pa(X_{i})}^{(r)})\right]p(\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta}),$
(14)

quantifying the model uncertainty.

Then, we develop a Gibbs sampler to generate the posterior samples from (14)
quantifying the model uncertainty. We derive the conditional posterior for
each parameter in $(\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta})$.
Let $\boldsymbol{\mu}_{-i}$, $\boldsymbol{v}_{-i}^{2}$ and
$\boldsymbol{\beta}_{-ij}$ denote the collection of parameters
$\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta}$ excluding the $i$-th
or ${(i,j)}$-th element. Let ${S}(X_{i})$ denote the set of direct succeeding
or child nodes of node $X_{i}$. We first derive the conditional posterior for
the coefficient $\beta_{ij}$,

$\displaystyle
p(\beta_{ij}|\mathcal{X},\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta}_{-ij})=\mathcal{N}(\theta_{ij}^{(R)},\tau_{ij}^{(R)2}),$
(15)

where
$\theta_{ij}^{(R)}=\dfrac{\tau_{ij}^{(0)2}\sum_{r=1}^{R}\alpha_{i}^{(r)}m_{ij}^{(r)}+v_{j}^{2}\theta_{ij}^{(0)}}{\tau_{ij}^{(0)2}\sum_{r=1}^{R}\alpha_{i}^{(r)2}+v_{j}^{2}}\quad\textrm{and}\quad\tau_{ij}^{(R)2}=\dfrac{\tau_{ij}^{(0)2}v_{j}^{2}}{\tau_{ij}^{(0)2}\sum_{r=1}^{R}\alpha_{i}^{(r)2}+v_{j}^{2}}$
with $\alpha_{i}^{(r)}=x_{i}^{(r)}-\mu_{i}$ and
$m_{ij}^{(r)}=(x_{j}^{(r)}-\mu_{j})-\sum_{X_{k}\in
Pa(X_{j})/\\{X_{i}\\}}\beta_{kj}(x_{k}^{(r)}-\mu_{k}).$ Then, we derive the
conditional posterior for $v_{i}^{2}=\mbox{Var}[X_{i}|Pa(X_{i})]$ with
$i=1,2,\ldots,m+1$,

$\displaystyle
p(v_{i}^{2}|\mathcal{X},\boldsymbol{\mu},\boldsymbol{v}_{-i}^{2},\boldsymbol{\beta})=\mbox{Inv-}\Gamma\left(\dfrac{\kappa_{i}^{(R)}}{2},\dfrac{\lambda_{i}^{(R)}}{2}\right),$
(16)

where $\kappa_{i}^{(R)}=\kappa_{i}^{(0)}+R$,
$\lambda_{i}^{(R)}=\lambda_{i}^{(0)}+\sum_{r=1}^{R}u_{i}^{(r)2}$ and
$u_{i}^{(r)}=(x_{i}^{(r)}-\mu_{i})-\sum_{X_{k}\in
Pa(X_{i})}\beta_{ki}(x_{k}^{(r)}-\mu_{k}).$ After that, we derive the
conditional posterior for the mean parameter $\mu_{i}$ with $i=1,2,\ldots,m+1$
for any CPP/CQA,

$\displaystyle
p(\mu_{i}|\mathcal{X},\boldsymbol{\mu}_{-i},\boldsymbol{v}^{2},\boldsymbol{\beta})\propto
p(\mu_{i})\prod_{r=1}^{R}\left[p(x_{i}^{(r)}|x_{Pa(X_{i})}^{(r)})\prod_{j\in\mathcal{S}(X_{i})}p(x_{j}^{(r)}|x_{Pa(X_{j})}^{(r)})\right]=\mathcal{N}(\mu_{i}^{(R)},\sigma_{i}^{(R)2}),$
(17)

where
$\mu_{i}^{(R)}=\sigma_{i}^{(R)2}\left[\dfrac{\mu_{i}^{(0)}}{\sigma_{i}^{(0)2}}+\sum_{r=1}^{R}\dfrac{a_{i}^{(r)}}{v_{i}^{2}}+\sum_{r=1}^{R}\sum_{X_{j}\in
S(X_{i})}\dfrac{\beta_{ij}c_{ij}^{(r)}}{v_{j}^{2}}\right]$ and
$\dfrac{1}{\sigma_{i}^{(R)2}}=\dfrac{1}{\sigma_{i}^{(0)2}}+\dfrac{R}{v_{i}^{2}}+\sum_{X_{j}\in
S(X_{i})}\dfrac{R\beta_{ij}^{2}}{v_{j}^{2}}$ with
$a_{i}^{(r)}=x_{i}^{(r)}-\sum_{X_{k}\in
Pa(X_{i})}\beta_{kj}(x_{k}^{(r)}-\mu_{k})$ and
$c_{ij}^{(r)}=\beta_{ij}x_{i}^{(r)}-(x_{j}^{(r)}-\mu_{j})+\sum_{X_{k}\in
Pa(X_{j})/\\{X_{i}\\}}\beta_{kj}(x_{k}^{(r)}-\mu_{k}).$ The Gibbs sampler
iteratively draws the posterior samples of
$(\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta})$ by applying the
conditional posterior distributions given in (15), (16), and (17) until
convergence [12].

Besides the case with complete production data, we often have additional
incomplete batch data. Since the lead time for biopharmaceutical production is
lengthy [34], we can have some batches in the middle of production. In
addition, the bio-drug quality requirements are restricted, especially for
human drugs. Following the quality control, we could discard some batches
after main fermentation or even in the middle of downstream purification.
Thus, we provide the Gibbs sampler (see Algorithm 3) for both cases with
complete or mixing data in Appendix D.2.

Next, we study the impact model uncertainty on the bioprocess risk and
sensitivity analyses and CPPs/CQAs criticality assessment. Based on Section 5,
the contribution from any random input factor $W_{k}$ to the output variance
$\mbox{Var}(X_{m+1})$ is measured by the Shapley value,
$\mbox{Sh}_{W_{k},X_{m+1}}(\boldsymbol{\theta}^{c})$. The unknown parameters
$\boldsymbol{\theta}^{c}$ specifying the underlying process probabilistic
model are estimated by using limited real-world data $\mathcal{X}$. Thus, the
estimation uncertainty of the contribution from factor $W_{k}$ can be
quantified by the posterior distribution,
$\mbox{Sh}_{W_{k},X_{m+1}}(\widetilde{\boldsymbol{\theta}})$ with
$\widetilde{\boldsymbol{\theta}}\sim p(\boldsymbol{\theta}|\mathcal{X})$. We
can use the posterior mean to estimate the expected variance contribution and
criticality,
$\mbox{E}^{*}[\mbox{Sh}_{W_{k},X_{m+1}}|\mathcal{X}]\equiv\mbox{E}^{*}_{p(\boldsymbol{\theta}|\mathcal{X})}[\mbox{Sh}_{W_{k},X_{m+1}}(\widetilde{\boldsymbol{\theta}})|\mathcal{X}]$
and
$\mbox{E}^{*}[p_{W_{k},X_{m+1}}|\mathcal{X}]\equiv\mbox{E}^{*}_{p(\boldsymbol{\theta}|\mathcal{X})}[p_{W_{k},X_{m+1}}(\widetilde{\boldsymbol{\theta}})|\mathcal{X}]$,
where
$p_{W_{k},X_{m+1}}(\widetilde{\boldsymbol{\theta}})={\mbox{Sh}_{W_{k},X_{m+1}}(\widetilde{\boldsymbol{\theta}})}/{\mbox{Var}(X_{m+1}|\widetilde{\boldsymbol{\theta}})}$.
The posterior variance is used to quantify the overall estimation uncertainty
induced by model uncertainty,
$\mbox{Var}^{*}[\mbox{Sh}_{W_{k},X_{m+1}}|\mathcal{X}]\equiv\mbox{Var}^{*}_{p(\boldsymbol{\theta}|\mathcal{X})}[\mbox{Sh}_{W_{k},X_{m+1}}(\widetilde{\boldsymbol{\theta}})|\mathcal{X}]$
and
$\mbox{Var}^{*}[p_{W_{k},X_{m+1}}|\mathcal{X}]\equiv\mbox{Var}^{*}_{p(\boldsymbol{\theta}|\mathcal{X})}[p_{W_{k},X_{m+1}}(\widetilde{\boldsymbol{\theta}})|\mathcal{X}].$

Since we do not have the closed form solutions, we can estimate the posterior
mean and variance of $\mbox{Sh}(W_{k})$ and $p_{W_{k},X_{m+1}}$ through the
sampling approach. By applying the Gibbs sampler in Appendix D, we can
generate posterior samples $\widetilde{\boldsymbol{\theta}}^{(b)}\sim
p(\boldsymbol{\theta}|\mathcal{X})$ with $b=1,2,\ldots,B$. At any
$\widetilde{\boldsymbol{\theta}}^{(b)}$, we can compute
$\mbox{Sh}_{W_{k},X_{m+1}}(\widetilde{\boldsymbol{\theta}}^{(b)})$ following
the description in Section 5. The expected contribution from $W_{k}$ to the
variance of $X_{m+1}$ is estimated by
$\widehat{\mbox{E}}^{*}[\mbox{Sh}_{W_{k},X_{m+1}}|\mathcal{X}]=\bar{\mbox{Sh}}_{W_{k},X_{m+1}}(\mathcal{X})=\dfrac{1}{B}\sum_{b=1}^{B}\mbox{Sh}_{W_{k},X_{m+1}}(\widetilde{\boldsymbol{\theta}}^{(b)}).$
And the overall estimation uncertainty can be estimated by sample variance,

$\widehat{\mbox{Var}}^{*}[\mbox{Sh}_{W_{k},X_{m+1}}|\mathcal{X}]=\dfrac{1}{B-1}\sum_{b=1}^{B}\left[\mbox{Sh}_{W_{k},X_{m+1}}(\widetilde{\boldsymbol{\theta}}^{(b)})-\bar{\mbox{Sh}}_{W_{k},X_{m+1}}(\mathcal{X})\right]^{2}.$
(18)

Similarly, we can estimate the expected criticality by
$\widehat{\mbox{E}}^{*}[p_{W_{k},X_{m+1}}|\mathcal{X}]=\bar{p}_{W_{k},X_{m+1}}=\dfrac{1}{B}\sum_{b=1}^{B}p_{W_{k},X_{m+1}}(\widetilde{\boldsymbol{\theta}}^{(b)})$
and estimate the overall estimation uncertainty by

$\widehat{\mbox{Var}}^{*}[p_{W_{k},X_{m+1}}|\mathcal{X}]=\dfrac{1}{B-1}\sum_{b=1}^{B}\left[p_{W_{k},X_{m+1}}(\widetilde{\boldsymbol{\theta}}^{(b)})-\bar{p}_{W_{k},X_{m+1}}\right]^{2}.$
(19)

###  6.2 Sensitivity Study for Model Uncertainty

Since there is often limited process data in biomanufacturing, model
uncertainty tends to be large. We propose the BN-SV-MU based sensitivity
analysis studying the effect of estimation uncertainty of each model
coefficient, which can guide the process monitoring and “most informative"
data collection. We provide the CPPs/CQAs criticality estimation uncertainty
quantification and BN-SV-MU based sensitivity analysis in Algorithm 2.
Specifically, Steps (1)–(3) evaluate
${\mbox{Var}}^{*}[\mbox{Sh}_{W_{k},X_{m+1}}|\mathcal{X}]$ quantifying the
overall estimation uncertainty of $\mbox{Sh}_{W_{k},X_{m+1}}$. Steps (4)–(13)
further study the impact from each model coefficient estimation uncertainty.

Here we use $\mbox{Var}^{*}[\mbox{Sh}_{W_{k},X_{m+1}}|\mathcal{X}]$ for
illustration and the similar procedure can be applied to CPPs/CQAs criticality
assessment $\mbox{Var}^{*}[p_{W_{k},X_{m+1}}|\mathcal{X}]$. Let
$\boldsymbol{\theta}(W_{k},X_{m+1})\subset\boldsymbol{\theta}$ represent the
subset of model coefficients that impacts on
$\mbox{Sh}(W_{k}|\boldsymbol{\theta})$ estimation. Notice that
$\boldsymbol{\mu}$ has no impact on
$\mbox{Sh}_{W_{k},X_{m+1}}(\boldsymbol{\theta})$. Since SV can account for the
probabilistic dependence of model coefficient estimation uncertainty,
characterized by the joint posterior distribution
$p(\boldsymbol{\theta}|\mathcal{X})$, and bioprocess structural interactions,
we can measure the contribution from any parameter
$\theta_{\ell}\in\boldsymbol{\theta}(W_{k},X_{m+1})$ through the posterior
variance decomposition,

$\mbox{Var}^{*}[\mbox{Sh}_{W_{k},X_{m+1}}|\mathcal{X}]=\sum_{{\theta_{\ell}\in\boldsymbol{\theta}(W_{k},X_{m+1})}}\mbox{Sh}^{*}_{\theta_{\ell}}\left[\left.\mbox{Sh}_{W_{k},X_{m+1}}\left(\widetilde{\boldsymbol{\theta}}\right)\right|\mathcal{X}\right]=\sum_{\theta_{\ell}\in\boldsymbol{\theta}(W_{k},X_{m+1})}\mbox{Sh}^{*}_{\theta_{\ell}}\left[\left.\mbox{Sh}_{W_{k},X_{m+1}}\right|\mathcal{X}\right].$

The proposed BN-SV-MU sensitivity analysis can provide the comprehensive and
interpretable understanding on how model uncertainty impacts on the process
risk analysis and identify those parameters $\theta_{\ell}$ contributing the
most on the estimation uncertainty of
$\mbox{Sh}_{W_{k},X_{m+1}}(\boldsymbol{\theta})$.

Input: BN structure $G(\mathbf{N}|\boldsymbol{\theta})$, data $\mathcal{X}$,
number of samples $N_{\pi}$, $B$, $B_{O}$ and $B_{I}$, index subset
$\mathcal{L}_{k}$.

Output: Return
$\widehat{\widehat{\mbox{Sh}}}^{*}_{\theta_{\ell}}\left[\left.\mbox{Sh}_{W_{k},X_{m+1}}\right|\mathcal{X}\right]$
and
$\widehat{\mbox{Sh}}_{\theta_{\ell}}^{*}\left[\left.p_{W_{k},X_{m+1}}\right|\mathcal{X}\right]$
for any
$W_{k}\in\boldsymbol{W}=\\{\mathbf{X}^{p}\cup\mathbf{X}^{a}\cup\mathbf{e}\\}$.

(1) Call Algorithm 3 in Appendix D.3 to obtain the posterior samples
$\widetilde{\boldsymbol{\theta}}^{(b)}=(\widetilde{\boldsymbol{\mu}}^{(b)},\widetilde{\boldsymbol{v}}^{(b)2},\widetilde{\boldsymbol{\beta}}^{(b)})$
with $b=1,2,\ldots,B$ for UQ and
$\widetilde{\boldsymbol{\theta}}^{(b_{O})}=(\widetilde{\boldsymbol{\mu}}^{(b_{O})},\widetilde{\boldsymbol{v}}^{(b_{O})2},\widetilde{\boldsymbol{\beta}}^{(b_{O})})$
with $b_{O}=1,2,\ldots,B_{O}$ for SA;

(2) Call Algorithm 1 to compute
$\mbox{Sh}_{W_{k},X_{m+1}}(\widetilde{\boldsymbol{\theta}}^{(b)})$ and
criticality $p_{W_{k},X_{m+1}}(\widetilde{\boldsymbol{\theta}}^{(b)})$ for
$b=1,2,\ldots,B$;

(3) Calculate the overall estimation uncertainty by using
$\widehat{\mbox{Var}}^{*}[\mbox{Sh}_{W_{k},X_{m+1}}|\mathcal{X}]$ and
$\widehat{\mbox{Var}}^{*}[p_{W_{k},X_{m+1}}|\mathcal{X}]$ in Equations (18)
and (19);

(4) Randomly generate $N_{\pi}$ permutations,
$\pi_{n}\sim\Pi(\mathcal{L}_{k})$ with $n=1,\ldots,N_{\pi}$;

for  _Each $\pi_{n}$_ do

(5) Set $\widehat{c}(P_{\pi_{n}(1)}(\pi_{n}))=0$;

for  _$\ell=1,\ldots,L_{k}$_ do

if  _$\ell <L_{k}$_ then

for  _$b_{O}=1,\ldots,B_{O}$_ do

(7) Set initial value
${\boldsymbol{\theta}}_{\mathcal{J}}^{(b_{O},0)}=\widetilde{\boldsymbol{\theta}}_{\mathcal{J}}^{(b_{O})}$
with ${\mathcal{J}}=P_{\pi_{n}(\ell+1)}(\pi_{n})$;

for  _$t=1,\ldots,T$_ do

(8) For each $\theta_{\mathcal{J}(\ell)}\in\boldsymbol{\theta}_{\mathcal{J}}$,
generate $\theta_{\mathcal{J}(\ell)}^{(b_{O},t)}\sim
p(\theta_{\mathcal{J}(\ell)}|\mathcal{X},\widetilde{\boldsymbol{\theta}}_{\mathcal{L}_{k}-\mathcal{J}}^{(b_{O})},\theta_{\mathcal{J}(1)}^{(b_{O},t)},\ldots,\theta_{\mathcal{J}(\ell-1)}^{(b_{O},t)},\theta_{\mathcal{J}(\ell+1)}^{(b_{O},t-1)},$
$\ldots,\theta_{\mathcal{J}(J)}^{(b_{O},t-1)})$ by applying Equations
(15)/(16)/(17) for the case with complete data or Equations (30)/(31)/(32) for
cases with mixing data (see Appendix D). Obtain the new sample
${\boldsymbol{\theta}}_{\mathcal{J}}^{(b_{O},t)}$;

(9) Set
${\widetilde{\boldsymbol{\theta}}}_{\mathcal{J}}^{(b_{O},b_{I})}={\boldsymbol{\theta}}_{\mathcal{J}}^{(b_{O},(b_{I}-1)h+1)}$
with some constant integer $h$ to reduce the correlation between consecutive
samples;

(10) Compute $\widehat{c}(P_{\pi_{n}(\ell+1)}(\pi_{n}))$ by Equations
(LABEL:eq.costfun2_h) and (LABEL:eq.costfun3_h);

else

(11) Set
$\widehat{c}(P_{\pi_{n}(\ell+1)}(\pi_{n}))=\widehat{\mbox{Var}}^{*}[\mbox{Sh}_{W_{k},X_{m+1}}|\mathcal{X}]$
and $\widehat{\mbox{Var}}^{*}[p_{W_{k},X_{m+1}}|\mathcal{X}]$;

(12) Compute
$\Delta_{\pi_{n}(\ell)}c(\pi_{n})=\widehat{c}(P_{\pi_{n}(\ell+1)}(\pi_{n}))-\widehat{c}(P_{\pi_{n}(\ell)}(\pi_{n}))$;

(13) Estimate
$\widehat{\widehat{\mbox{Sh}}}^{*}_{\theta_{\ell}}\left[\left.\mbox{Sh}_{W_{k},X_{m+1}}\right|\mathcal{X}\right]$
and
$\widehat{\mbox{Sh}}_{\theta_{\ell}}^{*}\left[\left.p_{W_{k},X_{m+1}}\right|\mathcal{X}\right]$
by using Equations (LABEL:eq.sensitivity_h2) and (LABEL:eq.sensitivity_h3).

Algorithm 2 Procedure for the BN-SV-MU Based UQ and SA

Then, we derive SV measuring the estimation uncertainty contribution from each
$\theta_{\ell}$,

$\mbox{Sh}^{*}_{\theta_{\ell}}\left[\left.\mbox{Sh}_{W_{k},X_{m+1}}\right|\mathcal{X}\right]=\sum_{\mathcal{J}\subset\mathcal{L}_{k}/\\{\ell\\}}\dfrac{(L_{k}-|\mathcal{J}|-1)!|\mathcal{J}|!}{L_{k}!}\\\
\left[c(\mathcal{J}\cup\\{\ell\\})-c(\mathcal{J})\right].$

Denote the size of relevant parameters by
$L_{k}=|\boldsymbol{\theta}(W_{k},X_{m+1})|$ and denote the index set by
$\mathcal{L}_{k}$,
$\boldsymbol{\theta}(W_{k},X_{m+1})=\boldsymbol{\theta}_{\mathcal{L}_{k}}$. We
further denote any subset by
$\boldsymbol{\theta}_{\mathcal{J}}\subset\boldsymbol{\theta}(W_{k},X_{m+1})$
with size $J=|\boldsymbol{\theta}_{\mathcal{J}}|$ and the corresponding index
set
$\mathcal{J}=\\{\mathcal{J}(1),\mathcal{J}(2),\ldots,\mathcal{J}(J)\\}\subset\mathcal{L}_{k}$.
For any $\mathcal{J}\subset\mathcal{L}_{k}$, the cost function is given as,

$c(\mathcal{J})=\mbox{E}^{*}_{p(\boldsymbol{\theta}_{\mathcal{L}_{k}-\mathcal{J}}|\mathcal{X})}[\mbox{Var}^{*}_{p(\boldsymbol{\theta}_{\mathcal{J}}|{\boldsymbol{\theta}}_{\mathcal{L}_{k}-\mathcal{J}},\mathcal{X})}[\mbox{Sh}_{W_{k},X_{m+1}}|\widetilde{\boldsymbol{\theta}}_{\mathcal{L}_{k}-\mathcal{J}}]],$
(20)

where
${\boldsymbol{\theta}}_{\mathcal{L}_{k}-\mathcal{J}}=\boldsymbol{\theta}_{\mathcal{L}_{k}\setminus\mathcal{J}}$.
Denote a permutation of $\mathcal{L}_{k}$ as $\pi$ and define the set
$P_{\ell}(\pi)$ as the index set preceding $\ell$ in $\pi$. The SV can be
rewritten as,

$\mbox{Sh}^{*}_{\theta_{\ell}}\left[\left.\mbox{Sh}_{W_{k},X_{m+1}}\right|\mathcal{X}\right]=\sum_{\pi\in\Pi(\mathcal{L}_{k})}\dfrac{1}{L_{k}!}\left[c(P_{\ell}(\pi)\cup\\{\ell\\})-c(P_{\ell}(\pi))\right],$
(21)

where $\Pi(\mathcal{L}_{k})$ denotes the set of all $L_{k}!$ permutations of
$\mathcal{L}_{k}$.

The number of all possible subsets $\mathcal{J}$ could grow exponentially as
$L_{k}$ increase. To address this computational issue, we use the Monte Carlo
sampling approximation, ApproShapley, suggested by [45, 5], which estimates
the Shapley value in (21) by

$\widehat{\mbox{Sh}}^{*}_{\theta_{\ell}}\left[\left.\mbox{Sh}_{W_{k},X_{m+1}}\right|\mathcal{X}\right]=\dfrac{1}{N_{\pi}}\sum_{n=1}^{N_{\pi}}\left[c(P_{\ell}(\pi_{n})\cup\\{\ell\\})-c(P_{\ell}(\pi_{n}))\right]\triangleq\dfrac{1}{N_{\pi}}\sum_{n=1}^{N_{\pi}}\Delta_{\ell}c(\pi_{n}),$
(22)

where $N_{\pi}$ denotes the number of permutations
$\pi_{1},\ldots,\pi_{N_{\pi}}$ randomly generated from $\Pi(\mathcal{L}_{k})$
and
$\Delta_{\ell}c(\pi_{n})=c(P_{\ell}(\pi_{n})\cup\\{\ell\\})-c(P_{\ell}(\pi_{n}))$
is the incremental posterior variance
$\mbox{Var}^{*}[\mbox{Sh}_{W_{k},X_{m+1}}|\mathcal{X}]$ induced by including
the $\ell$-th model parameter in $P_{\ell}(\pi_{n})$.

##  7 Empirical Study

To assess the performance of proposed risk and sensitivity analyses, we first
consider an integrated biomanufacturing process with simulated data in Section
7.1. Then, we study the performance by utilizing the real-world process data
collected from a cell culture process with multiple stages in Section 7.2.
Even though there are many factors impacting on the biomanufacturing outputs,
the amount of real-world bioprocess observations is often very limited.
Therefore, it is important to explore the causal relationships of the
biopharmaceutical production process, which can reduce the model uncertainty,
increase the interpretability for process sensitivity analysis, and guide the
decision making to improve the production stability.

###  7.1 Study the Performance of Proposed Framework with Simulation Data

We revisit the example described at the end of Section 3.2 in Fig. 3. In
total, the process graph model includes 20 nodes, consisting of 10 CPPs
($\mathbf{X}^{p}$) and 8 CQAs ($\mathbf{X}^{a}$) for intermediate product and
2 CQAs ($\mathbf{Y}$) for the final drug substance. The size of coefficients
$\boldsymbol{\theta}$ is 84, including 20 $\mu_{i}$’s, 20 $v_{i}$’s, and 44
$\beta_{ij}$’s coefficients. To study the performance of the proposed
framework, we generate the simulated production process data $\mathcal{X}$,
which mimics the real-world data collection. The BN-based probabilistic
knowledge graph with parameters
$\boldsymbol{\theta}^{c}=(\boldsymbol{\mu}^{c},(\boldsymbol{v}^{2})^{c},\boldsymbol{\beta}^{c})$,
characterizing the underlying bioprocess risk behaviors and CPPs/CQAs
interdependencies, is used for data generation, which is built on the
biomanufacturing domain knowledge; see the detailed setting in Appendix E. To
assess the performance of the proposed framework, we assume that the true
parameter values are unknown. We empirically study the convergence of process
model parameter inference in Appendix F. In Sections 7.1.1 and 7.1.3, we show
the capabilities of the proposed process risk and sensitivity analyses by
studying both process inherent stochastic uncertainty and model uncertainty.

####  7.1.1 Bioprocess Sensitivity Analysis and CPPs/CQAs Criticality
Assessment

We generate the data $\mathcal{X}$ with the number of batch $R=30$ to study
the performance of the proposed risk and sensitivity analyses. For any
intermediate and final product CQA output $X_{i}$ of interest, at each
posterior sample $\widetilde{\boldsymbol{\theta}}$, we follow Algorithm 1 to
assess the criticality of any input factor $W_{k}$ (i.e., CPPs/CQAs, residual
factors). Specifically, in the $h$-th macro-replication of simulation, we
first generate the “real-world" batch data $\mathcal{X}^{(h)}$ with
$h=1,2,\ldots,H$, which is used to mimic the process data collection.
Considering the criticality of input $W_{k}$ to the output variance
$p_{W_{k},X_{i}}(\widetilde{\boldsymbol{\theta}})={\mbox{Sh}_{W_{k},X_{i}}(\widetilde{\boldsymbol{\theta}})}/{\mbox{Var}(X_{i}|\widetilde{\boldsymbol{\theta}})}$,
we estimate the expected value $\mbox{E}[p_{W_{k},X_{i}}]=\iint
p_{W_{k},X_{i}}(\boldsymbol{\theta})dP(\boldsymbol{\theta}|\mathcal{X})dP(\mathcal{X}|\boldsymbol{\theta}^{c})\times
100\%$ by using
$\widehat{\mbox{E}}[p_{W_{k},X_{i}}]=\frac{1}{HB}\sum_{h=1}^{H}\sum_{b=1}^{B}p_{W_{k},X_{i}}(\widetilde{\boldsymbol{\theta}}^{(h,b)})\times
100\%$ with $\widetilde{\boldsymbol{\theta}}^{(h,b)}\sim
p(\boldsymbol{\theta}|\mathcal{X}^{(h)})$ for $h=1,\ldots,H$ and
$b=1,\ldots,B$, with $H=20$ and $B=1000$, and then record the results in terms
of percentage (%) in Tables 1 and 2. Each row records the criticality for each
CPP, CQA, or residual input factor $W_{k}$, and each column corresponds to an
intermediate or final product CQA output $X_{i}$.

Table 1: The estimated criticality level $\widehat{\mbox{E}}[p_{W_{k},X_{i}}]$
and standard deviation $\widehat{\mbox{SD}}[p_{W_{k},X_{i}}]$ (in %) of any
input CPP or other factor $W_{k}$ impacting on the variance of intermediate or
final product CQA $X_{i}$.

$\widehat{\mbox{E}}[p_{W_{k},X_{i}}]$ | $X_{i}=X_{5}$ | $X_{6}$ | $X_{7}$ | $X_{10}$ | $X_{11}$ | $X_{14}$ | $X_{15}$ | $X_{16}$ | $X_{19}$ | $X_{20}$
---|---|---|---|---|---|---|---|---|---|---
$W_{k}=X_{1}$ | 8.91(3.09) | 8.93(3.11) | 9.42(3.59) | 8.51(2.87) | 8.51(2.87) | 5.87(1.88) | 5.87(1.88) | 5.87(1.88) | 5.52(1.74) | 5.52(1.74)
$X_{2}$ | 0.82(0.38) | 0.76(0.35) | 0.96(0.75) | 0.75(0.32) | 0.75(0.32) | 0.52(0.21) | 0.52(0.21) | 0.52(0.2) | 0.49(0.19) | 0.49(0.19)
$X_{3}$ | 4.28(1.6) | 4.33(1.61) | 4.22(1.97) | 4.05(1.46) | 4(1.44) | 2.75(0.93) | 2.75(0.93) | 2.75(0.93) | 2.59(0.86) | 2.59(0.86)
$X_{4}$ | 85.75(4.02) | 85.73(4.03) | 83.29(4.84) | 81.52(4.55) | 81.6(4.53) | 58.2(7.32) | 58.22(7.32) | 58.2(7.32) | 55.09(7.21) | 55.09(7.21)
$e_{5}$ | 0.23(0.1) |  |  | 0.05(0.04) | 0.04(0.04) | 0.03(0.02) | 0.03(0.02) | 0.03(0.02) | 0.03(0.02) | 0.03(0.02)
$e_{6}$ |  | 0.24(0.1) |  | 0.04(0.03) | 0.05(0.04) | 0.03(0.02) | 0.03(0.02) | 0.03(0.02) | 0.03(0.02) | 0.03(0.02)
$e_{7}$ |  |  | 2.11(0.86) | 0.05(0.04) | 0.06(0.04) | 0.03(0.02) | 0.03(0.02) | 0.03(0.02) | 0.03(0.02) | 0.03(0.02)
$X_{8}$ |  |  |  | 1.02(0.41) | 1.01(0.41) | 0.68(0.24) | 0.68(0.25) | 0.69(0.25) | 0.64(0.23) | 0.64(0.23)
$X_{9}$ |  |  |  | 3.91(1.42) | 3.86(1.41) | 2.66(0.9) | 2.68(0.91) | 2.67(0.9) | 2.51(0.84) | 2.51(0.84)
$e_{10}$ |  |  |  | 0.1(0.04) |  | 0.02(0.01) | 0.02(0.01) | 0.02(0.01) | 0.02(0.01) | 0.02(0.01)
$e_{11}$ |  |  |  |  | 0.12(0.05) | 0.02(0.01) | 0.02(0.02) | 0.02(0.01) | 0.02(0.01) | 0.02(0.01)
$X_{12}$ |  |  |  |  |  | 1.86(0.63) | 1.9(0.66) | 1.86(0.63) | 1.76(0.59) | 1.76(0.59)
$X_{13}$ |  |  |  |  |  | 27.31(6.46) | 27.18(6.45) | 27.3(6.46) | 25.72(6.09) | 25.73(6.09)
$e_{14}$ |  |  |  |  |  | 0.02(0.01) |  |  |  $<$0.01($<$0.01) |  $<$0.01($<$0.01)
$e_{15}$ |  |  |  |  |  |  | 0.06(0.02) |  |  $<$0.01($<$0.01) |  $<$0.01($<$0.01)
$e_{16}$ |  |  |  |  |  |  |  | 0.02(0.01) |  $<$0.01($<$0.01) |  $<$0.01($<$0.01)
$X_{17}$ |  |  |  |  |  |  |  |  | 1.27(0.45) | 1.27(0.43)
$X_{18}$ |  |  |  |  |  |  |  |  | 4.23(1.39) | 4.26(1.39)
$e_{19}$ |  |  |  |  |  |  |  |  | 0.04(0.01) |
$e_{20}$ |  |  |  |  |  |  |  |  |  | 0.01($<$0.01)

Table 2: The estimated criticality level $\widehat{\mbox{E}}[p_{W_{k},X_{i}}]$
and standard deviation $\widehat{\mbox{SD}}[p_{W_{k},X_{i}}]$ (in %) of any
input CQA $W_{k}$ on the variance of intermediate or final product CQA
$X_{i}$.

$\widehat{\mbox{E}}[p_{W_{k},X_{i}}]$ | $X_{i}=X_{10}$ | $X_{11}$ | $X_{14}$ | $X_{15}$ | $X_{16}$ | $X_{19}$ | $X_{20}$
---|---|---|---|---|---|---|---
$W_{k}=X_{5}$ | 43.1(11.18) | 38.05(11.44) | 28.68(6.48) | 28.46(6.68) | 28.44(6.48) | 27.05(6.14) | 26.99(6.11)
$X_{6}$ | 37.97(10.79) | 42.44(11.27) | 28.64(6.32) | 28.8(6.5) | 28.88(6.28) | 27.13(5.98) | 27.18(5.95)
$X_{7}$ | 13.91(4.42) | 14.52(4.92) | 10.11(2.62) | 10.2(2.73) | 10.11(2.67) | 9.59(2.49) | 9.6(2.49)
$X_{10}$ |  |  | 37.17(6.55) | 33.62(9.62) | 34.59(6.7) | 33.52(5.69) | 33.01(5.16)
$X_{11}$ |  |  | 33.64(6.42) | 37.24(9.74) | 36.23(6.72) | 33.44(5.69) | 33.95(5.2)
$X_{14}$ |  |  |  |  |  | 32.65(12.69) | 31.91(7.8)
$X_{15}$ |  |  |  |  |  | 21.49(12.22) | 25.1(6.9)
$X_{16}$ |  |  |  |  |  | 40.31(14.51) | 37.45(7.64)

For each CQA output $X_{i}$, we record the criticality with the estimated mean
$\widehat{\mbox{E}}[p_{W_{k},X_{i}}]$ and standard deviation
$\widehat{\mbox{SD}}[p_{W_{k},X_{i}}]$ from any CPP or other factor $W_{k}$ in
Table 1. Under the example setting, we can see that the variations in $X_{4}$
(dissolved oxygen in main fermentation) and $X_{13}$ (temperature in
chromatography) have the dominant impact on both intermediate and final
product CQAs’ variance. Compared with main fermentation and chromatography,
the other two operation units (i.e., centrifuge and filtration) have
relatively small impact on the final product quality variation. Based on the
process risk and sensitivity analyses, we also provide the result
visualization; see for example Fig. 3.

By studying the subplots of the bioprocess probabilistic knowledge graph
illustrated in Fig. 3, we can study the contributions from the dependent CQAs
of intermediate products as inputs to the variance of final drug substance
CQAs outputs, i.e., nodes $\\{X_{19},X_{20}\\}$. We consider the subplots: (1)
starting from the end of main fermentation with $\\{X_{5},X_{6},X_{7}\\}$; (2)
starting from the end of centrifuge with $\\{X_{10},X_{11}\\}$; and (3)
starting from the end of chromatography with $\\{X_{14},X_{15},X_{16}\\}$. The
results of process sensitivity analysis are recorded in Table 2. The CQAs
after main fermentation, i.e., $\\{X_{5},X_{6},X_{7}\\}$, together account for
about 50% variance of final output $X_{19}$ or $X_{20}$; and CQAs after
chromatography, i.e., $\\{X_{14},X_{15},X_{16}\\}$ together account for about
90% of final output variation. Thus, the CQAs of intermediate product close to
the end of production process provides better explanation of the variation of
final drug substance CQAs and we can predict more accurate on its productivity
and quality. This information can be used to guide the production process
quality control and support the real-time release.

####  7.1.2 Criticality Assessment Estimation Performance Comparison

In this section, we use the same example studied in Section 7.1.1 to compare
the performance of criticality assessment obtained by the proposed BN-SV
approach (denoted by $p_{W_{k},X_{20}}^{BN-SV}$) with an existing approach,
which uses multiple linear regression and Morris sensitivity analysis
(represented by ML-M); see [14, 55, 15]. Basically, we first fit the multiple
linear regression to the random inputs (i.e., $W_{k}=X_{i}$ listed in the
first column of Table 1) and output $X_{20}$, and then use Morris sensitivity
analysis to measure the criticality of each input $W_{k}$. Here, we use the
same experiment setting with that used in Section 7.1.1. With the underlying
parameters setting
$\boldsymbol{\theta}^{c}=(\boldsymbol{\mu}^{c},(\boldsymbol{v}^{2})^{c},\boldsymbol{\beta}^{c})$
given in Appendix E, the true criticality of any input factor $W_{k}$ can be
calculated with
$p_{W_{k},X_{20}}^{c}={\mbox{Sh}_{W_{k},X_{20}}(\boldsymbol{\theta}^{c})}/{\mbox{Var}(X_{20}|\boldsymbol{\theta}^{c})}$,
where ${\mbox{Sh}_{W_{k},X_{20}}(\boldsymbol{\theta}^{c})}$ and
${\mbox{Var}(X_{20}|\boldsymbol{\theta}^{c})}$ are obtained by applying
Equations (11) and (12). Then, suppose the underlying process model
coefficients are unknown, and we can compare the criticality assessment
performance of both approaches. In Table 3, we record the mean and SD of
criticality estimates obtained from LM-M and proposed BN-SV approaches with
$H=30$ macro-replications and $R=30$ batches. The mean absolute error (MAE) is
calculated by,

$MAE(p_{W_{k},X_{20}}^{\gamma})=\frac{1}{HB}\sum_{h=1}^{H}\sum_{b=1}^{B}\left|p_{W_{k},X_{20}}^{\gamma}(\widetilde{\boldsymbol{\theta}}^{(h,b)})-p_{W_{k},X_{20}}^{c}\right|\times
100\%$ (23)

where $\gamma$ is ML-M or BN-SV. The results in Table 3 show that the proposed
BN-SV sensitivity analysis provides better criticality assessment of critical
inputs.

Table 3: The CPPs criticality estimation results obtained by BN-SV sensitivity analysis and existing multiple regression based sensitivity analysis.  Criticality (%) |  True Value $p_{W_{k},X_{20}}^{c}$ | $p_{W_{k},X_{20}}^{ML-M}$ | MAE | $p_{W_{k},X_{20}}^{BN-SV}$ | MAE
---|---|---|---|---|---
$W_{k}=X_{4}$ | 59.55 | 56.91 (14.94) | 11.14 | 55.09 (7.21) | 7.86
$X_{13}$ | 24.01 | 26.06 (9.47) | 6.41 | 25.73 (6.09) | 5.19
$X_{1}$ | 4.67 | 5.40 (1.63) | 1.24 | 5.52 (1.74) | 1.41
$X_{18}$ | 3.66 | 4.25 (1.64) | 1.21 | 4.26 (1.39) | 1.13
$X_{3}$ | 2.38 | 2.55 (0.77) | 0.61 | 2.59 (0.86) | 0.47
$X_{9}$ | 2.16 | 2.36 (0.77) | 0.64 | 2.51 (0.84) | 0.68
$X_{12}$ | 1.5 | 1.73 (0.52) | 0.45 | 1.76 (0.59) | 0.42
$X_{17}$ | 1.04 | 1.25 (0.59) | 0.43 | 1.27 (0.43) | 0.35

####  7.1.3 Sensitivity Analysis for Model Uncertainty

Here we consider the product protein content $X_{20}$ in Fig. 3 as the output
to study the performance of sensitivity analysis for model uncertainty. Based
on the results in Table 1, the CPPs $X_{4}$ and $X_{13}$ have the dominant
contributions to the variance of output $X_{20}$, and the estimates of
$p_{W_{k},X_{i}}$ also have the high estimation uncertainty. Thus, we conduct
the BN-SV-MU sensitivity analysis to study how the estimation uncertainty of
each model coefficient impacts on the criticality assessment for
$p_{X_{4},X_{20}}$ and $p_{X_{13},X_{20}}$.

Given the data $\mathcal{X}$, we provide the posterior variance decomposition
studying the criticality estimation uncertainty induced by the MU,
$\mbox{Var}^{*}_{p(\boldsymbol{\theta}|\mathcal{X})}[p_{W_{k},X_{i}}(\widetilde{\boldsymbol{\theta}})|\mathcal{X}]=\sum_{\theta_{\ell}\in\boldsymbol{\theta}(W_{k},X_{i})}\mbox{Sh}_{\theta_{\ell}}^{*}\left[\left.p_{W_{k},X_{i}}(\widetilde{\boldsymbol{\theta}})\right|\mathcal{X}\right].$
Then, we can estimate the expected relative contribution from each model
coefficient $\theta_{\ell}\in\boldsymbol{\theta}(W_{k},X_{i})$ with
$\mbox{EP}_{\theta_{\ell}}(p_{W_{k},X_{i}})\equiv\mbox{E}\left[\frac{{{\mbox{Sh}}}^{*}_{\theta_{\ell}}\left[\left.p_{W_{k},X_{i}}(\widetilde{\boldsymbol{\theta}})\right|\mathcal{X}\right]}{{\mbox{Var}^{*}_{p(\boldsymbol{\theta}|\mathcal{X})}\left[\left.p_{W_{k},X_{i}}(\widetilde{\boldsymbol{\theta}})\right|\mathcal{X}\right]}}\right].$
In the $h$-th macro-replication, given the data $\mathcal{X}^{(h)}$, we can
estimate the contribution from each $\theta_{\ell}$ by using
$\widehat{\widehat{\mbox{Sh}}}^{*}_{\theta_{\ell}}\left[\left.p_{W_{k},X_{i}}(\widetilde{\boldsymbol{\theta}})\right|\mathcal{X}^{(h)}\right]$
and
$\widehat{\mbox{Var}}^{*}_{p(\boldsymbol{\theta}|\mathcal{X})}\left[\left.p_{W_{k},X_{i}}(\widetilde{\boldsymbol{\theta}})\right|\mathcal{X}^{(h)}\right]$
following Equations (LABEL:eq.sensitivity_h2) and (19), which is estimated by
using $N_{\pi}=500$, $B_{O}=5$ and $B_{I}=20$; see [45] for the selection of
sampling parameter setting. Thus, we have the estimation uncertainty
proportion
$\widehat{\mbox{EP}}_{\theta_{\ell}}(p_{W_{k},X_{i}})\equiv\frac{1}{H}\sum_{h=1}^{H}\frac{\widehat{\widehat{\mbox{Sh}}}^{*}_{\theta_{\ell}}\left[\left.p_{W_{k},X_{i}}(\widetilde{\boldsymbol{\theta}})\right|\mathcal{X}^{(h)}\right]}{{\widehat{\mbox{Var}}^{*}_{p(\boldsymbol{\theta}|\mathcal{X})}\left[\left.p_{W_{k},X_{i}}(\widetilde{\boldsymbol{\theta}})\right|\mathcal{X}^{(h)}\right]}}$
with $H=20$.

The coefficients contributing to the estimation of $\mbox{Sh}_{X_{4},X_{20}}$
include $v_{4}^{2}$ and 18 linear coefficients $\boldsymbol{\beta}$ on the
paths from node $X_{4}$ to node $X_{20}$. The coefficients contributing to the
estimation of $\mbox{Sh}_{X_{13},X_{20}}$ include $v_{13}^{2}$ and 6 linear
coefficients $\boldsymbol{\beta}$ located on the paths from $X_{13}$ to
$X_{20}$. Due to the space limit, we only present the top five coefficients
contributing most to the estimation uncertainty of criticality
$p_{X_{4},X_{20}}$ and $p_{X_{13},X_{20}}$, and aggregate the results for
remaining coefficients. The sensitivity analysis results,
$\widehat{\mbox{EP}}_{\theta_{\ell}}(p_{W_{k},X_{i}})\pm\mbox{SE}\Big{[}\widehat{\mbox{EP}}_{\theta_{\ell}}(p_{W_{k},X_{i}})\Big{]}$,
for $p_{X_{4},X_{20}}$ and $p_{X_{13},X_{20}}$ are shown in Table 4, where SE
stands for the standard error (SE). Notice that the coefficients that
contribute the most to the estimation error of the criticality
$p_{X_{4},X_{20}}$ and $p_{X_{13},X_{20}}$ are the variance coefficients of
CPPs ($v_{4}^{2}$ and $v_{13}^{2}$). The estimation uncertainty of linear
coefficients have similar and relatively small contributions. Similar
information can be presented through the sample visualization of the
integrated bioprocess sensitivity analysis in Fig. 3. The darkness of directed
edges and circle boundaries indicates the seriousness of model uncertainty
from corresponding model coefficients $\boldsymbol{\beta}$ and
$\boldsymbol{v}$. Thus, this information can guide the process monitoring and
data collection to efficiently reduce the impact of model uncertainty and
facilitate learning and systematic risk control for integrated
biomanufacturing system.

Table 4: The estimated relative contribution of each BN parameter estimation uncertainty (in terms of %) on criticality assessment $\widehat{\mbox{EP}}_{\theta_{\ell}}(p_{W_{k},X_{i}})\pm\mbox{SE}\Big{[}\widehat{\mbox{EP}}_{\theta_{\ell}}(p_{W_{k},X_{i}})\Big{]}$ for $p_{X_{4},X_{20}}$ and $p_{X_{13},X_{20}}$.  $\theta_{\ell}\in\boldsymbol{\theta}(X_{4},X_{20})$ | $v_{4}^{2}$ | $\beta_{11,15}$ | $\beta_{10,14}$ | $\beta_{15,20}$ | $\beta_{14,20}$ | rest
---|---|---|---|---|---|---
$\widehat{\mbox{EP}}_{\theta_{\ell}}(p_{X_{4},X_{20}})$ |  73.75$\pm$1.96 |  1.59$\pm$0.13 |  1.58$\pm$0.19 |  1.57$\pm$0.21 |  1.56$\pm$0.13 |  19.95$\pm$4.70
$\theta_{\ell}\in\boldsymbol{\theta}(X_{13},X_{20})$ | $v_{13}^{2}$ | $\beta_{13,15}$ | $\beta_{16,20}$ | $\beta_{15,20}$ | $\beta_{13,16}$ | rest
$\widehat{\mbox{EP}}_{\theta_{\ell}}(p_{X_{13},X_{20}})$ |  68.16$\pm$6.40 |  5.57$\pm$0.50 |  5.54$\pm$0.49 |  5.42$\pm$0.44 |  5.23$\pm$0.46 |  10.08$\pm$8.87

###  7.2 Real Case Study for Multiple Phase Cell Culture Process Risk and
Sensitivity Analyses

To study the performance of proposed bioprocess risk and sensitivity analyses,
in this section, we consider the fed-batch fermentation process of Yarrowia
lipolytica yeast for citrate or citric acid (CA) production. This multiple-
phase cell culture process includes seed culture, cell growth and production
processes. In the seed culture, the thawed seed vial solution of Y. lipolytica
strain is transferred to a shake flask containing seed culture medium, and
then grown at $None$C and 280 rpm until cell concentration reaches around 2–5
in OD600 (optical density measured at a wavelength of 600 nm), which usually
takes 18–24h (hours). In the Fed-Batch Fermentation, the seed culture (50 mL)
is first transferred to the bioreactor, which contains the initial
fermentation medium (600 mL) and initial substrate (here we use 35g/L soybean
oil). The feeding starts when the substrate concentration decreases below
20g/L, while the rate is adjusted to maintain the concentration of substrate
about 20 g/L. During the fermentation, the dissolved oxygen level, denoted by
pO2, is set around 30% of air saturation by cascade controls of agitation
speed between 500 and 1,400 rpm, and the aeration rate is fixed at 0.3 L/min.
The pH is controlled at 6.0 during 0–12h, then increased to 7.0 in 6 hours,
and maintained at 7.0 in the remainder of run by feeding KOH (i.e., feed of
base). The temperature is maintained at $None$C for the entire run. At several
middle points of each run, the bioreactor state is estimated by using pH/pO2
probes and off-line sample measurement for residual substrate, which can guide
the adjustment of operation decisions (i.e., feed rate).

In this real case study, we focus on the critical CPPs during the fed-batch
fermentation, including cell concentration after seed culture process, feed
rate, dissolved oxygen (pO2), and residual oil. We consider main CQAs related
to cell (i.e., total cell biomass) and productivity (i.e., total CA
production). Experiments are conducted in Dr. Dongming Xie’s Lab to generate
process data generate the data $\mathcal{X}$ with $R=8$ batches during 140
hours; see the data in Fig. 6. We want to study how the CPPs at different time
contribute to the variation of intermediate and final CQAs outputs, while
evaluating the impact from model uncertainty.

Based on the interactions of CPPs/CQAs, we develop the BN-based bioprocess
probabilistic model with 62 nodes; see the illustration in Fig. 7. We first
estimate the expected criticality $\mbox{E}[p_{W_{k},X_{i}}]$ by using $B$
posterior samples of model coefficients,
$\widehat{\mbox{E}}[p_{W_{k},X_{i}}]=\frac{1}{B}\sum_{b=1}^{B}p_{W_{k},X_{i}}(\widetilde{\boldsymbol{\theta}}^{(b)})\times
100\%$ with $\widetilde{\boldsymbol{\theta}}^{(b)}\sim
p(\boldsymbol{\theta}|\mathcal{X})$ for $b=1,\ldots,B$, with $B=1000$. We
record the results in terms of percentage (%) in Tables 5 and 6 for cell
biomass and CA production respectively with each row and column corresponding
to random input $W_{k}$ and output $X_{i}$. In addition, the overall impact of
model uncertainty on the CPPs/CQAs criticality assessment can be quantified by
the posterior standard deviation (SD), which can be estimated by,
$\widehat{\mbox{SD}}[p_{W_{k},X_{i}}]=\sqrt{\dfrac{1}{(B-1)}\sum_{b=1}^{B}\left[p_{W_{k},X_{i}}\left(\widetilde{\boldsymbol{\theta}}^{(b)}\right)-\bar{p}_{W_{k},X_{i}}\right]^{2}}\times
100\%$, where
$\bar{p}_{W_{k},X_{i}}=\frac{1}{B}\sum_{b=1}^{B}p_{W_{k},X_{i}}(\widetilde{\boldsymbol{\theta}}^{(b)})$.
The results of SD are recorded in the bracket in Tables 5 and 6. Due to the
space limit, we only provide the dominant (high criticality level) part of
time points. We also study the subplots and assess the impact from
intermediate CQAs (biomass and CA amount) as inputs on the following output
variation in Table 7.

Figure 6: Data of citric acid fed-batch fermentation case study. Figure 7: BN
model for citric acid fed-batch fermentation case study. Table 5: The
estimated criticality level $\widehat{\mbox{E}}[p_{W_{k},X_{i}}]$ and standard
error $\widehat{\mbox{SD}}[p_{W_{k},X_{i}}]$ (in %) of any input CPP or other
factor $W_{k}$ impacting on the variance of intermediate or final biomass
$X_{i}$.

$\widehat{\mbox{E}}[p_{W_{k},X_{i}}]$ | $X_{i}=\mbox{BM\\_10}$ | BM_23 | BM_34 | BM_55 | BM_80 | BM_102 | BM_140
---|---|---|---|---|---|---|---
$W_{k}=\mbox{Cell\\_0}$ | 1.88(4.12) | 2.92(5.57) | 1.07(2.76) | 0.76(2.35) | 0.68(2.29) | 0.65(2.22) | 0.56(2.03)
Feed_5 | 15.62(15.79) | 4.29(7.15) | 0.84(1.97) | 0.66(1.81) | 0.57(1.47) | 0.52(1.43) | 0.5(1.51)
pO2_5 | 34.23(22.29) | 14.77(16.21) | 2.83(5.82) | 1.98(4.63) | 1.69(4.52) | 1.5(3.94) | 1.41(3.84)
rOil_5 | 16.07(16.33) | 15.31(15.16) | 5.62(9.25) | 4.02(7.3) | 3.33(6.54) | 3.07(6.41) | 2.81(5.95)
e(BM_5) | 13.6(15.38) | 19.87(18.02) | 7.46(11.46) | 5.35(9.5) | 4.54(8.95) | 4.23(8.65) | 3.85(8.11)
Feed_10 |  | 3.09(5.72) | 2(3.4) | 1.17(2.23) | 1.03(2.27) | 0.98(2.16) | 0.91(2.05)
pO2_10 |  | 8.4(10.58) | 6.34(8.82) | 4.11(6.13) | 3.51(5.65) | 3.29(5.44) | 3.04(5.15)
rOil_10 |  | 9.02(11.96) | 1.15(2.66) | 0.71(1.86) | 0.63(1.78) | 0.58(1.68) | 0.53(1.68)
e(BM_10) | 18.6(13.88) | 5.43(7.15) | 0.74(1.52) | 0.59(1.32) | 0.48(1.15) | 0.43(1.09) | 0.41(1.05)
e(CA_10) |  | 4.62(6.61) | 1.55(3.04) | 0.97(2.03) | 0.82(1.8) | 0.77(1.73) | 0.72(1.64)
Feed_23 |  |  | 21.69(21.63) | 7.23(12.14) | 6.19(11.11) | 5.71(10.56) | 5.17(9.8)
pO2_23 |  |  | 4.41(6.16) | 1.64(2.89) | 1.4(2.63) | 1.27(2.44) | 1.15(2.31)
rOil_23 |  |  | 12.52(15.96) | 5.93(8.72) | 5.22(8.27) | 4.77(8) | 4.31(7.31)
e(BM_23) |  | 12.28(12.71) | 1.47(2.84) | 0.8(1.62) | 0.7(1.54) | 0.63(1.36) | 0.56(1.24)
e(CA_23) |  |  | 0.63(1.39) | 0.4(0.83) | 0.35(1.01) | 0.33(1.03) | 0.31(0.93)
Feed_28 |  |  | 7.06(8.28) | 2.59(4.54) | 2.13(4.16) | 1.97(3.76) | 1.75(3.36)
pO2_28 |  |  | 1.25(2.27) | 1.07(1.91) | 0.85(1.47) | 0.77(1.38) | 0.71(1.35)
rOil_28 |  |  | 13.68(15.45) | 14.07(15.83) | 11.88(14.56) | 10.83(13.98) | 9.84(13.28)
e(BM_28) |  |  | 5.56(5.64) | 1.56(2.44) | 1.34(2.29) | 1.24(2.1) | 1.15(1.99)
e(CA_28) |  |  | 0.04(0.07) | 0.05(0.1) | 0.05(0.1) | 0.04(0.09) | 0.04(0.08)
Feed_34 |  |  |  | 4(6.98) | 3.2(6.14) | 2.67(5.43) | 2.27(5.33)
pO2_34 |  |  |  | 3.02(5.28) | 2.54(4.85) | 2.18(4.27) | 1.96(4.04)
rOil_34 |  |  |  | 3.47(6.19) | 2.76(5.38) | 2.55(5.2) | 2.32(5.04)
e(BM_34) |  |  | 2.09(4.34) | 0.45(1.17) | 0.38(1.08) | 0.32(0.94) | 0.29(0.71)
e(CA_34) |  |  |  | 1.23(2.44) | 1.01(2.3) | 0.88(2.11) | 0.77(1.87)
Feed_55 |  |  |  |  | 2.1(4.45) | 1.68(3.61) | 1.23(2.75)
pO2_55 |  |  |  |  | 3.99(7.64) | 3.07(6.24) | 2.16(4.77)
rOil_55 |  |  |  |  | 7.34(13.28) | 5.86(10.8) | 4.42(8.74)
e(BM_55) |  |  |  | 9.88(10.83) | 5.31(7.39) | 4.4(6.26) | 3.33(5.29)
e(CA_55) |  |  |  |  | 0.05(0.15) | 0.08(0.4) | 0.08(0.54)
Feed_80 |  |  |  |  |  | 0.94(2.56) | 0.67(2.05)
pO2_80 |  |  |  |  |  | 1.12(2.99) | 0.78(2.57)
rOil_80 |  |  |  |  |  | 6.2(11.46) | 4.09(7.93)
e(BM_80) |  |  |  |  | 0.15(0.47) | 0.1(0.27) | 0.07(0.19)
e(CA_80) |  |  |  |  |  | 0.02(0.06) | 0.03(0.08)
Feed_102 |  |  |  |  |  |  | 0.47(1.33)
pO2_102 |  |  |  |  |  |  | 0.4(1.07)
rOil_102 |  |  |  |  |  |  | 7.77(14.9)
e(BM_102) |  |  |  |  |  | 0.62(1.76) | 0.4(1.25)
e(CA_102) |  |  |  |  |  |  | 0(0.01)
Feed_120 |  |  |  |  |  |  | 1.97(3.28)
pO2_120 |  |  |  |  |  |  | 1.19(2.15)
rOil_120 |  |  |  |  |  |  | 4.33(7.71)

Table 6: The estimated criticality level $\widehat{\mbox{E}}[p_{W_{k},X_{i}}]$
and standard error $\widehat{\mbox{SD}}[p_{W_{k},X_{i}}]$ (in %) of any input
CPP or other factor $W_{k}$ impacting on the variance of intermediate or final
CA amount $X_{i}$.

$\widehat{\mbox{E}}[p_{W_{k},X_{i}}]$ | $X_{i}=\mbox{CA\\_10}$ | CA_23 | CA_34 | CA_55 | CA_80 | CA_102 | CA_140
---|---|---|---|---|---|---|---
$W_{k}=\mbox{Cell\\_0}$ | 5.22(7.79) | 2.41(4.68) | 1.22(3.11) | 0.96(2.83) | 0.89(2.69) | 0.82(2.39) | 0.79(2.46)
Feed_5 | 7.88(10.99) | 3.27(5.04) | 1.25(2.2) | 0.82(1.81) | 0.75(1.72) | 0.68(1.62) | 0.65(1.59)
pO2_5 | 5.5(8.81) | 4.98(9.24) | 3.05(6.24) | 2.41(5.46) | 2.21(5.09) | 2.09(4.95) | 2.01(4.74)
rOil_5 | 30.63(18.78) | 13.61(13.6) | 6.37(9.37) | 4.9(8.2) | 4.63(8.18) | 4.36(7.89) | 4.24(7.87)
e(BM_5) | 38.42(18.88) | 16.6(15.47) | 8.09(11.72) | 6.33(10.22) | 6(10.14) | 5.59(9.63) | 5.39(9.51)
Feed_10 |  | 12.64(12.02) | 2.49(3.9) | 1.61(2.76) | 1.45(2.53) | 1.35(2.45) | 1.3(2.46)
pO2_10 |  | 34.18(21.2) | 7.97(9.38) | 5.23(7.39) | 4.79(6.94) | 4.42(6.62) | 4.2(6.39)
rOil_10 |  | 1.52(3.62) | 1.2(2.53) | 0.97(2.58) | 0.9(2.5) | 0.84(2.34) | 0.77(2.16)
e(BM_10) |  | 1.6(3.54) | 1.06(2.14) | 0.71(1.46) | 0.66(1.45) | 0.62(1.41) | 0.6(1.38)
e(CA_10) | 12.36(13.95) | 3.69(5.19) | 1.52(2.98) | 1.18(2.24) | 1.12(2.2) | 1.03(2.08) | 0.99(2.07)
Feed_23 |  |  | 8.07(10.91) | 9.44(13.3) | 8.23(12.12) | 7.7(11.7) | 7.26(11.28)
pO2_23 |  |  | 2.71(3.52) | 2.17(3.5) | 1.91(3.08) | 1.82(3.03) | 1.74(2.96)
rOil_23 |  |  | 11.22(11.85) | 7.8(10.56) | 6.99(9.71) | 6.63(9.55) | 6.32(9.31)
e(BM_23) |  |  | 1.15(1.87) | 1.04(1.95) | 0.92(1.63) | 0.84(1.52) | 0.81(1.5)
e(CA_23) |  | 5.51(7.08) | 0.91(1.58) | 0.56(1.31) | 0.51(1.24) | 0.47(1.18) | 0.45(1.12)
Feed_28 |  |  | 3.5(7.99) | 3.91(6.66) | 3.6(6.27) | 3.38(5.89) | 3.19(5.48)
pO2_28 |  |  | 2.43(3.6) | 1.35(2.29) | 1.25(2.19) | 1.17(2.09) | 1.12(2.09)
rOil_28 |  |  | 29.32(20.85) | 17.37(17.85) | 16.13(17.22) | 14.9(16.39) | 14.35(16.17)
e(BM_28) |  |  | 1.17(2.25) | 2.28(3.33) | 2(2.97) | 1.88(2.83) | 1.77(2.73)
e(CA_28) |  |  | 0.15(0.24) | 0.07(0.12) | 0.06(0.11) | 0.06(0.11) | 0.05(0.1)
Feed_34 |  |  |  | 6.47(11.6) | 4.79(8.82) | 4.22(7.96) | 3.92(7.42)
pO2_34 |  |  |  | 6.22(9.95) | 4.72(8.05) | 4.18(7.17) | 3.84(6.77)
rOil_34 |  |  |  | 4.85(7.94) | 3.97(7.04) | 3.63(6.56) | 3.36(6.3)
e(BM_34) |  |  |  | 0.71(1.79) | 0.63(1.58) | 0.56(1.41) | 0.52(1.31)
e(CA_34) |  |  | 5.16(8.58) | 1.74(3.99) | 1.52(3.31) | 1.42(3.25) | 1.32(3.08)
Feed_55 |  |  |  |  | 0.5(1.55) | 0.54(1.59) | 0.57(1.57)
pO2_55 |  |  |  |  | 0.62(1.86) | 0.7(2.02) | 0.7(1.99)
rOil_55 |  |  |  |  | 2.97(5.6) | 2.76(5.43) | 2.82(6.04)
e(BM_55) |  |  |  |  | 1.58(3.7) | 1.72(3.62) | 1.78(3.58)
e(CA_55) |  |  |  | 0.5(1.34) | 0.46(1.68) | 0.38(1.23) | 0.35(1.21)
Feed_80 |  |  |  |  |  | 0.72(3.18) | 0.66(2.92)
pO2_80 |  |  |  |  |  | 0.7(2.98) | 0.67(2.9)
rOil_80 |  |  |  |  |  | 2.75(6.66) | 2.58(6.22)
e(BM_80) |  |  |  |  |  | 0(0.02) | 0.01(0.03)
e(CA_80) |  |  |  |  | 0.34(1.12) | 0.29(0.96) | 0.26(0.84)
Feed_102 |  |  |  |  |  |  | 0.16(0.49)
pO2_102 |  |  |  |  |  |  | 0.04(0.18)
rOil_102 |  |  |  |  |  |  | 1.77(4.04)
e(BM_102) |  |  |  |  |  |  | 0.02(0.16)
e(CA_102) |  |  |  |  |  | 0.06(0.22) | 0.05(0.21)
Feed_120 |  |  |  |  |  |  | 0.42(1.97)
pO2_120 |  |  |  |  |  |  | 0.32(1.48)
rOil_120 |  |  |  |  |  |  | 1.03(3.77)

Table 7: The estimated criticality level $\widehat{\mbox{E}}[p_{W_{k},X_{i}}]$
and standard error $\widehat{\mbox{SD}}[p_{W_{k},X_{i}}]$ (in %) of any input
CQA $W_{k}$ impacting on the variance of intermediate or final biomass or CA
amount $X_{i}$.

$\widehat{\mbox{E}}[p_{W_{k},X_{i}}]$ | $X_{i}=\mbox{BM\\_55}$ | CA_55 | BM_80 | CA_80 | BM_102 | CA_102 | BM_140 | CA_140
---|---|---|---|---|---|---|---|---
$W_{k}=\mbox{BM\\_10}$ | 3.59(7.62) | 4.29(8.09) | 2.99(6.7) | 3.95(7.58) | 2.67(6.23) | 3.74(7.34) | 2.45(5.83) | 3.61(7.15)
CA_10 | 11.09(16.78) | 13.45(18.7) | 9.41(15.62) | 12.72(18.59) | 8.74(14.99) | 11.82(17.68) | 8.04(14.08) | 11.44(17.56)
BM_23 | 11.29(18.41) | 13.97(20.33) | 9.59(16.82) | 12.93(19.42) | 8.71(15.77) | 12.01(18.61) | 7.99(14.9) | 11.53(18.18)
CA_23 | 11.02(15.02) | 13.72(17.42) | 9.43(14.17) | 12.81(16.92) | 8.88(13.8) | 11.94(16.22) | 8.16(13.01) | 11.46(16.07)
BM_28 | 20.66(25.26) | 27.56(27.27) | 17.94(23.75) | 24.68(26.25) | 16.92(22.87) | 23.42(25.67) | 15.56(21.97) | 22.26(25.2)
CA_28 | 17.94(19.6) | 21.7(20.84) | 15.18(18.26) | 20.08(20.67) | 13.63(17.48) | 18.47(19.75) | 12.34(16.23) | 17.74(19.64)
BM_34 | 25.11(29.47) | 35.52(33.18) | 21.67(26.95) | 31.32(31.48) | 20.52(26.31) | 29.67(30.59) | 18.67(25.32) | 27.91(29.73)
CA_34 | 35.03(33.43) | 40.71(35.87) | 29.32(31.31) | 38.45(35.08) | 26.44(29.72) | 35.59(33.99) | 24.04(28.29) | 34.43(33.78)
BM_47.5 | 37.75(35.28) | 5.79(13.43) | 25.9(30.21) | 12.71(20.88) | 22.94(28.07) | 12.96(20.94) | 19.66(26.42) | 13.03(20.92)
CA_47.5 | 36.49(34.21) | 93.23(14.45) | 36.68(32.59) | 75.43(28.48) | 34.07(31.48) | 68.86(30.46) | 32.04(31.25) | 64.62(31.61)
BM_55 |  |  | 58.85(34.64) | 12.7(22.51) | 49.22(35.71) | 14.46(24.33) | 39.99(34.47) | 15.33(25.24)
CA_55 |  |  | 19.49(25.99) | 82.28(25.87) | 22.22(28.09) | 75.43(29.57) | 23.04(28.06) | 70.83(31.2)
BM_72 |  |  | 95.2(9.58) | 3.25(9.46) | 73.79(29.1) | 7.78(15.83) | 57.99(33.23) | 9.38(18.1)
CA_72 |  |  | 3.03(8.26) | 95.77(10.77) | 13.71(22.8) | 86.34(20.9) | 16.95(24.58) | 81.4(24.37)
BM_80 |  |  |  |  | 76.49(27.65) | 6.85(14.48) | 59.42(32.6) | 9.15(17.69)
CA_80 |  |  |  |  | 12.56(21.31) | 88.21(19.51) | 16.54(23.74) | 82.39(23.6)
BM_95 |  |  |  |  | 94.44(10.84) | 1.48(3.17) | 70.6(28.14) | 5.15(11.38)
CA_95 |  |  |  |  | 4.39(9.63) | 98.43(3.24) | 12.25(18.9) | 91.22(15.09)
BM_102 |  |  |  |  |  |  | 72.79(27.28) | 5.28(11.54)
CA_102 |  |  |  |  |  |  | 10.93(17.17) | 91.31(15.3)
BM_120 |  |  |  |  |  |  | 84.91(18.97) | 3.16(8.28)
CA_120 |  |  |  |  |  |  | 7.81(14.19) | 95.19(11.36)

Differing with the simulation study in Section 7.1, there is no macro-
replication in the real case study. Notice that the posterior standard
deviation (SD) can measure the overall model uncertainty, i.e., the variation
of criticality estimates cross different posterior samples characterizing the
model coefficient estimation uncertainty. Based on the sample average of $B$
posterior samples
$\bar{p}_{W_{k},X_{i}}=\frac{1}{B}\sum_{b=1}^{B}p_{W_{k},X_{i}}(\widetilde{\boldsymbol{\theta}}^{(b)})\times
100\%$, the estimation accuracy of criticality $\bar{p}_{W_{k},X_{i}}$ is
measured by the standard error (SE) with
$\mbox{SE}(\bar{p}_{W_{k},X_{i}})=SD(\bar{p}_{W_{k},X_{i}})/\sqrt{B}$.

The results in Tables 5 and 6 show that the variations of residual oil and
feed rate in the cell growth phase (about from time $t=23h$ to $28h$) have
dominant impact on both intermediate and final cell biomass and CA
productivity. As fermentation time further increases, the criticality level of
input factors $W_{k}$ on the output $X_{i}$, CA production, tends to decrease.
It matches well with the data in Fig. 6, the cell growth and production both
become slower and more stable. This observation suggests that controlling the
CPPs (i.e., feed rate and residual oil) to ensure the good cell growth stage
is more important in order to improve the process stability. For the cell
total biomass output in Table 5, since the residual oil generates the
scattering particles impacting on OD600 and cell biomass measurement accuracy,
this effect becomes larger as the residual oil increases, which explains the
high contribution of residual oil at the end of process (i.e., rOil_102) to
the final cell biomass measurement variation.

By studying the subplots, we study the impact of middle step CQAs (i.e., cell
biomass and CA amount in the cell growth and production phases) on the final
output variation. We record the results in Table 7. As the fermentation time
$t$ increases, the explained variations of final biomass and CA by current
values increase. They reach to around 70% for biomass and 90% for CA at time
$t=95h$. This observation is consistent with the data in Fig. 6 and the growth
of biomass/CA is relative slow in the periods after it. However, compared with
CA, biomass has relatively larger prediction variation even in the later
stages of production, which can be explained by the measurement errors of Cell
OD induced by large amount of residual oil. In terms of biomass impacting on
final CA (or CA impacting on biomass), the most critical part is biomass at
time $t=34h$ (or CA amount at $t=47.5h$). Since cell growth needs nitrogen,
the production phase usually starts when nitrogen concentration becomes small.
During those periods (around 30h to 50h), nitrogen from the initial medium is
consumed and both intracellular lipid (which becomes part of biomass)
accumulation and extracellular CA production are induced by nitrogen
limitation. It can be also observed in Fig. 6, where the sudden increase of CA
total slopes happens around 30h to 50h, whose variations have critical
contribution to final CA output uncertainty.

We also conduct the sensitivity analysis studying the impact of model
uncertainty on the CPPs/CQAs criticality assessment. Here we focus on
criticality assessment estimation of $p_{rOil\\_{28},CA\\_{140}}$ and
$p_{rOil\\_{102},BM\\_{140}}$, which have high criticality and overall model
uncertainty; see Tables 5 and 6. The model coefficients contributing to the
estimation of $p_{rOil\\_{28},CA\\_{140}}$ include $v_{rOil\\_{28}}^{2}$ and
32 linear coefficients $\boldsymbol{\beta}$ on the paths from node
$rOil\\_{28}$ to node $CA\\_{140}$, whereas coefficients contributing to the
estimation of $p_{Feed\\_{23},CA\\_{140}}$ include $v_{Feed\\_{23}}^{2}$ and
36 linear coefficients $\boldsymbol{\beta}$. We present the top five
coefficients contributors to the estimation uncertainty of criticality
$p_{rOil\\_{28},CA\\_{140}}$ and $p_{Feed\\_{23},CA\\_{140}}$, and aggregate
the results for remaining coefficients in Table 8. From the results, the
estimation uncertainty of variance coefficients of CPPs ($v^{2}_{rOil\\_28}$
and $v^{2}_{Feed\\_23}$) have the largest contribution to the estimation
uncertainty of the criticality $p_{rOil\\_{28},CA\\_{140}}$ and
$p_{Feed\\_{23},CA\\_{140}}$. The estimation uncertainty of coefficients
$\boldsymbol{\beta}$ in both sets
$\boldsymbol{\theta}(rOil\\_{28},CA\\_{140})$ and
$\boldsymbol{\theta}(Feed\\_{23},CA\\_{140})$ have similar and relative lower
contributions.

Table 8: The estimated relative contribution of each BN parameter estimation uncertainty (in terms of %) on criticality assessment $\widehat{\mbox{EP}}_{\theta_{\ell}}(p_{W_{k},X_{i}})$ for $p_{Feed\\_{20},CA\\_{140}}$.  $\theta_{\ell}\in\boldsymbol{\theta}(rOil\\_{28},CA\\_{140})$ | $v^{2}_{rOil\\_28}$ | $\beta_{CA\\_72,BM\\_80}$ | $\beta_{CA\\_55,CA\\_72}$ | $\beta_{rOil\\_28,CA\\_34}$ | $\beta_{CA\\_120,CA\\_140}$ | rest
---|---|---|---|---|---|---
$\widehat{\mbox{EP}}_{\theta_{\ell}}(p_{rOil\\_{28},CA\\_{140}})$ | 23.38 | 4.00 | 4.00 | 3.86 | 3.55 | 61.21
$\theta_{\ell}\in\boldsymbol{\theta}(Feed\\_{23},CA\\_{140})$ | $v^{2}_{Feed\\_23}$ | $\beta_{BM\\_72,BM\\_80}$ | $\beta_{CA\\_72,BM\\_80}$ | $\beta_{BM\\_95,CA\\_102}$ | $\beta_{BM\\_72,CA\\_80}$ | rest
$\widehat{\mbox{EP}}_{\theta_{\ell}}(p_{Feed\\_{23},CA\\_{140}})$ | 19.52 | 5.10 | 4.10 | 3.97 | 3.95 | 63.36

##  8 Conclusions

Driven by the critical challenges in biomanufacturing, we create an integrated
bioprocess knowledge graph and propose interpretable risk and sensitivity
analyses, which can provide the production process risk- and science-based
understanding, guide the CPPs/CQAs specifications and production stability
control, and facilitate the process development. Since hundreds of factors can
impact on the product quality and productivity, and also the amount of process
observations is often very limited, we explore the process interactions and
causal relationships, and then develop a Bayesian network (BN) based
probabilistic knowledge graph characterizing the causal interdependencies of
production process CPPs/CQAs. Building on the knowledge graph, we propose the
BN-SV based sensitivity analysis to assess the criticality of each random
input factor on the variance of intermediate/final product quality attributes
by using the Shapley value (SV), which can correctly account for input
interdependencies and process structural interactions. We further introduce
the BN-SV-MU sensitivity analysis, which can provide the comprehensive
understanding on how the estimation uncertainty of each part of process model
coefficients impacts on the production risk analysis and CPPs/CQAs criticality
assessment. It can guide bioprocess sensor monitoring and “most informative"
data collection to facilitate bioprocess learning and model uncertainty
reduction. Both simulation and real case studies are used to demonstrate the
promising performance of proposed bioprocess risk and sensitivity analyses.

## ACKNOWLEDGEMENTS

The authors are grateful for constructive comments from Dr. Barry Nelson
(Northeastern University), Peter Baker (Green Mountain Quality Assurance,
LLC), help from Hua Zheng (NEU) on the development of bioprocess knowledge
graph visualization, and help from Na Liu (UMass Lowell) conducting lab
experiments.

## References

  * Ayech et al., [2012] Ayech, N., Chakour, C., and HARKAT, M.-F. (2012).  New adaptive moving window pca for process monitoring.  IFAC Proceedings Volumes, 45(20):606–611.
  * Borchert et al., [2019] Borchert, D., Suarez-Zuluaga, D. A., Sagmeister, P., Thomassen, Y. E., and Herwig, C. (2019).  Comparison of data science workflows for root cause analysis of bioprocesses.  Bioprocess and biosystems engineering, 42(2):245–256.
  * Borgonovo and Plischke, [2016] Borgonovo, E. and Plischke, E. (2016).  Sensitivity analysis: a review of recent advances.  European Journal of Operational Research, 248(3):869–887.
  * Castillo et al., [1997] Castillo, E., Gutiérrez, J. M., and Hadi, A. S. (1997).  Sensitivity analysis in discrete bayesian networks.  IEEE Transactions on Systems, Man, and Cybernetics-Part A: Systems and Humans, 27(4):412–423.
  * Castro et al., [2009] Castro, J., Gómez, D., and Tejada, J. (2009).  Polynomial calculation of the shapley value based on sampling.  Computers & Operations Research, 36(5):1726–1730.
  * Coleman and Block, [2006] Coleman, M. C. and Block, D. E. (2006).  Retrospective optimization of time-dependent fermentation control strategies using time-independent historical data.  Biotechnology and bioengineering, 95(3):412–423.
  * De Lira et al., [2010] De Lira, L. d. F. B., De Vasconcelos, F. V. C., Pereira, C. F., Paim, A. P. S., Stragevitch, L., and Pimentel, M. F. (2010).  Prediction of properties of diesel/biodiesel blends by infrared spectroscopy and multivariate calibration.  Fuel, 89(2):405–409.
  * Doran, [1995] Doran, P. M. (1995).  Bioprocess engineering principles.  Elsevier.
  * Doran, [2013] Doran, P. M. (2013).  Bioprocess Engineering Principles.  Elsevier.
  * Feng et al., [2014] Feng, N., Wang, H. J., and Li, M. (2014).  A security risk analysis model for information systems: Causal relationships of risk factors and vulnerability propagation analysis.  Information sciences, 256:57–73.
  * Fleischhacker and Zhao, [2011] Fleischhacker, A. J. and Zhao, Y. (2011).  Planning for demand failure: A dynamic lot size model for clinical trial supply chains.  European Journal of Operational Research, 211(3):496–506.
  * Gelman et al., [2004] Gelman, A., Carlin, J. B., Stern, H. S., and Rubin, D. B. (2004).  Bayesian Data Analysis.  Taylor and Francis Group, LLC, New York, 2nd edition.
  * Guideline et al., [2009] Guideline, I. H. T. et al. (2009).  Pharmaceutical development.  Q8 (R2) Current Step, 4.
  * Hassan et al., [2013] Hassan, S. S., Farhan, M., Mangayil, R., Huttunen, H., and Aho, T. (2013).  Bioprocess data mining using regularized regression and random forests.  BMC Systems Biology.
  * Helton, [1993] Helton, J. C. (1993).  Uncertainty and sensitivity analysis techniques for use in performance assessment for radioactive waste disposal.  Reliability Engineering & System Safety, 42(2-3):327–367.
  * Kabir et al., [2015] Kabir, G., Tesfamariam, S., Francisque, A., and Sadiq, R. (2015).  Evaluating risk of water mains failure using a bayesian belief network model.  European Journal of Operational Research, 240(1):220–234.
  * Kaminsky and Wang, [2015] Kaminsky, P. and Wang, Y. (2015).  Analytical models for biopharmaceutical operations and supply chain management: A survery of research literature, pharmaceutical bioprocess.  Pharmaceutical Bioprocess, 2:61–73.
  * Koller and Friedman, [2009] Koller, D. and Friedman, N. (2009).  Probabilistic graphical models: principles and techniques.  MIT press.
  * Kulkarni, [2015] Kulkarni, N. S. (2015).  A modular approach for modeling active pharmaceutical ingredient manufacturing plant: a case study.  In Proceedings of the 2015 Winter Simulation Conference, pages 2260–2271. IEEE Press.
  * Kyriakopoulos et al., [2018] Kyriakopoulos, S., Ang, K. S., Lakshmanan, M., Huang, Z., Yoon, S., Gunawan, R., and Lee, D.-Y. (2018).  Kinetic modeling of mammalian cell culture bioprocessing: The quest to advance biomanufacturing.  Biotechnology Journal, 13(3):1700229.
  * Lakhdar et al., [2007] Lakhdar, K., Savery, J., Papageorgiou, L., and Farid, S. (2007).  Multiobjective long-term planning of biopharmaceutical manufacturing facilities.  Biotechnology progress, 23(6):1383–1393.
  * Leachman et al., [2014] Leachman, R. C., Johnston, L., Li, S., and Shen, Z.-J. (2014).  An automated planning engine for biopharmaceutical production.  European Journal of Operational Research, 238(1):327–338.
  * Li and Yuan, [2006] Li, Y. and Yuan, J. (2006).  Prediction of key state variables using support vector machines in bioprocesses.  Chemical Engineering & Technology: Industrial Chemistry-Plant Equipment-Process Engineering-Biotechnology, 29(3):313–319.
  * Li and Venkatasubramanian, [2018] Li, Y. F. and Venkatasubramanian, V. (2018).  Neural network to understand process capability and process intermediates acceptance criteria in monoclonal antibody production process.  Journal of Pharmaceutical Innovation, 13(1):36–47.
  * Lim et al., [2004] Lim, A. C., Zhou, Y., Washbrook, J., Titchener-Hooker, N. J., and Farid, S. (2004).  A decisional-support tool to model the impact of regulatory compliance activities in the biomanufacturing industry.  Computers & chemical engineering, 28(5):727–735.
  * Lu et al., [2018] Lu, Q., Jiang, B., Gopaluni, R. B., Loewen, P. D., and Braatz, R. D. (2018).  Sparse canonical variate analysis approach for process monitoring.  Journal of Process Control, 71:90–102.
  * Lundberg and Lee, [2017] Lundberg, S. M. and Lee, S.-I. (2017).  A unified approach to interpreting model predictions.  In Advances in neural information processing systems, pages 4765–4774.
  * Martagan et al., [2018] Martagan, T., Krishnamurthy, A., and Leland, P. (2018).  Managing trade-offs in protein manufacturing: how much to waste?  Manufacturing & Service Operations Management.
  * Martagan et al., [2017] Martagan, T., Krishnamurthy, A., Leland, P. A., and Maravelias, C. T. (2017).  Performance guarantees and optimal purification decisions for engineered proteins.  Operations Research, 66(1):18–41.
  * Martagan et al., [2016] Martagan, T., Krishnamurthy, A., and Maravelias, C. T. (2016).  Optimal condition-based harvesting policies for biomanufacturing operations with failure risks.  IIE Transactions, 48(5):440–461.
  * Mitchell, [2013] Mitchell, M. (2013).  Determining criticality-process parameters and quality attributes part i: criticality as a continuum.  BioPharm International, 26(12).
  * Moullec et al., [2013] Moullec, M.-L., Bouissou, M., Jankovic, M., Bocquet, J.-C., Réquillard, F., Maas, O., and Forgeot, O. (2013).  Toward system architecture generation and performances assessment under uncertainty using bayesian networks.  Journal of Mechanical Design, 135(4):041002.
  * Ojha et al., [2018] Ojha, R., Ghadge, A., Tiwari, M. K., and Bititci, U. S. (2018).  Bayesian network modelling for supply chain risk propagation.  International Journal of Production Research, 56(17):5795–5819.
  * Otto et al., [2014] Otto, R., Santagostino, A., and Schrader, U. (2014).  From science to operations: Questions, choices, and strategies for success in biopharma.  McKinsey & Company.
  * Owen, [2014] Owen, A. B. (2014).  Sobol’indices and shapley value.  SIAM/ASA Journal on Uncertainty Quantification, 2(1):245–251.
  * Pharmaceutical Current Good Manufacturing Practices , 2004 [CGMPs] Pharmaceutical Current Good Manufacturing Practices (CGMPs) (2004).  Guidance for industry: Pat—a framework for innovative pharmaceutical development, manufacturing, and quality assurance.  Technical report.
  * Prinsloo et al., [2008] Prinsloo, N. M., Engelbrecht, J. P., Mashapa, T. N., and Strauss, M. J. (2008).  Acetone to mibk process optimization through multidisciplinary chemometrics and in-line nir spectroscopy.  Applied Catalysis A: General, 344(1-2):20–29.
  * Rabitz and Aliş, [1999] Rabitz, H. and Aliş, Ö. F. (1999).  General foundations of high-dimensional model representations.  Journal of Mathematical Chemistry, 25(2-3):197–233.
  * Rader and Langer, [2019] Rader, R. A. and Langer, E. S. (2019).  Single-use technologies in biopharmaceutical manufacturing: A 10-year review of trends and the future.  Single-Use Technology in Biopharmaceutical Manufacture, pages 193–200.
  * Rathore et al., [2010] Rathore, A., Bhambure, R., and Ghare, V. (2010).  Process analytical technology (pat) for biopharmaceutical products.  Analytical and bioanalytical chemistry, 398(1):137–154.
  * Severson et al., [2015] Severson, K., VanAntwerp, J. G., Natarajan, V., Antoniou, C., Thömmes, J., and Braatz, R. D. (2015).  Elastic net with monte carlo sampling for data-based modeling in biopharmaceutical manufacturing facilities.  Computers & Chemical Engineering, 80:30–36.
  * Shapley, [1953] Shapley, L. S. (1953).  A value for n-person games.  Contributions to the Theory of Games, 2(28):307–317.
  * Sobol, [1993] Sobol, I. M. (1993).  Sensitivity estimates for nonlinear mathematical models.  Mathematical modelling and computational experiments, 1(4):407–414.
  * Sokolov et al., [2018] Sokolov, M., Morbidelli, M., Butté, A., Souquet, J., and Broly, H. (2018).  Sequential multivariate cell culture modeling at multiple scales supports systematic shaping of a monoclonal antibody toward a quality target.  Biotechnology Journal, 13(4):1700461.
  * Song et al., [2016] Song, E., Nelson, B. L., and Staum, J. (2016).  Shapley effects for global sensitivity analysis: Theory and computation.  SIAM/ASA Journal on Uncertainty Quantification, 4(1):1060–1083.
  * Steinwandter et al., [2019] Steinwandter, V., Borchert, D., and Herwig, C. (2019).  Data science tools and applications on the way to pharma 4.0.  Drug discovery today, 24(9):1795–1805.
  * Telenko and Seepersad, [2014] Telenko, C. and Seepersad, C. C. (2014).  Probabilistic graphical modeling of use stage energy consumption: a lightweight vehicle example.  Journal of Mechanical Design, 136(10):101403.
  * Troyanskaya et al., [2003] Troyanskaya, O. G., Dolinski, K., Owen, A. B., Altman, R. B., and Botstein, D. (2003).  A bayesian framework for combining heterogeneous data sources for gene function prediction (in saccharomyces cerevisiae).  Proceedings of the National Academy of Sciences, 100(14):8348–8353.
  * Van der Gaag et al., [2007] Van der Gaag, L. C., Renooij, S., and Coupé, V. M. (2007).  Sensitivity analysis of probabilistic networks.  In Advances in probabilistic graphical models, pages 103–124. Springer.
  * Wagner, [1995] Wagner, H. M. (1995).  Global sensitivity analysis.  Operations Research, 43(6):948–969.
  * Walsh, [2013] Walsh, G. (2013).  Pharmaceutical biotechnology: concepts and applications.  John Wiley & Sons.
  * Wang et al., [2018] Wang, Y., Blache, R., Zheng, P., and Xu, X. (2018).  A knowledge management system to support design for additive manufacturing using bayesian networks.  Journal of Mechanical Design, 140(5):051701.
  * Wechselberger et al., [2012] Wechselberger, P., Sagmeister, P., Engelking, H., Schmidt, T., Wenger, J., and Herwig, C. (2012).  Efficient feeding profile optimization for recombinant protein production using physiological information.  Bioprocess and biosystems engineering, 35(9):1637–1649.
  * Zhai et al., [2014] Zhai, Q., Yang, J., Xie, M., and Zhao, Y. (2014).  Generalized moment-independent importance measures based on minkowski distance.  European Journal of Operational Research, 239(2):449–455.
  * Zi, [2011] Zi, Z. (2011).  Sensitivity analysis approaches applied to systems biology models.  The Institution of Engineering and Technology, 5:336–346.

##  Appendix A Ontology based Data and Process Integration

By exploring the causal relationships and interactions in the production
processes, we introduce bioprocess ontology-based data integration, which can
connect all distributed and heterogeneous data collected from bioprocess. This
relational graph can enable the connectivity of end-to-end process from drug
development to patient response; see Fig. 8 for a simplified illustration of
integrated biopharmaceutical manufacturing supply chain. Nodes represent
factors (i.e., CPPs/CQAs, media feed, bioreactor operating conditions, other
uncontrolled factors) impacting the process outputs, and the directed edges
model the causal relationships. Each dashed block could represent a module,
which can be each phase or each unit operation. In this relational graph, the
shaded nodes represent the variables with real-world observations, including
the testing and sensor monitoring data of CPPs/CQAs for raw materials,
operation conditions, and intermediate/final drug products. The unshaded and
dashed nodes represent variables without observations and residuals, including
the complete quality status of intermediate and final drug products, and other
uncontrollable factors (e.g., contamination) introduced during the process
unit operations.

Figure 8: Biopharmaceutical production process ontology based causal
relationships.

##  Appendix B Detailed Derivation of Equation (8)

In order to show Equation (8), we consider more general results as following,

$X_{n}=\mu_{n}+\sum_{k=1}^{m^{p}}\gamma_{k,n}(X_{k}-\mu_{k})+\sum_{k=m^{p}+1}^{n}\gamma_{k,n}e_{k},$
(24)

for $n=m^{p}+1,\ldots,m+1$, where $\gamma_{k,n}$ is given as Equations (9) and
(10). Notice according to linear Gaussian model (6), we can write
$X_{m^{p}+1}=\mu_{m^{p}+1}+\sum_{k=1}^{m^{p}}\beta_{k,m^{p}+1}(X_{k}-\mu_{k})+e_{m^{p}+1}$,
where $\beta_{k,m^{p}+1}=0$ for $k\notin Pa(X_{m^{p}+1})$. Suppose Equation
(24) holds for all $n=m^{p}+1,\ldots,n_{0}$. For $n=n_{0}+1$, by applying
linear Gaussian model, we have

$\displaystyle
X_{n_{0}+1}=\mu_{n_{0}+1}+\sum_{k=1}^{n_{0}}\beta_{k,n_{0}+1}(X_{k}-\mu_{k})+e_{n_{0}+1},$
$\displaystyle=\mu_{n_{0}+1}+\sum_{k=1}^{m^{p}}\beta_{k,n_{0}+1}(X_{k}-\mu_{k})+\sum_{\ell=m^{p}+1}^{n_{0}}\beta_{\ell,n_{0}+1}\left[\sum_{k=1}^{m^{p}}\gamma_{k,\ell}(X_{k}-\mu_{k})+\sum_{k=m^{p}+1}^{\ell}\gamma_{k,\ell}e_{k}\right]+e_{n_{0}+1},$
(25)
$\displaystyle=\mu_{n_{0}+1}+\sum_{k=1}^{m^{p}}\left[\beta_{k,n_{0}+1}+\sum_{\ell=m^{p}+1}^{n_{0}}\gamma_{k,\ell}\beta_{\ell,n_{0}+1}\right](X_{k}-\mu_{k})+\sum_{k=m^{p}+1}^{n_{0}}\left[\sum_{\ell=k}^{n_{0}}\gamma_{k,\ell}\beta_{\ell,n_{0}+1}\right]e_{k}+e_{n_{0}+1}$
$\displaystyle=\mu_{n_{0}+1}+\sum_{k=1}^{m^{p}}\gamma_{k,n_{0}+1}(X_{k}-\mu_{k})+\sum_{k=m^{p}+1}^{n_{0}+1}\gamma_{k,n_{0}+1}e_{k}.$
(26)

Step (25) follows by applying (24). Step (26) follows by applying Equations
(9) and (10). By mathematical induction, we can conclude that Equation (24)
holds for all $n=m^{p}+1,\ldots,m+1$.

##  Appendix C Detailed Derivation of Equation (11)

We consider $W_{k}$ and $\mathcal{J}\subset\mathcal{K}/\\{k\\}$. For
$\mathcal{J}=\emptyset$, we have

$\dfrac{(m-|\mathcal{J}|)!|\mathcal{J}|!}{(m+1)!}[c(\mathcal{J}\cup\\{k\\})-c(\mathcal{J})]=\dfrac{1}{m+1}\gamma_{k,m+1}^{2}\mbox{Var}(W_{k}).$

For $|\mathcal{J}|=m^{\prime}$ with $m^{\prime}=1,\ldots,m$, we have

$\displaystyle\sum_{\left\\{\mathcal{J}:|\mathcal{J}|=m^{\prime}\right\\}}\dfrac{(m-|\mathcal{J}|)!|\mathcal{J}|!}{(m+1)!}[c(\mathcal{J}\cup\\{k\\})-c(\mathcal{J})]$
$\displaystyle=\sum_{\left\\{\mathcal{J}:|\mathcal{J}|=m^{\prime}\right\\}}\dfrac{(m-m^{\prime})!m^{\prime}!}{(m+1)!}\left[\gamma_{k,m+1}^{2}\mbox{Var}(W_{k})+2\sum_{\ell\in\mathcal{J}}\gamma_{k,m+1}\gamma_{\ell,m+1}\mbox{Cov}(W_{k},W_{\ell})\right]$
$\displaystyle=\dfrac{(m-m^{\prime})!m^{\prime}!}{(m+1)!}\Bigg{\\{}{{m}\choose{m^{\prime}}}\gamma_{k,m+1}^{2}\mbox{Var}(W_{k})$
(27)
$\displaystyle~{}~{}~{}~{}+2\sum_{\ell\in\mathcal{K}/\\{k\\}}\left[\sum_{\left\\{\mathcal{J}:|\mathcal{J}|=m^{\prime}\mbox{
and
}\ell\in\mathcal{J}\right\\}}\gamma_{k,m+1}\gamma_{\ell,m+1}\mbox{Cov}(W_{k},W_{\ell})\right]\Bigg{\\}}$
(28)
$\displaystyle=\dfrac{1}{m+1}\gamma_{k,m+1}^{2}\mbox{Var}(W_{k})+2\dfrac{(m-m^{\prime})!m^{\prime}!}{(m+1)!}{{m-1}\choose{m^{\prime}-1}}\sum_{\ell\in\mathcal{K}/\\{k\\}}\gamma_{k,m+1}\gamma_{\ell,m+1}\mbox{Cov}(W_{k},W_{\ell})$
(29)
$\displaystyle=\dfrac{1}{m+1}\gamma_{ki}^{2}\mbox{Var}(W_{k})+\dfrac{2m^{\prime}}{m(m+1)}\sum_{\ell\in\mathcal{K}/\\{k\\}}\gamma_{k,m+1}\gamma_{\ell,m+1}\mbox{Cov}(W_{k},W_{\ell}).$

Step (27) holds because the number of all subsets $\mathcal{J}$ with size
$m^{\prime}$ is ${{m}\choose{m^{\prime}}}$. In Step (28), we shift the order
of sums over $\mathcal{J}$ and $\ell$. Then, Step (29) holds because given
$W_{\ell}$, the number of subset
$\left\\{\mathcal{J}:|\mathcal{J}|=m^{\prime}\mbox{ and
}\ell\in\mathcal{J}\right\\}$ is ${{m-1}\choose{m^{\prime}-1}}$. So, we get
the Shapley value,

$\displaystyle\mbox{Sh}_{W_{k},X_{m+1}}(\boldsymbol{\theta})=\sum_{\mathcal{J}\subset\mathcal{K}/\\{k\\}}\dfrac{(m-|\mathcal{J}|)!|\mathcal{J}|!}{(m+1)!}[c(\mathcal{J}\cup\\{k\\})-c(\mathcal{J})]$
$\displaystyle=$
$\displaystyle\sum_{m^{\prime}=0}^{m}\left[\dfrac{1}{m+1}\gamma_{ki}^{2}\mbox{Var}(W_{k})+\dfrac{2m^{\prime}}{m(m+1)}\sum_{\ell\in\mathcal{K}/\\{k\\}}\gamma_{k,m+1}\gamma_{\ell,m+1}\mbox{Cov}(W_{k},W_{\ell})\right]$
$\displaystyle=$
$\displaystyle\gamma_{k,m+1}^{2}\mbox{Var}(W_{k})+\sum_{\ell\in\mathcal{K}/\\{k\\}}\gamma_{k,m+1}\gamma_{\ell,m+1}\mbox{Cov}(W_{k},W_{\ell}).$

##  Appendix D Derivation and Procedure for BN Learning and Gibbs Sampler

We derive the posterior distribution of BN model parameters
$p(\boldsymbol{\theta}|\mathcal{X})$ and introduce a Gibbs sampling approach
to generate the posterior samples, $\widetilde{\boldsymbol{\theta}}^{(b)}\sim
p(\boldsymbol{\theta}|\mathcal{X})$ with $b=1,2,\ldots,B$ quantifying the
model uncertainty. In Section D.1, we first provide the derivation for
conditional posterior distribution with complete production process data
described in Section 6.1. Considering the situations where we could have some
additional incomplete batch data (e.g., batches in the middle of production or
thrown away at certain production step based on the quality control strategy),
we further extend the Bayesian learning approach to cases with mixing data in
Section D.2. Then, we provide the Gibbs sampling procedure to generate the
posterior samples $\widetilde{\boldsymbol{\theta}}^{(b)}$ with
$b=1,2,\ldots,B$ in Section D.3.

###  D.1 Knowledge Learning for Cases with Complete Production Process Data

Following Section 6.1, we first derive the conditional posterior distribution
for the weight coefficient $\beta_{ij}$,

$\displaystyle
p(\beta_{ij}|\mathcal{X},\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta}_{-ij})\propto\left[\prod_{r=1}^{R}p(x_{j}^{(r)}|x_{Pa(X_{j})}^{(r)})\right]p(\beta_{ij}),$
$\displaystyle\propto\exp\left\\{-\sum_{r=1}^{R}\dfrac{1}{2v_{j}^{2}}\left[(x_{j}^{(r)}-\mu_{j})-\beta_{ij}(x_{i}^{(r)}-\mu_{i})-\sum_{k\in
Pa(j)/\\{i\\}}\beta_{kj}(x_{k}^{(r)}-\mu_{k})\right]^{2}\right.$
$\displaystyle~{}~{}~{}~{}\left.-\dfrac{1}{2\tau_{ij}^{(0)2}}\left(\beta_{ij}-\theta_{ij}^{(0)}\right)^{2}\right\\},$
$\displaystyle\propto\exp\left\\{-\dfrac{1}{2v_{j}^{2}}\sum_{r=1}^{R}\left(\alpha_{i}^{(r)}\beta_{ij}-m_{ij}^{(r)}\right)^{2}-\dfrac{1}{2\tau_{ij}^{(0)2}}\left(\beta_{ij}-\theta_{ij}^{(0)}\right)^{2}\right\\},$
$\displaystyle\propto\exp\left\\{-\dfrac{\beta_{ij}^{2}}{2}\left(\sum_{r=1}^{R}\dfrac{\alpha_{i}^{(r)2}}{v_{j}^{2}}+\dfrac{1}{\tau_{ij}^{(0)2}}\right)+\beta_{ij}\left(\sum_{r=1}^{R}\dfrac{\alpha_{i}^{(r)}m_{ij}^{(r)}}{v_{j}^{2}}+\dfrac{\theta_{ij}^{(0)}}{\tau_{ij}^{(0)2}}\right)\right\\}=\mathcal{N}(\theta_{ij}^{(R)},\tau_{ij}^{(R)2}),$

where
$\theta_{ij}^{(R)}=\dfrac{\tau_{ij}^{(0)2}\sum_{r=1}^{R}\alpha_{i}^{(r)}m_{ij}^{(r)}+v_{j}^{2}\theta_{ij}^{(0)}}{\tau_{ij}^{(0)2}\sum_{r=1}^{R}\alpha_{i}^{(r)2}+v_{j}^{2}}$
and
$\tau_{ij}^{(R)2}=\dfrac{\tau_{ij}^{(0)2}v_{j}^{2}}{\tau_{ij}^{(0)2}\sum_{r=1}^{R}\alpha_{i}^{(r)2}+v_{j}^{2}}$
with $\alpha_{i}^{(r)}=x_{i}^{(r)}-\mu_{i},\quad\textrm{and}\quad
m_{ij}^{(r)}=(x_{j}^{(r)}-\mu_{j})-\sum_{X_{k}\in
Pa(X_{j})/\\{X_{i}\\}}\beta_{kj}(x_{k}^{(r)}-\mu_{k}).$

Second, we derive the conditional posterior distribution for the variance
coefficient $v_{i}^{2}=\mbox{Var}[X_{i}|Pa(X_{i})]$ with $i=1,2,\ldots,m+1$,

$\displaystyle
p(v_{i}^{2}|\mathcal{X},\boldsymbol{\mu},\boldsymbol{v}_{-i}^{2},\boldsymbol{\beta})\propto\left[\prod_{r=1}^{R}p(x_{i}^{(r)}|x_{Pa(X_{i})}^{(r)})\right]p(v_{i}^{2})$
$\displaystyle\propto$