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

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.

Introduction

Continuous-time finite Markov chains have been largely used to model biological and medical phenomena involving multiple discrete states (Bharucha-Reid, 1997). 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 (Feller, 1966). 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 (Neuts, 1975), (see Fig. 1).

Phase-type distributions were used in biology and medicine for modeling transcription and translation bursting (Kumar et al., 2015; Soltani et al., 2016; Kumar and Kulkarni, 2019; Bokes et al., 2020; Tantale et al., 2021; Pimmett et al., 2021; Douaihy et al., 2023; Zhang et al., 2024), stochastic enzymatic reactions (Moffitt and Bustamante, 2014), ion channel dynamics (LaMar et al., 2011), drug kinetics (Faddy, 1993), population genetics (Hobolth et al., 2019; Hössjer et al., 2018), and public health issues (Fackrell, 2009; Liquet et al., 2012; Asanjarani et al., 2021; Stone et al., 2022).

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 (Asanjarani et al., 2021; Stone et al., 2022).

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 (Dufresne, 2007; Douaihy et al., 2023). 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 (Maier, 1991; Commault and Chemla, 1996), 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 (Dufourt et al., 2018; Tantale et al., 2021; Pimmett et al., 2021; Bellec et al., 2022; Douaihy et al., 2023; Damour et al., 2023).

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.

[See PDF for image]

Fig. 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

Phase-Type Distribution, Direct Problem

Statistical Definition

Consider the successive observations of a state, as illustrated in Fig. 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 (Commault and Mocanu, 2003). Alternatively, it can be defined as having a Laplace transform that is a rational function (see Theorem 1 in (Commault and Mocanu, 2003)). 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:

1

$$\begin{array}{c}\hfill S(t)=\mathbb{P}\left[T_s,>,t\right]=\sum _{i=1}^N{A}_{i}exp({\lambda }_{i}t),\end{array}$$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^{\prime }(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.

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 (Bharucha-Reid, 1997; Meyn and Tweedie, 2012) 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 (Commault and Mocanu, 2003; Bladt, 2005).

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 $\mathit{Q}$, an $(N+1)\times (N+1)$ matrix with entries $Q_{\mathit{ij}}>0$ representing the transition rate from the state i to the state j for $j\ne i$, and $Q_{\mathit{ii}}=-{\sum }_{j\ne i}{Q}_{\mathit{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 (Feller, 1966) 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 $\mathit{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 Sect. 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 Sect. 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 Fig. 2 is

2

$$\begin{array}{c}\hfill \mathit{Q}=\left(\begin{array}{cccc}-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{array}\right).\end{array}$$For the calculation of the distribution of $T_s$ it is convenient to introduce the state probabilities $X_i(t)=\mathbb{P}\left[M,(,t,),=,i,,|,,M,(,0,),=,s\right]$ for $1\le i\le N+1$. The variables $X_i(t)$ satisfy the following system of linear differential equations (the master equation):

3

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

4

$$\begin{array}{c}\hfill \stackrel{\mathbf{~}}{\mathit{Q}}=\left(\begin{array}{ccc}-k_1& 0& k_3\\ 0& -k_2& k_4\\ k_1& k_2& -(k_3+k_4+k_5)\end{array}\right).\end{array}$$Then $\stackrel{\mathbf{~}}{\mathit{X}}(t)={(X_1(t),\dots ,X_N(t))}^T$ satisfies

5

$$\begin{array}{c}\hfill \frac{d\stackrel{\mathbf{~}}{\mathit{X}}}{\mathit{dt}}=\stackrel{\mathbf{~}}{\mathit{Q}}\stackrel{\mathbf{~}}{\mathit{X}}(t),\end{array}$$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

6

$$\begin{array}{c}\hfill \frac{d{X}_{N+1}}{\mathit{dt}}=Q_{N,N+1}{X}_N(t),\end{array}$$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 $\stackrel{\mathbf{~}}{\mathit{Q}}$ satisfy the following non-degeneracy condition

7

$$\begin{array}{c}\hfill {\lambda }_i\ne {\lambda }_j,\text{for all}1\le i\ne j\le N,\text{and}\underset{1\le i\le N}{max}({\lambda }_i)<0.\end{array}$$

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 $\mathit{Q}$, ${\stackrel{\mathbf{~}}{\mathit{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

8

$$\begin{array}{c}\hfill \stackrel{\mathbf{~}}{\mathit{X}}(t)=\sum _{i=1}^N{C}_i{\mathit{u}}_i{e}^{{\lambda }_{i}t}\text{and}{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),\end{array}$$where ${\lambda }_i$ and ${\mathit{u}}_i={(u_{i1},\dots ,u_{i,N})}^T$ for $1\le i\le N$ are eigenvalues and eigenvectors of $\stackrel{\mathbf{~}}{\mathit{Q}}$, respectively.

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

9

$$\begin{array}{c}\hfill S(t)=\mathbb{P}\left[T_s,>,t\right]=1-X_{N+1}(t)=\sum _{i=1}^N{A}_i{e}^{{\lambda }_{i}t},\end{array}$$where

10

$$\begin{array}{c}\hfill A_i=-\frac{Q_{N,N+1}{C}_i{u}_{i,N}}{{\lambda }_i}.\end{array}$$The constants $A_i$ satisfy

11

$$\begin{array}{c}\hfill \sum _{i=1}^N{A}_i=1,\end{array}$$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}$.

[See PDF for image]

Fig. 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

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 (Dufresne, 2007; Tantale et al., 2021; Douaihy et al., 2023; Liu, 2022). 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 (Commault and Mocanu, 2003):

Proposition 1

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 $\mathit{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 $\mathit{Q}$ has only $2N-1$ non-zero, non-diagonal elements. We range these elements in the $2N-1$ dimensional vector $\mathit{k}$. After reordering, we place the element $Q_{N,N+1}$ in the last position of $\mathit{k}$, namely$$k_{2N-1}=Q_{N,N+1}.$$We show below that the inverse problem consists in solving $2N-1$ polynomial equations for $\mathit{k}$.

Vieta’s Formulas

The eigenvalues ${\lambda }_1,\dots ,{\lambda }_N$ are the roots of the the characteristic polynomial of $\stackrel{\mathbf{~}}{\mathit{Q}}(\mathit{k})$, defined as

12

$$\begin{array}{c}\hfill P(\lambda )=\text{det}(\lambda \mathit{I}-\stackrel{\mathbf{~}}{\mathit{Q}}(\mathit{k}))={\lambda }^N+a_{N-1}(\mathit{k}){\lambda }^{N-1}+\dots +a_1(\mathit{k})\lambda +a_0(\mathit{k}).\end{array}$$The coefficients $a_i(\mathit{k})$ of the characteristic polynomial are polynomials with integer coefficients on the transition rates $\mathit{k}$.

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

13

$$\begin{array}{cc}\hfill L_1=& \sum _{i=1}^N{\lambda }_i=-a_{N-1}(\mathit{k}),\hfill \\ \hfill L_2=& \sum _{i<j}{\lambda }_i{\lambda }_j=a_{N-2}(\mathit{k}),\hfill \\ \hfill & ⋮\hfill & \hfill \\ \hfill L_N=& {\lambda }_1{\lambda }_2\cdots {\lambda }_N={(-1)}^N{a}_0(\mathit{k}).\hfill \end{array}$$

Eigenvector Equations

The amplitude parameters $A_1,A_2,\dots ,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 $\mathit{k}$.

First, we solve the eigenvector equation

14

$$\begin{array}{c}\hfill (\lambda \mathit{I}-\stackrel{\mathbf{~}}{\mathit{Q}}(\mathit{k}))\mathit{u}=0.\end{array}$$We look for solutions of (14) of the form$$\mathit{u}(\lambda ,\mathit{k})=(u_1(\lambda ,\mathit{k}),\cdots ,u_N(\lambda ,\mathit{k})),$$with $u_s(\lambda ,\mathit{k})=1$. In the Sect. 3.6 we will see for which models this choice is possible.

Because equations (14) have integer coefficients in $\lambda ,\mathit{k}$, the eigenvector components $u_i(\lambda ,\mathit{k})$ are rational functions of $\lambda$ and $\mathit{k}$.

The initial conditions satisfied by the variables $X_i$ and (8) provide a linear system of equations for the constants $C_i$:

15

$$\begin{array}{c}\hfill \sum _{j=1}^N{u}_i({\lambda }_j,\mathit{k})C_j={\delta }_{i,s},1\le i\le N.\end{array}$$Because of the non-degeneracy condition (7), the system (15) has a unique solution $C_i(\mathit{\lambda },\mathit{k}),1\le i\le N$, where $\mathit{\lambda }=({\lambda }_1,\dots ,{\lambda }_N)$. The solutions $C_i(\mathit{\lambda },\mathit{k})$ are rational functions of $\mathit{\lambda }$, and $\mathit{k}$.

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

16

$$\begin{array}{c}\hfill -k_{2N-1}{u}_N({\lambda }_i,\mathit{k})C_i(\mathit{\lambda },\mathit{k})=A_i{\lambda }_i,1\le i\le N-1,\end{array}$$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 $\mathit{k}$ can be expressed as functions of the survival function parameters ${\lambda }_i$ and $A_i$, $1\le i\le N$.

Illustrative Example

As an example let us consider again the model M9 with a return state $s=3$, represented in Fig. 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_1{k}_2+k_1{k}_4+k_2{k}_3+k_1{k}_5+k_2{k}_5)\lambda +k_1{k}_2{k}_5.$$The system (14) has the solution:$$\mathit{u}(\lambda ,\mathit{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

17

$$\begin{array}{cc}\hfill 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,\hfill \\ \hfill 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,\hfill \\ \hfill C_1+C_2+C_3=& 1.\hfill \end{array}$$It follows that

18

$$\begin{array}{cc}\hfill C_1=& \frac{{\lambda }_1^2+(k_1+k_2){\lambda }_1+k_1{k}_2}{({\lambda }_1-{\lambda }_2)({\lambda }_1-{\lambda }_3)},\hfill \\ \hfill C_2=& \frac{(k_1+{\lambda }_2)(k_2+{\lambda }_2)}{({\lambda }_2-{\lambda }_1)({\lambda }_2-{\lambda }_3)},\hfill \\ \hfill C_3=& \frac{(k_1+{\lambda }_3)(k_2+{\lambda }_3)}{({\lambda }_3-{\lambda }_1)({\lambda }_3-{\lambda }_2)}.\hfill \end{array}$$The inverse problem for the model M9 consists of solving the following system of five equations:

19

$$\begin{array}{cc}\hfill -(k_1+k_2+k_3+k_4+k_5)=& L_1,\hfill \end{array}$$

20

$$\begin{array}{cc}\hfill k_1{k}_2+k_1{k}_4+k_2{k}_3+k_1{k}_5+k_2{k}_5=& L_2,\hfill \end{array}$$

21

$$\begin{array}{cc}\hfill -k_1{k}_2{k}_5=& L_3,\hfill \end{array}$$

22

$$\begin{array}{cc}\hfill -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,\hfill \end{array}$$

23

$$\begin{array}{cc}\hfill -k_5\frac{(k_1+{\lambda }_2)(k_2+{\lambda }_2)}{({\lambda }_2-{\lambda }_1)({\lambda }_2-{\lambda }_3)}=& A_2{\lambda }_2.\hfill \end{array}$$

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 (Fulton and Harris, 1991)) that are valid for any distinct N numbers ${\lambda }_i$ ($1\le i\le N$):

24

$$\begin{array}{c}\hfill \sum _{i=1}^N{\lambda }_i^k\prod _{\stackrel{j=1}{j\ne i}}^N\frac{1}{{\lambda }_i-{\lambda }_j}=\left\{\begin{array}{cc}0,\hfill & \text{if}k<N-1\hfill \\ 1,\hfill & \text{if}k=N-1\hfill \\ h_{k-N+1}(\mathit{\lambda }),\hfill & \text{if}k>N-1\hfill \end{array}\right)\end{array}$$where $h_{k-N+1}(\mathit{\lambda })$ is the complete symmetric polynomial of degree $k-N+1$ in N variables and $\mathit{\lambda }=({\lambda }_1,\dots ,{\lambda }_N)$. More precisely, the complete symmetric polynomials in the N variables ${\lambda }_1,\dots ,{\lambda }_N$ have the form

25

$$\begin{array}{cc}\hfill h_1(\mathit{\lambda })=& \sum _{1\le i_1\le N}^N{\lambda }_{i_1},h_2(\mathit{\lambda })=\sum _{1\le i_1\le i_2\le N}^N{\lambda }_{i_1}{\lambda }_{i_2},\cdots ,\hfill \\ \hfill h_m(\mathit{\lambda })=& \sum _{1\le i_1\le \cdots \le i_m\le N}{\lambda }_{i_1}{\lambda }_{i_2}\cdots {\lambda }_{i_m},\cdots \hfill \end{array}$$In order to relate (16) and (24) we first relate the coefficients $C_i$ to eigenvalues.

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

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

26

$$\begin{array}{ccc}\hfill D_{\mathit{ss}}(\lambda ,\mathit{k})& =:\hfill & \hfill {\lambda }^{N-1}+c_{N-2}(\mathit{k}){\lambda }^{N-2}+\dots +c_1(\mathit{k})\lambda +c_0(\mathit{k}),\end{array}$$

27

$$\begin{array}{ccc}\hfill D_{\mathit{si}}(\lambda ,\mathit{k})& =:\hfill & \hfill d_{i,N-2}(\mathit{k}){\lambda }^{N-2}+d_{i,N-3}(\mathit{k}){\lambda }^{N-3}+\dots +d_{i,1}(\mathit{k})\lambda +d_{i,0}(\mathit{k}),\\ \hfill & \text{if}i\ne s.\hfill \\ \hfill \end{array}$$

Lemma 1

Assume that ${\lambda }_i\ne {\lambda }_j$ for $1\le i\ne j\le N$ and that $D_{\mathit{ss}}({\lambda }_j,\mathit{k})\ne 0$ and $u_s(\lambda ,\mathit{k})=1$ for $1\le j\le N$. Then

28

$$\begin{array}{c}\hfill u_i(\lambda ,\mathit{k})={(-1)}^{|s-i|}\frac{D_{\mathit{si}}(\lambda ,\mathit{k})}{D_{\mathit{ss}}(\lambda ,\mathit{k})},\end{array}$$and

29

$$\begin{array}{c}\hfill C_j=\frac{D_{\mathit{ss}}({\lambda }_j,\mathit{k})}{{\prod }_{1\le l\ne j\le N}({\lambda }_j-{\lambda }_l)},1\le j\le N,\end{array}$$where $u_i(\lambda ,\mathit{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,\mathit{k})C_j=0,1\le i\ne s\le N.$$

Proof

If $D_{\mathit{ss}}({\lambda }_j,\mathit{k})\ne 0$, then, using Cramer’s rule for the system (14) with $u_s(\lambda ,\mathit{k})=1$, we find$$\mathit{u}(\lambda ,\mathit{k})=({(-1)}^{|s-1|}{D}_{s1}(\lambda ,\mathit{k}),\cdots ,{(-1)}^{|s-N|}{D}_{\mathit{sN}}(\lambda ,\mathit{k}))/D_{\mathit{ss}}(\lambda ,\mathit{k}).$$Therefore, it holds$$\sum _{1\le j\le N}{C}_j=\sum _{1\le j\le N}\frac{D_{\mathit{ss}}({\lambda }_j,\mathit{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,\mathit{k})C_j=\sum _{1\le j\le N}\frac{{(-1)}^{|s-i|}{D}_{\mathit{si}}({\lambda }_j,\mathit{k})}{{\prod }_{1\le l\ne j\le N}({\lambda }_j-{\lambda }_l)}=0,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_{\mathit{si}}(\lambda ,\mathit{k})$ is $N-1$ and is smaller than $N-1$ for $i=s$ and $i\ne s$, respectively.

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,\cdots ,A_N{\lambda }_N^k$. Using (16),(28) and (29) we obtain for this sum the formula

30

$$\begin{array}{c}\hfill 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}_{\mathit{sN}}({\lambda }_i,\mathit{k})}{{\prod }_{1\le j\ne i\le N}({\lambda }_i-{\lambda }_j)}.\end{array}$$We can distinguish two cases:

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

31

$$\begin{array}{cc}\hfill S_1=& \sum _{i=1}^N{A}_i{\lambda }_i=-k_{2N-1},\hfill \\ \hfill S_2=& \sum _{i=1}^N{A}_i{\lambda }_i^2=-k_{2N-1}(h_1(\mathit{\lambda })+c_{N-2}(\mathit{k})),\hfill \\ \hfill & ⋮\hfill & \hfill \\ \hfill S_{N-1}=& \sum _{i=1}^N{A}_i{\lambda }_i^{N-1}=-k_{2N-1}(h_{N-2}(\mathit{\lambda })\hfill \\ \hfill & +c_{N-2}(\mathit{k})h_{N-3}(\mathit{\lambda })+\dots +c_2(\mathit{k})h_1(\mathit{\lambda })+c_1(\mathit{k})).\hfill \end{array}$$The inverse problem is thus equivalent to solving the symmetrized equations in (13) and (31) and computing the transition rates $\mathit{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(\mathit{\lambda }),\dots ,h_{N-2}(\mathit{\lambda })$ can be expressed using the elementary symmetric polynomials $L_1,\dots ,L_{N-2}$, that is $h_1(\mathit{\lambda })=L_1$, $h_2(\mathit{\lambda })=L_1^2-L_2$ and so on.

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

32

$$\begin{array}{c}\hfill S_1=\sum _{i=1}^N{A}_i{\lambda }_i=0.\end{array}$$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 $\mathit{k}$ uniquely. There remain $2N-2$ equations that can be used to constrain the $2N-1$ transition rate parameters

33

$$\begin{array}{cc}\hfill S_2=& {(-1)}^{|s-N|+1}{k}_{2N-1}{d}_{N,N-2}(\mathit{k}),\hfill \\ \hfill S_3=& {(-1)}^{|s-N|+1}{k}_{2N-1}(d_{N,N-2}(\mathit{k})h_1(\mathit{\lambda })+d_{N,N-3}(\mathit{k})),\hfill \\ \hfill & ⋮\hfill & \hfill \\ \hfill S_k=& {(-1)}^{|s-N|+1}{k}_{2N-1}(d_{N,N-2}(\mathit{k})h_{k-2}(\mathit{\lambda })\hfill \\ \hfill & +d_{N,N-3}(\mathit{k})h_{k-3}(\mathit{\lambda })+\dots +d_{N,N-k}(\mathit{k})),\hfill \\ \hfill & ⋮\hfill & \hfill \\ \hfill S_{N-1}=& {(-1)}^{|s-N|+1}{k}_{2N-1}(d_{N,N-2}(\mathit{k})h_{N-3}(\mathit{\lambda })+\dots +d_{N,2}(\mathit{k})h_1(\mathit{\lambda })+d_{N,1}(\mathit{k})).\hfill \\ \hfill \end{array}$$In this case the symmetrized system made of (13),(33) can be used to compute the transition rates $\mathit{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\ne 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.

Symmetrized System for the Illustrative Example

For the illustrative example of the model M9, the matrix $R(\lambda ,\mathit{k})$ reads$$R(\lambda ,\mathit{k})=\left(\begin{array}{ccc}\lambda +k_1& 0& -k_3\\ 0& \lambda +k_2& -k_4\\ -k_1& -k_2& \lambda +(k_3+k_4+k_5)\end{array}\right).$$Considering that the return state is $s=3$ we have$$D_{33}(\lambda ,\mathit{k})=det\left(\begin{array}{cc}\lambda +k_1& 0\\ 0& \lambda +k_2\end{array}\right)=(\lambda +k_1)(\lambda +k_2),$$and the symmetrized equations$$\begin{array}{cc}\hfill S_1=& -k_5,\hfill \\ \hfill S_2=& -k_5(L_1+k_1+k_2).\hfill \end{array}$$If the return state is $s=1$, we have$$D_{33}(\lambda ,\mathit{k})=det\left(\begin{array}{cc}0& \lambda +k_2\\ -k_1& -k_2\end{array}\right)=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 ,\mathit{k})=det\left(\begin{array}{cc}\lambda +k_1& 0\\ -k_1& -k_2\end{array}\right)=-k_2(\lambda +k_1),$$and the symmetrized equations$$S_1=0,$$$$S_2=k_5(L_1-k_2).$$

[See PDF for image]

Fig. 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

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)\mid i,j\in V,i\ne j,Q_{\mathit{ij}}\ne 0\}$.

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

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

2. 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 Sect. 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⪯w$ if there is a path in G from v to w. We say that v, w are equivalent if $v⪯w,w⪯v$. Equivalence classes of nodes are called strongly connected components (SCC). If G consists of a single SCC then we call Gstrongly 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 Fig. 3. In Fig. 3a), $s⪯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 Fig. 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 1

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

1. The transition graph is strongly connected.

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

3. 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 2

A solvable model necessarily satisfies $s=N$.

Proof

If $s\ne 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.

Remark 2

In Sect. 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 3

All ergodic models satisfy ${max}_{1\le i\le N}{\lambda }_i<0$, where ${\lambda }_i$ are the eigenvalues of $\stackrel{\mathbf{~}}{\mathit{Q}}(\mathit{k})$.

Proof

Suppose that some ${\lambda }_i=0$. Then we have$$\stackrel{\mathbf{~}}{\mathit{Q}}(\mathit{k})\mathit{u}({\lambda }_i,\mathit{k})=0.$$Summing all the equations of the above system term by term, we obtain $k_{2N-1}{u}_N({\lambda }_i,\mathit{k})=0.$ This implies that $u_N({\lambda }_i,\mathit{k})=0$. Therefore, $\mathit{u}({\lambda }_i,\mathit{k})$ also satisfies$${\stackrel{\mathbf{~}}{\mathit{Q}}}^{\text{red}}\mathit{u}({\lambda }_i,\mathit{k})=0,$$where ${\stackrel{\mathbf{~}}{\mathit{Q}}}^{\text{red}}$ is obtained from $\stackrel{\mathbf{~}}{\mathit{Q}}(\mathit{k})$ by setting $k_{2N-1}=0$.

Note that ${\stackrel{\mathbf{~}}{\mathit{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 $\mathit{u}({\lambda }_i,\mathit{k})/{\sum }_{i=1}^N{u}_i({\lambda }_i,\mathit{k})$ is a steady state probability distribution of the reduced Markov chain with states $\{1,\dots ,N\}$ and dual generator ${\stackrel{\mathbf{~}}{\mathit{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,\mathit{k})=0$.

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 4

Consider that G is strongly connected. Then the equation $(\lambda \mathit{I}-\stackrel{\mathbf{~}}{\mathit{Q}}(\mathit{k}))\mathit{u}=0$ has solutions $\mathit{u}(\lambda ,\mathit{k})=(u_1(\lambda ,\mathit{k}),\cdots ,u_N(\lambda ,\mathit{k}))$ with $u_i(\lambda ,\mathit{k})\ne 0$ for all $1\le i\le N$.

Proof

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

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.

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 $\mathit{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{d{X}_{N+1}}{\mathit{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}_{\mathit{sj}}({\lambda }_i,\mathit{Q})}{{\prod }_{1\le l\ne i\le N}({\lambda }_i-{\lambda }_l)},$$where $D_{\mathit{sj}}$ are defined as in Sect. 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 Sect. 3.4, we have two cases:

i) $s\notin \text{Pred}(N+1)$.

In this case $S_1=0$, and the argument in Proposition 2 holds: the model is not solvable. The remaining $N-2$ symmetrized equations resulting from eigenvectors read

34

$$\begin{array}{cc}\hfill S_2=& \sum _{j\in \text{Pred}(N+1)}{(-1)}^{|s-j|+1}{Q}_{j,N+1}{d}_{N,N-2}(\mathit{Q}),\hfill \\ \hfill S_3=& \sum _{j\in \text{Pred}(N+1)}{(-1)}^{|s-j|+1}{Q}_{j,N+1}(d_{N,N-2}(\mathit{Q})h_1(\mathit{\lambda })+d_{N,N-3}(\mathit{Q})),\hfill \\ \hfill & ⋮\hfill & \hfill \\ \hfill S_k=& \sum _{j\in \text{Pred}(N+1)}{(-1)}^{|s-j|+1}{Q}_{j,N+1}(d_{N,N-2}(\mathit{Q})h_{k-2}(\mathit{\lambda })+d_{N,N-3}(\mathit{Q})h_{k-3}(\mathit{\lambda })+\hfill \\ \hfill +& \dots +d_{N,N-k}(\mathit{Q})),\hfill \\ \hfill & ⋮\hfill & \hfill \\ \hfill S_{N-1}=& \sum _{j\in \text{Pred}(N+1)}{(-1)}^{|s-j|+1}{Q}_{j,N+1}(d_{N,N-2}(\mathit{Q})h_{N-3}(\mathit{\lambda })+\dots +d_{N,2}(\mathit{Q})h_1(\mathit{\lambda })+\hfill \\ \hfill +& d_{N,1}(\mathit{Q})),\hfill \end{array}$$where $d_{N,k}(\mathit{Q}),0\le k\le N-2$ are defined as in Sect. 3.4.

ii) $s\in \text{Pred}(N+1)$.

In this case we have $N-1$ equations resulting from the eigenvectors:

35

$$\begin{array}{cc}\hfill S_1=& -Q_{s,N+1},\hfill \\ \hfill S_2=& -Q_{s,N+1}(h_1(\mathit{\lambda })+c_{N-2}(\mathit{Q}))+\hfill \\ \hfill & +\sum _{j\in \text{Pred}(N+1),j\ne s}{(-1)}^{|s-j|+1}{Q}_{j,N+1}{d}_{j,N-2}(\mathit{Q}),\hfill \\ \hfill & ⋮\hfill & \hfill \\ \hfill S_{N-1}=& -Q_{s,N+1}(h_{N-2}(\mathit{\lambda })+c_{N-2}(\mathit{Q})h_{N-3}(\mathit{\lambda })+\dots +c_2(\mathit{Q})h_1(\mathit{\lambda })+c_1(\mathit{Q}))+\hfill \\ \hfill & +\sum _{j\in \text{Pred}(N+1),j\ne s}{(-1)}^{|s-j|+1}{Q}_{j,N+1}(d_{j,N-2}(\mathit{Q})h_{N-3}(\mathit{\lambda })+d_{j,N-3}(\mathit{Q})h_{N-4}(\mathit{\lambda })+\hfill \\ \hfill +& \dots +d_{j,1}(\mathit{Q})),\hfill \end{array}$$where $c_k(\mathit{Q}),0\le k\le N-2$ are defined as in Sect. 3.4.

Solution of the Inverse Problem for the Unbranched Chain Model

[See PDF for image]

Fig. 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

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 Fig. 4. For this model$$\begin{array}{cc}& {\tilde{Q}}_{i+1,i}=k_i^{+},{\tilde{Q}}_{i,i+1}=k_i^{-},1\le i<N,{\tilde{Q}}_{i,i}\hfill \\ \hfill & =-k_i^{+}-k_{i-1}^{-},1\le i<N,{\tilde{Q}}_{N,N}=-k_N-k_{N-1}^{-}.\hfill \end{array}$$This model is a generalization, of arbitrary length N, of the model M2 represented in Fig. 2 that has $N=3$. To ensure that the model is solvable, as stated in Proposition 2, we consider the return state to be $s=N$.

The solution of the direct problem presented in the Sect. 2 provides an algebraic map1

36

$$\begin{array}{cc}& f:{\mathbb{C}}^{2N-1}\to {\mathbb{C}}^{2N-1},f(k_1^{+},\cdots ,k_{N-1}^{+},k_1^{-},\cdots ,k_{N-1}^{-},k_N)\hfill \\ \hfill & =(A_1,\cdots ,A_{N-1},{\lambda }_1,\cdots ,{\lambda }_N).\hfill \end{array}$$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,\cdots ,C_N)}^T$ represent a vector, where $C_i$ are solutions of (15).

Note that (15), (31) imply that

37

$$\begin{array}{cc}\hfill \mathit{U}C=& {(0,\dots ,0,1)}^T,\hfill \end{array}$$

38

$$\begin{array}{cc}\hfill k_N=& -\sum _{1\le j\le N}{A}_j{\lambda }_j,\hfill \end{array}$$where $\mathit{U}$ is the $N\times N$ matrix with the columns $\mathit{u}({\lambda }_1,\mathit{k}),\cdots ,\mathit{u}({\lambda }_N,\mathit{k})$ and ${\lambda }_1,\dots ,{\lambda }_N$ are the eigenvalues of $\stackrel{\mathbf{~}}{\mathit{Q}}(\mathit{k})$.

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

Lemma 2

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,\cdots ,A_{N-1},k_{N-i+1}^{+},\cdots ,k_{N-1}^{+}{k}_{N-i}^{-},\cdots ,k_{N-1}^{-},k_N,{\lambda }_1,\cdots ,{\lambda }_N),$$for $1\le i\le N-1$.

Proof

For each eigenvalue $\lambda$ of the matrix $\tilde{\mathbf{Q}}(\mathit{k})$ it holds$$\begin{array}{cc}& k_{N-1}^{+}{u}_{N-1}(\lambda ,\mathit{k})=k_N+k_{N-1}^{-}+\lambda ,k_n^{+}{u}_n(\lambda ,\mathit{k})\hfill \\ \hfill & =(k_{n+1}^{+}+k_n^{-}+\lambda )u_{n+1}(\lambda ,\mathit{k})-k_{n+1}^{-}{u}_{n+2}(\lambda ,\mathit{k}),\hfill \end{array}$$for $0\le n\le N-2$, considering that $u_0(\lambda ,\mathit{k})=0$.

This provides by recursion on $N-n$ a rational function $g_n$ with rational coefficients such that

39

$$\begin{array}{c}\hfill k_n^{+}{u}_n(\lambda ,\mathit{k})=g_n(k_{n+1}^{+},\cdots ,k_{N-1}^{+},k_n^{-},\cdots ,k_{N-1}^{-},k_N,\lambda ).\end{array}$$Moreover, the denominator of $g_n$ equals $k_{n+1}^{+}\cdots k_{N-1}^{+}$.

Then for $i\ge 0$ it holds

40

$$\begin{array}{c}\hfill {\tilde{\mathbf{Q}}}^{i}UC=U·\text{diag}({\lambda }_1^i,\cdots ,{\lambda }_N^i)C={\tilde{\mathbf{Q}}}^i{(0,\cdots ,0,1)}^T,\end{array}$$where $\text{diag}({\lambda }_1^i,\cdots ,{\lambda }_N^i)$ denotes a diagonal matrix.

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

41

$$\begin{array}{c}\hfill \sum _{1\le j\le N}{u}_{N-i}({\lambda }_j,\mathit{k}){\lambda }_j^i(-A_j{\lambda }_j/k_N).\end{array}$$The same $(N-i)$-th coordinate of the right vector in (40) equals

42

$$\begin{array}{c}\hfill k_{N-i}^{-}\cdots k_{N-1}^{-}.\end{array}$$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

43

$$\begin{array}{cc}\hfill k_{N-i}^{+}{k}_{N-i}^{-}\cdots k_{N-1}^{-}=& G_{N-i,0}(A_1,\cdots ,A_{N-1},k_{N-i+1}^{+},\cdots ,k_{N-1}^{+},k_{N-i}^{-},\hfill \\ \hfill & \cdots \hfill & \hfill ,k_{N-1}^{-},k_N,{\lambda }_1,\cdots ,{\lambda }_N)\end{array}$$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).

Lemma 3

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

44

$$\begin{array}{c}\hfill k_{N-i}^{-}=P_{N-i}(A_1,\cdots ,A_{N-1},k_{N-i+1}^{+},\cdots ,k_{N-1}^{+},k_{N-i+1}^{-},\cdots ,k_{N-1}^{-},k_N,{\lambda }_1,\cdots ,{\lambda }_N).\end{array}$$

Proof

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

45

$$\begin{array}{c}\hfill \sum _{1\le j\le N}{u}_{N-i+1}({\lambda }_j,\mathit{k}){\lambda }_j^i(-A_j{\lambda }_j/k_N).\end{array}$$The same $(N-i+1)$-th coordinate of the right vector in (40) equals

46

$$\begin{array}{c}\hfill -k_{N-i}^{-}\cdots k_{N-1}^{-}+p_{N-i}(k_{N-i+1}^{+},\cdots ,k_{N-1}^{+},k_{N-i+1}^{-},\cdots ,k_{N-1}^{-},k_N)\end{array}$$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}^{+},\cdots ,k_{N-1}^{+},k_{N-i+1}^{-},\cdots ,k_{N-1}^{-},k_N)=$$$$P_{N-i,0}(A_1,\cdots ,A_{N-1},k_{N-i+2}^{+},\cdots ,k_{N-1}^{+},k_{N-i+1}^{-},\cdots ,k_{N-1}^{-},k_N,{\lambda }_1,\cdots ,{\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 (45) and of (46) (while the both latter multiplied by $k_{N-i+1}^{+}$).

Applying alternatingly Lemma 3 and Lemma 2 for $i=1,\cdots ,N-1$ consecutively, we conclude with the following main result of this section.

Theorem 1

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,\cdots ,A_{N-1},{\lambda }_1,\cdots ,{\lambda }_N),k_n^{-}=F_n^{-}(A_1,\cdots ,A_{N-1},{\lambda }_1,\cdots ,{\lambda }_N).$$

Remark 3

1. Together with (38) Theorem 1 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}^{+}·k_1^{-}\cdots k_{N-1}^{-}·k_N\ne 0$ and ${\lambda }_i\ne {\lambda }_j,1\le i<j\le N$;

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

Remark 4

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,\cdots ,A_{N-1},{\lambda }_1,\cdots ,{\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,\cdots ,A_{N-1},{\lambda }_1,\cdots ,{\lambda }_{N-1})\subset \mathbb{C}(k_1^{+},\cdots ,k_{n-1}^{+},k_1^{-},\cdots ,k_{n-1}^{-},k_N)$$where $\mathbb{C}(A_1,\cdots ,A_{N-1},{\lambda }_1,\cdots ,{\lambda }_{N-1})$ denotes the field generated by $A_1,\cdots ,A_{N-1},{\lambda }_1,\cdots ,{\lambda }_{N-1}$ over the field of complex numbers. It is known (see e.g. (Shafarevich, 1972), 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,\cdots ,A_{N-1}={\alpha }_{N-1},{\lambda }_1={\beta }_1,\cdots ,{\lambda }_N={\beta }_N$ at a generic point $({\alpha }_1,\cdots ,{\alpha }_{N-1},{\beta }_1,\cdots ,{\beta }_N)\in {\mathbb{C}}^{2N-1}$. Recall that the degree is defined as the dimension of the vector space $\mathbb{C}(k_1^{+},\cdots ,k_{n-1}^{+},k_1^{-},\cdots ,k_{n-1}^{-},k_N)$ over the field $\mathbb{C}(A_1,\cdots ,A_{N-1},{\lambda }_1,\cdots ,{\lambda }_{N-1})$. The degree can be infinite. Theorem 1 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.

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.

The Thomas Decomposition of an Algebraic System

In this section we introduce briefly the notion of the algebraic Thomas decomposition (see the appendix of (Lange-Hegermann et al., 2021) 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}[\mathit{x}]$ be a polynomial ring in n variables $\mathit{x}=(x_1,\cdots ,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

47

$$\begin{array}{c}\hfill \mathcal{S}=\{p_1(\mathit{x})=0,\cdots ,p_r(\mathit{x})=0,q_1(\mathit{x})\ne 0,\cdots ,q_s(\mathit{x})\ne 0\}\end{array}$$with polynomials $p_i(\mathit{x})$, $q_j(\mathit{x})$ in $\mathbb{C}[\mathit{x}]$ and integers r, $s\in {\mathbb{N}}_0$. The solution set $\text{Sol}(\mathcal{S})$ of the algebraic system (47) is defined as the set of all $\overline{\mathit{x}}=({\overline{x}}_1,\cdots ,{\overline{x}}_n)\in {\mathbb{C}}^n$ satisfying the equations and inequations of $\mathcal{S}$, that is as$$\begin{array}{c}\hfill \text{Sol}(\mathcal{S})=\{\overline{\mathit{x}}\in {\mathbb{C}}^n\mid p_i(\overline{\mathit{x}})=0,q_j(\overline{\mathit{x}})\ne 0\text{for}\text{all}1\le i\le r\text{and}1\le j\le s\}.\end{array}$$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$$\begin{array}{c}\hfill \{\overline{\mathit{x}}\in {\mathbb{C}}^n\mid p_1(\overline{\mathit{x}})=0,\cdots ,p_r(\overline{\mathit{x}})=0\}\text{and}\{\overline{\mathit{x}}\in {\mathbb{C}}^n\mid q_1(\overline{\mathit{x}})\cdots q_s(\overline{\mathit{x}})=0\}\end{array}$$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 $\mathit{x}=(x_1,\cdots ,x_n)$ of our polynomial ring $\mathbb{C}[\mathit{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(\mathit{x}))$ of a non-constant polynomial $p(\mathit{x})$ is defined as the greatest variable appearing in $p(\mathit{x})$. In case $p(\mathit{x})\in \mathbb{C}$ is a constant polynomial, we set $\text{ld}(p(\mathit{x}))=1$. If we consider every polynomial $p(\mathit{x})\in \mathbb{C}[\mathit{x}]$ as a univariate polynomial in its leader, say $\text{ld}(p(\mathit{x}))=x_k$, then the coefficients of $p(\mathit{x})$ as a polynomial in $x_k$ are polynomials in $\mathbb{C}[x_1,\cdots ,x_{k-1}]$. The coefficient of the highest power of $\text{ld}(p(\mathit{x}))$ in $p(\mathit{x})$ is called the initial of $p(\mathit{x})$ which we denote by $\text{init}(p(\mathit{x}))$. The separant $\text{sep}(p(\mathit{x}))$ of a polynomial $p(\mathit{x})$ is defined as the partial derivative of $p(\mathit{x})$ with respect to its leader.

Definition 2

Let $\mathcal{S}$ be the algebraic system of (47). 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(\mathit{x})),\cdots ,\text{ld}(p_r(\mathit{x})),\text{ld}(q_1(\mathit{x})),\cdots ,\text{ld}(q_s(\mathit{x}))\}\backslash \{1\})=r+s.$$ This property is called triangularity.

2. For every $p(\mathit{x})\in \{p_1(\mathit{x}),\cdots ,p_r(\mathit{x}),q_1(\mathit{x}),\cdots ,q_s(\mathit{x})\}$ the equation $\text{init}(p(\mathit{x}))=0$ has no solution in $\text{Sol}(\mathcal{S})$. We call this property non-vanishing initials.

3. For every $p(\mathit{x})\in \{p_1(\mathit{x}),\cdots ,p_r(\mathit{x}),q_1(\mathit{x}),\cdots ,q_s(\mathit{x})\}$ the equation $\text{sep}(p(\mathit{x}))=0$ has no solution in $\text{Sol}(\mathcal{S})$. This is called square-freeness.

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

Definition 3

A Thomas decomposition of an algebraic system $\mathcal{S}$ as in (47) consists of finitely many simple algebraic systems ${\mathcal{S}}_1,\cdots ,{\mathcal{S}}_m$ such that $\text{Sol}(\mathcal{S})$ is the disjoint union of $\text{Sol}({\mathcal{S}}_1),\cdots ,\text{Sol}({\mathcal{S}}_m)$.

It was proved by Thomas in (Thomas, 1962; Thomas, 1937) that any algebraic system has a Thomas decomposition which is in general not unique. A Thomas decomposition can be determined algorithmically (see (Bächler et al., 2012)) and there is an implementation in Maple. A description of the implementation can be found in (Bächler and Lange-Hegermann, 2008-2012; Gerdt et al., 2019).

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$$\begin{array}{c}\hfill \mathbb{C}[k_1,k_2,k_3,k_4,k_5,L_1,L_2,L_3,S_1,S_2].\end{array}$$To simplify notation we collect the variables into $\mathit{k}=(k_1,k_2,k_3,k_4,k_5)$ and $\mathit{v}=(L_1,L_2,L_3,S_1,S_2)$. Since we want to solve the symmetrized systems for the variables $\mathit{k}$, we rank them always higher than the collection of variables $\mathit{v}$. Among the variables $\mathit{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 $\mathit{k}>\mathit{v}$ implies that each simple algebraic system returned by the Thomas decomposition yields solutions for the variables $\mathit{k}$ which are only valid for solutions of the variables $\mathit{v}$ satisfying equation and inequation conditions in $\mathit{v}$ over $\mathbb{C}$.

A Thomas Decomposition for Model 2

The model studied here is a particular case of the one from Sect. 4 and Theorem 1. We are going to determine the solutions of the algebraic system of equations$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \mathcal{S}& =& \{-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,\hfill \\ \hfill & & 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\}.\hfill \end{array}\end{array}$$To this end we compute with Maple an algebraic Thomas decomposition of $\mathcal{S}$ with respect to the ranking$$k_1>k_3>k_2>k_4>k_5>\mathit{v}.$$Maple returns eleven simple systems ${\mathcal{S}}_1,\cdots ,{\mathcal{S}}_{11}$. We will only compute here the solutions in $\mathit{k}$ for the first simple systems, since it the most generic one, meaning that the solutions are valid for specializations of the variables $\mathit{v}$ satisfying only inequations, i.e the solutions for $\mathit{k}$ are valid over an Zariski open subset of ${\mathbb{C}}^5$. The first simple system consists of the equations and inequations$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill {\mathcal{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\hfill \\ \hfill & & +(-S_1^4+S_1^2{S}_2)L_2+(S_1^3-S_1{S}_2)L_3)k_1+(-2S_1^2{S}_2^2+(-S_1^4-S_1^2{S}_2)L_2\hfill \\ \hfill & & +(S_1^3+S_1{S}_2)L_3)L_1+(S_1^3{S}_2-2L_3{S}_1^2+S_1{S}_2^2)L_2+(S_1^4-3S_1^2{S}_2)L_3=0,\hfill \\ \hfill & & (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,\hfill \\ \hfill & & -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,\hfill \\ \hfill & & k_5+S_1=0,L_1{S}_1{S}_2-L_2{S}_1^2+L_3{S}_1-S_2^2\ne 0,S_1^3-S_1{S}_2\ne 0,S_2\ne 0\}.\hfill \end{array}\end{array}$$One easily checks that the last three inequations do not involve any variable $\mathit{k}$. Thus the inequations

48

$$\begin{array}{c}\hfill L_1{S}_1{S}_2-L_2{S}_1^2+L_3{S}_1-S_2^2\ne 0,S_1^3-S_1{S}_2\ne 0,S_2\ne 0\end{array}$$define the above described Zariski open subset. Solving now successively the remaining equations for the variables $\mathit{k}$ (they are all linear in their respective leaders), we obtain the solutions$$\begin{array}{cc}\hfill 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_1^2{S}_2)L_2+(S_1^3-S_1{S}_2)L_3)}\hfill \\ \hfill & (L_1^2{S}_1^3{S}_2+L_2^2{S}_1^3+S_1{S}_2^3+L_3^2{S}_1+(-2S_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\hfill \\ \hfill & +(S_1^3{S}_2-2L_3{S}_1^2+S_1{S}_2^2)L_2+(S_1^4-3S_1^2{S}_2)L_3),\hfill \\ \hfill 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)},\hfill \\ \hfill 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},\hfill \\ \hfill k_4& =\frac{S_1^2-S_2}{S_1},\hfill \\ \hfill k_5& =-S_1,\hfill \end{array}$$which are only valid for those $\overline{\mathit{v}}\in {\mathbb{C}}^5$ satisfying the inequations in (48). If one is now interested in a solutions for $\mathit{k}$ with respect to a $\overline{\mathit{v}}\in {\mathbb{C}}^5$ which does not satisfy the inequations of (48), then there is a simple system among ${\mathcal{S}}_2,\cdots ,{\mathcal{S}}_{11}$, which one can solve successively for $\mathit{k}$ as described in the previous subsection. The remaining simple systems ${\mathcal{S}}_2,\cdots ,{\mathcal{S}}_{11}$ can be found in Sect. A.1 of the appendix.

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

A Thomas Decomposition for Model 3

The symmetrized algebraic system for model 3 is$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \mathcal{S}& =& \{-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,\hfill \\ \hfill & & 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\}\hfill \end{array}\end{array}$$and we compute an algebraic Thomas decomposition for it with respect to the ranking$$k_1>k_2>k_3>k_4>k_5>\mathit{v}.$$Maple returns 10 simple systems ${\mathcal{S}}_1,\cdots ,{\mathcal{S}}_{10}$. One easily checks by comparing the number of equations only involving the variables $\mathit{v}$ (see appendix A.2 for the remaining simple systems ${\mathcal{S}}_2,\cdots ,{\mathcal{S}}_{10}$), that the first simple system$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill {\mathcal{S}}_1& =& \{k_1{k}_3{S}_1{S}_2+k_3^2{S}_1{S}_2+L_3{S}_2\hfill \\ \hfill & & +(L_2{S}_1^2-L_3{S}_1)k_3=0,k_2{k}_3{S}_1-L_3=0,k_3\ne 0,\hfill \\ \hfill & & k_4{S}_1-S_1^2+S_2=0,k_5+S_1=0,L_1{S}_1{S}_2\hfill \\ \hfill & & -L_2{S}_1^2+L_3{S}_1-S_2^2=0,S_1\ne 0,S_2\ne 0\}\hfill \end{array}\end{array}$$is the most generic one. Here, even in the most generic case, the solutions for $\mathit{k}$ are only valid for $\overline{\mathit{v}}\in {\mathbb{C}}^5$ lying on a Zariski locally closed subset defined by

49

$$\begin{array}{c}\hfill L_1{S}_1{S}_2-L_2{S}_1^2+L_3{S}_1-S_2^2=0,S_1\ne 0,S_2\ne 0.\end{array}$$Solving again the remaining equations successively for $k_5$, $k_4$, $k_3$, $k_2$, $k_1$, we obtain the solutions$$\begin{array}{cc}\hfill k_5& =-S_1\hfill \\ \hfill k_4& =\frac{S_1^2-S_2}{S_1}\hfill \\ \hfill k_3& \text{arbitrary}\text{in}\mathbb{C}\backslash \{0\},\hfill \\ \hfill k_2& =\frac{L_3}{k_3{S}_1},\hfill \\ \hfill 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),\hfill \end{array}$$where the expression for $k_2$ and $k_1$ depend on the choice made for $k_3$.

The conditions (49) 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 Sect. 3.6 this means that the model M3 is not solvable.

A Thomas Decomposition for Model 4

In case of model 4 the symmetrized algebraic system is$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \mathcal{S}& =& \{-S_1=k_5,-S_2=k_5(k_1+k_2+k_3+L_1),L_1=-k_1-k_2-k_3-k_4-k_5,\hfill \\ \hfill & & 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,\hfill \\ \hfill & & L_3=-k_1{k}_3{k}_5-k_2{k}_3{k}_5\}.\hfill \end{array}\end{array}$$Using the ranking$$k_1>k_2>k_3>k_4>k_5>\mathit{v}$$the Maple implementation of algebraic Thomas decomposition returns 21 simple systems (see Sect. A.3 of the appendix for the systems ${\mathcal{S}}_2,\cdots ,{\mathcal{S}}_{21}$). The first simple system$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill {\mathcal{S}}_1& =& \{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,\hfill \\ \hfill & & L_2{S}_1^2-L_1{S}_1{S}_2-L_3{S}_1+S_2^2\hfill \\ \hfill & & +(S_1^3-S_1{S}_2)k_2=0,k_3^2{S}_1+(L_1{S}_1-S_2)k_3+L_3=0,\hfill \\ \hfill & & k_4{S}_1-S_1^2+S_2=0,k_5+S_1=0,L_1^2{S}_1^2\hfill \\ \hfill & & -2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2\ne 0,L_3\ne 0,\hfill \\ \hfill & & S_1^3-S_1{S}_2\ne 0,S_2\ne 0\}.\hfill \end{array}\end{array}$$is the most generic one. Indeed, all inequations appearing in ${\mathcal{S}}_1$ involve only the variables $\mathit{v}$, namely$$\begin{array}{c}\hfill L_1^2{S}_1^2-2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2\ne 0,L_3\ne 0,S_1^3-S_1{S}_2\ne 0,S_2\ne 0,\end{array}$$and so the solutions for $\mathit{k}$ are valid for all $\overline{\mathit{v}}\in {\mathbb{C}}^5$ of the Zariski open subset defined by these inequations. We can now successively solve the remaining equations in ${\mathcal{S}}_1$ for the variables $k_5$, $k_4$, $k_3$, $k_2$ and $k_1$, where the equations for $k_5$, $k_4$, $k_2$, $k_1$ are all linear in their respective leaders except for the equation with leader $k_3$,$$k_3^2{S}_1+(L_1{S}_1-S_2)k_3+L_3=0,$$which is quadratic. Since we want to consider only real roots for $k_3$, we require that the discriminant is positive, that is we change the inequation $L_1^2{S}_1^2-2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2\ne 0$ into $L_1^2{S}_1^2-2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2>0$. Solving successively the equations with leader $k_5$, $k_4$, $k_3$, $k_2$, $k_1$ using for $k_3$ the quadratic formula, we obtain the solutions$$\begin{array}{cc}\hfill k_5& =-S_1,\hfill \\ \hfill k_4& =\frac{S_1^2-S_2}{S_1},\hfill \\ \hfill k_3^{1,2}& =-\frac{L_1{S}_1-S_2\pm \sqrt{L_1^2{S}_1^2-2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2}}{2S_1},\hfill \\ \hfill 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},\hfill \\ \hfill 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}\hfill \end{array}$$subject to the conditions

50

$$\begin{array}{c}\hfill L_1^2{S}_1^2-2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2>0,L_3\ne 0,S_1^3-S_1{S}_2\ne 0,S_2\ne 0.\end{array}$$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 (50). According to the Sect. 3.6 this means that the model M4 is solvable.

A Thomas Decomposition for Model 8

In case of model 8 the symmetrized algebraic system is$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \mathcal{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,\hfill \\ \hfill & & 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\}.\hfill \end{array}\end{array}$$We compute the algebraic Thomas decomposition of $\mathcal{S}$ with respect to the ranking$$k_1>k_3>k_4>k_2>k_5>\mathit{v}$$and obtain from Maple twelve simple systems ${\mathcal{S}}_1,\cdots ,{\mathcal{S}}_{12}$. One easily checks by comparing the number of equations which have leader one of the variables of $\mathit{v}$, that the first simple system$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill {\mathcal{S}}_1& =& \{S_1{k}_1{k}_2-L_3=0,k_3{k}_2{S}_1^2+L_1{S}_1{S}_2\hfill \\ \hfill & & -L_2{S}_1^2+L_3{S}_1-S_2^2+(-S_1^3+S_1{S}_2)k_2=0,\hfill \\ \hfill & & S_1^2{k}_2{k}_4-L_1{S}_1{S}_2+L_2{S}_1^2-L_3{S}_1\hfill \\ \hfill & & +S_2^2=0,k_2^2{S}_1+(L_1{S}_1-S_2)k_2+L_3=0,\hfill \\ \hfill & & k_5+S_1=0,L_1^2{S}_1^2-2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2\ne 0,L_3\ne 0,S_1\ne 0\}\hfill \end{array}\end{array}$$is the most generic one. Analogously as in case of model 4 one determines the solutions$$\begin{array}{cc}\hfill k_5& =-S_1,\hfill \\ \hfill k_2^{1,2}& =-\frac{L_1{S}_1-S_2\pm \sqrt{L_1^2{S}_1^2-2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2}}{2S_1},\hfill \\ \hfill 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},\hfill \\ \hfill 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},\hfill \\ \hfill k_1& =\frac{L_3}{S_1{k}_2^{1,2}}\hfill \end{array}$$of ${\mathcal{S}}_1$ subject to the conditions

51

$$\begin{array}{c}\hfill L_1^2{S}_1^2-2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2>0,L_3\ne 0,S_1\ne 0.\end{array}$$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-2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2\ne 0$ for the discriminant into $L_1^2{S}_1^2-2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2>0$ to guarantee that the roots are real. The remaining simple systems are presented in Sect. 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 (51). According to the Sect. 3.6 this means that the model M8 is solvable.

A Thomas Decomposition for Model 9

We compute an algebraic Thomas decomposition of the algebraic system

52

$$\begin{array}{cc}\hfill \mathcal{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,\hfill \\ \hfill & L_3=-k_1{k}_2{k}_5,L_2=k_1{k}_2+k_1{k}_4+k_1{k}_5+k_2{k}_3+k_2{k}_5\}\hfill \end{array}$$for model 9 with respect to the ranking$$\begin{array}{c}\hfill k_3>k_1>k_4>k_2>k_5>\mathit{v}.\end{array}$$The Maple implementation returns 15 simple systems ${\mathcal{S}}_1,\cdots ,{\mathcal{S}}_{15}$, where the most general simple system is$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill {\mathcal{S}}_1& =& \{L_1{S}_1{S}_2-L_2{S}_1^2+L_3{S}_1-S_2^2+(2k_2{S}_1^2+L_1{S}_1^2\hfill \\ \hfill & & -S_1{S}_2)k_3+(-S_1^3+S_1{S}_2)k_2=0,\hfill \\ \hfill & & -L_1{S}_1^2+L_2{S}_1+S_1{S}_2+(2k_2{S}_1+L_1{S}_1-S_2)\hfill \\ \hfill & & k_4+(-S_1^2+S_2)k_2-L_3=0,\hfill \\ \hfill & & 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,\hfill \\ \hfill & & k_5+S_1=0,L_1^2{S}_1^2-2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2\ne 0,L_3\ne 0,S_1\ne 0\}.\hfill \end{array}\end{array}$$Similar as for model 4 and 8 one determines the solutions

53

$$\begin{array}{cc}\hfill k_5=& -S_1,\hfill \\ \hfill k_2^{1,2}=& -\frac{L_1{S}_1-S_2\pm \sqrt{L_1^2{S}_1^2-2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2}}{2S_1},\hfill \\ \hfill 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}{2k_2^{1,2}{S}_1+L_1{S}_1-S_2},\hfill \\ \hfill k_1=& -\frac{k_2^{1,2}{S}_1+L_1{S}_1-S_2}{S_1},\hfill \\ \hfill 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(2k_2^{1,2}{S}_1+L_1{S}_1-S_2)}\hfill \end{array}$$subject to the conditions

54

$$\begin{array}{c}\hfill L_1^2{S}_1^2-2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2>0,L_3\ne 0,S_1\ne 0.\end{array}$$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-2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2\ne 0$ into $L_1^2{S}_1^2-2L_1{S}_1{S}_2-4L_3{S}_1+S_2^2>0$. The remaining 14 simple algebraic systems can be found in Sect. 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 (54). According to the Sect. 3.6 this means that the model M9 is solvable.

Model Degeneracy Effects

Several examples were used elsewhere (O’Cinneide, 1989) 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.

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 Fig. 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 (49). 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 Fig. 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 Sect. 7 dedicated to transcriptional bursting.

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 ${\stackrel{\mathbf{~}}{\mathit{Q}}}^{\text{red}}$ that is obtained from $\stackrel{\mathbf{~}}{\mathit{Q}}(\mathit{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:

55

$$\begin{array}{cc}\hfill {\stackrel{\mathbf{~}}{\mathit{Q}}}^{\text{red}}\left(\begin{array}{c}{p}_1(\mathit{k})\\ p_2(\mathit{k})\\ ⋮\\ p_N(\mathit{k})\end{array}\right)=& 0,\hfill \\ \hfill p_1(\mathit{k})+p_2(\mathit{k})+\dots +p_N(\mathit{k})=& 1.\hfill \end{array}$$Let ${\tau }_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 ${\stackrel{\mathbf{~}}{\mathit{Q}}}^{\text{red}}$:

56

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\tau }_i=-1/{\tilde{Q}}_{\mathit{ii}}^{\text{red}}=1/\sum _{j\ne i}{\tilde{Q}}_{\mathit{ij}}^{\text{red}}.\end{array}\end{array}$$The result of these calculations for $N=3$ is given in the Table 1.

Table 1. Reciprocal lifetimes ${\tau }_1^{-1}$ and steady state occupation probability $P_3$ for the solvable models M2, M4, M8, M9

Model	${\tau }_i^{-1}$	${\tau }_2^{-1}$	${\tau }_3^{-1}$	$P_3$
M2	$k_1+k_2$	$k_3$	$k_4$	$\frac{k_2{k}_3}{k_1{k}_4+k_2{k}_3+k_3{k}_4}$
M4	$k_1+k_2$	$k_3$	$k_4$	$\frac{k_1{k}_3+k_2{k}_3}{k_1{k}_3+k_1{k}_4+k_2{k}_3+k_3{k}_4}$
M8	$k_1$	$k_2$	$(k_3+k_4)$	$\frac{k_1{k}_2}{k_1{k}_2+k_1{k}_3+k_1{k}_4+k_2{k}_3}$
M9	$k_1$	$k_2$	$(k_3+k_4)$	$\frac{k_1{k}_2}{k_1{k}_2+k_1{k}_4+k_2{k}_3}$

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

Proposition 5

All solvable models that generate the same phase-type distribution have

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

2. the same value of the lifetime ${\tau }_N$ of the state N,

3. 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:$$\begin{array}{c}\hfill S_1-\frac{S_2}{S_1}=-k_{2N-1}+a_{N-1}(\mathit{k})-c_{N-2}(\mathit{k}).\end{array}$$From the definitions of $a_{N-1}(\mathit{k})$ and $c_{N-2}(\mathit{k})$ and elementary properties of the determinant we find$$-a_{N-1}(\mathit{k})=-c_{N-2}(\mathit{k})+{\tilde{Q}}_{\mathit{NN}}^{\text{red}}-k_{2N-1}.$$It follows$$S_1-\frac{S_2}{S_1}=-{\tilde{Q}}_{\mathit{NN}}^{\text{red}}={\tau }_N^{-1},$$which implies (ii).

Let $\overline{w}$ be the mean value of the phase-type distribution. The mean number of events on an interval [0, T] is $T/\overline{w}$. The same number is equal to $T{p}_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$$\overline{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 Fig. 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 (52) 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 (53) as

57

$$\begin{array}{cc}\hfill k_5=& -S_1,\hfill \\ \hfill k_1=& -\frac{L_1{S}_1-S_2\pm \sqrt{{(L_1{S}_1-S_2)}^2-4L_3{S}_1}}{2S_1},\hfill \\ \hfill k_2=& -\frac{L_1{S}_1-S_2\mp \sqrt{{(L_1{S}_1-S_2)}^2-4L_3{S}_1}}{2S_1},\hfill \\ \hfill k_3=& \frac{1}{2S_1}\left[S_1^2,-,S_2,\pm ,\frac{L_1{S}_1^3-S_1^2{S}_2-2L_2{S}_1^2+L_1{S}_1{S}_2+2L_3{S}_1-S_2^2}{\sqrt{L_1{S}_1-S_2{)}^2-4L_3{S}_1}}\right],\hfill \\ \hfill k_4=& \frac{1}{2S_1}\left[S_1^2,-,S_2,\mp ,\frac{L_1{S}_1^3-S_1^2{S}_2-2L_2{S}_1^2+L_1{S}_1{S}_2+2L_3{S}_1-S_2^2}{\sqrt{L_1{S}_1-S_2{)}^2-4L_3{S}_1}}\right].\hfill \end{array}$$(57) 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 6

There are bijective rational maps $f_{M,M^{\prime }}:D\subset {\mathbb{R}}^5\to D^{\prime }\subset {\mathbb{R}}^5$, relating parameters $\mathit{k}$ of model M to parameters ${\mathit{k}}^{\mathbf{\prime }}$ of the model $M^{\prime }$ (i.e. ${\mathit{k}}^{\mathbf{\prime }}=f_{M,M^{\prime }}(\mathit{k})$) where $M,M^{\prime }$ are any pair in the list $\{M4,M8,M9\}$, and $D,D^{\prime }$ are open domains containing all the real positive parameters resulting from the inverse problem of $M,M^{\prime }$. $f_{M,M^{\prime }}$ satisfy the following constraints

1. The parameter $k_5$ is the same for all the models.

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

3. Under steady state conditions, the probability to be in the state 3 is the same for all the models.

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

Proof

Consider that there are mappings relating parameters of solvable models that have the same phase-type distribution, in particular the same parameters of the phase type distributions. Then, (i) and (ii) follow from the Proposition 5. The same Proposition implies that ${\tau }_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 ${\tau }_1$ and ${\tau }_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 ${\tau }_1$ and ${\tau }_2$). Reciprocally, if we impose (i), (ii), (iii) one finds four equations relating $\mathit{k}$ and ${\mathit{k}}^{\mathbf{\prime }}$. These equations, if they have solutions, define the map $f_{M,M^{\prime }}$. As a matter of fact, solutions exist and are easy to compute. For instance if $k_i$ are the parameters of M9 corresponding to the parameters $k_i^{\prime }$ of M8, then using the Table 1 the conditions (ii), 6 read

58

$$\begin{array}{cc}\hfill k_1=& k_1^{\prime },\hfill \\ \hfill k_2=& k_2^{\prime },\hfill \\ \hfill k_3+k_4=& k_3^{\prime }+k_4^{\prime },\hfill \\ \hfill \frac{k_1{k}_2}{k_1{k}_2+k_1{k}_4+k_2{k}_3}=& \frac{k_1^{\prime }{k}_2^{\prime }}{k_1^{\prime }{k}_2^{\prime }+k_1^{\prime }{k}_3^{\prime }+k_1^{\prime }{k}_4^{\prime }+k_2^{\prime }{k}_3^{\prime }}.\hfill \end{array}$$The solutions of system (58) read

59

$$\begin{array}{cc}\hfill k_1^{\prime }=& k_1,\hfill \\ \hfill k_2^{\prime }=& k_2,\hfill \\ \hfill k_3^{\prime }=& \frac{k_3(k_2-k_1)}{k_2},\hfill \\ \hfill k_4^{\prime }=& \frac{k_1{k}_3+k_2{k}_4}{k_2},\hfill \end{array}$$and define the bijective rational map $f_{M9,M8}:\{k_1>k_2>0,k_3>0,k_4>0\}\to \{k_2^{\prime }>k_1^{\prime }>0,k_3^{\prime }>0,k_4^{\prime }>\frac{k_1^{\prime }{k}_3^{\prime }}{k_2^{\prime }-k_1^{\prime }}\}$.

Similarly, solving

60

$$\begin{array}{cc}\hfill k_1=& k_1^{\prime }+k_2^{\prime },\hfill \\ \hfill k_2=& k_3^{\prime },\hfill \\ \hfill k_3+k_4=& k_4^{\prime },\hfill \\ \hfill \frac{k_1{k}_2}{k_1{k}_2+k_1{k}_4+k_2{k}_3}=& \frac{k_1^{\prime }{k}_3^{\prime }+k_2^{\prime }{k}_3^{\prime }}{k_1^{\prime }{k}_3^{\prime }+k_1^{\prime }{k}_4^{\prime }+k_2^{\prime }{k}_3^{\prime }+k_3^{\prime }{k}_4^{\prime }},\hfill \end{array}$$with respect to ${\mathit{k}}^{\prime }$ leads to the equations

61

$$\begin{array}{c}\hfill k_1^{\prime }=(k_1{k}_4-k_2{k}_4)/(k_3+k_4),\\ \hfill k_2^{\prime }=(k_1{k}_3+k_2{k}_4)/(k_3+k_4),\\ \hfill k_3^{\prime }=k_2,\\ \hfill k_4^{\prime }=k_3+k_4,\end{array}$$that define the map $f_{M9,M4}:\{k_1>k_2>0,k_3>0,k_4>0\}\to \{k_1^{\prime }>k_2^{\prime }>0,k_2^{\prime }>k_3^{\prime }>0,k_4^{\prime }>0\}$.

Remark 5

The Proposition 6 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}(\mathit{k})$ and $f_{M9,M8}\circ {\sigma }_{M9}(\mathit{k})$ where $\mathit{k}$ is a solution of M9.

Remark 6

The Proposition 6 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 7

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 Sect. 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.

[See PDF for image]

Fig. 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 allow observation of mRNA synthesis. The observed signal is deconvolved to obtain the transcription initiation events. C The phase-type distribution is used to infer the Markov chain model

Application to Transcriptional Bursting Models

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.

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

Possible state types
Model	state 1	state 2	state 3	state 4
M2	OFF	OFF	ON	EL
M3	OFF	OFF	ON	EL
M4	PAUSE /OFF	OFF	ON	EL
M8	PAUSE /OFF	OFF	ON	EL
M9	PAUSE /OFF	OFF/PAUSE	ON	EL
Possible biological interpretation of state types
long OFF	Enhancer-Promoter disconnection / Chromatin closing
short OFF	Transcription Factor unbinding / Pre-initiation Complex disassembly
PAUSE	Promoter proximal pausing / Intrinsic pausing / Transcriptional roadblock

The promoter starts transcription in the productive ON state, 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 (Pichon et al., 2018). Using the tool BurstDeconv (Douaihy et al., 2023) 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 intiation can be abortive. Abortive transcription initiation events cannot be detected in our experimental settings. Therefore, in this paper, whenever we refer to transcription initiation, we consider the processive case.

We also consider that the state EL (denoted $N+1$ in our abstract Markov chain model) can be reached only from a state that can be either ON or PAUSE, and after reaching EL the promoter instantly switches to 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 Sect. 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.

Discriminating Between Variant Models Using State Probabilities and Lifetimes

As discussed in Sect. 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 experimantal data such as single molecule DNA footprinting assays (Krebs, 2021; Kleinendorst et al., 2021).

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

Additionally, Proposition 5 shows that for any N, the occupancy probability $p_N$ and the lifetime ${\tau }_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,{\tau }_1,{\tau }_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$, ${\tau }_1$, ${\tau }_2$ and the corresponding differences $\mathrm{\Delta }{p}_1$, $\mathrm{\Delta }log({\tau }_1)$, $\mathrm{\Delta }log({\tau }_2)$ between the maximum and minimum values among the four variant models. Finally, we use these differences to estimate the capacity of such marker values to discriminate between variants.

We found that $p_1$ can discriminate between variants in $89\%$ of the cases, whereas ${\tau }_1,{\tau }_2$ can discriminate between variants in $84\%$ of the cases, see Fig. 6. We have used random phase-type distribution parameters spanning domain values observed in studies of transcriptional bursting of mammalian and fruit-fly genes (Tantale et al., 2021; Pimmett et al., 2021). This suggests that $p_1$ is better for discriminating among variant models than ${\tau }_1$ and ${\tau }_2$. This is confirmed also for the particular case of a HIV gene promoter in mammalian cells studied in (Tantale et al., 2021), where ${\tau }_1$ and ${\tau }_2$ are practically the same for all experimental conditions for the 3 variants, see Fig. 7 below.

[See PDF for image]

Fig. 6

Testing the discrimination capacity of $p_1$, ${\tau }_1$, and ${\tau }_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$, ${\tau }_1$, and ${\tau }_2$ is measured by the variations $\mathrm{\Delta }p$, $\mathrm{\Delta }{log}_{10}({\tau }_1)$, and $\mathrm{\Delta }{log}_{10}({\tau }_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$, ${\tau }_1$, and ${\tau }_2$ fail to discriminate in $11\%$, $16\%$, and $16\%$ of the cases, respectively

[See PDF for image]

Fig. 7

Values of $p_i$ and ${\tau }_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 ${\tau }_1$, ${\tau }_2$ of the OFF states (see (Tantale et al., 2021)). 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 ${\tau }_1$,${\tau }_2$ cannot

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\ne N$) must be combined with equations (55) and (56) 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 (Douaihy et al., 2023).

We applied our approach to the study of transcriptional bursting in two different biological context, human cultured cells (Tantale et al., 2021) and Drosophila developing embryos (Pimmett et al., 2021).

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 (Douaihy et al., 2023). 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 M2 and M9 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 (Tanenbaum et al., 2014; Basyuk et al., 2021; Blake et al., 2024), 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 (Liquet et al., 2012).

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.

Data Availability

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

1

although the rate parameters are real positive numbers, for formal reasons we define this map on complex numbers.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

Asanjarani, A; Liquet, B; Nazarathy, Y. Estimation of semi-Markov multi-state models: a comparison of the sojourn times and transition intensities approaches. Int J Biostat; 2021; 18, 1 pp. 243-262. [DOI: https://dx.doi.org/10.1515/ijb-2020-0083]

Bokes, P; Borri, A; Palumbo, P; Singh, A. Mixture distributions in a stochastic gene expression model with delayed feedback: a WKB approximation approach. J Math Biol; 2020; 81, 1 pp. 343-367.4123899 [DOI: https://dx.doi.org/10.1007/s00285-020-01512-y]

Bellec, M; Dufourt, J; Hunt, G; Lenden-Hasse, H; Trullo, A; ZineEl Aabidine, A; Lamarque, M; Gaskill, MM; Faure-Gautron, H; Mannervik, M; Harrison, MM; Andrau, J-C; Favard, C; Radulescu, O; Lagha, M. The control of transcriptional memory by stable mitotic bookmarking. Nat Commun; 2022; 13, 1 1176. [DOI: https://dx.doi.org/10.1038/s41467-022-28855-y]

Blake LA, De La Cruz A, Wu B: Imaging spatiotemporal translation regulation in vivo. In: Seminars in cell and developmental biology, vol. 154, pp. 155–164 (2024). Elsevier

Bächler, T; Gerdt, VP; Lange-Hegermann, M; Robertz, D. Algorithmic Thomas decomposition of algebraic and differential systems. J Symbolic Comput; 2012; 47, 10 pp. 1233-1266.2926124 [DOI: https://dx.doi.org/10.1016/j.jsc.2011.12.043]

Bladt, M. A review on phase-type distributions and their use in risk theory. ASTIN Bull J IAA; 2005; 35, 1 pp. 145-161.2142688 [DOI: https://dx.doi.org/10.2143/AST.35.1.583170]

Bächler T, Lange-Hegermann M (2008) AlgebraicThomas and DifferentialThomas: Thomas decomposition for algebraic and differential systems. (2008–2012). (https://www.art.rwth-aachen.de/cms/MATHB/Forschung/Mathematische-Software/~lqnwi/)

Bharucha-Reid, AT. Elements of the theory of markov processes and their applications; 1997; Mineola, NY, Dover:

Basyuk E, Rage F, Bertrand E (2021) RNA transport from transcription to localized translation: a single molecule perspective. RNA Biol 18(9):1221–1237

Commault C, Chemla JP (1996) An invariant of representations of phase-type distributions and some applications. J Appl Prob 33(2):368–381. https://doi.org/10.2307/3215060

Commault C, Mocanu S (2003) Phase-type distributions and representations: some results and open problems for system theory. Int J Control 76(6):566–580. https://doi.org/10.1080/0020717031000114986. Publisher: Taylor & Francis _eprint:10.1080/0020717031000114986

Damour, A; Slaninova, V; Radulescu, O; Bertrand, E; Basyuk, E. lTranscriptional stochasticity as a key aspect of HIV-1 latency. Viruses; 2023; 15, 9 1969. [DOI: https://dx.doi.org/10.3390/v15091969]

Dufourt, J; Trullo, A; Hunter, J; Fernandez, C; Lazaro, J; Dejean, M; Morales, L; Nait-Amer, S; Schulz, KN; Harrison, MM; Favard, C; Radulescu, O; Lagha, M. Temporal control of gene expression by the pioneer factor Zelda through transient interactions in hubs. Nat Commun; 2018; 9, 1 5194. [DOI: https://dx.doi.org/10.1038/s41467-018-07613-z]

Douaihy, M; Topno, R; Lagha, M; Bertrand, E; Radulescu, O. BurstDECONV: a signal deconvolution method to uncover mechanisms of transcriptional bursting in live cells. Nucleic Acids Res; 2023; 51, 16 pp. 88-88. [DOI: https://dx.doi.org/10.1093/nar/gkad629]

Dufresne, D. Fitting combinations of exponentials to probability distributions. Appl Stochast Mod Bus Ind; 2007; 23, 1 pp. 23-48.2344603 [DOI: https://dx.doi.org/10.1002/asmb.635]

Fackrell, M. Modelling healthcare systems with phase-type distributions. Health Care Manage Sci; 2009; 12, pp. 11-26. [DOI: https://dx.doi.org/10.1007/s10729-008-9070-y]

Faddy, M. (1993) A structured compartmental model for drug kinetics. Biometrics, pp. 243–248

Feller, W. An introduction to probability theory and its applications; 1966; Hoboken, NJ, John Wiley & Sons: 2

Fulton, W; Harris, J. Representation theory a first course; 1991; New York, Springer: 129

Gerdt, VP; Lange-Hegermann, M; Robertz, D. The maple package TDDS for computing thomas decompositions of systems of nonlinear PDEs. Comp Phys Commun; 2019; 234, pp. 202-215.3864422 [DOI: https://dx.doi.org/10.1016/j.cpc.2018.07.025]

Hössjer O, Bechly G, Gauger A (2018) Phase-type distribution approximations of the waiting time until coordinated mutations get fixed in a population. In: Silvestrov S, Malyarenko A, Rančić M (eds) Stochastic processes and applications, pp 245–313. Springer, New York. https://doi.org/10.1007/978-3-030-02825-1_12

Hobolth, A; Siri-Jégousse, A; Bladt, M. Phase-type distributions in population genetics. Theor Popul Biol; 2019; 127, pp. 16-32. [DOI: https://dx.doi.org/10.1016/j.tpb.2019.02.001]

Kleinendorst, RW; Barzaghi, G; Smith, ML; Zaugg, JB; Krebs, AR. Genome-wide quantification of transcription factor binding at single-DNA-molecule resolution using methyl-transferase footprinting. Nat Protocols; 2021; 16, 12 pp. 5673-5706. [DOI: https://dx.doi.org/10.1038/s41596-021-00630-1]

Kumar N, Kulkarni RV (2019) Constraining the complexity of promoter dynamics using fluctuations in gene expression. Phys Biol 17(1):015001. https://doi.org/10.1088/1478-3975/ab4e57.Publisher:IOPPublishing

Krebs, AR. Studying transcription factor function in the genome at molecular resolution. Trends Genet; 2021; 37, 9 pp. 798-806. [DOI: https://dx.doi.org/10.1016/j.tig.2021.03.008]

Kumar N, Singh A, Kulkarni RV (2015) Transcriptional bursting in gene expression: analytical results for general stochastic models. PLOS Comput Biol 11(10):1004292. https://doi.org/10.1371/journal.pcbi.1004292.Publisher:PublicLibraryofScience

Lange-Hegermann, M; Robertz, D; Seiler, WM; Seiß, M. Singularities of algebraic differential equations. Adv Appl Math; 2021; 131, 4305883 [DOI: https://dx.doi.org/10.1016/j.aam.2021.102266] 102266.

Liu Y (2022) Non-parametric bayesian inference with application to system biology. PhD thesis, Department of Statistics, University of Chicago. https://doi.org/10.6082/uchicago.3931

LaMar, DM; Kemper, P; Smith, GD. Reduction of calcium release site models via moment fitting of phase-type distributions. Phys Biol; 2011; 8, 2 [DOI: https://dx.doi.org/10.1088/1478-3975/8/2/026015] 026015.

Liquet B, Timsit J-F, Rondeau V (2012) Investigating hospital heterogeneity with a multi-state frailty model: application to nosocomial pneumonia disease in intensive care units. BMC Med Res Methodol 12(1):1–14

Maier RS (1991) The algebraic construction of phase-type distributions. Commun Stat Stochast Mod 7(4):573–602. https://doi.org/10.1080/15326349108807207. Publisher: Taylor & Francis _eprint:10.1080/15326349108807207

Moffitt JR, Bustamante C (2014) Extracting signal from noise: kinetic mechanisms from a Michaelis-Menten-like expression for enzymatic fluctuations. FEBS J 281(2):498–517. https://doi.org/10.1111/febs.12545. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/febs.12545

Meyn, SP; Tweedie, RL. Markov chains and stochastic stability; 2012; New York, Springer:

Neuts MF (1975) Probability distributions of phase type. Liber Amicorum Prof. Emeritus H, Florin

O’Cinneide CA (1989) On non-uniqueness of representations of phase-type distributions. Commun Stat Stochast Mod 5(2):247–259. https://doi.org/10.1080/15326348908807108. Publisher: Taylor & Francis _eprint:10.1080/15326348908807108

Pimmett, VL; Dejean, M; Fernandez, C; Trullo, A; Bertrand, E; Radulescu, O; Lagha, M. Quantitative imaging of transcription in living Drosophila embryos reveals the impact of core promoter motifs on promoter state dynamics. Nat Commun; 2021; 12, 1 4504. [DOI: https://dx.doi.org/10.1038/s41467-021-24461-6]

Pichon, X; Lagha, M; Mueller, F; Bertrand, E. A growing toolbox to image gene expression in single cells: sensitive approaches for demanding challenges. Mol Cell; 2018; 71, 3 pp. 468-480. [DOI: https://dx.doi.org/10.1016/j.molcel.2018.07.022]

Shafarevich, I. Basic algebraic geometry 1; 1972; Heidelberg, Springer:

Soltani M, Vargas-Garcia CA, Antunes D, Singh A (2016) Intercellular variability in protein levels from stochastic expression and noisy cell cycle processes. PLOS Comput Biol 12(8):1004972. https://doi.org/10.1371/journal.pcbi.1004972. Publisher: Public Library of Science

Stone, K; Zwiggelaar, R; Jones, P; Mac Parthaláin, N. A systematic review of the prediction of hospital length of stay: towards a unified framework. PLOS Digital Health; 2022; 1, 4 0000017. [DOI: https://dx.doi.org/10.1371/journal.pdig.0000017]

Tantale, K; Garcia-Oliver, E; Robert, M-C; L’Hostis, A; Yang, Y; Tsanov, N; Topno, R; Gostan, T; Kozulic-Pirher, A; Basu-Shrivastava, M; Mukherjee, K; Slaninova, V; Andrau, J-C; Mueller, F; Basyuk, E; Radulescu, O; Bertrand, E. Stochastic pausing at latent HIV-1 promoters generates transcriptional bursting. Nat Commun; 2021; 12, 1 4503. [DOI: https://dx.doi.org/10.1038/s41467-021-24462-5]

Tanenbaum, ME; Gilbert, LA; Qi, LS; Weissman, JS; Vale, RD. A protein-tagging system for signal amplification in gene expression and fluorescence imaging. Cell; 2014; 159, 3 pp. 635-646. [DOI: https://dx.doi.org/10.1016/j.cell.2014.09.039]

Thomas, JM. Differential systems; 1937; New York, Colloquium Publications XXI. American Mathematical Society: [DOI: https://dx.doi.org/10.1090/coll/021]

Thomas, JM. Systems and roots; 1962; Richmond, VA, W. Byrd Press:

Zhang J, Chen A, Qiu H, Zhang J, Tian T, Zhou T (2024) Exact results for gene-expression models with general waiting-time distributions. Phys Rev E 109(2):024119. https://doi.org/10.1103/PhysRevE.109.024119.Publisher:AmericanPhysicalSociety