Integrability of Stochastic Birth-Death processes via Differential Galois Theory

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.


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 ∈ N ∪ {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 ∆t. Then birth and death rates satisfy Pr{N + 1, t + ∆t | N, t} = B N ∆t, Pr{N − 1, t + ∆t | N, t} = D N ∆t, (1.1) where Pr{N + 1, t + ∆t | N, t} is the conditional probability that the system undergoes a birth event at time t + ∆t given that there were N individuals at time t. Multiple births and deaths are usually ignored in the limit ∆t → 0 because their probability would be proportional to (∆t) 2 .
If the population is formed by N individuals at time t, at time t + ∆t the population can be composed by: Thus, using the theorem of total probability we can write an expression for the probability of observing N individuals at time t + ∆t in terms of the probabilities at time t: 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 ∆t → 0 to get the so-called master equation: The system is therefore fully described by a coupled system of infinitely many ordinary differential equations, given by equation (1.4) for N ≥ 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 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 ∈ N∪{0}. Then the initial probability distribution is precisely equal to P N (0) = δ N 0 N , where δ 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 i.e., the discrete variable N is transformed into to continuous variable z, 0 ≤ z ≤ 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 = βN b and D N = δN d , where b, d are natural exponents and β, δ 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) ∈ {(1, 1), (1, 2)}. For the combination (b, d) = (1, 1), the master equation is equivalent to the following PDE (see Section 3): with boundary conditions 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, (1.9) 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. (1.

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 ∂ x := ∂ ∂x (also ′ := ∂ ∂x ) throughout this section.

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 ∂ x a derivation, that is, a map ∂ x : for all a, b ∈ K. By C we denote the field of constants of K, 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 ∂ x = ′ .
Up to special considerations, we analyze second order linear homogeneous differential equations, that is, equations in the form L := y ′′ + ay ′ + by = 0, a, b ∈ K.
(2.1) Definition 2.2 (Picard-Vessiot Extension). Suppose that y 1 , y 2 is a basis of solutions of 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 C of these two. Let L = K y 1 , y 2 = K(y 1 , y 2 , y ′ 1 , y ′ 2 ) the differential extension of K such that 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 L. Definition 2.3 (Differential Galois Groups). Assume K, L and 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 Gal(L/K). This means that for σ ∈ Gal(L/K), σ(a) = (σ(a)) ′ for all a ∈ L and for all a ∈ K, σ(a) = a.
Assume that {y 1 , y 2 } is a fundamental system (basis) of solutions of L. If σ ∈ Gal(L/K) then {σy 1 , σy 2 } is another fundamental system of L. Hence there exists a matrix In a natural way, we can extend this to systems: This defines a faithful representation Gal(L/K) → GL(2, C) and it is possible to consider Gal(L/K) as a subgroup of GL(2, 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.1. The differential Galois group Gal(L/K) is an algebraic subgroup of GL(2, C). Definition 2.4 (Integrability). Consider the linear differential equation L such as in equation (2.1). We say that L is integrable if the Picard-Vessiot extension L ⊃ K is obtained as a tower of differential fields We remark that the usual terminology in differential algebra for integrable equations is that the corresponding Picard-Vessiot extensions are called Liouvillian.
Consider the differential equation We recall that equation (2.2) can be obtained from equation (2.1) through the change of variable 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 = ∂ x ζ/ζ we get the associated Riccati equation to equation (2.2), 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 L given by equation (2.2), it is very well known that Gal K (L) is an algebraic subgroup of SL(2, C). The well known classification of subgroups of SL(2, C) is the following. Theorem 2.3. Let G be an algebraic subgroup of SL(2, C). Then, up to conjugation, one of the following cases occurs.
1. G ⊆ B and then G is reducible and triangularizable.

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.4 (Kovacic). There are precisely four cases that can occur for equation (2.2): Case 1 It has a solution of the form e ω where ω ∈ C(x).
Case 2 It has a solution of the form e ω where ω is algebraic over C(x) of degree 2, and case 1 does not hold.
Case 3 All solutions of (2.2) are algebraic over 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) = ∞ N =0 P N (t)z N to transform the discrete variable N into the continuous variable z, 0 ≤ z ≤ 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 Recall that P N (t) := 0 for N < 0. We now use the following straightforward identities: to get the first-order PDE (1.7), The initial condition P N (0) = δ N 0 N reduces to g(z, 0) = z N 0 . Normalization of the probability distribution at any time, ∞ N =0 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.

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 β and δ 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 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, parameterized by the complex parameter s. We will say that equation 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 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 = C(x)indeed, along the rest of the paper we will assume that the coefficient field is the set of rational functions C(x). Hence, we introduce the Laplace transform acting over the time dependence of the generating function as which transforms (1.7) into the following first-order ODE, 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 .
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 = ∞ is an ordinary point. Moreover, when β = δ the points z = 1 and z = δ β are regular singular points, while when β = δ 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 (β = δ and β = δ) separately. For β = δ, the homogeneous equation can be solved immediately, .
Assume that δ > β (the calculations for the δ < β case are simple extensions of those provided here and are therefore left to Appendix B). Then equation (3.9) yields i.e., 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), 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, (3.13) Using (3.11) and (3.13) together, the Laplace transform of the generating function is expressed as (3.14) In terms of the new variable w(u) := α 1−u δ−βu with α := δ−βz 1−z , the integral above can be written as After a second change of variable, w(t) := e (β−δ)t , we finally get which allows us to identify the generating function Integration of (1.7) for β < δ yields exactly the same expression (see Appendix B). Now, considering β = δ, equation (3.9) yields i.e., where C is again an integration constant. Variation of the constant gives again a first-order ODE for C(z), for which we impose C(1) = 0 to avoid divergences, as above. Therefore, the general solution of equation (3.6) is (3.21) The previous function can be obtained through iterated partial integration and, for N 0 ∈ Z + , the result belongs to the family of exponential integrals, denoted by Ei, which is valid for R(z) > 0 -as in our case because 0 ≤ z ≤ 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) = 1 is the sought solution of (1.7) for β = δ, 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 for β = δ, and for β = δ.

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, , (3.26) and then divide the equation by the term Differentiating both sides of the equation above yields a second-order, linear, homogeneous equation whose coefficients are rational functions of z: G(z, s) = 0. (3.28) 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.

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

Proof . A first integral of equation (3.31) is given by the linear first order equation
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 f dz . where G(z, s) = H(z, s)ψ(z, s) and ψ(z, s) satisfies the first-order ODE 2ψ ′ (z, s) + a(z, s)ψ(z, s) = 0. (3.36) Note also that, for β = δ, (3.38) Integration of (3.36) yields (iii) r(z, s) We study the existence of case 1 solutions in Kovacic' where b is the residue of r at the singularity c (Appendix A). For z = 0 we obtain α + 0 = 1 + N 0 2 and α − 0 = − N 0 2 . For z = 1, we get α ± 1 = 1 2 1 ± s δ−β . For z = δ/β we obtain α ± δ/β = α ± 1 because the residues associated to both singularities coincide. Finally, for z = ∞ we obtain α + ∞ = N 0 2 and α − 0 = 1 − N 0 2 . We now form the 2 4 possible permutations of signs for the four singularities and compute the quantity m = α According to the algorithm, we have to consider only those permutations that yield a nonnegative integer m (Appendix A). This discards, for example, the cases (ε(∞), ε(δ/β), ε(1), ε(0)) = (+, ±, ∓, +) and (−, ±, ∓, +) -recall that N 0 is a non-negative integer number. It is possible to find integrability for the (+, ±, ∓, −) and (−, ±, ∓, −) 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 (+, ±, ∓, −).
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 In our example this function reduces to we get, according to equation (3.37), The algorithm now searches for a monic polynomial P m (z) of degree m that satisfies the differential equation If such polynomial exists, then a solution of the form P m e ω exists. In our case m = 0 and, as it can be easily checked using (3.41), the function ω(z, s) defined in equation (3.45) satisfies identically the condition ω ′ + ω 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): Therefore, using (3.39), we get and we recover the solution (3.11) obtained from the first-order homogeneous ODE. Note that G(z, s) = δ−βz 1−z s δ−β = e 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.1 with f and h given by (3.8) to construct the relevant solution as 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 β = δ. Now, equation (3.41) becomes Applying the case 1 of Kovacic's algorithm we obtain that the solution of H ′′ = r(z, s)H is as in (3.19). We can recover the sought Laplace transform (3.21) using Lemma 3.1 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 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 β = δ equation (3.28) has 4 singular regular points at z = 0, z = 1, z = β/δ and z = ∞. Therefore, it corresponds exactly to the general Heun's differential equation in the independent variable z with parameters δ/β, sN 0 /β, 0, 1 − N 0 , −N 0 , and (β − δ + s)/(β − δ). On the other hand, when β = δ, we can observe that this equation has two regular singularities at z = 0 and z = ∞, while it has one irregular singularity at z = 1. We conclude that, with the changes of variables G → Gz −N 0 −1 /(1 − z) 2 and z → (z − 1)/z, the equation corresponds to the confluent Heun's differential equation with parameters s/δ, N 0 − 1, N 0 + 1, 0, (N 2 0 δ − sN 0 + δ)/(2δ). Moreover, we observe that in the non-homogeneous first order linear differential equation the points z = ∞ and z = 0 are ordinary points, but with the procedure to transform it into an homogeneous second order linear differential equation the points z = ∞ and z = 0 are regular singular points. The type of singularity of z = 1 and z = β/δ is preserved under such procedure for the cases β = δ and β = δ, 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. 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 = βN , and the mortality rate is a quadratic function of N , D N = δN 2 . The generating function satisfies where we have used that P N (t) := 0 for N < 0. The following identities hold: Therefore, we can express 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 where, again, primes denote derivatives with respect to z. Here we denote hence (4.4) can be expressed as . (4.6) In order to find the invariant normal form of (4.6), we write G(z, s) = H(z, s)ψ(z) and impose that ψ(z) satisfies the first-oder ODE Note that, in this case, ψ is independent of s. Integration yields ψ(z) = z −1/2 e βz/2δ . Hence (4.6) reduces to the following second-order, non-homogeneous ODE for function H(z, s): . (4.8) Equivalently, H(z, s) satisfies 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 = ∞ (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.

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 (ii) r(z, s) = − s δ(z−1) + . . . about z = 1.
We observe in equation (4.10) that the order of r at ∞ is •(∞) = 0. We analyze the different cases in the algorithm by Kovacic (see further details in Appendix A): . 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 ω with P m a polynomial. For z = ∞, since •(∞) = 0 < 2, then E ∞ := {0}.
We now find the positive combinations of the sum m = 1 2 e ∞ − c∈Γ ′ e c for e p ∈ E p , p ∈ Γ. The only combination is (e 0 , e 1 , e ∞ ) = (2, 4, 0), hence m = 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 •(∞) ≥ 2, see [18]. There are no solutions of this type since •(∞) = 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].

Case 2.
Step 1. For each c ∈ Γ ′ and ∞ compute non-empty sets E c ⊂ Z and E ∞ ⊂ Z defined as follows: (∞ 2 ) If • (r ∞ ) = 2, and r = · · · + bx 2 + · · · , then Step 2. Find D = ∅ defined by If D = ∅, then we should start the case 3. Now, if Card(D) > 0, then for each n ∈ D we search a rational function θ defined by Step 3. For each n ∈ D, search for a monic polynomial P n of degree n, such that If P n does not exist, then case 2 cannot hold. If such a polynomial is found, set φ = θ + P ′ n /P n and let ω be a solution of Then ζ 1 = e ω is a solution of the differential equation.
Step 2. Find D = ∅ defined by In this case we start with m = 4 to obtain the solution, afterwards m = 6 and finally m = 12.
If D = ∅, then the differential equation is not integrable because it falls in case 4. Now, if Card(D) > 0, then for each n ∈ D with its respective m, search for a rational function Step 3. For each n ∈ 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: P m = −P, , for i ∈ {m, m − 1, . . . , 1, 0}, P −1 = 0. 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 ω such that m i=0 S i P (m − i)! ω i = 0, then a solution of the differential equation is given by where ω 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 δ < β. Here we consider two separate cases:  and the generating function is expressed as which exactly coincides with the expression obtained in Section 3.1.
(ii) δ β ≤ z ≤ 1: in this case we can write