# Identifying Markov Chain Models from Time-to-Event Data: An Algebraic Approach

Ovidiu Radulescu, Dima Grigoriev, Matthias Seiss, Maria Douaihy, Mounia Lagha, Edouard Bertrand

► To cite this version:

Ovidiu Radulescu, Dima Grigoriev, Matthias Seiss, Maria Douaihy, Mounia Lagha, et al. Identifying Markov Chain Models from Time-to-Event Data: An Algebraic Approach. Bulletin of Mathematical Biology, 2025, 87, pp.11. 10.1007/s11538-024-01385-y. hal-04958556

HAL Id: hal-04958556

https://hal.science/hal-04958556v1

Submitted on 27 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Identifying Markov chain models from time-to-event data: an algebraic approach

Ovidiu Radulescu¹,*, Dima Grigoriev², Matthias Seiss³, Maria Douaihy¹,⁴, Mounia Lagha⁴, Edouard Bertrand⁵,

¹ LPHI, University of Montpellier and CNRS, Place Eugène Bataillon, Montpellier 34095, France

² CNRS, Mathématiques, Université de Lille, 59655, Villeneuve d'Ascq, France

³ University of Kassel, Kassel, Germany

⁴ IGMM, University of Montpellier and CNRS, 1919 Rte de Mende, Montpellier 34090, France

⁵ IGH, University of Montpellier and CNRS, 141 Rue de la Cardonille, Montpellier 34094, France

November 19, 2024

## Abstract

Many biological and medical questions can be modeled using time-to-event data in finite-state Markov chains, with the phase-type distribution describing intervals between events. We solve the inverse problem: given a phase-type distribution, can we identify the transition rate parameters of the underlying Markov chain? For a specific class of *solvable* Markov models, we show this problem has a unique solution up to finite symmetry transformations, and we outline a recursive method for computing symbolic solutions for these models across any number of states. Using the Thomas decomposition technique from computer algebra, we further provide symbolic solutions for any model. Interestingly, different models with the same state count but distinct transition graphs can yield identical phase-type distributions. To distinguish among these, we propose additional properties beyond just the time to the next event. We demonstrate the method's applicability by inferring transcriptional regulation models from single-cell transcription imaging data.

**Keywords:** Markov chains, phase-type distribution, inverse problem, symmetric polynomials, Thomas decomposition, transcriptional bursting.
# 1 Introduction

Continuous-time finite Markov chains have been largely used to model biological and medical phenomena involving multiple discrete states [6]. For instance, the progression of a disease may have multiple stages, including varying severity grades, recovery after treatment, and potential relapse. Similarly, gene expression functions intermittently, and often requires several prior non-productive and productive stages to prepare and initiate transcription or translation, or to pause these processes. In both examples, the number of states and the timescales of transitions between them are usually unknown.

A common dataset for such models consists of a sequence of specific events, such as the occurrence of an observed state. In renewal processes, relevant to many applications, the times separating successive events are independent, identically distributed variables, and the sequence is characterized by the distribution of these intervals [18]. In simple cases, such as Poisson processes, this distribution is exponential. However, when multiple states are involved, this distribution deviates from the exponential. The renewal process scenario applies only when there is one observable state and the process resets to the same unobserved state after each observation. These scenarios are modeled by phase-type distributions, introduced by [34], (see Figure 1).

Phase-type distributions were used in biology and medicine for modeling transcription and translation bursting [26, 39, 25, 9, 42, 37, 13, 45], stochastic enzymatic reactions [33], ion channel dynamics [27], drug kinetics [17], population genetics [21, 22], and public health issues [16, 29, 1, 40].

In this paper, our aim is to determine the underlying Markov chain, the number of states, and the transition rates between states from the phase-type distribution. This inverse problem has previously been tackled using maximum likelihood approaches, where the data likelihood was computed for an arbitrarily chosen model and its maximization yielded optimal model parameters [1, 40].

In contrast to this direct inference approach, our method decomposes the inference problem into two parts. The first part involves regressing a parametric multi-exponential representation of the phase-type distribution [15, 13]. The second part of our approach is entirely algebraic. We formulate the inverse problem as a system of polynomial equations in the transition rate parameters and solve this system for different numbers of states and transition topologies. In doing so, we obtain explicit formulas relating transition rate parameters to phase-type distribution parameters for all Markov chains that generate the same multi-exponential distribution.

Other approaches to finding Markov chain representations of a phase-type distribution, using generating functions or the Laplace transform [31, 10], propose some, but not all, Markov chain representations and do not always offer explicit formulas for the transition rate parameters.

We provide a code implementation of our algorithms that formulate the inverse problem of any phase-type distribution as a system of polynomial equations and solve it using Thomas decomposition, a computer algebra method. Our approach also determines whether the model is solvable,
meaning whether the solution to the inverse problem is unique or unique
up to transformations by finite symmetries.

To demonstrate the relevance of our approach, we apply it to the
context of transcriptional bursting, a well-documented phenomenon in
gene expression research [14, 42, 37, 5, 13, 12].

The structure of the paper is as follows: In Section 2, we discuss the phase-type distribution problem and find its solution in the non-degenerate case, specifically when the generator of the Markov chain has distinct eigenvalues. Section 3 introduces the symbolic inverse problem, which entails finding the transition rate parameters from the phase-type distribution parameters. We demonstrate that solving this inverse problem involves a system of polynomial equations symmetric in the phase-type distribution parameters. Additionally, we discuss solvable models in Section 3, where the inverse problem has a unique solution. In Section 4, we present an example of a solvable model with an arbitrary number of states. Section 5 utilizes Thomas decomposition to solve the inverse problem. Section 6 addresses model degeneracy issues. Finally, in Section 7, we apply our findings to analyze transcriptional bursting data.

<figure><img src="image_2.png" /><figcaption>Figure 1: Phase-type distribution for a three-state Markov chain model with state 4 as the observable state. A typical dataset records the moments when 4 is reached. The chain returns instantly from 4 to 3 and restarts. The interval between successive observations, also defined as the time $T_3$ to reach 4 from 3, follows a phase-type distribution.</figcaption></figure>
# 2 Phase-type distribution, direct problem

## 2.1 Statistical definition

Consider the successive observations of a state, as illustrated in Figure 1. Following each observation, the system returns to a state $s$, so the interval between successive observations corresponds to the time required to reach the observed state from $s$, denoted by $T_s$. The distribution of $T_s$ is of the phase-type if it is derived from convolutions and mixtures (convex combinations) of exponential distributions [11]. Alternatively, it can be defined as having a Laplace transform that is a rational function (see Theorem 1 in [11]). Phase-type distributions generalize exponential distributions, including mixtures of exponential or Erlang distributions.

This paper considers multi-exponential distributions, characterized by a complementary cumulative distribution function (survival function) of the following type:

$$S(t) = \mathbb{P}[T_s > t] = \sum_{i=1}^{N} A_i \exp(\lambda_i t), \quad (1)$$

where $\lambda_i < 0$ are all negative, and $\sum_{i=1}^{N} A_i = 1$. $A_i$ are not necessarily all positive, but are constrained by the conditions $S(t) > 0, S'(t) < 0$ for all $t \ge 0$.

This family of distributions includes mixtures of exponentials when all $A_i$ are positive but is not restricted to this case (some $A_i$ can be negative). The memoryless exponential distribution corresponds to the case $N = 1$, while cases with $N > 1$ correspond to phase-type distributions with memory. Focusing on multi-exponential distributions is not restrictive, as any phase-type distribution can be approximated with arbitrary precision by distributions of this type. Furthermore, as will be shown in the next section, finite Markov chains generically lead to multi-exponential distributions.

## 2.2 Finite state continuous time Markov chains

In this subsection, we introduce phase-type distributions using Markov chain representations and demonstrate the conditions under which we obtain the multi-exponential distributions described by (1). We assume that the reader is familiar with the basic theory of finite-state continuous-time Markov chains and refer to [6, 32] for excellent presentations of this subject. We also use the well known result from operational research that any phase-type distribution can be defined as the distribution of the absorption time of a finite-state, continuous-time Markov chain [11, 7].

To formalize the problem, we consider a continuous time Markov chain $M(t)$, where the first $N$ states are unobservable and the state $N + 1$ is the observed state.

The chain is defined by its generator $Q$, an $(N + 1) \times (N + 1)$ matrix with entries $Q_{ij} > 0$ representing the transition rate from the state $i$ to the state $j$ for $j \neq i$, and $Q_{ii} = - \sum_{j \neq i} Q_{ij}$ (zero row sum).
We are interested in the distribution of $T_s$, the first passage time to state $N+1$ starting from the state $s$ with $1 \le s \le 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 similar models. In the *model with return*, the chain returns instantly to state $s$ once reaching $N+1$ ($Q_{N+1,s} = \infty$) and the process restarts. Since $s$ is the unique return state, the observation timings form a renewal process with inter-event times following the distribution of $T_s$ (see [18] for an introduction to renewal processes). In the *model with absorption*, $N+1$ is absorbing and $Q_{N+1,i} = 0$ for all $i \ne N+1$. In this case, $T_s$ is the time to absorption starting from $s$.

The distribution of $T_s$ is the same for both models but suited to different applications: the absorption model applies when observations stop after the event, like in survival analysis or gene fixation, while the return model fits recurring events, such as transcription bursting or epidemiology. Even when the return model is contextually appropriate, we use the equivalent absorption model for distribution calculations due to technical considerations.

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

As examples of Markov models we have generated all the models having $N$ unobservable states and satisfying the two conditions:

* The matrix $Q$ has $2N-1$ independent, non-diagonal elements (kinetic parameters). This assumption is necessary to establish a one-to-one relation between the parameters of the $T_s$ distribution and the model's kinetic parameters, ensuring a unique solution to the inverse problem.

* Only the state $N$ leads to $N+1$, a simplifying assumption that fits the application discussed and aids in presentation. General results without this assumption are presented in the Subsection 3.7.

For $N=3$, up to permutations of the states, there are nine models that satisfy these assumptions. We label them as M1 to M9. Among these models, only five are ergodic, i.e. any state can be reached from any other state through transitions of the Markov chain. We will see in Section 3.6 that ergodic models are interesting because they can yield a unique solution to the inverse problem.

With the absorption assumption, the generator of the Markov chain model M9 that has $N=3$, represented in the Figure 2 is

$$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}. \quad (2)$$

For the calculation of the distribution of $T_s$ it is convenient to introduce the state probabilities $X_i(t) = P[M(t) = i | M(0) = s]$ for $1 \le i \le N+1$
The variables $X_i(t)$ satisfy the following system of linear differential equations (the master equation):

$$ \frac{d\mathbf{X}}{dt} = \mathbf{Q}^T \mathbf{X}(t), \qquad (3) $$

with the initial conditions $X_n(0) = \delta_{n,s}$, where $\delta_{n,s}$ is the Kronecker symbol and $\mathbf{Q}^T$ (called dual generator) is the transpose of the generator matrix $\mathbf{Q}$, and $\mathbf{X} = (X_1, \dots, X_N, X_{N+1})^T$. Because the last column of $\mathbf{Q}^T$ is zero ($N+1$ is absorbing), the variables $X_1, \dots, X_N$ satisfy an autonomous ODE system. Indeed, let $\tilde{\mathbf{Q}}$ be the $N \times N$ matrix obtained by eliminating the last line and the last column of $\mathbf{Q}^T$. For the example considered, we have

$$ \tilde{\mathbf{Q}} = \begin{pmatrix} -k_1 & 0 & k_3 \\ 0 & -k_2 & k_4 \\ k_1 & k_2 & -(k_3 + k_4 + k_5) \end{pmatrix}. \qquad (4) $$

Then $\tilde{\mathbf{X}}(t) = (X_1(t), \dots, X_N(t))^T$ satisfies

$$ \frac{d\tilde{\mathbf{X}}}{dt} = \tilde{\mathbf{Q}}\tilde{\mathbf{X}}(t), \qquad (5) $$

with initial conditions $X_i(0) = \delta_{i,s}$, where $\delta_{i,s}$ is the Kronecker delta.

Unless stated otherwise, we assume that the state $N+1$ can be reached only from state $N$. Thus, the remaining variable $X_{N+1}$ satisfies

$$ \frac{dX_{N+1}}{dt} = Q_{N,N+1}X_N(t), \qquad (6) $$

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

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

$$ \lambda_i \neq \lambda_j, \text{ for all } 1 \le i \neq j \le N, \text{ and } \max_{1 \le i \le N} (\lambda_i) < 0. \qquad (7) $$

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

Considering that the condition (7) is satisfied, the solutions of (5), (6) read

$$ \tilde{\mathbf{X}}(t) = \sum_{i=1}^{N} C_i \mathbf{u}_i e^{\lambda_i t} \quad \text{and} \quad X_{N+1}(t) = Q_{N,N+1} \sum_{i=1}^{N} \frac{C_i u_{i,N}}{\lambda_i} (e^{\lambda_i t} - 1), \qquad (8) $$

where $\lambda_i$ and $\mathbf{u}_i = (u_{i1}, \dots, u_{i,N})^T$ for $1 \le i \le N$ are eigenvalues and eigenvectors of $\tilde{\mathbf{Q}}$, respectively.
Because $N + 1$ is absorbing,

$$X_{N+1}(t) = \mathbb{P}[M(t) = N + 1 | M(0) = s] = \mathbb{P}[T_s \le t]$$

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}[T_s > t] = 1 - X_{N+1}(t) = \sum_{i=1}^{N} A_i e^{\lambda_i t}, \quad (9)$$

where

$$A_i = - \frac{Q_{N,N+1} C_i u_{i,N}}{\lambda_i}. \quad (10)$$

The constants $A_i$ satisfy

$$\sum_{i=1}^{N} A_i = 1, \quad (11)$$

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

This shows that in the non-degenerate case the distribution of $T_s$ is multi-exponential with $2N - 1$ independent parameters $\lambda_1, \dots, \lambda_N$ and $A_1, \dots, A_{N-1}$.

<figure><img src="image_3.png" /><figcaption>Figure 2: Models with $N = 3$ and return to the state 3. We show only those that are ergodic (every state can be reached from any other state through a directed path) and not related by symmetries (permutations of states 1, 2, 3). All these models can generate exactly the same phase-type distribution.</figcaption></figure>

## 3 Inverse problem

The inverse problem involves assuming that the survival function, and thus the parameters $\lambda_1, \dots, \lambda_N$, and $A_1, \dots, A_{N-1}$, are known. These can be estimated from data using least squares or maximum likelihood [15, 42, 13, 30]. The goal is to identify which Markov chains can represent this phase-type distribution. A number of general results are known for this problem. We recall here the following, important result, that follows from Theorem 12 in [11]:

**Proposition 2.** A multi-exponential phase-type distribution with $N$ exponentials can be represented by a Markov chain model of order $N$ (the order is defined as the number of unobservable states).
In this section we introduce the problem of computing the transition rate parameters, which are the non-zero, non-diagonal elements of the matrix Q of a Markov chain of order N, from the 2N − 1 parameters of the phase-type distribution. We are interested in models where this problem could have a unique solution, therefore the matrix Q has only 2N − 1 non-zero, non-diagonal elements. We range these elements in the 2N − 1 dimensional vector k. After reordering, we place the element Q$_{N,N}$$_{+1}$ in the last position of k, namely

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

We show below that the inverse problem consists in solving 2N − 1 polynomial equations for k.

3.1 Vieta’s formulas

The eigenvalues λ$_{1}$, . . . , λ$_{N}$ are the roots of the the characteristic polyno- mial of ˜ Q ( k ) , defined as

$$
P(\lambda) = \det(\lambda \mathbf{I} - \tilde{Q}(\mathbf{k})) = \lambda^N + a_{N-1}(\mathbf{k})\lambda^{N-1} + \dots + a_1(\mathbf{k})\lambda + a_0(\mathbf{k}). \quad (12)
$$

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

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

$$
\begin{align*}
L_1 &= \sum_{i=1}^{N} \lambda_i = -a_{N-1}(\mathbf{k}), \\
L_2 &= \sum_{i<j} \lambda_i \lambda_j = a_{N-2}(\mathbf{k}), \\
&\vdots \\
L_N &= \lambda_1 \lambda_2 \dots \lambda_N = (-1)^N a_0(\mathbf{k}). \tag{13}
\end{align*}
$$

3.2 Eigenvector equations

The amplitude parameters A$_{1}$, A$_{2}$, . . . , A$_{N}$ of the survival function occur in (10) together with eigenvector components u$_{i,N}$ and solution coefficients C$_{i}$ . We need to relate the latter to the transition rate parameters k .

First, we solve the eigenvector equation

$$
(\lambda \mathbf{I} - \tilde{Q}(\mathbf{k})) \mathbf{u} = 0. \tag{14}
$$

We look for solutions of (14) of the form

$$
\mathbf{u}(\lambda, \mathbf{k}) = (u_1(\lambda, \mathbf{k}), \dots, u_N(\lambda, \mathbf{k})),
$$

with u$_{s}$ ( λ, k ) = 1. In the subsection 3.6 we will see for which models this choice is possible.

Because equations (14) have integer coefficients in λ, k , the eigenvector components u$_{i}$ ( λ, k ) are rational functions of λ and k .
The initial conditions satisfied by the variables $X_i$ and (8) provide a linear system of equations for the constants $C_i$:

$$ \sum_{j=1}^{N} u_i(\lambda_j, \mathbf{k}) C_j = \delta_{i,s}, \quad 1 \le i \le N. \qquad (15) $$

Because of the non-degeneracy condition (7), the system (15) has a unique solution $C_i(\lambda, \mathbf{k})$, $1 \le i \le N$, where $\lambda = (\lambda_1, \dots, \lambda_N)$. The solutions $C_i(\lambda, \mathbf{k})$ are rational functions of $\lambda$, and $\mathbf{k}$.

From (10) we obtain $N-1$ independent equations for the transition rates $\mathbf{k}$:

$$ -k_{2N-1} u_N(\lambda_i, \mathbf{k}) C_i(\lambda, \mathbf{k}) = A_i \lambda_i, \quad 1 \le i \le N-1, \qquad (16) $$

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

The inverse problem is defined by the system of $2N-1$ equations formed by (13) and (16). When this system has solutions, the transition rates $\mathbf{k}$ can be expressed as functions of the survival function parameters $\lambda_i$ and $A_i$, $1 \le i \le N$.

### 3.3 Illustrative example

As an example let us consider again the model M9 with a return state $s = 3$, represented in Figure 2. For this model the characteristic polynomial is

$$ P(\lambda) = \lambda^3 + (k_1+k_2+k_3+k_4+k_5)\lambda^2 + (k_1k_2+k_1k_4+k_2k_3+k_1k_5+k_2k_5)\lambda + k_1k_2k_5. $$

The system (14) has the solution:

$$ \mathbf{u}(\lambda, \mathbf{k}) = \left( \frac{k_3}{k_1 + \lambda}, \frac{k_3}{k_2 + \lambda}, 1 \right). $$

The coefficients $(C_1, C_2, C_3)$ satisfy

$$ \begin{aligned} C_1 \frac{k_3}{k_1 + \lambda_1} + C_2 \frac{k_3}{k_1 + \lambda_2} + C_3 \frac{k_3}{k_1 + \lambda_3} &= 0, \\ C_1 \frac{k_3}{k_2 + \lambda_1} + C_2 \frac{k_3}{k_2 + \lambda_2} + C_3 \frac{k_3}{k_2 + \lambda_3} &= 0, \\ C_1 + C_2 + C_3 &= 1. \end{aligned} \qquad (17) $$

It follows that

$$ \begin{aligned} C_1 &= \frac{\lambda_1^2 + (k_1 + k_2)\lambda_1 + k_1 k_2}{(\lambda_1 - \lambda_2)(\lambda_1 - \lambda_3)}, \\ C_2 &= \frac{(k_1 + \lambda_2)(k_2 + \lambda_2)}{(\lambda_2 - \lambda_1)(\lambda_2 - \lambda_3)}, \\ C_3 &= \frac{(k_1 + \lambda_3)(k_2 + \lambda_3)}{(\lambda_3 - \lambda_1)(\lambda_3 - \lambda_2)}. \end{aligned} \qquad (18) $$

The inverse problem for the model M9 consists of solving the following system of five equations:

$$ -(k_1 + k_2 + k_3 + k_4 + k_5) = L_1, \qquad (19) $$
$$k_1 k_2 + k_1 k_4 + k_2 k_3 + k_1 k_5 + k_2 k_5 = L_2, \quad (20)$$

$$-k_1 k_2 k_5 = L_3, \quad (21)$$

$$-k_5 \frac{\lambda_1^2 + (k_1 + k_2)\lambda_1 + k_1 k_2}{(\lambda_1 - \lambda_2)(\lambda_1 - \lambda_3)} = A_1 \lambda_1, \quad (22)$$

$$-k_5 \frac{(k_1 + \lambda_2)(k_2 + \lambda_2)}{(\lambda_2 - \lambda_1)(\lambda_2 - \lambda_3)} = A_2 \lambda_2. \quad (23)$$

## 3.4 Symmetrized systems

There is a difference between systems (13) and (16). System (13) is entirely expressed using elementary symmetric polynomials in $\lambda_i$, whereas there is no obvious symmetry in system (16). This difference is clearly visible in the example of model M9: equations (19), (20), and (21) are symmetrized, while equations (22) and (23) are not symmetric in $\lambda_i$ and $A_i$.

We show here that the system formed by (13) and (16) 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 [19]) that are valid for any distinct $N$ numbers $\lambda_i$ ($1 \le i \le N$):

$$\sum_{i=1}^{N} \lambda_i^k \prod_{\substack{j=1 \\ j \neq i}}^{N} \frac{1}{\lambda_i - \lambda_j} = \begin{cases} 0, & \text{if } k < N-1 \\ 1, & \text{if } k = N-1 \\ h_{k-N+1}(\lambda), & \text{if } k > N-1 \end{cases} \quad (24)$$

where $h_{k-N+1}(\lambda)$ is the complete symmetric polynomial of degree $k-N+1$ in $N$ variables and $\lambda = (\lambda_1, \dots, \lambda_N)$. More precisely, the complete symmetric polynomials in the $N$ variables $\lambda_1, \dots, \lambda_N$ have the form

$$h_1(\lambda) = \sum_{1 \le i_1 \le N}^{N} \lambda_{i_1}, \quad h_2(\lambda) = \sum_{1 \le i_1 \le i_2 \le N}^{N} \lambda_{i_1} \lambda_{i_2}, \dots, \\ h_m(\lambda) = \sum_{1 \le i_1 \le \dots \le i_m \le N} \lambda_{i_1} \lambda_{i_2} \dots \lambda_{i_m}, \dots \quad (25)$$

In order to relate (16) and (24) we first relate the coefficients $C_i$ to eigenvalues.

Denote by $D_{si}(\lambda, \mathbf{k})$ the determinant of the $(N-1) \times (N-1)$ sub-matrix of the matrix $R(\lambda, \mathbf{k}) = \lambda \mathbf{I} - \tilde{Q}(\mathbf{k})$ obtained by deleting its $s$-th row and $i$-th column, for $1 \le i \le N$. Because elements of $R(\lambda, \mathbf{k})$ are linear combinations with integer coefficients of $\lambda$ and components of $\mathbf{k}$, $D_{si}(\lambda, \mathbf{k}) \in \mathbb{Z}[\lambda, \mathbf{k}]$, where $\mathbb{Z}[\lambda, \mathbf{k}]$ is the ring of polynomials with integer coefficients in $\lambda$, and $\mathbf{k}$.

From the definition of $D_{si}(\lambda, \mathbf{k})$, it follows for $i = s$ that $\deg_{\lambda}(D_{ss}(\lambda, \mathbf{k})) = N-1$ and that the (leading) coefficient of $D_{ss}(\lambda, \mathbf{k})$ at the monomial $\lambda^{N-1}$ equals 1, while for $1 \le i \neq s \le N$ we have $\deg_{\lambda}(D_{si}(\lambda, \mathbf{k})) \le N-2$. Thus, we have

$$D_{ss}(\lambda, \mathbf{k}) =: \lambda^{N-1} + c_{N-2}(\mathbf{k})\lambda^{N-2} + \dots + c_1(\mathbf{k})\lambda + c_0(\mathbf{k}), \quad (26)$$
$$ D_{si}(\lambda, \mathbf{k}) =: d_{i,N-2}(\mathbf{k})\lambda^{N-2} + d_{i,N-3}(\mathbf{k})\lambda^{N-3} + \dots + d_{i,1}(\mathbf{k})\lambda + d_{i,0}(\mathbf{k}), \quad (27) $$
if $i \neq s$.

**Lemma 3.** Assume that $\lambda_i \neq \lambda_j$ for $1 \le i \neq j \le N$ and that $D_{ss}(\lambda_j, \mathbf{k}) \neq 0$ and $u_s(\lambda, \mathbf{k}) = 1$ for $1 \le j \le N$. Then

$$ u_i(\lambda, \mathbf{k}) = (-1)^{|s-i|} \frac{D_{si}(\lambda, \mathbf{k})}{D_{ss}(\lambda, \mathbf{k})}, \quad (28) $$

and

$$ C_j = \frac{D_{ss}(\lambda_j, \mathbf{k})}{\prod_{1 \le l \ne j \le N} (\lambda_j - \lambda_l)}, \quad 1 \le j \le N, \quad (29) $$

where $u_i(\lambda, \mathbf{k})$ are the eigenvector components, and $C_j$ form the unique solution of the system

$$ \sum_{1 \le j \le N} C_j = 1, \\ \sum_{1 \le j \le N} u_i(\lambda_j, \mathbf{k}) C_j = 0, \quad 1 \le i \ne s \le N. $$

*Proof.* If $D_{ss}(\lambda_j, \mathbf{k}) \neq 0$, then, using Cramer's rule for the system (14) with $u_s(\lambda, \mathbf{k}) = 1$, we find

$$ \mathbf{u}(\lambda, \mathbf{k}) = ((-1)^{|s-1|} D_{s1}(\lambda, \mathbf{k}), \dots, (-1)^{|s-N|} D_{sN}(\lambda, \mathbf{k})) / D_{ss}(\lambda, \mathbf{k}). $$

Therefore, it holds

$$ \sum_{1 \le j \le N} C_j = \sum_{1 \le j \le N} \frac{D_{ss}(\lambda_j, \mathbf{k})}{\prod_{1 \le l \ne j \le N} (\lambda_j - \lambda_l)} = 1, $$

while

$$ \sum_{1 \le j \le N} u_i(\lambda_j, \mathbf{k}) C_j = \sum_{1 \le j \le N} \frac{(-1)^{|s-i|} D_{si}(\lambda_j, \mathbf{k})}{\prod_{1 \le l \ne j \le N} (\lambda_j - \lambda_l)} = 0, \quad 1 \le i \ne s \le N, $$

where $C_j, 1 \le j \le N$ are taken from (29)) due to (24) and to the fact that the degree in $\lambda$ of $D_{si}(\lambda, \mathbf{k})$ is $N-1$ and is smaller than $N-1$ for $i=s$ and $i \ne s$, respectively. $\square$

The structure of the equations (16) suggests a natural way to symmetrize them. For $1 \le k \le N-1$ we define $S_k$ to be the sum of $A_1\lambda_1^k, \dots, A_N\lambda_N^k$. Using (16), (28) and (29) 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{(-1)^{|s-N|+1} D_{sN}(\lambda_i, \mathbf{k})}{\prod_{1 \le j \ne i \le N} (\lambda_i - \lambda_j)}. \quad (30) $$

We can distinguish two cases:
i) $s = N$. In this case, according to (26) $D_{sN}$ has degree $N-1$ in $\lambda$. Using (30) and the Jacobi-Trudi identities (24) we find the following $N-1$ symmetrized equations:

$$S_1 = \sum_{i=1}^{N} A_i \lambda_i = -k_{2N-1},$$

$$S_2 = \sum_{i=1}^{N} A_i \lambda_i^2 = -k_{2N-1}(h_1(\lambda) + c_{N-2}(k)),$$

$$S_{N-1} = \sum_{i=1}^{N} A_i \lambda_i^{N-1} = -k_{2N-1}(h_{N-2}(\lambda) + c_{N-2}(k)h_{N-3}(\lambda) + \dots + c_2(k)h_1(\lambda) + c_1(k)). \quad (31)$$

The inverse problem is thus equivalent to solving the symmetrized equations in (13) and (31) and computing the transition rates $k$ as functions of the symmetric polynomials $L_1, \dots, L_N$ and $S_1, \dots, S_{N-1}$. We should note that the complete symmetric polynomials $h_1(\lambda), \dots, h_{N-2}(\lambda)$ can be expressed using the elementary symmetric polynomials $L_1, \dots, L_{N-2}$, that is $h_1(\lambda) = L_1, h_2(\lambda) = L_1^2 - L_2$ and so on.

ii) $s < N$. In this case, according to (27) $D_{sN}$ has degree at most $N-2$. Using (30) and (24) it follows that

$$S_1 = \sum_{i=1}^{N} A_i \lambda_i = 0. \quad (32)$$

This means that the inverse problem is not well posed in this case. The $2N-1$ parameters $\lambda_1, \dots, \lambda_N, A_1, \dots, A_{N-1}$ are no longer independent, they need to satisfy the constraint (32). Furthermore, we can no longer determine all the transition rate parameters $k$ uniquely. There remain $2N-2$ equations that can be used to constrain the $2N-1$ transition rate parameters

$$S_2 = (-1)^{|s-N|+1} k_{2N-1} d_{N,N-2}(k),$$

$$S_3 = (-1)^{|s-N|+1} k_{2N-1} (d_{N,N-2}(k)h_1(\lambda) + d_{N,N-3}(k)),$$

$$S_k = (-1)^{|s-N|+1} k_{2N-1} (d_{N,N-2}(k)h_{k-2}(\lambda) + d_{N,N-3}(k)h_{k-3}(\lambda) + \dots + d_{N,N-k}(k)),$$

$$S_{N-1} = (-1)^{|s-N|+1} k_{2N-1} (d_{N,N-2}(k)h_{N-3}(\lambda) + \dots + d_{N,2}(k)h_1(\lambda) + d_{N,1}(k)). \quad (33)$$

In this case the symmetrized system made of (13), (33) can be used to compute the transition rates $k$ as functions of the symmetric polynomials $L_1, \dots, L_N$ and $S_2, \dots, S_{N-1}$ and of one or several indeterminate transition rates.
While the examples in this paper assume $N = s$, the case $N \neq s$ is also relevant in practice. In this case, additional data such as state occupancy probabilities and lifetimes help to constrain the kinetic parameters by eliminating the remaining indeterminate variables. Symmetrized equations provide further constraints that must be solved alongside those from the supplementary data.

## 3.5 Symmetrized system for the illustrative example

For the illustrative example of the model M9, the matrix $R(\lambda, k)$ reads

$$R(\lambda, k) = \begin{pmatrix} \lambda + k_1 & 0 & -k_3 \\ 0 & \lambda + k_2 & -k_4 \\ -k_1 & -k_2 & \lambda + (k_3 + k_4 + k_5) \end{pmatrix}.$$

Considering that the return state is $s = 3$ we have

$$D_{33}(\lambda, k) = \det \begin{pmatrix} \lambda + k_1 & 0 \\ 0 & \lambda + k_2 \end{pmatrix} = (\lambda + k_1)(\lambda + k_2),$$

and the symmetrized equations

$$S_1 = -k_5,$$

$$S_2 = -k_5(L_1 + k_1 + k_2).$$

If the return state is $s = 1$, we have

$$D_{33}(\lambda, k) = \det \begin{pmatrix} 0 & \lambda + k_2 \\ -k_1 & -k_2 \end{pmatrix} = k_1(\lambda + k_2),$$

and the symmetrized equations

$$S_1 = 0,$$

$$S_2 = -k_5(L_1 + k_1).$$

Finally, if the return state is $s = 2$, we have

$$D_{33}(\lambda, k) = \det \begin{pmatrix} \lambda + k_1 & 0 \\ -k_1 & -k_2 \end{pmatrix} = -k_2(\lambda + k_1),$$

and the symmetrized equations

$$S_1 = 0,$$

$$S_2 = k_5(L_1 - k_2).$$
<figure><img src="image_4.png" /><figcaption>Figure 3: Non ergodic models with return. In both examples the state 2 is not in the strongly connected component of the return state 3. a) the state 2 is absorbing, so the return time can be infinite with non-zero probability, meaning that the distribution is not of the phase-type. b) the state 2 cannot be reached from the return state 3, allowing the ergodic chain obtained by eliminating state the state 2 to be used to compute the phase-type distribution.</figcaption></figure>

## 3.6 Solvable models

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

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

(C1) $N+1$ is reachable only from $N$.

(C2) From $N+1$ there is only one return arc, towards $s$, where $1 \le s \le N$.

C1 is a simplifying assumption that is satisfied in a number of applications, including the biological application discussed in this paper. All the results in this section remain valid, with slight modifications, even without this condition (see Subsection 3.7). The condition C2 is needed in order to have a phase-type distribution.

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 if $v \preceq w, w \preceq v$. Equivalence classes of nodes are called strongly connected components (SCC). If $G$ consists of a single SCC then we call *G strongly connected*. A Markov chain model is called *ergodic* if and only if its transition graph is strongly connected, meaning that every node can be reached from any other node through a directed path. Let us consider the SCC of the return state $s$, i.e. the smallest strongly connected sub-graph containing $s$. Non ergodic models contain nodes $w$ outside this SCC. Two situations may arise for such nodes, as represented in Figure 3. In Figure 3a), $s \preceq w$ but the converse is not true. With non-zero probability, the chain starting in $s$ never returns to $s$. This situation does not lead to a phase-
type distribution and will not be considered. In Figure 3b), $w$ can not be reached from $s$. In such a case, an ergodic model with fewer states, limited to the SCC of $s$, generates the phase-type distribution. For these reasons we can restrict our analysis to ergodic models only. We would like to know for which ergodic models the inverse problem has a unique solution, eventually up to transformation by discrete symmetries.

**Definition 4.** We say that an model is *solvable* if the following solutions are satisfied:

(i) The transition graph is strongly connected.

(ii) The model has $2N - 1$ transition rate parameters (non-zero, non-diagonal elements of $\tilde{Q}$).

(iii) The system made by the equations (13) and (16) or equivalently the system made by the equations (13) and (31) or the system (13),(33) 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 necessary condition.

**Proposition 5.** A *solvable model necessarily satisfies* $s = N$.

*Proof.* If $s \neq N$, then (32) holds, making the survival function parameters independent. The inverse problem then involves at most $2N - 2$ independent equations for the $2N - 1$ kinetic parameters. $\square$

**Remark 6.** In Subsection 3.7 we relax this necessary condition to $s \in \text{Pred}(N+1)$, where $\text{Pred}(N+1)$ is the set of predecessors of $N+1$ in the transition graph.

Ergodic models also satisfy the following property that has been used in condition (7) and is needed for writing the solution (8).

**Proposition 7.** All ergodic models satisfy $\max_{1 \le i \le N} \lambda_i < 0$, where $\lambda_i$ are the eigenvalues of $\tilde{Q}(k)$.

*Proof.* Suppose that some $\lambda_i = 0$. Then we have

$$ \tilde{Q}(k)u(\lambda_i, k) = 0. $$

Summing all the equations of the above system term by term, we obtain $k_{2N-1}u_N(\lambda_i, k) = 0$. This implies that $u_N(\lambda_i, k) = 0$. Therefore, $u(\lambda_i, k)$ also satisfies

$$ \tilde{Q}^{\text{red}} u(\lambda_i, k) = 0, $$

where $\tilde{Q}^{\text{red}}$ is obtained from $\tilde{Q}(k)$ by setting $k_{2N-1} = 0$.

Note that $\tilde{Q}^{\text{red}}$ satisfies the zero column-wise sum and is the dual generator of a reduced Markov chain obtained by deleting the state $N+1$. Therefore $u(\lambda_i, k)/\sum_{i=1}^N u_i(\lambda_i, k)$ is a steady state probability distribution of the reduced Markov chain with states $\{1, \dots, N\}$ and dual generator $\tilde{Q}^{\text{red}}$.

The full model, and consequently the reduced chain, are ergodic. Or, any ergodic Markov chain has a unique steady state distribution in which all states have non-zero probabilities, which contradicts $u_N(\lambda_i, k) = 0$. $\square$
The following Proposition guarantees that the system of equations (16), that are used to define the inverse problem, can be obtained for all ergodic models.

**Proposition 8.** Consider that $G$ is strongly connected. Then the equation $(\lambda I - \tilde{Q}(k))u = 0$ has solutions $u(\lambda, k) = (u_1(\lambda, k), \dots, u_N(\lambda, k))$ with $u_i(\lambda, k) \neq 0$ for all $1 \le i \le N$.

*Proof.* With the notations of the Subsection 3.2, the polynomial $D_{ss}(\lambda, k)$ has leading term $\lambda^{N-1}$ and therefore is not identically zero (see (26)). Lemma 3 implies that there are eigenvector solutions with $u_s(\lambda, k) = 1$. Furthermore, for $i \neq s$, $u_i(\lambda, k)$ are not identically zero because, otherwise, $u_i(0, k)$ would also be zero. However, $u(0, k)$ is an eigenvector of $\tilde{Q}(k)$ for the eigenvalue $\lambda = 0$, representing a steady state of the Markov chain. Following the same argument as in the proof of Proposition 7, for ergodic chains, $u_i(0, k) \neq 0$ for all $1 \le i \le N$. $\square$

Figure 2 shows all the ergodic models with $N = 3$ satisfying the conditions C1 and C2 and having $2N - 1$ non-zero transition rate parameters. The numbering of the models stems from the fact that there are nine models that satisfy the conditions, and among them, only five are ergodic.

## 3.7 Symmetrized equations and solvability without the condition C1.

In this subsection, we lift the simplifying assumption that $N + 1$ is reachable only from $N$ and allow multiple states to lead to $N$. The $2N - 1$ independent kinetic parameters of the model are non-zero, non-diagonal elements of $Q$.

Let $\text{Pred}(N+1)$ be the set of states that lead to $N+1$ (the predecessors of $N+1$ in the transition graph).

Then, (6), (10), and (30) become

$$ \frac{dX_{N+1}}{dt} = \sum_{j \in \text{Pred}(N+1)} Q_{j,N+1} X_j(t), $$

$$ A_i = -\frac{C_i}{\lambda_i} \sum_{j \in \text{Pred}(N+1)} Q_{j,N+1} u_{i,j}, $$

$$ S_k = \sum_{i=1}^{N} \frac{\sum_{j \in \text{Pred}(N+1)} (-1)^{|s-j|+1} Q_{j,N+1} \lambda_i^{k-1} D_{sj}(\lambda_i, Q)}{\prod_{1 \le l \ne i \le N} (\lambda_i - \lambda_l)}, $$

where $D_{sj}$ are defined as in Subsection 3.4.

The Vieta's formulas (13), which provide the first $N$ equations of the inverse problem, remain unchanged. The remaining $N - 1$ symmetrized equations resulting from the eigenvectors are now more complex.

Like in Subsection 3.4, we have two cases:
i) $s \notin \text{Pred}(N+1)$. In this case $S_1 = 0$, and the argument in Proposition 5 holds: the model is not solvable. The remaining $N-2$ symmetrized equations resulting from eigenvectors read

$$S_2 = \sum_{j \in \text{Pred}(N+1)} (-1)^{|s-j|+1} Q_{j,N+1} d_{N,N-2}(\mathbf{Q}),$$

$$S_3 = \sum_{j \in \text{Pred}(N+1)} (-1)^{|s-j|+1} Q_{j,N+1} (d_{N,N-2}(\mathbf{Q})h_1(\boldsymbol{\lambda}) + d_{N,N-3}(\mathbf{Q})),$$

$\vdots$

$$S_k = \sum_{j \in \text{Pred}(N+1)} (-1)^{|s-j|+1} Q_{j,N+1} (d_{N,N-2}(\mathbf{Q})h_{k-2}(\boldsymbol{\lambda}) + d_{N,N-3}(\mathbf{Q})h_{k-3}(\boldsymbol{\lambda}) + \dots + d_{N,N-k}(\mathbf{Q})),$$

$\vdots$

$$S_{N-1} = \sum_{j \in \text{Pred}(N+1)} (-1)^{|s-j|+1} Q_{j,N+1} (d_{N,N-2}(\mathbf{Q})h_{N-3}(\boldsymbol{\lambda}) + \dots + d_{N,2}(\mathbf{Q})h_1(\boldsymbol{\lambda}) + \dots + d_{N,1}(\mathbf{Q})), \quad (34)$$

where $d_{N,k}(\mathbf{Q})$, $0 \le k \le N-2$ are defined as in Subsection 3.4.

ii) $s \in \text{Pred}(N+1)$. In this case we have $N-1$ equations resulting from the eigenvectors:

$$S_1 = -Q_{s,N+1},$$

$$S_2 = -Q_{s,N+1}(h_1(\boldsymbol{\lambda}) + c_{N-2}(\mathbf{Q})) + \sum_{j \in \text{Pred}(N+1), j \neq s} (-1)^{|s-j|+1} Q_{j,N+1} d_{j,N-2}(\mathbf{Q}),$$

$\vdots$

$$S_{N-1} = -Q_{s,N+1}(h_{N-2}(\boldsymbol{\lambda}) + c_{N-2}(\mathbf{Q})h_{N-3}(\boldsymbol{\lambda}) + \dots + c_2(\mathbf{Q})h_1(\boldsymbol{\lambda}) + c_1(\mathbf{Q})) + \sum_{j \in \text{Pred}(N+1), j \neq s} (-1)^{|s-j|+1} Q_{j,N+1} (d_{j,N-2}(\mathbf{Q})h_{N-3}(\boldsymbol{\lambda}) + d_{j,N-3}(\mathbf{Q})h_{N-4}(\boldsymbol{\lambda}) + \dots + d_{j,1}(\mathbf{Q})), \quad (35)$$

where $c_k(\mathbf{Q})$, $0 \le k \le N-2$ are defined as in Subsection 3.4.

# 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 and propose a method to compute the solutions of the inverse problem.

Consider now a Markov chain with $N+1$ states arranged in a line and reversibly connected, as illustrated in Figure 4. For this model

$$\tilde{Q}_{i+1,i} = k_i^+, \tilde{Q}_{i,i+1} = k_i^-, 1 \le i < N, \tilde{Q}_{i,i} = -k_i^+ - k_{i-1}^-, 1 \le i < N, \tilde{Q}_{N,N} = -k_N - k_{N-1}^+.$$
<figure><img src="image_5.png" /><figcaption>Figure 4: A solvable model with an arbitrarily large number of states : the unbranched chain model. For this model, the solution of the inverse problem can be computed symbolically for any number of states $N$.</figcaption></figure>

This model is a generalization, of arbitrary length $N$, of the model $M2$ represented in Figure 2 that has $N = 3$. To ensure that the model is solvable, as stated in Proposition 5, we consider the return state to be $s = N$.

The solution of the direct problem presented in the Section 2 provides an algebraic map<fnref n="1" />

$$f: \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). \tag{36}$$

The goal of this section is to prove that $f$ is invertible and that its inverse is a rational map. This also means that the unbranched chain model is solvable for any $N$. However, as we will see later, the inverse map is not a rational function for all solvable models; for instance inversion formulas may involve radicals. From this perspective, the unbranched chain model is special and, in some ways, simpler.

Let $C := (C_1, \dots, C_N)^T$ represent a vector, where $C_i$ are solutions of (15).

Note that (15), (31) imply that

$$UC = (0, \dots, 0, 1)^T, \tag{37}$$

$$k_N = - \sum_{1 \le j \le N} A_j \lambda_j, \tag{38}$$

where $U$ is the $N \times N$ matrix with the columns $u(\lambda_1, k), \dots, u(\lambda_N, k)$ and $\lambda_1, \dots, \lambda_N$ are the eigenvalues of $\tilde{Q}(k)$.

The following lemmas will provide an algorithmic approach to constructing the solution to the inverse problem.

**Lemma 9.** For a suitable rational function $G_{N-i}$ with rational coefficients and with a denominator $k_{N-i+1}^+ \cdots k_{N-1}^+ k_{N-i}^- \cdots k_{N-1}^- k_N$ it holds

$$k_{N-i}^{+} = G_{N-i}(A_1, \dots, A_{N-1}, k_{N-i+1}^{+}, \dots, k_{N-1}^{+} k_{N-i}^{-}, \dots, k_{N-1}^{-}, k_N, \lambda_1, \dots, \lambda_N), \\ \text{for } 1 \le i \le N-1.$$

*Proof.* For each eigenvalue $\lambda$ of the matrix $\tilde{Q}(k)$ it holds

$$k_{N-1}^{+} u_{N-1}(\lambda, k) = k_N + k_{N-1}^{-} + \lambda, \quad k_n^{+} u_n(\lambda, k) = (k_{n+1}^{+} + k_n^{-} + \lambda) u_{n+1}(\lambda, k) - k_{n+1}^{-} u_{n+2}(\lambda, k),$$

<footnote n="1">¹although the rate parameters are real positive numbers, for formal reasons we define this map on complex numbers.</footnote>
for $0 \le n \le N - 2$, considering that $u_0(\lambda, \mathbf{k}) = 0$. This provides by recursion on $N - n$ a rational function $g_n$ with rational coefficients such that

$$k_n^+ u_n(\lambda, \mathbf{k}) = g_n(k_{n+1}^+, \dots, k_{N-1}^+, k_n^-, \dots, k_{N-1}^-, k_N, \lambda). \quad (39)$$

Moreover, the denominator of $g_n$ equals $k_{n+1}^+ \cdots k_{N-1}^+$.
Then for $i \ge 0$ it holds

$$\tilde{Q}^i UC = U \cdot \operatorname{diag}(\lambda_1^i, \dots, \lambda_N^i) C = \tilde{Q}^i (0, \dots, 0, 1)^T, \quad (40)$$

where $\operatorname{diag}(\lambda_1^i, \dots, \lambda_N^i)$ denotes a diagonal matrix.

The $(N-i)$-th coordinate of the middle vector in (40) equals

$$\sum_{1 \le j \le N} u_{N-i}(\lambda_j, \mathbf{k}) \lambda_j^i (-A_j \lambda_j / k_N). \quad (41)$$

The same $(N-i)$-th coordinate of the right vector in (40) equals

$$k_{N-i}^{-} \cdots k_{N-1}^{-} \quad (42)$$

Observe that the $j$-th coordinate of the right vector in (40) vanishes for $1 \le j < N - i$.

Multiplying both sides of (39) for $n = N - i$, $\lambda = \lambda_j$ by $\lambda_j^i (-A_j \lambda_j / k_N)$ and summing them up over $1 \le j \le N$, we obtain that

$$k_{N-i}^{+} k_{N-i}^{-} \cdots k_{N-1}^{-} = G_{N-i,0} (A_1, \dots, A_{N-1}, k_{N-i+1}^{+}, \dots, k_{N-1}^{+}, k_{N-i}^{-}, \dots, k_{N-1}^{-}, k_N, \lambda_1, \dots, \lambda_N) \quad (43)$$

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

$$k_{N-i+1}^{+} \cdots k_{N-1}^{+} k_N$$

taking into account the equality of (41) and of (42). $\square$

**Lemma 10.** For an appropriate rational function $P_{N-i}$ with rational coefficients and with a denominator $k_{N-i+1}^{+} \cdots k_{N-1}^{+} k_{N-i+1}^{-} \cdots k_{N-1}^{-} k_N$ it holds

$$k_{N-i}^{-} = P_{N-i}(A_1, \dots, A_{N-1}, k_{N-i+1}^{+}, \dots, k_{N-1}^{+}, k_{N-i+1}^{-}, \dots, k_{N-1}^{-}, k_N, \lambda_1, \dots, \lambda_N).$$

*Proof.* The $(N - i + 1)$-th coordinate of the middle vector in (40) equals

$$\sum_{1 \le j \le N} u_{N-i+1}(\lambda_j, \mathbf{k}) \lambda_j^i (-A_j \lambda_j / k_N). \quad (44)$$

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

$$-k_{N-i}^{-} \cdots k_{N-1}^{-} + p_{N-i}(k_{N-i+1}^{+}, \dots, k_{N-1}^{+}, k_{N-i+1}^{-}, \dots, k_{N-1}^{-}, k_N) \quad (45)$$

for an appropriate polynomial $p_{N-i}$ with integer coefficients.

Multiplying both sides of (39) for $n = N - i + 1$, $\lambda = \lambda_j$ by $\lambda_j^i (-A_j \lambda_j / k_N)$ and summing them up over $1 \le j \le N$, we obtain that

$$-k_{N-i+1}^{+} k_{N-i}^{-} \cdots k_{N-1}^{-} + h_{N-i}(k_{N-i+1}^{+}, \dots, k_{N-1}^{+}, k_{N-i+1}^{-}, \dots, k_{N-1}^{-}, k_N) = P_{N-i,0}(A_1, \dots, A_{N-1}, k_{N-i+2}^{+}, \dots, k_{N-1}^{+}, k_{N-i+1}^{-}, \dots, k_{N-1}^{-}, k_N, \lambda_1, \dots, \lambda_N)$$

for a suitable rational function $P_{N-i,0}$ with rational coefficients and with a denominator $k_{N-i+2}^{+} \cdots k_{N-1}^{+} k_N$ taking into account the equality of (44) and of (45) (while the both latter multiplied by $k_{N-i+1}^{+}$). $\square$
Applying alternatingly Lemma 10 and Lemma 9 for $i = 1, \dots, N - 1$ consecutively, we conclude with the following main result of this section.

**Theorem 11.** There is an algorithm which produces rational functions $F_n^+, F_n^-$, $1 \le n \le 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 12.** (i) Together with (38) Theorem 11 assures the inverse rational map to $f$ (see (36)) 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 \le i < j \le N$;

(ii) the proof of Theorem 11 provides an algorithm which produces explicitly rational functions $F_n^+, F_n^-$ by recursion on $N - n$.

**Remark 13.** In fact, for any model (not necessarily, the unbranched chain model elaborated in this section) one can consider a rational map similar to the one defined by (36). Thus, $A_1, \dots, A_{N-1}, \lambda_1, \dots, \lambda_{N-1}$ are rational functions in $k_n^+, k_n^-, 1 \le n < N, k_N$. Therefore, there is a field 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)$$

where $\mathbb{C}(A_1, \dots, A_{N-1}, \lambda_1, \dots, \lambda_{N-1})$ denotes the field generated by $A_1, \dots, A_{N-1}, \lambda_1, \dots, \lambda_{N-1}$ over the field of complex numbers. It is known (see e.g. [38], 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 11 states that for the unbranched chain model the degree equals 1. This provides a rigorous algebraic reformulation of the statement that the inverse problem is in some ways simpler for the unbranched chain model.

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 computer algebra algorithm that decomposes systems of polynomial equations and inequations into simpler systems. These can be more easily tackled by iteratively solving univariate polynomial equations. We apply this technique to compute the solutions of the inverse problem. The method is illustrated by solving the inverse problem for all ergodic models with $N = 3$ but can also be applied to larger $N$. Additionally, the method helps determine whether a model is solvable.

## 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 [28] 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\} \quad (46) $$

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 $\text{Sol}(\mathcal{S})$ of the algebraic system (46) is defined as the set of all $\bar{\boldsymbol{x}} = (\bar{x}_1, \dots, \bar{x}_n) \in \mathbb{C}^n$ satisfying the equations and inequations of $\mathcal{S}$, that is as

$$ \text{Sol}(\mathcal{S}) = \{\bar{\boldsymbol{x}} \in \mathbb{C}^n \mid p_i(\bar{\boldsymbol{x}}) = 0, q_j(\bar{\boldsymbol{x}}) \neq 0 \text{ for all } 1 \le i \le r \text{ and } 1 \le j \le s\}. $$

To fix ideas, we will use terminology from algebraic geometry, but this can be skipped without any loss by non-specialist readers. In algebraic geometry, Zariski closed sets, also called varieties, are sets of solutions of systems of polynomial equations. Zariski open sets are complements of closed sets (thus sets of solutions of systems of polynomial inequations) and Zariski locally closed sets are intersections between open and closed sets.

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

$$ \{\bar{\boldsymbol{x}} \in \mathbb{C}^n \mid p_1(\bar{\boldsymbol{x}}) = 0, \dots, p_r(\bar{\boldsymbol{x}}) = 0\} \text{ and } \{\bar{\boldsymbol{x}} \in \mathbb{C}^n \mid q_1(\bar{\boldsymbol{x}}) \cdots q_s(\bar{\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 $\text{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 $\text{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 $\text{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 $\text{ld}(p(\boldsymbol{x}))$ in $p(\boldsymbol{x})$ is called the initial of $p(\boldsymbol{x})$ which we denote by $\text{init}(p(\boldsymbol{x}))$. The separant $\text{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 14.** Let $\mathcal{S}$ be the algebraic system of (46). Then $\mathcal{S}$ is called a simple algebraic system with respect to a ranking, if the following conditions are satisfied:

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

$$ \text{card}(\{\text{ld}(p_1(\boldsymbol{x})), \dots, \text{ld}(p_r(\boldsymbol{x})), \text{ld}(q_1(\boldsymbol{x})), \dots, \text{ld}(q_s(\boldsymbol{x}))\} \setminus \{1\}) = r+s. $$

This property is called triangularity.

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 $\text{init}(p(\boldsymbol{x})) = 0$ has no solution in $\text{Sol}(\mathcal{S})$. We call this property non-vanishing initials.
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 $\text{sep}(p(\boldsymbol{x})) = 0$ has no solution in $\text{Sol}(S)$. This is called square-freeness.

The advantage of a simple algebraic system $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 $\bar{x}_1 \in \mathbb{C}$ of $p(\boldsymbol{x})$ can be chosen as the first coordinate of a solution $\bar{x}$ of $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 $\bar{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 $\bar{x}_1$ for $x_1$ in this equation or inequation and obtain so a univariate polynomial in $x_2$. Condition 2 of Definition 14 guarantees that the degree of the so obtained polynomial is independent of the choice of $\bar{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 $\bar{x}_2 \in \mathbb{C}$ for the second coordinate of a solution of $S$. An iteration of this process yields successively a solution $\bar{x} = (\bar{x}_1, \dots, \bar{x}_n) \in \mathbb{C}^n$ of $S$. Moreover, any solution of the simple algebraic system $S$ can be obtained by this process.

**Definition 15.** A Thomas decomposition of an algebraic system $S$ as in (46) consists of finitely many simple algebraic systems $S_1, \dots, S_m$ such that $\text{Sol}(S)$ is the disjoint union of $\text{Sol}(S_1), \dots, \text{Sol}(S_m)$.

It was proved by Thomas in [44, 43] that any algebraic system has a Thomas decomposition which is in general not unique. A Thomas decomposition can be determined algorithmically (see [2]) and there is an implementation in MAPLE. A description of the implementation can be found in [3, 20].

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 $k > v$ implies that each simple algebraic system returned by the Thomas decomposition yields solutions for the variables $k$ which are only valid for solutions of the variables $v$ satisfying equation and inequation conditions in $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 11. We are going to determine the solutions of the algebraic system of equations

$$ S = \begin{cases} -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, \\ 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 \end{cases} $$

To this end we compute with MAPLE an algebraic Thomas decomposition of $S$ with respect to the ranking

$$ k_1 > k_3 > k_2 > k_4 > k_5 > v. $$

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

$$ \begin{aligned} S_1 = & \{ 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 \\ & + (-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 \\ & + (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, \\ & (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, \\ & -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, \\ & 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 \}. \end{aligned} $$

One easily checks that the last three inequations do not involve any variable $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 \quad (47) $$

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

$$ \begin{aligned} k_1 = & \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_2^2 S_2) L_2 + (S_1^3 - S_1 S_2) L_3)} \\ & (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 \\ & + (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), \end{aligned} $$

$$ k_3 = \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)}, $$
$$k_2 = \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},$$

$$k_4 = \frac{S_1^2 - S_2}{S_1},$$

$$k_5 = -S_1,$$

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

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

## 5.3 A Thomas decomposition for model 3

The symmetrized algebraic system for model 3 is

$$S = \begin{cases} -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, \\ 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 \end{cases}$$

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

$$k_1 > k_2 > k_3 > k_4 > k_5 > v.$$

MAPLE returns 10 simple systems $S_1, \dots, S_{10}$. One easily checks by comparing the number of equations only involving the variables $v$ (see appendix A.2 for the remaining simple systems $S_2, \dots, S_{10}$), that the first simple system

$$S_1 = \begin{cases} 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, \\ 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 \end{cases}$$

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

$$L_1 S_1 S_2 - L_2 S_1^2 + L_3 S_1 - S_2^2 = 0, \quad S_1 \neq 0, \quad S_2 \neq 0. \quad (48)$$

Solving again the remaining equations successively for $k_5, k_4, k_3, k_2, k_1$, we obtain the solutions

$$k_5 = -S_1$$

$$k_4 = \frac{S_1^2 - S_2}{S_1}$$

$$k_3 \text{ arbitrary in } \mathbb{C} \setminus \{0\},$$

$$k_2 = \frac{L_3}{k_3 S_1},$$
$$k_1 = \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 (48) 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.6 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

$$S = \begin{cases} -S_1 = k_5, & -S_2 = k_5 (k_1 + k_2 + k_3 + L_1), \quad L_1 = -k_1 - k_2 - k_3 - k_4 - k_5, \\ 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, \quad L_3 = -k_1 k_3 k_5 - k_2 k_3 k_5 \end{cases}.$$

Using the ranking

$$k_1 > k_2 > k_3 > k_4 > k_5 > v$$

the MAPLE implementation of algebraic Thomas decomposition returns 21 simple systems (see subsection A.3 of the appendix for the systems $S_2, \dots, S_{21}$). The first simple system

$$S_1 = \begin{cases} 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, \\ 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, \quad k_3^2 S_1 + (L_1 S_1 - S_2) k_3 + L_3 = 0, \\ k_4 S_1 - S_1^2 + S_2 = 0, \quad k_5 + S_1 = 0, \quad L_1^2 S_1^2 - 2 L_1 S_1 S_2 - 4 L_3 S_1 + S_2^2 \neq 0, \quad L_3 \neq 0, \\ S_1^3 - S_1 S_2 \neq 0, \quad S_2 \neq 0 \end{cases}.$$

is the most generic one. Indeed, all inequations appearing in $S_1$ involve only the variables $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, \quad L_3 \neq 0, \quad S_1^3 - S_1 S_2 \neq 0, \quad S_2 \neq 0,$$

and so the solutions for $k$ are valid for all $\bar{v} \in \mathbb{C}^5$ of the Zariski open subset defined by these inequations. We can now successively solve the remaining equations in $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 - 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$. 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

$$k_5 = -S_1,$$

$$k_4 = \frac{S_1^2 - S_2}{S_1},$$

$$k_3^{1,2} = - \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},$$
$$k_2 = \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},$$

$$k_1 = \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, \quad L_3 \neq 0, \quad S_1^3 - S_1 S_2 \neq 0, \quad S_2 \neq 0. \tag{49}$$

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 (49). According to the Section 3.6 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

$$S = \begin{cases} -S_1 = k_5, & -S_2 = k_5(k_1 + k_2 + L_1), \quad L_1 = -k_1 - k_2 - k_3 - k_4 - k_5, \\ 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 \end{cases}$$

We compute the algebraic Thomas decomposition of $S$ with respect to the ranking

$$k_1 > k_3 > k_4 > k_2 > k_5 > v$$

and obtain from MAPLE twelve simple systems $S_1, \dots, S_{12}$. One easily checks by comparing the number of equations which have leader one of the variables of $v$, that the first simple system

$$S_1 = \begin{cases} 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, \\ 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, \\ 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 \end{cases}$$

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

$$k_5 = -S_1,$$

$$k_2^{1,2} = - \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},$$

$$k_4 = \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},$$

$$k_3 = - \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},$$

$$k_1 = \frac{L_3}{S_1 k_2^{1,2}}$$

of $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, \quad L_3 \neq 0, \quad S_1 \neq 0. \tag{50}$$
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 (50). According to the Section 3.6 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

$$ S = \{-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, \\ L_3 = -k_1k_2k_5, L_2 = k_1k_2 + k_1k_4 + k_1k_5 + k_2k_3 + k_2k_5\} \quad (51) $$

for model 9 with respect to the ranking

$$ k_3 > k_1 > k_4 > k_2 > k_5 > v. $$

The MAPLE implementation returns 15 simple systems $S_1, \dots, S_{15}$, where the most general simple system is

$$ S_1 = \{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, \\ -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, \\ 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, \\ 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\}. $$

Similar as for model 4 and 8 one determines the solutions

$$ \begin{align*} k_5 &= -S_1, \\ k_2^{1,2} &= -\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}, \\ k_4 &= \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} \end{align*} $$

$$ k_1 = -\frac{k_2^{1,2} S_1 + L_1 S_1 - S_2}{S_1}, $$

$$ k_3 = -\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)} \quad (52) $$

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, \quad L_3 \neq 0, \quad S_1 \neq 0. \quad (53) $$

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 (53). According to the Section 3.6 this means that the model M9 is solvable.
# 6 Model degeneracy effects

Several examples were used elsewhere [35] to illustrate the model degeneracy of the phase-type distribution. Indeed, given a phase-type distribution, one can have two or more distinct representations in terms of the same chain but with different rate parameters, and a Markov chain may provide representations of all distributions represented by a different chain. In this section, we explore the degeneracy of multi-exponential phase-type distributions and illustrate this phenomenon with examples significant in biological applications.

## 6.1 Discrimination between variant models

For a fixed number of states $N$, there are several variant models that predict exactly the same phase-type distribution and fit the data 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. Among these models, $M2$, $M4$, $M8$, and $M9$ are solvable. One can not discriminate between any two of the four solvable models by using only time-to-event data. The non-solvable model $M3$ does not generate any multi-exponential phase-type distribution, but only those distributions whose parameters satisfy (48). Additionally, the inverse problem has an infinite set of solutions for the model $M3$. 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 a 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.

To identify model variants experimentally, we aim to evaluate the potential of additional data, beyond time-to-event measurements, to differentiate between models. Two natural candidates for these data are the steady-state occupation probabilities of the Markov chain states and the lifetime of each state, defined as the average time spent in a state before transitioning. Depending on the application, experimental methods are available to estimate these quantities. We will discuss some of these methods in the Section 7 dedicated to transcriptional bursting.

## 6.2 Symmetries and classification of solvable models

Solvable models that generate the same phase-type distribution must satisfy certain common conditions, which we outline in this section. These
conditions allow us to identify mappings that relate the parameters of one
solvable model to those of another,

Let us consider the matrix $\tilde{Q}^{\text{red}}$ that is obtained from $\tilde{Q}(\mathbf{k})$ by setting $k_{2N+1}$ to zero. This matrix is the dual generator of a reduced Markov chain describing the transitions between the states $1, \dots, N$. Then the steady state probabilities of the states $1, \dots, N$ can be computed as solutions of the system:

$$ 
\begin{gathered}
\tilde{Q}^{\text{red}} \begin{pmatrix} p_1(\mathbf{k}) \\ p_2(\mathbf{k}) \\ \vdots \\ p_N(\mathbf{k}) \end{pmatrix} = 0, \\
p_1(\mathbf{k}) + p_2(\mathbf{k}) + \dots + p_N(\mathbf{k}) = 1.
\end{gathered}
\quad (54) $$

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

$$ T_i = -1/\tilde{Q}_{ii}^{\text{red}} = 1/\sum_{j \ne i} \tilde{Q}_{ij}^{\text{red}}. \quad (55) $$

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

<table><caption>Table 1: Reciprocal lifetimes $T_i^{-1}$ and steady state occupation probability $P_3$ for the solvable models $M2$, $M4$, $M8$, $M9$.</caption><thead><tr><th>Model</th><th>$T_1^{-1}$</th><th>$T_2^{-1}$</th><th>$T_3^{-1}$</th><th>$P_3$</th></tr></thead><tbody><tr><td>$M2$</td><td>$k_1 + k_2$</td><td>$k_3$</td><td>$k_4$</td><td>$\frac{k_2 k_3}{k_1 k_4 + k_2 k_3 + k_3 k_4}$</td></tr><tr><td>$M4$</td><td>$k_1 + k_2$</td><td>$k_3$</td><td>$k_4$</td><td>$\frac{k_1 k_3 + k_1 k_4 + k_2 k_3 + k_3 k_4}{k_1 k_2}$</td></tr><tr><td>$M8$</td><td>$k_1$</td><td>$k_2$</td><td>$(k_3 + k_4)$</td><td>$\frac{k_1 k_2 + k_1 k_3 + k_1 k_4 + k_2 k_3}{k_1 k_2 + k_1 k_4 + k_2 k_3}$</td></tr><tr><td>$M9$</td><td>$k_1$</td><td>$k_2$</td><td>$(k_3 + k_4)$</td><td>$\frac{k_1 k_2 + k_1 k_3 + k_1 k_4 + k_2 k_3}{k_1 k_2 + k_1 k_4 + k_2 k_3}$</td></tr></tbody></table>

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

**Proposition 16.** All solvable models that generate the same phase-type distribution have

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

(ii) the same value of the lifetime $T_N$ of the state $N$,

(iii) the same value of the steady state probability $p_N$ of the state $N$.

*Proof.* (i) is a direct consequence of (31). From (31) and (13) it follows:

$$ S_1 - \frac{S_2}{S_1} = -k_{2N-1} + a_{N-1}(\mathbf{k}) - c_{N-2}(\mathbf{k}). $$

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

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

$$S_1 - \frac{S_2}{S_1} = -\tilde{Q}_{NN}^{\text{red}} = 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 (51) invariant and therefore transforms any solution of the problem into another solution. As $\sigma_{M9}^2 = e$, where $e$ is the identity, the symmetry group is isomorphic to $\mathbb{Z}_2$. Solutions that are not symmetric are then two-fold degenerated.

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

$$ \begin{align} k_5 &= -S_1, \notag \\ k_1 &= -\frac{L_1 S_1 - S_2 \pm \sqrt{(L_1 S_1 - S_2)^2 - 4 L_3 S_1}}{2 S_1}, \notag \\ k_2 &= -\frac{L_1 S_1 - S_2 \mp \sqrt{(L_1 S_1 - S_2)^2 - 4 L_3 S_1}}{2 S_1}, \notag \\ k_3 &= \frac{1}{2 S_1} \left[ S_1^2 - S_2 \pm \frac{L_1 S_1^3 - S_1^2 S_2 - 2 L_2 S_1^2 + L_1 S_1 S_2 + 2 L_3 S_1 - S_2^2}{\sqrt{(L_1 S_1 - S_2)^2 - 4 L_3 S_1}}, \notag \\ k_4 &= \frac{1}{2 S_1} \left[ S_1^2 - S_2 \mp \frac{L_1 S_1^3 - S_1^2 S_2 - 2 L_2 S_1^2 + L_1 S_1 S_2 + 2 L_3 S_1 - S_2^2}{\sqrt{(L_1 S_1 - S_2)^2 - 4 L_3 S_1}} \right]. \tag{56} \end{align} $$

(56) 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 17.** There are bijective rational maps $f_{M,M'}: D \subset \mathbb{R}^5 \to D' \subset \mathbb{R}^5$, relating parameters $k$ of model M to parameters $k'$ of the model $M'$ (i.e. $k' = f_{M,M'}(k)$) where $M, M'$ are any pair in the list $\{M4, M8, M9\}$, and $D, D'$ are open domains containing all the real positive parameters resulting from the inverse problem of $M, M'$. $f_{M,M'}$ satisfy the following constraints

(i) The parameter $k_5$ is the same for all the models.

(ii) The lifetimes of the states 1, 2, 3 is the same for all the models.

(iii) Under steady state conditions, the probability to be in the state 3 is the same for all the models.
Model M with parameters $k$ has exactly the same phase-type distribution as model $M'$ with phase-type $k'$. 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 Proposition 16. The same Proposition 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 $k$ and $k'$. These equations, if they have solutions, define the map $f_{M,M'}$. 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'_i$ of M8, then using the Table 1 the conditions (ii), (iii) read

$$ 
\begin{aligned}
k_1 &= k'_1, \\
k_2 &= k'_2, \\
k_3 + k_4 &= k'_3 + k'_4, \\
\frac{k_1 k_2}{k_1 k_2 + k_1 k_4 + k_2 k_3} &= \frac{k'_1 k'_2}{k'_1 k'_2 + k'_1 k'_3 + k'_1 k'_4 + k'_2 k'_3}.
\end{aligned} 
\quad (57) $$

The solutions of system (57) read

$$ 
\begin{aligned}
k'_1 &= k_1, \\
k'_2 &= k_2, \\
k'_3 &= \frac{k_3(k_2 - k_1)}{k_2}, \\
k'_4 &= \frac{k_1 k_3 + k_2 k_4}{k_2},
\end{aligned} 
\quad (58) $$

and define the bijective rational map $f_{M9,M8} : \{k_1 > k_2 > 0, k_3 > 0, k_4 > 0\} \to \{k'_2 > k'_1 > 0, k'_3 > 0, k'_4 > \frac{k'_1 k'_3}{k'_2 - k'_1}\}$.

Similarly, solving

$$ 
\begin{aligned}
k_1 &= k'_1 + k'_2, \\
k_2 &= k'_3, \\
k_3 + k_4 &= k'_4, \\
\frac{k_1 k_2}{k_1 k_2 + k_1 k_4 + k_2 k_3} &= \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},
\end{aligned} 
\quad (59) $$

with respect to $k'$ leads to the equations

$$ 
\begin{aligned}
k'_1 &= (k_1 k_4 - k_2 k_4)/(k_3 + k_4), \\
k'_2 &= (k_1 k_3 + k_2 k_4)/(k_3 + k_4), \\
k'_3 &= k_2, \\
k'_4 &= k_3 + k_4,
\end{aligned} 
\quad (60) $$
that define the map $f_{M9,M4} : \{k_1 > k_2 > 0, k_3 > 0, k_4 > 0\} \to \{k'_1 > k'_2 > 0, k'_2 > k'_3 > 0, k'_4 > 0\}$. $\square$

**Remark 18.** The Proposition 17 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}(\mathbf{k})$ and $f_{M9,M8} \circ \sigma_{M9}(\mathbf{k})$ where $\mathbf{k}$ is a solution of M9.

**Remark 19.** The Proposition 17 allows to compute the symmetries of the models M8 and M4. The symmetry groups are isomorphic to $\mathbb{Z}_2$ and generated by $f_{M9,M8} \circ \sigma_{M9} \circ f_{M9,M8}^{-1}$, and $f_{M9,M4} \circ \sigma_{M9} \circ f_{M9,M4}^{-1}$, respectively.

**Remark 20.** Of course, there are also mappings that relate the parameters of models M4, M8, or M9 to those of model M2, resulting in the same phase-type distribution; however, these mappings are many-to-one. This is consistent with the remark in Section 4 that the degree of the field extension from the field generated by survival function parameters to the field generated by kinetic parameters is one for M2 and larger than one for the other models.
<figure><img src="image_6.png" /><figcaption>Figure 5: 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 allows 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.</figcaption></figure>
# 7 Application to transcriptional bursting models

## 7.1 Biological considerations and assumptions

Transcription is the fundamental molecular process that copies the DNA into RNA. The synthesis of RNA requires a dedicated molecular machinery, including an 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 5 A and Table 2). We classify these states in three categories: productive ON states that can initiate transcription, non-productive OFF states that are unable to initiate or resume transcription, and paused states PAUSE in which initiated transcription halts and can resume later or abort. The promoter starts transcription in the productive ON state,

<table><caption>Table 2: Possible biological interpretation of the states in the models *M2*, *M3*, *M4*, *M8*, *M9*. For *M4*, *M8*, *M9* pausing is always abortive.</caption><thead><tr><th colspan="5">Possible state types</th></tr><tr><th>Model</th><th>state 1</th><th>state 2</th><th>state 3</th><th>state 4</th></tr></thead><tbody><tr><td>M2</td><td>OFF</td><td>OFF</td><td>ON</td><td>EL</td></tr><tr><td>M3</td><td>OFF</td><td>OFF</td><td>ON</td><td>EL</td></tr><tr><td>M4</td><td>PAUSE /OFF</td><td>OFF</td><td>ON</td><td>EL</td></tr><tr><td>M8</td><td>PAUSE /OFF</td><td>OFF</td><td>ON</td><td>EL</td></tr><tr><td>M9</td><td>PAUSE /OFF</td><td>OFF/PAUSE</td><td>ON</td><td>EL</td></tr></tbody><thead><tr><th colspan="5">Possible biological interpretation of state types</th></tr></thead><tbody><tr><td>long OFF</td><td colspan="4">Enhancer-Promoter disconnection / Chromatin closing</td></tr><tr><td>short OFF</td><td colspan="4">Transcription Factor unbinding / Pre-initiation Complex disassembly</td></tr><tr><td>PAUSE</td><td colspan="4">Promoter proximal pausing / Intrinsic pausing / Transcriptional roadblock</td></tr></tbody></table>

when it can trigger several departures of RNAP molecules. After ON, the RNAP can eventually stop in PAUSE or commit to irreversible productive elongation 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 [36]. Using the tool BurstDeconv [13] we can deconvolve the live imaging signal and identify the time of each elongation event. This means that for each
transcription site we observe the sequence of EL states (Figure 5B). The other states of the promoter are unobservable.

Strictly speaking, transcription initiation is not synonymous with the start of processive elongation, as transcription initiation 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 a unique return state $s$, that can be either productive or non-productive. The second assumption is needed for applying the phase-type distribution formalism. The first assumption, not applying to the case of multiple ON states, was lifted in Subsection 3.7. Using this extension of our approach we can also address scenarios in which transcription occurs simultaneously on duplicate DNA templates, as for example newly replicated sister chromatids.

## 7.2 Discriminating between variant models using state probabilities and lifetimes.

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.

In this subsection we test the possibility to discriminate between solvable variant models using measurements of state lifetimes or probabilities. The probability that a given state is occupied can be estimated from experimental data such as single molecule DNA footprinting assays [24, 23].

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

Additionally, Proposition 16 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
<figure><img src="image_7.png" /><figcaption>Figure 6: Testing the discrimination capacity of $p_1$, $T_1$, and $T_2$: We generated 100,000 random parameter sets $(\lambda_1, \lambda_2, \lambda_3, A_1, A_2, A_3)$, where $A_1$ and $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]$. From these, 61,411 sets were retained for which the inverse problem yielded positive rate parameters $k_1, k_2, k_3, k_4, k_5$ in at least one model (M2, M4, M8, or M9). Discrimination capacity for $p_1$, $T_1$, and $T_2$ is measured by the variations $\Delta p$, $\Delta \log_{10}(T_1)$, and $\Delta \log_{10}(T_2)$ across models with the same phase-type distribution parameters, where larger variations indicate better discrimination. The figures show frequency histograms of these variations, with no discrimination if the variation is zero. Results show that $p_1$, $T_1$, and $T_2$ fail to discriminate in 11%, 16%, and 16% of the cases, respectively.</figcaption></figure>

variants.

We found that $p_1$ 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 6. We have used random phase-type distribution parameters spanning domain values observed in studies of transcriptional bursting of mammalian and fruit-fly genes [42, 37]. This suggests that $p_1$ 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 [42], where $T_1$ and $T_2$ are practically the same for all experimental conditions for the 3 variants, see Figure 7 below.
<figure><img src="image_8.png" /><figcaption>Figure 7: 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 [42]). 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.</figcaption></figure>

# 8 Conclusion and perspectives.

In contrast to statistical inference, which treats one model at a time, using symbolic solutions for the inverse phase-type distribution problem allows us to simultaneously infer all solvable models. Once the parameters of the phase-type distributions are estimated, our formulas provide the transition rate parameters for all solvable models. Thomas decomposition can be used to solve the inverse problem symbolically 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 inverse problem for phase-type distributions leads us to a classification of Markov chain models compatible with the data. Some models are solvable, while others are not. Among the solvable three-state models, we discovered relationships stemming from the discrete symmetries of the inverse problem. Extending these results to an arbitrary number of states requires different mathematical tools, which we will address in future work.

We identified examples of solvable models with an arbitrary number of states and established necessary conditions for their solvability. While sufficient conditions for solvability remain to be found, this task presents an exciting challenge for future research.

In practical applications, Thomas decomposition provides formal solutions of the inverse problem. However, if a model has a solution set with a dimension greater than zero, or if multiple models have solutions, degeneracy arises, preventing unique identification from the experimental phase-type distribution. In such cases, additional data, such as state occupancy probabilities or state lifetimes, is needed. To achieve unique solutions, equations (13) and (31) (for $s = N$ or (13) and (33) (for $s \neq N$) must be combined with equations (54) and (55) that relate state occupancies and lifetimes to kinetic parameters.

However, limitations may arise from the complexity of the problem. Although models and parameters can theoretically be distinguished, many combinations may produce similar results, leading to practical degeneracy and overfitting. Some practical recipes to penalize complexity and avoid overfitting have been proposed in [13].

We applied our approach to the study of transcriptional bursting in two different biological context, human cultured cells [42] and *Drosophila* developing embryos [37].

By mapping the parameters of the phase-type distribution to kinetic model parameters, we can identify mechanisms impacting transcriptional bursting from imaging datasets. This concept led to the development of BurstDeconv, a tool for analyzing transcriptional bursting data [13]. Utilizing Thomas decomposition, the computations presented here have enabled us to extend the range of models that BurstDeconv can analyze. Currently, BurstDeconv supports only the two-state random telegraph and the three-state *M*2 and *M*9 models. Additionally, we demonstrated that variant models indistinguishable by time-to-event data can be discriminated using state occupancy probabilities.

In contrast to BurstDeconv, which only addresses transcription models with a single ON state, the results in this paper also apply to models with multiple ON states. This broader approach enables the analysis of phenomena such as sister chromatids (homologous DNA sequences from replication that can undergo transcription) where transcription sites are considered to have multiple active ON states. In this case the symmetrized equations to solve are (35) (when the return state is an ON state) and (34) (when the return state is not an ON state).

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 [41, 4, 8], 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 [29].

## Availability

The implementation of the Thomas decomposition in Maple is available at https://github.com/Computational-Systems-Biology-LPHI/Thomas_Decomposition/.

## 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, and John Reinitz for critical reading the manuscript.

## References

[1] 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(1):243–262, 2021.

[2] T. Bächler, V.P. Gerdt, M. Lange-Hegermann, and D. Robertz. Algorithmic Thomas decomposition of algebraic and differential systems. *J. Symbolic Comput.*, 47(10):1233–1266, 2012.

[3] T. Bächler and M. Lange-Hegermann. *Algebraic Thomas and Differential Thomas: Thomas Decomposition for algebraic and differential systems*, 2008-2012. (https://www.art.rwth-aachen.de/cms/MATHB/Forschung/Mathematische-Software/~lqnwi/).

[4] Eugenia Basyuk, Florence Rage, and Edouard Bertrand. RNA transport from transcription to localized translation: a single molecule perspective. *RNA Biology*, 18(9):1221–1237, 2021.

[5] 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, Heloise Faure-Gautron, Mattias Mannervik, Melissa M. Harrison, Jean-Christopheandrau, Cyril Favard, Ovidiu Radulescu, and Mounia Lagha. The control of transcriptional memory by stable mi-rotic bookmarking. *Nature Communications*, 13(1):1176, Mar 2022.
[6] Albert T Bharucha-Reid. *Elements of the Theory of Markov Processes and their Applications*. Dover, Mineola, NY, 1997.

[7] Mogens Bladt. A Review on Phase-type Distributions and their Use in Risk Theory. *ASTIN Bulletin: The Journal of the IAA*, 35(1):145–161, 2005.

[8] Lauren A Blake, Ana De La Cruz, and Bin Wu. Imaging spatiotemporal translation regulation in vivo. In *Seminars in Cell & Developmental Biology*, volume 154, pages 155–164. Elsevier, 2024.

[9] Pavol Bokes, Alessandro Borri, Pasquale Palumbo, and Abhyudai Singh. Mixture distributions in a stochastic gene expression model with delayed feedback: a WKB approximation approach. *Journal of Mathematical Biology*, 81(1):343–367, 2020.

[10] C. Commault and J. P. Chemla. An invariant of representations of phase-type distributions and some applications. *Journal of Applied Probability*, 33(2):368–381, 1996.

[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(6):566–580, 2003. Publisher: Taylor & Francis eprint: https://doi.org/10.1080/0020717031000114986.

[12] Alexia Damour, Vera Slaninova, Ovidiu Radulescu, Edouard Bertrand, and Eugenia Basyuk. Transcriptional Stochasticity as a Key Aspect of HIV-1 Latency. *Viruses*, 15(9):1969, 2023.

[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(16):e88–e88, 2023.

[14] 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(1):5194, Dec 2018.

[15] Daniel Dufresne. Fitting combinations of exponentials to probability distributions. *Applied Stochastic Models in Business and Industry*, 23(1):23–48, 2007.

[16] Mark Fackrell. Modelling healthcare systems with phase-type distributions. *Health Care Management Science*, 12:11–26, 2009.

[17] MJ Faddy. A structured compartmental model for drug kinetics. *Biometrics*, pages 243–248, 1993.

[18] William Feller. *An introduction to probability theory and its applications, Volume 2*. John Wiley & Sons, Hoboken, NJ, 1966.

[19] William Fulton and Joe Harris. *Representation theory. A first course*, volume 129. Springer, New York, 1991.
[20] 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:202-215, 2019.

[21] Asger Hobolth, Arno Siri-Jégousse, and Mogens Bladt. Phase-type distributions in population genetics. *Theoretical Population Biology*, 127:16-32, 2019.

[22] Ola Hössjer, Günter Bechly, and Ann Gauger. Phase-Type Distribution Approximations of the Waiting Time Until Coordinated Mutations Get Fixed in a Population. In Sergei Silvestrov, Anatoliy Malyarenko, and Milica Rančić, editors, *Stochastic Processes and Applications*, pages 245-313, New York, 2018. Springer International Publishing.

[23] 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(12):5673-5706, 2021.

[24] Arnaud R Krebs. Studying transcription factor function in the genome at molecular resolution. *Trends in Genetics*, 37(9):798-806, 2021.

[25] Niraj Kumar and Rahul V. Kulkarni. Constraining the complexity of promoter dynamics using fluctuations in gene expression. *Physical Biology*, 17(1):015001, 2019. Publisher: IOP Publishing.

[26] Niraj Kumar, Abhyudai Singh, and Rahul V. Kulkarni. Transcriptional Bursting in Gene Expression: Analytical Results for General Stochastic Models. *PLOS Computational Biology*, 11(10):e1004292, 2015. Publisher: Public Library of Science.

[27] Drew M. LaMar, Peter Kemper, and Gregory D. Smith. Reduction of calcium release site models via moment fitting of phase-type distributions. *Physical Biology*, 8(2):026015, 2011.

[28] Markus Lange-Hegermann, Daniel Robertz, Werner M. Seiler, and Matthias Seiß. Singularities of algebraic differential equations. *Advances in Applied Mathematics*, 131:102266, 2021.

[29] 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(1):1-14, 2012.

[30] Yi Liu. *Non-Parametric Bayesian Inference with Application to System Biology*. PhD thesis, Department of Statistics, University of Chicago, 2022.

[31] Robert S. Maier. The algebraic construction of phase-type distributions. *Communications in Statistics. Stochastic Models*, 7(4):573-602, 1991. Publisher: Taylor & Francis eprint: https://doi.org/10.1080/15326349108807207.

[32] Sean P Meyn and Richard L Tweedie. *Markov chains and stochastic stability*. Springer Science & Business Media, New York, 2012.
[33] Jeffrey R. Moffitt and Carlos Bustamante. Extracting signal from noise: kinetic mechanisms from a Michaelis–Menten-like expression for enzymatic fluctuations. *The FEBS Journal*, 281(2):498–517, 2014. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/febs.12545.

[34] Marcel F Neuts. Probability distributions of phase type. *Liber Amicorum Prof. Emeritus H. Florin*, 1975.

[35] Colm Art O'Cinneide. On non-uniqueness of representations of phase-type distributions. *Communications in Statistics. Stochastic Models*, 5(2):247–259, 1989. Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/15326348908807108.

[36] Xavier Pichon, Mounia Lagha, Florian Mueller, and Edouard Bertrand. A growing toolbox to image gene expression in single cells: sensitive approaches for demanding challenges. *Molecular Cell*, 71(3):468–480, 2018.

[37] 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(1):4504, Jul 2021.

[38] Igor Shafarevich. *Basic Algebraic Geometry 1*. Springer, Heidelberg, 1972.

[39] Mohammad Soltani, Cesar A. Vargas-Garcia, Duarte Antunes, and Abhyudai Singh. Intercellular Variability in Protein Levels from Stochastic Expression and Noisy Cell Cycle Processes. *PLOS Computational Biology*, 12(8):e1004972. Publisher: Public Library of Science.

[40] 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(4):e0000017, 2022.

[41] 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(3):635–646, 2014.

[42] 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(1):4503, Jul 2021.

[43] J.M. Thomas. *Differential Systems*. Colloquium Publications XXI. American Mathematical Society, New York, 1937.

[44] J.M. Thomas. *Systems and Roots*. W. Byrd Press, Richmond, VA, 1962.
[45] Jinqiang Zhang, Aimin Chen, Huahai Qiu, Jiajun Zhang, Tian-hai Tian, and Tianshou Zhou. Exact results for gene-expression models with general waiting-time distributions. *Physical Review E*, 109(2):024119, 2024. Publisher: American Physical Society.

## A Appendix: List of simple Thomas decomposition systems

In the main text of the paper we have used the most generic simple Thomas decomposition systems to solve the inverse problem. These most generic systems contain a minimum number of equations in $v = (L_1, L_2, L_3, S_1, S_2)$ (no equation in $v$ for solvable models). The remaining Thomas decomposition systems contain more equations in $v$ and eventually some equations in $k$ are replaced by inequations. Systems that contain only equations in $v$ and also some equations in $v$, lead to simpler formulas for $k$, valid for special values of $v$ satisfying equations (belonging to Zariski closed or locally closed subsets). Systems that contain inequations in $k$ lead to formulas where not all $k_i$, $1 \le i \le N$ are determined.

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

$$ S_2 = \begin{cases} -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, & 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, \\ L_2 S_1 - L_3 \neq 0, S_1 \neq 0, S_2 = 0 \end{cases} $$

$$ S_3 = \begin{cases} 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, \\ 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 \end{cases} $$

$$ S_4 = \begin{cases} 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, \\ 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 \end{cases} $$

$$ S_5 = \begin{cases} 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, \\ S_1 \neq 0, S_2 = 0 \end{cases} $$

$$ S_6 = \begin{cases} 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, \\ k_2 \neq 0, k_5 = 0, L_3 = 0, S_1 = 0, S_2 = 0 \end{cases} $$

$$ S_7 = \begin{cases} 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, \\ L_2 \neq 0, L_3 = 0, S_1 = 0, S_2 = 0 \end{cases} $$

$$ S_8 = \begin{cases} 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 \end{cases} $$

$$ S_9 = \begin{cases} 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 \end{cases} $$

$$ S_{10} = \begin{cases} 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, \end{cases} $$
$S_1 = 0, S_2 = 0\},$

$$
S_{11} = \{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\}.
$$

A.2 The Remaining Simple Systems for Model 3

$$
S_2 = \{ 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, \\
k_5 + S_1 = 0, L_2 S_1 - L_3 = 0, S_1 \neq 0, S_2 = 0 \},
$$

$$
S_3 = \{ 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, \\
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 \},
$$

$$
S_4 = \{ 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, \\
S_1 \neq 0, S_2 = 0 \},
$$

$$
S_5 = \{ 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, \\
k_3 \neq 0, k_5 = 0, L_3 = 0, S_1 = 0, S_2 = 0 \},
$$

$$
S_6 = \{ 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, \\
L_2 \neq 0, L_3 = 0, S_1 = 0, S_2 = 0 \},
$$

$$
S_7 = \{ 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, \\
S_2 = 0 \},
$$

$$
S_8 = \{ 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 \},
$$

$$
S_9 = \{ 2k_1 + 2k_2 + L_1 = 0, k_3 = 0, 2k_4 + L_1 = 0, k_5 = 0, L_1^2 - 4L_2 = 0, L_2 \neq 0, \\
L_3 = 0, S_1 = 0, S_2 = 0 \},
$$

$$
S_{10} = \{ 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 \}.
$$

A.3 The Remaining Simple Systems for Model 4

$$
S_2 = \{ 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, \\
-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, \\
L_1 S_1 - S_2 \neq 0, L_3 = 0, S_1^3 - S_1 S_2 \neq 0, S_2 \neq 0 \},
$$

$$
S_3 = \{ 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, \\
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 \},
$$

$$
S_4 = \{ -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,
$$
$$ 
\begin{aligned}
& -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, \quad 2 k_3 S_1 + L_1 S_1 - S_2 = 0, \\
& k_4 S_1 - S_1^2 + S_2 = 0, \quad k_5 + S_1 = 0, \quad 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, \quad S_1^3 - S_1 S_2 \neq 0, \quad S_2 \neq 0
\end{aligned}
$$

$$ S_5 = \{ -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, \\ 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 \}, $$

$$ S_6 = \{ 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, \\ 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 \}, $$

$$ S_7 = \{ 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, \\ S_1 \neq 0, S_2 = 0 \}, $$

$$ S_8 = \{ 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, \\ L_3 = 0, S_1 \neq 0, S_2 = 0 \}, $$

$$ S_9 = \{ 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, \\ 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 \}, $$

$$ S_{10} = \{ 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, \\ S_1 \neq 0, S_2 = 0 \}, $$

$$ S_{11} = \{ 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, \\ 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, \\ 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 \}, $$

$$ S_{12} = \{ 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, \\ L_3 = 0, S_1^2 - S_2 = 0, S_2 \neq 0 \}, $$

$$ S_{13} = \{ 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, \\ L_3 = 0, S_1^2 - S_2 = 0, S_2 \neq 0 \}, $$

$$ S_{14} = \{ 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, \\ 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, \\ L_3 \neq 0, S_1^2 - S_2 = 0, S_2 \neq 0 \}, $$

$$ S_{15} = \{ 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, \\ S_1^2 - S_2 = 0, S_2 \neq 0 \}, $$

$$ S_{16} = \{ 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, $$
$$k_5 = 0, L_3 = 0, S_1 = 0, S_2 = 0\},$$

$$S_{17} = \{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, L_3 = 0, S_1 = 0, S_2 = 0\},$$

$$S_{18} = \{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\},$$

$$S_{19} = \{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\},$$

$$S_{20} = \{2k_1 + 2k_2 + L_1 = 0, 2k_3 + L_1 = 0, k_4 = 0, k_5 = 0, L_1^2 - 4L_2 = 0, L_2 \neq 0, L_3 = 0, S_1 = 0, S_2 = 0\},$$

$$S_{21} = \{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\}.$$

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

$$S_2 = \{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,$$

$$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,$$

$$L_1 S_1 - S_2 \neq 0, L_3 = 0, S_1 \neq 0\},$$

$$S_3 = \{(L_1 S_1 - S_2) k_1 + 2 L_3 = 0,$$

$$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,$$

$$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,$$

$$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\},$$

$$S_4 = \{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,$$

$$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\},$$

$$S_5 = \{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\},$$

$$S_6 = \{k_1 + k_3 + k_4 + k_2 + L_1 = 0,$$

$$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,$$

$$-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\},$$

$$S_7 = \{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,$$

$$L_1^2 - 4 L_2 \neq 0, L_2 \neq 0, L_3 = 0, S_1 = 0, S_2 = 0\},$$

$$S_8 = \{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\},$$

$$S_9 = \{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\},$$

$$S_{10} = \{2 k_1 + k_2 + L_1 = 0, L_1^2 - k_2^2 + 4 k_2 k_3 - 4 L_2 = 0,$$

$$-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\},$$
$$
\mathcal{S}_{11} = \begin{cases} 2k_1 + L_1 = 0, & 2k_3 + 2k_4 + L_1 = 0, & k_2 = 0, & k_5 = 0, & L_1^2 - 4L_2 = 0, & L_2 \neq 0, \\ L_3 = 0, & S_1 = 0, & S_2 = 0 \end{cases},
$$

$$
\mathcal{S}_{12} = \{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\}.
$$

A.5 The Remaining Simple Systems for Model 9

$$
\begin{align*}
S_2 &= \{-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, \\
    &\phantom{=} 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, \\
    &\phantom{=} L_1 S_1 - S_2 \neq 0, L_3 = 0, S_1 \neq 0\},
\end{align*}
$$

$$
\begin{align*}
S_3 &= \{ 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, \\
    &\quad - 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, \\
    &\quad L_3 = 0, S_1 \neq 0 \},
\end{align*}
$$

$$
\begin{align*}
S_4 &= \{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, \\
    &\quad 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, \\
    &\quad 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\},
\end{align*}
$$

$$
\begin{align*}
S_5 &= \{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, \\
    &\quad L_3 = 0, S_1^3 + S_1 S_2 \neq 0, S_2 \neq 0\},
\end{align*}
$$

$$
\begin{align*}
S_6 &= \{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, \\
    &\quad L_1 S_1^2 - 2 L_2 S_1 + 2 L_3 = 0, L_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\},
\end{align*}
$$

$$
S_7 = \{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, \\
S_2 = 0\},
$$

$$
\begin{align*}
S_8 &= \{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, \\
    &\quad 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\},
\end{align*}
$$

$$
S_9 = \{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, \\
S_1^2 + S_2 = 0, S_2 \neq 0\},
$$

$$
\mathcal{S}_{10} = \{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, \\
k_4 \neq 0, k_5 = 0, L_3 = 0, S_1 = 0, S_2 = 0\},
$$

$$
\mathcal{S}_{11} = \{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, \\
L_3 = 0, S_1 = 0, S_2 = 0\},
$$

$$
\mathcal{S}_{12} = \{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\},
$$
$$S_{13} = \{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\},$$

$$S_{14} = \begin{cases} 2k_3 + 2k_1 + L_1 = 0, & k_4 = 0, 2k_2 + L_1 = 0, k_5 = 0, L_1^2 - 4L_2 = 0, L_2 \neq 0, \\ L_3 = 0, & S_1 = 0, S_2 = 0 \end{cases},$$

$$S_{15} = \{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\}.$$