Sanjay Pant
email: sanjay.pant@inria.fr

Damiano Lombardi
email: damiano.lombardi@inria.fr

An information-theoretic approach to assess practical identifiability of parametric dynamical systems

Keywords: identifiability, entropy, mutual information, conditional mutual information, uncertainty quantification, non-parametric estimation, dynamical systems, computational mechanics

A new approach for assessing parameter identifiability of dynamical systems in a Bayesian setting is presented. The concept of Shannon entropy is employed to measure the inherent uncertainty in the parameters. The expected reduction in this uncertainty is seen as the amount of information one expects to gain about the parameters due to the availability of noisy measurements of the dynamical system. Such expected information gain is interpreted in terms of the variance of a hypothetical measurement device that can measure the parameters directly, and is related to practical identifiability of the parameters. If the individual parameters are unidentifiable, correlation between parameter combinations is assessed through conditional mutual information to determine which sets of parameters can be identified together. The information theoretic quantities of entropy and information are evaluated numerically through a combination of Monte Carlo and k-nearest neighbour methods in a non-parametric fashion. Unlike many methods to evaluate identifiability proposed in the literature, the proposed approach takes the measurement-noise into account and is not restricted to any particular noise-structure. Whilst computationally intensive for large dynamical systems, it is easily parallelisable and is non-intrusive as it does not necessitate re-writing of the numerical solvers of the dynamical system. The application of such an approach is presented for a variety of dynamical systemsranging from systems governed by ordinary differential equations to partial differential equations -and, where possible, validated against results previously published in the literature.

Introduction

This work is devoted to the setting up of a semi-empirical framework to assess identifiability in parameter estimation problems.

The identifiability problem was first stated by Bellman in [1] and it is a key issue in inverse problems (see [2]). It is still an open problem from a mathematical point of view, when a generic system described by a Partial Differential Equation (PDE) is investigated. The identifiability question can be summarised as follows: "given a system described by a set of equations and parametrised by certain scalar quantities and a set of measurements of the system, called observables, is it possible to determine the values of the parameters that account for the observables?"

For Ordinary Differential Equation (ODE) systems, several approaches are proposed in the literature to study identifiability (see [3] for a comprehensive review). Among them, the methods based on an algebraic approach, such as power expansion ( [4]), differential algebra (see for instance [5]) and Gröbner basis (see [6,7]) proved to be a valuable tool. The analyses based on differential algebra investigate the parameters-toobservable by assuming that the measurements are not affected by noise.

Other approaches, as the one proposed in the present work, are based on the construction of databases of pre-computed solutions or on the application of numerical methods in order to assess pragmatical identifiability (see [8]).

The proposed approach is based on a Bayesian framework. Its application in the context of inverse problems is commented in detail in [START_REF] Kaipio | Statistical and Computational inverse problems[END_REF]. Recently, it has also been exploited for model selection problems [START_REF] Oden | Mathematical Models and Methods in Applied Sciences[END_REF] and for experimental optimal design (see, for instance, [START_REF] Huan | [END_REF]12]). A work on the use of Fisher information to assess parameter identifiability is proposed in [13].

The present work aims at setting up a generic framework that could be useful for Computational Mechanics and Biophysical applications. A database of precomputed simulations is computed. In Uncertainty Quantification (UQ) [14][START_REF] Fishman | Monte Carlo: concepts, algorithms and applications[END_REF][START_REF] Xiu | [END_REF], the moments of the model outputs are computed when a structured uncertainty in the parameters is assumed. Instead, in the present work, the parameters-to-solution map is studied by means of entropies. These are estimated by using a Monte Carlo (MC) approach.

The use of entropies have been studied extensively in Information Theory since the pioneering work in [17]. The main goal is to set up a framework that allows to assess the identifiability from a pragmatical point of view, by taking into account the a priori knowledge about the parameters and the noise level in the measurements. The resulting approach is quite general and can be applied to a large class of systems. Although the approach can be computationally expensive (depending on the computational complexity of the dynamical system under consideration), it is easily parallelisable.

The structure of the work is as follows. In sections 2, 3, 4, and 5, the definition of the entropies and a general discussions about the key quantities involved is presented. Then, the numerical method, its implementation, and its scalability, are detailed in section 6. In the last part several numerical experiments of increasing complexity are presented. First, a Windkessel model, a commonly employed boundary condition in haemodynamics simulations, is investigated. Then, a non-linear epidemiological model is studied and the results are compared to those presented in [6]. The last two examples deal with PDE systems. The first example concerns the identification of potential in a system of harmonic waves, and the second example deals with parameter estimation in an Advection-Diffusion system.

Quantifying information gain through entropy

In the context of parameter estimation of dynamical systems a reasonable question to ask is the following: "How much information does one expect to gain about the parameters of interest through noisy measurements of an observable (or a set of observables) at a given set of discrete time instants?" A prerequisite to answer this question is to establish what is meant by information about a parameter. In a Bayesian context one may treat all the parameters as random variables and a priori, i.e. before any measurements are taken, express the beliefs about these parameters in the form of a prior probability distribution. Following the acquisition of measurements, one typically updates such beliefs using Bayes' rule to obtain the posterior probability distribution. This posterior reflects updated beliefs about the parameters based on both prior information and the measurements. If one can quantify the uncertainty in both the prior and posterior probability distributions then the difference between them can provide a reasonable quantification of the amount of information gained by the measurements. A precursory requirement, hence, is to quantify the uncertainty of a random variable with a given probability density function.

Shannon [17] introduced a measure of uncertainty or 'missing information' for discrete probability distributions that is widely used in the field of information theory. His measure of uncertainty, referred as Shannon entropy, is based on three intuitive notions about uncertainty: continuity, monotonic increase with increasing uncertainty, and the composition law (see [17] for details). For a discrete set of probabilities P 1 , P 2 , . . . , P n Shannon entropy, denoted by H, is defined as

S = n i P i log 1 P i . (1)
The units of Shannon entropy depend on the base of the logarithm used; for base 2, entropy is measured in bits, while if natural base is used then entropy is measured in nats. For a continuous random variable X with a probability density function p X (x), the analogue of Shannon entropy, often referred as differential entropy is defined as

H S (X) = H S (p X (x)) = X p X (x) log 1 p X (x) dx. ( 2
)
where X is the support of p X (x). The continuous version does not have the intuitive and desirable properties of the discrete version. The two notable drawbacks associated with the notion of uncertainty as defined by equation (2) are: first, that if x is not dimensionless then p X (x) has dimensions and taking log of a dimensional quantity presents problems; and second, that the uncertainty measure is not invariant under a change of variables. Jaynes [START_REF] Jaynes | Information Theory and Statistical Mechanics[END_REF] showed that a renormalised continuous limit of equation ( 1) that overcomes the drawbacks of equation ( 2) is

H J (X) = H J (p X (x)) = X p X (x) log m(x) p X (x) dx, (3)
where m(x) is an 'invariant measure' function introduced by Jaynes. Kullback and Leibler [19] (see also [START_REF] Kullback | Information theory and statistics[END_REF]) proposed that a more fundamental measure of information, referred as the discrimination information [START_REF] Ghosh | [END_REF] or generally the Kullback-Leibler (KL) distance/divergence, between two probability distributions p X (x) and q X (x) is

D (p X (x)||q X (x)) = X p X (x) log p X (x) q X (x) dx. (4)
The above expression represents the mean information, i.e. averaged over all outcomes x from p X (x), for discriminating between p X (x) and q X (x). Similar to the expression by Jaynes, equation (3), the expression by Kullback and Leibler, equation ( 4), is invariant under transformations of the random variable under consideration. The KL-distance can also be viewed as the loss of information if q X (x) is used to approximate p X (x). Alternatively, it can be viewed as the information gained about the random variable when its probability distribution changes from q X (x) to p X (x) [22]. Hobson [23] proved that the KL-distance of equation ( 4) satisfies all the intuitively reasonable properties desired in an information measure including those originally proposed by Shannon in [17]. Hobson and Cheng [22] argued that the uncertainty (missing information) in a probability distribution can be obtained by the KL-distance as follows: if p m X (x) denotes the distribution with maximum information content (subject to the constraints of the random variable X), and p o X (x) denotes the distribution with minimum information content, then the missing information in a probability distribution p X (x) can be seen as the difference between the information gain when the probability distribution changes from p o X (x) to p m X (x) and the information gain when the probability distribution changes from

p o X (x) to p X (x) H H (X) = H H (p X (x)) = D(p m X (x)||p o X (x)) -D(p X (x)||p o X (x)) (5) = X p m X (x) log p m X (x) p o X (x) dx - X p X (x) log p X (x) p o X (x) dx. (6)
They also showed that Shannon's measure of equation ( 2), Jaynes' measure of equation (3), and the measure of equation ( 5) based on the KL-distance, are identical (except for different scales), when the prior probability distribution p o X (x) is uniform and viewed as the invariant measure function m(x) of Jaynes.

While the quantities H S , H J , and H H can be used to quantify the inherent uncertainty of a random variable, the quantity of interest in this article is the gain in information (reduction in uncertainty). Consider a random variable X and a related observable Y (note that Y , too, is a random variable). Prior to any measurement of Y , let the probability distribution of X be p X (x), the prior probability distribution. Through a measurement Y = y one can employ Bayes' theorem to update the distribution of X from p X (x) to p X|Y (x|y), the posterior distribution. Note that this posterior is a distribution of X, and while writing it as p X|Y =y (x) or p X (x|Y = y) is more informative, it is denoted as p X|Y (x|y) for notational simplicity. Through Bayes' theorem one may write p X|Y (x|y) as

p X|Y (x|y) = likelihood p Y |X (y|x) prior p X (x) p Y (y) evidence . (7)
Consequently, the gain in information about X through the measurement Y = y can be written as

G X|Y =y (y) = H (p X (x)) -H p X|Y (x|y) H ∈ {H S , H J , H H } (8)
where G X|Y =y (y) denotes the aforementioned gain in information (measured through corresponding H: H S , H J , or H H ) and is a function of the measurement y. For an a priori analysis, i.e. without any real measurement of Y available, the average/expected gain in information is considered. This expected gain in information, I X|Y , can be obtained by integrating G X|Y =y (y) over all values of Y

I X|Y = Y G X|Y =y (y) p Y (y) dy, ( 9
)
where Y is the support of Y . Although different measures H (H S , H J , H H ) lead to different quantification of the inherent uncertainty in a random variable, all these measures, including the measure of KL-distance, lead to identical quantification of expected gain in information when the probability distribution of X changes from p X (x) to p X|Y (x|y). Proof of this proposition is presented in appendix Appendix A. Since the continuous form of Shannon entropy, the differential entropy of equation ( 2), has the simplest form, it is chosen to measure the gain in information. In what follows, the symbol H always denotes this continuous form of Shannon entropy. Hence, H(X) denotes the entropy of the marginal/prior probability distribution of X and is calculated as

H(X) = H (p X (x)) = X p X (x) log 1 p X (x) dx. (10)
H(X|Y ) = Y H (X|Y = y) p Y (y) dy = Y X p X|Y (x|y) log 1 p X|Y (x|y) p Y (y) dx dy H(X|Y ) = Y X p X,Y (x, y) log 1 p X|Y (x|y) dx dy, (11)
where p X,Y (x, y) denotes the joint distribution of the random variables X and Y . By equations ( 9), [START_REF] Oden | Mathematical Models and Methods in Applied Sciences[END_REF], and [START_REF] Huan | [END_REF], the expected reduction in uncertainty, ∆H(X|Y ), or equivalently the gain in information, I X|Y , about X by observing Y can be written as

∆H(X|Y ) = I X|Y = H(X) -H(X|Y ) (12)
Remark 2.1. The quantity ∆H(X|Y ) or I X|Y is also known as mutual information.

The mutual information between two random variables X and Y is a familiar concept in information theory and is interpreted as the average amount of information X carries about Y , and vice versa. In this manuscript it is denoted as M (X; Y ). It is usually defined in the following form, which is easy to derive from equations ( 10), [START_REF] Huan | [END_REF], and ( 12)

M (X; Y ) = Y X p X,Y (x, y) log p X,Y (x, y) p X (x) p Y (y) dx dy. (13)
Hence, the following three quantities are equivalent

∆H(X|Y ) = I X|Y = M (X; Y ) ( 14
)
Remark 2.2. Mutual information is also the KL-distance, see equation ( 4), between the joint distribution p X,Y (x, y) and the joint distribution of X and Y if they were independent. In the latter case the joint distribution factorises into the product of the marginal distributions, i.e. p X,Y (x, y) = p X (x) p Y (y) if X and Y are independent.

Remark 2.3. Mutual information is non-negative, M (X; Y ) ≥ 0 (see [START_REF] Cover | Elements of information theory[END_REF] for proof).

Remark 2.4. Mutual information is invariant under homeomorphic transformations of X and Y . If X = P (X) and Y = Q(Y ), where P and Q are smooth and uniquely invertible maps then M (X; Y ) = M (X ; Y ) (see [START_REF] Kraskov | [END_REF] for proof).

Remark 2.5. While the above theory is presented for scalar random variables X and Y , it also extends to vector random variables. For example, the continuous shannon entropy of a vector random variable X with probability distribution function p X (x) is

H(X) = -p X (x) log (p X (x)) dx.
Given the above background, the concept of expected entropy decrease (expected information gain) is applied to parameter estimation of dynamical systems in the next section.

Parameter Information gain in dynamical systems through observations

Consider a dynamical system of the following form

ẋ(t) = F (x(t), t, ϑ) , (15)
where x ∈ R Nx is the state of the system, F ∈ R Nx is a non-linear Lipschitz vector function, and ϑ ∈ R N ϑ are the parameters. Equation ( 15) can be considered, for instance, as a space discretization of a non-linear PDE. An observable for the system is a vector z ∈ R Nz such that:

z(t) = G(x(t)) + n(t), ( 16
)
where G is a non-linear function, often known as the observation operator, and n(t) ∈ R Nz is the noise component. The goal of an inverse problem for this system is to determine the parameter vector ϑ by measurements of z(t).

In equation ( 16) the noise is considered to be additive and i.i.d at all times. This, however, as will be presented later, is not a requirement of the presented method. The noise can also present itself in the forward model of equation [START_REF] Fishman | Monte Carlo: concepts, algorithms and applications[END_REF]. Furthermore, the presented method is not restricted to Gaussian nature of noise; any statistical process of noise, as long as one can draw independent samples from it, can be used.

The random variable for the parameter vector is denoted by Θ and its (prior) probability distribution by p Θ (ϑ). By virtue of the parameters and initial conditions for the forward model, equation [START_REF] Fishman | Monte Carlo: concepts, algorithms and applications[END_REF], being random variables the state vector and the observation vectors too are random variables. The random variables for the state vector and the observation vector at time t are denoted by X(t) and Z(t), respectively. Consequently, in equations ( 15) and ( 16), the variables x(t), ϑ, z(t), and n(t) are realisations of the corresponding random variables X(t), Θ, Z(t), and N(t) (the noise process), respectively. The question of interest in this article is to estimate how much information about the parameters, Θ, is contained in a series of observations of Z(t).

Let Z j denote the random observation vector at time t j , and z j ∈ R Nz denote the measurement of this random vector (one realisation). Similarly, let Z 1:n denote the set of n observation vectors at discrete time instants labelled t 1 to t n , and z 1:n ∈ R Nz×n denote the corresponding realisations. Starting from a prior probability distribution of the parameters p Θ (ϑ), the posterior probability distribution of the parameters evolves as follows by successive observations Z j = z j p Θ (ϑ) → p Θ|Z 1 (ϑ|z 1 ) → p Θ|Z 1:2 (ϑ|z 1:2 ) → . . . → p Θ|Z 1:j (ϑ|z 1:j ) → . . . → p Θ|Z 1:n (ϑ|z 1:n ) For the above probability distributions, the differential entropy evolves as

H(Θ) → H(Θ|Z 1 ) → H(Θ|Z 1:2 ) → . . . → H(Θ|Z 1:j ) → . . . → H(Θ|Z 1:n )
where, as before for an a priori analysis, the conditional differential entropies H(Θ|Z 1:j ) are averaged quantities, see equation [START_REF] Huan | [END_REF], over all realisations of Z 1:j

H(Θ|Z 1:j ) = Z 1:j T p Θ|Z 1:j (ϑ|z 1:j ) log 1 p Θ|Z 1:j (ϑ|z 1:j ) p Z 1:j (z 1:j ) dϑ dz 1:j , (17)
where Z 1:j and T denote the support of p Z 1:j (z 1:j ) and p Θ (ϑ), respectively. Following this, the quantity of interest -the expected gain in information-evolves as

0 → ∆H(Θ|Z 1 ) → ∆H(Θ|Z 1:2 ) → . . . → ∆H(Θ|Z 1:j ) → . . . → ∆H(Θ|Z 1:n )
where each ∆H(Θ|Z 1:j ) = I Θ|Z 1:j = M (Θ; Z 1:j ) is the difference

∆H(Θ|Z 1:j ) = I Θ|Z 1:j = M (Θ; Z 1:j ) = H(Θ) -H(Θ|Z 1:j ). ( 18
)
While each of the above expected gains in information is non-negative (see remark 2.3), the evolution of expected gains in information as successive observations are added is also non-negative.

Proposition 3.1. H(Θ|Z 1:j ) is non-increasing for increasing values of j.

Proof. Consider the difference between two successive conditional entropies H(Θ|Z 1:j+1 ) and H(Θ|Z 1:j )

D = H(Θ|Z 1:j+1 ) -H(Θ|Z 1:j ) (19)
From equation (17), D can be written as

D = Z 1:j+1 T p Θ,Z 1:j+1 (ϑ, z 1:j+1 ) log p Θ|Z 1:j (ϑ|z 1:j ) p Θ|Z 1:j+1 (ϑ|z 1:j+1 ) dϑ dz 1:j+1 (20)
D = Z 1:j+1 T p Θ,Z 1:j+1 (ϑ, z 1:j+1 ) log p Θ|Z 1:j (ϑ|z 1:j ) p Z 1:j+1 (z 1:j+1 ) p Θ,Z 1:j+1 (ϑ, z 1:j+1 ) dϑ dz 1:j+1 (21) Since log(x) is a concave function, the inequality E[log(x)] ≤ log(E[x]), where E[•]
denotes expectation, holds due to Jensen [26]. Applying this inequality to the expression for D in the above equation yields

D ≤ log Z 1:j+1 T p Θ,Z 1:j+1 (ϑ, z 1:j+1 ) p Θ|Z 1:j (ϑ|z 1:j ) p Z 1:j+1 (z 1:j+1 ) p Θ,Z 1:j+1 (ϑ, z 1:j+1 ) dϑ dz 1:j+1 (22) D ≤ log Z 1:j+1 T
p Θ|Z 1:j (ϑ|z 1:j ) p Z 1:j+1 (z 1:j+1 ) dϑ dz 1:j+1

D ≤ log Z 1:j T p Θ,Z 1:j (ϑ, z 1:j ) dϑ dz 1:j D ≤ log (1) = 0
A direct consequence of the above proposition is that the expected gain in information represented by equation ( 18) is non-decreasing as more (increasing j) observations are available. Intuitively, this is reasonable as the addition of more observations should not result in a decrease in the information that has already been gained. It is also reasonable that D in equation ( 19) should be zero if and only if the observation Z j+1 is totally uncorrelated to the parameter Θ given the observations Z 1:j . Alternatively, once Z 1:j are observed, if there remains no dependence between the random variables Z j+1 and Θ, then no information can be gained about the parameter by additionally observing Z j+1 .

Remark 3.1. The equality in equation (22) occurs when the numerator and denominator in the log term of equation ( 21) are identical, i.e. when p Θ,Z 1:j+1 (ϑ, z 1:j+1 ) = p Θ|Z 1:j (ϑ|z 1:j ) p Z 1:j+1 (z 1:j+1 )

= p Θ|Z 1:j (ϑ|z

1:j ) p Z 1:j |Z j+1 (z 1:j |z j+1 ) p Z j+1 (z j+1 ) (24) = p Θ (ϑ) p Z 1:j |Θ (z 1:j |ϑ) p Z j+1 |Z 1:j (z j+1 |z 1:j ). ( 25
)
Equations ( 24) and [START_REF] Kraskov | [END_REF] show the familiar factorisation of the joint probability distribution of a Markov chain. Hence, if a Markov chain of the form Θ → Z 1:j → Z j+1 is formed between the random variables Θ, Z 1:j , and Z j+1 , then there is no information gain about Θ by observing Z j+1 when Z 1:j have already been observed. One can also view this result in the form of the Markov property that conditioned on Z 1:j the variables Θ and Z j+1 are independent and hence carry no information about each other.

Remark 3.2. The quantity of interest presented in equation ( 18) can also be written as

∆H(Θ|Z 1:j ) = M (Θ; Z 1:j ) = H(Θ) + H(Z 1:j ) -H(Θ, Z 1:j ). ( 26
)
where

H(Θ, Z 1:j ) = Z 1:j T p Θ,Z 1:j (ϑ, z 1:j ) log 1 p Θ,Z 1:j (ϑ, z 1:j ) dϑ dz 1:j (27)
is the joint entropy of the parameters and the observables, and

H(Z 1:j ) = Z 1:j p Z 1:j (z 1:j ) log 1 p Z 1:j (z 1:j ) p Z 1:j (z 1:j ) dz 1:j (28)
and

H(Θ) = T p Θ (ϑ) log 1 p Θ (ϑ p Θ (ϑ) dϑ (29)
are the marginal entropies of the observations and the parameters, respectively. The structure of equation ( 26) is useful as it does not contain any conditional entropy terms and all the quantities in the RHS have the same functional form, namely that they are pure continuous Shannon entropies of the joint random variable [Θ, Z 1:j ] and the marginal random variables Θ and Z 1:j . Consequently, if one has an estimator for continuous Shannon entropy of a random variable, then the expected gain in information can be evaluated using equation (26).

Remark 3.3. Although Θ in equation (26) represents the random vector corresponding to all the parameters, Θ = [Θ (1) , Θ (2) , . . . , Θ (N ϑ ) ], of the dynamical system represented by equation [START_REF] Fishman | Monte Carlo: concepts, algorithms and applications[END_REF], it can be used to calculate the expected gains in information for any combinations of parameters. For example, expected information gains for an individual parameters Θ (l) is given by equation ( 26) when Θ is replaced by Θ (l) . Similarly, expected information gain of two parameters Θ (l) and Θ (m) considered together is given by ∆H(Θ (l) , Θ (m) |Z 1:j ) = H(Θ (l) , Θ (m) ) + H(Z 1:j )) -H(Θ (l) , Θ (m) , Z 1:j ).

From expected information gain to assessing identifiability

Consider a single parameter Θ (l) of the dynamical system. From section 3, the expected gain in information about this parameter by observing Z 1:j is ∆H(Θ (l) |Z 1:j ), where 1 ≤ j ≤ n represents the time when the observation Z j is taken. This section deals with interpretation of the expected gain in information in relation to identifiability of the parameter Θ (l) . It is assumed that one can numerically estimate this information gain (such methods are presented in section 6). An intuitive interpretation of the expected gain in information relies on its magnitude. One would expect high ∆H(Θ (l) |Z 1:j ) to imply that a substantial amount of information about the parameter Θ (l) is contained in the observations Z 1:n , and hence Θ (l) to be identifiable. One would also expect the unidentifiable parameters to show little to no expected gain in information. Furthermore, the time-intervals t a:b (corresponding to Z a:b ) where ∆H(Θ (l) |Z 1:j ) shows most increase should be the time-intervals where the measurements of Z a:b are most informative about the parameter. This would hint that addition of more observations in this time-interval can lead to a better estimate of the parameter.

While the evolution of ∆H(Θ (l) |Z 1:j ) certainly conveys how easy or difficult (based on its magnitude) it is to identify a particular parameter Θ (l) from observations of Z, it is relatively hard to make physical sense of a statement such as "by observing Z 1:n at times t 1:n , the expected gain in information for a particular parameter is q nats". Unless q takes extremely high (or extremely low) values relative to the prior entropy H(Θ (l) ), thereby implying certain identifiability (or unidentifiability), the magnitude of q nats in itself has little physical interpretation in terms of whether the parameter is identifiable or not. Indeed, if two parameters have similar prior entropies and resulted in q and r nats of expected information gain, the relative magnitudes of q and r can be used to compare the identifiability of one parameter against the other, i.e. to ascertain how easy/difficult it is to identify one parameter with respect to the other.

In order to assign a physical interpretation to the expected information gain in relation to identifiability, a concept based on the intuitive and physically interpretable argument of direct observation is introduced. The following hypothetical question is posed:

Suppose there existed a hypothetical device that could directly measure the parameter Θ, i.e. it measures the parameter corrupted by an additive Gaussian noise

Z = Θ + Λ Λ ∼ N (0, σ 2 e ) ( 30
)
where Z represents the hypothetical measurement and σ 2 e represents the variance of zero-mean noise Λ in the hypothetical device. Given the gain in information that we expect from our complex forward and observation processes, what is the σ 2 e that would result in precisely the same expected gain in information for a single measurement through the hypothetical device? .

The hypothetical noise variance σ 2

e is referred as the entropy equivalent variance of a single direct observation (EEV) and represents the precision of the hypothetical instrument that is used to directly measure the individual parameter Θ. It converts the hard-to-interpret nats of information to the variance of a measurement device with Gaussian error statistics. Since, it converts the complex forward and observation operators -that hide the manner in which information is gained about the parameters from the observables -to a simple method of directly measuring the parameter only once, it is readily interpretable.

EEV is meaningful as it offers an alternate/equivalent measure for interpreting nats of expected information gain. For instance, if the initial (prior) parameter variance is 10 -1 square-units and σ 2 e for this parameter is 10 2 square-units, then the interpretation is that if one were to measure this parameter directly with a hypothetical measurement device, the error in such measurement would be ±20.0 (± two standard deviations) units. Comparing this with the error in prior knowledge of ± 0.6 units, one may conclude that such an instrument is clearly futile to identify the parameter. Consequently, by the equivalence of entropy between the real measurements and the hypothetical measurements, the parameter can be deemed unidentifiable. On the other hand if σ 2 e is 10 -3 square-units then the hypothetical instrument's measurement precision is ±0.06 units and in this case the parameter can be deemed identifiable. In what follows the formulation of σ 2 e is presented. In this formulation, the entropies of the hypothetical devices are represented by H and the entropies of the real observation process are represented by H. Consider a parameter ϑ, the corresponding random variable Θ, and related observables Z 1:n measured at times t 1:n . Assume that the real expected gain in information, ∆H(Θ|Z 1:n ), has been estimated (method presented in section 6). Let the prior entropy Θ be H(Θ). If the prior probability distribution of this parameter were Gaussian with mean µ 0 and variance

σ 2 0 p Θ (ϑ) = N (ϑ : µ 0 , σ 2 0 ) = 1 σ 0 √ 2π exp - (ϑ -µ 0 ) 2 2σ 2 0 , (31)
then its entropy is given by H(Θ) = 1 2 log (2 π e σ 2 0 ). This can be equated to H(Θ) in order to obtain an entropy equivalent σ 2 0

H(Θ) = 1 2 log 2 π e σ 2 0 = H(Θ) (32)
Consider the direct observation process for Θ presented in equation (30). By Bayes' rule the posterior probability distribution of Θ by observing

Z = z is p Θ|Z (ϑ|z) = p Z|Θ (z|ϑ) p Θ (ϑ) p Z (z) . ( 33
)
From equations ( 30) and ( 31) it is easy to see that the above conditional distribution is Gaussian

p Θ|Z (ϑ|z) = N Θ : µ 1 , σ 2 1 , (34)
where

µ 1 = σ 2 e σ 2 0 + σ 2 e µ 0 + σ 2 0 σ 2 0 + σ 2 e z, (35)
and

1 σ 2 1 = 1 σ 2 0 + 1 σ 2 e . ( 36
)
The entropy of this posterior distribution is

H (Θ| Z = z) = 1 2 log 2 π e σ 2 1 . (37)
Since the above conditional entropy does not depend on the value z of the observation, the average conditional entropy (expectation over all values of Z) is identical to the above

H (Θ| Z) = Z H (Θ| Z = z) p Z (z) dz = 1 2 log 2 π e σ 2 1 . (38)
Hence, the decrease in entropy (increase in information) from equations ( 32) and ( 38) is

∆H(Θ|Z) = H(Θ) -H(Θ|Z) = 1 2 log σ 2 0 σ 2 1 ( 39
)
Substituting σ 2 1 from equation (36) in the above equation yields

∆H(Θ|Z) = 1 2 log 1 + σ 2 0 σ 2 e ( 40
)
In order to obtain the same entropy decrease in the hypothetical measurement and the real measurements, ∆H(Θ|Z) must be equal to ∆H(Θ|Z 1:n ), which yields

σ 2 e = σ 2 0 exp {2 ∆H(Θ|Z 1:n )} -1 (41)
Equation ( 41) depends on the expected gain in information and the equivalent prior variance that depends on the prior entropy, H(Θ), of the parameter through equation (32). It is interesting to see from the above equation that the following ratio

σ 2 u = σ 2 e σ 2 0 = 1 exp {2 ∆H(Θ|Z 1:n )} -1 (42)
depends only on the information gain by measuring Z 1:n . When the prior variance σ 2 0 is unity then σ 2 u = σ 2 e , and is therefore termed the entropy equivalent variance of a single direct measurement with respect to unit prior variance (EEVU). This is a standardised/normalised measure which enables easy comparison of equivalent variances between different parameters irrespective of their prior entropy magnitudes. It is meaningful due to the property of invariance of the mutual information (also the expected gain in information) under homeomorphic transformations. Following remark 2.4, since the expected gain in information remains invariant under homeomorphic transformations, one can imagine that for each parameter Θ there exists such a transformation under which the prior distribution is transformed to a distribution with unit σ 2 0 . While the existence of such a transformation is a hypothesis, the interpretation of σ 2 u under such a transformation standardises the concept of equivalent variances. To conclude this concept, irrespective of the prior distributions chosen for the parameterswhich might differ both in form and location -σ 2 u represents the error of the hypothetical measurement device under an assumption that all the prior distributions have been transformed to distributions with unit σ 2 0 .

Remark 4.1. It should be noted that σ 2 e in the above is the variance of the Gaussianerror measurement, and not the posterior expected variance of the parameters. The latter, denoted by σ 2 f , is given by σ 2 f = σ 2 e σ 2 0 /(σ 2 e + σ 2 0 ), see equation (36).

Noise process of choice

In the above section the choice of Gaussian noise for the hypothetical measurement device is purely a result of the universal ease of interpretability of the Gaussian distribution. In cases where both the prior and posterior distributions are strongly non-Gaussian, the measurement error of the hypothetical device can be described by user's choice of noise process. Let the hypothetical observation process be

Z = Θ + Λ(φ) ( 43
)
where φ are the parameters of the noise process that one wishes to infer. In the above, additivity of noise process is not necessary; any transformation of Θ to Z through noise Λ(φ) can be used. The argument is same as before: what are the parameters φ of the observation process Λ(φ) that result in the same gain in information as by measurement of the real observables Z 1:n . While analytically intractable in most cases, this problem can be solved in an optimisation (minimisation) framework. The idea is that the numerical method to estimate the expected gain in information (presented later in section 6) relies only on samples of the random variables. Given a guess of the noise parameters φ, one can obtain samples of the noise process Λ(φ) and use equation ( 43) to obtain samples of Z from the samples of Θ. From the samples of Θ and Z the expected information gain ∆H(Θ|Z) can be obtained by the same method (section 6) that is used to estimate the expected information gain ∆H(Θ|Z 1:n ) of the real observation process. Our goal is to find φ such that ∆H(Θ|Z(φ)) = ∆H(Θ|Z 1:n ). The optimal noise parameter vector can hence be found by minimising the L-2 norm of the difference between the real information gain ∆H(Θ|Z 1:n ) and the information gain obtained through the hypothetical noise process ∆H(Θ|Z(φ))

φ = argmin φ || ∆H(Θ|Z(φ)) -∆H(Θ|Z 1:n ) || 2 (44)
The above minimisation problem can be solved by any optimisation algorithm of user's choice, and is always in 2-dimensions for calculation of the expected information gain ∆H(Θ|Z(φ)) as both Θ and Λ(φ) are scalar random variables. This makes the minimisation problem relatively less expensive when compared to estimating the real ∆H(Θ|Z 1:n ) in the first place (which is in n + 1 dimensions), and offers the users to interpret the results consistent with their choice of priors and noise process.

Mutual information and correlation

While the analysis presented in sections 3 and 4 can be used to assess identifiability of individual parameters, it does not reveal if two or more parameter estimates are correlated. High correlation among two or more parameters implies that these parameters can only be estimated together (not separately) [3]. This section deals with assessment of such situations. While only the case of correlation between two parameters is considered here, its extension to three or more groups of parameters is not difficult. Suppose two parameters, Θ (l) and Θ (m) are individually unidentifiable based on the analysis presented in sections 3 and 4. A reasonable question then to ask is the following: "After observing Z 1:j , suppose one had an estimate -perhaps the true value or a reasonable estimate obtained through other means -of one of the parameters, say Θ (l) . Given that Θ (l) is now known how much more information does this knowledge carry about the parameter Θ (m) , and how significant is this information regarding the identifiability of parameter Θ (m) ". This question is answered by means of correlation between Θ (l) and Θ (m) and the quantity conditional mutual information answers precisely this question. In section 3, for individual identifiability analysis, the mutual information between each parameter and the observations is considered , see equation [START_REF] Jaynes | Information Theory and Statistical Mechanics[END_REF]. On the contrary, in reference to the aforementioned question, the quantity of interest in this section is the conditional mutual information, i.e. mutual information between Θ (l) and Θ (m) given Z 1:j . This quantity, denoted as M (Θ (l) ; Θ (m) )|Z 1:j ), reflects the sum of two correlations: the prior correlation between Θ (l) and Θ (m) , and the expected correlation between Θ (l) and Θ (m) that is developed due to the observations of Z 1:j . Alternatively, having observed Z 1:j , it quantifies the average amount of information that would be gained about one of the parameters by observing the other.

From equation (13), for the two random variables Θ (l) and Θ (m) mutual information, M (Θ (l) ; Θ (m) ), is defined as follows

M (Θ (l) ; Θ (m) ) = Θ (m) Θ (l) p Θ (l) ,Θ (m) (ϑ (l) , ϑ (m) ) log p Θ (l) ,Θ (m) (ϑ (l) , ϑ (m) ) p Θ (l) (ϑ (l) ) p Θ (m) (ϑ (m) ) dϑ (l) dϑ (m) (45)
From the above, the conditional mutual information for the observation

Z 1:j = z 1:j is M (Θ (l) ; Θ (m) |Z 1:j = z 1:j ) = Θ (m) Θ (l) p Θ (l) ,Θ (m) |Z 1:j (ϑ (l) , ϑ (m) |z 1:j ) log p Θ (l) ,Θ (m) |Z 1:j (ϑ (l) , ϑ (m) )|z 1:j ) p Θ (l) |Z 1:j (ϑ (l) |z 1:j ) p Θ (m) |Z 1:j (ϑ (m) |z 1:j ) dϑ (l) dϑ (m) . ( 46
)
and its expected value over all all possible observations, the conditional mutual information, M (Θ (l) ; Θ (m) )|Z 1:j ) is

M (Θ (l) ; Θ (m) |Z 1:j ) = Z 1:j M (Θ (l) ; Θ (m) |Z 1:j = z 1:j ) p Z 1:j (z 1:j ) dz 1:j (47)
By methods similar to those presented in section 3, it is easy to prove that M (Θ (l) ; Θ (m) |Z 1:j ) ≥ 0. In the discussion that follows, it is assumed that the prior correlation between the parameters is zero and hence M (Θ (l) ; Θ (m) |Z 1:j ) reflects purely the correlation developed between the parameters due to the observations. Consequently, if the observations of Z 1:j induce any correlation between the parameters then M (Θ (l) ; Θ (m) |Z 1:j ) > 0, while if no correlation is induced then M (Θ (l) ; Θ (m) |Z 1:j ) = 0. The former happens when the parameters Θ (l) or Θ (m) form a common effect for the observations Z 1:j ; for example, if the observable equation in terms of the parameters has a term containing a combination of the two parameters. For instance if z(t) = g(t, ϑ) + f (t)/(ϑ l ϑ m ) + (t), and if the second term in the RHS is dominant then it is difficult to differentiate between θ l and θ m by measuring z 1:j and consequently to identify them separately (their product, however, can be possibly identified together). Such a case is presented in section 7.1. For increasing j, unlike the quantity ∆H(Θ (l) |Z 1:j ) = M (Θ (l) ; Z 1:j ) which is monotonically non-decreasing, the evolution of M (Θ (l) ; Θ (m) |Z 1:j ), is not necessarily monotonic. The difference between two successive conditional mutual informations is given by

D M = M (Θ (l) ; Θ (m) |Z 1:j+1 ) -M (Θ (l) ; Θ (m) |Z 1:j ) (48) = M (Θ (l) ; Θ (m) |{Z 1:j , Z j+1 }) -M (Θ (l) ; Θ (m) |Z 1:j ), (49)
and reflects how much more information about Θ (l) can be obtained if Z j+1 was also observed after having observed both Z 1:j and Θ (m) . This difference can be positive, negative, or zero, and is referred as interaction information in literature.

The evolution of M (Θ (l) ; Θ (m) |Z 1:j ) can reveal interesting features of the problem under consideration. The regions of time where M (Θ (l) ; Θ (m) |Z 1:j ) is increasing reflect the time intervals where the observations of Z j induce a correlation between the parameters, and the regions where M (Θ (l) ; Θ (m) |Z 1:j ) decreases reflect the time intervals where the observations of Z j destroy the hitherto built correlation between the parameters. The former happens when the observations are a common effect of the parameters while the latter happens when the observations are a common cause. When the common effect is fixed (known/observed) a correlation is induced (referred as explaining away in the Bayesian literature), and when the common cause if fixed, the correlation is destroyed as the observation partly explains or accounts for the hitherto built correlation. The latter can be thought of as the situation where the observation of Z j+1 makes part of the shared information between Θ (l) and Θ (m) (developed due to the observations of Z 1:j ) redundant.

Remark 5.1. In terms of the differential entropies, the conditional mutual information can be written as

M (Θ (l) ; Θ (m) |Z 1:j ) = H(Θ (l) , Z 1:j ) + H(Θ (m) , Z 1:j ) -H(Z 1:j ) -H(Θ (l) , Θ (m) , Z 1:j ). (50)
In light of remark 3.2, the above form is useful as it, too, contains only terms of continuous Shannon (differential) entropy. Consequently, if an estimator of differential entropy is available, then it can be used to estimate both mutual information and conditional mutual information through equations ( 26) and (50), respectively.

Numerical methods to estimate (conditional) mutual information

From sections 3 and 5, the quantities of interest to be evaluated are the mutual information, M (Θ; Z 1:j ) = ∆H(Θ|Z 1:j ), and the conditional mutual information M (Θ (l) ; Θ (m) |Z 1:j ). Furthermore, through equations ( 26) and (50) both these quantities can be expressed as a combination of the continuous Shannon (differential) entropies. In light of remark 5.1, in this section, numerical methods to estimate differential entropies are presented.

To reiterate, the differential entropy of a vector random variable X is given by

H(X) = X p X (x) log 1 p X (x) dx, (51)
where p X (x) is the probability density function of X and has a support X . If the analytical form of the probability density p X (x) is known then one might employ numerical integration techniques to estimate the differential entropy. However, the probability distributions of interest are generally not tractable due to high non-linearity in the forward and observation operators. Furthermore, when X is Z 1:j then as j increases the dimension of the vector random variable increases, which makes analytical tractability impractical in most cases. Nonetheless, given the forward and observation models, i.e. equations ( 15) and ( 16), and prior probability distributions on the parameters and initial conditions of the forward model, one may be able to simulate a large, yet finite, number of particles (samples from the prior distribution) forward in time to obtain samples from the probability distribution of Z 1:j . Hence, methods to estimate H(X) only from samples of X are sought. An overview of non-parametric estimation of entrapy estimation methods is presented in [27]. In what follows, two such concepts are briefly described. Entropy estimation methods can be broadly classified into the following categories (i) Kernel-density estimators based: In such methods p X (x) is first approximated by a kernel-density estimate (KDE). Such KDE is constructed through the samples of x, and is used in equation (51) to evaluate the integral (either by numerical integration or Monte-Carlo estimation). See references [28,29] for details.

(ii) Nearest neighbours based: Such estimates are based on approximating the probability density p X (x) at each sample point through the distances to its nearest neighbours, and then using equation ( 51) to obtain a Monte-Carlo estimate of H(X). Kozachenko and Leonenko [30] proposed such an entropy estimator estimator based on the first nearest neighbour. This estimator was extended to more robust k-nearest neighbours based estimators by Singh et. al, [31] and Kraskov et. al. [START_REF] Kraskov | [END_REF].

In this article, the estimator proposed by Kraskov et. al. [START_REF] Kraskov | [END_REF] after a minor modification to aid the comparison of entropies in increasing dimensions is employed. These estimators are presented next. If one has an approximation p(x i ) of the probability density p(x i ) then the entropy in equation ( 51) can be estimated using a Monte Carlo approximation as follows

Ĥ(X) = 1 N ens Nens i=1 log 1 p(x i ) = - 1 N ens Nens i=1 log p(x i ) , (52)
where x i represent the N ens number of samples from the probability density p(x). In k-nearest neighbour based entropy estimators the estimator p(x) is constructed using the k-nearest neighbours of x i . Kraskov et. al. [START_REF] Kraskov | [END_REF] proposed the following estimator for log (p(x i ))

log p(x i ) = ψ(k) -ψ(N ens ) -N x E [log( )] -log(c d ), (53)
where 2 Nx Γ(1+Nx/2) , where Γ represents the gamma function . Based on this, the estimator for H(X) can be written from equation (52) as

ψ(•) is the digamma function, E[•] is the expectation operator, N x is the dimension of x (x ∈ R Nx ),
Ĥ(X) = ψ(N ens ) -ψ(k) + log(c d ) + N x N ens Nens i=1 log ( (i)) (54)
where (i) represents twice the distance from the i th sample to its k th neighbour.

Given the above differential entropy estimator, the mutual information between two random variables X ∈ R Nx and Y ∈ R Ny can be calculated from the samples (x i , y i ) is M (X; Y) = H(X) + H(Y) -H(X, Y). However, Kraskov et. al. pointed out that such an approach leads to biased results as the k-nearest neighbour distances are calculated in different spaces -the joint and the marginal -which use different distance scales. To remedy this, they suggested that while equation (54) can be used to estimate the entropy in the joint space, the distance scale in the joint space can be used in the marginal spaces too in order to cancel the biases. This can be achieved by calculating the k th -neighbour distances (i) in the joint space and then finding the number of nearest neighbours within this distance in the marginal spaces for marginal entropy estimation. Thus, one uses a variable k in the marginal spaces based on the distance calculated in the joint space. The entropy estimators for the joint and the marginal spaces can then be written as

Ĥ(X, Y) = ψ(N ens ) -ψ(k) + log(c dx c dy ) + N x + N y N ens Nens i=1 log J (i) (55)
and

Ĥ(X) = ψ(N ens ) - 1 N ens Nens i=1 ψ (k x (i) + 1) + log(c dx ) + N x N ens Nens i=1 log J (i) (56) Ĥ(Y) = ψ(N ens ) - 1 N ens Nens i=1 ψ (k y (i) + 1) + log(c dy ) + N y N ens Nens i=1 log J (i) (57)
where for each sample point (x i , y i ) the distance J (i) is twice the distance to the k-th nearest neighbour in the joint space, k x (i) and k y (i) represent the number of points in the marginal spaces of X and Y that are at distances less than J (i)/2 from the points x i and y i , respectively, and c d (•) represents the volume of a unit ball in the dimension of (•). The mutual information estimator then simplifies to

M (X, Y) = ψ(N ens ) + ψ(k) - 1 N ens Nens i=1 ψ (k x (i) + 1) - 1 N ens Nens i=1 ψ (k y (i) + 1) (58)

Modification of the estimator

As mentioned in section 3 the quantity of interest in order to quantify the information gain about the parameters is the evolution of M (Θ; Z 1:j ), see equation ( 26) for increasing j. Equations ( 55), ( 56) and ( 57) can be used to estimate this mutual information, equation (58), without incorporating additional biases due to different distance scales in the joint and marginal spaces. However, a comparison of M (Θ; Z 1:j ) with M (Θ; Z 1:j+1 ) is unfair as the space (Θ, Z 1:j ) is of lower dimension than that of (Θ; Z 1:j+1 ). This is particularly noticeable (see Figure 4 in section 7.1) in regions of time where no information is available about the parameters, i.e. regions where M (Θ; Z 1:j ) remains constant for changes in j. In such regions the the evolution of M (Θ; Z 1:j ) can decrease for increasing j, which is contrary to the non-decreasing nature of M (Θ; Z 1:j ), see proposition 3.1. To avoid the above problem a minor modification in the entropy estimator is proposed. Once the sampling plan is defined, i.e. once the times at which observations are available is fixed as {t j }, 1 ≤ j ≤ n, the largest joint space is of (Θ, Z 1:n ). In this largest space, for each sample point (ϑ i , z i 1:n ), let M (i) represent twice the distance to its k-th nearest neighbour. It is proposed that to calculate the mutual information estimates of M (Θ; Z 1:j ) for all 1 ≤ j ≤ n, the distances used be M (i). Following equation (58), these estimates can be written as

M (Θ; Z 1:j ) = + 1 N ens Nens i=1 ψ k ϑ,z 1:j (i) + 1 - 1 N ens Nens i=1 ψ k z 1:j (i) + 1 - 1 N ens Nens i=1 ψ (k ϑ (i) + 1) (59)
where for each sample point (ϑ i , z i 1:j ), k ϑ,z 1:j (i) represents the number of points at distances less than M (i)/2 in the joint space of (Θ, Z 1:j ). Similarly, in the marginal spaces of Θ and Z 1:j , k ϑ (i) and k z 1:j (i) represent the number of points at distances less than M (i)/2 from the sample points ϑ i and z i 1:j , respectively. Since all the entropy estimators use the same distance scale M , the successive mutual information estimates, M (Θ; Z 1:j ), can be compared to each other. In what follows, this modified mutual information estimator is used to study the evolution of gain in parameter information as more observations are available.

Lastly, the conditional mutual information M (Θ (l) ; Θ (m) |Z 1:j ) can be estimated through equation ( 50) by estimating all the differential entropies on the RHS of equation (50) through the estimator of Kraskov et. al. Here too, for comparison of M (Θ (l) ; Θ (m) |Z 1:j ) for different j, the k th -th neighbour distances can be calculated in the largest space of (Θ (l) , Θ (m) , Z 1:n ) and all terms on the RHS of equation ( 50) can be estimated by finding the number of neighbours within this distance in the marginal spaces.

The numerical method

Given a dynamical system and its numerical discretisation the proposed approach consists of the following steps:

(i) generate N ens samples from the prior distributions of the parameters and initial conditions.

(ii) for each sample compute the numerical solutions.

(iii) compute the corresponding observations by generating samples from the observation noise process.

(iv) compute the (conditional) mutual informations using the modified estimator presented above.

(v) compute corresponding equivalent variances.

The pseudo-code for such computation is shown in algorithm 1. One difficulty associated with the above presented estimators is that if one of the parameter magnitudes (or one of the observables magnitude) is significantly larger than others then this parameter dominates M . Consequently, in the marginal spaces, k (•) can encompass a large number of sample points -sometimes even the entire ensemble -which leads to erroneous results. In order to alleviate this difficulty the following procedure is performed: since the gain in information is invariant under transformations, see remark 2.4, , scaling and centering (to zero mean and unit variance) of all the components of the parameter and observation vectors is performed before calculating M (i) and number of points inside the -balls. The computational features of the estimator, the convergence with respect to the number of samples and neighbours, as well as its complexity are detailed in [START_REF] Kraskov | [END_REF]. The advantage of using a Monte Carlo approach is that despite its slow convergence rate, the convergence is independent of the dimension of the stochastic space considered. This is of the utmost relevance when considering the present application, in which the dimension is growing. Moreover, it is non-intrusive and it does not require for the direct solvers to be re-written in order to take the stochasticity into account.

Algorithm 1: Algorithm to calculate gain in information for the parameters

Input:

• The forward model: ẋ(t) = F (x(t), t, ϑ);

x ∈ R Nx , ϑ ∈ R N ϑ • The observation model: z(t) = G(x(t)) + n(t); z, n ∈ R Nz
• Prior probability distribution on the parameter vector: p Θ (ϑ).

• Prior probability distribution on the initial state vector: p X (x).

• Number of samples to be simulated: N ens .

• List of measurements times: t j , 1 ≤ j ≤ n.

• Number of nearest neighbours to used for mutual information estimation: k

• The set S containing parameter combinations for which ∆ Ĥ0→1:j is estimated. For example, S = {ϑ (1) , ϑ (2) , (ϑ (1) , ϑ (2) )}

Output: ∆ Ĥ0→1:j , 1 ≤ j ≤ n, for all elements of S. In this section the numerical experiments are commented, of increasing complexity, aiming at validating the proposed approach.

The first example deals with an RCR circuit, which is an example of 0-D model, often used as boundary condition for simulation of cardiovascular flows, and whose identification is a key step to the setting up of realistic simulations.

Then, an epidemiological SIR model is investigated aiming at comparing the results with those obtained in the literature by means of more complex differential algebra methods. This is still an ODE model, but it is non-linear, with a large number of parameters.

A first example with PDEs is on a potential identification problem for an hyperbolic system in 1 + 1 dimensions.

In the last section, the identification of sources for a problem of advection diffusion in 2D is commented.

For all the testcases presented, the number of direct simulation computed for the MC estimation of the entropies were empirically evaluated by considering the convergence of the results. The experiments presented uses a number of particles such that the estimated quantities vary less than 0.01 when the number of particles is doubled.

The number of k-neighbours has been adjusted in an analogous way.

Three-element Windkessel model

Parameter estimation in a three-element Windkessel model [32], a commonly imposed boundary condition to study haemodynamics [33], is considered. This model, as shown in Figure 1a, consists of three elements: a proximal resistance, R p , representing the sum of large vessels; a capacitance, C, representing the elastance of large vessels; and a distal resistance, R d , representing the resistance of smaller vessels in the microcirculation. The behaviour of the Windkessel model of Figure 1a is governed by the following differential algebraic equations

P i (t) -P m (t) = R p q i (t), ( 60
)
P m (t) -P ven (t) = R d q o (t), ( 61
)
q i (t) -q o (t) = C d dt P m -P ext (t) , (62)
where q i and q o represent the inlet and outlet flow-rates, P i and P m represent the inlet and the mid-Windkessel pressure, and P ven and P ext represent the venous and external pressures, respectively. P ven and P ext are both assumed to be zero. The inlet flow-rate, the control curve, is known and is shown in Figure 1b. The goal of this analysis is to determine the identifiability of the Windkessel parameters -R p , C, and R d -by discrete measurements of P i (t) where z(t j ) represents the measurement of P i (t j ) at time t j , and ˜ (t j ) represents the associated measurement noise. To maintain positivity of the Windkessel parameters and pressure the following transformation is used for these variables

z(t j ) = P i (t j ) + ˜ (t j ), ( 63
) P i R p q i R d q o C q i -q o P ven P ext P m (a)
Ψ (•) = Ψ (•) 0 2 ψ (•) where (•) ∈ {P i , R p , C, R d } (64)
where

Ψ (•)
o represent a reference value, Ψ (•) represents the real values of P, R p , C, R d which can lie between 0 and ∞, and ψ represents the transformed values which can take any value between -∞ and ∞. The advantage of such transformation is that the probability distributions on {P i , R p , C, R d } in the ψ-space can be assumed to have support over the entire real line while ensuring that the model is physically correct in terms of positivity of these parameters. Moreover, while the assumption of a Normal process for (t i ) in equation ( 63) would have been wrong by leaving a finite probability that the pressure measurement, z(t i ), could be negative, the assumption of Normal noise in the measurement of the transformed variable ψ P i is not unreasonable

z(t j ) = ψ P i (t j ) + (t j ),
where n , are considered. These noise levels transform to variances of approximately 2.0, 18.0, 36.0, and 184.0 for the measurement of pressure in the real Ψ-space (a log-normal distribution when a normal distribution is considered the ψ-space).

(t i ) = N (0, σ 2 n ). ( 65
The expected information gains for individual parameters, pairs of parameters considered together, and the triplet consisting of all the three parameters, are shown in Figure 2 The plot for P i represents the mean of the observations generated through the propagated samples. σ 2 n corresponds to the transformed variable ψ P i (equation ( 65)). For σ 2 n = {0.0005, 0.005, 0.01, 0.05}, the corresponding variances of the noise process for the real parameter P i (untransformed, see equation ( 63 1: Real parameters are in the Ψ-space and a parameterisation of the form Ψ = 2 ψ is used to ensure positivity of the parameters and pressure.

magnitudes of expected information gains for the individual parameters it can be concluded that R d is most easily identifiable (largest gain in information) while R p is most difficult to identify (lowest gain in information). Furthermore, it is observed that most of the information about the parameter R p is contained in the sharp systolic phase of pressure (t ∈ [0.4, 0.5]). This is consistent with the results obtained through generalised sensitivity function (GSF) analysis presented in [34,35]. In fact, if one normalises the expected gain in information for the individual parameters to lie between 0 and 1 then the normalised plots look similar to GSFs (see Figure 10 of [35]). While GSFs can oscillate in the presence of parameter correlations, thereby making them relatively hard to interpret, the expected information is always non-decreasing. This property makes the expected information gain easier to interpret. Furthermore, in contrast to GSFs the magnitudes of expected information gain provide a quantitative measure of not only which regions of time contain most information about a particular parameter but also how much.

The effect of increasing noise in Figure 2 is consistent with the intuition that a higher noise should result in a lesser amount of information gain. The EEV values for the parameters for the four levels of noise considered are shown in Table 1. For the lowest level of noise corresponding to noise variance of 1.8 mmHg 2 for P i , it is observed that the σ 2 e for R d is 0.1, implying that relative to the prior variance of 1.0 one obtains a sharp measurement of variance 0.1 through the hypothetical measurement device. One may conclude with these statistics that the parameter R d is easily identifiable. Similar arguments hold for the parameter C which has an σ 2 e = 0.28, a low value compared to the prior unit variance. However, for R p the σ 2 e = 0.97, a value of roughly the same magnitude as the prior unit variance. This implies that even with a low level of noise, since the hypothetical measurement device has an error variance of the same amount as the prior variance, the variance of the final estimate is only halved, see equation (36). This suggests that even in the presence of low noise, the parameter R p is hard to estimate. Furthermore, for the highest level of noise the parameter R p has an σ 2 e = 3.4, implying that virtually no information is contained in the measurements for this level of noise and that the parameter cannot be identified. On the other hand, the parameter R d has an σ 2 e of 0.91 for the highest level of noise, which is lower than the σ 2 e of R p for the lowest level of noise, implying that even with the highest level of noise considered, the uncertainty (judging by the variance) can be halved for R d . Such an analysis can be used in determining the level of acceptable noise in the real measurement device when used for parameter identification. For example, consider two pressure measurement devices of variances 1.8 and 17.8 mmHg 2 . If the goal of buying this measurement device was solely to estimate the Windkessel parameters, then one may conclude that the latter device, even though significantly more error-prone, will do only marginally worse than the former device when it comes to estimating R p , C, and R d . Consequently, in a costbenefit analysis under the assumption that the former device is considerably expensive than the latter, the latter device should be preferred.

Conditional mutual information:

In Figure 2 the expected gains in information when any two parameters are considered together, and when all the three parameters are considered together, are also presented. These are presented for sake of completeness; a utility of these plots lies in cases where the expected gain in information of the pair of parameters considered together is significantly larger than the sum of individual expected information gains. Such a case would imply that one can learn a lot more about the joint parameter distribution compared to the individual marginal distributions. This implies that the parameter estimates are correlated. However, since kNN estimators are used -making it unwise to quantitatively add or subtract the plots shown in Figure 2 since different length scales are used in each plot -it is proposed that that correlations be studied through the conditional mutual information estimated by equation (50) where the same distance scale is used for calculating all the RHS terms. The pair-wise conditional mutual informations are plotted in Figure 3. The first observation is that the behaviour with respect to noise is the same as for the individual information gains. This is expected. The second observation is that the conditional mutual informations for the pairs {R p , C} and {R p , R d } are negligible, implying that these pairs of parameters are not correlated. The evolution of conditional mutual information for the pair {R d , C}, however, shows an interesting pattern. It shows an increasing behaviour in the diastolic phase (t ∈ [0.0, 0.4] ∪ [0.55, 0.7]), and a decreasing behaviour in the systolic phase (t ∈ [0.4, 0.55]). From the discussion presented in section 5, regions of increasing conditional mutual information are regions where the observations build up a correlation between the parameters. This implies that the observations are largely a result of a common effect of the two parameters (hinting that a functional term containing the two parameters drives the observations) thereby building a correlation between the parameters. This makes sense if one looks at the solution to the system presented in equations ( 60) -(62)

M (nats) M (R p ; R d |Z 1:j ) σ 2 n = 0.0005 σ 2 n = 0.005 σ 2 n = 0.01 σ 2 n = 0.05
P i (t) = term-1 (P i (0) -R p q i (0)) e -t/(R d C) + R p q i (t) term-2 + term-3 e -t/(R d C) t 0 e t/(R d C) C q i ( t)d t term-4 (66)
It can be seen that solution consists of many parts where the product of the parameters R d and C comes together as a rate of exponential decay. Their combined effect is much more dominant in the diastolic phase where term-2 and term-4 remain approximately constant (since q i (t) shows little relative variation in these regions, see Figure 1b). Consequently, in these regions, a large part of the solution, P i (t), is determined by the combined effect of these two parameters, which results in an induced correlation by measurements of P i (t). It is interesting to note that even if the solution in equation ( 66) was not known, as typically would be the case for complex dynamical systems, the conditional mutual information plots hint at which functional terms of parameter combinations appear together in the solution or influence the solution together. Furthermore, they isolate regions of time where this combined effect is most dominant. On the other hand, regions of decreasing conditional mutual information are regions where given the observations some of the information that one parameter contains about the other becomes redundant. For example in the systolic phase of pressure (t ∈ [0.4, 0.55]) the built-up correlation between the R d and C is destroyed as observing P i (t) in this region partly accounts for the correlation. After the cycles (if any) of the increase and decrease of conditional mutual information, the value of the conditional mutual information at the end, i.e. when all measurements have been taken, is of high interest. This quantity reflects the final correlation that is induced between the parameters after all the observations have been taken into account. From Figure 3, this final value for the pair {R d , C} is approximately 0.1 nats for the lowest level of noise. This implies that at the end of the experiment, if one fixes/knows the value of one of the parameters, then 0.1 nats of information could be gained about the other parameter on top of the already gained information shown in Figure 2. Since this is only 12.5% more information for the parameter C and 0.8% more information for R d , it is safe to conclude that the correlations are not significant to affect identifiability.

Modified Windkessel model:

In order to demonstrate an example where correlations affect identifiability, a simple modification of the Windkessel model is considered where an extra distal resistance is added (see Figure 5). In this model, R d 1 + R d 2 behaves as R d of the previous model and hence individually, we expect Table 2: Windkessel model with two distal resistances: σ 2 u before and after considering mutual information between the two distal resistances that R d 1 and R d 2 should have lower expected gain in information, but their mutual information should be high. This is exactly what is observed in the numerical results as shown in Figure 6. Compared to the previous example, where the individual gain in information was 1.2 nats, here R d 1 and R d 2 show an expected information gain of only 0.35 nats. On the other hand, the conditional mutual information between R d 1 and R d 2 is now 0.22 nats, which is comparable to 0.35 nats in magnitude, implying that large correlations between these parameters exist. The results of EEV are shown in Table 2. It is observed that these parameters now depict an σ 2 e of roughly the same magnitude as the prior variance. Furthermore, if one adds the conditional mutual information to the original individual gain (implying that one of the two parameters was known/fixed), the σ 2 e drops to half of the original value -a significant increase in identifiability. This further confirms that correlation between R d 1 and R d 2 cannot be ignored.

P i R p q i R d1 R d2 q o C q i -q o P ven=0 P ext=0
Lastly, in Figure 4, the differences between the original estimator by Kraskov et. al. [START_REF] Kraskov | [END_REF] and its modified version, where the distances in the largest space are used, are shown for the parameter R p . It is observed that while the original estimator performs quite well, it shows an undesirable negative slope (decrease in entropy) in regions of time, t ∈ [0.45, 0.5], where entropy remains constant. The modified estimator respects the non-decreasing nature of conditional entropy better than the original estimator.

SIR model.

In this section an application to a non-linear ODE model used in epidemiology (Susceptible-Infected-Recovered, SIR) is presented. This system has been studied in [6] by using a Gröbner basis approach.

The model is a non-linear compact SIR model of the form:

Ṡ = µN -µS - η N SI, (67)
İ = -(µ + ω)I + η N SI, (68)
where µ, N, ω, η are scalar parameters, S and I are the number of susceptible and the infected individuals, respectively. The total number of imdividuals N = S + I + R is a constant parameter for this model, where R is the number of recovered individuals. The state is x = (S, I) and the observation is a fraction of the infected individuals:

z = κI, (69)
where κ ∈ [0, 1] is an unknown parameter. The initial conditions S 0 , I 0 are unknown and can be considered as part of the parameters. Therefore, the parameters vector is: ϑ = (µ, N, η, ω, κ, S 0 , I 0 ). The particularity of this model is that the parameter κ does not affect the dynamics, but only the observations. The level of noise models the number of false positive and negative arising in the diagnosis of the infected individuals that actually go to the medical doctor.

The choice for the prior is discussed hereafter. The total number of individuals has to be positive, so a lognormal prior is assumed of the form: N = N 0 • 2 ϑ N , where N 0 = 100, ϑ N is normally distributed with mean ϑ N = 0 and variance σ 2 N = 0.05. The initial conditions for the susceptible individuals (namely S 0 ) is chosen analogously. Let S 0 = γN . A prior for the fraction γ is chosen to be γ 0 • 2 ϑγ , where and γ 0 = 0.3 and the exponent is a random variable normally distributed, with mean ϑ γ = 0 and variance σ 2 γ = 0.05. The prior for the infected individuals at initial time is chosen to be I 0 = αN , where the prior for α is lognormal, α = α 0 • 2 ϑα , with α 0 = 0.5 and the exponent being normally distributed, with the same statistics as γ. The prior for κ is assumed to be gaussian, with mean κ = 0.4, and variance σ 2 k = 0.1. For the other variables, normal priors are assumed: for µ the mean value is µ = 2.0 and variance σ 2 µ = 0.4; for η: mean η = 1.75, and variance σ 2 η = 0.25; and for the parameter ω, mean ω = 0.6 and variance σ 2 ω = 0.1. The above choice of the priors is made to be reasonably consistent with respect to each other. For instance, the number of individuals in each of the susceptible, infected, and recovered, classes should be less than the total number of individuals.

A number of N ens = 10000 samples is generated and 25 neighbours are used to approximate the entropies. The model is integrated in time by using a second order Crank-Nicolson scheme and 36 observations are retained, uniformly distributed between initial and final time (T = 6).

The expected information gain for the individual parameters are shown in Fig. 7. From the plot it can be inferred that the expected information gains for the parameter N , S 0 , I 0 , and κ are negligible when compared to the other remaining parameters, hinting that they are not identifiable.

In Table 3 the results are shown in terms of EEVs, confirming that the parameters N, κ and the initial conditions for the population are not identifiable. The results obtained are in good agreement with those obtained in [6], showing that for this epidemiologic model, the parameters governing the disease dynamics (namely µ, η, ω) 2). In the first column, the parameters, in the second and third column, the average and the variance of the prior, in the last column the variances of the EEVs. can be identified, while the population characteristic sizes cannot.

Remark 7.1. It should be noted that the choice of priors can influence the outcome of identifiability. This is particularly true for this example where the parameter κ does not affect the dynamics but only the observations. From equation (69) it can be seen that if the prior on I 0 is strong, and the observations at the beginning of the experiment have low noise, then the posterior on κ is relatively precise irrespective of its prior. In such a case, which might not be physically reasonable in terms of the chosen priors, the parameter κ is identifiable.

Potential identification for harmonic waves.

An example where the system dynamics is governed by a PDE is presented in this section. The system is a 1 + 1-dimensional hyperbolic PDE defined for (ξ, t) ∈ [0, 1] × [0, T ], where T = 4 is the final time. The system describes harmonic waves in the presence of an unknown potential field. The system equation reads:

∂ 2 t u = ∂ 2 ξ u + V (ξ)u, ( 70
)
where u is the displacement and homogeneous Dirichlet boundary conditions are imposed u(0, t) = u(1, t) = 0. The system has a deterministic initial condition u(ξ, 0) = sin(πξ). The potential is parametrised by a harmonic expansion as follows:

V (ξ) = A N h j vj sin(πjξ), (71)
where N h = 4 is the maximum number of harmonics considered, A = 20 is a scaling coefficient.

The identifiability of different space wave numbers vj is investigated. The measurements are the displacements u at several points in the domain. The direct simulation is performed by discretising in space by means of piecewise linear Finite Elements (N x = 256 degrees of freedom), in time by an implicit second order Crank-Nicolson scheme, with ∆t = 8 • 10 -3 . The observation consists in measuring the position u in N z = 12 points, uniformly distributed (excluding the extrema), with a time sampling frequency of ν = 12.5 Hz (corresponding to 1 over 10 simulated time instants). By adopting this discretisation, the state x = (u(ξ 1 ), ..., u(ξ N x )) is the value of the solution in the points of the uniform mesh, the observation z = (u(ξ 0 1 ), ..., u(ξ 0 N z )) is the value of the solution in the observation points, corrupted by a gaussian noise with zero mean and variance σ 2 n = 10 -3 . The parameter vector is ϑ = (v 1 , . . . , vN h ). The joint prior for the parameters is assumed to be gaussian, with mean ϑ and covariance matrix S, defined as:

ϑ = (0, 0, 0, 0), (72)
S = diag(1, 1, 1, 1). ( 73
)
In Fig. 8 the expected information gains for the different harmonics are shown. For every parameter, most of the expected information gain occurs in the first part of the evolution, i.e. for t < T w /4, where T w is the time period of the wave. After T w /4, the expected information gain stabilises and there is a continuous learning throughout the dynamical system evolution. This provides a useful information for the setting up of the measuring procedure. As it is expected, low space frequency harmonics are well identifiable, while high space frequencies are not. Indeed, given the observation process, there exists a cut-off on the space frequencies, which is determined by the number of space points monitored. Frequencies higher than 1/N z are not identifiable by construction.

The results are summarised in 3). In the first column, the parameters, in the second and third column, the mean and the variance of the prior, in the last columns the mean value of the prior and the variances of the EEVs.

An analysis of the mutual information between the parameters is performed, showing that all the mutual informations are negligible relative to the values of the expected information gains for the single parameters. The interpretation is that the unidentifiability of high space frequencies does not depend on the correlation between the parameters, but is an intrinsic physical property of the system and the observation process adopted.

Advection-Difffusion equation in 2D.

The last example is a system governed by a 2 + 1 dimensional PDE, that is aimed to mimic a source detection inverse problem for an advection-diffusion equation. A physical point domain is (ξ, t) = (ξ 1 , ξ 2 , t) ∈ Ω = [0, 1] 2 × [0, 0.1], where ξ 1 and ξ 2 are the horizontal and vertical coordinates respectively. The system unknown is a passive scalar c that diffuses and is simultaneously advected by a velocity v. The equation governing the dynamics reads:

∂ t c + v • ∇c = µ∇ 2 c, (74) c
(ξ 1 | (0,1) , ξ 2 , t) = 0 on ∂Ω lat , (75) c
(ξ 1 , 0, t) = S(ξ 1 ) on ∂Ω bot , (76)
c(ξ, 0) = 0. ( 77
)
where S(ξ 1 ) is the source function, defined by:

S(ξ 1 ) = s 1 ξ 1 (1/2 -ξ 1 ) ξ 1 ≤ 1/2 s 2 (ξ 1 -1/2)(1 -ξ 1 ) 1/2 < ξ 1 ≤ 1 (78)
The intensities of the parabolic sources are positive scalar parameters (s 1 , s 2 ). The diffusivity of the passive scalar c is denoted by µ. The velocity is mostly aligned with the ξ 2 axis, so that on the upper boundary characteristics are exiting and no boundary conditions need to be imposed. The probes are located at a constant value ξ 2 = 0.6, uniformly distributed along the horizontal coordinate. A number of N s = 20 probes are considered. This setting is shown in Fig. 9. The equations are discretised by means of standard finite element methods, considering N dof = 10 4 degrees of freedom in space. The equations are integrated in time by using a Crank-Nicolson scheme and taking 800 time steps. The goal of the study is the identifiability of the parameters s 1 , s 2 , µ, v 1 , and v 2 (where v 1 and v 2 are the components of the drift velocity) by measuring the concentration at the probe locations. Two different scenarios are investigated: a drift regime and a diffusive regime. It should be noted that these scenarios are implemented by simply changing the prior. The differences in terms of identifiability are assessed.

Drift regime:

In this section, the drift dominated regime is investigated. Intuitively, owing to low diffusion, the information about the sources is transported downstream relatively unaffected to the probes. The prior is assumed as follows. The source intensities (s 1 , s 2 ) and the diffusivity (µ) have a lognormal prior. For the sources, given s 0 = 25, the following parameterisation for priors for s 1 and s 2 are chosen: s 1 = s 0 2 ϑs 1 and s 2 = s 0 2 ϑs 2 . Here, for the prior distribution, both ϑ s 1 and ϑ s 2 are normally distributed random variables with zero mean and variance 0.2. Similarly, µ = µ 0 2 ϑµ , where µ 0 = 0.1, and for the prior distribution ϑ µ is a normally distributed random variable with zero mean and variance 5×10 -2 . The priors for v 1 and v 2 are assumed normal with means 0.0 and 10.0, and variances 0.05 and 1.0, respectively.

The results in terms of expected information gains are shown in Fig. 10. The source In the first column, the parameters, in the second and third column, the mean and the variance of the prior, in the last columns the mean value of the prior and the variances of the EEVs.

intensities are well identifiable. The less identifiable parameters are the diffusivity µ (which is quite small in the drift-dominated regime) and the horizontal velocity v 1 .

It should be noted that in Figure 10 the expected gain in information concerning the vertical component of the velocity v 2 increases before that related to the source intensities s 1 and s 2 . This is due to the fact that the front propagates from the bottom to the top, and, when it reaches the probes location, the change in the measurements can be first related to the advection speed and then to the intensity of the sources. The results are summarised in Table 5: the EEV analysis confirms that the source intensities and the drift-velocity v 2 are well identifiable.

Even if in both cases the source intensities are mathematically identifiable in a noiseless setting, there is a difference from a pragmatical standpoint. In the presence of noise, when different regimes of motion are considered, the information learned about the parameters can vary significantly. This is an important contribution of this kind of analysis.

Remark 7.2 (Limitations of the k-nearest neighbour mutual information estimator). While the estimator for mutual information and conditional mutual information works well for the above presented examples, it suffers from a drawback. In particular, in equation ( 26), if the joint probability distribution of Θ and Z 1:j lie on a lowerdimensional manifold, then the estimator yields erroneous results due to the incorrect assumption of constant uniform density in the -balls; see [START_REF] Kraskov | [END_REF] for this assumption in the formulation of the estimator and [36] for a discussion on this drawback. To ensure that such a drawback does not effect the above presented results, the results of the k-nearest neighbour based estimator are validated against those obtained by other methods, such as kernel density estimation.

Conclusion and Perspectives.

In this work a semi-empirical method to assess practical identifiability of parametric dynamical systems is proposed . It is set in a Bayesian framework and consists of analysing several information-theoretic measures. In particular, decrease in uncertainty of the system parameters due to the availability of noisy system measurements is computed via Shannon entropy. This decrease is viewed as the expected information gain and is related to identifiability of the parameters. Furthermore, the concept of a hypothetical measurement device that measures the parameters directly is proposed to facilitate interpretation of uncertainty reduction. In order to determine which subsets of parameters can only be identified collectively (i.e. individual parameters in this subset are unidentifiable), conditional mutual information is employed.

Four advantages of the proposed method are particularly relevant: first, that the method can take the measurement-noise of any structure into account; second, that it can be applied to generic dynamical systems including those governed by PDEs; third, that it is non-intrusive and does not require any modification of the existing numerical codes to solved the dynamical systems; and finally, that it has the ability to differentiate between parameter identifiability in different regions of the parametric space. The application of this approach is presented in a range of dynamical systems (including non-linear ODEs and PDEs) and, where possible, results are validated against those independently published in the literature.

There are two main drawbacks of the proposed method. First, that it is computationally intensive; however, this drawback is partly mitigated by the easily parallelisable nature of the method. Second, that the employed mutual information estimator presents problems if the joint distribution of the parameters and the observables lies in a low-dimensional manifold. To avoid the latter problem, future work Since Y p X,Y (x, y) dy = p X (x), the last two terms in the above equation cancel and hence

I H X|Y = X p X (x) log 1 p X (x) dx - Y X

6. 1 .

 1 k-nearest neighbour based estimator byKraskov et. al.

  is twice the distance between x i and its k th nearest neighbour, and c d is the volume of the N x -dimensional unit ball. The value of c d depends on the kind of norm used for calculating the distances between neighbours: for max norm c d = 1 and for Euclidean norm c d = π Nx/2

Figure 1 :

 1 Figure 1: The three-element Windkessel model and the flow-rate curve

) 7 . 1 . 1 . 0 = 9 .

 71109 Expected information gain: For the calculation of expected information gains, 150 measurements of z(t j ), for t j uniformly distributed in [0.0, 0.75] are considered. For the transformation the reference values are: Ψ P i 0 = 78.64 mmHg; Ψ Rp 0 = 0.838 mmHg-s/cm 3 ; Ψ C 0 = 0.0424 cm 3 /mmHg; and Ψ R d 109 mmHg-s/cm 3 . In the ψspace, Normal priors of zero mean and variance 1.0 are considered for P i , R p , C, and R d . The sample size is N = 3000 and the number of nearest neighbours is k = 10. Four successively increasing values of the noise variance, σ 2

Figure 2 :

 2 Figure 2: Gain in information for various parameter combinations in the Windkessel model.

Figure 3 :

 3 Figure 3: Mutual information between all pairs of the parameters in the Windkessel model. The plot for P i represents the mean of the observations generated through the propagated samples. σ 2 n corresponds to the transformed variable ψ P i (equation (65)). For σ 2 n = {0.0005, 0.005, 0.01, 0.05}, the corresponding variances of the noise process for the real parameter P i (untransformed, see equation (63)) are {1.8, 17.8, 35.7, 183.7} mmHg 2 .

Figure 4 :

 4 Figure 4: Comparison of the original Kraskov estimator and the modified estimator. The modified estimator handles the regions of constant entropy, t ∈ [0.45, 0.5], better than the original estimator which shows a decrease in such regions.

Figure 5 :d 1 , R d 2 |Z 1 M (R d 1 ;Figure 6 : 2 Parameter prior variance σ 2 e

 51211622 Figure 5: Schematic of a modified three-element Windkessel model

Figure 7 :

 7 Figure 7: Expected information gains in (nats) for the SIR model (section 7.2) against time in (s) for all the parameters.

Figure 8 :

 8 Figure 8: Information gains in (nats), for the 1D wave model (section 7.3) against time in (s) for all the harmonics of the potential.

Figure 9 :

 9 Figure 9: Setting for the PDE problem of advection-diffusion, see Sec.7.4. The domain is the unit square, the probes are located at y = a = 0.6. The capital letter D corresponds to Dirichlet boundary conditions while T are transparent boundary conditions. The velocity is mostly aligned with the ξ 2 axis. The function that mimics the source is characterised by two parameters S 1,2 . Superimposed, the contours of the solution at t = 5 10 -2 normalised, 40 values between the maximum and the minimum. The parameters to obtain it are: S 1 , S 2 = 25, v 1 = 0, v 2 = 10, µ = 0.1.

Figure 10 :Table 5 :

 105 Figure 10: Expected information gains in (nats), for the advection-diffusion equation in the drift-dominated regime of motion (section 7.4.1) against time in (s) for all the parameters.

  )) are {1.8, 17.8, 35.7, 183.7} mmHg 2 .

	Θ		Prior			Observation Noise Expected Posterior
	ψ-space		Real (Ψ-space)		ψ-space	Ψ-space
	µ ϑ σ 2 ϑ	Arith. mean	Arith. var.	Geom. mean	Geom. var.	σ 2 n	σ 2 n	σ 2 e	σ 2 u
	P 0.0 1.0 99.99 6167 78.64 1.62	-	-	-	-
						0.0005	1.8	0.976	0.976
	R p 0.0 1.0 1.065	0.7	0.838 1.62	0.005 0.01	17.8 35.7	1.066 1.285	1.066 1.285
						0.05	183.7	3.426	3.426
						0.0005	1.8	0.283	0.283
	C 0.0 1.0 0.054 0.0018 0.0424 1.62	0.005 0.01	17.8 35.7	0.392 0.557	0.392 0.557
						0.05	183.7	1.747	1.747
						0.0005	1.8	0.100	0.101
	R d 0.0 1.0 11.58 82.7 9.109 1.62	0.005 0.01	17.8 35.7	0.188 0.257	0.188 0.257
						0.05	183.7	0.912	0.912
	Table

Table 3 :

 3 Results for the SIR model (seeSection 7.

	ϑ ϑ (prior mean) σ 2 ϑ (prior variance)	σ 2 e	σ 2 u
	µ	2.0	0.4	0.101 0.265
	ϑ N	0	0.05	8.50 120.07
	η	1.75	0.25	0.054 0.213
	ω	0.6	0.1	0.14	1.24
	κ	0.4	0.1	0.801 20.82
	ϑ α	0	0.05	2.09 27.68
	ϑ γ	0	0.05	1.05 14.75

Table 4 .

 4 The EEV values clearly indicate that the identification is more and more difficult as the space frequency increases.

	ϑ ϑ (prior mean) σ 2 ϑ (prior variance)	σ 2 e	σ 2 u
	v1	0	1.0	0.094 0.094
	v2	0	1.0	0.247 0.247
	v3	0	1.0	0.417 0.417
	v4	0	1.0	0.967 0.967

Table 4 :

 4 Results for the Wave model (seeSection 7.

includes assessment and development of non-parametric entropy estimation methods in manifolds (see for example [37]).

). In the first column, the parameters, in the second and third column, the mean and the variance of the prior, in the last columns the mean value of the prior and the variances of the EEVs.

Diffusive regime:

In this section a diffusive regime, where the passive scalar dynamics is dominated by diffusion, is considered. All the parameters are kept constant as in the previous section except those for diffusivity and drift velocity. For µ the parameterisation µ = µ 0 2 ϑµ is considered, where µ 0 = 1.5 and ϑ µ is a normally distributed random variable with zero mean and variance 10 -2 . For v 2 a normal distribution with mean 7.5 and variance 1.0 is considered. The expected information gains for these settings are shown in Fig. 11. The sources are still individually identifiable. This is due to the fact that the drift velocity is still high. However, the gain in information about the source intensities is much less when compared to the drift-dominated case (0.4 nats over 1.4). These results are confirmed by the values reported in Table 6. The values of the EEVs are about 3 times larger.

Appendix A. Different measures of uncertainty lead to the same expected information gain

This appendix presents details of the argument that irrespective of which measure of inherent uncertainty is used, the expected gain in information when the probability distribution changes from p X (x) tp p X|Y (x|y) remains the same.

Appendix A.1. Shannon uncertainty

If equation ( 2) is seen as an inherent measure of uncertainty then the gain in information G X|Y =y (y) from equation ( 8) is

Consequently the expected gain in information, I S X|Y , from equation ( 9) is

p Y (y) dx dy. (A.2)

where p X|Y (x, y) represents the joint distribution of X and Y .

Appendix A.2. Jaynes' uncertainty

If equation ( 3) is seen as an inherent measure of uncertainty then the gain in information G X|Y =y (y) from equation ( 8) is

Since m(x) is an invariant measure of the random space and, it remains constant. Consequently, I J X|Y , can be written from equation ( 9) as

p Y (y) dx dy (A.5)

Since Y p X,Y (x, y) dy = p X (x), the last two terms in the above equation cancel and hence

According to the interpretation of the Kullback-Leibler distance that D(p X (x)||q X (x)) represents the gain in information when the probability distribution changes from q X (x) to p X (x). Hence, G X|Y =y (y) from equation ( 4) is

Consequently, I K X|Y can be written from equation [START_REF] Kaipio | Statistical and Computational inverse problems[END_REF] as If equation ( 5) is seen as an inherent measure of uncertainty then the gain in information G X|Y =y (y) from equation ( 8) is

Consequently, I H X|Y can be written from equation [START_REF] Kaipio | Statistical and Computational inverse problems[END_REF] as