Quality Assurance/Control (QA/QC): QA and QC are performed to ensure the quality of the raw material selection, the production process, and the final bio-drug product.

In sum, Step (1) belongs to upstream fermentation and drug substance synthesis, Steps (2)–(4) belong to downstream purification, and Steps (5)–(7) are for finished drug filling/formulation, freeze drying, and product quality control testing.
There are interactions of hundreds of factors at different productions steps impacting drug quality, yield and production cycle time.
These factors can be divided into critical process parameters (CPPs) and critical quality attributes (CQAs) in general; see the definitions of CPPs/CQAs in ICH-Q8R2 [1]}.
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.

KG hybrid model for Biomanufacturing Stochastic Decision Process
The probabilistic KG hybrid model proposed in [2]} and [3]} can provide the risk- and science-based understanding of underlying stochastic decision process (SDP) mechanisms for controlled production processes. The input-output relationship in each step is modeled by a hybrid (“mechanistic+statistical") model that can leverage the prior knowledge on biophysicochemical mechanisms from existing mechanistic models and further advance scientific learning from process data.
Specifically, at any time \(t\) , the process state \(s_t\)  (such as CQAs) composed of observable and latent state variables, i.e., \(s_t=(x_t,z_t)\)  (e.g., particle size distribution, RNA sequence and integrity level), and CPPs action \(a_t\)  (e.g., temperature, mixing flow rate, pH) interactively influence on the dynamics and variations of output trajectories (e.g., mRNA-LNP formulation and self-assembling processes).
Given the existing nonlinear ODE/PDE-based mechanistic model (such as biomolecular dynamics, thermodynamics, molecular interactions of mRNA and lipids which can affect the RNA integrity), represented by
\({\mbox{d}s}/{\mbox{d}t} = f\left(s,a; \beta \right),\) 
by applying the finite difference approximations on derivatives, we construct the hybrid model for state transition,
\(x_{t+1} &=& x_t + \Delta t \cdot f_x(x_t,z_t,a_t; \beta _t) + e^{x}_{t+1},\nonumber \\z_{t+1} &=& z_t + \Delta t \cdot f_z(x_t,z_t,a_t; \beta _t) + e^{z}_{t+1}, \nonumber \) 
with unknown random kinetic coefficients \(\beta _t\in \mathbb {R}^{d_\beta }\)  (e.g., particle clustering rates) accounting for batch-to-batch variation. The function structures of \(f_x(\cdot )\)  and \(f_z(\cdot )\)  are derived from \(f(\cdot )\)  in the mechanistic models.
The residual terms are modeled by multivariate Gaussian distributions \(e_{t+1}^{x} \sim \mathcal {N}(0,V^{x}_{t+1})\)  and \(e_{t+1}^{z} \sim \mathcal {N}(0,V^{z}_{t+1})\)  with zero means and covariance matrices \(V^{x}_{t+1}\)  and \(V^{z}_{t+1}\)  by applying CLT.
The kinetic coefficients \(\beta _t\)  can change cross different phases of bioprocess accounting for the fact that the process dynamics can be time-varying.
The statistical residual terms \(e_t=(e_t^x,e_t^z)\)  allow us to account for the impact from bioprocess noises, raw material variations, ignored CPPs, sensor measurement errors, and other uncontrollable factors (e.g., contamination) occurring at any time step \(t\) .
The probabilistic KG
of integrated biomanufacturing process can be visualized by a directed network as shown in Figure REF .
The observed state variables \(x_t\)  and the latent state variables \(z_t\)  are represented by solid and shaded nodes respectively. The directed edges represent causal interactions.
At any time period \(t+1\) , the process state output node
\(s_{t+1}=(x_{t+1},z_{t+1})\)  depends on its parent nodes: \(s_{t+1}=f(Pa(s_{t+1});{\theta }_t)\) 
with \(Pa(s_{t+1})=(s_t,a_t,e_{t+1})\)  and model parameters denoted by \(\theta _t\) . To represent the underlying controlled SDP, we create a policy augmented KG network by including additional edges: 1) connecting state \(s_t\)  to action \(a_t\)  representing the causal effect of the policy, \(a_t=\pi _t (s_t|\phi )\)  specified by parameters \(\phi \) ; and 2) connecting actions and states to the immediate reward \(r_t(s_t,a_t)\)  (e.g., cost and RNA production).
This KG network models how the effect of current state and action, \(\lbrace s_t,a_t\rbrace \) , propagates through mechanism pathways impacting on the output trajectory and the accumulated reward.
<FIGURE>
Bayesian KG representing risk- and science-based mechanism understanding. Leveraging the information from existing mechanistic models and heterogeneous online/offline measurements, the probabilistic KG hybrid model represents the understanding of bioprocess mechanisms. It allows us to inference the latent state (e.g., RNA and nanoparticle binding strength, mRNA-LNP folding structure),
which can support biomanufacturing online monitoring and real-time release.
Correctly quantifying all sources of uncertainty can facilitate optimal learning, guide risk reduction, and support robust control.
Therefore, the Bayesian KG, accounting for inherent stochasticity and model uncertainty, can be created and used to support integrated bioprocess risk, sensitivity, and predictive analyses.
Given finite real-world data of the partially observed bioprocess trajectory with size \(m\) , denoted by \(\mathcal {D}_m=\lbrace \tau _x^{(i)}: i=1,2,\ldots ,m\rbrace \)  with \(\tau _x\equiv (x_1,a_1,\ldots ,x_H,a_H,x_{H+1})\) , the model uncertainty is quantified by a posterior distribution,
\(p(\theta |\mathcal {D}_m) \propto p(\theta )p(\mathcal {D}_m|\theta )=p(\theta )\prod _{i=1}^m p\left(\left.\tau _x^{(i)} \right|\theta  \right)\nonumber \) 
where \(p(\theta )\)  represents the prior distribution. Since the likelihood evaluation of the KG hybrid model, with high fidelity to capture the critical properties of bioprocessing, is intractable, i.e., \(p(\tau _x| \theta )= \int \cdots \int p(\tau |\theta ) d z_1 \cdots d z_{H+1}\) , approximate Bayesian computation sampling with Sequential Monte Carlo (ABC-SMC) is developed to approximate the posterior distribution [2]}, [5]}.
In the naive ABC implementation, we draw a candidate sample from the prior \(\theta  \sim p(\theta )\)  and then generate a simulation dataset \(\mathcal {D}^\star \)  from the hybrid model.
If the simulated dataset \(\mathcal {D}^\star \)  is “close" to the observed real-world observations \(\mathcal {D}_m\) , we accept the sample \(\theta \) ; otherwise reject it.
The accept rate can be very low since it is very computationally expensive to match random trajectories from complex bioprocesses especially under the situations with high stochastic and model uncertainties.
The ABC-sequential Monte Carlo (ABC-SMC) methods [6]}, [7]} can improve the sampling efficiency through: (1) generating candidate samples from updated posterior approximates by using sequential importance sampling (SIS); and (2) matching “well-designed" summary statistics, denoted by \({\eta }(\mathcal {D})\) , instead trajectory observations.
Following the spirit of the auxiliary likelihood-based ABC [7]},
we create a linear Gaussian dynamic Bayesian network (LG-DBN) auxiliary model and derive summary statistics for ABC-SMC that can accelerate online Bayesian inference on KG hybrid models [5]}.
This simple LG-DBN auxiliary model, in conjunction with SIS, can capture critical biophysicochemical interactions and variations of bioprocess trajectory, ensure the computational efficiency, and enable high quality of inference. Therefore, the proposed LG-DBN auxiliary likelihood-based ABC-SMC approach can support process soft sensor monitoring, facilitate mechanism online learning, and guide robust process control.

Interpretable prediction and sensitivity analysis.
Given process model parameters and policy parameters, denoted by \((\theta ,\phi )\) , the spatial-temporal interdependencies of the bioprocess SDP trajectory \(\tau =(s_1,a_1,\ldots ,s_H,a_H,s_{H+1})\)  is quantified by the joint distribution,
\(p(\tau |\theta ,\phi )=p(s_1)\prod ^{H}_{t=1}p(s_{t+1}|s_t,a_t;\theta )\pi _{\phi }(a_t|s_t)\) ,
which depends on underlying bioprocess mechanisms, sensor network design, and data collection strategies.
For each batch of production, given inputs denoted by \(\mathbf {X}\)  (e.g., mRNA sequence, lipid design), we can predict any intermediate or final outputs, denoted by \(\mathbf {Y}\)  (e.g., CQAs of mRNA-LNP) by using Bayesian KG. The prediction risk can be quantified by the posterior predictive distribution,
\(P(\mathbf {Y}| \mathbf {X}) = \int P(\mathbf {Y}| \mathbf {X}, \theta ) p(\theta |\mathcal {D}_m) d \theta \) , accounting for both stochastic and model uncertainties.
We create a Shapley value (SV)-based sensitivity analysis scheme on the Bayesian KG, called “KG-SV", to support backward root cause analysis and forward interpretable predictive analysis [3]}, [11]} .
Since the proposed Bayesian KG-SV can faithfully account for bioprocess causal interdependencies and biophysicochemical interactions, it can correctly assess the effect from each input variation (such as RNA virus mutation), as well as the impact of each source of stochastic and model uncertainties on the prediction risks.
The criticality assessment of input factors is
based on the estimated values and estimation uncertainties of interpretable KG hybrid model parameters – such as mechanism pathways in the bioprocess KG from inputs to output (i.e., biomolecular reaction rates, mRNA lipid nanoparticle clustering kinetic parameters).
Since model uncertainty can be efficiently reduced by most “informative" data collection and SDP inherent stochasticity can be controlled by decision making, the Bayesian KG based risk, sensitivity, and predictive analyses can identify bottlenecks, guide optimal learning based data collection, and enhance biomanufacturing process CPPs/CQAs specifications
for QbD.

Reinforcement Learning for Bioprocess Design and Control
The proposed Bayesian KG built in conjunction with reinforcement learning (RL) can support long-term prediction and guide interpretable, robust, and optimal decision making.
Given any feasible policy specified by parameters \(\phi  \in \mathbb {C}\) , i.e., \(a_t=\pi _t (s_t|\phi )\) , the optimization of the policy \(\pi \)  is to maximize the expected accumulated reward,
\(J(\phi ) \equiv \mbox{E}_{\theta \sim p\left(\theta |\mathcal {D}\right)} \left[ \mbox{E}_{\tau \sim p(\tau |{\theta })}\left[\left.\sum _{t=1}^{H+1} r_t(s_t,a_t)\right|\pi ,{{\theta }}\right] \right],\) 
accounting for bioprocess inherent stochasticity and model uncertainty, where \(\mathbb {C}\)  is a feasible region.
At any \(k\) -th iteration, we can use the policy gradient to solve the optimization,
\(\phi _{k+1} = \Pi _{\mathbb {C}}\left(\phi _k + \eta _k \nabla {J}\left(\phi _k\right)\right),\) 
where \(\eta _k\)  is a suitable
stepsize and \(\Pi _{\mathbb {C}}\)  is a projection onto \(\mathbb {C}\) .

Reinforcement learning scheme on Bayesian KG.
We propose model-based RL scheme on the Bayesian KG [11]}, which can provide an insightful prediction on how the effect of input factor propagates through bioprocess mechanism pathways and impacts on the outputs.
It can find control policies that are interpretable and robust against heterogeneous model uncertainty, and overcome the key challenges of biopharmaceutical manufacturing, i.e., high complexity, high uncertainty, and very limited process data.
To support real-time control for complex biomanufacturing processes, we provide a provably convergent
stochastic policy gradient optimization and it is computationally efficient through reusing computations associated with similar input-output mechanism pathways.

Hybrid Model Likelihood Ratio based Historical Observation Reuse.
Since each experiment run is very computationally expensive especially for multi-scale bioprocess hybrid model, we propose KG assisted multiple important sampling (“KG-MIS") to accelerate policy gradient optimization [13]}.
Basically, we can select and reuse the “most relevant" historical trajectories, improve policy gradient estimation, and accelerate the search for the optimal robust policy.
For high dimensional SDP,
this study can selectively reuse historical trajectories
having similar underlying distributions with that of target SDP and improve the estimation of policy gradient.
In classical policy gradient (PG) approach, at any \(k\) -th iteration,
the sample average approximation (SAA) is used to estimate the gradient based on \(n\)  new trajectories generated,
\(\nabla \widehat{J}^{PG}(\theta _k)=\frac{1}{n}\sum ^{n}_{j=1}g\left(\tau ^{(k,j)}|\theta ^{(k,j)},\phi _k \right)\)  with the scenario gradient \(g(\tau |\theta ,\phi )=\nabla _{\theta }\mbox{E}_{\tau }\left[\left.\sum _{t=1}^H r_t(s_t,a_t)\right|\theta ,\phi \right]\) ,
where \(\theta ^{(k,j)} \sim p(\theta |\mathcal {D}_{m})\)  and \(\tau ^{(k,j)} \sim p(\tau |\theta ^{(k,j)},\phi _k)\) .
The target SDP mixture distribution
\(p_k(\tau ) = \frac{1}{n} \sum _{j=1}^np(\tau |\theta ^{(k,j)},\phi _k)\)  accounts for both process stochastic and model uncertainties.
Motivated by the studies on multiple important sampling (MIS) [14]}, [15]},
we create a KG-MIS policy gradient unbiased estimator,
\(\nabla \widehat{J}^{MIS}(\phi _k)=\frac{1}{n|U_k|}\sum _{i\in U_k}\sum ^n_{j=1}f_{k}\left(\tau ^{(i,j)}\right)g\left.\left(\tau ^{(i,j)}\right|\theta ^{(k,j)}, \phi _k\right)~~~\mbox{with}~~~ f_{k}(\tau ) = \frac{p_k(\tau )}{\sum _{i \in U_k} p(\tau |\theta _i,\phi _i)/|U_k|}.\) 
Since an inappropriate selection of reuse set \(U_k\)  can lead to the inflated estimation variance of policy gradient, we propose a variance reduction based experience replay criteria [13]},
which can automatically select the most relevant historical trajectories generated under different decision policies \(\phi \)  and model parameters \(\theta \)  from different posterior distributions resulting from online learning and process control.
We prove that the proposed approach is asymptotically convergent and show it significantly outperforms classical policy gradient approach. Furthermore, we extend the proposed KG-MIS framework so that it can select and reuse the most relevant partial trajectories from historical observations [17]}, i.e., the reuse unit is defined based on state-action transition \((s,a,s^\prime )\) . This study can allow us to flexibly integrate and leverage the relevant information collected from different production lines and facilitate personalized bio-drug manufacturing.

CHALLENGES AND OPPORTUNITIES FOR OPERATIONS RESEARCH
Increased flexibility, that is necessary to achieve personalized products and manufacturing, should be considered early on as integral part of product design. In fact, for achieving a “full circle”, not only the manufacturing technology needs to be flexible, but also the drug design and the process control need to support it. Novel operations research approaches and simulation platforms can substantially improve the performance of CMO allowing for larger variability of products with potentially small volume per variant capitalizing upon single use/disposable technologies.
Drug discovery is positively impacted by optimization methods. These should embed scarce data and low fidelity physical models characterizing the existing understanding of bioprocess mechanisms. In fact, gray-box search methods are a very active field of research and we believe drug design represents a leading opportunity for further development. The expert-based framework is an example of such approaches, but more efforts are necessary.
Novel simulation methodologies are necessary for analyzing end-to-end biomanufacturing processes and supporting interoperability. Among new upcoming techniques are hybrid modeling, data integration, risk management, and interpretable robust process control. However these are an example and more development in the ares is required.

Acknowledgements
We acknowledge Hua Zheng and Keqi Wang (Northeastern University) for their help on the paper preparation. This work was partially supported by the grants NSF#2046588-2007861 (PI Pedrielli) and NIIMBL PC4.1-206 (PI Xie).

AUTHOR BIOGRAPHIES
WEI XIE is an assistant professor in MIE at Northeastern University.
Her research
interests include interpretable Artificial Intelligence (AI), machine learning (ML),
computer simulation, digital twin, data analytics, and stochastic optimization for cyber-physical system risk management, learning, and automation.
Her email address is w.xie@northeastern.edu. Her website is http://www1.coe.neu.edu/\(\sim \) wxie/

GIULIA PEDRIELLI is Assistant Professor for the School of Computing and Augmented Intelligence at Arizona State University. Her expertise is in the design and analysis of sampling based algorithms, Monte Carlo methods, reinforcement learning, and dynamic simulation. Applications range from molecular design, process design and control, testing and verification of Cyber Physical Systems. Her email address is giulia.pedrielli@asu.edu, her webpage can be found at https://www.gpedriel.com/.
