# Some Identities in Quantum Torus Arising from Ringel-Hall Algebras

Jiuzhao Hua

###### Abstract

We define two classes of representations of quivers over arbitrary fields,
called monomorphic representations and epimorphic representations. We show
that every representation has a unique maximal nilpotent subrepresentation and
the associated quotient is always monomorphic, and every representation has a
unique maximal epimorphic subrepresentation and the associated quotient is
always nilpotent. The uniquenesses of such subrepresenations imply two
identities in the Ringel-Hall algebra. By applying Reineke’s integration map,
we obtain two identities in the corresponding quantum torus.

##  1 Introduction

Let $\mathbb{N}$ be the set of all non-negative integers, $\mathbb{Z}$ the
ring of integers and $\mathbb{Q}$ the field of rational numbers. Let
$Q=(Q_{0},Q_{1})$ be a finite quiver and assume that $Q$ is connected and
$Q_{0}=\\{1,2,\cdots,n\\}$ with $n\geq 1$. For any arrow $h\in Q_{1}$, let
$h^{\prime}$ and $h^{\prime\prime}$ be its source vertex and target vertex
respectively, depicted as
${}_{h^{\prime}}\bullet\overset{h}{\longrightarrow}\bullet_{h^{\prime\prime}}$.
The Euler form associated with $Q$ is defined as follows:

$\langle v,w\rangle:=\sum_{i\in Q_{0}}v_{i}w_{i}-\sum_{h\in
Q_{1}}v_{h^{\prime}}w_{h^{\prime\prime}}\,\,\text{ for }v,w\in\mathbb{Z}^{n}.$
(1)

A representation of $Q$ over a field $k$ is a collection of finite dimensional
$k$-vector spaces and $k$-linear maps $M=(V_{i},f_{h})_{i\in Q_{0},h\in
Q_{1}}$ such that $f_{h}$ is a $k$-linear map from $V_{h^{\prime}}$ to
$V_{h^{\prime\prime}}$.
$(\dim_{k}V_{1},\cdots,\dim_{k}V_{n})\in\mathbb{N}^{n}$ is called the
dimension vector of $M$, denoted by $\dim M$. A morphism from $M$ to another
representation $N=(W_{i},g_{h})_{i\in Q_{0},h\in Q_{1}}$ is a collection of
$k$-linear maps $(\phi_{i}:V_{i}\to W_{i})_{i\in Q_{0}}$ such that
$\phi_{h^{\prime\prime}}f_{h}=g_{h}\phi_{h^{\prime}}$ for all $h\in Q_{1}$. If
$\phi_{i}$ is non-singular for all $i\in Q_{0}$, then $M$ and $N$ are called
isomorphic. $N$ is called a subrepresentation of $M$ if $W_{i}\subset V_{i}$
for all $i\in Q_{0}$ and $g_{h}(x)=f_{h}(x)$ for all $h\in Q_{1}$ and $x\in
W_{h^{\prime}}$. In this case we write $g_{h}=f_{h}|_{W_{h^{\prime}}}$.

A vertex $i\in Q_{0}$ is called a source vertex if it is not a target of any
arrow in $Q_{1}$, and $i$ is called a sink vertex if it is not a source of any
arrow in $Q_{1}$. We assume that $Q$ has no sink vertices and no source
vertices. If $Q$ does have sink vertices or source vertices, we will work over
an extended quiver $\overline{Q}$, which is constructed as follows:

  1. 1.

add a new vertex labelled as $n+1$,

  2. 2.

add a new arrow from vertex $n+1$ to each source vertex,

  3. 3.

add a new arrow from each sink vertex to vertex $n+1$.

In this setting, every representation of $Q$ can be extended to a
representation of $\overline{Q}$ by adding the zero vector space to vertex
$n+1$ and the zero map to every arrow connected to vertex $n+1$.

Let $N=(W_{i},f_{h}|W_{h^{\prime}})_{i\in Q_{0},h\in Q_{1}}$ be a
subrepresentation of $M=(V_{i},f_{h})_{i\in Q_{0},h\in Q_{1}}$, we define an
operator $im^{-}$ as follows:

$im^{-}(N):=(U_{i},f_{h}|{U_{h^{\prime}}})_{i\in Q_{0},h\in Q_{1}}\text{ where
}U_{i}=\bigcap_{h^{\prime}=i}f^{-1}_{h}(W_{h^{\prime\prime}})\subset
V_{i}\text{ for }i\in Q_{0}.$

Moreover, we let

$im^{-1}(N)=im^{-}(N)\text{ and }im^{-i}(N)=im^{-}(im^{-(i-1)}(N))\text{ for
}i>1.$

Note that $im^{-}(N)$ is always a subrepresentation of $M$ and $N\subset
im^{-}(N)$. Thus we have the following chain of subrepresentations of $M$:

$0\subset im^{-1}(0)\subset im^{-2}(0)\subset\cdots\subset
im^{-s}(0)\subset\cdots$ (2)

$M$ is called a nilpotent representation if $im^{-s}(M)=M$ for some
$s\in\mathbb{N}$ and $M$ is called a monomorphic representation if
$im^{-1}(0)=0$. Note that $M$ is nilpotent if and only there exists
$s\in\mathbb{N}$ such that each vector space $V_{i}$ (for $i\in Q_{0}$) is
mapped to the zero vector space by any path starting from vertex $i$ as long
as the length of the path is greater $s$. Thus, the definition of nilpotency
here is equivalent to the one given by Bozec, Schiffmann & Vasserot [1]. Note
that if $Q$ has no oriented cycles, then the only nilpotent representation of
$Q$ is the zero representation.

We define another operator $im^{+}$ as follows:

$im^{+}(N):=(U_{i},f_{h}|U_{h^{\prime}})_{i\in Q_{0},h\in Q_{1}}\text{ where
}U_{i}=\sum_{h^{\prime\prime}=i}f_{h}(W_{h^{\prime}})\subset W_{i}\text{ for
}i\in Q_{0}.$

Moreover, we let

$im^{+1}(N)=im^{+}(N)\text{ and }im^{+i}(N)=im^{+}(im^{+(i-1)}(N))\text{ for
}i>1.$

Note that $im^{+}(N)$ is always a subrepresentation of $M$ and
$im^{+}(N)\subset N$. Thus we have the following chain of subrepresentations
of $M$:

$M\supset im^{+1}(M)\supset im^{+2}(M)\supset\cdots\supset
im^{+s}(M)\supset\cdots$ (3)

$M$ is called an epimorphic represenation if $im^{+1}(M)=M$. Note that $M$ is
nilpotent if there exists $s\in\mathbb{N}$ such that $im^{+s}(M)=0$.

Let $\mathrm{mod}_{k}(Q)$ be the category of representations of $Q$ over $k$
and let $\mathrm{mod}^{n}_{k}(Q)$ (resp. $\mathrm{mod}^{m}_{k}(Q)$,
$\mathrm{mod}^{e}_{k}(Q)$) be the subcategory consisting of all nilpotent
(resp. monomorphic, epimorphic) representations. All three subcategories are
closed under direct summands and extensions. Let
$\underline{\mathrm{mod}}_{k}(Q)$ (resp.
$\underline{\mathrm{mod}}^{n}_{k}(Q)$, $\underline{\mathrm{mod}}^{m}_{k}(Q)$,
$\underline{\mathrm{mod}}^{e}_{k}(Q)$) denote the set of isomorphism classes
in $\mathrm{mod}_{k}(Q)$(resp. $\mathrm{mod}^{n}_{k}(Q)$,
$\mathrm{mod}^{m}_{k}(Q)$, $\mathrm{mod}^{e}_{k}(Q)$). The isomorphism class
of an object $M$ in a category is denoted by $[M]$ and the cardinality of a
set $X$ is denoted by $|X|$.

For any finite field $k$, Ringel [10] defines an associative algebra
$\mathcal{H}(Q)$ known as the Ringel-Hall algebra associated with $Q$. As a
$\mathbb{Q}$-vector space, $\mathcal{H}(Q)$ has a basis
$\\{[M]\,|\,[M]\in\underline{\mathrm{mod}}_{k}(Q)\\}$, and the multiplication
of two basis elements is given by:

$[M]\circ[N]:=\sum_{[X]}g^{X}_{MN}[X],$

where $g^{X}_{MN}$ is the number of subrepresentations $U$ of $X$ such that
$U\cong N$ and $X/U\cong M$.

Ringel [10] proved that when $Q$ is of finite representation type (precisely
when the underlying graph of $Q$ is a Dynkin diagram by Gabriel [2]) the
twisted Ringel-Hall algebra is isomorphic to the quantized enveloping algebra
of the positive part of the semisimple Lie algebra associated with $Q$. For
general quivers, Green [3] proved that the composition subalgebra of the
twisted Ringel-Hall algebra is isomorphic to the quantized enveloping algebra
of the positive part of the Kac-Moody algebra associated with $Q$.

Let $X_{1},\cdots,X_{n}$ be $n$ non-commuting indeterminates, and let
$X^{v}=\prod_{i=1}^{n}X_{i}^{v_{i}}$ for any vector
$v=(v_{1},\cdots,v_{n})\in\mathbb{N}^{n}$. Let $\mathcal{T}(Q)$ be the torus
algebra associated with $Q$. Thus as a $\mathbb{Q}$-vector space,
$\mathcal{T}(Q)$ has a basis $\\{X^{v}|\,v\in\mathbb{N}^{n}\\}$, and the
multiplication of two basis elements is defined by the following rule:

$X^{v}\circ X^{w}=q^{-\langle v,w\rangle}X^{v+w}\,\text{ for all
}v,w\in\mathbb{N}^{n}.$ (4)

Reineke’s integration map is a $\mathbb{Q}$-linear map from $\mathcal{H}(Q)$
to $\mathcal{T}(Q)$ defined as follows:

$\delta:[M]\,\,\to\,\,\frac{1}{|\mathrm{Aut}(M)|}X^{\dim M},$ (5)

where $\mathrm{Aut}(M)$ is the automorphism group of $M$. The Riedtmann-Peng
formula ([9][6]) can be stated as:

$g^{X}_{MN}=\frac{|\mathrm{Aut}(X)|\cdot|\mathrm{Ext}^{1}(M,N)_{X}|}{|\mathrm{Hom}(M,N)|\cdot|\mathrm{Aut}(M)|\cdot|\mathrm{Aut}(N)|},$

where $\mathrm{Ext}^{1}(M,N)_{X}$ is the set of extension classes
corresponding to short exact sequences with middle term isomorphic to $X$. The
Riedtmann-Peng formula implies that Reineke’s integration map $\delta$ is a
homomorphism of algebras (Lemma 3.3 of Reineke [8]). Reineke’s integration map
can be naturally extended to a homomorphism from the completion algebra of
$\mathcal{H}(Q)$ to the completion algebra of $\mathcal{T}(Q)$. In these two
completion algebras, the sum and product of two infinite series are well
defined because both algebras are graded by dimension vectors.

Let $\mathbb{Q}(t)$ be the field of rational functions in $t$ over
$\mathbb{Q}$ and let $\mathcal{T}_{t}(Q)$ be the quantum torus over
$\mathbb{Q}(t)$ associated with $Q$. Thus as a $\mathbb{Q}(t)$-vector space,
$\mathcal{T}_{t}(Q)$ has a basis $\\{X^{v}|\,v\in\mathbb{N}^{n}\\}$ and the
multiplication of two basis elements is defined by the following rule:

$X^{v}\circ X^{w}=t^{-\langle v,w\rangle}X^{v+w}\,\text{ for all
}v,w\in\mathbb{N}^{n}.$ (6)

In the next section, we show that each representation of $Q$ is a unique
extension of a monomorphic representation by a nilpotent representation, and
it is also a unique extension of a nilpotent representation by an epimorphic
representation. The uniquenesses of those extensions yield two identities in
the Ringel-Hall algebra involving the numbers of representations of $Q$ over
finite fields. These identities are then translated into identities in the
quantum torus by Reineke’s integration map. In the last section we present two
identities which count the isomorphism classes of absolutely simple (resp.
indecomposable) conservative representations of $Q$ over finite fields.

##  2 Two Identities in Quantum Torus

Each representation $M=(V_{i},f_{h})_{i\in Q_{0},h\in Q_{1}}$ gives rise to
two linear maps $\sigma_{i}$ and $\tau_{i}$ for each $i\in Q_{0}$ as follows:

$\displaystyle\sigma_{i}:V_{i}$
$\displaystyle\,\,\to\,\,\bigoplus_{h^{\prime}=i}V_{h^{\prime\prime}},$
$\displaystyle{}\\!\\!\\!\\!\\!\\!\tau_{i}:\bigoplus_{h^{\prime\prime}=i}V_{h^{\prime}}$
$\displaystyle\,\,\to\,\,V_{i},$ $\displaystyle v_{i}$
$\displaystyle\,\,\mapsto\,\,(f_{h}(v_{i}))_{h^{\prime}=i},$
$\displaystyle{}(v_{h^{\prime}})$
$\displaystyle\,\,\mapsto\,\,\sum_{h^{\prime\prime}=i}f_{h}(v_{h^{\prime}}).$

The following lemma is a direct consequence of the definitions of monomorphic
and epimorphic representations.

######  Lemma 2.1.

Let $M=(V_{i},f_{h})_{i\in Q_{0},h\in Q_{1}}$ be a representation of $Q$. Then
we have

  1. 1.

$M$ is monomorphic if and only $\sigma_{i}$ is injective for all $i\in Q_{0}$.

  2. 2.

$M$ is epimorphic if and only $\tau_{i}$ is surjective for all $i\in Q_{0}$,

######  Lemma 2.2.

For any representation $M$ of $Q$ over a field $k$, there exists a unique
maximal nilpotent subrepresentation $N\subset M$. Moreover, $M/N$ is
monomorphic.

###### Proof.

Since $M$ is finite dimensional, chain (2) must be stable in finite steps,
i.e., there exists an integer $s$ such that $im^{-t}(0)=im^{-s}(0)$ for all
$t>s$. We claim that $im^{-s}(0)$ is the unique maximal nilpotent
subrepresentation of $M$.

$im^{-s}(0)$ is obviously nilpotent. Let
$N=(W_{i},f_{h}|_{W_{h^{\prime}}})_{i\in Q_{0},h\in Q_{1}}$ be any nilpotent
subrepresentation of $M$, we claim that $N\subset im^{-s}(0)$. Since $N$ is
nilpotent, there exists $r\in\mathbb{N}$ such that each vector space $W_{i}$
(for $i\in Q_{0}$) is mapped to the zero vector space by any path starting
from vertex $i$ as long as the length of the path is greater than $r$. This
implies that $N\subset im^{-(s+r)}(0)$. Since $im^{-(s+t)}(0)=im^{-s}(0)$, we
have $N\subset im^{-s}(0)$.

To show that $M/im^{-s}(0)$ is monomorphic, we only need to show the
$\sigma_{i}$ is injective for each $i\in Q_{0}$ when $\sigma_{i}$ is acting on
vector spaces in $M/im^{-s}(0)$, which is an easy consequence of the
definition of $im^{-}$.

∎

######  Lemma 2.3.

For any representation $M$ of $Q$ over a field $k$, there exists a unique
maximal epimorphic subrepresentation $E\subset M$. Moreover, $M/E$ is
nilpotent.

###### Proof.

Since $M$ is finite dimensional, chain (3) must be stable in finite steps,
i.e., there exists an integer $s$ such that $im^{+t}(M)=im^{+s}(M)$ for all
$t>s$. We claim that $im^{+s}(M)$ is the unique maximal epimorphic
subrepresentation of $M$.

$im^{+s}(M)$ is obviously epimorphic. Let $E$ be any epimorphic
subrepresentation of $M$, then $im^{+1}(E)=E$ and hence $im^{+s}(E)=E$.
$E\subset M$ implies that $im^{+s}(E)\subset im^{+s}(M)$ and hence $E\subset
im^{+s}(M)$. Thus $im^{+s}(M)$ is the unique maximal epimorphic
subrepresentation. $M/im^{+s}(M)$ is nilpotent because
$im^{+i}(M)/im^{+(i+1)}(M)$ is nilpotent for all $i\geq 1$. ∎

In what follows, let $k$ be the finite field $\mathbb{F}_{q}$ with $q$
elements, where $q$ is a prime power. Lemma 2.2 and Lemma 2.3 implies the
following theorem.

######  Theorem 2.1.

The following identities hold in the completion algebra of the Ringel-Hall
algebra $\mathcal{H}(Q)$:

$\displaystyle\sum_{[M]\in\underline{\mathrm{mod}}_{k}(Q)}[M]$
$\displaystyle=\Bigg{(}\sum_{[M]\in\underline{\mathrm{mod}}_{k}^{m}(Q)}[M]\Bigg{)}\circ\Bigg{(}\sum_{[M]\in\underline{\mathrm{mod}}_{k}^{n}(Q)}[M]\Bigg{)},$
$\displaystyle\sum_{[M]\in\underline{\mathrm{mod}}_{k}(Q)}[M]$
$\displaystyle=\Bigg{(}\sum_{[M]\in\underline{\mathrm{mod}}_{k}^{n}(Q)}[M]\Bigg{)}\circ\Bigg{(}\sum_{[M]\in\underline{\mathrm{mod}}_{k}^{e}(Q)}[M]\Bigg{)}.$

Applying Reineke’s integration map $\delta$ on both sides of the identities in
Theorem 2.1, it yields the following results.

######  Corollary 2.1.

The following identities hold in the completion algebra of the torus algebra
$\mathcal{T}(Q)$:

$\displaystyle\sum_{v\in\mathbb{N}^{n}}\frac{|\mathrm{R}(v,\mathbb{F}_{q})|}{|\mathrm{GL}(v,\mathbb{F}_{q})|}X^{v}$
$\displaystyle=\left(\sum_{v\in\mathbb{N}^{n}}\frac{|\mathrm{M}(v,\mathbb{F}_{q})|}{|\mathrm{GL}(v,\mathbb{F}_{q})|}X^{v}\right)\circ\left(\sum_{v\in\mathbb{N}^{n}}\frac{|\mathrm{N}(v,\mathbb{F}_{q})|}{|\mathrm{GL}(v,\mathbb{F}_{q})|}X^{v}\right),$
$\displaystyle\sum_{v\in\mathbb{N}^{n}}\frac{|\mathrm{R}(v,\mathbb{F}_{q})|}{|\mathrm{GL}(v,\mathbb{F}_{q})|}X^{v}$
$\displaystyle=\left(\sum_{v\in\mathbb{N}^{n}}\frac{|\mathrm{N}(v,\mathbb{F}_{q})|}{|\mathrm{GL}(v,\mathbb{F}_{q})|}X^{v}\right)\circ\left(\sum_{v\in\mathbb{N}^{n}}\frac{|\mathrm{E}(v,\mathbb{F}_{q})|}{|\mathrm{GL}(v,\mathbb{F}_{q})|}X^{v}\right).$

If $Q$ has source vertices or sink vertices, then Corollary 2.1 holds for
$\overline{Q}$. By taking $X^{n+1}=0$ we recover those two identities for $Q$.

For any dimension vector $v=(v_{1},\cdots,v_{n})\in\mathbb{N}^{n}$, let
$v\cdot v=\sum_{i=1}^{n}v_{i}^{2}\in\mathbb{N}$,
$v_{\bullet}=(v_{\bullet}^{i})\in\mathbb{N}^{n}$ with
$v_{\bullet}^{i}=\sum_{h^{\prime}=i}v_{h^{\prime\prime}}$,
${}_{\bullet}v=(_{\bullet}v^{i})\in\mathbb{N}^{n}$ with
${}_{\bullet}v^{i}=\sum_{h^{\prime\prime}=i}v_{h^{\prime}}$, and
$\mathrm{GL}(v,\mathbb{F}_{q})=\prod_{i=1}^{n}\mathrm{GL}(v_{i},\mathbb{F}_{q})$,
where $\mathrm{GL}(v_{i},\mathbb{F}_{q})$ is the General Linear Group of order
$v_{i}$ over $\mathbb{F}_{q}$. Let
$\mathit{gl}(s,t)=\prod_{i=0}^{s-1}(t^{s}-t^{i})\in\mathbb{Q}(t)$ for
$s\in\mathbb{N}$ and $\mathit{gl}(v,t)=\prod_{i=1}^{n}gl(v_{i},t)$, then we
have $|\mathrm{GL}(v,\mathbb{F}_{q})|=\mathit{gl}(v,q)$.

For any given dimension vector $v\in\mathbb{N}^{n}$, let
$\mathrm{R}(v,\mathbb{F}_{q})$ (resp. $\mathrm{N}(v,\mathbb{F}_{q})$,
$\mathrm{M}(v,\mathbb{F}_{q})$, $\mathrm{E}(v,\mathbb{F}_{q})$) denote the set
of all representations (resp. nilpotent, monomorphic and epimorphic
representations) of $Q$ over $\mathbb{F}_{q}$ with dimension vector $v$.

Define the following rational functions in $\mathbb{Q}(t)$:

$\displaystyle\mathit{r}(v,t)=t^{v\cdot v-\langle v,v\rangle},$
$\displaystyle\mathit{m}(v,t)=t^{\langle
v,v\rangle}\frac{\mathit{gl}(v_{\bullet},t)}{\mathit{gl}(v_{\bullet}-v,t)}\,\,\text{
if }v_{\bullet}\geq v\text{, else }\mathit{m}(v,t)=0,$
$\displaystyle\mathit{e}(v,t)=t^{\langle
v,v\rangle}\frac{\mathit{gl}(_{\bullet}v,t)}{\mathit{gl}(_{\bullet}v-v,t)}\,\,\text{
if }_{\bullet}v\geq v\text{, else }\mathit{e}(v,t)=0.$

Then we have $|\mathrm{R}(v,\mathbb{F}_{q})|=\mathit{r}(v,q)$,
$|\mathrm{M}(v,\mathbb{F}_{q})|=\mathit{m}(v,q)$ and
$|\mathrm{E}(v,\mathbb{F}_{q})|=\mathit{e}(v,q)$.

Define the rational function $\mathit{n}(v,t)$ as follows:

$\frac{\mathit{n}(v,t)}{\mathit{gl}(v,t)}=\sum_{v^{*}}t^{-\sum_{k<l}\langle
v^{(k)},v^{(l)}\rangle}\prod_{k\geq
1}\frac{H(v^{(k)},v^{(k+1)},t^{-1})}{\mathit{gl}(v^{(k)},t)},$ (7)

where the sum ranges over all tuples $v^{*}=(v^{(1)},v^{(2)},\dots)$ of non-
zero dimension vectors such that $\sum_{k\geq 1}v^{(k)}=v$ and
${}_{\bullet}v^{(k)}\geq v^{(k+1)}$ for $k\geq 1$, and where the function $H$
is defined by:

$H(v,w,t^{-1})=t^{(_{\bullet}v-w)\cdot(_{\bullet}v-w)-_{\bullet}v\cdot_{\bullet}v}\frac{\mathit{gl}(_{\bullet}v,t)}{\mathit{gl}(_{\bullet}v-w,t)}.$
(8)

Then we have $|\mathrm{N}(v,\mathbb{F}_{q})|=\mathit{n}(v,q)$, which is due to
Bozec, Schiffmann & Vasserot [1].

For the Jordan quiver i.e., the quiver with one vertex and one loop, the
identities in Corollary 2.1 are reduced to the following identity which is
trivial:

$\sum_{n=0}^{\infty}\frac{q^{n^{2}}X^{n}}{\prod_{i=0}^{n-1}(q^{n}-q^{i})}=\left(\sum_{n=0}^{\infty}\frac{q^{n^{2}-n}X^{n}}{\prod_{i=0}^{n-1}(q^{n}-q^{i})}\right)\cdot\left(\sum_{n=0}^{\infty}X^{n}\right).$
(9)

For the wild quiver which has one vertex and 2 loops, the first few terms of
the identities in Corollary 2.1 are as follows:

$\displaystyle R(1,q)=q^{2},$ $\displaystyle R(2,q)=q^{8},$ $\displaystyle
R(3,q)=q^{18},$ $\displaystyle N(1,q)=1,$ $\displaystyle
N(2,q)=q^{3}+q^{2}-q,$ $\displaystyle N(3,q)=q^{9}+2q^{8}-q^{6}-2q^{5}+q^{3},$
$\displaystyle M(1,q)=q^{2}-1,$ $\displaystyle M(2,q)=q^{8}-q^{5}-q^{4}+q,$
$\displaystyle M(3,q)=q^{18}-q^{14}-q^{13}-q^{12}+q^{9}+q^{8}+q^{7}-q^{3},$
$\displaystyle E(1,q)=q^{2}-1,$ $\displaystyle E(2,q)=q^{8}-q^{5}-q^{4}+q,$
$\displaystyle E(3,q)=q^{18}-q^{14}-q^{13}-q^{12}+q^{9}+q^{8}+q^{7}-q^{3}.$

Because $|\mathrm{R}(v,\mathbb{F}_{q})|$, $|\mathrm{N}(v,\mathbb{F}_{q})|$,
$|\mathrm{E}(v,\mathbb{F}_{q})|$, $|\mathrm{M}(v,\mathbb{F}_{q})|$ and
$|\mathrm{GL}(v,\mathbb{F}_{q})|$ are all polynomial functions in $q$ with
rational coefficients and there are infinitely many prime numbers, the two
identities in Corollary 2.1 must also hold in the quantum torus.

######  Corollary 2.2.

The following identities hold in the completion algebra of the quantum torus
$\mathcal{T}_{t}(Q)$:

$\displaystyle\sum_{v\in\mathbb{N}^{n}}\frac{r(v,t)}{{gl}(v,t)}X^{v}$
$\displaystyle=\left(\sum_{v\in\mathbb{N}^{n}}\frac{m(v,t)}{{gl}(v,t)}X^{v}\right)\circ\left(\sum_{v\in\mathbb{N}^{n}}\frac{n(v,t)}{{gl}(v,t)}X^{v}\right),$
$\displaystyle\sum_{v\in\mathbb{N}^{n}}\frac{r(v,t)}{{gl}(v,t)}X^{v}$
$\displaystyle=\left(\sum_{v\in\mathbb{N}^{n}}\frac{n(v,t)}{{gl}(v,t)}X^{v}\right)\circ\left(\sum_{v\in\mathbb{N}^{n}}\frac{e(v,t)}{{gl}(v,t)}X^{v}\right).$

##  3 Numbers of absolutely indecomposable conservative representations

A representation of $Q$ is called conservative if it is monomorphic and
epimorphic at the same time. The subcategory of conservative representations
is denoted by $\mathrm{mod}_{k}^{c}(Q)$. It is closed under direct summands
and extensions. A representation of $Q$ over a field is called absolutely
indecomposable if it is still indecomposable when the ground field is extended
to its algebraic closure. Let $c(v,q)$ be the number of conservative
representations of $Q$ over $\mathbb{F}_{q}$ with dimension
$v\in\mathbb{N}^{n}$ and let $a(v,q)$ (resp. $s(v,q)$) be the number of
isomorphism classes of absolutely indecomposable conservative (resp.
absolutely simple conservative) representations of $Q$ over $\mathbb{F}_{q}$
with dimension $v$.

Following the methods of Mozgovoy and Reineke [8][7], we have the following
result.

######  Corollary 3.1.

Assuming that $c(v,q)$ is a polynomial function in $q$ with rational
coefficients for all $v\in\mathbb{N}^{n}$, we have the following identity in
the completion algebra of the quantum torus $\mathcal{T}_{q}(Q)$:

$\Bigg{(}\sum_{v\in\mathbb{N}^{n}}\frac{c(v,q)}{\mathit{gl}(v,q)}X^{v}\Bigg{)}\\!\circ\mathrm{Exp}\Bigg{(}\sum_{v\in\mathbb{N}^{n}\backslash\\{0\\}}\\!\frac{s(v,q)}{1-q}X^{v}\Bigg{)}=1,\\\
$ (10)

where $\mathrm{Exp}$ is the plethystic exponential map (see [7] for details).

Following the methods of Hua [4], we have the following result.

######  Corollary 3.2.

Assuming that $c(v,q)$ is a polynomial function in $q$ with rational
coefficients for all $v\in\mathbb{N}^{n}$, we have the following identity in
the formal power series ring $\mathbb{Q}(q)[[X_{1},\cdots,X_{n}]]$:

$1+\sum_{v^{*}}\\!\Bigg{(}\prod_{s\geq 1}\frac{q^{\langle
v^{s}\\!,\,v^{s}\rangle}}{q^{\langle\beta^{s}\\!,\,\beta^{s}\rangle}}\frac{c(v^{s},q)}{\mathit{gl}(v^{s},q)}X^{sv^{s}}\Bigg{)}=\,\mathrm{Exp}\Bigg{(}\sum_{v\in\mathbb{N}^{n}\backslash\\{0\\}}\\!\frac{a(v,q)}{q-1}X^{v}\Bigg{)},$
(11)

where the sum runs over all tuples of dimension vectors
$v^{*}=(v^{1},v^{2},\cdots,v^{r})$ such that $r\geq 1$,
$v^{i}\in\mathbb{N}^{n}$ for $1\leq i\leq r$ and $v^{r}\neq 0$, and
$\beta^{s}=\sum_{i\geq s}v^{i}$ for $1\leq s\leq r$.

Consequently, $s(v,q)$ and $a(v,q)$ are polynomials in $q$ with rational
coefficients. Elementary calculations suggest that their coefficients are
always integers. Using the classification of indecomposable representations of
the Kronecker quiver, one can verify that the assumptions in Corollary 3.1 and
3.2 are ture for the wild quiver which has one vertex and 2 loops.

Acknowledgments

I would like to thank Yingbo Zhang, Bangming Deng and Xueqing Chen for their
valuable discussions and comments.

## References

  * [1] T. Bozec, O. Schiffmann & E. Vasserot, _On the numbers of points of nilpotent quiver varieties over finite fields_ , Annales scientifiques de l’Ecole normale supérieure, 2020, ff10.24033/asens.2452ff. ffhal-01709055.
  * [2] P. Gabriel, _Unzerlegbare Darstellungen. I_ , Manuscripta Math. 6 (1972) 71–103; correction, ibid. 6, 309 (1972).
  * [3] J. A. Green, _Hall algebras, hereditary algebras and quantum groups_ , Invent. Math. 120, 361-377 (1995).
  * [4] J. Hua, _On Numbers of Semistable Representations of Quivers over Finite Fields_ , arXiv:2109.08905 (2021).
  * [5] V. G. Kac, _Root Systems, Representations of Quivers and Invariant Theory_ , Lecture Notes in Mathematics 996, pp. 74–108, Springer-Verlag, Berlin (1983).
  * [6] L. Peng, _Lie algebras determined by finite Auslander-Reiten quivers_ , Comm. in Alg., 26, no.9, 2711–2725 (1998).
  * [7] S. Mozgovoy, & M. Reineke, _On the number of stable quiver representations over finite fields_ , Journal of Pure and Applied Algebra 213, 430-439 (2009)
  * [8] M. Reineke, _Counting Rational Points of Quiver Moduli_ , International Mathematical Research Notices 2006, pp. 1-19 (2006).
  * [9] C. Riedtmann, _Lie algebras generated by indecomposables_ , Journal of Algebra 170, no. 2, pp. 526–546 (1994).
  * [10] C. M. Ringel, _Hall algebras and quantum groups_ , Inv. Math. 101, pp. 583-591 (1990).

Email address<EMAIL_ADDRESS>


Uncertainty and Sensitivity Analyses Methods for Agent-Based Mathematical
Models: An Introductory Review

Sara Hamis1,2, Stanislav Stratiev3, Gibin G Powathil2.

1 School of Mathematics and Statistics, University of St Andrews, St Andrews
KY16 9SS, Scotland.

2 Department of Mathematics, College of Science, Swansea University, Swansea,
SA1 8EN, United Kingdom.

3 Department of Physics, College of Science, Swansea University, Swansea, SA2
8PP, United Kingdom.

###### Abstract

Multiscale, agent-based mathematical models of biological systems are often
associated with model uncertainty and sensitivity to parameter perturbations.
Here, three uncertainty and sensitivity analyses methods, that are suitable to
use when working with agent-based models, are discussed. These methods are
namely Consistency Analysis, Robustness Analysis and Latin Hypercube Analysis.
This introductory review discusses origins, conventions, implementation and
result interpretation of the aforementioned methods. Information on how to
implement the discussed methods in MATLAB is included.

##  1 Introduction

Mathematical models of biological systems are abstractions of highly complex
reality. It follows that parameters used in such models often are associated
with some degree of uncertainty, where the uncertainty can be derived from
various origins. Epistemic uncertainty refers to uncertainty resulting from
limited knowledge about the biological system at hand, whilst aleatory
uncertainty stems from naturally occurring stochasticity, intrinsic to
biological systems [23, 1, 15]. Model parameters may thus be naturally
stochastic, theoretically unknown, and unfeasible or impossible to measure
precisely (or at all). Further magnifying the contributions of uncertainty in
mathematical models of biological systems, in particular, is the fact that one
parameter in the mathematical model may correspond to a multitude of
underlying biological mechanisms and factors in the real, biological system.
This is especially true for minimal parameter models, i.e. mathematical models
that aspire to be as non-complex as possible whilst still capturing all
biological details of interest [8]. In this review, we focus our attention on
uncertainty and sensitivity analyses techniques that are suitable for use with
agent-based models. A mathematical, agent-based model comprises several
distinct agents that may interact with each other and their environment. In an
agent-based tumour model, for example, an agent typically corresponds to one
tumour cell or a group of tumour cells. This naturally allows for
heterogeneity amongst tumour cells, which is useful as tumour heterogeneity is
associated with many complications involved in modelling (and treating) solid
tumours. Accordingly, many modellers in the field of mathematical oncology
choose to work with agent-based models [24].

There already exist multiple method papers that describe how to perform
uncertainty and sensitivity analyses when working with agent-based models,
authors Alden et al. even provide a free R-based software package (Spartan
[1]) that enables the user to perform different such methods, including the
three methods discussed in this review. However, as these methods have been
developed across multiple research fields, both inside and outside of the
natural sciences, it is difficult to find one comprehensive review that
discusses not only how to perform these methods, but also where these methods
come from, and why certain conventions are proposed and/or used. To this end,
we have in this review gathered such information for three uncertainty and
sensitivity analyses techniques, namely Consistency Analysis, Robustness
Analysis and Latin Hypercube Analysis. Our aim is that this will allow the
reader to better evaluate uncertainty and sensitivity analyses presented by
other authors, and encourage the reader to consider performing these methods
when suitable.

In order to understand the impact that parameter uncertainty and parameter
perturbations have on results produced by a mathematical model, uncertainty
and sensitivity analyses can be used. A mathematical model that comprises a
set of uncertain model parameters (or inputs), is able to produce a range of
possible responses (or outputs). Uncertainty analysis assesses the range of
these outputs overall, and provides information regarding how certain (or
uncertain) we should be with our model results, and the conclusions that we
draw from them [2]. Sensitivity analysis describes the relationship between
uncertainty in inputs and uncertainty in outputs. It can be used to identify
which sources of input uncertainty (i.e. which model parameters) significantly
influence the uncertainty in the output and, equally importantly, which do not
[2]. Assessing how sensitive the output is to small input perturbations is a
healthy way to scrutinise our mathematical models [14]. Moreover, for a well-
formulated model, knowledge regarding how input uncertainty influences output
uncertainty can yield insight into the biological system that has not yet been
empirically observed [1]. Furthermore, if the uncertainty in some input
parameter is shown to not affect output uncertainty, the modeller may consider
fixing that parameter, and thus reducing model complexity in accordance with a
minimal-parameter modelling approach. In local sensitivity analysis
techniques, model parameters (inputs) are perturbed one at a time whilst other
parameters remain fixed at their calibrated value. In global sensitivity
analysis techniques, all model parameters are simultaneously perturbed [4].

There exist several sensitivity and uncertainty analyses techniques, but here
we will focus on three such techniques that are suitable to use in conjunction
with agent-based mathematical models. These techniques are namely Consistency
Analysis, Robustness Analysis and Latin Hypercube Analysis, which all answer
important, and complementary, questions about mathematical models and their
corresponding in silico responses [23, 1].

Note that Consistency Analysis is only meaningful when analysing models with
stochastic variables.

The statistical techniques described in this review have been developed and
applied across multiple academic disciplines, both inside and outside of the
natural sciences. Consequently, terminology and notations vary in the
literature. The aim of this review is to combine pertinent literature from
various academic fields whilst keeping terminology and mathematical notations
consistent, unambiguous and tailored towards a mathematical and scientific
audience. Therefore, when needed, certain algorithms from the literature are
here reformulated into expressions that a mathematician would consider to be
conventional. This review is intended to provide gentle, yet comprehensive,
instructions to the modeller wanting to perform uncertainty and sensitivity
analyses on agent-based models. Thorough directions on how to perform
Consistency Analysis (Section 3), Robustness Analysis (Section 4) and Latin
Hypercube Sampling and Analysis (Section 5) are provided. Consistency Analysis
utilises the A-measure of stochastic superiority, which is therefore discussed
in Section 2. Throughout this review, we have included some historical
information that elucidates why certain statistical conventions are used. Each
section also contains pictorial, step-by-step instructions on how to perform
the aforementioned techniques. Worked examples of all methods discussed in
this review are provided in Section 6. These worked examples use in silico
data produced in a previous agent-based, multiscale, mathematical oncology
study [9].

#### Methods outside the scope of this review

Note that there exist other uncertainty and sensitivity analyses techniques,
suitable for agent-based models, that are outside the scope of this review.
For example, Bayesian inference is a statistical technique that uses Bayes
theorem, prior beliefs and observed data to infer input parameter values and
their associated uncertainties [12]. When this inference is difficult to
calculate, computational methods such as Approximate Bayesian Computation
(ABC) [13] or, if the model is highly computationally expensive, Approximate
Approximate Bayesian Computation (AABC) [3] can be used to perform inference
on input model parameters. The Sobol method is a variance-based sensitivity
analysis method that, using Monte-Carlo multidimensional integration, allows
for the evaluation of each individual parameter’s fractional contribution to
the output variance [31]. Another global, variance based sensitivity analysis
method that enables quantification of the input parameter’s fractional
contributions to the output variance is Fourier amplitude sensitivity testing
(FAST) which uses an underlying algorithm that involves Fourier decomposition
[26]. Although the Sobol method and FAST use different underlying techniques,
they both allow us to say that “input parameter $r_{i}$ contributes $R_{i}$
percent to the output variance”. Newer sensitivity analysis methods and tools
tend to take advantage of the current abundance of computing power. The
recently introduced MASSIVE (Massively parallel Agent-based Simulations and
Subsequent Interactive Visualization-based Exploration) methodology, for
example, combines parallel computing and interactive data visualisation to
produce a graphical user interface that provides an overview of input-output
relations for a broad input parameter range [22].

##  2 The A-measure of stochastic superiority

###  2.1 The Common Language Statistics

In 1992, McGraw and Wong introduced the common language statistics ($CL$) as
an intuitive way to compare two distributions of data [18]. The $CL$ was
initially introduced as a tool to compare data from normal distributions, but
was later on approximated for use on any continuous distributions. The $CL$
describes the probability that a random data sample from one of the
distributions is greater than a random data sample from the other
distribution. For example, if we have two continuous data distributions $B$
and $C$, and we are comparing the distributions with respect to some variable
$X$, then the $CL$ is simply given by

$CL_{BC}(X)=P(X_{B}>X_{C}),$ (1)

where standard probability notations have been used so that $P(X_{B}>X_{C})$
denotes the probability that a random data sample $X_{B}$ from distribution
$B$ is greater than a random data sample $X_{C}$ from distribution $C$ [18].
Thus the subscript of $X$ here signifies the distribution from which the data
sample $X$ was taken.

###  2.2 The A-measure of stochastic superiority

The $CL$ was developed to compare continuous data distributions, but Vargha
and Delaney [29] introduced the A-measure of stochastic superiority (or
A-measure for short) as a generalisation of the $CL$ that can directly be
applied to compare both continuous and discrete distributions of variables
that are at least ordinally scaled. When comparing two distributions $B$ and
$C$, with respect to the variable $X$, the A-measure $A_{BC}(X)$ is given by

$A_{BC}(X)=P(X_{B}>X_{C})+0.5P(X_{B}=X_{C}),$ (2)

where $P(X_{B}=X_{C})$ denotes the probability that a random data sample from
distribution $B$ is equal to a random data sample from distribution $C$. By
comparing Equations 1 and 2, it is clear that in the continuous case, where
$P(X_{B}=X_{C})=0$, the A-measure reduces to the $CL$.

If two distributions that are identical with respect to the variable $X$ are
compared, then $P(X_{B}>X_{C})=P(X_{C}>X_{B})$ and we say that the
distributions $B$ and $C$ are stochastically equal with respect to the
variable $X$. On the other hand, if $P(X_{B}>X_{C})>P(X_{C}>X_{B})$, then we
say that the distribution $B$ is stochastically greater than distribution $C$,
and accordingly, that distribution $C$ is stochastically smaller than
distribution $B$ [29]. If distribution $B$ is stochastically greater than
distribution $C$ with respect to the variable $X$, it simply occurs more often
that the sample $X_{B}$ is greater than the sample $X_{C}$ when two random
samples $X_{B}$ and $X_{C}$ are compared. These definitions of stochastic
relationships (stochastically equal to, stochastically greater than,
stochastically smaller than), used by Vargha and Delaney [29], amongst others,
are weaker than definitions used by some other authors, but sufficient and
appropriate for our current purposes: comparing distributions of discrete data
samples produced by in silico simulations based on stochastic, agent-based
mathematical models.

When comparing two samples $X_{B}$ and $X_{C}$, the possible outcomes are (i)
that $X_{B}$ is greater than $X_{C}$, (ii) that $X_{B}$ is equal to $X_{C}$
and (iii) that $X_{B}$ is smaller than $X_{C}$. These three possible outcomes
must sum up to one so that,

$P(X_{B}>X_{C})+P(X_{B}=X_{C})+P(X_{C}>X_{B})=1.$ (3)

In the continuous case, $P(X_{B}=X_{C})=0$ as previously stated, and thus it
follows that

$P(X_{C}>X_{B})=1-P(X_{B}>X_{C}),\quad\text{for continuous distributions,}$
(4)

and thus it suffices to know only one of the values $P(X_{B}>X_{C})$ or
$P(X_{C}>X_{B})$ in order to determine the stochastic relationship between the
distributions $B$ and $C$ with respect to $X$.

  * $\blacktriangleright$

For example: if $P(X_{B}>X_{C})=0.4$, then it is clear that
$P(X_{C}>X_{B})=0.6$ and thus that $P(X_{B}>X_{C})<P(X_{C}>X_{B})$, or
equivalently, that distribution $B$ is stochastically smaller than
distribution $C$.

However, in the discrete case, $P(X_{B}=X_{C})$ is not generally equal to zero
and therefore,

$P(X_{C}>X_{B})=1-P(X_{B}>X_{C})-P(X_{B}=X_{C})\quad\text{for discrete
distributions.}$ (5)

Consequently, one single value $P(X_{B}>X_{C})$ or $P(X_{C}>X_{B})$ alone can
generally not be used to determine the stochastic relationship between the
distributions $B$ and $C$.

  * $\blacktriangleright$

For example: if, again, $P(X_{B}>X_{C})=0.4$, it follows that
$P(X_{C}>X_{B})=0.6-P(X_{B}=X_{C})$. This does not give us enough information
to determine the stochastic relationship between the two distributions $B$ and
$C$.

In order to proceed to compare the distributions $B$ and $C$ in this case, the
stochastic difference $\delta$ is introduced, where $\delta$ is given by

$\delta=P(X_{B}>X_{C})-P(X_{C}>X_{B}),\quad\delta\in[-1,1].$ (6)

Via a linear transformation, the transformed stochastic difference,
$\delta^{\prime}\in[0,1]$, can be obtained using Equation 5 so that

$\delta^{\prime}=\frac{\delta+1}{2}=\frac{P(X_{B}>X_{C})-P(X_{C}>X_{B})+1}{2}=\\\
=\frac{P(X_{B}>X_{C})-\big{(}1-P(X_{B}>X_{C})-P(X_{B}=X_{C})\big{)}+1}{2}=\\\
=P(X_{B}>X_{C})+0.5P(X_{B}=X_{C})=A_{BC}(X),$ (7)

from which we can see that the A-measure, $A_{BC}(X)$ (Equation 2), measures
the stochastic difference between $P(X_{B}>X_{C})$ and $P(X_{C}>X_{B})$ under
a linear transformation [29].

In order to estimate the A-measure using samples from two distributions, the
point estimate of the A-measure, here denoted the $\hat{A}$-measure (with a
hat), is used. (In the Spartan package [1], this is referred to as the A test
score). For example, if we want to compare two discrete distributions $B$ and
$C$, where $B$ comprises $m$ data samples (of some variable $X$) so that
$B=\\{b_{1},b_{2},..,b_{m}\\}$ and $C$ comprises $n$ data samples (of some
variable $X$) so that $C=\\{c_{1},c_{2},..,c_{n}\\}$ then

$\hat{A}_{BC}(X)=\frac{\mathcal{\\#}(b_{i}>c_{j})}{mn}+0.5\thinspace\frac{\mathcal{\\#}(b_{i}=c_{j})}{mn},$
(8)

where $i=1,2,..,m$ and $j=1,2,..,n$ and $\mathcal{\\#}(\text{event})$ is the
‘counting function’ that simply denotes the number of times that a certain
event occurs when comparing all possible pairs of data samples
$(b_{i},c_{j})$. For clarity, Figure 1 provides an example of how the
$\hat{A}$-measure can be computed by simply counting events.

Figure 1: Using Equation 8 to compute the point estimate of the A-measure,
i.e. the $\hat{A}$-measure or $\hat{A}_{BC}$, of the two distributions of data
samples $B$ and $C$ of sizes $m$ and $n$ respectively.

Using more conventional mathematical notation, the $\hat{A}$-measure is given
by

$\hat{A}_{BC}(X)=\frac{1}{mn}{\sum}_{i=1}^{m}{\sum}_{j=1}^{n}H(b_{i}-c_{j}),$
(9)

where $H(x)$ is the Heaviside step function such that

$H(x)=\left\\{\begin{array}[]{ll}1\quad\text{ for $x>0$},\\\
\frac{1}{2}\quad\text{ for $x=0$},\\\ 0\quad\text{ for
$x<0$}.\end{array}\right.$ (10)

If $\hat{A}_{BC}(X)=0.5$, then the distributions $B$ and $C$ are
stochastically equal with respect to the variable $X$. The $\hat{A}$-measure
can thus be used to measure ‘how equal’ two discrete distributions $B$ and $C$
are, by assessing how much the $\hat{A}$-measure $(\in[0,1])$ deviates from
equality, i.e. the value $0.5$. The closer the $\hat{A}$-measure is to $0.5$,
the ‘more equal’ the two compared distributions are [29]. In many
applications, we are only interested in ‘how equal’ two distributions $B$ and
$C$ are, and it is not important which distribution is the stochastically
greater one. In such cases we are only interested in how much the
$\hat{A}$-measure deviates from stochastic equality (i.e. the value $0.5$) but
the direction is not important. Or in mathematical terms: the magnitude of the
difference between the $\hat{A}$-measure and stochastic equality is important
but the sign is not. The magnitudal (or scaled) $\hat{A}$-measure (or
$\hat{A}$-value), here denoted $\underline{\hat{A}}$ with an underscore,
ignores the sign of deviation from equality and is given by

$\underline{\hat{A}}=\begin{cases}\hat{A}_{BC}(X)&\text{ if
$\hat{A}_{BC}(X)\geq 0.5$,}\\\ 1-\hat{A}_{BC}(X)&\text{ if
$\hat{A}_{BC}(X)<0.5$}.\end{cases}$ (11)

The statistical significance is used to describe the effect of the stochastic
difference between two distributions $B$ and $C$. If two distributions $B$ and
$C$ are ‘fairly equal’ (i.e. if they yield an
$\underline{\hat{A}}_{BC}$-measure close to 0.5) then the statistical
significance is classified as small. The statistical significance is
classified using the magnitudal $\underline{\hat{A}}$-measure and, using
guidelines from Vargha and Delaney [29], the statistical significance is
classified to be small, medium or large with respect to $X$ according to the
following threshold values for $\underline{\hat{A}}_{BC}(X)$,

$\text{Statistical Significance}=\begin{cases}\text{small}&\text{if
$\underline{\hat{A}}_{BC}(X)\in[0.5,0.56]$},\\\ \text{medium}&\text{if
$\underline{\hat{A}}_{BC}(X)\in(0.56,0.64]$},\\\ \text{large}&\text{if
$\underline{\hat{A}}_{BC}(X)\in(0.64,0.71]$}.\end{cases}$ (12)

These threshold values (that might appear somewhat arbitrary) were first
introduced by psychologist and statistician Jacob Cohen [6, 7] in the 1960s
when comparing normal distributions, but then in terms of another statistical
measurement: the effect size (Cohen’s) d where

$\textbf{d}=\frac{\big{|}\text{(mean of population $B$) - (mean of population
$C$)}\big{|}}{\sigma},$ (13)

and $\sigma$ is the standard deviation of either $B$ or $C$ (as B and C here
are assumed to have the same standard deviation) [7, 25]. Omitting details
from statistics, a small d-value essentially corresponds to a big overlap
between distributions $B$ and $C$, whilst a large d-value corresponds to a
small overlap between distributions $B$ and $C$, as is illustrated in Figure
2. Cohen decided to use the threshold d-values for describing ‘small’,
‘medium’ and ‘large’ effect sizes to be 0.2, 0.5 and 0.8 respectively [7]. If
we hold on to the assumption that $B$ and $C$ are two normal distributions
with the same variability, and furthermore say that they contain the same
number of data samples, we can use measures of overlap to get a further ‘feel’
for the previously discussed effect sizes, as illustrated in Figure 2. Cohen’s
d value can also be converted into ‘the probability that a random data sample
$X_{B}$ from (normal) distribution $B$ is larger than a random data sample
$X_{C}$ from (normal) distribution $C$ [18], but that is exactly what the
$\hat{A}$-measure $\hat{A}_{BC}(X)$ measures! So this is where the threshold
values for the descriptors ‘small’, ‘medium’ and ‘large’ statistical
differences, listed in Equation 12, come from.

Figure 2: The small (left), medium (centre) and large (right) threshold values
for the scaled A-measure of stochastic superiority
($\underline{\hat{A}}_{BC}$) are based on Cohen’s d-values comparing two
normal distributions $B$ and $C$ with the same variance. The higher the
overlap between $B$ and $C$, the smaller the d-value, and the smaller the
$\underline{\hat{A}}_{BC}$-measure ($\underline{\hat{A}}_{BC}\in[0.5,1]$).

Now, Cohen motivated his choice of the d-value thresholds using a blend of
intuitive ‘everyday’ examples and mathematical reasoning [7]. However, he did
issue a warning regarding the fact that the threshold values should be
determined based on the research methodology at hand. Thus the modeller should
not blindly use Cohen’s suggested thresholds, but instead reason what
constitutes a small enough statistical significance in the study at hand. The
modeller must also decide how fine the data samples in the data distributions
should be before performing Consistency Analysis. In many applications, it is
likely that the amount of data samples required in order to achieve a small
statistical significance increases with the fineness of the data. Nonetheless,
scientific conventions are useful and thus in the remainder of this review we
will use the threshold values suggested by Cohen, as is done in other
mathematical biology studies [1].

##  3 Consistency Analysis

In silico simulations based on mathematical models with built-in stochasticity
will not produce the same output data every simulation run. Consistency
Analysis (also called aleatory analysis) is a stochastic technique that
answers the question: how many data samples do we need to produce in order to
mitigate uncertainty originating from intrinsic model stochasticity? In our
case, one data sample is the product of one in silico simulation, so an
equivalent question is: how many in silico simulations should we run before
describing our results in terms of, for example, average values, standard
deviations or similar?

Let us say that one in silico simulation produces one data sample of some
output response $X$. This data sample can for example correspond to ‘the
population size at time point $T$’, or something similar. It is up to the
modeller to identify and decide what the meaningful output response(s) should
be, and Consistency Analysis can be performed on multiple output responses at
multiple time steps, for comprehensiveness. Before we begin, note that when
performing Consistency Analysis, we always use the calibrated model
parameters.

The first step involved in performing Consistency Analysis is to produce
multiple distributions of data of various sizes. We say that a distribution
with $n$ data samples has a distribution size $n$, and the goal of Consistency
Analysis is to find the smallest $n$ value (here denoted $n^{*}$) that yields
a small stochastic significance. To do this, we create various distribution
groups that all contain 20 distributions each of some distribution size $n$,
as is shown in Step 1 in Section 3.1. Following the methodology described in
previous work by Alden et al., and the Spartan package that they developed
[1], we create one distribution group that contains 20 distributions of size
$n=1$, one distribution group that contains 20 distributions of size $n=5$ and
so on. Here, the $n$ values 1, 5, 50, 100 and 300 are evaluated [1] and thus
we must produce a total of $20\times(1+5+50+100+300)=9120$ in silico runs.
(Note that if the desired accuracy is not achieved for the highest
investigated $n$ value, here $n=300$, higher values of $n$ can be explored).

We here let a distribution $D_{n,k}$ denote the $k$th distribution of
distribution size $n$ so that

$D_{n,k}=\\{d_{n,k}^{1},d_{n,k}^{2},\>..\>,d_{n,k}^{n}\\}$ (14)

where $d_{n,k}^{h}$ is the the $h$th data sample in distribution $D_{n,k}$ and
$h=1,2,..,n$. The $\hat{A}$-measure resulting from comparing two distributions
$D_{n,k}$ and $D_{n,k^{\prime}}$ with respect to the variable $X$ is denoted
by $\hat{A}^{n}_{k,k^{\prime}}(X).$

Now, within every distribution-group, we compare the first distribution
($k=1$) to all other distributions ($k^{\prime}=2,3,..,20$) using the
$\hat{A}$-measure. This yields 19 $\hat{A}$-measures per distribution-group,
as is shown in Step 2 in Section 3.1. The maximum scaled $\hat{A}$-measure
with respect to $X$, occurring in a distribution-group $g$ that contains
distributions of size $n_{g}$, is denoted
$\underline{\hat{A}}^{n_{g}}_{max}(X)$. The smallest value $n_{g}$ for which
$\underline{\hat{A}}^{n_{g}}_{max}(X)\leq 0.56$ is denoted $n^{*}$. In other
words: $n^{*}$ corresponds to the smallest distribution size for which all of
the 19 computed $\hat{A}$-measures yield a small stochastic significance, as
is shown in Step 3 in Section 3.1. This answers the question that we set out
to answer via Consistency Analysis: $n^{*}$ data samples (or in silico runs)
are needed in order to mitigate uncertainty originating from intrinsic model
stochasticity. The procedure on how to perform Consistency Analysis is
outlined Section 3.1.

###  3.1 Quick Guide: Consistency Analysis

Here follows a quick guide for how to perform Consistency Analysis.

##  4 Robustness Analysis

Robustness Analysis answers the question: how robust are model responses to
local parameter perturbations? Robustness Analysis investigates if, and how,
perturbing the value of one input parameter significantly changes an output
$X$. Using the ${\hat{A}}$-measure, data distributions containing output data
produced by perturbed input parameters, are compared to a data distribution
containing output data produced by the calibrated input parameters. All
perturbed data distributions are here of size $n^{*}$, where $n^{*}$ is
decided in the Consistency Analysis process, previously described in Section
3, when analysing stochastic models.

Before commencing the Robustness Analysis, we must identify the uncertain
model parameters that we want to investigate the robustness of. We denote
these parameters $p^{i}$, where $i=1,2,..,q$, and thus we have a total of $q$
parameters whose robustness we will investigate. Now, as illustrated in Step 1
in Section 4.1, we let each such parameter $p^{i}$ be investigated at
$r(p^{i})$ different parameter values (including the calibrated value), and
thus we need to generate a total of $P$ distributions of sample size $n^{*}$
where

$P=\sum_{i=1}^{q}r(p^{i}).$ (15)

Note that the number of investigated parameter values, $r(p^{i})$, need not be
the same for every input parameter $p^{i}$. Investigated distributions of
sample size $n^{*}$ are here denoted $D_{n^{*},p^{i}_{j}}$, where $i=1,2,..,q$
denotes which parameter is being perturbed and $j=1,2,..,r(p^{i})$ denotes the
specific perturbation of parameter $p^{i}$. For some perturbation $j=C$, the
parameter value $p^{i}_{j}$ equals the calibrated value for input parameter
$p^{i}$. For each parameter that we are investigating, the $\hat{A}$-measure
is used to compare the calibrated distribution $D_{n^{*},p^{i}_{C}}$ to all
distributions $D_{n^{*},p^{i}_{j}}$. Note that when $j=C$, the calibrated
distribution is compared to itself and thus the $\hat{A}$-measure equals 0.5.
These $\hat{A}$-measures provide information regarding the statistical
significance, specifically if it can be described to be small, medium or large
under parameter perturbations. Plotting the corresponding $\hat{A}$-measure
over the parameter value $p^{i}_{j}$ for each parameter $p^{i}$, paints an
informative picture of local parameter robustness, as shown in Step 2, in
Section 4.1. Another descriptive way to demonstrate the influence that
parameter values $p^{i}_{j}$ have on some output response $X$ is to use
boxplots. As is illustrated in Step 3 in Section 4.1, boxplots can be used to
clearly show the median, different percentiles, and outliers of some data
distribution $D_{n^{*},p^{i}_{j}}$ as a function of the parameter value
$p^{i}_{j}$. The methodology to perform Robustness Analysis is outlined in
Section 4.1. Note that Robustness Analysis does not pick up on any non-linear
effects between an input parameter $p^{i}$ and an output $X$, that occur when
more than one model parameter is simultaneously perturbed [4]. Such effects
can however be identified using a global sensitivity analysis technique, such
as Latin Hypercube Analysis, as described in Section 5.

###  4.1 Quick Guide: Robustness Analysis

##  5 Latin Hypercube Sampling and Analysis

Latin Hypercube Analysis answers the question: how robust are model responses
to global parameter perturbations? Latin Hypercube Analysis is a type of
global sensitivity analysis that investigates the relationship between input
parameters and output responses when all input parameters are simultaneously
perturbed. The parameters that we want to perturb are (as in Section 4)
denoted $p^{i}$, where $i=1,2,..,q$. Thus the parameters
$p^{1},p^{2},...,p^{q}$ together span a parameter space of dimension $q$. It
is impossible to test every possible combination of input parameter values if
they are picked from continuous ranges. In fact, even if we select a finite
number of parameter values $r(p^{i})$ to test for each parameter $p^{i}$, or
if we pick discrete parameter values, comparing every possible combination of
parameter values may require us to produce an impractically large number of
simulation runs. Thus performing in silico simulations for all possible
combinations of input parameters will in many cases be at worst impossible,
and at best impractical. In order to circumvent this issue, Latin Hypercube
Sampling can be used [1]. It is a sampling technique that ensures
comprehensive testing coverage over the parameter space whilst keeping the
number of tested parameter combinations low enough to be applicable in
practice [19, 20]. After Latin Hypercube Sampling (Section 5.1), Latin
Hypercube Analysis (Section 5.2) is used in order to assess global
sensitivity.

###  5.1 Latin Hypercube Sampling

In the two-dimensional case, a Latin Square is an $\ell\times\ell$ square grid
containing $\ell$ (traditionally Latin, hence the name) different symbols such
that each symbol occurs exactly once in every row and exactly once in every
column [28], as illustrated in Figure 3. Analogously, in the Latin Hypercube
Sampling framework, consider two parameters $p^{1}$ and $p^{2}$, spanning a
parameter space of dimension $q=2$, where both $p^{1}$ and $p^{2}$ are
sectioned into $\ell$ intervals. We then pick $\ell$ combinations of input
parameter values (or sampling points) $(p^{1}_{j},p^{2}_{j})$, where
$j=1,2,...,\ell$, such that every $p^{1}$-interval is sampled from exactly
once and every $p^{2}$-interval is sampled from exactly once. Within the
parameter range of an interval, the sampled parameter value $p^{i}_{j}$ is
randomly selected (unless of course the interval contains only one possible
value $p^{i}_{j}$). Note that the $j$ index denotes the coordinate combination
that $p^{i}_{j}$ belongs to, not the interval from which the parameter value
$p^{i}_{j}$ was taken. Thus there is no condition demanding that the values
$p^{i}_{j}$ are ordered in a way such that
$p^{i}_{1}<p^{i}_{2}<...<p^{i}_{\ell}$.

The analogy between a Latin Square and Latin Hypercube Sampling from a two-
dimensional parameter space is illustrated in Figure 3. The Latin Square can
be extended to higher dimensions to form a Latin Cube (dimension = 3) or a
Latin Hypercube (dimension > 3) and, analogously, the two-dimensional sampling
space illustrated in Figure 3 can be extended to $q$ dimensions, spanned by
the input parameters $p^{1},p^{2},..,p^{q}$ [28].

Figure 3: Left: An $\ell\times\ell$ Latin Square in which each Latin symbol
occurs $\ell$ times, exactly once in each row and exactly once in each column.
Right (analogously): A two-dimensional parameter space spanned by the input
parameters $p^{1}$ and $p^{2}$ that are both sectioned into $\ell$ intervals.
Using Latin Hypercube sampling, $\ell$ parameter combinations
$(p^{1}_{j},p^{2}_{j})$ are sampled where $j=1,2,..,\ell$ and each
$p^{1}$-interval is sampled from exactly once and each $p^{2}$-interval is
sampled from exactly once.

.

For each parameter $p^{i}$, the total investigated parameter range is
$[min(p^{i}),\;max(p^{i})]$, where $min(p^{i})$ and $max(p^{i})$ respectively
denote the minimum and maximum values of $p^{i}$ to be investigated. Now each
parameter range $[min(p^{i}),\;max(p^{i})]$ is sectioned into $N$ intervals,
and we denote these intervals by
$u_{p^{i}}^{1},u_{p^{i}}^{2},....,u_{p^{i}}^{N}$. Note that all input
parameters $p^{i}$ must be sectioned into the same number of intervals. If the
intervals are of equal size, then the size of an interval, $w(p^{i})$, is

$w(p^{i})=\frac{max(p^{i})-min(p^{i})}{N}$ (16)

and the $r$th interval $u_{p^{i}}^{r}$ has a parameter range such that

$u_{p^{i}}^{r}=\big{[}min(p^{i})+w\cdot(r-1),\;min(p^{i})+w\cdot r\big{]}$
(17)

where $r=1,2,...,N$.

Note that there are more than one way to populate Latin symbols in a Latin
Square, this can be realised by regarding Figure 3 and noticing that the
A-symbols and the B-symbols cover the Latin Square in different ways.
Analogously, and by extension, there are multiples ways to populate sampling
coordinates in a Latin Hypercube Sampling framework. Some of these ways
provide better coverage of the parameter space than do others [28], but
details regarding such sampling-optimisation are outside the scope of this
review. Here, we use the built-in MATLAB function lhsdesign [17] to select
which parameter combinations to use according to a Latin Hypercube Sampling
approach, details about the implementation are available in the Appendix. Note
that, in our case, all $N$ intervals
$u_{p^{i}}^{1},u_{p^{i}}^{2},....,u_{p^{i}}^{N}$ for a parameter $p^{i}$ are
uniformly spaced, but the choice of spacing can be adjusted to the specific
application at hand [17].

Now let us address the choice of intervals $N$, as this is not
straightforward. Using the Latin Hypercube Sampling framework, every parameter
$p^{i}$, where $i=1,2,..,q$, is partitioned into $N$ intervals and,
consequently, $N$ combinations comprising $q$ parameter values are sampled and
tested. Compared to a small $N$ value, a large value of $N$ will provide more
data to use, and draw conclusions from, in the Latin Hypercube Analysis stage,
however, it will also increase the computational cost in the Latin Hypercube
Sampling stage. There is no strict rule for how to choose $N$, but suggested
values for $N$ in the literature are $N=2q$ for large values of $q$ (i.e.
high-dimensional parameter spaces) or $N=4q/3$ which has been described to be
‘usually satisfactory’ [16, 10]. Authors of the Spartan package use a lot
larger numbers in their provided examples [1]. In this example study, we
decide to use $N=100$ uniform intervals. At the end of the day, the choice of
$N$ is up to the modeller, who must outweigh the (computational) cost of
producing a large number of data samples, with the advantage of having a vast
amount of data, and thus plentiful information, in the analysis stage. Details
regarding quantitative choices of $N$ are outside the scope of this review.

###  5.2 Latin Hypercube Analysis

During the Latin Hypercube Sampling process, $N$ different points in the
$q$-dimensional parameter space spanned by the input parameters
$p^{1},p^{2},...,p^{q}$ are selected as sampling points, as shown in Step 1 in
Section 5.3. One such sampling point, $C_{j}$, can be described by its
coordinates in the parameter space so that
$C_{j}=(p^{1}_{j},p^{2}_{j},...,p^{q}_{j})$. Each sampling point $C_{j}$ is
used to generate $n^{*}$ output responses $X(C_{j})$, where $n^{*}$ is
determined using Consistency Analysis. Subsequently, the median output value,
here denoted $\underaccent{\tilde}{X}(C_{j})$, is computed for every $C_{j}$.
Now, our overall aim is to investigate the relationship between an input
parameter $p^{i}$ and an output response ${X}$. We investigate this input-
output relationship in two steps, one of which is qualitative and one of which
is quantitative. In the first and qualitative step, we produce two-dimensional
scatterplots in which median output data,

$\underaccent{\tilde}{X}(C_{1}),\underaccent{\tilde}{X}(C_{2}),...,\underaccent{\tilde}{X}(C_{N})=\underaccent{\tilde}{X}(p^{1}_{1},p^{2}_{1},...,p^{q}_{1}),\underaccent{\tilde}{X}(p^{1}_{2},p^{2}_{2},...,p^{q}_{2}),...,\underaccent{\tilde}{X}(p^{1}_{N},p^{2}_{N},...,p^{q}_{N}),$

are plotted over parameter values

$p^{i}_{1},p^{i}_{2},...,p^{i}_{N},$

for one of the input parameters $p^{i}$. We do this for every input parameter
$i=1,2,...,q$ and thus $q$ scatterplots are created. By simply visually
analysing the data in the scatterplots, we are able to make qualitative
observations regarding the relationship between the input and the output.
Examples of such observations are provided in Step 2 in Section 5.3

As a second step, we use a quantitative measure, such as the Pearson Product
Moment Correlation Coefficient (or the correlation coefficient for short), to
quantitatively describe the correlation between input parameters and output
responses, as done in Step 3 in Section 5.3. The correlation coefficient is
denoted $r$, where $r\in[-1,+1]$. It describes the linear association between
the input parameter and the output response in terms of both magnitude and
direction. A positive (linear) correlation between $p^{i}$ and
$\underaccent{\tilde}{X}(C_{j})$ means that if either the input value or the
output value increases, so does the other one, and thus $r$ is positive.
Conversely, a negative correlation means that if either $p^{i}$ or
$\underaccent{\tilde}{X}(C_{j})$ increases, the other one decreases, and thus
$r$ is negative. The magnitude of $r$ describes the strength of the
correlation, where a magnitude of $1$ corresponds to a strong linear
association, and a small magnitude corresponds to a weak correlation. An $r$
value of approximately zero indicates that there is no linear correlation
between the two investigated variables. Note that the Pearson Product Moment
Correlation Coefficient picks up linear associations only, thus there may
exist other, non-linear correlations that are not captured by the correlation
coefficient $r$. Therefore it is important to, not only quantitatively compute
input-output correlations, but to also qualitatively assess the relationships
between inputs and outputs via data visualisation in scatterplots as
previously described.

The correlation coefficient, $r^{i}$, describing the correlation between an
input parameter $p^{i}$, and an output response $X$ (in median form) is given
by [21],

$r^{i}=\frac{{{{\sum}}}_{j=1}^{N}\bigg{(}p^{i}_{j}-\bar{p}^{i}\bigg{)}\bigg{(}\underaccent{\tilde}{X}(C_{j})-\underaccent{\tilde}{\bar{X}}\bigg{)}}{\sqrt{\bigg{(}{{{\sum}}}_{j=1}^{N}\big{(}p^{i}_{j}-\bar{p}^{i}\big{)}^{2}\bigg{)}\bigg{(}{{{\sum}}}_{j=1}^{N}\big{(}\underaccent{\tilde}{X}(C_{j})-\underaccent{\tilde}{\bar{X}}\big{)}^{2}\bigg{)}}},$
(18)

where a bar denotes the mean value.

When it comes to interpreting quantitative input-output relationships based on
the correlation coefficient $r$, there are no all-encompassing threshold
values to use for descriptors such as ‘weak’, ‘moderate’, ‘strong’ [21, 27,
11]. Relationships quantified by correlation coefficient values close to the
extrema 0 or 1 may be easy to describe as ‘negligible’ or ‘strong’,
respectively. However, correlation coefficient values in the middle of the
[0,1] range are more difficult to label. Various ‘rules of thumb’ have been
suggested in the literature but, at the end of the day, it is up to the
modeller to appropriately judge what constitutes a ‘weak’, ‘moderate’ or
‘strong’ input-output relationship in the specific (modelling) application at
hand, taking into account the research area, the number of data samples, and
the range of investigated input values [27]. However, even without rigid
descriptor threshold values, we can compare the correlation coefficient values
for all input-output pairs and see which input values are the most influential
within the ranges of regarded input values. As a guide, suggested correlation
coefficient descriptor threshold values presented in the literature are listed
in Table 1. The methodology to perform Latin Hypercube Sampling and Analysis
is outlined in Section 5.3.

reference descriptor |  negligible |  weak |  moderate |  strong |  very strong
---|---|---|---|---|---
Mukaka [21] |  [0,0.3) |  [0.3,0.5) |  [0.5,0.7) |  [0.7,0.9) |  [0.9,1]
Schober et al. [27] |  [0,0.1) |  [0.1,0.4) |  [0.4,0.7) |  [0.7,0.9) |  [0.9,1]
Krehbiel [11] | “A linear relationship exists if $|r|\geq 2/\sqrt{\text{number of samples}}$.”
Table 1: Suggested descriptor threshold values for the magnitude of the
correlation coefficient, $|r|$, reported in the literature.

###  5.3 Quick Guide: Latin Hypercube Sampling and Analysis

##  6 A worked example: Analysing a mathematical cancer model

In this section we will perform Consistency Analysis, Robustness Analysis and
Latin Hypercube Analysis on an agent-based mathematical model that describes a
population of cancer cells (in vitro) that are subjected to an anti-cancer
drug (AZD6738) that may inhibit DNA damage repair in cells and, by extension,
cause cell death. Full details of this model are available in one of our
recent research papers [9], but a pictorial model summary is provided in
Figure 4. This summary contains sufficient information for our current
purposes: performing uncertainty and sensitivity analysis through a worked
example. In order to do this, we need to specify a set of model inputs and
outputs. The full model includes seven input parameters $p^{i}$, $i=1,..,7$,
three of which are mentioned in the model summary and will be investigated in
this review. These input parameters are $p^{1}$: the probability $\Pi_{D-S}$
that a cell enters the damaged S phase in the cell cycle, $p^{2}$: the drug’s
EC50 value and $p^{3}$: the Hill-exponent ($\gamma$) used to compute cellular
drug responses. Furthermore, we consider two in silico measurements as
outputs, specifically $X^{1}$: the percentage of DNA-damaged (i.e.
$\gamma$-H2AX positive) cells at the end of the simulation and $X^{2}$: the
cell count at the end of the simulation.

Figure 4: A summary of an agent-based model that is used to simulate human
colon carcinoma cells in vitro subjected to a drug that targets cellular DNA
damage responses. In this review, we perform uncertainty and sensitivity
analyses on this model as a worked example. Full details are available in the
research paper in which this model was first introduced [9].

###  6.1 Worked example: Consistency Analysis

In order to perform Consistency Analysis, we follow steps 1, 2 and 3 outlined
in the quick guide in Section 3.1.

Step 1: Using the calibrated model parameters, we run our in silico experiment
$20\times(1+5+50+100+300)=9120$ times in order to produce 9120 data samples.
Note that, in this case, we have two output responses of interest, and thus
one data sample consists of an output-pair ($X^{1},X^{2}$). Post the in silico
production of data, we organise our data samples into five groups of
distributions, where each group consists of 20 distributions of data samples.
In the first group, each one of the 20 distributions includes only one data
sample. In the second, third, forth and fifth group, each distribution
respectively includes 5, 50, 100 or 300 data samples.

Step 2: In each distribution group, we compute and plot the $\hat{A}$-measure
(in original and scaled form) for each of its 20 distributions, as is done in
Figures 5 through to 9. These figures clearly demonstrate that the statistical
significance decreases with increasing distribution size $n$. Note that the
two output responses of interest are computed and plotted independently of
each other, as we are aiming to find a distribution size that yields a small
statistical significance for both $X^{1}$ and $X^{2}$.

Step 3: The largest scaled $\hat{A}$-measure in each distribution group is
computed and plotted over the group’s distribution size. The smallest
distribution size for which the statistical significance is small (i.e. $\leq$
0.56) for both $X^{1}$ and $X^{2}$ is denoted $n^{*}$. By regarding Table 2
and Figure 10, we can see that, in this case, $n^{*}=100$. Accordingly, we
determine that $100$ simulation runs are sufficient to mitigate uncertainty
originating from intrinsic model stochasticity. Thus when talking about, for
example, average values and standard deviations produced by this model, we
should base these measures on data samples from $100$ in silico runs.

output distribution size |  n=1 |  n=5 |  n=50 |  n=100 |  n=300
---|---|---|---|---|---
$X^{1}$ |  1 |  0.92 |  0.61 |  0.55 |  0.54
$X^{2}$ |  1 |  0.84 |  0.59 |  0.56 |  0.54
Table 2: Maximal scaled $\hat{A}$-values for various distribution sizes $n$.
The output responses are $X^{1}$: the percentage of $\gamma$H2AX-positive
(i.e. DNA-damaged) cells, and $X^{2}$: the cell count at the end of the
simulation. Figure 5: Consistency Analysis, $\hat{A}$-values in initial (top)
and scaled (bottom) form for distribution size $n=1$.  Figure 6: Consistency
Analysis, $\hat{A}$-values in initial (top) and scaled (bottom) form for
distribution size $n=5$.  Figure 7: Consistency Analysis, $\hat{A}$-values in
initial (top) and scaled (bottom) form for distribution size $n=50$.  Figure
8: Consistency Analysis, $\hat{A}$-values in initial (top) and scaled (bottom)
form for distribution size $n=100$.  Figure 9: Consistency Analysis,
$\hat{A}$-values in initial (top) and scaled (bottom) form for distribution
size $n=300$.  Figure 10: Consistency Analysis, maximal scaled
$\hat{A}$-values for various distribution sizes $n$.

###  6.2 Worked example: Robustness Analysis

We now set out to perform Robustness Analysis using steps 1, 2 and 3 described
in the quick guide in Section 4.1. Here, our goal is to investigate how
sensitive model output responses are to local (one at a time) parameter
perturbations of the inputs $\Pi_{D-S}$, $EC_{50}$ and $\gamma$.

Step 1: The first thing we must do here is to decide appropriate parameter
ranges to investigate for each input parameter. This decision should ideally
be guided by information from experimental data or biological/mathematical
knowledge about the modelling scenario at hand. The mathematical model
summarised in Figure 4 is driven by in vitro data, and from this data suitable
parameter ranges can be deduced and established (full details are available in
the original model paper [9]).

For example, recall that the parameter $EC_{50}$ corresponds to the drug
concentration $C$ that achieves half of $E_{max}$, the maximal drug effect.
Now the in vitro data shows that $C=1\mu$M yields $\approx 0.5E_{max}$ and
that $C=3\mu$M (or higher) yields $\approx E_{max}$, whilst a drug
concentration of $C=0.3\mu$M evokes roughly the same drug response as the
control case (when no drug at all is applied and thus $C=0$) [9, 5]. From
this, we can reason that it is appropriate to investigate $EC_{50}$ parameter
values in the range $(1\pm 0.75)\mu$M. Similarly, we may decide to investigate
$\Pi_{D-S}$ values in the parameter range $(75\pm 10)$%, and $\gamma$ values
in $(2\pm 1)$.

Now that we have established our parameter ranges of interest, it remains to
decide how many (here evenly spaced) parameter values in these ranges to test.
The more parameter values we test, the more detailed information we get, but
recall that for every investigated parameter we need to create $n^{*}$ data
samples (where $n^{*}$ is determined in the Consistency Analysis). The
decision of how many parameter values to include in the Robustness Analysis
should be informed by the fineness of available in vitro/in vivo data,
computational costs, scientific reasoning and model application. If a model is
being used in a pharmaceutical setting, for example, a detailed analysis may
be of extra importance. Here, we respectively choose to investigate 9,7 and 9
evenly spaced parameter values (including the calibrated parameter values) in
the parameter ranges for $\Pi_{D-S}$, $EC_{50}$ and $\gamma$. Thus we need to
produce $n^{*}\times(9+7+9)=100\times(25)$ in silico data samples.

Steps 2 and 3: Post in silico simulations, the $\hat{A}$-measures are computed
and plotted over all parameter values $p^{i}_{j}$ for each input parameter
$p^{i}$. Such plots are here shown for input parameters $\Pi_{D-S}$ (Figure
11), $EC_{50}$ (Figure 12) and $\gamma$ (Figure 13).

The data samples produced in the in silico experiments are also represented
using boxplots in these figures. The $\hat{A}$-measures and the boxplots
demonstrate the effect that local perturbations of the input parameters have
on both output responses $X^{1}$ and $X^{2}$.

Figure 11 illustrates that increasing the probability $\Pi_{D-S}$ that a cell
enters the damaged S state increases the percentage of $\gamma$H2AX-positive
cells and decreases the cell count, as is to be expected. Further, Figure 12
demonstrates that the model output is highly sensitive to perturbations of
$EC_{50}$. Increasing $EC_{50}$ results in a higher percentage of
$\gamma$H2AX-positive cells and a lower cell count. Finally, Figure 13 shows
that the regarded output responses are less sensitive to small perturbations
of the Hill-exponent $\gamma$ than to small perturbations of $\Pi_{D-S}$ and
$EC_{50}$.

Figure 11: Robustness Analysis, Left: $\hat{A}$-values resulting from
comparisons between distributions of data samples produced with perturbed
$\Pi_{D-S}$ values, and the distribution produced with the calibrated
(unperturbed) $\Pi_{D-S}$ value. Right: Output responses, in terms of
percentage of $\gamma$H2AX-positive (i.e. damaged) cells, and cell count as a
result of perturbations to the input variable $\Pi_{D-S}$.

Figure 12: Robustness Analysis, Left: $\hat{A}$-values resulting from
comparisons between distributions of data samples produced with perturbed
$EC_{50}$ values, and the distribution produced with the calibrated
(unperturbed) $EC_{50}$ value. Right: Output responses, in terms of percentage
of $\gamma$H2AX-positive (i.e. damaged) cells, and cell count as a result of
perturbations to the input variable $EC_{50}$.

Figure 13: Robustness Analysis, Left: $\hat{A}$-values resulting from
comparisons between distributions of data samples produced with perturbed
$\gamma$ values, and the distribution produced with the calibrated
(unperturbed) $\gamma$ value. Right: Output responses, in terms of percentage
of $\gamma$H2AX-positive (i.e. damaged) cells, and cell count as a result of
perturbations to the input variable $\gamma$.

###  6.3 Worked example: Latin Hypercube Analysis

Following steps 1, 2 and 3, as described in the quick guide in Section 5.3,
Latin Hypercube Analysis is here performed in order to investigate how
sensitive output responses are to global parameter perturbations. We here
investigate parameter values within parameter ranges that we consider to be
‘plausible’ post Robustness Analysis.

Step 1: For each input parameter, we decide to split the investigated
parameter range into $N=100$ intervals. Why 100 you might ask? Well, in the
original paper where we first introduced this model [9], the model took seven
input parameters and thus $q=7$. Accordingly, we tried using $N=4q/3\approx
10$ and $N=2q=14$ intervals at first (following the suggestions discussed in
Section 5) but neither of these options produced enough data samples to yield
meaningful information in steps 2 and 3 below. Therefore, we decided to use
$N=100$ instead, as this choice covered a larger range of the input parameter
space whilst coming at a feasible computational cost. Now, we can use the
built-in MATLAB function lhsdesign [17] to create combinations of input
parameter values (represented by a point $C_{p}$ in input parameter space)
that shall be used to produce the in silico data samples needed for Latin
Hypercube Analysis.

Step 2: For each $C_{p}$, median output responses ($X^{1}$ and $X^{2}$) of
$n^{*}=100$ in silico runs are computed. These are plotted in ’output-over-
input’ scatterplots for each investigated input parameter in Figures 14, 15
and 16. From these figures we can make some qualitative remarks: Figure 14
indicates that the relationships between the input variable $\Pi_{D-S}$ and
the output responses $X^{1}$ and $X^{2}$ are, respectively, positively and
negatively correlated. This agrees with the intuitive notion that if the
probability that a cell enters the D-S state increases, so does the percentage
of damaged cells ($X^{1}$) whilst the cell count ($X^{2}$) decreases as more
cells will be susceptible to the drug and potentially die. The scatterplots in
Figure 15 demonstrate that the input variable $EC_{50}$ impacts the output
responses more than do other investigated input parameters (within the
regarded ranges). $EC_{50}$ is negatively, linearly correlated with $X^{1}$
and positively, linearly correlated with $X^{2}$. In Figure 16, however, there
is no visually apparent correlation between the input parameter $\gamma$ and
the output.

Step 3: After making some qualitative remarks in Step 2, we now compute the
Pearson Product Moment Correlation Coefficients between the various input-
output pairs for a quantitative analyses. These correlation coefficients are
listed in Table 3. To decide threshold values for correlation coefficient
descriptors, we here compromise between threshold values suggested by other
authors (listed in table 1) whilst taking into account that we are only
regarding parameter values within ‘plausible’ ranges. With this as a guide, we
here decide to refer to the linear input-output relationship as being
‘negligible’ for $\gamma$, where the obtained correlation coefficients are
0.05 and 0.12 for $X^{1}$ and $X^{2}$ respectively. We further say that for
$\Pi_{D-S}$ the linear input-output relationship is positively/negatively weak
for $X^{1}$ and $X^{2}$ respectively. For $EC_{50}$, however, the linear
input-output relationship is moderately negative for $X^{1}$ and strongly
positive for $X^{2}$. Clearly, we must be careful when choosing our $EC_{50}$
value in the model, as this highly influences the model output!

output input |  $\Pi_{D-S}$ |  $EC_{50}$ |  $\gamma$
---|---|---|---
$X^{1}$ |  0.19 |  -0.59 |  0.05
$X^{2}$ |  -0.24 |  0.84 |  0.12
Table 3: The Pearson Product Moment Correlation Coefficients, computed in the
Latin Hypercube Analysis, quantitatively describe input-output relationships.
Figure 14: Latin Hypercube Analysis. Output responses in terms of
$\gamma$H2AX-positive cells (left) and cell count (right) when global
parameter perturbations are performed. The scatterplots show the correlation
between outputs and the input value of $\Pi_{D-S}$.  Figure 15: Latin
Hypercube Analysis. Output responses in terms of $\gamma$H2AX-positive cells
(left) and cell count (right) when global parameter perturbations are
performed. The scatterplots show the correlation between outputs and the input
value of $EC_{50}$.  Figure 16: Latin Hypercube Analysis. Output responses in
terms of $\gamma$H2AX-positive cells (left) and cell count (right) when global
parameter perturbations are performed. The scatterplots show the correlation
between outputs and the input value of $\gamma$.

##  7 Conclusion

This review is intended as a gentle, introductory review to three uncertainty
and sensitivity analyses methods, namely, Consistency Analysis, Robustness
Analysis and Latin Hypercube Analysis. Information on how to implement these
methods in MATLAB are available in the Appendix. Alternatively, all methods
discussed in this review can be implemented using the R-based software package
Spartan, developed by Alden et al. [1]. In fact, many of the proceedings and
conventions used in this review follow those suggested by Alden et al. in
order to allow the reader to, as easily as possible, use Spartan if desired.
Scrutinising mathematical models using uncertainty and sensitivity analyses
methods is an important part in model development. In many applications,
knowledge about a model’s robustness is crucial [30]. In the context of
quantitative pharmacology, for example, a mathematical model may be used to
guide preclinical or, ultimately, clinical proceedings. In such cases,
understanding how confident we can be with model results, and how sensitive a
model is to parameter perturbations, is of the utmost importance.

## Acknowledgements

SH was supported by the Medical Research Council [grant code MR/R017506/1] and
Swansea University PhD Research Studentship. SS was supported by an STFC
studentship under the DTP grant ST/N504464/1.

## Appendix – MATLAB code snippets

### Computing measure of stochastic superiority

We here list two different MATLAB functions that can be used in order to
compute the point estimate of the A-measure of stochastic superiority in the
original form, $\hat{A}\in[0,1]$, and in the scaled form,
$\hat{\underline{A}}\in[0.5,1]$. The function getA_measure_naive, listed
below, uses direct implementations of Equations 9 and 11 to compute and return
values for $\hat{A}_{x_{0},x_{1}}$ and $\underline{\hat{A}}_{x_{0},x_{1}}$,
given two input vectors $x_{0}$ and $x_{1}$. The function getA_measure uses
the built-in MATLAB function ranksum to do the same.

⬇

function [A_measure, scaled_A_measure] = getA_measure(x0, x1)

[p,h,stats] = ranksum(x0,x1);

% Compute the A measure

A_measure=(stats.ranksum/length(x0) - (length(x0)+1)/2)/length(x1);

% Compute the scaled A measure

scaled_A_measure=0.5+abs(0.5 -A_measure);

end

⬇

function [A_measure, scaled_A_measure] = getA_measure_naive(x0, x1)

% Compute the A measure

A_measure = 0;

for i = 1:length(x0)

for j = 1:length(x1)

if(x0(i)>x1(j))

A_measure = A_measure + 1;

elseif(x0(i)==x1(j))

A_measure = A_measure + 0.5;

elseif(x0(i)<x1(j))

A_measure = A_measure + 0;

end

end

end

A_measure = A_measure/(length(x0)*length(x1));

% Compute the scaled A measure

if(A_measure>=0.5)

scaled_A_measure = A_measure;

else

scaled_A_measure = 1-A_measure;

end

end

### Creating boxplots

The MATLAB function boxplot can be used to create boxplots. The input data in
one column is represented by one box in the boxplot. For details regarding
labeling and style alternatives, please see the MATLAB documentation [17].

⬇

boxplot(M);

### Choosing Latin Hypercube Sampling points

A Latin Hypercube Sampling matrix can be created using the MATLAB function
lhsdesign, which returns a matrix of size $n\times q$, where $n$ denotes the
number of samples to be tested, and $q$ denotes the number of input parameters
to investigate (and thus perturb).

⬇

LHC_Matrix=lhsdesign(n,q)

Each row $i$, in the created matrix (here denoted LHC_Matrix), corresponds to
the $i$th sampling point. Each element ($i,j$) corresponds to the parameter
value of the $j$th input parameter in sampling point $i$, where each parameter
ranges between 0 and 1. For different criteria on how to chose the specific
parameter values within each sampled interval, please refer to the MATLAB
documentation [17]. Sampling points can, for example, be chosen in a way that
maximises the distance between sampling points in the $q$-dimensional sampling
space.

### Qualitative and Quantitative Latin Hypercube Sampling Analysis

In order to qualitatively asses the correlation between an input parameter
$p$, and an output response $X$, one can use the MATLAB function scatter. In
the below listings, p and X are two data vectors.

⬇

scatter(p,X)

Further, to quantify the linear correlation between $p$ and $X$, the MATLAB
function corrcoef can be used to compute correlation coefficients.

⬇

R=corrcoef(p, X);

## References

  * [1] K. Alden, M. Read, J. Timmis, P. S. Andrews, H. Veiga-Fernandes, and M. Coles.  Spartan: a comprehensive tool for understanding uncertainty in simulations of biological systems.  PLoS Comput. Biol., 9(2):e1002916, 2013.
  * [2] S.M. Blower and Hadi Dowlatabadi.  Sensitivity and uncertainty analysis of complex models of disease transmission: An hiv model, as an example.  International Statistical Review, 62, 08 1994.
  * [3] E. O. Buzbas and N. A. Rosenberg.  AABC: approximate approximate Bayesian computation for inference in population-genetic models.  Theor Popul Biol, 99:31–42, Feb 2015.
  * [4] A. Charzyńska, A. Nałęcz, M. Rybiński, and A. Gambin.  Sensitivity analysis of mathematical models of signaling pathways.  BioTechnologia., 93 (3):291–308, 2012.
  * [5] S. Checkley, L. MacCallum, J. Yates, P. Jasper, H. Luo, J. Tolsma, and C. Bendtsen.  Bridging the gap between in vitro and in vivo: Dose and schedule predictions for the ATR inhibitor AZD6738.  Sci Rep, 5:13545, Aug 2015.
  * [6] J. Cohen.  The statistical power of abnormal-social psychological research: a review.  J Abnorm Soc Psychol, 65:145–153, Sep 1962.
  * [7] O. Cohen.  Statistical Power Analysis for the Behavioral Sciences (Second Edition) .  Lawrence Erlbaum Associates, 1988.
  * [8] J. L. Gevertz and J. R. Wares.  Developing a Minimally Structured Mathematical Model of Cancer Treatment with Oncolytic Viruses and Dendritic Cell Injections.  Comput Math Methods Med, 2018:8760371, 2018.
  * [9] S. Hamis, J. Yates, M.A.J Chaplain, and G.G. Powathil.  Bridging in vitro and in vivo research via an agent-based modelling approach: predicting tumour responses to an atr-inhibiting drug.  Preprint: bioRxiv, doi: 10.1101/841270.
  * [10] R.L. Iman and J.C. Helton.  Comparison of uncertainty and sensitivity analysis techniques for computer models.  Report NUREGICR-3904, SAND 84-1461, Sandia National Laboratories, Albuquerque, New Mexico, 3 1985.
  * [11] Timothy Krehbiel.  Correlation coefficient rule of thumb.  Decision Sciences Journal of Innovative Education, 2:97–100, 01 2004.
  * [12] B. Lambert, A. L. MacLean, A. G. Fletcher, A. N. Combes, M. H. Little, and H. M. Byrne.  Bayesian inference of agent-based models: a tool for studying kidney branching morphogenesis.  J Math Biol, 76(7):1673–1697, 06 2018.
  * [13] J. Liepe, P. Kirk, S. Filippi, T. Toni, C. P. Barnes, and M. P. Stumpf.  A framework for parameter estimation and model selection from experimental data in systems biology using approximate Bayesian computation.  Nat Protoc, 9(2):439–456, Feb 2014.
  * [14] A. Ligmann-Zielinska, D. B. Kramer, K. Spence Cheruvelil, and P. A. Soranno.  Using uncertainty and sensitivity analyses in socioecological agent-based models to improve their analytical performance and policy relevance.  PLoS ONE, 9(10):e109779, 2014.
  * [15] Shenglin Lin, Wei Li, Xiaochao Qian, Ping Ma, and Ming Yang.  A Simulation Model Validation and Calibration Platform.  pages 687–693, 12 2018.
  * [16] G. Manache and C. Melching.  Sensitivity of Latin Hypercube Sampling to sample size and distributional assumptions.  07 2007.
  * [17] MATLAB.  version 1.8.0_202 (R2019n).  The MathWorks Inc., Natick, Massachusetts, 2019.
  * [18] K. O. McGraw and S. P. Wong.  A common language effect size statistic.  Psychological Bulletin, 111(2):361–365, 1992.
  * [19] M. D. McKay, R. J. Beckman, and W. J. Conover.  Comparison of three methods for selecting values of input variables in the analysis of output from a computer code.  Technometrics, 21(2):239–245, 1979.
  * [20] Michael D. McKay.  Latin hypercube sampling as a tool in uncertainty analysis of computer models.  In Proceedings of the 24th Conference on Winter Simulation, WSC ’92, pages 557–564, New York, NY, USA, 1992. ACM.
  * [21] M. M. Mukaka.  Statistics corner: A guide to appropriate use of correlation coefficient in medical research.  Malawi Med J, 24(3):69–71, Sep 2012.
  * [22] A. Niida, T. Hasegawa, and S. Miyano.  Sensitivity analysis of agent-based simulation utilizing massively parallel computation and interactive data visualization.  PLoS ONE, 14(3):e0210678, 2019.
  * [23] M. Read, P. S. Andrews, J. Timmis, and V. Kumar.  Techniques for grounding agent-based simulations in the real domain: a case study in experimental autoimmune encephalomyelitis.  Mathematical and Computer Modelling of Dynamical Systems, 18(1):67–86, 2012.
  * [24] K. A. Rejniak and A. R. Anderson.  Hybrid models of tumor growth.  Wiley Interdiscip Rev Syst Biol Med, 3:115–125, 2011.
  * [25] J. Ruscio and T. Mullen.  Confidence Intervals for the Probability of Superiority Effect Size Measure and the Area Under a Receiver Operating Characteristic Curve.  Multivariate Behavioral Research, 47(2):201–223, 2012.
  * [26] Andrea Saltelli and Ricardo Bolado.  An alternative way to compute fourier amplitude sensitivity test (fast).  Comput. Stat. Data Anal., 26(4):445–460, February 1998.
  * [27] P. Schober, C. Boer, and L. A. Schwarte.  Correlation Coefficients: Appropriate Use and Interpretation.  Anesth. Analg., 126(5):1763–1768, 05 2018.
  * [28] Razi Sheikholeslami and Saman Razavi.  Progressive latin hypercube sampling: An efficient approach for robust sampling-based analysis of environmental models.  Environmental Modelling and Software, 93:109–126, 07 2017.
  * [29] A. Vargha and H. D. Delaney.  A Critique and Improvement of the CL Common Language Effect Size Statistics of McGraw and Wong.  Journal of Educational and Behavioral Statistics, 25(2):101–132, 2000.
  * [30] S. A. Visser, D. P. de Alwis, T. Kerbusch, J. A. Stone, and S. R. Allerheiligen.  Implementation of quantitative and systems pharmacology in large pharma.  CPT Pharmacometrics Syst Pharmacol, 3:e142, Oct 2014.
  * [31] X. Y. Zhang, M. N. Trame, L. J. Lesko, and S. Schmidt.  Sobol Sensitivity Analysis: A Tool to Guide the Development and Evaluation of Systems Pharmacology Models.  CPT Pharmacometrics Syst Pharmacol, 4(2):69–79, 02 2015.