Integrability of Stochastic Birth-Death processes via Differential Galois Theory

TL;DR

This paper investigates the mathematical integrability of stochastic birth-death processes using differential Galois theory, revealing that only linear rate cases are integrable, which impacts the solvability of these models.

Contribution

It applies differential Galois theory to analyze the integrability of birth-death processes with polynomial rates, providing new insights into their solvability.

Findings

PDE is non-integrable for polynomial rates

Linear rate cases are integrable

Differential Galois theory effectively analyzes stochastic process integrability

Abstract

Stochastic birth-death processes are described as continuous-time Markov processes in models of population dynamics. A system of infinite, coupled ordinary differential equations (the so-called master equation) describes the time-dependence of the probability of each system state. Using a generating function, the master equation can be transformed into a partial differential equation. In this contribution we analyze the integrability of two types of stochastic birth-death processes (with polynomial birth and death rates) using standard differential Galois theory. We discuss the integrability of the PDE via a Laplace transform acting over the temporal variable. We show that the PDE is not integrable except for the (trivial) case in which rates are linear functions of the number of individuals.

Tables

$\epsilon \left(\infty \right)$	$\epsilon \left(\delta /\beta \right)$	$\epsilon \left(1\right)$	$\epsilon \left(0\right)$	$m={\alpha }_{\infty }^{\epsilon \left(\infty \right)}-{\sum }_{c\in {\Gamma }^{\prime }}{\alpha }_c^{\epsilon \left(c\right)}$
$+$	$\pm$	$\pm$	$+$	$-2\mp \hat{s}$
$+$	$\pm$	$\pm$	$-$	$N_0-1\mp \hat{s}$
$-$	$\pm$	$\pm$	$+$	$-N_0-1\mp \hat{s}$
$-$	$\pm$	$\pm$	$-$	$\mp \hat{s}$
$+$	$\pm$	$\mp$	$+$	$-2$
$+$	$\pm$	$\mp$	$-$	$N_0-1$
$-$	$\pm$	$\mp$	$+$	$-N_0-1$
$-$	$\pm$	$\mp$	$-$	$0$

$\epsilon \left(\infty \right)$	$\epsilon \left(0\right)$	$\epsilon \left(1\right)$	$m={\alpha }_{\infty }^{\epsilon \left(\infty \right)}-{\alpha }_0^{\epsilon \left(0\right)}-{\alpha }_1^{\epsilon \left(1\right)}$
$+$	$\pm$	$\pm$	$-\frac{1}{2}-\frac{1}{2}-1=-2$
$-$	$\pm$	$\pm$	$\frac{1}{2}-\frac{1}{2}-1=-1$
$+$	$\mp$	$\mp$	$-\frac{1}{2}-\frac{1}{2}-1=-2$
$-$	$\mp$	$\mp$	$\frac{1}{2}-\frac{1}{2}-1=-1$

Equations

$$\text{Pr}\{N+1,t+\Delta t|N,t\}=B_N\Delta t,$$

$$\text{Pr}\{N-1,t+\Delta t|N,t\}=D_N\Delta t,$$

$$\text{Pr}\{N,t+\Delta t|N-1,t\}=B_{N-1}\Delta t,$$

$$\text{Pr}\{N,t+\Delta t|N+1,t\}=D_{N+1}\Delta t,$$

$$\text{Pr}\{N,t+\Delta t|N,t\}=1-(B_N+D_N)\Delta t.$$

$$P_N(t+\Delta t)=P_{N-1}(t)B_{N-1}\Delta t+P_{N+1}(t)D_{N+1}\Delta t+P_N(t)[1-(B_N+D_N)\Delta t].$$

$$P_N^{\prime }(t)=B_{N-1}{P}_{N-1}(t)+D_{N+1}{P}_{N+1}(t)-(B_N+D_N)P_N(t).$$

$$P_0^{\prime }(t)=D_1{P}_1(t)-B_0{P}_0(t).$$

$$g(z,t)=\sum _{N=0}^{\infty }{P}_N(t)z^N,$$

$$\frac{\partial g(z,t)}{\partial t}=(1-z)(\delta -\beta z)\frac{\partial g(z,t)}{\partial z}$$

$$g(z,0)=z^{N_0},g(1,t)=1.$$

$$\frac{\partial g(z,t)}{\partial t}=(1-z)\left[(\delta -\beta z)\frac{\partial g(z,t)}{\partial z}+\delta z\frac{{\partial }^{2}g(z,t)}{\partial z^2}\right].$$

$$g(z,t)={\left[\frac{\delta -\beta z-(1-z)\delta e^{(\beta -\delta )t}}{\delta -\beta z-(1-z)\beta e^{(\beta -\delta )t}}\right]}^{N_0}.$$

$$g(z,t)={\left[\frac{\delta (1-z)t+z}{\delta (1-z)t+1}\right]}^{N_0}.$$

$$\mathcal{C}=\{c\in K|c^{\prime }=0\},$$

$$\mathcal{L}:=y^{\prime \prime }+a{y}^{\prime }+by=0,a,b\in K.$$

$$A_{\sigma }=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\in \mathrm{GL}(2,\mathbb{C}),$$

$$\sigma \left(\begin{array}{c}{y}_1\\ y_2\end{array}\right)=\left(\begin{array}{c}\sigma (y_1)\\ \sigma (y_2)\end{array}\right)=\left(\begin{array}{cc}{y}_1& y_2\end{array}\right)A_{\sigma }.$$

$$\sigma \left(\begin{array}{cc}{y}_1& y_2\\ y_1^{\prime }& y_2^{\prime }\end{array}\right)=\left(\begin{array}{cc}\sigma (y_1)& \sigma (y_2)\\ \sigma (y_1^{\prime })& \sigma (y_2^{\prime })\end{array}\right)=\left(\begin{array}{cc}{y}_1& y_2\\ y_1^{\prime }& y_2^{\prime }\end{array}\right)A_{\sigma }.$$

$$\mathcal{L}:={\zeta }^{\prime \prime }=r\zeta ,r\in K.$$

$$y=e^{-\frac{1}{2}\int a}\zeta ,r=\frac{a^2}{4}+\frac{a^{\prime }}{2}-b$$

$${\partial }_{x}v=r-v^2,v=\frac{{\zeta }^{\prime }}{\zeta },$$

$$\frac{\partial g(z,t)}{\partial t}=\sum _{N=0}^{\infty }\left\{\beta (N-1)P_{N-1}(t)z^N+\delta (N+1)P_{N+1}(t)z^N-(\beta +\delta )N{P}_N(t)z^N\right\}.$$

$$\frac{\partial g(z,t)}{\partial t}=\left\{\beta z^2-(\beta +\delta )z+\delta \right\}\frac{\partial g(z,t)}{\partial z}.$$

$$\frac{\partial g}{\partial t}=Mg$$

$$MG=sG+g_0(z),$$

$$MG=sG,$$

$$G(z,s)={\int }_0^{\infty }g(z,t)e^{-st}dt,$$

$$(1-z)(\delta -\beta z)G^{\prime }(z,s)=sG(z,s)-z^{N_0},$$

$$G^{\prime }(z,s)=f(z,s)G(z,s)+h(z)$$

$$f(z,s)=\frac{s}{(1-z)(\delta -\beta z)},$$

$$h(z)=-\frac{z^{N_0}}{(1-z)(\delta -\beta z)}.$$

Full text

\FirstPageHeading

\ShortArticleName

Integrability of Stochastic Birth-Death processes via Differential Galois Theory

\ArticleName

Integrability of Stochastic Birth-Death processes via Differential Galois Theory

\Author

Primitivo B. ACOSTA-HUMÁNEZ {}^{{\dagger}^{1}}$${}^{{\dagger}^{2}}, José A. CAPITÁN {}^{{\dagger}^{3}}$${}^{{\dagger}^{4}} and Juan J. MORALES-RUIZ ${}^{{\dagger}^{3}}$

\AuthorNameForHeading

P.B. Acosta-Humánez, J.A. Capitán and J.J. Morales-Ruiz

\Address

${}^{{\dagger}^{1}}$ Instituto Superior de Formación Docente Salomé Ureña - ISFODOSU,

Recinto Emilio Prud’Homme, Calle R.C. Tolentino # 51, esquina 16 de Agosto,

Los Pepines, Santiago de los Caballeros, República Dominicana \EmailDprimitivo.acosta-humanez@isfodosu.edu.do \Address${}^{{\dagger}^{2}}$ Facultad de Ciencias Básicas y Biomédicas, Universidad Simón Bolívar,

Sede 3, Carrera 59 No. 58-135, Barranquilla, Colombia

\Address

${}^{{\dagger}^{3}}$ Depto. de Matemática Aplicada. E.T.S. Edificación. Avda. Juan de Herrera 6.

Universidad Politécnica de Madrid. 28040, Madrid, Spain \EmailDja.capitan@upm.es \EmailDjuan.morales-ruiz@upm.es \Address${}^{{\dagger}^{4}}$ Grupo de Sistemas Complejos. Universidad Politécnica de Madrid. 28040, Madrid, Spain

\Abstract

Stochastic birth-death processes are described as continuous-time Markov processes in models of population dynamics. A system of infinite, coupled ordinary differential equations (the so-called master equation) describes the time-dependence of the probability of each system state. Using a generating function, the master equation can be transformed into a partial differential equation. In this contribution we analyze the integrability of two types of stochastic birth-death processes (with polynomial birth and death rates) using standard differential Galois theory. We discuss the integrability of the PDE via a Laplace transform acting over the temporal variable. We show that the PDE is not integrable except for the (trivial) case in which rates are linear functions of the number of individuals.

\Keywords

differential Galois theory; stochastic processes; population dynamics; Laplace transform.

\Classification

12H05; 35A22; 35C05; 92D25.

1 Introduction

Stochastic birth-death processes [27, 10, 16] are widely used in the mathematical modeling of interacting populations. They are a special case of continuous-time Markov processes [15] for which transitions between states are either births (which increase the state variable in one unity) and deaths (which decrease the state variable by one). Birth-death processes have been used in different fields of applied science, with many applications in ecology [21, 13, 5], queueing theory [25], epidemiology [6] and population genetics [22], to mention just a few. In contrast to deterministic models, these kind of processes make the assumption that population changes take place in discrete numbers, and this fact introduces variability and noise when compared to deterministic dynamics [7, 8]. In the limit of infinite system size, these models are the counterpart of deterministic dynamics that usually appear in demography and population dynamics [19, 26].

Only few instances of birth-death processes are analytically tractable in mathematical terms. Most of the results are related to the probability distributions observed at stationarity [15, 12]. Little is known, however, about how these probability distributions change over time before reaching the equilibrium state. The existence of closed-form, analytical solutions for certain families of birth-death processes, even when certain restrictions on the parameters are forced to ensure the existence of analytical solutions, would be a powerful way to get valuable insights on how probability distributions behave over time and reach the steady state in such processes. In this contribution we focus on the existence of closed-form analytical solutions for two widely-used stochastic birth-death models with non-constant birth and death rates.

In order to properly contextualize the problem, we start by describing the mathematical framework used in the theory of stochastic birth-death processes. The central quantity used to characterize quantitatively a population formed by a number of individuals is precisely $N$, the number of individuals observed, where $N\in\mathbb{N}\cup\{0\}$. Let $P_{N}(t)$ be the probability that there are $N$ individuals at time $t$ —the latter variable is regarded as a continuous time— in the system. Given a particular way of how new individuals enter the system (births events) or leave the population (death events), the goal of the theory is to describe mathematically this probability. The stochastic process is fully described once the probability rates $B_{N}$ (births) and $D_{N}$ (deaths) are defined. For simplicity we assume here that rates are time-independent.

As such, birth and death rates, $B_{N}$ and $D_{N}$, are regarded as probabilities per unit time that a birth occurs (hence the populations moves from $N$ to $N+1$ individuals) or, correspondingly, that a death event occurs (and then the system changes from having $N$ to $N-1$ individuals). Consider an infinitesimal time interval $\Delta t$. Then birth and death rates satisfy

[TABLE]

where $\text{Pr}\{N+1,t+\Delta t\,|\,N,t\}$ is the conditional probability that the system undergoes a birth event at time $t+\Delta t$ given that there were $N$ individuals at time $t$. Multiple births and deaths are usually ignored in the limit $\Delta t\to 0$ because their probability would be proportional to $(\Delta t)^{2}$.

If the population is formed by $N$ individuals at time $t$, at time $t+\Delta t$ the population can be composed by: (a) $N+1$ individuals with probability $B_{N}\Delta t$; (b) $N-1$ individuals with probability $D_{N}\Delta t$; (c) $N$ individuals with probability $1-\{B_{N}+D_{N}\}\Delta t$. Therefore, the conditional probabilities that end up with a population formed exactly by $N$ individuals are:

[TABLE]

Thus, using the theorem of total probability we can write an expression for the probability of observing $N$ individuals at time $t+\Delta t$ in terms of the probabilities at time $t$:

[TABLE]

Here we are assuming a Markovian hypothesis, according to which the state of the system at a given time is conditioned only by the potential states at previous times but infinitely close to the current time. Now we subtract $P_{N}(t)$ from both sides of (1.3) and take the limit $\Delta t\to 0$ to get the so-called master equation:

[TABLE]

The system is therefore fully described by a coupled system of infinitely many ordinary differential equations, given by equation (1.4) for $N\geq 1$. For $N=0$, since the number of individuals has to remain non-negative, we have to impose that $D_{0}=0$ and $B_{-1}=0$. In this case, the corresponding equation reduces to

[TABLE]

Therefore, the central problem in the theory of stochastic birth-death processes for a population of individuals is to solve the master equation (1.4)–(1.5) as a problem of initial value. To be precise, if we know the probability distribution $P_{N}(t_{0})$ at some initial time $t_{0}$, then the system of differential equations allows for obtaining the probability distribution at any time $t$, $P_{N}(t)$. In what follows we will assume, without loss of generality, that there are $N_{0}$ individuals at time $t=0$, $N_{0}\in\mathbb{N}\cup\{0\}$. Then the initial probability distribution is precisely equal to $P_{N}(0)=\delta_{N_{0}N}$, where $\delta_{ij}$ stands for the usual Kronecker delta symbol.

To make the master equation tractable, in some cases it can be transformed into a partial differential equation by using a generating function defined as

[TABLE]

i.e., the discrete variable $N$ is transformed into to continuous variable $z$, $0\leq z\leq 1$. In this contribution we are interested in the conditions under which we can find a closed-form, analytical solution for the generating function $g(z,t)$. Knowledge of the generating function allows the calculation of important properties of the stochastic processes —for example, the average number of individuals, the variance of the population, or even the probability of extinction of the system at time $t$, $P_{0}(t)=g(0,t)$. We will follow a Laplace transform strategy to solve the corresponding partial differential equation, and we will analyze ecologically meaningful examples for the birth and death rates, which yield useful insights about the integrability of these kind of systems by considering one fully integrable case and a fully non-integrable case (the sense in which we use the term ‘integrable’ will be precisely defined in Section 2).

To be more specific, from now on we focus on the birth-death process defined by the rates (as mentioned, the term ‘rate’ stands for probability per unit time) $B_{N}=\beta N^{b}$ and $D_{N}=\delta N^{d}$, where $b$, $d$ are natural exponents and $\beta$, $\delta$ are positive real numbers. Usually death rates are taken as a quadratic function ($d=2$) since it is commonly assumed that two individuals compete with each other in death events, whereas birth processes (asexual reproduction) are described as linear functions of $N$ ($b=1$, i.e., the probability of a birth event is proportional to the number of individuals in the population). In this contribution we will consider two combinations of exponents: $(b,d)\in\{(1,1),(1,2)\}$. For the combination $(b,d)=(1,1)$, the master equation is equivalent to the following PDE (see Section 3):

[TABLE]

with boundary conditions

[TABLE]

This equation turns out to be integrable (in a sense specified below) via a Laplace transform technique. Roughly speaking, integrability here means solvability in closed form.

If death rates are quadratic functions of $N$, in Section 4 we show that the generating function satisfies the following PDE,

[TABLE]

We will refer to this PDE as the $(b,d)=(1,2)$ case. As before, the generating function has to satisfy the conditions (1.8). This case of quadratic death rates, which is the more relevant one in biological terms, remains as non-integrable, as we will show in Section 4 using results from Differential Galois Theory.

Proposition 1.1**.**

The PDE given by equation (1.7) with boundary conditions (1.8) is integrable and the solution is given by:

For $\beta\neq\delta$,

[TABLE] 2. 2.

For $\beta=\delta$,

[TABLE]

Proposition 1.2**.**

The PDE given by equation (1.9) is non-integrable.

Proposition 1.2 is a non-integrability result and tell us that any search for a closed-form, analytical solution for equation (1.9) is doomed to failure.

2 Differential Galois Theory

Differential Galois Theory, also known as Picard-Vessiot theory, is the Galois theory of linear differential equations. In classical Galois theory, the main object is a group of permutations of the polynomial’s roots, whereas in the Picard-Vessiot theory it is a linear algebraic group. For polynomial equations we look for solutions in terms of radicals. According to classical Galois theory, this form of the solution will exist whenever the Galois group is a solvable group. An analogous situation holds for linear homogeneous differential equations.

As a notational convention we will use $\partial_{x}:=\frac{\partial}{\partial x}$ (also ${}^{\prime}:=\frac{\partial}{\partial x})$ throughout this section.

2.1 Definitions and Known Results

The following theoretical background can be found in the references [9, 20, 23]. We recall that although differential Galois theory is more general, here we just summarize results from theory for second order differential equations.

Definition 2.1** (Differential Fields).**

Let $K$ (depending on $x$) be a commutative field of characteristic zero, and $\partial_{x}$ a derivation, that is, a map $\partial_{x}:K\rightarrow K$ satisfying $\partial_{x}(a+b)=\partial_{x}a+\partial_{x}b$ and $\partial_{x}(ab)=\partial_{x}a\cdot b+a\cdot\partial_{x}b$ for all $a,b\in K$. By $\mathcal{C}$ we denote the field of constants of $K$,

[TABLE]

which is also of characteristic zero and will be assumed algebraically closed. In this terms, we say that $K$ is a differential field with the derivation $\partial_{x}=\,\,^{\prime}$.

Up to special considerations, we analyze second order linear homogeneous differential equations, that is, equations in the form

[TABLE]

Definition 2.2** (Picard-Vessiot Extension).**

Suppose that $y_{1},y_{2}$ is a basis of solutions of $\mathcal{L}$ given in equation (2.1), i.e., $y_{1},y_{2}$ are linearly independent over $K$ and every solution is a linear combination over $\mathcal{C}$ of these two. Let $L=K\langle y_{1},y_{2}\rangle=K(y_{1},y_{2},y_{1}^{\prime},y_{2}^{\prime})$ the differential extension of $K$ such that $\mathcal{C}$ is the field of constants for $K$ and $L$. In this terms, we say that $L$, the smallest differential field containing $K$ and $\{y_{1},y_{2}\}$, is the Picard-Vessiot extension of $K$ for $\mathcal{L}$.

Definition 2.3** (Differential Galois Groups).**

Assume $K$, $L$ and $\mathcal{L}$ as in the previous definition. The group of all differential automorphisms (automorphisms that commute with derivation) of $L$ over $K$ is called the differential Galois group of $L$ over $K$ and is denoted by ${\rm Gal}(L/K)$. This means that for $\sigma\in\mathrm{Gal}(L/K)$, $\sigma(a)=(\sigma(a))^{\prime}$ for all $a\in L$ and for all $a\in K,$ $\sigma(a)=a$.

Assume that $\{y_{1},y_{2}\}$ is a fundamental system (basis) of solutions of $\mathcal{L}$. If $\sigma\in\mathrm{Gal}(L/K)$ then $\{\sigma y_{1},\sigma y_{2}\}$ is another fundamental system of $\mathcal{L}$. Hence there exists a matrix

[TABLE]

such that

[TABLE]

In a natural way, we can extend this to systems:

[TABLE]

This defines a faithful representation $\mathrm{Gal}(L/K)\to\mathrm{GL}(2,\mathbb{C})$ and it is possible to consider $\mathrm{Gal}(L/K)$ as a subgroup of $\mathrm{GL}(2,\mathbb{C})$. It depends on the choice of the fundamental system $\{y_{1},y_{2}\}$, but only up to conjugacy.

One of the fundamental results of the Picard-Vessiot theory is the following theorem (see [14, 17]).

Theorem 2.4**.**

The differential Galois group $\mathrm{Gal}(L/K)$ is an algebraic subgroup of $\mathrm{GL}(2,\mathbb{C})$.

Definition 2.5** (Integrability).**

Consider the linear differential equation $\mathcal{L}$ such as in equation (2.1). We say that $\mathcal{L}$ is integrable if the Picard-Vessiot extension $L\supset K$ is obtained as a tower of differential fields $K=L_{0}\subset L_{1}\subset\cdots\subset L_{m}=L$ such that $L_{i}=L_{i-1}(\eta)$ for $i=1,\ldots,m$, where either

$\eta$* is algebraic over $L_{i-1}$, that is $\eta$ satisfies a polynomial equation with coefficients in $L_{i-1}$.* 2. 2.

$\eta$* is primitive over $L_{i-1}$, that is $\partial_{x}\eta\in L_{i-1}$.* 3. 3.

$\eta$* is exponential over $L_{i-1}$, that is $\partial_{x}\eta/\eta\in L_{i-1}$.*

We remark that the usual terminology in differential algebra for integrable equations is that the corresponding Picard-Vessiot extensions are called Liouvillian.

Theorem 2.6** (Kolchin).**

The equation $\mathcal{L}$ given in (2.1) is integrable if and only if $(\mathrm{Gal}(L/K))^{0}$ is solvable.

Consider the differential equation

[TABLE]

We recall that equation (2.2) can be obtained from equation (2.1) through the change of variable

[TABLE]

and equation (2.2) is called the reduced form (also known as invariant normal form) of equation (2.1).

On the other hand, introducing the change of variable $v=\partial_{x}\zeta/\zeta$ we get the associated Riccati equation to equation (2.2),

[TABLE]

where $r$ is given by equation (2.3). Moreover, the Riccatti equation (2.4) has one algebraic solution over the differential field $K$ if and only if the differential equation (2.2) is integrable.

For $\mathcal{L}$ given by equation (2.2), it is very well known that ${\rm Gal}_{K}(\mathcal{L})$ is an algebraic subgroup of ${\rm SL}(2,\mathbb{C})$. The well known classification of subgroups of $\mathrm{SL}(2,\mathbb{C})$ is the following.

Theorem 2.7**.**

Let $G$ be an algebraic subgroup of ${\rm SL}(2,\mathbb{C})$. Then, up to conjugation, one of the following cases occurs.

$G\subseteq\mathbb{B}$* and then $G$ is reducible and triangularizable.* 2. 2.

$G\nsubseteq\mathbb{B}$, $G\subseteq\mathbb{D}_{\infty}$ and then $G$ is imprimitive. 3. 3.

$G\in\{A_{4}^{\mathrm{SL}_{2}},S_{4}^{\mathrm{SL}_{2}},A_{5}^{\mathrm{SL}_{2}}\}$* and then $G$ is primitive (finite)* 4. 4.

$G={\rm SL}(2,\mathbb{C})$* and then $G$ is primitive (infinite).*

2.2 Kovacic’s Algorithm

In 1986, Kovacic ([18]) introduced an algorithm to solve the differential equation (2.2), where $K=\mathbb{C}(x)$, showing that (2.2) is integrable if and only if the solution of the Riccati equation (2.4) is a rational function (case 1), is a root of polynomial of degree two (case 2) or is a root of polynomial of degree 4, 6, or 12 (case 3). We leave the details of the algorithm to Appendix A. We summarize here the main result by Kovacic as the following theorem.

Theorem 2.8** (Kovacic).**

There are precisely four cases that can occur for equation (2.2):

Case 1

It has a solution of the form $e^{\int\omega}$ where $\omega\in\mathbb{C}(x)$.

Case 2

It has a solution of the form $e^{\int\omega}$ where $\omega$ is algebraic over $\mathbb{C}(x)$ of degree 2, and case 1 does not hold.

Case 3

All solutions of (2.2) are algebraic over $\mathbb{C}(x)$ of degree 4, 6 or 12 and cases 1 and 2 do not hold.

Case 4

The differential equation (2.2) has no Liouvillian solution.

In the following sections we will apply the algorithm to equations (1.7) and (1.9) using a Laplace transform acting over the time variable.

3 Linear rates: $(b,d)=(1,1)$

As mentioned in Section 1, we introduce the generating function $g(z,t)=\sum_{N=0}^{\infty}P_{N}(t)z^{N}$ to transform the discrete variable $N$ into the continuous variable $z$, $0\leq z\leq 1$. This converts the master equation into a PDE: if we multiply both sides of the master equation (1.4) by $z^{N}$ and sum over $N$, we get

[TABLE]

Recall that $P_{N}(t):=0$ for $N<0$. We now use the following straightforward identities:

• (i)

$\frac{\partial g}{\partial z}=\sum_{N=0}^{\infty}(N+1)P_{N+1}(t)z^{N}=\sum_{N=1}^{\infty}NP_{N}(t)z^{N-1}$,

• (ii)

$z\frac{\partial g}{\partial z}=\sum_{N=1}^{\infty}NP_{N}(t)z^{N}$,

• (iii)

$z^{2}\frac{\partial g}{\partial z}=\sum_{N=1}^{\infty}(N-1)P_{N-1}(t)z^{N}$,

• (iv)

$\sum_{N=0}^{\infty}P_{N-1}(t)z^{N}=\sum_{N=1}^{\infty}P_{N-1}(t)z^{N}=\sum_{N=0}^{\infty}P_{N}(t)z^{N+1}=zg(z,t)$,

to get the first-order PDE (1.7),

[TABLE]

The initial condition $P_{N}(0)=\delta_{N_{0}N}$ reduces to $g(z,0)=z^{N_{0}}$. Normalization of the probability distribution at any time, $\sum_{N=0}^{\infty}P_{N}(t)=1$, implies that $g(1,t)=1$. Then we have to solve (3.2) with the boundary conditions $g(z,0)=z^{N_{0}}$ and $g(1,t)=1$.

3.1 Solution via a Laplace transform

Although this case leads to a first-order PDE, which can be solved via the method of characteristics (see reference [16] for the application of this method to the $(1,1)$ case in a more general setting in which the coefficients $\beta$ and $\delta$ are functions of time), we calculate here the solution explicitly via the Laplace transform method to illustrate our methodology. Previously, however, it is convenient to clarify what kind of integrability we are considering in this work.

Let

[TABLE]

be a partial differential equation, $M$ being a linear differential operator in the one-dimensional spatial variable $z$. Then we can state the following

Problem. *Solve the PDE (3.3) subject to suitable boundary conditions, including the initial Cauchy problem $g(z,0)=g_{0}(z)$. *

Applying the Laplace transform with respect to time to (3.3), we obtain a family of linear ODE equations,

[TABLE]

parameterized by the complex parameter $s$. We will say that equation (3.4) is integrable if the homogeneous equation

[TABLE]

is integrable in the sense of Picard-Vessiot theory. This is natural, because from the general solution of the homogeneous equation we obtain the general solution of (3.4) by quadratures. Another approach to the integrability of (3.4) is by transforming it to an homogeneous equation: later we will point out an explicit example of this point of view. Of course, we are assuming here that the coefficients of $M$, and the function $g_{0}(z)$ belong to a suitable differential field $K$, for instance, the set of complex rational functions. Then

Definition 3.1**.**

We say that the equation (3.3) is integrable if the family of linear ODE equations (3.4) is integrable in the sense of the Picard-Vessiot theory for almost any complex $s$.

We remark that despite it is usually assumed that the Laplace transformed function of the variable $s$ is defined in some half plane of the complex variable $s$, we are assuming here that this function can be prolongated analytically to other values of $s$.

Now focus on the PDE (1.7). It is clearly integrable, according to definition 3.1, because the associated linear ODE (3.4) is a first order ODE, being the coefficient field $K=\mathbb{C}(x)$ —indeed, along the rest of the paper we will assume that the coefficient field is the set of rational functions $\mathbb{C}(x)$. Hence, we introduce the Laplace transform acting over the time dependence of the generating function as

[TABLE]

which transforms (1.7) into the following first-order ODE,

[TABLE]

where we regard $s$ as a parameter and primes denote derivatives with respect to $z$. We now focus on solving equation (3.6) for arbitrary values of $s$. Note also that equation (3.6) can be expressed as

[TABLE]

with

[TABLE]

This form of the ODE will be convenient later in our computations. We observe that, in the homogeneous part of equation (3.7), the point $z=\infty$ is an ordinary point. Moreover, when $\beta\neq\delta$ the points $z=1$ and $z=\frac{\delta}{\beta}$ are regular singular points, while when $\beta=\delta$ the point $z=1$ is a singularity of irregular type, see [3] for a detailed explanation about differential Galois theory of non-homogeneous equations.

From now on we shall consider these two cases ($\beta\neq\delta$ and $\beta=\delta$) separately. For $\beta\neq\delta$, the homogeneous equation can be solved immediately,

[TABLE]

Assume that $\delta>\beta$ (the calculations for the $\delta<\beta$ case are simple extensions of those provided here and are therefore left to Appendix B). Then equation (3.9) yields

[TABLE]

i.e.,

[TABLE]

where $C$ is an integration constant. Variation of the constant in equation (3.6) yields a first-order ODE for the unknown function $C(z)$,

[TABLE]

for which we impose $C(1)=0$ to avoid possible divergences in the generating function $g(z,t)$ at $z=1$ —recall that $g(z,t)$ has to be an analytic function of $z$ because the probability distribution $P_{N}(t)$ is to be determined through a series expansion of $g(z,t)$ about $z=0$, see equation (1.6). Therefore,

[TABLE]

Using (3.11) and (3.13) together, the Laplace transform of the generating function is expressed as

[TABLE]

In terms of the new variable $w(u):=\alpha\left(\frac{1-u}{\delta-\beta u}\right)$ with $\alpha:=\frac{\delta-\beta z}{1-z}$, the integral above can be written as

[TABLE]

After a second change of variable, $w(t):=e^{(\beta-\delta)t}$, we finally get

[TABLE]

which allows us to identify the generating function

[TABLE]

Integration of (1.7) for $\beta<\delta$ yields exactly the same expression (see Appendix B).

Now, considering $\beta=\delta$, equation (3.9) yields

[TABLE]

i.e.,

[TABLE]

where $C$ is again an integration constant. Variation of the constant gives again a first-order ODE for $C(z)$,

[TABLE]

for which we impose $C(1)=0$ to avoid divergences, as above. Therefore, the general solution of equation (3.6) is

[TABLE]

The previous function can be obtained through iterated partial integration and, for $N_{0}\in\mathbb{Z^{+}}$, the result belongs to the family of exponential integrals, denoted by Ei, which is valid for $\mathfrak{R}(z)>0$ —as in our case because $0\leq z\leq 1$. Ei functions are not elementary functions, see [1] for further details. But in fact, we are interested here in the inverse-Laplace transformed function, $g(z,t)$, that becomes an elementary function. So, by means of the change t(u)=\frac{1}{\delta}\big{(}\frac{1}{1-u}-\frac{1}{1-z}\big{)}, we obtain

[TABLE]

Then

[TABLE]

is the sought solution of (1.7) for $\beta=\delta$, satisfying the boundary conditions $g(z,0)=z^{N_{0}}$ and $g(1,t)=1$.

In summary, proposition 1.1 has been proved. We consider the $(b,d)=(1,1)$ case as completely solved since the probability distribution $P_{N}(t)$ could eventually be obtained through a series expansion of the generating function. In particular, useful expressions for the mean and the variance of the distribution (or even any moment) can be computed for arbitrary values of $N$ and $t$. In addition, the probability of extinction at time $t$ is given by

[TABLE]

for $\beta\neq\delta$, and

[TABLE]

for $\beta=\delta$.

3.2 Solution via Kovacic’s algorithm

For quadratic death rates we find a second-order PDE for the generating function, see equation (1.9) and Section 4. The integrability of this case can be analyzed using Kovacic’s algorithm [18] since the Laplace transform yields a second-order, linear ODE whose coefficients are rational functions. As we anticipated in Proposition 1.2, the $(b,d)=(1,2)$ PDE is not integrable. However, we have just shown that, for the linear-rate case $(b,d)=(1,1)$, the problem is integrable. Kovacic’s algorithm usually restricts the values of the parameters in the differential equation in order to yield integrability. In both cases, the Laplace transform method introduces a new parameter in the equations —the parameter $s$ associated to the time dependence.

In this section we apply the algorithm by Kovacic to the $(b,d)=(1,1)$ case in order to gain some insight about integrability of the PDE via the Laplace transform: obviously, we have to recover the solution (3.11) with no restrictions imposed by the algorithm on the Laplace transform parameter $s$, in agreement with our definition of integrability.

We can apply Kovacic’s algorithm to a second-order, linear ODE whose coefficients are rational functions. In order to apply Kovacic’s algorithm to the inhomogeneous, first-order ODE (3.6), we transform the equation as follows: first eliminate the first derivative,

[TABLE]

and then divide the equation by the term $\frac{z^{N_{0}}}{(1-z)(\delta-\beta z)}$ to get

[TABLE]

Differentiating both sides of the equation above yields a second-order, linear, homogeneous equation whose coefficients are rational functions of $z$:

[TABLE]

Now it is convenient to clarify the relation between the solutions of the linear equation (3.26) and of the second order (3.28), that we write as a lemma for future reference.

Lemma 3.2**.**

Consider a first order linear ODE,

[TABLE]

with general solution

[TABLE]

Then the general solution of the associated second order, linear ODE obtained by derivation over equation (3.29) divided by $h$,

[TABLE]

is given by

[TABLE]

Proof 3.3**.**

A first integral of equation (3.31) is given by the linear first order equation

[TABLE]

which coincides with (3.29) for $C_{2}=1$. Then solving equation (3.33), we obtain (3.32).

In other words, a fundamental system of solutions of (3.31) is given by a non-trivial solution of the homogeneous part of (3.29) and by any of the particular solutions of (3.29) (like the one obtained by variation of constants). In particular, (3.31) has always a solution given by the exponential of an integral: $e^{\int fdz}$.

Now we focus on the solutions of the second-order ODE (3.28) yielded by Kovacic’s algorithm. For that purpose we normalize (3.28) to write it in the form $H^{\prime\prime}-r(z,s)H=0$ for a new function $H(z,s)$. If we define

[TABLE]

then the invariant normal form of (3.28) is obtained using equation (2.3):

[TABLE]

where $G(z,s)=H(z,s)\psi(z,s)$ and $\psi(z,s)$ satisfies the first-order ODE

[TABLE]

Note also that, for $\beta\neq\delta$,

[TABLE]

while for $\beta=\delta$,

[TABLE]

Integration of (3.36) yields

[TABLE]

and

[TABLE]

Now we apply Kovacic’s algorithm to (3.28) to check the integrability of this equation. Together with (3.39), we will construct solutions for the Laplace transform of the generating function as $G(z,s)=H(z,s)\psi(z,s)$. We recall that, by Lemma 3.2, the equation (3.28) has always a solution given by the exponential of an integral of a rational function,

[TABLE]

Then equation (3.28) has a solution given by the exponential of an integral in $K$,

[TABLE]

This implies that case 1 of Kovacic’s algorithm always holds for equation (3.35).

The computations that lead to the closed-form solution of equation (3.35) go as follows (see Appendix A for details on how the algorithm proceeds in a general setup). As can be easily checked, the rational function

[TABLE]

has three finite singularities at $z=0$, $z=1$ and $z=\delta/\beta$ if $\beta\neq\delta$. The algorithm is based on the orders of the poles of $r(z,s)$ in the complex plane, considering the singularity $z=\infty$ as well. Let $\Gamma^{\prime}$ be the set of finite poles of $r(z,s)$ in the complex plane. Let $\Gamma=\Gamma^{\prime}\cup\{\infty\}$. The method is based on the Laurent series expansions of $r(z,s)$ about the singularities in $\Gamma$ (see Appendix A). Let $\circ(c)$ denote the order of the pole $c$ in the Laurent series expansion. In this example, the following series expansions hold:

• (i)

$r(z,s)=\frac{N_{0}(N_{0}+2)}{4z^{2}}+\dots$ about $z=0$.

• (ii)

$r(z,s)=\frac{1}{4}\left(-1+\frac{s^{2}}{(\delta-\beta)^{2}}\right)\frac{1}{(z-1)^{2}}+\dots$ about $z=1$.

• (iii)

$r(z,s)=\frac{1}{4}\left(-1+\frac{s^{2}}{(\delta-\beta)^{2}}\right)\frac{1}{(z-\delta/\beta)^{2}}+\dots$ about $z=\frac{\delta}{\beta}$.

• (iv)

$r(z,s)=\frac{N_{0}(N_{0}-2)}{4z^{2}}+\dots$ about $z=\infty$.

We study the existence of case 1 solutions in Kovacic’s algorithm: all the poles have order 2, hence $\circ(c)=2$ for all $c\in\Gamma$. Therefore we write the Laurent series expansion of $\sqrt{r}$ about $c$, $[\sqrt{r}]_{c}$, as $[\sqrt{r}]_{c}=0$ for all $c\in\Gamma$ (see details about the general notation used in case 1 in Appendix A). Then we compute

[TABLE]

where $b$ is the residue of $r$ at the singularity $c$ (Appendix A). For $z=0$ we obtain $\alpha_{0}^{+}=1+\frac{N_{0}}{2}$ and $\alpha_{0}^{-}=-\frac{N_{0}}{2}$. For $z=1$, we get $\alpha_{1}^{\pm}=\frac{1}{2}\left(1\pm\frac{s}{\delta-\beta}\right)$. For $z=\delta/\beta$ we obtain $\alpha_{\delta/\beta}^{\pm}=\alpha_{1}^{\pm}$ because the residues associated to both singularities coincide. Finally, for $z=\infty$ we obtain $\alpha_{\infty}^{+}=\frac{N_{0}}{2}$ and $\alpha_{0}^{-}=1-\frac{N_{0}}{2}$.

We now form the $2^{4}$ possible permutations of signs for the four singularities and compute the quantity $m=\alpha_{\infty}^{\varepsilon(\infty)}-\sum_{c\in\Gamma^{\prime}}\alpha_{c}^{\varepsilon(c)}$. Let $\hat{s}:=\frac{s}{\delta-\beta}$. Then

According to the algorithm, we have to consider only those permutations that yield a non-negative integer $m$ (Appendix A). This discards, for example, the cases $(\varepsilon(\infty),\varepsilon(\delta/\beta),\varepsilon(1),\varepsilon(0))=(+,\pm,\mp,+)$ and $(-,\pm,\mp,+)$ —recall that $N_{0}$ is a non-negative integer number. It is possible to find integrability for the $(+,\pm,\mp,-)$ and $(-,\pm,\mp,-)$ cases. The remaining $8$ cases depend explicitly on $s$, hence imposing that $m$ is a non-negative integer would restrict the possible values of $s$ yielding closed-form solutions. In principle we are interested in finding a solution valid for any value of $s$, so we do not enter in the discussion of these cases in the main text. We leave the analysis of some of them for Appendix C, including computations for the sign combinations $(+,\pm,\mp,-)$.

Here we will discuss only one of the potential $4$ cases that can yield a solution: consider the permutation $(-,+,-,-)$, which corresponds to $m=0$. Then the algorithm proceeds by considering the rational function

[TABLE]

In our example this function reduces to

[TABLE]

Given that $\alpha_{0}^{-}=-\frac{N_{0}}{2}$, $\alpha_{1}^{-}=\frac{1}{2}\left(1-\frac{s}{\delta-\beta}\right)$ and $\alpha_{\delta/\beta}^{+}=\frac{1}{2}\left(1+\frac{s}{\delta-\beta}\right)$ we get, according to equation (3.37),

[TABLE]

The algorithm now searches for a monic polynomial $P_{m}(z)$ of degree $m$ that satisfies the differential equation

[TABLE]

If such polynomial exists, then a solution of the form $P_{m}e^{\int\omega}$ exists. In our case $m=0$ and, as it can be easily checked using (3.41), the function $\omega(z,s)$ defined in equation (3.45) satisfies identically the condition $\omega^{\prime}+\omega^{2}-r=0$. Therefore equation (3.46) is satisfied by the constant monic polynomial $P_{0}=1$ and we find the following closed-form solution for (3.35):

[TABLE]

Therefore, using (3.39), we get

[TABLE]

and we recover the solution (3.11) obtained from the first-order homogeneous ODE. Note that $G(z,s)=\left(\frac{\delta-\beta z}{1-z}\right)^{\frac{s}{\delta-\beta}}=e^{\int f(z,s)dz}$, with $f(z,s)$ given by (3.8). But this solution is not the Laplace transform $G(z,s)$ —equation (3.14)— of the generating function $g(z,t)$ we are looking for. However, we can use (3.48) and Lemma 3.2 with $f$ and $h$ given by (3.8) to construct the relevant solution as

[TABLE]

which is the exact same solution obtained in (3.14). Alternatively, it would be also possible to obtain this solution by applying to equation (3.28) the D’Alambert order reduction of a linear equation when a particular solution is known —we, however, skip the details here.

In a similar way we apply Kovacic’s algorithm for $\beta=\delta$. Now, equation (3.41) becomes

[TABLE]

Applying the case 1 of Kovacic’s algorithm we obtain that the solution of $H^{\prime\prime}=r(z,s)H$ is

[TABLE]

Now, using equation (3.40), we conclude that

[TABLE]

as in (3.19). We can recover the sought Laplace transform (3.21) using Lemma 3.2 as presented above.

An important insight that we infer thanks to the analysis of the first-order equation via Kovacic’s algorithm is the following conjecture: if we were to obtain integrability of the corresponding PDE via a Laplace transform, we conjecture that a necessary condition to obtain solutions of the form of Kovacic’s first case is that the combination

[TABLE]

remains independent of $s$, as our definition of integrability for equation (3.3) requires integrability of the linear ODE (3.4) for any value of the parameter $s$.

Remark: we observe that for $\beta\neq\delta$ equation (3.28) has 4 singular regular points at $z=0$, $z=1$, $z=\beta/\delta$ and $z=\infty$. Therefore, it corresponds exactly to the general Heun’s differential equation in the independent variable $z$ with parameters $\delta/\beta$, $sN_{0}/\beta$, [math], $1-N_{0}$, $-N_{0}$, and $(\beta-\delta+s)/(\beta-\delta)$. On the other hand, when $\beta=\delta$, we can observe that this equation has two regular singularities at $z=0$ and $z=\infty$, while it has one irregular singularity at $z=1$. We conclude that, with the changes of variables $G\mapsto Gz^{-N_{0}-1}/(1-z)^{2}$ and $z\mapsto(z-1)/z$, the equation corresponds to the confluent Heun’s differential equation with parameters $s/\delta$, $N_{0}-1$, $N_{0}+1$, [math], $(N_{0}^{2}\delta-sN_{0}+\delta)/(2\delta)$. Moreover, we observe that in the non-homogeneous first order linear differential equation the points $z=\infty$ and $z=0$ are ordinary points, but with the procedure to transform it into an homogeneous second order linear differential equation the points $z=\infty$ and $z=0$ are regular singular points. The type of singularity of $z=1$ and $z=\beta/\delta$ is preserved under such procedure for the cases $\beta=\delta$ and $\beta\neq\delta$, though. For further details about Heun’s differential equations, we refer the reader to reference [24]. We remark that a complete characterization of the integrability of Heun’s equations is today an open problem. Here it was possible to solve the integrability problem because the equations correspond to very special subfamilies of Heun’s general families.

4 Mixed rates: $(b,d)=(1,2)$

As we have shown in the previous section, Kovacic’s algorithm turns out to be a powerful tool to analyze the integrability of PDE associated to birth-death processes via a Laplace transform. In biological terms, a relevant case arises when mortality processes involve pairs of individuals, i.e., when the death rate is a quadratic function of the number of individuals. In this section we apply the same technology to analyze the integrability of the PDE associated to this situation.

As in the case of linear birth and death rates, we start by finding the PDE satisfied by the generating function when the birth rate is linear, $B_{N}=\beta N$, and the mortality rate is a quadratic function of $N$, $D_{N}=\delta N^{2}$. The generating function satisfies

[TABLE]

where we have used that $P_{N}(t):=0$ for $N<0$. The following identities hold:

• (i)

$\frac{\partial g}{\partial z}+z\frac{\partial^{2}g}{\partial z^{2}}=\sum_{N=1}^{\infty}N^{2}P_{N}(t)z^{N-1}=\sum_{N=0}^{\infty}(N+1)^{2}P_{N+1}(t)z^{N}$,

• (ii)

$z\frac{\partial g}{\partial z}=\sum_{N=1}^{\infty}NP_{N}(t)z^{N}$,

• (iii)

$z^{2}\frac{\partial g}{\partial z}=\sum_{N=1}^{\infty}(N-1)P_{N-1}(t)z^{N}$.

Therefore, we can express

[TABLE]

and

[TABLE]

As a consequence, we obtain a second-order PDE to be satisfied by the generating function, see equation (1.9). Similarly, we impose here the initial condition $g(z,0)=z^{N_{0}}$ and the normalization condition $g(1,t)=1$. In order to find solutions of equation (1.9), we follow the same procedure as for the $(1,1)$ case: we introduce the Laplace transform $G(z,s)$ of the generating function and try to solve the parametric ODE satisfied by $G(z,s)$ for arbitrary values of $s$. In terms of $G(z,s)$, the ODE reads

[TABLE]

where, again, primes denote derivatives with respect to $z$. Here we denote

[TABLE]

hence (4.4) can be expressed as

[TABLE]

In order to find the invariant normal form of (4.6), we write $G(z,s)=H(z,s)\psi(z)$ and impose that $\psi(z)$ satisfies the first-oder ODE

[TABLE]

Note that, in this case, $\psi$ is independent of $s$. Integration yields $\psi(z)=z^{-1/2}e^{\beta z/2\delta}$. Hence (4.6) reduces to the following second-order, non-homogeneous ODE for function $H(z,s)$:

[TABLE]

Equivalently, $H(z,s)$ satisfies

[TABLE]

This is the second-order normal invariant form of the original ODE. We want to see whether we can find closed-form solutions for this ODE for any value of the parameter $s$.

Remark: this equation has 3 singular points (as shown below) at $z=0$, $z=1$ (regular ones), and $z=\infty$ (irregular). Therefore, it belongs to the family of Heun’s confluent equations [24]. Thus, the $(b,d)=(1,2)$ case also belongs to Heun’s families.

4.1 Solution via Kovacic’s algorithm

In line with our definition of integrability (definition 3.1), we look for closed-form solutions of the homogeneous part of equation (4.9). For that purpose we define

[TABLE]

and apply Kovacic’s algorithm to search for closed-form solutions. In our case, $\Gamma^{\prime}=\{0,1\}$ with orders $\circ(0)=2$ and $\circ(1)=1$. The following series expansions for $r(z,s)$ about the three elements in $\Gamma$ hold:

• (i)

$r(z,s)=-\frac{1}{4z^{2}}+\dots$ about $z=0$.

• (ii)

$r(z,s)=-\frac{s}{\delta(z-1)}+\dots$ about $z=1$.

• (iii)

$r(z,s)=\frac{\beta^{2}}{4\delta^{2}}-\frac{\beta}{2\delta z}+\dots$ about $z=\infty$.

We observe in equation (4.10) that the order of $r$ at $\infty$ is $\circ(\infty)=0$. We analyze the different cases in the algorithm by Kovacic (see further details in Appendix A):

Case 1

Since $\circ(0)=2$, we set $[\sqrt{r}]_{0}=0$ and obtain $\alpha_{0}^{\pm}=\frac{1}{2}\pm\frac{1}{2}\sqrt{1+4b}=\frac{1}{2}$ because the residue at $z=0$ is $b=-\frac{1}{4}$.

For $z=1$, since $\circ(1)=1$, we set $[\sqrt{r}]_{1}=0$ and $\alpha_{1}^{\pm}=1$.

For $z=\infty$, since $\circ(\infty)=0=-2\nu$ and $r(z,s)=q^{2}+b/z+\dots$, with $q=\frac{\beta}{2\delta}$ and $b=-\frac{\beta}{2\delta}$, we set $[\sqrt{r}]_{\infty}=q=\frac{\beta}{2\delta}$ and $\alpha_{\infty}^{\pm}=\left(\pm\frac{b}{q}-\nu\right)/2=\mp\frac{1}{2}$.

We now consider all the possible combinations of signs for the three points:

Since all the values of $m$ are negative integers, Kovacic’s algorithm does not find solutions of the form $P_{m}e^{\int\omega}$ with $P_{m}$ a polynomial.

Case 2

Given the orders of the singularities of $r(z,s)$, we define the following subsets of $\mathbb{Z}$ (see a full description of how the algorithm proceeds in this case in Appendix A):

For $z=0$, since $\circ(0)=2$ and the residue at $z=0$ is $b=-\frac{1}{4}$, we have $E_{0}:=\{2+k\sqrt{1+4b},k=0,\pm 2\}=\{2\}$.

For $z=1$, since $\circ(1)=1$, we define $E_{1}:=\{4\}$.

For $z=\infty$, since $\circ(\infty)=0<2$, then $E_{\infty}:=\{0\}$.

We now find the positive combinations of the sum $m=\frac{1}{2}\left(e_{\infty}-\sum_{c\in\Gamma^{\prime}}e_{c}\right)$ for $e_{p}\in E_{p}$, $p\in\Gamma$. The only combination is $(e_{0},e_{1},e_{\infty})=(2,4,0)$, hence $m=\frac{1}{2}(0-2-4)<0$. Therefore the set of positive $m$ is empty and there are no solutions in this case.

Case 3

A necessary condition for this case to work is that $\circ(\infty)\geq 2$, see [18]. There are no solutions of this type since $\circ(\infty)=0$.

Therefore, we conclude that equation (4.9) is non-integrable for any value of $s$. Hence, the homogeneous part of equation (4.4) is also non-integrable and, as a consequence, the PDE (1.9) becomes non-integrable as well. This proves proposition 1.2.

Appendix A

This Appendix describes Kovacic’s algorithm in detail. In our presentation here, we follow the original version given by Kovacic in reference [18] with an adapted version presented in [2, 4].

Each case in Kovacic’s algorithm is related with each one of the algebraic subgroups of ${\rm SL}(2,\mathbb{C})$ and the associated Riccatti equation

[TABLE]

According to Theorem 2.7, there are four cases in Kovacic’s algorithm. Only for cases 1, 2 and 3 we can solve the differential equation, but for the case 4 the differential equation is not integrable. It is possible that Kovacic’s algorithm can provide us only one solution ($\zeta_{1}$), so that we can obtain the second solution ($\zeta_{2}$) through

[TABLE]

Notations. For the differential equation given by

[TABLE]

we use the following notations.

Denote by $\Gamma^{\prime}$ be the set of (finite) poles of $r$, $\Gamma^{\prime}=\left\{c\in\mathbb{C}:t(c)=0\right\}$. 2. 2.

Denote by $\Gamma=\Gamma^{\prime}\cup\{\infty\}$. 3. 3.

By the order of $r$ at $c\in\Gamma^{\prime}$, $\circ(r_{c})$, we mean the multiplicity of $c$ as a pole of $r$. 4. 4.

By the order of $r$ at $\infty$, $\circ\left(r_{\infty}\right),$ we mean the order of $\infty$ as a zero of $r$. That is $\circ\left(r_{\infty}\right)=\mathrm{deg}(t)-\mathrm{deg}(s)$.

**The four cases

**

Case 1. In this case $\left[\sqrt{r}\right]_{c}$ and $\left[\sqrt{r}\right]_{\infty}$ stand for the Laurent series of $\sqrt{r}$ at $c$ and the Laurent series of $\sqrt{r}$ at $\infty$ respectively. Furthermore, we define $\varepsilon(p)$ as follows: if $p\in\Gamma,$ then $\varepsilon\left(p\right)\in\{+,-\}.$ Finally, the complex numbers $\alpha_{c}^{+},\alpha_{c}^{-},\alpha_{\infty}^{+},\alpha_{\infty}^{-}$ will be defined in the first step. If the differential equation has no poles it only can fall in this case.

Step 1. For each $c\in\Gamma^{\prime}$ and for $\infty$ consider the following possibilities:

$(c_{0})$

If $\circ\left(r_{c}\right)=0$, then

[TABLE]

$(c_{1})$

If $\circ\left(r_{c}\right)=1$, then

[TABLE]

$(c_{2})$

If $\circ\left(r_{c}\right)=2,$ and

[TABLE]

[TABLE]

$(c_{3})$

If $\circ\left(r_{c}\right)=2v\geq 4$, and

[TABLE]

[TABLE]

$(\infty_{1})$

If $\circ\left(r_{\infty}\right)>2$, then

[TABLE]

$(\infty_{2})$

If $\circ\left(r_{\infty}\right)=2,$ and $r=\cdots+bx^{2}+\cdots$, then

[TABLE]

$(\infty_{3})$

If $\circ\left(r_{\infty}\right)=-2v\leq 0$, and

[TABLE]

[TABLE]

Step 2. Find $D\neq\emptyset$ defined by

[TABLE]

If $D=\emptyset$, then we should start with the case 2. Now, if $\mathrm{Card}(D)>0$, then for each $n\in D$ we search $\omega$ $\in\mathbb{C}(x)$ such that

[TABLE]

Step 3. For each $n\in D$, search for a monic polynomial $P_{n}$ of degree $n$ with

[TABLE]

If success is achieved then $\zeta_{1}=P_{n}e^{\int\omega}$ is a solution of the differential equation. Otherwise, case 1 cannot hold.

Case 2.

Step 1. For each $c\in\Gamma^{\prime}$ and $\infty$ compute non-empty sets $E_{c}\subset\mathbb{Z}$ and $E_{\infty}\subset\mathbb{Z}$ defined as follows:

($c_{1}$)

If $\circ\left(r_{c}\right)=1$, then $E_{c}=\{4\}$.

($c_{2}$)

If $\circ\left(r_{c}\right)=2,$ and r=\cdots+b(x-c)^{-2}+\cdots,\ then

[TABLE]

($c_{3}$)

If $\circ\left(r_{c}\right)=v>2$, then $E_{c}=\{v\}$.

$(\infty_{1})$

If $\circ\left(r_{\infty}\right)>2$, then $E_{\infty}=\{0,2,4\}$.

$(\infty_{2})$

If $\circ\left(r_{\infty}\right)=2,$ and $r=\cdots+bx^{2}+\cdots$, then

[TABLE]

$(\infty_{3})$

If $\circ\left(r_{\infty}\right)=v<2$, then $E_{\infty}=\{v\}$.

Step 2. Find $D\neq\emptyset$ defined by

[TABLE]

If $D=\emptyset,$ then we should start the case 3. Now, if $\mathrm{Card}(D)>0,$ then for each $n\in D$ we search a rational function $\theta$ defined by

[TABLE]

Step 3. For each $n\in D,$ search for a monic polynomial $P_{n}$ of degree $n$, such that

[TABLE]

If $P_{n}$ does not exist, then case 2 cannot hold. If such a polynomial is found, set $\phi=\theta+P_{n}^{\prime}/P_{n}$ and let $\omega$ be a solution of

[TABLE]

Then $\zeta_{1}=e^{\int\omega}$ is a solution of the differential equation.

Case 3.

Step 1. For each $c\in\Gamma^{\prime}$ and $\infty$ compute non-empty sets $E_{c}\subset\mathbb{Z}$ and $E_{\infty}\subset\mathbb{Z}$ defined as follows:

$(c_{1})$

If $\circ\left(r_{c}\right)=1$, then $E_{c}=\{12\}$.

$(c_{2})$

If $\circ\left(r_{c}\right)=2,$ and $r=\cdots+b(x-c)^{-2}+\cdots$, then

[TABLE]

$(\infty)$

If $\circ\left(r_{\infty}\right)=v\geq 2,$ and $r=\cdots+bx^{2}+\cdots$, then

[TABLE]

Step 2. Find $D\neq\emptyset$ defined by

[TABLE]

In this case we start with $m=4$ to obtain the solution, afterwards $m=6$ and finally $m=12$. If $D=\emptyset$, then the differential equation is not integrable because it falls in case 4. Now, if $\mathrm{Card}(D)>0,$ then for each $n\in D$ with its respective $m$, search for a rational function

[TABLE]

and a polynomial $S$ defined as

[TABLE]

Step 3. For each $n\in D$, with its respective $m$, search for a monic polynomial $P_{n}=P$ of degree $n,$ such that $P$ can be determined by the following polynomial recursion:

[TABLE]

This can be done by using undetermined coefficients for $P$. If $P$ does not exist, then the differential equation is not integrable because it falls in case 4. Now, if $P$ exists search $\omega$ such that

[TABLE]

then a solution of the differential equation is given by

[TABLE]

where $\omega$ is solution of the previous polynomial equation of degree $m$.

Appendix B

In this appendix we show that the generating function for the case of linear death rates is given also by equation (1.10) if we assume $\delta<\beta$. Here we consider two separate cases:

• (i)

$0\leq z\leq\frac{\delta}{\beta}$: here (3.9) reduces to

[TABLE]

hence

[TABLE]

$C(z)$ being the constant obtained after integration. Variation of the constant in equation (3.6) yields the first-order ODE for $C(z)$,

[TABLE]

We impose the condition $C(\delta/\beta)=0$ for $G(z,s)$ to be non-singular at $z=\frac{\delta}{\beta}<1$. Hence

[TABLE]

Then the Laplace transform of the generating function can be written as

[TABLE]

We change variable $u$ to $w(u):=\alpha\left(\frac{\delta-\beta u}{1-u}\right)$ with $\alpha:=\frac{1-z}{\delta-\beta z}$ and obtain

[TABLE]

Finally we introduce a second change of variable, $w(t):=e^{-(\beta-\delta)t}$, which yields

[TABLE]

and the generating function is expressed as

[TABLE]

which exactly coincides with the expression obtained in Section 3.1.

• (ii)

$\frac{\delta}{\beta}\leq z\leq 1$: in this case we can write

[TABLE]

i.e.,

[TABLE]

Variation of the constants implies

[TABLE]

which can be integrated as

[TABLE]

(notice the condition $C(\delta/\beta)=0$ for $G(z,s)$ to be finite at $z=\frac{\delta}{\beta}<1$). We can write

[TABLE]

We change variables to $w(u):=\alpha\left(\frac{\beta u-\delta}{1-u}\right)$ with $\alpha:=\frac{1-z}{\beta z-\delta}$,

[TABLE]

and after a second change of variable, $w(t):=e^{-(\beta-\delta)t}$, we finally obtain

[TABLE]

The generating function, in this case, is

[TABLE]

which coincides with (1.10).

Appendix C

In this Appendix we analyze in detail some cases of Kovacic’s algorithm applied to equation (3.35). We are particularly interested in those cases that impose a restriction on the values that $s$ can take to yield a closed-form solution, as well as in the cases where $m=N_{0}-1$ is a nonnegative integer. We start considering the combination of signs $(\varepsilon(\infty),\varepsilon(\delta/\beta),\varepsilon(1),\varepsilon(0))=(-,-,-,-)$. In this case, $\hat{s}=\frac{s}{\delta-\beta}$ has to be a non-negative integer. Let $m=\hat{s}$ be a non-negative integer. Then we form the rational function

[TABLE]

Hence it can be checked that

[TABLE]

We search for a polynomial of degree $m$ that satisfies $P_{m}^{\prime\prime}+2\omega P^{\prime}_{m}+(\omega^{\prime}+\omega^{2}-r)P_{m}=0$. It turns out that $P_{m}(z)=(\delta-\beta z)^{m}$ satisfies the equation and there is a solution of the form

[TABLE]

where $m=\frac{s}{\delta-\beta}$ is a non-negative integer. Using (3.39) we find the solution

[TABLE]

i.e., we obtain the solution given by equation (3.11) but specialized to non-negative integer values of the variable $m=\frac{s}{\delta-\beta}$.

Similarly, consider the combination of signs $(\varepsilon(\infty),\varepsilon(\delta/\beta),\varepsilon(1),\varepsilon(0))=(-,+,+,-)$. Here $-\hat{s}=-\frac{s}{\delta-\beta}$ has to be a non-negative integer. We define $m=-\hat{s}$ as a non-negative integer and we find, as before,

[TABLE]

Hence we obtain

[TABLE]

In this case the polynomial $P_{m}(z)=(1-z)^{m}$ satisfies (3.46) and there is a solution of the form

[TABLE]

Using (3.39) we finally get

[TABLE]

which coincides with (3.11) because $m=-\frac{s}{\delta-\beta}$ in this case. Again, we recover the same solution specialized to non-negative integer values of $m$.

Now we consider the combination of signs $(\varepsilon(\infty),\varepsilon(\delta/\beta),\varepsilon(1),\varepsilon(0))=(+,+,-,-)$. In this case $m=N_{0}-1$ has to be a non-negative integer; thus, we replace $N_{0}$ by $m+1$ in the following computations. Then we form the rational function

[TABLE]

Hence it can be checked that

[TABLE]

We search for a polynomial of degree $m$ that satisfies $P_{m}^{\prime\prime}+2\omega P^{\prime}_{m}+(\omega^{\prime}+\omega^{2}-r)P_{m}=0$ and we observe that for $m=0$, the polynomial $P_{0}=1$ does not satisfies such algebraic equation. Moreover, for $m>0$ the polynomial could exist with algebraic coefficients depending on $\beta$, $\delta$ and $s$. For example, the polynomial of degree $m=1$ is $P_{1}(z)=z-\delta/s$ only for $\delta=s-\beta$.

Consider finally the combination of signs $(\varepsilon(\infty),\varepsilon(\delta/\beta),\varepsilon(1),\varepsilon(0))=(+,-,+,-)$. Then we form the rational function

[TABLE]

Then it holds

[TABLE]

We search for a polynomial of degree $m$ that satisfies $P_{m}^{\prime\prime}+2\omega P^{\prime}_{m}+(\omega^{\prime}+\omega^{2}-r)P_{m}=0$ and we observe that for $m=0$, the polynomial $P_{0}=1$ satisfies such algebraic equation. In addition, for $m>0$ the polynomial exists with algebraic coefficients on $\beta$, $\delta$ and $s$. For example, we observe that the polynomial of degree $m=1$ is $P_{1}(z)=z$ when $\delta=0$ and $s=\beta-1$.

Remark: although $N_{0}$ must be non-negative for its meaning in the birth-death processes, we can apply Kovacic algorithm when $m=-N_{0}-1$ is the degree of the polynomial $P_{m}$ assuming that $-N_{0}\in\mathbb{Z}^{+}$. Thus, we obtain the same results as above because the combinations of signs $(\varepsilon(\infty),\varepsilon(\delta/\beta),\varepsilon(1),\varepsilon(0))=(+,+,-,-)$ and $(\varepsilon(\infty),\varepsilon(\delta/\beta),\varepsilon,\varepsilon(0))=(-,+,-,-)$ coincide with the combinations of signs

[TABLE]

and

[TABLE]

respectively.

Acknowledgements

The authors kindly thank to the members of our Integrability Madrid Seminar for many fruitful discussions: Rafael Hernández-Heredero, Sonia Jiménez-Verdugo, Alvaro Pérez-Raposo, José Rojo-Montijano, Sonia L. Rueda, Raquel Sánchez-Cauce and Maria A. Zurro.