A second, but especially challenging task would be to attempt the scoring of active vs. inactive compounds. This task is rarely accomplished due to the interplay of several factors that are external to the nano-particle.
We focus on the search algorithms that have been designed to efficiently predict the docking pose. The process of docking a target and a small molecule falls into the class of NP-hard problems due to the non countability of the number of possible poses. Hence, search becomes the approach to identify candidate solution and improving on those. In the literature, search methods can be classified into deterministic (also referred to as systematic) and stochastic [1]}. Systematic search methods sample within the binding molecule search space at predefined intervals and are deterministic. Within this class, we can still differentiate between exhaustive, fragmentation
or conformational ensemble methods. The main difference between them is in the approach they take to deal with the binding molecule flexibility. In exhaustive search methods, for example, the docking is performed by systematically rotating all possible rotatable bonds in the binding molecule at a given interval. The drawback of this family of approaches is computational in nature as the number of possible combinations to consider goes with the number of rotatable bonds in the the binding molecule. An common exhaustive sampling method is Glide presented in [2]}, [3]}. Fragmentation represents an attempt to improve on the computational efficiency by incrementally forming binding over fragments that the binding molecule is divided into. An approach that relies on fragmentation is FlexX [4]}. Finally, in conformational ensemble methods, the binding molecule flexibility is represented by rigidly docking an ensemble of pre-generated conformations, thus improving the approach efficiency since using this approach removes the computational cost due to the exploration of the conformational space.
In stochastic algorithms, the binding molecule orientations and conformations are sampled by making changes to the molecule that are informed by random score values iteratively generated by a random algorithm. The orientation, conformation change is then treated as a incumbent that is accepted or rejected according to an algorithm-dependent criterion. The advantage of stochastic algorithms is that they can generate large ensembles of molecular conformations and explore a broader range of the
energy landscape increasing the probability of finding a global energy minimum.
However, this comes at computational cost. Examples are genetic algorithm, Monte
Carlo, ant colony optimization (ACO) and tabu search methods. GOLD in [5]} uses a genetic algorithm, DockVision [6]} is an example of docking program using Monte Carlo stochastic method where the probability to accept a random change is calculated by using the Boltzmann probability function. An example of ACO-based approach is PLANTS [7]}, while PSI-DOCK uses a tabu search [8]}.

Data Sets.
Search algorithms usually rely on data sets to retrieve potential components of the nano-particle to be assembled. In these datasets, the crystal structures of the complexes are specified using different microscopy technologies (with potentially different resolutions). Prior to any docking study, the virtual compounds database that will be screened must be carefully selected and prepared. This compound collection, often referred to as the virtual library, can encompass up to millions of compounds. Already prepared virtual libraries can be used, but users can also generate their own. Commercial databases represent an important source of compounds for virtual screening (often containing more than 1 million molecules). Suppliers commonly offer free access to the files containing molecules structures in several formats (2D or 3D). These databases undergo frequent
updates, new products being added while other being either removed or out of stock. Virtual screening libraries constructed from these databases should ideally be prepared when the whole virtual screening protocol is settled and ready to be used.
A particular category of databases are bioactivity databases providing knowledge on biological targets and their modulation mechanisms. These data are frequently used for data-set benchmarking in the context of docking protocols design prior to prospective virtual screenings and
to construct predictive models of activity [9]}, [10]}, [11]}. Examples of databases in this category are ZINC [12]} or PubChem [13]}.
Finally, natural products data bases are available. Natural products were the first drugs ever used and have always been a source of drugs. In a recent retrospective study, it has been reported that, between January 1, 1981 and September 30, 2019, 23.5% of all new approved drugs and 33.6% of new approved small molecules drugs were natural products or derivatives of natural products [14]}. The chemical space covered by natural products is quite dissimilar to the one occupied by synthetic drug-like compounds [15]} and natural products are believed to constitute promising starting point for drug discovery [16]}. Hence, these data sets can represent a source for virtual screening libraries. Numerous databases of this type are available; please refer to [17]} and [18]}).

Scoring Functions.
There are three important applications of scoring functions in molecular docking [19]}.
The first of these is
the determination of the binding mode and site of small molecule to its target.
Specifically, once a target has been defined, molecular docking generates hundreds of thousands of possible binding orientations/conformations for the small molecule (e.g., ligand, peptide) at the active site around the target molecule. The scoring function is used to rank such small molecule orientations/conformations by evaluating the binding tightness of each of the candidate complexes. Ideally, the scoring function would rank the highest the experimentally determined binding mode. Given the determined binding mode, scientists would be able to gain a deep understanding of the molecular binding mechanism and to further design an efficient drug by modifying either the target or the small molecule. It is important to highlight that, instead of scoring functions, other computational methodologies based on molecular dynamics or Monte Carlo simulations could be used to model the dynamics of the binding process thus resulting in a more accurate prediction of binding affinity. However, these models result in computationally prohibitive free energy calculations, ultimately resulting impractical for the evaluation of large numbers of molecular complexes. As a result the application of high fidelity simulation techniques is generally reduced to predicting binding affinity in small nano-particles [20]}.
The second application of a scoring function, related to the first, is the prediction of the absolute binding affinity between the active compounds. This is particularly important in lead optimization, i.e., the process to improve the tightness of binding for low-affinity hits or lead compounds that have been identified. In this phase, an accurate scoring function can greatly increase the optimization efficiency and save costs by computationally predicting the binding affinities between the conformed small molecule and the target before the much more expensive synthesis and experimental steps.
The third application, is related to the foundational task of structure-based design, that fundamentally attempts to identify the potential drug hits/leads for a given target by searching a large compound database, this is commonly referred to as virtual database screening. A reliable scoring function should be able to associate higher rank to known binders following their binding scores during database screening. In fact, due to the expensive cost of experimental screening and sometimes unavailability of high-throughput assays, virtual database screening has played an increasingly important role in drug discovery.
Classical scoring functions can be classified into three groups: forcefield, knowledge-based and empirical [21]}, [22]}, [23]}, [24]}, [25]}. An alternative to the classical approach to the design of scoring functions, a non-parametric machine-learning approach can be taken to implicitly capture the binding interactions that are hard to explicitly model by classical approaches in a computationally efficient way. By not imposing a particular functional form for the function, the collective effect of intermolecular inter-actions in binding can be directly inferred from experimental data, which should lead to scoring mechanisms that are characterized by greater generality and, consequently, prediction accuracy.
BIOMANUFACTURING HYBRID MODELING AND ANALYTICS
Driven by the critical challenges and needs of biopharmaceutical manufacturing, we propose the probabilistic knowledge graph (KG) hybrid model characterizing the spatial-temporal causal interdependencies of critical process parameters (CPPs) and critical quality attributes (CQAs)  [26]}, [27]}.
To faithfully represent underlying bioprocessing mechanisms, this hybrid model can capture the important properties, including
nonlinear reactions, partially observed state, and nonstationary dynamics.
It has time-varying kinetic coefficients (such as molecular reaction rates) with uncertainty representing batch-to-batch variation.
To avoid the evaluation of intractable likelihood, we further investigate a computational Bayesian inference approach, called Approximate Bayesian Computation (ABC), to efficiently approximate the posterior distribution.
Then, assisted by the Bayesian KG, accounting for inherent stochasticity and model uncertainty, we present interpretable risk, sensitivity, and predictive analyses.
This study can support: (1) biomanufacturing stochastic decision process (SDP) mechanism online learning; (2) bioprocess soft sensor monitoring (such as tracking latent metabolic state associated with therapeutic cell function and product quality); and (3) optimal and robust process control.

Illustration Example: mRNA lipid nanoparticle formulation process.
Here we use mRNA lipid nanoparticle (mRNA-LPN) formulation process as an illustration example; see Figure REF .
Lipid-based formulations have rapidly emerged as carriers in nucleic-acid therapies (i.e., DNA and RNA)
to
increase cellular uptake, support RNA delivery, and improve drug stability. The mRNA lipid nanoparticle (LNP) formulation utilizes microfluidic or T-junction mixing to rapidly combine an ethanol phase containing hydrophobic lipids and an aqueous phase containing mRNA in a buffer. Then, the self-assembly of LNPs with mRNA is driven by the hydrophobic and electrostatic force field that is influenced by the design of lipids, the selection of LNP formulation and CPPs, as well as the phases of mRNA-LNP complex. The pH-responsive ionizable cationic lipids have the surface charge modulated, controlling the efficient binding with the oppositely charged polynucleotides, which will support self-assembly, influence mRNA-LNP thermodynamic stability, prolong the circulation time of mRNA-LNP complexes, facilitate endosomal escape and mRNA release, and increase their therapeutic.
Therefore, the dynamics and variations of mRNA-LNP formulation process output trajectory depends on complex interactions of (1) the design of lipids and (2) CPPs (e.g., flow rate ratio, total flow rate, temperature, lipid choices, buffer choices), which impacts on the quality of LNP delivery system specified CQAs including a) RNA integrity level; b) species composition/concentrations/phases, particle size distribution, zeta potential, surface charge; and c) mRNA-LNP, bound/unbound mRNA.
In Section REF , we present the macro-phases or operation units involved in the manufacturing of a bioproduct, while Sections REF  and REF  illustrate novel methods developed for the efficient modeling, analytics, and control of such process when targeting novel synthetic products (e.g., variations of the mRNA-LPN).
Operation Unites
of a Biomanufacturing process
Biopharmaceutical manufacturing process is crucially important to determine product quality and productivity.
It typically includes the main unit operations listed below.
(1)
Fermentation and Drug Substance Synthesis: Living organisms (e.g., cells, yeasts) are mixed with appropriate medium and enzymes under carefully controlled conditions to grow and synthesize the target drug substance.
The byproducts or unwanted impurities are also generated at the meantime, which impacts downstream purification operation and cost.
During different growth and production phases of cell and yeast life cycle, different media compositions and feeding strategies are used to improve the productivity and reduce the waste generation.
(2)
Centrifugation(s): The centrifuge device is used for separation of particles, e.g., cells, subcellular organelles, viruses, large molecules such as proteins, from a solution according to their size, shape, density, viscosity of the medium and rotor speed, during which bulk of impurities would be removed.
(3)
Chromatography and Purification: As the mixture of solutes flows through a packed resin bed, specific solutes are separated as they are bound. Chromatography serves as the most critical part for purification, and it usually determines the purity of product. Since removing more impurities often results in removing more drug substance during chromatography step, there is often a trade-off between productivity and purity.
(4)
Filtration: It is applied at several stages for capture (i.e., concentrate the product), intermediate purification, and polishing (i.e., eliminates trace contaminants and impurities) purpose.
(5)
Formulation and Filling: To maintain the safety and efficacy of the drug substance during the storage, delivery, and facilitate the patient absorption of active pharmaceutical ingredient (API), the purified drug substance is usually formulated with carefully selected excipients and nanoparticle carriers into stable drug products and filled into dose containers [28]}.
(6)
Freeze Drying: It is used to stabilize bio-drugs through removing water or other solvents from the frozen matrix and converting the water directly from solid phase to vapor phase through sublimation. Freeze drying is critical for immobilizing the bio-drug product in storage and delivery, as the kinetics of most chemical and physical degradation reactions are significantly decreased.
(7)
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 [29]}.
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 [26]} and [27]} 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 [26]}, [33]}.
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 [34]}, [35]} 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 [35]},
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 [33]}.
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 [27]}, [39]} .
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 [39]}, 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 [41]}.
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) [42]}, [43]},
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 [41]},
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 [45]}, 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/.
