# Identifying Markov chain models from time-to-event data: an algebraic
approach

Ovidiu Radulescu1,∗, Dima Grigoriev2, Matthias Seiss3,
Maria Douaihy1,4, Mounia Lagha4, Edouard Bertrand5

( 1 LPHI, University of Montpellier and CNRS, Place Eugène Bataillon,
Montpellier 34095, France
2 CNRS, Mathématiques, Université de Lille, 59655, Villeneuve d’Ascq, France
3 University of Kassel, Kassel, Germany
4 IGMM, University of Montpellier and CNRS, 1919 Rte de Mende, Montpellier
34090, France
5 IGH, University of Montpellier and CNRS, 141 Rue de la Cardonille,
Montpellier 34094, France
)

###### Abstract

In many fields, including biology, medicine, physics, chemistry, economy and
actuarial science, data can be represented as the time-to-event of a finite
state Markov chain model. The distribution of time intervals between
successive recorded events is known as a phase-type distribution. We
demonstrate that in cases where the eigenvalues of the Markov chain transition
matrix are distinct (non-degenerate), the phase-type distribution is multi-
exponential. We then pose and solve an inverse problem: given knowledge of the
phase-type distribution, can we determine the transition rate parameters of
the underlying Markov chain? To tackle this challenge, we initially convert
the inverse problem into a computer algebraic task, involving the solution of
a system of polynomial equations. These equations are symmetrized with respect
to the parameters of the phase-type distribution. For a specific subset of
Markov models that we refer to as ”solvable,” the inverse problem yields a
unique solution, up to transformations by finite symmetries. We outline a
recursive approach to compute these solutions for specific families of models,
regardless of their number of states. Additionally, we use the Thomas
decomposition technique to calculate solutions for models possessing 2, 3, or
4 states. Interestingly, we show that models having the same number of states
but different transition graphs can yield identical phase-distributions. In
such ”Rashomon effect” scenarios, time-to-event data permits the formulation
of multiple models and interpretations, all of which are consistent with the
observed experimental data. To differentiate among these models, we propose
additional distinguishing properties beyond just the time until the next
event. In order to prove its applicability, we discuss how this method can be
used to infer models of transcription regulation from transcriptional bursting
data. Keywords: Markov chains, phase-type distribution, inverse problem,
symmetric polynomials, Thomas decomposition, transcriptional bursting.

##  1 Introduction

Many problems in biology, medicine, physics, chemistry, economy, actuarial
science can be modeled by using continuous time Markov chains [8]. A common
dataset produced by such models comprises a sequence of states that denote the
incidence of specific events within the system. When all states within the
model are observed, the time until the occurrence of the next event follows
the memoryless exponential distribution. However, in many practical
applications, not all the states of the model are observed. In this case the
time to the next event is no longer exponentially distributed.

A particular case is when only one state of the model can be observed. In this
case the time to the next event distribution is referred to as phase-type
distribution. Phase-type distributions were introduced by Neuts [26] and are
defined as the distribution of the time to absorption of a finite Markov chain
in the presence of only a single absorbing state. For a single observed state,
computing the time-to-event distribution is equivalent to computing a phase-
type distribution (see Figure 1). Phase-type distributions have been used in
various domains such as queueing theory [30, 29], drug kinetics [18], public
health problems [17, 25, 2, 32]. In spite of their importance for applications
it is surprising how little is known about the general properties of phase-
type distributions. Algorithmic calculations of phase-type distributions can
be found in [3]. For some general results and open questions one could also
see [27, 11, 9].

In this paper we aim at determining the Markov chain from the phase-type
distribution. This inverse problem has already been addressed using maximum
likelihood approaches in which data likelihood was computed for a chosen model
and its maximization provided the optimal model parameters [2, 32]. In
contrast to this direct inference approach, here we decompose the inference
problem into two parts. The first part consists in the regression of a
parametric multi-exponential representation of the phase-type distribution.
This part has been reviewed elsewhere [16, 13] and will not be developed here.
The second part of our approach consists in solving symbolically the inverse
phase-type distribution problem. More precisely, we compute the transition
rate parameters of all Markov chains that generate the same multi-exponential
distribution from the parameters of the distribution. This approach has
several advantages. It allows us to compute simultaneously all the Markov
chain models, having different transition parameters and diagrams, but exactly
the same likelihood. It is thus a precious tool for analyzing parametric and
structural uncertainty in inference.

Here we demonstrate the relevance of our approach by selecting the context of
transcriptional bursting, a highly regarded phenomenon in the gene expression
research community [15, 34, 28, 7, 13, 12]. However, our approach can be
applied to other fields involving phase-type distributions such as in
reliability, traffic and telecommunication engineering, actuarial science,
epidemiology, and medical research [1, 17, 10].

The structure of the paper is the following. In Section 2 we discuss the
phase-type distribution problem and find its solution in the non-degenerate
case, when the generator of the Markov chain has distinct eigenvalues. In
Section 3 we introduce the symbolic inverse problem consisting in finding the
transition rate parameters of the model from the parameters of the phase-type
distribution. We show that the inverse problem consists in solving a system of
polynomial equations, symmetric in the phase-type distribution parameters. We
also discuss solvable models, i.e. for which the inverse problem has a unique
solution. In section 4 we discuss an example of solvable model with arbitrary
number of states. In section 5 we use Thomas decomposition to solve the
inverse problem. In section 6 we discuss non-identifiability issues that we
regroup under the name Rashomon effect. In section 7 we present the
application of these findings in the analysis of transcriptional bursting
data.

$k_{1}$$k_{3}$$k_{2}$$k_{4}$$k_{5}$

Figure 1:  Phase-type distribution for a three state Markov chain model. 4
represents the unique observed state. A typical dataset is represented by the
moments when the state 4 is observed. One does not know what’s going on
between two successive occurrences of 4 but can measure the time $T_{3}$
separating them. The probability distribution of this time is called phase-
type distribution. The Markov model is represented by its transition diagram
(arcs connecting states) and by the transition rate parameters. It is
considered that from 4 the chain returns to 3 instantly. For this reason,
$T_{3}$ is the first hitting time of 4 starting from 3. The timechart
represents repeated observation of the state 4 for the same system. It could
be equivalently seen as a succession of replicas of a system that stops when
it reaches the state 4. The latter interpretation corresponds to the
traditional definition of the phase-type distribution.

$k_{1}$$k_{3}$$k_{2}$$k_{4}$$k_{5}$

$k_{1}$$k_{3}$$k_{2}$$k_{4}$$k_{5}$

$k_{1}$$k_{3}$$k_{2}$$k_{4}$$k_{5}$

$k_{2}$$k_{4}$$k_{1}$$k_{3}$$k_{5}$

$k_{1}$$k_{3}$$k_{2}$$k_{4}$$k_{5}$

$M2$ $M3$ $M4$ $M8$ $M9$

Figure 2:  All strongly ergodic models with $N=3$. We show only those models
that are not related to one another by permutations of the states
$1,\ldots,N$. All these models can generate exactly the same phase-type
distribution. Model $M9$ is symmetric and generates the same phase-type
distribution if its parameters are transformed by the permutation
$k_{1}\leftrightarrow k_{2}$, $k_{3}\leftrightarrow k_{4}$ (the inverse
problem has two solutions). Similarly, the inverse problem has two solutions
for $M4$ and $M8$, whereas the solution is unique for $M2$. Models
$M2,M4,M8,M9$ are solvable (the inverse problem has a unique or a finite
number of solutions) and the model $M3$ is not solvable (the inverse problem
has infinitely many solutions). Models generating the same phase-type
distributions spend the same proportion of the total time in state 3. However,
they could be discriminated by the proportion of the total time spent in the
states 1 and 2 (any comparison between $M4$,$M8$,$M9$) or by the lifetimes of
the states 1 and 2 or by both types of times ($M2$ compared to any other).

##  2 Phase-type distribution, direct problem

###  2.1 Finite state continuous time Markov chain model

We consider continuous time Markov chains with N+1 states
$M(t)\in\\{1,\ldots,N+1\\}$. The first $N$ states are not observed. A
supplementary state $N+1$ designates the event that we observe. For
simplicity, we consider that $N+1$ can only be reached from $N$. After $N+1$
there is instantaneous transition to the state $s\in\\{1,\dots,N\\}$.

The chain is defined by its generator $\boldsymbol{Q}$, a $(N+1)\times(N+1)$
matrix with entries $Q_{ij}$ such that $Q_{ij}>0$ represents the probability
per unit time to jump from the state $A_{i}$ to the state $A_{j}$ for $j\neq
i$ and $Q_{ii}=-\sum_{j\neq i}Q_{ij}$ (zero row sum rule).

We are interested in the distribution of $T_{s}$, the first time to reach
$N+1$ starting from the state $s$ with $1\leq s\leq N$:

$T_{s}=\inf\\{t\,|\,t>0,\,M(t)=N+1,\,M(0)=s\\},$

where $M(t)$ is the state of the chain at time $t$.

Let us consider two models. In the model with return, once in $N+1$, the
Markov chain returns to $s$ instantly, $Q_{N+1,s}=\infty$ and then the chain
starts again. Because the reset state $s$ is unique, the observation times in
the model with return form a renewal process (the inter-event times are
positive, independent, identically distributed random variables [19]). Because
after each event, the model is instantaneously reset to $s$ the inter-event
time is always distributed as $T_{s}$.

In the model with absorption, $N+1$ is absorbing and $Q_{N+1,i}=0$ for all
$i\neq N+1$. In this case $T_{s}$ is the time to absorption starting from $s$.

The distribution of $T_{s}$ is the same for the two models, with return and
with absorption. However, the two models are relevant to different
applications. The model with absorption is relevant in scenarios where
observations cease after the event, as seen in cases like failure analysis or
survival analysis during clinical trials. The model with return is applicable
in situations where the event of interest is observed repeatedly, as is the
case in phenomena like transcriptional bursting or epidemiology.

If several return states are possible, then the model with return is no longer
a renewal process. The distribution of the inter-event time depends on the
return state and is no longer a phase-type distribution.

While the return model may not precisely fit the context of a specific
application, we will employ the equivalent absorption model for distribution
calculations due to technical considerations.

For example, with the absorption assumption, the generator of the model $M9$
represented in the Figure 2 is

$\boldsymbol{Q}=\begin{pmatrix}-k_{1}&0&k_{1}&0\\\ 0&-k_{2}&k_{2}&0\\\
k_{3}&k_{4}&-(k_{3}+k_{4}+k_{5})&k_{5}\\\ 0&0&0&0\end{pmatrix}.$ (1)

For the calculation of the distribution of $T_{s}$ it is convenient to
introduce for $1\leq i\leq N+1$ the state probabilities
$X_{i}={\mathbb{P}}\\!\left[M(t)=i\,|\,M(0)=s\right]$.

The variables $X_{i}(t)$ satisfy the following system of linear differential
equations (the master equation):

$\frac{d\boldsymbol{X}}{dt}=\boldsymbol{Q}^{T}\boldsymbol{X},$ (2)

with the initial conditions $X_{n}(0)=\delta_{n,s}$, where $\delta_{n,s}$ is
the Kronecker symbol and $\boldsymbol{Q}^{T}$ is the transpose of the
generator matrix $\boldsymbol{Q}$.

Because the last column of $\boldsymbol{Q}^{T}$ is zero ($N+1$ is absorbing),
the variables $X_{1},\ldots,X_{N}$ satisfy an autonomous ODE system. Indeed,
let $\boldsymbol{\tilde{Q}}$ be the $N\times N$ matrix obtained by eliminating
the last line and the last column of $\boldsymbol{Q}^{T}$. Then
$\boldsymbol{\tilde{X}}=(X_{1},\ldots,X_{N})$ is the solution of

$\frac{d\boldsymbol{\tilde{X}}}{dt}=\boldsymbol{\tilde{Q}}\boldsymbol{\tilde{X}},$
(3)

with initial conditions $X_{i}(0)=\delta_{i,s}$.

The remaining variable $X_{N+1}$ satisfies

$\frac{dX_{N+1}}{dt}=Q_{N,N+1}X_{N},$ (4)

with the initial condition $X_{N+1}(0)=0$.

Let us consider that the eigenvalues $\lambda_{i},1\leq i\leq N$ of the matrix
$\boldsymbol{\tilde{Q}}$ satisfy the following non-degeneracy condition

$\lambda_{i}\neq\lambda_{j},\text{ for all }1\leq i\neq j\leq N,\text{ and
}\max_{1\leq i\leq N}(\lambda_{i})<0.$ (5)

######  Remark 1.

Using properties of solutions of (3) it can be shown that $\max_{1\leq i\leq
N}(\lambda_{i})\leq 0$ is always satisfied. Indeed, $X_{i}(t),1\leq i\leq N$
are probabilities, therefore remain bounded for all $t$, which means that one
cannot have $\lambda_{i}>0$ for some $1\leq i\leq N$. If some $\lambda_{i}>0$,
then there are solutions of (3) that increase without bound when $t$
increases. We will show later that strong connectedness of the transition
graph reduced to the vertices $1,\ldots,N$ implies that none of the
eigenvalues $\lambda_{i}$ can be zero. One should note that, unlike
$\boldsymbol{Q}$, $\boldsymbol{\tilde{Q}}^{T}$ is not a continuous time Markov
chain generator and does not satisfy the zero row sum rule.

Considering that the condition (5) is satisfied, we have

$\boldsymbol{\tilde{X}}=\sum_{i=1}^{N}C_{i}\boldsymbol{u}_{i}e^{\lambda_{i}t}\quad\mathrm{and}\quad
X_{N+1}=Q_{N,N+1}\sum_{i=1}^{N}\frac{C_{i}u_{i,N}}{\lambda_{i}}(e^{\lambda_{i}t}-1),$
(6)

where $\lambda_{i}$ and $\boldsymbol{u}_{i}=(u_{i1},\ldots,u_{i,N})^{T}$ for
$1\leq i\leq N$ are eigenvalues and eigenvectors of $\boldsymbol{\tilde{Q}}$,
respectively.

Because $N+1$ is absorbing,
$X_{N+1}={\mathbb{P}}\\!\left[M(t)={N+1}\,|\,M(0)=s\right]={\mathbb{P}}\\!\left[T_{s}\leq
t\right]$ is the cumulative distribution function of the time $T_{s}$. We also
define the survival function of $T_{s}$ as follows

$S(t)={\mathbb{P}}\\!\left[T_{s}>t\right]=1-X_{N+1}(t)=\sum_{i=1}^{N}A_{i}e^{\lambda_{i}t},$
(7)

where

$A_{i}=-\frac{Q_{N,N+1}C_{i}u_{i,N}}{\lambda_{i}}.$ (8)

The constants $A_{i}$ are not independent. They satisfy

$\sum_{i=1}^{N}A_{i}=1,$ (9)

which follows from $s(0)=1$.

This shows that in the non-degenerate case the distribution of the first
hitting time $T_{s}$ is multi-exponentional with $2N-1$ independent parameters
$\lambda_{1},\ldots,\lambda_{N}$ and $A_{1},\ldots,A_{N-1}$.

##  3 Inverse problem

The inverse problem consists in considering that the survival function and
therefore the parameters $\lambda_{1},\ldots,\lambda_{N}$ and
$A_{1},\ldots,A_{N-1}$ are known. These can be inferred from data using least
squares or maximum likelihood regression.

We want to find the transition rate parameters, which are the non-diagonal
elements of the matrix $\boldsymbol{Q}$. We are interested in models where
this problem could have a unique solution, therefore the matrix
$\boldsymbol{Q}$ has only $2N-1$ non-zero elements. We denote these transition
rates by the $2N-1$ dimensional vector $\boldsymbol{k}$.

We show below that the inverse problem consists in solving $2N-1$ polynomial
equations for the transition rates $\boldsymbol{k}$.

###  3.1 Vieta’s formulas

The eigenvalues $\lambda_{1},\ldots,\lambda_{N}$ are the roots of the the
characteristic polynomial of $\boldsymbol{\tilde{Q}}$, defined as

$P(\lambda)=\mathrm{det}(\boldsymbol{\tilde{Q}}-\lambda\boldsymbol{I})=(-1)^{N}\lambda^{N}+a_{N-1}(\boldsymbol{k})\lambda^{N-1}+\ldots+a_{1}(\boldsymbol{k})\lambda+a_{0}(\boldsymbol{k}).$
(10)

The coefficients $a_{i}(\boldsymbol{k})$ of the characteristic polynomial are
polynomials with integer coefficients on the transition rates
$\boldsymbol{k}$.

A first set of equations relating eigenvalues to the kinetic equations results
from the Vieta’s formulas

$\displaystyle L_{1}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\lambda_{i}=(-1)^{N-1}a_{N-1}(\boldsymbol{k}),$
$\displaystyle L_{2}$ $\displaystyle=$
$\displaystyle\sum_{i<j}\lambda_{i}\lambda_{j}=(-1)^{N-2}a_{N-2}(\boldsymbol{k}),$
$\displaystyle\vdots$ $\displaystyle L_{N}$ $\displaystyle=$
$\displaystyle\lambda_{1}\lambda_{2}\dots\lambda_{N}=a_{0}(\boldsymbol{k}).$
(11)

###  3.2 Eigenvector equations

The amplitude parameters $A_{1},A_{2},\ldots,A_{N}$ of the survival function
occur in (8) together with eigenvector components $u_{i,N}$ and solution
coefficients $C_{i}$. We need to relate the latter to the transition rate
parameters $\boldsymbol{k}$.

First, we need to solve the eigenvector equation

$(\boldsymbol{\tilde{Q}}-\lambda\boldsymbol{I})\boldsymbol{u}=0.$ (12)

From (12) it follows that the eigenvector components are rational functions of
$\lambda$ and $\boldsymbol{k}$.

We consider that the eigenvectors can be chosen of the form
$\boldsymbol{u}(\lambda,\boldsymbol{k})=(u_{1}(\lambda,\boldsymbol{k}),\dots,u_{N}(\lambda,\boldsymbol{k}))$
with $u_{s}(\lambda,\boldsymbol{k})=1$. In the subsection 3.4 we will see for
which models this choice is possible. It follows that

$u_{i}(\lambda,\boldsymbol{k})=\frac{b_{i}(\lambda,\boldsymbol{k})}{b_{s}(\lambda,\boldsymbol{k})},\quad
1\leq i\leq N$

where $b_{i}(\lambda,\boldsymbol{k})$ are polynomials. We also choose
$b_{s}(\lambda,\boldsymbol{k})$ prime relatively to
$b_{i}(\lambda,\boldsymbol{k})$ for all $i\neq s$ in
$\mathbb{Z}[\boldsymbol{k},\lambda]$.

The initial conditions satisfied by the variables $X_{i}$ provide a linear
system of equations for the constants $C_{i}$:

$\sum_{j=1}^{N}u_{i}(\lambda_{j},\boldsymbol{k})C_{j}=\delta_{i,s},\quad 1\leq
i\leq N.$ (13)

Because of the non-degeneracy condition (5), if
$\boldsymbol{u}(\lambda,\boldsymbol{k})$ is not identically zero, then
$\boldsymbol{u}(\lambda_{i},\boldsymbol{k})=\boldsymbol{u}_{i}$ with $1\leq
i\leq N$ are independent. Therefore (13) has a unique solution
$C_{i}(\boldsymbol{\lambda},\boldsymbol{k}),\,1\leq i\leq N$.

From (8) we obtain $N-1$ equations for the transition rates $\boldsymbol{k}$:

$k_{2N-1}u_{N}(\lambda_{i},\boldsymbol{k})C_{i}(\boldsymbol{\lambda},\boldsymbol{k})=-A_{i}\lambda_{i},\quad
1\leq i\leq N-1,$ (14)

where $k_{2N-1}=Q_{N,N+1}$.

The solution of the inverse problem is the solution of the system formed by
(11) and (14). When this solution exists, the transition rates
$\boldsymbol{k}$ can be expressed as functions of the survival function
parameters $\lambda_{i}$ and $A_{i}$ for $1\leq i\leq N$.

###  3.3 Symmetrized systems

There is a difference between (11) and (14). (11) is entirely expressed using
elementary symmetric polynomials in $\lambda_{i}$, whereas there is no obvious
symmetry in (14). We show here that the system formed by (11) and (14) is
equivalent to a system symmetrized in both $\lambda_{i}$ and $A_{i}$. The
advantage of a symmetrized system over a non-symmetrized one is that it
handles simpler formulas, decreasing the computational burden of the symbolic
tools.

To this aim we use the following identities (due to Jacobi-Trudi [20]) that
are valid for any distinct $N$ numbers $\lambda_{i}$ ($1\leq i\leq N$):

$\sum_{i=1}^{N}\lambda_{i}^{k}\prod_{\overset{j=1}{j\neq
i}}^{N}\frac{1}{\lambda_{i}-\lambda_{j}}=\left\\{\begin{array}[]{ll}0,&\text{if
}k<N-1\\\ 1,&\text{if }k=N-1\\\ h_{k-N+1}(\boldsymbol{\lambda}),&\text{if
}k>N-1\end{array}\right.$ (15)

where $h_{k-N+1}(\boldsymbol{\lambda})$ is the complete symmetric polynomial
of degree $k-N+1$ in $N$ variables and
$\boldsymbol{\lambda}=(\lambda_{1},\ldots,\lambda_{N})$. More precisely, the
complete symmetric polynomials in the $N$ variables
$\lambda_{1},\ldots,\lambda_{N}$ have shape

$h_{1}(\boldsymbol{\lambda})=\sum_{1\leq i_{1}\leq N}^{N}\lambda_{i_{1}},\
h_{2}(\boldsymbol{\lambda})=\sum_{1\leq i_{1}\leq i_{2}\leq
N}^{N}\lambda_{i_{1}}\lambda_{i_{2}},\ \dots,\
h_{m}(\boldsymbol{\lambda})=\sum_{1\leq i_{1}\leq\dots\leq i_{m}\leq
N}\lambda_{i_{1}}\lambda_{i_{2}}\cdots\lambda_{i_{m}},\ \dots\ .$ (16)

In order to relate (14) and (15) we first related the coefficients $C_{i}$ to
eigenvalues. This can be done using Cramer’s rule.

Denote by $D_{i}(\lambda,\boldsymbol{k})\in\mathbb{Z}[\lambda,\boldsymbol{k}]$
for $1\leq i\leq N$ the determinant of $(N-1)\times(N-1)$ submatrix of the
matrix $\lambda\boldsymbol{I}-\boldsymbol{\tilde{Q}}$ obtained by deleting its
$s$-th row and $i$-th column.

######  Lemma 2.

Assume that $\lambda_{i}\neq\lambda_{j}$ for $1\leq i\neq j\leq N$ and that
$D_{s}(\lambda_{j},\boldsymbol{k})\neq 0$ and
$u_{s}(\lambda_{j},\boldsymbol{k})=1$ for $1\leq j\leq N$. Then

$C_{j}=\frac{D_{s}(\lambda_{j},\boldsymbol{k})}{\prod_{1\leq l\neq j\leq
N}(\lambda_{j}-\lambda_{l})},\quad 1\leq j\leq N,$ (17)

where $C_{j}$ form the unique solution of the system

$\sum_{1\leq j\leq N}u_{s}(\lambda_{j},\boldsymbol{k})C_{j}=1,$ $\sum_{1\leq
j\leq N}u_{i}(\lambda_{j},\boldsymbol{k})C_{j}=0,\quad 1\leq i\neq s\leq N.$

###### Proof.

If $D_{s}(\lambda_{j},\boldsymbol{k})\neq 0$, then, using Cramer’s rule, we
have

$\boldsymbol{u}(\lambda_{j},\boldsymbol{k})=(D_{1}(\lambda_{j},\boldsymbol{k}),\dots,D_{N}(\lambda_{j},\boldsymbol{k}))/D_{s}(\lambda_{j},\boldsymbol{k}).$

Therefore, it holds

$\sum_{1\leq j\leq N}u_{s}(\lambda_{j},\boldsymbol{k})C_{j}=\sum_{1\leq j\leq
N}\frac{D_{s}(\lambda_{j},\boldsymbol{k})}{\prod_{1\leq l\neq j\leq
N}(\lambda_{j}-\lambda_{l})}=1,$

while

$\sum_{1\leq j\leq N}u_{i}(\lambda_{j},\boldsymbol{k})C_{j}=\sum_{1\leq j\leq
N}\frac{D_{i}(\lambda_{j},\boldsymbol{k})}{\prod_{1\leq l\neq j\leq
N}(\lambda_{j}-\lambda_{l})}=0,\quad 1\leq i\neq s\leq N$

(where $C_{j},1\leq j\leq N$ are taken from (17)) due to (15). ∎

For $1\leq k\leq N-1$ we define $S_{k}$ to be the sum of
$A_{1}\lambda_{1}^{k},\dots,A_{N}\lambda_{N}^{k}$. Using (14) and (17) we
obtain for this sum the formula

$S_{k}=\sum_{i=1}^{N}A_{i}\lambda_{i}^{k}=-k_{2N-1}\sum_{i=1}^{N}\lambda_{i}^{k-1}\frac{D_{s}(\lambda_{i},\boldsymbol{k})}{\prod_{1\leq
j\neq i\leq N}(\lambda_{i}-\lambda_{j})}.$ (18)

From the definition of $D_{i}(\lambda,\boldsymbol{k})$, it follows for $i=s$
that $\mathrm{deg}_{\lambda}(D_{s}(\lambda,\boldsymbol{k}))=N-1$ and that the
(leading) coefficient of $D_{s}(\lambda,\boldsymbol{k}))$ at the monomial
$\lambda^{N-1}$ equals 1, while for $1\leq i\neq s\leq N$ we have
$\mathrm{deg}_{\lambda}(D_{i}(\lambda,\boldsymbol{k}))\leq N-2$. Thus, we have

$\displaystyle D_{s}(\lambda,\boldsymbol{k})$ $\displaystyle=:$
$\displaystyle\lambda^{N-1}+c_{N-2}(\boldsymbol{k})\lambda^{N-2}+\ldots+c_{1}(\boldsymbol{k})\lambda+c_{0}(\boldsymbol{k}),$
(19) $\displaystyle D_{i}(\lambda,\boldsymbol{k})$ $\displaystyle=:$
$\displaystyle
d_{N-2}(\boldsymbol{k})\lambda^{N-2}+d_{N-3}(\boldsymbol{k})\lambda^{N-3}+\ldots+d_{1}(\boldsymbol{k})\lambda+d_{0}(\boldsymbol{k}),\text{
if }i\neq s.$ (20)

We can distinguish two cases:

$\mathbf{s}=N$

Then, according to (19) $D_{N}$ has degree $N-1$. Using (18) and (15) we find
the following $N-1$ symmetrized equations:

$\displaystyle S_{1}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}A_{i}\lambda_{i}=-k_{2N-1},$ $\displaystyle S_{2}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{N}A_{i}\lambda_{i}^{2}=-k_{2N-1}(h_{1}(\boldsymbol{\lambda})+c_{N-2}(\boldsymbol{k})),$
$\displaystyle\vdots$ $\displaystyle S_{N-1}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}A_{i}\lambda_{i}^{N-1}=-k_{2N-1}(h_{N-2}(\boldsymbol{\lambda})+c_{N-2}(\boldsymbol{k})h_{N-3}(\boldsymbol{\lambda})+\ldots+c_{2}(\boldsymbol{k})h_{1}(\boldsymbol{\lambda})+c_{1}(\boldsymbol{k})).$

The symmetrized system is consists of the equations in (11) and (3.3) and is
used to compute the transition rates $\boldsymbol{k}$ as functions of the
symmetric polynomials $L_{1},\ldots,L_{N}$ and $S_{1},\ldots,S_{N-1}$. We
should note that the complete symmetric polynomials
$h_{1}(\boldsymbol{\lambda}),\ldots,h_{N-2}(\boldsymbol{\lambda})$ can be
expressed using the elementary symmetric polynomials $L_{1},\ldots,L_{N-2}$,
that is $h_{1}(\boldsymbol{\lambda})=L_{1}$,
$h_{2}(\boldsymbol{\lambda})=L_{1}^{2}-L_{2}$ and so on.

$\mathbf{s}<N$

Then, according to (20) $D_{N}$ has degree at most $N-2$. Using (18) and (15)
it follows that

$S_{1}=\sum_{i=1}^{N}A_{i}\lambda_{i}=0.$ (22)

This means that the inverse problem is not well posed in this case. The $2N-1$
parameters $\lambda_{1},\ldots,\lambda_{N}$, $A_{1},\ldots,A_{N-1}$ are no
longer independent and they do not determine $\boldsymbol{k}$ uniquely.
Although we can no longer obtain $N-1$ independent symmetrized equations we
can obtain a smaller number of equations

$\displaystyle S_{2}$ $\displaystyle=$ $\displaystyle-
k_{2N-1}d_{N-2}(\boldsymbol{k}),$ $\displaystyle S_{3}$ $\displaystyle=$
$\displaystyle-
k_{2N-1}(d_{N-2}(\boldsymbol{k})h_{1}(\boldsymbol{\lambda})+d_{N-3}(\boldsymbol{k})),$
$\displaystyle\vdots$ $\displaystyle S_{k}$ $\displaystyle=$ $\displaystyle-
k_{2N-1}(d_{N-2}(\boldsymbol{k})h_{k-2}(\boldsymbol{\lambda})+d_{N-3}(\boldsymbol{k})h_{k-3}(\boldsymbol{\lambda})+\ldots+d_{N-k}(\boldsymbol{k})),$
$\displaystyle\vdots$ $\displaystyle S_{N-1}$ $\displaystyle=$ $\displaystyle-
k_{2N-1}(d_{N-2}(\boldsymbol{k})h_{N-3}(\boldsymbol{\lambda})+\ldots+d_{2}(\boldsymbol{k})h_{1}(\boldsymbol{\lambda})+d_{1}(\boldsymbol{k})).$
(23)

In this case the symmetrized system made of (11),(3.3) can be used to compute
the transition rates $\boldsymbol{k}$ as functions of the symmetric
polynomials $L_{1},\ldots,L_{N}$ and $S_{2},\ldots,S_{N-1}$ and of one or
several indeterminate transition rates.

###  3.4 Solvable models

Different Markov chain models are distinguished by their transition graph
defined as the directed graph $G=(V,A)$ with vertices $V=\\{1,\ldots,N+1\\}$
and arcs $A=\\{(i,j)\mid i\neq j,\,Q_{ij}\neq 0\\}$.

Let us consider that this graph satisfies some conditions related to the
underlying modeled process:

  1. (C1)

$N+1$ is reachable only from $N$.

  2. (C2)

From $N+1$ there is only one outgoing arc, towards $s$, where $1\leq s\leq N$.

Conditions C1 and C2 where used for setting the problem and writing the
equations defining the inverse problem. C2 is needed in order to have a phase-
type distribution. C1 is a simplifying assumption that is satisfied in a
number of applications, including the biological application discussed in this
paper.

We would like to know for which models satisfying C1 and C2 the inverse
problem has a unique solution, eventually up to transformation by discrete
symmetries.

######  Definition 3.

We say that a model is solvable if the following solutions are satisfied:

  1. (i)

The equation $(\boldsymbol{\tilde{Q}}-\lambda\boldsymbol{I})\boldsymbol{u}=0$
has solutions
$\boldsymbol{u}(\lambda,\boldsymbol{k})=(u_{1}(\lambda,\boldsymbol{k}),\dots,u_{N}(\lambda,\boldsymbol{k}))$
with $u_{s}(\lambda,\boldsymbol{k})$ not identically zero where $s$ is the
return state.

  2. (ii)

The model has $2N-1$ transition rate parameters.

  3. (iii)

The system made by the equations (11) and (14) or equivalently the system made
by the equations (11) and (3.3) has unique solutions up to transformations by
discrete symmetries, on a open domain of dimension $2N-1$.

It is very difficult to obtain general sufficient conditions for solvability
but we can state a number of necessary conditions.

The following proposition follows from the proof of Lemma 2 in the subsection
3.3.

######  Proposition 4.

A solvable model necessarily satisfies $s=N$.

For a directed graph $G$ and its vertices $v,w$ we denote $v\preceq w$ if
there is a path in $G$ from $v$ to $w$. We say that $v,w$ are equivalent is
$v\preceq w,w\preceq v$. If $G$ consists of a single equivalence class then we
call $G$ strongly connected. A Markov chain model is called strongly ergodic
if and only if its transition graph is strongly connected.

The following condition is another necessary condition of solvability.

######  Proposition 5.

A solvable model is strongly ergodic.

###### Proof.

We will prove that a model with a non strongly connected transition graph
contradicts (i) of Definition 3.

Denote by $G$ the graph of the Markov chain. Suppose the $G$ is not ergodic.
There exists its vertex $v$ such that not for every vertex $w$ of $G$ it holds
$w\preceq v$. Reorder the vertices of $G$ in such a way that first are listed
all $r$ vertices $w$ for which $w\preceq v$ holds, and $v$ is listed as the
last one among them (thus, $v$ is numbered by $r$ in the ordering). The rest
of $N-r$ vertices are listed in an arbitrary order. The resulting matrix of
the Markov chain we still denote by $\boldsymbol{\tilde{Q}}$ (slightly abusing
the notations).

Making use of the notations from Lemma 2 we claim that the polynomial
$D_{r+1}(\lambda,\boldsymbol{k})=0$ vanishes identically. Indeed, in the
matrix obtained from $\boldsymbol{\tilde{Q}}$ by means of deleting its $r$-th
row and $(r+1)$-th column, all the entries in places $(i,j)$ where $i\geq
r+1,j\leq r$ vanish, which justifies the claim.

On the other hand $\mathrm{deg}_{\lambda}(D_{r}(\lambda,\boldsymbol{k}))=N-1$,
and the coefficient of $D_{r}(\lambda,\boldsymbol{k})$ at $\lambda^{N-1}$
equals 1 (cf. the proof of Lemma 2). Hence
$D_{r}(\lambda_{t},\boldsymbol{k})\neq 0$ for a suitable $1\leq t\leq N$.
Therefore
$u_{r+1}(\lambda_{t},\boldsymbol{k})=(D_{r+1}/D_{r})(\lambda_{t},\boldsymbol{k})=0$
(see again the proof of Lemma 2). ∎

Reciprocally, we have

######  Proposition 6.

A strongly ergodic model satisfying the conditions $\mathrm{C1}$,
$\mathrm{C2}$ and $s=N$ satisfies the condition (i) of Definition 3, in other
words the symmetrized equations (3.3) can be written down for such models.

###### Proof.

If $s=N$, then
$D_{N}(\lambda,\boldsymbol{k})\in\mathbb{Z}[\lambda,\boldsymbol{k}]$ is the
determinant of $(N-1)\times(N-1)$ submatrix of the matrix
$\lambda\boldsymbol{I}-\boldsymbol{\tilde{Q}}$ obtained by deleting its $N$-th
row and $N$-th column. This is a polynomial in $\lambda$, whose leading term
is $\lambda^{N-1}$ and therefore not identically zero. In order to prove (i)
it is then enough to follow the proof of Lemma 2. ∎

Strongly ergodic models also satisfy the following property that has been used
in condition (5) and is needed for writing the solution (6).

######  Proposition 7.

All strongly ergodic models satisfy $\max_{1\leq i\leq N}\lambda_{i}<0$, where
$\lambda_{i}$ are the eigenvalues of $\boldsymbol{\tilde{Q}}$.

###### Proof.

Suppose that some $\lambda_{i}=0$. Then the corresponding eigenvector
satisfies
$\frac{d\boldsymbol{X}_{i}}{dt}=\boldsymbol{\tilde{Q}}\boldsymbol{X}_{i}=0$.
Using the condition C1 we find that
$\frac{d\left(\sum_{j=1}^{N}X_{ij}\right)}{dt}=-Q_{N,N+1}X_{iN}$. It follows
that $X_{iN}=0$. Also, $\boldsymbol{X}_{i}$ is a steady state distribution of
the Markov chain having the generator $\boldsymbol{\tilde{Q}}^{T}$ in which
$Q_{N,N+1}$ is set to zero. This chain is strongly ergodic. Or, any strongly
ergodic Markov chain has a unique steady state distribution in which all
states have non-zero probabilities, which contradicts $X_{iN}=0$. ∎

Figure 2 shows all strongly ergodic models with $N=3$ satisfying the
conditions C1 and C2.

##  4 Solution of the inverse problem for the unbranched chain model

In this section we present an example of a solvable model with an arbitrary
number of states. For such model we propose a method to compute the solutions
of the inverse problem.

Consider now a Markov chain with $N+1$ states located on a line and reversibly
connected one to next such that

${\tilde{Q}}_{n+1,n}=k_{n}^{+},{\tilde{Q}}_{n,n+1}=k_{n}^{-},1\leq
n<N,{\tilde{Q}}_{n,n}=-k_{n}^{+}-k_{n-1}^{-},1\leq
n<N,{\tilde{Q}}_{N,N}=-k_{N}-k_{N-1}^{-}.$

This model is a generalization, of arbitrary length $N$, of the model $M2$
represented in Figure 2 that has $N=3$.

We have an algebraic map

$f:=f_{N}:\mathbb{C}^{2N-1}\to\mathbb{C}^{2N-1},\,f(k_{1}^{+},\dots,k_{N-1}^{+},k_{1}^{-},\dots,k_{N-1}^{-},k_{N})=(A_{1},\dots,A_{N-1},\lambda_{1},\dots,\lambda_{N}).$
(24)

The goal of this section is to prove that $f$ is invertible and that its
inverse is a rational function. In other words, we will show that the
unbranched chain model is solvable for any $N$.

Note that (13), (14) imply that

$C_{1}+\cdots+C_{N}=A_{1}+\cdots+A_{N}=1,\,k_{N}=-\sum_{1\leq j\leq
N}A_{j}\lambda_{j}.$ (25)

Denote a vector $C:=(C_{1},\dots,C_{N})^{T}$.

For each eigenvalue $\lambda$ of the matrix $\bf\tilde{Q}$ it holds

$k_{N-1}^{+}u_{N-1}(\lambda,\boldsymbol{k})=k_{N}+k_{N-1}^{-}+\lambda,\,k_{n}^{+}u_{n}(\lambda,\boldsymbol{k})=(k_{n+1}^{+}+k_{n}^{-}+\lambda)u_{n+1}(\lambda,\boldsymbol{k})-k_{n+1}^{-}u_{n+2}(\lambda,\boldsymbol{k})$

for $0\leq n\leq N-2$. This provides by recursion on $N-n$ a rational function
$g_{n}$ with rational coefficients such that

$k_{n}^{+}u_{n}(\lambda,\boldsymbol{k})=g_{n}(k_{n+1}^{+},\dots,k_{N-1}^{+},k_{n}^{-},\dots,k_{N-1}^{-},k_{N},\lambda).$
(26)

Moreover, the denominator of $g_{n}$ equals $k_{n+1}^{+}\cdots k_{N-1}^{+}$.

Consider $N\times N$ matrix $U$ with the columns
$\boldsymbol{u}(\lambda_{1},\boldsymbol{k}),\dots,\boldsymbol{u}(\lambda_{N},\boldsymbol{k})$.
Then for $s\geq 0$ it holds

${\bf\tilde{Q}}^{s}UC=U\cdot\mathrm{diag}(\lambda_{1}^{s},\dots,\lambda_{N}^{s})C={\bf\tilde{Q}}^{s}(0,\dots,0,1)^{T}$
(27)

where $\mathrm{diag}(\lambda_{1}^{s},\dots,\lambda_{N}^{s})$ denotes a
diagonal matrix.

The $(N-s)$-th coordinate of the middle vector in (27) equals

$\sum_{1\leq j\leq
N}u_{N-s}(\lambda_{j},\boldsymbol{k})\lambda_{j}^{s}(-A_{j}\lambda_{j}/k_{N}).$
(28)

The same $(N-s)$-th coordinate of the right vector in (27) equals

$k_{N-s}^{-}\cdots k_{N-1}^{-}.$ (29)

Observe that the $j$-th coordinate of the right vector in (27) vanishes for
$1\leq j<N-s$.

Multiplying the both sides of (26) for $n=N-s,\lambda=\lambda_{j}$ by
$\lambda_{j}^{s}(-A_{j}\lambda_{j}/k_{N})$ and summing them up over $1\leq
j\leq N$, we obtain that

$k_{N-s}^{+}k_{N-s}^{-}\cdots
k_{N-1}^{-}=G_{N-s,0}(A_{1},\dots,A_{N-1},k_{N-s+1}^{+},\dots,k_{N-1}^{+},k_{N-s}^{-},\dots,k_{N-1}^{-},k_{N},\lambda_{1},\dots,\lambda_{N})$
(30)

for a suitable rational function $G_{N-s,0}$ with rational coefficients and
with a denominator

$k_{N-s+1}^{+}\cdots k_{N-1}^{+}k_{N}$

taking into account the equality of (28) and of (29). Summarizing, we have
established the following statement.

######  Lemma 8.

For a suitable rational function $G_{N-s}$ with rational coefficients and with
a denominator $k_{N-s+1}^{+}\cdots k_{N-1}^{+}k_{N-s}^{-}\cdots
k_{N-1}^{-}k_{N}$ it holds

$k_{N-s}^{+}=G_{N-s}(A_{1},\dots,A_{N-1},k_{N-s+1}^{+},\dots,k_{N-1}^{+}k_{N-s}^{-},\dots,k_{N-1}^{-},k_{N},\lambda_{1},\dots,\lambda_{N}).$

The $(N-s+1)$-th coordinate of the middle vector in (27) equals

$\sum_{1\leq j\leq
N}u_{N-s+1}(\lambda_{j},\boldsymbol{k})\lambda_{j}^{s}(-A_{j}\lambda_{j}/k_{N}).$
(31)

The same $(N-s+1)$-th coordinate of the right vector in (27) equals

$-k_{N-s}^{-}\cdots
k_{N-1}^{-}+d_{N-s}(k_{N-s+1}^{+},\dots,k_{N-1}^{+},k_{N-s+1}^{-},\dots,k_{N-1}^{-},k_{N})$
(32)

for an appropriate polynomial $h_{N-s}$ with integer coefficients.

Multiplying the both sides of (26) for $n=N-s+1,\lambda=\lambda_{j}$ by
$\lambda_{j}^{s}(-A_{j}\lambda_{j}/k_{N})$ and summing them up over $1\leq
j\leq N$, we obtain that

$-k_{N-s+1}^{+}k_{N-s}^{-}\cdots
k_{N-1}^{-}+d_{N-s}(k_{N-s+1}^{+},\dots,k_{N-1}^{+},k_{N-s+1}^{-},\dots,k_{N-1}^{-},k_{N})=$
$D_{N-s,0}(A_{1},\dots,A_{N-1},k_{N-s+2}^{+},\dots,k_{N-1}^{+},k_{N-s+1}^{-},\dots,k_{N-1}^{-},k_{N},\lambda_{1},\dots,\lambda_{N})$

for a suitable rational function $D_{N-s,0}$ with rational coefficients and
with a denominator $k_{N-s+2}^{+}\cdots k_{N-1}^{+}k_{N}$ taking into account
the equality of (31) and of (32) (while the both latter multiplied by
$k_{N-s+1}^{+}$). Summarizing, we have established the following statement.

######  Lemma 9.

For an appropriate rational function $D_{N-s}$ with rational coefficients and
with a denominator $k_{N-s+1}^{+}\cdots k_{N-1}^{+}k_{N-s+1}^{-}\cdots
k_{N-1}^{-}k_{N}$ it holds

$k_{N-s}^{-}=D_{N-s}(A_{1},\dots,A_{N-1},k_{N-s+1}^{+},\dots,k_{N-1}^{+},k_{N-s+1}^{-},\dots,k_{N-1}^{-},k_{N},\lambda_{1},\dots,\lambda_{N}).$

Applying alternatingly Lemma 9 and Lemma 8 for $s=1,\dots,N-1$ consecutively,
we conclude with the following main result of this section.

######  Theorem 10.

One can produce rational functions $F_{n}^{+},F_{n}^{-},1\leq n\leq N-1$ with
rational coefficients such that

$k_{n}^{+}=F_{n}^{+}(A_{1},\dots,A_{N-1},\lambda_{1},\dots,\lambda_{N}),k_{n}^{-}=F_{n}^{-}(A_{1},\dots,A_{N-1},\lambda_{1},\dots,\lambda_{N}).$

######  Remark 11.

  1. (i)

Together with (25) Theorem 10 assures the inverse rational map to $f$ (see
(24)) on the open dense subset of $\mathbb{C}^{2N-1}$ determined by conditions
$k_{1}^{+}\cdots k_{N-1}^{+}\cdot k_{1}^{-}\cdots k_{N-1}^{-}\cdot k_{N}\neq
0$ and $\lambda_{i}\neq\lambda_{j},1\leq i<j\leq N$;

  2. (ii)

the proof of Theorem 10 provides an algorithm which produces explicitly
rational functions $F_{n}^{+},F_{n}^{-}$ by recursion on $N-n$.

######  Remark 12.

In fact, for any model (not necessary, the unbranched chain model elaborated
in this section) one can consider a similar to (24) rational map. Thus,
$A_{1},\dots,A_{N-1},\lambda_{1},\dots,\lambda_{N-1}$ are rational functions
in $k_{n}^{+},k_{n}^{-},1\leq n<N,k_{N}$. Therefore, there is a fields
extension

$\mathbb{C}(A_{1},\dots,A_{N-1},\lambda_{1},\dots,\lambda_{N-1})\subset\mathbb{C}(k_{1}^{+},\dots,k_{n-1}^{+},k_{1}^{-},\dots,k_{n-1}^{-},k_{N}).$

It is known (see e.g. [31], Ch. 1) that the degree of this extension equals
the number of solutions in $\mathbb{C}^{2N-1}$ of the system of rational
equations
$A_{1}=\alpha_{1},\dots,A_{N-1}=\alpha_{N-1},\lambda_{1}=\beta_{1},\dots,\lambda_{N}=\beta_{N}$
at a generic point
$(\alpha_{1},\dots,\alpha_{N-1},\beta_{1},\dots,\beta_{N})\in\mathbb{C}^{2N-1}$.
Recall that the degree is defined as the dimension of the vector space
$\mathbb{C}(k_{1}^{+},\dots,k_{n-1}^{+},k_{1}^{-},\dots,k_{n-1}^{-},k_{N})$
over the field
$\mathbb{C}(A_{1},\dots,A_{N-1},\lambda_{1},\dots,\lambda_{N-1})$. The degree
can be infinite. Theorem 10 states that for the unbranched chain model the
degree equals 1.

We conjecture that for all other models the degree is greater than 1.

##  5 Using the Thomas decomposition for solving the inverse problem

Thomas decomposition is a computational algebra algorithm allowing to
decompose systems of polynomial equations and inequations into simple systems
whose solutions can be more easily found by iteratively solving unvariate
polynomial equations. We apply this technique to compute the solutions of the
inverse problem for all the strongly ergodic models with $N=3$ and $N=4$.

###  5.1 The Thomas decomposition of an algebraic system

In this section we introduce briefly the notion of the algebraic Thomas
decomposition (see the appendix of [24] for a similar introduction). We will
then apply the algebraic Thomas decomposition in the subsequent sections to
our symmetrized systems to determine their solutions.

Let $\mathbb{C}[\boldsymbol{x}]$ be a polynomial ring in $n$ variables
$\boldsymbol{x}=(x_{1},\dots,x_{n})$ over the complex numbers $\mathbb{C}$. An
algebraic system $\mathcal{S}$ is defined as a finite set of polynomial
equations and inequations, that is as the set

$\mathcal{S}=\\{p_{1}(\boldsymbol{x})=0,\,\dots,\,p_{r}(\boldsymbol{x})=0,\,q_{1}(\boldsymbol{x})\neq
0,\,\dots,\,q_{s}(\boldsymbol{x})\neq 0\\}$ (33)

with polynomials $p_{i}(\boldsymbol{x})$, $q_{j}(\boldsymbol{x})$ in
$\mathbb{C}[\boldsymbol{x}]$ and integers $r$, $s\in\mathbb{N}_{0}$. The
solution set $\mathrm{Sol}(\mathcal{S})$ of the algebraic system (33) is
defined as the set of all
$\overline{\boldsymbol{x}}=(\overline{x}_{1},\dots,\overline{x}_{n})\in\mathbb{C}^{n}$
satisfying the equations and inequations of $\mathcal{S}$, that is as

$\mathrm{Sol}(\mathcal{S})=\\{\overline{\boldsymbol{x}}\in\mathbb{C}^{n}\mid
p_{i}(\overline{\boldsymbol{x}})=0,\,q_{j}(\overline{\boldsymbol{x}})\neq 0\
\mathrm{for}\ \mathrm{all}\ 1\leq i\leq r\ \mathrm{and}\ 1\leq j\leq s\\}.$

Geometrically, $\mathrm{Sol}(\mathcal{S})$ is the difference of the two
varieties

$\\{\overline{\boldsymbol{x}}\in\mathbb{C}^{n}\mid
p_{1}(\overline{\boldsymbol{x}})=0,\,\dots,\,p_{r}(\overline{\boldsymbol{x}})=0\\}\
\mathrm{and}\ \\{\overline{\boldsymbol{x}}\in\mathbb{C}^{n}\mid
q_{1}(\overline{\boldsymbol{x}})\cdots q_{s}(\overline{\boldsymbol{x}})=0\\}$

and so it is a locally Zariski closed subset of $\mathbb{C}^{n}$.

In order to introduce the notion of an algebraic Thomas decomposition of a
system $\mathcal{S}$ we make the following definitions. On the variables
$\boldsymbol{x}=(x_{1},\dots,x_{n})$ of our polynomial ring
$\mathbb{C}[\boldsymbol{x}]$ we define a total ordering (sometimes also called
a ranking) by setting $x_{i}<x_{j}$ for $i<j$. With respect to this ranking
the leader $\mathrm{ld}(p(\boldsymbol{x}))$ of a non-constant polynomial
$p(\boldsymbol{x})$ is defined as the greatest variable appearing in
$p(\boldsymbol{x})$. In case $p(\boldsymbol{x})\in\mathbb{C}$ is a constant
polynomial, we set $\mathrm{ld}(p(\boldsymbol{x}))=1$. If we consider every
polynomial $p(\boldsymbol{x})\in\mathbb{C}[\boldsymbol{x}]$ as a univariate
polynomial in its leader, say $\mathrm{ld}(p(\boldsymbol{x}))=x_{k}$, then the
coefficients of $p(\boldsymbol{x})$ as a polynomial in $x_{k}$ are polynomials
in $\mathbb{C}[x_{1},\dots,x_{k-1}]$. The coefficient of the highest power of
$\mathrm{ld}(p(\boldsymbol{x}))$ in $p(\boldsymbol{x})$ is called the initial
of $p(\boldsymbol{x})$ which we denote by $\mathrm{init}(p(\boldsymbol{x}))$.
The separant $\mathrm{sep}(p(\boldsymbol{x}))$ of a polynomial
$p(\boldsymbol{x})$ is defined as the partial derivative of
$p(\boldsymbol{x})$ with respect to its leader.

######  Definition 13.

Let $\mathcal{S}$ be the algebraic system of (33). Then $\mathcal{S}$ is
called a simple algebraic system with respect to a ranking, if the following
conditions are satisfied:

  1. 1.

The leaders of all equations and inequations are pairwise different, i.e. we
have

$\mathrm{card}\big{(}\\{\mathrm{ld}(p_{1}(\boldsymbol{x})),\,\dots,\mathrm{ld}(p_{r}(\boldsymbol{x})),\,\mathrm{ld}(q_{1}(\boldsymbol{x})),\,\dots,\,\mathrm{ld}(q_{s}(\boldsymbol{x}))\\}\setminus\\{1\\}\big{)}=r+s.$

This property is called triangularity.

  2. 2.

For every
$p(\boldsymbol{x})\in\\{p_{1}(\boldsymbol{x}),\,\dots,\,p_{r}(\boldsymbol{x}),\,q_{1}(\boldsymbol{x}),\,\dots,\,q_{s}(\boldsymbol{x})\\}$
the equation $\mathrm{init}(p(\boldsymbol{x}))=0$ has no solution in
$\mathrm{Sol}(\mathcal{S})$. We call this property non-vanishing initials.

  3. 3.

For every
$p(\boldsymbol{x})\in\\{p_{1}(\boldsymbol{x}),\,\dots,\,p_{r}(\boldsymbol{x}),\,q_{1}(\boldsymbol{x}),\,\dots,\,q_{s}(\boldsymbol{x})\\}$
the equation $\mathrm{sep}(p(\boldsymbol{x}))=0$ has no solution in
$\mathrm{Sol}(\mathcal{S})$. This is called square-freeness.

The advantage of a simple algebraic system $\mathcal{S}$ is that one can
obtain its solution set by iteratively solving univariate polynomials. This is
a consequence of the triangularity of a simple algebraic system. Indeed, the
triangularity implies that there is at most one equation $p(\boldsymbol{x})=0$
with leader $x_{1}$ or at most one inequation $q(\boldsymbol{x})\neq 0$ with
leader $x_{1}$. Note that $p(\boldsymbol{x})$ or respectively
$q(\boldsymbol{x})$ is a univariate polynomial in $x_{1}$. The square-freeness
implies that the number of zeros in $\mathbb{C}$ of $p(\boldsymbol{x})$
(respectively $q(\boldsymbol{x})$) is equal the degree of $p(\boldsymbol{x})$
(respectively $q(\boldsymbol{x})$). In case of the equation
$p(\boldsymbol{x})=0$, any root $\overline{x}_{1}\in\mathbb{C}$ of
$p(\boldsymbol{x})$ can be chosen as the first coordinate of a solution
$\overline{\boldsymbol{x}}$ of $\mathcal{S}$. In case of the inequation
$q(\boldsymbol{x})\neq 0$, all elements of $\mathbb{C}$ except for the roots
of $q(\boldsymbol{x})$ can here be chosen as the first coordinate of a
solution. If there is no equation or inequation with leader $x_{1}$, then the
first coordinate is free, that is $\overline{x}_{1}$ can be chosen arbitrary
in $\mathbb{C}$. Now we make the first iteration step. Again by triangularity
there is at most one equation or inequation with leader $x_{2}$. If there is
an equation or inequation with leader $x_{2}$, then we substitute
$\overline{x}_{1}$ for $x_{1}$ in this equation or inequation and obtain so an
univariate polynomial in $x_{2}$. Condition 2 of Definition 13 guarantees that
the degree of the so obtained polynomial is independent of the choice of
$\overline{x}_{1}$ and the square-freeness implies that the number of roots of
this polynomial is equal to its degree. According to the three possible cases,
that is there is an equation, inequation or neither of them, we determine as
described above $\overline{x}_{2}\in\mathbb{C}$ for the second coordinate of a
solution of $\mathcal{S}$. An iteration of this process yields successively a
solution
$\overline{\boldsymbol{x}}=(\overline{x}_{1},\dots,\overline{x}_{n})\in\mathbb{C}^{n}$
of $\mathcal{S}$. Moreover, any solution of the simple algebraic system
$\mathcal{S}$ can be obtained by this process.

######  Definition 14.

A Thomas decomposition of an algebraic system $\mathcal{S}$ as in (33)
consists of finitely many simple algebraic systems
$\mathcal{S}_{1},\dots,\mathcal{S}_{m}$ such that $\mathrm{Sol}(\mathcal{S})$
is the disjoint union of
$\mathrm{Sol}(\mathcal{S}_{1}),\dots,\mathrm{Sol}(\mathcal{S}_{m})$.

It was proved by Thomas in [36, 35] that any algebraic system has a Thomas
decomposition which is in general not unique. A Thomas decomposition can be
determined algorithmically (see [4]) and there is an implementation in Maple.
A description of the implementation can be found in [5, 21].

We will compute in the subsequent sections a Thomas decomposition of the
symmetrized systems $M2$, $M3$, $M4$, $M8$ and $M9$ ($N=3$) by applying the
Maple implementation to them. To this end we need to define a ranking on the
polynomial ring

$\mathbb{C}[k_{1},k_{2},k_{3},k_{4},k_{5},L_{1},L_{2},L_{3},S_{1},S_{2}].$

To simplify notation we collect the variables into
$\boldsymbol{k}=(k_{1},k_{2},k_{3},k_{4},k_{5})$ and
$\boldsymbol{v}=(L_{1},L_{2},L_{3},S_{1},S_{2})$. Since we want to solve the
symmetrized systems for the variables $\boldsymbol{k}$, we rank them always
higher than the collection of variables $\boldsymbol{v}$. Among the variables
$\boldsymbol{k}$ we sometimes change the ranking for the different symmetrized
systems. This is done to minimize the output of the Thomas decomposition and
has no deeper meeting. The ranking of the variables
$\boldsymbol{k}>\boldsymbol{v}$ implies that each simple algebraic system
returned by the Thomas decomposition yields solutions for the variables
$\boldsymbol{k}$ which are only valid for solutions of the variables
$\boldsymbol{v}$ satisfying equation and inequation conditions in
$\boldsymbol{v}$ over $\mathbb{C}$.

###  5.2 A Thomas decomposition for model 2

The model studied here is a particular case of the one from Section 4 and
Theorem 10. We are going to determine the solutions of the algebraic system of
equations

$\begin{array}[]{rcl}\mathcal{S}&=&\big{\\{}-S_{1}=k_{5},\
-S_{2}=k_{5}\,(L_{1}+k_{1}+k_{2}+k_{3}),\
L_{1}=-k_{1}-k_{2}-k_{3}-k_{4}-k_{5},\\\\[5.0pt] &&\
L_{2}=k_{1}\,k_{4}+k_{1}\,k_{5}+k_{2}\,k_{3}+k_{2}\,k_{5}+k_{3}\,k_{4}+k_{3}\,k_{5},\
L_{3}=-k_{2}\,k_{3}\,k_{5}\big{\\}}.\end{array}$

To this end we compute with Maple an algebraic Thomas decomposition of
$\mathcal{S}$ with respect to the ranking

$k_{1}>k_{3}>k_{2}>k_{4}>k_{5}>\boldsymbol{v}.$

Maple returns eleven simple systems $\mathcal{S}_{1},\dots,\mathcal{S}_{11}$.
We will only compute here the solutions in $\boldsymbol{k}$ for the first
simple systems, since it the most generic one, meaning that the solutions are
valid for specializations of the variables $\boldsymbol{v}$ satisfying only
inequations, i.e the solutions for $\boldsymbol{k}$ are valid over an Zariski
open subset of $\mathbb{C}^{5}$. The first simple system consists of the
equations and inequations

$\begin{array}[]{rcl}\mathcal{S}_{1}&=&\big{\\{}L_{1}^{2}\,S_{1}^{3}\,S_{2}+L_{2}^{2}\,S_{1}^{3}+S_{1}\,S_{2}^{3}+L_{3}^{2}\,S_{1}+(-S_{1}^{2}\,S_{2}^{2}+S_{2}^{3}+(S_{1}^{3}\,S_{2}-S_{1}\,S_{2}^{2})\,L_{1}\\\\[5.0pt]
&&\
+(-S_{1}^{4}+S_{1}^{2}\,S_{2})\,L_{2}+(S_{1}^{3}-S_{1}\,S_{2})\,L_{3})\,k_{1}+(-2\,S_{1}^{2}\,S_{2}^{2}+(-S_{1}^{4}-S_{1}^{2}\,S_{2})\,L_{2}\\\\[5.0pt]
&&\
+(S_{1}^{3}+S_{1}\,S_{2})\,L_{3})\,L_{1}+(S_{1}^{3}\,S_{2}-2\,L_{3}\,S_{1}^{2}+S_{1}\,S_{2}^{2})\,L_{2}+(S_{1}^{4}-3\,S_{1}^{2}\,S_{2})\,L_{3}=0,\\\\[5.0pt]
&&\
(L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}+L_{3}\,S_{1}-S_{2}^{2})\,k_{3}+(-S_{1}^{2}+S_{2})\,L_{3}=0,\\\\[5.0pt]
&&\
-L_{1}\,S_{1}\,S_{2}+L_{2}\,S_{1}^{2}-L_{3}\,S_{1}+S_{2}^{2}+(S_{1}^{3}-S_{1}\,S_{2})\,k_{2}=0,\
k_{4}\,S_{1}-S_{1}^{2}+S_{2}=0,\\\\[5.0pt] &&\ k_{5}+S_{1}=0,\
L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}+L_{3}\,S_{1}-S_{2}^{2}\neq 0,\
S_{1}^{3}-S_{1}\,S_{2}\neq 0,\ S_{2}\neq 0\big{\\}}.\end{array}$

One easily checks that the last three inequations do not involve any variable
$\boldsymbol{k}$. Thus the inequations

$L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}+L_{3}\,S_{1}-S_{2}^{2}\neq 0,\quad
S_{1}^{3}-S_{1}\,S_{2}\neq 0,\quad S_{2}\neq 0$ (34)

define the above described Zariski open subset. Solving now successively the
remaining equations for the variables $\boldsymbol{k}$ (they are all linear in
their respective leaders), we obtain the solutions

$\displaystyle k_{1}$
$\displaystyle=\frac{-1}{(-S_{1}^{2}\,S_{2}^{2}+S_{2}^{3}+(S_{1}^{3}\,S_{2}-S_{1}\,S_{2}^{2})\,L_{1}+(-S_{1}^{4}+S_{1}^{2}\,S_{2})\,L_{2}+(S_{1}^{3}-S_{1}\,S_{2})\,L_{3})}$
$\displaystyle(L_{1}^{2}\,S_{1}^{3}\,S_{2}+L_{2}^{2}\,S_{1}^{3}+S_{1}\,S_{2}^{3}+L_{3}^{2}\,S_{1}+(-2\,S_{1}^{2}\,S_{2}^{2}+(-S_{1}^{4}-S_{1}^{2}\,S_{2})\,L_{2}+(S_{1}^{3}+S_{1}\,S_{2})\,L_{3})\,L_{1}$
$\displaystyle+(S_{1}^{3}\,S_{2}-2\,L_{3}\,S_{1}^{2}+S_{1}\,S_{2}^{2})\,L_{2}+(S_{1}^{4}-3\,S_{1}^{2}\,S_{2})\,L_{3}),$
$\displaystyle k_{3}$
$\displaystyle=\frac{(S_{1}^{2}-S_{2})\,L_{3}}{(L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}+L_{3}\,S_{1}-S_{2}^{2})},$
$\displaystyle k_{2}$
$\displaystyle=\frac{L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}+L_{3}\,S_{1}-S_{2}^{2}}{S_{1}^{3}-S_{1}\,S_{2}},$
$\displaystyle k_{4}$ $\displaystyle=\frac{S_{1}^{2}-S_{2}}{S_{1}},$
$\displaystyle k_{5}$ $\displaystyle=-S_{1},$

which are only valid for those $\overline{\boldsymbol{v}}\in\mathbb{C}^{5}$
satisfying the inequations in (34). If one is now interested in a solutions
for $\boldsymbol{k}$ with respect to a
$\overline{\boldsymbol{v}}\in\mathbb{C}^{5}$ which does not satisfy the
inequations of (34), then there is a simple system among
$\mathcal{S}_{2},\dots,\mathcal{S}_{11}$, which one can solve successively for
$\boldsymbol{k}$ as described in the previous subsection. The remaining simple
systems $\mathcal{S}_{2},\dots,\mathcal{S}_{11}$ can be found in subsection
A.1 of the appendix.

The conditions (34) do not contain equalities, therefore the inverse problem
has solutions on an open domain of dimension 5. According to the Section 3.4
this means that the model $M2$ is solvable.

###  5.3 A Thomas decomposition for model 3

The symmetrized algebraic system for model 3 is

$\begin{array}[]{rcl}\mathcal{S}&=&\big{\\{}-S_{1}=k_{5},\
-S_{2}=k_{5}\,(L_{1}+k_{1}+k_{2}+k_{3}),\
L_{1}=-k_{1}-k_{2}-k_{3}-k_{4}-k_{5},\\\\[5.0pt] &&\
L_{3}=-k_{2}\,k_{3}\,k_{5},\
L_{2}=k_{1}\,k_{4}+k_{1}\,k_{5}+k_{2}\,k_{3}+k_{2}\,k_{4}+k_{2}\,k_{5}+k_{3}\,k_{4}+k_{3}\,k_{5}\big{\\}}\end{array}$

and we compute an algebraic Thomas decomposition for it with respect to the
ranking

$k_{1}>k_{2}>k_{3}>k_{4}>k_{5}>\boldsymbol{v}.$

Maple returns $10$ simple systems $\mathcal{S}_{1},\dots,\mathcal{S}_{10}$.
One easily checks by comparing the number of equations only involving the
variables $\boldsymbol{v}$ (see appendix A.2 for the remaining simple systems
$\mathcal{S}_{2},\dots,\mathcal{S}_{10}$), that the first simple sytem

$\begin{array}[]{rcl}\mathcal{S}_{1}&=&\big{\\{}k_{1}\,k_{3}\,S_{1}\,S_{2}+k_{3}^{2}\,S_{1}\,S_{2}+L_{3}\,S_{2}+(L_{2}\,S_{1}^{2}-L_{3}\,S_{1})\,k_{3}=0,\
k_{2}\,k_{3}\,S_{1}-L_{3}=0,\ k_{3}\neq 0,\\\\[5.0pt] &&\
k_{4}\,S_{1}-S_{1}^{2}+S_{2}=0,\ k_{5}+S_{1}=0,\
L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}+L_{3}\,S_{1}-S_{2}^{2}=0,\ S_{1}\neq 0,\
S_{2}\neq 0\big{\\}}\end{array}$

is the most generic one. Here, even in the most generic case, the solutions
for $\boldsymbol{k}$ are only valid for
$\overline{\boldsymbol{v}}\in\mathbb{C}^{5}$ lying on a Zariski open subset of
a hypersurface defined by

$L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}+L_{3}\,S_{1}-S_{2}^{2}=0,\ S_{1}\neq
0,\,S_{2}\neq 0.$ (35)

Solving again the remaining equations successively for $k_{5}$, $k_{4}$,
$k_{3}$, $k_{2}$, $k_{1}$, we obtain the solutions

$\displaystyle k_{5}$ $\displaystyle=-S_{1}$ $\displaystyle k_{4}$
$\displaystyle=\frac{S_{1}^{2}-S_{2}}{S_{1}}$ $\displaystyle k_{3}$
$\displaystyle\quad\mathrm{arbitrary}\,\mathrm{in}\,\mathbb{C}\setminus\\{0\\},$
$\displaystyle k_{2}$ $\displaystyle=\frac{L_{3}}{k_{3}\,S_{1}},$
$\displaystyle k_{1}$
$\displaystyle=\frac{-1}{k_{3}\,S_{1}\,S_{2}}(k_{3}^{2}\,S_{1}\,S_{2}+L_{3}\,S_{2}+(L_{2}\,S_{1}^{2}-L_{3}\,S_{1})\,k_{3}),$

where the expression for $k_{2}$ and $k_{1}$ depend on the choice made for
$k_{3}$.

The conditions (35) contain one equality, therefore the inverse problem has
solutions on an open domain of dimension 4. Furthermore, the set of solutions
is infinite. According to the Section 3.4 this means that the model $M3$ is
not solvable.

###  5.4 A Thomas decomposition for model 4

In case of model 4 the symmetrized algebraic system is

$\begin{array}[]{rcl}\mathcal{S}&=&\big{\\{}-S_{1}=k_{5},\
-S_{2}=k_{5}\,(k_{1}+k_{2}+k_{3}+L_{1}),\
L_{1}=-k_{1}-k_{2}-k_{3}-k_{4}-k_{5},\\\\[5.0pt] &&\
L_{2}=k_{1}\,k_{3}+k_{1}\,k_{4}+k_{1}\,k_{5}+k_{2}\,k_{3}+k_{2}\,k_{5}+k_{3}\,k_{4}+k_{3}\,k_{5},\
L_{3}=-k_{1}\,k_{3}\,k_{5}-k_{2}\,k_{3}\,k_{5}\big{\\}}.\end{array}$

Using the ranking

$k_{1}>k_{2}>k_{3}>k_{4}>k_{5}>\boldsymbol{v}$

the Maple implementation of algebraic Thomas decomposition returns 21 simple
systems (see subsection A.3 of the appendix for the systems
$\mathcal{S}_{2},\dots,\mathcal{S}_{21}$). The first simple system

$\begin{array}[]{rcl}\mathcal{S}_{1}&=&\big{\\{}L_{1}\,S_{1}^{2}-L_{2}\,S_{1}-S_{1}\,S_{2}+(S_{1}^{2}-S_{2})\,k_{1}+(S_{1}^{2}-S_{2})\,k_{3}+L_{3}=0,\\\\[5.0pt]
&&\
L_{2}\,S_{1}^{2}-L_{1}\,S_{1}\,S_{2}-L_{3}\,S_{1}+S_{2}^{2}+(S_{1}^{3}-S_{1}\,S_{2})\,k_{2}=0,\
k_{3}^{2}\,S_{1}+(L_{1}\,S_{1}-S_{2})\,k_{3}+L_{3}=0,\\\\[5.0pt] &&\
k_{4}\,S_{1}-S_{1}^{2}+S_{2}=0,\ k_{5}+S_{1}=0,\
L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}\neq 0,\
L_{3}\neq 0,\\\\[5.0pt] &&\ S_{1}^{3}-S_{1}\,S_{2}\neq 0,\ S_{2}\neq
0\big{\\}}.\end{array}$

is the most generic one. Indeed, all inequations appearing in
$\mathcal{S}_{1}$ involve only the variables $\boldsymbol{v}$, namely

$L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}\neq 0,\
L_{3}\neq 0,\ S_{1}^{3}-S_{1}\,S_{2}\neq 0,\ S_{2}\neq 0,$

and so the solutions for $\boldsymbol{k}$ are valid for all
$\overline{\boldsymbol{v}}\in\mathbb{C}^{5}$ of the Zariski open subset
defined by these inequations. We can now successively solve the remaining
equations in $\mathcal{S}_{1}$ for the variables $k_{5}$, $k_{4}$, $k_{3}$,
$k_{2}$ and $k_{1}$, where the equations for $k_{5}$, $k_{4}$, $k_{2}$,
$k_{1}$ are all linear in their respective leaders except for the equation
with leader $k_{3}$,

$k_{3}^{2}\,S_{1}+(L_{1}\,S_{1}-S_{2})\,k_{3}+L_{3}=0,$

which is quadratic. Since we want to consider only real roots for $k_{3}$, we
require that the discriminant is positive, that is we change the inequation
$L_{1}^{2}S_{1}^{2}-2L_{1}S_{1}S_{2}-4L_{3}S_{1}+S_{2}^{2}\neq 0$ into
$L_{1}^{2}S_{1}^{2}-2L_{1}S_{1}S_{2}-4L_{3}S_{1}+S_{2}^{2}>0$. Solving
successively the equations with leader $k_{5}$, $k_{4}$, $k_{3}$, $k_{2}$,
$k_{1}$ using for $k_{3}$ the quadratic formula, we obtain the solutions

$\displaystyle k_{5}$ $\displaystyle=-S_{1},$ $\displaystyle k_{4}$
$\displaystyle=\frac{S_{1}^{2}-S_{2}}{S_{1}},$ $\displaystyle k_{3}^{1,2}$
$\displaystyle=-\frac{L_{1}\,S_{1}-S_{2}\pm\sqrt{L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}}}{2\,S_{1}},$
$\displaystyle k_{2}$
$\displaystyle=\frac{L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}+L_{3}\,S_{1}-S_{2}^{2}}{S_{1}^{3}-S_{1}\,S_{2}},$
$\displaystyle k_{1}$
$\displaystyle=\frac{-L_{1}\,S_{1}^{2}+L_{2}\,S_{1}+S_{1}\,S_{2}-(S_{1}^{2}-S_{2})\,k_{3}^{1,2}-L_{3}}{S_{1}^{2}-S_{2}}$

subject to the conditions

$L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}>0,\
L_{3}\neq 0,\ S_{1}^{3}-S_{1}\,S_{2}\neq 0,\ S_{2}\neq 0.$ (36)

Note that the solution for $k_{1}$ depends on the choice of the root
$k_{3}^{1,2}$.

The generic simple system provides a finite set of solutions of the inverse
problem on the open domain of dimension 5 defined by (36). According to the
Section 3.4 this means that the model $M4$ is solvable.

###  5.5 A Thomas decomposition for model 8

In case of model 8 the symmetrized algebraic system is

$\begin{array}[]{rcl}\mathcal{S}&=&\big{\\{}-S_{1}=k_{5},\
-S_{2}=k_{5}\,(k_{1}+k_{2}+L_{1}),\
L_{1}=-k_{1}-k_{2}-k_{3}-k_{4}-k_{5},\\\\[5.0pt] &&\
L_{3}=-k_{1}\,k_{2}\,k_{5},\
L_{2}=k_{1}\,k_{2}+k_{1}\,k_{3}+k_{1}\,k_{4}+k_{1}\,k_{5}+k_{2}\,k_{3}+k_{2}\,k_{5}\big{\\}}.\end{array}$

We compute the algebraic Thomas decomposition of $\mathcal{S}$ with respect to
the ranking

$k_{1}>k_{3}>k_{4}>k_{2}>k_{5}>\boldsymbol{v}$

and obtain from Maple twelve simple systems
$\mathcal{S}_{1},\dots,\mathcal{S}_{12}$. One easily checks by comparing the
number of equations which have leader one of the variables of
$\boldsymbol{v}$, that the first simple system

$\begin{array}[]{rcl}\mathcal{S}_{1}&=&\big{\\{}S_{1}\,k_{1}\,k_{2}-L_{3}=0,\
k_{3}\,k_{2}\,S_{1}^{2}+L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}+L_{3}\,S_{1}-S_{2}^{2}+(-S_{1}^{3}+S_{1}\,S_{2})\,k_{2}=0,\\\\[5.0pt]
&&\
S_{1}^{2}\,k_{2}\,k_{4}-L_{1}\,S_{1}\,S_{2}+L_{2}\,S_{1}^{2}-L_{3}\,S_{1}+S_{2}^{2}=0,\
k_{2}^{2}\,S_{1}+(L_{1}\,S_{1}-S_{2})\,k_{2}+L_{3}=0,\\\\[5.0pt] &&\
k_{5}+S_{1}=0,\
L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}\neq 0,\
L_{3}\neq 0,\ S_{1}\neq 0\big{\\}}\end{array}$

is the most generic one. Analogously as in case of model 4 one determines the
solutions

$\displaystyle k_{5}$ $\displaystyle=-S_{1},$ $\displaystyle k_{2}^{1,2}$
$\displaystyle=-\frac{L_{1}\,S_{1}-S_{2}\pm\sqrt{L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}}}{2\,S_{1}},$
$\displaystyle k_{4}$
$\displaystyle=\frac{L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}+L_{3}\,S_{1}-S_{2}^{2}}{k_{2}^{1,2}\,S_{1}^{2}},$
$\displaystyle k_{3}$
$\displaystyle=-\frac{-S_{1}^{3}\,k_{2}^{1,2}+L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}+S_{1}\,S_{2}\,k_{2}^{1,2}+L_{3}\,S_{1}-S_{2}^{2}}{k_{2}^{1,2}\,S_{1}^{2}},$
$\displaystyle k_{1}$ $\displaystyle=\frac{L_{3}}{S_{1}\,k_{2}^{1,2}}$

of $\mathcal{S}_{1}$ subject to the conditions

$L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}>0,\
L_{3}\neq 0,\ S_{1}\neq 0.$ (37)

Note that here the solution for $k_{4}$, $k_{3}$ and $k_{1}$ depend on the
choice of the root $k_{2}^{1,2}$ and that we changed the inequation
$L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}\neq 0$
for the discriminant into
$L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}>0$ to
guarantee that the roots are real. The remaining simple systems are presented
in subsection A.4 of the appendix.

The generic simple system provides a finite set of solutions of the inverse
problem on the open domain of dimension 5 defined by (37). According to the
Section 3.4 this means that the model $M8$ is solvable.

###  5.6 A Thomas decomposition for model 9

We compute an algebraic Thomas decomposition of the algebraic system

$\displaystyle\mathcal{S}$ $\displaystyle=$
$\displaystyle\big{\\{}-S_{1}=k_{5},\ -S_{2}=k_{5}\,(k_{1}+k_{2}+L_{1}),\
L_{1}=-k_{1}-k_{2}-k_{3}-k_{4}-k_{5},$ (38) $\displaystyle\
L_{3}=-k_{1}k_{2}k_{5},\
L_{2}=k_{1}\,k_{2}+k_{1}\,k_{4}+k_{1}\,k_{5}+k_{2}\,k_{3}+k_{2}\,k_{5}\big{\\}}$

for model 9 with respect to the ranking

$k_{3}>k_{1}>k_{4}>k_{2}>k_{5}>\boldsymbol{v}.$

The Maple implementation returns $15$ simple systems
$\mathcal{S}_{1},\dots,\mathcal{S}_{15}$, where the most general simple system
is

$\begin{array}[]{rcl}\mathcal{S}_{1}&=&\big{\\{}L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}+L_{3}\,S_{1}-S_{2}^{2}+(2\,k_{2}\,S_{1}^{2}+L_{1}\,S_{1}^{2}-S_{1}\,S_{2})\,k_{3}+(-S_{1}^{3}+S_{1}\,S_{2})\,k_{2}=0,\\\\[5.0pt]
&&\
-L_{1}\,S_{1}^{2}+L_{2}\,S_{1}+S_{1}\,S_{2}+(2\,k_{2}\,S_{1}+L_{1}\,S_{1}-S_{2})\,k_{4}+(-S_{1}^{2}+S_{2})\,k_{2}-L_{3}=0,\\\\[5.0pt]
&&\ k_{1}\,S_{1}+k_{2}\,S_{1}+L_{1}\,S_{1}-S_{2}=0,\
k_{2}^{2}\,S_{1}+(L_{1}\,S_{1}-S_{2})\,k_{2}+L_{3}=0,\\\\[5.0pt] &&\
k_{5}+S_{1}=0,\
L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}\neq 0,\
L_{3}\neq 0,\ S_{1}\neq 0\big{\\}}.\end{array}$

Similar as for model 4 and 8 one determines the solutions

$\displaystyle k_{5}$ $\displaystyle=$ $\displaystyle-S_{1},$ $\displaystyle
k_{2}^{1,2}$ $\displaystyle=$
$\displaystyle-\frac{L_{1}\,S_{1}-S_{2}\pm\sqrt{L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}}}{2\,S_{1}},$
$\displaystyle k_{4}$ $\displaystyle=$
$\displaystyle\frac{L_{1}\,S_{1}^{2}+S_{1}^{2}\,k_{2}^{1,2}-L_{2}\,S_{1}-S_{1}\,S_{2}-S_{2}\,k_{2}^{1,2}+L_{3}}{2\,k_{2}^{1,2}\,S_{1}+L_{1}\,S_{1}-S_{2}},$
$\displaystyle k_{1}$ $\displaystyle=$
$\displaystyle-\frac{k_{2}^{1,2}\,S_{1}+L_{1}\,S_{1}-S_{2}}{S_{1}},$
$\displaystyle k_{3}$ $\displaystyle=$
$\displaystyle-\frac{-S_{1}^{3}\,k_{2}^{1,2}+L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}+S_{1}\,S_{2}\,k_{2}^{1,2}+L_{3}\,S_{1}-S_{2}^{2}}{S_{1}\,(2\,k_{2}^{1,2}\,S_{1}+L_{1}\,S_{1}-S_{2})}$
(39)

subject to the conditions

$L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}>0,\
L_{3}\neq 0,\ S_{1}\neq 0.$ (40)

Note that the solutions for $k_{4}$, $k_{1}$, $k_{3}$ depend on the choice of
the root $k_{2}^{1,2}$ and to guarantee real roots we changed the inequation
representing the discriminant
$L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}\neq 0$
into
$L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}>0$. The
remaining 14 simple algebraic systems can be found in subsection A.5 of the
appendix.

The generic simple system provides a finite set of solutions of the inverse
problem on the open domain of dimension 5 defined by (40). According to the
Section 3.4 this means that the model $M9$ is solvable.

##  6 Rashomon effect

###  6.1 Discrimination between variant models

For a fixed number of states $N$, one has several variant models that predict
exactly the same phase-type distribution and fit to data exactly equally well.
The variant models can be generated automatically. For instance in Figure 2 we
list, up to permutation symmetries of the transition graphs, all the models
that produce three exponential ($N=3$) phase-type distributions. One can not
discriminate between any two such variant models by using only data leading to
phase-type distributions, i.e. from the sequence of successive times-to-event.

Of course, there are other properties that can be used to distinguish between
variants. One of these properties, the solvability of the inverse problem, has
been discussed above. Among the variant models in Figure 2 models $M2$, $M4$,
$M8$, and $M9$ are solvable and the inverse problem has a finite number of
solutions. Model $M3$ is not solvable and the inverse problem has an infinite
set of solutions. Although in practice we favor solvable models because in
this case the inverse problem is well posed, ultimately only experimental
results should be used to exclude non-solvable variants.

Mathematically, it is interesting to classify solvable models. One
classification can be performed, based on the number of solutions. For
instance in Figure 2 the inverse problem has unique solution for model $M2$,
whereas it has two distinct solutions for models $M4$, $M8$, $M9$. Model $M9$
is obviously symmetric with respect to permutation of the vertices 1 and 2,
which explains the two-fold degeneracy of the solutions, but $M4$ and $M8$
have no obvious permutation symmetries. It is therefore interesting to
identify possible algebraic relations relating these three models together.
These relations will be discussed in the next subsection.

For experimental identification of the variants we would like to evaluate the
capacity of new data, other than the time-to-event, to discriminate between
variant models. Among possible data there are two natural candidates. The
first one is the steady state occupation probabilities of different states of
the Markov chain. Another candidate is the lifetime of each state, defined as
the mean time spent in the state before switching to another state. We should
emphasize that the lifetimes are independent on the state probabilities
because one can have the same probabilities with short and long lifetimes.
Reciprocally, the same state lifetimes can be arranged in different sequences
to obtain different cumulative occupation times and therefore different state
probabilities.

###  6.2 Symmetries and classification of solvable models

Solvable models are related one to another by some constraints that we
formulate in this section.

Let us consider the matrix $\boldsymbol{{\tilde{Q}}}^{red}(\boldsymbol{k})$
that is obtained from $\boldsymbol{\tilde{Q}}$ by setting $k_{2N+1}$ to zero.
$\left(\boldsymbol{{\tilde{Q}}}^{red}\right)^{T}$ generates a reduced Markov
chain describing the transitions between the states $1,\ldots,N$ when the
return state $s$ is $N$ (these are not affected by the existence of the state
$N+1$ because the chain gets back to $N$ instantly after a transition to
$N+1$, so one can ignore this transition).

Then the steady state probabilities of the states $1,\ldots,N$ can be computed
as solutions of the system:

$\displaystyle\boldsymbol{{\tilde{Q}}}^{red}(\boldsymbol{k})\begin{pmatrix}p_{1}(\boldsymbol{k})\\\
p_{2}(\boldsymbol{k})\\\ \vdots\\\ p_{N}(\boldsymbol{k})\end{pmatrix}$
$\displaystyle=$ $\displaystyle 0$ $\displaystyle
p_{1}(\boldsymbol{k})+p_{2}(\boldsymbol{k})+\ldots+p_{N}(\boldsymbol{k})$
$\displaystyle=$ $\displaystyle 1$ (41)

Let $T_{i}$ be the lifetime of the state $i$, where $1\leq i\leq N$. This can
be calculated from the diagonal elements of the matrix
$\boldsymbol{{\tilde{Q}}}^{red}$:

$T_{i}=-1/{\tilde{Q}}^{red}_{ii}=1/\sum_{j\neq i}{\tilde{Q}^{red}_{ij}}.$ (42)

The result of these calculations for $N=3$ is given in the Table 1.

The steady state probabilities and lifetimes allow us to formulate some
general constraints on solvable models as follows

######  Lemma 15.

All models that satisfy condition (i) of Definition 3 and have the same phase-
type distribution have

  1. (i)

the same value of $k_{2N-1}$,

  2. (ii)

the same value of the lifetime $T_{N}$ of the state $N$,

  3. (iii)

the same value of the steady state probability $p_{N}$ of the state $N$.

###### Proof.

If the first condition of Definition 3 is satisfied, then one can write the
symmetrized system of equations (3.3).

(i) is a direct consequence of (3.3).

From (3.3) it also follows:

$S_{1}-\frac{S_{2}}{S_{1}}=-k_{2N-1}+L_{1}+c_{N-2}(\boldsymbol{k}).$

Using the Vieta formulas (11) we find

$\frac{S_{2}}{S_{1}}=(-1)^{N-1}a_{N-1}(\boldsymbol{k})+c_{N-2}(\boldsymbol{k}).$

From the definitions of $a_{N-1}(\boldsymbol{k})$ and
$c_{N-2}(\boldsymbol{k})$ and elementary properties of the determinant we find

$(-1)^{N-1}a_{N-1}(\boldsymbol{k})=-c_{N-2}(\boldsymbol{k})+{\tilde{Q}}^{red}_{NN}-k_{2N-1}.$

It follows

$S_{1}-\frac{S_{2}}{S_{1}}=-{\tilde{Q}}^{red}_{NN}=T_{N}^{-1},$

which implies (ii).

Let $\bar{w}$ be the mean value of the phase-type distribution. The mean
number of events on an interval $[0,T]$ is $T/\bar{w}$. The same number is
equal to $Tp_{N}k_{2N-1}$ because the model goes to $N+1$ only from the state
$N$ with an intensity $k_{2N-1}$. It follows that

$\bar{w}=\frac{1}{p_{N}k_{2N-1}},$

for all models.

We already know that $k_{2N-1}$ is the same for all models, therefore (iii)
follows. ∎

As can be seen in Figure 2, $M9$ is a symmetrical model. The transformation
$\sigma_{M9}:\mathbb{R}^{5}\to\mathbb{R}^{5}$ defined by
$\sigma_{M9}(k_{1},k_{2},k_{3},k_{4},k_{5})=(k_{2},k_{1},k_{4},k_{2},k_{5})$,
that is equivalent to a permutation of the vertex labels $1$ and $2$, leaves
the inverse problem (38) invariant and therefore transforms any solution of
the problem into another solution. Solutions that are not symmetric are then
two-fold degenerated.

In order to see the symmetry of the solutions we rewrite (39) as

$\displaystyle k_{5}$ $\displaystyle=-S_{1},$ $\displaystyle k_{1}$
$\displaystyle=-\frac{L_{1}\,S_{1}-S_{2}\pm\sqrt{(L_{1}\,S_{1}-S_{2})^{2}-4\,L_{3}\,S_{1}}}{2\,S_{1}},$
$\displaystyle k_{2}$
$\displaystyle=-\frac{L_{1}\,S_{1}-S_{2}\mp\sqrt{(L_{1}\,S_{1}-S_{2})^{2}-4\,L_{3}\,S_{1}}}{2\,S_{1}},$
$\displaystyle k_{3}$
$\displaystyle=\frac{1}{2S_{1}}\left[S_{1}^{2}-S_{2}\pm\frac{L_{1}S_{1}^{3}-S_{1}^{2}S_{2}-2L_{2}S_{1}^{2}+L_{1}S_{1}S_{2}+2L_{3}S_{1}-S_{2}^{2}}{\sqrt{L_{1}S_{1}-S_{2})^{2}-4L_{3}S_{1}}}\right],$
$\displaystyle k_{4}$
$\displaystyle=\frac{1}{2S_{1}}\left[S_{1}^{2}-S_{2}\mp\frac{L_{1}S_{1}^{3}-S_{1}^{2}S_{2}-2L_{2}S_{1}^{2}+L_{1}S_{1}S_{2}+2L_{3}S_{1}-S_{2}^{2}}{\sqrt{L_{1}S_{1}-S_{2})^{2}-4L_{3}S_{1}}}\right].$
(43)

(43) shows that each of the two solutions of the inverse problem is not
symmetric, whereas the set of two solutions is symmetric with respect to
$\sigma_{M9}$ that transforms one solution into the other.

The models $M4$ and $M8$ have no symmetry of the transition graphs and no
obvious algebraic symmetry of the inverse problem. However, the following
property shows that these models are tightly related to $M9$.

######  Proposition 16.

There are bijective transformations
$f_{M,M^{\prime}}:D\subset\mathbb{R}^{5}\to D^{\prime}\subset\mathbb{R}^{5}$,
relating parameters $\boldsymbol{k}$ of model $M$ to parameters
$\boldsymbol{k^{\prime}}$ of the model $M^{\prime}$ (i.e.
$\boldsymbol{k^{\prime}}=f_{M,M^{\prime}}(\boldsymbol{k})$) where
$M,M^{\prime}$ are any pair in the list $\\{M4,M8,M9\\}$, and $D,D^{\prime}$
are open domains containing all the parameters for which the inverse problem
of $M,M^{\prime}$ has a solution. $f_{M,M^{\prime}}$ satisfy the following
constraints

  1. (i)

The parameter $k_{5}$ should be the same for all the models.

  2. (ii)

The lifetimes of the states $1,2,3$ considering the reduced $3$ states Markov
chain whose generator is $\left(\boldsymbol{\tilde{Q}}^{red}\right)^{T}$, i.e.
$\boldsymbol{\tilde{Q}}^{T}$ with $k_{5}$ set to zero, is the same for all the
models.

  3. (iii)

The probability to be in the state $3$ in the steady state of the same reduced
Markov chain, is the same for all the models.

Model $M$ with parameters $\boldsymbol{k}$ has exactly the same phase-type
distribution as model $M^{\prime}$ with phase-type $\boldsymbol{k^{\prime}}$.
The property does not hold if the model $M2$ is included in the list.

###### Proof.

Consider that there are mappings relating parameters of solvable models that
have the same phase-type distribution, in particular the same parameters of
the phase type distributions.

Then, (i) and (iii) follow from the Lemma 15. The same Lemma implies that
$T_{3}$ is the same for all models.

Longer, but elementary calculations using the solutions of the inverse problem
obtained in the section allow us to check that $T_{1}$ and $T_{2}$ are the
same for all the models in the list $M4$,$M8$,$M9$. However, this is no true
if we include $M2$ in the list (this model has different values of $T_{1}$ and
$T_{2}$).

Reciprocally, if we impose (i), (ii), (iii) one finds four equations relating
$\boldsymbol{k}$ and $\boldsymbol{k^{\prime}}$. These equations, if they have
solutions, define the mapping $f_{M,M^{\prime}}$. As a matter of fact,
solutions exist and are easy to compute.

For instance if $k_{i}$ are the parameters of $M9$ corresponding to the
parameters $k^{\prime}_{i}$ of $M8$, then using the Table 1 the conditions
(ii), (iii) read

$\displaystyle k_{1}$ $\displaystyle=$ $\displaystyle k^{\prime}_{1}$
$\displaystyle k_{2}$ $\displaystyle=$ $\displaystyle k^{\prime}_{2}$
$\displaystyle k_{3}+k_{4}$ $\displaystyle=$ $\displaystyle
k^{\prime}_{3}+k^{\prime}_{4}$
$\displaystyle\frac{k_{1}k_{2}}{k_{1}k_{2}+k_{1}k_{4}+k_{2}k_{3}}$
$\displaystyle=$
$\displaystyle\frac{k^{\prime}_{1}k^{\prime}_{2}}{k^{\prime}_{1}k^{\prime}_{2}+k^{\prime}_{1}k^{\prime}_{3}+k^{\prime}_{1}k^{\prime}_{4}+k^{\prime}_{2}k^{\prime}_{3}}$
(44)

The solutions of system (44) read

$\displaystyle k^{\prime}_{1}$ $\displaystyle=$ $\displaystyle k_{1}$
$\displaystyle k^{\prime}_{2}$ $\displaystyle=$ $\displaystyle k_{2}$
$\displaystyle k^{\prime}_{3}$ $\displaystyle=$
$\displaystyle\frac{k_{3}(k_{2}-k_{1})}{k_{2}}$ $\displaystyle k^{\prime}_{4}$
$\displaystyle=$ $\displaystyle\frac{k_{1}k_{3}+k_{2}k_{4}}{k_{2}},$ (45)

and define the bijective mapping
$f_{M9,M8}:\\{k_{1}>k_{2}>0,k_{3}>0,k_{4}>0\\}\to\\{k^{\prime}_{2}>k^{\prime}_{1}>0,k^{\prime}_{3}>0,k^{\prime}_{4}>\frac{k^{\prime}_{1}k^{\prime}_{3}}{k^{\prime}_{2}-k^{\prime}_{1}}\\}$.

Similarly, solving

$\displaystyle k_{1}$ $\displaystyle=$ $\displaystyle
k^{\prime}_{1}+k^{\prime}_{2}$ $\displaystyle k_{2}$ $\displaystyle=$
$\displaystyle k^{\prime}_{3}$ $\displaystyle k_{3}+k_{4}$ $\displaystyle=$
$\displaystyle k^{\prime}_{4}$
$\displaystyle\frac{k_{1}k_{2}}{k_{1}k_{2}+k_{1}k_{4}+k_{2}k_{3}}$
$\displaystyle=$
$\displaystyle\frac{k^{\prime}_{1}k^{\prime}_{3}+k^{\prime}_{2}k^{\prime}_{3}}{k^{\prime}_{1}k^{\prime}_{3}+k^{\prime}_{1}k^{\prime}_{4}+k^{\prime}_{2}k^{\prime}_{3}+k^{\prime}_{3}k^{\prime}_{4}},$
(46)

leads to

$\displaystyle k^{\prime}_{1}=(k_{1}k_{4}-k_{2}k_{4})/(k_{3}+k_{4}),$
$\displaystyle k^{\prime}_{2}=(k_{1}k_{3}+k_{2}k_{4})/(k_{3}+k_{4}),$
$\displaystyle k^{\prime}_{3}=k_{2},$ $\displaystyle
k^{\prime}_{4}=k_{3}+k_{4},$ (47)

that defines the map
$f_{M9,M4}:\\{k_{1}>k_{2}>0,k_{3}>0,k_{4}>0\\}\to\\{k^{\prime}_{1}>k^{\prime}_{2}>0,k^{\prime}_{2}>k^{\prime}_{3}>0,k^{\prime}_{4}>0\\}$.
∎

######  Remark 17.

The Proposition 16 explains the two-fold degeneracy of the solutions of the
inverse problem for models $M4$,$M8$. Indeed the solutions of $M4$ read as
$f_{M9,M8}(\boldsymbol{k})$ and $f_{M9,M8}\circ\sigma_{M9}(\boldsymbol{k})$
where $\boldsymbol{k}$ is a solution of $M9$.

######  Remark 18.

The Proposition 16 allows to compute the symmetries of the models $M8$ and
$M4$ as $f_{M9,M8}\circ\sigma_{M9}\circ f^{-1}_{M9,M8}$, and
$f_{M9,M4}\circ\sigma_{M9}\circ f^{-1}_{M9,M4}$, respectively.

Model | $T_{1}^{-1}$ | $T_{2}^{-1}$ | $T_{3}^{-1}$ | $P_{3}$
---|---|---|---|---
$M2$ | $k_{1}+k_{2}$ | $k_{3}$ | $k_{4}$ | $\frac{k_{2}k_{3}}{k_{1}k_{4}+k_{2}k_{3}+k_{3}k_{4}}$
$M4$ | $k_{1}+k_{2}$ | $k_{3}$ | $k_{4}$ | $\frac{k_{1}k_{3}+k_{2}k_{3}}{k_{1}k_{3}+k_{1}k_{4}+k_{2}k_{3}+k_{3}k_{4}}$
$M8$ | $k_{1}$ | $k_{2}$ | $(k_{3}+k_{4})$ | $\frac{k_{1}k_{2}}{k_{1}k_{2}+k_{1}k_{3}+k_{1}k_{4}+k_{2}k_{3}}$
$M9$ | $k_{1}$ | $k_{2}$ | $(k_{3}+k_{4})$ | $\frac{k_{1}k_{2}}{k_{1}k_{2}+k_{1}k_{4}+k_{2}k_{3}}$
Table 1: Reciprocal lifetimes $T_{i}^{-1}$ and steady state occupation
probability $P_{3}$ for the solvable models $M2$, $M4$, $M8$, $M9$.

##  7 Application to transcriptional bursting models

###  7.1 Biological considerations

Transcription is the fundamental molecular process that copies the DNA into
RNA. The synthesis of RNA requires a dedicated molecular machinery, including
a RNA polymerase (RNAP) enzyme that binds the DNA in specific regulatory
regions able to activate this process, called promoter regions. Transcription
requires multiple preceding steps such as chromatin opening and formation of
transcription complexes needed for activation of the polymerase. Here we
consider that transcription has a small number of limiting steps that can be
modeled as transition between discrete states (Figure 3 A). We classify these
states in three categories: productive ON states that can initiate
transcription, non-productive OFF states that can not initiate nor resume
transcription, and paused states PAUSE in which initiated transcription stops
and can resume later or abort.

The promoter starts transcription in the productive ON state, when it can
trigger several departures of RNAP molecules along DNA. After ON, the RNAP can
eventually stop in PAUSE or commit to irreversible productive elongation whose
start is modeled by the state EL. Elongating polymerases move away from the
promoter (to proceed into the gene body), thereby freeing the promoter that
instantly enters either a new productive ON state or a non-productive state
OFF.

Live imaging of transcription allows the detection of the synthesis of new
mRNA molecules in real time and for each transcription site. Using the tool
BurstDeconv [13] we can deconvolve the live imaging signal and identify the
time of each RNA synthesis event. This means that for each transcription site
we observe the sequence of EL states (Figure 3B). The other states of the
promoter are not observed.

Strictly speaking, transcription initiation is not synonymous with the start
of processive elongation, as transcription intiation can be abortive. Abortive
transcription initiation events cannot be detected in our experimental
settings. Therefore, in this paper, whenever we refer to transcription
initiation, we consider the processive case. We also consider that the state
EL (denoted $N+1$ in our abstract Markov chain model) can be reached only from
a state that can be either ON or PAUSE, and after reaching EL the promoter
instantly switches to an unique return state $s$, that can be either
productive or non-productive. This important assumption simplifies the problem
for two, closely related reasons. First, because elongation is instantly
followed by the transition to the $s$ state, therefore the observation of the
elongation is equivalent to the observation of the $s$ state. Second, the
times of successive observations of elongations form a renewal process. For
these reasons, the distribution of the waiting time between successive
observed events is a phase-type distribution.

The uniqueness of the return state rules out the presence of multiple ON
states. Due to this, our paper does not address scenarios in which
transcription occurs simultaneously on duplicate DNA copies formed during
recombination, a phenomenon referred to as sister chromatids. This specific
scenario, which could be effectively modeled using mixtures of phase-type
distributions, will be addressed elsewhere.

For a number of models, although the direct phase-type distribution problem
has a meaning, the inverse problem is not well posed. For instance,
Proposition 4 excludes models with obligatory pausing, because in this case
the state $N$ would be PAUSE, different from the return state $s$ which is ON.

Figure 3: Experimental setting for identifying phase-type distributions and
Markov chain models of transcriptional bursting. A) Transcription requires
multiple preceding steps such as chromatin opening and formation of
transcription complexes needed for activation of the polymerase. These steps
can be modeled as transitions between discrete steps. B) Fluorescent tagging
of mRNA structures allow observation of mRNA synthesis. The observed signal is
deconvolved to obtain the transcription initiation events. C) The phase-type
distribution is used to infer the Markov chain model.

###  7.2 Using alternative datasets to discriminate between variant models

As discussed in Section 6, several different models can lead to the same
phase-type distribution. For instance, three exponential phase-type
distribution data can have the same likelihood for models $M2$, $M4$, $M8$,
$M9$. Furthermore, the states of these models can be interpreted in various
ways. A possible interpretation is provided in the Table 2.

Possible state types
---
Model |  state $1$ |  state $2$ |  state $3$ |  state $4$
$M2$ |  OFF |  OFF |  ON |  EL
$M3$ |  OFF |  OFF |  ON |  EL
$M4$ |  PAUSE /OFF |  OFF |  ON |  EL
$M8$ |  PAUSE /OFF |  OFF |  ON |  EL
$M9$ |  PAUSE /OFF |  OFF/PAUSE |  ON |  EL

Possible biological interpretation of state types
---
long OFF |  Enhancer-Promoter disconnection / Chromatin closing
short OFF |  Transcription Factor unbinding / Pre-initiation Complex disassembly
PAUSE |  Promoter proximal pausing / Intrinsic pausing / Transcriptional roadblock

Table 2: Possible biological interpretation of the states in the models $M2$,
$M3$, $M4$, $M8$, $M9$. For $M4,M8,M9$ pausing is always abortive.

In this subsection we test the possibility to discriminate between solvable
variant models using measurements of state lifetimes or probabilities. State
occupation probabilities can be measured using various techniques, for
instance single molecule methylation footprinting [23, 22].

Proposition 16 shows that state lifetimes can not discriminate between models
$M4$,$M8$,$M9$.

The Lemma 15 shows that for any $N$, the occupancy probability $p_{N}$ and the
lifetime $T_{N}$ are the same for all models with $N$ states and the same
phase-type distribution.

After these considerations, for $N=3$ we remain with three independent
parameters $p_{1},T_{1},T_{2}$ as potential markers, discriminating between
variant solvable models, except for models $M4$,$M8$,$M9$ where the only
discriminating parameter could be $p_{1}$.

To evaluate the capacity of these markers to distinguish between variant
models we generate a dataset consisting of random parameters of the phase-type
distribution $S_{1},S_{2},L_{1},L_{2},L_{3}$. For each set of phase-type
distribution parameters we use the solutions of the inverse problem to compute
the transition rate parameters of the models $M2$,$M4$,$M8$,$M9$, excluding
cases leading to negative rate parameters. From the rate parameters we compute
the values of $p_{1}$, $T_{1}$, $T_{2}$ and the corresponding differences
$\Delta p_{1}$, $\Delta log(T_{1})$, $\Delta\log(T_{2})$ between the maximum
and minimum values among the four variant models. Finally, we use these
differences to estimate the capacity of such marker values to discriminate
between variants.

We found that $p_{3}$ can discriminate between variants in $89\%$ of the
cases, whereas $T_{1},T_{2}$ can discriminate between variants in $84\%$ of
the cases, see Figure 4. We have used random phase-type distribution
parameters spanning domain values observed in studies of transcriptional
bursting of mammalian and fruit-fly genes. This suggests that $p_{3}$ is
better for discriminating among variant models than $T_{1}$ and $T_{2}$. This
is confirmed also for the particular case of a HIV gene promoter in mammalian
cells studied in [34], where $T_{1}$ and $T_{2}$ are practically the same for
all experimental conditions for the 3 variants, see Figure 5.

Figure 4: Testing discrimination capacity of $p_{3}$, $T_{2}$ and $T_{1}$. We
have generated 100000 random parameter sets
$(\lambda_{1},\lambda_{2},\lambda_{3},A_{1},A_{2},A3)$ such that $A_{1}$,
$A_{2}$ are uniformly distributed in $[0,1]$, $A_{3}=1-A_{1}-A_{2}$,
$\lambda_{i}<0$ and $|\lambda_{i}|$ are log-uniformly distributed in
$[10^{-4},1]$ (4 decades of variation) for $1\leq i\leq 3$. From these
parameter sets we have kept only 61411 for which the solution of the inverse
problem resulted into real positive rate parameters
$k_{1},k_{2},k_{3},k_{4},k_{5}$ for at least one of the models
$M2$,$M4$,$M8$,$M9$. The discrimination capacity of $p_{3}$, $T_{2}$ and
$T_{1}$ is given by the variations $\Delta p$, $\Delta log_{10}(T_{1})$, or
$\Delta\log_{10}(T_{2})$ over models $M2$,$M4$,$M8$,$M9$ having the same
phase-type distribution parameters, the higher the better. The figures are
frequency histograms of these variations. There is no discrimination if the
variations are zero. We can see that $p_{3}$, $T_{1}$, $T_{2}$ do not
discriminate at all in $11\%$, $16\%$, and $16\%$ of the cases, respectively.

Figure 5: Values of $p_{i}$ and $T_{i}$ for three phase-type distributions
corresponding to three experiments on cell lines infected with HIV. Increasing
the level of the protein tat corresponds to decreasing the durations $T_{1}$,
$T_{2}$ of the OFF states (see [34]). For these cells, under no-tat, low-tat,
and high-tat conditions, the probabilities $p_{1}$, $p_{2}$ can discriminate
between variant models, but the lifetimes $T_{1}$,$T_{2}$ cannot.

##  8 Conclusion and perspectives.

In contrast to statistical inference that treats one model at a time, using
symbolic solutions for the inverse phase-type distribution problem allows us
to infer simultaneously all solvable models.

Thomas decomposition can be used to solve the inverse problem for models with
a small number of states (we provide results for $N=3$ but we tested the
method also for $N=4$).

The algebraic structure of the phase-type distribution inverse problem lead us
a to a classification of the Markov chain models compatible with data. Some
models are solvable and others are not. Among the three state solvable models
we found relations resulting from the discrete symmetry of the inverse
problem. Its generalization asks for different mathematical tools and we leave
it for future work.

We found examples of solvable models with an arbitrary number of states and
necessary solvability conditions but we did no find sufficient solvability
conditions. The task of finding sufficient solvability conditions, as well as
classifying all solvable models, can be very challenging.

Solvable models have exactly the same phase-type distribution and can not be
discriminated by time-to-event data. In this case, other quantities such as
state occupancy probabilities, could further discriminate between these
models.

We discussed a biological application of our approach. Mapping the parameters
of the phase-type distribution to kinetic model parameters allows to identify
the mechanisms impacting transcriptional bursting, from imaging datasets. This
idea led to a tool, BurstDeconv used for analyzing transcriptional bursting
data [13]. By employing Thomas decomposition, the computations presented here
have enabled us to expand the range of models amenable to analysis using
BurstDeconv (the current implementation of this tool covers only the two state
random-telegraph and the three states $M2$ and $M9$ models).

We showed here that variant models that can not be discriminated by time-to-
event data could be discriminated by state occupation probabilities that for
the transcriptional bursting application can be given for instance by single
molecule methylation footprinting.

Phase-type distributions arise when the sequence of observed transcription
events resulting from transcription imaging data is a renewal process.
However, we can consider other situations when this property is no longer
valid. In the context of transcriptional bursting, one illustrative scenario
involves the presence of sister chromatids, i.e. homologous DNA sequences
resulting from replication and capable of undergoing transcription. In this
case, the transcription site should be considered as having multiple active ON
states. This case has to be treated using different methods and the
corresponding problem will be addressed elsewhere.

While our paper primarily centers on gene transcription applications, it’s
important to note that the phase-type distribution models discussed herein
possess a broader range of potential uses. The most immediate of them concern
the next steps of gene expression. Monitoring RNA translation initiation
events as time-to-event data is feasible by using techniques such as SunTag
[33, 6, 14], and the method presented here can also be applied to reverse
engineer translation regulation mechanisms. Applying our multi-exponential
models in epidemiology and medical research could also prove intriguing, as it
may help discern disease stages and the trajectories of progression [25].

## Acknowledgements

OR and ED are founded by ANRS (ANRS0068). MD is supported by a CNRS Prime80
grant to ML and OR. ML is sponsored by CNRS and ERC (SyncDev and
LightRNA2Prot). We thank Werner Seiler for engaging and insightful
discussions.

## References

  * [1] Odd O Aalen, _Phase type distributions in survival analysis_ , Scandinavian journal of statistics (1995), 447–463.
  * [2] Azam Asanjarani, Benoit Liquet, and Yoni Nazarathy, _Estimation of semi-markov multi-state models: a comparison of the sojourn times and transition intensities approaches_ , The International Journal of Biostatistics 18 (2021), no. 1, 243–262.
  * [3] Søren Asmussen and Mogens Bladt, _Renewal theory and queueing algorithms for matrix-exponential distributions_ , Lecture notes in pure and applied mathematics (1996), 313–342.
  * [4] T. Bächler, V.P. Gerdt, M. Lange-Hegermann, and D. Robertz, _Algorithmic Thomas decomposition of algebraic and differential systems_ , J. Symbolic Comput. 47 (2012), no. 10, 1233–1266.
  * [5] T. Bächler and M. Lange-Hegermann, _AlgebraicThomas and DifferentialThomas: Thomas decomposition for algebraic and differential systems_, 2008-2012, (https://www.art.rwth-aachen.de/cms/MATHB/Forschung/Mathematische-Software/~lqnwi/).
  * [6] Eugenia Basyuk, Florence Rage, and Edouard Bertrand, _Rna transport from transcription to localized translation: a single molecule perspective_ , RNA biology 18 (2021), no. 9, 1221–1237.
  * [7] Maëlle Bellec, Jérémy Dufourt, George Hunt, Hélène Lenden-Hasse, Antonio Trullo, Amal Zine El Aabidine, Marie Lamarque, Marissa M. Gaskill, Heloïse Faure-Gautron, Mattias Mannervik, Melissa M. Harrison, Jean-Christophe Andrau, Cyril Favard, Ovidiu Radulescu, and Mounia Lagha, _The control of transcriptional memory by stable mitotic bookmarking_ , Nature Communications 13 (2022), no. 1, 1176.
  * [8] Albert T Bharucha-Reid, _Elements of the theory of markov processes and their applications_ , Courier Corporation, 1997.
  * [9] Mogens Bladt, _A review on phase-type distributions and their use in risk theory_ , ASTIN Bulletin: The Journal of the IAA 35 (2005), no. 1, 145–161.
  * [10] Peter Buchholz, Jan Kriege, and Iryna Felko, _Input modeling with phase-type distributions and markov models: theory and applications_ , Springer, 2014.
  * [11] Christian Commault and Stéphane Mocanu, _Phase-type distributions and representations: some results and open problems for system theory_ , International Journal of Control 76 (2003), no. 6, 566–580.
  * [12] Alexia Damour, Vera Slaninova, Ovidiu Radulescu, Edouard Bertrand, and Eugenia Basyuk, _Transcriptional stochasticity as a key aspect of hiv-1 latency_ , Viruses 15 (2023), no. 9, 1969.
  * [13] Maria Douaihy, Rachel Topno, Mounia Lagha, Edouard Bertrand, and Ovidiu Radulescu, _Burstdeconv: a signal deconvolution method to uncover mechanisms of transcriptional bursting in live cells_ , Nucleic Acids Research 51 (2023), no. 16, e88–e88.
  * [14] Jeremy Dufourt, Maelle Bellec, Antonio Trullo, Matthieu Dejean, Sylvain De Rossi, Cyril Favard, and Mounia Lagha, _Imaging translation dynamics in live embryos reveals spatial heterogeneities_ , Science 372 (2021), no. 6544, 840–844.
  * [15] Jeremy Dufourt, Antonio Trullo, Jennifer Hunter, Carola Fernandez, Jorge Lazaro, Matthieu Dejean, Lucas Morales, Saida Nait-Amer, Katharine N. Schulz, Melissa M. Harrison, Cyril Favard, Ovidiu Radulescu, and Mounia Lagha, _Temporal control of gene expression by the pioneer factor zelda through transient interactions in hubs_ , Nature Communications 9 (2018), no. 1, 5194.
  * [16] Daniel Dufresne, _Fitting combinations of exponentials to probability distributions_ , Applied Stochastic Models in Business and Industry 23 (2007), no. 1, 23–48.
  * [17] Mark Fackrell, _Modelling healthcare systems with phase-type distributions_ , Health care management science 12 (2009), 11–26.
  * [18] MJ Faddy, _A structured compartmental model for drug kinetics_ , Biometrics (1993), 243–248.
  * [19] William Feller, _An introduction to probability theory and its applications, volume 2_ , John Wiley & Sons, 1966.
  * [20] William Fulton and Joe Harris, _Representation theory. a first course_ , vol. 129, Springer, 1991.
  * [21] V.P. Gerdt, M. Lange-Hegermann, and D. Robertz, _The Maple package TDDS for computing Thomas decompositions of systems of nonlinear PDEs_ , Comp. Phys. Comm. 234 (2019), 202–215.
  * [22] Rozemarijn WD Kleinendorst, Guido Barzaghi, Mike L Smith, Judith B Zaugg, and Arnaud R Krebs, _Genome-wide quantification of transcription factor binding at single-dna-molecule resolution using methyl-transferase footprinting_ , Nature protocols 16 (2021), no. 12, 5673–5706.
  * [23] Arnaud R Krebs, _Studying transcription factor function in the genome at molecular resolution_ , Trends in Genetics 37 (2021), no. 9, 798–806.
  * [24] Markus Lange-Hegermann, Daniel Robertz, Werner M. Seiler, and Matthias Seiß, _Singularities of algebraic differential equations_ , Advances in Applied Mathematics 131 (2021), 102266.
  * [25] Benoit Liquet, Jean-François Timsit, and Virginie Rondeau, _Investigating hospital heterogeneity with a multi-state frailty model: application to nosocomial pneumonia disease in intensive care units_ , BMC medical research methodology 12 (2012), no. 1, 1–14.
  * [26] Marcel F Neuts, _Probability distributions of phase type_ , Liber Amicorum Prof. Emeritus H. Florin (1975).
  * [27] Colm Art O’Cinneide, _Phase-type distributions: open problems and a few properties_ , Stochastic Models 15 (1999), no. 4, 731–757.
  * [28] Virginia L. Pimmett, Matthieu Dejean, Carola Fernandez, Antonio Trullo, Edouard Bertrand, Ovidiu Radulescu, and Mounia Lagha, _Quantitative imaging of transcription in living drosophila embryos reveals the impact of core promoter motifs on promoter state dynamics_ , Nature Communications 12 (2021), no. 1, 4504.
  * [29] V Ramaswami and David M Lucantoni, _Algorithms for the multi-server queue with phase type service_ , Stochastic Models 1 (1985), no. 3, 393–417.
  * [30] V Ramaswami and Marcel F Neuts, _Some explicit formulas and computational methods for infinite-server queues with phase-type arrivals_ , Journal of Applied Probability 17 (1980), no. 2, 498–514.
  * [31] Igor Shafarevich, _Basic algebraic geometry 1_ , Springer, 1972.
  * [32] Kieran Stone, Reyer Zwiggelaar, Phil Jones, and Neil Mac Parthaláin, _A systematic review of the prediction of hospital length of stay: Towards a unified framework_ , PLOS Digital Health 1 (2022), no. 4, e0000017.
  * [33] Marvin E Tanenbaum, Luke A Gilbert, Lei S Qi, Jonathan S Weissman, and Ronald D Vale, _A protein-tagging system for signal amplification in gene expression and fluorescence imaging_ , Cell 159 (2014), no. 3, 635–646.
  * [34] Katjana Tantale, Encar Garcia-Oliver, Marie-Cécile Robert, Adèle L’Hostis, Yueyuxiao Yang, Nikolay Tsanov, Rachel Topno, Thierry Gostan, Alja Kozulic-Pirher, Meenakshi Basu-Shrivastava, Kamalika Mukherjee, Vera Slaninova, Jean-Christophe Andrau, Florian Mueller, Eugenia Basyuk, Ovidiu Radulescu, and Edouard Bertrand, _Stochastic pausing at latent hiv-1 promoters generates transcriptional bursting_ , Nature Communications 12 (2021), no. 1, 4503.
  * [35] J.M. Thomas, _Differential systems_ , Colloquium Publications XXI, American Mathematical Society, New York, 1937.
  * [36] , _Systems and roots_ , W. Byrd Press, 1962.

##  Appendix A Appendix: List of simple Thomas decomposition systems

###  A.1 The Remaining Simple Systems for Model 2

$\displaystyle\mathcal{S}_{2}$
$\displaystyle=\big{\\{}-L_{2}^{2}\,S_{1}^{2}-L_{3}\,S_{1}^{3}+2\,L_{2}\,L_{3}\,S_{1}-L_{3}^{2}+(L_{2}\,S_{1}^{3}-L_{3}\,S_{1}^{2})\,k_{1}+(L_{2}\,S_{1}^{3}-L_{3}\,S_{1}^{2})\,L_{1}=0,$
$\displaystyle\ L_{3}\,S_{1}+(L_{2}\,S_{1}-L_{3})\,k_{3}=0,\
k_{2}\,S_{1}^{2}+L_{2}\,S_{1}-L_{3}=0,\ k_{4}-S_{1}=0,\ k_{5}+S_{1}=0,$
$\displaystyle\ L_{2}\,S_{1}-L_{3}\neq 0,\ S_{1}\neq 0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{3}$
$\displaystyle=\big{\\{}k_{1}\,k_{2}\,S_{1}\,S_{2}+k_{2}^{2}\,S_{1}\,S_{2}+L_{3}\,S_{2}+(L_{2}\,S_{2}-L_{3}\,S_{1})\,k_{2}=0,\
k_{3}\,k_{2}\,S_{1}-L_{3}=0,\ k_{2}\neq 0,$ $\displaystyle\ k_{4}=0,\
k_{5}+S_{1}=0,\ L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{2}+L_{3}\,S_{1}-S_{2}^{2}=0,\
S_{1}^{2}-S_{2}=0,\ S_{2}\neq 0\big{\\}},$ $\displaystyle\mathcal{S}_{4}$
$\displaystyle=\big{\\{}k_{1}\,S_{2}+k_{3}\,S_{2}+L_{2}\,S_{1}=0,k_{2}=0,\
k_{4}\,S_{1}-S_{1}^{2}+S_{2}=0,\ k_{5}+S_{1}=0,$ $\displaystyle\
L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}-S_{2}^{2}=0,\ L_{3}=0,\ S_{1}\neq 0,\
S_{2}\neq 0\big{\\}},$ $\displaystyle\mathcal{S}_{5}$
$\displaystyle=\big{\\{}k_{1}+k_{3}+L_{1}=0,\ k_{2}=0,\ k_{4}-S_{1}=0,\
k_{5}+S_{1}=0,\ L_{2}=0,\ L_{3}=0,$ $\displaystyle\ S_{1}\neq 0,\
S_{2}=0\big{\\}},$ $\displaystyle\mathcal{S}_{6}$
$\displaystyle=\big{\\{}k_{1}\,k_{2}+k_{2}^{2}+k_{4}^{2}+k_{4}\,L_{1}+(2\,k_{4}+L_{1})\,k_{2}+L_{2}=0,\
k_{3}\,k_{2}-k_{2}\,k_{4}-k_{4}^{2}-k_{4}\,L_{1}-L_{2}=0,$ $\displaystyle\
k_{2}\neq 0,\ k_{5}=0,\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{7}$
$\displaystyle=\big{\\{}k_{1}+k_{3}+k_{4}+L_{1}=0,\ k_{2}=0,\
k_{4}^{2}+k_{4}\,L_{1}+L_{2}=0,\ k_{5}=0,\ L_{1}^{2}-4\,L_{2}\neq 0,$
$\displaystyle\ L_{2}\neq 0,\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{8}$ $\displaystyle=\big{\\{}k_{1}+k_{3}=0,\
k_{2}=0,\ k_{4}+L_{1}=0,\ k_{5}=0,\ L_{1}\neq 0,\ L_{2}=0,\ L_{3}=0,\
S_{1}=0,\ S_{2}=0\big{\\}},$ $\displaystyle\mathcal{S}_{9}$
$\displaystyle=\big{\\{}k_{1}+k_{3}+L_{1}=0,\ k_{2}=0,\ k_{4}=0,\ k_{5}=0,\
L_{1}\neq 0,\ L_{2}=0,\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{10}$
$\displaystyle=\big{\\{}2\,k_{1}+2\,k_{3}+L_{1}=0,\ k_{2}=0,\
2\,k_{4}+L_{1}=0,\ k_{5}=0,\ L_{1}^{2}-4\,L_{2}=0,\ L_{2}\neq 0,\ L_{3}=0,$
$\displaystyle\ S_{1}=0,\ S_{2}=0\big{\\}},$ $\displaystyle\mathcal{S}_{11}$
$\displaystyle=\big{\\{}k_{1}+k_{3}=0,\ k_{2}=0,\ k_{4}=0,\ k_{5}=0,\
L_{1}=0,\ L_{2}=0,\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}}.$

###  A.2 The Remaining Simple Systems for Model 3

$\displaystyle\mathcal{S}_{2}$
$\displaystyle=\big{\\{}k_{1}\,k_{3}\,S_{1}+k_{3}^{2}\,S_{1}+k_{3}\,L_{1}\,S_{1}+L_{3}=0,\
k_{2}\,k_{3}\,S_{1}-L_{3}=0,\ k_{3}\neq 0,\ k_{4}-S_{1}=0,$ $\displaystyle\
k_{5}+S_{1}=0,\ L_{2}\,S_{1}-L_{3}=0,\ S_{1}\neq 0,\ S_{2}=0\\},$
$\displaystyle\mathcal{S}_{3}$
$\displaystyle=\big{\\{}k_{1}\,S_{2}+k_{2}\,S_{2}+L_{2}\,S_{1}=0,\ k_{3}=0,\
k_{4}\,S_{1}-S_{1}^{2}+S_{2}=0,\ k_{5}+S_{1}=0,$ $\displaystyle\
L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}-S_{2}^{2}=0,\ L_{3}=0,\ S_{1}\neq 0,\
S_{2}\neq 0\big{\\}},$ $\displaystyle\mathcal{S}_{4}$
$\displaystyle=\big{\\{}k_{1}+k_{2}+L_{1}=0,\ k_{3}=0,\ k_{4}-S_{1}=0,\
k_{5}+S_{1}=0,\ L_{2}=0,\ L_{3}=0,\ $ $\displaystyle\ S_{1}\neq 0,\
S_{2}=0\big{\\}},$ $\displaystyle\mathcal{S}_{5}$
$\displaystyle=\big{\\{}k_{1}\,k_{3}+k_{3}^{2}+k_{4}^{2}+k_{4}\,L_{1}+(k_{4}+L_{1})\,k_{3}+L_{2}=0,\
k_{2}\,k_{3}-k_{4}^{2}-k_{4}\,L_{1}-L_{2}=0,$ $\displaystyle\ k_{3}\neq 0,\
k_{5}=0,\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{6}$
$\displaystyle=\big{\\{}k_{1}+k_{2}+k_{4}+L_{1}=0,\ k_{3}=0,\
k_{4}^{2}+k_{4}\,L_{1}+L_{2}=0,\ k_{5}=0,\ L_{1}^{2}-4\,L_{2}\neq 0,$
$\displaystyle\ L_{2}\neq 0,\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{7}$ $\displaystyle=\big{\\{}k_{1}+k_{2}=0,\
k_{3}=0,\ k_{4}+L_{1}=0,\ k_{5}=0,\ L_{1}\neq 0,\ L_{2}=0,\ L_{3}=0,\
S_{1}=0,$ $\displaystyle\ S_{2}=0\big{\\}},$ $\displaystyle\mathcal{S}_{8}$
$\displaystyle=\big{\\{}k_{1}+k_{2}+L_{1}=0,\ k_{3}=0,\ k_{4}=0,\ k_{5}=0,\
L_{1}\neq 0,\ L_{2}=0,\ L_{3}=0,\ S_{1}=0,$ $\displaystyle\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{9}$
$\displaystyle=\big{\\{}2\,k_{1}+2\,k_{2}+L_{1}=0,\ k_{3}=0,\
2\,k_{4}+L_{1}=0,\ k_{5}=0,\ L_{1}^{2}-4\,L_{2}=0,\ L_{2}\neq 0,$
$\displaystyle\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{10}$ $\displaystyle=\big{\\{}k_{1}+k_{2}=0,\
k_{3}=0,\ k_{4}=0,\ k_{5}=0,\ L_{1}=0,\ L_{2}=0,\ L_{3}=0,\ S_{1}=0,\
S_{2}=0\big{\\}}.$

###  A.3 The Remaining Simple Systems for Model 4

$\displaystyle\mathcal{S}_{2}$
$\displaystyle=\big{\\{}L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}-S_{2}^{2}+(S_{1}^{3}-S_{1}\,S_{2})\,k_{1}=0,\
k_{3}\,S_{1}+L_{1}\,S_{1}-S_{2}=0,$ $\displaystyle\
-L_{1}\,S_{1}\,S_{2}+L_{2}\,S_{1}^{2}+S_{2}^{2}+(S_{1}^{3}-S_{1}\,S_{2})\,k_{2}=0,\
k_{4}\,S_{1}-S_{1}^{2}+S_{2}=0,\ k_{5}+S_{1}=0,$ $\displaystyle\
L_{1}\,S_{1}-S_{2}\neq 0,\ L_{3}=0,\ S_{1}^{3}-S_{1}\,S_{2}\neq 0,\ S_{2}\neq
0\big{\\}},$ $\displaystyle\mathcal{S}_{3}$
$\displaystyle=\big{\\{}L_{1}\,S_{1}^{2}-L_{2}\,S_{1}-S_{1}\,S_{2}+(S_{1}^{2}-S_{2})\,k_{1}=0,\
-L_{1}\,S_{1}\,S_{2}+L_{2}\,S_{1}^{2}+S_{2}^{2}+(S_{1}^{3}-S_{1}\,S_{2})\,k_{2}=0,$
$\displaystyle\ k_{3}=0,\ k_{4}\,S_{1}-S_{1}^{2}+S_{2}=0,\ k_{5}+S_{1}=0,\
L_{1}\,S_{1}-S_{2}\neq 0,\ L_{3}=0,\ S_{1}^{3}-S_{1}\,S_{2}\neq 0,\ S_{2}\neq
0\big{\\}},$ $\displaystyle\mathcal{S}_{4}$
$\displaystyle=\big{\\{}-2\,L_{2}\,S_{1}^{2}-S_{1}^{2}\,S_{2}+2\,L_{3}\,S_{1}-S_{2}^{2}+(2\,S_{1}^{3}-2\,S_{1}\,S_{2})\,k_{1}+(S_{1}^{3}+S_{1}\,S_{2})\,L_{1}=0,$
$\displaystyle\
-L_{1}\,S_{1}\,S_{2}+L_{2}\,S_{1}^{2}-L_{3}\,S_{1}+S_{2}^{2}+(S_{1}^{3}-S_{1}\,S_{2})\,k_{2}=0,\
2\,k_{3}\,S_{1}+L_{1}\,S_{1}-S_{2}=0,$ $\displaystyle\
k_{4}\,S_{1}-S_{1}^{2}+S_{2}=0,\ k_{5}+S_{1}=0,\
L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}=0,$
$\displaystyle\ L_{3}\neq 0,\ S_{1}^{3}-S_{1}\,S_{2}\neq 0,\ S_{2}\neq
0\big{\\}},$ $\displaystyle\mathcal{S}_{5}$
$\displaystyle=\big{\\{}-L_{2}\,S_{1}+(S_{1}^{2}-S_{2})\,k_{1}=0,\
L_{2}\,S_{1}+(S_{1}^{2}-S_{2})\,k_{2}=0,\ k_{3}=0,\
k_{4}\,S_{1}-S_{1}^{2}+S_{2}=0,$ $\displaystyle\ k_{5}+S_{1}=0,\
L_{1}\,S_{1}-S_{2}=0,\ L_{3}=0,\ S_{1}^{3}-S_{1}\,S_{2}\neq 0,\ S_{2}\neq
0\big{\\}},$ $\displaystyle\mathcal{S}_{6}$
$\displaystyle=\big{\\{}k_{1}\,S_{1}^{2}+k_{3}\,S_{1}^{2}+L_{1}\,S_{1}^{2}-L_{2}\,S_{1}+L_{3}=0,\
k_{2}\,S_{1}^{2}+L_{2}\,S_{1}-L_{3}=0,\ k_{4}-S_{1}=0,$ $\displaystyle\
k_{3}^{2}\,S_{1}+k_{3}\,L_{1}\,S_{1}+L_{3}=0,\ k_{5}+S_{1}=0,\
L_{1}^{2}\,S_{1}-4\,L_{3}\neq 0,\ L_{3}\neq 0,\ S_{1}\neq 0,\
S_{2}=0\big{\\}},$ $\displaystyle\mathcal{S}_{7}$
$\displaystyle=\big{\\{}k_{1}\,S_{1}-L_{2}=0,\ k_{2}\,S_{1}+L_{2}=0,\
k_{3}+L_{1}=0,\ k_{4}-S_{1}=0,\ k_{5}+S_{1}=0,\ L_{1}\neq 0,\ L_{3}=0,$
$\displaystyle\ S_{1}\neq 0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{8}$
$\displaystyle=\big{\\{}k_{1}\,S_{1}+L_{1}\,S_{1}-L_{2}=0,\
k_{2}\,S_{1}+L_{2}=0,\ k_{3}=0,\ k_{4}-S_{1}=0,\ k_{5}+S_{1}=0,\ L_{1}\neq 0,$
$\displaystyle\ L_{3}=0,\ S_{1}\neq 0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{9}$
$\displaystyle=\big{\\{}2\,k_{1}\,S_{1}^{2}+L_{1}\,S_{1}^{2}-2\,L_{2}\,S_{1}+2\,L_{3}=0,\
k_{2}\,S_{1}^{2}+L_{2}\,S_{1}-L_{3}=0,\ 2\,k_{3}+L_{1}=0,\ k_{4}-S_{1}=0,$
$\displaystyle\ k_{5}+S_{1}=0,\ L_{1}^{2}\,S_{1}-4\,L_{3}=0,\ L_{3}\neq 0,\
S_{1}\neq 0,\ S_{2}=0\big{\\}},$ $\displaystyle\mathcal{S}_{10}$
$\displaystyle=\big{\\{}k_{1}\,S_{1}-L_{2}=0,\ k_{2}\,S_{1}+L_{2}=0,\
k_{3}=0,\ k_{4}-S_{1}=0,\ k_{5}+S_{1}=0,\ L_{1}=0,\ L_{3}=0,$ $\displaystyle\
S_{1}\neq 0,\ S_{2}=0\big{\\}},$ $\displaystyle\mathcal{S}_{11}$
$\displaystyle=\big{\\{}k_{1}\,S_{1}\,S_{2}+k_{2}\,S_{1}\,S_{2}+k_{3}\,S_{1}\,S_{2}+L_{2}\,S_{2}-L_{3}\,S_{1}=0,\
k_{4}=0,\ k_{5}+S_{1}=0,$ $\displaystyle\
k_{3}^{2}\,S_{1}\,S_{2}+L_{3}\,S_{2}+(L_{2}\,S_{2}-L_{3}\,S_{1})\,k_{3}=0,\
L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{2}+L_{3}\,S_{1}-S_{2}^{2}=0,$ $\displaystyle\
L_{2}^{2}\,S_{2}-2\,L_{2}\,L_{3}\,S_{1}-4\,L_{3}\,S_{1}\,S_{2}+L_{3}^{2}\neq
0,\ L_{3}\neq 0,\ S_{1}^{2}-S_{2}=0,\ S_{2}\neq 0\big{\\}},$
$\displaystyle\mathcal{S}_{12}$ $\displaystyle=\big{\\{}k_{1}+k_{2}=0,\
k_{3}\,S_{1}+L_{2}=0,\ k_{4}=0,\ k_{5}+S_{1}=0,\ L_{1}\,S_{1}-L_{2}-S_{2}=0,\
L_{2}\neq 0,$ $\displaystyle\ L_{3}=0,\ S_{1}^{2}-S_{2}=0,\ S_{2}\neq
0\big{\\}},$ $\displaystyle\mathcal{S}_{13}$
$\displaystyle=\big{\\{}k_{1}\,S_{1}+k_{2}\,S_{1}+L_{2}=0,\ k_{3}=0,\
k_{4}=0,\ k_{5}+S_{1}=0,\ L_{1}\,S_{1}-L_{2}-S_{2}=0,\ L_{2}\neq 0,$
$\displaystyle\ L_{3}=0,\ S_{1}^{2}-S_{2}=0,\ S_{2}\neq 0\big{\\}},$
$\displaystyle\mathcal{S}_{14}$
$\displaystyle=\big{\\{}2\,k_{1}\,S_{1}\,S_{2}+2\,k_{2}\,S_{1}\,S_{2}+L_{2}\,S_{2}-L_{3}\,S_{1}=0,\
2\,k_{3}\,S_{1}\,S_{2}+L_{2}\,S_{2}-L_{3}\,S_{1}=0,\ k_{4}=0,$ $\displaystyle\
k_{5}+S_{1}=0,\ L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{2}+L_{3}\,S_{1}-S_{2}^{2}=0,\
L_{2}^{2}\,S_{2}-2\,L_{2}\,L_{3}\,S_{1}-4\,L_{3}\,S_{1}\,S_{2}+L_{3}^{2}=0,$
$\displaystyle\ L_{3}\neq 0,\ S_{1}^{2}-S_{2}=0,\ S_{2}\neq 0\\},$
$\displaystyle\mathcal{S}_{15}$ $\displaystyle=\big{\\{}k_{1}+k_{2}=0,\
k_{3}=0,\ k_{4}=0,\ k_{5}+S_{1}=0,\ L_{1}\,S_{1}-S_{2}=0,\ L_{2}=0,\ L_{3}=0,$
$\displaystyle\ S_{1}^{2}-S_{2}=0,\ S_{2}\neq 0\big{\\}},$
$\displaystyle\mathcal{S}_{16}$
$\displaystyle=\big{\\{}k_{1}\,k_{4}-k_{3}^{2}-k_{3}\,L_{1}-L_{2}=0,\
k_{2}\,k_{4}+k_{3}^{2}+k_{4}^{2}+k_{4}\,L_{1}+(k_{4}+L_{1})\,k_{3}+L_{2}=0,\
k_{4}\neq 0,$ $\displaystyle\ k_{5}=0,\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{17}$
$\displaystyle=\big{\\{}k_{1}+k_{2}+k_{3}+L_{1}=0,\
k_{3}^{2}+k_{3}\,L_{1}+L_{2}=0,\ k_{4}=0,\ k_{5}=0,\ L_{1}^{2}-4\,L_{2}\neq
0,\ L_{2}\neq 0,$ $\displaystyle\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{18}$ $\displaystyle=\big{\\{}k_{1}+k_{2}=0,\
k_{3}+L_{1}=0,\ k_{4}=0,\ k_{5}=0,\ L_{1}\neq 0,\ L_{2}=0,\ L_{3}=0,\
S_{1}=0,\ S_{2}=0\big{\\}},$ $\displaystyle\mathcal{S}_{19}$
$\displaystyle=\big{\\{}k_{1}+k_{2}+L_{1}=0,\ k_{3}=0,\ k_{4}=0,\ k_{5}=0,\
L_{1}\neq 0,\ L_{2}=0,\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{20}$
$\displaystyle=\big{\\{}2\,k_{1}+2\,k_{2}+L_{1}=0,\ 2\,k_{3}+L_{1}=0,\
k_{4}=0,\ k_{5}=0,\ L_{1}^{2}-4\,L_{2}=0,\ L_{2}\neq 0,$ $\displaystyle\
L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$ $\displaystyle\mathcal{S}_{21}$
$\displaystyle=\big{\\{}k_{1}+k_{2}=0,\ k_{3}=0,\ k_{4}=0,\ k_{5}=0,\
L_{1}=0,\ L_{2}=0,\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}}.$

###  A.4 The Remaining Simple Systems for Model 8

$\displaystyle\mathcal{S}_{2}$ $\displaystyle=\big{\\{}k_{1}=0,\
-L_{1}\,S_{1}^{2}+L_{2}\,S_{1}+S_{1}\,S_{2}+(L_{1}\,S_{1}-S_{2})\,k_{3}=0,$
$\displaystyle\
L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}-S_{2}^{2}+(L_{1}\,S_{1}^{2}-S_{1}\,S_{2})\,k_{4}=0,\
L_{1}\,S_{1}+S_{1}\,k_{2}-S_{2}=0,\ k_{5}+S_{1}=0,$ $\displaystyle\
L_{1}\,S_{1}-S_{2}\neq 0,\ L_{3}=0,\ S_{1}\neq 0\big{\\}},$
$\displaystyle\mathcal{S}_{3}$
$\displaystyle=\big{\\{}(L_{1}\,S_{1}-S_{2})\,k_{1}+2\,L_{3}=0,$
$\displaystyle\
2\,L_{2}\,S_{1}^{2}+S_{1}^{2}\,S_{2}-2\,L_{3}\,S_{1}+S_{2}^{2}+(L_{1}\,S_{1}^{2}-S_{1}\,S_{2})\,k_{3}+(-S_{1}^{3}-S_{1}\,S_{2})\,L_{1}=0,$
$\displaystyle\
2\,L_{1}\,S_{1}\,S_{2}-2\,L_{2}\,S_{1}^{2}+2\,L_{3}\,S_{1}-2\,S_{2}^{2}+(L_{1}\,S_{1}^{2}-S_{1}\,S_{2})\,k_{4}=0,\
L_{1}\,S_{1}+2\,S_{1}\,k_{2}-S_{2}=0,$ $\displaystyle\ k_{5}+S_{1}=0,\
L_{1}^{2}\,S_{1}^{2}-2\,L_{1}\,S_{1}\,S_{2}-4\,L_{3}\,S_{1}+S_{2}^{2}=0,\
L_{3}\neq 0,\ S_{1}\neq 0\big{\\}},$ $\displaystyle\mathcal{S}_{4}$
$\displaystyle=\big{\\{}L_{2}\,S_{1}+S_{2}\,k_{1}=0,\
-S_{1}^{2}+S_{1}\,k_{3}+S_{1}\,k_{4}+S_{2}=0,\ k_{2}=0,\ k_{5}+S_{1}=0,$
$\displaystyle\ L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}-S_{2}^{2}=0,\ L_{3}=0,\
S_{1}\neq 0,\ S_{2}\neq 0\big{\\}},$ $\displaystyle\mathcal{S}_{5}$
$\displaystyle=\big{\\{}k_{1}+L_{1}=0,\ k_{3}+k_{4}-S_{1}=0,\ k_{2}=0,\
k_{5}+S_{1}=0,\ L_{2}=0,\ L_{3}=0,\ S_{1}\neq 0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{6}$
$\displaystyle=\big{\\{}k_{1}+k_{3}+k_{4}+k_{2}+L_{1}=0,$ $\displaystyle\
k_{3}^{2}+k_{4}^{2}+k_{2}^{2}+k_{2}\,L_{1}+(2\,k_{4}+k_{2}+L_{1})\,k_{3}+(2\,k_{2}+L_{1})\,k_{4}+L_{2}=0,$
$\displaystyle\
-L_{1}^{2}+2\,L_{1}\,k_{2}+3\,k_{2}^{2}+4\,k_{2}\,k_{4}+4\,L_{2}\neq 0,\
k_{2}\neq 0,\ k_{5}=0,\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{7}$
$\displaystyle=\big{\\{}k_{1}+k_{3}+k_{4}+L_{1}=0,\
k_{3}^{2}+k_{4}^{2}+k_{4}\,L_{1}+(2\,k_{4}+L_{1})\,k_{3}+L_{2}=0,\ k_{2}=0,\
k_{5}=0,$ $\displaystyle\ L_{1}^{2}-4\,L_{2}\neq 0,\ L_{2}\neq 0,\ L_{3}=0,\
S_{1}=0,\ S_{2}=0\big{\\}},$ $\displaystyle\mathcal{S}_{8}$
$\displaystyle=\big{\\{}k_{1}=0,\ k_{3}+k_{4}+L_{1}=0,\ k_{2}=0,\ k_{5}=0,\
L_{1}\neq 0,\ L_{2}=0,\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{9}$ $\displaystyle=\big{\\{}k_{1}+L_{1}=0,\
k_{3}+k_{4}=0,\ k_{2}=0,\ k_{5}=0,\ L_{1}\neq 0,\ L_{2}=0,\ L_{3}=0,\
S_{1}=0,\ S_{2}=0\big{\\}},$ $\displaystyle\mathcal{S}_{10}$
$\displaystyle=\big{\\{}2\,k_{1}+k_{2}+L_{1}=0,\
L_{1}^{2}-k_{2}^{2}+4\,k_{2}\,k_{3}-4\,L_{2}=0,$ $\displaystyle\
-L_{1}^{2}+2\,L_{1}\,k_{2}+3\,k_{2}^{2}+4\,k_{2}\,k_{4}+4\,L_{2}=0,\ k_{2}\neq
0,\ k_{5}=0,\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{11}$ $\displaystyle=\big{\\{}2\,k_{1}+L_{1}=0,\
2\,k_{3}+2\,k_{4}+L_{1}=0,\ k_{2}=0,\ k_{5}=0,\ L_{1}^{2}-4\,L_{2}=0,\
L_{2}\neq 0,$ $\displaystyle\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{12}$ $\displaystyle=\big{\\{}k_{1}=0,\
k_{3}+k_{4}=0,\ k_{2}=0,\ k_{5}=0,\ L_{1}=0,\ L_{2}=0,\ L_{3}=0,\ S_{1}=0,\
S_{2}=0\big{\\}}.$

###  A.5 The Remaining Simple Systems for Model 9

$\displaystyle\mathcal{S}_{2}$
$\displaystyle=\big{\\{}-L_{1}\,S_{1}^{2}+L_{2}\,S_{1}+S_{1}\,S_{2}+(L_{1}\,S_{1}-S_{2})\,k_{3}=0,\
k_{1}=0,$ $\displaystyle\
L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}-S_{2}^{2}+(L_{1}\,S_{1}^{2}-S_{1}\,S_{2})\,k_{4}=0,\
k_{2}\,S_{1}+L_{1}\,S_{1}-S_{2}=0,\ k_{5}+S_{1}=0,$ $\displaystyle\
L_{1}\,S_{1}-S_{2}\neq 0,\ L_{3}=0,\ S_{1}\neq 0\big{\\}},$
$\displaystyle\mathcal{S}_{3}$
$\displaystyle=\big{\\{}L_{1}\,S_{1}\,S_{2}-L_{2}\,S_{1}^{2}-S_{2}^{2}+(L_{1}\,S_{1}^{2}-S_{1}\,S_{2})\,k_{3}=0,\
k_{1}\,S_{1}+L_{1}\,S_{1}-S_{2}=0,$ $\displaystyle\
-L_{1}\,S_{1}^{2}+L_{2}\,S_{1}+S_{1}\,S_{2}+(L_{1}\,S_{1}-S_{2})\,k_{4}=0,\
k_{2}=0,\ k_{5}+S_{1}=0,\ L_{1}\,S_{1}-S_{2}\neq 0,$ $\displaystyle\ L_{3}=0,\
S_{1}\neq 0\big{\\}},$ $\displaystyle\mathcal{S}_{4}$
$\displaystyle=\big{\\{}k_{3}\,S_{1}+k_{4}\,S_{1}-S_{1}^{2}+S_{2}=0,\
L_{2}\,S_{1}+(S_{1}^{2}+S_{2})\,k_{1}-L_{3}=0,\ k_{5}+S_{1}=0,$
$\displaystyle\ L_{2}\,S_{1}+(S_{1}^{2}+S_{2})\,k_{2}-L_{3}=0,\
-2\,L_{2}\,S_{1}^{2}-S_{1}^{2}\,S_{2}+2\,L_{3}\,S_{1}-S_{2}^{2}+(S_{1}^{3}+S_{1}\,S_{2})\,L_{1}=0,$
$\displaystyle\
L_{2}^{2}\,S_{1}^{3}-2\,L_{2}\,L_{3}\,S_{1}^{2}+L_{3}^{2}\,S_{1}+(-S_{1}^{4}-2\,S_{1}^{2}\,S_{2}-S_{2}^{2})\,L_{3}=0,\
L_{3}\neq 0,\ S_{1}^{3}+S_{1}\,S_{2}\neq 0,\ S_{2}\neq 0\big{\\}},$
$\displaystyle\mathcal{S}_{5}$
$\displaystyle=\big{\\{}k_{3}\,S_{1}+k_{4}\,S_{1}-S_{1}^{2}+S_{2}=0,\
k_{1}=0,\ k_{2}=0,\ k_{5}+S_{1}=0,\ L_{1}\,S_{1}-S_{2}=0,\ L_{2}=0,$
$\displaystyle\ L_{3}=0,\ S_{1}^{3}+S_{1}\,S_{2}\neq 0,\ S_{2}\neq
0\big{\\}},$ $\displaystyle\mathcal{S}_{6}$
$\displaystyle=\big{\\{}k_{3}+k_{4}-S_{1}=0,\
k_{1}\,S_{1}^{2}+L_{2}\,S_{1}-L_{3}=0,\
k_{2}\,S_{1}^{2}+L_{2}\,S_{1}-L_{3}=0,\ k_{5}+S_{1}=0,$ $\displaystyle\
L_{1}\,S_{1}^{2}-2\,L_{2}\,S_{1}+2\,L_{3}=0,\
L_{2}^{2}\,S_{1}^{2}-L_{3}\,S_{1}^{3}-2\,L_{2}\,L_{3}\,S_{1}+L_{3}^{2}=0,\
L_{3}\neq 0,\ S_{1}\neq 0,\ S_{2}=0\big{\\}},$ $\displaystyle\mathcal{S}_{7}$
$\displaystyle=\big{\\{}k_{3}+k_{4}-S_{1}=0,\ k_{1}=0,\ k_{2}=0,\
k_{5}+S_{1}=0,\ L_{1}=0,\ L_{2}=0,\ L_{3}=0,\ S_{1}\neq 0,$ $\displaystyle\
S_{2}=0\big{\\}},$ $\displaystyle\mathcal{S}_{8}$
$\displaystyle=\big{\\{}k_{3}\,S_{1}+k_{4}\,S_{1}+2\,S_{2}=0,\
2\,k_{1}\,S_{1}+L_{1}\,S_{1}-S_{2}=0,\ 2\,k_{2}\,S_{1}+L_{1}\,S_{1}-S_{2}=0,\
k_{5}+S_{1}=0,$ $\displaystyle\
L_{1}^{2}\,S_{2}+2\,L_{1}\,S_{1}\,S_{2}+4\,L_{3}\,S_{1}-S_{2}^{2}=0,\
L_{2}\,S_{2}+L_{3}\,S_{1}=0,\ L_{3}\neq 0,\ S_{1}^{2}+S_{2}=0,\ S_{2}\neq
0\big{\\}},$ $\displaystyle\mathcal{S}_{9}$
$\displaystyle=\big{\\{}k_{3}\,S_{1}+k_{4}\,S_{1}+2\,S_{2}=0,\ k_{1}=0,\
k_{2}=0,\ k_{5}+S_{1}=0,\ L_{1}+S_{1}=0,\ L_{2}=0,\ L_{3}=0,$ $\displaystyle\
S_{1}^{2}+S_{2}=0,\ S_{2}\neq 0\big{\\}},$ $\displaystyle\mathcal{S}_{10}$
$\displaystyle=\big{\\{}k_{3}\,k_{4}+k_{4}^{2}+k_{2}^{2}+k_{2}\,L_{1}+(2\,k_{2}+L_{1})\,k_{4}+L_{2}=0,\
k_{1}\,k_{4}-k_{4}\,k_{2}-k_{2}^{2}-k_{2}\,L_{1}-L_{2}=0,$ $\displaystyle\
k_{4}\neq 0,\ k_{5}=0,\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{11}$
$\displaystyle=\big{\\{}k_{3}+k_{1}+k_{2}+L_{1}=0,\ k_{4}=0,\
k_{2}^{2}+k_{2}\,L_{1}+L_{2}=0,\ k_{5}=0,\ L_{1}^{2}-4\,L_{2}\neq 0,\
L_{2}\neq 0,$ $\displaystyle\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{12}$ $\displaystyle=\big{\\{}k_{3}+k_{1}=0,\
k_{4}=0,\ k_{2}+L_{1}=0,\ k_{5}=0,\ L_{1}\neq 0,\ L_{2}=0,\ L_{3}=0,\
S_{1}=0,\ S_{2}=0\big{\\}},$ $\displaystyle\mathcal{S}_{13}$
$\displaystyle=\big{\\{}k_{3}+k_{1}+L_{1}=0,\ k_{4}=0,\ k_{2}=0,\ k_{5}=0,\
L_{1}\neq 0,\ L_{2}=0,\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{14}$
$\displaystyle=\big{\\{}2\,k_{3}+2\,k_{1}+L_{1}=0,\ k_{4}=0,\
2\,k_{2}+L_{1}=0,\ k_{5}=0,\ L_{1}^{2}-4\,L_{2}=0,\ L_{2}\neq 0,$
$\displaystyle\ L_{3}=0,\ S_{1}=0,\ S_{2}=0\big{\\}},$
$\displaystyle\mathcal{S}_{15}$ $\displaystyle=\big{\\{}k_{3}+k_{1}=0,\
k_{4}=0,\ k_{2}=0,\ k_{5}=0,\ L_{1}=0,\ L_{2}=0,\ L_{3}=0,\ S_{1}=0,\
S_{2}=0\big{\\}}.$