# Learning Physics between Digital Twins with Low-Fidelity Models and Physics-
Informed Gaussian Processes

Michail Spitieris1,2, Ingelin Steinsland1
1Department of Mathematical Sciences
2Department of Public Health and Nursing
Norwegian University of Science and Technology (NTNU)
{michail.spitieris<EMAIL_ADDRESS>
###### Abstract

A digital twin is a computer model that represents an individual, for example
a component, a patient or a process. In many situations we want to gain
knowledge about an individual from its data while incorporating imperfect
physical knowledge and also learn from data from other individuals. In this
paper we introduce and demonstrate a fully Bayesian methodology for learning
between digital twins in a setting where the physical parameters of each
individual are of interest. For each individual the methodology is based on
Bayesian calibration with model discrepancy. Through the discrepancy, modeled
as a Gaussian process, the imperfect low-fidelity physical model is accounted
for. Using ideas from Bayesian hierarchical models, a joint probabilistic
model of digital twins is constructed by connecting them through a new level
in the hierarchy. For the physical parameters the methodology can be seen as
using a prior distribution in the individual model that is the posterior of
the corresponding hyperparameter in the joint model. For learning the
imperfect physics between individuals two approaches are introduced, one that
assumes the same discrepancy for all individuals and one that can be seen as
using a prior learned from all individuals for the parameters of the Gaussian
processes representing the discrepancies. Based on recent advances related to
physics informed priors, Hamiltonian Monte Carlo methods and using these for
inverse problems we set up an inference methodology that allows our approach
to be computational feasible also for physical models based on partial
differential equations and individual data that are not aligned. The
methodology is demonstrated in two synthetic case studies, a toy example
previously used in the literature extended to more individuals and an example
based on a cardiovascular differential equation model relevant for treatment
of hypertension. The case studies show that 1) models not accounting for
imperfect physical models are biased and over-confident, 2) the models
accounting for imperfect physical models are more uncertain, but covers the
truth, 3) the joint models have less uncertainty than the corresponding
independent individual models, but are not over-confident.

##  1 Introduction

A definition of a digital twin is a virtual representation of a physical asset
enabled through data and simulators [15]. Simulators refer to mathematical
models and physical model based on first principles are often preferred both
to incorporate system knowledge and to gain explainability. A known challenge
of deployment of digital twins is scalability, i.e. the ability to accessible
robust digital twin implementations at scale [15, 10].

The ability to make better decisions is the motivation for digital twins. In
this paper, we consider the case that better decisions can be made with
knowledge about the physical parameters of an imperfect physical model. In one
of our case studies, a low fidelity model for the cardiovascular system, the
Windkessel model [21] is used. This model is a differential equation linking
blood flow and blood pressure and has two physically interpretable parameters
that can direct the treatment of hypertension for patients. The parameters
will vary between patients, and we aim to estimate these using noisy flow and
pressure observations. It is known that low-fidelity models produce biased
parameter estimates if we don’t account for model discrepancy [5].

Recently, Spitieris and Steinsland [20] introduced and demonstrated a
methodology for Bayesian analysis of imperfect physical models using physics
informed Gaussian process priors. They combine a model formulation that
accounts for discrepancy introduced by Kennedy and O’Hagan [11] with Physics-
informed Gaussian process priors for physical models described by differential
equations [14]. This gives computational efficiency and uncertainty
quantification of the physical parameters of interest through a fully Bayesian
approach.

For each digital twin accounting for model discrepancy combined with noisy
data can result in too uncertain estimates to be of practical use. Further,
there can be identification issues related to the discrepancy. The working
hypothesis is that the uncertainty of the parameters for the individual
digital twins can be reduced by using information from other digital twins. In
this paper, we introduce a framework for how individual digital twins can
share data to reduce uncertainty. The physical models are cast into a Bayesian
hierarchical model framework which is extended with an extra level to allow
both the physical parameters and the discrepancy to gain information from the
other individuals. The code to replicate all the results in the paper is
available at https://github.com/MiSpitieris/Learning-Physics-between-DTs.

##  2 Background

###  2.1 Accounting for model discrepancy

Let $\mathbf{x}=(x_{1},\ldots,x_{k})$ be the $k$ input variables, and
$\boldsymbol{\phi}=(\phi_{1},\ldots,\phi_{p})$ the vector of unknown
parameters for the physical model $\eta(\bf{x},\boldsymbol{\phi})$. This model
is typically a simplified representation of the reality $\mathcal{R},$ and
does not fit the observed data well resulting in biased physical parameters
$\boldsymbol{\phi}$ estimates [5]. To account for this model form uncertainty,
Kennedy and O’Hagan [11] (KOH) suggested a Bayesian calibration framework
which incorporates a functional model discrepancy in the model formulation.
More specifically, the noise corrupted observed data
$\mathbf{y}^{\mathcal{R}}$ are described by the physical model $\eta$ and the
systematic model discrepancy $\delta$ as follows

$y^{R}(\mathbf{x})=\eta(\mathbf{x},\boldsymbol{\phi})+\delta(\mathbf{x})+\varepsilon,$
(1)

where $\varepsilon$ is the noise term. Since the discrepancy is an unknown
function of the inputs $\mathbf{x},$ KOH used a flexible Gaussian process
prior [16] to model the discrepancy, $\delta(\mathbf{x})\sim
GP(0,K_{\delta}(\mathbf{x},\mathbf{x}^{\prime}\mid\boldsymbol{\omega})),$
where $K_{\delta}$ is the covariance function and $\boldsymbol{\omega}$ are
kernel hyperparameters. The physical model is usually computational expensive
and thus KOH replaced the physical model numerical simulator with an emulator
which is a GP model trained on $N$ simulator output runs trained on an
experimental designed on the input and parameter space. This results to a GP
model which utilizes the two sources of information, $N$ simulator runs data
and $n$ observed data and has a computational complexity of
$\mathcal{O}((N+n)^{3}),$ where typically $N\gg n.$ If the deterministic model
is fast to evaluate, the computational complexity reduces. If an i.i.d.
Gaussian error term, $\varepsilon\sim N(0,\sigma^{2})$ is assumed we get the
following model formulation [8]

$y^{\mathcal{R}}\sim
GP(\eta(\mathbf{x},\boldsymbol{\phi}),K_{\delta}(\mathbf{x},\mathbf{x^{\prime}}\mid\boldsymbol{\omega})+\sigma^{2}),$
(2)

and this formulation will be used in Section 4. To infer model parameters
($\boldsymbol{\omega},\boldsymbol{\phi},\sigma$) in equation 2, Bayesian
inference is used by assigning priors to the model unknowns and specifically
priors that reflect underlying knowledge about $\boldsymbol{\phi},$ and the
posterior distribution is sampled using MCMC.

###  2.2 Physics-informed (PI) priors

We focus on the construction of physics-informed Gaussian process priors for
physical models that are described by linear differential equations
$\mathcal{L}_{x}^{\boldsymbol{\phi}}u(x)=f(x),$ where $\mathcal{L}$ is the
linear differential operator and $\boldsymbol{\phi}$ the vector of physical
parameters. To simplify the notation we take $x$ to be univariate, however the
following results apply to higher dimensions (e.g. see examples in [14, 20]).

Suppose we have observations of the functions $u$ and $f$ at potentially
different locations $x_{u1},\ldots,x_{un_{u}}$ and $x_{f1},\ldots,x_{fn_{f}},$
where $n_{u}$ and $n_{f}$ are the number of data for the functions $u$ and
$f$, respectively. The corresponding observed data are
$\mathbf{y}_{u}=(y_{u1},\ldots,y_{un_{u}})$ and
$\mathbf{y}_{f}=(y_{f1},\ldots,y_{fn_{f}}).$ The observed data are functional
and usually $u(x)$ is a smooth function. Thus it is often reasonable to assume
a GP prior, $u(x)\sim GP(0,K_{uu}(x,x^{\prime}\mid\boldsymbol{\theta}))$ to
describe the $u(x)$. The derivative of a GP is also a GP [1] and more
specifically $\textrm{Cov}\left(u(x),\frac{\partial u(x^{\prime})}{\partial
x^{\prime}}\right)=\frac{\partial
K_{uu}(x,x^{\prime}\mid\boldsymbol{\theta})}{\partial x^{\prime}}$ and
$\textrm{Cov}\left(\frac{\partial u(x)}{\partial x},\frac{\partial
u(x^{\prime})}{\partial
x^{\prime}}\right)=\frac{\partial^{2}K_{uu}(x,x^{\prime}\mid\boldsymbol{\theta})}{\partial
x\partial x^{\prime}},$ if the kernel is differentiable. Thus a convenient
choice (but not the only) is the squared exponential kernel since it is
infinitely differentiable.

Raissi et al. [14] utilized this result to build physics-informed priors for
$\mathcal{L}_{x}^{\boldsymbol{\phi}}u(x)=f(x)$ by assuming a GP prior on
$u(x).$ Then we have that $f(t)\sim
GP(0,K_{ff}(x,x^{\prime}\mid\boldsymbol{\theta},\boldsymbol{\phi})),$ where
$K_{ff}(x,x^{\prime}\mid\boldsymbol{\theta},\boldsymbol{\phi})=\mathcal{L}_{x}^{\boldsymbol{\phi}}\mathcal{L}_{x^{\prime}}^{\boldsymbol{\phi}}K_{uu}(x,x^{\prime}\mid\boldsymbol{\theta})$,
and also
$K_{uf}(x,x^{\prime}\mid\boldsymbol{\theta},\boldsymbol{\phi})=\mathcal{L}_{x^{\prime}}^{\boldsymbol{\phi}}K_{uu}(x,x^{\prime}\mid\boldsymbol{\theta})$
and
$K_{fu}(x,x^{\prime}\mid\boldsymbol{\theta},\boldsymbol{\phi})=\mathcal{L}_{x}^{\boldsymbol{\phi}}K_{uu}(x,x^{\prime}\mid\boldsymbol{\theta}).$
This results in a multi-output (of $u$ and $f$) GP model that satisfies the
differential equation. It is also convenient that the physical parameters
$\boldsymbol{\phi}$ have become hyperparameters of the physics-inspired
covariance matrix and therefore standard tools for inference can be used,
where Raissi et al. [14] used maximum likelihood and obtained point estimates.

This approach has been extended in a Bayesian framework [20] which also
accounts for model discrepancy in a model formulation similar to equation 1.
More specifically, for noise corrupted data $y_{u}$ and $y_{f}$, we assume a
discrepancy function on $u(x),$ and we have the following model formulation

$\begin{split}y_{u}&=u(x_{u})+\delta_{u}(x_{u})+\varepsilon_{u},\text{ where
}\delta_{u}(x)\sim GP(0,K_{\delta}(x,x^{\prime}\mid\boldsymbol{\omega}))\\\
y_{f}&=f(x_{f})+\varepsilon_{f},\end{split}$

where $\varepsilon_{u}\sim N(0,\sigma_{u}^{2})$ and $\varepsilon_{f}\sim
N(0,\sigma_{f}^{2})$ are the error terms. By assuming that $u(x)\sim
GP(\mu_{u}(x_{u}\mid\boldsymbol{\beta}),K_{uu}(x,x^{\prime}\mid\theta)),$
where $\mu_{u}$ is a mean function with parameters $\boldsymbol{\beta},$ and
consequently
$\mu_{f}(x_{f}\mid\boldsymbol{\beta},\boldsymbol{\phi})=\mathcal{L}_{x}^{\boldsymbol{\phi}}\mu_{u}(x_{u}\mid\boldsymbol{\beta})$
we get the following mutli-output GP that accounts for model discrepancy

$p(\mathbf{y}\mid\boldsymbol{\beta},\boldsymbol{\theta},\boldsymbol{\omega},\boldsymbol{\phi},\sigma_{u},\sigma_{f})=\mathcal{N}(\boldsymbol{\mu},\mathbf{K}_{\text{disc}}+\mathbf{S})$
(3)

where $\bf{y}=\begin{bmatrix}\bf{y}_{u}\\\ \bf{y}_{f}\end{bmatrix}$,
$\mathbf{K}_{\text{disc}}=\begin{bmatrix}K_{uu}(\mathbf{X}_{u},\mathbf{X}_{u}\mid\boldsymbol{\theta})+K_{\delta}(\mathbf{X}_{u},\mathbf{X}_{u}\mid\boldsymbol{\omega})&K_{uf}(\mathbf{X}_{u},\mathbf{X}_{f}\mid\boldsymbol{\theta},\boldsymbol{\phi})\\\
K_{fu}(\mathbf{X}_{f},\mathbf{X}_{u}\mid\boldsymbol{\theta},\boldsymbol{\phi})&K_{ff}(\mathbf{X}_{f},\mathbf{X}_{f}\mid\boldsymbol{\theta},\boldsymbol{\phi})\end{bmatrix}$,
$\mathbf{S}=\begin{bmatrix}\sigma_{u}^{2}I_{u}&0\\\
0&\sigma_{f}^{2}I_{f}\end{bmatrix}$ and
$\boldsymbol{\mu}=\begin{bmatrix}\boldsymbol{\mu}_{u}(\mathbf{X}_{u}\mid\boldsymbol{\beta})\\\
\boldsymbol{\mu}_{f}(\mathbf{X}_{f}\mid\boldsymbol{\beta,\phi})\end{bmatrix}.$
Finally, to infer model hyperparameters Hamiltonian Monte Carlo (HMC) [12] is
used and more specifically the NUTS variation [9].

For the rest of this paper we denote
$\mathbf{X}=(\mathbf{X}_{u},\mathbf{X}_{f}),$ where $\mathbf{X}$ can be
multivariate and eq. (3) can be written as
$p(\mathbf{y}\mid\boldsymbol{\beta},\boldsymbol{\theta},\boldsymbol{\omega},\boldsymbol{\phi},\sigma_{u},\sigma_{f})=\mathcal{N}(\boldsymbol{\mu}(\mathbf{X}\mid\boldsymbol{\beta,\phi}),\mathbf{K}_{\text{disc}}(\mathbf{X},\mathbf{X}\mid\boldsymbol{\theta},\boldsymbol{\omega},\boldsymbol{\phi},\sigma_{u},\sigma_{f})+\mathbf{S}).$

##  3 Hierarchical physical models accounting for model discrepancy

In this section, the model and inference method for sharing information
between individuals about parameters in physical models and discrepancy are
introduced. We assume that we have data from M individuals, where
$m=\\{1,\ldots,M\\}$ denotes the individual id, and $\mathbf{X}_{m}$ and
$\mathbf{y}_{m}$ the corresponding matrix of observed input data and vector
output for individual $m.$

Note that in the case of physics informed priors, the vectors $\mathbf{X}_{m}$
and $\mathbf{y}_{m}$ are joint for both the part related to $u$ and $f$ for
individual $m$ as described in Section 2.2. Our aim is to infer the physical
parameters for each individual, denoted $\boldsymbol{\phi}_{m}.$

###  3.1 Hierarchical physical model, shared global parameters

We now set up our most general hierarchical model where the individuals are
connected through distributions of the conditional priors of the individual
parameters. The model is illustrated in Figure 1 and mathematically formulated
below:

Figure 1: The hierarchical model detailed in Section 3 represented as directed
acyclic graph (DAG). It presents the general case where the personalized
physical model for the first and M-th individual is denoted as
$\boldsymbol{\eta}_{1}=\eta(\mathbf{X}_{1},\boldsymbol{\phi}_{1})$ and
$\boldsymbol{\eta}_{M}=\eta(\mathbf{X}_{M},\boldsymbol{\phi}_{M})$
respectively and we are interested in learning the individual physical
parameters $\boldsymbol{\phi}_{1},\ldots,\boldsymbol{\phi}_{M}.$ In the
individual latent field row, the red line represents the posterior mean of the
(low fidelity) latent physical model and the black dashed line the posterior
mean of the latent discrepancy model. At the bottom we have the posterior mean
of the bias corrected (low fidelity) model, where the black dots are the
observed data.

Individual likelihood: $\displaystyle\mathbf{y}_{m}$
$\displaystyle\mid\boldsymbol{\eta}_{m},\boldsymbol{\delta}_{m}\sim
P(\mathbf{y}_{m}\mid\boldsymbol{\eta}_{m},\boldsymbol{\delta}_{m})$ (4)
Individual latent field:
$\displaystyle[\boldsymbol{\eta}_{m},\boldsymbol{\delta}_{m}]$
$\displaystyle\mid\boldsymbol{\psi}_{m}\sim
GP(\boldsymbol{\mu}_{m},K_{m}\mid\boldsymbol{\zeta}_{m})$ (5) Individual
parameter priors: $\displaystyle\boldsymbol{\zeta}_{m}$
$\displaystyle\mid\boldsymbol{\psi}\sim
P(\boldsymbol{\zeta}_{m}\mid\boldsymbol{\psi})$ (6) Global parameters priors:
$\displaystyle\boldsymbol{\psi}$ $\displaystyle\sim P(\boldsymbol{\psi})$ (7)

We start by describing the model component for each individual $m$,
demonstrated for $m=1$ and $m=M$ in Figure 1. Given the individual latent
field $[\boldsymbol{\eta}_{m},\boldsymbol{\delta}_{m}]$ the observations are
assumed independent. There are possible parameters in the likelihood function,
and we will use independent Gaussian noise $N(0,\sigma^{2})$. The parameter
$\sigma^{2}$ is suppressed from the notation in this section due to
readability. The latent field consists of the physical model
$\boldsymbol{\eta}_{m}$ and the discrepancy $\boldsymbol{\delta_{m}}$.
Examples of these are given in Figure 1. We define
$\boldsymbol{g}_{m}=[\boldsymbol{\eta}_{m},\boldsymbol{\delta}_{m}]$. If we
have a model formulated as in Equation 2 we have $\mathbf{g}_{m}\sim
GP(\eta(\mathbf{X}_{m},\boldsymbol{\phi}_{m}),K_{\delta}(\mathbf{X}_{m},\mathbf{X}_{m}\mid\boldsymbol{\omega}_{m}))$
with $K_{\delta}$ and parameters as described in Section 2.1. If the PI prior
formulation in Equation (3) is used we have $\mathbf{g}_{m}\sim
GP(\boldsymbol{\mu}(\mathbf{X}_{m}\mid\boldsymbol{\beta}_{m},\boldsymbol{\phi}_{m}),\mathbf{K}_{\text{disc}}(\mathbf{X}_{m},\mathbf{X}_{m}\mid\boldsymbol{\theta}_{m},\boldsymbol{\phi}_{m},\boldsymbol{\omega}_{m})),$
with $\mathbf{K}_{\text{disc}}$ and parameters as described in Section 2.2.
The vector of all individual parameters for individual $m$ are denoted
$\boldsymbol{\zeta}_{m}$. If the model is formulated as in Equation 2 the
vector $\boldsymbol{\zeta}_{m}$ consists of the physical parameters
$\boldsymbol{\phi}_{m}$ and the discrepancy kernel hyperparameters
$\boldsymbol{\omega_{m}},$
$\boldsymbol{\zeta}_{m}=(\boldsymbol{\phi}_{m},\boldsymbol{\omega_{m})}$. When
the physics informed priors formulation in Equations 3 is used the vector
$\boldsymbol{\zeta}_{m}$ in addition includes the kernel and mean parameters
$(\boldsymbol{\theta}_{m},\boldsymbol{\beta}_{m}),$ and consequently
$\boldsymbol{\zeta}_{m}=(\boldsymbol{\phi}_{m},\boldsymbol{\theta}_{m},\boldsymbol{\beta}_{m},\boldsymbol{\omega}_{m}).$

The individual parameters are assumed to be conditionally independent given
the global parameters $\boldsymbol{\psi}$. Further, the global parameter is
given a prior distribution. Hence, all these $M$ individual models are
connected through the global parameters.

Due to the conditional independence in the hierarchical model, the joint
density of the model decomposes as follows

$P(\mathbf{y}_{1},\ldots,\mathbf{y}_{M},\mathbf{g}_{1},\ldots,\mathbf{g}_{M},\boldsymbol{\zeta}_{1},\ldots,\boldsymbol{\zeta}_{M},\boldsymbol{\psi})=P(\boldsymbol{\psi})\cdot\prod_{m=1}^{M}P(\mathbf{y}_{m}\mid\mathbf{g}_{m})\cdot
P(\mathbf{g}_{m}\mid\boldsymbol{\zeta}_{m})\cdot
P(\boldsymbol{\zeta}_{m}\mid\boldsymbol{\psi}),$

which means that the individual models are conditionally independent given the
global parameters $\boldsymbol{\psi}$.

###  3.2 Hierarchical physical model, common discrepancy and shared global
parameters

In some cases the discrepancy can be (almost) similar for all individuals, and
hence we want to assume common discrepancy parameters for all individuals. The
$m$ discrepancies are modeled as Gaussian processes with identical
hyperparameters. This corresponds to $\delta_{m}(\mathbf{X}_{m})\sim
GP(0,K(\mathbf{X}_{m},\mathbf{X}_{m}\mid\boldsymbol{\omega})),$ with
discrepancy hyperparameters $\boldsymbol{\omega}_{m}=\boldsymbol{\omega}$ that
are identical for all the individuals.

###  3.3 Inference

Suppose that we have observations from $M$ individuals where,
$\mathbf{y}_{\text{pop}}=(\mathbf{y}_{1},\ldots,\mathbf{y}_{M})$ and the
corresponding inputs are
$\mathbf{X}_{\text{pop}}=(\mathbf{X}_{1},\ldots,\mathbf{X}_{M}).$ The
posterior distribution of the unknown parameters
$\boldsymbol{\zeta}=(\boldsymbol{\zeta}_{1},\ldots,\boldsymbol{\zeta}_{M}),$
$\mathbf{g}=(\mathbf{g}_{1},\ldots,\mathbf{g}_{m}),$ $\boldsymbol{\psi}$ and
$\sigma$ is given by the following equation

$\displaystyle
P(\mathbf{g},\boldsymbol{\zeta},\boldsymbol{\psi}\mid\mathbf{y}_{\text{pop}},\mathbf{X}_{\text{pop}})=P(\boldsymbol{\psi})\cdot\prod_{m=1}^{M}P(\mathbf{y}_{m}\mid\mathbf{g}_{m},\mathbf{X}_{m})\cdot
P(\mathbf{g}_{m}\mid\boldsymbol{\zeta}_{m},\mathbf{X}_{m})\cdot
P(\boldsymbol{\zeta}_{m}\mid\boldsymbol{\psi}).$ (8)

This posterior distribution is analytically intractable, and we rely on
sampling methods. Since the dimension of the posterior is considerably large
traditional MCMC methods, for example a Metropolis within Gibbs implementation
can have slow mixing and fail to converge in practice. Hamiltonian Monte Carlo
[12, 3] provides an efficient alternative for sampling high dimensional spaces
(see for example [13]). The complex funnel-shape geometry of the posterior of
hierarchical models can be hard to sample, while non-center parametrization
can alleviate this problem [4] and it is used in this paper. More
specifically, we use the NUTS [9] variation of HMC, implemented in the
probabilistic language STAN [6].

If we assume for simplicity that each of the $M$ individual has $n$
observations, the computational complexity of the model is $\mathcal{O}(M\cdot
n^{3}).$ Where the cubic complexity of the Gaussian process is scalable, since
$n$ is typically relatively small, and the complexity increases just linearly
with the number of individuals $M.$

##  4 Toy example

We consider a conceptually simple model with one input parameter $x$ and one
physical parameter $u$ which has been used in the literature [2]. The model
represents the reality $\mathcal{R}$ and the misspecified model, $\eta$ are
the following exponential models

$\displaystyle y^{R}(x)$ $\displaystyle=3.5\cdot\exp(-u\cdot x)+b+\varepsilon$
$\displaystyle\eta(x,u)$ $\displaystyle=5\cdot\exp(-u\cdot x),$

where $\varepsilon$ is a Gaussian noise term. For noise free data and for a
given value of the physical parameter, $u=u_{0},$ the discrepancy between the
two models is $-1.5\cdot\exp(-u_{0}\cdot x)+b.$ However, we do not know the
functional form in practice and thus assume a zero mean Gaussian process prior
to describe the discrepancy function, $\delta(x)\sim
GP(0,K_{\delta}(x,x^{\prime}\mid\boldsymbol{\omega})).$

Figure 2: Toy example: Simulated noisy data for four of the ten individuals in
the data set. The blue line is the true model ($5\cdot\exp(-u_{m}\cdot
x_{m})+b_{m}$), and the dots are the corresponding noisy observed data
$y^{\mathcal{R}}$.

For the simulation study we assume there are $M=10$ individuals. For each
individual, data are simulated according to model $y^{R}(x)$. The true
physical parameter $u_{m}$ is set to be $0.7+0.1\cdot m$, and hence
$u_{1}=0.8$ and $u_{10}=1.7$. The individual offsets $b_{m}$ are sampled
randomly from a uniform distribution on the interval $[0.5,5].$ For all
individuals $\varepsilon\sim N(0,0.3^{2})$. In Figure 2 the true model and
simulated data are plotted for four individuals.

Figure 3: Toy example: $95\%$ credibility intervals for the physical
parameters based on posterior distributions for the four different models. The
dots represent the true values of the parameter $u_{m}$ for each of the
individual.

We fit four different models to the simulated data. The first model (no-
without delta in Figure 3), is the model $\eta$ with Gaussian noise without
assuming any model discrepancy, $\delta$ and therefore is the regression model
$y(x)=3.5\cdot\exp(-u\cdot x)+\varepsilon,$ where $\varepsilon\sim
N(0,\sigma^{2}).$ The second model (no-with delta in Figure 3) accounts for
model discrepancy and is given by equation 2. Both models do not share
individual information and are fitted for each of the $m$ participants
independently. The third model (yes-common delta in Figure 3) shares
information between individuals through a common discrepancy and inclusion of
a global level parameter as described in Section 3.2. The fourth model (yes-
shared delta in Figure 3) shares information between the individuals through
global parameters for both the discrepancy and the physical parameters as
described in Section 3.1. It allows the discrepancies to differ between
individuals, $\delta_{m}$ and also share information through the parameters.
Furthermore, for the models that account for discrepancy we assume that
$\delta_{m}(x_{m})\sim GP(0,K_{\delta}(x_{m},x_{m}^{\prime})),$ and we use the
squared exponential kernel
$K_{\delta}(x_{m},x_{m}^{\prime})=\alpha_{m}^{2}\exp\left(-\frac{(x_{m}-x_{m}^{\prime})^{2}}{2\rho_{m}^{2}}\right).$

In Figure 3, we see the 95% credible intervals (CIs) of the physical parameter
$u_{m}$ for the four different model fits for each individual. First, we
observe that if we do not account for model discrepancy (no-without delta) the
model produces biased and over-confident parameter estimates which comes in
line with Brynjarsdóttir and O’Hagan [5]. The posterior credible intervals
(CIs) of the independent model which accounts for model discrepancy (no-with
delta) cover the true values of the parameter. However, the uncertainties are
quite large and this can be impractical for decision purposes. The posterior
CIs of the models that share information (yes/common and shared delta), cover
the true parameter value and they drastically reduce the posterior
uncertainties. The model with individual discrepancies has the smallest
uncertainties. This can be explained by its flexibility that allows individual
discrepancies to share individual information between them without assuming a
common discrepancy parameters. More information about the models, priors,
predictions and implementation can be found in the Appendix.

##  5 Cardiovascular model example

###  5.1 Models

The Windkessel (WK) models [21] are linear differential equations that
describe the relationship between the blood pressure, $P(t)$ and blood flow,
$Q(t)$ in terms of physically interpretable parameters. The two parameters WK
model (WK2) is the basis for more complex models and is given by the following
time-dependent linear differential equation

$Q(t)=\frac{1}{R}P(t)+C\frac{dP(t)}{dt},$ (9)

where $R$ is the total vascular resistance and $C$ is the arterial compliance.
These hemodynamical parameters are the physical parameters of interest. The
three parameters WK model (WK3) is described by the following differential
equation,
$\frac{dP(t)}{dt}+\frac{P(t)}{R_{2}C}=\frac{Q(t)}{C}\left(1+\frac{R_{1}}{R_{2}}\right)+R_{1}\frac{dQ(t)}{dt}.$
The addition of the third parameter $R_{1}$ can increase model flexibility and
might fit the observed data better, though it overestimates the total arterial
compliance $C$ [19]. In practice, $R_{1}$ controls the amplitude of the blood
pressure waveform (see Figure 4) and therefore controls the discrepancy
between the two models. The WK2 model is a special case of the WK3 when
$R_{1}=0$ and also we have that
$R^{\text{WK2}}=R_{1}^{\text{WK3}}+R_{2}^{\text{WK3}}$ [21], and this is an
important connection for the simulation study in the following section.

Figure 4: Blood pressure (three cardiac cycles/heart beats) generated from WK2
model (red) and for a range of range of $R_{1}$ values $[0.02,0.12]$ from WK3
model (grey). The flow and $C$ values are identical for both models. The
amplitude of WK3-generated curve decreases linearly with $R_{1},$ while the
models become equivalent for $R_{1}=0.$

###  5.2 Simulation study

To validate our approach we need a model that is more complex than our
modelling choice, but also with known connections between their physical
parameters. Therefore, we use the WK2 model as a modelling choice, ($\eta$ in
Section 2.1) and simulate data from the WK3 model (see Figure 5), which we
consider as the true model ($\mathbf{y}^{\mathcal{R}}$ in Section 2.1). We
simulate data for $M=9$ individuals. For each individual $m=1,2,\ldots,9$ we
use an observed flow and given individual parameters $R_{1m},$ $R_{2m}$ and
$C_{m}$ to simulate pressure observations. The individual parameters are
chosen such that we have individuals with all the 9 possible combinations of
the values $R_{2}=(1,1.15,1.3)\text{ and }C=(0.95,1.1,1.25).$ The $R_{1m}$
parameter which controls the discrepancy between WK2 and WK3 is sampled
randomly from a uniform distribution on the interval $[0.02,0.1].$ More
specifically, for the blood flow $Q(t),$ we simulate the individual pressure
$P_{m}(t)=\text{WK3}(Q(t),R_{1m},R_{2m},C_{m}),$ and we use $\mathbf{t}_{P}$
temporal locations for the pressure observations and $\mathbf{t}_{Q}$ temporal
locations for the flow observations. The pressure observations and flow
observations are not required to be aligned. We add i.i.d Gaussian noise to
get flow and pressure observations as follows

$\displaystyle\mathbf{y}_{P_{m}}$
$\displaystyle=P_{m}(\mathbf{t}_{P})+\varepsilon_{P},\varepsilon_{P}\sim
N(0,4^{2})$ $\displaystyle\mathbf{y}_{Q_{m}}$
$\displaystyle=P_{m}(\mathbf{t}_{Q})+\varepsilon_{Q},\varepsilon_{Q}\sim
N(0,10^{2}).$

Figure 5: Cardiovascular simulation study: Simulated noisy pressure data for
four of the nine individuals in the data set. The red line is the true model
(WK3) and the dots are the observed blood pressure data.

Our modelling choice, the WK2 model (9) is a linear differential equation
which can be written as $\mathcal{L}_{t}^{\boldsymbol{\phi}}P(t)=Q(x),$ where
$\mathcal{L}_{t}^{\boldsymbol{\phi}}=\frac{1}{R}+C\frac{d}{dt}.$ By assuming
that $P(t)\sim GP(\mu_{P},K_{PP}(t,t^{\prime}\mid\boldsymbol{\theta})),$ where
$K_{PP}(t,t^{\prime})=\alpha_{\text{WK2}}^{2}\exp\left(-\frac{(t_{P}-t^{\prime}_{P})^{2}}{2\rho_{\text{WK2}}^{2}}\right),$
we construct a PI prior that accounts for model discrepancy as described in
Section 2.2, given by the formulation (3) and we fit four different models as
in the Section 4. The first model (no-without delta in Figure 6), is the PI
prior without the model discrepancy, which is equivalent to model (3) but
without the term $K_{\delta}$ in the first element of the covariance matrix.
The second model (no-with delta in Figure 6) accounts for model discrepancy
and is given by equation (3). Both models that do not share individual
information and are fitted for each participant $m$ independently. The third
model (yes-common delta in Figure 6) shares information between individuals by
assuming a common discrepancy and including global level parameters as
described in Section 3.2. The fourth model (yes-shared delta in Figure 6)
shares information between individuals through global parameters for both the
model and physical parameters. For the models that account for discrepancy we
assume that $\delta_{m}(t_{P})\sim
GP(0,K_{\delta_{m}}(t_{P},t^{\prime}_{P})),$ and we use the squared
exponential kernel,
$K_{\delta}(t_{P},t^{\prime}_{P})=\alpha^{2}\exp\left(-\frac{(t_{P}-t^{\prime}_{P})^{2}}{2\rho^{2}}\right).$

Figure 6: Cardiovascular simulation study: $95\%$ credibility intervals for
the physical parameters based on the posterior distributions for the four
different models for each individual in the . The dots represent the true
individual parameter values.

In Figure 6, we see the $95\%$ credible intervals for the posterior
distributions of the physical parameters $R_{m}$ and $C_{m},$ with the
corresponding true values represented by the dots. Note, that the true $R_{m}$
values are equal to $R_{1m}+R_{2m}.$ First, we observe that the model which
does not account for model discrepancy (no-without delta), systematically
underestimates both physical parameters $R$ and $C$ which comes in line with
[5] the results for the toy example in Section 4. For the $C$ parameter the
true value is within the $95\%$ CI, though it is in the upper tail of the
distribution. However, in practice the real observed data are different from
the WK3 simulated data and the fit without accounting for discrepancy might
result in larger biases. The independent models that account for model
discrepancy (no-with delta), cover the true values for all individuals, though
in some cases the uncertainty is quite large. For example, for the individuals
7-9 the posterior covers almost all the prior range, which is $R\sim
Unif(0.5,3),$ and thus the result can have low practical value. The models
that share information between individuals (yes/common and shared delta)
reduce the uncertainty substantially and still cover the true parameter
values. However, compared to the toy example in Section 4, the model that
shares discrepancy parameter information has similar performance to the model
which assumes the same discrepancy for all individuals. This can be understood
by observing that the discrepancies between the WK3 (true model) and the WK2
(modelling choice) are very similar, thus the correlation length scale and the
marginal variance ($\rho_{m}$ and $\alpha_{m}^{2}$ parameters in
$K_{\delta}$), are similar for all individuals. Consequently, there is no need
for considering different discrepancy parameters for each individual, and the
more parsimonious parametrization of the model is sufficient. More details
about the models, priors predictions and implementation can be found in the
Appendix.

##  6 Discussion

We have presented a fully probabilistic modelling framework for learning
physics between low fidelity physical models. To correct for model bias, we
have included a functional discrepancy in the model formulation. The physical
model and the model discrepancy have been incorporated in a Gaussian process
prior, and our aim is to learn the physical parameters and account for their
uncertainty. Therefore we introduced a Bayesian Hierarchical model where the
individual physical and discrepancy parameters are connected through global
parameters. We demonstrated our approach in two case studies, a toy example
from the literature and a cardiovascular model where the hemodynamical
parameters are of interest. We simulated data from a more complex model than
our modelling choice, and we show that we can learn the physical parameters of
the more complex model but also reduce the posterior uncertainty compared to
fitting the models independently for each individual.

The proposed modelling framework is scalable since the computational
complexity increases linearly with the number of individuals and is cubic for
each individual, where usually the number of individual data is typically
small. In cases where the individual data is large, GP methods for big data
can be used (see e.g. [7, 18]). A known problem is the identifiability between
the discrepancy and physical parameters in cases where prior information is
vague. A remedy in this case is to express shape constraints in the
discrepancy function through the derivative of the GP [5, 17].

## References

  * Adler [2010] R. J. Adler.  _The geometry of random fields_.  SIAM, 2010.
  * Bayarri et al. [2007] M. J. Bayarri, J. O. Berger, R. Paulo, J. Sacks, J. A. Cafeo, J. Cavendish, C.-H. Lin, and J. Tu.  A framework for validation of computer models.  _Technometrics_ , 49(2):138–154, 2007.
  * Betancourt [2017] M. Betancourt.  A conceptual introduction to Hamiltonian Monte Carlo.  _arXiv preprint arXiv:1701.02434_ , 2017.
  * Betancourt and Girolami [2015] M. Betancourt and M. Girolami.  Hamiltonian Monte Carlo for hierarchical models.  _Current trends in Bayesian methodology with applications_ , 79(30):2–4, 2015.
  * Brynjarsdóttir and O’Hagan [2014] J. Brynjarsdóttir and A. O’Hagan.  Learning about physical parameters: The importance of model discrepancy.  _Inverse problems_ , 30(11):114007, 2014.
  * Carpenter et al. [2017] B. Carpenter, A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Betancourt, M. Brubaker, J. Guo, P. Li, and A. Riddell.  Stan: A probabilistic programming language.  _Journal of statistical software_ , 76(1), 2017.
  * Hensman et al. [2015] J. Hensman, A. G. Matthews, M. Filippone, and Z. Ghahramani.  Mcmc for variationally sparse Gaussian processes.  _Advances in Neural Information Processing Systems_ , 28, 2015.
  * Higdon et al. [2004] D. Higdon, M. Kennedy, J. C. Cavendish, J. A. Cafeo, and R. D. Ryne.  Combining field data and computer simulations for calibration and prediction.  _SIAM Journal on Scientific Computing_ , 26(2):448–466, 2004.
  * Hoffman et al. [2014] M. D. Hoffman, A. Gelman, et al.  The no-u-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo.  _J. Mach. Learn. Res._ , 15(1):1593–1623, 2014\.
  * Kapteyn et al. [2021] M. G. Kapteyn, J. V. Pretorius, and K. E. Willcox.  A probabilistic graphical model foundation for enabling predictive digital twins at scale.  _Nature Computational Science_ , 1(5):337–347, 2021.
  * Kennedy and O’Hagan [2001] M. C. Kennedy and A. O’Hagan.  Bayesian calibration of computer models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 63(3):425–464, 2001.
  * Neal et al. [2011] R. M. Neal et al.  Mcmc using Hamiltonian dynamics.  _Handbook of markov chain monte carlo_ , 2(11):2, 2011.
  * Piironen and Vehtari [2017] J. Piironen and A. Vehtari.  Sparsity information and regularization in the horseshoe and other shrinkage priors.  _Electronic Journal of Statistics_ , 11(2):5018–5051, 2017.
  * Raissi et al. [2017] M. Raissi, P. Perdikaris, and G. E. Karniadakis.  Machine learning of linear differential equations using Gaussian processes.  _Journal of Computational Physics_ , 348:683–693, 2017\.
  * Rasheed et al. [2020] A. Rasheed, O. San, and T. Kvamsdal.  Digital twin: Values, challenges and enablers from a modeling perspective.  _Ieee Access_ , 8:21980–22012, 2020.
  * Rasmussen [2003] C. E. Rasmussen.  Gaussian processes in machine learning.  In _Summer School on Machine Learning_ , pages 63–71. Springer, 2003\.
  * Riihimäki and Vehtari [2010] J. Riihimäki and A. Vehtari.  Gaussian processes with monotonicity information.  In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , pages 645–652. JMLR Workshop and Conference Proceedings, 2010.
  * Rossi et al. [2021] S. Rossi, M. Heinonen, E. Bonilla, Z. Shen, and M. Filippone.  Sparse Gaussian processes revisited: Bayesian approaches to inducing-variable approximations.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1837–1845. PMLR, 2021.
  * Segers et al. [2008] P. Segers, E. Rietzschel, M. De Buyzere, N. Stergiopulos, N. Westerhof, L. Van Bortel, T. Gillebert, and P. Verdonck.  Three-and four-element windkessel models: assessment of their fitting performance in a large cohort of healthy middle-aged individuals.  _Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine_ , 222(4):417–428, 2008\.
  * Spitieris and Steinsland [2022] M. Spitieris and I. Steinsland.  Bayesian calibration of imperfect computer models using physics-informed priors.  _arXiv preprint arXiv:2201.06463_ , 2022.
  * Westerhof et al. [2009] N. Westerhof, J.-W. Lankhaar, and B. E. Westerhof.  The arterial windkessel.  _Medical & biological engineering & computing_, 47(2):131–141, 2009.

## Appendix

##  Appendix A Hierarchical physical model, common discrepancy and shared
global parameters (Details on Section 3.2)

This is the case that inputs are observed at the same domain and discrepancy
is expected to have similar characteristics. For example, in 1D and for the
squared exponential kernel,
$K_{\delta}(x_{m},x_{m}^{\prime})=\alpha_{m}^{2}\exp\left(-\frac{(x_{m}-x_{m}^{\prime})^{2}}{2\rho_{m}^{2}}\right)$
it means that the correlation decay characterized by $\rho_{m}$ and the
marginal variance $\alpha_{m}^{2}$ is similar for all individuals, thus
$\rho_{m}=\rho$ and $\alpha_{m}^{2}=\alpha^{2}$ for $m=1,\ldots,M.$ Therefore,
$\rho$ and $\alpha$ are not controlled by global parameters and their priors
distributions have fixed parameter values.

For simplicity, lets consider the case that the physical model $\eta$ is fast
to evaluate and the input and the physical parameters are univariate and we
have the following formulation
$y^{\mathcal{R}}(x)=\eta(x,\phi)+\delta(x)+\varepsilon$. We follow the
notation of Section 3.1 and for individual input data $\mathbf{X}_{m}$ we have
the latent field $\mathbf{g}_{m}\sim
GP(\eta(\mathbf{X}_{m},\boldsymbol{\phi}_{m}),K_{\delta}(\mathbf{X}_{m},\mathbf{X}_{m}\mid\boldsymbol{\omega}_{m})).$
In this case the vector of individual parameters $\boldsymbol{\zeta}_{m}$
consist of the individual physical parameter $\phi_{m}$ and the kernel
hyperparameters $\boldsymbol{\omega}_{m},$
$\boldsymbol{\zeta}_{m}=(\phi_{m},\boldsymbol{\omega}_{m}).$ In Section 3.1
both have priors $\phi_{m}\sim P(\phi_{m}\mid
a_{\phi},b_{\phi}),\boldsymbol{\omega}_{m}\sim
P(\boldsymbol{\omega}_{m}\mid\mathbf{a}_{\boldsymbol{\omega}},\mathbf{b}_{\boldsymbol{\omega}})$
that depend on the global parameters
$a_{\phi},b_{\phi},\mathbf{a}_{\boldsymbol{\omega}},\mathbf{b}_{\boldsymbol{\omega}},$
where $a_{\phi}\sim P(a_{\phi}),b_{\phi}\sim
P(b_{\phi}),\mathbf{a}_{\boldsymbol{\omega}}\sim
P(\mathbf{a}_{\boldsymbol{\omega}})\text{ and
}\mathbf{b}_{\boldsymbol{\omega}}\sim P(\mathbf{b}_{\boldsymbol{\omega}}).$
Now if we assume the same discrepancy across the individuals $m=1,\ldots,M,$
we have again that $\phi_{m}\sim P(\phi_{m}\mid a_{\phi},b_{\phi})$, but
$\mathbf{b}\sim P(\mathbf{b}).$ More specifically if we use the squared
exponential kernel for the discrepancy GP model, we have that
$\boldsymbol{\omega}=(\rho,\alpha)$ and the joint density from Section 3.1
decomposes as follows

$P(a_{\phi})P(b_{\phi})\cdot\prod_{m=1}^{M}P(\mathbf{y}_{m}\mid\mathbf{g}_{m})\cdot
P(\mathbf{g}_{m}\mid\phi_{m},\rho,\alpha)\cdot P(\phi_{m}\mid
a_{\phi},b_{\phi}),$

where $\rho\sim P(\rho),\alpha\sim P(\alpha).$

##  Appendix B Prior specification

In this Section we specify the non-center parameterization for the general
case where the model is fast to evaluate and in the case of the PI GP prior.

###  B.1 General formulation:
$y^{\mathcal{R}}(\mathbf{x})=\eta(\mathbf{x},\boldsymbol{\phi})+\delta(\mathbf{x})+\varepsilon$

For simplicity, lets assume that the input $x$ and the physical parameter
$\phi$ are both univariate. We take a Gaussian prior on the physical parameter
$\phi,$ $\phi\sim N(\mu_{\phi},\sigma^{2}_{\phi})),$ and $\mu_{\phi},$
$\sigma_{\phi}$ are the global level parameters, where $\mu_{\phi}\sim
P(\mu_{\phi})$ and $\sigma_{\phi}\sim P(\sigma_{\phi}).$ Then we use the
following non-center parameterization

$\begin{split}\tilde{\nu}_{\phi}&\sim N(0,1)\\\
\phi&=\mu_{\phi}+\sigma_{\phi}\cdot\tilde{\nu}_{\phi}\sim
N(\mu_{\phi},\sigma^{2}_{\phi}).\end{split}$ (10)

For the GP prior on the discrepancy $\delta(x)\sim
GP(0,K_{\delta}(x,x^{\prime}\mid\boldsymbol{\omega})),$ we use the squared-
exponential kernel,
$K_{\delta}(x,x^{\prime})=\alpha^{2}\exp\left(-\frac{(x-x^{\prime})^{2}}{2\rho^{2}}\right)$
and $\boldsymbol{\omega}=(\alpha,\rho).$ For both parameters we use
$\mathrm{Log\text{-}normal}$ priors. More specifically,
$\rho\sim\mathrm{Log\text{-}normal}(\mu_{\rho},\sigma^{2}_{\rho}),$ where the
distribution can be equivalently parameterized by the median
$\textrm{m}_{\rho},$
$\rho\sim\mathrm{Log\text{-}normal}(\log(\textrm{m}_{\rho}),\sigma^{2}_{\rho}).$
Then we use the following non-center parameterization

$\begin{split}\tilde{\nu}_{\rho}&\sim N(0,1)\\\
\rho&=\exp(\log(\textrm{m}_{\rho})+\sigma_{\rho}\cdot\tilde{\nu}_{\rho})\sim\mathrm{Log\text{-}normal}(\log(\textrm{m}_{\rho}),\sigma^{2}_{\rho})\end{split}$
(11)

and similarly for $\alpha$ we have that

$\begin{split}\tilde{\nu}_{\alpha}&\sim N(0,1)\\\
\alpha&=\exp(\log(\textrm{m}_{\alpha})+\sigma_{\alpha}\cdot\tilde{\nu}_{\alpha})\sim\mathrm{Log\text{-}normal}(\log(\textrm{m}_{\alpha}),\sigma^{2}_{\alpha}),\end{split}$
(12)

where $\textrm{m}_{\rho}\sim P(\textrm{m}_{\rho}),\sigma_{\rho}\sim
P(\sigma_{\rho}),\textrm{m}_{\alpha}\sim
P(\textrm{m}_{\alpha}),\sigma_{\alpha}\sim P(\sigma_{\alpha}).$

This can be generalized in cases where the physical parameters
$\boldsymbol{\phi}$ and inputs $\mathbf{x}$ have higher dimension.
Furthermore, the global level parameters the priors are chosen in a way that
reflects the population level properties.

###  B.2 Physics-informed Gaussian process

The main difference when we use the PI GP prior compared to B.1, is that it
involves the mean and kernel hyperparameters,
$(\boldsymbol{\beta},\boldsymbol{\theta})$ in addition to the physical
parameters $\boldsymbol{\phi}$. For the mean parameters, $\boldsymbol{\beta}$
we can use the same parameterization with the physical parameter $\phi$ as in
B.1, while for the $\boldsymbol{\theta}$ parameters we can use the same
parameterization we used for the discrepancy kernel in B.1.

##  Appendix C Prediction equations

We provide prediction equations for the two modelling cases discussed in the
paper. The first case is when the physical model
$\eta(\mathbf{x},\boldsymbol{\phi})$ is fast to run (e.g. a fast numerical
solver). The second case is when we can use the physics-informed Gaussian
process prior, and therefore no numerical discretization is needed.

###  C.1 General case

Standard GP predictions formulas can be used in this case, where the mean of
the GP prior is the output of the physical model
$\eta(\mathbf{x},\boldsymbol{\phi})$ (eq. (2) in Section 2.1). More
specifically, for the observed inputs $\mathbf{X},$ if $\mathbf{f}\sim
GP(\eta(\mathbf{X},\boldsymbol{\phi}),K(\mathbf{X,X^{\prime}})),$ at new
points $\mathbf{X_{*}}$ the joint distribution of the noise corrupted data
$\mathbf{y}=f(\mathbf{X})+\varepsilon,\varepsilon\sim N(0,\sigma^{2}I)$ and
$f(\mathbf{X_{*}})=\mathbf{f_{*}}$ is expressed as

$\begin{bmatrix}\mathbf{y}\\\
\mathbf{f}_{*}\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}\eta(\mathbf{X},\boldsymbol{\phi})\\\
\eta(\mathbf{X}_{*},\boldsymbol{\phi})\end{bmatrix},\begin{bmatrix}\mathbf{K}+\sigma^{2}I&\mathbf{K}_{*}\\\
\mathbf{K}_{*}^{T}&\mathbf{K}_{**}\end{bmatrix}\right),$ (13)

where $\mathbf{K}=K(\mathbf{X,X}),$ $\mathbf{K_{*}}=K(\mathbf{X,X_{*}})$ and
$\mathbf{K_{**}}=K(\mathbf{X_{*},X_{*}}).$ The conditional distribution o
$p(\mathbf{f}_{*}\mid\mathbf{X}_{*},\mathbf{X},\mathbf{y})$ is also
multivariate normal and more specifically

$\begin{split}p(\mathbf{f}_{*}\mid\mathbf{X}_{*},\mathbf{X},\mathbf{y})&=\mathcal{N}(\boldsymbol{\mu}_{*},\boldsymbol{\Sigma}_{*})\\\
\text{where
}\boldsymbol{\mu_{*}}&={\eta(\mathbf{X}_{*},\boldsymbol{\phi})}+\mathbf{K}_{*}^{T}(\mathbf{K}+\sigma^{2}I)^{-1}(\mathbf{y}-\eta(\mathbf{X},\boldsymbol{\phi}))\\\
\text{and
}\boldsymbol{\Sigma_{*}}&=\mathbf{K}_{**}-\mathbf{K}_{*}^{T}(\mathbf{K}+\sigma^{2}I)^{-1}\mathbf{K}_{*}\,.\end{split}$
(14)

###  C.2 Physics-informed Gaussian Process case

We present the prediction equations for the models formulated as PI GP priors
in Section 3.1 and 3.2. Suppose want to predict both functional outputs, $u$
and $f$ at new inputs $\mathbf{X}^{*}_{u}$ and $\mathbf{X}^{*}_{f}$
respectively. Let $\mathbf{u}^{*}=u(\mathbf{X}^{*}_{u})$ to be the predictions
for $u.$ The conditional distribution
$p(\mathbf{u}^{*}\mid\mathbf{X}^{*}_{u},\mathbf{X},\mathbf{y},\boldsymbol{\zeta})$
is multivariate Gaussian and more specifically

$\begin{split}p(\mathbf{u}_{*}\mid\mathbf{X}_{u}^{*},\mathbf{X},\mathbf{y},\boldsymbol{\zeta}_{\delta})&=\mathcal{N}(\boldsymbol{\mu}_{u}^{*},\boldsymbol{\Sigma}_{u}^{*})\\\
\boldsymbol{\mu_{u}^{*}}&=\mu_{u}\mathbf{(X_{u}^{*})}+\mathbf{V}_{u}^{*}{{}^{T}}(\mathbf{K_{\text{disc}}}+\mathbf{S})^{-1}(\mathbf{y}-\boldsymbol{\mu})\\\
\boldsymbol{\Sigma_{u}^{*}}&=K_{uu}(\mathbf{X}_{u}^{*},\mathbf{X}_{u}^{*})+K_{\delta}(\mathbf{X}_{u}^{*},\mathbf{X}_{u}^{*})-\mathbf{V}_{u}^{*}{{}^{T}}(\mathbf{K_{\text{disc}}}+\mathbf{S})^{-1}\mathbf{V}_{u}^{*},\end{split}$
(15)

where
$\mathbf{V}_{u}^{*}{{}^{T}}=\begin{bmatrix}K_{uu}(\mathbf{X}_{u}^{*},\mathbf{X}_{u})+K_{\delta}(\mathbf{X}_{u}^{*},\mathbf{X}_{u})&K_{uf}(\mathbf{X}_{u}^{*},\mathbf{X}_{f})\end{bmatrix}.$
The conditional distribution
$p(\mathbf{f}_{*}\mid\mathbf{X}_{f}^{*},\mathbf{X},\mathbf{y},\boldsymbol{\zeta})$
is multivariate Gaussian and more specifically

$\begin{split}p(\mathbf{f}_{*}\mid\mathbf{X}_{f}^{*},\mathbf{X},\mathbf{y},\boldsymbol{\zeta}_{\delta})&=\mathcal{N}(\boldsymbol{\mu}_{f}^{*},\boldsymbol{\Sigma}_{f}^{*})\\\
\boldsymbol{\mu_{f}^{*}}&=\mu_{f}\mathbf{(X_{f}^{*})}+\mathbf{V}_{f}^{*}{{}^{T}}(\mathbf{K_{\text{disc}}}+\mathbf{S})^{-1}(\mathbf{y}-\boldsymbol{\mu})\\\
\boldsymbol{\Sigma_{f}^{*}}&=K_{ff}(\mathbf{X}_{f}^{*},\mathbf{X}_{f}^{*})-\mathbf{V}_{f}^{*}{{}^{T}}(\mathbf{K_{\text{disc}}}+\mathbf{S})^{-1}\mathbf{V}_{f}^{*},\end{split}$
(16)

where
$\mathbf{V}_{f}^{*}{{}^{T}}=\begin{bmatrix}K_{fu}(\mathbf{X}_{f}^{*},\mathbf{X}_{u})&K_{ff}(\mathbf{X}_{f}^{*},\mathbf{X}_{f})\end{bmatrix}.$

##  Appendix D Toy example (Section 4)

In this Section additional results for the toy example in Section 4 are
presented.

###  D.1 Posteriors for individual parameters

In Figure 7, we plot the posterior distributions for the physical parameter of
interest u for all $M$ individuals (ID) and the four different approaches.
First, we observe that if we do not account for model discrepancy (no-without
delta), the posterior distributions for all individuals are biased and do not
cover the true value. The other three models which account for model
discrepancy cover the true value; however, the model that does not share
information (no-with delta) has too large uncertainty. In contrast, the other
two models share information on the physical parameter with common discrepancy
(yes/common delta) or shared information on the discrepancy (yes/ shared
delta), reducing the posterior uncertainty. The latter is more flexible,
allows for different discrepancies, and has the smallest uncertainty.

Figure 7: Toy example, posteriors.

###  D.2 Predictions

We now present the predictions obtained for all individuals as described in C.
In Figure 8, we plot the predictions for all four approaches. We have two
regions of predictions, on the left of the vertical dashed line where we have
observed data (interpolation) and on the right where we have not observed data
(extrapolation). First, we observe that the model which does not account for
model discrepancy (no-without delta) does not fit the observed data well, and
the prediction uncertainty is quite large. The other three approaches which
account for model discrepancy perform similarly in the regions where we have
observed data. While all three models perform well in regions without observed
data, the model that shares information about both the physical parameter $u$
and the discrepancy $\delta,$ has the smallest uncertainty.

Figure 8: Toy example, predictions.

##  Appendix E Cardiovascular model (Section 5)

###  E.1 Physics-informed GP model

We provide details about the construction of the physics-informed prior WK
model. The model is formulated as in Section 2.2, eq.(3) as follows

$\begin{split}y_{P}&=P^{\text{WK2}}(t_{P})+\delta(t_{P})+\varepsilon_{P}\\\
y_{Q}&=Q^{\text{WK2}}(t_{Q})+\varepsilon_{Q},\end{split}$ (17)

where $P^{\text{WK2}}(t_{P})\sim
GP(\mu_{P},K(t_{P},t_{P}^{\prime})),\varepsilon_{P}\sim N(0,\sigma_{P}^{2})$
and $\varepsilon_{Q}\sim N(0,\sigma_{Q}^{2})$. In addition we assume a GP
prior for the model discrepancy $\delta(t_{P}),$ $\delta(t_{P})\sim
GP(\mu_{P},K_{\delta}(t_{P},t_{P}^{\prime})),$ resulting in the following
multi-output GP prior

$p(\mathbf{y}\mid\boldsymbol{\theta},\boldsymbol{\phi},\sigma_{P},\sigma_{Q})=\mathcal{N}(\boldsymbol{\mu},\mathbf{K})$
(18)

where

$\begin{split}\bf{y}&=\begin{bmatrix}\bf{y}_{P}\\\ \bf{y}_{Q}\end{bmatrix},\\\
\boldsymbol{\mu}&=\begin{bmatrix}\boldsymbol{\mu}_{P}\\\
R^{-1}\cdot\boldsymbol{\mu}_{P}\end{bmatrix}\\\
\mathbf{K}&=\begin{bmatrix}K_{PP}(\mathbf{t}_{P},\mathbf{t}_{P}\mid\boldsymbol{\theta})+K_{\delta}(\mathbf{t}_{P},\mathbf{t}_{P}\mid\boldsymbol{\omega})+\sigma_{P}^{2}I_{P}&K_{PQ}(\mathbf{t}_{P},\mathbf{t}_{Q}\mid\boldsymbol{\theta},\boldsymbol{\phi})\\\
K_{QP}(\mathbf{t}_{Q},\mathbf{t}_{P}\mid\boldsymbol{\theta},\boldsymbol{\phi})&K_{QQ}(t\mathbf{t}_{Q},\mathbf{t}_{Q}\mid\boldsymbol{\theta},\boldsymbol{\phi})+\sigma_{Q}^{2}I_{Q}\end{bmatrix}\end{split}$
(19)

and

$\begin{split}K_{PQ}(t,t^{\prime})&=R^{-1}K_{PP}(t,t^{\prime})+C\frac{\partial
K_{PP}(t,t^{\prime})}{\partial t^{\prime}}\\\
K_{QP}(t,t^{\prime})&=R^{-1}K_{PP}(t,t^{\prime})+C\frac{\partial
K_{PP}(t,t^{\prime})}{\partial t}\\\
K_{QQ}(t,t^{\prime})&=R^{-2}K_{PP}(t,t^{\prime})+C^{2}\frac{\partial^{2}K_{PP}(t,t^{\prime})}{\partial
t\partial t^{\prime}}.\end{split}$ (20)

The PI GP model which does not account for model discrepancy is the same with
the difference that from the first element of the covariance matrix
$\mathbf{K},$ we remove the discrepancy kernel $K_{\delta}.$

###  E.2 Posteriors for individual parameters $R$ and $C$

In Figures 9 and 10, we plot the posterior distributions for $R$ and $C$ for
all four approaches. In both Figures, we observe that for the model without
discrepancy (no-without delta), the posteriors are over-confident and
underestimate both $R$ and $C.$ The model which accounts for model discrepancy
but does not share information (no-with delta) produces more reasonable
estimates of the physical parameters, though in some cases, it can be too
uncertain. For example, in Figure 10 the posterior can cover the whole range
of possible values. The models that share information (yes/common delta and
yes/shared delta) have reduced posterior uncertainty while covering the true
values of $R$ and $C.$ Furthermore, posterior densities for the two models are
very similar, suggesting that the more parsimonious model (yes/common
discrepancy) is sufficient.

Figure 9: Cardiovascular model, posterior distribution, R. Figure 10:
Cardiovascular model, posterior distribution, C.

###  E.3 Predictions

In Figures 11 and 12, the predictions for all four models for the two model
outputs (blood pressure and blood flow) are plotted. The model which does not
account for discrepancy produces more uncertain predictions than the other
three models. We see that this uncertainty increases with the values of the
parameter $R_{1},$ which controls the discrepancy between the two models (WK2
and WK3). Even if the posterior uncertainty for $R$ and $C$ is larger for the
model that does not share information (no-with delta), its prediction
uncertainty is similar to the two models that share information (yes/common
delta and yes/shared delta). In Figure 12, the blood flow predictions are
plotted. All models perform similarly, but the model that does not account for
discrepancy (no-without delta) has slightly larger prediction uncertainty.

Figure 11: Cardiovascular model, pressure predictions. Figure 12:
Cardiovascular model, flow predictions.