# Sums of three Fibonacci numbers close to a power of 2

Bibhu Prasad Tripathy and Bijan Kumar Patel

###### Abstract

In this paper, we find all the sums of three Fibonacci numbers which are close
to a power of 2. This paper continues and extends the previous work of
Hasanalizade [8].

Keywords: Fibonacci number, linear forms in logarithms, reduction method.
2010 Mathematics Subject Classification: 11B39; 11J86.

##  1 Introduction

The Fibonacci sequence $\\{F_{n}\\}_{n\geq 0}$ is the binary recurrence
sequence given by

$F_{n+2}=F_{n+1}+F_{n}~{}~{}{\rm for}~{}n\geq 0$

with initials $F_{0}=0$ and $F_{1}=1$. The Fibonacci numbers are famous for
possessing wonderful and amazing properties. The problem of determining all
integer solutions to Diophantine equations with Fibonacci numbers has piqued
the curiosity of mathematicians and there is a wide literature on this
subject. For instance, Bugeaud et al. [4] studied the Fibonacci numbers, that
is, perfect powers and found that 1, 2 and 8 are the only powers of 2 in the
Fibonacci sequence. Later, Bravo and Luca [3] solved the Diophantine equation
$F_{n}+F_{m}=2^{a}$ in non-negative integers $n,m$ and $a$ with $n\geq m$.
Following that Bravo and Bravo [1] found all solutions of the Diophantine
equation $F_{n}+F_{m}+F_{l}=2^{a}$ in non-negative integers $n,m,l$ and $a$
with $n\geq m\geq l$. Consequently, Chim and Ziegler [6] determined all
solutions of the Diophantine equations

$F_{n_{1}}+F_{n_{2}}=2^{a_{1}}+2^{a_{2}}+2^{a_{3}}$

and

$F_{m_{1}}+F_{m_{2}}+F_{m_{3}}=2^{t_{1}}+2^{t_{2}}$

in non-negative integers
$n_{1},n_{2},m_{1},m_{2},m_{3},a_{1},a_{2},a_{3},t_{1},t_{2}$ with $n_{1}\geq
n_{2},m_{1}\geq m_{2}\geq m_{3},a_{1}\geq a_{2}\geq a_{3}$ and $t_{1}\geq
t_{2}$.

An integer $n$ is said to be close to a positive integer $m$ if it satisfies

$|n-m|<\sqrt{m}.$

By using the above concept Chern and Cui [5] found all the Fibonacci numbers
which are close to the power of 2. More precisely they found all solutions of
the inequality

$|F_{n}-2^{m}|<2^{m/2}.$

Later, Bravo et al. [2] extended the previous work of [5] by considering the
$k$-generalized Fibonacci sequence $\\{F_{n}^{(k)}\\}$ and solved the
Diophantine equation

$|F_{n}^{(k)}-2^{m}|<2^{m/2}$

in non-negative integer $n,k,m$ with $k\geq 2$ and $n\geq 1$.

Recently, Hasanlizade [8] extended the previous work of [5] by considering the
sum of Fibonacci numbers and studied the sum of two Fibonacci numbers close to
a power of 2. In particular, he solved

$|F_{n}+F_{m}-2^{a}|<2^{a/2},$

in positive integers $n,m$ and $a$ with $n\geq m$.

Motivated by the above literature, we extend the work of [5] and [8] search
for the sum of three Fibonacci numbers which are close to the power of 2. More
specifically, we study the Diophantine inequality

$|F_{n}+F_{m}+F_{l}-2^{a}|<2^{a/2},$ (1.1)

in positive integers $n,m,l$ and $a$ with $n\geq m\geq l$.

In particular, our main result concerning (1.1) is the following.

######  Theorem 1.1.

There are exactly 214 solutions $(n,m,l,a)$ to Diophantine inequality (1.1).
All the solutions satisfy $n\leq 42$ and $a\leq 28$. A list of solutions is
given in the appendix.

##  2 Auxiliary results

In this part, we will review several well-known results which will be used
after. The Binet formula for Fibonacci sequence is

$F_{n}=\frac{\alpha^{n}-\beta^{n}}{\alpha-\beta}\quad(n\geq 0),$ (2.2)

where $\alpha=\frac{1+\sqrt{5}}{2}$ and $\beta=\frac{-1}{\alpha}$ are the
roots of the characteristic equation $x^{2}-x-1=0$.

######  Lemma 2.1.

For every positive integer $n\geq 1$, we have

$\alpha^{n-2}\leq F_{n}\leq\alpha^{n-1}.$ (2.3)

######  Lemma 2.2.

([10, Eq. 5]) For every positive integer $n>1$, we have

$0.38\alpha^{n}\leq F_{n}\leq 0.48\alpha^{n}.$

To prove our main result, we use a few times a Baker-type lower bound for non-
zero linear forms in logarithms of algebraic numbers. We state a result of
Matveev [9] about the general lower bound for linear forms in logarithms.
Applying a version of the Baker-Davenport reduction method we reduce the large
bound, but first, recall some basic notations from algebraic number theory.

Let $\eta$ be an algebraic number of degree $d$ with minimal primitive
polynomial

$f(X):=a_{0}X^{d}+a_{1}X^{d-1}+\cdots+a_{d}=a_{0}\prod_{i=1}^{d}(X-\eta^{(i)})\in\mathbb{Z}[X],$

where the $a_{i}$’s are relatively prime integers, $a_{0}>0$, and the
$\eta^{(i)}$’s are conjugates of $\eta$. Then

$h(\eta)=\frac{1}{d}\left(\log
a_{0}+\sum_{i=1}^{d}\log\left(\max\\{|\eta^{(i)}|,1\\}\right)\right)$ (2.4)

is called the _logarithmic height_ of $\eta$. In particular, if $\eta=p/q$ is
a rational number with gcd$(p,q)=1$ and $q>0$, then
$h(\eta)=\log\max\\{|p|,q\\}$.

With these established notations, Matveev (see [9] or [4, Theorem 9.4]),
proved the ensuing result.

######  Theorem 2.3.

Let $\mathbb{K}$ be a number field of degree $D$ over $\mathbb{Q}$,
$\gamma_{1},\gamma_{2},\dots,\gamma_{t}$ be positive real numbers of
$\mathbb{K}$, and $b_{1},...,b_{t}$ rational integers. Put

$B\geq max\\{|b_{1}|,|b_{2}|,\dots,|b_{t}|\\},$

and

$\Lambda:=\gamma_{1}^{b_{1}}\dots\gamma_{t}^{b_{t}}-1.$

Let $A_{1},...,A_{t}$ be real numbers such that

$A_{i}\geq max\\{Dh(\gamma_{i}),|\log\gamma_{i}|,0.16\\},\quad i=1,\dots,t.$

Then, assuming that $\Lambda\neq 0$, we have

$|\Lambda|>exp\left(-1.4\times 30^{t+3}\times t^{4.5}\times D^{2}(1+\log
D)(1+logB)A_{1}\cdots A_{t}\right).$

The following criterion of Legendre, a well-known result from the theory of
Diophantine approximation, is used to reduce the upper bounds on variables
that are too large.

######  Lemma 2.4.

Let $\tau$ be an irrational number,
$\frac{p_{0}}{q_{0}},\frac{p_{1}}{q_{1}},\frac{p_{2}}{q_{2}},\dots$ be all the
convergents of the continued fraction of $\tau$, and $M$ be a positive
integer. Let $N$ be a non-negative integer such that $q_{N}>M$. Then putting
$a(M):=\max\\{a_{i}:i=0,1,2,\dots,N\\}$, the inequality

$\Bigm{|}\tau-\frac{r}{s}\Bigm{|}>\frac{1}{(a(M)+2)s^{2}},$

holds for all pairs $(r,s)$ of positive integers with $0<s<M$.

Another result which will play an important role in our proof is due to
Dujella and Pethö [7, Lemma 5 (a)].

######  Lemma 2.5.

Let $M$ be a positive integer, let $p/q$ be a convergent of the continued
fraction of the irrational $\gamma$ such that $q>6M$, and let $A,B,\mu$ be
some real numbers with $A>0$ and $B>1$. Let $\epsilon:=||\mu q||-M||\gamma
q||$, where $||\cdot||$ denotes the distance from the nearest integer. If
$\epsilon>0$, then there exists no solution to the inequality

$0<|u\gamma-v+\mu|<AB^{-w},$

in positive integers $u$, $v$ and $w$ with

$u\leq M\quad\text{and}\quad w\geq\frac{\log(Aq/\epsilon)}{\log B}.$

######  Lemma 2.6.

For any non-zero real number $x$, we have the following

(a) $0<x<e^{x}-1$.

(b) If $x<0$ and $|e^{x}-1|<1/2$, then $|x|<2|e^{x}-1|$.

##  3 Proof of Theorem 1.1

###  3.1 Upper bound for $n$.

First of all observe that, $(n,m,l,a)$ is a solution of (1.1) for $m=l=1$ and
$m=l=2$. So from now on, we assume that $m\geq 2$ and $l\geq 2$. Since
$F_{n}+F_{n-1}=F_{n+1}$, we can assume that $n>m+1$ and $n>l+1$. In particular
$n-m\geq 2$ and $n-l\geq 2$.

######  Proposition 3.1.

There are exactly $214$ solutions $(n,m,l,a)\in\mathbb{N}^{4}$ to (1.1) with
$n\leq 550$. All solutions satisfy $n\leq 42$ and $a\leq 28$. The list of
solutions is given in the appendix.

###### Proof.

The solutions were found by a brute force search with Mathematica. ∎

Due to Proposition 3.1, for the rest of paper, we assume that $n>550$. Using
Lemma 2.2 for $n>m>l>1$, we have

$0.38\alpha^{n}<F_{n}<F_{n}+F_{m}+F_{l}<0.48\alpha^{n}+0.48\alpha^{n-1}+0.48\alpha^{n-2}<0.97\alpha^{n}.$
(3.5)

Let us now get a relationship between $n$ and $a$. Without affecting
generality, we may now assume $a\geq 2$. Combining (1.1) with the right
inequality of (3.5), we get that

$2^{a-1}\leq
2^{a}-2^{\frac{a}{2}}<F_{n}+F_{m}+F_{l}<0.97\alpha^{n}<\alpha^{n}.$ (3.6)

When we combine (1.1) with the left inequality of (3.5), we get

$0.38\alpha^{n}<F_{n}+F_{m}+F_{l}<2^{a}+2^{\frac{a}{2}}<2^{a+1}.$

Thus

$n\frac{\log\alpha}{\log 2}+\frac{\log 0.38}{\log
2}-1<a<n\frac{\log\alpha}{\log 2}+1,$ (3.7)

where $\frac{\log\alpha}{\log 2}=0.6942\dots$. Hence, we have $a<n$. Using
(2.2) in (1.1), we get

$\frac{\alpha^{n}+\alpha^{m}+\alpha^{l}}{\sqrt{5}}-F_{n}+F_{m}+F_{l}=\frac{\beta^{n}+\beta^{m}+\beta^{l}}{\sqrt{5}}.$
(3.8)

Taking the absolute values on both sides in (3.8), we obtain

$\left|\frac{\alpha^{n}(1+\alpha^{m-n}+\alpha^{l-n})}{\sqrt{5}}-2^{a}\right|<\frac{|\beta|^{n}+|\beta|^{m}+|\beta|^{l}}{\sqrt{5}}+2^{\frac{a}{2}}<\frac{3}{5}+2^{\frac{a}{2}}<2^{\frac{a}{2}+1},$

for all $n\geq 4$ and $m,l\geq 2$.
Dividing both sides of the above relation by $2^{a}$, we get

$\left|\frac{\alpha^{n}(1+\alpha^{m-n}+\alpha^{l-n})}{2^{a}\sqrt{5}}-1\right|<2^{-\frac{a}{2}+1}.$
(3.9)

For the left-hand side, we apply Theorem 2.3 with the following data. Set

$\Lambda_{1}:=2^{-a}\alpha^{n}\frac{(1+\alpha^{m-n}+\alpha^{l-n})}{\sqrt{5}}-1.$

Note that $\Lambda_{1}\neq 0$. If $\Lambda_{1}=0$, we would get the relation

$2^{a}\sqrt{5}=\alpha^{n}+\alpha^{m}+\alpha^{l}.$ (3.10)

Conjugating the above relation in $\mathbb{Q}(\sqrt{5})$, we get

$-2^{a}\sqrt{5}=\beta^{n}+\beta^{m}+\beta^{l}.$ (3.11)

Equations (3.10) and (3.11) lead to

$\alpha^{n}<\alpha^{n}+\alpha^{m}+\alpha^{l}=|\beta^{n}+\beta^{m}+\beta^{l}|<|\beta^{n}|+|\beta^{m}|+|\beta^{l}|<1,$

which is impossible for positive integer $n$. Hence $\Lambda_{1}\neq 0$. We
take $t:=3$,

$\gamma_{1}:=2,\quad\gamma_{2}:=\alpha,\quad\gamma_{3}:=\frac{(1+\alpha^{m-n}+\alpha^{l-n})}{\sqrt{5}},$

and

$b_{1}:=-a,\quad b_{2}:=n,\quad b_{3}:=1.$

Note that $\mathbb{K}:=\mathbb{Q}(\sqrt{5})$ contains
$\gamma_{1},\gamma_{2},\gamma_{3}$ and has $D:=[\mathbb{K}:\mathbb{Q}]=2$.
Since $a<n$, we can take $B:=\max\\{|B_{1}|,|B_{2}|,|B_{3}|\\}=n$. Since
$h(\gamma_{1})=\log 2$ and $h(\gamma_{2})=\frac{\log\alpha}{2}$, it follows
that we can take $A_{1}:=1.4>Dh(\gamma_{1})$ and $A_{2}:=0.5>Dh(\gamma_{2})$.
To estimate $h(\gamma_{3})$, we begin by observing that

$\gamma_{3}=\frac{(1+\alpha^{m-n}+\alpha^{l-n})}{\sqrt{5}}<\frac{3}{\sqrt{5}}\quad\text{and}\quad\gamma_{3}^{-1}=\frac{\sqrt{5}}{(1+\alpha^{m-n}+\alpha^{l-n})}<\sqrt{5},$

it follows that $|\log\gamma_{3}|<1$. Using the properties of logarithmic
height, we get

$\displaystyle h(\gamma_{3})$
$\displaystyle\leq\log\sqrt{5}+|m-n|\left(\frac{\log\alpha}{2}\right)+|l-n|\left(\frac{\log\alpha}{2}\right)+2\log
2$
$\displaystyle=\log(4\sqrt{5})+(n-m)\left(\frac{\log\alpha}{2}\right)+(n-l)\left(\frac{\log\alpha}{2}\right).$

Thus, we can take
$A_{3}:=5+(2n-m-l)\log\alpha>\max\\{2h(\gamma_{3}),|\log\gamma_{3}|,0.16\\}$.
Then by using Theorem 2.3, we have

$\displaystyle\log|\Lambda_{1}|$ $\displaystyle>-1.4\times 30^{6}\times
3^{4.5}\times 2^{2}\times(1+\log 2)\times(1+\log n)\times 1.4\times
0.5\times(5+(2n-m-l)\log\alpha)$ $\displaystyle>-1.4\times 10^{12}\times\log
n\times(5+(2n-m-l)\log\alpha),$

where $1+\log n<2\log n$ holds for all $n\geq 3$. By comparing the above
inequality with (3.9), we get

$\left(\frac{a}{2}-1\right)\log 2<1.4\times 10^{12}\times\log
n\times(5+(2n-m-l)\log\alpha).$ (3.12)

Now consider the second linear form in logarithms by rewriting (3.8)
differently. Using (2.2), we get that

$\frac{\alpha^{n}}{\sqrt{5}}-(F_{n}+F_{m}+F_{l})=\frac{\beta^{n}}{\sqrt{5}}-F_{m}+F_{l}.$

Again, combining the above relation with (1.1), we get that

$\left|\frac{\alpha^{n}}{\sqrt{5}}-2^{a}\right|<2^{\frac{a}{2}}+\frac{|\beta|^{n}}{\sqrt{5}}+F_{m}+F_{l},$
(3.13)

where $|\beta|^{n}<1/2$ for all $n\geq 2$. Dividing both sides of inequality
(3.13) by $\frac{\alpha^{n}}{\sqrt{5}}$ and taking account that
$\alpha>\sqrt{2}$ and $n>m>l$, we obtain

$\displaystyle|1-2^{a}\alpha^{-n}\sqrt{5}|$
$\displaystyle<\frac{2^{\frac{a}{2}}\sqrt{5}}{\alpha^{n}}+\frac{\sqrt{5}}{2\alpha^{n}}+\frac{\sqrt{5}}{\alpha^{n-m}}+\frac{\sqrt{5}}{\alpha^{n-l}}$
$\displaystyle<2\sqrt{5}\max\\{\alpha^{m-n},\alpha^{l-n},\alpha^{a-n}\\}.$
(3.14)

For the left-hand side, we apply Theorem 2.3 with the following data. Set

$\Lambda_{2}:=2^{a}\alpha^{-n}\sqrt{5}-1.$

Note that $\Lambda_{2}\neq 0$. If $\Lambda_{2}=0$, then we have that
$2^{a}=\frac{\alpha^{n}}{\sqrt{5}}$. So $\alpha^{2n}\in\mathbb{Z}$, which is
not possible. Therefore $\Lambda_{2}\neq 0$. We take $t:=3$,

$\gamma_{1}:=2,\quad\gamma_{2}:=\alpha,\quad\gamma_{3}:=\sqrt{5},$

and

$b_{1}:=a,\quad b_{2}:=-n,\quad b_{3}:=1.$

Note that $\mathbb{K}:=\mathbb{Q}(\sqrt{5})$ contains
$\gamma_{1},\gamma_{2},\gamma_{3}$ and has $D:=2$. Since $a<n$, we deduce that
$B:=\max\\{|b_{1}|,|b_{2}|,|b_{3}|\\}=n.$ The logarithmic heights for
$\gamma_{1}$, $\gamma_{2}$ and $\gamma_{3}$ are calculated as follows:

$h(\gamma_{1})=\log 2,\quad
h(\gamma_{2})=\frac{\log\alpha}{2}\quad\text{and}\quad
h(\gamma_{3})=\log\sqrt{5}.$

Thus, we can take

$A_{1}:=1.4,\quad A_{2}:=0.5\quad\text{and}\quad A_{3}:=1.7.$

As before, applying Theorem 2.3, we have

$\displaystyle\log|\Lambda_{2}|$ $\displaystyle>-1.4\times 30^{6}\times
3^{4.5}\times 2^{2}\times(1+\log 2)\times(1+\log n)\times 1.4\times 0.5\times
1.7$ $\displaystyle>-2.31\times 10^{12}\times\log n.$

Comparing the above inequality with (3.1) implies that

$\min\\{(n-m)\log\alpha,(n-l)\log\alpha,(n-a)\log\alpha\\}<2.4\times
10^{12}\log n.$

Now the argument splits into three cases.
Case 1.
$\min\\{(n-m)\log\alpha,(n-l)\log\alpha,(n-a)\log\alpha\\}=(n-m)\log\alpha$.
In this case, we have

$(n-m)\log\alpha<2.4\times 10^{12}\log n.$

Case 2.
$\min\\{(n-m)\log\alpha,(n-l)\log\alpha,(n-a)\log\alpha\\}=(n-l)\log\alpha.$
In this case, we have

$(n-l)\log\alpha<2.4\times 10^{12}\log n.$

Now by using Case 1 and Case 2 in (3.12), we have

$\displaystyle\left(\frac{a}{2}-1\right)\log 2$ $\displaystyle<1.4\times
10^{12}\times\log n\times(5+(2n-m-l)\log\alpha)$ $\displaystyle<1.4\times
10^{12}\times\log n\times(5+4.8\times 10^{12}\log n)$
$\displaystyle<6.86\times 10^{24}\log^{2}n.$

Combining it with the left inequality of (3.7) and by a calculation in
Mathematica, we get

$a<9\times 10^{28}\quad\text{and}\quad n<1.87\times 10^{29}.$

Case 3.
$\min\\{(n-m)\log\alpha,(n-l)\log\alpha,(n-a)\log\alpha\\}=(n-a)\log\alpha.$
In this case, we have

$(n-a)\log\alpha<2.4\times 10^{12}\log n.$ (3.15)

Note that the right inequality of (3.7) yields

$n-a>n\left(1-\frac{\log\alpha}{\log 2}\right)-1.$ (3.16)

From (3.15) and (3.16), we obtain

$a<3.84\times 10^{14}\quad\text{and}\quad n<5.53\times 10^{14}.$

Thus, in all the three cases, we have

$a<9\times 10^{28}\quad\text{and}\quad n<1.87\times 10^{29}.$

Now we need to reduce the bound of $n$.

###  3.2 Reducing the bound on $n$

Let us assume that $n-m>550$, $n-l>550$ and $n-a>550$. In order to apply Lemma
2.5, we put

$z_{1}:=a\log 2-n\log\alpha+\log\sqrt{5}.$

Since we have assumed that $\min\\{n-m,n-l,n-a\\}>550$ therefore
$|e^{z_{1}}-1|<1/2$. Thus, by Lemma 2.6, we have that
$|z_{1}|<2|e^{z_{1}}-1|$. Replacing $z_{1}$ by its formula and by (3.1), we
get that

$|a\log
2-n\log\alpha+\log\sqrt{5}|<\frac{4\sqrt{5}}{\alpha^{\min\\{n-m,n-l,n-a\\}}}.$

Dividing both sides of the above inequality by $\log\alpha$, we can conclude

$0<\left|a\frac{\log
2}{\log\alpha}-n+\frac{\log\sqrt{5}}{\log\alpha}\right|<\frac{4\sqrt{5}}{\log\alpha}\cdot\alpha^{-w},$
(3.17)

where $w=\min\\{n-m,n-l,n-a\\}$. Putting now

$\gamma:=\frac{\log
2}{\log\alpha},\quad\mu:=\frac{\log\sqrt{5}}{\log\alpha},\quad
A:=\frac{4\sqrt{5}}{\log\alpha}\quad\text{and}\quad B:=\alpha,$

the above inequality (3.17) implies

$0<|a\gamma-n+\mu|<AB^{-w}.$

It is clear that $\gamma$ is an irrational number. We also put $M:=9\times
10^{28}$ which is an upper bound for $a$. We find that the convergent
$\frac{p}{q}=\frac{p_{69}}{q_{69}}$ is such that $q=q_{69}>6M$. By using this
we have that $\epsilon=||\mu q_{69}||-M||\gamma q_{69}||>0$. Therefore

$n-m<\frac{\log\left(\frac{4\sqrt{5}q}{\epsilon\
\log\alpha}\right)}{\log\alpha}<157\quad\text{or}\quad
n-l<\frac{\log\left(\frac{4\sqrt{5}q}{\epsilon\
\log\alpha}\right)}{\log\alpha}<157\quad\text{or}\quad
n-a<\frac{\log\left(\frac{4\sqrt{5}q}{\epsilon\
\log\alpha}\right)}{\log\alpha}<157.$

Thus, we have that either $n-m<157$ or $n<517$. The latter case contradicts
our assumption that $n>550$. Inserting the upper bound for $n-m$ into (3.12),
we get that $a<3.93\times 10^{15}$.

Finally, we shall use inequality (3.9) to improve the bound on $n$, where we
put

$z_{2}:=a\log 2-n\log\alpha+\log\psi(n-m,n-l),$

where $\psi(t,s):=\sqrt{5}(1+\alpha^{-t}+\alpha^{-s})^{-1}$. Therefore, (3.9)
implies that

$|1-e^{z_{2}}|<\frac{2}{2^{a/2}}.$ (3.18)

Note that $z_{2}\neq 0$. Thus, we distinguish the following cases. If
$z_{2}>0$, then

$0<z_{2}<e^{z_{2}}-1\leq\frac{2}{2^{a/2}}.$

Suppose $z_{2}<0$. Then, from (3.18) we have that $|1-e^{z_{2}}|<1/2$ and
therefore $e^{|z_{2}|}<2.$ Since $z_{2}<0$, we have

$0<|z_{2}|\leq e^{|z_{2}|}-1=e^{|z_{2}|}|e^{z_{2}}-1|<\frac{4}{2^{a/2}}.$

In any case, we have that

$0<|z_{2}|\leq\frac{4}{2^{a/2}}.$

Replacing $z_{2}$ in the above inequality and dividing both sides by
$\log\alpha$, we conclude that

$0<\left|a\frac{\log
2}{\log\alpha}-n+\frac{\log\psi(n-m,n-l)}{\log\alpha}\right|<\frac{4}{\log\alpha}\cdot
2^{-a/2}.$ (3.19)

Here we take $M:=3.93\times 10^{15}$ which is an upper bound for $a$. By
virtue of Lemma 2.5 with

$\gamma:=\frac{\log
2}{\log\alpha},\quad\mu:=\frac{\log\psi(n-m,n-l)}{\log\alpha},\quad
A:=\frac{4}{\log\alpha},B:=\sqrt{2}$

to inequality (3.19) for all choices $n-m\in\\{1,\dots,157\\}$ and
$n-l\in\\{1,\dots,157\\}$ except when

$(n-m,n-l)\in\\{(0,3),(1,1),(1,5),(3,0),(3,4),(4,3),(5,1),(7,8),(8,7)\\}.$

Furthermore, with the help of Mathematica, we find that if $(n,m,l,a)$ is a
possible solution of (1.1), excluding those cases mentioned previously, then
$n\leq 550$. This is false because we have assumed $n>550$.
Finally, we deal with the special cases where

$(n-m,n-l)\in\\{(0,3),(1,1),(1,5),(3,0),(3,4),(4,3),(5,1),(7,8),(8,7)\\}.$

As we are seeing, there are certain symmetric cases and there is no
significant difference in their study. Therefore, we deal with the cases
$(n-m,n-l)\in\\{(1,1),(3,0),(4,3),(5,1),(8,7)\\}.$ For these special cases, we
have that

$\frac{\log\psi(t,s)}{\log\alpha}=\begin{cases}0&\text{if}~{}~{}(t,s)=(1,1);\\\
0&\text{if}~{}~{}(t,s)=(3,0);\\\ 1&\text{if}~{}~{}(t,s)=(4,3);\\\ 2-\frac{\log
2}{\log\alpha}&\text{if}~{}~{}(t,s)=(5,1);\\\ 3-\frac{\log
2}{\log\alpha}&\text{if}~{}~{}(t,s)=(8,7).\end{cases}$

Note that, when we apply Lemma 2.5 to the inequality (3.19), the parameters
$\gamma$ and $\mu$ appearing in Lemma 2.5 are linearly dependent and the
corresponding value of $\epsilon$ is always negative. When $n-m=1$ and $n-l=1$
from (3.19), we have

$0<|a\gamma-n|<\frac{4}{\log\alpha}\cdot
2^{-\frac{1}{2}(n\frac{\log\alpha}{\log 2}+\frac{\log 0.38}{\log 2}-1)}.$
(3.20)

We recall that $a<3.93\times 10^{15}$. Let
$[a_{0},a_{1},a_{2},a_{3},a_{4},\dots]=[1,2,3,1,2,\dots]$ be the continued
fraction of $\gamma$. A quick search using Mathematica reveals that

$q_{34}<3.93\times 10^{15}<q_{35}.$

Furthermore, $a_{M}:=\max\\{a_{i}:i=0,1,\dots,35\\}=a_{17}=134$. So by Lemma
2.4 we have

$|a\gamma-n|>\frac{1}{(a_{M}+2)a}.$

Comparing estimates (3.19) and (3.20) we get that $n<112$. Similarly, one can
get $n<112$ in all other cases. This again contradicts our assumption that
$n>550$. Hence, the result is proved.

##  4 Appendix

The solutions of Diophantine inequality (1.1) are given in Table 1.

Table 1 $|F_{2}+F_{2}+F_{2}-2|<\sqrt{2}$ |  $|F_{11}+F_{8}+F_{i}-2^{7}|<2^{7/2},i=6,7,8$
---|---
$|F_{2}+F_{2}+F_{2}-2^{2}|<2$ |  $|F_{11}+F_{9}+F_{i}-2^{7}|<2^{7/2}$, $i=2,3,\dots,8$
$|F_{3}+F_{i}+F_{2}-2^{2}|<2,i=2,3$ |  $|F_{11}+F_{11}+F_{11}-2^{8}|<2^{4}$
$|F_{4}+F_{2}+F_{2}-2^{2}|<2$ |  $|F_{12}+F_{10}+F_{10}-2^{8}|<2^{4}$
$|F_{4}+F_{3}+F_{i}-2^{3}|<2^{3/2},i=2,3$ |  $|F_{12}+F_{11}+F_{i}-2^{8}|<2^{4},i=6,7,8,9$
$|F_{4}+F_{4}+F_{i}-2^{3}|<2^{3/2},i=2,3,4$ |  $|F_{13}+F_{5}+F_{i}-2^{8}|<2^{4},i=4,5$
$|F_{5}+F_{2}+F_{2}-2^{3}|<2^{3/2}$ |  $|F_{13}+F_{6}+F_{i}-2^{8}|<2^{4},i=2,3,4,5,6$
$|F_{5}+F_{3}+F_{i}-2^{3}|<2^{3/2},i=2,3$ |  $|F_{13}+F_{7}+F_{i}-2^{8}|<2^{4},i=2,3,4,5,6,7$
$|F_{5}+F_{4}+F_{i}-2^{3}|<2^{3/2},i=2,3$ |  $|F_{13}+F_{8}+F_{i}-2^{8}|<2^{4},i=2,3,4,5,6,7$
$|F_{6}+F_{2}+F_{2}-2^{3}|<2^{3/2}$ |  $|F_{13}+F_{9}+F_{i}-2^{8}|<2^{4},i=2,3,4$
$|F_{5}+F_{5}+F_{i}-2^{4}|<2^{2},i=4,5$ |  $|F_{13}+F_{12}+F_{12}-2^{9}|<2^{9/2}$
$|F_{6}+F_{4}+F_{i}-2^{4}|<2^{2},i=3,4$ |  $|F_{13}+F_{13}+F_{i}-2^{9}|<2^{9/2},i=9,10$
$|F_{6}+F_{5}+F_{i}-2^{4}|<2^{2},i=2,3,4,5$ |  $|F_{14}+F_{11}+F_{i}-2^{9}|<2^{9/2},i=9,10$
$|F_{6}+F_{6}+F_{i}-2^{4}|<2^{2},i=2,3,4$ |  $|F_{14}+F_{12}+F_{i}-2^{9}|<2^{9/2},i=2,3,4,5,6,7$
$|F_{7}+F_{2}+F_{2}-2^{4}|<2^{2}$ |  $|F_{15}+F_{14}+F_{i}-2^{10}|<2^{5},i=6,7,8,9,10$
$|F_{7}+F_{3}+F_{i}-2^{4}|<2^{2},i=2,3$ |  $|F_{16}+F_{4}+F_{4}-2^{10}|<2^{5}$
$|F_{7}+F_{4}+F_{i}-2^{4}|<2^{2},i=2,3,4$ |  $|F_{16}+F_{5}+F_{i}-2^{10}|<2^{5},i=2,3,4,5$
$|F_{7}+F_{5}+F_{2}-2^{4}|<2^{2}$ |  $|F_{16}+F_{6}+F_{i}-2^{10}|<2^{5},i=2,3,4,5,6$
$|F_{7}+F_{6}+F_{6}-2^{5}|<2^{5/2}$ |  $|F_{16}+F_{7}+F_{i}-2^{10}|<2^{5},i=2,3,4,5,6,7$
$|F_{7}+F_{7}+F_{i}-2^{5}|<2^{5/2},i=2,3,4,5,6$ |  $|F_{16}+F_{8}+F_{i}-2^{10}|<2^{5},i=2,3,4,\dots,8$
$|F_{8}+F_{4}+F_{4}-2^{5}|<2^{5/2}$ |  $|F_{16}+F_{9}+F_{i}-2^{10}|<2^{5},i=2,3,\dots,9$
$|F_{8}+F_{5}+F_{i}-2^{5}|<2^{5/2},i=2,3,4,5$ |  $|F_{16}+F_{10}+F_{i}-2^{10}|<2^{5},i=2,3,4,5,6,7$
$|F_{8}+F_{6}+F_{i}-2^{5}|<2^{5/2},i=2,3,4,5,6$ |  $|F_{16}+F_{16}+F_{i}-2^{10}|<2^{5},i=9,10,11,13$
$|F_{8}+F_{7}+F_{i}-2^{5}|<2^{5/2},i=2,3,4$ |  $|F_{17}+F_{13}+F_{13}-2^{11}|<2^{11/2}$
$|F_{9}+F_{i}+F_{2}-2^{5}|<2^{5/2},i=2,3$ |  $|F_{17}+F_{14}+F_{i}-2^{11}|<2^{11/2},i=9,10,11$
$|F_{8}+F_{8}+F_{8}-2^{6}|<2^{3}$ |  $|F_{20}+F_{16}+F_{14}-2^{13}|<2^{13/2}$
$|F_{9}+F_{7}+F_{7}-2^{6}|<2^{3}$ |  $|F_{21}+F_{21}+F_{21}-2^{15}|<2^{15/2}$
$|F_{9}+F_{8}+F_{i}-2^{6}|<2^{3},i=3,4,5,6,7$ |  $|F_{22}+F_{21}+F_{19}-2^{15}|<2^{15/2}$
$|F_{9}+F_{9}+F_{i}-2^{6}|<2^{3},i=2,3,4$ |  $|F_{23}+F_{18}+F_{17}-2^{15}|<2^{15/2}$
$|F_{10}+F_{i}+F_{2}-2^{6}|<2^{3},i=2,3$ |  $|F_{23}+F_{19}+F_{i}-2^{15}|<2^{15/2}$,$i=2,3,\dots,11$
$|F_{10}+F_{i}+F_{3}-2^{6}|<2^{3},i=3,4$ |  $|F_{24}+F_{22}+F_{17}-2^{16}|<2^{8}$
$|F_{10}+F_{4}+F_{i}-2^{6}|<2^{3},i=2,4$ |  $|F_{26}+F_{20}+F_{18}-2^{17}|<2^{17/2}$
$|F_{10}+F_{5}+F_{i}-2^{6}|<2^{3},i=2,3,4,5$ |  $|F_{29}+F_{20}+F_{18}-2^{19}|<2^{19/2}$
$|F_{10}+F_{6}+F_{i}-2^{6}|<2^{3},i=2,3,4,5,6$ |  $|F_{41}+F_{40}+F_{29}-2^{28}|<2^{14}$
$|F_{10}+F_{7}+F_{i}-2^{6}|<2^{3},i=2,3,4$ |  $|F_{42}+F_{28}+F_{27}-2^{28}|<2^{14}$
$|F_{10}+F_{9}+F_{9}-2^{7}|<2^{7/2}$ |  $|F_{42}+F_{29}+F_{i}-2^{28}|<2^{14}$,$i=2,3,\dots,22$
$|F_{10}+F_{10}+F_{i}-2^{7}|<2^{7/2},i=6,7,8$ |

## References

  * [1] E. F. Bravo and J. J. Bravo, Powers of 2 as sums of three Fibonacci numbers, Lith. Math. J. 55(3) (2015), 301–311.
  * [2] J. J. Bravo, C. A. Gomez and J. L. Herrera, $k$-Fibonacci numbers close to a power of 2, Quaest. Math. 44(12) (2021), 1681–1690.
  * [3] J. J. Bravo and F. Luca, On the Diophantine equation $F_{n}+F_{m}=2^{a}$, Quaest. Math. 39(3) (2016), 391–400.
  * [4] Y. Bugeaud, M. Mignotte and S. Siksek, Classical and modular approaches to exponential Diophantine equations. I. Fibonacci and Lucas perfect powers, Ann. of Math. 163(3) (2006), 969–1018.
  * [5] S. Chern and A. Cui, Fibonacci numbers close to a power of 2, Fibonacci Quart. 52(4) (2014), 344–348.
  * [6] K. C. Chim and V. Ziegler, On Diophantine equations involving sums of Fibonacci numbers and powers of 2, Integers. 18(30) (2018), A99, 30pp.
  * [7] A. Dujella and A. Pethő, A generalization of a theorem of Baker and Davenport, Quart. J. Math. Oxford Ser. 49(195) (1998), 291–306.
  * [8] E. Hasanalizade, Sums of Fibonacci numbers close to a power of 2, (2022), arXiv:2206.10356.
  * [9] E. M. Matveev, An explicit lower bound for a homogeneous rational linear form in the logarithms of algebraic numbers, II, Izv. Ross. Akad. Nauk Ser. Mat., 64(6) (2000), 125–180.
  * [10] V. Ziegler, Sums of Fibonacci numbers that are perfect powers, (2021), arXiv:2103.14487.

Department of Mathematics,
KIIT University, Bhubaneswar,
Bhubaneswar 751024, Odisha, India.
Email<EMAIL_ADDRESS>
Department of Mathematics,
KIIT University, Bhubaneswar,
Bhubaneswar 751024, Odisha, India.
Email<EMAIL_ADDRESS>


# Synergistic Signal Denoising for Multimodal Time Series of Structure
Vibration

Yang Yu  Han Chen

###### Abstract

Structural Health Monitoring (SHM) plays an indispensable role in ensuring the
longevity and safety of infrastructure. With the rapid growth of sensor
technology, the volume of data generated from various structures has seen an
unprecedented surge, bringing forth challenges in efficient analysis and
interpretation. This paper introduces a novel deep learning algorithm tailored
for the complexities inherent in multimodal vibration signals prevalent in
SHM. By amalgamating convolutional and recurrent architectures, the algorithm
adeptly captures both localized and prolonged structural behaviors. The
pivotal integration of attention mechanisms further enhances the model’s
capability, allowing it to discern and prioritize salient structural responses
from extraneous noise. Our results showcase significant improvements in
predictive accuracy, early damage detection, and adaptability across multiple
SHM scenarios. In light of the critical nature of SHM, the proposed approach
not only offers a robust analytical tool but also paves the way for more
transparent and interpretable AI-driven SHM solutions. Future prospects
include real-time processing, integration with external environmental factors,
and a deeper emphasis on model interpretability.

##  1 Introduction

Structural health monitoring (SHM) has emerged as a vital field of research,
geared towards preserving the longevity and safety of civil infrastructure
[1]. A critical component of SHM is the analysis of vibration time series
data, which offers insights into the behavior, health, and performance of
structures [2]. As infrastructure, especially in urban regions, is subject to
a myriad of dynamic forces—ranging from wind to traffic loads - it becomes
pivotal to extract clear and meaningful data from the complex vibration
signatures that these forces induce. However, one of the significant
challenges plaguing SHM practitioners is the interference of noise in these
vibration signals, which can distort interpretations and lead to unreliable
conclusions.

The dynamic response of structures is often manifested as multimodal
vibrations, meaning multiple modes or patterns of vibration coexist. These
modes, each characterized by its frequency and shape, provide a fingerprint of
the structure’s health and dynamic properties. For example, a bridge might
have one mode that represents a swaying motion, while another might represent
a bouncing motion [3]. Distinguishing between these modes and their individual
characteristics is a nuanced task, made increasingly difficult when the
signals are clouded by noise.

Noise in the context of vibration time series can arise from multiple sources:
instrumental noise from sensors, environmental interferences, or other
unpredictable external factors [4, 5]. Traditional denoising methods, while
effective to a degree, sometimes struggle to cater to the unique intricacies
of multimodal structural vibrations [6]. Such challenges underscore the
necessity for novel denoising techniques tailored to the specific nature of
these signals [7]. Enter the realm of synergistic signal denoising—a cutting-
edge approach that promises a more holistic treatment of multimodal vibration
data [8]. Unlike conventional methods that treat each mode of vibration as a
separate entity, the synergistic approach recognizes the interconnectedness of
these modes [9]. By leveraging the relationships between different modal
responses, this methodology aims to achieve superior denoising outcomes,
preserving the integrity of the structural response while effectively
mitigating the noise. The importance of an effective denoising technique
cannot be understated [10, 11]. Clear, noise-free signals allow for accurate
modal identification, which in turn can be pivotal in detecting anomalies or
damage in structures. Furthermore, with the increasing integration of machine
learning and artificial intelligence in SHM, the quality of input data, i.e.,
the vibration signals—directly influences the quality of predictions and
insights [12].

This paper ventures into the depths of synergistic signal denoising for
multimodal structure vibration time series. We explore its foundational
principles, contrast its performance with traditional techniques, and
demonstrate its potential in real-world applications. The ultimate objective
is to illuminate the path forward in harnessing clean and precise vibration
data, a cornerstone in the edifice of modern structural health monitoring.

##  2 Prior Arts and Methods

Vibration-based Structural Health Monitoring (SHM) has been extensively
studied over the past few decades, primarily due to its potential to detect
anomalies and ensure structural safety. Central to SHM’s efficacy is the
clarity of the acquired vibration signals. However, these signals,
particularly from civil structures, are often contaminated with noise from
various sources, complicating their analysis [13, 9]. Traditional denoising
methods like wavelet decomposition and empirical mode decomposition have been
applied to structure-borne signals with varying degrees of success. While
these techniques can reduce noise, they sometimes struggle to retain the
integrity of multimodal signals, leading to potential information loss [14].

A shift in the paradigm is noted in recent years towards more integrated
denoising techniques. A synergistic approach, treating multimodal vibrations
as interrelated rather than isolated events, suggest an improved retention of
signal characteristics post-denoising compared to conventional methods.
Despite its promise, synergistic denoising in the context of SHM remains an
underexplored territory, warranting further research to validate its full
potential and broader applicability [15]. Structural responses to dynamic
loads often manifest in the form of vibration signals that are complex and
multifaceted. Analyzing these signals is essential to understand the
underlying structural behavior, particularly when it is governed by multiple
modes of vibration. In this section, we delve into the fundamentals of
multimodal vibrations and discuss common noise sources and characteristics
inherent to these signals.

Every structure possesses a set of natural frequencies at which it tends to
vibrate when subjected to external stimuli. These frequencies correspond to
different modes of vibration, with each mode showcasing a distinct deformation
shape or pattern. This phenomenon is widely recognized in the context of
structural dynamics as ’modal analysis’. For instance, a simple beam might
primarily deflect in a singular arc-like shape at its fundamental frequency, a
mode often termed as the first mode. However, when vibrated at higher
frequencies, the same beam might exhibit complex deformation patterns like
forming two or more arcs—these represent the higher modes of vibration.

In real-world scenarios, structures don’t vibrate in a singular mode. Instead,
they display a combination of these modes, giving rise to what is termed as
’multimodal vibrations’. The presence of multiple vibration modes, especially
in large and complex structures, makes signal analysis more intricate. Each
mode carries unique information about the structure’s health and behavior.
Disentangling these concurrent modes from vibration time series and accurately
identifying their individual frequencies and shapes is of paramount importance
in structural health monitoring. Multimodal vibration signals, while rich in
information, are often marred by noise, which can impede accurate signal
interpretation [16, 17]. The noise sources contaminating these signals can be
broadly categorized as [18]:

Instrumental Noise: Every measurement tool, be it an accelerometer, a strain
gauge, or any other sensor, has an inherent level of noise. Factors like the
sensor’s sensitivity, its electronic components, and even its placement can
introduce unwanted signals, often drowning the subtle nuances of certain
vibration modes.

Environmental Interferences: Structures are not isolated entities. They
interact with their surroundings, which means their recorded vibration signals
are a cumulative effect of their intrinsic response and the myriad of
environmental factors. Wind-induced vibrations, seismic activities, or even
nearby vehicular traffic can introduce extraneous noise components into the
signals (Doebling et al., 1996).

Operational Loads: For operational structures, like bridges or buildings, the
regular loads they bear—whether from vehicles, occupants, or machinery—can
mask or distort their natural vibration patterns, complicating the analysis.

The characteristics of noise in multimodal signals are often unpredictable.
They might manifest as random spikes in the time domain, irregular
fluctuations in amplitude, or unexpected frequency components in the frequency
domain. Separating this noise, especially without compromising the integrity
of the original vibration modes, remains a significant challenge in the realm
of structural health monitoring.

Figure 1: Mechanical motor shaft vibrational signals. Three noisy signals, all
originating from the same vibration source, are displayed.

##  3 Results and Discussion

Designing a deep learning algorithm tailored for multiple time series,
especially one that incorporates attention mechanisms, necessitates a thorough
approach that ensures the model can effectively capture underlying temporal
patterns and dependencies among series. When dealing with multiple time
series, the inherent relationships and dependencies between different
sequences can be intricate. Traditional analytical methods sometimes fall
short in accurately deciphering these complexities. This is where modern deep
learning, fortified by attention mechanisms, offers promising avenues.

Before diving into the architecture, it’s worth emphasizing the importance of
preprocessing. Like with most data-centric models, the quality and structure
of the input can profoundly impact the performance. Time series data, given
its temporal nature, often needs to be normalized, usually between a scale of
0 and 1, using methods like Min-Max scaling. This step ensures that all time
series are on a consistent scale, an essential aspect for neural networks to
function optimally. After normalization, the time series data should be
divided into overlapping windows. This segmentation aids in feeding the
network with structured sequences of data. If, for example, a time series has
a length of T and you choose a window size of W, the model will predict the
subsequent point by taking into account the preceding W points.

With the data preprocessed, the architectural design of the deep learning
model becomes the next focus. The starting point is an input layer designed to
accommodate segments of the time series. Its size is dictated by the chosen
window size, W. While traditional time series methods, such as ARIMA or
Exponential Smoothing, operate on predefined statistical rules, deep learning
offers a more adaptive, data-driven approach. A series of convolutional layers
can follow the input layer, although this is optional. Convolutional Neural
Networks (CNNs) have proven their merit in detecting local patterns within a
sequence. By deploying one-dimensional convolutions, these layers can sieve
out salient features from the input series, providing a richer representation
for the subsequent layers.

Recurrent layers, particularly those using Long Short-Term Memory (LSTM) cells
or Gated Recurrent Units (GRU), are integral for models dealing with temporal
data. Their design inherently allows them to understand sequences, making them
apt for time series data. These layers are adept at capturing dependencies and
patterns over extended durations, crucial when dealing with multiple time
series. Therefore, incorporating two or more LSTM or GRU layers that
sequentially pass data can provide depth to the model.

The pivotal aspect of this architecture is the attention mechanism. When
dealing with multiple sequences, not all segments or points within a series
are equally significant. The attention mechanism facilitates the model’s
ability to focus on specific sequences, allowing it to allocate different
attention scores to various segments. These scores are calculated based on the
hidden states of the recurrent layers, determining the relevance of each time
step in the series. By employing a softmax function, these scores are then
normalized to fit between 0 and 1, which subsequently determine the weights of
each step in the sequence. The culmination of this process results in the
derivation of a context vector. This vector is essentially a weighted sum of
the recurrent layers’ hidden states and the attention scores, summarizing the
most influential segments of the time series.

Post attention layer, dense or fully connected layers can be layered into the
model. The context vector obtained from the attention layer, infused with the
distilled essence of the most pertinent parts of the series, is processed
through these dense layers, enabling the model to further learn from the data
[19]. The tail end of this architecture concludes with an output layer.
Depending on the task, if the aim is forecasting just one step ahead, a
singular neuron is used. If multiple steps are to be predicted, multiple
neurons corresponding to the number of steps are integrated.

The final stages involve training the model. It’s essential to choose an
appropriate loss function, with Mean Squared Error (MSE) being apt for
regression-based tasks. The choice of optimizer is also critical. The Adam
optimizer, with its adaptive learning capabilities, is often a favorable
choice for deep learning models. To counteract the pitfalls of overfitting,
early stopping can be employed. By consistently monitoring the validation loss
and halting training when there’s no discernible improvement, the model can be
kept in check. Lastly, the model’s performance should be assessed on unseen
data using metrics like the Root Mean Squared Error (RMSE) or the Mean
Absolute Error (MAE) to gauge the model’s predictive accuracy.

In conclusion, designing a deep learning model embedded with attention
mechanisms for multiple time series signals a paradigm shift in time series
analysis. The ability to prioritize certain segments over others allows for a
nuanced understanding of the data, leading to more robust and accurate
predictions.

The designed deep learning algorithm for multiple time series, fortified with
attention mechanisms, brings forth several significant advantages that cater
to the complexities inherent to temporal data. Firstly, by going beyond
traditional statistical methods, it offers a data-driven approach that is
highly adaptive, learning directly from the unique patterns and intricacies
embedded within the data. This adaptability ensures that the model isn’t
confined by pre-defined rules, allowing it to evolve with the nuances of
different datasets.

Secondly, the incorporation of attention mechanisms marks a revolutionary
stride in handling multiple sequences. In most time series, certain segments
or points are inherently more significant than others, and the attention
mechanism ensures the model recognizes and prioritizes these pivotal
sequences. By assigning different attention scores to different segments, the
model can focus on crucial sequences while simultaneously downplaying less
relevant segments. This focused approach not only enhances prediction accuracy
but also aids in interpreting which segments of the data the model deems most
influential, offering valuable insights into the underlying patterns of the
series.

Additionally, the blend of convolutional and recurrent layers ensures a
holistic understanding of the data. While convolutional layers effectively
detect local patterns, recurrent layers, especially LSTM or GRU cells, retain
memory of past sequences, capturing long-term dependencies, a crucial aspect
in time series analysis. This fusion guarantees that the model is well-
equipped to understand both immediate and historical patterns.

Moreover, the model’s modular architecture ensures flexibility. Depending on
the data’s nature, one can adjust the depth of convolutional or recurrent
layers or even refine the attention mechanism. This flexibility means that the
algorithm can be tailored to suit various datasets with varying complexities.
The training methodologies, such as the use of adaptive optimizers and early
stopping, further enhance the model’s robustness, ensuring it doesn’t overfit
and generalizes well on unseen data.

In essence, this designed algorithm represents a culmination of some of the
best practices in deep learning, tailored explicitly for multiple time series.
Its ability to learn from data adaptively, prioritize significant sequences,
and capture both local and global patterns, positions it as a formidable tool
in the realm of time series analysis.

##  4 Conclusion and Future Work

In the realm of Structural Health Monitoring (SHM), our proposed deep learning
algorithm, complemented by attention mechanisms, offers a leap forward in
efficiently analyzing and interpreting multimodal vibration signals. The
challenges presented by the vast and intricate temporal data inherent in SHM
demand solutions that can adapt, learn, and most importantly, prioritize
salient structural responses over extraneous noise [20]. This algorithm, with
its blend of convolutional and recurrent architectures, adeptly captures both
instantaneous and prolonged structural behaviors, crucial for early damage
detection and predictive maintenance.

Moving forward, there’s substantial potential in refining this algorithm for
SHM-specific challenges. Given the criticality of timely damage detection in
infrastructure, real-time processing and analysis can be explored further,
possibly integrating edge computing for immediate on-site assessments. Another
avenue lies in fusing SHM data with external environmental factors like
weather conditions, which often play a role in structural health. As the world
of SHM continually advances, there will also be an increasing need to ensure
our deep learning models are not just accurate but also transparent.
Interpretability, especially in safety-critical applications like SHM, is
paramount. Future work will need to delve deeper into making these complex
models more understandable, ensuring engineers and stakeholders can trust and
act upon the insights provided.

## References

  * [1] Jingzhou Xin et al.  “Bridge deformation prediction based on SHM data using improved VMD and conditional KDE”  In _Engineering Structures_ 261 Elsevier, 2022, pp. 114285
  * [2] GT Webb, Paul J Vardanega and Campbell R Middleton  “Categories of SHM deployments: technologies and capabilities”  In _Journal of Bridge Engineering_ 20.11 American Society of Civil Engineers, 2015, pp. 04014118
  * [3] C Scuro et al.  “Internet of Things (IoT) for masonry structural health monitoring (SHM): Overview and examples of innovative systems”  In _Construction and Building Materials_ 290 Elsevier, 2021, pp. 123092
  * [4] Yue Hu, Fucai Li, Hongguang Li and Chengliang Liu  “An enhanced empirical wavelet transform for noisy and non-stationary signal processing”  In _Digital signal processing_ 60 Elsevier, 2017, pp. 220–229
  * [5] Yue Hu et al.  “An adaptive and tacholess order analysis method based on enhanced empirical wavelet transform for fault detection of bearings with varying speeds”  In _Journal of Sound and Vibration_ 409 Elsevier, 2017, pp. 241–255
  * [6] Youzhi Liang, Wen Liang and Jianguo Jia  “Structural Vibration Signal Denoising Using Stacking Ensemble of Hybrid CNN-RNN”  In _arXiv e-prints_ , 2023, pp. arXiv–2303
  * [7] Youzhi Liang  “Structural Vibration Signal Denoising Using Stacking Ensemble of Hybrid CNN-RNN. Advances in Artificial Intelligence and Machine Learning. 2022; 3 (2): 65”, 2006
  * [8] Dongming Feng and Maria Q Feng  “Experimental validation of cost-effective vision-based structural health monitoring”  In _Mechanical Systems and Signal Processing_ 88 Elsevier, 2017, pp. 199–211
  * [9] Dongming Feng and Maria Q Feng  “Computer vision for SHM of civil infrastructure: From dynamic response measurement to damage detection–A review”  In _Engineering Structures_ 156 Elsevier, 2018, pp. 105–117
  * [10] N Dervilis, K Worden and EJ Cross  “On robust regression analysis as a means of exploring environmental and operational conditions for SHM data”  In _Journal of Sound and Vibration_ 347 Elsevier, 2015, pp. 279–296
  * [11] Chao Che et al.  “Constrained transformer network for ECG signal processing and arrhythmia classification”  In _BMC Medical Informatics and Decision Making_ 21.1 BioMed Central, 2021, pp. 1–13
  * [12] Lili Ge et al.  “Study on the scattered sound modulation with a programmable chessboard device”  In _AIP Advances_ 13.4 AIP Publishing, 2023
  * [13] Dongming Feng and Maria Q Feng  “Vision-based multipoint displacement measurement for structural health monitoring”  In _Structural Control and Health Monitoring_ 23.5 Wiley Online Library, 2016, pp. 876–890
  * [14] Hong-Nan Li, Dong-Sheng Li and Gang-Bing Song  “Recent applications of fiber optic sensors to health monitoring in civil engineering”  In _Engineering structures_ 26.11 Elsevier, 2004, pp. 1647–1657
  * [15] G Song et al.  “Concrete structural health monitoring using embedded piezoceramic transducers”  In _Smart Materials and Structures_ 16.4 IOP Publishing, 2007, pp. 959
  * [16] Cole S Hamey, Wahyu Lestari, Pizhong Qiao and Gangbing Song  “Experimental damage identification of carbon/epoxy composite beams using curvature mode shapes”  In _Structural Health Monitoring_ 3.4 Sage Publications Sage CA: Thousand Oaks, CA, 2004, pp. 333–353
  * [17] G Song, YL Mo, K Otero and H Gu  “Health monitoring and rehabilitation of a concrete structure using intelligent materials”  In _Smart materials and structures_ 15.2 IOP Publishing, 2006, pp. 309
  * [18] Enrique Garcia-Macias et al.  “Micromechanics modeling of the uniaxial strain-sensing property of carbon nanotube cement-matrix composites for SHM applications”  In _Composite Structures_ 163 Elsevier, 2017, pp. 195–215
  * [19] Youzhi Liang and Wen Liang  “ResWCAE: Biometric Pattern Image Denoising Using Residual Wavelet-Conditioned Autoencoder”  In _arXiv preprint arXiv:2307.12255_ , 2023
  * [20] Duc-Duy Ho et al.  “Solar-powered multi-scale sensor node on lmote 2 platform for hybrid SHM in cable-stayed bridge”  In _Smart Structures and Systems_ 9.2 Techno-Press, P. O. Box 33 Yusong Taejon 305-600 Korea, Republic of, 2012, pp. 145–164


11institutetext:  1Saarbrücken Graduate School of Computer Science, Saarland
University, Saarland Informatics Campus E1 3, Saarbrücken, Germany
🖂 11email<EMAIL_ADDRESS>
2Univeristy of Trieste, Trieste, Italy
3 Saarland University, Saarland Informatics Campus E1 3, Saarbrücken, Germany

# Analysis of Markov Jump Processes under Terminal Constraints

Michael Backenköhler,🖂 11 Luca Bortolussi 2233 Gerrit Großmann 11 Verena Wolf
1133

###### Abstract

Many probabilistic inference problems such as stochastic filtering or the
computation of rare event probabilities require model analysis under initial
and terminal constraints. We propose a solution to this _bridging problem_ for
the widely used class of population-structured Markov jump processes. The
method is based on a state-space lumping scheme that aggregates states in a
grid structure. The resulting approximate bridging distribution is used to
iteratively refine relevant and truncate irrelevant parts of the state-space.
This way, the algorithm learns a well-justified finite-state projection
yielding guaranteed lower bounds for the system behavior under endpoint
constraints. We demonstrate the method’s applicability to a wide range of
problems such as Bayesian inference and the analysis of rare events.

###### Keywords:

Bayesian Inference Bridging problem Smoothing Lumping Rare Events.

##  1 Introduction

Discrete-valued continuous-time Markov Jump Processes (MJP) are widely used to
model the time evolution of complex discrete phenomena in continuous time.
Such problems naturally occur in a wide range of areas such as chemistry [16],
systems biology [49, 46], epidemiology [36] as well as queuing systems [10]
and finance [39]. In many applications, an MJP describes the stochastic
interaction of populations of agents. The state variables are counts of
individual entities of different populations.

Many tasks, such as the analysis of rare events or the inference of agent
counts under partial observations naturally introduce terminal constraints on
the system. In these cases, the system’s initial state is known, as well as
the system’s (partial) state at a later time-point. The probabilities
corresponding to this so-called _bridging problem_ are often referred to as
_bridging probabilities_ [17, 19]. For instance, if the exact, full state of
the process $X_{t}$ has been observed at time $0$ and $T$, the bridging
distribution is given by

$\Pr(X_{t}=x\mid X_{0}=x_{0},X_{T}=x_{g})$

for all states $x$ and times $t\in[0,T]$. Often, the condition is more
complex, such that in addition to an initial distribution, a terminal
distribution is present. Such problems typically arise in a Bayesian setting,
where the a priori behavior of a system is filtered such that the posterior
behavior is compatible with noisy, partial observations [11, 25]. For example,
time-series data of protein levels is available while the mRNA concentration
is not [1, 25]. In such a scenario our method can be used to identify a good
truncation to analyze the probabilities of mRNA levels.

Bridging probabilities also appear in the context of rare events. Here, the
rare event is the terminal constraint because we are only interested in paths
containing the event. Typically researchers have to resort to Monte-carlo
simulations in combination with variance reduction techniques in such cases
[14, 26].

Efficient numerical approaches that are not based on sampling or ad-hoc
approximations have rarely been developed.

Here, we combine state-of-the-art truncation strategies based on a forward
analysis [28, 4] with a refinement approach that starts from an abstract MJP
with lumped states. We base this lumping on a grid-like partitioning of the
state-space. Throughout a lumped state, we assume a uniform distribution that
gives an efficient and convenient abstraction of the original MJP. Note that
the lumping does not follow the classical paradigm of Markov chain lumpability
[12] or its variants [15]. Instead of an approximate block structure of the
transition-matrix used in that context, we base our partitioning on a
segmentation of the molecule counts. Moreover, during the iterative refinement
of our abstraction, we identify those regions of the state-space that
contribute most to the bridging distribution. In particular, we refine those
lumped states that have a bridging probability above a certain threshold
$\delta$ and truncate all other macro-states. This way, the algorithm learns a
truncation capturing most of the bridging probabilities. This truncation
provides guaranteed lower bounds because it is at the granularity of the
original model.

In the rest of the paper, after presenting related work (Section 2) and
background (Section 3), we discuss the method (Section 4) and several
applications, including the computation of rare event probabilities as well as
Bayesian smoothing and filtering (Section 5).

##  2 Related Work

The problem of endpoint constrained analysis occurs in the context of Bayesian
estimation [41]. For population-structured MJPs, this problem has been
addressed by Huang et al. [25] using moment closure approximations and by
Wildner and Köppl [48] further employing variational inference. Golightly and
Sherlock modified stochastic simulation algorithms to approximatively augment
generated trajectories [17]. Since a statistically exact augmentation is only
possible for few simple cases, diffusion approximations [18] and moment
approximations [35] have been employed. Such approximations, however, do not
give any guarantees on the approximation error and may suffer from numerical
instabilities [43].

The bridging problem also arises during the estimation of first passage times
and rare event analysis. Approaches for first-passage times are often of
heuristic nature [42, 22, 8]. Rigorous approaches yielding guaranteed bounds
are currently limited by the performance of state-of-the-art optimization
software [6]. In biological applications, rare events of interest are
typically related to the reachability of certain thresholds on molecule counts
or mode switching [45]. Most methods for the estimation of rare event
probabilities rely on importance sampling [26, 14]. For other queries,
alternative variance reduction techniques such as control variates are
available [5]. Apart from sampling-based approaches, dynamic finite-state
projections have been employed by Mikeev et al. [34], but are lacking
automated truncation schemes.

The analysis of countably infinite state-spaces is often handled by a pre-
defined truncation [27]. Sophisticated state-space truncations for the
(unconditioned) forward analysis have been developed to give lower bounds and
rely on a trade-off between computational load and tightness of the bound [37,
28, 4, 24, 31].

Reachability analysis, which is relevant in the context of probabilistic
verification [8, 38], is a bridging problem where the endpoint constraint is
the visit of a set of goal states. Backward probabilities are commonly used to
compute reachability likelihoods [2, 50]. Approximate techniques for
reachability, based on moment closure and stochastic approximation, have also
been developed in [8, 9], but lack error guarantees. There is also a
conceptual similarity between computing bridging probabilities and the
forward-backward algorithm for computing state-wise posterior marginals in
hidden Markov models (HMMs) [40]. Like MJPs, HMMs are a generative model that
can be conditioned on observations. We only consider two observations (initial
and terminal state) that are not necessarily noisy but the forward and
backward probabilities admit the same meaning.

##  3 Preliminaries

###  3.1 Markov Jump Processes with Population Structure

A population-structured Markov jump process (MJP) describes the stochastic
interactions among agents of distinct types in a well-stirred reactor. The
assumption of all agents being equally distributed in space, allows to only
keep track of the overall copy number of agents for each type. Therefore the
state-space is $\mathcal{S}\subseteq\mathbb{N}^{n_{S}}$ where $n_{S}$ denotes
the number of agent types or populations. Interactions between agents are
expressed as _reactions_. These reactions have associated gains and losses of
agents, given by non-negative integer vectors ${v}_{j}^{-}$ and ${v}_{j}^{+}$
for reaction $j$, respectively. The overall effect is given by
$v_{j}={v}_{j}^{+}-{v}_{j}^{-}$. A reaction between agents of types
$S_{1},\dots,S_{n_{S}}$ is specified in the following form:

$\sum_{\ell=1}^{n_{S}}v_{j\ell}^{-}S_{\ell}\xrightarrow{\alpha_{j}(x)}\sum_{\ell=1}^{n_{S}}v_{j\ell}^{+}S_{\ell}\,.$
(1)

The propensity function $\alpha_{j}$ gives the rate of the exponentially
distributed firing time of the reaction as a function of the current system
state $x\in\mathcal{S}$. In population models, _mass-action_ propensities are
most common. In this case the firing rate is given by the product of the
number of reactant combinations in $x$ and a _rate constant_ $c_{j}$, i.e.

$\alpha_{j}({x})\coloneqq
c_{j}\prod_{\ell=1}^{n_{S}}\binom{x_{\ell}}{v_{j\ell}^{-}}\,.$ (2)

In this case, we give the rate constant in (1) instead of the function
$\alpha_{j}$. For a given set of $n_{R}$ reactions, we define a stochastic
process $\\{{{X}}_{t}\\}_{t\geq 0}$ describing the evolution of the population
sizes over time $t$. Due to the assumption of exponentially distributed firing
times, $X$ is a continuous-time Markov chain (CTMC) on $\mathcal{S}$ with
infinitesimal generator matrix $Q$, where the entries of $Q$ are

$Q_{x,y}=\begin{cases}\sum_{j:x+v_{j}=y}\alpha_{j}(x)\,,&\text{if}\;x\neq
y,\\\\[4.30554pt]
-\sum_{j=1}^{n_{R}}\alpha_{j}(x)\,,&\text{otherwise.}\end{cases}$ (3)

The probability distribution over time can be analyzed as an initial value
problem. Given an initial state $x_{0}$, the distribution111In the sequel,
$x_{i}$ denotes a state with index $i$ instead of its $i$-th component.

$\pi(x_{i},t)=\Pr(X_{t}=x_{i}\mid X_{0}=x_{0}),\quad t\geq 0$ (4)

evolves according to the Kolmogorov forward equation

$\frac{d}{dt}\pi(t)=\pi(t)Q\,,$ (5)

where $\pi(t)$ is an arbitrary vectorization
$(\pi(x_{1},t),\pi(x_{2},t),\dots,\pi(x_{|\mathcal{S}|},t))$ of the states.

Let $x_{g}\in\mathcal{S}$ be a fixed goal state. Given the terminal constraint
$\Pr(X_{T}=x_{g})$ for some $T\geq 0$, we are interested in the so-called
backward probabilities

$\beta(x_{i},t)=\Pr(X_{T}=x_{g}\mid X_{t}=x_{i}),\quad t\leq T\,.$ (6)

Note that $\beta(\cdot,t)$ is a function of the conditional event and thus is
no probability distribution over the state-space. Instead $\beta(\cdot,t)$
gives the reaching probabilities for all states over the time span of $[t,T]$.
To compute these probabilities, we can employ the Kolmogorov backward equation

$\frac{d}{dt}\beta(t)=Q\beta(t)^{\top}\,,$ (7)

where we use the same vectorization to construct $\beta(t)$ as we used for
$\pi(t)$. The above equation is integrated backwards in time and yields the
reachability probability for each state $x_{i}$ and time $t<T$ of ending up in
$x_{g}$ at time $T$.

The state-space of many MJPs with population structure, even simple ones, is
countably infinite. In this case, we have to truncate the state-space to a
_reasonable_ finite subset. The choice of this truncation heavily depends on
the goal of the analysis. If one is interested in the most “common” behavior,
for example, a dynamic mass-based truncation scheme is most appropriate [32].
Such a scheme truncates states with small probability during the numerical
integration. However, common mass-based truncation schemes are not as useful
for the bridging problem. This is because trajectories that meet the specific
terminal constraints can be far off the main bulk of the probability mass. We
solve this problem by a state-space lumping in connection with an iterative
refinement scheme.

Consider as an example a birth-death process. This model can be used to model
a wide variety of phenomena and often constitutes a sub-module of larger
models. For example, it can be interpreted as an M/M/1 queue with service
rates being linearly dependent on the queue length. Note, that even for this
simple model, the state-space is countably infinite.

######  Model 3.1 (Birth-Death Process)

The model consists of exponentially distributed arrivals and service times
proportional to queue length. It can be expressed using two mass-action
reactions:

$\varnothing\xrightarrow{10}X\qquad\text{and}\qquad
X\xrightarrow{.1}\varnothing\,.$

The initial condition $X_{0}=0$ holds with probability one.

###  3.2 Bridging Distribution

The process’ probability distribution given both initial and terminal
constraints is formally described by the conditional probabilities

$\gamma(x_{i},t)=\Pr(X_{t}=x_{i}\mid X_{0}=x_{0},X_{T}=x_{g}),\quad 0\leq
t\leq T$ (8)

for fixed initial state $x_{0}$ and terminal state $x_{g}$. We call these
probabilities the _bridging probabilities_. It is straight-forward to see that
$\gamma$ admits the factorization

$\gamma(x_{i},t)=\pi(x_{i},t)\beta(x_{i},t)/\pi(x_{g},T)$ (9)

due to the Markov property. The normalization factor, given by the
reachability probability $\pi(x_{g},T)=\beta(x_{0},0)$, ensures that
$\gamma(\cdot,t)$ is a distribution for all time points $t\in[0,T]$. We call
each $\gamma(\cdot,t)$ a _bridging distribution_. From the Kolmogorov
equations (5) and (7) we can obtain both the forward probabilities
$\pi(\cdot,t)$ and the backward probabilities $\beta(\cdot,t)$ for $t<T$.

We can easily extend this procedure to deal with hitting times constrained by
a finite time-horizon by making the goal state $x_{g}$ absorbing.

In Figure 1 we plot the forward, backward, and bridging probabilities for
Model 3.1. The probabilities are computed on a $[0,100]$ state-space
truncation. The approximate forward solution $\hat{\pi}$ shows how the
probability mass drifts upwards towards the stationary distribution
$\text{Poisson}(100)$. The backward probabilities are highest for states below
the goal state $x_{g}=40$. This is expected because upwards drift makes
reaching $x_{g}$ more probable for “lower” states. Finally, the approximate
bridging distribution $\hat{\gamma}$ can be recognized to be proportional to
the product of forward $\hat{\pi}$ and backward probabilities $\hat{\beta}$.

Figure 1: Forward, backward, and bridging probabilities for Model 3.1 with
initial constraint $X_{0}=0$ and terminal constraint $X_{10}=40$ on a
truncated state-space. Probabilities over $0.1$ in $\hat{\pi}$ and
$\hat{\beta}$ are given full intensity for visual clarity. The lightly shaded
area ($\geq 60$) indicates a region being more relevant for the forward than
for the bridging probabilities.

##  4 Bridge Truncation via Lumping Approximations

We first discuss the truncation of countably infinite state-spaces to analyze
backward and forward probabilities (Section 4.1). To identify effective
truncations we employ a lumping scheme. In Section 4.2, we explain the
construction of macro-states and assumptions made, as well as the efficient
calculation of transition rates between them. Finally, in Section 4.3 we
present an iterative refinement algorithm yielding a suitable truncation for
the bridging problem.

###  4.1 Finite State Projection

Even in simple models such as a birth-death Process (Model 3.1), the reachable
state-space is countably infinite. Direct analyzes of backward (6) and forward
equations (4) are often infeasible. Instead, the integration of these
differential equations requires working with a finite subset of the infinite
state-space [37]. If states are truncated, their incoming transitions from
states that are not truncated can be re-directed to a _sink state_. The
accumulated probability in this sink state is then used as an error estimate
for the forward integration scheme. Consequently, many truncation schemes,
such as dynamic truncations [4], aim to minimize the amount of “lost mass” of
the forward probability. We use the same truncation method but base the
truncation on bridging probabilities rather than the forward probabilities.

###  4.2 State-Space Lumping

When dealing with bridging problems, the most likely trajectories from the
initial to the terminal state are typically not known a priori. Especially if
the event in question is rare, obtaining a state-space truncation adapted to
its constraints is difficult. We devise a lumping scheme that groups nearby
states, i.e. molecule counts, into larger _macro-states_. A macro-state is a
collection of states treated as one state in a lumped model, which can be seen
as an abstraction of the original model. These macro-states form a
partitioning of the state-space. In this lumped model, we assume a uniform
distribution over the constituent micro-states inside each macro-state. Thus,
given that the system is in a particular macro-state, all of its micro-states
are equally likely. This partitioning allows us to analyze significant regions
of the state-space efficiently albeit under a rough approximation of the
dynamics. Iterative refinement of the state-space after each analysis moves
the dynamics closer to the original model. In the final step of the iteration,
the considered system states are at the granularity of the original model such
that no approximation error is introduced by assumptions of the lumping
scheme. Computational efficiency is retained by truncating in each iteration
step those states that contribute little probability mass to the
(approximated) bridging distributions.

We choose a lumping scheme based on a grid of hypercube macro-states whose
endpoints belong to a predefined grid. This topology makes the computation of
transition rates between macro-states particularly convenient. Mass-action
reaction rates, for example, can be given in a closed-form due to the
Faulhaber formulae. More complicated rate functions such as Hill functions can
often be handled as well by taking appropriate integrals.

Our choice is a scheme that uses $n_{S}$-dimensional hypercubes. A macro-state
$\bar{x}_{i}(\ell^{(i)},u^{(i)})$ (denoted by $\bar{x}_{i}$ for notational
ease) can therefore be described by two vectors $\ell^{(i)}$ and $u^{(i)}$.
The vector $\ell^{(i)}$ gives the corner closest to the origin, while
$u^{(i)}$ gives the corner farthest from the origin. Formally,

$\bar{x}_{i}=\bar{x}_{i}(\ell^{(i)},u^{(i)})=\\{x\in\mathbb{N}^{n_{S}}\mid\ell^{(i)}\leq
x\leq u^{(i)}\\},$ (10)

where ’$\leq$’ stands for the element-wise comparison. This choice of topology
makes the computation of transition rates between macro-states particularly
convenient: Suppose we are interested in the set of micro-states in macro-
state $\bar{x}_{i}$ that can transition to macro-state $\bar{x}_{k}$ via
reaction $j$. It is easy to see that this set is itself an interval-defined
macro-state $\bar{x}_{i\xrightarrow{j}k}$. To compute this macro-state we can
simply shift $\bar{x}_{i}$ by $v_{j}$, take the intersection with
$\bar{x}_{k}$ and project this set back. Formally,

$\bar{x}_{i\xrightarrow{j}k}=((\bar{x}_{i}+v_{j})\cap\bar{x}_{k})-v_{j}\,,$
(11)

where the additions are applied element-wise to all states making up the
macro-states. For the correct handling of the truncation it is useful to
define a general exit state

$\bar{x}_{i\xrightarrow{j}}=((\bar{x}_{i}+v_{j})\setminus\bar{x}_{i})-v_{j}.$
(12)

This state captures all micro-states inside $\bar{x}_{i}$ that can leave the
state via reaction $j$. Note that all operations preserve the structure of a
macro-state as defined in (10). Since a macro-state is based on intervals the
computation of the transition rate is often straight-forward. Under the
assumption of polynomial rates, as it is the case for mass-action systems, we
can compute the sum of rates over this transition set efficiently using
Faulhaber’s formula. We define the lumped transition function

$\bar{\alpha}_{j}({\bar{x}})=\sum_{x\in\bar{x}}\alpha_{j}(x)$ (13)

for macro-state $\bar{x}$ and reaction $j$. As an example consider the
following mass-action reaction $2X\xrightarrow{c}\varnothing\,.$ For macro-
state $\bar{x}=\\{0,\dots,n\\}$ we can compute the corresponding lumped
transition rate

$\bar{\alpha}(\bar{x})=\frac{c}{2}\sum_{i=1}^{n}i(i-1)=\frac{c}{2}\sum_{i=1}^{n}(i^{2}-i)=\frac{c}{2}\left(\frac{2n^{3}+3n^{2}+n}{6}-\frac{n^{2}+n}{2}\right)$

eliminating the explicit summation in the lumped propensity function.

For polynomial propensity functions $\alpha$ such formulae are easily obtained
automatically. For non-polynomial propensity functions, we can use the
continuous integral as an approximation. This is demonstrated on a case study
in Section 5.2.2.

Figure 2: A lumping approximation of Model 3.1 on the state-space truncation
to $[0,200]$ on $t\in[0,50]$. On the left-hand side solutions of a regular
truncation approximation and a lumped truncation (macro-state size is 5) are
given. On the right-hand side the respective terminal distributions
$\Pr(X_{50}=x_{i})$ are contrasted.

Using the transition set computation (11) and the lumped propensity function
(13) we can populate the $Q$-matrix of the finite lumping approximation:

$\bar{Q}_{\bar{x}_{i},\bar{x}_{k}}=\begin{cases}\sum_{j=1}^{n_{R}}{\bar{\alpha}}_{j}\left(\bar{x}_{i\xrightarrow{j}k}\right)/{\mathrm{vol}\left(\bar{x}_{i}\right)}\,,&\text{if}\;\bar{x}_{i}\neq\bar{x}_{k}\\\\[4.30554pt]
-\sum_{j=1}^{n_{R}}{\bar{\alpha}}_{j}\left(\bar{x}_{i\xrightarrow{j}}\right)/{\mathrm{vol}\left(\bar{x}_{i}\right)}\,,&\text{otherwise}\end{cases}$
(14)

In addition to the lumped rate function over the transition state
$\bar{x}_{i\xrightarrow{j}k}$, we need to divide by the total volume of the
lumped state $\bar{x}_{i}$. This is due to the assumption of a uniform
distribution inside the macro-states. Using this $Q$-matrix, we can compute
the forward and backward solution using the respective Kolmogorov equations
(5) and (7).

Interestingly, the lumped distribution tends to be less concentrated. This is
due to the assumption of a uniform distribution inside macro-states. This
effect is illustrated by the example of a birth-death process in Figure 2. Due
to this effect, an iterative refinement typically keeps an over-approximation
in terms of state-space area. This is a desirable feature since relevant
regions are less likely to be pruned due to lumping approximations.

###  4.3 Iterative Refinement Algorithm

The iterative refinement algorithm (Alg. 1) starts with a set of large macro-
states that are iteratively refined, based on approximate solutions to the
bridging problem. We start by constructing square macro-states of size $2^{m}$
in each dimension for some $m\in\mathbb{N}$ such that they form a large-scale
grid $\mathcal{S}^{(0)}$. Hence, each initial macro-state has a volume of
${\left(2^{m}\right)}^{n_{S}}$. This choice of grid size is convenient because
we can halve states in each dimension. Moreover, this choice ensures that all
states have equal volume and we end up with states of volume $2^{0}=1$ which
is equivalent to a truncation of the original non-lumped state-space.
An iteration of the state-space refinement starts by computing both the
forward and backward probabilities (lines 1 and 1) via integration of (5) and
(7), respectively, using the lumped $\hat{Q}$-matrix. Based on the resulting
approximate forward and backward probabilities, we compute an approximation of
the bridging distributions (line 1). This is done for each time-point in an
equispaced grid on $[0,T]$. The time grid granularity is a hyper-parameter of
the algorithm. If the grid is too fine, the memory overhead of storing
backward $\hat{\beta}^{(i)}$ and forward solutions $\hat{\pi}^{(i)}$
increases.222We denote the approximations with a hat (e.g. $\hat{\pi}$) rather
than a bar (e.g. $\bar{\pi}$) to indicate that not only the lumping
approximation but also a truncation is applied and similarly for the
$Q$-matrix. If, on the other hand, the granularity is too low, too much of the
state-space might be truncated. Based on a threshold parameter $\delta>0$
states are either removed or split (line 1), depending on the mass assigned to
them by the approximate bridging probabilities $\hat{\gamma}^{(i)}_{t}$. A
state can be split by the split-function which halves the state in each
dimension. Otherwise, it is removed. Thus, each macro-state is either split
into $2^{n_{S}}$ new states or removed entirely. The result forms the next
lumped state-space $\mathcal{S}^{(i+1)}$. The $Q$-matrix is adjusted (line 1)
such that transition rates for $\mathcal{S}^{(i+1)}$ are calculated according
to (14). Entries of truncated states are removed from the transition matrix.
Transitions leading to them are re-directed to a sink state (see Section 4.1).
After $m$ iterations (we started with states of side lengths $2^{m}$) we have
a standard finite state projection scheme on the original model tailored to
computing an approximation of the bridging distribution.

input : Initial partitioning $\mathcal{S}^{(0)}$, truncation threshold
$\delta$

output : approximate bridging distribution $\hat{\gamma}$

1 for  _$i=1,\dots,m$_ do

2 ${\hat{\pi}}^{(i)}_{t}\leftarrow$ solve approximate forward equation on
$\mathcal{S}^{(i)}$;

3 ${\hat{\beta}}^{(i)}_{t}\leftarrow$ solve approximate backward equation on
$\mathcal{S}^{(i)}$;

${\hat{\gamma}}^{(i)}_{t}\leftarrow{\hat{\beta}}^{(i)}{\hat{\pi}}^{(i)}/\hat{\pi}(x_{g},T)$;

/* approximate bridging distribution */

4 $\mathcal{S}^{(i+1)}\leftarrow\emptyset$;

5 foreach  _$\bar{x}\in\mathcal{S}^{(i)}$_ do

_ _ if  _ $\exists t.{\hat{\gamma}}^{(i)}_{t}(\bar{x})\geq\delta$;_

__ /* refine based on bridging probabilities */

6  __ then

7
$\mathcal{S}^{(i+1)}\leftarrow\mathcal{S}^{(i+1)}\cup\textnormal{{split}}(\bar{x})$;

8

9

10 update $\hat{Q}$-matrix;

11

12return ${\hat{\gamma}}^{(i)}$;

Algorithm 1 Iterative refinement for the bridging problem

In Figure 3 we give a demonstration of how Algorithm 1 works to refine the
state-space iteratively. Starting with an initial lumped state-space
$\mathcal{S}^{(0)}$ covering a large area of the state-space, repeated
evaluations of the bridging distributions are performed. After five iterations
the remaining truncation includes all states that significantly contribute to
the bridging probabilities over the times $[0,T]$.

Figure 3: The state-space refinement algorithm on two parallel unit-rate
arrival processes. The bridging problem from $(0,0)$ to $(64,64)$ and $T=10$
and truncation threshold $\delta=\text{5e-3}$. States with a bridging
probability below $\delta$ are light grey. The macro-state containing the goal
state is marked in black. The initial macro-states are of size $16\times 16$.

It is important to realize that determining the most relevant states is _the_
main challenge. The above algorithm solves this problem by considering only
those parts of the state-space that contribute most to the bridging
probabilities. The truncation is tailored to this condition and might ignore
regions that are likely in the unconditioned case. For instance, in Fig. 1 the
bridging probabilities mostly remain below a population threshold of $\\#X=60$
(as indicated by the lighter/darker coloring), while the forward probabilities
mostly exceed this bound. Hence, in this example a significant portion of the
forward probabilities $\hat{\pi}_{t}^{(i)}$ is captured by the sink state.
However, the condition in line 1 in Algorithm 1 ensures that states
contributing significantly to $\hat{\gamma}_{t}^{(i)}$ will be kept and
refined in the next iteration.

##  5 Results

We present four examples in this section to evaluate our proposed method. A
prototype was implemented in Python 3.8. For numerical integration we used the
Scipy implementation [47] of the implicit method based on backward-
differentiation formulas [13]. The analysis as a Jupyter notebook is made
available online333https://www.github.com/mbackenkoehler/mjp˙bridging.

###  5.1 Bounding Rare Event Probabilities

We consider a simple model of two parallel Poisson processes describing the
production of two types of agents. The corresponding probability distribution
has Poisson product form at all time points $t\geq 0$ and hence we can compare
the accuracy of our numerical results with the exact analytic solution. We use
the proposed approach to compute lower bounds for rare event probabilities.
444These bounds are rigorous up to the approximation error of the numerical
integration scheme. However, the forward solution could be replaced by an
adaptive uniformization approach [3] for a more rigorous integration error
control.

######  Model 5.1 (Parallel Poisson Processes)

The model consists of two parallel independent Poisson processes with unit
rates.

$\varnothing\xrightarrow{1}A\qquad\text{and}\qquad\varnothing\xrightarrow{1}B\,.$

The initial condition $X_{0}=(0,0)$ holds with probability one. After $t$ time
units each species abundance is Poisson distributed with rate $\lambda=t$.

We consider the final constraint of reaching a state where both processes
exceed a threshold of 64 at time 20. Without prior knowledge, a reasonable
truncation would have been $160\times 160$. But our analysis shows that just
20% of the states are necessary to capture over 99.6% of the probability mass
reaching the target event (cf. Table 1). Decreasing the threshold $\delta$
leads to a larger set of states retained after truncation as more of the
bridging distribution is included (cf. Figure 4). We observe an increase in
truncation size that is approximately logarithmic in $\delta$, which, in this
example, indicates robustness of the method with respect to the choice of
$\delta$.

Figure 4: State-space truncation for varying values of the threshold parameter $\delta$: Two parallel Poisson processes under terminal constraints $X_{20}^{(A)}\geq 64$ and $X_{20}^{(B)}\geq 64$. The initial macro-states are $16\times 16$ such that the final states are regular micro states. threshold $\delta$ | 1e-2 | 1e-3 | 1e-4 | 1e-5
---|---|---|---|---
truncation size | 1154 | 2354 | 3170 | 3898
overall states | 2074 | 3546 | 4586 | 5450
estimate | 8.8851e-30 | 1.8557e-29 | 1.8625e-29 | 1.8625e-29
rel. error | 5.2297e-01 | 3.6667e-03 | 3.7423e-05 | 9.5259e-08
Table 1: Estimated reachability probabilities based on varying truncation
thresholds $\delta$: The true probability is 1.8625e-29. We also report the
size of the final truncation and the accumulated size of all truncations
during refinement iterations (overall states).

##### Comparison to other methods

The truncation approach that we apply is similar to the one used by Mikeev et
al. [34] for rare event estimation. However, they used a given linearly biased
MJP model to obtain a truncation. A general strategy to compute an appropriate
biasing was not proposed. It is possible to adapt our truncation approach to
the dynamic scheme in Ref. [34] where states are removed in an on-the-fly
fashion during numerical integration.

A finite state-space truncation covering the same area as the initial lumping
approximation would contain $25,\\!600$ states.555Here, the goal is not
treated as a single state. Otherwise, it consists of $24,\\!130$ states. The
standard approach would be to build up the entire state-space for such a model
[27]. Even using a conservative truncation threshold $\delta=\text{1e-5}$, our
method yields an accurate estimate using only about a fifth (5450) of this
accumulated over all intermediate lumped approximations.

###  5.2 Mode Switching

Mode switching occurs in models exhibiting _multi-modal_ behavior [44] when a
trajectory traverses a potential barrier from one mode to another. Often, mode
switching is a rare event and occurs in the context of gene regulatory
networks where a mode is characterized by the set of genes being currently
active [30]. Similar dynamics also commonly occur in queuing models where a
system may for example switch its operating behavior stochastically if traffic
increases above or decreases below certain thresholds. Using the presented
method, we can get both a qualitative and quantitative understanding of
switching behavior without resorting to Monte-Carlo methods such as importance
sampling.

####  5.2.1 Exclusive Switch

The exclusive switch [7] has three different modes of operation, depending on
the DNA state, i.e. on whether a protein of type one or two is bound to the
DNA.

######  Model 5.2 (Exclusive Switch)

The exclusive switch model consists of a promoter region that can express both
proteins $P_{1}$ and $P_{2}$. Both can bind to the region, suppressing the
expression of the other protein. For certain parameterizations, this leads to
a bi-modal or even tri-modal behavior.

$D\xrightarrow{\rho}D+P_{1}\qquad D\xrightarrow{\rho}D+P_{2}\qquad
P_{1}\xrightarrow{\lambda}\varnothing\qquad
P_{2}\xrightarrow{\lambda}\varnothing$
$D+P_{1}\xrightarrow{\beta}D.P_{1}\qquad D.P1\xrightarrow{\gamma}D+P_{1}\qquad
D.P_{1}\xrightarrow{\alpha}D.P_{1}+P_{1}$
$D+P_{2}\xrightarrow{\beta}D.P_{2}\qquad D.P2\xrightarrow{\gamma}D+P_{2}\qquad
D.P_{2}\xrightarrow{\alpha}D.P_{2}+P_{2}$

The parameter values are $\rho=\text{1e-1}$, $\lambda=\text{1e-3}$,
$\beta=\text{1e-2}$, $\gamma=\text{8e-3}$, and $\alpha=\text{1e-1}$.

Since we know a priori of the three distinct operating modes, we adjust the
method slightly: The state-space for the DNA states is not lumped. Instead we
“stack” lumped approximations of the $P_{1}$-$P_{2}$ phase space upon each
other. Special treatment of DNA states is common for such models [28].

To analyze the switching, we choose the transition from (variable order:
$P_{1}$, $P_{2}$, $D$, $D.P_{1}$, $D.P_{2}$) $x_{1}=(32,0,0,0,1)$ to
$x_{2}=(0,32,0,1,0)$ over the time interval $t\in[0,10]$. The initial lumping
scheme covers up to 80 molecules of $P_{1}$ and $P_{2}$ for each mode. Macro-
states have size $8\times 8$ and the truncation threshold is
$\delta=\text{1e-4}$.

In the analysis of biological switches, not only the switching probability but
also the switching dynamics is a central part of understanding the underlying
biological mechanisms. In Figure 5 (left), we therefore plot the time-varying
probabilities of the gene state conditioned on the mode. We observe a rapid
unbinding of $P_{2}$, followed by a slow increase of the binding probability
for $P_{1}$. These dynamics are already qualitatively captured by the first
lumped approximation (dashed lines).

Figure 5: (left) Mode probabilities of the exclusive switch bridging problem
over time for the first lumped approximation (dashed lines) and the final
approximation (solid lines) with constraints $X_{0}=(32,0,0,1,0)$ and
$X_{10}=(0,32,0,0,1)$. (right) The expected occupation time (excluding initial
and terminal states) for the switching problem of the toggle switch using
Hill-type functions. The bridging problem is from initial $(0,120)$ to a first
passage of $(120,0)$ in $t\in[0,10]$.

####  5.2.2 Toggle Switch

Next, we apply our method to a toggle switch model exhibiting non-polynomial
rate functions. This well-known model considers two proteins $A$ and $B$
inhibiting the production of the respective other protein [29].

###### Model 5.3

Toggle Switch (Hill functions) We have population types $A$ and $B$ with the
following reactions and reaction rates.

$\varnothing\xrightarrow{\alpha_{1}(\cdot)}A\,,\quad\text{where}\quad\alpha_{1}(x)=\frac{\rho}{1+x_{B}},\qquad
A\xrightarrow{\lambda}\varnothing$
$\varnothing\xrightarrow{\alpha_{1}(\cdot)}B\,,\quad\text{where}\quad\alpha_{1}(x)=\frac{\rho}{1+x_{A}},\qquad
B\xrightarrow{\lambda}\varnothing$

The parameterization is $\rho=10$, $\lambda=0.1$.

Due to the non-polynomial rate functions $\alpha_{1}$ and $\alpha_{2}$, the
transition rates between macro-states are approximated by using the continuous
integral

$\bar{\alpha}_{1}(\bar{x})\approx\int_{a-0.5}^{b+0.5}\frac{\rho}{1+x}\,dx=\rho\left(\log{\left(b+1.5\right)}-\log{\left(a+0.5\right)}\right)$

for a macro-state $\bar{x}=\\{a,\dots,b\\}$.

We analyze the switching scenario from $(0,120)$ to the first visit of state
$(120,0)$ up to time $T=10$. The initial lumping scheme covers up to 352
molecules of $A$ and $B$ and macro-states have size $32\times 32$. The
truncation threshold is $\delta=\text{1e-4}$. The resulting truncation is
shown in Figure 5 (right). It also illustrates the kind of insights that can
be obtained from the bridging distributions. For an overview of the switching
dynamics, we look at the expected occupation time under the terminal
constraint of having entered state $(120,0)$. Letting the corresponding
hitting time be $\tau=\inf\\{t\geq 0\mid X_{t}=(120,0)\\}$, the expected
occupation time for some state $x$ is
$E\left(\int_{0}^{\tau}1_{=x}(X_{t})\,dt\mid\tau\leq 10\right)$. We observe
that in this example the switching behavior seems to be asymmetrical. The main
mass seems to pass through an area where initially a small number of $A$
molecules is produced followed by a total decay of $B$ molecules.

###  5.3 Recursive Bayesian Estimation

We now turn to the method’s application in recursive Bayesian estimation. This
is the problem of estimating the system’s past, present, and future behavior
under given observations. Thus, the MJP becomes a hidden Markov model (HMM).
The observations in such models are usually noisy, meaning that we cannot
infer the system state with certainty.

This estimation problem entails more general distributional constraints on
terminal $\beta(\cdot,T)$ and initial $\pi(\cdot,0)$ distributions than the
point mass distributions considered up until now. We can easily extend the
forward and backward probabilities to more general initial distributions and
terminal distributions $\beta(T)$. For the forward probabilities we get

$\pi(x_{i},t)=\sum_{j}\Pr(X_{t}=x_{i}\mid X_{0}=x_{j})\pi(x_{j},0),$ (15)

and similarly the backward probabilities are given by

$\beta(x_{i},t)=\sum_{j}\Pr(X_{T}=x_{j}\mid X_{t}=x_{i})\beta_{T}(x_{j})\,.$
(16)

We apply our method to an SEIR (susceptible-exposed-infected-removed) model.
This is widely used to describe the spreading of an epidemic such as the
current COVID-19 outbreak [23, 20]. Temporal snapshots of the epidemic spread
are mostly only available for a subset of the population and suffer from
inaccuracies of diagnostic tests. Bayesian estimation can then be used to
infer the spreading dynamics given uncertain temporal snapshots.

######  Model 5.4 (Epidemics Model)

A population of susceptible individuals can contract a disease from infected
agents. In this case, they are exposed, meaning they will become infected but
cannot yet infect others. After being infected, individuals change to the
removed state. The mass-action reactions are as follows.

$S+I\xrightarrow{\lambda}E+I\qquad E\xrightarrow{\mu}I\qquad
I\xrightarrow{\rho}R$

The parameter values are $\lambda=0.5$, $\mu=3$, $\rho=3$. Due to the
stoichiometric invariant
$X_{t}^{(S)}+X_{t}^{(E)}+X_{t}^{(I)}+X_{t}^{(R)}=\mathrm{const.}$, we can
eliminate $R$ from the system.

We consider the following scenario: We know that initially ($t=0$) one
individual is infected and the rest is susceptible. At time $t=0.3$ all
individuals are tested for the disease. The test, however, only identifies
infected individuals with probability $0.99$. Moreover, the probability of a
false positive is 0.05. We like to identify the distribution given both the
initial state and the measurement at time $t=0.3$. In particular, we want to
infer the distribution over the latent counts of $S$ and $E$ by _recursive
Bayesian estimation_.

The posterior for $n_{I}$ infected individuals at time $t$, given measurement
$Y_{t}=\hat{n}_{I}$ can be computed using Bayes’ rule

$\Pr(X_{t}^{(I)}=n_{I}\mid Y_{t}=\hat{n}_{I})\propto\Pr(Y_{t}=\hat{n}_{I}\mid
X_{t}^{(I)}=n_{I})\Pr(X_{t}^{(I)}=n_{I})\,.$ (17)

This problem is an extension of the bridging problem discussed up until now.
The difference is that the terminal posterior is estimated it using the result
of the lumped forward equation and the measurement distribution using (17).
Based on this estimated terminal posterior, we compute the bridging
probabilities and refine the truncation tailored to the location of the
posterior distribution. In Figure 6 (left), we illustrate the bridging
distribution between the terminal posterior and initial distribution. In the
context of filtering problems this is commonly referred to as smoothing. Using
the learned truncation, we can obtain the posterior distribution for the
number of infected individuals at $t=0.3$ (Figure 6 (middle)). Moreover, can
we infer a distribution over the unknown number of susceptible and exposed
individuals (Figure 6 (right)).

Figure 6: (left) A comparison of the prior dynamics and the posterior
smoothing (bridging) dynamics. (middle) The prior, likelihood, and posterior
of the number of infected individuals $n_{I}$ at time $t=0.3$ given the
measurement $\hat{n}_{I}=30$. (right) The prior and posterior distribution
over the latent types $E$ and $S$.

##  6 Conclusion

The analysis of Markov Jump processes with constraints on the initial and
terminal behavior is an important part of many probabilistic inference tasks
such as parameter estimation using Bayesian or maximum likelihood estimation,
inference of latent system behavior, the estimation of rare event
probabilities, and reachability analysis for the verification of temporal
properties. If endpoint constraints correspond to atypical system behaviors,
standard analysis methods fail as they have no strategy to identify those
parts of the state-space relevant for meeting the terminal constraint.

Here, we proposed a method that is not based on stochastic sampling and
statistical estimation but provides a direct numerical approach. It starts
with an abstract lumped model, which is iteratively refined such that only
those parts of the model are considered that contribute to the probabilities
of interest. In the final step of the iteration, we operate at the granularity
of the original model and compute lower bounds for these bridging
probabilities that are rigorous up to the error of the numerical integration
scheme.

Our method exploits the population structure of the model, which is present in
many important application fields of MJPs. Based on experience with other work
based on truncation, the approach can be expected to scale up to at least a
few million states [33]. Compared to previous work, our method neither relies
on approximations of unknown accuracy nor additional information such as a
suitable change of measure in the case of importance sampling. It only
requires a truncation threshold and an initial choice for the macro-state
sizes.

In future work, we plan to extend our method to hybrid approaches, in which a
moment representation is employed for large populations while discrete counts
are maintained for small populations. Moreover, we will apply our method to
model checking where constraints are described by some temporal logic [21].

####  6.0.1 Acknowledgements

This project was supported by the DFG project MULTIMODE and Italian PRIN
project SEDUCE.

## References

  * [1] Adan, A., Alizada, G., Kiraz, Y., Baran, Y., Nalbant, A.: Flow cytometry: basic principles and applications. Critical reviews in biotechnology 37(2), 163–176 (2017)
  * [2] Amparore, E.G., Donatelli, S.: Backward solution of Markov chains and Markov regenerative processes: Formalization and applications. Electron. Notes Theor. Comput. Sci. 296, 7–26 (2013)
  * [3] Andreychenko, A., Crouzen, P., Mikeev, L., Wolf, V.: On-the-fly uniformization of time-inhomogeneous infinite Markov population models. arXiv preprint arXiv:1006.4425 (2010)
  * [4] Andreychenko, A., Mikeev, L., Spieler, D., Wolf, V.: Parameter identification for Markov models of biochemical reactions. In: International Conference on Computer Aided Verification. pp. 83–98. Springer (2011)
  * [5] Backenköhler, M., Bortolussi, L., Wolf, V.: Control variates for stochastic simulation of chemical reaction networks. In: Bortolussi, L., Sanguinetti, G. (eds.) Computational Methods in Systems Biology. pp. 42–59. Springer, Cham (2019)
  * [6] Backenköhler, M., Bortolussi, L., Wolf, V.: Bounding mean first passage times in population continuous-time Markov chains. To appear in Proc. of QEST’20 (2020)
  * [7] Barzel, B., Biham, O.: Calculation of switching times in the genetic toggle switch and other bistable systems. Physical Review E 78(4), 041919 (2008)
  * [8] Bortolussi, L., Lanciani, R.: Stochastic approximation of global reachability probabilities of Markov population models. In: Computer Performance Engineering - 11th European Workshop, EPEW 2014, Florence, Italy, September 11-12, 2014. Proceedings. pp. 224–239 (2014)
  * [9] Bortolussi, L., Lanciani, R., Nenzi, L.: Model checking markov population models by stochastic approximations. Inf. Comput. 262, 189–220 (2018)
  * [10] Breuer, L.: From Markov jump processes to spatial queues. Springer Science & Business Media (2003)
  * [11] Broemeling, L.D.: Bayesian Inference for Stochastic Processes. CRC Press (2017)
  * [12] Buchholz, P.: Exact and ordinary lumpability in finite Markov chains. Journal of applied probability pp. 59–75 (1994)
  * [13] Byrne, G.D., Hindmarsh, A.C.: A polyalgorithm for the numerical solution of ordinary differential equations. ACM Transactions on Mathematical Software (TOMS) 1(1), 71–96 (1975)
  * [14] Daigle Jr, B.J., Roh, M.K., Gillespie, D.T., Petzold, L.R.: Automated estimation of rare event probabilities in biochemical systems. The Journal of Chemical Physics 134(4), 01B628 (2011)
  * [15] Dayar, T., Stewart, W.J.: Quasi lumpability, lower-bounding coupling matrices, and nearly completely decomposable Markov chains. SIAM Journal on Matrix Analysis and Applications 18(2), 482–498 (1997)
  * [16] Gillespie, D.T.: Exact stochastic simulation of coupled chemical reactions. The journal of physical chemistry 81(25), 2340–2361 (1977)
  * [17] Golightly, A., Sherlock, C.: Efficient sampling of conditioned Markov jump processes. Statistics and Computing 29(5), 1149–1163 (2019)
  * [18] Golightly, A., Wilkinson, D.J.: Bayesian inference for stochastic kinetic models using a diffusion approximation. Biometrics 61(3), 781–788 (2005)
  * [19] Golightly, A., Wilkinson, D.J.: Bayesian parameter inference for stochastic biochemical network models using particle Markov chain monte carlo. Interface focus 1(6), 807–820 (2011)
  * [20] Grossmann, G., Backenköhler, M., Wolf, V.: Importance of interaction structure and stochasticity for epidemic spreading: A COVID-19 case study. In: Seventeenth international conference on the quantitative evaluation of systems (QEST 2020). IEEE (2020)
  * [21] Hajnal, M., Nouvian, M., Šafránek, D., Petrov, T.: Data-informed parameter synthesis for population markov chains. In: International Workshop on Hybrid Systems Biology. pp. 147–164. Springer (2019)
  * [22] Hayden, R.A., Stefanek, A., Bradley, J.T.: Fluid computation of passage-time distributions in large Markov models. Theoretical Computer Science 413(1), 106–141 (2012)
  * [23] He, S., Peng, Y., Sun, K.: SEIR modeling of the COVID-19 and its dynamics. Nonlinear Dynamics pp. 1–14 (2020)
  * [24] Henzinger, T.A., Mateescu, M., Wolf, V.: Sliding window abstraction for infinite Markov chains. In: International Conference on Computer Aided Verification. pp. 337–352. Springer (2009)
  * [25] Huang, L., Pauleve, L., Zechner, C., Unger, M., Hansen, A.S., Koeppl, H.: Reconstructing dynamic molecular states from single-cell time series. Journal of The Royal Society Interface 13(122), 20160533 (2016)
  * [26] Kuwahara, H., Mura, I.: An efficient and exact stochastic simulation method to analyze rare events in biochemical systems. The Journal of chemical physics 129(16), 10B619 (2008)
  * [27] Kwiatkowska, M., Norman, G., Parker, D.: Prism 4.0: Verification of probabilistic real-time systems. In: International conference on computer aided verification. pp. 585–591. Springer (2011)
  * [28] Lapin, M., Mikeev, L., Wolf, V.: SHAVE: stochastic hybrid analysis of Markov population models. In: Proceedings of the 14th international conference on Hybrid systems: computation and control. pp. 311–312 (2011)
  * [29] Lipshtat, A., Loinger, A., Balaban, N.Q., Biham, O.: Genetic toggle switch without cooperative binding. Physical review letters 96(18), 188101 (2006)
  * [30] Loinger, A., Lipshtat, A., Balaban, N.Q., Biham, O.: Stochastic simulations of genetic switch systems. Physical Review E 75(2), 021904 (2007)
  * [31] Mikeev, L., Neuhäußer, M.R., Spieler, D., Wolf, V.: On-the-fly verification and optimization of DTA-properties for large Markov chains. Formal Methods in System Design 43(2), 313–337 (2013)
  * [32] Mikeev, L., Sandmann, W.: Approximate numerical integration of the chemical master equation for stochastic reaction networks. arXiv preprint arXiv:1907.10245 (2019)
  * [33] Mikeev, L., Sandmann, W., Wolf, V.: Efficient calculation of rare event probabilities in Markovian queueing networks. In: Proceedings of the 5th International ICST Conference on Performance Evaluation Methodologies and Tools. pp. 186–196 (2011)
  * [34] Mikeev, L., Sandmann, W., Wolf, V.: Numerical approximation of rare event probabilities in biochemically reacting systems. In: International Conference on Computational Methods in Systems Biology. pp. 5–18. Springer (2013)
  * [35] Milner, P., Gillespie, C.S., Wilkinson, D.J.: Moment closure based parameter inference of stochastic kinetic models. Statistics and Computing 23(2), 287–295 (2013)
  * [36] Mode, C.J., Sleeman, C.K.: Stochastic processes in epidemiology: HIV/AIDS, other infectious diseases, and computers. World Scientific (2000)
  * [37] Munsky, B., Khammash, M.: The finite state projection algorithm for the solution of the chemical master equation. The Journal of chemical physics 124(4), 044104 (2006)
  * [38] Neupane, T., Myers, C.J., Madsen, C., Zheng, H., Zhang, Z.: Stamina: Stochastic approximate model-checker for infinite-state analysis. In: International Conference on Computer Aided Verification. pp. 540–549. Springer (2019)
  * [39] Pardoux, E.: Markov processes and applications: algorithms, networks, genome and finance, vol. 796. John Wiley & Sons (2008)
  * [40] Rabiner, L., Juang, B.: An introduction to hidden Markov models. IEEE ASSP Magazine 3(1), 4–16 (1986)
  * [41] Särkkä, S.: Bayesian filtering and smoothing, vol. 3. Cambridge University Press (2013)
  * [42] Schnoerr, D., Cseke, B., Grima, R., Sanguinetti, G.: Efficient low-order approximation of first-passage time distributions. Phys. Rev. Lett. 119, 210601 (Nov 2017). https://doi.org/10.1103/PhysRevLett.119.210601
  * [43] Schnoerr, D., Sanguinetti, G., Grima, R.: Validity conditions for moment closure approximations in stochastic chemical kinetics. The Journal of chemical physics 141(8), 08B616_1 (2014)
  * [44] Siegal-Gaskins, D., Mejia-Guerra, M.K., Smith, G.D., Grotewold, E.: Emergence of switch-like behavior in a large family of simple biochemical networks. PLoS Comput Biol 7(5), e1002039 (2011)
  * [45] Strasser, M., Theis, F.J., Marr, C.: Stability and multiattractor dynamics of a toggle switch based on a two-stage model of stochastic gene expression. Biophysical journal 102(1), 19–29 (2012)
  * [46] Ullah, M., Wolkenhauer, O.: Stochastic approaches for systems biology. Springer Science & Business Media (2011)
  * [47] Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S.J., Brett, M., Wilson, J., Jarrod Millman, K., Mayorov, N., Nelson, A.R.J., Jones, E., Kern, R., Larson, E., Carey, C., Polat, İ., Feng, Y., Moore, E.W., Vand erPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E.A., Harris, C.R., Archibald, A.M., Ribeiro, A.H., Pedregosa, F., van Mulbregt, P., Contributors, S…: SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods 17, 261–272 (2020). https://doi.org/https://doi.org/10.1038/s41592-019-0686-2
  * [48] Wildner, C., Koeppl, H.: Moment-based variational inference for Markov jump processes. arXiv preprint arXiv:1905.05451 (2019)
  * [49] Wilkinson, D.J.: Stochastic modelling for systems biology. CRC press (2018)
  * [50] Zapreev, I., Katoen, J.P.: Safe on-the-fly steady-state detection for time-bounded reachability. In: Third International Conference on the Quantitative Evaluation of Systems-(QEST’06). pp. 301–310. IEEE (2006)