Improvement of the k-nn Entropy Estimator with Applications in Systems Biology

Abstract

In this paper, we investigate efficient estimation of differential entropy for multivariate random variables. We propose bias correction for the nearest neighbor estimator, which yields more accurate results in higher dimensions. In order to demonstrate the accuracy of the improvement, we calculated the corrected estimator for several families of random variables. For multivariate distributions, we considered the case of independent marginals and the dependence structure between the marginal distributions described by Gaussian copula. The presented solution may be particularly useful for high dimensional data, like those analyzed in the systems biology field. To illustrate such an application, we exploit differential entropy to define the robustness of biochemical kinetic models.

1. Introduction

Determining an entropy of random variable has been broadly studied since the nineteenth century. Originally, the entropy concept crystallized in thermodynamics and quantum mechanics [1,2]. In information theory, the entropy was introduced by Claude Shannon [3] in 1948 for random variables with a discrete probability space. A natural extension of a discrete entropy is a differential entropy [4] defined for a continuous random variable with probability density function (pdf) $p\left(x\right)$ by:

$$H\left(X\right):=\mathbb{E}[-logp\left(x\right)]={\int }_X-log\left(p\left(x\right)\right)p\left(x\right)dx$$

It is worth mentioning that both notions, although analogous, have different properties, e.g., differential entropy is not necessarily non-negative.

The notion of entropy has many applications in various research fields. As an example, we mention here the emerging field of systems biology. It combines knowledge from various disciplines, including molecular biology, engineering, mathematics and physics, to model and analyze cellular processes. One of the challenging tasks in systems biology is to study interactions of regulatory, signaling and metabolism processes by means of biochemical kinetic models [5]. To check the robustness of such models, adequate sensitivity analysis procedures have to be applied.

Efficient entropy estimation is crucial for sensitivity analysis methods based on the mutual information of random variable measurements. An interesting approach for multi-variables system has been proposed recently by Lüdtke et al. [6]. However, it is based on a discrete entropy estimator and requires a computationally-inefficient variable discretization procedure. Having this concerns in mind, we propose here a more efficient k-nn differential entropy estimator, which can be used to calculate sensitivity indices in lieu of its discrete equivalent. Its application to sensitivity analysis for the model of the p53-Mdm2 feedback loop is presented in Section 5. In our approach, there is no need for variable discretization, and more importantly, an improved entropy estimator can be successfully used for highly multidimensional variables.

The estimation of differential entropy from the sample is an important problem with many applications; consequently, many methods have been proposed. Furthermore, several estimators based on the k-nn approach were analyzed, and their mathematical properties were investigated; see the discussion in the next section. We would like to emphasize that our main motivation here was very practical. Therefore, we focused on efficient implementation and bias correction for the well-established k-nn estimator with already proven consistency. The theory justifying our bias correction algorithm is based on the consistency analysis by Kozachnko and Leonenko [7]. Moreover, we performed extensive simulations for random samples from many multidimensional distributions recognized as most adequate in systems biology models.

Our improved k-nn entropy estimator explores the idea of correcting the density function evaluation near the boundary of random variable support. This smart idea has been proposed previously in [8] for several distributions. However, the novelty and usefulness of our approach lies in efficient bias correction of the classical k-nn estimator, which is applicable for any data sample (i.e., the a priori knowledge of the underlying distribution is not required).

2. k-nn Estimators of Differential Entropy

There are several approaches to entropy estimation from data samples; some are based on random variable discretization by partitioning [9], while others aim to estimate the probability mass according to the pdf [8,10]. The most efficient methods for probability mass estimation are based on k-th nearest neighbor (k-nn) estimation. The first k-nn approach to entropy estimation is dated to 1987 with the entropy estimator provided by Kozachenko and Leonenko [7] together with the proof of convergence. In this section, we recall definition and main properties of the k-nn entropy estimator by Kozachenko and Leonenko and we compare it to other approaches proposed in the literature.

Denote by $x_i$ for $i=1,\dots n$ i.i.d. samples of continuous random variable $X\in {\mathbb{R}}^d$. Let $r_d\left(x_i\right)$ be the d-dimensional distance of point $x_i$ to its nearest neighbor in the sample and $V_d$ be the volume of a d-dimensional unit ball. The estimator is defined for d-dimensional distributions as follows [7]:

$$\widehat{H}\left(X\right):=\frac{1}{n}\sum _{i=1}^n[-log\widehat{p}\left(x_i\right)]+\gamma$$

(1)

$$\widehat{p}\left(x_i\right)={[(n-1)·r_d{\left(x_i\right)}^d·V_d]}^{-1}$$

white where $\gamma =exp\left(-{\int }_0^{\infty }{e}^{-t}logtdt\right)\simeq 0.5772\dots$ is the Euler–Mascheroni constant.

For completeness, we recall here the further study conducted by Kraskov et al. [11], who proposed a variant of the entropy estimator for the k-th nearest neighbor:

$${\widehat{H}}_k\left(X\right):=-\psi \left(k\right)+\psi \left(n\right)+log\left(V_d\right)+\frac{d}{n}\sum _{i=1}^{n}log{r}_d\left(x_i\right)$$

(2)

where $r_d\left(x_i\right)$ is the d-dimensional distance of point $x_i$ to its k-th nearest neighbor in the sample and $\psi (.)$ is the digamma function with $\psi \left(1\right)=-\gamma$ and $\psi \left(n\right)\sim log(n-1)$.

Notice that for $k=1$, the Kraskov estimator given by Equation (2) performs almost equally to the Kozachenko–Leonenko estimator defined in Equation (1) with only a small unnoticeable difference stemming from: $\psi \left(n\right)-log(n-1)\simeq 0.$

According to analyses presented by Kozachenko et al. (1987) and Kraskov et al. (2004) in [7,11], k-nn estimators defined by Equations (1) and (2) converge in probability towards the real entropy value with increasing sample size:

$${\widehat{H}}_k\left(X\right)\stackrel{P}{\to }H\left(X\right)$$

Both estimators were shown to be consistent, i.e., asymptotically unbiased with the variance vanishing as the sample size increases. Moreover, for one-dimensional random variables, the estimator can be considered unbiased, even for a relatively small sample size, whereas the increase of dimension results in significant estimator bias even for a moderate sample size; see the exemplary case in Figure 1. The proof presented in [7] of the Kozachenko–Leonenko estimator in Equation (1) requires the following condition on the pdf $p\left(x\right)$ for the asymptotic unbiasedness of the estimator:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathbb{E}|logp\left(x\right)|}^{1+ϵ}<\infty ,\\ \hfill {\mathbb{E}|log\rho (x,y)|}^{1+ϵ}<\infty \end{array}\end{array}$$

(3)

and for the estimator consistency:

$$\begin{array}{c}\hfill {\mathbb{E}|logp\left(x\right)|}^{2+ϵ}<\infty ,\\ \hfill {\mathbb{E}|log\rho (x,y)|}^{2+ϵ}<\infty \end{array}$$

(4)

where $\rho (x,y)$ denote a metric in ${\mathbb{R}}^d$ and the expectations are over independent variables x and y.

Figure 1. Performance of the k-nn estimator Equation (1) for increasing sample size: box-plots of estimated entropy values for samples from a uniform distribution: on interval [0,1] (top); on hypercube ${[0,1]}^5$ (center); on hypercube ${[0,1]}^{15}$ (bottom); with real entropy value denoted by the red line.

Figure 1. Performance of the k-nn estimator Equation (1) for increasing sample size: box-plots of estimated entropy values for samples from a uniform distribution: on interval [0,1] (top); on hypercube ${[0,1]}^5$ (center); on hypercube ${[0,1]}^{15}$ (bottom); with real entropy value denoted by the red line.

Pérez-Cruz (2008) in [12] requires a somehow weaker assumption that the pdf must exist, i.e., the probability measure must be absolutely continuous. Consequently, the almost sure convergence of the entropy estimator ${\widehat{H}}_k\left(X\right)\stackrel{a.s.}{\to }H\left(X\right)$ is stated by the author.

The k-nn estimator for the related notion of Shannon entropy was proposed by Leonenko et al. (2008) in [13]. It can be shown that it is equivalent to those defined by Equation (2) with component $\psi \left(n\right)$ replaced by $log(n-1)$. Similarly, the Shannon entropy estimator presented by Penrose et al. (2013) in [14] can be expressed by Equation (1) with the component $log(n-1)$ replaced by $log\left(n\right)$. Both estimators were proven to be convergent to the real entropy value in ${\mathbb{L}}^2$ with some restrictions for density $p\left(x\right)$. Namely, in [13], the pdf must be bounded, and $\mathbb{E}\left\{p^{q-1}\left(x\right)\right\}={\int }_{{\mathbb{R}}^d}{p}^q\left(x\right)dx<\infty$ for some $q<1$. The boundness condition has been relaxed in [15] to an unbounded pdf with additional constraint $sup\{r>0:\mathbb{E}\parallel X{\parallel }^r<\infty \}>0$; white while in [14], the pdf must be bounded or the pdf must be defined on ${\mathbb{R}}^d$ with $sup\{r⩾0:\mathbb{E}\parallel X{\parallel }^r<\infty \}>0.$

It is worth mentioning that also, estimators for the Kullback–Leibler divergence were proposed and analyzed by Pérez-Cruz and Wang in [12,16]. Nevertheless, these estimators were based on the same idea of probability mass estimation by the distance to the k-nn in a sample; consequently, when used for estimating mutual information (cf. Definition 2), they suffer from the same problems as the above-mentioned entropy estimators.

A slightly different concept based on order statistics generalized to Voronoi regions in multiple dimensions was explored by Miller [17] to estimate the differential entropy and applied in [18] for independent component analysis. Unfortunately, the complexity of calculating the volume of Voronoi regions appears to be exponential in the dimension, and therefore, the method does not seem to be promising for our applications.

3. Bias Correction of the k-nn Entropy Estimator

In higher dimensions, k-nn estimators for many important distributions converge slowly and with noticeable bias; cf. Figure 1. Moreover, the entropy value is often overestimated; see Figure 1 for a d-dimensional uniform distribution and Appendix Figure A1 for a d-dimensional distribution with independent marginals sampled from an exponential distribution. On the other hand, for some distributions, like multivariate Gaussian ones, the estimator performs well, even in higher dimensions; cf. Appendix Figure A1 (top panel). As argued in [8], the overestimation problem affects random variables with bounded support. Before giving the description of our bias correction procedure, we briefly present the density mass estimation problem.

Recall that the k-nn estimation for continuous random variables requires probability mass estimation in small neighborhoods of sampled points. Consequently, when the sample points are located densely near the boundary of a given distribution, the calculations might be overestimated, as the samples’ neighborhoods might exceed the boundary of random variable support.

Therefore, the entropy estimator bias stems from overestimation of the pdf that occurs when the sample point is located near the boundaries of a random variable support. The curse of dimensionality implies that the boundary area increases with the variable dimension. Thus, in higher dimensions, neighborhoods of sampled points more frequently exceed the boundaries of a variable support.

Here, we focus on correcting the bias of the k-nn estimator. We restrict our further analysis to the Kozachenko and Leonenko estimator in Equation (1), but our procedure may be easily adapted to reduce the bias of the estimator defined by Equation (2) and other estimators discussed in the previous section. The precise reason for the biased performance of the k-nn entropy estimator may be exposed by analyzing the proof of its convergence [7].

Let ${\left\{x_i\right\}}_{i=1\dots n}$ be a set of independent random samples of a variable $X\in {\mathbb{R}}^d$ with probability density function $p\left(x\right)$. Let $\rho (x,y)$ be a metric in ${\mathbb{R}}^d$ with a neighborhood of radius r centered at point x defined by:

$$\nu (x,r)=\{y\in {\mathbb{R}}^d:\rho (x,y)<r\}$$

Then, the volume of the neighborhood is equal to $\left|\nu (x,r)\right|=r^d·V_d$, where $V_d$ is the unit ball volume according to the metric ρ.

For differential entropy estimation, the probability mass of a random variable X in a point $x_i$ is estimated using the following Lebesgue theorem with a small volume neighborhoods centered at sampled points $x_i$.

Theorem 1 (Lebesgue). Let $p\left(x\right)\in L_1\left({\mathbb{R}}^d\right)$; then, for almost all $x\in {\mathbb{R}}^d$ and open balls with radius $r_n\to 0$:

$$\underset{n\to \infty }{lim}\frac{1}{\left|\nu \right(x,r_n\left)\right|}{\int }_{\nu (x,r_n)}p\left(y\right)dy=p\left(x\right)$$

By analyzing the proof of convergence of the k-nn estimator (see the Appendix), we noticed that the application of Lebesgue’s theorem for distributions with support $supp\left(X\right)\ne {\mathbb{R}}^d$ can introduce a systematic error to the estimated value. The formal analysis reduces to a quotation of the main steps of the proof of convergence. The key point of the proof in which the systematic error is to be noticed is presented in Equation (B1) in the Appendix. The bias stems from the relatively slow convergence with the sample size to the exponent of the averaged pdf within neighborhoods of vanishing volumes that can exceed variable support.

To illustrate this issue, let us consider the exemplary case of the uniform distribution X on d-dimensional unit cube ${[0,1]}^d$. Every point in a cube is equally probable and the pdf $p\left(x\right)\equiv 1$ for points $x\in {[0,1]}^d$ and $p\left(x\right)\equiv 0$ for points $x\notin {[0,1]}^d$ outside the the cube. Lebesgue’s theorem applied for X states:

$$\underset{n\to \infty }{lim}\frac{1}{\left|\nu \right(x,r_n\left)\right|}{\int }_{\nu (x,r_n)}{\mathbb{1}}_{{[0,1]}^d}\left(y\right)dy=1$$

and it introduces the following estimation error for points x near the boundaries, i.e., with neighborhood $\nu (x,r_n)$ partially outside the cube ${[0,1]}^d$:

$$-log\left(\frac{p\left(x\right)}{\widehat{p}\left(x\right)}\right)\simeq -log\left(\frac{|\nu (x,r_n)\cap {[0,1]}^d|}{\left|\nu \right(x,r_n\left)\right|}\right)$$

where $r_n\left(x\right)=\sqrt[d]{\frac{r_d\left(x\right)}{V_d\gamma (n-1)}}$. Our correction to the entropy estimator is based on the logarithm of the proportions of neighborhood volumes $\nu (x_i,r_n)$ inside the support of X to the volume of the whole neighborhood $\nu (x_i,r_n)$:

$$\mathrm{CORRECTION}\left({\left\{x_i\right\}}_{i=1\dots n}\right)=\frac{1}{n}\sum _{i=1}^{n}log\left(\frac{|\nu (x_i,r_n)\cap supp\left(X\right)|}{\left|\nu \right(x_i,r_n\left)\right|}\right)$$

(5)

Finally, for $\widehat{H}\left(X\right)$ defined by Equation (1), the corrected entropy estimator is defined as follows:

$${\widehat{H}}_{corr}\left(X\right)=\widehat{H}\left(X\right)+\mathrm{CORRECTION}\left({\left\{x_i\right\}}_{i=1\dots n}\right)$$

(6)

To ensure the efficient implementation of our approach, we use the maximum metric in ${\mathbb{R}}^d$:

$$\rho (\mathbf{x},\mathbf{y})=max\left\{\right|x\left(1\right)-y\left(1\right)|,\dots ,|x\left(d\right)-y\left(d\right)\left|\right\}$$

for points $\mathbf{x}=\left[x\right(1\left)\dots x\right(d\left)\right]$ and $\mathbf{y}=\left[y\right(1\left)\dots y\right(d\left)\right]$. In this metric, the neighborhoods are d-dimensional hypercubes, and the unit ball volume equals $V_d=2^d$.

The support boundaries are estimated based on the sample and are determined simply as the leftmost and the rightmost values for each dimension. Determined boundary values (in an amount of at most $2·d$) are removed from the sample in the preprocessing phase, before entropy estimation.

Algorithm 1 presents the pseudocode of the corrected k-nn entropy estimator. The bias correction presented on the uniform example can be applied for any multivariate distribution. However, due to the source of bias, our correction is much more efficient for random variables with bounded support. Therefore, in the next section, we focus on the usefulness of Algorithm 1 for multidimensional uniform and beta distributions with both independent and dependent marginals. The exemplary dependency structure is modeled by a Gaussian copula.

The issue of estimator convergence for a one-dimensional variable has been analyzed by van der Meulen and Tsybakov [19], and the bias was proven to be of order $\mathcal{O}\left(\frac{1}{\sqrt{n}}\right)$, where n is the sample size.

Furthermore, the bounded support of a random variable as a source of k-nn estimator bias was already identified in papers by Sricharan et al. [8,20]. The authors estimated the bias using Taylor expansion of pdf $p\left(x\right)$, as:

$$\mathbb{E}\left[\widehat{H}\left(X\right)\right]-H\left(X\right)=\mathbb{E}{\left[p^{-(d+2)/d}\left(X\right)c\left(X\right)\right]}^{d/2}+1/2k+o(1/k+{(k/n)}^{2/d})$$

(7)

where $c\left(X\right)={\Gamma }^{2/d}\left(\frac{d+2}{2}\right)tr\left[{▽}^{2}p\left(x\right)\right]$, under the assumption that p is two-times continuously differentiable and is bounded away from zero. Compare also further work by Sricharan et al. [21,22].

It should be emphasized that the usefulness of analytical bias estimations is limited to variables with known pdf $p\left(x\right)$. In contrast to previous approaches, we did not aim at proving asymptotic bounds for differential entropy estimator bias. Instead, we propose the efficient heuristic procedure for bias correction. Moreover, in the applications, e.g., in systems biology, we have no a priori knowledge about the probability distributions of a random sample. The analytical bounds as in Equation (7) is of limited value without knowledge of pdf. To compare our bias correction procedure (Algorithm 1) to the estimation given by Equation (7), we calculated both corrections for a d-dimensional uniform distribution (see Figure 2). It occurred that compared to our method , the Equation (7) significantly underestimates the bias for higher dimensions.

Figure 2. Comparison of the k-nn entropy bias estimate by Sricharan et al. [8] (blue line) with the real bias based on the k-nn entropy estimation from sampled points (green line) and the k-nn entropy bias estimate obtained by our method (red line).

Figure 2. Comparison of the k-nn entropy bias estimate by Sricharan et al. [8] (blue line) with the real bias based on the k-nn entropy estimation from sampled points (green line) and the k-nn entropy bias estimate obtained by our method (red line).

4. k-nn Estimator Performance for Different Distributions

As argued above, the heavy bias for the k-nn differential entropy estimator characterizes multidimensional distributions supported on a bounded interval. If the sampled points appear frequently near interval boundaries, the probability mass determined with the use of the Lebesgue’s theorem is overestimated, due to the fact that the radius of an estimated neighborhood exceeds the distribution support. The reason for bias of sample based k-nn density estimators can be found in [8].

For the highly multivariate variables, improved estimator defined in Equation (6) performed significantly better with much smaller bias than the original k-nn estimator in Equation (1). Here, we demonstrate the efficiency of the proposed correction for two families of marginal distributions defined on the [0,1] interval: uniform and Beta. To study the estimator behavior, we sampled random points from multivariate distributions with both correlated and independent marginals.

4.1. Independent Marginals Case

Novel bias correction yields spectacular improvement of estimator accuracy for a multidimensional uniform distribution on unit hypercube (i.e., a multivariate distribution with independent uniform marginals on interval [0,1]). Analogous outcomes were achieved for multidimensional distributions with independent beta marginals (with supports on interval [0,1]); see Figure 3. Notice that the application of our bias correction is necessary, as comparable results could not be obtained by the increase of sample size; cf. Figure 4.

For completeness, we also tested the performance of estimator correction for multidimensional distributions with independent marginals that have unbounded support, e.g., Gaussian distribution, and one side bounded support, e.g., exponential distribution; see Appendix Figure A1. For a multivariate Gaussian distribution, the correction slightly underestimates the real value, but not significantly, especially for a bigger sample size. While for a multivariate distribution with independent exponential marginals, the corrected estimator introduces a small improvement; cf. Figure A2 in the Appendix.

Figure 3. Bias of the entropy estimator for growing dimensions for original k-nn entropy estimator in Equation (1) (red) and corrected entropy estimator in Equation (6) (blue) for multivariate random variables with independent marginals sampled from a uniform distribution on interval [0,1] (left); and from the Beta(3,1) distribution (right); for two different sample sizes (top and bottom).

Figure 3. Bias of the entropy estimator for growing dimensions for original k-nn entropy estimator in Equation (1) (red) and corrected entropy estimator in Equation (6) (blue) for multivariate random variables with independent marginals sampled from a uniform distribution on interval [0,1] (left); and from the Beta(3,1) distribution (right); for two different sample sizes (top and bottom).

Figure 4. Box plots for the growing sample size of estimated entropy values with k-nn entropy estimator in Equation (1) (top); and corrected entropy estimator in Equation (6) (center); for four-dimensional random variables with independent marginals sampled from a uniform distribution on interval [0,1] (left); and the Beta(3,1) distribution (right); w.r.t. the real entropy value denoted by the red line. The bottom panels demonstrate the histograms of marginal distributions.

Figure 4. Box plots for the growing sample size of estimated entropy values with k-nn entropy estimator in Equation (1) (top); and corrected entropy estimator in Equation (6) (center); for four-dimensional random variables with independent marginals sampled from a uniform distribution on interval [0,1] (left); and the Beta(3,1) distribution (right); w.r.t. the real entropy value denoted by the red line. The bottom panels demonstrate the histograms of marginal distributions.

4.2. Dependent Marginals Case

We aim here to test the accuracy of the corrected k-nn estimator for samples from multidimensional distributions with dependent marginals. The exemplary structure of dependency is provided by the Gaussian copula:

$$C_{\Sigma }={\mathsf{\Phi }}_{\Sigma }\left({\mathsf{\Phi }}^{-1}\left(u_1\right),\dots ,{\mathsf{\Phi }}^{-1}\left(u_d\right)\right)$$

where ${\mathsf{\Phi }}_{\Sigma }$ is a joint cumulative distribution function of the multivariate Gaussian distribution with mean vector zero and the covariance matrix equal to the correlation matrix Σ, and ${\mathsf{\Phi }}^{-1}$ is the inverse cumulative distribution function of the standard normal distribution. See Figure A3 in the Appendix for a copula probability density function and a copula function in the two-dimensional case.

The dependent marginals case is much more realistic in many applications, e.g., for sensitivity analysis in systems biology. We decided that the dependency structure between marginal distributions by means of an appropriate copula is flexible enough to model various dependencies between model parameters. Moreover, the controlled dependence structure allows one to calculate the real estimated entropy value and to estimate the bias of an estimator [23]. It occurred that in the case of dependent variables, the accuracy gained by the corrected estimator is as significant as in the case of independent marginals. Improvement of the k-nn estimator was tested for samples from multivariate distributions with marginals supported on a bounded interval, i.e., uniformly distributed on interval [0,1] and non-uniformly distributed on interval [0,1] (beta marginals); cf. Figure 5 and Figure 6.

Figure 5. Bivariate distributions with a Gaussian copula dependence structure with correlation coefficient $\rho =0.5$ and marginals sampled from a uniform distribution on interval [0,1] (left); vs. marginals sampled from the Beta(3,1) distribution (right).

Figure 5. Bivariate distributions with a Gaussian copula dependence structure with correlation coefficient $\rho =0.5$ and marginals sampled from a uniform distribution on interval [0,1] (left); vs. marginals sampled from the Beta(3,1) distribution (right).

Figure 6. Bias of the entropy estimator for growing dimensions for original k-nn entropy estimator in Equation (1) (red) and corrected entropy estimator in Equation (6) (blue) for multivariate random variables with dependent marginals sampled from a uniform distribution on interval [0,1] (left); and the Beta(3,1) distribution (right); for two different sample sizes (top-bottom). The dependence structure is given by the Gaussian copula with correlation coefficients among marginals $\rho =0.5$.

Figure 6. Bias of the entropy estimator for growing dimensions for original k-nn entropy estimator in Equation (1) (red) and corrected entropy estimator in Equation (6) (blue) for multivariate random variables with dependent marginals sampled from a uniform distribution on interval [0,1] (left); and the Beta(3,1) distribution (right); for two different sample sizes (top-bottom). The dependence structure is given by the Gaussian copula with correlation coefficients among marginals $\rho =0.5$.

5. Sensitivity Indices Based on the k-nn Entropy Estimator

Efficient entropy estimation is crucial in many applications, especially in sensitivity analysis based on mutual information [6]. Particularly, in systems biology, aiming to formally describe complex natural phenomena, an in-depth sensitivity analysis is required to better understand the modeled system behavior. Although, the method proposed by Lüdtke et al. [6] provides an approach to sensitivity analysis for multi-variables system, it is based on entropy estimation requiring variable discretization and, consequently is computationally inefficient. Moreover, it may provide biased or inadequate results for continuous random variables.

Therefore, here, we propose to modify the approach from [6] by applying our corrected k-nn differential entropy estimation with no need for variable discretization. The main concept in sensitivity analysis, taken from information theory, is mutual information (MI) between random variables. Mutual information indicates the amount of shared information between two random variables, i.e., input variable (parameter or explanatory variable) and output variable (explained variable). When used for samples of parameters and samples of corresponding model response, estimated MI can indicate the sensitivity of a model to the parameters.

Definition 2. The mutual information between continuous random variables $X\sim p\left(x\right)$ and $Y\sim p\left(y\right)$ is defined by:

$$I(X;Y):={\mathbb{E}}_{X,Y}\left[log\frac{p(x,y)}{p\left(x\right)p\left(y\right)}\right]=H\left(Y\right)+H\left(X\right)-H(X,Y)$$

where $(X,Y)\sim p(x,y)$ is the joint distribution of X and Y and $p\left(x\right)$, $p\left(y\right)$ and $p(x,y)$ are probability density functions.

Definition 3. The conditional entropy of a random variable $X\sim p\left(x\right)$ given random variable $Y\sim p\left(y\right)$ is defined as:

$$H\left(X\right|Y):={\mathbb{E}}_{X,Y}[-logp\left(x\right|y)]={\int }_{X,Y}-log\left(p\left(x\right|y)\right)p(x,y)dxdy$$

where $p(x,y)$ and $p\left(x\right|y)$ are respectively the joint and the conditional probabilities densities functions of variables X and Y.

There exist the following analogies between the information theoretic measures and set theory operations (see Figure B1 in the Appendix):

information theory	set theory
$H(X,Y)$	$X\cup Y$
$I(X;Y)$	$X\cap Y$
$H\left(X\right|Y)$	$X\setminus Y$

Inspired by [6], we redefine sensitivity indices for parameters (considered as continuous random variables) using the notion of differential entropy:

Definition 4. white Assume that $X_i$ are the parameters of the model and Y is the model output, then single sensitivity indices are defined as:

$$I(X_i;Y):=H\left(Y\right)+H\left(X_i\right)-H(X_i,Y)=H\left(Y\right)-H\left(Y\right|X_i)$$

white Analogously, group sensitivity indices for pairs of parameters are defined as:

$$I(X_i,X_j;Y):=H\left(Y\right)+H(X_i,X_j)-H(Y,X_i,X_j)=H\left(Y\right)-H\left(Y\right|X_i,X_j)$$

Clearly, Definition 4 can be extended for any subset of parameters.

The sensitivity indices reflect the impact of parameters on the model output. It may happen that the group sensitivity index for a pair of parameters had a high value, indicating the significant influence of these parameters, while two single sensitivity indices for these two parameters had a low value. We interpret such a case as the opposite (negative) interaction between this pair of parameters. To quantify such interactions, we use the following notion:

Definition 5. If $X_i$ are parameters of the model and Y is model output, then interactions indices within a pair of parameters are defined by:

$$\begin{array}{cc}\hfill I(X_i;X_j;Y)& ={\mathbb{E}}_{X_i,X_j,Y}\left[-log\frac{p\left(x_i\right)p\left(x_j\right)p\left(y\right)p(x_i,x_j,y)}{p(x_i,x_j)p(x_i,y)p(x_j,y)}\right]\hfill \\ & =H\left(X_i\right)+H\left(X_j\right)+H\left(Y\right)-H(X_i,Y)-H(X_j,Y)-H(X_i,X_j)\hfill \\ & +H(X_i,X_j,Y)\hfill \end{array}$$

An analogous definition of interactions was introduced by McGill in [24] for discrete random variables. The following corollary expresses the interaction indices by the group sensitivity indices; see also the Venn diagrams depicted in Figure B1 in the Appendix.

Corollary 6. Interaction index within the pair of parameters $X_i$ and $X_j$ can be expressed as the sum of single sensitivity indices of parameter $X_i$ and parameter $X_j$ to the output variable Y minus the group sensitivity index for a pair of parameters:

$$I(X_i;X_j;Y)=I(X_i;Y)+I(X_j;Y)-I(X_i,X_j;Y)$$

5.1. Case Study: Model of the p53-Mdm2 Feedback Loop

To validate the new approach to the global sensitivity analysis based on the entropy measure, we tested the method on a well-known and widely-studied example of the negative feedback loop of the p53 protein and Mdm2 ligase.

Tumor suppressor p53 protein, also known as TP53 transcription protein 53, is a transcription factor determining the fate of a cell in the case of DNA damage; p53 indirectly, via activation of the transcription of the encoding p21 gene, can block the cell cycle to repair DNA or activate a process of programmed cell death, called apoptosis. The main regulator of the concentration of p53 protein is ligase Mdm2/Hdm2 (double minute 2 mouse/human double minute 2), which, through ubiquitination, leads to degradation of p53 in the proteasome. In more than half of the cases of human cancers, p53 is inactivated or absent, which allows the mutated tumor cells to replicate and determines their immortality. Consequently, this protein is under investigation due to its property of leading to self-destruction of cancer cells, which could be successfully used as therapy in many types of cancer [25]. The modelled system (Figure 7 left panel) can be formally described with the use of ordinary differential equations (ODE) as follows:

$$\begin{array}{cc}\hfill \frac{d{Y}_1}{dt}=& X_1-X_6{Y}_1-X_2{Y}_3\frac{Y_1}{Y_1+k}\hfill \\ \hfill \frac{d{Y}_2}{dt}=& X_3{Y}_1-X_4{Y}_2\hfill \\ \hfill \frac{d{Y}_3}{dt}=& X_4{Y}_2-X_5{Y}_3\hfill \end{array}$$

where the variables of the modeled system $Y_i$ correspond to the concentration of species p53, mRna Mdm2 and Mdm2, respectively. The initial values of the variables are given by species concentrations at starting time point $t_0$:

$Y_1\left(t_0\right)=0$	p53 protein;
$Y_2\left(t_0\right)=0.8$	Mdm2 ligase;
$Y_3\left(t_0\right)=0.1$	mRna Mdm2 precursor of Mdm2 ligase.

The parameters of the model denoted by $X_i$ correspond to reaction rates with the following base values $p_i=\mathbb{E}\left(X_i\right)$ and their biological interpretations:

$p_1=0.9$	p53 production;	$p_4=0.8$	Mdm2 transcription;
$p_2=1.7$	Mdm2-dependent p53 degradation;	$p_5=0.8$	Mdm2 degradation;
$p_3=1.1$	p53-dependent Mdm2 production;	$p_6=0$	independent p53 degradation;
$k=0.0001$	p53 threshold for degradation by Mdm2.

To investigate the sensitivity of the model to the parameter changes, the parameters ($X_i$) were perturbed around the base value $p_i$ according to Gamma distributions:

$$X_i\sim \frac{x^{\kappa -1}{e}^{-x/\theta }}{\Gamma \left(\kappa \right){\theta }^{\kappa }}$$

with shape parameters $\kappa =3$ and $\theta =p_i/\kappa$ and mean values equal to parameter base values $\mathbb{E}\left(X_i\right)=p_i$. The negative regulation of p53 protein by Mdm2 ligase results in an oscillatory behavior in time (see Figure 7 right panel). A parameter responsible for negative regulation is Mdm2-dependent p53 degradation rate denoted by $X_2$.

Figure 7. Graphical scheme of the modelled system (left); and the oscillatory behavior for the considered parameters values (right).

Figure 7. Graphical scheme of the modelled system (left); and the oscillatory behavior for the considered parameters values (right).

Different possible system trajectories resulting from the perturbation of parameters are depicted in Figure B2 in the Appendix. It is worth noting that oscillatory behavior is sensitive to specific parameter values.

The output of the model Y is a three-dimensional variable with components denoting the species concentration that are averaged with respect to time:

$$Y=(Y_1,Y_2,Y_3)$$

To interpret novel sensitivity indices based on differential entropy, we compared them to standard local sensitivity indices, which are partial derivatives of model output w.r.t. the parameters of interest. As illustrated by the heat maps in Figure 8, both methods indicated the model inflow and outflow parameters ($X_1$ and $X_6$, respectively) as crucial for the systems’ behavior. Interestingly, the standard method (LSA) did not capture the parameter responsible for the negative feedback loop, whereas the method based on mutual information showed high sensitivity to the parameter $X_2$ crucial for the oscillatory behavior of the model.

Figure 8. Sensitivity indices based on mutual information (MI) (left); local sensitivity analysis based on derivatives of variables w.r.t. parameters averaged in time (right).

Figure 8. Sensitivity indices based on mutual information (MI) (left); local sensitivity analysis based on derivatives of variables w.r.t. parameters averaged in time (right).

The global sensitivity analysis aims at identifying the influence of some group of parameters on the system behavior. There exist several well-established approaches to this task, like the Sobol method or multi-parameter sensitivity analysis; see, e.g., [26]. However, determining the nature of mutual interactions as compatible or opposite seems a novelty, which was made possible thanks to defining the interaction indices.

In our case study, we determined the pair of parameters $X_1$ (model inflow) and $X_2$ (negative feedback) to have the highest sensitivity index among all pairs of parameters. Moreover, these parameters play the opposite roles in the model, which was illustrated as a negative interaction index for this pair of parameters. Analogously, parameters $X_4$ and $X_5$ responsible for Mdm2 transcription and degradation have a negative interaction and high common impact on the model behavior; see Figure 9.

Figure 9. Sensitivity indices for pairs of parameters (left); interactions within pairs of parameters, respectively, for the model output (right).

Figure 9. Sensitivity indices for pairs of parameters (left); interactions within pairs of parameters, respectively, for the model output (right).

6. Conclusions

In this work, we study the issue of efficient estimation of the differential entropy. By analyzing the proof of convergence of the estimator proposed in [7], we designed the efficient bias correction procedure that can by applied to any data sample (no a priori knowledge of the probability density function is required).

The source of systematic bias that impairs all existing sample based k-nn estimators has been identified previously in [8]. The problem affects all distributions defined on a finite support, and due to the curse of dimensionality, the bias increases for multivariate distributions.

We have compared the bias correction proposed by Sricharan et al. [8] to our procedure, which has been exhaustively tested for many families of distributions. It turned out that the correction spectacularly proposed in this paper improves the accuracy of the estimator in tested multivariate distributions (we illustrated the accuracy gain for multivariate uniform and beta distributions with both independent and dependent marginals). It is worth mentioning that, in the case of the distribution with support bounded at one side (e.g., exponential), the proposed amendment is relatively small. The correction does not introduce improvement for a multivariate Gaussian distribution with unbounded support. However, for the latter family of distributions, the original estimator bias is not observed.

Summarizing, our results can be successfully used in many applications, the most prominent for us being the sensitivity analysis for kinetic biochemical models, where one often assumes parameters (e.g., reaction rates or initial conditions) originate from distributions with bounded support. The applicability of an efficient entropy estimator for sensitivity analysis has been presented on the model of the p53-Mdm2 feedback loop.