# Efficient Analysis of Chemical Reaction Networks Dynamics based on Input-
Output Monotonicity

Lucia Nasti Department of Computer Science - University of Pisa
Largo Bruno Pontecorvo, 3,
56127, Pisa, Italy  Roberta Gori Department of Computer Science - University
of Pisa
Largo Bruno Pontecorvo, 3,
56127, Pisa, Italy  Paolo Milazzo Department of Computer Science - University
of Pisa
Largo Bruno Pontecorvo, 3,
56127, Pisa, Italy  Federico Poloni Department of Computer Science -
University of Pisa
Largo Bruno Pontecorvo, 3,
56127, Pisa, Italy

###### Abstract

Motivation: A Chemical Reaction Network (CRN) is a set of chemical reactions,
which can be very complex and difficult to analyze. Indeed, dynamical
properties of CRNs can be described by a set of non-linear differential
equations that rarely can be solved in closed-form, but that can instead be
used to reason on the system dynamics. In this context, one of the possible
approaches is to perform numerical simulations, which may require a high
computational effort. In particular, in order to investigate some dynamical
properties, such as robustness or global sensitivity, many simulations have to
be performed by varying the initial concentration of chemical species.
Results: In order to reduce the computational effort required when many
simulations are needed to assess a property, we exploit a new notion of
monotonicity of the output of the system (the concentration of a target
chemical species at the steady state) with respect to the input (the initial
concentration of another chemical species). To assess such monotonicity
behaviour, we propose a new graphical approach that allows us to state
sufficient conditions for ensuring that the monotonicity property holds. Our
sufficient conditions allow us to efficiently verify the monotonicity property
by exploring a graph constructed on the basis of the reactions involved in the
network. Once established, our monotonicity property allows us to drastically
reduce the number of simulations required to assess some dynamical properties
of the CRN.

Chemical Reaction Networks — Monotonicity — Simulations — ODEs

##  1 Introduction

From the discovery of DNA structure, in 1953, there has been a growing
interest in understanding the morphological and functional organization of
living cells (Kitano, 2002). Cells are very complex to analyze because they
consist of many components that interact with each other, through multiple
sequences of chemical reactions, Chemical Raction Networks (CRNs), which
regulate the overall behavior. Besides, several fluctuations can alter the
cell functionalities, such as internal errors propagation and variation in the
concentrations of chemical species.

Bioinformatics and systems biology emerge as powerful tools to investigate CRN
dynamics merging computational methods and real data. Through simulations, for
instance, it is possible to mimic the internal dynamics of a natural system
and, therefore, to predict its behaviour. Moreover, model-based analysis
techniques can be used to interpret some less intuitive features of the
system.

In this context, computer scientists developed many formalisms to study
systems of interacting components, which can be applied to model and describe
CRNs and, in general, biological systems. Among these formalisms, those that
have been applied in systems biology (Bernini et al., 2018) include Petri nets
(Behinaein et al., 2014; Murata, 1989; Koch, 2010), Hybrid systems (Alur et
al., 1992; Henzinger, 2000; Li et al., 2017), process calculi such as the
$\pi$-calculus (Regev et al., 2000) as well as many ad-hoc biologically
inspired calculi (Danos et al., 2008) and rule-based systems (Barbuti et al.,
2011).

The development of models that can help in predicting the system behaviour
requires precise and detailed information about the set of initial conditions
of the biological system under study. In many cases, obtaining such precise
informations is unfortunately very challenging (or even impossible) because of
the noisy nature of biological data (Dresch et al., 2010). Moreover, some
parameters may be affected by fluctuations that alter ordinary system
behavior. Finally, in many cases such precise informations simply cannot be
measured.

In order to determine which model parameters are more critical in case of
perturbations and approximate measurements, it is possible to apply
sensitivity analysis methods, which give a measure of the behavioral change of
the system under perturbation of one of its parameters. As underlined in
(Iooss & Lemaître, 2015; Zi, 2011), there exist two fundamental approaches to
sensitivity analysis, local and global. While the local sensitivity analysis
investigates the effects of small perturbations, the global one studies the
effects of large perturbations and determines also the most (or the least)
influencing parameters.

Sensitivity analysis methods (in particular, the global ones) typically
require performing many simulations, by varying the system parameters one by
one. These methods are, in general, quite expensive, because of the number of
parameters and the large range of values to be tested.

For these reasons, alternative dynamical properties of CRNs, such as
monotonicity (Angeli et al., 2006; Gori et al., 2019) and steady-state
reachability (Feinberg, 1987), have been studied. Establishing such
properties, indeed, provides information on the CRN dynamics without the need
of performing several numerical simulations (Nasti, 2020).

Monotonicity, in particular, is a property stating that a given measurable
aspect of the system dynamics increases (or decreases) with the increase of a
given system parameter. Many more specific definitions of monotonicity exist
(Angeli et al., 2006; Gori et al., 2019) and in some cases they can be tested
simply by inspecting the structure of the system models. It is worth noting
that the validation of many different biological properties, such as
robustness (Kitano, 2002; Rizk et al., 2011; Shinar & Feinberg, 2010),
persistence (Angeli et al., 2007), and adaptation (Shinar et al., 2009),
greatly benefits from the assessment of monotonicity properties of the
network, since this typically allows reducing the number of cases to be
analyzed through numerical simulations.

Consider, for example, the robustness property. It is observed in many
biological systems and it expresses the ability of the system to preserve its
functions despite the presence of perturbations (Kitano, 2007). Without any
assumption on monotonicity, verifying robustness would require, in general, to
consider all possible perturbations, usually expressed as different initial
states of the system. In particular, regarding the CRNs, it would be necessary
to test the system behaviour by examining all the possible combinations of
initial concentrations of chemical species and, in practice, this would
require a huge (in principle, infinite) number of simulations (Nasti et al.,
2018).

The same reasoning applies also to the case in which the initial
concentrations of a biological network are simply unknown. Without any
assumption on monotonicity, understanding the qualitative system behaviour
would require, in principle, to consider all possible initial concentrations.

In this paper we propose a sufficient condition for CRNs that, if satisfied,
ensures that a form of monotonicity, called Input-Output monotonicity, holds.
Given a species considered as the _input_ of a CRN, and another species
considered as the _output_ , we say that input and output are in a
monotonicity relation if the concentration of the output species always
increases (or decreases) in response to an increase in the initial
concentration of the input. If this monotonicity relationship holds, than it
is possible to consider an interval of initial concentrations for the input
species and obtain the corresponding concentrations interval for the output
species by performing only two simulations, one for each extreme value of the
input initial concentration (Gori et al., 2019).

The sufficient condition we propose is based on a condition on the structure
of the CRN that can be efficiently evaluated, without the need of performing
any simulation. Following the lines of (Angeli et al., 2006), our condition is
based on a graph representation of the CRN enriched with information about
cooperation and competition among reactions and it is expressed as a set of
constraints on the graph structure.

The paper is organized as follows. In the rest of this section we recall some
background definitions of CRNs and related concepts. In Section 2 we define
Input-Output monotonicity and state Theorem 5 expressing our sufficient
condition. In Section 3, we apply our methodology to study the case of the ERK
signaling pathway. In Section 5 we draw our conclusions and discuss future
work. Finally, in Section 4, we give the proof of Theorem 5.

###  1.1 Chemical Reaction Networks

A Chemical Reaction Network (CRN) is a set of reactions that can be formally
defined as follows.

######  Definition 1 (Chemical reaction network).

Given an indexed set of chemical species
$\mathcal{S}=(\mathcal{S}_{1},\mathcal{S}_{2},\dots,\mathcal{S}_{s})$, a
_Chemical Reaction Network_ is an indexed set of reactions
$\mathcal{R}=(\mathcal{R}_{1},\mathcal{R}_{2},\dots,\mathcal{R}_{r})$
involving such species. Each chemical reaction is denoted as

$\mathcal{R}_{i}:\sum_{j=1}^{s}\alpha_{ij}\mathcal{S}_{j}\rightarrow\sum_{j=1}^{s}\beta_{ij}\mathcal{S}_{j}\qquad\mbox{
or
}\qquad\mathcal{R}_{i}:\sum_{j=1}^{s}\alpha_{ij}\mathcal{S}_{j}\rightleftharpoons\sum_{j=1}^{s}\beta_{ij}\mathcal{S}_{j},$

where $\alpha_{ij}$ and $\beta_{ij}$ are non-negative integers called the
stoichiometric coefficients.

The arrow is used to indicate the direction in which a chemical reaction takes
place (from reactants to products). When we have a single arrow
($\rightarrow$), the reaction is irreversible, namely that a reaction
transforming products back into reactants cannot take place. When there is a
double arrow ($\rightleftharpoons$) it means that it is possible to have both
the forward and the backward transformation, and then the reaction is
reversible. In a CRN, a reversible reaction can be equivalently represented as
a pair of irreversible reactions with opposite directions.

The _stoichiometric coefficient matrix_ $\Gamma\in\mathbb{R}^{s\times r}$ is
defined as

$\Gamma_{ji}=\beta_{ij}-\alpha_{ij}.$

Note that we allow in our networks the presence of _promoters_ : a promoter is
a species $\mathcal{S}_{j}$ (such as an enzyme) that affects the rate of a
reaction $\mathcal{R}_{i}$ but does not get produced nor consumed by it, i.e.,
$\Gamma_{ji}=0$.

The _rate_ of a reaction is expressed as a function of the concentrations of
the species in the network. The vector of species concentrations at a given
time is denoted by $S=(S_{1},S_{2},...,S_{s})^{\prime}\in\mathbb{R}^{s}$. The
vector of reaction rates (which are functions of $S$) is denoted by
$R(S)=(R_{1}(S),R_{2}(S),...,R_{r}(S))^{\prime}\in\mathbb{R}^{r}$.

Following (Angeli et al., 2010), we assume that all entries of the Jacobian of
$R(S)$ (which we denote by $DR\in\mathbb{R}^{r\times s}$) have a well-defined
constant sign that does not depend on $S$ (although they may become zero for
certain values of $S$), and that

$DR_{ij}\Gamma_{ji}\leq 0\quad\text{for all $i,j$}.$ (1)

The assumption (1) is a natural one because in most cases the rate of a
reaction increases with the quantity of reactants. This assumption covers, in
particular, the most common case, described by the well-established law of
_mass action_ , in which the rate of a reaction is proportional to the
concentration of each of its reactants.

Following a deterministic approach, the evolution of the concentrations in
time is usually described by a system of differential equations:

$\frac{dS}{dt}=\Gamma\cdot{R(S)},\quad S\in\mathbb{R}_{\geq 0}^{s}.$

######  Example 1 (Enzyme kinetics).

We consider the simple enzymatic reaction network

$\ch{E+S<=>[k_{1}][k_{-1}]ES->[k_{2}]E+P}.$

In this CRN, the set of chemical species is $\mathcal{S}$ = (E, S, ES, P). The
enzyme E, binding the substrate S, forms a complex ES, which releases the
product P and the original enzyme E. According to _mass-action kinetics_ , the
rate of each reaction is directly proportional to the concentrations of its
reactants, via a coefficient marked next to each arrow ($k_{1}$,$k_{-1}$ and
$k_{2}$, respectively). In our framework, this system is represented with a
reversible reaction $\mathcal{R}_{1}$ and an irreversible one
$\mathcal{R}_{2}$, with rate vector

$R(S)=\begin{bmatrix}k_{1}[E][S]-k_{-1}[ES]\\\ k_{2}[ES]\end{bmatrix}.$

Here, we have used the symbol $[E]$ to denote the concentration of E, and so
on.

The differential equations describing the behavior of the CRN are

$\begin{cases}&\frac{d[E]}{dt}=-k_{1}[E][S]+k_{-1}[ES]+k_{2}[ES],\\\
&\frac{d[S]}{dt}=-k_{1}[E][S]+k_{-1}[ES],\\\
&\frac{d[ES]}{dt}=+k_{1}[E][S]-k_{-1}[ES]-k_{2}[ES],\\\
&\frac{d[P]}{dt}=+k_{2}[ES].\\\ \end{cases}$

As we can notice, a CRN is usually described by a set of non-linear
differential equations, which makes difficult the analysis of the dynamics of
the system. Indeed, different initial concentrations of the species involved
in the reactions can affect the internal and external behavior of the network.
Besides, the exact values of parameters are often unknown. Then, in order to
investigate the behaviour of the system, we need to study it by performing
many simulations under all possible combinations of chemical species
concentrations.

To reduce the computational effort, one of the possible approaches is to study
the qualitative behaviour of CRNs, without making assumptions on the structure
of the dynamical equations involved. In this context, establishing some kind
of monotonicity property can help to answer questions concerning the network
asymptotic dynamics, such as which are the functionalities of specific
chemical pathways or how parameter variations influence the network. Indeed,
monotonicity describes the capacity of a system to respond in a natural way to
perturbation on its components.

According to Angeli et al. (Angeli et al., 2006, 2010) a system is monotonic
if its forward flow preserves some order defined on the state space. The
system dynamics is expressed in terms of reaction coordinates, for which all
the reaction processes are taken into account as a unique flux. As a
consequence, with this approach, while it is possible to study how internal or
external perturbations influence the fluxes generated by the reactions, it is
not possible to address (and therefore understand) how different
concentrations of the chemical species involved in the network influence the
overall behaviour. The main advantage of Angeli’s approach is that the authors
provide efficiently verifiable sufficient conditions to establish if a system
is globally monotonic, based on a graphical representations of reactions. In
particular, they investigate the dynamics of the system using a particular
kind of graph, the so called R-graph. Their sufficient condition for global
monotonicity is the following: the network is globally monotonic if each
closed path of the R-graph contains an even number of negative edges (positive
loop property). In the next section, we shall see a formal definition of this
graph and comment on this condition.

##  2 Results

### Input-Output monotonicity definition

Global monotonicity, proposed in (Angeli et al., 2006), is a very strong
property, since it is based on an unique ordering on the whole reaction
network. Unfortunately, such a strong property does not hold on most realistic
chemical reaction networks.

For this reason, our goal is to assess monotonicity properties between species
concentrations that would allow us to infer the behaviour of the network under
different initial conditions and to define sufficient conditions for chemical
reactions networks that are easy to test and guarantee that such monotonicity
properties hold.

The monotonicity property we are interested in can be summarised as follows.
Given two species, that we call input and output species of the network, we
say that there is a monotonicity relation if the concentration of output
species at any time either increases or decreases due to an increase in the
initial concentration of the input species. If established, our monotonicity
property would allow us to substantially reduce the number of simulations
required to study the system. Indeed, if two species are in a monotonicity
relation and we are interested in studying the dynamics of the output when the
input varies, we can avoid to simulate the chemical reaction network for all
possible values of the initial concentrations of the input species. To
formally define the previous intuition, we give a new definition of
monotonicity, namely the Input-Output monotonicity, and then, following the
approach in (Angeli et al., 2010), we give sufficient conditions that
guarantee the monotonicity relation of the input and output species.

The following two definitions describe the concept of monotonicity we are
interested in: it describes whether the output species reacts in a monotonic
way to the increase of the input concentration.

We consider two initial states $S^{0},\overline{S^{0}}$ such that
$\overline{S^{0}_{I}}>S^{0}_{I}$ for one particular species $I$ (the _input
species_), and $S^{0}_{k}=\overline{S^{0}_{k}}$ for all other species $k\neq
I$. With $S_{i}(t)$ we indicate the solution of the ODEs for the species
$S_{i}$ with initial value $S(0)=S^{0}$, and with $\overline{S_{i}}(t)$ the
solution with initial value $\overline{S}(0)=\overline{S^{0}}$.

######  Definition 2 (Positive Input-Output Monotonicity).

Given a set of reactions $\mathcal{R}$, species
$\mathcal{S}_{O}\neq\mathcal{S}_{I}$ is _positively monotonic_ with respect to
$\mathcal{S}_{I}$ in $\mathcal{R}$ if, for any two initial states
$S^{0},\overline{S^{0}}$ as above, $\overline{S_{O}}(t)\geq S_{O}(t)$, for
every time $t\in\mathbb{R}_{\geq 0}$.

######  Definition 3 (Negative Input-Output Monotonicity).

Given a set of reactions $\mathcal{R}$, species
$\mathcal{S}_{O}\neq\mathcal{S}_{I}$ is _negatively monotonic_ with respect to
$\mathcal{S}_{I}$ in $\mathcal{R}$ if, for any two initial states
$S^{0},\overline{S^{0}}$ as above, $\overline{S_{O}}(t)\leq S_{O}(t)$, for
every time $t\in\mathbb{R}_{\geq 0}$.

### Directed SR-graph and R-graph

We define two graphs associated to a reaction network. The directed SR-graph
models the interplay between species and reactions.

######  Definition 4 (Directed SR-graph).

Given a finite set of reactions $\mathcal{R}$ over a set of species
$\mathcal{S}$, the _directed SR-graph_ is the directed graph
$\langle\mathcal{S}\cup\mathcal{R},E\rangle$, where
$E\subseteq(\mathcal{S}\cup\mathcal{R})\times(\mathcal{S}\cup\mathcal{R})$ is
defined as follows.

  * •

If $\Gamma_{jk}\neq 0$, i.e., $\mathcal{R}_{k}$ affects the concentration of
$\mathcal{S}_{j}$, then we include the edge
$(\mathcal{R}_{k},\mathcal{S}_{j})$ in $E$.

  * •

If $DR_{ij}\neq 0$, i.e., $\mathcal{S}_{j}$ affects the speed of
$\mathcal{R}_{i}$, then we include the edge
$(\mathcal{S}_{j},\mathcal{R}_{i})$ in $E$.

No other edges exist in $E$ (in particular, there no edges in
$\mathcal{S}\times\mathcal{S}$ nor in $\mathcal{R}\times\mathcal{R}$.

In particular, for each pair $(\mathcal{R}_{i},\mathcal{S}_{j})$ so that
$\mathcal{S}_{j}$ is involved in $\mathcal{R}_{i}$, one of these three cases
applies.

  * C1

if $\mathcal{S}_{j}$ is a promoter for $\mathcal{R}_{i}$ that does not get
consumed, i.e., $DR_{ij}\neq 0$ but $\Gamma_{ji}=0$, then the edge
$\mathcal{S}_{j}\to\mathcal{R}_{i}$ exists only in this direction;

  * C2

if $\mathcal{S}_{j}$ is produced or consumed by $\mathcal{R}_{i}$ but does not
affect its rate (for instance the product of an irreversible reaction), i.e.,
$DR_{ij}=0$ but $\Gamma_{ji}\neq 0$, then the edge
$\mathcal{R}_{i}\to\mathcal{S}_{j}$ exists only in this direction;

  * C3

in all other cases in which $\mathcal{S}_{j}$ is involved $\mathcal{R}_{i}$,
the edge $\mathcal{R}_{i}\leftrightarrow\mathcal{S}_{j}$ is bidirectional.

######  Definition 5 (Consistent labeling).

A _consistent labeling_ of the set of reactions $\mathcal{R}$ is a map
$\sigma:\mathcal{R}\to\\{+1,-1\\}$ such that

$\sigma(i)DR_{ij}\Gamma_{jk}\sigma(k)\geq 0,\quad\forall i\neq k,\,\forall j.$
(2)

The second graph that we define is a signed but undirected graph that
represents the constraints for the existence of a consistent labeling.

######  Definition 6 (R-graph).

Given a finite set of reactions $\mathcal{R}$ over a set of species
$\mathcal{S}$, the R-graph of $\mathcal{R}$ is the signed graph
$\langle\mathcal{R},E_{+},E_{-}\rangle$, where
$E_{+}\subseteq(\mathcal{R}\times\mathcal{R})$ and
$E_{-}\subseteq(\mathcal{R}\times\mathcal{R})$ are defined as follows:

  * •

$(\mathcal{R}_{i},\mathcal{R}_{k})\in E_{+}$ if $i\neq k$ and there is a
species $\mathcal{S}_{j}$ such that $DR_{ij}\Gamma_{jk}>0$.

  * •

$(\mathcal{R}_{i},\mathcal{R}_{k})\in E_{-}$ if $i\neq k$ and there is a
species $\mathcal{S}_{j}$ such that $DR_{ij}\Gamma_{jk}<0$.

Note that an edge in $E_{+}\cup E_{-}$ exists if and only if there is a path
$\mathcal{R}_{i}\leftarrow\mathcal{S}_{j}\leftarrow\mathcal{R}_{k}$ in the
directed SR-graph, so the R-graph is the graph obtained from the SR-graph by
removing the vertices corresponding to species and connecting reactions
directly.

In biological terms, edges of the R-graph describe the relationship between
reactions behavior. If the reactions cooperate, helping each other, the
reaction nodes are linked by a positive edge. This occurs when the chemical
species produced by a reaction is the reactant of an another one. Otherwise,
if the reactions compete, for example when they share the same reactants, the
reaction nodes are linked by a negative edge.

Note tht these definitions are very similar, but subtly different, from those
in (Angeli et al., 2006, 2010): indeed, in those papers the case of enzymes
(C1) is excluded _a priori_ , and the SR-graph is constructed in an undirected
version based only on $\alpha_{ij}$ and $\beta_{ij}$, neglecting the
difference between cases C2 and C3.

######  Example 2.

Consider the network

$\mathcal{R}_{1}:A\to C,\quad\mathcal{R}_{2}:B\to C,\quad\mathcal{R}_{3}:C\to
D,$ (3)

where all reactions are irreversible. According to our definition, there is no
edge between $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ in the R-graph. Under the
definition in (Angeli et al., 2006, 2010), there is a negative edge between
them, instead.

The following graph-theoretical characterization appears in (Angeli et al.,
2006).

######  Lemma 1.

Given a finite set of reactions $\mathcal{R}$ over a set of species
$\mathcal{S}$, a consistent labeling of $\mathcal{R}$ exists if and only if
every loop in the R-graph includes an even number of edges in $E_{-}$ (the
_positive-loop property_).

Indeed, to construct a consistent labeling, we can assign any value to
$\sigma(\mathcal{R}_{1})$ and then traverse the graph assigning signs so each
encountered vertex so that (2) holds; if the graph has the positive-loop
property we will never encounter a contradiction.

In addition, we can find a necessary condition for this labeling to exist.

######  Lemma 2 (Rule of 2).

For each species $\mathcal{S}_{j}$ in a reaction network, determine a number
$n(\mathcal{S}_{j})$ as follows. Consider all reactions in which
$\mathcal{S}_{j}$ is involved, and determine which of the cases C1, C2, C3
above holds.

  * •

Take the number of reactions in case C3.

  * •

Add 1 if there are one or more reactions that fall in case C1.

  * •

Add 1 if there are one or more reactions that fall in case C2.

Call this number $n(\mathcal{S}_{j})$ If, for some species $\mathcal{S}_{j}$,
$n(\mathcal{S}_{j})>2$, then $\mathcal{R}$ does _not_ admit a consistent
labeling.

###### Proof.

Suppose by contradiction that a consistent labeling exists, but that
$n(\mathcal{S}_{j})>2$ for a certain species $\mathcal{S}_{j}$. Then, we can
find reactions $i,k,\ell$ so that $\mathcal{R}_{i}$ falls in case C1 or C3,
$\mathcal{R}_{k}$ falls in case C2 or C3, and $\mathcal{R}_{\ell}$ falls in
case C3. Thanks to the definition of the three cases and to (2),

$\displaystyle\sigma(i)DR_{ij}\Gamma_{j\ell}\sigma(\ell)$ $\displaystyle>0,$
$\displaystyle\sigma(i)DR_{ij}\Gamma_{jk}\sigma(k)$ $\displaystyle>0,$
$\displaystyle\sigma(\ell)DR_{\ell j}\Gamma_{j\ell}\sigma(\ell)$
$\displaystyle>0.$

Taking the product of these three inequalities, one obtains $DR_{\ell
j}\Gamma_{j\ell}>0$, which contradicts our assumption (1). ∎

######  Example (2, continued).

In the network (3), the only species involved in more than one reaction is C.
We have $n(C)=2$, as it is involved in $\mathcal{R}_{3}$ under case C3 (1),
and in $\mathcal{R}_{1},\mathcal{R}_{2}$ under case C2 (+1). Thus Lemma 2 does
_not_ prohibit the existence of a consistent labeling. Indeed, one can verify
that
$\sigma(\mathcal{R}_{1})=\sigma(\mathcal{R}_{2})=\sigma(\mathcal{R}_{3})=1$ is
a consistent labeling. The same network (3), does not admit a consistent
labeling with the definitions in (Angeli et al., 2006, 2010), instead, as
Condition 2 in (Angeli et al., 2006, Proposition 1) fails with their
definition.

### Monotonicity in reaction coordinates

Let $x(t)=\int_{0}^{t}R(S(\tau))d\tau$ be the vector such that $x_{i}$ is the
extent of the $i$th reaction. The vector $x(t)$ solves, with initial condition
$x(0)=0$, the system of differential equations

$\frac{d}{dt}x(t)=R(S^{0}+\Gamma x(t)),$ (4)

where $S^{0}$ is the vector of initial concentrations.

The following result appears in (Angeli et al., 2006), and continues to hold
with our definition of R-graph (as proved in Section 4).

######  Theorem 3.

If the R-graph has the positive-loop property, then the system (4) is
_orthant-cooperative_ (or _orthant-monotone_), i.e., for a certain diagonal
matrix $\Sigma$, given two initial conditions $x^{0}$ and $\overline{x^{0}}$
such that $\Sigma(x^{0}-\overline{x^{0}})\geq 0$, then the solutions
$x(t),\overline{x}(t)$ of (4) with those initial conditions satisfy
$\Sigma(x(t)-\overline{x}(t))\geq 0$ for each time $t\geq 0$.

An explicit choice for this matrix $\Sigma$ is obtained by taking
$\Sigma_{ii}=\sigma(i)$, where $\sigma$ is a consistent labeling of
$\mathcal{R}$ (which exists thanks to the positive-loop property).

Note that the physical meaning of Theorem 3 is not apparent, since $x^{0}=0$
(i.e., empty history of reactions occurred before the initial time) is the
only initial condition that has a direct relevance to the application to CRNs.
Rather, this theorem is used in (Angeli et al., 2006, 2010) to obtain
statements about the asymptotic behavior of the reaction network.

### Handling reversible reactions

As in (Angeli et al., 2006), in the case of reversible reactions, we construct
the graph by considering only one of the two possible orientations. The choice
of reaction orientation does not influence Angeli’s monotonicity result that
we will exploit in the proof of our theorems. Moreover, even if this choice
changes the way we process the R-graph, as we will discuss later, it does not
influence the final result of our method.

Note that, as we already pointed out, the choice of the orientation of
reversible reactions does not influence the existence of a consistent
labeling.

######  Lemma 4.

Given a set of reactions $\mathcal{R}$ and a reaction
$\mathcal{R}_{i}\in\mathcal{R}$, let $\mathcal{R}^{rev}_{i}$ be the reaction
obtained from $\mathcal{R}_{i}$ by swapping reactants with products. There
exists a consistent labeling for the R-graph of $\mathcal{R}$ if and only if
there exists a consistent labeling for the R-graph of
$\mathcal{R}^{\prime}=\mathcal{R}\setminus\\{\mathcal{R}_{i}\\}\cup\\{\mathcal{R}^{rev}_{i}\\}$.

###### Proof.

Let $\langle\mathcal{R},E_{+},E_{-}\rangle$ be the R-graph of $\mathcal{R}$.
The R-graph of $\mathcal{R}^{\prime}$ can be obtained from the previous one as
follows: $\langle\mathcal{R}^{\prime},E^{\prime}_{+},E^{\prime}_{-}\rangle$
where $E^{\prime}_{+}=E_{+}\setminus\\{(\mathcal{R}_{i},\mathcal{R}_{j})\in
E_{+}\\}\cup\\{(\mathcal{R}^{rev}_{i},\mathcal{R}_{j})\mid(\mathcal{R}_{i},\mathcal{R}_{j})\in
E_{-}\\}$ and
$E^{\prime}_{-}=E_{-}\setminus\\{(\mathcal{R}_{i},\mathcal{R}_{j})\in
E_{-}\\}\cup\\{(\mathcal{R}^{rev}_{i},\mathcal{R}_{j})\mid(\mathcal{R}_{i},\mathcal{R}_{j})\in
E_{+}\\}$. Hence, the R-graph of $\mathcal{R}^{\prime}$ is the same as the one
of $\mathcal{R}$, but in which edges that were connected to $R_{i}$ are now
connected to $R^{rev}_{i}$ with opposite sign. Given a consistent labeling for
the R-graph of $\mathcal{R}$ it is possible to derive a consistent labeling
for $\mathcal{R}^{\prime}$ by simply assigning to $\mathcal{R}^{rev}_{i}$ an
opposite label with respect to $\mathcal{R}$, and vice-versa. ∎

### Sufficient condition for input-output monotonicity

Our main result is the following, proved in Section 4.

######  Theorem 5.

Let a set of chemical reactions $\mathcal{R}$ be given, with two distinct
species $\mathcal{S}_{I}$ and $\mathcal{S}_{O}$ designated as the input and
the output species.

Augment the network by adding two dummy reactions: an irreversible reaction
with output $\mathcal{S}_{I}$

$\mathcal{R}_{IN}:\emptyset\to\mathcal{S}_{I},$

and a reaction which has $\mathcal{S}_{O}$ as promoter

$\mathcal{R}_{OUT}:\mathcal{S}_{O}\to\mathcal{S}_{O}.$

If $\mathcal{R}\cup\\{\mathcal{R}_{IN},\mathcal{R}_{OUT}\\}$ has a consistent
labeling (i.e., if its R-graph has the positive-loop property), then
$\mathcal{S}_{O}$ is monotonic with respect to $\mathcal{S}_{I}$. In
particular, it is positively monotonic if $\sigma(IN)\sigma(OUT)=+1$, and
negatively monotonic if $\sigma(IN)\sigma(OUT)=-1$.

To show how Theorem 5 works, consider again our running Example 1, where we
want to study the monotonicity between the species S and the species P, which
we consider as the input and the output of the CRN.

We build the stoichiometric matrix of the network, considering only one
direction for the reversible reaction $\mathcal{R}_{1}$

$\Gamma=\bordermatrix{~{}&\mathcal{R}_{1}&\mathcal{R}_{2}\cr E&-1&+1\cr
S&-1&0\cr ES&+1&-1\cr P&0&+1\cr},$

and the matrix $DR$

$DR=\bordermatrix{~{}&E&S&ES&P\cr\mathcal{R}_{1}&+k_{1}[S]&+k_{1}[E]&-k_{1}[ES]&0\cr\mathcal{R}_{2}&0&0&+k_{2}&0\cr}.$

Figure 1: Directed SR-graph of Example 1, representing Michaelis-Menten
kinetics.

Then, following Definition 4, we build the directed SR-Graph, represented in
Figure 1, from which we derive the R-graph, which has one positive edge since
we obtain $DR_{ij}\Gamma_{jk}>0$ for each species $j$. We represent the
R-graph in Figure 2.

Figure 2: The R-graph of Example 1, representing Michaelis-Menten kinetics.

Then, we augment the network adding the two following dummy reactions:

$\displaystyle\begin{split}&\ch{0->S}\\\ &\ch{P->P},\\\ \end{split}$

and we verify trivially that the new network has a consistent labeling, which
means that output is monotonic with respect to the input. Moreover, since
$\sigma(IN)\sigma(OUT)=+1$ we can affirm that the output of this CRN, the
chemical species P, is positively monotonic with respect to the species S,
which means that if we increase the initial concentration of the input S, then
the concentration of the output P increases at any time.

As expected, this result is confirmed by simulations. Once we have fixed the
chemical rates $(k_{1}=0.1,k_{2}=1000,k_{3}=0.3)$, we vary the initial
concentration of the input S in the range $[100,2000]$. In Figure 3, we show
that by increasing the initial concentration of S (denoted $[S]_{O}$), the
concentration of P at the steady state (denoted $[P]_{ss}$) increases as well.
This actually does not hold only for the steady state: the concentration of P
is increased at any time point.

Figure 3: Simulation results of Example 1, representing Michaelis-Menten
kinetics. To show how the concentration of the species P is positively
monotonic with respect to the species S, we plot on the horizontal axis the
initial concentration of S, in a range $[100,2000]$ and on the vertical axis
the concentration of P at the steady state.

##  3 Case study

We now apply the result stated in Theorem 5 to the more complex network of ERK
signalling pathway.

### Signalling pathway example.

A signalling pathway usually consists of enzymatic cascades, having a starting
species that triggers the other connected processes. An initial stimulus,
perceived by a transductor (a sort of the first messenger), activates the
cascade amplifying the signal for the next enzymatic reaction. Many
biochemical processes are associated with signalling pathways as protein
activation, repression, and expression of genes: anomalies in these processes
could give rise to diseases like cancer, diabetes, and others.

One of the most important examples of such processes is the ERK pathway, which
is involved in growth, survival, proliferation, and differentiation of cells.
We consider the mathematical model of the ERK pathway implemented by Schilling
et al. (Schilling et al., 2009) and available to the public on the BioModels
Database (BIOMD0000000270). It consists of many fast phosphorylation
reactions, which spread the signal along the enzymatic cascade. For
simplicity, let us focus on a particular portion of the entire pathway, which
we will denote as ERK$*$. We indicate the species and the kinetics rates as
originally denoted in the model in (Schilling et al., 2009). For simplicity,
we refer to the reaction using the notation $R_{i}$, where $i$ is the kinetics
rate index. The reactions involved are the following:

$\displaystyle\begin{split}&\ch{Raf<=>[$k_{18}$][$k_{19}$]PRaf}\\\
&\ch{Mek1<=>[$k_{21}$[PRaf]][$k_{27}$]PMek1}\\\
&\ch{PMek1<=>[$k_{23}$][$k_{25}$]PPMek1},\\\ \end{split}$ (5)

in Table 1 we reported the coefficient rates and the initial conditions of
ERK$*$ system.

Table 1: The initial concentrations and the rates of ERK* system. Initial concentrations | Rates
---|---
Raf = 10 | $k_{18}=0.1445$
Praf = 0 | $k_{19}=0.37$
Mek1= 1 | $k_{21}=0.02$
PMek1 = 0 | $k_{23}=667.957$
PPMek1 = 0 | $k_{25}=0.13$
| $k_{27}=0.07$

The species PRaf, in the reaction $R_{21}$, is involved as catalyst promoter,
which means that its concentration positively influences the production of the
species PMek1. In our framework, the rate vector for this network is

$R(S)=\begin{bmatrix}k_{18}[Raf]-k_{19}[PRaf]\\\
k_{21}[Mek1][PRaf]-k_{27}[PMek1]\\\
k_{23}[PMek1]-k_{25}[PPMek1]\end{bmatrix}.$

We can apply Theorem 5 on the network, considering Raf and PPMek1,
respectively, as the input and the output. We augment the network adding the
two following dummy reactions:

$\displaystyle\begin{split}&\ch{0->Raf}\\\ &\ch{PPMek1->PPMek1}.\\\
\end{split}$

We build the R-graph and we can easily verify that it has a consistent
labeling, as we show in Figure 4. Moreover, since $\sigma(IN)\sigma(OUT)=+1$,
we can conclude that PPMek1 is positively monotonic with respect to Raf. This
can be confirmed with simulations, as shown in Figure 5.

Figure 4: R-graph representation of the network LABEL:erk_reaction1, augmented
by the two dummy reactions. As we can notice it has a consistent labeling.

Figure 5: Simulation results of CRN LABEL:erk_reaction1, representing ERK
signalling pathway. To show how the concentration of the species PPMek1 is
positively monotonic with respect to the species Raf, we plot on the
horizontal axis the initial concentration of Raf, in a range $[1,100]$ and on
the vertical axis the concentration of PPMek1 at the steady state.

##  4 Methods

In this section, we illustrate the proofs of the two main results.

###### Proof of Theorem 3.

This proof is basically the same one that appears in (Angeli et al., 2006),
but adapted to our slightly different definitions.

If the R-graph has the positive-loop property, then there is a consistent
labelling $\sigma$; let $\Sigma\in\mathbb{R}^{r\times r}$ be the diagonal
matrix with $\Sigma_{ii}=\sigma(i)$. The Jacobian of (4) is $J=DR\cdot\Gamma$.
The off-diagonal entries of $\Sigma J\Sigma$ are given by

$\Sigma_{ii}J_{ik}\Sigma_{kk}=\sum_{j}\sigma(i)DR_{ij}\Gamma_{jk}\sigma(k),\quad
i\neq k,$

where all summands are non-negative by the definition of consistent labeling
(2). The fact that $\Sigma J\Sigma$ has non-negative diagonal elements is a
sufficient condition for the system to be orthant-monotone (see (Smith,
1995)). ∎

###### Proof of Theorem 5.

Let $e_{I}\in\mathbb{R}^{s}$ (resp. $e_{O}$) be the vector that has an entry
$1$ in the corresponding to $\mathcal{S}_{I}$ (resp. $\mathcal{S}_{O}$) and
all other entries equal to $0$.

The rate and stoichiometry matrix of the augmented network are

$\widehat{DR}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left[\kern
0.0pt\kern-2.5pt\kern-5.55557pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\\\1\dots r$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle DR\\\IN$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0\\\OUT$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle ke_{O}^{T}$\hfil\kern
5.0pt\crcr}}}}\right]$}},\quad\widehat{\Gamma}=\hbox{}\vbox{\kern
0.86108pt\hbox{$\kern 0.0pt\kern 2.5pt\kern-5.0pt\left[\kern
0.0pt\kern-2.5pt\kern-5.55557pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
1\dots r$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle IN$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle OUT\\\$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\Gamma$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
e_{I}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern
5.0pt\crcr}}}}\right]$}}$

and its Jacobian is

$\widehat{J}=\widehat{DR}\cdot\widehat{\Gamma}=\hbox{}\vbox{\kern
0.86108pt\hbox{$\kern 0.0pt\kern 2.5pt\kern-5.0pt\left[\kern
0.0pt\kern-2.5pt\kern-5.55557pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
1\dots r$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle IN$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle OUT\\\1\dots r$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle J$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
DRe_{I}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0\\\IN$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0\\\OUT$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle e_{O}^{T}\Gamma$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt\crcr}}}}\right]$}}.$

Let $\sigma$ be a consistent labeling for the augmented network, and assume
without loss of generality (up to replacing $\sigma$ with $-\sigma$) that
$\sigma(IN)=1$.

Hence $\Sigma DRe_{I}\geq 0$, $\sigma(OUT)e_{O}^{T}\Gamma\Sigma\geq 0$, and
$\Sigma J\Sigma$ has non-negative diagonal entries.

We now take two initial values $S^{0}$ and $\overline{S^{0}}=S^{0}+\delta
e_{I}$ with $\delta>0$, and aim to prove that
$\sigma(OUT)(\overline{S}_{O}(t)-S_{O}(t))\geq 0$ for each $t\geq 0$; indeed,
this is the input-output monotonicity statement that we need to prove.

Define the dynamical system

$\begin{cases}\frac{d}{dt}y(t)=R(S^{0}+ze_{I}+\Gamma x(t)),\\\
\frac{d}{dt}z(t)=0.\end{cases}$ (6)

Let

$\check{J}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left[\kern
0.0pt\kern-2.5pt\kern-5.55557pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
y$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle z\\\y$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle J$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
DRe_{I}\\\z$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern
5.0pt\crcr}}}}\right]$}},\quad\check{\Sigma}=\hbox{}\vbox{\kern
0.86108pt\hbox{$\kern 0.0pt\kern 2.5pt\kern-5.0pt\left[\kern
0.0pt\kern-2.5pt\kern-5.55557pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
y$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle z\\\y$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\Sigma$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0\\\z$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 1$\hfil\kern 5.0pt\crcr}}}}\right]$}}.$

It is simple to verify that $\check{J}$ is the Jacobian of (6), and that
$\check{\Sigma}\check{J}\check{\Sigma}$ has non-negative off-diagonal, hence
the system (6) is orthant-monotone by the same result in (Smith, 1995) that we
have used in the proof of Theorem 3.

Direct verification shows that the solution of this system with initial value
$y(0)=0,z(0)=0$ is $y(t)=x(t),z(t)=0$, whereas the solution
$\overline{y}(t),\overline{z}(t)$ with initial value
$\overline{y}(0)=0,\overline{z}(0)=\delta$ is
$\overline{y}(t)=\overline{x}(t),\overline{z}(t)=\delta$, where the quantity
$\overline{x}(t)$ is defined analogously to (4) but with initial value
$\overline{S^{0}}$.

We have

$\begin{bmatrix}0\\\
\delta\end{bmatrix}=\check{\Sigma}\begin{bmatrix}\overline{y}(0)\\\
\overline{z}(0)\end{bmatrix}\geq\check{\Sigma}\begin{bmatrix}y(0)\\\
z(0)\end{bmatrix}=\begin{bmatrix}0\\\ 0\end{bmatrix}.$

By the orthant-monotonicity of (6), this implies

$\check{\Sigma}\begin{bmatrix}\overline{y}(t)\\\
\overline{z}(t)\end{bmatrix}\geq\check{\Sigma}\begin{bmatrix}y(t)\\\
z(t)\end{bmatrix}\quad\text{for all $t\geq 0$},$

and this means in particular that $\Sigma\overline{x}(t)\geq\Sigma x(t)$.

The concentrations of the output species $S_{O}(t),\overline{S}_{O}(t)$ at
time $t$ with the two initial conditions $S^{0},\overline{S^{0}}$ are given by

$S_{O}(t)=S^{0}_{O}+e_{O}^{T}\cdot\Gamma
x(t),\quad\overline{S}_{O}(t)=S^{0}_{O}+e_{O}^{T}\cdot\Gamma\overline{x}(t),$

respectively ($S^{0}_{O}=\overline{S^{0}}_{O}$ because
$\mathcal{S}_{O}\neq\mathcal{S}_{I}$), hence

$\sigma(OUT)(\overline{S}_{O}(t)-S_{O}(t))=\underbrace{\sigma(OUT)e_{O}^{T}\Gamma\Sigma}_{\geq
0}\underbrace{\Sigma(\overline{x}(t)-x(t))}_{\geq 0}\geq 0,$

which completes the proof. ∎

##  5 Discussion and Conclusions

In this paper, we proposed a new notion of monotonicity, namely the Input-
Output monotonicity, for which two species, considered as input and output of
the network, are monotonic if the variation of the initial concentration of
the input implies a monotonic variation in the concentration of the output.
Moreover, we established a new sufficient condition based on the R-graph that
guarantees that such monotonicity property holds. We showed how this new
notion of monotonicity can have a great practical impact. Indeed, monotonicity
assessment can drastically reduce the number of simulations necessary to study
the dynamical behaviour of a chemical reaction network under uncertain initial
conditions. This can be very useful in several different cases, for example
when initial concentrations are not exactly known or when the focus is on
studying the effects of perturbations of initial concentrations (as in the
case of robustness). The proposed method allows us to study complex chemical
reaction networks, performing a preliminary analysis of the system dynamics,
in order to investigate different biological properties, such as robustness or
sensitivity to perturbations.

We have shown the application of our approach to the small example of
Michaelis-Menten kinetics and to the quite complex model of the ERK signaling
pathway (Kwang-Hyun et al., 2003).

In order to apply our sufficient condition on larger and more complex network,
different approaches of model reduction can be applied. For instance, in
(Küken et al., 2021) the authors show how to remove particular nodes of the
network preserving the steady state fluxes of the system. Moreover, for
analysing large-scale biochemical models, we can also simplify the model using
the common approach of the separation of timescale (Ingalls, 2013), which
allows us to consider the reaction networks as divided into processes having
different timescales leading to an approximate – but accurate – version of the
original model.

In addition, in (Bove et al., 2020), the authors propose an innovative
approach based on machine learning on graphs to predict whether a biological
system is robust studying its topological features. As future work, we intend
to apply a similar method to predict the monotonicity of a system.

## Author’s contributions

Conceptualization & Formal analysis: R. Gori, P. Milazzo, L. Nasti, F. Poloni.
Simulations & Experiments: L. Nasti. Writing original draft: R. Gori, P.
Milazzo, L. Nasti, F. Poloni. Writing - review & editing: R. Gori, P. Milazzo,
L. Nasti, F. Poloni.

## Funding

R. Gori, P. Milazzo, L. Nasti are partially supported by the University of
Pisa’s project PRA_2020_26 “Metodi informatici integrati per la biomedica”. F.
Poloni is partially supported by the University of Pisa’s project PRA_2020_61
“Analisi di reti complesse: dalla teoria alle applicazioni” and by GNCS/INDAM
(Istituto Nazionale di Alta Matematica).

## References

  * Alur et al. (1992) Alur, R., Courcoubetis, C., Henzinger, T. A., & Ho, P.-H. 1992, in Hybrid systems, ed. R. L. Grossman, A. Nerode, A. P. Ravn, & H. Rischel, 209–229
  * Angeli et al. (2010) Angeli, D., De Leenheer, P., & Sontag, E. 2010, Journal of mathematical biology, 61, 581
  * Angeli et al. (2006) Angeli, D., De Leenheer, P., & Sontag, E. D. 2006, in Decision and Control, 2006 45th IEEE Conference, IEEE, 7–12
  * Angeli et al. (2007) Angeli, D., De Leenheer, P., & Sontag, E. D. 2007, Mathematical biosciences, 210, 598
  * Barbuti et al. (2011) Barbuti, R., Maggiolo-Schettini, A., Milazzo, P., Pardini, G., & Tesei, L. 2011, Natural Computing, 10, 3
  * Behinaein et al. (2014) Behinaein, B., Rudie, K., & Sangrar, W. 2014, in Electrical and Computer Engineering (CCECE), 2014 IEEE 27th Canadian Conference on, ed. IEEE
  * Bernini et al. (2018) Bernini, A., Brodo, L., Degano, P., Falaschi, M., & Hermith, D. 2018, Natural Computing, 17, 345
  * Bove et al. (2020) Bove, P., Micheli, A., Milazzo, P., & Podda, M. 2020, in 13th International conference on bioinformatics models, methods and algorithms (BIOINFORMATICS 2020), SCITEPRESS, 32–43
  * Danos et al. (2008) Danos, V., Feret, J., Fontana, W., & Krivine, J. 2008, in International Workshop on Verification, Model Checking, and Abstract Interpretation, ed. Springer, 83–97
  * Dresch et al. (2010) Dresch, J. M., Liu, X., Arnosti, D. N., & Ay, A. 2010, BMC systems biology, 4, 142
  * Feinberg (1987) Feinberg, M. 1987, Chemical Engineering Science, 42, 2229
  * Gori et al. (2019) Gori, R., Milazzo, P., & Nasti, L. 2019, in 10th International conference on bioinformatics models, methods and algorithms (BIOINFORMATICS 2019), ed. SciTePress, 250–257
  * Henzinger (2000) Henzinger, T. A. 2000, in Verification of Digital and Hybrid Systems, ed. Springer, Vol. 736, 265–292
  * Ingalls (2013) Ingalls, B. P. 2013, Mathematical modeling in systems biology: an introduction (MIT press)
  * Iooss & Lemaître (2015) Iooss, B., & Lemaître, P. 2015, in Uncertainty management in simulation-optimization of complex systems, ed. Springer, 101–122
  * Kitano (2002) Kitano, H. 2002, Systems biology: towards systems-level understanding of biological systems. Kitano H ed, Foundations of Systems Biology, MIT Press, Cambridge
  * Kitano (2007) —. 2007, Molecular systems biology, 3
  * Koch (2010) Koch, I. 2010, Molecular Informatics, 29, 838
  * Küken et al. (2021) Küken, A., Wendering, P., Langary, D., & Nikoloski, Z. 2021, bioRxiv, doi: 10.1101/2021.03.17.435785
  * Kwang-Hyun et al. (2003) Kwang-Hyun, C., Sung-Young, S., Hyun-Woo, K., et al. 2003, in International Conference on Computational Methods in Systems Biology, ed. Springer, 127–141
  * Li et al. (2017) Li, X., Omotere, O., Qian, L., & Dougherty, E. R. 2017, EURASIP J. on Bioinformatics and Systems Biology, 2017, 8
  * Murata (1989) Murata, T. 1989, Proceedings of the IEEE, 77, 541
  * Nasti (2020) Nasti, L. 2020, PhD thesis, Universitità di Pisa, Dipartimento di Informatica
  * Nasti et al. (2018) Nasti, L., Gori, R., & Milazzo, P. 2018, in Federation of International Conferences on Software Technologies: Applications and Foundations, Springer, 81–97
  * Regev et al. (2000) Regev, A., Silverman, W., & Shapiro, E. 2000, in Biocomputing 2001, ed. W. Scientific, 459–470
  * Rizk et al. (2011) Rizk, A., Batt, G., Fages, F., & Soliman, S. 2011, Theoretical Computer Science, 412, 2827
  * Schilling et al. (2009) Schilling, M., Maiwald, T., Hengl, S., et al. 2009, Molecular systems biology, 5
  * Shinar et al. (2009) Shinar, G., Alon, U., & Feinberg, M. 2009, SIAM Journal on Applied Mathematics, 69, 977
  * Shinar & Feinberg (2010) Shinar, G., & Feinberg, M. 2010, Science, 327, 1389
  * Smith (1995) Smith, H. L. 1995, Mathematical Surveys and Monographs, Vol. 41, Monotone dynamical systems (American Mathematical Society, Providence, RI), x+174
  * Zi (2011) Zi, Z. 2011, IET systems biology, 5, 336


# A Standard Model Explanation for the “ATOMKI Anomaly”

A. Aleksejevs  Grenfell Campus, Memorial University of Newfoundland, Corner
Brook, NL, Canada  S. Barkanova  Grenfell Campus, Memorial University of
Newfoundland, Corner Brook, NL, Canada  Yu.G. Kolomensky  Department of
Physics, University of California, Berkeley, CA, USA  Nuclear Science
Division, Lawrence Berkeley National Laboratory, Berkeley, CA, USA  B. Sheff
Department of Physics, University of Michigan, Ann Arbor, MI, USA

###### Abstract

Using the $e^{+}e^{-}$ pair spectrometer at the 5 MV Van de Graaff accelerator
at the Institute for Nuclear Research, Hungarian Academy of Sciences (ATOMKI),
Krasznahorkay et al. have claimed a 6.8$\sigma$ excess at high $e^{+}e^{-}$
opening angles in the internal pair creation isoscalar transition
${}^{8}\mathrm{Be}(18.15)\to\,^{8}\mathrm{Be}\,e^{+}e^{-}$. A hypothetical
gauge boson with the mass circa $17$ MeV, “X17”, has been proposed as an
explanation for the excess. We show that the observed experimental structure
can be reproduced within the Standard Model by adding the full set of second-
order corrections and the interference terms to the Born-level decay
amplitudes considered by Krasznahorkay et al. We implement a detailed model of
the ATOMKI detector, and also show how experimental selection and acceptance
bias exacerbate the apparent difference between the experimental data and the
Born-level prediction.

## Introduction

The Standard Model (SM) of particle physics has survived many challenges over
the past fifty years. Occasionally, experimental observations disagree with
the SM predictions; such deviations are often reconciled after additional
theoretical or experimental scrutiny. One such deviation is claimed by an
experiment at the Institute for Nuclear Research, Hungarian Academy of
Sciences (ATOMKI). In this experiment, a 5 MV Van De Graaf-accelerator was
used to produce an 18.15 MeV excited state of ${}^{8}\mathrm{Be}(18.15)$
(called ${}^{8}\mathrm{Be}^{*}$ henceforth) with a subsequent internal pair
conversion (IPC) to the ground state of ${}^{8}\mathrm{Be}$:
${}^{7}\mathrm{Li}(p,e^{+}e^{-})^{8}\mathrm{Be}$ Announcement . This channel
is of particular interest, as the creation of $e^{+}e^{-}$ pairs provides a
potential avenue for detection of low-mass boson candidates that are not
present in the Standard Model Theory . A pair spectrometer was constructed to
focus on detecting the IPC process with a fine resolution in the opening angle
of $e^{+}e^{-}$ pairs $\theta_{+-}$ Detector .

Recent results from this experiment indicate an excess of events at large
$\theta_{+-}$ over the leading-order theoretical predictions, consistent with
a new boson with a mass of $16.7$ MeV Announcement . This experimental
evidence was found to be inconsistent with a dark photon hypothesis Theory .
Feng et al. proposed a new protophobic, fifth force gauge boson, which has
sparked a flurry of additional model-building work and interest in popular
press. In this paper, we discuss our effort to critically analyze the idea
that a non-SM, resonant process is needed to produce the observed excess in
$e^{+}e^{-}$ production at high $\theta_{+-}$. In order to fully understand
the interplay between the theoretical and experimental effects, we construct a
Monte Carlo (MC) model of the detector used at ATOMKI. We use this model to
demonstrate that an interference between the Born-level IPC amplitude and a
subleading nonresonant component (an amplitude with a broad phase-space
structure) could produce effects observed by ATOMKI SheffThesis . Similar
ideas have been explored by other authors Kalman:2020meg ; Zhang:2017zap . The
most basic subleading contribution to the IPC process arises from the second
order electromagnetic corrections, i.e. contributions beyond the Born
approximation: two-photon box diagrams, vertex corrections, etc. While the
calculation of the higher-order corrections is technically challenging, one
would naively expect such contributions to be suppressed by a factor of
$Z\alpha$ compared to the Born term. However, the interference between the box
and tree-level diagrams and the nontrivial structure of the box diagrams,
including kinematic singularities and the presence of excited nuclear states
in the box, can produce unexpected effects. In this paper, we report the full
second-order calculation of the
${}^{8}\mathrm{Be}^{*}\to\,^{8}\mathrm{Be}\,e^{+}e^{-}$ process, and compare
the results to the experimental distributions reported by ATOMKI.

## NLO QED Model of ${}^{8}\mathrm{Be}^{*}$ Decay

In the rest frame of ${}^{8}\mathrm{Be}^{*}$, the doubly-differential decay
rate comprises from a phase-space component and a square of the matrix
element, which has a leading-order (LO) contribution $\left|M_{LO}\right|^{2}$
and an interference term between LO and next-to-the-leading (NLO) matrix
element, $2\Re\left[M_{LO}M_{NLO}^{*}\right]$:

$\frac{d^{2}\Gamma}{d\theta_{+}d\theta_{-}}=\left(\left|M_{LO}\right|^{2}+2\Re\left[M_{LO}M_{NLO}^{*}\right]\right)\Phi,$
(1)

where details on phase space element $\Phi$ are given in supplemental part of
the paper.

At the energy scale of the ATOMKI experiment, an effective operator approach
can be employed, with ${}^{8}\mathrm{Be}^{*}$ and ${}^{8}\mathrm{Be}$ as the
fundamental degrees of freedom. Specifically, for the transition
${}^{8}\mathrm{Be}^{*}\rightarrow{}^{8}\mathrm{Be}+\gamma$, we can write:

$\displaystyle\mathfrak{L}^{(1)}$
$\displaystyle=\frac{e}{\Lambda_{\gamma}}\epsilon^{\mu\nu\alpha\beta}B_{\mu\nu}^{*}F_{\alpha\beta}B,$
(2)

where $B_{\mu\nu}^{*}=\partial_{\mu}B_{\nu}^{*}-\partial_{\nu}B_{\mu}^{*}$ and
$F_{\alpha\beta}=\partial_{\alpha}A_{\beta}-\partial_{\beta}A_{\alpha}$ are
the usual field strength tensors for ${}^{8}\text{Be}^{*}$ and electromagnetic
fields, respectively. The constant $\Lambda_{\gamma}$ is related to the
transition nuclear matrix element and is used as a scaling parameter. However,
Eq. (2) is only describing the tree-level transition, which is insufficient to
explain a small peak at $\sim 16.7$ MeV of electron-positron invariant mass.
Thus, we need to extend our analysis from the Born (LO) matrix element of
order $\alpha$ to the next-to-leading-order (NLO) QED contributions, which
requires one-loop calculations. The four categories of Feynman graphs
corresponding to $\alpha^{2}$ transitions are shown on Fig. 1.

Figure 1: Higher order QED contributions in decay of ${}^{8}\mathrm{Be}^{*}$.
Grey bubble corresponds to spin $1\rightarrow 0$ transition, black bubble for
spin $0\rightarrow 0$ and white bubble for spin $1\rightarrow 1$ transitions.

Group (a) in Fig. 1 is the vacuum polarization contributions, which we treat
by splitting the leptonic and hadronic contributions and using effective
masses of light quarks to calculate the hadronic contributions. Group (b) in
Fig. 1) is the electron current vertex corrections. Groups (c) and (d) are
represented by boxes with ${}^{8}\mathrm{Be}$ and ${}^{8}\mathrm{Be}^{*}$ in
the loop, respectively. The diagrams (b)-(d) are infrared-divergent. To treat
infrared (IR) divergences in the vertex correction graph (c), we use the soft-
photon approximation in bremsstrahlung diagrams and introduce a cut on energy
of soft photons given by the threshold conditions in the electron-positron
pair production. For the boxes, we take into account only the infrared-finite
part. Since ATOMKI does not distinguish between the electrons and positrons,
we sum over events in which either electron or positron are observed in the
same detector element in the MC simulation. Due to the CP symmetry of the
underlying QED lagrangian, this cancels out the IR-divergence in the box
diagrams.

Besides the QED transition
${}^{8}\mathrm{Be}^{*}\rightarrow^{8}\mathrm{Be}+\gamma$, in the box diagrams,
we also use

$\displaystyle\mathfrak{L}^{(2)}=-ieZ\left[B^{*\mu\nu}B_{\mu}^{*\dagger}A_{\nu}-B^{*\mu\nu\dagger}B_{\mu}^{*}A_{\nu}+B_{\mu}^{*}B_{\nu}^{*\dagger}F^{\mu\nu}\right]$
(3)

to describe a coupling
${}^{8}\mathrm{Be}^{*}\rightarrow^{8}\mathrm{Be}^{*}+\gamma$, and

$\displaystyle\mathfrak{L}^{(3)}=-ieZ\left[B\partial_{\mu}B^{\dagger}-\partial_{\mu}BB^{\dagger}\right]A^{\mu}$
(4)

to introduce coupling ${}^{8}\mathrm{Be}\rightarrow^{8}\mathrm{Be}+\gamma$. In
both lagrangians, $Z$ is equal to four. Calculations are done in several
stages. Based on the lagrangian densities in Eqs. (2-4) we develop the model
file, which we used in FeynArts feynarts  to generate one-loop topologies and
produce matrix elements, which we evaluate using Passarino-Veltman basis in
the FormCalc frmcloop  package. Numerical calculations are carried out with
the help of LoopTools frmcloop  and confirmed by Collier collier  packages.

Figure 2: The interference term in Eq. (1)
($2\Re\left[M_{LO}M_{NLO}^{*}\right]$) as a function of $e^{+}e^{-}$ invariant
mass $\omega$ for electron (top) and positron (bottom) angles fixed to
$\theta_{-}(\theta_{+})=40^{\circ}$. Dashed line (yellow) represents vacuum
polarization graphs (a), dot-dashed line (green) is the result of vertex
correction graph (b) with soft-photon bremsstrahlung treatment. Dotted graph
(blue) is the IR-finite part of the boxes (c)-(d).

The results of the higher-order QED corrections for the cases when one of the
angles (either electron or positron) is fixed are given in Fig. 2. Self-
energies and vertex corrections are completely symmetric with respect to
electron-positron flip. On the contrary, boxes show asymmetrical behavior and
almost cancel out each other when we combine the electron and positron events.
The look-up tables for the MC simulations are available upon request.

## Monte Carlo Simulation

The momenta and angles of the final state particles ${}^{8}\mathrm{Be}$,
$e^{+}$, and $e^{-}$ are generated over the full three-body phase space with
the probability proportional to the differential decay rate in Eq. (1). The
kinematics are taken in the rest frame of ${}^{8}\mathrm{Be}^{*}$, neglecting
its boost due to the initial proton momentum. $e^{+}$ and $e^{-}$ are assumed
to be indistinguishable in the detector, and we ignore the 4-vector of
${}^{8}\mathrm{Be}$ in the subsequent steps. The final momenta of $e^{+}e^{-}$
are then fed into the detector simulation and weighted accordingly. We then
histogram over the opening angle and over the invariant mass of the
$e^{+}e^{-}$ pair and compare against the experimental data.

The ATOMKI detector model was constructed based on description of Detector .
We assume each detector block to be a rectangle in $\theta\phi$ space with
perfect efficiency. The center of each rectangle is placed at
$\phi=0,\pi/2,5\pi/6,7\pi/6,3\pi/2$ Detector . We tune the angular acceptance
of each block by simulating a uniform angular distribution for each particle
in $e^{+}e^{-}$ pairs, and matching the angular acceptance of the $e^{+}e^{-}$
pairs to the released detector response curve published Detector , as shown in
Fig. 3. The model represents the measured detector acceptance very well for
angles larger than $40^{\circ}$, with comparable variance from data to the
published expected detector response. The final block sizes are set to
$\Delta\theta=0.59$ rad and $\delta\phi=0.58$ rad, with the uncertainty of
roughly $0.02$ rad in each dimension. These values are then applied to the MC
simulation by applying a 0 weight to any $e^{+}e^{-}$ pair that includes at
least one particle not at an angle that crosses one of the detectors.

Figure 3: Comparison of simulated and measured detector responses to uniformly
distributed particle pairs, along with the published expected response.
$\chi^{2}$ calculated with all standard errors assumed to be 1 for comparison
of relative variance of new model and ATOMKI’s model

Based on Announcement , an additional cut is applied to the detected particle
pair, removing particularly asymmetric $e^{+}e^{-}$ pairs. The cut was applied
to only allow $|y|<0.5$, where disparity
$y=\left(E_{e^{-}}-E_{e^{+}}\right)/\left(E_{e^{-}}+E_{e^{+}}\right)$

We note that while $e^{+}$ and $e^{-}$ are indistinguishable in the ATOMKI
detector, the response of the plastic scintillator to electrons and positrons
is different eePlastic . It is not clear if this difference was incorporated
in the original calibration of the ATOMKI apparatus Detector . We approximate
the difference in response by adding an additional 1 MeV of visible energy to
the positron signal before applying the disparity cut. This introduces a small
asymmetry in the detector response.

## Analysis

We fit the simulated distributions of $e^{+}e^{-}$ opening angles to the
ATOMKI data using a binned $\chi^{2}$ fit by varying two parameters: the
overall normalization of the Born contribution, and the coefficient of the
interference term between the Born and the second-order diagrams. The sign and
the magnitude of the interference term are allowed to float, allowing for a
possible phase shift between the two terms. The best fit reproduces the ATOMKI
spectrum with a $\chi^{2}=38$ for $7$ degrees of freedom. We note that the
largest $\chi^{2}$ contribution comes from the lowest $\theta_{+-}$ bin, where
our assumption of the uniform detector efficiency may be incorrect. The best-
fit value for the interference coefficient is $\epsilon=57.72\pm 0.06$,
indicating that the box diagram propagator may have contributions from more
than one $J^{P}=1^{+}$ state (e.g. 18.15 MeV and 17.64 MeV states). This fit
is shown in Fig. 4.

Figure 4:  Results of a Monte Carlo simulation using underlying particle
distributions from the calculations above and applying detector acceptance and
disparity limiting cuts. Shown are the (orange dashed) tree level calculation,
(green dot-dashed) loop-level, (blue solid) their fitted sum, as compared to
(purple dotted) the published result Announcement  and (red points) the data
collected in the same source.

### Mass Distribution

With a reasonable agreement between the SM fit and the experimental data in
the X17 ”signal” region at the large values of $\theta_{+-}$, we compute the
$e^{+}e^{-}$ invariant mass distribution. First, for each of the two angular
distributions used in the fit, a histogram of the $e^{+}e^{-}$ invariant mass
is produced. We use the X17 contribution reported in Announcement  as the
detector mass resolution at $\bar{m}=16.7$ MeV, and scale that resolution with
mass according to

$\text{g}_{m_{0}}(m)=\text{f}\left(\frac{m-m_{0}}{m_{0}}\bar{m}+\bar{m}\right)$
(5)

where $m_{0}$ is the true $e^{+}e^{-}$ mass and $m$ is the measured mass. The
theoretically predicted mass distribution is then smeared with this resolution
function.

The resulting mass distribution is shown in Fig. 5. We do not attempt to make
a quantitative comparison with the ATOMKI mass distribution Announcement ,
since it depends critically on the detector energy resolution function,
detection efficiency as a function of $e^{-}$ and $e^{+}$ energies, and
details of the detector calibration. Such details are not provided in the
experimental publications Announcement ; Detector . However, we note a
striking qualitative feature: a bump-like structure in the predicted invariant
mass distribution that does not require a new boson. This shows how non-
resonant SM processes can combine to produce a false bump in the detector.

Figure 5:  Distribution of electron positron pair creation over the invariant
mass. Significant peak structure shown arising from interference between loop
and tree level effects. Relative magnitudes of the components are taken from
the fit in Fig. 4.

### Angular Distribution Bias

Of interest in this investigation are the origins of the bump-like structure
in the angular and mass distributions in Fig. 4-5. Much of this effect appears
to arise from non-uniformities in the angular acceptance of the ATOMKI
spectrometer. Running the simulation with complete $4\pi$ acceptance, and
setting the interference coefficient to be the same as in the above fit, we
arrive at the results shown in Fig. 6-7. As shown, the bump structure is
highly suppressed, implying that it is a result of an interplay between
systematic biases in the detector and the effects of loop diagrams.

Figure 6:  Angular separation distribution with detector set to full $4\pi$
acceptance. Energy disparity cuts still included. Shown are the (orange
dashed) tree level calculation, (green dot-dashed) loop-level, and (blue
solid) their sum. Relative magnitude of each component is taken from the fit
found in Fig. 4. Figure 7:  Mass Distribution with the detector set to full
$4\pi$ acceptance. Energy disparity cuts still included. Shown are the (orange
dashed) tree level calculation, (green dot-dashed) loop-level, and (blue
solid) their sum. Relative magnitude of each component is taken from the fit
found in Fig. 4

## Conclusions

We show that the ATOMKI data Announcement  can be reasonably described by
taking into account the second-order contributions to the IPC process, and
accounting for detector and analysis bias. As these are non-resonant, Standard
Model processes, they do not require the addition of a new boson. Clearly, an
element of caution must be taken before confidence can be placed in an
entirely new particle having been discovered.

An independent measurement of the IPC spectra in the ${}^{8}\mathrm{Be}$ and
other isoscalar nuclear systems is clearly called for. Such measurements
should aim to measure the IPC distributions with large acceptance, minimal
detector bias, and ideally better invariant mass resolution, in order to be
able to discriminate between a resonant and a non-resonant contributions.
Tagging the initial nuclear state in order to eliminate initial-state
interference effects Kalman:2020meg  would also be of value.

## Acknowledgments

We would like to thank the members of UC Berkeley and Lawrence Berkeley
National Lab (LBNL) Weak Interactions Group for their support in this work,
and Gerald Miller and Xilin Zhou for stimulating discussions. This work was
supported in part by the Natural Sciences and Engineering Research Council of
Canada (NSERC), by the National Science Foundation (NSF), and by the US
Department of Energy (DOE) Office of Science. B. Sheff thanks UC Berkeley
Institute for International Studies Merit Scholarship Program and UC Berkeley
Physics department for their support. A. Aleksejevs and S. Barkanova would
also like to thank UC Berkeley and LBNL for hospitality and support.

## References

  * [1] A.J. Krasznahorkay et al., Phys. Rev. Lett., 116, 042501 (2016).
  * [2] J.L. Feng et al., Phys. Rev. Lett. 117, 071803 (2016).
  * [3] J. Gulyás et al., Nucl. Instrum. Methods Phys. Res. A 808, 21 (2016).
  * [4] B. Sheff, Honors Thesis, UC Berkeley, 2017 (unpublished).
  * [5] P. Kálmán and T. Keszthelyi, Eur. Phys. J. A 56, 205 (2020).
  * [6] X. Zhang and G. A. Miller, Phys. Lett. B 773, 159 (2017).
  * [7] T. Hahn, Comput. Phys. Commun. 140, 418 (2001) [arXiv:hep-ph/0012260 [hep-ph]].
  * [8] T. Hahn and M. Perez-Victoria, Comput. Phys. Commun. 118, 153 (1999).
  * [9] A. Denner, S. Dittmaier and L. Hofer, Comput. Phys. Commun. 212, 220 (2017).
  * [10] W.J. Meiring, J. van Klinken, and V.A. Wichers, Phys. Rev. A 44, 2960 (1991).

## Appendix

The phase space for doubly differential decay rate in Eq. (1) can be written
as follows:

$\Phi=\frac{1}{\left(2\pi\right)^{3}}\frac{J}{8m_{Be^{*}}}\frac{E_{-}^{2}}{E_{Be}}\frac{\sin\theta_{-}}{\sin\theta_{+}}$

and

$\displaystyle
J=\frac{\sin^{2}\left(\theta_{+}+\theta_{-}\right)}{2R\sin^{2}\theta_{-}}\left[\frac{m_{Be*}}{\left(\cos\left(\theta_{+}+\theta_{-}\right)-1\right)}\left(\frac{m_{Be^{*}}\sin\left(\theta_{+}+\theta_{-}\right)}{\tan\theta_{-}}-R\right)-m_{Be}^{2}\right],$
$\displaystyle
R=\left\\{m_{Be}^{2}+\frac{\sin\theta_{+}}{\sin\theta_{-}}\left(2m_{Be^{*}}^{2}+m_{Be}^{2}\frac{\sin\theta_{+}}{\sin\theta_{-}}\right)+\left(m_{Be^{*}}^{2}-m_{Be}^{2}\right)\frac{\sin\left(2\theta_{+}+\theta_{-}\right)}{\sin\theta_{-}}\right\\}^{1/2},$

where subscripts $+/-$ correspond to the positron and the electron,
respectively.

Figure 8: Phase space of three body decay for
${}^{8}\mathrm{Be}^{*}\rightarrow^{8}\mathrm{Be}e^{+}e^{-}$. Figure 9:
Separation angle between electron and positron as a function of invariant mass
of the $e^{+}e^{-}$ system and the $e^{+}~{}^{8}\mathrm{Be}$ system

The phase space of ${}^{8}\mathrm{Be}$ and the electron-positron pair in the
final state is shown on Fig. 8, demonstrating the characteristic three-body
decay behavior. The correlation between the invariant masses in the Dalitz
plot plane and the $e^{+}e^{-}$ opening angle is shown in Fig. 9. The matrix
elements projected on the 3-body phase space Dalitz plot are shown in Fig. 10.

Figure 10:  $\alpha^{3}$ order matrix element contribution as a function of
invariant mass of the $e^{+}e^{-}$ system and the $e^{+}~{}^{8}\mathrm{Be}$
system on linear (left) and log (right) scales