CELL DIVISION AND THE PANTOGRAPH EQUATION

Simple models for size structured cell populations undergoing growth and division produce a class of functional ordinary differential equations, called pantograph equations, that describe the long time asymptotics of the cell number density. Pantograph equations arise in a number of applications outside this model and, as a result, have been studied heavily over the last five decades. In this paper we review and survey the rôle of the pantograph equation in the context of cell division. In addition, for a simple case we present a method of solution based on the Mellin transform and establish uniqueness directly from the transform equation.


Introduction
A simple model for a size structured population of cells undergoing growth and division into α > 1 daughter cells of equal size is given by the functional partial differential equation n t (x, t) + (G(x)n(x, t)) x + B(x)n(x, t) = α 2 B(αx)n(αx, t), where n(x, t) is the number density of cells of size x at time t, and G and B are non negative functions that model the growth and division rates respectively.Here, "size" means either mass or DNA content, and α is generally 2 in applications.Solutions to equation (1) are required to satisfy an initial condition n(x, 0) = n 0 (x), (2) and the boundary conditions The model was first presented by Sinko and Streifer [33] in the 1960's and revisited by several researchers (e.g.[13,21]).Equation (1) along with conditions (2) and ( 3) is an example of an initial boundary value problem that features a partial differential equation containing a non local term.There are few analytical solution techniques for such problems, although an analytical solution has been found for the special case when the growth and division rates are constant [44].In fact, it is the long time behaviour of solutions that is of central interest in applications and this is the focus of most of the research on this model.The model has been placed in the framework of a fragmentation coagulation model for which there is a formidable array of results regarding the long time asymptotics of solutions (cf.[27], [28], [31], [14] and references therein).Briefly, it has been shown for constant growth rate and a large class of division rates that there exists an eigenvalue and a corresponding non negative eigenfunction y such that lim t→∞ e −λt n(x, t) − y(x) = 0, where • is a weighted norm.The details can be found in [31] and [32].The above results use a relative entropy method.We note that a related result was given by [9] using a different approach.We note also that the original motivation behind the work of Hall and Wake [21] to seek long time asymptotic solutions was experimental evidence in plant tissues, which showed that, regardless of the initial distribution n 0 , the long time shape of the distribution was the same [20].
The eigenfunction y satisfies the functional differential equation where denotes d/dx, along with the conditions It is further required that y be a non negative function and that (i.e.y is a probability density function).Equation ( 4) corresponds to the equation for y derived under the assumption of a separable solution n(x, t) = N (t)y(x).The remainder of this paper will focus primarily on the solution of equation ( 4) subject to conditions (5) and ( 6).We will refer to this as Problem A.

The pantograph equation
Suppose that y is a solution to Problem A and that By ∈ L 1 [0, ∞).Then integration of both sides of equation (4) from 0 to ∞ gives If, in addition, xy ∈ L 1 [0, ∞) and Gy ∈ L 1 [0, ∞), then multiplying both sides of (4) by x and integrating from 0 to ∞ yields The above calculations highlight two special cases when λ can be determined explicitly: (1) Constant division rate: If B(x) = b = const., then equation (7) shows that (2) Linear growth rate: If G(x) = gx, where g is a constant, then These calculations also highlight an obvious problem: suppose that we have a linear growth and a constant division rate, then Since g, b and α are independent parameters that depend on the biology of the cell, equation (11), in general, is not satisfied.The simplest case is when G and B are both constant functions.Let G(x) = g and B(x) = b.Then equation (4) reduces to y (x) + κy(x) = ακy(αx), where κ = bα/g.Equation ( 12) is an example of a pantograph equation.This equation has been studied extensively owing to the numerous applications in which it appears.In addition to the cell division model, this equation arises in models for light absorption in the Milky Way galaxy [2], the collection of current in an electric locomotive [30] and in a ruin problem [16] among other applications.Certainly one reason that the equation appears in many applications is that it models the limiting distribution of a Markov process.Derfel et al. [12] show this explicitly for the cell division model.In fact, this relationship is a special case of a more general connection between a class of functional equations and Markov processes ( [8,10]).Key features of the pantograph equation are that it contains a non local term of the form y(Λx), where Λ is a constant other than 1, and that the differential equation is linear in y, y and the non local term.Equation ( 4) is regarded as a generalization.The (original) pantograph equation is of the form where A, B and Λ are constants.Equation ( 12) is an example of an advanced functional differential equation, since The retarded version (|Λ| < 1) was studied much earlier by Flamant [18] in 1924 and Adams [1] in 1931.Starcher [34] used variants of this equation to prove identities among theta functions.It is also noted that this equation arises in an asymptotic estimate of the number of partitions in the Mahler partition problem [26].
The advanced and retarded versions of the pantograph equation are related through the Laplace transform, yet their solution behaviour is markedly different.Kato and McLeod [24] and later Iserles [23] studied the solutions to the pantograph equation subject to the initial condition y(0) = 1.If the equation is retarded, then the initial value problem has a unique solution.In fact this solution is an entire function.In contrast, if the equation is advanced then it is easy to show that there are no solutions holomorphic at x = 0.In fact, it can be shown that the imaginary axis forms a natural boundary for holomorphic continuation [19,38].More strikingly, it can be shown that the problem is ill posed in that there are an infinite number of solutions to the initial value problem.The asymptotic behaviour of solutions to the advanced problem as x → ∞ is very sensitive to the relative values of A and B. It is only when some sort of asymptotic behaviour is imposed on solutions that there is any link to uniqueness for the advanced problem.Note that equation (12) along with condition (6) can be converted into the initial value problem studied by Kato and McLeod by the transformation

The solution
Despite the bewildering possibilities for solutions in the general case, the problem arising from cell division is quite tractable.We have the following result Theorem 3.1.There exists a unique non negative solution to equation ( 12) that satisfies conditions ( 5) and (6).Moreover, this solution is unimodal.
The proof of existence and uniqueness can be found in [21]; a proof that the solution is unimodal can be found in [9].The condition that y be a probability density function is sufficient to ensure the uniqueness of a solution.The solution can be found explicitly in a number of ways.Kato and McLeod [24] presented a solution along with a host of other solutions for more general coefficients that reflect various assumptions on the asymptotic behaviour of solutions as x → ∞.Hall and Wake [21] used Laplace transforms to construct a solution.A related but distinct method of solution was also devised by Derfel et al. [12] using a sequence of Laplace convolutions.In practice, this last method is perhaps a less efficient method than simply using Laplace transforms for deriving the solution, but it does offer a greater insight into nature of the solution as a limiting distribution to a Markov process.Yet another method to derive the solution is based on the Mellin transform, and it is this method that we will outline as it has some novel features.
Recall that the Mellin transform of a function y is The Mellin transform has already been recognised as a useful tool for solving more general equations of the pantograph type [39].This transform is particularly useful when the growth and division rates are polynomials.
It also has the useful property that Once the transform is obtained, the asymptotic behaviour of the solution y as x → 0 + and as x → ∞ can be gleaned by position of the poles that bound the strip of holomorphy.For succinctness, let and condition (6) yields A specific solution M 0 (s) to ( 13) can be found as follows: let M 0 (s) = F (s)Q(s), where Note that F (s) is simply the Mellin transform to a nontrivial solution of the differential equation y h (x) + κy h (x) = 0, so that we can use Substituting this form of M 0 into equation (13) gives which has a solution We thus have the particular solution Since α > 1, the infinite product Q converges for all s ∈ C. The zeros for Q are of the form s = −m + 2πki/ ln α, where m is a non negative integer and k is an integer.These zeros are all of order 1.The gamma function has simple poles at s = 0, −1, −2, . .., where Q has a zero and we thus see that M 0 is an entire function.
It is possible at this stage to find the inverse Mellin transform of M 0 and then use equation ( 12) directly to show that, suitably normalised, M 0 gives the solution to the problem.We will invert this transform below, but it is of interest to examine the uniqueness of solutions to (13) directly.Determining the inverse transform of M 0 relies on a partition identity which is usually not available for more general problems.The problem that we consider is what if we solve a Mellin transform equation such as (13), and we can confirm that the solution has an inverse transform (e.g. through a Payley Wiener theorem).Knowledge of the transform gives us direct access to the asymptotics of the inverse transform and is thus useful in its own right.The problem is that difference equations such as (13) have any number of solutions.How do we deduce that we have the correct transform?
Let M (s) = M 0 (s)R(s).Substituting this form of M into (13) yields Any function R that satisfies the above equation can be used to get another solution to (13), but the solution must: (1) satisfy condition ( 14); (2) have an inverse transform; and (3) the inverse transfom must satisfy (5).It turns out that the above restrictions force R to be a constant function whose value is determined by ( 14), provided we restrict our attention to meromorphic transforms.The proof of this will be sketched.
Let M be a meromorphic function that satisfies (13) that has an inverse transform.Let H(Ω) denote the set of functions holomorphic in Ω ⊂ C and υ, Υ denote the strip {s : υ < Re(s) < Υ}.Condition (14) indicates that the strip of holomorphy for M must include the line Re(s) = 1 so that there are numbers υ < 1 and Υ > 1 such that M ∈ H( υ, Υ ).We will first show that there is a number υ 0 < 0 such that M ∈ H( υ 0 , Υ ).
Suppose that there is a υ, 0 ≤ υ < Υ such that M ∈ H( υ, Υ ) but M has at least one pole on the line Re(s) = υ.The poles to the left of the line Re(s) = 1 and closest to this line provide the leading order terms for the asymptotic expansion of y as x → 0 + .Suppose there is a pole of order k at s 0 = υ + iτ 0 with a principal part of the form Then it can be shown [17] that the contribution to the leading order term as x → 0 + is where M < υ.Condition (5), however, requires that lim x→0 + y(x) = 0. Suppose there is another singularity s 1 that cancels part or all of these leading terms for all x > 0 close to 0. Evidently such a pole must lie on the line Re(s) = υ, so if this is possible, s 0 is not on the real axis.Let s 1 = υ + iτ 1 .Suppose that this singularity cancels the contribution from say d j = 0 Then, there must be a number D j such that for all x in some non empty interval (0, ).Since τ 1 = τ 0 this is not possible.In this manner we see that M ∈ H( υ 0 , Υ ), where υ 0 < 0 and it follows immediately from ( 13) that M can be holomorphically continued in C. Any solution M must thus be an entire function.Now, so that R is holomorphic for all s except perhaps where M 0 (s) = 0. We know that the zeros of M 0 are of the form s = −m + 2kπi/ ln α where m is a non negative integer and k is a non-zero integer.In particular, R must be holomorphic in the half plane Π 0 = {s : Re(s) > 0}.Condition ( 16) can be used to holomorphically continue R to the other half plane so that R must also be an entire function.
Let σ ∈ R. The inverse transform of M is given by and the Riemann-Lebesgue lemma shows that as |τ | → ∞.We can use equation (17) to glean a growth condition on R.
Stirling's formula can be used to discern the growth of the gamma function.In particular, in any strip υ 1 , Υ 1 , where υ 1 > 0 and Υ 1 − υ 1 > 1 Stirling's formula indicates that, for s = σ + iτ , |Γ(s)| cannot decay faster that Ce − π 2 |s| as |τ | → ∞, where is C is some positive constant.In such a strip the function Q is bounded away from 0, and equation ( 17) thus shows that as |τ | → ∞.The periodicity of R shows that equation ( 18) is valid for all σ ∈ R and hence the growth condition is valid throughout C.
We now show that R must be a constant.To establish this we use the following result the proof of which can be found in [35], section 5.8: Theorem 3.2 (Carlson).Suppose f ∈ H(C) and f (s) = O(e k|s| ) as s → ∞ , where k < π, for Re(s) ≥ 0. If f (s) = 0 for s = 0, 1, 2, 3 . .., then f (s) = 0 identically.Let f (s) = R(s) − R(0).The function f is entire and satisfies the same growth condition as R; hence, the function f is of exponential type k for some k < π.Equation ( 16) implies that f (m) = 0 for all integers m and it follows by Carlson's theorem that R must be constant.
The above argument shows that the only meromorphic solution M that satisfies (14) and has an inverse transform that meets conditions (5) is M (s) = KM 0 (s), where .
The constant K can be evaluated as elliptic integrals or by use of the Euler pentagonal number theorem [21].It can also be expressed in terms of theta functions or the Dedekind eta function (cf.[29] for a full study).Note that at this stage, without knowing the inverse transform, the fact that M ∈ H(C) implies that the inverse transform y must have derivatives of all orders, that all the derivatives of y at x = 0 vanish (easy to show directly from the pantograph equation), and that the solution must decay faster than x −p for any positive number p.The inverse transform of M is straightforward to determine in this case owing to the Euler partition identity, which provides a means to convert the infinite product to an infinite series.The Euler identity is , which is valid for |q| < 1 and z ∈ C (cf. [3], pg.19).Let q = 1/α and z = −1/α s , then where .
The inverse transform is thus given by the Dirichlet series We stress that in this simple case all of the solution properties can be gleaned more efficiently by direct proofs from the pantograph equation (cf.[21]) coupled with a candidate solution that is known explicitly thanks to the Euler identity.The direct use of the Mellin transform equation, though clearly awkward in this simple case, nonetheless gives a glimmer of what types of tools are needed when we cannot invert the transform analytically.Although there may be simpler proofs that R must be constant, certainly the interplay between integrability requirements, the initial condition at x = 0 and the value distribution of R as a meromorphic function play a key rôle here.If, for instance, a more general problem allowed R to have poles, then tools from Nevanlinna theory would seem appropriate.

Some Generalizations
The pantograph equation has been generalized in many different directions.It has been studied in the complex plane [11].Connections have been made with complex dynamics [25,38], second order singular Sturm Liouville problems [37], and eigenfunction expansions even in the first order case [41,42].In terms of probability theory, the pantograph equation has been interpreted as a special case of a much broader class of equations [8].
Returning to the basic cell division model there have been at least three distinct generalizations apart from choices of variable growth and division rates such as found in [22,40].These generalizations are: (1) Systems of equations (at least one of which is the pantograph type) that model cell division taking into account the cell cycle; (2) Models of asymmetric cell division; and (3) The extension of the model to incorporate dispersion.The first generalization is motivated by tumour growth and problems relating to the response of cells to chemotherapy.The reader is referred to [4] as an example.The second generalization models cells of size x that divide into two daughter cells of size x/α and x/β, where α > 1 and β > 1 are numbers such that 1/α + 1/β = 1.The governing equation here is and the corresponding pantograph equation is (G(x)y(x)) + (B(x) + λ) y(x) = αB(αx)y(αx) + βB(βx)y(βx).
This equation was studied for the special case of constant growth and division rates: G(x) = g, B(x) = b, and it is shown that the solution is a double Dirichlet series of the form [36,45]).
The third generalization extends the basic first order model to include a dispersion term.Such a term arises, for instance, when the growth rate is stochastic.The governing equation in this case is where D is a non negative function.The boundary conditions are The corresponding pantograph equation has been studied for the simplest case when D, G, B are constant functions [43] and for certain other non constant coefficient cases where the solution can be found analytically [5,39].Of note here is the work of Begg et al. [6] directly on equation (19).The relative entropy methods developed by Mischler et al. [27,28,31] were extended in this work for a broad class of division rates (cf.also [32]).The results in [6] were targeted for particular rates B and G that reflect a tumour growth model, yet the authors set up machinery that can readily be adapted to more general problems.

Conclusions
The pantograph equation and its generalizations play a key rôle in the study of long time asymptotic solutions to the simple cell division/growth model.These equations appear naturally not only through the separable solutions to the differential equation but also in the "adjoint equation" used in the relative entropy methods, where the equation is retarded [31].Even the simplest versions of this equation produce analytical surprises: they are very sensitive to whether the argument is advanced or retarded, the division and growth rate coefficients, and the asymptotic requirements on solutions as x → 0 + and x → ∞.
In this paper we have surveyed the current rôle of the pantograph equation as it applies to the cell division model.We have also considered the problem of extracting solution properties directly from the Mellin transform equation.Although it may seem a Pyhrric victory for the simple case studied here, extracting solution information from a transform equation is illuminating because it illustrates the mathematical tools that must play a part in models where we do not have the luxury of an analytical inverse transform.
Although the pantograph equation plays an important part when G is constant, the situation becomes murky when G is variable.In particular, when G is a linear function, then there is no dominant eigenvalue [13], and for a broad class of division rates B we know that the long time solutions oscillate in time with asymptotic period ln α [7].Finally, there remain anomalous cases such as when the growth rate is linear and the division rate is constant.These and more general critical cases have been studied in [15].