# A simple algorithm for global sensitivity analysis with Shapley effects

Takashi Goda  School of Engineering, University of Tokyo, 7-3-1 Hongo, Bunkyo-
ku, Tokyo 113-8656, Japan.

###### Abstract

Global sensitivity analysis aims at measuring the relative importance of
different variables or groups of variables for the variability of a quantity
of interest. Among several sensitivity indices, so-called Shapley effects have
recently gained popularity mainly because the Shapley effects for all the
individual variables are summed up to the overall variance, which gives a
better interpretability than the classical sensitivity indices called main
effects and total effects. In this paper, assuming that all the input
variables are independent, we introduce a quite simple Monte Carlo algorithm
to estimate the Shapley effects for all the individual variables
simultaneously, which drastically simplifies the existing algorithms proposed
in the literature. We present a short Matlab implementation of our algorithm
and show some numerical results. A possible extension to the case where the
input variables are dependent is also discussed.

###### keywords:

Global sensitivity analysis , Shapley effect , Monte Carlo methods

##  1 Introduction

Global sensitivity analysis provides an indispensable framework in measuring
the relative importance of different variables or groups of variables for the
variability of a quantity of interest [1, 2], and has long been considered one
of the central problems in reliability engineering [3, 4, 5]. In particular,
since his pioneering work by Sobol’ [6, 7], variance-based sensitivity
analysis has been studied intensively and applied to a variety of subjects in
science and engineering. Two classical but still major sensitivity indices are
called _main effect_ and _total effect_. The main effect, also called the
first-order effect, measures the variance explained only by an input variable
or a group of input variables, whereas the total effect is given by the
overall variance minus the variance explained only by the complement
variables. When looking at these effects for the individual variables, the sum
of the main effects is always less than or equal to the overall variance,
while the sum of the total effects is always larger than or equal to the
overall variance. In this way the classical sensitivity indices have a
difficulty in normalization, and may cause some trouble when judging whether
one input variable is more important than another.

Recently the connection between variance-based global sensitivity analysis and
Shapley value from game theory has been studied by Owen [8]. The resulting
sensitivity index defined for individual variables is called _Shapley effect_.
The Shapley effect takes its value between the main effect and the total
effect, and importantly, the sum of the Shapley effects over all the
individual variables is exactly equal to the overall variance. This
normalization property gives us a better interpretability in determining
whether one input variable is more important than another. Moreover, the
Shapley effect has been proven to remove the conceptual problem inherit to
both the main and total effects when the input variables are correlated [9,
10].

In this paper we mainly focus on the case where the input variables are
independent and study Monte Carlo algorithms to estimate the Shapley effects
efficiently. Building upon the previous works by Song et al. [9] and Broto et
al. [11], which themselves build upon the idea from [12], in conjunction with
the famous pick-freeze scheme [7, 13, 14, 15, 16, 17], we introduce a quite
simple Monte Carlo algorithm to estimate the Shapley effects for all the
individual variables simultaneously at the cost of $(d+1)N$, where $d$ denotes
the number of variables and $N$ denotes the sample size. Our algorithm offers
the following advantages:

  1. 1.

Each Shapley effect is estimated unbiasedly.

  2. 2.

The variance of the estimator is also estimated unbiasedly and decays at the
canonical $1/N$ rate under some assumption.

  3. 3.

The above advantages lead to an approximate confidence interval for the
Shapley effect without requiring independent Monte Carlo trials or bootstrap
resampling.

  4. 4.

The sum of the estimates for all the Shapley effects is an unbiased estimator
of the overall variance.

We also discuss how to extend our algorithm to the case where the input
variables are dependent. Although we are particularly interested in Monte
Carlo sampling-based approaches in this paper, note that there are numerous
different approaches to carry out a global sensitivity analysis proposed in
the literature, such as random balance designs [18], polynomial chaos
expansion [19], fast Fourier transform [20], low-rank tensor approximation
[21] and random forests [22]. We refer the reader to [23] for a recent review
on various estimation algorithms for classical sensitivity indices.

The rest of this paper is organized as follows. In Section 2, we give an
overview of variance-based global sensitivity analysis and introduce several
sensitivity indices. In Section 3 we present our simple Monte Carlo algorithm
to estimate the Shapley effects for all the input variables simultaneously. We
also show some good properties of our algorithm as explained above. Numerical
experiments in Section 4 confirm the effectiveness of our algorithm for three
test cases. We conclude the paper in Section 5 with discussing two possible
extensions of our algorithm, one for an additional cost saving and the other
for the case where the input variables are dependent. In A, we provide a short
Matlab implementation of our algorithm, which we use for one test case.

##  2 Variance-based sensitivity analysis

###  2.1 ANOVA decomposition

Let
$\boldsymbol{x}=(x_{1},\ldots,x_{d})\in\Omega_{1}\times\cdots\times\Omega_{d}=:\Omega\subseteq\mathbb{R}^{d}$
be a vector of input random variables. Each variable $x_{j}$ follows a
probability distribution with density $\rho_{j}$ defined over the interval
$\Omega_{j}\subseteq\mathbb{R}$. Throughout this paper except Section 5.2, we
assume that all the random variables are independent with each other.
Moreover, for simplicity of notation, we write $[1:d]:=\\{1,\ldots,d\\}$. For
a subset $u\subseteq[1:d]$, we denote the complement of $u$ by
$-u=[1:d]\setminus u$ and denote the cardinality of $u$ by $|u|$. For vectors
$\boldsymbol{x},\boldsymbol{y}\in\Omega$, we write
$\boldsymbol{x}_{u}=(x_{j})_{j\in u}$ and
$(\boldsymbol{x}_{u},\boldsymbol{y}_{-u})=\boldsymbol{z}$ with $z_{j}=x_{j}$
if $j\in u$ and $z_{j}=y_{j}$ otherwise. The Cartesian product $\prod_{j\in
u}\Omega_{j}$ is denoted by $\Omega_{u}$, and the product of density functions
$\prod_{j\in u}\rho_{j}(x_{j})$ is denoted by $\rho_{u}(\boldsymbol{x}_{u})$
with an exception for $u=[1:d]$ in which case we simply write
$\rho(\boldsymbol{x})$.

Now let $f\colon\Omega\to\mathbb{R}$ be a function which outputs a quantity of
our interest. If the variance of $f$ with respect to $\boldsymbol{x}$ is
finite, $f$ can be decomposed as

$\displaystyle
f(\boldsymbol{x})=\sum_{u\subseteq[1:d]}f_{u}(\boldsymbol{x}_{u}),$ (1)

where each summand is recursively defined by

$f_{\emptyset}=\int_{\Omega}f(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}=:\mu$

and

$f_{u}(\boldsymbol{x}_{u})=\int_{\Omega_{-u}}f(\boldsymbol{x})\rho_{-u}(\boldsymbol{x}_{-u})\,\mathrm{d}\boldsymbol{x}_{-u}-\sum_{v\subsetneq
u}f_{v}(\boldsymbol{x}_{v})$

for any non-empty subset $u$. The following lemma shows important properties
of this decomposition.

######  Lemma 1.

With the notation above, the following holds true.

  1. 1.

For any non-empty subset $u$ and coordinate $j\in u$, we have

$\int_{\Omega_{j}}f_{u}(\boldsymbol{x}_{u})\rho_{j}(x_{j})\,\mathrm{d}x_{j}=0.$

  2. 2.

For any $u,v\subseteq[1:d]$, we have

$\int_{\Omega}f_{u}(\boldsymbol{x}_{u})f_{v}(\boldsymbol{x}_{v})\rho(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}=\begin{cases}\displaystyle\sigma_{u}^{2}:=\int_{\Omega_{u}}\left(f_{u}(\boldsymbol{x}_{u})\right)^{2}\rho_{u}(\boldsymbol{x}_{u})\,\mathrm{d}\boldsymbol{x}_{u}&\text{if
$u=v$,}\\\ 0&\text{otherwise.}\end{cases}$

###### Proof.

We refer the reader to [24, Appendix A.3] for the proof of the special case
where $\rho_{j}$ is the uniform distribution over the unit interval $[0,1]$,
which can be easily generalized to a proof of this lemma. ∎

It follows from the second assertion of Lemma 1 that the variance of $f$ can
be decomposed as

$\displaystyle\sigma^{2}$
$\displaystyle:=\int_{\Omega}\left(f(\boldsymbol{x})-\mu\right)^{2}\rho(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}=\int_{\Omega}\left(\sum_{\emptyset\neq
u\subseteq[1:d]}f_{u}(\boldsymbol{x}_{u})\right)^{2}\rho(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}$
$\displaystyle\>=\sum_{\emptyset\neq
u,v\subseteq[1:d]}\int_{\Omega}f_{u}(\boldsymbol{x}_{u})f_{v}(\boldsymbol{x}_{v})\rho(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}=\sum_{\emptyset\neq
u\subseteq[1:d]}\sigma_{u}^{2}.$

This way, the _overall variance_ $\sigma^{2}$ is decomposed into $2^{d}-1$
terms with each term $\sigma_{u}^{2}$ being the variance of a lower-
dimensional function $f_{u}$. This is why we call the decomposition (1) the
_analysis of variance (ANOVA) decomposition_ of $f$.

###  2.2 Sensitivity indices

For a non-empty subset $u\subseteq[1:d]$, the main effect and the total effect
for a group of input variables $\boldsymbol{x}_{u}$ are defined by

$\underline{\tau}^{2}_{u}:=\sum_{\emptyset\neq v\subseteq
u}\sigma_{v}^{2}\quad\text{and}\quad\overline{\tau}^{2}_{u}:=\sum_{\begin{subarray}{c}\emptyset\neq
v\subseteq[1:d]\\\ v\cap
u\neq\emptyset\end{subarray}}\sigma_{v}^{2}=\sigma^{2}-\underline{\tau}^{2}_{-u},$

respectively. We see from the definition that the main effect measures the
variance explained by the variables $\boldsymbol{x}_{u}$. On the other hand,
the total effect takes into account all possible interactions with the
variables $\boldsymbol{x}_{u}$, so that we always have
$\underline{\tau}^{2}_{u}\leq\overline{\tau}^{2}_{u}$. We see from the second
equality that the total effect measures the difference between the overall
variance and the variance explained by the complement variables
$\boldsymbol{x}_{-u}$.

The following identities on these effects are well-known:

$\underline{\tau}^{2}_{u}=\int_{\Omega}\int_{\Omega}f(\boldsymbol{x})f(\boldsymbol{x}_{u},\boldsymbol{y}_{-u})\rho(\boldsymbol{x})\rho(\boldsymbol{y})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{y}-\mu^{2},$

and

$\displaystyle\overline{\tau}^{2}_{u}=\frac{1}{2}\int_{\Omega}\int_{\Omega}\left(f(\boldsymbol{x})-f(\boldsymbol{y}_{u},\boldsymbol{x}_{-u})\right)^{2}\rho(\boldsymbol{x})\rho(\boldsymbol{y})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{y},$
(2)

see, e.g., [24, Appendix A.6]. These identities have been used for
constructing so-called _pick-freeze_ Monte Carlo estimators of the main and
total effects, respectively. We refer the reader to [7, 14, 15, 16] among many
others.

Looking at these effects for an individual variable, we have

$\sum_{j=1}^{d}\underline{\tau}^{2}_{\\{j\\}}=\sum_{j=1}^{d}\sigma^{2}_{\\{j\\}}\leq\sigma^{2},$

and

$\sum_{j=1}^{d}\overline{\tau}^{2}_{\\{j\\}}=\sum_{\emptyset\neq
u\subseteq[1:d]}|u|\sigma_{u}^{2}\geq\sigma^{2}.$

Therefore, the _rescaled_ versions of these effects,
$\underline{\tau}^{2}_{\\{j\\}}/\sigma^{2}$ and
$\overline{\tau}^{2}_{\\{j\\}}/\sigma^{2}$, are not summed up to 1. This
normalization issue can be circumvented by employing a new class of
sensitivity indices called Shapley effects [8]. For an individual variable
$x_{j}$, the Shapley effect is defined as follows:

$\displaystyle\phi_{j}:=\frac{1}{d}\sum_{u\subseteq-\\{j\\}}\binom{d-1}{|u|}^{-1}\left(\underline{\tau}^{2}_{u+j}-\underline{\tau}^{2}_{u}\right)=\frac{1}{d}\sum_{u\subseteq-\\{j\\}}\binom{d-1}{|u|}^{-1}\left(\overline{\tau}^{2}_{u+j}-\overline{\tau}^{2}_{u}\right),$
(3)

where we write
$\underline{\tau}^{2}_{\emptyset}=\overline{\tau}^{2}_{\emptyset}=0$. Here the
second equality is proven in [9, Theorem 1]. As pointed out in [8],

$\sum_{j=1}^{d}\phi_{j}=\sigma^{2}$

holds, so that the Shapley effect enables an easy interpretation when
measuring the relative importance of individual variables. Moreover, it is
known from [8, Theorem 1] that

$\phi_{j}=\sum_{\begin{subarray}{c}\emptyset\neq u\subseteq[1:d]\\\ j\in
u\end{subarray}}\frac{\sigma_{u}^{2}}{|u|},$

which leads to an inequality
$\underline{\tau}^{2}_{\\{j\\}}\leq\phi_{j}\leq\overline{\tau}^{2}_{\\{j\\}}$.

##  3 Monte Carlo estimator for Shapley effect

###  3.1 Algorithm

In this section, we introduce a simple unbiased Monte Carlo algorithm to
estimate the Shapley effects $\phi_{j}$ for all $1\leq j\leq d$. Although we
shall use the second expression of $\phi_{j}$ in (3), it is also possible to
construct a similar Monte Carlo algorithm based on the first expression. For a
start, let us consider a single value $\phi_{j}$. Let
$\ell\in\\{0,1,\ldots,d-1\\}$ be a uniformly distributed discrete random
variable, and given $\ell$, let $\mathcal{U}_{j}(\ell)\subseteq-\\{j\\}$ be a
uniformly distributed random subset with fixed cardinality $\ell$. Then it
follows from (3) and (2) that the Shapley effect for a variable $x_{j}$ is
given by

$\displaystyle\phi_{j}$
$\displaystyle=\frac{1}{d}\sum_{\ell=0}^{d-1}\binom{d-1}{\ell}^{-1}\sum_{\begin{subarray}{c}u\subseteq-\\{j\\}\\\
|u|=\ell\end{subarray}}\left(\overline{\tau}^{2}_{u+j}-\overline{\tau}^{2}_{u}\right)$
$\displaystyle=\mathbb{E}_{\ell}\mathbb{E}_{\mathcal{U}_{j}(\ell)}\left[\overline{\tau}^{2}_{\mathcal{U}_{j}(\ell)+j}-\overline{\tau}^{2}_{\mathcal{U}_{j}(\ell)}\right]$
$\displaystyle=\mathbb{E}_{\ell}\mathbb{E}_{\mathcal{U}_{j}(\ell)}\left[\int_{\Omega}\int_{\Omega}g_{\mathcal{U}_{j}(\ell)}(\boldsymbol{x},\boldsymbol{y})\rho(\boldsymbol{x})\rho(\boldsymbol{y})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{y}\right]$
$\displaystyle=\int_{\Omega}\int_{\Omega}\mathbb{E}_{\ell}\mathbb{E}_{\mathcal{U}_{j}(\ell)}\left[g_{\mathcal{U}_{j}(\ell)}(\boldsymbol{x},\boldsymbol{y})\right]\rho(\boldsymbol{x})\rho(\boldsymbol{y})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{y},$

wherein we have defined

$\displaystyle g_{\mathcal{U}_{j}(\ell)}(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=\frac{1}{2}\left(f(\boldsymbol{x})-f(\boldsymbol{y}_{\mathcal{U}_{j}(\ell)+j},\boldsymbol{x}_{-(\mathcal{U}_{j}(\ell)+j)})\right)^{2}-\frac{1}{2}\left(f(\boldsymbol{x})-f(\boldsymbol{y}_{\mathcal{U}_{j}(\ell)},\boldsymbol{x}_{-\mathcal{U}_{j}(\ell)})\right)^{2}$
$\displaystyle=\left(f(\boldsymbol{x})-\frac{f(\boldsymbol{y}_{\mathcal{U}_{j}(\ell)},\boldsymbol{x}_{-\mathcal{U}_{j}(\ell)})+f(\boldsymbol{y}_{\mathcal{U}_{j}(\ell)+j},\boldsymbol{x}_{-(\mathcal{U}_{j}(\ell)+j)})}{2}\right)$
$\displaystyle\quad\times\left(f(\boldsymbol{y}_{\mathcal{U}_{j}(\ell)},\boldsymbol{x}_{-\mathcal{U}_{j}(\ell)})-f(\boldsymbol{y}_{\mathcal{U}_{j}(\ell)+j},\boldsymbol{x}_{-(\mathcal{U}_{j}(\ell)+j)})\right).$

This representation naturally leads to the following unbiased Monte Carlo
estimator of $\phi_{j}$:

$\widehat{\phi}_{j,N}=\frac{1}{N}\sum_{n=1}^{N}g_{\mathcal{U}^{(n)}_{j}(\ell^{(n)})}(\boldsymbol{x}^{(n)},\boldsymbol{y}^{(n)}),$

where $\boldsymbol{x}^{(1)},\ldots,\boldsymbol{x}^{(N)}$ and
$\boldsymbol{y}^{(1)},\ldots,\boldsymbol{y}^{(N)}$ are i.i.d. samples
generated from the density $\rho$, $\ell^{(1)},\ldots,\ell^{(N)}$ are i.i.d.
samples of $\ell\in\\{0,1,\ldots,d-1\\}$, and for each $n$,
$\mathcal{U}^{(n)}_{j}(\ell^{(n)})$ denotes a random sample of
$\mathcal{U}_{j}(\ell^{(n)})$ conditional on $\ell^{(n)}$.

The key difference from the algorithms proposed in [9] and [11] is that we
directly apply the pick-freeze scheme to estimate the difference
$\overline{\tau}^{2}_{u+j}-\overline{\tau}^{2}_{u}$ instead of estimating two
terms independently. This idea not only gives a tight coupling in estimating
the difference, but also exploits a sensible recycling idea from [9] well
further. In fact, if each of $\phi_{1},\ldots,\phi_{d}$ is estimated
independently, we need $3dN$ total function evaluations since computing each
of
$g_{\mathcal{U}^{(n)}_{j}(\ell^{(n)})}(\boldsymbol{x}^{(n)},\boldsymbol{y}^{(n)})$
requires evaluating the function $f$ three times, whereas we can reduce this
cost to $(d+1)N$ through the recycling technique by estimating all
$\phi_{1},\ldots,\phi_{d}$ simultaneously.

Let $\boldsymbol{\pi}=(\pi(1),\ldots,\pi(d))$ denote a random permutation of
$[1:d]$. For a fixed $1\leq j\leq d$, generating $\ell\in\\{0,1,\ldots,d-1\\}$
and $\mathcal{U}_{j}(\ell)$ randomly is equivalent to finding $\ell$ such that
$\pi(\ell+1)=j$ and set $\mathcal{U}_{j}(\ell)=\\{\pi(1),\ldots,\pi(\ell)\\}$.
Note that we set $\mathcal{U}_{j}(0)$ to the empty set. Because of this
equivalence, by writing $\boldsymbol{\pi}_{j}=\\{\pi(1),\ldots,\pi(\ell)\\}$
with $\ell$ satisfying $\pi(\ell+1)=j$, we can rewrite our estimator as

$\widehat{\phi}_{j,N}=\frac{1}{N}\sum_{n=1}^{N}g_{\boldsymbol{\pi}^{(n)}_{j}}(\boldsymbol{x}^{(n)},\boldsymbol{y}^{(n)}),$

for i.i.d. permutations $\boldsymbol{\pi}^{(1)},\ldots,\boldsymbol{\pi}^{(N)}$
of $[1:d]$. The recycling technique from [9] shares the same
$\boldsymbol{\pi}^{(1)},\ldots,\boldsymbol{\pi}^{(N)}$ for computing all of
$\widehat{\phi}_{1,N},\ldots,\widehat{\phi}_{d,N}$, which still ensures the
unbiaseness of estimators. After some rearrangements, our algorithm can be
summarized as follows.

######  Algorithm 1 (Monte Carlo estimation of Shapley effects for all input
variables).

Let $d$ be the number of input variables and $N$ be the sample size.
Initialize $\widehat{\phi}_{1,N}=\cdots=\widehat{\phi}_{d,N}=0$. For $1\leq
n\leq N$, do the following:

  1. 1.

Generate $\boldsymbol{x}^{(n)},\boldsymbol{y}^{(n)}\in\Omega$ and
$\boldsymbol{\pi}^{(n)}$ randomly.

  2. 2.

Let $\ell=1$ and compute
$\mathcal{F}_{n}=\mathcal{F}^{-}_{n}=f(\boldsymbol{x}^{(n)})$.

  3. 3.

    1. (a)

Compute
$\mathcal{F}^{+}_{n}=f(\boldsymbol{y}^{(n)}_{\\{\pi(1),\ldots,\pi(\ell)\\}},\boldsymbol{x}^{(n)}_{-\\{\pi(1),\ldots,\pi(\ell)\\}})$.

    2. (b)

Update

$\widehat{\phi}_{\pi(\ell),N}=\widehat{\phi}_{\pi(\ell),N}+\frac{1}{N}\left(\mathcal{F}_{n}-\frac{\mathcal{F}^{-}_{n}+\mathcal{F}^{+}_{n}}{2}\right)\left(\mathcal{F}^{-}_{n}-\mathcal{F}^{+}_{n}\right).$

    3. (c)

Let $\ell=\ell+1$. If $\ell\leq d$, let
$\mathcal{F}^{-}_{n}=\mathcal{F}^{+}_{n}$ and go to Step (a).

It is obvious that we evaluate function values only $d+1$ times for each $n$,
leading to the total computational cost of $(d+1)N$.

######  Remark 1.

In Algorithm 1, instead of generating
$\boldsymbol{x}^{(n)},\boldsymbol{y}^{(n)}$ randomly for all $n$, one can use
Latin hypercube sampling or (randomized) quasi-Monte Carlo sampling. However,
improving the convergence of the estimated Shapley effects seems quite hard.
The reason is that we also need to generate a discrete object
$\boldsymbol{\pi}^{(n)}$ and the integrand
$g_{\boldsymbol{\pi}_{j}}(\boldsymbol{x},\boldsymbol{y})$ changes depending on
the samples $\boldsymbol{x}^{(n)}$ and $\boldsymbol{y}^{(n)}$, so that the
standard theory of quasi-Monte Carlo sampling which requires a smoothness of
integrands does not apply. This problem does not occur for the main and total
effects.

###  3.2 Some properties

Because of the linearity of expectation, the first property on the
unbiasedness of our Monte Carlo estimator, mentioned in Section 1, can be
easily shown as follows:

$\mathbb{E}\left[\widehat{\phi}_{j,N}\right]=\frac{1}{N}\sum_{n=1}^{N}\mathbb{E}\left[g_{\mathcal{U}^{(n)}_{j}(\ell^{(n)})}(\boldsymbol{x}^{(n)},\boldsymbol{y}^{(n)})\right]=\mathbb{E}_{\ell}\mathbb{E}_{\mathcal{U}_{j}(\ell)}\mathbb{E}_{\boldsymbol{x},\boldsymbol{y}\sim\rho}\left[g_{\mathcal{U}_{j}(\ell)}(\boldsymbol{x},\boldsymbol{y})\right]=\phi_{j}.$

This leads to the fourth property

$\mathbb{E}\left[\sum_{j=1}^{d}\widehat{\phi}_{j,N}\right]=\sum_{j=1}^{d}\mathbb{E}\left[\widehat{\phi}_{j,N}\right]=\sum_{j=1}^{d}\phi_{j}=\sigma^{2}.$

We now show the second property of the estimator. The third property
immediately follows from the central limit theorem. Let us make the second
property more explicit.

######  Theorem 1.

Assume that

$M_{4}(f):=\int_{\Omega}\left(f(\boldsymbol{x})\right)^{4}\rho(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}<\infty.$

Then the following holds true.

  1. 1.

The variance of $\widehat{\phi}_{j,N}$ is finite and decays at the rate of
$1/N$ for all $j$.

  2. 2.

The variance of $\widehat{\phi}_{j,N}$ is estimated unbiasedly by

$\frac{1}{N(N-1)}\sum_{n=1}^{N}\left(g_{\mathcal{U}^{(n)}_{j}(\ell^{(n)})}(\boldsymbol{x}^{(n)},\boldsymbol{y}^{(n)})-\widehat{\phi}_{j,N}\right)^{2}.$

###### Proof.

Since the second assertion is well-known [24, Chapter 2], we only give a proof
for the first assertion. Because of the independence between different
samples, we have

$\displaystyle\mathbb{V}\left[\widehat{\phi}_{j,N}\right]=\frac{\mathbb{V}\left[g_{\mathcal{U}_{j}(\ell)}(\boldsymbol{x},\boldsymbol{y})\right]}{N},$

where the variance on the right-hand side is taken with respect to
$\boldsymbol{x},\boldsymbol{y}\sim\rho$, $\ell\in\\{0,1,\ldots,d-1\\}$ and
$\mathcal{U}_{j}(\ell)\subseteq-\\{j\\}$. Thus it suffices to prove that
$\mathbb{V}\left[g_{\mathcal{U}_{j}(\ell)}(\boldsymbol{x},\boldsymbol{y})\right]$
is finite under the assumption $M_{4}(f)<\infty$. In fact, applying Jensen’s
inequality twice, we have

$\displaystyle\mathbb{V}\left[g_{\mathcal{U}_{j}(\ell)}(\boldsymbol{x},\boldsymbol{y})\right]\leq\mathbb{E}\left[\left(g_{\mathcal{U}_{j}(\ell)}(\boldsymbol{x},\boldsymbol{y})\right)^{2}\right]$
$\displaystyle=\mathbb{E}_{\ell}\mathbb{E}_{\mathcal{U}_{j}(\ell)}\left[\int_{\Omega}\int_{\Omega}\left(g_{\mathcal{U}_{j}(\ell)}(\boldsymbol{x},\boldsymbol{y})\right)^{2}\rho(\boldsymbol{x})\rho(\boldsymbol{y})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{y}\right]$
$\displaystyle=\frac{1}{d}\sum_{\ell=0}^{d-1}\binom{d-1}{\ell}^{-1}\sum_{\begin{subarray}{c}u\subseteq-\\{j\\}\\\
|u|=\ell\end{subarray}}\int_{\Omega}\int_{\Omega}\rho(\boldsymbol{x})\rho(\boldsymbol{y})$
$\displaystyle\quad\quad\times\left(\frac{1}{2}\left(f(\boldsymbol{x})-f(\boldsymbol{y}_{u+j},\boldsymbol{x}_{-(u+j)})\right)^{2}-\frac{1}{2}\left(f(\boldsymbol{x})-f(\boldsymbol{y}_{u},\boldsymbol{x}_{-u})\right)^{2}\right)^{2}\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{y}$
$\displaystyle\leq\frac{1}{d}\sum_{\ell=0}^{d-1}\binom{d-1}{\ell}^{-1}\sum_{\begin{subarray}{c}u\subseteq-\\{j\\}\\\
|u|=\ell\end{subarray}}\int_{\Omega}\int_{\Omega}\rho(\boldsymbol{x})\rho(\boldsymbol{y})$
$\displaystyle\quad\quad\times\left(\frac{1}{2}\left(f(\boldsymbol{x})-f(\boldsymbol{y}_{u+j},\boldsymbol{x}_{-(u+j)})\right)^{4}+\frac{1}{2}\left(f(\boldsymbol{x})-f(\boldsymbol{y}_{u},\boldsymbol{x}_{-u})\right)^{4}\right)\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{y}$
$\displaystyle\leq\frac{1}{d}\sum_{\ell=0}^{d-1}\binom{d-1}{\ell}^{-1}\sum_{\begin{subarray}{c}u\subseteq-\\{j\\}\\\
|u|=\ell\end{subarray}}\int_{\Omega}\int_{\Omega}\rho(\boldsymbol{x})\rho(\boldsymbol{y})$
$\displaystyle\quad\quad\times\left(8\left(f(\boldsymbol{x})\right)^{4}+4\left(f(\boldsymbol{y}_{u+j},\boldsymbol{x}_{-(u+j)})\right)^{4}+4\left(f(\boldsymbol{y}_{u},\boldsymbol{x}_{-u})\right)^{4}\right)\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{y}$
$\displaystyle=16M_{4}(f)<\infty.$

Thus we are done. ∎

##  4 Numerical experiments

We conduct some numerical experiments to confirm the effectiveness of our
Shapley effect estimator and compare the result with the shapleyPermRand
estimator proposed in [9], which is implemented in the sensitivity package of
the R software [25]. For the first test case, we consider Ishigami function
which is a function of only three variables but is interesting in terms of
interaction and non-linearity. For the second test case, we consider Sobol’
$g$ function which has been often used as a standard test problem in the
context of global sensitivity analysis in high dimensions [1]. The third
example is taken from [26, Section 7] which involves the assessment of
structural component strength. A short implementation of our algorithm in
Matlab for the second test case is available in A.

###  4.1 Ishigami function

Let us consider the case where $\rho$ is the uniform distribution over the
domain $\Omega=[-\pi,\pi]^{3}$. Ishigami function is defined by

$f(x_{1},x_{2},x_{3})=\left(1+bx_{3}^{4}\right)\sin x_{1}+a\left(\sin
x_{2}\right)^{2},$

with some coefficients $a,b>0$. For this function, we see that

$\displaystyle\underline{\tau}^{2}_{\\{1\\}}=\frac{1}{2}\left(1+\frac{\pi^{4}b}{5}\right)^{2},\quad\underline{\tau}^{2}_{\\{2\\}}=\frac{a^{2}}{8},\quad\underline{\tau}^{2}_{\\{3\\}}=0,$
$\displaystyle\overline{\tau}^{2}_{\\{1\\}}=\underline{\tau}^{2}_{\\{1\\}}+\frac{8\pi^{8}b^{2}}{225},\quad\overline{\tau}^{2}_{\\{2\\}}=\underline{\tau}^{2}_{\\{2\\}},\quad\overline{\tau}^{2}_{\\{3\\}}=\frac{8\pi^{8}b^{2}}{225},$
$\displaystyle\phi_{1}=\underline{\tau}^{2}_{\\{1\\}}+\frac{4\pi^{8}b^{2}}{225},\quad\phi_{2}=\underline{\tau}^{2}_{\\{2\\}},\quad\phi_{3}=\frac{4\pi^{8}b^{2}}{225}.$

This means that the second variable $x_{2}$ does not have any interaction with
the other variables and that the third variable $x_{3}$ affects the variation
of $f$ only through an interaction with the first variable $x_{1}$.

In what follows, we choose $a=7$ and $b=0.1$. The Shapley effects for all the
input variables are estimated according to Algorithm 1. As a comparison, they
are also estimated by the shapleyPermRand estimator with the choice of
$N_{v}=2N$, $m=\lceil N/3\rceil$, $N_{o}=1$ and $N_{i}=3$ such that the number
of the total function evaluations is almost equal to that of our proposed
algorithm, which is $(d+1)N$. Although we refer the reader to [25] for what
each input parameter means, the choice ($m=\lceil N/3\rceil$, $N_{o}=1$ and
$N_{i}=3$) is nearly optimal according to the theoretical analysis given in
[9, Claim 2].

Figure 1 shows the result for the case $N=2^{14}$. Not only the values of
Shapley effects estimated by the two algorithms but also the exact values of
the overall variance, main effects, Shapley effects and total effects are
plotted. Besides, the 95% confidence intervals for the Shapley effects are
shown by error bars. It is clear that our proposed algorithm provides better
estimates of the Shapley effects than the compared algorithm. In particular,
the compared algorithm underestimates the Shapley effect for the second
variable $x_{2}$ and the upper confidence limit is slightly larger than the
exact value. Moreover, the confidence intervals for our proposed algorithm are
narrower than those for the compared algorithm for all the input variables.

Figure 1: The estimated Shapley effects of individual variables for Ishigami
function with $N=2^{14}$. The results obtained by our proposed algorithm are
shown in yellow together with the confidence intervals, while those by the
shapleyPermRand estimator are in purple. The exact overall variance, main
effects, Shapley effects and total effects are also shown for comparison.

In order to evaluate the errors from the exact values quantitatively, we run
$R=10$ independent trials with a fixed $N$ and estimate the expected sum of
squared errors (SSE) by

$\displaystyle\mathbb{E}\left[\sum_{j=1}^{d}\left(\widehat{\phi}_{j,N}-\phi_{j}\right)^{2}\right]\approx\frac{1}{R}\sum_{r=1}^{R}\sum_{j=1}^{d}\left(\widehat{\phi}_{j,N}^{(r)}-\phi_{j}\right)^{2},$
(4)

where $\widehat{\phi}_{j,N}^{(r)}$ denotes the estimated value of $\phi_{j}$
for the $r$-th trial. Figure 2 compares the expected SSEs computed for the two
algorithms with various values of $N$. The expected SSE for our proposed
algorithm is smaller than that for the compared algorithm, which is consistent
with the observation of Figure 1 that the confidence interval for our proposed
algorithm is narrower. The expected SSE decays at the rate almost of $1/N$,
which supports the theoretical result given in Theorem 1.

Figure 2: The expected sum of squared errors for Ishigami function as a
function of $N$.

###  4.2 Sobol’ $g$ function

As a high-dimensional example, let us consider the case where $\rho$ is the
uniform distribution over the domain $\Omega=[0,1]^{d}$ for large $d$. With
its simplest form, Sobol’ $g$ function is defined by

$f(\boldsymbol{x})=\prod_{j=1}^{d}\frac{|4x_{j}-2|+a_{j}}{1+a_{j}},$

with non-negative weight parameters $a_{1},\ldots,a_{d}$. For this function,
we can easily obtain

$\sigma^{2}=\prod_{j=1}^{d}\left[1+\frac{1}{3(1+a_{j})^{2}}\right]-1\quad\text{and}\quad\sigma_{u}^{2}=\prod_{j\in
u}\frac{1}{3(1+a_{j})^{2}},$

for any non-empty $u\subseteq[1:d]$. Thus, for an individual variable $x_{j}$,
the main and total effects are given by

$\underline{\tau}^{2}_{\\{j\\}}=\frac{1}{3(1+a_{j})^{2}}\quad\text{and}\quad\overline{\tau}^{2}_{\\{j\\}}=\frac{1}{3(1+a_{j})^{2}}\times\prod_{\begin{subarray}{c}\ell=1\\\
\ell\neq j\end{subarray}}^{d}\left[1+\frac{1}{3(1+a_{\ell})^{2}}\right]$

respectively. The Shapley effect is

$\phi_{j}=\sum_{\begin{subarray}{c}\emptyset\neq u\subseteq[1:d]\\\ j\in
u\end{subarray}}\frac{1}{|u|}\prod_{\ell\in u}\frac{1}{3(1+a_{\ell})^{2}},$

which seems hard to simplify further.

Figure 3: The estimated Shapley effects of individual variables for Sobol’ $g$
function with $N=2^{14}$. The results obtained by our proposed algorithm are
shown in yellow together with the confidence intervals, while those by the
shapleyPermRand estimator are in purple. The exact overall variance, main
effects, Shapley effects and total effects are also shown for comparison.

In what follows, we set $d=10$ and $a_{j}=j-1$. With this choice of $a_{j}$,
the relative importance of individual variables is given by the ascending
order, i.e., $x_{1}$ is the most important, $x_{2}$ is the second most
important, and so on. The Shapley effects for all the input variables are
estimated according to Algorithm 1 and the shapleyPermRand estimator with the
same choice of $N_{v},m,N_{o},N_{i}$ as considered in the first test case. The
result with $N=2^{14}$ is shown in Figure 3. Similarly to Figure 1, we also
plot the exact values of the overall variance, main effects, Shapley effects
and total effects. We can see that all of the exact sensitivity effects
decrease as the index increases. Here again, our proposed algorithm provides
better estimates with narrower confidence intervals of the Shapley effects
than the compared algorithm. The Shapley effect for the variable $x_{10}$, the
least important variable, is estimated to be negative by the compared
algorithm, although it must take a non-negative value by definition. Note that
our proposed algorithm does not ensure non-negativity, so that such an erratic
result may be obtained probabilistically particularly for small $N$. However,
it can be confirmed from our result with $N=2^{14}$ that the estimated Shapley
effects take the values between the corresponding main and total effects,
which agrees with the theory.

An interesting observation is that the confidence interval for our proposed
algorithm gets narrower as the index increases, whereas it stays almost the
same width for the compared algorithm. The crucial difference between the two
algorithms is that our proposed algorithm directly estimates the difference
$\overline{\tau}^{2}_{u+j}-\overline{\tau}^{2}_{u}$ instead of estimating
$\overline{\tau}^{2}_{u}$ and $\overline{\tau}^{2}_{u+j}$ independently. If
the variance of the function
$[(f(\boldsymbol{x})-f(\boldsymbol{y}_{u+j},\boldsymbol{x}_{-(u+j)}))^{2}-(f(\boldsymbol{x})-f(\boldsymbol{y}_{u},\boldsymbol{x}_{-u}))^{2}]/2$
is smaller than those of
$(f(\boldsymbol{x})-f(\boldsymbol{y}_{u+j},\boldsymbol{x}_{-(u+j)}))^{2}/2$
and $(f(\boldsymbol{x})-f(\boldsymbol{y}_{u},\boldsymbol{x}_{-u}))^{2}/2$, a
direct estimation of the difference
$\overline{\tau}^{2}_{u+j}-\overline{\tau}^{2}_{u}$ should be more accurate.
This way our proposed algorithm avoids an unnecessary increment of the
variance of the estimator, particularly for the input variables with small
Shapley effects.

The expected SSEs are estimated by (4) with $R=10$ for the two algorithms.
Figure 4 compares the results with various values of $N$. The expected SSE for
our proposed algorithm is one order of magnitude smaller than that for the
compared algorithm and decays at the rate almost of $1/N$.

Figure 4: The expected sum of squared errors for Sobol’ $g$ function as a
function of $N$.

###  4.3 Plate buckling

Table 1: List of input variables for plate buckling example. CV denotes the coefficient of variation. variable | description | mean | CV | distribution type
---|---|---|---|---
$x_{1}$ | width | 23.808 | 0.028 | normal
$x_{2}$ | thickness | 0.525 | 0.044 | log-normal
$x_{3}$ | yield stress | 44.2 | 0.1235 | log-normal
$x_{4}$ | elastic modulus | 28623 | 0.076 | normal
$x_{5}$ | initial deflection | 0.35 | 0.05 | normal
$x_{6}$ | residual stress | 5.25 | 0.07 | normal

Here we consider a more realistic example from structural engineering. As
explained in [26], let us consider the buckling strength of a rectangular
plate that is supported on all four edges subjected to uniaxial compression.
As described in Table 1, we have $d=6$ input random variables
$x_{1},\ldots,x_{6}$ which are related to the material, geometry and
imperfection of the plate. The buckling strength is defined as a function of
these variables, and is explicitly given by

$f(\boldsymbol{x})=\left(\frac{2.1}{\lambda}-\frac{0.9}{\lambda^{2}}\right)\left(1-\frac{0.75x_{5}}{\lambda}\right)\left(1-\frac{2x_{2}x_{6}}{x_{1}}\right)\qquad\text{with}\quad\lambda=\frac{x_{1}}{x_{2}}\sqrt{\frac{x_{3}}{x_{4}}}.$

Since closed formulas for main and total effects are no longer available for
this example, we estimate them for all the input variables by using the
estimators proposed in the literature. We use the one from [15] for the main
effects, whereas the standard one based on the identity (2) is used for the
total effects, see, e.g., [14]. The sample size is set large enough
($N=2^{24}$) to make sure that the estimates are converged sufficiently. The
Shapley effects for all the input variables are estimated according to
Algorithm 1 and the shapleyPermRand estimator with the same choice of
$N_{v},m,N_{o},N_{i}$ as before. The results are shown in Figure 5.
Interestingly, there is no clear difference between the main and total effects
for all the input variables. This means that each variable does not have
significant interactions with other variables for the buckling strength, so
that the buckling strength can be effectively written as a sum of $d$ one-
dimensional functions $f_{\\{1\\}},\ldots,f_{\\{d\\}}$.

Figure 5: The estimated Shapley effects of individual variables for the plate
buckling example with $N=2^{14}$. The results obtained by our proposed
algorithm are shown in yellow together with the confidence intervals, while
those by the shapleyPermRand estimator are in purple. The estimated overall
variance, main effects, and total effects are also shown for comparison.

As already pointed out, the Shapley effect takes its value between the
corresponding main and the total effects for each input variable. Since the
difference between the latter two values is quite small in this example, it
might be challenging to estimate the Shapley effect in a way that the
resulting value is bracketed by them. However, our proposed algorithm with
$N=2^{14}$ already provides good estimates for all the variables. Similarly to
the second test case, the confidence interval for our proposed algorithm is
narrower for a variable with smaller Shapley effect, while it is not the case
for the compared algorithm.

For this example, the expected SSE is estimated by

$\frac{1}{R-1}\sum_{r=1}^{R}\sum_{j=1}^{d}\left(\widehat{\phi}_{j,N}^{(r)}-\overline{\phi}_{j}\right)^{2}\quad\text{with}\quad\overline{\phi}_{j}=\frac{1}{R}\sum_{r=1}^{R}\widehat{\phi}_{j,N}^{(r)}.$

Again we choose $R=10$. The results with various values of $N$ are compared in
Figure 6. We can confirm that the expected SSE for our proposed algorithm is
always smaller than that for the compared algorithm and decays at the rate
almost of $1/N$.

Figure 6: The expected sum of squared errors for the plate buckling example as
a function of $N$.

##  5 Possible extensions

Finally we discuss two possible extensions of our Monte Carlo estimator. The
first extension allows for an additional cost reduction from $(d+1)N$ to
$dN+1$ or $dN$, which might be useful when $d$ is small but $N$ is large. The
second extension considers the case where the input variables are dependent.

###  5.1 Additional cost saving

In Algorithm 1, we generate two independent sequences for $\boldsymbol{x}$ and
$\boldsymbol{y}$, and consider a pairing
$(\boldsymbol{x}^{(n)},\boldsymbol{y}^{(n)})$ for $1\leq n\leq N$. Because of
the linearity of expectation, however, the unbiasedness of our Monte Carlo
estimator still holds even by generating only one random sequence for
$\boldsymbol{x}$ and then considering a pairing
$(\boldsymbol{x}^{(n)},\boldsymbol{x}^{(n+1)})$ for $1\leq n\leq N$ instead.
This strategy is called the _winding star_ sampling [27] and enables to reduce
the necessary cost from $(d+1)N$ to $dN+1$. The cost of $dN$ is possible by
setting $\boldsymbol{x}^{(N+1)}=\boldsymbol{x}^{(1)}$, although the algorithm
then becomes not extensible in $N$.

The variance of this lower-cost estimator is finite and decays at the $1/N$
rate under the same assumption given in Theorem 1. However, because of the
dependence between consecutive samples, an unbiased estimation of the variance
of the estimator is not possible in the same way as the second item of Theorem
1.

###  5.2 Dependent input case

For $\boldsymbol{x},\boldsymbol{y}\sim\rho$ and any subset $u\subseteq[1:d]$,
it holds that $(\boldsymbol{x}_{u},\boldsymbol{y}_{-u})\sim\rho$ as long as
the input variables are independent. This property plays a fundamental role in
establishing the simple concatenation approach for two independent samples
$\boldsymbol{x}^{(n)},\boldsymbol{y}^{(n)}$ in Algorithm 1. In fact, if the
input variables are dependent, we generally have
$(\boldsymbol{x}_{u},\boldsymbol{y}_{-u})\not\sim\rho$, which prevents us from
extending Algorithm 1 to the dependent input case in a straightforward manner.

In order to adapt Algorithm 1 to the dependent input variables, we need the
following modifications:

  1. 1.

In the first step, generate only $\boldsymbol{x}^{(n)}\in\Omega$ and
$\boldsymbol{\pi}^{(n)}$ randomly.

  2. 2.

In the third step, generate
$\boldsymbol{y}^{(n)}_{\\{\pi(1),\ldots,\pi(\ell)\\}}$ conditional only on
$\boldsymbol{x}^{(n)}_{-\\{\pi(1),\ldots,\pi(\ell)\\}}$ and independently of
the existing ones
$\boldsymbol{y}^{(n)}_{\\{\pi(1)\\}},\ldots,\boldsymbol{y}^{(n)}_{\\{\pi(1),\ldots,\pi(\ell-1)\\}}$
for all $1\leq\ell\leq d$.

This way, we still have an unbiased Monte Carlo estimator for the Shapley
effects with the cost of $(d+1)N$. However, the independence between
$\boldsymbol{y}^{(n)}_{\\{\pi(1),\ldots,\pi(\ell-1)\\}}$ and
$\boldsymbol{y}^{(n)}_{\\{\pi(1),\ldots,\pi(\ell)\\}}$ may lead to a large
variance in estimating the difference
$\overline{\tau}^{2}_{\\{\pi(1),\ldots,\pi(\ell)\\}}-\overline{\tau}^{2}_{\\{\pi(1),\ldots,\pi(\ell-1)\\}}$.
Further detailed investigation is needed to address this issue.

## Funding

The work of T.G. was supported by JSPS KAKENHI Grant Number 20K03744.

## Acknowledgments

The author would like to thank the associate editor and the referees for
helpful suggestions and remarks. He also thanks Elmar Plischke for pointing
out the relevant literature and suggesting a way to simplify the Matlab code.

## References

  * [1] A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana and S. Tarantola, _Global Sensitivity Analysis. The Primer_ , John Wiley & Sons, Ltd, 2008.
  * [2] B. Iooss and P. Lemaître, A review on global sensitivity analysis methods, In: G. Dellino and C. Meloni (eds.) Uncertainty Management in Simulation-Optimization of Complex Systems. Springer, Boston, 2015, pp. 101–122.
  * [3] T. Homma and A. Saltelli, Importance measures in global sensitivity analysis of nonlinear models, Reliab. Eng. Syst. Safe. 52 (1996), 1–17.
  * [4] E. Borgonovo, G. E. Apostolakis, S. Tarantola and A. Saltelli, Comparison of global sensitivity analysis techniques and importance measures in PSA, Reliab. Eng. Syst. Safe. 79 (2003), 175–185.
  * [5] J. C. Helton, J. D. Johnson, C. J. Sallaberry and C. B. Storlie, Survey of sampling-based methods for uncertainty and sensitivity analysis, Reliab. Eng. Syst. Safe. 91 (2006), 1175–1209.
  * [6] I. M. Sobol’, Sensitivity estimates for nonlinear mathematical models, Math. Model. Comput. Exper. 1 (1993), 407–414.
  * [7] I. M. Sobol’, Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates, Math. Comput. Simulation 55 (2001), 271–280.
  * [8] A. B. Owen, Sobol’ indices and Shapley values, SIAM/ASA J Uncertainty Quantification 2 (2014), 245–251.
  * [9] E. Song, B. L. Nelson and J. Staum, Shapley effects for global sensitivity analysis: theory and computation, SIAM/ASA J Uncertainty Quantification 4 (2016), 1060–1083.
  * [10] A. B. Owen and C. Prieur, On Shapley value for measuring importance of dependent inputs, SIAM/ASA J Uncertainty Quantification 5 (2017), 986–1002.
  * [11] B. Broto, F. Bachoc and M. Depecker, Variance reduction for estimation of Shapley effects and adaptation to unknown input distribution, SIAM/ASA J Uncertainty Quantification 8 (2020), 693–716.
  * [12] J. Castro, D. Gómez and J. Tejada, Polynomial calculation of the Shapley value based on sampling, Comput. Oper. Res. 36 (2009), 1726–1730.
  * [13] A. Saltelli, Making best use of model evaluations to compute sensitivity indices, Comput. Phys. Comm. 145 (2002), 280–297.
  * [14] A. Saltelli, P. Annoni, I. Azzini, F. Campolongo, M. Ratto and S. Tarantola, Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index, Comput. Phys. Comm. 181 (2010), 259–270.
  * [15] A. B. Owen, Better estimation of small Sobol’ sensitivity indices, ACM Trans. Model. Comput. Simul. 23 (2013), Article 11.
  * [16] A. Janon, T. Klein, A. Lagnoux, M. Nodet and C. Prieur, Asymptotic normality and efficiency of two Sobol index estimators, ESAIM: PS 18 (2014), 342–364.
  * [17] T. Goda, Computing the variance of a conditional expectation via non-nested Monte Carlo, Oper. Res. Lett. 45 (2017), 63–67.
  * [18] S. Tarantola, D. Gatelli and T. A. Mara, Random balance designs for the estimation of first order global sensitivity indices, Reliab. Eng. Syst. Safe. 91 (2006), 717–727.
  * [19] B. Sudret, Global sensitivity analysis using polynomial chaos expansions, Reliab. Eng. Syst. Safe. 93 (2008), 964–979.
  * [20] E. Plischke, An effective algorithm for computing global sensitivity indices (EASI), Reliab. Eng. Syst. Safe. 95 (2010), 354–360.
  * [21] K. Konakli and B. Sudret, Global sensitivity analysis using low-rank tensor approximations, Reliab. Eng. Syst. Safe. 156 (2016), 64–83.
  * [22] A. Antoniadis, S. Lambert-Lacroix and J.-M. Poggi, Random forests for global sensitivity analysis: A selective review, Reliab. Eng. Syst. Safe. 206 (2021), 107312.
  * [23] C. Prieur and S. Tarantola, Variance-based sensitivity analysis: Theory and estimation algorithms, In: R. Ghanem, D. Higdon and H. Owhadi (eds.) Handbook of Uncertainty Quantification. Springer, Switzerland, 2017, pp. 1217–1239.
  * [24] A. B. Owen, _Monte Carlo theory, methods and examples_ , http://statweb.stanford.edu/~owen/mc/ Last accessed 1 September 2020.
  * [25] B. Iooss et al., Package ‘sensitivity’, version 1.23.1, 2020.
  * [26] M. D. Shields and J. Zhang, The generalization of Latin hypercube sampling, Reliab. Eng. Syst. Safe. 148 (2016), 96–108.
  * [27] M. J. W. Jansen, W. A. H. Rossing and R. A. Daamen, Monte Carlo estimation of uncertainty contributions from several independent multivariate sources, In: J. Grasman and G. van Straten (eds.) Predictability and Nonlinear Modelling in Natural Sciences and Economics. Springer, Dordrecht, 1994, pp. 334–343.

##  Appendix A Matlab implementation

    n = 2^14;                    % sample size
    d = 10;                      % dimension
    x = rand(n,d);               % Monte Carlo samples for x
    y = rand(n,d);               % Monte Carlo samples for y
    [~, pm] = sort(rand(n,d),2); % random permutation matrix

    func = @(x,a)prod((abs(4.*x-2)+a)./(1+a),2); % Sobol g function
    a = (0:d-1);

    z = x;
    fz1 = func(z,a);
    fx = fz1;

    phi1 = zeros(1,d); % initialization
    phi2 = zeros(1,d); % initialization
    for j=1:d
        ind = bsxfun(@eq,pm(:,j),1:d); % compare j-th column with 1:d
        z(ind) = y(ind);
        fz2 = func(z,a);
        fmarg = ((fx-fz1/2-fz2/2).*(fz1-fz2))’;
        phi1 = phi1 + fmarg*ind/n;
        phi2 = phi2 + fmarg.^2*ind/n;
        fz1 = fz2;
    end

    s_all = sum(phi1);            % variance of function
    phi2 = (phi2-phi1.^2)./(n-1); % variance of Shapley estimates

    disp(phi1);                  % Shapley estimates
    disp(phi1-1.96*sqrt(phi2));  % lower confidence limit
    disp(phi1+1.96*sqrt(phi2));  % upper confidence limit


# Thermodynamics and SARS-CoV-2: neurological effects in post-Covid 19
syndrome

Umberto Lucia 1,a, Giulia Grisolia 1,b & Thomas S. Deisboeck 2,c
1 Dipartimento Energia “Galileo Ferraris”,
Corso Duca degli Abruzzi 24, 10129 Torino, Italy
2 Department of Radiology,
Harvard-MIT Martinos Center for Biomedical Imaging,
Massachusetts General Hospital and
Harvard Medical School, Charlestown, MA 02129, USA
a<EMAIL_ADDRESS>
b<EMAIL_ADDRESS>
c<EMAIL_ADDRESS>
###### Abstract

There is increasing evidence that infection with SARS-CoV-2 can cause a
spectrum of neurological symptoms. In this paper, we develop a theoretical
concept underlying such neurological COVID-19 consequences by employing a non-
equilibrium thermodynamic approach that allows linking the neuronal electric
potential with a virus-induced pH variation. Our theoretical findings support
further experimental work on therapeutically correcting electrolyte
imbalances, such as Na+ and K+, to attenuate the neurological effects of SARS-
CoV-2.

Keyword: Irreversible thermodynamics; Membrane potential; Neurological
effects; Post-Covid Syndrome; SARS-CoV-2.

##  1 Introduction

Infection with the new severe acute respiratory syndrome coronavirus 2, i.e.,
SARS-CoV-2, can lead to a variety of clinical symptoms ranging from
respiratory and circulatory effects to neurological ones; some of these
symptoms require immediate therapeutic intervention to try to stabilize
vulnerable patients, while others are feared to potentially cause long term
morbidity in some. (Gu et al., 2020).

In this theoretical paper, we focus our analysis on the neurological
consequences of SARS-CoV-2. Since the virus has seldom been found in the
cerebrospinal fluid of the patients, the damage caused by the anti-virus
immune response-mediated damage seems to be the culprit (Solomon, 2021).
Severe neurological manifestations, albeit rare by comparison, run the gamut
from Guillain-Barre (demyelinating polyneuritis) syndrome, and ischemic stroke
(Berlit et al., 2020; Kandimalla et al., 2020) to encephalitis (Ellul et al.,
2020) while milder symptoms include temporary memory loss, altered mental
state or ‘brain fog’ as well as olfactory and gustatory dysfunctions (Agyeman
et al., 2020). In fact, partial or complete loss of smell, dysosmia or
anosmia, as well as dysgeusia (loss of taste) are common symptoms during a
SARS-CoV-2 infection, even in the absence of any other symptoms (Solomon,
2021; Kanberg et al., 2020a, b; Matschke et al., 2020a; Paniz-Mondolfi et al.,
2020; Meinhardt et al., 2021); as to the loss of smell, the virus seems to
infect the olfactory epithelium, rather than the sensory neurons themselves
(Kanberg et al., 2020a, b; Matschke et al., 2020a; Paniz-Mondolfi et al.,
2020; Meinhardt et al., 2021). This notwithstanding, in autopsy studies of
patients dying of COVID viral RNA transcripts were found in brain tissue and
viral proteins in the endothelial cells within the olfactory bulb (Song et
al., 2021); in another study, inflammatory changes appear particularly
extensive in the brainstem in proximity to cranial nerve origins (Matschke et
al., 2020b).

At present, the lasting consequences of this neuro-invasion are not fully
understood, neither is the optimal therapeutic strategy as, for instance,
targeting an over-active immune response with corticosteroids would be
potentially dangerous in the presence of virus (Marshall, 2020). In this
paper, we therefore develop a non-equilibrium thermodynamic analysis of this
mechanism, in order to suggest a possible new viewpoint that may hold promise
for designing future therapies.

##  2 Materials and Methods

Protein phosphorylation is a fundamental biochemical mechanism regulating of
the cell functions, due to its ability to activate and deactivate some enzymes
and receptors (Rudolph et al., 2006; Strong, 2002). In this context the
relation with kinases is of particular interest because kinases are related to
cellular transduction signalling (Ardito et al., 2017).

Ions actively cross the cell membrane against its electrochemical potential by
deriving the required energy from the hydrolysis of ATP, where the H+-ATPase
plays a fundamental role; this is related to the movement of positive charges
into the cell, by generating large membrane voltage (inside negative and
outside positive) and a pH gradient (Nakanishi-Matsui et al., 2010; Stevens &
Forgac, 1997; Tuszynski & Kurzynski, 2003):

$\text{ATP}+\text{H}_{2}\text{O}\rightarrow\text{ADP}+\text{P}$ (1)
$\text{H}^{+}_{out}\rightarrow\text{H}^{+}_{in}$ (2)

with a subsequent variation of the pH because (Ashrafuzzaman & Tuszynski,
2013):

$\Delta\text{pH}=\frac{F}{2.3RT}\,\big{(}\Delta\phi_{m}-\Delta
G_{\text{H}^{+}}\big{)}$ (3)

where $G$ depicts the Gibbs potential. The phosphorylation potential,
$\Delta\bar{g}_{p}$ [kJ mol-1], is well known and described by the following
equation (Tuszynski & Kurzynski, 2003; Grabe et al., 2000; Lucia, 2015b, a;
Lucia & Grisolia, 2018a; Lucia et al., 2014, 2018; Lucia & Grisolia, 2017;
Lucia et al., 2017, 2016; Lucia & Grisolia, 2018b):

$\Delta\bar{g}_{p}=-nF\Delta\phi$ (4)

where $n$ is the number of moles of ions per ATP synthesized, $F=96.485\times
10^{3}$ A s mol-1 is the Faraday constant, and $\Delta\phi$ stands for the
membrane potential.

The movement of the ions can be analysed by introducing the Onsager general
phenomenological relationships as they pertain to both the electrochemical
potential and the heat flux (Yourgrau et al., 1982; Callen, 1960; Lucia &
Grisolia, 2020a, b, c; Goupil, 2011):

$\left\\{\begin{array}[]{l}\mathbf{J}_{e}=-L_{11}\,\dfrac{\nabla\mu_{e}}{T}-L_{12}\,\dfrac{\nabla
T}{T^{2}}\\\ \\\
\mathbf{J}_{Q}=-L_{21}\,\dfrac{\nabla\mu_{e}}{T}-L_{22}\,\dfrac{\nabla
T}{T^{2}}\end{array}\right.$ (5)

where $\textbf{J}_{e}$ is the current density [A m-2], $\textbf{J}_{Q}$
denotes the heat flux [W m-2], $\mu_{e}=\mu+ze\phi$ is the electrochemical
potential [J mol-1], with $\mu$ the chemical potential [J mol-1], $ze$ the
electric charge [A s mol-1], and $\phi$ the membrane potential [V],
respectively; $T$ is the living cell temperature and $L_{ij}$ represent the
phenomenological coefficients, such that (Katchalsky & Currant, 1965)
$L_{12}(\mathbf{B})=L_{21}(-\mathbf{B})$ (Onsager-Casimir relation (Degroot &
Mazur, 1962)), and $L_{11}\geq 0$ and $L_{22}\geq 0$, and (Katchalsky &
Currant, 1965) $L_{11}L_{22}-L_{12}L_{21}>0$.

The result consists of a model of the life cycle based on two related
processes (Lucia & Grisolia, 2020a, b):

  * •

A continuous energy generation (metabolism), due to ion fluxes: The ion and
metabolite fluxes can be described by imposing $\mathbf{J}_{e}\neq\mathbf{0}$
and $\mathbf{J}_{Q}=\mathbf{0}$;

  * •

A continuous heat flux from the cell to its microenvironment: The heat
exchange towards the environment can be described by imposing
$\mathbf{J}_{e}=\mathbf{0}$ and $\mathbf{J}_{Q}\neq\mathbf{0}$.

In this way, we can split the life cycle into two thermodynamic processes, as
it is usually done in irreversible thermodynamics for any complex process
(Callen, 1960).

Now, if ion and metabolite fluxes occur, $\mathbf{J}_{e}\neq\mathbf{0}$ and
$\mathbf{J}_{Q}=\mathbf{0}$, it follows (Callen, 1960; Yourgrau et al., 1982;
Lucia & Grisolia, 2020a)

$\frac{d\mu_{e}}{dT}=-\frac{L_{21}}{L_{11}}\,\frac{1}{T}$ (6)

with a related heat flux (Callen, 1960; Yourgrau et al., 1982):

$\frac{du}{dt}=-\nabla\cdot\mathbf{J}_{Q}$ (7)

where $u$ is the internal energy density [W m-3].

Living cells exchange heat power towards their environment by convection, and
so, we can write (Lucia & Grisolia, 2020d)

$\frac{du}{dt}\,dV=\delta\dot{Q}=-\alpha\,(T-T_{0})\,dA$ (8)

where $\alpha\approx 0.023Re^{0.8}Pr^{0.35}\lambda/\langle R\rangle$ is the
coefficient of convection, $A$ the area of the external surface of the cell
membrane, $V$ is the cell volume, $T$ depicts the mean temperature of the
external surface of the cell’s membrane, and $T_{0}$ is the temperature of the
cell environment.

So, considering Equations (7) and (8), and the Divergence Theorem (Apostol,
1969), the heat flux can be written as:

$J_{Q}=\alpha\,(T-T_{0})$ (9)

and the related power flux yields:

$\dot{Q}=\int_{A}\mathbf{J}_{Q}\cdot\hat{\mathbf{n}}dA=\alpha\,(T-T_{0})A$
(10)

Furthermore, considering Equation (5), together with the second hypothesis of
our modelling ($\mathbf{J}_{e}=\mathbf{0}$, $\mathbf{J}_{Q}\neq\mathbf{0}$),
it follows (Lucia & Grisolia, 2020a):

$\frac{d\mu_{e}}{d\ell}=\frac{T\,J_{Q}}{\Bigg{(}L_{22}\frac{L_{11}}{L_{12}}-L_{21}\Bigg{)}}=-\frac{\alpha\,T(T-T_{0})}{\Bigg{(}L_{22}\frac{L_{11}}{L_{12}}-L_{21}\Bigg{)}}$
(11)

where $\ell$ is the length of a cell membrane and $|\nabla\mu_{e}|\approx
d\mu_{e}/d\ell$. This relation is the link between the cell membrane’s
electric potential and the temperature of the cell itself.

Equations (9) and (11) allow us to obtain:

$J_{Q}=\alpha\,(T-T_{0})=-\frac{1}{T}\Bigg{(}L_{22}\frac{L_{11}}{L_{12}}-L_{21}\Bigg{)}\,\frac{d\mu_{e}}{d\ell}$
(12)

where:

$\Bigg{(}L_{22}-L_{21}\frac{L_{12}}{L_{11}}\Bigg{)}=K_{J}T^{2}$ (13)

with $K_{J}$ being the Thomson coefficient. Consequently, it follows:

$\frac{\partial\mu_{e}}{\partial\ell}=\frac{\partial\mu_{e}}{\partial
T}\,\frac{\alpha}{K_{J}}\big{(}T_{surf}-T_{0}\big{)}$ (14)

from which, taking into account that $\mu_{e}=\mu+ze\phi$, becomes:

$\frac{\partial\mu}{\partial\ell}=-ze\frac{d\phi}{d\ell}+\frac{\partial\mu_{e}}{\partial
T}\,\frac{\alpha}{K_{J}}\big{(}T_{surf}-T_{0}\big{)}$ (15)

Now, considering Equation (3) we can obtain:

$\frac{\partial\mu_{e}}{\partial
T}=\frac{K_{J}}{\alpha}\frac{F+ze}{T_{surf}-T_{0}}\,\frac{d\phi}{d\ell}-\frac{K_{J}}{\alpha}\frac{2.3R\,T_{0}}{T_{surf}-T_{0}}\,\frac{d\text{pH}}{d\ell}$
(16)

which links the electrochemical potential to the pH.

In order to understand the effect of SARS-CoV-2 on the brain, we use a simple
model of information coding, introduced previously in our thermodynamic
analysis of Alzheimer’s disease (Lucia et al., 2020b).

Specifically, a brain cell requires a Na+-inflow, and a countering flow of
K+-outflow to develop the functionality of processing signals (Bustamante et
al., 2004). During this function, the consumption of one molecule of ATP
(adenosine triphosphate) requires that the membrane pump extrudes 3 Na+-ions
and imports 1 K+-ion. At the stationary state, a neuron maintains its pump
current:

$I_{p}=\frac{\Delta\phi_{\text{Na}^{+}}-\Delta\phi_{m}}{R_{\text{Na}^{+}}}+\frac{\Delta\phi_{\text{K}^{+}}-\Delta\phi_{m}}{R_{\text{K}^{+}}}$
(17)

where $R_{i}$ ($i=$ [Na+] or [K+]) stands for the electric resistance of the
ion considered during its current flux through the membrane, and
$\Delta\phi_{m}$ is the membrane electric potential variation,
$\Delta\phi_{\text{Na}^{+}}$ is the electric potential variation due to the
Na+-flux (Attwell & Laughlin, 2001; Goldman, 1943; Ashrafuzzaman & Tuszynski,
2013):

$\Delta\phi_{\text{Na}^{+}}=-\frac{RT}{F}\,\ln\Bigg{(}\frac{[\text{Na}^{+}]_{f}}{[\text{Na}^{+}]_{i}}\Bigg{)}$
(18)

$\Delta\phi_{\text{K}^{+}}$ is the electric potential variation due to the
K+-flux (Attwell & Laughlin, 2001; Goldman, 1943; Ashrafuzzaman & Tuszynski,
2013):

$\Delta\phi_{\text{K}^{+}}=-\frac{RT}{F}\,\ln\Bigg{(}\frac{[\text{K}^{+}]_{f}}{[\text{K}^{+}]_{i}}\Bigg{)}$
(19)

where $R=8314$ J mol-1 K-1 denotes the constant of the ideal gasses,
$F=96,485$ C mol-1 is the Faraday constant, $f$ and $i$ means final and
initial respectively, and they are referred to the initial and finale state of
the neuronal signalling process, and $T$ is the temperature, and

$R_{in}=\frac{1}{\dfrac{1}{R_{\text{Na}^{+}}}+\dfrac{1}{R_{\text{K}^{+}}}}$
(20)

under the biochemical constraint:

$\frac{d[\text{Na}^{+}]}{dt}=-\frac{d[\text{K}^{+}]}{dt}$ (21)

where [A] is the concentration of the A-ion (Na+/K+). In order to maintain a
normal membrane potential of around (Attwell & Laughlin, 2001) -70 mV a neuron
($R_{in}=$ 200 M$\Omega$ of input resistance) requires an influx of around
$1.02\times 10^{9}$ Na+-K+ ions s-1 ($\Delta\phi_{Na^{+}}=-50$ mV and
$\Delta\phi_{K^{+}}=-100$ mV) which necessitates $3.42\times 10^{8}$
hydrolysed ATP molecules s-1 (Attwell & Laughlin, 2001), consumed at a rate of
$I_{p}/F$: it generates a pump current $I_{p}$ of $1.63\times 10^{-10}$ A.

Next, considering Equations (21), (18) and (19) it follows that (19) becomes:

$\displaystyle\frac{\partial\mu_{e,\text{Na}^{+}}}{\partial
T}=-\frac{K_{J}RT_{0}}{F\alpha}\frac{F+ze}{T_{surf}-T_{0}}\,\frac{1}{[\text{Na}^{+}]}\,\frac{d[\text{Na}^{+}]}{d\ell}-\frac{K_{J}}{\alpha}\frac{2.3R\,T_{0}}{T_{surf}-T_{0}}\,\frac{d\text{pH}}{d\ell}$
(22) $\displaystyle\frac{\partial\mu_{e,\text{K}^{+}}}{\partial
T}=-\frac{K_{J}RT_{0}}{F\alpha}\frac{F+ze}{T_{surf}-T_{0}}\,\frac{1}{[\text{K}^{+}]}\,\frac{d[\text{K}^{+}]}{d\ell}-\frac{K_{J}}{\alpha}\frac{2.3R\,T_{0}}{T_{surf}-T_{0}}\,\frac{d\text{pH}}{d\ell}$

which points out that, in order to maintain a normal chemical potential the
cell, or neuron in this case, must actively change its concentration of ions
if a change in the pH occurs as a result of the viral infection.

As a consequence of the previous steps, a density entropy rate due to
irreversibility (dissipation function (Yourgrau et al., 1982)) is generated
(Lucia & Grisolia, 2017):

$\sigma=-\frac{1}{T_{0}}\,\sum_{i=1}^{N}\mathbf{J}_{i}\cdot\nabla\mu_{i}\geq
0$ (23)

where $T_{0}$ represents the environmental temperature,
$\sum_{i=1}^{N}\mu_{i}\,\mathbf{J}_{i}$ is the contribution of the inflows and
outflows, and $\mu$ denotes the chemical potential. So, using the previous
relations we obtain:

$\displaystyle\sigma$
$\displaystyle\approx-\frac{R\,J_{\text{Na${}^{+}$}}\,(ze)_{\text{Na${}^{+}$}}}{F\ell}\,\ln\Bigg{(}\frac{[\text{Na}^{+}]_{f}}{[\text{Na}^{+}]_{i}}\Bigg{)}+\frac{J_{\text{Na${}^{+}$}}}{T_{0}}\,\frac{\partial\mu_{e,\text{Na${}^{+}$}}}{\partial
T}\,\frac{\alpha}{K_{J}}\big{(}T_{surf}-T_{0}\big{)}+$ (24)
$\displaystyle-\frac{R\,J_{\text{K${}^{+}$}}\,(ze)_{\text{K${}^{+}$}}}{F\ell}\,\ln\Bigg{(}\frac{[\text{K}^{+}]_{f}}{[\text{K}^{+}]_{i}}\Bigg{)}+\frac{J_{\text{K${}^{+}$}}}{T_{0}}\,\frac{\partial\mu_{e,\text{K${}^{+}$}}}{\partial
T}\,\frac{\alpha}{K_{J}}\big{(}T_{surf}-T_{0}\big{)}$

Starting from this last equation, and considering the previous condition of
non-negative entropy density (Katchalsky & Currant, 1965), we arrive at the
following condition:

$\displaystyle\frac{R\,J_{\text{Na${}^{+}$}}\,(ze)_{\text{Na${}^{+}$}}}{F\ell}$
$\displaystyle\ln\Bigg{(}\frac{[\text{Na}^{+}]_{f}}{[\text{Na}^{+}]_{i}}\Bigg{)}-\frac{J_{\text{Na${}^{+}$}}}{T_{0}}\,\frac{\partial\mu_{e,\text{Na${}^{+}$}}}{\partial
T}\,\frac{\alpha}{K_{J}}\big{(}T_{surf}-T_{0}\big{)}\leq$ (25)
$\displaystyle-\frac{R\,J_{\text{K${}^{+}$}}\,(ze)_{\text{K${}^{+}$}}}{F\ell}\,\ln\Bigg{(}\frac{[\text{K}^{+}]_{f}}{[\text{K}^{+}]_{i}}\Bigg{)}+\frac{J_{\text{K${}^{+}$}}}{T_{0}}\,\frac{\partial\mu_{e,\text{K${}^{+}$}}}{\partial
T}\,\frac{\alpha}{K_{J}}\big{(}T_{surf}-T_{0}\big{)}$

which suggests that to maintain stability of the neuronal cell system the
effect of sodium fluxes is less pronounced than that of potassium.

##  3 Results

Our conjecture yielded Equations (3), (15) and (17). These equations highlight
the link between the neuronal signalling process and the neurons’ membrane
transport.

Specifically, Equation (3) states that a change in pH determines a related
variation in membrane potential and in proton flux, related to Gibbs energy
and chemical potential (Equation (15)). As such, a variation in the neuronal
ion current pump occurs, with the consequence of modifying the membrane
potential related to Na+ and K+. In turn, this change determines a modulation
in the concentrations of these chemical species with a symmetry breaking in
the stationary condition for the neurons. Consequently, the signalling process
changes which may offer an explanation for the neurological consequences of
infection with SARS-CoV-2. Indeed, we proved (Lucia et al., 2020a) that SARS-
CoV-2 leads to changes in pH-homeostasis due to modifications of H+-fluxes.
Taken together, this suggests that a promising therapeutic strategy would seek
to control the neurons’ membrane potential through manipulation of the ions
responsible for it.

##  4 Discussion and Conclusions

In patients COVID-19 is characterized by a wide variety of symptoms, some of
them neurological with at times very significant morbidity (Marshall, 2020).
To shed more light onto this, we have developed our non-equilibrium
thermodynamics approach that focuses on the effect of SARS-CoV-2-induced pH
variation on the ion and thermal fluxes across the cell membrane. As
emphasized in Eq. 16 and 22, this ‘pH-ion flux’ link can potentially support
the development of new therapeutic strategies to combat SARS-CoV-2’s neuronal
effects.

Electrolyte imbalances are common in Covid-19 patients, including hypokalemia
which has been found to be an independent predictor of the requirement for
mechanical ventilation (Moreno-P et al., 2020) and which may predispose to
cardiac complications and necessitate potassium supplement therapy (Chen et
al., 2020), and hypocalcemia (Zhou et al., 2020) which can lead to muscle
twitches and tremors. Furthermore, hyponatremia i.e. a serum sodium [Na+]
level of below 135 mmol L-1 can cause serious central nervous system symptoms
ranging from headaches, lethargy, and cramps to seizures, coma, and
respiratory arrest (Giuliani & Peri, 2014). Interestingly, hyponatremia has
indeed been noted as an early sign in Covid-19 infections (Gheorghe et al.,
2021). Clinical management of this common, multifactorial sodium imbalance in
Covid-19 patients depends on the exact etiology and therefore can involve
electrolyte replacement therapy in patients suffering from hyponatremia
primarily due to fluid losses or infusion of hypertonic saline in cases of
insufficient antidiuretic hormone secretion, seen as a consequence of the
systemic inflammation triggered by the virus (Gheorghe et al., 2021). It is
intriguing in this context that a recent article, available so far only as
preprint, describes in vitro experiments where hypertonic saline (1.5% NaCl)
inhibits SARS-Cov-2 replication completely, presumed to be achieved by cell
plasma membrane depolarization and intracellular energy deprivation (Machado
et al., 2020). If these experimental findings can be confirmed, our
theoretical conjecture supports advancing this concept from in vitro to in
vivo studies, in an effort to gain insights if hypertonic saline could help
mitigate SARS-CoV-2 effects on the central nervous system. Finally, in light
of these very recent experimental findings on sodium, Eq. 25 may deserve
attention as it indicates that the impact of correcting potassium levels may
be even more pronounced. Still, it must be noted that in clinics addressing
electrolyte imbalances is highly non-trivial and extensive, experimental work
would be required to properly assess risk vs. benefit of serum adjustments of
Na+ and/or K+.

In summary, our thermodynamics approach conceptualizes that SARS-CoV-2-induced
pH fluctuations trigger ion-flux changes across the neuron cell membrane which
in turn alters signalling throughout the system. Cautiously extrapolated, this
may explain some of the neurological symptoms seen in COVID patients and it
supports further research in therapeutically addressing electrolyte
imbalances.

## Authors contributions

Conceptualization, U.L and T.S.D.; methodology, U.L., G.G. and T.S.D.;
software, G.G.; validation, U.L., T.S.D. and G.G.; formal analysis, U.L.;
investigation, G.G.; resources, U.L.; data curation, G.G.; writing-original
draft preparation, U.L., G.G. and T.S.D.; writing—review and editing, U.L.,
G.G. and T.S.D.; visualization, G.G.; supervision, U.L., T.S.D.; project
administration, U.L.; funding acquisition, U.L. All authors have read and
agreed to the published version of the manuscript.

## References

  * Agyeman et al. (2020) Agyeman, A. A., Chin, K. L., Landersdorfer, C. B., Liew, D., , & Ofori-Asenso, R. (2020).  Smell and taste dysfunction in patients withcovid-19: A systematic reviewand meta-analysis.  Mayo Clinic Proceedings, 95, 1621–1631.
  * Apostol (1969) Apostol, T. S. (1969).  Calculus. Volume 2: Multi-variable calculus and linear algebra with applications to differential equations and probability.  Hoboken: Wiley.
  * Ardito et al. (2017) Ardito, F., Giuliani, M., Perrone, D., Troiano, G., & Muzio, L. L. (2017).  The crucial role of protein phosphorylation in cell signaling and its use as targeted therapy.  International Journal of Molecular Medicine, 40, 271–280.
  * Ashrafuzzaman & Tuszynski (2013) Ashrafuzzaman, M., & Tuszynski, J. (2013).  Membrane Biophysics.  Berlin: Springer.
  * Attwell & Laughlin (2001) Attwell, D., & Laughlin, S. B. (2001).  An energy budget for signaling in the grey matter of the brain.  Journal of Cerebral Blood Flow and Metabolism, 21, 1133–1145.
  * Berlit et al. (2020) Berlit, P., Bösel, J., Gahn, G., Isenmann, S., Meuth, S. G., Nolte, C. H., Pawlitzki, M., Rosenow, F., Schoser, B., Thomalla, G., & Hummel, T. (2020).  “neurological manifestations of covid-19” - guideline of the german society of neurology.  Neurological Research and Practice, 2, 51.
  * Bustamante et al. (2004) Bustamante, C., Chemla, Y. R., Forde, N. R., & Izhaky, D. (2004).  Mechanical processes in biochemistry.  Annual Review of Biochemistry, 73, 705–748.
  * Callen (1960) Callen, H. B. (1960).  Thermodynamics.  New York: Wiley.
  * Chen et al. (2020) Chen, D., Li, X., Song, Q., Hu, C., Su, F., Dai, J., Ye, Y., Huang, J., & Zhang, X. (2020).  Assessment of hypokalemia and clinical characteristics in patients with coronavirus disease 2019 in wenzhou, china.  JAMA Network Open, 3(6), e2011122–e2011122.
  * Degroot & Mazur (1962) Degroot, S. R., & Mazur, P. (1962).  Non-Equilibrium Thermodynamics.  Amsterdam: North-Holland Publishing Company.
  * Ellul et al. (2020) Ellul, M. A., Benjamin, L., Singh, B., Lant, S., Michael, B. D., Easton, A., Kneen, R., Defres, S., Sejvar, J., & Solomon, T. (2020).  Neurological associations of covid-19.  The Lancet Neurology, (pp. 767–783).
  * Gheorghe et al. (2021) Gheorghe, G., Ilie, M., Bungau, S., Stoian, A. M. P., Bacalbasa, N., & Diaconu, C. C. (2021).  Is there a relationship between covid-19 and hyponatremia?  Medicina, 57, 55.
  * Giuliani & Peri (2014) Giuliani, C., & Peri, A. (2014).  Effects of hyponatremia on the brain.  Journal of Clinical Medicine, 3, 1163–1177.
  * Goldman (1943) Goldman, D. E. (1943).  Potential impedance, and rectification in membranes.  Journal of General Physiology, 27, 37–60.
  * Goupil (2011) Goupil, C. (2011).  Thermodynamics of thermoelectricity.  In T. Mizutani (Ed.) Thermodynamics. Shanghai: IntechOpen.
  * Grabe et al. (2000) Grabe, M., Wang, H., & Oster, G. (2000).  The mechanochemistry of v-atpase proton pumps.  Biophysical Journal, 78, 2798–2813.
  * Gu et al. (2020) Gu, J., Han, B., & Wang, J. (2020).  Covid-19: Gastrointestinal manifestations and potential fecal–oral transmission.  Gastroenterology, March 3, in press.
  * Kanberg et al. (2020a) Kanberg, N., Ashton, N. J., Andersson, L.-M., Yilmaz, A., Lindh, M., Nilsson, S., Price, R. W., Blennow, K., Zetterberg, H., & Gisslén, M. (2020a).  Neurochemical evidence of astrocytic and neuronal injury commonly found in covid-19.  Neurology, 95.
  * Kanberg et al. (2020b) Kanberg, N., Ashton, N. J., Andersson, L.-M., Yilmaz, A., Lindh, M., Nilsson, S., Price, R. W., Blennow, K., Zetterberg, H., & Gisslén, M. (2020b).  Neurochemical evidence of astrocytic and neuronal injury commonly found in covid-19.  Neurology, 95.
  * Kandimalla et al. (2020) Kandimalla, R., John, A., Abburi, C., Vallamkondu, J., & Reddy, P. H. (2020).  Current status of multiple drug molecules, and vaccines: An update in sars-cov-2 therapeutics.  Molecular Neurobiology, 15 July 2020, 1481–1517.
  * Katchalsky & Currant (1965) Katchalsky, A., & Currant, P. F. (1965).  Nonequilibrium Thermodynamics in Biophysics.  Boston: Harvard University Press.
  * Lucia (2015a) Lucia, U. (2015a).  Bioengineering thermodynamics: An engineering science for thermodynamics of biosystems.  International Journal of Thermodynamics, 18, 254–265.
  * Lucia (2015b) Lucia, U. (2015b).  Bioengineering thermodynamics of biological cells.  Theoretical Biology and Medical Modelling, 12, 29\.
  * Lucia & Grisolia (2017) Lucia, U., & Grisolia, G. (2017).  Second law efficiency for living cells.  Frontiers in Bioscience, 9, 270–275.
  * Lucia & Grisolia (2018a) Lucia, U., & Grisolia, G. (2018a).  Constructal law and ion transfer in normal and cancer cells.  Proceedings of the Romanian Academy Series A, 19, 213–218.
  * Lucia & Grisolia (2018b) Lucia, U., & Grisolia, G. (2018b).  Cyanobacteria and microalgae: Thermoeconomic considerations in biofuel production.  Energies, 11, 156.
  * Lucia & Grisolia (2020a) Lucia, U., & Grisolia, G. (2020a).  How life works — a continuous seebeck-peltier transition in cell membrane?  Entropy, 22, 960.
  * Lucia & Grisolia (2020b) Lucia, U., & Grisolia, G. (2020b).  Non-equilibrium thermodynamic approach to ca2+-fluxes in cancer.  Applied Sciences, 10, 6737.
  * Lucia & Grisolia (2020c) Lucia, U., & Grisolia, G. (2020c).  Thermal physics and glaucoma: From thermodynamic to biophysical considerations to design future therapies.  Applied Sciences, accepted, in printing.
  * Lucia & Grisolia (2020d) Lucia, U., & Grisolia, G. (2020d).  Thermal resonance and cell behavior.  Entropy, 22, 774.
  * Lucia et al. (2017) Lucia, U., Grisolia, G., & Astori, M. R. (2017).  Constructal law analysis of Cl- transport in eyes aqueous humor.  Scientific Reports, 7, 6856.
  * Lucia et al. (2020a) Lucia, U., Grisolia, G., & Deisboeck, T. S. (2020a).  Seebeck-like effect in SARS-CoV-2 bio-thermodynamics.  Atti della Accademia Peloritana dei Pericolanti, 98, A6.
  * Lucia et al. (2020b) Lucia, U., Grisolia, G., & Deisboek, T. S. (2020b).  Alzheimer’s disease: A thermodynamic perspective.  Applied Sciences, 10, 7562.
  * Lucia et al. (2016) Lucia, U., Grisolia, G., Dolcino, D., Astori, M. R., Massa, E., & Ponzetto, A. (2016).  Constructal approach to bio-engineering: The ocular anterior chamber temperature.  Scientific Reports, 6, 31099.
  * Lucia et al. (2018) Lucia, U., Grisolia, G., Ponzetto, A., & Deisboeck, T. S. (2018).  Thermodynamic considerations on the role of heat and mass transfer in biochemical causes of carcinogenesis.  Physica A, 490, 1164–1170.
  * Lucia et al. (2014) Lucia, U., Ponzetto, A., & Deisboeck, T. S. (2014).  A thermo-physical analysis of the proton pump vacuolar-ATPase: The constructal approach.  Scientific Reports, 4, 1.
  * Machado et al. (2020) Machado, R. R. G., Glaser, T., Araujo, D. B., Petiz, L. L., Oliveira, D. B. L., Durigon, G. S., Leal, A. L., Pinho, J. R. R., Ferreira, L. C. S., Ulrich, H., Durigon, E. L., & Guzzo, C. R. (2020).  Hypertonic saline solution inhibits sars-cov-2 in vitro assay.  bioRxiv.
URL https://www.biorxiv.org/content/10.1101/2020.08.11.244996v3

  * Marshall (2020) Marshall, M. (2020).  How covid-19 can damage the brain.  Nature, 585, 342–343.
  * Matschke et al. (2020a) Matschke, J., Lütgehetmann, M., Hagel, C., Sperhake, J. P., Schröder, A. S., Edler, C., Mushumba, H., Fitzek, A., Allweiss, L., Dandri, M., Dottermusch, M., Heinemann, A., Pfefferle, S., Schwabenland, M., Magruder, D. S., Bonn, S., Prinz, M., Gerloff, C., Püschel, K., Krasemann, S., Aepfelbacher, M., & Glatzel, M. (2020a).  Neuropathology of patients with covid-19 in germany: a post- mortem case series.  The Lancet Neurology, 19.
  * Matschke et al. (2020b) Matschke, J., Lütgehetmann, M., Hagel, C., Sperhake, J. P., Schröder, A. S., Edler, C., Mushumba, H., Fitzek, A., Allweiss, L., Dandri, M., Dottermusch, M., Heinemann, A., Pfefferle, S., Schwabenland, M., Magruder, D. S., Bonn, S., Prinz, M., Gerloff, C., Püschel, K., Krasemann, S., Aepfelbacher, M., & Glatzel, M. (2020b).  Neuropathology of patients with covid-19 in germany: a post-mortem case series.  The Lancet Neurology, 19, 919––929.
  * Meinhardt et al. (2021) Meinhardt, J., Radke, J., Dittmayer, C., Franz, J., Thomas, C., Mothes, R., Laue, M., Schneider, J., Brünink, S., Greuel, S., Lehmann, M., Hassan, O., Aschman, T., Schumann, E., Chua, R. L., Conrad, C., Eils, R., Stenzel, W., Windgassen, M., Rößler, L., Goebel, H.-H., Gelderblom, H. R., Martin, H., Nitsche, A., Schulz-Schaeffer, W. J., Hakroush, S., Winkler, M. S., Tampe, B., Scheibe, F., Körtvélyessy, P., Reinhold, D., Siegmund, B., Kühl, A. A., Elezkurtaj, S., Horst, D., Oesterhelweg, L., Tsokos, M., Ingold-Heppner, B., Stadelmann, C., Drosten, C., Corman, V. M., Radbruch, H., & Heppner, F. L. (2021).  Olfactory transmucosal sars-cov-2 invasion as a port of central nervous system entry in individuals with covid-19.  Nature Nauroscience, 24.
  * Moreno-P et al. (2020) Moreno-P, O., Leon-Ramirez, J.-M., Fuertes-Kenneally, L., Andres, M. P. M., Garcia-Navarro, M., Ruiz-Torregrosa, P., Boix, V., Gil, J., & Merino, E. (2020).  Hypokalemia as a sensitive biomarker of disease severity and the requirement for invasive mechanical ventilation requirement in covid-19 pneumonia: A case series of 306 mediterranean patients.  International Journal of Infectious Diseases, 100, 449–454.
  * Nakanishi-Matsui et al. (2010) Nakanishi-Matsui, M., Sekiya, M., & Futai, R. K. N. M. (2010).  The mechanism of rotating proton pumping atpases.  Biochimica et Biophysica Acta - Bioenergetics, 1797, 1343–1352.
  * Paniz-Mondolfi et al. (2020) Paniz-Mondolfi, A., Bryce, C., Grimes, Z., Gordon, R. E., Reidy, J., Lednicky, J., Sordillo, E. M., & Fowkes, M. (2020).  Central nervous system involvement by severe acute respiratory syndrome coronavirus-2 (sars-cov-2).  Journal of Medical Virology, 92.
  * Rudolph et al. (2006) Rudolph, M. G., Stanfield, R. L., & Wilson, I. A. (2006).  How TCRs bind MHCs, peptides, and coreceptors.  Annual Review of Immunology, 24, 419–466.
  * Solomon (2021) Solomon, T. (2021).  Neurological infection with sars- cov-2 – the story so far.  Nature Review Neurology.
  * Song et al. (2021) Song, E., Zhang, C., Israelow, B., Lu-Culligan, A., Prado, A. V., Skriabine, S., Lu, P., Weizman, O.-E., Liu, F., Dai, Y., Szigeti-Buck, K., Yasumoto, Y., Wang, G., Castaldi, C., Heltke, J., Ng, E., Wheeler, J., Alfajaro, M. M., Levavasseur, E., Fontes, B., Ravindra, N. G., Dijk, D. V., Mane, S., Gunel, M., Ring, A., Kazmi, S. A. J., Zhang, K., Wilen, C. B., Horvath, T. L., Plu, I., Haik, S., Thomas, J.-L., Louvi, A., Farhadian, S. F., Huttner, A., Seilhean, D., Renier, N., Bilguvar, K., & Iwasaki, A. (2021).  Neuroinvasion of sars-cov-2 in human andmouse brain.  Journal of Experimental Medicine, 218, e2020213.
  * Stevens & Forgac (1997) Stevens, T. H., & Forgac, M. (1997).  Structure, function and regulation of the vacuolar (H+)-atpase.  Annual Review of Cell and Developmental Biology, 13, 779–808.
  * Strong (2002) Strong, R. K. (2002).  Asymmetric ligand recognition by the activating natural killer cell receptor NKG2D, a symmetric homodimer.  Molecular Immunology, 38, 1029–1037.
  * Tuszynski & Kurzynski (2003) Tuszynski, J. A., & Kurzynski, M. (2003).  Introduction to Molecular Biophysics.  Boca Raton: CRC Press.
  * Yourgrau et al. (1982) Yourgrau, W., van der Merwe, A., & Raw, G. (1982).  Treatise on Irreversible and Statistical Thermophysics.  New York: Dover.
  * Zhou et al. (2020) Zhou, X., Chen, D., Wang, L., Zhao, Y., Wei, L., Chen, Z., & Yang, B. (2020).  Low serum calcium: a new, important indicator of covid-19 patients from mild/moderate to severe/critical.  Bioscience Reports, 40, BSR20202690.