current address: ]Max Planck Institute for the Physics of Complex Systems,
Nöthnitzer Straße 38, 01187 Dresden, Germany

# Turing’s diffusive threshold in random reaction-diffusion systems

Pierre A. Haas<EMAIL_ADDRESS>[  Mathematical Institute, University of
Oxford, Woodstock Road, Oxford OX2 6GG, United Kingdom  Raymond E. Goldstein
<EMAIL_ADDRESS>Department of Applied Mathematics and
Theoretical Physics, Centre for Mathematical Sciences,
University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom

###### Abstract

Turing instabilities of reaction-diffusion systems can only arise if the
diffusivities of the chemical species are sufficiently different. This
threshold is unphysical in most systems with $N=2$ diffusing species, forcing
experimental realizations of the instability to rely on fluctuations or
additional nondiffusing species. Here we ask whether this diffusive threshold
lowers for $N>2$ to allow “true” Turing instabilities. Inspired by May’s
analysis of the stability of random ecological communities, we analyze the
probability distribution of the diffusive threshold in reaction-diffusion
systems defined by random matrices describing linearized dynamics near a
homogeneous fixed point. In the numerically tractable cases $N\leqslant 6$, we
find that the diffusive threshold becomes more likely to be smaller and
physical as $N$ increases and that most of these many-species instabilities
cannot be described by reduced models with fewer species.

In 1952, Turing described the pattern-forming instability that now bears his
name [1]: diffusion can destabilize a fixed point of a system of reactions
that is stable in well-mixed conditions. Nigh on threescore and ten years on,
the contribution of Turing’s mechanism to chemical and biological
morphogenesis remains debated, not least because of the _diffusive threshold_
inherent in the mechanism: chemical species in reaction systems are expected
to have roughly equal diffusivities, yet Turing instabilities cannot arise at
equal diffusivities [2, 3]. It remains an open problem to determine how much
of a diffusivity difference is required for generic systems to undergo this
instability, yet this diffusive threshold has been recognized at least since
reduced models of the Belousov–Zhabotinsky reaction [4, 5] only produced
Turing patterns at unphysically large diffusivity differences.

For this reason, the first experimental realizations of Turing instabilities
[6, 7, 8] obviated the threshold by using gel reactors that greatly reduced
the effective diffusivity of one species [9, 10]. (Analogously, biological
membrane transport dynamics can increase the effective diffusivity difference
[11].) Later work showed how binding to an immobile substrate, or more
generally, a third, nondiffusing species, can allow Turing instabilities even
if the $N=2$ diffusing species have equal diffusivities [12, 13, 14]. Such
nondiffusing species continue to permeate more recent work on the network
topology of Turing systems [15, 16].

Moreover, Turing instabilities need not be deterministic: fluctuation-driven
instabilities in reaction-diffusion systems have noise-amplifying properties
that allow their pattern amplitude to be comparable to that of deterministic
Turing patterns [17], with a lower diffusive threshold than the deterministic
one [18, 19, 20, 21]. A synthetic bacterial population with $N=2$ species that
exhibits patterns in agreement with such a stochastic Turing instability, but
does not satisfy the conditions for a deterministic instability [22], was
reported recently.

These experimental instabilities relying on fluctuations or the dynamics of
additional nondiffusing species and the nonlinear instabilities arising from
finite-amplitude perturbations [2] are not however instabilities in Turing’s
own image. Can such instabilities be realized instead in systems with $N>2$
diffusing species? Equivalently, is the diffusive threshold lower in such
systems? These questions have remained unanswered, perhaps because, in
contrast to the textbook case $N=2$ and the concomitant picture of an
“inhibitor” out-diffusing an “activator” [23, *[][Vol.I, AppendixB.1, pp.507--
509andVol.II, Chap.2, pp.71--140, 3rded.]murray], the complicated instability
conditions for $N>2$ [25] do not lend themselves to analytical progress.

Here, we analyze the diffusive threshold for Turing instabilities with
$2\leqslant N\leqslant 6$ diffusing species. Inspired by May’s work on the
stability of random ecological communities [26], we analyze _random Turing
instabilities_ by sampling random matrices that represent the linearized
reaction dynamics of otherwise unspecified reaction-diffusion systems. A
semianalytic approach shows that the diffusive threshold is more likely to be
smaller and physical for $N=3$ compared to $N=2$, and that two of the three
diffusivities are equal at the transition to instability. We extend these
results to the remaining numerically tractable cases of reaction-diffusion
systems with $4\leqslant N\leqslant 6$ and two different diffusivities: their
Turing instabilities are still more likely to have a smaller and physical
diffusive threshold, but most of them cannot be described by reduced models
with fewer species.

We begin with the simplest case, $N=2$, in which species $u$, $v$ obey

$\displaystyle\dot{u}=f(u,v)+d_{u}\nabla^{2}u,$
$\displaystyle\dot{v}=g(u,v)+d_{v}\nabla^{2}v.$ (1)

The conditions for Turing instability in this system [24] only depend on the
four entries of the Jacobian

$\displaystyle\mathsf{J}=\left(\begin{array}[]{cc}f_{u}&f_{v}\\\
g_{u}&g_{v}\end{array}\right),$ (4)

the partial derivatives of the reaction system at a fixed point
$(u_{\ast},v_{\ast})$ of the homogeneous system. This fixed point is stable to
homogeneous perturbations iff ${J\equiv\det{\mathsf{J}}>0}$ and
${I_{1}\equiv\operatorname{tr}{\mathsf{J}}<0}$. A stable fixed point of this
kind is unstable to a Turing instability only if $p\equiv-f_{u}g_{v}>0$ [24].
Defining the diffusion coefficient ratio
$D_{2}=\max{\\{d_{u}/d_{v},d_{v}/d_{u}\\}}\geqslant 1$, a Turing instability
occurs iff these conditions hold along with 111See Supplemental Material at
[url to be inserted], which includes Refs. [3, 24, 28, 29, 30, 31, 32], for
details of calculations for $N=2$, the derivation of the semianalytic approach
for $N>2$ and a discussion of its numerical implementation, statistics of the
wavenumber of Turing instabilities, a discussion of “slow” species, the proof
of an asymptotic result, and for python3 code.

$\displaystyle D_{2}\geqslant
D_{2}^{\ast}\equiv\left(\dfrac{\sqrt{J\vphantom{p}}+\sqrt{J+p}}{\min{\\{|f_{u}|,|g_{v}|\\}}}\right)^{2}.$
(5)

Figure 1: Turing’s diffusive threshold for $N=2$. (a) Cartoon of the diffusive
threshold and the fine-tuning (FT) problem for $R\approx 1$ and $R\gg 1$. The
diffusivity difference required mathematically is unphysical in the hatched
region ${\mathcal{D}\leqslant D_{2}^{\ast}\leqslant D_{2}^{\max}}$. (b)
Distribution $P(D_{2}^{\ast})$, supported on the (scaled) interval
$[1,D_{2}^{\max}(R)]$, estimated for different $R$. (c) Plot of
$\mathbb{P}(D_{2}^{\ast}<\mathcal{D})$ [shaded areas in panels (a) and (b)]
against $R$, revealing the diffusive threshold. Markers: estimates from panel
(b); solid line: exact result [27] for ${R>\mathcal{D}}$ [33].

This diffusivity difference $D_{2}^{\ast}$, which is required mathematically
for instability, is unphysical [Fig. 1(a)] if it exceeds the diffusivity
difference $\mathcal{D}>1$ of the physical system: $\mathcal{D}\approx 1$ for
similarly sized molecules in solution, but, e.g., $\mathcal{D}\approx 20$ for
the stochastic Turing instability observed in Ref. [22]. Hereinafter we take
$\mathcal{D}=5$ arbitrarily (but have checked that the value of $\mathcal{D}$
does not affect results qualitatively). To quantify $D_{2}^{\ast}$, we
introduce the _range_ $R$ of kinetic parameters,

$\displaystyle
R\equiv\frac{\max{\\{|f_{u}|,|f_{v}|,|g_{u}|,|g_{v}|\\}}}{\min{\\{|f_{u}|,|f_{v}|,|g_{u}|,|g_{v}|\\}}}.$
(6)

Equivalently, $f_{u},f_{v},g_{u},g_{v}\in I\equiv[-R,-1]\cup[1,R]$ up to
scaling, with one parameter equal to $\pm 1$ and one equal to $\pm R$. One
deduces [27] that

$\displaystyle D_{2}^{\ast}\leqslant
D_{2}^{\max}(R)\equiv\left(R+\sqrt{R^{2}-1}\right)^{2},$ (7)

as shown in Fig. 1(a). Note that $D_{2}^{\max}\\!\to\\!1$ as $R\\!\to\\!1$; in
this limit, there is no diffusive threshold: $R\approx 1$ is a particular
instance of the converse _fine-tuning problem_ for the reaction kinetics that
allows Turing instabilities at nearly equal diffusivities more generally [3].
If $R\gg 1$, then $D_{2}^{\max}=O\bigl{(}R^{2}\bigr{)}$. This does not imply
the existence of a threshold, for it does not preclude most systems with range
$R$ having $D_{2}^{\ast}\ll D_{2}^{\max}$. The existence of a diffusive
threshold therefore relates to the distribution of $D_{2}^{\ast}$ for systems
with range $R$.

To understand this distribution, we draw inspiration from May’s statistical
analysis of the stability of ecological communities [26], which studies random
Jacobians, corresponding to equilibria of otherwise unspecified population
dynamics. By analogy, we study random Turing instabilities, sampling uniformly
and independently random Jacobians corresponding to equilibria of otherwise
unspecified reaction kinetics, and analyze the criteria for them to be Turing
unstable. There is of course no more reason to expect the kinetic parameters
to be independent or uniformly distributed than there is reason to expect the
linearized population dynamics in May’s analysis [26] to be independent or
normally distributed. Yet, in the absence of experimental understanding of
what the distributions of these parameters should be (in either context), the
potential of the random matrix approach to reveal stability principles has
been amply demonstrated in population dynamics [34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44].

We sample the kinetic parameters in Eq. (4) independently and uniformly from
$I$, set one of them equal to $\pm 1$ and one equal to $\pm R$, and thus
estimate the probability distribution $P(D_{2}^{\ast})$ for fixed $R$ [Fig.
1(b)]. The threshold is quantified by the probability of a Turing instability
being physical,

${\mathbb{P}(D_{2}^{\ast}<\mathcal{D})}=\int_{1}^{\mathcal{D}}{P(D_{2}^{\ast})\,\mathrm{d}D_{2}^{\ast}}.$
(8)

Both from the estimates in Fig. 1(b) and by evaluating the integral in closed
form [27], we find that ${\mathbb{P}(D_{2}^{\ast}<\mathcal{D})}$ is tiny [Fig.
1(c)], except if $R$ is small, which is the fine-tuning problem again. In
other words, the required diffusivity difference is very likely to be
unphysical. This expresses Turing’s diffusive threshold for $N=2$.

To investigate how this threshold changes with $N$ we consider next the $N=3$
system

$\displaystyle\dot{u}$ $\displaystyle=f(u,v,w)+d_{u}\nabla^{2}u,$ (9a)
$\displaystyle\dot{v}$ $\displaystyle=g(u,v,w)+d_{v}\nabla^{2}v,$ (9b)
$\displaystyle\dot{w}$ $\displaystyle=h(u,v,w)+\nabla^{2}w,$ (9c)

where we have rescaled space to set $d_{w}=1$. We introduce the matrix of
diffusivities and the reaction Jacobian,

$\displaystyle\mathsf{D}=\left(\begin{array}[]{ccc}d_{u}&0&0\\\ 0&d_{v}&0\\\
0&0&1\end{array}\right),$
$\displaystyle\mathsf{J}=\left(\begin{array}[]{ccc}f_{u}&f_{v}&f_{w}\\\
g_{u}&g_{v}&g_{w}\\\ h_{u}&h_{v}&h_{w}\end{array}\right),$ (16)

in which the entries of $\mathsf{J}$ are again the partial derivatives
evaluated at a fixed point $(u_{\ast},v_{\ast},w_{\ast})$ of the homogeneous
system. This fixed point is unstable to a Turing instability if it is stable
but, for some eigenvalue $-k^{2}<0$ of the Laplacian,
$\smash{\overline{\mathsf{J}}}\bigl{(}k^{2}\bigr{)}=\mathsf{J}-k^{2}\mathsf{D}$
is unstable [3], i.e. has an eigenvalue $\lambda$ such that
$\operatorname{Re}(\lambda)<0$. More precisely, a Turing instability arises
when a real eigenvalue of $\smash{\overline{\mathsf{J}}\bigl{(}k^{2}\bigr{)}}$
crosses zero, i.e. when
$\mathcal{J}\smash{\bigl{(}k^{2}\bigr{)}}\equiv\det{\smash{\overline{\mathsf{J}}\bigl{(}k^{2}\bigr{)}}}=0$,
and therefore arises first at a wavenumber $k=k_{\ast}$ with
$\mathcal{J}\bigl{(}k_{\ast}^{2}\bigr{)}=\partial\mathcal{J}/\partial
k^{2}\bigl{(}k_{\ast}^{2}\bigr{)}=0$ [3]. Hence $\mathcal{J}$, a cubic
polynomial in $k^{2}$, has a double root at $k^{2}=k_{\ast}^{2}>0$, so its
discriminant [29] vanishes. This discriminant, $\Delta(d_{u},d_{v})$, is a
polynomial in $d_{u},d_{v}$. We denote by $K(d_{u},d_{v})$ the double root of
$\mathcal{J}$ corresponding to a point $(d_{u},d_{v})$ on the curve
$\Delta(d_{u},d_{v})=0$.

Figure 2: Results for $N=3$. (a) Contours of $D_{3}(d_{u},d_{v})$ in the
positive $(d_{u},d_{v})$ quadrant. (b) Smoothed distribution
$P(D_{3}^{\ast})$, estimated for different $R$. Inset: same plot, scaled to
$[1,D_{2}^{\max}(R)]$ for comparison to $N=2$ in Fig. 1(a). (c)
$\mathbb{P}(D_{N}^{\ast}<\mathcal{D})$ against $R$ for ${N\in\\{2,3\\}}$: the
diffusive threshold lowers for $N=3$ compared to $N=2$. (d) Proportion
$\phi_{N}(\mathcal{D})$ of random Jacobians that have a physical Turing
instability, plotted against $R$, for $N\in\\{2,3\\}$. Inset: proportion
$\tau_{N}$ of random Jacobians that have a (physical or unphysical) Turing
instability, averaged over $R$, for $N\in\\{2,3\\}$ [33].

Determining the diffusive threshold for Turing instability in Eqs. (9) thus
requires solving the problem

$\displaystyle\text{minimize }D_{3}(d_{u},d_{v})\quad\text{subject
to}\left\\{\begin{array}[]{ll}\Delta(d_{u},d_{v})=0,\\\
K(d_{u},d_{v})>0,\end{array}\right.$ (19)

in which the diffusion coefficient ratio is

$\displaystyle
D_{3}(d_{u},d_{v})\\!=\\!\max{\\!\\{d_{u},1/d_{u},d_{v},1/d_{v},d_{u}/d_{v},d_{v}/d_{u}\\}}.$
(20)

With the aim in mind of obtaining statistics for the minimal value
$D_{3}^{\ast}$, direct numerical solution of this constrained optimization
problem is obviously not a feasible approach. In the Supplemental Material
[27], we therefore show how to reduce solving problem (19) to polynomial root
finding. This semianalytic approach reveals a particular class of minima,
attained at the vertices of the contours of $D_{3}(d_{u},d_{v})$ [Fig. 2(a)],
i.e. at $d_{u}=1$, $d_{v}=1$, or $d_{u}=d_{v}$. In these cases,
$\Delta(d_{u},d_{v})=0$ is a (sextic) polynomial in the single variable
$d_{v}$, $d_{u}$, or $d=d_{u}=d_{v}$, respectively. We call these minima
“binary”, since the corresponding systems have only two different
diffusivities. We implement this approach numerically [27], and sample random
systems similarly to the case $N=2$, drawing the entries of $\mathsf{J}$ in
Eq. (16) uniformly and independently at fixed range $R$.

Remarkably, all global minima we found numerically were binary [27]. This
means that the minimizing systems come in two flavors: those with two “fast”
diffusers and one “slow” diffuser, and those with one “fast” diffuser and two
“slow” diffusers. Systems with a nondiffusing species are a limit of the
former; this point will be discussed below. The latter arise in models of
scale pattern formation in fish and lizards [45, 46], in which short-range
pigments respectively activate and inhibit a long-range factor.

The distribution of $D_{3}^{\ast}$, shown for different values of $R$ in Fig.
2(b), has a different shape from that of $D_{2}^{\ast}$ [Figs. 1(a) and 2(b),
inset]. While the support of the distribution of $D_{3}^{\ast}$ does not
appear to be bounded, Fig. 2(c) shows that
${\mathbb{P}(D_{3}^{\ast}<\mathcal{D})>\mathbb{P}(D_{2}^{\ast}<\mathcal{D})}$.
Hence the diffusivity difference is more likely to be physical for $N=3$ than
for $N=2$: the diffusive threshold is lowered.

The proportion $\tau_{N}$ of random kinetic Jacobians that have a Turing
instability (be it physical or unphysical) is smaller for $N=3$ than for $N=2$
[Fig. 2(d), inset]. This is not surprising, because a random Jacobian is less
likely to correspond to a stable fixed point (which, we recall, is a necessary
condition for Turing instability) for $N=3$ than for $N=2$, essentially
because its entries have to satisfy more conditions for stability if $N=3$. It
is therefore striking that the threshold is reduced sufficiently for $N=3$
compared to $N=2$ for the proportion
$\phi_{N}(\mathcal{D})=\tau_{N}\mathbb{P}(D_{N}^{\ast}<\mathcal{D})$ of random
Jacobians that have a physical Turing instability to be larger for $N=3$ than
for $N=2$ [Fig. 2(d)], even though a Turing instability of any kind is more
likely if $N=2$.

To extend these results to $N>3$ diffusing species, we consider the
(linearized) reaction-diffusion system

$\boldsymbol{\dot{u}}=\mathsf{J}\cdot\boldsymbol{u}+\mathsf{D}\cdot\nabla^{2}\boldsymbol{u},$
(21)

where $\mathsf{J}$ is a random kinetic Jacobian, and $\mathsf{D}$ is a
diagonal matrix of diffusivities. Even with our semianalytic approach, this
cannot be analyzed for general $\mathsf{D}$: not even for $N=4$ were we able
to obtain closed forms of the required polynomials. To make further progress,
we therefore restrict to binary $\mathsf{D}$ in which the $N$ diffusivities
take two different values only, since we showed above that $D_{3}^{\ast}$ is
attained for such binary $\mathsf{D}$. As in the case $N=3$, this reduces the
discriminant condition ${\Delta(\mathsf{D})=0}$ to polynomial equations in one
variable that determine the minimum diffusivity difference $D_{N}^{\ast}$ for
Turing instability in these binary systems [27].

Figure 3: Results for “binary” systems with $4\leqslant N\leqslant 6$. (a)
${\mathbb{P}(D_{N}^{\ast}<\mathcal{D})}$ against $R$ for ${3\leqslant
N\leqslant 6}$, revealing further lowering of the diffusive threshold compared
to the case $N=3$. (b) Proportion $\sigma_{N}$ of random stable kinetic
Jacobian that have a (binary, if $N>3$) Turing instability, averaged over $R$,
and plotted against $N$. (c) Proportion $\phi_{N}(\mathcal{D})$ of random
Jacobians that have a physical Turing instability plotted against $R$, for
$3\leqslant N\leqslant 6$. Inset: proportion $\tau_{N}$ of random Jacobians
that have a (physical or unphysical) Turing instability, averaged over $R$,
for $3\leqslant N\leqslant 6$ [33].

Figure 3(a) shows that the diffusive threshold lowers further for $4\leqslant
N\leqslant 6$ in these systems. At the same time, the fact that most stable
random kinetic Jacobians undergo such a binary Turing instability [Fig. 3(b)]
suggests that these provide a useful picture of the diffusive threshold.
However, $\tau_{N}$ decreases further for $N\geqslant 4$ [Fig. 3(c), inset],
and the widening of the bottleneck is not sufficient to prevent
$\phi_{N}(\mathcal{D})$ from decreasing for $N\geqslant 4$. Nonetheless, since
both $\mathbb{P}(D_{N}^{\ast}<\mathcal{D})$ and the proportion $\sigma_{N}$ of
stable random Jacobians that are Turing unstable increase [Figs. 3(a) and
3(b)], so does the proportion of stable random Jacobians that have a physical
Turing instability.

How then to realize “true” Turing instabilities experimentally? Our analysis
shows that the diffusive threshold of a Turing instability is more likely to
be physical the more species there are, but how to find an experimental Turing
instability in the first place? Turing instabilities remain rare in random
reaction systems even as the number of species is increased, but the above
shows that this rareness mainly results from the rareness of stable equilibria
in such systems. The proverbial search for the needle in a haystack can
therefore be avoided by exploring biochemical systems that admit a stable
equilibrium, and evolving them towards a “true” Turing instability.

This analysis does not however reveal whether these instabilities lead to
patterns that are observable at the physical scale of the system. Analysis of
the wavenumber at which the linear instability first arises [27] suggests that
we can extend our conclusions: Turing instabilities with more species are more
likely to have physical diffusivity differences and to be observable. However,
our statistical, linearized analysis cannot fully answer this question of
observability, because it fundamentally depends on the system through details
of the nonlinearities of its reaction kinetics, which set the precise nature
and scale of the Turing patterns that develop beyond onset of the instability;
this is why we have relegated this discussion to the Supplemental Material
[27].

The different species in the systems with $3\leqslant N\leqslant 6$ analyzed
above separate into “fast” and “slow” diffusers. The diffusion of these “slow”
species is often ignored in the analysis of systems of many chemical reactions
[30], such as the full Belousov–Zhabotinsky reaction [47]. Corresponding
reduced models are obtained by substituting the steady-state kinetics of the
“slow” species into the remaining equations, thereby eliminating them from the
system [30]. The conditions for Turing instability in these reduced models are
(almost) equivalent to those for the full model with nondiffusing “slow”
species [30]. However, the diffusion of the “slow” species cannot in general
be ignored: up to reordering species and rescaling space,

$\displaystyle\mathsf{D}=\left(\begin{array}[]{c|c}\mathsf{I}&\mathsf{0}\\\
\hline\cr\mathsf{0}&d\mathsf{I}\end{array}\right),$
$\displaystyle\mathsf{J}=\left(\begin{array}[]{c|c}\mathsf{J_{11}}&\mathsf{J_{12}}\\\
\hline\cr\mathsf{J_{21}}&\mathsf{J_{22}}\end{array}\right),$ (26)

where $d<1$ is the common diffusivity of the slow diffusers. Results of Ref.
[30] imply that there is a Turing instability with nondiffusing “slow”
species, i.e. with $d=0$, only if
$\mathsf{J_{11}}-\mathsf{J_{12}}\mathsf{J}^{-1}_{\mathsf{22}}\mathsf{J_{21}}$
has a positive (real) eigenvalue [27]. Although the proportion of Turing
unstable systems that have $n\geqslant 2$ fast diffusers (and hence could _a
priori_ still undergo a Turing instability with $d=0$) is large [Fig. 4(a)],
the proportion of systems that do undergo such an instability is small, even
if we restrict to those systems with physical diffusivity differences [Fig.
4(b)]. Hence most of these Turing instabilities with $N>2$ species require all
species to diffuse.

Figure 4: “Slow” diffusers in binary Turing instabilities with $3\leqslant
N\leqslant 6$. (a) Proportion $f_{N}$ of Turing unstable systems with
$n\geqslant 2$ “fast” diffusers plotted against $N$, averaged over $R$. (b)
Proportion $r_{N}$ of systems that remain Turing unstable at $d=0$, plotted
against $N$, averaged over $R$. Closed markers: all Turing systems; open
markers: physical Turing systems with $D_{N}^{\ast}<\mathcal{D}$ [33].

These many-species Turing instabilities, although binary, are thus more
general than the instabilities of systems with nondiffusing species realized
experimentally in gel reactors [6, 7, 8] and analyzed theoretically in Ref.
[30]. In particular, this shows that reduced models give but an incomplete
picture of Turing instabilities. Together with our main result, that the
diffusive threshold lowers as $N$ increases, this implies that the failure of
a reduced model to produce a physical Turing instability cannot be taken as an
indication that a Turing instability cannot exist in the full system that the
reduced model seeks to describe.

In this Letter, we have analyzed random Turing instabilities to show how the
diffusive threshold that has hampered experimental efforts to generate “true”
Turing instabilities in systems of $N=2$ diffusing species lowers for systems
with $N\geqslant 3$, most of whose instabilities cannot be described by
reduced models with fewer species. All of this does not, however, explain the
existence of a “large” threshold in the first place: even though Turing
instabilities at equal diffusivities are impossible [2, 3], this does not mean
that the threshold needs to be “large”. In this context, we prove an
asymptotic result in the Supplemental Material [27]: for a Jacobian
$\mathsf{J}$ to allow a Turing instability at almost equal diffusivities
$\mathsf{D}\approx\mathsf{I}$, $\mathsf{J}$ must be even closer to a singular
matrix $\mathsf{J_{0}}$, i.e.
$\mathsf{J}-\mathsf{J_{0}}\ll\mathsf{D}-\mathsf{I}$. In this sense, the
threshold $\mathsf{D}-\mathsf{I}$ is asymptotically “large”. Understanding how
a large threshold arises more generally outside this asymptotic regime and
lowers as $N$ increases remains an open problem, as do extending the present
analysis to include the nonlocal interactions [48, 49] that arise for example
in vegetation patterns [50] and extending previous work [16, 51] on the
robustness of Turing patterns to $N\geqslant 3$. The latter in particular may
help to identify those chemical or biological pattern forming systems with
$N\geqslant 3$ in which the “true” Turing instabilities discussed here can be
realized experimentally.

###### Acknowledgements.

We thank N. Goldenfeld, A. Krause, and P. K. Maini for discussions. This work
was supported in part by a Nevile Research Fellowship from Magdalene College,
Cambridge, and a Hooke Research Fellowship (P.A.H.), Established Career
Fellowship EP/M017982/1 from the Engineering and Physical Sciences Research
Council and Grant 7523 from the Marine Microbiology Initiative of the Gordon
and Betty Moore Foundation (R.E.G.).

## References

  * Turing [1952] A. M. Turing, The chemical basis of morphogenesis, Phil. Trans. R. Soc. B 237, 37 (1952).
  * Vastano  _et al._ [1987] J. A. Vastano, J. E. Pearson, W. Horsthemke, and H. L. Swinney, Chemical pattern formation with equal diffusion coefficients, Phys. Lett. A 124, 320 (1987).
  * Pearson and Horsthemke [1989] J. E. Pearson and W. Horsthemke, Turing instabilities with nearly equal diffusion coefficients, J. Chem. Phys. 90, 1588 (1989).
  * Becker and Field [1985] P. K. Becker and R. J. Field, Stationary concentration patterns in the oregonator model of the Belousov–Zhabotinskii reaction, J. Phys. Chem. 89, 118 (1985).
  * Rovinskii [1987] A. B. Rovinskii, Turing bifurcation and stationary patterns in the ferroin-catalyzed Belousov–Zhabotinskii reaction, J. Phys. Chem. 91, 4606 (1987).
  * Castets  _et al._ [1990] V. Castets, E. Dulos, J. Boissonade, and P. De Kepper, Experimental evidence of a sustained standing Turing-type nonequilibrium chemical pattern, Phys. Rev. Lett. 64, 2953 (1990).
  * De Kepper  _et al._ [1991] P. De Kepper, V. Castets, E. Dulos, and J. Boissonade, Turing-type chemical patterns in the chlorite-iodide-malonic acid reaction, Physica D 49, 161 (1991).
  * Ouyang and Swinney [1991] Q. Ouyang and H. L. Swinney, Transition from a uniform state to hexagonal and striped Turing patterns, Nature (London) 352, 610 (1991).
  * Lengyel and Epstein [1991] I. Lengyel and I. R. Epstein, Modeling of Turing structures in the chloride-iodid-malonic acid-starch reaction system, Science 251, 650 (1991).
  * Dulos  _et al._ [1996] E. Dulos, J. Boissonade, J. J. Perraud, B. Rudovics, and P. De Kepper, Chemical morphogenesis: Turing patterns in an experimental chemical system, Acta Biother. 44, 249 (1996).
  * Recho  _et al._ [2019] P. Recho, A. Hallou, and E. Hannezo, Theory of mechanochemical patterning in biphasic biological tissues, Proc. Natl. Acad. Sci. USA 116, 5344 (2019).
  * Lengyel and Epstein [1992] I. Lengyel and I. R. Epstein, A chemical approach to designing Turing patterns in reaction-diffusion systems, Proc. Natl. Acad. Sci. USA 89, 3977 (1992).
  * Pearson [1992] J. E. Pearson, Pattern formation in a $(2+1)$-species activator-inhibitor-immobilizer system, Physica A 188, 178 (1992).
  * Korvasová  _et al._ [2015] K. Korvasová, E. A. Gaffney, P. K. Maini, M. A. Ferreira, and V. Klika, Investigating the Turing conditions for diffusion-driven instability in the presence of a binding immobile substrate, J. Theor. Biol. 367, 286 (2015).
  * Marcon  _et al._ [2016] L. Marcon, X. Diego, J. Sharpe, and P. Müller, High-throughput mathematical analysis identifies Turing networks for patterning with equally diffusing signals, eLife 5, e14022 (2016).
  * Diego  _et al._ [2018] X. Diego, L. Marcon, P. Müller, and J. Sharpe, Key features of Turing systems are determined purely by network topology, Phys. Rev. X 8, 021071 (2018).
  * Biancalani  _et al._ [2017] T. Biancalani, F. Jafarpour, and N. Goldenfeld, Giant amplification of noise in fluctuation-induced pattern formation, Phys. Rev. Lett. 118, 018101 (2017).
  * Butler and Goldenfeld [2009] T. Butler and N. Goldenfeld, Robust ecological pattern formation induced by demographic noise, Phys. Rev. E 80, 030902(R) (2009).
  * Biancalani  _et al._ [2010] T. Biancalani, D. Fanelli, and F. Di Patti, Stochastic Turing patterns in the Brusselator model, Phys. Rev. E 81, 046215 (2010).
  * Butler and Goldenfeld [2011] T. Butler and N. Goldenfeld, Fluctuation-driven Turing patterns, Phys. Rev. E 84, 011112 (2011).
  * Di Patti  _et al._ [2018] F. Di Patti, L. Lavacchi, R. Arbel-Goren, L. Schein-Lubomirsky, D. Fanelli, and J. Stavans, Robust stochastic turing patterns in the development of a one-dimensional cyanobacterial organism, PLoS Biol. 16, e2004877 (2018).
  * Karig  _et al._ [2018] D. Karig, K. M. Martini, T. Lu, N. A. DeLateur, N. Goldenfeld, and R. Weiss, Stochastic Turing patterns in a synthetic bacterial population, Proc. Natl. Acad. Sci. USA 115, 6572 (2018).
  * Koch and Meinhardt [1994] A. J. Koch and H. Meinhardt, Biological pattern formation: from basic mechanisms to complex structures, Rev. Mod. Phys. 66, 1481 (1994).
  * Murray [2002] J. D. Murray, in  _Mathematical Biology_ (Springer, Berlin, Germany, 2002).
  * Satnoianu  _et al._ [2000] R. A. Satnoianu, M. Menzinger, and P. K. Maini, Turing instabilities in general systems, J. Math. Biol. 41, 493 (2000).
  * May [1972] R. M. May, Will a large complex ecosystem be stable?, Nature (London) 238, 413 (1972).
  * Note [1] See Supplemental Material at [url to be inserted], which includes Refs. [3, 24, 28, 29, 30, 31, 32], for details of calculations for $N=2$, the derivation of the semianalytic approach for $N>2$ and a discussion of its numerical implementation, statistics of the wavenumber of Turing instabilities, a discussion of “slow” species, the proof of an asymptotic result, and for python3 code.
  * Johansson [2018] F. Johansson,  _mpmath : a python library for arbitrary-precision floating-point arithmetic_ (version 1.1.0, 2018).
  * Yap [2000] C.-K. Yap, in  _Fundamental Problems in Algorithmic Algebra_ (Oxford University Press, Oxford, England, 2000) Chap. 6 and Chap. 7, pp. 141–218.
  * Smith and Dalchau [2018] S. Smith and N. Dalchau, Model reduction enables Turing instability analysis of large reaction–diffusion models, J. R. Soc. Interface 15, 20170805 (2018).
  * Hinch [1991] E. J. Hinch, in  _Perturbation Methods_ (Cambridge University Press, Cambridge, UK, 1991) Chap. 1.6, pp. 15–18.
  * Lidskii [1966] V. B. Lidskii, Perturbation theory of non-conjugate operators, USSR Comput. Math. Math. Phys. 6, 73 (1966).
  * [33] All results shown are for $\mathcal{D}=5$. The asymmetric error bars in Figs. 3 and 4 show $95\%$ confidence intervals larger than the plot markers, corrected for systems for which the numerics failed.
  * Allesina and Tang [2012] S. Allesina and S. Tang, Stability criteria for complex ecosystems, Nature (London) 483, 205 (2012).
  * Mougi and Kondoh [2012] A. Mougi and M. Kondoh, Diversity of interaction types and ecological community stability, Science 337, 349 (2012).
  * Coyte  _et al._ [2015] K. Z. Coyte, J. Schluter, and K. R. Foster, The ecology of the microbiome: Networks, competition, and stability, Science 350, 663 (2015).
  * Grilli  _et al._ [2016] J. Grilli, T. Rogers, and S. Allesina, Modularity and stability in ecological communities, Nat. Commun. 7, 12031 (2016).
  * Gibbs  _et al._ [2018] T. Gibbs, J. Grilli, T. Rogers, and S. Allesina, Effect of population abundances on the stability of large random ecosystems, Phys. Rev. E 98, 022410 (2018).
  * Serván  _et al._ [2018] C. A. Serván, J. A. Capitán, J. Grilli, K. E. Morrison, and S. Allesina, Coexistence of many species in random ecosystems, Nat. Ecol. Evol. 2, 1237 (2018).
  * Butler and O’Dwyer [2018] S. Butler and J. P. O’Dwyer, Stability criteria for complex microbial communities, Nat. Commun. 9, 2970 (2018).
  * Stone [2018] L. Stone, The feasibility and stability of large complex biological networks: a random matrix approach, Sci. Rep. 8, 8246 (2018).
  * Maynard  _et al._ [2019] D. S. Maynard, C. A. Serván, J. A. Capitán, and S. Allesina, Phenotypic variability promotes diversity and stability in competitive communities, Ecol. Lett. 22, 1776 (2019).
  * Haas  _et al._ [2020] P. A. Haas, N. M. Oliveira, and R. E. Goldstein, Subpopulations and stability in microbial communities, Phys. Rev. Research 2, 022036(R) (2020).
  * Barron and Galla [2020] J. W. Barron and T. Galla, Dispersal-induced instability in complex ecosystems, Nat. Commun. 11, 6032 (2020).
  * Nakamasu  _et al._ [2009] A. Nakamasu, G. Takahashi, A. Kanbe, and S. Kondo, Interactions between zebrafish pigment cells responsible for the generation of Turing patterns, Proc. Natl. Acad. Sci. USA 106, 8429 (2009).
  * Manukyan  _et al._ [2017] L. Manukyan, S. A. Montandon, A. Fofonjka, S. Smirnov, and M. C. Milinkovitch, A living mesoscopic cellular automaton made of skin scales, Nature (London) 544, 173 (2017).
  * Györgi  _et al._ [1990] L. Györgi, T. Turányi, and R. J. Field, Mechanistic details of the oscillator Belousov–Zhabotinskii reaction, J. Phys. Chem. 94, 7162 (1990).
  * Ohta  _et al._ [1989] T. Ohta, M. Mimura, and R. Kobayashi, Higher-dimensional localized patterns in excitable media, Physica D 34, 115 (1989).
  * Goldstein  _et al._ [1996] R. E. Goldstein, D. J. Muraki, and D. M. Petrich, Interface proliferation and the growth of labyrinths in a reaction-diffusion system, Phys. Rev. E 53, 3933 (1996).
  * Meron [2019] E. Meron, Vegetation pattern formation: The mechanisms behind the forms, Physics Today 72, 30 (2019).
  * Scholes  _et al._ [2019] N. S. Scholes, D. Schnoerr, M. Isalan, and M. P. H. Stumpf, A comprehensive network atlas reveals that Turing patterns are common but not robust, Cell Systems 9, 243 (2019).


# A general framework for probabilistic sensitivity analysis with respect to
distribution parameters

Jiannan Yang
Department of Engineering
University of Cambridge
Trumpington Street, Cambridge CB2 1PZ, UK
<EMAIL_ADDRESS>

(February 1, 2023)

###### Abstract

Probabilistic sensitivity analysis identifies the influential uncertain input
to guide decision-making. We propose a general sensitivity framework with
respect to the input distribution parameters that unifies a wide range of
sensitivity measures, including information theoretical metrics such as the
Fisher information. The framework is derived analytically via a constrained
maximization and the sensitivity analysis is reformulated into an eigenvalue
problem. There are only two main steps to implement the sensitivity framework
utilising the likelihood ratio/score function method, a Monte Carlo type
sampling followed by solving an eigenvalue equation. The resulting
eigenvectors then provide the directions for simultaneous variations of the
input parameters and guide the focus to perturb uncertainty the most. Not only
is it conceptually simple, but numerical examples demonstrate that the
proposed framework also provides new sensitivity insights, such as the
combined sensitivity of multiple correlated uncertainty metrics, robust
sensitivity analysis with an entropic constraint, and approximation of
deterministic sensitivities. Three different examples, ranging from a simple
cantilever beam to an offshore marine riser, are used to demonstrate the
potential applications of the proposed sensitivity framework to applied
mechanics problems.

_Keywords_ sensitivity matrix, parametric sensitivity, combined sensitivity,
information theoretical sensitivity, decision under uncertainty

##  1 Introduction

The use of mathematical models to simulate real-world phenomena is firmly
established in many areas of science and technology. The input data for the
models are often uncertain, as they could be from multiple sources and of
different levels of relevance. The uncertain inputs of a mathematical model
induce uncertainties in the output and sensitivity analysis identifies the
influential inputs to guide decision-making. A broad range of approaches can
be found in the literature, but in practice, the input uncertainties are
commonly quantified by a joint probability distribution. The analysis of the
input and output relationship in this probabilistic setting is called
probabilistic sensitivity analysis [1].

A suitable measure can be used to summarise the induced output uncertainties.
Commonly used metrics are the (central) moment functions of the uncertain
output, such as the mean and variance, and the probability of failure, i.e.,
the probability that the random output would exceed a certain threshold. In
addition, the average uncertainty or information content can be measured using
entropy that is based on the entire distribution function of the random output
[2]. The probabilistic sensitivity analysis then examines the relationship
between the uncertain input and the induced uncertainty of the output. In
particular, we are interested in identifying which input parameters would
impact the output metrics the most, i.e., the largest output change for the
same input variation, to guide decision-making.

In this setting, the sensitivity of the point estimates, such as the moment
functions and the failure probabilities, can be obtained using the partial
derivatives of the metrics with respect to (w.r.t) the input distribution
parameters. Although the general application of the derivative-based
sensitivity analysis can be limited by the difficulty of computing the
derivatives, the derivatives w.r.t the input distribution parameters can be
more easily evaluated by differentiation inside the expectation operator (c.f.
eqs. 1, 2, 3 and 4 in Section 2 ). This is possible because the individual
samples of the random output are not directly dependent on the input
distribution parameters. As a result, the partial derivative operation is only
evaluated w.r.t the joint probability density function (PDF) of the input, and
this approach is called the likelihood ratio/score function method (LR/SF) [3,
4].

As described, the LR/SF method is merely a mathematical trick. Nevertheless,
if used together with a sampling method, it is efficient as the uncertainty
metric and its sensitivity can be evaluated in a single simulation run (c.f.
Section 3.4). The LR/SF method has been applied to general objective functions
in stochastic optimization [3], the failure probability in reliability
engineering [5] and some distribution-free properties of the LR/SF method are
discussed in [6].

The sensitivity of entropy, on the other hand, cannot be directly evaluated
using the LR/SF method. Instead, sensitivity related to entropy is often
analysed using the Kullback–Leibler (K-L) divergence (aka relative entropy),
by measuring the divergence between two PDFs (probability density functions)
corresponding to two different cases. This approach is studied in [7] for
safety assessment to explore the impact on risk profile due to input
uncertainties and in [8] for engineering design before and after uncertainty
reduction of the random variables of interest. A similar approach using the
mutual information between the input and the output has also been studied for
sensitivity analysis [9]. The mutual information can be regarded as a special
form of the K-L divergence except that it requires the use of the joint PDF.
As the K-L divergence is not a metric, alternative distance measures such as
the Hellinger distance has been proposed to quantify the difference between
two PDFs and the corresponding sensitivities [10].

It should be noted although the relative entropy is not a metric, its
infinitesimal form is directly linked to the Fisher information [11] which is
a metric tensor and this link has been explored in [12] for probabilistic
sensitivity analysis using the Fisher information matrix (FIM). The LR/SF
method can then be used to compute the FIM efficiently for sensitivity
analysis of the output entropy [12].

In this paper, we propose a new sensitivity matrix $\bf r$ that unifies the
sensitivity of a wide range of commonly used uncertainty metrics, from moments
of the uncertain output to the entropy of the entire distributions, in a
single framework. This is made possible by the likelihood ratio/score function
method (LR/SF) where the sensitivity to the input distribution parameters of
different metrics can be expressed in the same form (c.f. Eq 4). The 2nd
moment of the sensitivity matrix,
$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]$, arises
naturally when the impact of input perturbation on the output is examined.
Moreover, the maximization of the perturbation of the output uncertainty
metrics leads to an eigenvalue problem of the matrix
$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]$. The
eigenvalues represent the magnitudes of the sensitivities with respect to
(w.r.t) simultaneous variations of the input distribution parameters
$\mathbf{b}$, and the relative magnitudes and directions of the variations are
given by the corresponding eigenvectors. Therefore, the eigenvectors
corresponding to the largest eigenvalues are the most sensitive directions to
guide decision-making.

The sensitivity matrix $\bf r$ can be seen as a counterpart of the
deterministic sensitivity matrix (Jacobian matrix) as the elements of $\bf r$
are the normalised partial derivatives of the output uncertainty metrics w.r.t
to the distribution parameters of the uncertain input (c.f. Eq 6). The
resulting eigenvectors, therefore, have direct sensitivity interpretation. It
should be noted that although the sensitivity matrix $\bf r$ is formulated and
estimated using the LR/SF method, the use of the 2nd moment matrix and its
eigenvectors for sensitivity analysis additionally captures the interactions
of the sensitivities of different metrics.

In addition, the current work is motivated by a recent study [13] where a
special case of the proposed sensitivity matrix has been applied successfully
to the combined sensitivity analysis of multiple failure modes. We are going
to show that, not only does
$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]$ capture the
combined perturbation effect of multiple metrics, e.g., multiple failure modes
or multiple moment functions, but also include the Fisher information matrix
(FIM) as a special case. Application of the FIM for sensitivity analysis can
be found in many areas of science and engineering. For example, the Fisher
Information Matrix (FIM) has been applied to the parametric sensitivity study
of stochastic biological systems [14], to assess the most sensitive directions
for climate change given a model for the present climate [15] and as one of
the process-tailored sensitivity metrics for engineering design [12].

It should be noted that there are two main differences between the proposed
framework and the commonly used variance-based sensitivity analysis [16].
First, variance-based approaches study how the variance of the output can be
decomposed into contributions from uncertain inputs. It ranks the factors
based on the assumption that the factor can be fixed to its true value, i.e.,
complete reduction of the uncertainties, which is rarely possible in practice
[1]. In contrast, the proposed framework uses partial derivatives to examine
the perturbation of the output metrics due to a variation of the input
distribution parameters. As the distribution parameters are often based on
data, it is equivalent to asking which uncertain dataset the decision-makers
should focus on to change the output the most. And this is particularly
pertinent to data-driven applications like digital twins [12]. Second, the
output sensitivity measure from the variance-based methods is the percentage
contribution, of each factor or the interactions between factors, to the
output variance. The proposed framework, on the other hand, outputs the
eigenvectors of the sensitivity moment matrix
$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]$ as the
principal sensitivity directions for simultaneous variations of the input
distribution parameters. This is based on a more pragmatic view that given a
finite budget to change the parameters, maximizing the impact on the output
follows the principal sensitivity directions, which tend to be a simultaneous
variation of the parameters because their effects on the output are likely to
be correlated. More discussions on the budget constraint can be found in
Section 3.2 with a generalization to the generalised eigenvalue problem.

It should be noted that despite the differences, for some cases, the
aggregated index for individual parameters given in Eq 25 can be used to
compare the Fisher sensitivity results against the variance-based main and
total sensitivity indices. This has been done in [17] for a 15-dimensional
problem, and in that case, the dominant first eigenvector of the FIM seems to
correspond to the main effects from the variance-based sensitivity analysis.

In what follows, the general sensitivity framework is introduced in Section 2
where the sensitivity analysis is reformulated as a standard eigenvalue
problem. In Section 3, we discuss various properties of the proposed
framework, including the link to the Fisher information matrix and the
possible extension to a generalised eigenvalue problem for robust sensitivity
analysis. A benchmark study, using two commonly used functions for sensitivity
analysis, is conducted in Section 4 for the Fisher sensitivity. Three
different examples are considered in Section 5, ranging from a simple
cantilever beam to an offshore marine riser, to demonstrate the potential
applications of the proposed sensitivity framework. Concluding remarks are
given in Section 6.

##  2 Sensitivity framework

Consider a general function $\bf y=h(x)$, the probabilistic sensitivity
analysis characterises the uncertainties of the outputs $\bf y$ that are
induced by the random inputs $\bf x$. It is assumed that the uncertainties of
$\bf x$ can be described by parametric probability distributions, i.e., $x\sim
p(\mathbf{x|b})$, where $\bf b$ are the distribution parameters.

One commonly used summary statistic is the (central) moment function of the
uncertain output, such as the mean and variance. More generally, the moment
function is taken with respect to a function of the uncertain output
$g(\mathbf{y})$. This might arise when there is a stochastic process present,
such as the random forces considered in some of the examples in Section 5, and
the $g(\cdot)$ function could represent max, min or root mean square (r.m.s).
In this setting, the $q^{\text{th}}$ moment function and its partial
derivative w.r.t the input distribution parameters can be expressed as:

$\displaystyle m_{q}$
$\displaystyle=\mathbb{E}_{X}\left[g^{q}(\mathbf{h(x)})\right]=\int
g^{q}(\mathbf{h(x)})p(\mathbf{x|b})\mathrm{d}\mathbf{x}$ (1a)
$\displaystyle\frac{\partial m_{q}}{\partial\mathbf{b}}$ $\displaystyle=\int
g^{q}(\mathbf{h(x)})\frac{\partial
p(\mathbf{x|b})}{\partial\mathbf{b}}\mathrm{d}\mathbf{x}$ (1b)

where it has been assumed that the differential and integral operators are
commutative, i.e. the order of the two operations can be exchanged under
regularity conditions of continuous and bounded functions.

Another metric is the probability of failure and its gradient:

$\displaystyle P_{f}$
$\displaystyle=\mathbb{E}_{X}\left[\mathrm{H}\left[g\left(\mathbf{h(x)}\right)-z\right]\right]=\int\mathrm{H}\left[g\left(\mathbf{h(x)}\right)-z\right]p(\mathbf{x|b})\mathrm{d}\mathbf{x}$
(2a) $\displaystyle\frac{\partial P_{f}}{\partial\mathbf{b}}$
$\displaystyle=\int\mathrm{H}\left[g\left(\mathbf{h(x)}\right)-z\right]\frac{\partial
p(\mathbf{x|b})}{\partial\mathbf{b}}\mathrm{d}\mathbf{x}$ (2b)

where $\mathrm{H}(\cdot)$ is the Heaviside step function and $z$ represents
the failure threshold. It is noted in passing that the application of failure
probability is not limited to reliability engineering. For example, the
probability of cost-effectiveness in health economics [18] and the probability
of acceptability in design [19] can both be formulated in the same way as Eq
2a.

When the quantity of interest is the underlying distribution function of the
uncertain outputs, the density function and its gradient w.r.t the
distribution parameters can be expressed as [12]:

$\displaystyle p(\mathbf{y})$
$\displaystyle=\mathbb{E}_{X}\left[\prod_{n}\delta\left[y_{n}-h_{n}(\mathbf{x})\right]\right]=\int\prod_{n}\delta\left[y_{n}-h_{n}(\mathbf{x})\right]p(\mathbf{x|b})\mathrm{d}\mathbf{x}$
(3a) $\displaystyle\frac{\partial p(\mathbf{y})}{\partial\mathbf{b}}$
$\displaystyle=\int\prod_{n}\delta\left[y_{n}-h_{n}(\mathbf{x})\right]\frac{\partial
p(\mathbf{x|b})}{\partial\mathbf{b}}\mathrm{d}\mathbf{x}$ (3b)

where $\delta(\cdot)$ is the Dirac delta function.

Although the aforementioned diverse metrics measure different aspects of the
uncertain output, it is clear that all of them can be more compactly described
using a general utility function:

$\displaystyle U$
$\displaystyle=\mathbb{E}_{X}\left[u\left(\mathbf{x}\right)\right]=\int
u\left(\mathbf{x}\right)p(\mathbf{x|b})\mathrm{d}\mathbf{x}$ (4a)
$\displaystyle\frac{\partial U}{\partial\mathbf{b}}$ $\displaystyle=\int
u\left(\mathbf{x}\right)\frac{\partial
p(\mathbf{x|b})}{\partial\mathbf{b}}\mathrm{d}\mathbf{x}$ (4b)

where the utility function $u(\mathbf{x})$ represents the $g^{q}(\cdot)$ in Eq
1, $\mathrm{H}(\cdot)$ in Eq 2 and $\delta(\cdot)$ in Eq 3. It should be noted
that the utility function could also depend on other variables, such as the
failure threshold $z$ for the case of failure probability. However,
$u(\mathbf{x})$ is not directly dependent on the parameters $\bf b$, as
$\mathbf{b}\rightarrow\mathbf{x}\rightarrow u(\mathbf{x})$ forms a Markov
chain. As a result, it is possible to differentiate the joint PDF
$p(\mathbf{x|b})$ within the integral in Eq 4b. And that is the same for eqs.
1, 2 and 3. As mentioned in the introduction, this approach is sometimes
called the likelihood ratio/score function method (LR/SF). An advantage of
this approach is that, if used together with a sampling method such as the
Monte Carlo method, the uncertainty quantification and sensitivity analysis
can be conducted in a single simulation run and more details are given in
Section 3.4.

The purpose of our sensitivity analysis is to identify the most important
uncertain parameters, i.e., which set of parameters would perturb the output
of interest the most. This perturbation can be quantified as $\Omega=(\Delta
U/U)^{2}$, where the normalisation leads to percentage perturbation and the
square operation quantifies the absolute value. If a first-order perturbation
is assumed, a general form of the normalised perturbation is:

$\begin{split}\Omega&=\mathbb{E}\left[\sum_{k}\left(\frac{\Delta
U_{k}}{U_{k}}\right)^{2}\right]\\\
&=\mathbb{E}\left[\sum_{k}\left(\frac{1}{U_{k}}\frac{\partial
U_{k}}{\partial\mathbf{b}}\Delta\mathbf{b}\right)^{2}\right]\\\
&=\sum_{i}\sum_{j}\Delta b_{i}\Delta
b_{j}\mathbb{E}\left[\sum_{k}r_{ik}r_{jk}\right]\\\
&=\Delta\mathbf{b}^{\operatorname{\mathsf{T}}}\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]\Delta\mathbf{b}\end{split}$
(5)

where the $jk^{\text{th}}$ entry of the matrix $\mathbf{r}$ is defined
accordingly as:

$r_{jk}=\frac{1}{U_{k}}\frac{\partial U_{k}}{\partial b_{j}}$ (6)

The matrix $\bf r$ can be seen as a counterpart of the deterministic
sensitivity matrix (Jacobian matrix) and therefore called sensitivity matrix
in this paper. It is interesting to note that the 2nd moment of the
sensitivity matrix,
$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]$, arises
naturally from the perturbation analysis. As it is in the form of a Gram
matrix, $\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]$ is
symmetric positive semi-definite (also evident from the quadratic form of Eq
5).

The general form of the perturbation in Eq 5 considers the combined effect of
multiple utilities. For example, there could be multiple failure modes where
$U_{k}=P_{f}^{(k)}$ denotes the $k^{\text{th}}$ failure mode; it is also often
of interest to consider the combined sensitivity of multiple responses or
moments of the same uncertain output where $U_{k}=m_{k}$ denotes the
$k^{\text{th}}$ moment. It is noted in passing that a weighting could be added
to each $U_{k}$ and that would result in a weighting of $r_{jk}$ in Eq 6. The
weighted scenario is not considered further in this paper as the weighting is
strongly case-dependent but will not alter the general form of Eq 5. The
expectation operation $\mathbb{E}[\cdot]$ in Eq 5 takes account of any
additional uncertainties that might arise in different cases. For example, the
failure threshold $z$ could be uncertain in Eq 2; for the case of the joint
density function, where $U=p(\mathbf{y})$, the gradient of the log utility
described in Eq 6 is uncertain due to randomness of the output $\bf y$.

Using the general perturbation function described in Eq 5, the sensitivity
analysis can be formulated as a constrained optimization problem:

$\begin{split}\text{max}&\quad\frac{1}{2}\Omega=\frac{1}{2}\Delta\mathbf{b}^{\operatorname{\mathsf{T}}}\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]\Delta\mathbf{b}\\\
\text{s.t.}&\quad\Delta\mathbf{b}^{\operatorname{\mathsf{T}}}\Delta\mathbf{b}=\epsilon\end{split}$
(7)

where the method of Lagrange Multiplier can be used:

$\displaystyle L$
$\displaystyle=\Delta\mathbf{b}^{\operatorname{\mathsf{T}}}\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]\Delta\mathbf{b}-\lambda(\Delta\mathbf{b}^{\operatorname{\mathsf{T}}}\Delta\mathbf{b}-\epsilon)$
(8a) $\displaystyle\frac{\partial L}{\partial\Delta\mathbf{b}}$
$\displaystyle=\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]\Delta\mathbf{b}-\lambda\Delta\mathbf{b}$
(8b)

where $\lambda$ is the Lagrange multiplier. Setting the first order optimality
condition for the Lagrangian, Eq 8b then leads to the following standard
eigenvalue problem:

$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]\mathbf{q}=\lambda\mathbf{q}$
(9)

The eigenvalues represent the magnitudes of the sensitivities with respect to
(w.r.t) simultaneous variations of the parameters $\mathbf{b}$, and the
relative magnitudes and directions of the variations are given by the
corresponding eigenvectors. As the solution to a maximization problem, the
eigenvectors corresponding to the largest eigenvalues then provide the most
perturbation of $\Omega$ in Eq 5.

##  3 Discussion

###  3.1 Information theoretical metrics as a special case

When the utility $U$ corresponds to a probably or probability density, the
expression for $\bf r$ in Eq 6 can be more compactly written as:

$r_{jk}=\frac{\partial\log U_{k}}{\partial b_{j}}$ (10)

As the log probability can be seen as the information content of a random
event in information theory, the expression in Eq 10 reveals the information
link of the sensitivity framework.

More concretely, for the case described in Eq 3, the utility corresponds to
$p(\mathbf{y})$ and the resulted perturbation is (using eqs. 5 and 6):

$\begin{split}\Omega&=\mathbb{E}_{Y}\left[\left(\frac{\Delta
p(\mathbf{y})}{p(\mathbf{y})}\right)^{2}\right]\\\
&=\Delta\mathbf{b}^{\operatorname{\mathsf{T}}}\mathbb{E}_{Y}\left[\frac{\partial\log
p(\mathbf{y})}{\partial\mathbf{b}}{\frac{\partial\log
p(\mathbf{y})}{\partial\mathbf{b}}}^{\operatorname{\mathsf{T}}}\right]\Delta\mathbf{b}\end{split}$
(11)

where the expectation part in the right hand side of the equation can be more
explicitly written as:

$\begin{split}\mathbb{E}_{Y}\left[\frac{\partial\log
p(\mathbf{y})}{\partial\mathbf{b}}{\frac{\partial\log
p(\mathbf{y})}{\partial\mathbf{b}}}^{\operatorname{\mathsf{T}}}\right]&=\int\frac{\partial\log
p(\mathbf{y})}{\partial\mathbf{b}}{\frac{\partial\log
p(\mathbf{y})}{\partial\mathbf{b}}}^{\operatorname{\mathsf{T}}}p(\mathbf{y})\mathrm{d}\mathbf{y}\\\
&=\int\frac{\partial p(\mathbf{y})}{\partial\mathbf{b}}{\frac{\partial
p(\mathbf{y})}{\partial\mathbf{b}}}^{\operatorname{\mathsf{T}}}\frac{1}{p(\mathbf{y})}\mathrm{d}\mathbf{y}\end{split}$
(12)

and this is the Fisher information matrix (FIM) [2] and is denoted as
$\mathbf{F}$. Therefore, the perturbation in Eq 11 can be rewritten as:

$\Omega=\Delta\mathbf{b}^{\operatorname{\mathsf{T}}}\mathbf{F}\Delta\mathbf{b}\approx
2\text{KL}\left[p(\mathbf{y|b})||p(\mathbf{y|b}+\Delta\mathbf{b})\right]$ (13)

where $\text{KL}[\cdot]$ indicates Kullback-Leibler (K-L) divergence and it is
also called relative entropy [2]. The approximation in Eq 13 can be found via
a Taylor expansion of the perturbed PDF $p(\mathbf{y|b}+\Delta\mathbf{b})$
from the relative entropy expression, with the third and higher order terms
ignored. It should be noted that Eq 13 has been derived in [12] to link the
relative entropy and the FIM. However, in this paper, we have extended that
link to a general perturbation metric $\Omega$ which unifies several other
metrics as well.

In Section 4, a benchmark study is conducted for the Fisher sensitivity. As
mentioned in the introduction, the application of the Fisher information for
sensitivity analysis can be found in many areas of science and engineering. In
the numerical examples given in Section 5, we will demonstrate that the FIM
can be utilised in different ways using the proposed framework.

###  3.2 Extension to generalised eigenvalue problems

The constraint in Eq 7 controls the potential change of the parameters. This
limit can be seen as a result of the finiteness of resources which is true for
all physical systems. A more general decision-oriented constraint can be
written as:

$\Delta\mathbf{b}^{\operatorname{\mathsf{T}}}\mathbf{W}\Delta\mathbf{b}=\epsilon$
(14)

where $\mathbf{W}$ is a weighting matrix which is symmetric. Instead of the
standard eigenvalue equation in Eq 9, the weighted constraint from Eq 14 leads
to a generalised eigenvalue problem:

$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]\mathbf{q}=\lambda\mathbf{W}\mathbf{q}$
(15)

Consider a scenario where the interest is to understand the sensitivity of the
failure probability. Eq 9 then leads us to the standard eigenvalue analysis of
the following matrix:

$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]=\mathbb{E}\left[\frac{\partial\log
P_{f}}{\partial\mathbf{b}}{\frac{\partial\log
P_{f}}{\partial\mathbf{b}}}^{\operatorname{\mathsf{T}}}\right]$ (16)

where the expectation is w.r.t the potential uncertain failure threshold. Eq
16 is rank-1 for a deterministic threshold. However, the resulting sensitivity
directions for $\Delta\mathbf{b}$ might impact both the safe and failure
regions. One way to mitigate the unwanted perturbation is to set an
uncertainty constraint using the relative entropy:

$\begin{split}\text{max}&\quad\frac{1}{2}\Omega=\frac{1}{2}\Delta\mathbf{b}^{\operatorname{\mathsf{T}}}\mathbb{E}\left[\frac{\partial\log
P_{f}}{\partial\mathbf{b}}{\frac{\partial\log
P_{f}}{\partial\mathbf{b}}}^{\operatorname{\mathsf{T}}}\right]\Delta\mathbf{b}\\\
\text{s.t.}&\quad\text{KL}\left[p(\mathbf{y|b})||p(\mathbf{y|b}+\Delta\mathbf{b})\right]=\epsilon\end{split}$
(17)

and this leads to the following generalised eigenvalue problem (with the
substitution of Eq 13 for the approximation of the relative entropy):

$\mathbb{E}\left[\frac{\partial\log
P_{f}}{\partial\mathbf{b}}{\frac{\partial\log
P_{f}}{\partial\mathbf{b}}}^{\operatorname{\mathsf{T}}}\right]\mathbf{q}=\lambda\mathbf{F}\mathbf{q}$
(18)

where the Fisher information matrix $\bf F$ now takes the position of $\bf W$
in Eq 15. Eq 18 will be utilised in one of the numerical examples to
demonstrate the application for constrained perturbation of a failure
probability.

###  3.3 Reparameterization and normalisation

The matrix $\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]$
depends on the parametrization used. Suppose
$b_{j}=\phi_{j}(\theta_{i}),\>i=1,2,\dots,s$, then the $k^{\text{th}}$ column
of the sensitivity matrix $\mathbf{r}$ with respect to $\bm{\uptheta}$ is:

$\mathbf{r}_{k}=\frac{1}{U_{k}}\frac{\partial
U_{k}}{\partial\bm{\uptheta}}=\mathbb{J}^{\operatorname{\mathsf{T}}}\frac{1}{U_{k}}\frac{\partial
U_{k}}{\partial\mathbf{b}}$ (19)

where $\mathbb{J}$ is the Jacobian matrix with $\mathbb{J}_{ji}={\partial
b_{j}}/{\partial\theta_{i}}$. As a result, the matrix
$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]$ w.r.t the
parameters $\bm{\uptheta}$ can be found via a reparameterization:

$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]_{\bm{\uptheta}}=\mathbb{J}^{\operatorname{\mathsf{T}}}\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]_{\mathbf{b}}\mathbb{J}$
(20)

Normalization w.r.t the parameters $\bf b$ is equivalent to a
reparametrization. In the case of proportional normalization, where the
sensitivity matrix becomes dimensionless, the change of parameter is
$b_{j}=\bar{b}_{j}\theta_{j}$ with $\bar{b}_{j}$ the nominal value for
normalization, and the Jacobian matrix in Eq 20 is just a diagonal matrix with
$\bar{b}_{j}$ on the diagonal. Due to its simplicity for sensitivity analysis,
proportional normalization is applied to the numerical cases considered in
Section 5.

###  3.4 Numerical considerations

A general function $\bf y=h(x)$ can rarely be solved analytically. The unique
mathematical form of eqs. 1, 2 and 3, and more generally Eq 4, allows
efficient computation of the gradient if a sampling method such as Monte Carlo
(MC) method is used. The Monte Carlo approximation of the integrals in Eq 4
results:

$\displaystyle\begin{split}U&=\int
u\left(\mathbf{x}\right)p(\mathbf{x|b})\mathrm{d}\mathbf{x}\\\
&\approx\frac{1}{N}\sum_{i}u\left(\mathbf{x}_{i}\right)\end{split}$ (21a)
$\displaystyle\begin{split}\frac{\partial U}{\partial\mathbf{b}}&=\int
u\left(\mathbf{x}\right)\frac{\partial\log
p(\mathbf{x|b})}{\partial\mathbf{b}}p(\mathbf{x|b})\mathrm{d}\mathbf{x}\\\
&\approx\frac{1}{N}\sum_{i}u\left(\mathbf{x}_{i}\right)\frac{\partial\log
p(\mathbf{x}_{i}|\mathbf{b})}{\partial\mathbf{b}}\end{split}$ (21b)

where $\mathbf{x}_{i}$ is a MC realisation of the random variable $\bf x$ and
$N$ MC simulations are considered.
For many commonly used distributions, analytical closed-form expressions can
be obtained for the partial derivatives with respect to (w.r.t) a distribution
parameter. For example, for a Gaussian distribution:

$\displaystyle\frac{\partial\log p(x|\mu,\sigma)}{\partial\mu}$
$\displaystyle=\frac{x-\mu}{\sigma}$ (22a) $\displaystyle\frac{\partial\log
p(x|\mu,\sigma)}{\partial\sigma}$
$\displaystyle=\frac{(x-\mu)^{2}-\sigma^{2}}{\sigma^{3}}$ (22b)

where $\mu$ and $\sigma$ are the mean and standard deviation of the Gaussian
PDF. Therefore, with the analytical expressions available, the utility of
interest and its gradient in Eq 21 can be obtained in a single computational
run. This is one of the main numerical advantages of the adopted likelihood
ratio/score function method (LR/SF) as described in the introduction. The
application of Eq 21 has been validated in [13] for a failure probability
sensitivity, where the perturbation approximated using Eq 21b agrees well with
the exact results from a direct perturbation. In [12], the Fisher information
estimated from Eq 21 also identifies the influential parameters for a marine
riser as expected. Not only is it efficient, but it is also clear that Eq 21
is independent of the parameter dimensions and thus the sensitivity approach
is not limited by the dimension of the input variables. It should be noted
that although the sensitivity framework introduced applies to dependent
inputs, it is assumed for simplicity in the numerical implementation below
that the components of $\bf x$ are independent.

As the numerical computation is based on the LR/SF method, the sensitivity
estimation is independent of the dimension of the input parameters [4]. This
is in contrast to the variance-based methods where the computational cost is
proportional to the input dimension [16]. However, a large number of sampling
points might be required to maintain a low variance from the LR/SF estimation.
As our main purpose is to present the new sensitivity framework, only the
standard Monte Carlo (MC) method is used. Advanced simulation methods, such as
importance sampling and subset simulation, can be used to improve the sampling
efficiency.

Based on the LR/SF method estimation, the proposed framework requires an
eigen-decomposition of the
$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]$ matrix. As
the size of the 2nd moment matrix depends on the dimension of the input
parameters, for high dimension problems, an iterative approach might be needed
to solve the eigenvalue decomposition. Nevertheless, as the 2nd moment matrix
is symmetric and it is often the first few largest eigenvalues that are of
interest for sensitivity analysis, additional cost for the eigenvalue problem
should be minimal compared to the MC sampling. As a result, the computational
cost of our method is similar to the LR/SF method in general.

##  4 Benchmark study

The Fisher information has been introduced in [12] for sensitivity analysis
with respect to the input distribution parameters. The examples given in [12]
have focused on the demonstration for engineering design applications, but
lack benchmarking against commonly used functions for sensitivity analysis. As
the Fisher sensitivity will be used throughout the examples in Section 5, in
this section, we present benchmark studies against two of the commonly used
variable screening benchmark functions. For screening purposes, the input
variables for both functions are assumed to be from a standard Gaussian
distribution, i.e., $x\sim\mathbb{N}(0,1)$.

###  4.1 Decreasing coefficient function

The first benchmark function is a linear function with decreasing
coefficients:

$f(\mathbf{x})=0.2x_{1}+\frac{0.2}{2}x_{2}+\frac{0.2}{4}x_{3}+\frac{0.2}{8}x_{4}+\frac{0.2}{16}x_{5}+\frac{0.2}{32}x_{6}+\frac{0.2}{64}x_{7}+\frac{0.2}{128}x_{8}+e$
(23)

where $e\sim\mathbb{N}(0,0.05^{2})$ is a small error term. This function has
been used for input variable screening in [20], where it concluded that the
first three variables, $x_{1},\ x_{2},\ \text{and}\ x_{3}$, are more important
than the rest.

The sensitivity analysis using the Fisher information matrix (FIM) from Eq 12
is applied to the function above and the results are presented in Figure 1. In
contrast to variance-based sensitivity approaches, as explained in the
introduction, the proposed framework provides the principal sensitivity
directions w.r.t the distribution parameters, e.g., the mean and standard
deviation (Std Dev) in Figure 1.

The sensitivity results in Figure 1 display the first two eigenvectors of the
FIM, as the rest of the eigenvalues have negligible amplitudes as shown in
Figure 2. The results in Figure 1 clearly show that the importance of the
variables decreases as the subscript index increases. It is clear that the
first three variables are more important than the rest as seen from the 1st
eigenvector, and that agrees with the conclusion from [20]. Moreover, for the
1st eigenvector shown in Figure 1(a), the ratio between the sensitivities to
the standard deviations of the first three variables are
$[x_{1}:x_{2}:x_{3}]=[1:0.236:0.065]$. This ratio is very close to the ratio
of the square of the coefficients, i.e., $[1:1/4:1/16]$. Similarly, the
dominant sensitivities to the mean of the first three variables, as given in
Figure 1(b), have the ratio of $[1:0.496:0.260]$. This suggests that not only
does the FIM sensitivity index give the correct ranking of the variables, but
it has also accurately captured the relative strength of the linear
coefficients of Eq 23.

In Figure 2, the spectrum of the FIM eigenvalues is given, for different
numbers of Monte Carlo samples (c.f. Section 3.4). This serves as a
convergence test for the FIM estimation. The same type of convergence tests
has also been performed for the numerical examples in Section 5 but not
presented for conciseness.

Figure 1: Fisher sensitivity results for the decreasing coefficient function
described in Eq 23.  Figure 2: The eigenvalue spectrum of the Fisher
information matrix, for the decreasing coefficient function described in Eq
23.

###  4.2 Roos & Arnold Function

The second benchmark function is a product function:

$f(\mathbf{x})=\prod_{i=1}^{d}|4x_{i}-2|$ (24)

This function occurs as an integrand multiple times in the literature. It has
been classified as a Type C function in [21], meaning that it has dominant
high-order interaction terms. It is clear from Eq 24 that there is no
difference between the input variables. The results from the FIM, similar to
the previous example, are shown in Figure 3 with five input variables. In this
case, the first eigenvalue is much bigger than the rest as shown in Figure 4
(which also shows the convergence results as previous example). The
sensitivity results in Figure 3 agree well with expectations that the input
variables are of similar importance. It should be noted that the amplitudes of
the five variables shown in Figure 3 are not exactly the same from the Fisher
sensitivity. This is because the FIM is based on a sampling approach and the
samples for the five input variables will not be exactly the same as they are
independently sampled.

Figure 3: Fisher sensitivity results for the Roos & Arnold function described
in Eq 24.  Figure 4: The eigenvalue spectrum of the Fisher information matrix,
for the Roos & Arnold function described in Eq 24.

##  5 Application examples

(a) Case 1

(b) Case 2

(c) Case 3

Figure 5: Three demonstrating examples. (a) Case 1, a cantilever beam with
white noise excitation; (b) Case 2, an offshore marine riser subject to random
wave loading; (c) Case 3, a model floating column in a wave tank.

###  5.1 Combined sensitivity for multiple responses

In this section, we use a simple cantilever beam as an example to demonstrate
the combined sensitivity analysis for multiple responses using the proposed
sensitivity framework. The cantilever beam, case 1 in Figure 5, is subject to
a white noise excitation of unit amplitude. The frequency response functions
for both acceleration and strain responses, at different positions along the
beam, are obtained via modal summation. To keep it simple, only the first mode
is considered with a constant modal damping of 0.1. The linear vibration
equation of a cantilever beam can be found in many vibration/mechanics
textbooks, see e.g., [22], and it is therefore not repeated here for
conciseness. The code for this example study can be found in the address given
in Data availability statement.

The quantities of interest for our sensitivity analysis are the peak r.m.s
acceleration and strain responses. ‘peak’ indicates the maximum response along
the beam for each sample of the random inputs. The $g(\cdot)$ function in Eq 1
is thus the composition of the r.m.s and maximum functions in this case. The
two types of response are normalised by the maximum values across the ensemble
of the random samples, i.e., the peak r.m.s results are between 0 and 1.

Table 1: Mean ($\mu$) and Coefficient of Variation (CoV) for the random
variables.

Young’s Modulus Density Length Width Thickness $E[Pa]$ $\rho[kg/m^{3}]$ $L[m]$
$w[m]$ $t[m]$ Mean 69e9 2700 0.45 2e-2 2e-3 Scenario-A Normal distribution
with CoV = 0.1 Scenario-B Gamma distribution with CoV = 0.5

Figure 6: Variable ranking for multiple responses of the cantilever beam in
Figure 5a. $\mathbf{r_{1}}$, $\mathbf{r_{2}}$ indicate the sensitivity of the
$U_{1}$ (acceleration) and $U_{2}$ (strain) respectively. $\Omega$ is for the
combined sensitivity results for $U_{1}$ and $U_{2}$ using the metric proposed
in this paper, and ‘Fisher’ denotes the sensitivity results from the FIM using
the jPDF of the responses. Error bars indicate the standard deviation from 10
repetitions of estimation with 20000 samples for each run. (a) Scenario-A with
Normal random input; (b) Scenario-B with Gamma random input.

The parameters for the random variable are listed in Table 1. Two scenarios
are considered with the same mean value but different types of distributions,
one with Normal (Gaussian) distribution and the other one with Gamma
distribution, of which the coefficients of variation (CoV) are different. The
sensitivity results are shown in Figure 6 for the relative ranking of the five
input random variables. One measure for the relative importance of the
$j^{\text{th}}$ variable can be obtained as:

$s_{j}^{2}=\sum_{i}\lambda_{i}q_{ji}^{2}$ (25)

where $\lambda$ and $\bf q$ are the eigenvalue and eigenvector of the
sensitivity matrix from Eq 9.

Eq 25 takes a Pythagorean view to estimate the contributions, in analogy to
principal component analysis, within each eigenvector and across different
eigenvectors using the eigenvalue amplitudes. However, this summary index from
Eq 25 assumes that different principal sensitivity directions, represented by
the eigenvectors, can be varied at the same time for the parameters. This
aggregated view of importance measure essentially only considers the diagonal
entries of the moment matrix
$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]$, as shown in
Appendix A. It therefore ignores the interactions between the parameters (off-
diagonal terms). Furthermore, the phase information of the sensitivity
vectors, i.e., increase or decrease, are also lost using this summary index
from Eq 25. This summary sensitivity index is used only in this case to give a
better indication that the proposed sensitivity metric does account for the
combined effect, as a direct comparison between eigenvectors is difficult
between two different metrics. As the purpose of the sensitivity framework is
to find the influential set of distribution parameters, as discussed in the
introduction, sensitivity results from the eigenvectors will be used for the
next two examples.

In this case study, the utility function to form $\bf r$ in Eq 6 is the mean
value of the uncertain response (the peak r.m.s response for each sample of
the random input), i.e., $U_{1}$ and $U_{2}$ are the averaged peak r.m.s
acceleration and strain responses respectively. $\mathbf{r}_{1}$ and
$\mathbf{r}_{2}$ are the normalised LR/SF sensitivity vectors of the two
corresponding utility functions as in Eq 6.

In Figure 6, the combined sensitivity result is indicated by $\Omega$. This is
computed using Eq 25 with the eigenvalues/eigenvectors of the
$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]$ matrix from
Eq 9, where $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$ are the intermediate results
for forming the
$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]$ matrix in
this case.

In comparison, the sensitivity results using the Fisher information matrix
(FIM) form Eq 12 is also presented in Figure 6. It can be seen in Figure 6
that $\mathbf{r}_{1}$ ranks $L$ as more important than $t$, while $\rho$ and
$L$ are of little importance from $\mathbf{r}_{2}$. In contrast, although the
absolute agreement varies, the resulting relative ranking of the random
variables are very similar between the combined analysis ($\Omega$) and the
FIM. The conclusion is similar between scenario-A and scenario-B.

The FIM results are based on the joint probability density function of the
acceleration and strain responses. This is in contrast to the combined
sensitivity analysis where the two types of responses are assumed to be
independent. Nevertheless, the comparison in Figure 6 demonstrates that the
combined analysis takes account of the combined perturbation of multiple
responses, as indicated by the derivation in Eq 5. This is in line with the
findings in [13] for multiple correlated failure modes, which can be seen as a
special case of the proposed sensitivity framework in this paper.

###  5.2 Robust failure sensitivity

Table 2: Mean ($\mu$) and Coefficient of Variation (CoV) for the random
variables of Case 2

Morison’s added Morison’s Riser Young’s Oil Top Material S-N curve
coefficients $C_{a}$ [-] $C_{d}$ [-] $\rho$ [kg/m3] $E$ [GPa] $\rho_{0}$
[kg/m3] $T_{0}$ [kN] $\alpha$ [GPa] $\delta$ [-] Mean 1.5 1.1 7840 200 920
4905 199 3 CoV 0.20 0.20 0.05 0.05 0.10 0.10 0.10 0.10 In this section, the
generalised eigenvalue problem introduced in Section 3.2 is considered. In
particular, we use Eq 18 to demonstrate the possibility to perturb a failure
probability with an entropic constraint. This might arise in practice to avoid
unwanted perturbation of the system responses, e.g., no perturbation of the
response in the safe region, and it can be regarded as a robust failure
sensitivity analysis. The example system considered is an offshore marine rise
shown in Figure 5b (Case 2) that is subject to a random wave excitation and
the nonlinear wave structure interaction is included in the model. The
simulation model for this case study, and Case 3 below, has been developed
using the CHAOS hydrodynamic code [23] which uses the semi-empirical Morison’s
equation [24] to estimate wave forces. This case study is taken from [12]
where the details of this model and its sensitivities can be found. Different
from [12] where the failure sensitivity is compared to the Fisher sensitivity
results, an entropic constraint is applied to the failure sensitivity in this
paper. The code for this example study can be found in the address given in
Data availability statement.

Figure 7: Constrained sensitivity of the failure probability due to fatigue
for Case 2. In comparison, the standard sensitivity results of the fatigue
failure is also shown.

The parameters for the uncertain variables of Case 2 are listed in Table 5.2
and the corresponding failure sensitivity results are given in Figure 7. The
results marked as ‘Constrained’ are obtained from solving the generalised
eigenvalue equation in Eq 18. The corresponding results of the standard
failure sensitivity vector, which is equivalent to substituting an identity
matrix for the matrix $\bf F$ in Eq 18, is shown in comparison and is denoted
as ‘Standard’. Note that the standard results are the same as in [12].

As can be seen in Figure 7, the entropy-constrained sensitivity is completely
dominated by the S-N coefficients $\alpha$ and $\delta$. This is because the
FIM is not dependent on the S-N coefficients that are specific to the fatigue
failure analysis. Although a relatively extreme example, this case study
demonstrates the potential application of the entropy-constrained sensitivity
analysis. It also shows the general applicability of the proposed sensitivity
framework once formulated as a generalised eigenvalue problem.

###  5.3 Limiting approximation for deterministic inputs

Figure 8: Sensitivity of the two natural frequencies, with a comparison
between the dominant Fisher eigenvector and the sensitivity obtained from the
analytical approach. The nominal values of the input parameters are used for
the deterministic case, and as the mean values for the random case. (a) 1st
natural frequency $\omega_{1}$; (b) 2nd natural frequency $\omega_{2}$ Table
3: Mean ($\mu$) and Coefficient of Variation (CoV) for the random variables of
Case 3

Material Water Column Tether Ballast Column Column Ballast Mass coefficient
$\rho$ [kg/m3] $\rho_{f}$ [kg/m3] $L$ [m] $L_{S}$ [m] $L_{b}$ [m] $r$ [m] $t$
[m] $m_{b}$ [kg] $C_{a}$ [-] Mean 1180 1025 1 0.2 0.15 4.5e-2 3.5e-3 3 1 CoV
1e-4 ($\delta$-approximation with a small CoV) In this section, we demonstrate
the potential application of the proposed sensitivity framework for
deterministic inputs, where a $\delta$-type distribution is assumed for
approximation. In this example, the natural frequency sensitivity of a model
floating column in a wave tank environment is analysed. This two-degrees-of-
freedom floating system, as seen in Figure 5c, is tethered to the bottom of
the wave tank and a ballast mass is added to adjust the centre of gravity.
This example is chosen here mainly to represent the commonly encountered
natural frequency design problem. Furthermore, a closed-form sensitivity
analysis for the natural frequencies, as given in Appendix B, can be obtained
straightforwardly as a benchmark. The example has been modelled using this
hydrodynamic code [23] where the added inertia effect of the water is
considered. The partial derivatives of the mass and stiffness matrices, which
are required for the analytical sensitivity results as seen in Appendix B, are
obtained using the symbolic differentiation module within Matlab. The code for
this example study can be found in the address given in Data availability
statement.

The sensitivity results for the two natural frequencies, namely $\omega_{1}$
and $\omega_{2}$, are displayed in Figure 8 for the parameters listed in Table
5.3. The ‘Deterministic sensitivity’ shows the results from the analytical
analysis, as given in Appendix B, using the mean values of the input variables
as the evaluation values. To use the proposed sensitivity framework, the input
variables are assumed to be uncertain with a Gaussian distribution and the
mean and CoV listed in Table 5.3. The Fisher information matrix (FIM), based
on Eq 12, is formed for the two natural frequencies and the eigenvectors of
the FIM with large eigenvalues then provide us with the sensitivity
information. In order to compare with the analytical deterministic
sensitivities, the FIM is calculated for each natural frequency separately,
i.e., one-dimensional PDFs of each natural frequency rather than the joint
PDF. In this case, only the 1st eigenvectors of FIM are shown as only the 1st
eigenvalue dominates. Very good agreement can be observed between the FIM
results and the analytical sensitivities. In particular, the relative phase of
the parameter $L_{b}$ is captured by the FIM despite its relatively low
sensitivity.

##  6 Conclusions

A sensitivity matrix $\bf r$ is proposed as a new probabilistic sensitivity
metric with respect to the input distribution parameters. The sensitivity of a
wide range of commonly used uncertainty metrics, from moments of the uncertain
output to the entropy of the entire distribution, can be formulated as an
eigenvalue problem of the 2nd moment of the proposed sensitivity matrix. The
resulting framework has a solid mathematical underpinning, is numerically
efficient, and unifies the sensitivity analysis in a general but conceptually
simple framework. And that is the main contribution of the present work.

The proposed framework is derived analytically via a constrained maximization
of the perturbation of the output uncertainties. On top of this mathematical
foundation is the framework’s conceptual simplicity, where its implementation
only consists of two main steps, a Monte Carlo type sampling followed by
solving an eigenvalue equation. Through the numerical examples, it is
demonstrated that the sensitivity framework can be applied for the combined
sensitivity analysis of multiple responses, even if the degree of correlation
between the responses is unknown. This is in line with the findings for
multiple correlated failure modes studied in [13]. In addition, using the
failure sensitivity as an example, it is shown that a robust sensitivity
analysis can be formed using entropy as a constraint and solved as a
generalised eigenvalue problem. Furthermore, the Fisher information matrix, a
special case of the proposed sensitivity metric, is shown to approximate
deterministic sensitivities very closely using $\delta$-type distribution
inputs.

A key element of the framework is the sensitivity matrix $\bf r$. As its
elements are the normalised partial derivatives of the expected utility of
interest, the resulting eigenvectors based on the sensitivity matrix have
direct sensitivity interpretations. Utilising the likelihood ratio/score
function method, when a sampling approach is used, the expected utility and
its derivatives can be obtained in a single simulation run. This allows
efficient computation of the sensitivity matrix $\bf r$ and the corresponding
$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]$ matrix.

The stochastic aspects of the output of interest have been considered
implicitly in this study. For example, the stochastic responses along the beam
in case 1 have been considered by setting the $g(\cdot)$ function as the r.m.s
operation. Future work will consider a stochastic variant of the proposed
framework, including sensitivities with respect to time-dependent inputs and
treating the expected utility as stochastic by updating the expectation
operation in Eq 5.

## Acknowledgment

This work has been funded by the Engineering and Physical Sciences Research
Council through the award of a Programme Grant “Digital Twins for Improved
Dynamic Design”, Grant No. EP/R006768. For the purpose of open access, the
author has applied a Creative Commons Attribution (CC BY) licence to any
Author Accepted Manuscript version arising. The author is grateful to
Professor Robin Langley, University of Cambridge, for the support to publish
this work.

## Data availability statement

The datasets generated during and/or analysed during the current study are
available in the GitHub repository:

https://github.com/longitude-jyang/Probabilistic-sensitivity-framework

## References

  * [1] Jeremy E. Oakley and Anthony O’Hagan.  Probabilistic sensitivity analysis of complex models: a Bayesian approach.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 66(3):751–769, 2004.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1467-9868.2004.05304.x.
  * [2] Thomas M. Cover and Joy A. Thomas.  Elements of Information Theory.  John Wiley & Sons, Incorporated, Hoboken, UNITED STATES, 2006.
  * [3] James C. Spall.  Introduction to stochastic search and optimization: estimation, simulation, and control.  Wiley-Interscience series in discrete mathematics and optimization. Wiley-Interscience, Hoboken, N.J, 2003.
  * [4] Reuven Y. Rubinstein and Dirk P. Kroese.  Simulation and the Monte Carlo Method.  Wiley Publishing, 3rd edition, 2016.
  * [5] Jinghui Li, Ali Mosleh, and Rui Kang.  Likelihood ratio gradient estimation for dynamic reliability applications.  Reliability Engineering & System Safety, 96(12):1667–1679, December 2011.
  * [6] Harry Millwater.  Universal properties of kernel functions for probabilistic sensitivity analysis.  Probabilistic Engineering Mechanics, 24(1):89–99, January 2009\.
  * [7] Chang K. Park and Kwang-Il Ahn.  A new approach for measuring uncertainty importance and distributional sensitivity in probabilistic safety assessment.  Reliability Engineering System Safety, 46(3):253–261, 1994.
  * [8] Huibin Liu, Wei Chen, and Agus Sudjianto.  Relative Entropy Based Method for Probabilistic Sensitivity Analysis in Engineering Design.  Journal of Mechanical Design, 128(2):326–336, April 2005.
  * [9] Niklas Lüdtke, Stefano Panzeri, Martin Brown, David S Broomhead, Joshua Knowles, Marcelo A Montemurro, and Douglas B Kell.  Information-theoretic sensitivity analysis: a general method for credit assignment in complex networks.  Journal of The Royal Society Interface, 5(19):223–235, February 2008.  Publisher: Royal Society.
  * [10] Gaofeng Jia and Alexandros A. Taflanidis.  Sample-based evaluation of global probabilistic sensitivity measures.  Computers & Structures, 144:103–118, 2014.
  * [11] R. A. Fisher and Edward John Russell.  On the mathematical foundations of theoretical statistics.  Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 222(594-604):309–368, 1 1922.  publisher: Royal Society.
  * [12] Jiannan Yang, Robin S. Langley, and Luis Andrade.  Digital twins for design in the presence of uncertainties.  Mechanical Systems and Signal Processing, 179:109338, November 2022\.
  * [13] Jiannan Yang, Arnau Clot, and Robin S. Langley.  Combined sensitivity analysis for multiple failure modes.  Computer Methods in Applied Mechanics and Engineering, 395:115030, May 2022.
  * [14] Rudiyanto Gunawan, Yang Cao, Linda Petzold, and Francis J. Doyle.  Sensitivity Analysis of Discrete Stochastic Systems.  Biophysical Journal, 88(4):2530–2540, April 2005.  Publisher: Elsevier.
  * [15] Andrew J. Majda and Boris Gershgorin.  Quantifying uncertainty in climate change science through empirical information theory.  Proceedings of the National Academy of Sciences, 107(34):14958–14963, August 2010.  Publisher: Proceedings of the National Academy of Sciences.
  * [16] A. Saltelli, editor.  Global sensitivity analysis: the primer.  John Wiley, Chichester, England ; Hoboken, NJ, 2008.  OCLC: ocn180852094.
  * [17] Jiannan Yang.  Decision-oriented two-parameter fisher information sensitivity using symplectic decomposition, 2022.
  * [18] Gianluca Baio and A Philip Dawid.  Probabilistic sensitivity analysis in health economics.  Statistical Methods in Medical Research, 24(6):615–634, December 2015.  Publisher: SAGE Publications Ltd STM.
  * [19] David R Wallace, Mark J Jakiela, and Woodie C Flowers.  Design search under probabilistic specifications using genetic algorithms.  Computer-Aided Design, 28(5):405–421, May 1996.
  * [20] Crystal Linkletter, Derek Bingham, Nicholas Hengartner, David Higdon, and Kenny Q Ye.  Variable selection for gaussian process models in computer experiments.  Technometrics, 48(4):478–490, 2006.
  * [21] Sergei Kucherenko, Balazs Feil, Nilay Shah, and Wolfgang Mauntz.  The identification of model effective dimensions using global sensitivity analysis.  Reliability Engineering & System Safety, 96(4):440–449, 2011.
  * [22] L. Meirovitch.  Elements of Vibration Analysis.  McGraw-Hill, 1986.
  * [23] Jiannan Yang.  Code for Hydrodynamic Analysis of Offshore Structures (CHAOS), January 2022.  10.5281/zenodo.5831833.
  * [24] Turgut ’Sarp’ Sarpkaya.  Wave Forces on Offshore Structures.  Cambridge University Press, 2010.

##  Appendix A The summary sensitivity index

The eigenvalue equation from Eq 9 can be rewritten as a matrix decomposition
(as $\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]$ is
symmetric):

$\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]\mathbf{=}\bf
Q\Lambda Q^{\operatorname{\mathsf{T}}}$ (A.1)

where $\bf Q$ is the eigenvector matrix with $\mathbf{q}_{i}$ as its $i_{th}$
column and $\Lambda$ is the diagonal eigenvalue matrix. Eq A.1 can be further
expressed as a summation over the eigenvalues:

$\begin{split}\mathbb{E}\left[\mathbf{rr}^{\operatorname{\mathsf{T}}}\right]&=\sum_{i}\lambda_{i}\mathbf{q}_{i}\mathbf{q}_{i}^{\operatorname{\mathsf{T}}}\\\
&=\sum_{i}\lambda_{i}\begin{bmatrix}q_{1i}^{2}&\dots&q_{1i}q_{ji}&\dots\\\
&\ddots&&\\\ \text{sym}&&q_{ji}^{2}&\\\ &&&\ddots\end{bmatrix}\end{split}$
(A.2)

Therefore, the sensitivity summary index in Eq 25, where
$s_{j}^{2}=\sum_{i}\lambda_{i}q_{ji}^{2}$, is essentially the diagonal entries
of the moment matrix.

##  Appendix B Analytical derivation for natural frequency sensitivity

The natural frequency $\omega$ of a linear discrete vibration system can be
found from the following eigenvalue equation:

$\mathbf{K\bm{\upphi}}=\omega^{2}\mathbf{M\bm{\upphi}}$ (B.1)

where $\bf K$ and $\bf M$ are the stiffness and mass matrices, $\bm{\upphi}$
is the eigenvector and $\omega^{2}$ is the corresponding eigenvalue. Without
loss of generality, the eigenvectors are assumed to be mass normalised, i.e.,
$\bm{\upphi}^{\operatorname{\mathsf{T}}}\mathbf{M}\bm{\upphi}=1$. Therefore,
the partial derivatives of the squared natural frequency with respect to the
parameters are:

$\begin{split}\frac{\partial\omega^{2}}{\partial
b_{j}}&=\frac{\partial}{\partial
b_{j}}\left(\bm{\upphi}^{\operatorname{\mathsf{T}}}\mathbf{K}\bm{\upphi}\right)\\\
&=\omega^{2}\left(\bm{\upphi}_{b_{j}}^{\operatorname{\mathsf{T}}}\mathbf{M}\bm{\upphi}+\bm{\upphi}^{\operatorname{\mathsf{T}}}\mathbf{M}\bm{\upphi}_{b_{j}}\right)+\bm{\upphi}^{\operatorname{\mathsf{T}}}\mathbf{K}_{b_{j}}\bm{\upphi}\\\
&=\omega^{2}\left(-\bm{\upphi}^{\operatorname{\mathsf{T}}}\mathbf{M}_{b_{j}}\bm{\upphi}\right)+\bm{\upphi}^{\operatorname{\mathsf{T}}}\mathbf{K}_{b_{j}}\bm{\upphi}\\\
&=\bm{\upphi}^{\operatorname{\mathsf{T}}}\left(\mathbf{K}_{b_{j}}-\omega^{2}\mathbf{M}_{b_{j}}\right)\bm{\upphi}\end{split}$
(B.2)

where the subscript $b_{j}$ denotes the partial derivative with respect to the
parameter $b_{j}$. In our case study with random input variables, the
parameter $b_{j}$ is taken as the mean value of the $j_{th}$ variable. The 2nd
to 3rd step of Eq B.2 makes use of the fact that $\frac{\partial}{\partial
b_{j}}\left(\bm{\upphi}^{\operatorname{\mathsf{T}}}\mathbf{M}\bm{\upphi}\right)=0$.
The partial derivatives of the natural frequency are then:

$\frac{\partial\omega}{\partial
b_{j}}=\frac{1}{2\omega}\bm{\upphi}^{\operatorname{\mathsf{T}}}\left[\frac{\partial\mathbf{K}}{\partial
b_{j}}-\omega^{2}\frac{\partial\mathbf{M}}{\partial b_{j}}\right]\bm{\upphi}$
(B.3)

and the normalised sensitivity is:

$r_{\omega}^{j}=\frac{\partial\omega}{\partial b_{j}}\frac{b_{j}}{\omega}$
(B.4)