# Digital Twin Calibration for Biological System-of-Systems:
Cell Culture Manufacturing Process

Fuqiang Cheng
Northeastern University

& Wei Xie
Northeastern University
& Hua Zheng
Northeastern University
Corresponding author. Email<EMAIL_ADDRESS>
###### Abstract

Biomanufacturing innovation relies on an efficient design of experiments (DoE)
to optimize processes and product quality. Traditional DoE methods, ignoring
the underlying bioprocessing mechanisms, often suffer from a lack of
interpretability and sample efficiency. This limitation motivates us to create
a new optimal learning approach that can guide a sequential DoEs for digital
twin model calibration. In this study, we consider a multi-scale mechanistic
model for cell culture process, also known as Biological Systems-of-Systems
(Bio-SoS), as our digital twin. This model with modular design, composed of
sub-models, allows us to integrate data across various production processes.
To calibrate the Bio-SoS digital twin, we evaluate the mean squared error of
model prediction and develop a computational approach to quantify the impact
of parameter estimation error of individual sub-models on the prediction
accuracy of digital twin, which can guide sample-efficient and interpretable
DoEs.

_Keywords_ Digital Twin Calibration $\cdot$ Sequential Design of Experiments
$\cdot$ Cell Culture $\cdot$ Biomanufacturing Process $\cdot$ Biological
System-of-Systems $\cdot$ Linear Noise Approximation

##  1 INTRODUCTION

To support interpretable predictions and optimal control of biomanfuacturing
processes, in this paper, we develop a digital twin calibration approach to
guide a sequential design of experiments (DoEs). It can sample efficiently
improve model fidelity of multi-scale bioprocess mechanistic model for
Biological System-of-Systems (Bio-SoS) (Zheng et al., 2024) that characterizes
the fundamental bioprocessing mechanisms and causal interdependence from
molecular- to cellular- to macro-kinetics. Even though this study is motivated
by cell culture process, the proposed optimal learning approach can be
extended to calibrate the digital twin of general Bio-SoS with modular design.
In specific, the dynamics and variations of cell culture process is
characterized by a multi-scale mechanistic model composed of sub-modules or
modules: (1) a single cell mechanistic model characterizing each living cell
behaviors and their interactions with environment; (2) a metabolic shift model
characterizing the change of cell metabolic phase and behaviors as a response
to culture condition changes and cell age; and (3) macro-kinetic model of a
bioreactor system composed of many living cells under different metabolic
phases.

The benefits of considering the Bio-SoS mechanistic model with modular design
include: a) support flexible manufacturing through assembling a system of
modules to account for biomanufacturing processes under different conditions
and with various inputs; and b) facilitate the integration of experimental
data from different production systems, such as 2D culture and 3D culture for
Induced pluripotent stem cell (iPSCs). By accounting for the structure
property of the Bio-SoS mechanistic model, we can quantify how the model
uncertainties or estimation errors of different modules interact with each
other and propagate through the reaction mechanism pathways to outputs, which
can support the interpretable design of experiments to enhance model fidelity
and improve predictions.

Before describing the proposed calibration approach, we review the key aspects
for digital twin calibration. First, the model uncertainty quantification
approaches for digital twins and simulators can be divided into two main
categories: Bayesian and frequentist approaches (Corlu et al., 2020). Bayesian
approaches treat unknown model parameters as random variables and quantify our
belief by using posterior distributions. It involves specifying prior
distributions for model parameters and updating these distributions based on
the information collected from observed data by applying Bayes’ theorem. On
the other hand, frequentist approaches estimate model parameters using
traditional statistical methods, such as maximum likelihood estimation, to
find parameter values that maximize the likelihood of observed data. The model
estimation uncertainty can be quantified by asymptotic normal approximation or
bootstrap approach. These frequentist estimation approaches could be easier
and faster compared with Bayesian inference approaches since there may not
exist conjugate priors and posterior sampling update can be computationally
expensive.

Second, for the propagation of input and model uncertainty to simulation
outputs, to save computational budget and time, a metamodel is often used to
approximate the response surface, especially for complex stochastic systems,
since easy experiment can be expensive. In the classical calibration approach
(Kennedy and O’Hagan, 2001), the Gaussian processes (GPs) are used to
characterize the mean response and model discrepancies. Although GP metamodel
is often used to propagate the input to output and further support the digital
twin calibration, it’s hard to interpret and difficult to leverage prior
knowledge of the real biomanufacturing system mechanisms. To enhance
interpretability and sample efficiency, the mechanistic models can be employed
to construct the response surface in our study.

Third, the selection of calibration criteria plays a critical role to guide
the experiments and most informative data collection to improve model fidelity
and support digital twin calibration. The existing calibration approaches
often focus on system design. Mean response and Mean Squared Error (MSE) are
widely used as a criteria for assessing model prediction accuracy. Tuo and
Jeff Wu (2016) shows Kennedy’s method, which models the mean output as a GP,
may lead to asymptotically $L_{2}$-inconsistent. They modify the criteria into
$L_{2}$ norm of the discrepancy between real and digital system outputs and
propose a $L_{2}$-consistent calibration method, which can ensure an optimal
convergence rate.

To support process control and facilitate biomanufacturing process automation,
in this paper, we aim to calibrate the Bio-SoS mechanistic model to guide a
sequential DoE and improve model fidelity and prediction accuracy of digital
twin. We first employ the Maximum Likelihood Estimation (MLE) method to
estimate model parameters and utilize bootstrap techniques to quantify
estimation uncertainty across different modules. Then, to optimize the digital
twin calibration policy that will be used to sequentially guide the
experiments and the most informative data collection, we derive the gradient
of prediction MSE and use the steepest descent to efficiently search for the
optimal policy parameters, which makes the learning process more
interpretable. Since it is challenging to find a closed form solution of the
mechanistic model with stochastic differential equations (SDEs) form, we
utilize the linear noise approximation (LNA) for uncertainty propagation and
constructs a surrogate model for MSE. Thus, in an online setting, a
calibration policy is iteratively updated by using a stochastic gradient
method where the gradient of MSE with respect to the calibration policy
parameter follows the backward direction of the uncertainty propagation. The
proposed calibration approach accounts for the Bio-SoS mechanism structure.
That means it quantifies Bio-SoS module error interaction and their
propagation through mechanistic pathways to output, which guides the
sequential design of experiments to efficiently improve model fidelity and
prediction accuracy.

In sum, we develop an interpretable and sample efficient optimal learning
approach to support digital twin calibration for a multi-scale bio-process
mechanistic model so that the digital twin of the cell culture process can
improve process prediction and support optimal control. Even though this paper
focuses on cell culture, the proposed calibration approach can be extendable
to general biomanufacturing systems. The key contributions of the proposed
calibration approach and the benefits are summarized as follows.

  * •

We proposed a new DoE method for calibrating a multi-scale bioprocess
mechanistic model by updating the parameters using closed-form gradients that
incorporate the mechanistic information.

  * •

We assess the MSE of the process output prediction and use a LNA-based
metamodel along with Euler’s method to estimate how model uncertainty
propagates through Bio-SoS mechanism pathways and impacts on the output
prediction accuracy. This approach can advance our understanding on how the
errors in individual module parameters affect the overall prediction accuracy
of digital twin model.

  * •

Built on the LNA, we further develop a gradient-based policy optimization to
guide most informative experiments and support sample-efficient optimal
learning.

The organization of the paper is as follows. In Section 2, we describe the
multi-scale mechanistic model. Then, we propose the novel calibration method
for the multi-scale model in Sections 3, which includes two main procedures:
model inference and policy update. A stochastic gradient method is developed
for the model inference in Section 4. The calibration policy gradient
estimation and update method is further developed in Section 5. Then we
empirically validate the proposed digital twin calibration approach in Section
6 and conclude this paper in Section 7.

##  2 Problem Description and Bio-SoS Mechanistic Model

In this study, we explore a multi-scale bioprocess mechanistic model designed
to elucidate the interdependencies that exist across molecular, cellular, and
macroscopic levels within cell culture processes. This comprehensive model
integrates the complexities of biological interactions within and between
these scales and extends applicably to a broader range of biological systems,
referred to as Biological System-of-Systems (Bio-SoS). A detailed summary of
the variable notations used within our model is provided in Table 1.

The core of our model is based on the structured Markov chain, capturing both
the dynamics and variability inherent in cellular processes through a state
transition probabilistic model. Conceptually, each individual cell operates as
a complex system, and collectively, numerous cells at various metabolic phases
within the bioreactor form an intricate system of systems. This interaction is
depicted in Figure 1, where, for clarity, only three distinct metabolic phases
are illustrated. Cells interact with each other by altering their environment,
notably through nutrient uptake and the production of metabolic wastes, and in
turn, respond to these environmental changes.

$\tilde{\mathbf{s}}$ |  Single cell state |  $\boldsymbol{s}$ |  Macro-state
---|---|---|---
$r$ |  Reaction indices ($r=1,2,\ldots,R$) |  $Z_{t}$ |  Cell phases for single cell ($Z_{t}=i,i=0,1,\ldots,I$)
$X_{i,t}$ |  Cell density in $i$-th phase ($i=0,1,\ldots,I$) |  $m$ |  State component indices ($m=1,2,\ldots,M$)
$t$ |  Time ($t=0,1,2,\ldots,T$) |  $k$ |  Calibration iteration number ($k=1,2\ldots,K$)
$\mathbb{P}_{ii^{\prime}}$ |  Transition probability from $i$\- to $i^{\prime}$-th phase |  $\boldsymbol{\beta}$ |  Parameters for phase shift model
$\boldsymbol{\alpha}_{i}$ |  Model parameters for cells in phase $i$ |  $\boldsymbol{\theta}^{c}$ |  True model parameters
$\widehat{\boldsymbol{\theta}}_{k}$ |  Estimated model parameters at iteration $k$ |  $\boldsymbol{\omega}_{k}$ |  Policy parameters
$\lambda$ |  Learning rate for updating parameter |  $\gamma_{k}$ |  Learning rate for updating policy
$\boldsymbol{u}_{k}$ |  Virtual model update function |  $\tilde{\boldsymbol{\theta}}_{k}$ |  bootstrapping estimation for $\widehat{\boldsymbol{\theta}}_{k}$
Table 1: Summary of key variables and their descriptions

At any time $t$, for each single cell in the $i$-th metabolic phase, denoted
by $Z_{t}=i$ with $i\in\\{0,1,\ldots,I\\}$ (i.e., growth, stationary, death
phases, etc), the dynamic behaviors of its metabolic network can be
characterized by a stochastic model specified by parameters
$\boldsymbol{\alpha}_{i}$. The regulation mechanism and reaction rate of
metabolic network respond to the change of environmental condition is denoted
by $\boldsymbol{\nu}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})$. The
metabolic phase shift model with
$\mathbb{P}_{ij}(\boldsymbol{s}_{t};\boldsymbol{\beta})$ represents the
transition probability from $Z_{t}=i$ to $Z_{t}=j$. Therefore, in this paper,
we focus on calibrating the mechanistic model of the Bio-SoS specified by the
calibration parameters
$\boldsymbol{\theta}\equiv\\{\boldsymbol{\alpha}_{i},\boldsymbol{\beta}\\}$:
(1) $\boldsymbol{\alpha}_{i}$ characterizing cell dynamics in $i$-th metabolic
phase or class; and (2) $\boldsymbol{\beta}$ characterizing the phase shift
probability that impacts the percentage of cells in each phase or class.

Figure 1: An illustration of the multi-scale mechanistic model for cell
culture process and Bio-SoS.

###  2.1 Single Cell Model

This section briefly describes a stochastic metabolic model characterizing the
dynamics of a single cell across different metabolic phases: growth ($Z=0$),
production ($Z=1$) and so on. The model includes two main components: the
Stochastic Modular Reaction Networks (SMRN) and the metabolic phase shift
dynamics. Following the study (Anderson and Kurtz, 2011), the SMRN can be
constructed

$\tilde{\mathbf{s}}^{i}_{t+1}=\tilde{\mathbf{s}}^{i}_{t}+\mathbf{N}\cdot\boldsymbol{R}_{s}^{i}(t),$
(1)

where $\tilde{\mathbf{s}}^{i}_{t}$ is the state for cells in $i$-th phase,
$\mathbf{N}$ represents the stoichiometric matrix characterizing the structure
of molecular reaction network, and $\boldsymbol{R}_{s}^{i}(t)$ follows a
multivariate Poisson process with molecular reaction rate vector
$\boldsymbol{\nu}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})$. This reaction
rate depends on macro-environment $\boldsymbol{s}_{t}$, typically modeled by
using a Michaelis-Menten (MM) kinetic equation (Kyriakopoulos et al., 2018;
Wang et al., 2024).

The second component on metabolic shift model addresses the transitions
between metabolic phases in response to environmental changes. At any time
$t$, we have the metabolic phase-shift probability matrix
$\mathbb{P}(\boldsymbol{s}_{t};\boldsymbol{\beta})$ depending on the
environmental condition $\boldsymbol{s}_{t}$, where each element represents

$\mathbb{P}_{ii^{\prime}}(\boldsymbol{s}_{t})\equiv\mathbb{P}[Z_{t+1}=i^{\prime}|Z_{t}=i]~{}~{}~{}\mbox{for}~{}~{}~{}i,i^{\prime}=0,1,\ldots,I.$

###  2.2 Macro-Kinetic State Transition

Suppose the environmental condition $\boldsymbol{s}_{t}$ in bioreactor is
homogeneous over space. The cell population dynamics are characterized by the
evolution of cell densities in each metabolic phase $i$, i.e.

$\frac{\text{d}X_{i,t}}{\text{d}t}=\mu_{i}{(\boldsymbol{s}_{t})}X_{i,t}+\sum_{i^{\prime}\neq
i}\mathbb{P}_{i^{\prime}i,t}({\boldsymbol{s}_{t}})X_{i^{\prime},t},$ (2)

where $\mu_{i}{(\boldsymbol{s}_{t})}$ represents growth rate characterizing
the growth of the cell population and the second term counts the cells
metabolic shifting from $i^{\prime}$-th phase to $i$-th phase.

Let $\boldsymbol{s}_{t}$ denote metabolite concentrations (i.e., the number of
molecules per unit of volume) in the system at time $t$. Let
$\Delta\tilde{\mathbf{s}}_{t+1}^{(i,n)}=\tilde{\mathbf{s}}^{(i,n)}_{t+1}-\tilde{\mathbf{s}}^{(i,n)}_{t}$
represent the change in metabolic concentration for the $n$-th cell in phase
$i$. Then at any time $t$, given the density $X_{i,t}$ of cells located in
phase $i$, the overall change of metabolite concentration is the sum of the
contributions from individual cells,

$\boldsymbol{s}_{t+1}-\boldsymbol{s}_{t}=\sum_{i=0}^{I}\sum^{X_{i,t}}_{n=1}\Delta\tilde{\mathbf{s}}^{(i,n)}_{t+1}=\sum_{i=0}^{I}\mathbf{N}\boldsymbol{R}^{i}(t),$
(3)

where $\boldsymbol{R}^{i}(t)$ is a multivariate Poisson process with reaction
rate vector
$X_{i,t}\boldsymbol{v}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})$; as
detailed in Eq. (10) from Zheng et al. (2024). Without losing generality, in
the following sections, we consider a two-phase mechanistic model with
parameters
$\boldsymbol{\theta}=[\boldsymbol{\alpha}_{0},\boldsymbol{\alpha}_{1},\boldsymbol{\beta}]$.

###  2.3 State Transition Probability

In an infinitesimal time change $\mathrm{d}t$, the change of state
$\boldsymbol{s}_{t}$ becomes,

$\mathrm{d}\boldsymbol{s}_{t}=\boldsymbol{s}_{t+\text{d}t}-\boldsymbol{s}_{t}=\sum_{i=0}^{1}\mathbf{N}\text{d}\boldsymbol{R}^{i}(t)=\mathbf{N}\text{d}\boldsymbol{R}(t),$

where $\boldsymbol{R}(t)$ is a multivariate Poisson process with parameter
$\sum_{i=0}^{1}X_{i,t}\boldsymbol{\nu}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})$.
Then
${\mathbb{E}}(\text{d}\boldsymbol{R}\mid\boldsymbol{s}_{t})=\sum_{i=0}^{1}X_{i,t}\boldsymbol{\nu}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})$
and ${\rm
Var}(\text{d}\boldsymbol{R}\mid\boldsymbol{s}_{t})\approx\text{diag}\\{\sum_{i=0}^{1}X_{i,t}\boldsymbol{\nu}^{r}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\\}$,
with $r=1,2,\cdots,R$ representing the component indices of
$\boldsymbol{\nu}$. The approximation is made under the assumption that the
time interval is small, so the state doesn’t change too much and different
reactions can be treated independent. By forming as a Stochastic Differential
Equations (SDEs), we have

$\displaystyle\mathbf{N}\text{d}\boldsymbol{R}(t)\mid\boldsymbol{s}_{t}$
$\displaystyle=\mathbf{N}{\mathbb{E}}(\text{d}\boldsymbol{R}\mid\boldsymbol{s}_{t})\text{d}t+\mathbf{N}{\rm
Var}(\text{d}\boldsymbol{R}\mid\boldsymbol{s}_{t})\text{d}\boldsymbol{B}_{t}$
$\displaystyle\approx\mathbf{N}\sum_{i=0}^{1}X_{i,t}\boldsymbol{\nu}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\text{d}t+\left\\{\mathbf{N}\text{diag}\\{\sum_{i=0}^{1}X_{i,t}\boldsymbol{\nu}^{r}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\\}\mathbf{N}^{\top}\right\\}^{\frac{1}{2}}\text{d}\boldsymbol{B}_{t},$
(4)

where $\boldsymbol{B}_{t}$ denotes standard Brownian motion. It is hard to
find the analytic solution to (4), so we take an approximation in a short time
interval $(t,t+\Delta t]$, i.e.,

$\displaystyle\mathbf{N}\Delta\boldsymbol{R}(t)\mid\boldsymbol{s}_{t}$
$\displaystyle\approx\mathbf{N}\sum_{i=0}^{1}X_{i,t}\boldsymbol{\nu}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\Delta
t+\left\\{\mathbf{N}\text{diag}\\{\sum_{i=0}^{1}X_{i,t}\boldsymbol{\nu}^{r}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\\}\mathbf{N}^{\top}\right\\}^{\frac{1}{2}}\Delta\boldsymbol{B}_{t}$
$\displaystyle\sim\mathcal{N}\left\\{\mathbf{N}\sum_{i=0}^{1}X_{i,t}\boldsymbol{\nu}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\Delta
t,\mathbf{N}\text{diag}\\{\sum_{i=0}^{1}X_{i,t}\boldsymbol{\nu}^{r}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\\}\mathbf{N}^{\top}\text{diag}\\{\Delta
t\\}\right\\},$

where $\Delta\boldsymbol{B}_{t}$ is a Gaussian random vector with mean zero
and covariance matrix $\text{diag}\\{\Delta t\\}$. Finally, we approximate
conditional distribution of $\boldsymbol{s}_{t+1}$ as

$\displaystyle(\boldsymbol{s}_{t+1}-\boldsymbol{s}_{t})|\boldsymbol{s}_{t},\boldsymbol{\theta}\sim\mathcal{N}\\{\boldsymbol{\mu}(\boldsymbol{\theta}),\Sigma(\boldsymbol{\theta})\\},$
(5)

where
$\boldsymbol{\mu}(\boldsymbol{\theta})=\mathbf{N}\sum_{i=0}^{1}X_{i,t}\boldsymbol{\nu}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\Delta
t$,
$\Sigma(\boldsymbol{\theta})=\mathbf{N}\text{diag}\\{\sum_{i=0}^{1}X_{i,t}\boldsymbol{\nu}^{r}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\\}\mathbf{N}^{\top}\text{diag}\\{\Delta
t\\}\\}$.

##  3 Digital Twin Calibration for Bio-SoS Mechanistic Model

To facilitate the mechanistic model digital twin calibration, in this paper,
we consider the physical system as a finite horizon stochastic process. Its
dynamics can be characterized by a Bio-SoS mechanistic model specified by Eq.
(2) and (3) with underlying true parameters, denoted by
$\boldsymbol{\theta}^{c}=[\boldsymbol{\alpha}_{0}^{c},\boldsymbol{\alpha}_{1}^{c},\boldsymbol{\beta}^{c}]$.
Our goal is to develop an optimal learning calibration approach that can
efficiently guide the sequential design of experiments and most informative
data collection to improve model fidelity and the prediction accuracy of the
digital twin. For simplification, we consider batch-based cell culture
experiments with the selection of initial concentration $\boldsymbol{s}_{0}$
as the decision variable.

The trajectory
$\tau=(\boldsymbol{s}_{0},\boldsymbol{s}_{1},\ldots,\boldsymbol{s}_{T})$ over
a horizon of $T$ time steps has the joint probability density function

$p(\boldsymbol{\tau}|\boldsymbol{\theta},\boldsymbol{s}_{0})=\prod_{t=0}^{T-1}p(\boldsymbol{s}_{t+1}|\boldsymbol{s}_{t};\boldsymbol{\theta})\quad\text{
with $\boldsymbol{s}_{0}=\pi(\boldsymbol{\omega})$,}$ (6)

where $\pi(\boldsymbol{\omega})$ is a deterministic design or calibration
policy parameterized by $\boldsymbol{\omega}$. Each transition probability
$p(\boldsymbol{s}_{t+1}|\boldsymbol{s}_{t};\boldsymbol{\theta})$ follows a
normal distribution as defined in Eq. (5).

At each $k$-th calibration iteration, we select one design following the
latest policy, i.e. $\boldsymbol{s}_{0}^{(k)}=\pi(\boldsymbol{\omega}_{k})$,
run an experiment to collect one sample path from the physical system, and
then update the parameter $\widehat{\boldsymbol{\theta}}_{k}$. Specifically,
given the historical samples
$\mathcal{D}_{k}=\\{\boldsymbol{\tau}_{n}\\}_{n=1}^{k}$, we consider
maximizing the log-likelihood for model parameters $\boldsymbol{\theta}$,
i.e.,

$\widehat{\boldsymbol{\theta}}_{k}\equiv\mathop{\rm
arg\,max}_{\boldsymbol{\theta}}L(\mathcal{D}_{k};\boldsymbol{\theta})=\mathop{\rm
arg\,max}_{\boldsymbol{\theta}}\sum_{n=1}^{k}\ell\left(\boldsymbol{\tau}_{n}|\boldsymbol{\theta},\boldsymbol{s}_{0}^{(n)}\right),$
(7)

where $L(\mathcal{D}_{k},\boldsymbol{\theta})$ denotes the sum of log-
likelihoods and
$\ell(\tau_{n}|\boldsymbol{\theta},\boldsymbol{s}_{0}^{(n)})=\log
p(\boldsymbol{\tau}_{n}|\boldsymbol{\theta},\boldsymbol{s}_{0}^{(n)})$ is the
log-likelihood of trajectory $\tau_{n}$ with initial state
$\boldsymbol{s}_{0}^{(n)}=\pi(\boldsymbol{\omega}_{n})$. The detailed
derivation can be found in Section 4.

For digital twin calibration, we select the design of new experiment at the
beginning of $k$-th iteration. This selection is based on the evaluation of
our current model on a pre-specified set of evaluation and prediction points
$\mathcal{D}_{test}=\\{\boldsymbol{s}_{0}^{h}\\}_{h=1}^{H}$. The objective is
to minimize the MSE of digital twin prediction,

$\boldsymbol{\omega}^{\star}_{k}\triangleq\mathop{\rm
arg\,min}_{\boldsymbol{\omega}}J_{k}(\boldsymbol{\omega})\quad\mbox{with}\quad
J_{k}(\boldsymbol{\omega})\equiv\sum_{h=1}^{H}\mathbb{E}\left[Y^{\text{p}}\left(\boldsymbol{s}^{h}_{0},\boldsymbol{\theta}^{c}\right)-Y^{\text{d}}\left(\boldsymbol{s}^{h}_{0},{\boldsymbol{u}}_{k}(\widehat{\boldsymbol{\theta}}_{k-1},\boldsymbol{\omega})\right)\right]^{2},$
(8)

where
$\boldsymbol{u}_{k}(\widehat{\boldsymbol{\theta}}_{k-1},\boldsymbol{\omega})$
denote the virtual model update function with the simulated trajectory
$\boldsymbol{\tau}^{\text{d}}_{k}$ (see Eq. (10)) and
$Y=\phi(\boldsymbol{s}_{T})$ represents the system outputs of interest, such
as yield and product quality attributes. The superscript "p" indicates the
output from the physical system and "d" for digital twin.

To sample efficiently guide digital twin calibration, we propose a gradient-
based optimal learning approach with the procedure illustration as shown in
Figure 2 and develop a calibration algorithm summarized in Algorithm 1. In
specific, for the optimization problem (7) , we use the gradient-based method
to search for the MLE, i.e.,

$\widehat{\boldsymbol{\theta}}_{k}\leftarrow\widehat{\boldsymbol{\theta}}_{k-1}+\text{Adam}_{k}\left(\nabla_{\boldsymbol{\theta}}\ell(\boldsymbol{\tau}_{k};{\boldsymbol{\theta}},\boldsymbol{s}_{0}^{(k)})\mid_{\boldsymbol{\theta}=\widehat{\boldsymbol{\theta}}_{k-1}}\right),$
(9)

where “Adam” represents Adam optimizer (Kingma and Ba, 2014). The details of
the optimization steps will be discussed in section 4. Now we can define
virtual model update function

$\boldsymbol{u}_{k}(\widehat{\boldsymbol{\theta}}_{k-1},\boldsymbol{\omega})=\widehat{\boldsymbol{\theta}}_{k-1}+\text{Adam}_{k}\left(\nabla_{\boldsymbol{\theta}}\ell(\boldsymbol{\tau}_{k}^{\text{d}};\boldsymbol{\theta},\boldsymbol{\omega})\mid_{\boldsymbol{\theta}=\widehat{\boldsymbol{\theta}}_{k-1}}\right)$
(10)

with a simulated trajectory $\boldsymbol{\tau}^{\text{d}}_{k}\sim
p(\tau|\widehat{\boldsymbol{\theta}}_{k-1},\boldsymbol{s}_{0})$ and
$\boldsymbol{s}_{0}=\pi(\boldsymbol{\omega})$. Then to guide informative
experiment and improve digital twin prediction, we solve the optimization (8)
by using the stochastic gradient descent

$\boldsymbol{\omega}_{k}=\boldsymbol{\omega}_{k-1}-\gamma_{k}\nabla_{\boldsymbol{\omega}}J_{k}(\boldsymbol{\omega})\mid_{\boldsymbol{\omega}=\boldsymbol{\omega}_{k-1}}.$
(11)

Figure 2: The procedure illustration of the proposed digital twin calibration
approach.

Inspired by Figure 2, we compute the gradient with respect to
$\boldsymbol{\omega}$ using the chain rule. Specifically, we take the gradient
with respect to $\boldsymbol{u}_{k}$, and then take the gradient to
$\boldsymbol{u}_{k}$ with respect to $\boldsymbol{\omega}$, i.e.,

$\boldsymbol{\omega}_{k}=\boldsymbol{\omega}_{k-1}-\gamma_{k}\nabla_{\boldsymbol{u}_{k}}\sum_{h=1}^{H}\mathbb{E}\left[Y^{\text{p}}\left(\boldsymbol{s}^{h}_{0},\boldsymbol{\theta}^{c}\right)-Y^{\text{d}}\left(\boldsymbol{s}^{h}_{0},\boldsymbol{u}_{k}\right)\right]^{2}\Bigm{|}_{\boldsymbol{u}_{k}=(\widehat{\boldsymbol{\theta}}_{k-1},\boldsymbol{\omega}_{k-1})}\times\nabla_{\boldsymbol{\omega}}\boldsymbol{u}_{k}\left(\widehat{\boldsymbol{\theta}}_{k-1},\boldsymbol{\omega}\right)\Bigm{|}_{\boldsymbol{\omega}=\boldsymbol{\omega}_{k-1}}.$
(12)

Since $\boldsymbol{\theta}^{c}$ is unknown, we use the bootstrapping to
quantify the model parameter estimation uncertainty in Section 5. Following
the upper part of Figure 2, the prediction error propagates from
$\boldsymbol{\omega}$ to model parameter update $\boldsymbol{u_{k}}$, and then
to our objective $J_{k}(\boldsymbol{\omega})$. By taking the gradient of
$J_{k}(\boldsymbol{\omega})$, we can discern the influence of the experimental
design specified by $\boldsymbol{\omega}$ on the final prediction error, which
assists in selecting a more effective design. Since deriving an exact
expression for this gradient is complex, we develop a linear noise
approximation closed-form representation, which will be detailed in Section 5.
This approximation provides an efficient method for improving calibration
policy.

Input:initial policy parameter $\boldsymbol{\omega}_{0}$, initial model
parameter $\widehat{\boldsymbol{\theta}}_{0}$, step size $\lambda$ for
parameter update, initial step size $\gamma_{0}$ for design policy update

For $k=1:K$

  1. 1.

Compute
$\boldsymbol{u}_{k}(\widehat{\boldsymbol{\theta}}_{k-1},\boldsymbol{\omega})$
by Eq. (10)

* 2.

Update design $\boldsymbol{\omega}_{k}$ by Eq. (19) 3. Do experiments in
physical system with design $\pi{(\boldsymbol{\omega}_{k})}$ to get new data
$\boldsymbol{\tau}_{k}$,
$\mathcal{D}_{k}=\mathcal{D}_{k-1}\cup\\{\boldsymbol{\tau}_{k}\\}$ 4. Update
model parameters $\widehat{\boldsymbol{\theta}}_{k}$ by Eq. (13) End

Algorithm 1 Gradient-based Optimal Learning

##  4 Model Inference and Parameter Update

The gradient of the log-likelihood is calculated as

$\displaystyle\nabla_{\boldsymbol{\theta}}L(\mathcal{D}_{k};\boldsymbol{\theta})=\nabla_{\boldsymbol{\theta}}\left[\sum_{n=1}^{k}\log
p(\boldsymbol{\tau}_{k}|\boldsymbol{\theta},\boldsymbol{s}_{0}^{(n)})\right]=\sum_{n=1}^{k}\sum_{t=0}^{T-1}\left[\nabla_{\boldsymbol{\theta}}\log
p(\boldsymbol{s}^{(n)}_{t+1}|\boldsymbol{s}^{(n)}_{t},\boldsymbol{\theta})\right].$

Each term $\log
p(\boldsymbol{s}_{t+1}^{(n)}|\boldsymbol{s}_{t}^{(n)},\boldsymbol{\theta})$ in
this summation involves gradients that are detailed in the following lemma.

######  Lemma 1.

Conditioning on the current state $\boldsymbol{s}_{t}$ and model parameter
$\boldsymbol{\theta}$, the state change during in any small time interval
$(t,t+\Delta t]$ is
$\Delta\boldsymbol{s}_{t+1}=\boldsymbol{s}_{t+1}-\boldsymbol{s}_{t}$ following
a multi-variate normal distribution
$\mathcal{N}\\{\boldsymbol{\mu}(\boldsymbol{\theta}),\Sigma(\boldsymbol{\theta})\\}$.
Then the gradient of log-likelihood $\log
p(\boldsymbol{s}_{t+1}|\boldsymbol{s}_{t},\boldsymbol{\theta})$ with respect
to $\boldsymbol{\theta}$ is given by

$\displaystyle\nabla_{\boldsymbol{\theta}}\log
p(\boldsymbol{s}_{t+1}|\boldsymbol{s}_{t},\boldsymbol{\theta})=$
$\displaystyle-\frac{1}{2}\text{tr}\left(\Sigma(\boldsymbol{\theta})^{-1}\nabla_{\theta}\Sigma(\boldsymbol{\theta})\right)+\nabla_{\theta}{\boldsymbol{\mu}}(\boldsymbol{\theta})^{\top}\Sigma(\boldsymbol{\theta})^{-1}(\boldsymbol{s}_{t+1}-\boldsymbol{s}_{t}-\boldsymbol{\mu}(\boldsymbol{\theta}))$
$\displaystyle-\frac{1}{2}(\boldsymbol{s}_{t+1}-\boldsymbol{s}_{t}-\boldsymbol{\mu}(\boldsymbol{\theta}))^{\top}\left(-\Sigma(\boldsymbol{\theta})^{-1}{\nabla_{\boldsymbol{\theta}}\Sigma(\boldsymbol{\theta})}\Sigma(\boldsymbol{\theta})^{-1}\right)(\boldsymbol{s}_{t+1}-\boldsymbol{s}_{t}-\boldsymbol{\mu}(\boldsymbol{\theta})).$

###### Proof.

By the definition of multi-variant normal distribution, we have

$p(\boldsymbol{s}_{t+1}|\boldsymbol{s}_{t},\boldsymbol{\theta})=\frac{1}{(2\pi)^{M/2}|\Sigma(\boldsymbol{\theta})|^{1/2}}\exp\left(-\frac{1}{2}(\boldsymbol{s}_{t+1}-\boldsymbol{s}_{t}-\boldsymbol{\mu}(\boldsymbol{\theta}))^{\top}\Sigma(\boldsymbol{\theta})^{-1}(\boldsymbol{s}_{t+1}-\boldsymbol{s}_{t}-\boldsymbol{\mu}(\boldsymbol{\theta}))\right).$

Taking the logarithm on both side, we have

$\log
p(\boldsymbol{s}_{t+1}|\boldsymbol{s}_{t},\boldsymbol{\theta})=-\frac{1}{2}\log|\Sigma(\boldsymbol{\theta})|-\frac{1}{2}(\boldsymbol{s}_{t+1}-\boldsymbol{s}_{t}-\boldsymbol{\mu}(\boldsymbol{\theta}))^{\top}\Sigma(\boldsymbol{\theta})^{-1}(\boldsymbol{s}_{t+1}-\boldsymbol{s}_{t}-\boldsymbol{\mu}(\boldsymbol{\theta}))-\frac{M}{2}\log(2\pi).$

The gradient for the first term is
$\nabla_{\boldsymbol{\theta}}\log|\Sigma(\boldsymbol{\theta})|=\text{tr}\left(\Sigma(\boldsymbol{\theta})^{-1}\nabla_{\boldsymbol{\theta}}\Sigma(\boldsymbol{\theta})\right),$
where $\nabla_{\boldsymbol{\theta}}\Sigma(\boldsymbol{\theta})$ is a tensor of
$M\cdot M\cdot q$ with $q$ as the length of $\boldsymbol{\theta}$, and "tr"
represents the trace for first two dimensions. The gradient for the second
term is

$\displaystyle{\nabla_{\boldsymbol{\theta}}}(\boldsymbol{s}_{t+1}-\boldsymbol{s}_{t}-\boldsymbol{\mu}(\boldsymbol{\theta}))^{\top}\Sigma(\boldsymbol{\theta})^{-1}(\boldsymbol{s}_{t+1}-\boldsymbol{s}_{t}-\boldsymbol{\mu}(\boldsymbol{\theta}))$
$\displaystyle=$
$\displaystyle-2\nabla_{\boldsymbol{\theta}}{\boldsymbol{\mu}}(\boldsymbol{\theta})^{\top}\Sigma(\boldsymbol{\theta})^{-1}(\boldsymbol{s}_{t+1}-\boldsymbol{s}_{t}-\boldsymbol{\mu}(\boldsymbol{\theta}))+(\boldsymbol{s}_{t+1}-\boldsymbol{s}_{t}-\boldsymbol{\mu}(\boldsymbol{\theta}))^{\top}\left(-\Sigma(\boldsymbol{\theta})^{-1}{\nabla_{\boldsymbol{\theta}}\Sigma(\boldsymbol{\theta})}\Sigma(\boldsymbol{\theta})^{-1}\right)(\boldsymbol{s}_{t+1}-\boldsymbol{s}_{t}-\boldsymbol{\mu}(\boldsymbol{\theta})).$

∎

Given the calculated gradient, we can use the Adam optimizer Kingma and Ba
(2014) for updating $\widehat{\boldsymbol{\theta}}$. Following the algorithm
in their paper, for the $k$-th iteration, we have

$\widehat{\boldsymbol{\theta}}_{k}=\widehat{\boldsymbol{\theta}}_{k-1}+\lambda\cdot\frac{\left(1-\xi_{1}\right)\sum_{l=1}^{k}\xi_{1}^{k-l}G_{l}}{\left(1-\xi_{1}^{k}\right)\sqrt{\frac{\left(1-\xi_{2}\right)\sum_{l=1}^{k}\xi_{2}^{k-l}G_{l}^{2}}{\left(1-\xi_{2}^{k}\right)}+\epsilon}},$
(13)

where $\lambda$ is the learning rate,
$G_{l}=\nabla_{\boldsymbol{\theta}}\ell(\boldsymbol{\tau}_{l};\boldsymbol{\theta},\boldsymbol{s}_{0}^{(l)})\mid_{\boldsymbol{\theta}=\widehat{\boldsymbol{\theta}}_{l-1}}$
is the gradient of the likelihood function at $l$-th iteration, $g_{l}^{2}$ is
the element-wise square, $\xi_{1}=0.9,\xi_{2}=0.999$ are hyper-parameters
related to exponential decay rates, $\epsilon=10^{-8}$ is a small number to
avoid zero in the denominator.

##  5 Policy Gradient Estimation and Optimal Learning

In this section, we turn our focus on learning the calibration policy by
minimizing the prediction MSE in (8). To develop a surrogate model of this
MSE, we employ the linear noise approximation in Section 5.1 to approximate
the stochastic dynamics of the system and the solution of SDE-based
mechanistic model, reflecting its sensitivity to initial conditions and
parameter values. This approach allows us to describe the system state at any
given time as a normally distributed random variable, facilitating the
derivation of analytical expressions for the gradient estimator. To keep the
mechanism information, in Section 5.2, a first-order Euler’s Method is
utilized to derive the closed-form solutions of the SDEs based mechanisms
obtained from linear noise approximation.

###  5.1 Linear Noise Approximation on Bio-SoS Dynamics

We use linear noise approximation (LNA) (Fearnhead et al., 2014) to derive the
solution of SDEs as shown in Lemma 2 that allows us to analyze the estimation
error propagation from ${\boldsymbol{\theta}}$ to the output prediction and
develop a surrogate model of the MSE (8). Conditional on $\boldsymbol{s}_{0}$
and $\boldsymbol{\theta}$, we denote the state of systems at time $t$ as
$\boldsymbol{s}_{t}(\boldsymbol{\theta})$. Considering the SDEs characterizing
the Bio-SoS dynamics as shown in (4), we have

$\displaystyle\mathrm{d}\boldsymbol{s}_{t}(\boldsymbol{\theta})$
$\displaystyle=\boldsymbol{b}(\boldsymbol{s}_{t}(\boldsymbol{\theta}))\text{d}t+\mathbf{Q}(\boldsymbol{s}_{t}(\boldsymbol{\theta}))\text{d}\boldsymbol{B}_{t},$
(14)

where
$\boldsymbol{b}=\mathbf{N}\sum_{i=0}^{1}X_{i,t+\text{d}t}\boldsymbol{\nu}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})$
and
$\mathbf{Q}=\left\\{\mathbf{N}\text{diag}\\{\sum_{i=0}^{1}X_{i,t+\text{d}t}\boldsymbol{\nu}^{k}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\\}\mathbf{N}^{\top}\right\\}^{\frac{1}{2}}.$

######  Lemma 2.

(Linear Noise Approximation) Suppose $\boldsymbol{s}_{0}\sim
N\left(\boldsymbol{\mu}_{0}^{*},\boldsymbol{\Sigma}_{0}^{*}\right)$,
$\boldsymbol{\eta}_{0}=\boldsymbol{\mu}_{0}^{*}$ and
$\boldsymbol{\Psi}_{0}=\boldsymbol{\Sigma}_{0}^{*}$. Then the solution of SDE
(14) can be approximated as

$\boldsymbol{s}_{t}(\boldsymbol{\theta})\sim
N\left(\boldsymbol{\eta}_{t},\boldsymbol{\Psi}_{t}\right),$

where $\boldsymbol{\eta}_{t},\boldsymbol{\Psi}_{t}$ are the solution of
ordinary differential equations (ODEs)

$\frac{\text{d}\boldsymbol{\eta}}{\text{d}t}=\boldsymbol{b}(\boldsymbol{\eta}),\quad\frac{\mathrm{d}\boldsymbol{\Psi}}{\mathrm{d}t}=\boldsymbol{\Psi}\mathbf{F}^{\top}+\mathbf{F}\boldsymbol{\Psi}+\mathbf{QQ}^{\top},$
(15)

where we assume $\boldsymbol{s}_{t}=\boldsymbol{\eta}_{t}+\mathbf{M}_{t}$ with
$\boldsymbol{\eta}_{t}$ representing the deterministic path and
$\mathbf{M}_{t}$ representing stochastic path,
$\boldsymbol{\Psi}_{t}\triangleq\operatorname{Var}[\mathbf{M}_{t}],$
$\mathbf{F}_{t}$ is the $M\times M$ matrix with components

$\mathbf{F}_{t,ij}=\left.\frac{\partial\boldsymbol{b}_{i}}{\partial\boldsymbol{s}_{t,j}}\right|_{\boldsymbol{\eta}_{t}}\text{
and }\mathbf{Q}_{t}=\mathbf{Q}(\boldsymbol{\eta}_{t}).$

###### Proof.

Proof can be find in Fearnhead et al. (2014) with
$\mathbf{m}_{t}\triangleq\mathbb{E}[\mathbf{M}_{t}]=0$. ∎

Obtaining an analytical solution for the series of ODEs described in (15) is
challenging. Therefore, we use a numerical solution approach, as described in
Section 5.2. This method involves applying the solution iteratively at each
time step $t$, allowing us to approximate the system’s behavior over time. For
the physical system and digital twin, we denote

$\boldsymbol{s}_{t}(\boldsymbol{\theta}^{c})\sim
N\left(\boldsymbol{\eta}_{t},\boldsymbol{\Psi}_{t}\right),\quad\boldsymbol{s}_{t}({\boldsymbol{u}}_{k})\sim
N\left(\widehat{\boldsymbol{\eta}}_{t},\widehat{\boldsymbol{\Psi}}_{t}\right).$

Then we have
$\boldsymbol{s}_{T}(\boldsymbol{\theta}^{c})-\mathbf{s}_{T}({\boldsymbol{u}}_{k})\sim
N\left(\boldsymbol{\eta}_{T}-\widehat{\boldsymbol{\eta}}_{T},\boldsymbol{\Psi}_{T}+\widehat{\boldsymbol{\Psi}}_{T}\right).$
For any $m$-th component $Y$, we have a surrogate model of the MSE (8)

$\displaystyle{\mathbb{E}}\left[\left(Y^{\text{p}}(\boldsymbol{\theta}^{c})-Y^{\text{d}}({\boldsymbol{u}}_{k})\right)^{2}\right]$
$\displaystyle=\left[{\mathbb{E}}\left(Y^{\text{p}}(\boldsymbol{\theta}^{c})-Y^{\text{d}}({\boldsymbol{u}}_{k})\right)\right]^{2}+{\rm
Var}\left(Y^{\text{p}}(\boldsymbol{\theta}^{c})-Y^{\text{d}}({\boldsymbol{u}}_{k})\right)$
$\displaystyle=\left[\boldsymbol{\eta}_{T,m}-\widehat{\boldsymbol{\eta}}_{T,m}\right]^{2}+\boldsymbol{\Psi}_{T,mm}+\widehat{\boldsymbol{\Psi}}_{T,mm}.$
(16)

###  5.2 First-Order Euler’s Method for Solving LNA

To get the solution of $\boldsymbol{s}_{t}$, we solve the ODEs outlined in
(15) numerically by using first-order Euler’s method. For the term
$\boldsymbol{\eta}_{t}$, the update equation is given by

$\displaystyle\boldsymbol{\eta_{T}}$
$\displaystyle=\boldsymbol{\eta}_{T-1}+b(\boldsymbol{\eta_{T-1}})\Delta
t=\boldsymbol{\eta}_{0}+\sum_{t=0}^{T-1}b(\boldsymbol{\eta_{t}})\Delta
t=\boldsymbol{\eta}_{0}+\sum_{t=0}^{T-1}\sum_{i=0}^{1}X_{i,t}\mathbf{N}\boldsymbol{\nu}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\Delta
t\triangleq\boldsymbol{\eta}_{0}+\sum_{t=0}^{T-1}f_{t}{(\boldsymbol{\theta})}.$
(17)

For the term $\boldsymbol{\Psi}_{t}$, the update equation is given by

$\displaystyle\boldsymbol{\Psi}_{T}$
$\displaystyle=\boldsymbol{\Psi}_{T-1}+\left[\boldsymbol{\Psi}_{0}\mathbf{F}_{T-1}^{\top}+\mathbf{F}_{T-1}\boldsymbol{\Psi}_{0}+\mathbf{Q_{T-1}Q_{T-1}}^{\top}\right]\Delta
t+O((\Delta t)^{2})$
$\displaystyle=\boldsymbol{\Psi}_{0}+\sum_{t=0}^{T-1}\left[\boldsymbol{\Psi}_{0}\mathbf{F}_{t}^{\top}+\mathbf{F}_{t}\boldsymbol{\Psi}_{0}+\mathbf{Q_{t}Q_{t}}^{\top}\right]\Delta
t+O((\Delta t)^{2})$
$\displaystyle=\boldsymbol{\Psi}_{0}+\boldsymbol{\Psi}_{0}\left[\sum_{t=0}^{T-1}\sum_{i=0}^{1}\left[\mathbf{N}\boldsymbol{\nu}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\nabla_{\boldsymbol{s}_{t}}X_{i,t+\text{d}t}\right]\right]^{\top}\Delta
t+\left[\sum_{t=0}^{T-1}\sum_{i=0}^{1}\left[\mathbf{N}\boldsymbol{\nu}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\nabla_{\boldsymbol{s}_{t}}X_{i,t+\text{d}t}\right]\right]\boldsymbol{\Psi}_{0}\Delta
t$
$\displaystyle~{}~{}~{}~{}+\boldsymbol{\Psi}_{0}\left[\sum_{t=0}^{T-1}\sum_{i=0}^{1}\left[X_{i,t+\text{d}t}\mathbf{N}\nabla_{\boldsymbol{s}_{t}}\boldsymbol{\nu}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\right]\right]^{\top}\Delta
t+\left[\sum_{t=0}^{T-1}\sum_{i=0}^{1}\left[X_{i,t+\text{d}t}\mathbf{N}\nabla_{\boldsymbol{s}_{t}}\boldsymbol{\nu}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\right]\right]\boldsymbol{\Psi}_{0}\Delta
t$
$\displaystyle~{}~{}~{}~{}+\sum_{t=0}^{T-1}\left[\mathbf{N}\text{diag}\\{\sum_{i=0}^{1}X_{i,t+\text{d}t}\boldsymbol{\nu}^{r}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\\}\mathbf{N}^{\top}\right]\Delta
t+O((\Delta t)^{2})$
$\displaystyle\triangleq\boldsymbol{\Psi}_{0}+\sum_{t=0}^{T-1}g_{t}(\boldsymbol{\theta})+O((\Delta
t)^{2}).$ (18)

Then by (16)-(18), the final output discrepancy can be written as

$\displaystyle{\mathbb{E}}\left[\left(Y^{\text{p}}(\boldsymbol{\theta}^{c})-Y^{\text{d}}({\boldsymbol{u}}_{k})\right)^{2}\right]=\left(\sum_{t=0}^{T-1}[f_{t,m}(\boldsymbol{\theta}^{c})-f_{t,m}({\boldsymbol{u}}_{k})]\right)^{2}+2\Psi_{0,mm}+\sum_{t=0}^{T-1}[g_{t,mm}(\boldsymbol{\theta}^{c})+g_{t,mm}({\boldsymbol{u}}_{k})].$

Then conditional on $\widehat{\boldsymbol{\theta}}_{k-1}$, by applying a chain
rule, we can update the policy parameters by

$\displaystyle\boldsymbol{\omega}_{k}\approx$
$\displaystyle\boldsymbol{\omega}_{k-1}-\gamma_{k}\nabla_{{\boldsymbol{u}}_{k}}\sum_{h=1}^{H}\left[\left(\sum_{t=0}^{T-1}[f_{t,m}(\boldsymbol{\theta}^{c})-f_{t,m}({\boldsymbol{u}}_{k})]\right)^{2}+\sum_{t=0}^{T-1}[g_{t,mm}(\boldsymbol{\theta}^{c})+g_{t,mm}({\boldsymbol{u}}_{k})]\middle]\right|_{\boldsymbol{u}_{k}=(\widehat{\boldsymbol{\theta}}_{k-1},\boldsymbol{\omega}_{k-1})}$
$\displaystyle\times\lambda\nabla_{\boldsymbol{\omega}}\frac{\left(1-\xi_{1}\right)\sum_{t=1}^{k}\xi_{1}^{k-t}\tilde{G}_{t}}{\left(1-\xi_{1}^{k}\right)\sqrt{\frac{\left(1-\xi_{2}\right)\sum_{t=1}^{k}\xi_{2}^{k-t}\tilde{G}_{t}^{2}}{\left(1-\xi_{2}^{k}\right)}+\epsilon}}\Biggm{|}_{\boldsymbol{\omega}=\boldsymbol{\omega}_{k-1}},$
(19)

where
$\tilde{G}_{k}=\nabla_{\boldsymbol{\theta}}\ell(\boldsymbol{\tau}_{k}^{\text{d}};\boldsymbol{\theta},\boldsymbol{s}_{0}^{(k)})\mid_{\boldsymbol{\theta}=\widehat{\boldsymbol{\theta}}_{k-1}}$
and $\tilde{G}_{t}=G_{t},\text{for }t=1,\ldots,k-1$. Here $\tilde{G}_{k}$ is a
gradient based on simulation trajectory $\boldsymbol{\tau}_{k}^{\text{d}}$
since $\boldsymbol{\omega}_{k}$ is updated at the beginning of $k$-th
iteration and we can not get physical data $\boldsymbol{\tau}_{k}$.

To quantify the estimation error of underlying true model parameters
$\boldsymbol{\theta}^{c}$, we use the parametric bootstrapping. Based on MLE
and data $\boldsymbol{\tau}_{1},\ldots,\boldsymbol{\tau}_{k-1}$, we have
$\widehat{\boldsymbol{\theta}}_{k-1}$ as the estimator for
$\boldsymbol{\theta}^{c}$. Then we can generate $l$ bootstrap samples, denoted
as $\tilde{\boldsymbol{\theta}}^{1}_{k-1}$,
$\tilde{\boldsymbol{\theta}}^{2}_{k-1},\ldots,$
$\tilde{\boldsymbol{\theta}}^{l}_{k-1}$. Then we can take the following
estimation

$f_{t,m}(\boldsymbol{\theta}^{c})\approx\frac{1}{l}\sum_{n=1}^{l}f_{t,m}(\tilde{\boldsymbol{\theta}}_{k-1}^{n}),\quad
g_{t,mm}(\boldsymbol{\theta}^{c})\approx\frac{1}{l}\sum_{n=1}^{l}g_{t,mm}(\tilde{\boldsymbol{\theta}}_{k-1}^{n}).$

##  6 Empirical Study

To assess the performance of the proposed calibration approach, we employ a
simpler version of the Chinese Hamster Ovary (CHO) cell culture model as
described by Ghorbaniaghdam et al. (2014) as our SMRN model. We use this model
with a state shift model to simulate cell culture dynamics and compare our
gradient-based policy against a random policy for a DoE problem. We detail the
model and its validation in Section 6.1, and present the comparative results
in Section 6.2.

###  6.1 Cell Culture Model and Validation

CHO cells have become the most commonly used mammalian hosts for industrial
production of monoclonal antibodies (mAbs) and recombinant protein
therapeutics. This is due to their capacity to correctly fold, assemble, and
modify human proteins, ensuring biological activity and stability. The
robustness of CHO cells, coupled with their genetic and metabolic
malleability, makes them ideal candidates for genetic engineering and
optimization for specific production processes. The cell growth model is
described by a function characterizing different phases of cell culture

$\displaystyle v_{growth0}$
$\displaystyle=v_{maxg0}\frac{GLC}{K_{mGLC_{0}}+GLC}X_{0},\quad\frac{\mathrm{d}X_{0}}{\text{d}t}=v_{growth0}+\beta_{1}X_{1}-\beta_{0}X_{0},\quad\frac{\mathrm{d}X_{1}}{\text{d}t}=\beta_{0}X_{0}-\beta_{1}X_{1},$

where $v_{growth0}$ represents the growth rate that depends on glucose
concentration in the growth phase (phase 0) and where $\beta_{0}$ and
$\beta_{1}$ are phase transition rates, reflecting the dynamic interplay
between growth and production phases. A simplified SMRN focuses on three key
metabolisms: glucose consumption, lactate production, and monoclonal antibody
synthesis with chemical equations

$\displaystyle GLC$ $\displaystyle\rightarrow 2LAC,\quad 0.02GLC\rightarrow
X,\quad 0.01GLC\rightarrow mAb$

where GLC, LAC and mAb represent glucose, lactate and monoclonal antibodies.
The state of the system at any time $t$ is denoted by
$\boldsymbol{s}_{t}=\left[[GLC]_{t},[LAC]_{t},[mAb]_{t}\right]$ representing
the concentrations of glucose, lactate, and mAb, respectively. The flux rates
$\boldsymbol{v}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})=(-v_{HK},2v_{HK},v_{mAb})$,
where $-v_{HK}$ denotes the rate of glucose consumption, $2v_{HK}$ represents
the stoichiometric production of lactate , and $v_{mAb}$ quantifies the mAb
production rate. The kinetic equations of two phases $i=0,1$ are described as

$\displaystyle v_{HK_{i}}$
$\displaystyle=v_{maxHK_{i}}\frac{GLC}{K_{mGLC_{i}}+GLC}X_{i,t},\quad
v_{HK}=v_{HK_{0}}+v_{HK_{1}}$ $\displaystyle v_{mAb_{i}}$
$\displaystyle=v_{maxmAb_{i}}\frac{GLC}{K_{mGLC_{i}}+GLC}X_{i,t},\quad
v_{mAb}=v_{mAb_{0}}+v_{mAb_{1}}$

The metabolic reactions are formulated as a SDE system (14) to account for the
inherent randomness and fluctuations in biochemical reactions. The general
form of this SDE is given by

$\mathrm{d}\boldsymbol{s}_{t}(\boldsymbol{\theta})=\boldsymbol{b}(\boldsymbol{s}_{t}(\boldsymbol{\theta}))\text{d}t+\mathbf{Q}(\boldsymbol{s}_{t}(\boldsymbol{\theta}))\text{d}\boldsymbol{B}_{t},$

where the mean function
$\boldsymbol{b}=\mathbf{N}\sum_{i=0}^{1}X_{i,t+\text{d}t}\boldsymbol{\nu}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})$
is modeled by the stoichiometric matrix

$\mathbf{N}=\begin{bmatrix}-1&-0.02&-0.01\\\ 2&0&0\\\ 0&0&1\end{bmatrix}$

and reaction rates
$\boldsymbol{v}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})$. The system
variability is modeled by
$\mathbf{Q}=\left\\{\mathbf{N}\text{diag}\\{\sum_{i=0}^{1}X_{i,t+\text{d}t}\boldsymbol{\nu}^{r}(\boldsymbol{s}_{t};\boldsymbol{\alpha}_{i})\\}\mathbf{N}^{\top}\right\\}^{\frac{1}{2}}$.

We fit our model by minimizing
$\text{MSE}=\frac{1}{n}\sum_{j=1}^{n}{[Y^{s}_{j}-Y^{e}_{j}]^{2}}$ between
experiment data $Y^{e}$ from Ghorbaniaghdam et al. (2014) and simulation model
prediction $Y^{s}$. The results of our model validation are depicted in Figure
3. The figure presents experimental data as red dots and the simulation
results as blue lines, which illustrate both the mean and the 95% confidence
interval (CI) of the model predictions. The CI is calculated by the formula
$[\bar{x}-z{\sigma}/{\sqrt{n}},\bar{x}+z{\sigma}/{\sqrt{n}}]$, where $\bar{x}$
is the sample mean, $z$ is the z-score corresponding to the 95% confidence
level, which is approximately 1.96, $\sigma$ is the standard deviation of the
sample, $n$ is the number of observations in the sample. This result shows the
model’s effectiveness in capturing the dynamics of the CHO cell culture
system.

Figure 3: The cell culture mechanistic model was validated by using the
experimental data from the literature (Ghorbaniaghdam et al., 2014). Red dots
represent the experimental data from the literature and the simulation
predictions by our model is demonstrated in blue lines with mean and 95% CI
computed from $n=20$ replications.

In the empirical study assessing the performance of the proposed calibration
approach, we use the fitted parameters as the true model parameters of the
physical system, i.e.,

$\displaystyle\boldsymbol{\alpha}_{0}$
$\displaystyle=[v_{maxHK0},v_{maxmAb0},K_{mGLC0},v_{maxg0}]=[4.15\times
10^{-5},3.5\times 10^{-8},4.2,0.126]$ (20)
$\displaystyle\boldsymbol{\alpha}_{1}$
$\displaystyle=[v_{maxHK1},v_{maxmAb1},K_{mGLC1}]=[1.95\times
10^{-5},7.5\times 10^{-8},3.5]$ (21) $\displaystyle\boldsymbol{\beta}$
$\displaystyle=[\beta_{0},\beta_{1}]=[0.3,0.1].$ (22)

and assume the digital twin has the same model structure with unknown
parameters.

###  6.2 Digital Twin Calibration

With the true parameters defined in (20) to (22), our objective is to
calibrate the digital twin model by strategically designing CHO cell culture
experiments. The aim is to minimize the model prediction error by choosing the
initial glucose concentration as design decision variable. In each experiment,
we adjust the initial glucose concentration and then collect a new batch of
data to sequentially refine our model.

For the calibration experiments setting, we assume time course data of the
states is collected every 12h until 72h, so each batch includes 7-time
sequence data. We focus on calibrating parameters
$\boldsymbol{\alpha}_{0},\boldsymbol{\alpha}_{1}$. The design space for the
initial concentrations of glucose is set within the range $[18,38]$ mM, and
the initial concentration of lactate and mAb are $[LAC]_{0}=0.8$ mM,
$[mAb]_{0}=0$ mM, the initial cell densities in Phase 0 and 1 are set to be
$X_{0,0}=180\times 10^{6}$ cells/L, $X_{1,0}=10\times 10^{6}$ cells/L,
respectively. The initial parameters for the digital twin are estimated based
on the first batch of data. Since the model parameters are in different
scales, we set a learning rate $\lambda=[1\times 10^{-5},1\times
10^{-7},1,0.1,1\times 10^{-5},1\times 10^{-7},0.1]$ for updating the
calibrating parameters $(\boldsymbol{\alpha}_{0},\boldsymbol{\alpha}_{1})$ and
$\gamma_{k}=5$ for updating the design policy by applying Eq. (19). The fixed
set of test data $\mathcal{D}_{test}$ is composed of initial states
$\boldsymbol{s}_{0}^{1}=[28,0.8,0]$ and $\boldsymbol{s}_{0}^{2}=[38,0.8,0]$.

We conducted 15 macro-replications, each consisting of $K=200$ iterations
under varying initial conditions, with initial glucose levels randomly chosen
by applying the uniform distribution over the interval $[18,38]$. We compare
the proposed gradient-based calibration approach with a random design
approach, where the design is uniformly sampled from $[8,38]$. To ensure
consistency across the two approaches, we used the same random seeds for these
two approaches. The mean prediction performance and the 95% confidence
interval, calculated across macro-replications, represented by the Mean
Relative Errors
$\text{MRE}=\frac{1}{7}\sum_{j=1}^{7}{|\widehat{\theta}_{(j)}-\theta_{(j)}^{c}|}/{|\theta_{(j)}^{c}|}$
with script $(j)$ representing the $j$-th component, are presented in Figure 5
for mAb prediction and in Figure 5 for parameter estimation, respectively. For
both figures, blue lines represent the performance of our approach, while
yellow lines represent the random design performance.

Figure 4: Yield Prediction Error

Figure 5: Parameter Estimation Error

The results demonstrate that the proposed Bio-SoS model calibration approach
is significantly more sample-efficient than a random design. After 125
iterations, MRE for mAb production predictions using our calibration approach
is reduced to approximately 60%, in contrast to roughly 200% for predictions
made through random design. Furthermore, the MRE for parameters using our
approach decreases to around 30% after 125 iterations, whereas the MRE for the
random policy stays as high as 400%.

##  7 CONCLUSION

This study develops a robust digital twin calibration approach for a Bio-SoS
mechanistic model in the context of cell culture processes. By strategically
guiding more informative data collection through the proposed gradient-based
calibration approach, we significantly enhance the model fidelity and
predictive accuracy of the digital twin. The empirical validation using the
CHO cell culture model underscores the superiority of our approach over the
traditional random design, showcasing its potential applicability across
various biological systems. Overall, this study not only contributes to the
theoretical advancements in bioprocess digital twin development but also holds
promise for practical implementations that could improve the efficiency and
effectiveness of biomanufacturing in the pharmaceutical industry.

## References

  * Zheng et al. [2024] Hua Zheng, Sarah Harcum, Jinxiang Pei, and Wei Xie.  Stochastic biological system-of-systems modelling for ipsc culture.  _Communications Biology_ , 7, 01 2024.  doi:10.1038/s42003-023-05653-w.
  * Corlu et al. [2020] Canan G. Corlu, Alp Akcay, and Wei Xie.  Stochastic simulation under input uncertainty: A review.  _Operations Research Perspectives_ , 7:100162, 2020.  ISSN 2214-7160.  doi:https://doi.org/10.1016/j.orp.2020.100162.  URL https://www.sciencedirect.com/science/article/pii/S221471602030052X.
  * Kennedy and O’Hagan [2001] Marc C. Kennedy and Anthony O’Hagan.  Bayesian calibration of computer models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 63(3):425–464, 2001.  doi:https://doi.org/10.1111/1467-9868.00294.  URL https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.00294.
  * Tuo and Jeff Wu [2016] Rui Tuo and C. F. Jeff Wu.  A theoretical framework for calibration in computer models: Parametrization, estimation and convergence properties.  _SIAM/ASA Journal on Uncertainty Quantification_ , 4(1):767–795, 2016.  doi:10.1137/151005841.  URL https://doi.org/10.1137/151005841.
  * Anderson and Kurtz [2011] David F. Anderson and Thomas G. Kurtz.  _Continuous Time Markov Chain Models for Chemical Reaction Networks_ , pages 3–42.  Springer New York, New York, NY, 2011.  ISBN 978-1-4419-6766-4.  doi:10.1007/978-1-4419-6766-4_1.  URL https://doi.org/10.1007/978-1-4419-6766-4_1.
  * Kyriakopoulos et al. [2018] Sarantos Kyriakopoulos, Kok Siong Ang, Meiyappan Lakshmanan, Zhuangrong Huang, Seongkyu Yoon, Rudiyanto Gunawan, and Dong-Yup Lee.  Kinetic modeling of mammalian cell culture bioprocessing: The quest to advance biomanufacturing.  _Biotechnology Journal_ , 13(3):1700229, 2018.  doi:https://doi.org/10.1002/biot.201700229.  URL https://analyticalsciencejournals.onlinelibrary.wiley.com/doi/abs/10.1002/biot.201700229.
  * Wang et al. [2024] Keqi Wang, Wei Xie, and Sarah W. Harcum.  Metabolic regulatory network kinetic modeling with multiple isotopic tracers for ipscs.  _Biotechnology and Bioengineering_ , 121(4):1335–1353, 2024.  doi:https://doi.org/10.1002/bit.28609.  URL https://analyticalsciencejournals.onlinelibrary.wiley.com/doi/abs/10.1002/bit.28609.
  * Kingma and Ba [2014] Diederik Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _International Conference on Learning Representations_ , 2014.
  * Fearnhead et al. [2014] Paul Fearnhead, Vasilieos Giagos, and Chris Sherlock.  Inference for reaction networks using the linear noise approximation.  _Biometrics_ , 70(2):457–466, 2014.  doi:https://doi.org/10.1111/biom.12152.  URL https://onlinelibrary.wiley.com/doi/abs/10.1111/biom.12152.
  * Ghorbaniaghdam et al. [2014] Atefeh Ghorbaniaghdam, Jingkui Chen, Olivier Henry, and Mario Jolicoeur.  Analyzing clonal variation of monoclonal antibody-producing cho cell lines using an in silico metabolomic platform.  _PLOS ONE_ , 9, 2014.  URL https://api.semanticscholar.org/CorpusID:5913733.