Copyright is owned by the Author of the thesis. Permission is given for a copy to be downloaded by an individual for the purpose of research and private study only. The thesis may not be reproduced elsewhere without the permission of the Author. Steady Size Distributions in Cell Populations A thesis presented in partial fuJflJment of the requirements for the degree of Doctor of Philosophy in Mathematics at Massey University Alistair John Hall 1991 Abstract In any population of cells, individual cells grow for some period of time and then divide into two or more parts, called daughters. 1b describe this process mathematically, we need to specify functions describing the growth rate, size at division, and proportions into which each cell divides. In this thesis, it is assumed that the growth rate of a cell can be determined precisely from its size, but that both its size at division and the proportions into which it divides may be described stochastically, by probability density functions whose parameters are dependent on cell size and age (or birth-size). Special cases are also considered where all cells with the same birth-size divide at the same size, or where all cells divide exactly in half. We consider a population of cells growing and dividing steadily, such that the total cell population is increasing, but the proportion of cells in any size class remains constant. In Chapter 1, equations are derived which need to be solved in order to deduce the shape of the steady size distribution (or steady sizelage or sizelbirth-size distributions) from any given growth rate and probability distributions describing the division rate and division proportions. In the general case, a Fredholm-type integral equation is obtained, but if the probability of cell division depends on cell size only (i.e. not age or birth-size), and all cells divide into equal-sized daughters, then we obtain a functional differential equation. In two special cases, the resulting equations simplify considerably, and it is these cases which are explored further in this thesis. The first is where the probability of a cell dividing in any instant of time is a constant, independent of cell age or size. In Chapter 2, the functional differential equation resulting when cells divide into equal-sized daughters is solved for the special case where the growth rate is constant, and in an appendix the case where the growth rate is described by a power law is dealt with. The second case which simplifies is where the time-independent part of the growth rate of a cell is proportional to cell size. This case is particularly important, as it is a good first-order approximation to the real cell growth rate in some structured tissues, and in some bacteria. The special case in which this leads to a functional difl'erential equation is discussed in Chapter 3,-and the integral equation arising in the general case is dealt with in Chapter 4. Finally, the conditions under which the integral operator in Chapter 4 will be both square-integrable and non-factorable are discussed in Chapter 5. It is shown that if these conditions are satisfied then a unique, stable, steady size distribution will exist. i Acknowledgements My thanks to the New Zealand Department of Scientific and Industrial Research (DSIR) for allowing me to spend half of my time on this work while still on full salary. Without this support, firstly from Plant Physiology Division then latterly (thanks to the current enthusiasm for restructuring) from DSm Fruit and Trees, this work could not have been carried out. I would particularly like to thank the following for their help during the preparation of this thesis: • Prof. G.C.Wake (Massey University), who has taught me a great deal during the course of this work, including not only a great deal of mathematics but also an approach to applied mathematics problems in general. The material for much of the thesis has come from the application of one or both of two principles which seem to characterise his approach to difficult or apparently insoluble problems: , - Try something even if the chances of success look poor, or - Solve a simpler problem then look for ways to generalise! I have also very much appreciated his willingness, however busy, to take time to discuss my mathematical problems. • Dr. P.W.Gandar (DSm Fruit and '!rees), who provided the initial motivation and enthusiasm for this area of research, and who from a biological perspective has provided numerous ideas and insights which have assisted the mathematical development. • My wife Anne, son Michael, and daughter Joanne, who have with good grace put up with my odd hours of work and variable moods! Thanks are also due to Dr. M.R.Carter (Massey University) for a number of helpful suggestions, including the expansion used in (2.49); Dr. J.F.Harper (Victoria University of Wellington) for suggesting the transformation used in (2.14); An unknown referee for a number of constructive suggestions in Chapter 2, including the use of the Mellin transform in section 2.2.4 and the connection with Jacobian elliptic functions in section 2.4.1 .; and Prof. O.Diekmann (Centre for Mathematics and Computer Science, Amsterdam) who as a referee for the Journal of Mathematical Biology made a number ofhelpfu} corrections and suggestions on Chapter 4. ii Contents 1 Introduction and Model Formulation 1 1 .1 Unstructured Population Models . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 .2 Age-Structured Populations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1 .3 Populations Structured on Many Properties . . . . . . . . . . . . . . . . . . . 4 1 .3.1 Models with Age. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1 .3.2 Models without Age. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1 .4 Population Models with Non-Deterministic Growth Rates . . . . . . . . . . . 8 1 .5 Cell Population Models with Size Structure Only . . . . . . . . . . . . . . . . 11 1 .6 Steady Size Distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1 .7 Cell Populations with age and size structure . . . . . . . . . . . . . . . . . . . 17 1 .8 Steady Size/Age Distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.9 The Simplest Case: A Constant Probability of Division . . . . . . . . . . . . . 21 1 .10 A Special Case: Cells in One-Dimensional Structured Tissues. . . . . . . . . 22 2 The Equationy'(x) = -ay(x) + aay(tu) ill 27 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.2 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.2.1 y(O)andy(oo) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.2.2 Restrictions ony. «. and a . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.2.3 Uniqueness of the Solution . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.2.4 Moments of the SSD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.3 Solution of the SSD Equation . . . . . . . . . . . . . . . . . . . . . . . . . . " 33 2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . " 35 2.4.1 'I'he Constant K . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 35 2.4.2 Properties of the Solution . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.4.3 An Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.4.4 'I'he Solution for Extreme Values of a . . . . . . . . . . . . . . . . . " 38 2.4.5 'I'he Form of the Time-dependent Number Density n(x, t) . . . . . . " 42 3 Size-Structured Models for Cells Growing Exponentially 43 3.1 In.troduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 43 3.2 PreIimjnsries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.3 'I'he case b(x) = bxt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.3.1 'I'he Case k :F m-I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.3.2 'I'he Case k = m-I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.3.3 Discussion and Examples . . . . . . . . . . . . . . . . . . . . . . . . . . 52 iv 3.4 A Minimum Size for Division: The Case b(x) = bxkH(x - Xl) • • • • • • • • • • 56 3.4.1 Region 1 , X � Xl • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • , 56 3.4.2 Region 2, � � X � Xl • • • • • • • • • • • • • • • • • • • • • • • • • • • . , 57 3.4.3 Region 3, X � � • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 58 3.4.4 Discussion and Examples . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.5 A Minimum and Maximum Size for Cell Division . . . . . . . . . . . . . . . . 59 3.5.1 The Solution Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.5.2 An example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4 Cell Populations with Age and Size Structure 64 4.1 Part I: An Integral Equation for the Birth-size Distribution . . . . . . . . . . 65 4.1.1 Introduction: Size/Age and SizelBirth-size Distributions . . . . . . . . 65 4.1 .2 Steady Size/birth-size Distributions . . . . . . . . . . . . . . . . . . . . 67 4.1 .3 Growth Rate Proportional to Cell Size . . . . . . . . . . . . . . . . . .. 69 4.1 .4 Size, size/age, and birth-size distributions. . . . . . . . . . . . . . . . . 71 4.1 .5 A Fredholm Integral Equation for ;(s) and its Unique Solution . . . . 73 4.1 .6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.2 Part 11: Solving the Integral Equation . . . . . . . . . . . . . . . . . . . . . . . 77 4.2.1 Constant cell division size . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.2.2 Hazard Rate Independent of Initial Cell Size . . . . . . . . . . . . . . . 79 4.2.3 Cell Division Size Precisely Determined by Birth-size . . . . . . . . . 83 v 4.2.4 Numerical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5 Existence, Uniqueness, and Stability of SSD. 5.1 Existence and Uniqueness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.1.1 Square-integrability of the Kernel . . . . . . . . . . . . . . . . . . . . . 96 5.1 .2 Non-factorability of the Kernel . . . . . . . . . . . . . . . . . . . . . . . 97 5.2 The Time-dependent Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.2.1 Formulation and Simplification . . . . . . . . . . . . . . . . . . . . . . . 99 5.2.2 Solution using Characteristics . . . . . . . . . . . . . . . . . . . . . . . 101 5.2.3 The Region Influenced Directly by the Initial Condition . . . . . . . . 101 5.2.4 The Region Influenced by the Boundary Condition . . . . . . . . . . . . 102 5.3 Stability of the SizelBirtb-size Distributions . . . . . . . . . . . . . . . . . . . 104 5.3.1 The Laplace Transform of the Operator . . . . . . . . . . . . . . . . . . 104 5.3.2 Eigenvalues of the Operator K(,t) . . . . . . . . . . . . . . . . . . . . . . 106 5.3.3 The Inverse Laplace Transform . . . . . . . . . . . . . . . . . . . . . . . 108 5.4 Examples and Counter-examples . . . . . . . . . . . . . . . . . . . . . . . . . . 112 5.4.1 Some Stable Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 5.4.2 Some Unstable Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 A Moments of a Limiting Distribution 120 vi A.l Formulation as a Combinatoric Limit . . . . . . . . . . . . . . . . . . . . . . . 120 A.2 Formulation Involving Polynomial Coefficients . . . . . . . . . . . . . . . . . . 122 A.3 An expression for rP1f(P) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 A.4 Proof of the Main Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 A.5 Proof By Induction of the Expression for Sf . . . . . . . . . . . . . . . . . . . . 127 B The case b(x) = b with g(x) = gxk 180 C The General Equationy'(x) = -p(x»'(x) + q(x)y(ax) 182 D Generalised Solutions of the Functional Differential Equation 186 Bibliography 189 vii Chapter! Introduction and Model ForlDulation 1.1 Unstructured Population Models Historically, the earliest population models developed were unstructured - they ignored the properties of individuals and considered only the total number of individuals in the population. In the eighteenth century, Malthus [35] analysed what we now term the human "population explosion", using the unstructured model of populations which we now write as dN �(t) == r}{(t), (1 .1 ) where }{(t) is the total population at time t. The population's intrinsic growth rate or Malthusian parameter, r, is given by r == f3 - p, (1.2) where f3 is the population's birth rate per unit time (per head of population), and p is the corresponding death rate. If r is independent of time, and r > 0, then the population will grow exponentially for all time. Note that in this model, and throughout this thesis, it is assumed that there is no immigration or emigration of individuals. An extension of this idea is to allow the Malthusian parameter, r == P - p, to vary with time. For example, r may be assumed to depend on the environment in some way, with 1 the population growing at a given instant if r > 0 and contracting if r < O. It is possible to look back at historical records and deduce the Malthusian parameter for any given periods in the history of a population - not surprisingly, this turns out to be positive during the "good times", and negative during wars and famines! For a more predictive model, it is often assumed that r is directly dependent on the total population N(t), and therefore only indirectly dependent on t. For example, we may decide that r should decrease as N(t) increases (possibly because of crowding effects or competition for resources), in which case a differential equation of the form tiN = 1()(1 _ ( N(t) )4)N(t) dt Nmax (1 .3) may be appropriate, with A. > 0 so the intrinsic growth rate is negative if N(t) > Nmax• The intrinsic growth rate for small N(t) is just 1(). If A. = 1 , then we have the usual form of the logistic equation first discussed by Verhulst [65] in 1838. The solution of (1 .3) with A. = 1 , in conjunction with the initial condition N(O) = No, is where NllJIAX N(t) = 1 -Ae7Ol' Nmax A=l - - . No (1.4) (1 .5) The various other extensions of the logistic and alternative unstructured continuous-time population models that have appeared during the last two centuries can generally also be regarded as variations on (1 .1) in which r is dependent either on time t explicitly or on the total population N(t). One alternative which should be mentioned here because of the application it has found in biology generally, and in the study of tumor growth particularly, (Laird [32]; Simpson-Herren and Uoyd [51 ]), is the Gompertz growth law, proposed by Gompertz [1 2] in 1825. We can regard this model as one where the Malthusian parameter is proportional to the log of the ratio of the current population N(t) to some maximum population size Nmza., so ':: = (-aIOg(:� )) N(t). If we once again let No be the population at time t = 0, then the solution is (N )[1-'-4 N(t) =No ;: . (1 .6) (1 .7) Unstructured population models, based on total population alone, are useful where the individuals in a population are indistinguishable or where the observed differences do not 2 appear to influence the processes of birth and death. However, this is not the case for the cell populations discussed later in this thesis, nor is it the case for human populations. Because of the interest in human population growth, where the death rates and birth rates are most obviously dependent on the age of individuals, the first structured population models to appear were those structured on age. 1.2 Age-Structured PopulatioDS The obvious deficiency of unstructured population models is that they ignore the differences between individuals. For example, in animal and human populations it is clearly a better approximation to make the birth rate fJ dependent on the number of females of child­ bearing age rather than on the total population. In order to do this, we need a model structured on age at least (the sex structure is often allowed for by ignoring the males altogether, and just modelling the female population!). Early work on age-structured populations was carried out by Euler in the 18th century (for an English translation see Smith and Keyfitz [54]), and in the early twentieth century by Lotka [34] and Sharpe and Lotka [50], with some mathematical rigour added in 1941 by Feller [11 ]. The problem was first formulated as a partial differential equation by M'Kendrick [36] in 1926, but it was not until the work ofScherbaum and Rasch [49] in 1957 and von Foerster [66] in 1959 that this approach was taken up again by mathematical demographers. In this approach, we consider a continuous number density n(a,t), such that at any time t, 1:: n(a,t)da is the number of individuals with age a lying between a1 and a2. Note that this approach requires that the numbers involved are in some sense "large", 80 that the approximation of the population by a continuous number density is reasonable. We assume that the proportions of individuals in any age range dying or giving birth are independent of time. It is then possible to define a death rate pea) and a birth rate (or "procreation rate") pea) such that p(a)dt and p(a)dt give the proportions of individuals of age a which die or give birth respectively in a short time period tit. Then for t,a > 0 we have a partial differential equation, lJn on &t + oa = -p(a)n(a, t), (1 .8) which can be regarded as a statement of number conservation. Equation (1 .8) is known, somewhat unjustly, as "von Foerster's equation". The initial condition at time t = 0 is given by n(a, 0) = 1IO(a), 3 (1 .9) and the boundary condition or "renewal equation" at age a = 0 is given by n(O,t) = looo p(a)n(a, t)da. (1.10) The book by Metz and Diekmann [37] discusses a wide range of structured population models, and has one chapter devoted exclusively to age-structured models. Age-structured models and the von Foerster equation are also discussed in Oster [39], and applications to cell populations are discussed by Trucco ([57] , [58] ) and Rubinow [46] . A more general approach would be to allow Jl and p to depend on t as well (which caters for possible environmental changes), or on the total population N(t) (Gurtin and Macamy [15]). 1.3 Populations Structured on Many Properties Instead of considering a population structured on just one variable (such as age in the previous section), consider a population where the individuals are distinguished by say m properties, ! = {Xl ... X".}. One of these properties may be age, as in the previous section, and others may include such things as mass, height or chemical content. Each individual is born with a set of values for these properties, which may depend in some way on the properties of the parent, then these properties change as time goes on. In this thesis, attention is restricted to continuous properties of individuals only, by which we mean that between birth and death the properties of an individual change continuously. Following the ideas of the previous section, we define a continuous number density ne!, t) such that J ... ne!, t}4.x, tDab = {!:;a � ! � !b}, gives the number of individuals with properties in the range!a to !b. We will call the region in which it is possible for n to be non-zero o. Then given that there is no immigration or emigration, we can write a simple number balance at! E Cl symbolically as a . al(!' t) + v 1! = (birth rate) - (death rate), (1.11) where I! is the vector flux of individuals at ! in the m-dimensional property space. The birth and death rates are defined more clearly below. One way of characterising continuous structured population models is according to whether the processes are deterministic, in the sense that they are completely determined by the properties of the individUal! and the time t, or whether they are in some sense stochastic. Processes we may wish to consider in this light are: 4 • Growth rates, or rates of change of each property Xj. Section 1 .4 of this chapter looks briefly at the case where the growth of individuals is described as a stochastic process, but in the rest of this thesis it is assumed that the growth rate of the ith property of an individual is a determined entirely by the individual's properties; and the time I, so we may write the growth rate of the ith property as dxj di = gi(;, I). (1.12) • Birth Rates. Whether or not an adult will give birth will in general be dependent on its properties ; and time I. We will allow birth to be a stochastic process, by defining a birth rate P(;, I) such that P(;, 1)tIt gives the probability that an adult with properties; at time I will give birth in the next short time interval dt. We show later that the case of deterministic birth of individuals can be regarded as a special case of this general case. • Death Rates. We also treat death as a stochastic process, defining a death rate Jl(;, I) such that Jl(;, l)dt gives the probability that an adult with properties; at tune I will die in the next short time interval dl. Once again, deterministic death can be regarded as a special case of this general case. • Heredity. Each new-born individual must be given a set of properties to begin with. If heredity were to be modelled deterministically, then the properties of the newborn individual would be precisely given by some function of the properties of the parent. Instead, we describe heredity in a stochastic manner, defining a heredity function 71(;,�,I) such that 71(;,�,I)4t gives the expected number of offspring with properties in the range ; to ; + 4t produced by an individual with properties � giving birth at time I. This heredity function 71 is non-local in that it describes the properties of an individual at one point ; in the property space in terms of the properties of another point�. Note that this formulation also means that the number of newborn individuals produced when an adult gives birth need not be deterministic. We need to consider two possibilities with respect to newborn individuals: • Newborn individuals may appear on a boundary of the region Q. Usually this means that one of the variables is age, a, as in the previous section, and newborn individuals appear as a flux over the boundary at a = o. In Chapter 4 we use a slightly different formulation in terms of birth-size instead of age, where newborn individuals appear on a boundary of a different type. 5 • Newborn individuals may appear throughout the region n. Usually this occurs when there is no age-like variable in the structured population model. 1b simplify the following discussion we assume that all properties Xi are restricted to the region [0,00), so n = [0,00) x ... x [0,00). 1.3.1 Models with Age Consider the case where all our functions are dependent on age a, time t, and m other variables -! = (Xl, . . . Xm). Hence we write n(-!, a, t), gi(-!,a, t), p(-!,a, t), Jl(-!,a, t), and 77(-!,�, a, t). In this case, in the region n the vector flux I!. is given by I!. = �(-!, a, t)n(-!, a, t), so the number conservation law (1.11 ) is written as the partial differential equation an an � a -a + -a + � o: (gi(-!,a, t)n) = -Jl(-!,a, t)n. t a i=l u'X, (1 .13) This equation is discussed in Oster and Takahashi [40], and a detailed derivation for the special case of populations structured on age a and one other variable X is given in Sinko and Streifer [52]. For a brief derivation of the general case, see Chapter 3 ofRandolph and Larson [43]. Note that because birth occurs on the boundary of n, there is no source term on the right hand side of (1 .13), and that the term containing the partial derivative with respect to a is particularly simple because � = 1 . The initial condition corresponding to this PDE should be given in the form n(-!,a, O) = no(-!, a). The boundary condition (or "renewal equation") at a = ° may be written as n(-!, O, t) = kP(�, a, t)77(-!,�,a, t)n(�, a,t')d§da. (1.14) (1 .15) Fig. 1 .1 shows why this is so diagramatically for m = 2 - the number of individuals born in time dl with property 1 in the range Xl to Xl + dxl and property 2 in the range X2 to X2 + dx2 is simply the number in the lower box (which has height dl because � = 1), and this must equal the sum of the numbers of offspring from the individuals in all possible boxes like the upper one. If a boundary condition is needed at any of the other boundaries of n, we simply use the zero flux conditions gi(-!, a, t)n(-!, a, t)\x;=o = gi(-!, a, t)n(-!, a, t)\x;=oo = 0, i = 1, ... ,m. 6 (1 .16) a Figure 1 .1 : The number of new-born individuals appearing over the boundary a = 0 in time tit with properties in the range (Xl ,X2) to (Xl + dxl , X2 + dx2) must equal the sum of all contributions from cells dividing in the appropriate proportions in the same time tit. For large age, we use n(�,oo,t) = O. (1 .17) Note that in some population models, a physwlogical or developmental age is used instead of chronological age. Under this formulation, individuals are still all born with age zero, but some small changes need to be made to (1 .13) and (1 .15) because � 1: 1 for physiological or developmental age. Van Sickle [64], Rubinow [47] and Sundareshan and Fundakowski [56] discuss models of this type, while Weiss [68] allows stochastic variation in the rate of change of physiological age. L3.2 Models without Age If none of the relevant properties of individuals are age-like variables, then a source term is required on the right-hand side of the PDE corresponding to (1 .13), but no boundary condition is necessary. The PDE becomes (1 .18) The last term in (1 .18) comes from equating the number of newborn cells formed in time tin some small m-cube in m-dimensional property space to the sum of contributions from cells 7 ·····;;;··10 •...•... 1 - 'ch ' , , , I' , , , , , , , , tU2g /' Figure 1 .2: The number of new·born individuals appearing in some short time interval lit, with properties in the the range (Xl,X2) to (Xl + dxl,X2 + dx2), must equal the sum of all contributions from cells with properties (st ,$2) dividing in appropriate proportions in the same time dt. anywhere in the space dividing in appropriate proportions. Fig. 1 .2 shows the situation for m = 2, where the number of newborn cells formed in a short time lit with property 1 in the range Xl to Xl + dxl and property 2 in the range X2 + dx2 is equal to the sum of contributions from cells with any properties (st ,$2) dividing in time lit. Alternatively, (1 .18) may be obtained by considering a model with age incorporated, but where the functions jJ, p, f, and gi, i = 1 , ... ,m, are in fact independent of a. We can then integrate (1 .13) over all age a, so that :,<1000 nda) + n(�, oo, t) - n(�, O, t) + i� �/gi(�, t) 1000 nda) = -jJ(�,t) 1000 nda. (1 .19) Now we use (1 .15) and (1 .17) to replace n(�, O, t) and n(�, oo, t) respectively, then replace fo n(�,a,t)da by n(�, t) to obtain (1 .18). Once again, if boundary conditions are required they are given by the zero-flux conditions (1 .20) 1.4 Population Models with Non-Deterministic Growth Rates This thesis is concerned with steady size distributions arising in models of cell populations where birth, death, and heredity processes are stochastic, but growth rates are determine istic. However, it may be that there is experimental evidence showing significant variation 8 in the growth rates of individuals all with the same measured properties, and in this case the deterministic growth models discussed here would be inappropriate. In this section we take a brieflook at the situation where the growth rate of individuals may be regarded as a stochastic process. The changes in a property are assumed to be local in the sense that the change in some time interval tu must tend to zero as III - o. Consider the growth rate of just one property of individuals, x say, which we will call size. 'Ib simplify notation we will write just II(X, t). Then instead of writing the growth rate of x as g(x), we choose some short time period tu and for each x and t define a probability distribution r(x, t, e) such that r(x, t, e)de gives the probability that an individual with size x at time t will increase size x by between e and e + de in the next tu. Then for any x and t, ignoring birth and death terms for the moment, we have II(X,t + Ill) = 1: II(X - e,t)r(x - e,t,e)de. (1 .21) Expanding the expression inside the integral as a partial Taylor series about II(X, t)r(x, t, e) gives II(X,t + tu) = 1: [1I(X,t)r(x,t,e) - e !(II(x,t)r(x,t,e» + � ::2 (II(X,t)r(x,t, e» - .. .J de (1 .22) We now assume that the probability density function r is sharply peaked about its mean for each x and t, so that we may ignore terms of order higher than two. Changing the order of differentiation and integration then gives n(x,t + tu) = n(x,t) 1: r(x,t,e)de -! [ n(x,t) L: £1(x,t,e) 4£] + �:; [1I(X,t) 1: �r(x,t,e) de] (1 .23) Now to simplify notation we make the following substitutions, using EO to indicate expected values for fixed x and t: L: r(e)de = 1 (1 .24) f: £1(e)de = E[e} (1 .25) L: �r(e)de = E[�] (1 .26) If we now include mortality and birth as in the deterministic case, writing J.l(x, t), fJ(x, t), and 11(X, s, t), (1 .23) becomes o 1 � [ ] II(X,t + M) = n(x,t) - Ox [(n(x,t)E[e]] + 2 {Jx2 II(X,t)E[�] 9 so that - p(x, t)n(x, t)&I+ k P(S, t)T1(X, S, t)n(s, t)ds&l, (1 .27) n(x, t + &I) - n(x, t) _ &I - a [ E[E] ] 1 [j2 [ E[e2]j - Ox (n(x, t) &I + '2 Br n(x, t)--:it - p(x, t)n(x, t) + k P(S, t)T1(X, s, t)n(S, t)ds. (1 .28) Now we let &I - 0 and define so g(x, t) D(x, t) an a 1 Q2 at = - Bx [g(x, t)n(x, t)] + '2 or [D(x, t)n(x, t)] -p(x, t)n(x, t) + kP(S, t)T1(X,s, t)n(s, t)ds. (1 .29) (1 .30) (1 .31 ) This is in the form of a Fokker-Planck equation (Cox and Miller [8]), where D(x, t) is called the dispersion coefficient, and g(x, t) approximates the mean growth rate. This second order partial differential equation (1 .31) incorporates stochastic birth and death, as in (1 .18), as well as stochastic growth rates of individuals. The only difference between this model and the deterministic models we consider elsewhere in this thesis is the term involving the dispersion coefficient D(x, t). For an alternative but similar approach to deriving (1 .31 ) in the context of structured populations, using the limit as &I - 0 of short finite transitions, see Ricciardi [44] or Weiss [68]. The book by Cox and Miller [8] diBCUBseS both this approach and an approach utilising Wiener processes in a more general context, and also contains a wider di8CU88ion of stochastic processes with the property that both the mean and variance are proportional to &I, so the definitions (1 .29) and (1 .30 above make sense. For a more general model A involving non-deterministic growth and transition probabilities, and another approach to deriving the Fokker-Planck equation, see Rotenberg [45]. In a "dispersion-like" model as described above, the growth rate of an individual is stochastic in nature and the trajectories of individuals with the same initial x value would be expected to cross and re-cross in x - t space (Fig. 1.3(a». Such a model would be inappropriate if the tnijectories of individuals were more like those shown in Fig. 1.3(b). A more appropriate model in this case might be to hypothesise some extra property of 10 t t Figure 1.3: Sample tnijectories of individuals with size Xo at time t = O. In case (a) the "dispersion-like" model given by (1.31) may be appropriate, whereas in (b) it would not. individuals, previously unused, which affects the growth rate of x, and include this as an extra dimension in a deterministic model. However, additional information or hypotheses would then be required to describe how this new variable changes with time, and how it is passed onto newborn individuals, so this approach may not be feasible in many situations. Kirkpatrick [29] makes the same distinction, calling the case shown in Fig. 1.3(b) the assignment at birth growth model, and that shown in Fig. 1 .3(a) the variable growth model. 1.5 Cell Population Models with Size Structure Only In this thesis we will be concerned with populations of cells which reproduce by fission into two or more parts. During fission the original or parent cell disappears and is replaced by two or more newborn cells, called daughters. We will call a property of a cell a size if it is conserved during cell division, in the sense that the sum of the sizes of the new-born daughter cells must equal the size of the original parent cell. Properties which could be called the size of a cell include mass, volume, quantity of DNA, or in some cases lengtb.. A size may only take on positive values. We will restrict our attention to cell populations structured on a single size variable I only, with a strictly positive growth rate which we write as gl(/,t) > 0, VI. Our choice ofthe symbol I for the size variable here is motivated by the work in Chapter 4, where cell length is the size 11 variable chosen. We now write (1.18) as 8 8 tX> 8tn(/, t) + 8lg1(/, t)n(/, t) = -p(/, t)n(/, t) + 10 P(S, t)71(/, s, t)n(s, t)ds. (1.32) S ize-structured models of cell populations were given prominence by a number of re­ searchers in the 1960's, including Collins and Richmond [7], Koch and Schaechter [30], Powell [42] and Painter and Marr [41]. More recently, there has been the mini-review by Tyson [6P], the papers by Koppes et al. [31] and Grover et al. [13], papers by 1Yson and DiekmsDD [61] and 1Yson and Hannsgen [62], and the book by Metz and Diekmann [37] which contains a good discussion of the special case of (1.32) where all cells divide exactly in half. When a cell divides, it generally does not split into daughters all of the same size, although this is clearly one limiting case which should be considered (Metz and DiekmsDD [37]). Experimentally, the most important quantity in predicting the size of a daughter cell from the size of the parent is the ratio of the daughter size to the parent size. This suggests that instead of using the heredity function 71(I, s, t) for cells reproducing by fission, we should define a probability density function/(P,I, t) such thatl(p, I, t)dp gives the probability that a particular daughter lies between fraction p and p + dp of the parent cell size, given that the cell was derived from a parent which divided at time t and size I. Now because p is the ratio of the size of a daughter cell to its parent, and we do not wish to allow negative cell sizes or daughters larger than the parent, we demand that the support of I as a function of p be some subset of [0, 1] for all I and t. Further, I is a probability distribution so for all I and t we must havel(P, I, t) � 0, Vp and also fo1 /(P, I, t)dp = 1. (1.33) Finally, conservation of size properties requires that if the mean number of cells formed when an adult cell of size I divides at time t is a(/, t), then we must have 101 pl(p, I, t)dp = [a(/, t)r1 , (1.34) which requires a(/, t) > 1. Note that it is often assumed that cells divide into exactly two parts, in which case a is given the value 2 for alII and t. We note that 71(I,s,t)dlshould give the average number of daughter cells with size lto I +dl produced by a parent cell of size s at time t, which is terms ofl is I dl q(l,s, t)dl = a(s, t)f(- ,s, t)-. (1.35) s s 12 For simplicity, and because we are interested primarily in populations of cells which do not disappear or die by any mechanism other than fission, we assume that fission is the only way in which cells "die-. Hence the two separate functions P and p. as in Section 1.3 are unnecessary, and we simply define a function b(l, t) such that b(/, t)dt gives the probability that a cell of size I at time t will divide in the next short time interval dt. Using (1.35) and replacing P(l, t) and p.(l, t) by b(/, t), we write (1.32) as a a tJO I ds a,n(/, t) + 8lg1(/, t)n(/, t) = -b( I, t)n(/, t) + 10 b(s, t)a(s, ty"<-;,s, t)n(s, t)s· (1.36) The equation given by Sinko and Streifer [53] for the special case where cells always divide into two pieces whose sizes are in some fixed ratio, can be obtained by substitution of the appropriate combination of Dirac delta functions into (1.36). Equation (1.36) is the most general form of the equation for the number density n(/, t) in size-structured populations with deterministic growth and no cell death. However, noting as earlier that heredity is primarily dependent on the ratio of daughter to parent cell size, we assume throughout the remainder of this thesis that: • The fractional-size pdf,/, is dependent only on the proportionp, so we may wrltef(P) instead of/(P,I, t). It then follows from (1.34) that a is a constant. The PDE for a size-structured cell population with deterministic growth and stochastic division is thus a a lOO I ds a,n(l, t) + 8lg1(l, t)n(l, t) = -b(/, t)n(l, t) + a 10 b(s ,t)f(-;)n(s, t)s' which must be solved subject to an initial condition n(/, 0) = no(/) and the zero-flux boundary condition gIO , t)n(O, t) = #oo, t)n(oo, t) = o. 1.6 Steady Size Distributions (1.37) (1.38) (1.39) If we now define a hazard rate hi such that hl(l, t)dl gives the probability that a cell which reaches size I at time t will divide before reaching size I + dl. Then hl(/, t)dl = b(l, t)dt, (1.40) 13 so where we know g(l, t) > O. b(l, t) hl(/, t) = g(l, t) , (1.41) Now we are interested in steady size distributions (S S Ds) which can arise where the total number of cells in the population is increasing, but the proportion of the population in any given size range remains constant. This idea is £emiJiar from the demographic literature, where distributions of such quantities as age and height in a population are regularly tabulated and graphed, and the shape of "steady" age distributions, which develop in the presence of time-independent birth and death rates, are studied (Keyfitz [28]). Steady size distributions have been studied extensively in the bacteriological literature (e.g. Collins and Richmond [7], Lasota and Mackey [33], Koch and Schaechter [30], Tyson and Hannsgen [63]). Fig. (1.4) is an example of a experimentally measured S S D, taken from Harvey et al [22]. .6 .I 1 U lA Cell Volume epI) Figure 1.4: An experimentally measured steady size distribution (S S D) for Escherichia coli bacteria, taken from Harvey et al [22]. 'lb find the equations governing the form of S S Ds, we look for separable solutions of (1.37) such that n(l, t) = N(t)y(/), (1.42) where N(t) is the total cell population at time t, and y(/) is a probability density function with 1000 y(/)dl = 1, 14 (1.43) and y(/) � 0, VI. In order to find such separable solutions, we assume: • The hazard rate is time-independent, so we may write h,(/) instead of hl(l, t) • The growth rate is separable with respect to time, 80 we may replace g,(l, t) by g(l)r(t) say These assumptions mean we restrict our attention to situations where environmental changes may increase or decrease the growth rate of all cells by the same proportion r(t), but these changes have no effect on the size at which cells divide. Given these assumptions, and making the substitutions b(/, t) = h,(/)g(l)r(t) from (1.41) and n(l, t) = N(t)y(l) from (1.42), (1.37) becomes dN d [00 I ds dty(/) + r(t)N(t) dlg(I)y(/) = -r(t)N(t}h,(/)g(I)y(I) + r(t)N(t)a 10 h,(s)g(s)f(;)Y(s)s' (1.44) Separation of variables in this equation leads to � _ -�g(l)y(l) - h,(l)g(l)y(l) + a 10 hl(S)g(Sy(�)Y(s)� = Q (1.45) r(t)N(t) - y(l) - . An expression for the unknown Q may be obtained from the right hand equality here, d [00 I ds Qy(l) = - d/g(l)y(l» - h,(l)g(l)y(l) + a 10 h,(s)g(s)f(;)y(s)s' (1.46) We now integrate over all I and apply the zero flux boundary condition which from (1.39) is g(O)y(O) = g( oo)y( 00) = 0 (1.47) and the normalising condition (1.43) to obtain Q = (a - 1) looo h,( l)g( l)y( l)dl. (1.48) Also from (1.45), we see that the growth rate of the population as a whole is dependent on this functional Q, with dN � = Qr(t)N(t). (1.49) In the simplest case, where the growth rate of individuals is independent of time, we may set r(t) = 1 so Q becomes the exponential rate of growth of the population. 15 An alternative expression for Q may be obtained by considering the rate of increase of the total size of all the individuals in the population. This can be found either by integrating first over all sizes then differentiating with respect to t, or simply by integrating the growth rate over all individuals in the population. Equating these gives ! [1000 W(t)Y(l)dl] = 1000 N(t)y(l)g(l)dl, or ': [1000 ly(l)dl] = r(t)N(t) 1000 y(l)g(l)dl, then using (1 .49) we obtain Q = fooo g(/)y(l)dl fooo Iy( l)dl , (1 .50) (1 .51) (1 .52) so Q may also be regarded as ratio of the average growth rate of the individuals in the population to their average size. Rearranging (1 .46), the PDE for the SSD y(1) becomes d lOO I ds &g(/)y(/) = -(hl(/)g(1) + Q)y(/) + a hl(S)g(S)f( - )y(s)- , o s s (1 .53) which must be solved subject to the normalising condition (1 .43) and the boundary condition (1.47). Note that (1 .48) was deduced by integrating (1 .46) and applying the normalisation (1 .43) and the boundary condition (1 .47), so (1 .48) should be regarded as an alternative normalising condition rather than an independent equation. A non-zero normalisation condition of some sort is clearly necessary in order to obtain a non-trivial solution, because (1 .53) is homogeneous and therefore admits the solution y(/) == O. One special case is worth mentioning here, as it is used in both Chapter 2 and Chapter 3 and has been considered extensively in the literature (e.g. Sinko and Streifer [53], Koch and Schaechter [30]). If all cells divide into exactly a pieces each of the same size, we can write f(P) = 6(p - a-I ), (1 .54) where 6(p - a-I ) is a translated Dirac delta function. We thus obtain the functional differential equation d &g(/)y(/) = -(hl(l)g(/) + Q)y(l) + �hl(al)g(al)y(al). (1.55) A more direct derivation of (1 .55) is discussed in the book by Metz and Diekmann [37]. 16 In Chapter 2, we solve this equation for the simplest possible case, where both g and hi are constants, and in Chapter 3 the case where g(/) = gl is considered along with various forms of hi (I). In Appendix D, a general method of solution is presented for equations of the form y'(x) = -p(x)y(x) + q(x)y(ax), which could be used to solve (1 .55) in any situation where Q is known. 1.7 Cell Populations with age and size structure The case where an age variable is included along with a single size variable requires a similar set of definitions to the size-only case: • n(/, a, t) is the population number density, with n(1, a, t)dlda giving the number of cells in the size range I to I + dl and age range a to a + da at time t. • g(/, a, t) is the growth rate ofa cell of size I and age a at time t. • f(P) is a pdfsuch thatf(P)dp gives the probability that a particular daughter cell lies between fraction p and fraction p + dp of the parent cell size. • a, the average number of newborn cells formed when an adult divides, is then given by a = (fJ pf(P )dpJ-l . • bC/, a, t) is the probability per unit time that an adult divides, 80 that b(/, a, t)dt gives the probability that an adult cell of size I and age a at time t will divide in the next short time interval dt. Using this notation, the pde for the number denSity. n(/, a, t) becomes (from (1 .13» () () () ()tn(/,a, t) + ()a n(/, a, t) + mgl(l, a, t)n(/, a, t) = -b(/, a, t)n(/, a, t). (1 .56) This equation is discussed by Heijmans [24], as well as in a review by Oster [38], and a detailed derivation is found in Webb [67]. The initial condition corresponding to (1 .56) is given by n(/, a, O) = 1I()(l, a), (1 .57) and the zero flux boundary condition is gl(O, a, t)n(O, a, t) = gl(oo,a, t)n(oo, a, t) = 0. (1 .58) 17 'lb obtain a renewal condition at age a = ° corresponding to (1 .35) we substitute TI(I, s, a, t) = !If(�) into (1 .15), giving 1000 1000 I ds n(I, O , t) = a b(s, a, t)f(-)n(s, a, t)-da. o 0 s s (1 .59) Streifer [55] and Sinko and Streifer [52] contain good general discussions of size/age models, including forms of equations (1 .56) and (1 .59). Gyllenberg and Webb [17] discuss an elaboration of this model in which the cells are divided into populations of "quiescent" and "normal" cells. 1.8 Steady Size/Age Distributions Corresponding to SSDs in the case of size-structured populations, steady size I age distri­ butwns can arise in a population of cells with size and age structure, where the proportion of cells in any given size and age range remains steady while the total population grows. This idea was first given prominence in a series of papers by Anderson, Bell and others in the late 1960's (Bell and Anderson [6], Anderson and Peterson [4], Bell [5], and Anderson et al [3]). Following the approach in section (1 .6), we look for separable solutions of (1 .56) such that n(l, a, t) = N(t)y(l, a), (1 .60) and we define the following: • N(t) is the total population at time t. • y(l, a) is the steady size/age distribution such that y(l, a)dlda gives the proportion of cells in the size range I to I + dl and the age range a to a + da, 80 1000 1000 y(l, a)dlda = 1 . (1 .61) • r(t) and g(l, a) are, respectively, the time-varying and time-independent components of the separable growth rate gl(l, a, t), such that gl(l,a, t) = g(l, a)r(t). (1 .62) • h,(l, a) is the time-independent hazard rate such that h/(I, a)dl gives the probability that a parent cell of size I and age a will divide before reaching size 1 + dl, 80 h (I ) b(l, a, t) ( I , a = g(l,a)r(t) ' 1 .63) 18 Our assumption that hi is independent of time is equivalent to assuming that b is proportional to r( I). Note that these definitions assume the same separability and time-independence properties used in Section 1 .6. Substituting (1 .60), (1 .62), and (1 .63) into (1 .56), we obtain dN 8y f) (jiy(/, a) + N(/) f)a (/,a) + r(t)N(t) 8lg(/, a)y(/, a) = -r(t)N(t)hl(/,a)g(/, a)y(/,a) . (1 .64) The zero-flux boundary condition (1 .58) becomes g(O,a)y(O,a) = g(oo, a)y(oo, a) = 0, and the renewal condition (1 .59) gives 1000 1000 I ds y(/, O) = r(t)a h,(s, a)g(s, ay( -)y(s, a)-da. o 0 s s (1 .65) (1 .66) Separation of variables in (1 .64) is not immediately possible as it was in the case where the population was structured on size only (equation (1 .45). This is reasonable, in the sense that we cannot expect separable solutions for n(/, a, t) of the form given by (1 .60) if the growth rate varies with time, because ;. = 1 implies that any variation of growth in size with time must also affect age structure. We proceed here assuming that r(t) = 1 80 the growth rate ofindividuals is independent of time, but in later chapters we explore two more general ways of proceeding: • In Chapter 4, we reformulate the problem so that instead of looking for steady size/age distributions, we look for steady si.zelbirth-size distributions, which can exist in the presence of a separable time-dependent growth rate, then calculate the corresponding time-dependent size/age distribution from this . • In Chapter 5, we use a transformed time variable with respect to which steady size/age distributions are possible. Now setting r(/) == 1 , 19 (1 .67) separation of variables in (1 .64) leads to � _ -/ay(l, a) - &(g(l, a)y(l, a» - h/(l, a)g(l, a)y(l, a) _ Q � - �� - . (1.68) An expression for the unknown constant Q may be obtained by rewriting the right hand equality as () () Qy(l,a) = - ()ay(l, a ) - 81 (g(l, a)y(l, a » - h/(l, a)g(l, a)y(l, a). (1.69) . Proceeding as in Section 1.6, we integrate (1 .69) over all I and a and apply the normalisation condition (1 .61 ), the zero-flux condition (1 .65), and the renewal condition (1 .66) to obtain Q = (a - 1 ) f 1000 h/(l, a)g(l, a)y(l, a)dlda. The exponential growth rate of the total population is simply Q, because dN di = QN(t). (1 .70) (1.71) Once again, an alternative expression for Q may be obtained by considering the rate of increase of the total size of all the individuals in the population. Following the approach in Section 1 .6, we have or ! [1000 1000 LN(t)Y(I, a)dlda] = 1000 1000 N(t)y(l, a)g(l, a)dlda, ':: [loOO 1000 ly(l, a)dlda] = N(t) 1000 1000 y(l, a)g(l, a)dlda, then using (1 .71 ) we obtain Q _ 10 1000 g(l, a)y(l, a)dlda - fO 1000 ly(l, a)dlda , (1 .72) (1 .73) (1 .74) 80 once again Q may be regarded as ratio of the average growth rate of the individuals in the population to their average size. Rearranging (1 .69), the PDE for the SSD y(l) becomes () () 8ay(l,a) + 8lg(I, �)y(I, a) = - (h/(l, a)g(l, a) + Q)y(l, a), and substituting r(t) = 1 into (1 .66) the renewal condition at age a = 0 is 10 00 10 00 I ds y(I, O) = a h,(s, aY( - )g(s, a)y(s, a)-da. o 0 s s 20 (1 .75) (1 .76) This system must be solved subject to the normalising condition (1 .61), and the zero flux boundary condition (1 .65). In the special case where all cells divide exactly into a pieces each of the same size, we have f(P) = 6(p - a-I ), (1 .77) so the renewal condition may be written simply as y(/, O) = a2 1000 h,(al,a)g(al, a)y(al, a)da. (1 .78) 1.9 The Simplest Case: A Constant Probability of Division A major difficulty in dealing with (1 .53) or (1 .75) is the presence of the functional Q. In general, finding Q is diffi.cult, because it involves solving a transcendental equation, or some other iterative technique. For example, (1 .53), or (1 .75) in conjunction with (1.76), may be solved assuming a zero-flux boundary condition at 1 = 00 to obtain a solution in terms of the unknown Q. '1b find Q we then apply the zero-flux boundary condition at 1 = 0, which gives in general a transcendental equation known as a characteristic equation (Gyllenberg [17]). This terminology is similar to that used in demography, where a "characteristic equation" is used to determine the intrinsic rate of growth of age-structured populations (Keyfitz [28]). There appear to be just two situations in which Q can be known in advance. The first, simpler but less interesting case, arises when h,(l, a)g(/, a) = b, (1 .79) or, in the case of populations structured on size alone, h,(l)g(l) = b, (1 .80) where b is a constant. Note that this b is not the same as b(l, a, t) and b(/, t) defined in sections (1 .7) and (1 .5) respectively, as in this special case we have b(l, a, t) = r(t)b and b(/, t) = r(t)b. 21 Using (1 .79) or (1.80) and (1.48) or (1 .70), along with the appropriate normalisation condition for y, we then have simply Q = ( a - 1 )b. (1 .81 ) Chapter 2 and Appendix B deal with size-structured populatioDB satisfying (1 .80), in which each cell divides into a daughters all of the same size. In this case we substitute (1 .81) into (1 .55), giving the simple functional differential equation d dlg(l)y(l) = - bay(l) + ba2y(al). (1 .82) In Chapter 2, a series solution is obtained for this equation in the special case where the growth rate g is a constant, and some interesting properties of the solution are outlined. In Appendix B, a series solution is obtained for the case where gel) is given by a power law. 1.10 A Special Case: Cells in One-Dimensional Structured Tis­ sues The second case in which Q may be evaluated immediately is the case where the time­ independent component of the growth rate of cells is proportional to cell size, so that g(l, a) = gl. (1 .83) or, in the case of a population structured on size only, gel) = gl, (1 .84) where the g on the right hand side is a constant. The case where cells grow exponentially in size is a sub-case of this, with r(t) a constant. The constant g is arbitrary, as a given cell growth rate g, specifies what the product r(t)g must be, but neither g nor r(t) is specified uniquely. Indeed, we describe below a situation in which it is natural to choose g = 1 . 'Ib find Q with (1.83) or (1 .84), we substitute into (1 .74) or (1 .52), respectively, to obtain Q = g. (1 .85) One situation in which (1.83) or (1.84) holds is of particular interest to plant scientists and others studying structured tissues. Fig. 1 .5 shows a longitudinal section of part of the 22 • • direction of growth Figure 1 .5: A longitudinal section of part of the meristem of a Zea mays root. This reproduction is taken from the data collected by Erickson for the classical paper by Erickson and Sax [9]. Growth of the tissue involves expansion in the direction shown, with cell division occurring by insertion of new cross-walls at right angles to the direction of growth. 23 meristem of a Zea mays root, which is typical oflongitudinal sections from roots or leaves in a large number of plant species. Growth of the tissue occurs mainly by expansion in one direction, with cell division occurring by insertion of new cross-walls at right angles to this direction. For this reason we refer to tissues like that shown in Fig. (1.5) as one-dimensional structured tissues. Direction of growth .. . Figure 1 .6: An idealised section of plant tissue growing by elongation in one dimension. Cells are arranged in files and are locked together by common longitudinal walls. When cells divide, new cross-walls (dashed lines) form within the parent cell, and two new cells are formed (i.e. a = 2). We assume there are sufficient files at right angles to the direction of growth to allow use of continuous number density and probability density functions. Fig. 1 .6 shows an idealised section of a root or leaf, elongating in a spatially homogeneous manner such that any two points a unit distance apart separate at a time-varying rate r(t). The assumption of spatial homogeneity will usually be applicable only to small regions of the growth zone; in large regions the rate of expansion r(t) will vary spatially in the direction of growth shown. As can be seen from Fig. 1.5, different files of cells in a one-dimensional structured tissue may contain dift'erent types of cells. We restrict our attention here to files of cells of a single type, where the idealisation exemplified by Fig. 1 .6 is a good approximation. We take the size of a cell to be its length, I, in the direction of growth. The important characteristic of a structured tissue in this context is that it grows symplastically, i.e. with longitudinal files of cells locked together so that cells in adjacent files do not slip past each other. This means that all cells in a small region can be assumed to grow at the same relative growth rate r(t), so the growth rate for a cell of size I is gl(l, a, t) = Ir(t), 24 (1 .86) or, for size-only models, gJ(I, I) = lr(I). (1 .87) The time-independent part of the growth rate is then given by (1 .83) or (1 .84) with constant g = 1 . There are other situations where the time-independent component of the cell growth rate is proportional to cell size. In populations of bacterial cells, an exponential growth rate oftbe form (1 .83) or (1 .84) has been assumed by Tyson and Diekmann [61], while experimental evidence shows that exponential growth in cell size is a good approximation for some types of bacteria (Schaechter et al [48], Grover et al [13], Painter and Marr [41 ]), but not for others (Harvey et al [22], Anderson et al [3], Collins and Richmond [7], Grover et al [13], Painter and Marr [41]). In the case of a size-structured cell population, substitution of (1 .84) and (1 .85) into (1 .53) leads to (1 .88) or in the case where all cells divide into exactly a pieces each of the same size, we obtain from (1 .55) the functional differential equation : = -(hJ(l) + � )y(l) + rl hJ(al)y(al). (1 .89) In Chapter 3, (1 .89) is solved subject to the normalisation condition (1 .43), for a range of forms of the hazard rate hJ(I). For a population structured on size and age, we substitute (1 .83) and (1.85) into (1.75) to obtain the partial differential equation fiy fiy Ba + gl m = -(glhJ(l, a) + 2g)y(I,a), (1 .90) and from (1 .76) the renewal condition at a = ° is 1000 1000 I y(I, O) = ga hJ(s, a)f(-)y(s,a)dsda. o 0 s (1 .91) However, as pointed out in section (1 .8), the solution of these equations for the steady size/age distribution y(l, a) is valid only in the case where the growth rate of individuals is independent of time, so r(l) = 1 and the growth is strictly exponential. In Chapter 4 we reformulate the problem in terms of cell size and cell birth-size, instead of cell size and age, in order to obtain more widely applicable solutions. 25 In Chapter 5, we look at the question of whether the steady size, size/age, or sizelbirth-size distributions arising in this case are stable. In other words; given any initial distribution of cell size and age or birth-size, we look for conditions under which the size (or size/age or sizelbirth-size) distribution will approach the calculated steady distributions. 26 Cbapter 2 The Equation y' (x) -ay(x) + aay( ax) This chapter is a slight modification of a paper published in the Journal of the Australian Mathematics Society, series B (Hall and Wake [18]). The introduction has been abbreviated as most of this material is covered in Chapter 1 (particularly Sections 1 .5,1.6, and 1 .9), but some extra mathematical detail, which was omitted in the paper for the sake of brevity, has been inserted into other sections. The notation used is the same as that used in Chapter 1 except that we use x instead of I to represent cell size. 2.1 Introduction We consider here a cell population structured on one size variable x only, where each parent cell divides evenly to produce exactly a daughters all of the same size. Substituting f(P) = o(p - a-I ) into (1 .37) in Section 1 .5, the partial differential equation for the number density n(x, t) is {J {J {J/n(x, t» + ax (g(x, t)n(x, t» = - b(x, t)n(x, t) + a2b(ax, t)n(ax, t) . (2.1 ) We now look for a steady size distribution, y(x), which can develop under the conditions outlined in Section 1 .6. After following the development in that section, we replace h,( l)g(l) in (1 .55) by the single function b(x) to obtain a functional differential equation for y(x), d dx (g(x)y(x» = - [b(x) + Q]y(x) + �b(ax)y(ax), (2.2) 27 where Q is a functional on y(x) given from (1 .48) by Q = (a - 1 ) 1000 b(x)y(x)dt, and y is a probability density function 80 must satisfy the normalisation condition 1000 y(x)dt = 1 . (2.3) (2.4) The case we consider here is the simplest non-trivial case, with b(x) and g(x) both constants, say b and g respectively. As pointed out in Section 1 .9, in the case where b(x) is a constant the value of the functional Q can be calculated in advance from (2.3) to be Q = (a - 1 )b. (2.5) While this case where g and b are both constants is not particularly realistic, its solution does offer some ideas which can be further explored in more realistic examples, as in Chapter 3. The equation to solve then becomes, after simplification using (2.4), g: = -bay(x) + ba2y(ax), which is just (1 .82) with g( l) replaced by the constant g. (2.6) Finally, assnming g :f 0 we can substitute a = a:, forming the functional differential equation : = -ay(x) + aay( ax), (2.7) which we wish to solve subject to the normalising condition (2.4) and a zero-flux boundary condition, which from (1 .47) is written y(O) = y(oo) = O. (2.8) Kato and Mcleod [27] have studied the equation y'(x) = ay(Ax) + by(x) subject to the initial condition y(O) = Yo. Here we show how with the aid of the integral condition the solution can be obtained directly by application of the Laplace Transform technique, and that the application of the integral condition (2.4) (instead of the more usual initial condition y(O) = Yo) leads to a solution which can be regarded as a probability density function (pdi), which has a number of interesting properties not previously investigated. We do not 28 restrict our interest to the case a = 2 which is usually used for cell growth, but rather treat (2.7) as an equation in which a can take on any real value a > 1 . Fractional values, a = � say, could be regarded as representing a conbination of fusion processes, in which a number (q) of individuals fuse together to form one large entity, then this large entity divides into a number (p) of equal parts. We show that pure fusion processes, which would require a < 1 , cannot have a steady size distribution y(x) satisfying (2.7). 2.2 Preliminaries We now proceed to study the solution of (2.7) and (2.4) in the region x � 0, ignoring the manner in which these equations were derived. As a first step in this process, we note that we may restrict our attention to a > 0, as y(tu) is undefined for a < 0, and a = 0 leads trivially to the simple solution y = ae-II%. We may also assume a :F 1 , as if a = 1 then (2.7) reduces immediately to y(x) = 0 which has no solutions satisfying (2.4). We shall first derive necessary properties of the solution, assuming at this stage that it exists. 2.2.1 y(O) and y( 00 ) Condition (2.4) does not directly provide us with the initial condition y(O), but we can show that the realistic assumption y(O) = y( (0) = 0 can be derived from (2.7) and (2.4) alone. H we integrate (2.7) then take limits as x - 00 we get and so lim (y(x) - y(O » = lim (-a r'y(s)ds+ aa r' y(as)ds), %-+00 %-+00 Jo Jo %�(y(x» = -a + a + y(O), %�(y(x» = y(O) . (2.9) Now if y(O) :I 0 then y( (0) :I 0 so Jooo y(x)dx does not converge, contradicting (2.4). Hence we must have y(O) = 0 (2.10) 29 and lim (y(x» = O. %-00 (2.11 ) The first of these results gives us the appropriate initial condition, and the second ensures that y(x) is of exponential order so Laplace transforms can be applied. It is interesting to note that ify(O) = 0, all derivatives ofy are also zero at x = O. This can be shown iteratively by first substituting x = 0 into (2.7) to show y' = 0, then differentiating (2.7) and substituting x = 0 to show y't = 0, and so on. Hence the solution is "infinitely flat" and therefore not analytic at x = o. 2.2.2 Restrictions on y, a, and a Given y(O) = y( 00) = 0, and 10 y(x)dx = 1 > 0, it follows that y(x) is positive for some range of x values. It therefore has a maximum at some point X say, so from (2.7) we have o = -ay(x) + aay( ai), and so y( ax) = y(x) . a (2.12) Hence for a < 1 we would conclude y(ax) > y(x) contradicting the assumption that y has its maximum at x. Hence we must have a > 1 , (2.13) and no SSD can exist in the fu,sron case, where a < 1 . We can show y(x) > 0 (\Ix > 0) by considering the variable Z(x) = f%oo y(s)ds. Integrating equation (2.7) leads to 100 y(s)ds = -aZ(x) + aa 100 y(as)ds, so Z(x) = -aZ(x) + aZ(ax). (2.14) This is subject of course to the conditions Z(O) = 1 (from (2.4» and Z( 00) = o. 30 Assume for the moment that there exists a point Xl E (0, 00) at which Z'(Xl ) = -Y(Xl ) = 0, and Z(Xl ) i- O. Then Z( aXl ) = Z(Xl ) and there must be a point X2 � aXl such that Z'(X2) = 0 and IZ(X2)1 � lZ(xt } l· Without loss of generality, assume Z(Xl ) > O. Then one of the following must hold: • (a) Z'( aXl ) > 0 in which case such a X2 must exist, or • (b) Z'(axl ) = 0 so we can choose X2 = axl , or • (c) Z'(axt } < 0 so 3x* E (Xl , axt ) such that Z'(x*) = 0 and Z(x*) > Z(xt }, so Z(ax*) > Z(axl ) also. Then we have axl < ax* < 00 but Z(ax*) > Z(axl ) and Z(ax*) > Z(oo), 80 there must be some point X2 > axl where Z'(X2) = 0 and Z(X2) > Z(Xl ) as required. The argument for Z(Xl ) < 0 is almost identical. Hence we can find a sequence of points (XII) -- 00 as n -- 00, at which Z'(xlI) = 0, whose images under the function Z, (Z(lIt» f+ O. This contradicts Z( 00) = 0, so it follows that Z has no non-zero stationary points on (0, 00), so Z is monotonically decreasing with Z'(x) � 0 ("Ix > 0). As y(x) = -Z'(x) we have for X > 0 y(X) � 0, so y(x) is a probability density function as desired. Now equation (2.7) can be written as the integral equation y(X) = a 1«% y(s)ds, so as y(x) > 0, a necessary condition for the existence of an SSD must be a > O. 2.2.3 Uniqueness of the Solution (2.15) (2.16) (2.17) Let there be two solutions to (2.7) and (2.4), Yl (x) and Y2(X). Then the variable Z(x) = fzoo(Yl (s) - Y2(s»m must satisfy equation (2.14) but this time with the conditions Z(O) = 0 and Z( 00) = O. Now if Z(x) i- 0 for any x, then there must be a stationary point Xl of Z such 31 that Z(Xl ) :F 0, which by a similar argument to that given above leads to a contradiction. Hence we have Z(x) == 0 so Z'(x) == 0 and Y1 (x) == Y2(X), so there can be at most one solution to (2.7) and (2.4). 2.2.4 Moments of the SSD (2.18) Because of the integral condition (2.4), it is in fact easier to find the moments of the SSD, than the function y(x) itself. The nth moment ofy(x) about the origin is p" = 1000 x"y(x)d% so we multiply (2.7) by X' and integrate, giving 1000 i'''(x)dx = -ap" + aa r. x"y(ax)d% and (2.19) (2.20) (2.21 ) where we have used integration by parts on the left and the substitution z = ax on the right. This procedure is equivalent to taking the Mellin Transform. of each side if the equation. We now rearrange (2.21) to obtain an iterative equation for J.4., (2.22) with from (2.4), JIG = ! . (2.23) Now iterating n times, the nth moment ofy about the origin is given by n! (2.24) p" = a" 11::1 (1 - a-I) ' Hence the mean cell size, /J, is given by 1 (2.25) JJ = a(l _ a-I ) , 32 and the variance of the cell size, UJ, is The third moment about the mean, ma, can be shown from (2.24) to be 2 ma = Q8(1 _ a-i) ' 80 the coefficient of SkeWneB8 �fJt is given by 1fJJ. _ ma _ 2(a + 1 )�a2 - 1 v' - a:a - a2 + a + 1 Thus the coefficient of skewnes8 is dependent on a only. 2.3 Solution of the SSD Equation (2.26) (2.27) (2.28) We now establish that equations ( 2.7) and (2.4) have a solution, which by section (2.2.3) is unique. We first take Laplace transforms of each side of (2.7), remembering y(O) = 0, giving (2.29) so PY(P) = ay(�) - ay(P). (2.30) Thus we have an iterative equation with an initial value given by the integral condition (2.4), with Iterating n + 1 times, we get 1 p y(P) = 1 + eY(;)' tJ y(O) = 1 . y(P) = (1 + �)(1 + �) ... (1 + faal( a�+l )' 33 (2.31) (2.32) (2.33) Now as a > 1 we know that so we get �Y(;+I ) = y(O) = 1 , (P) _ lim TIN 1 Y - N-ooo,,:o (1 + a-"�r Expressing this product 88 a sum of partial fractions we can write N 1 N Q.( -aa") TI (1 + a-"! ) = L (1 + a-"! ) ' ,,:0 • tc:O • (2.34) (2.35) (2.36) where Q,,(P) is the Pt:Oduct on the left hand side of this equation excluding the factor involving a-". Thus, 1 Q,,(P) = (1 + �)(1 + a-1 !) ... (1 + a-(,,-I )!)(l + a-("+1 )�) .. . (l + a-N�) (2.37) and Q,,( -aa") = (1 _ ate)(l _ ate-I ) .. . (1 _ a)(l � a-I )(1 _ a-2) ... (1 _ a-(N-,,» · (2.38) Simplifying and taking limits, we get 1 00 1 1 y(P) = K ,,� (1 - a)(l - a2) . . . (1 - ate) (1 + a-"�) ' where the product in the denominator is given the value 1 when 11 = 0, and 00 K = TI (1 - a-"). 11=1 Now taking the inverse laplace transform term by term we get ( ) 1 � 1 a" -11th Y X = K II� (1 _ a)(l _ a2) . . . (1 _ a")a e , so noting that a > 1 , we write a 00 ( -1 )"a"e-a-x y(x) = K II� (a - 1 )(a2 - l ) . . . (ate - 1 ) " (2.39) (2.40) (2.41) (2.42) This solution is in the same general form as that given by Kato and Mcleod [27], but we have simplified the notation slightly by using the convention that when 11 = 0, the empty product in the denominator is given the value 1. 34 2.4 Discussion 2.4.1 The Constant K The infinite product in (2.40), which gives the value of K, has been studied in some depth by number theorists and others. Perhaps the best known result involving this product, and one which offers a rapidly convergent way of calculating K, is Euler's Pentagonal Number Theorem (Andrews [2]). This is 00 il (l - q") = f (�1 )"q",..-l) 11=1 where Iql < 1 , so we have 11=-00 = 1 - (q + q'A) + (t! + q7) _ (tl'l + q15) + ... K = 1 - (a-1 + a-2) + (a-5 + a-7) _ (a-12 + a-15) + .. . . (2.43) (2.44) Another interesting (and rapidly convergent) series for K can be deduced by substituting p = 0 in (2.39) and using (2.32), giving K = I, (-1)" , 11=0 (a - 1 )( a2 - 1 ) . . . ( a" - 1 ) which could also have been deduced by substitution of q = a-I and z = -q in 00 qill(II-1)z" 00 !o (1 - q)(1 _ if) ... (1 _ q") = !! (1 + zq"), (2.45) (2.46) which is also given by Andrews [2]. Once again, we are using the convention that the product in the denominator is given the value 1 when n = O. Another expre88ion for K, but one we do not actually use here, can be derived from a formula involving Jacobian elliptic functions. Formula 16.37.1 in Abramowitz and Stegun [1 ] leads to (2.47) which could be used to calculate K if we make the substitution q = a-i. This transfers the problem to relating the elliptic function parameter m to the nome q, which is possible via Landen's transformation. 35 2.4.2 Properties of the Solution From a cursory inspection of (2.42) alone, it is not at all obvious that y(x) satisfies all the properties that we deduced earlier it must have. Firstly, we calculate y(O) by substituting x = 0 in equation (2.42), which leads to a � (-1 )" a" y(O) = K !:o (a - I )(a2 - I ) . . . (a" - 1 ) " (2.48) Consider the finite form. of the sum on the right hand side. Expanding the a" in the numerator as (a" - 1 ) + 1 , we get f (-1 )" a" 11::0 (a - I )(a2 - I ) . . . (a" - 1 ) f ( -1 r( a" - 1 ) + f (-1 r (2.49) = ,.::() (a - 1 )(a2 - 1 ) . . . (a" - 1) Il=O (a - 1 )(a2 - 1 ) . . . (a" - 1) " We then drop the first term. in the first sum on the right (which is zero), which means that we can cancel (a" - 1) from numerator and denominator then change the summation variable to m = n - 1 so f (-I)"a" 11=0 (a - 1)(a2 - 1 ) . . . (eJ'I - 1) = _ Nf (-I f' + f (-I r m=O (a - 1 )(a2 - I ) . . . (aM - 1) Il=O (a - I )(a2 - I ) . . . (a" - 1) (-I f = (a - I )(a2 - I ) . . . (aA' - 1) " (2.50) Taking limits as N -+ 00, by substituting (2.50) into (2.48) we have y(O) = � lim (-If = 0, K N-oo (a - I )(a2 - I ) . .. (aN - 1) because a > 1, so the denominator grows without limit as N -+ 00. (2.51 ) 1b verify by substitution that (2.42) satisfies the integral condition (2.4), we interchange the integral and summation signs (which is valid in view of the uniform convergence of the series), so 1000 a 1000 00 ( -1 r et'rlla"% O y(x)d% - - L ��-:---:;.----::-�,--� - K 0 11=0 (a - I )(a2 - 1) .. . (all - 1) = !. f (-1)" K ,,::0 (a - 1 )( a2 - 1 ) . . . ( a" - 1) , 36 then apply (2.45), giving (2.52) as required. Verifying the expression for the nth moment (2.24) by direct substitution of (2.42) in Jl,t = fooo x"y(x)d% requires a similar interchange of order, Jl,t = looo x"y(x)d% = a � ( -1 'ra'" 100 x" -ot-%dx K".-:O (a - 1 )(a2 - 1 ) ... (tr"- 1 ) 10 e n! 00 ( -1)'" = Kd' lo (a - 1)( a2 - 1 ) ... ( aW'-- 1 )a- . We now apply (2.46) with z = _a-(II+I) and q = a-I then simplify further using (2.40), which gives n! = d'(1 - a-I )(1 - a-2) • • . (1 - a-II) , (2.53) as required. 2.4.3 An Example Consider a population of cells undergoing binary fission (a = 2) with a constant growth rate g of 0.2 size units per second, and a division rate b of 0 .1 per cell per second. Then and 00 1 11 ba a = - = 1 g K = 1£1 - 2 ) 11= " 1 1 1 1 1 1 = 1 - (2 + 4) + ( 32 + 128) - (4096 + 32768 ) + . . . ::::: .288788 . 37 The mean, variance, and coefficient of skewness of the SSD are obtained by substitution in (2.25), (2.26), and (2.28) respectively, giving JJ = 2, (2.54) a2 4 (2.55) = 3 � 1 .3333, and v'/Jt 6-/3 (2.56) = 7 � 1 .4846. .5 .4 .3 - )( - >- .2 . 1 6 7 Figure 2.1 : The solution for a = 1 and a = 2. The dashed line is the graph ofy(x) = ie-G. Figure 2.1 shows the shape of the SSD functiony(x) in this case as given by (2.42). The qualitative features of Figure 2.1 remain for all a and a. The parameter a acts as a simple scaling factor for x (the substitutions t = ax and z(t) = y(x) make equation (2.7) independent of a). Equation (2.28) shows that the coefficient of skewness is always positive, increasing from 0 to 2 as a increases from 1 to 00. 2.4.4 The Solution for Extreme Values of a In the light of the comments at the end of the example, it is worth considering the form of the solution when a - 00 and a - 1 +. The first is easy - if we write (2.42) in a form 38 similar to that given in Kato and Mcleod [27], namely a _Q 00 (-1 )1Ia"e-4I(a"-lF y(x) = Ke (1 + II� (a - 1 )(a2 - 1 ) ... (aII - 1 » ' (2.57) then it is clear that for all a, y(x) � fe-tlZ for large x, and that in particular for all x > 0 we have (2.58) Hence in the limit as a -- 00, y(x) approaches an exponential pdfwith parameter a. 1 .g .8 .7 .6 - N .5 "-' � .... . 3 .2 . 1 -4 3 4 z Figure 2.2: The normalised pM Y(z) = ay( az + Jl) for various values of a. The dashed line is the pdf of the standard normal distribution. The limit as a - 1 + is much more difficult. Inspection of (2.42) shows that each term in the series becomes infinite as a -- 1 +, and (2.24) shows that all the moments become infinite as well. This suggests that we "normalise" the equation using the transformation and x - Jl z =-­a Y(z) = oy(x), where Jl and a are given by (2.25) and (2.26). 39 (2.59) (2.60) Figure 2.2, which shows Y( z) plotted for various values of a, suggests that in the limit as a - I +, Y(z) approaches the pdfofthe standard normal distribution N(O, I ). Substituting (2.60) and (2.59) into (2.4), Y must satisfy the normalisation condition Now differentiating (2.60) and (2.59), we get .J( ) _ r(:) J x - -:;r-, C1 and making the appropriate substitutions in y(ax) = Y(z(ax)) we have ( ) _ Y(az + (a - l );) y ax - . C1 Now substituting (2.62) and (2.63) into (2.7) we obtain yI(z) = _a Y(z) + aa Y(az + (a - 1 )i) er- C1 C1 ' (2.61 ) (2.62) (2.63) (2.64) which after substitution for Jl and C1 using (2.25) and (2.26) gives "normalised" functional differential equation (2.65) Now to find out what happens as a - I +, we expand the Y( az + .J a2 - 1 ) as a Taylor series about Y(z) up to all terms of order (a - 1 ), giving 80 yea: + v' a2 - 1 ) = Y(z + [(a - 1): + v' a2 - I D, (2.66) yea: + Ja2 - 1 ) � Y(z) + [(a - 1 )z + .j a2 - IJY'(z) + a2 2 - 1 Y" (z). (2.67) Substituting this expression into (2.65) and simplifying, we get aY(z) + a2zy'(z) + [(a + 1 ).j a2 - 1] Y'(:) + a2(� + 1 ) Y" (:) = O. (2.68) Substitution of a = 1 now leads to the second order differential equation Y"(:) + zY'(:) + Y(z) = 0, 40 (2.69) with an integral condition given from (2.61 ) by 1:00 Y(z)dz = 1 , (2.70) because from (2.25) and (2.26) we have lima-+l (Jl/O') = 00. It is trivial to show that Y(z) = Yl (Z) = e-� (2.71) is a solution of (2.69). A second (independent) solution can then be found by substituting say Y(z) = Y2(Z) = U(Z)Yl (z), leading (after cancellation Of Yl (x» to U"(z) - zU'(z) = o. (2.72) One solution of this equation is U'(z) = e� , (2.73) so we may choose (2.74) The general solution of (2.69) is then (2.75) where Cl and C2 are arbitrary constants, and Yl (Z) and Y2(Z) are given by (2.71) and (2.74) respectively. Now applying (2.70), we note that Icoo Y2(Z)dz is infinite for all c, so we must have C2 = o. Hence Cl = Ai, and the only solution to (2.69) and (2.70) is 1 _.t Y(z) = rn=e T , v21r which is the pdf of the standard normal distribution as expected. (2.76) Consideration of the neglected terms in (a - 1) shows that we can write, for all a close to 1, (2.77) 41 Reverting to the original coordinates, we can then say that for a close to 1 , y(x) = " �e- (·;JI (1 + o(Ja - 1 ») , (2.78) where J.l and a2 are given by (2.25) and (2.26) respectively. Thus as a - 1 +, y(x) approaches the pM of a normal distribution. Clearly, if the limit as a - 1 + of the SSD is a normal distribution, the expression given in (2.24) for the nth moment of the distribution should tend to the expression for the corresponding nth moment of the normal distribution. In Appendix A, we show that this is true provided that lim [ (1 + q) 1 I,(-1t n! ] _ { 0 , n odd q_1- 1 - q t=O (n - k)!k!q - (n - 1 )! ! , n even (2.79) This result is then proved, as this has apparently not been done elsewhere in the combinatorial literature. 2.4.5 The Form of the Time-dependent Number Density n(x, t) Assuming that r( t) = 1 80 the growth rate is independent of time, then substituting b(x) = b and Q = (a - 1 )b) into (1 .49), gives so dN di = N(t)( a - 1 )b, (2.80) N(t) = Noe(a-1)bt, (2.81 ) where No = N(O). We can now substitute (2.42) and (2.81 ) into n(x, t) = N(t)y(x) to obtain the separable form of n(x, t) which applies when the size distribution is given by the SSD, n(x, t) = N(t)y(x) = aNo (a-1 )bt I, (-1 )"a"e-a-z K e lI=o (a - 1 )(a2 - 1 ) .. . (a" - 1 ) " 42 (2.82) Chapter 3 Size-Structured Models for Cells Growing Exponentially This chapter is a slight modification of a second paper published in the Journal of the Australian Mathematics Society, series B (Hall and Wake [19]). The introduction has once again been abbreviated, the first part being in fact identical to that of Chapter 2. The remainder of the paper has been slightly expanded, with some extra mathematical detail being inserted. The notation used is the same as that used in Chapter 1 except that we use x instead of I to represent cell size. 3.1 Introduction We consider here, as in Chapter 2, a cell population structured on one size variable x only, where each parent cell divides evenly to produce exactly a daughters all of the same size. Substituting/(p) = 6(p - a-I ) into (1 .37) in Section 1.5, the partial differential equation for the number density n(x, t) is {J {J {J,(n(x, t» + {Jx (g(x, t)n(x, t» = -b(x, t)n(x, t) + a2b(ax, t)n(ax, t). (3.1 ) We once again look for a steady size distribution, y(x), which can develop under the conditions outlined in Section 1 .6. After following the development in that section, we replace h,(l)g(l) in (1.55) by the single function b(x) to obtain a functional differential 43 equation for y(x), d dt(g(x)y(x» = -[b(x) + Q}y(x) + a2b(ax)y(ax), (3.2) where Q is a functional on y(x) given from (1 .48) by Q = (a - 1 ) 1000 b(x)y(x)d:t. (3.3) We note here that all cells must be of positive size, so we may set y(x) == 0 for x < 0, and as y(x) is a probability density function it must satisfy the conditions and 1000 y(x)dx = 1 , (3.4) y(x) � O. (3.5) It is difficult to see how to obtain closed-form solutions to (3.2) subject to (3.4) and (3.5) in the most general case, most obviously because of the presence of the unknown functional Q. However, as pointed out in Sections 1.9 and 1 .10, in two simple cases we can find Q readily. The first of these is if b(x) = b, a constant, so that from (3.3) Q = (a - 1)b. A solution for this case with g(x) also constant was given in Chapter 2. In Appendix B, a solution with b(x) constant and g(x) = gr-k (k > 0) is obtained by using the transformation Z(x) = g(x)y(x) to obtain an equation of the same form as equation (3.20) below. As discussed in Section 1.10, the case where g(x) = gx, (3.6) where the g on the right hand side is a constant, is of particular importance, both for cells in one-dimensional tissues and also for some types of bacterial cells. In this case, as pointed out in Section 1.10, (3.3) may be replaced by Q = g. (3.7) Ifr(t) = 1 , then g(x) = gx corresponds to the case where the cells grow in size exponentially during their life cycle. Substituting (3.6) and (3.7) into (3.2), we obtain a functional differential equation similar to (1.89), gxy'(x) = -[b(x) + 2g)y(x) + a2b(ax)y(ax). (3.8) 44 In this chapter we seek solutions to this equation, subject to the normalising and non­ negativity conditions (3.4) and (3.5). 3.2 Preliminaries In Appendix D, it is shown that y = 6(x), (3.9) where 6(x) is the Dirac delta function, is a generalised solution of(3.8) for all forms of b(x), and is in fact a solution of (3.2) for any g(x) such that g(O) = O. Further, it is shown that y(x) = 6(x) is the only non-negative1 generalised solution satisfying (3.4), and therefore the only generalised solution'which may be regarded as a probability distribution. The existence of this solution is reasonable in that it says that if we have a population all of size zero at any time, this situation will not change provided g(O) = o. We now look for classical solutions to (3.8) subject to conditions (3.4) and (3.5). Multiplying (3.8) by x and rearranging gives g[x2y'(x) + 2xy(x)] = -xb(x)y(x) + a2xb(ax)y(ax). This suggests the change of variable to give Setting simplifies (3.12) to Z(x) = .ry(x) gZ'(x) = _ b(x)Z(x) + b(ax)Z(ax). x x a(x) = b(x) gx Z'(x) = -a(x)Z(x) + aa(ax)Z(ax). (3.10) (3.11) (3.12) (3.13) (3.14) 1 A generalised function /(x) is defined as non-negative if for all appropriately -BJDooth- Don-negative teat functions ;(x) we have f�oo/(x);(x)dz � O. 45 This equation needs to be solved subject to Z(X) � 0 and either of the conditions or 1000 Z�:) dx = 1 , 100 a(x)Z(x) dx = 1 10 x a - I ' (3.15) (3.16) Here the second condition follows either from Q = g, or from integration by parts and substituting from (3.14). Equations (3.14) and (3.15) together are equivalent to (3.14) and (3.16) together. Note that equation (3.14) can be written in the integral form Z(x) = 1"" a(s)Z(s)dJ, (3.17) which despite superficial similarities is very different from the classical Volterra integral equation. It follows that Z(O) = 0, (3.18) then substitution into (3.14) and its derivatives shows that z{1I}(0) = 0 , "In � O. (3.19) Hence any solution to (3.14) must be non-analytic (-mfinitely flat") at the origin. In the next three sections we examine solutions of (3.14) subject to (3.15) or (3.16) for different forms of b(x). We find solutions involving series of exponentials for the two cases b(x) = bxl and b(x) = bxlH(x - Xl ), and then describe a general method of solution in cases where b(x) is zero below some mjnjmum size for division XI . and becomes infinite at some maximum cell size X2. A general method for solving functional differential equations of the form (3.14), involving series of multiple integrals, is outlined in Appendix C. / 46 3.3 The case b(x) = bxk In this case, we have a(x) = � = axl-1 say (where a = :), so that (3.14) becomes Z'(x) = xk-l (-aZ(x) + a�Z(ax» . We now make the substitutions z = � and Y(z) = Z(x), so Z'(x) = Y'(z): = xk-1Y'(z) and Z(ax) = Y« a;)k ) = Y(�z), then (3.20) becomes (after cancellation of xk-1 ) Y'(z) = -aY(z) + a�Y(�z). (3.20) (3.21) (3.22) (3.23) This equation is in the same form as that dealt with in Chapter 2, except that the normalisation required will be different. Using (2.42), we can therefore write the solution immediately, which in terms of the variables used here is 00 ( -1 )" ab' _01-, Y(z) = c ,,� (at - 1 ) . . . (ab' _ 1 )e , (3.24) where by convention the expression in the denominator is given the value 1 when n = o. This solution is unique apart form the multiplicative constant C. Now reverting to our Z and x notation we obtain _ 00 ( -1 )"ab' -G� Z(x) - c 1 (at - 1) ... (ab' _ 1 )e The probability density function for the SSD is then given (from (3.11» by Z(x) y(x) = -2 ' X (3.25) (3.26) with the constant C to be evaluated from condition (3.16), which can be written in this case as (00 �k-2 1 lo X"� Z(x)dt = a(a - 1 r 47 (3.27) Note that if k = 0 the transformation, and therefore the solution, is ill-defined, and for k < 0 the series does not converge. The solution is therefore valid only for k > O. (3.28) Obtaining the solution (3.25) was straightforward because we managed to transform the original problem into one for which the solution was known. However, using the integral condition (3.27) to find the constant C is not so easy, as it involves carrying out the the integration in (3.29) for all k > O. For 0 < k 5 1, care must be taken over interchanging the integral and summation signs in this expression, as Jooo xPe- %dx does not exist for p 5 -1 . It is necessary to treat two cases separately: k = m-I for some integer m, and k f; m-I for any such m. 3.3.1 The Case k =f:. m-I Firstly, consider the integral Jooo xPZ(x)dx, where p f; -1 . Integrating by parts, we get 1000 [ xP+I 1 00 1 1000 X'Z(x)dx = --lZ(x) - --1 X'+1Z'(x)dx. o p + 0 p + 0 (3.30) Now the square-bracket term is zero because Z is infinitely fiat at the origin, and l'Hopital's rule together with the zero-flux boundary condition can be used to show that lim%_ooxPZ(x) = 0 for all realp. We can then substitute from (3.20) for Z'(x) to give 1000 X'Z(x)dx = -p! 1 [-a 1000 x'+lZ(x)dx + , aa1 1000 X'+lZ(ax)dx] , (3.31) then change variable from x to ax, giving foo X'Z(x)dx = a(l - a-CP+ l » foo X'+lZ(x)dx. Jo p + 1 Jo (3.32) If we apply this iteration m times, starting with p = k - 2, (and remembering that km # 1 , 'im E [+) we get foo xl-2Z(x)dx = ""(1 - a l-k)(l - al-21) • • • (1 - ai-IRk) foo xCm+l )k-2Z(x)dx. (3.33) 10 (k - 1 )(2k - 1 ) . . . (mk - 1 ) Jo 48 Now for all k > 0, it is possible to choose m E 1+ such that k > "'�i - for k > 1 we would choose m = O. With such a choice of m, we can now substitute for Z(x) and interchange the summation and integral signs to give a"'(1 - a1-k)(1 - al-2A:) . . . (1 _ a1--) C Ck - 1 )(2k - 1 ) . . . (mk - 1 ) . I, (-1 ) lIaD (00 X<"'+1)A:-2e-a� dx. 11=0 (al - 1 ) ... (ab' - 1) 10 (3.34) Note that for k < "'�i the interchange would not be permitted as the resulting integral would not be finite. Making the natural substitution z = iaMxk we get and 80 r-1 = (;) "'-1 �(1-->z"'-1 -Ilk xk-1 dx = �dz, a = c (!)"'-l a"'-1 (1 - a1-k)(l - al-2A:) . • • (1 - a1--) a (k - 1 )(2k - 1 ) . . . (mk - 1) 00 (-1 )lIaII(l-_> (00 1 • II� (ak - 1 ) .. . ( ab' - 1 ) 10 z"'- e-'dz. In order to simplify this complicated-looking expression, we define 00 (-1 )11/1' K(a,fJ) = It;' (a - 1) . . . (� - 1 ) ' (3.35) (3.36) (3.37) (3.38) which may be regarded as a generalisation of the definition of the constant K in 2.40. The function K( a, fJ) defined in (3.38) has a number of useful properties. Firstly, if we take the well-known identity (e.g. Andrews [2], Goulden and Jackson [14]) 00 qlll(lI-l )r" 00 1: ( q2 = n(l + uf), 11::0 (1 - q) 1 - ) ... (1 -q") 11::0 and substitute q = a-k then we get 00 aD'zII 00 II� (at - 1 )(a2A: - l) . . . (aM - 1 ) = !l (1 + za-b) . 49 (3.39) (3.40) Substituting fJ = -akz gives f (-1 )"/1' = fi (1 - pa-b) ,,::0 (at - 1 ) ... (aM - 1 ) 11=1 so that using the definition of K given in (3.38) we obtain 00 K(a1,fJ) = n (1 - pa-ka). ,,=1 Secondly, substituting fJ = at" into (3.42) gives 00 K( a1, at") = n (1 _ a-1(,,-1 )+--t), ,,=1 so Now by definition, rep) = 1000 x,P-le-JCdx, so using (3.38) we can write (3.37) 88 1000 xt-2Z(x)th (3.41 ) (3.42) (3.43) (3.44) (3.45) = c (�)"'-l a"'-� (1 - a1-1)(1 - al-2A:) . . . (1 - al-lfIl) K(ak �-IfIl)r( 1 _ !.) a (k - 1 )(2k - 1 ) . . . (mk - 1 ) , m + A: = C (�) -1 (1 - a1-1)(! - al-2A:) • • • (1 - l a1--) K( a1, al--)r(m + 1 _ !. ) . (3.46) a a(1 - I )(2 - i') . . . (m - I) k Now we have 1 1 1 1 1 r(m + 1 - i) = (m - i ) ' . . (2 - i )(1 - i )r(1 - i) and, using (3.44) repeatedly, so (3.46) can be simplified to fOO (k) -1 1 1 10 xt-2Z(x)dx = C a ;K(a1 , a)r(l - i)' 50 (3.47) (3.49) So from the integral condition (3.16) we get (Ie) l 1 C = a (a - 1 )K(al, a)r(l - if (3.50) Note that using (3.42), we could write K(ak, a) in a somewhat neater form as the infinite product 3.3.2 The Case k = m-I 00 K(�, a) = TI(I - aI-a). 11=1 (3.51) For le = �, m E 1+, if we start with p = le - 2 we can only apply the iterative relationship (3.32) used in the previous section m - I times, giving tX> -2 _ a"'-l (1 - al-k)(1 - al-2k) • • • (1 - a1-(Ift-I).t) tx) -1 Jo xk Z(x)dx - (le - 1 )(2/c - 1 ) . . . «m - 1 )1e - 1 ) Jo x Z(x)dx. The integral on the right can be integrated by parts once more, giving loOO x-1Z(x)d% = - looo 10gXZ'(x)d% = a [loOO (logx�-IZ(x)dx - � looo (logX�-IZ(ax)dx] (3.52) = a [loOO (logx�-IZ(x)dx - looo (log;f - IOgaV-1Z(s)ds] , (3.53) where once again we have used (3.14) to replace Z'(x), so (3.54) Now substitute for Z(x) from (3.25) and interchange the integral and summation signs (which is now permitted) to give (3.55) Now fooo se-Ills = 1, so using (3.38) this simplifies to 1000 x-1Z(x)dx = CK(�,l ) log a. Thus (3.52) becomes /00 -2 a"'-1(1 - al-k)(l - a1-2k) . . . (1 -al-(",-1)k) 10 xl Z(x)dx = CK(�,l ) loga (t - 1 )(21 - 1 ) . . . «m - l)t - 1 ) . Applying the integral condition (3.27) leads to C _ (t - 1 )(2t - 1 ) . . . «m - 1 )t - 1) 1 - (1 - a1-k)(1 - a1-2k) . . . (1 - a1-(",-1)k) a"'(a - l )K(ak, l ) loga · The numerator in the first fraction can be simplified by taking out m factors of t, so (3.56) (3.57) (3.58) 1 1 1 1 1 (t - 1 )(21 - 1 ) . . . « m - 1 )t - 1 ) = (-1 j- t"'i( i - 1 )( i - 2) . . . (i - (m - 1», (3.59) then replacing i by m we have (k - 1 )(2k - 1 ) . . . « m - 1 )t - 1 ) = (-1 j-1k"'m(m - 1 )(m - 2) . . . (1 ) = (-1 )",-1t"'m! (3.60) The denominator of the first fraction can also be simplified a little, replacing the nwneral 1 in the exponents by mk, then simplifying so (1 _ a1-k)(1 _ �-2k) (1 _ �-("'-1 )k) = (_l j-1 (a(m-1)k _ 1 )(a(m-2)k - 1 ) . . . (� - 1) , (3.61) then from (3.58) we have • C _ (�) 1 m! 1 - a (ak - l )(a2k - 1 ) . . . (amk - l) K(ak, l ) log a ' where of course mk = 1 . 3.3.3 Discussion and Examples We have shown that the SSD in the case g(x) = gx and b(x) = bxl, C 00 (-l)lIak"e-� y(x) = x2 ! (ak - 1 ) ... (ab' - 1 ) ' where a = � and { (k) l _ 1 • (01-1 )(CIS -1 ) ... (11-'-1) K(OI,1) 101 a C = (k) l 1 • (a-1)K(at,a)r(1-t> , t = �, m E 1+ , t # �, m E 1+. 52 (3.62) (3.63) (3.M) Continuity of C In order to show that C is a continuous function of le, we need to prove that lim K(�,a)r(l _ �) = (ak - 1)(a2k - 1); . . (a(lII-l)k - 1 ) K(�, l ) loga. (3.65) 1 __ -1 � m. Now using (3.42) and separating out the one factor which becomes zero as le _ m-I , we write 00 K(ak, a) = (1 - ai-h) n (1 - ai-kit). II=I,llPt (3.66) Similarly, we separate out the factor in r(1 - i) which causes difficulties as le _ m-I , by applying the iterative relationship r(x - 1 ) = ��l m times, giving r 1 _ !. - [_1_] r(l - i + m) ( t ) - m - I (1 - i)(2 - i) . . . «m - 1) - i> " Substituting these expansions into the left hand side of(3.65) we have lim K(ak a)r(l - !.) - lim { r(l - i + m)�I,lle(l - ai-kit) [_l _-_al"7"-_h] } 1 __ -1 ' le - 1 __ -1 (1 - t )(2 - t) . . . «m - 1) -i) m - i (3.67) (3.68) Now the limit of the expression in square brackets can be found by a single application of l'Hopital's rule, (3.69) m Substituting le = m-I into the rest of (3.68) then gives Now lim K(ak, a)r(l - !.) = r (l )�I,ll#,(l - aI-� ) log a 1 ..... 111-1 le (-1 )",-1 (m - 1 )(m - 2) . . . 1 m (3.70) (3.71 ) so after expanding the left-hand product then changing the summation variable in the right-hand product to get it into the form required by (3.42) we get 00 n (1 - aI-� ) = (-1 j-l (� - 1 )(� - 1) . . . (a(lII-l)k - l )K(ak, l ) . (3.72) 1I=I,llPt Substitution ofthis expression into (3.70) now gives (3.65), so C is a continuous function of le as expected. 53 Moments of �he SSD When k is an integer, we can obtain the kth moment of the SSD about the origin, Ilk = fooo xky(x)dx, directly from the integral condition (3.27) as 1 Pt. = a(a - 1 ) " (3.73) An iterative expression for moments beyond the kth is obtainable by applying a Mellin Transform to (3.20): multiplying through by .%"'-1 , where m is an integer, and integrating over all x, remembering that y(x) = x-2Z(x) gives 1000 ,t"-lZ'(x)dx = -ap",+k + a 1000 ak,t"+A:[x-2Z(ax)]dt. (3.74) Integration by parts on the left and the substitution z = ax on the right then leads to (3.75) so we have the iterative equation (m - 1 ) (3.76) provided m > 1. This enables us to deduce all the moments of the SSD once we have found the kth moment by using (3.73) and all other moments up to the (t + 1 )th by integration of (3.63) directly. For example, in the particular case t = 1 , the mean is 1 (3.77) Jl = a(a - 1 ) and the second moment about the origin (Jl2 = fO'ry(x)dx = fooo Z(x)dt) can be obtained by integrating (3.25) term by term to give 1 Jl2 = a2(a - 1 ) log a · We can then apply (3.76) repeatedly to obtain (n - 2)! JI,a = ---........:..-��--� a"( a - 1 ) log a n;j(1 - a-A:) 54 (3.78) , n > 2. (3.79) 1 .2 ....... . 8 x - >- . 6 .4 .2 .5 2 .5 3 x Figure 3.1 : SSD probability density functions y(x) for g(x) = gx and hex) = hxk for a range of values of le. The variables a and a have been fixed at 1 and 2 respectively. An Example Figure 3.1 shows the shape of the SSD function y(x) with a = 1 and a = 2 for a range of values of /c. Consideration of the differential equation formed in the extreme cases shows that as le -- 0, the SSD y(x) approaches the Dirac delta function 6(x), and as /c __ 00 it approaches (Koch and Schaechter [30)) the function y(x) = (a-l)ii ' a < x <.Q { 11 11 o , OtherwISe. 55 (3.80) 3.4 A Minimum Size for Division: The Case b(x) = bxkH(x - Xl ) This case is the same as that discussed in the previous section, except that there is a minimum size Xl > 0 below which cells do not divide. Substituting b(x) = bxkH(x - Xl )' k > 0, (where H is the Heavyside or unit step function and b and Xl are constants) into equation (3.14) we get Z'(X) = xf-l ( -aH(x - XI )Z(X) + a�H(ax - XI)Z(ax» , (3.81 ) where as in Section 3.3 we have made the substitution a = � . In this case we choose to couple (3.81 ) with the integral condition (3.16) as this simplifies to a condition dependent only on Z(x) for X � Xl , namely 100 -2 1 ;q xf Z(x)dt = a(a - 1 r (3.82) It is convenient to solve the equation in three regions starting with the region X > Xl . We let { ZI (X),X � Xl (region 1 ) Z(x) = Z:a (x) , � � X � Xl (region 2) Zs(x),x � � (region 3) where continuity considerations require ZI (XI ) = Z:a(XI ) and Zs(�) = Z:a(�). 3.4.1 Region 1, x � Xl In this region, (3.81 ) becomes (3.20), so the solution is 00 ( -1 )" aD !!!{t ZI (X) = C,,� (a-t - l ) . . . (aG - 1/- 'lb find the constant C, we substitute Z(x) = Zl (x) into (3.82) above, giving c 100 [I ( -1 )"ab' xl-2 -/l�l dx _ 1 ;q ,,=0 (a-t - 1) .. . (aG - 1 ) e - a(a - 1 ) " (3.83) (3.84) (3.85) This time there is no difficulty in immediately interchanging the integral and summation operations as the integration starts at Xl > 0 and therefore avoids the singularity at X = 0 56 when le � 1 . Making the change of variable z = a� in each integral leads to • Hence we can write (le) 1 1 c - -a (a - l)S' where and r(c,x) = f%oo r-le-'ds is an incomplete gamma (unction. In the particular case le = 1 , (3.87) reduces to c = le a(a - 1)S where and El is the exponential integral. 8.4.2 Region 2, � :::; x < Xl (3.86) (3.87) (3.88) (3.89) (3.90) In this region we have H(x - Xl ) = 0 and H( c:u - Xl ) = 1 80 we may immediately integrate (3.81 ) to obtain (3.91 ) This integration can be carried out term by term, using the substitution z = a�, to give 00 (-1 'f [ _* _tl,,*l)A/j �(x) = C,.!O (at - 1) .. . (abt - 1) e - e I , where the constant C is the same as that in the expression for Zl (x) in (3.84). 57 (3.92) 3.4.3 Region 3, x < � In this region the (3.81 ) becomes �(X) = 0 so as we have Z3(�) = ZJ(�) = 0 from (3.92) the solution is simply Za(X) = 0 as expected. 3.4.4 Discussion and EXamples (3.93) (3.94) Combining the results above and noting Z(x) = ry(x), the SSD in the case g(x) = gx and b(x) = bxkH(x - Xl ) is given by (3.95) where (k) l 1 c - -a (a - l)S (3.96) and (3.97) Figure 3.2 shows this piecewise function for y(x) for a range of mjnjmum sizes for cell divisionxl t with k = 2, a = � = 1 , and a = 2. The discontinuity in the slope ofy(x) atx = Xl is a natural consequence of the step change in the probability of cell division b(x) (from 0 to b�) at this size. The form ofy(x) for small values of Xl approaches the form given in Fig. 2.1 for k = 2. 58 1 .2 ......... .8 x ......... >. .6 .4 .2 00 .5 1 .5 2 x Figure 3.2: SSD probability density functions y(x) for g(x) = gx and b(x) = bxkH(x - Xl ), with k = 2, a = 1 , a = 2, and a range of values of Xl 3.5 A Minimum and Maximum Size for Cell Division Consider the case where g(x) = gx and there exists a mjnjmum size for cell division Xl as in the previous section, and there is also a maximum cell size X2 (0 < Xl < X2) such that all cells divide before they reach size X2. This means we are considering non-negative forms of a(x) = W such that a(x) = O,x � Xl , a(x) -- oo, x -- Xi , (3.98) (3.99) and we will assume y(x) = 0 for X > X2. For reasons that will become apparent later we will in fact restrict our attention to forms of a(x) which are integrable over all closed regions [c,d] such that 0 < c < d < X2, but with a non-integrable singularity at X2, 80 1%t a(x)dx = 00 (3.100) for any c < X2. 59 3.5.1 The Solution Method The solution, Z(x), to (3.14) can be constructed as follows. Firstly, let m be the largest integer such that Xl < X2a-IJI. In typical cell populations, m is usually quite small - lets than 3. Tyson and Diekman [61] give the solution for the special case m = 0 and a = 2, and show how in this case solutions may be obtained for forms of g(x) other than g(x) = gx by treating Q (see introduction) as an unknown constant to be found numerically. The method outlined here could be similarly extended to more general forms of g(x). We split Z(x) into regions so that Z(x) = { Zo(X), X � X2 Zk (x) , X2a-k 5 X < X2a-k-1 , 1 < k 5 m + 1 . region: I I I I I m+ I I " • k I 0 I l(x)- �_.(x� . . . z,(x) l.(x) Iz.(X)-O ��F�! x X2), so it follows from (3.102) that (3.103) where C is a constant to be evaluated, andA'(x) = a(x). Ifwe now assume that the solution in region le is known, then (3.102) becomes a simple first-order ODE in Zl+l (X), and the solution is ... -.t l-aI-.t ,..-1 Zl+l (X) = ef. a(.r)drZl(X2a-1) - ca eJ• II(.r'}dI' a(s)Zl(s}df, (3.104) where the limits of the integration have been fixed by the continuity requirement Zl+l (X2a-1) = Zl(X2a-1). Hence we can deduce the solution Zl(X) for all le � 1 , then calculate the constant C from the either of the integral conditions (3.15) or (3.16). The desired SSD can then be simply calculated from y(x) = x-2Z(x). Two points are worth noting in this solution process. Firstly, from (3.103), continuity of Zl with Zo at x = X2 for non-zero C requires that A(X2) = 00, so X2 must be a non-integrable singularity of a(x) as pointed out earlier. The alternative would be to allow a discontinuity in Z(x) at X2, in which case y(x) would be discontinous at x = xa-1 for 0 � le � m + 1 , 80 the argument involving continuity used above would be invalid. Equation (3.103) must clearly be used in place of (3.104) for le = 0 because the first term in (3.104) becomes a product of infinity and zero. Secondly, once the solution has been found in the size region between Xl and X2 using the above process, the solution in the region x < Xl is most simply found by noting that a(x) == 0, so the integral form of the equation (3.17) can be used in the form Z(x) = L«Z a(s)Z(s)df , x < Xl . (3.105) It is clear from this formulation that we will have Z(x) = 0 for x < Xl as expected, with continuity ofZ(x) atx = Xl . 3.5.2 An example As an analytical example of this process, consider the simplest form of a(x) which satisfies the conditions (3.99), namely a(x) = { t-z (3.106) ,X < Xl , 61 which corresponds to a birth rate b(x) of { -'!... b(x) = �-z (3.107) where g is the constant in g(x) = gx. We will choose Xl = a-2x2 exactly, as this is the simplest choice which involves use of all of the equations (3.103),(3.104), and (3.105). In region 1 we then find, using (3.103), In region 2, using (3.104) we obtain x X -1 Z2(x) = CX2(1 - -)[1 - a 10g(1 - -) + a log(1 - a )] . X2 X2 Applying integral condition (3.16) we get f7JJa-1 Zt(x) dx + f7JJ Zl (X) dx = � lzsa-I X(X2 - x) lzsa-1 X(X2 - x) a - I which leads to • 1 C = [(a - l )[(2 + a log(l - a-1 » 10g a + a(S(a-1 ) - S(a-2»]] - , where Se") = I..�1 S· We now use (3.105) to find the solution in the region x < Xl as Zs(x) = CX2a [(1 + a log(1 - a-l » [� - a-a] X2 + (1 - a� )�og(1 - a�) - 1] - (1 - a-2)�og(1 - a-2) - 1 ]] . X2 X2 Now using y(x) = x-2Z(x) we can write the SSD y(x) as y(x) = o � a [(1 + a log(1 - a-I » [� - a-8] +(1 - a�)[log(l - a�) - 1] -(1 - a-2)�og(1 - a-2) - 1]] �(1 - �)[1 - a 10g(1 - �) + alog(l - a-I )] �(1 - �) o 62 ,X2a-2 < X � x2a-l ,X2a-1 < X � X2 ,X > X2· (3.108) (3.109) (3.110) (3.111) (3.112) (3.113) x Figure 3.4: The SSD probability density function y(x) for a(x) = � = (1 - x)-lH(x - 0 .25) and a = 2. where C is given by (3.111 ). Fig. 3.4 shows this function with binary fission (a = 2) and with X2 equal to one size unit. Substitution is straightforward apart from the calculation of S(0.5) and S(0.25) required to find C. Jolley [26] gives S(O 5) = tfJ _ (log2� . 12 2 ' (3.114) but the only way to find S(0.25) seems to be to substitute into the series, which actually converges very rapidly in this case. From the form of (3.113) it is clear that X2 isjust a scaling factor, so the solution for other values of %2 is simply obtained by relabelling the horizontal and vertical axes in Fig. (3.4) as � and X2Y(x) respectively. 63 Cbapter 4 Cell Populations with Age and Size Structure This chapter is a slight modification of a paper published in the Journal of Mathematical Biology (Hall, Wake, and Gandar [20]), with a few general comments having been moved to Chapter 1 and some extra mathematical detail added in places. The development provided in the first three sections of Part I of this chapter, based on the use of cell birth-size rather than cell age, parallels the development based on age in Sections 1 .7 and 1 .8 of Chapter 1 but is more general in that it allows for non-steady size/age distributions. In Part I, we first discuss the case where cells grow without structural constraints, then show how the theory of this general case can be simplified and developed when the cells are constrained to grow in one dimension within structured tissues, as in Section 1 .10 of Chapter 1 . Steady size/age and sizelbirth-size distributions of cens are expressed in terms of the steady cell birth-size distribution, which is shown to be the principal eigenfunction of a Fredholm integral operator. In Part II, we solve some special cases analytically and present a simple numerical method which can be used in analytically intractible cases. 4.1 Part I: An Integral Equation for the Birth-size Distribution 4.Ll Introduction: Size/Age and SizelBirth-size Distributions Consider a population of cells, growing and dividing without mortality so that the total cell population is always increasing. We consider the important attributes of cells to be age, a, and size, I, and assume that the cell population is sufficiently large to be adequately described by a continuous number density n.,(I, a, t) at any time t. Then n.,(I, a, t)dadJ gives the number of cells in the age range a to a + da and size range I to I + dl, and the total cell population at t is N(t) = fO Jo fIt,(l, a, t)dadl. This size/age approach is the one taken in Sinko and Streifer [52] and Dster [38] amongst others, and contrasts with the simpler size-based approach of Koch and Schaechter [30] and Collins and Richmond [7]. In this chapter, cell size may be taken as being synonymous with cell length. We assume that the growth of each individual cell is detenninistic in the sense that given its age a (time since it was formed by fission of a parent cell) and size I, then it is possible to determine its growth rate gG(I,a, t) precisely at any time t. Dispersive growth as discussed in Section 1 .4 is not considered. If cell growth rates are deterministic and positive, models in which all functions and number densities are defined in terms of initial cell size s and current size I are equivalent to size/age models. Using this sizelbirth-size approach, we let n(I, s, t) be the continuous cell number density such that n(I, s, t)dlds gives the number of cells that were born in the size range s to s + ds and lie in the size range I to I + dl at time t, and we write the cell growth rate as gl(l, s, t). Fig. 4.1 shows typical paths followed by a cell in l-a and l-s spaces. In this chapter we choose to work in l-s space because, as we shall show later, steady sizelbirth-size distributions may arise in some situations where the size/age distribution is non-steady. Some simplification is also achieved because individual cells follow straight paths parallel to the I axis in l-s space, whereas paths are generally curved in I-a space (Fig. 4.1). In the cell division process, we allow for possible variation both in the size at which a cell may divide, and in the proportions into which it divides. Firstly, we specify a "hazard rate­ hl(l, s, t) such that hl(l, s, t)dl gives the probability that a cell which was born at size s and which reaches size I at time t will divide before reaching size 1 + dl. Secondly, we introduce a "fractional-size" probability density function (pdf) /(p, l, s, t) such that /(p, l, s, t)dp gives 65 (a) l=s Size, 1 Size, 1 Figure 4.1 : Possible paths followed by a single cell in (a) sizelage space, and (b) sizelbirth­ size space. The solid lines show cell trajectories, the dashed lines show where cell division may take place, and the dotted lines show where new cells may appear. the probability that the birth-size of a daughter cell lies between fraction p and p + tip of the parent cell size, given that the parent cell was born at size s and divides at time t and size I. For all I, s, and t, we require/(p, I, s, t) = 0 for p > 1 andp < 0, and 10 1 /(p, I, s, t)dp = 1 . (4.1 ) If the mean number of cells formed when an adult cell born at size s and dividing at time t and size I is a(l,s, t), then conservation of total cell size during division means that for alI I, s, and t. 10 1 p/(P, I, s, t)dp = [a(l, s, t)t1 (4.2) The special case where all cells divide exactly in half and hi is independent of initial cell size s (which means we write hl(l) only and set/ equal to a Dirae delta function) has been discussed in Chapter 2 of this thesis and in Tyson and DiekmsDD [61], amongst others. In this chapter, we deal with uneven cell division, discussing the theory of the more general case in Part I, then, in Part 11, solving extreme cases where either the division size of a cell is precisely determined by its birth size s, or where hi is completely independent of s. Using the notation established above, number conservation arguments based on the equation (1 .11 ) in Chapter 1 show that the number density, n(l, s, t), must satisfy {J (J (J,n(I, s, t) + mgl(l,s, t)n(l,s, t) = -hl(l,s, t)gl(l,s, t)n(l,s, t). (4.3) 66 for s > 0 and 1 > s, along with an appropriate boundary condition on the line 1 = s. r + I1r ..................................... ········�; ID r ................................. . ...•........... - :11f: Size 1 + 111 Figure 4.2: Formation of new cells along the line 1 = s by division of cells in the region 1 > s. The parallelogram indicates the region containing all cells born in the size range s to s + & in time M. New cells come into existence along 1 = s as a result of cell division in the region 1 > s. Consider the number of new cells born in the small size range s to s + & in the short time interval III (Fig. 4.2). Summing the contributions from cells dividing all over the feasible region 1 > s, we obtain n(s, s , t)gl(s, s, t)M& = f t (hl(/, r, t)gl(/, r, t)fJ) (a(/, r, t)f(7 , I, r, t)�) (n(l, r, t�) . (4.4) 1::3,::0 The first of the bracketed group of terms on the right hand side is the probability that a cell of size 1 and birth-size r will divide in time M, the second is the mean number of daughter cells in the size range s to s + & produced by such a division, and the third is the number of cells in the size range 1 to 1 + l!J with birth-size r to r + M . Dividing through by Ill& and taking limits we obtain the boundary condition in the form of a renewal equation, [00 (I S dl gl(S,S, t)n(s, s, t) = I Jo hl(/, r, t)gll, r, t)a(/, r, tY{ I ' I, r, t)n(/, r, t)drT ' 4.1.2 Steady Sizelbirth-size Distributions We now make the following assumptions: 67 (4.5) • The hazard rate, hi, is time-independent, i.e. hi = hl(/,s). • The fractional-size pdf,/, is dependent only on the proportionp, i.e. / = /(p). It then follows from (4.2) that a is a constant. Then (4.1 ) and (4.2) become lol /(P)dp = 1 (4.6) and lol p/(P,I, s, t)dp = [a(/, s , t)r1 (4.7) respectively. • The growth rate is separable with respect to time, i.e. gl(/, s, t) = g(/,s)r(t) say. We are interested in the situation where the total number of cells in the population is increasing, but the proportion of the population currently in any given size range and born in a given size range remains constant. Under these conditions, the cell number density function n(l, s,t) is separable, and we may write n(l, s, t) = N(t)y(/, s), (4.8) where y(/, s) is a probability density function so loOO 100 y(/, s)dlds = 1 , (4.9) and of course we need y(/, s) � 0, VI,s. The proportion of cells in the size It to 12 which were born in the size range SI to S2 is given by Jkl J: y(/, s)dsdl. We will refer to y(l,s) as the steady sizeJbirth-size distribution. Substituting (4.8) into (4.3), we now obtain dN {) diy(l, s) + N( t)r(t) m(g(l, s)y(/, s)) = -N(t)r(t)hl(l, s)g(/, s)y(/, s), (4.10) and the renewal equation, (4.5), becomes fOO fl s dl g(s,s)y(s, s) = a JI Jo y(/,r)hl(/, r)g(l, r)f(I )drT · Integrating (4.10) over all s and I, substituting g(s,s)y(s, s) from (4.11 ) leads to dN = �(t)r(t) , dl 68 (4.11 ) (4.12) where Q is the functional on y defined by Q = (a - 1 ) 1000 100 y(l, s)h/(l, s)g(l, s)dlds. (4.13) Substituting the expression for dN /tit given in (4.12) into (4.10), we obtain a partial differential equation for y(l,s), 8g(I,�:<" s) = -(hI(l,s)g(l, s) + Q)y(l,.r). (4.14) Thus the problem of finding the steady sizelbirth-size distribution y(l,s) becomes one of solving (4.14) with its boundary condition, (4.11 ), subject to the normalising condition (4.9). Note that (4.13) can be obtained by integrating (4.14) over all I and s then using (4.9) and (4.11 ), so (4.13) is an alternative normalising condition rather than an independent equation. 4.L3 Growth Rate Proportional to Cell Size Direction of growth . .. Figure 4.3: An idealised section of plant tissue growing by elongation in one dimension. Cells are arranged in files and are locked together by common longitudinal walls. When cells divide, new cross-walls (dashed lines) form within the parent cell, and two new cells are formed (i.e. a = 2). We assume there are sufficient files at right angles to the direction of growth to allow use of continuous number density and probability density functions. In this chapter we are interested in steady cell size distributions in plant tissues growing in one dimension only. Fig. 4.3 reproduces Fig 1 .6 from Chapter 1 , showing an idealised 69 section of a root or leaf, elongating in a spatially homogeneous manner such that any two points a unit distance apart separate at a time-varying rate r(t). We take the size of a cell to be its length, I, in the direction of growth. Plants grow symplastically, i.e. with longitudinal files of cells locked together 80 that cells in adjacent files do not slip past each other. This means that all cells in a small region may be assumed to grow at the same relative growth rate r(t), so that g,(/, I, t) = lr(t). (4.15) This assumption is central to all that follows, as it leads to considerable simplification of the more general theory. Now the total length of all cells in a section of tissue as in Fig. 4.3 is lOO 1 to n(/, I, t)dsdl, and the rate of change caused by cell growth is 1000 J/, g,(/, I ,t)n(/, I, t)dsdl. Using (4.8) and (4.15) we then have ! (1000 110' Y(/,I)N(t)dsdl) = 1000 lo' lr(t)y(/, ')N(t)dsdl. (4.16) Interchanging the order of differentiation and integration and substituting for dN / tIt using (4.12) yields simply Q = 1 , (4.17) which is the same as (1 .85) with the constant g == 1 . Substituting this value of Q and g(/, I) = 1 into (4.14), we obtain IOy�t) = - (h,(/, 3)J + 2)y(/,3), (4.18) and the renewal equation (4.11) may be rewritten as SY(I, I) = a 100 t y(/, r)h,(/, r)f(j)drdl. (4.19) Since both (4.18) and (4.19) are homogeneous, any solution of this system will contain an arbitrary multiplicative constant. This constant can be evaluated either using the unit integral condition (4.9), or using(4.13) with Q = 1 , i.e., (a - 1) 1000 100 ly(/, I)h,(/,I)dlds = 1 . (4.20) Following standard reliability or demographic theory, we now introduce a survivor (unction, B,(/, 3), which is the probability that a cell initially of size 1 will survive without division to at least size I. In terms of the hazard rate h" BI(/, I) is given by B,(1, 3) = e-I: It, (tJ,I)da. (4.21 ) 70 If for a given birth size s, there is a maximum ceU size L(s) such that f; h,(CT,S)dCT is finite for I < L(s) but f,(6) h,( CT, S)dCT is infinite, then B(L(s) ,s) = 0, and we set y(/, s) = 0 outside the region s ::; I ::; L(s). The derivative of B" b (I ) = lJB,(/,s) I , s lJI ' (4.22) is a pdf such that b,(/, s)dl gives the probability that a cell initially of size s will divide in the size range I to 1 + dJ. From (4.21) and (4.22) it follows that (4.23) The characteristics of (4.18) are parallel to the I axis, so (4.18) can be integrated directly (cf. Sinko and Streifer [52] where the problem is formulated in terms of cell size and age) to give (4.24) where ; is a function yet to be found, and B,(/,s) is the survivor function given by (4.21 ). 4.L4 Size, size/age, and birth-size distributions We leave to one side for the moment the question of how ;(s) in (4.24) is to be evaluated, and look at ways in which various other distributions may be expressed in terms of ;(s). Integrating (4.24), we obtain the one-dimensional cell size pdf Y(l) = Fo y(/, s)ds), Y(/) = � t s;(s)B,(/, s)ds. (4.25) 1b describe a size/age distribution, we define a pdfy.{l, a, t) such that Y4l(/, a, t)dadl gives the proportion of cells whose ages lie in the range a to a + da and whose sizes lie between I and 1 + dl. In general, the size/age distribution may be non-steady when the sizelbirth-size distribution is steady. The growth rate of an individual cell of size I (and any initial size s) is given by (4.15), so we solve dl/tit = r(t)1 to obtain s = le - t_ ,(t'}dt' • 71 (4.26) Given this relationship between s and a , we must have y.(I, a, t)da = -y(l, s)ds where ds = - r(t - a)le- t_r(t')dt' da, so YIJ(I,a,t) = r(t - a)le- t_ r(t')dt' y(l, le- t_ r(t"�'). (4.27) Using (4.24), we can then obtainy.(I,a, t) in terms of ;(1) as Steady size/age distributions Y.( I, a) are only possible if the growth rate is steady, i.e. when r(t) is a constant, say g. In this case, we obtain (4.29) From (4.28), the size distribution of newborn cells is just ,.(I, O , t) = r(t);(l), (4.30) so we may interpret ;(1) as a multiple of the pdf describing the newborn cell size distribution. Substituting (4.24) into (4.9) and integrating by parts, we have 1000 ;(s)ds - 1000 s;(s) [f bl(�'S) dl] ds = 1 . (4.31 ) Similarly, substituting (4.24) into (4.20) gives (4.32) Adding (4.31) to (4.32), we then obtain the condition (4.33) which will be used to normalise ;(s) where required in Part IT of this chapter. If we now define a birth-size distribution yo(s) such that yo(s)ds gives the probability that a given newborn cell lies in the size range s to s + ds, then as Yo is proportional to ; and Iooo yo(s)ds = 1, we must have Yo(s) = (1 - a-I );(s) . (4.34) 72 4.LIS A Fredholm InteJral Equation for ;(s) and its Unique Solution An integral equation for ;(s) is obtained by substituting the solution (4.24) into the boundary condition (4.19), giving 100 r s 1 �s) = a , 10 r;(r)b,(I, r)f(l) pdrdl, (4.35) where b, is given by (4.22). We then change the order of integration, writing (4.35) in the form. of a homogeneous Fredholm integral equation, ;(/) = a 1000 K(/,s);(s)ds, where the kernel K(l,s) is given by (4.36) (4.37) and we note that!a) = 0 ifr < I. Changing the variable of integration top = f gives as an alternative form for the kernel, S {1 I K(/,s) = 1 10 !(P)b,(p 's)dp. (4.38) We assume that ! and b, are such that this kernel K(/,s) exists and is square-integrable on the appropriate subset of the 1 - s plane. We wish to show that (4.36), with a kernel given by (4.37) or (4.38), has just one independent, non-negative solution. Hochstadt [25] shows that every continuous, strictly positive kernel K(/, s) on a finite rectangle has a largest positive eigenvalue, and that associated with this principal eigenvalue there is only one independent eigenfunction. This eigenfunction is also strictly positive (or strictly negative). Zabreyko et al. [69] show that this result can be extended to non-negative kernels like K(/,s) given by (4.37) or (4.38), provided that the kernels are non-factorable in the appropriate subset of the 1 - s plane. That is, for any ; � O, � 0 and 1 > 0, we must have (4.39) for some N � 0, where K(N) is the Nth iterated kernel as defined in either Hochstadt [25] or Kabreyko et al. [69]. A non-negative kernel satisfying (4.39) must then have a largest positive eigenvalue, and corresponding to this eigenvalue there is just one independent eigenfunction, which is non-negative (or non-positive) but not identically zero. Biologically, (4.39) is equivalent to restricting! and b, to functions such that an ancestral cell of some 73 birth-size can produce, after a finite number of generations, a descendent cell with any other birth-size in the domain considered. A set of precise conditions onf and b" which are sufficient to ensure that the kernel K(I,s) given by (4.38) is both square-integrable and non-factorable, is presented in Section 5.1 . Let the largest positive eigenvalue be 4, and let _(I) be a multiple of the corresponding eigenfunction chosen such that _ is non-negative. Then 1000 K(I,s)_(s)ds = .a;(I), and, using (4.37), it follows that .t 1000 1_(1)dl = 1000 I (1000 K(I,s)_(s)ds) dl = 1000 '1000 s [00 f(�)bI(r,S) �dr;(S)dsdl. Changing the order of integration gives 1000 1000 100 1000 I I dl .t 1;(1)dl = s_(s) b,(r,s) f( -)--drds. o 0 , 0 r r r (4.40) (4.41 ) (4.42) Noting thatf(P) = 0 for p > 1 and using (4.7) and the fact that J,oo b,(r, s)dr = 1 , Vs, we obtain .t 1000 1_(1)dl = a-I 1000 s;(s)ds. As ;(s) is non-negative and not identically zero, we therefore have .t = a-I . (4.43) (4.44) Thus a-I is the principal eigenvalue of the kernel K(I, s) in (4.37) or (4.38), and it follows that the Fredholm integral equation (4.36) has just one independent solution, which is also non-negative (or non-positive) but not identically zero. Hence there will be just one function ;(1) satisfying both (4.36) and (4.33), and as the integral is positive this ;(1) will in fact be non-negative and not identically zero. 4.L6 Discussion We have shown in (4.24) that the steady sizelbirth-size pM, y(l, s), for cells in small regions within one-dimensional plant tissues, can be expressed in terms of ;(s), a multiple 74 of the cell birth-size pdI. The corresponding size/age distribution Ya(/, a, t) is given by (4.28) and the one-dimensional steady size distribution Y(/) given by (4.25). All of these solutions depend on ;(s), which is the eigenfunction of the Fredholm kernel (4.37) or (4.38), normalised by (4.33). In the Part II of this chapter we show that this eigenfunction can be found in some special cases using analytical techniques, and more generally by a simple numerical iteration. However, even without finding the birth-size distribution ;(s), it is possible to make some comments on the general shapes of the functions y(/,s), Y(/), and Ya(/, a, t), given realistic forms of the pdfs/(P) and b,(/,s). Biological considerations dictate that there will be limits on the proportions into which a cell can divide, so the support of/ must be rPt ,P2] with o < Pt < P2 < 1. Similarly, for any initial cell size s, there will be some mjnjmum size L*(s), below which a cell cannot divide, and some maximum size, L(s), by which it will definitely have divided, with s < L*(s) < L(s) (Fig. 4.4). Both L*(s) and L(s) should be 11Wnotonic non-decreasing functions of s, because cells which are born at a larger size should in general also divide at a larger size, but the functions L*(s)ls and L(s)ls should be monotonic decreasing in the sense that cells which are born at a smaller size in general need to grow to a greater multiple of their initial size before dividing. Ifwe let the support of ;(s) be [SmiD,smu], a cell of birth-size Smin must divide as an adult so that the smallest possible size for any of its offspring is also Smin. Thus Smin is the solution of SmiD = PlL *(Smin) , (4.45) and this will be unique because of the monotonicity of L*(s)ls. Similarly, Smu is the unique solution of (4.46) In Fig. 4.4 the regions in which y(/,s) can be non-zero are shaded. In region A there is no cell division, so B,( I, s) = l, and from (4.24) y( I, s) varies as liP along any line parallel to the I axis. Cell division occurs in region B, so y(/,s) drops faster than l iP. there, reaching zero on the line I = L(s). In some situations, it is possible for L*(smin) to be greater than Smu, in which case for L*(smin) > I > Smax, the integral in (4.25) is a constant, so Y(/) obeys an inverse-square law. This occurs irrespective of the precise forms of/(P) or b,(/, s), provided that the condition L *(Smin) > Smu is satisfied. This "inverse square" behaviour in a central region is apparent in published examples of steady cell size pdfs Y(/) where the relative growth rate of cells is roughly constant in this region. (e.g. Collins and Richmond [7]). 75 1 = & ,' 1 = L * (&) 1 = L(&) 8 &- Size, 1 Figure 4.4: Regions in which y(I,s) is non-zero. Cell division takes place in region B only, where I > L*(s). 76 4.2 Part 11: Solving the Integral Equation In Part 11, we solve (4.36) for various choices of/(P) and bl(l, s), using (4.33) to normalise the solution, ;(s). In each case, the sizelbirth-size distribution, y(l, s) and the overall size distribution, yell are then calculated from (4.24) and (4.25) respectively and an example is displayed in grapbical form. 4.2.1 Constant cell division size Consider the simplest possible case with uneven division, where all cells divide at the same size, L. We do not impose restrictions on the form of/(P), so we consider the feasible region for cells to be the triangle shaded in Fig. 4.5. Since all cells divide on the line I = L, bl(l, s) = 6(1 - L), (4.47) where 6 represents the Dirac delta function. Hence the kernel, (4.37), of the homogeneous Fredholm equation becomes and (4.36) becomes a 1 fL ;(1) = Li/(I) 10 s;(s)ds. The integral in (4.49) does not depend on 1, so we may write simply ;(1) = A/(L )' We now substitute this ;(1) into the integral condition, (4.33), to obtain A - 1 - (1 - a-I )L ' so that From (4.24) we then have 1 s s y(l, s) = (1 _ a-I )L ?P(I. )' 77 (4.48) (4.49) (4.50) (4.51) (4.52) (4.53) 1 = 8 Size, l Figure 4.5: The feasible region when all cells divide at the same size, L. A cell born on the line I = s divides on reaching the vertical line I = L. and the one-dimensional steady size distribution is found using (4.25) to be L rf Y(/) = (1 - a-I )P Jo sf(s)ds. (4.54) Fig. 4.6(a) shows y(/, s) in the case where f(P) is a normal distribution with mean 0.5 (so a = 2) and standard deviation 0.1 . Note that if the standard deviation of the normal distribution were chosen to be much larger, care would need to be taken in truncating it so that the support off is limited to [0,1]. The corresponding one-dimensional size distribution Y(/) is shown in Fig. 4.6(b), with the line Y = LIP. shown as a dashed line for comparison. The size/age distribution corresponding to (4.53) may be obtained from (4.28). In the particular case where the growth rate is steady with r(t) = g, (4.29) gives the steady size/age distribution as (4.55) 78 s, which is clearly not independent of s in general. It is useful, therefore, to define Bi(l) = e- 1011,(" )1 , with the corresponding pM Thus we can write bi(l) = -dB;l) . Bi(l) B,(l, s) = Bi(s) ' I > s, 79 (4.56) (4.57) (4.58) (4.59) so for any 1 and s, bj(l) bl(l, s) = Bj(s)H(l - s), (4.60) where H is the unit step function. Because the hazard rate is independent of cell birth­ size, the minimum and maximum sizes for cell division, L*(s) and L(s), (Fig. 4.4) become constants, L* and L respectively. All cells must then divide in the region (L* ,L) (Fig. 4.7(a) ). If the support of! is rPt ,P:a], with 0 < Pl < P2 < I , the feasible region for cells is that shaded in Fig. 4.7(a). With this notation, bj(/) may be specified as any pdfon (L* ,L), then bl(l, s) can be found using (4.60). For s < L*, bl(l, s) is identical to bj(I), 80 bj(l)dl may be interpreted the probability that a cell with birth size less than L * will divide between size 1 and size {at dl. 1=s (b) .. . M .. , � I t:Q 0 PIL" I " L· , L*. (b) The regions used in solving for the eigenfunction ,(1). See text for details. Given this form of bl( i, s), we may write the kernel (4.38) of the Fredholm integral equation as a 1P1 I I KCl,s) = lBj(a) PI !(P)bj(p)H(p - s)dp, (4.61 ) so we obtain from (4.36>, after a change of integration order, ;(/) = 7 [!(P)bj(;) [fo� ��i:�dfl dp, (4.62) where the upper limit of the integral with respect to s has come from Heaviside function in (4.61 ). We now define (I) - f(1) " - roo �d.s ' JO Hf\ij 80 (4.63) so (4.62) may be rewritten as (4.64) Now P < P2 < 1 , so in the inner integral s > I/P2 > I. Hence (4.64) expresses f'(I) in terms of f'(s) Is > I/P2, which suggests that we may solve iteratively for f'(I) in a finite number of steps provided that P2 is strictly less than 1 . 1b provide a starting point for this iteration, we observe from Fig. 4.7(a) that f'(I) = 0 for I > P2L and for 1 < PJ.L*. We divide the region PIL* < 1 � L (Fig. 4.7(b» into m sub-regions, where m is the smallest integer such that �L � PJ.L*, with the kth region RA; as RA; = {I I �� < 1 ��-IL}, 1 � ,, � m. (4.65) The mth region may be "truncated", having a lower bound of PJ.L * instead of IIrL. Then we denote the solution in region RA; by f'1(I). and the solution in the union of regions RI to RA; by f';(I), as shown in Fig. 4.7(b). Starting with region RI we have 9'l(I) = 0, 92(1) a lh 1 = T PI !(P)b;(;)dp, (4.66) (4.67) (4.68) The iteration stops at ",,(I), at which stage we know f'(1) in the entire region (p].L*,P2L). Hence, from (4.63), ;(1) is (4.69) where from the integral condition (4.33). (4.70) We can then use (4.24), (4.25), and (4.28) to find y(l, s), Y(I), and Ya(l, a, t) respectively. 81 The "no overlap" case In the particular case where there is no overlap in size between the size distributions of newborn cells and of parent cells (L* � p� in Fig. 4.7), a much simpler solution can be presented. In this case the inner integral in (4.62) is a constant, Al say, so that (4.62) may be written as a (PI I ;(I) = AI 7 JPI f(P)bj(p )dp. (4.71 ) The constant Al is evaluated using the integral condition (4.33). After simplification using (4.6) and (4.2) we obtain Substituting (4.71 ) into (4.24), we have Bj(l) 11'1 s y(l,s) = Al a p PI f(P)bj(p)dp, and then from (4.25) the steady size pdf is Y(l) = Al aB�l) [Pf(P)(l - Bj(;»)dp. (4.72) (4.73) (4.74) This expression for Y(l) is a particular case of a more general result given in Tyson and Diekmsnn [61]. The general form of the size-age distribution at time t may be found by substituting (4.71) into (4.28), or if the growth rate is steady with r(t) = g we use (4.29) to obtain B*(le-fG) 11'1 le-fG � y(l,a) = Al a I l PI 1(P)bj( p )dp.� (4.75) As an illustrative example which can be solved analytically, we choosef(P) to be the quartic with minima at {In , 0) and (P2, 0), scaled so that (4.6) is satisfied and truncated so that f(P) = ° for P > P2 and for P < Pt , f(P) - 30 CR2 - p) 2(P - Pt )2H(p - P)H(P - Pt ) (4.76) � _ pt )5 2 • For bj(l) we choose a quartic similarly truncated at minima at (L* , 0) and (L, O), bj(l) = 30 (L - 1'f(I - L*)2H(L _ l)H(I _ L*) (L - L*)" . 82 (4.77) (a) (b) Figure 4.8: (a) The steady sizelbirth-size pdfy(l,s) and (b) the steady size pdf Y(l) where J(P) is a quartic truncated at minima at P1 = 0.35 and P2 = 0.65, and bi is also a quartic, truncated at minima at L* = O.7L and L. The factor 30 comes from the requirement that J and bi are both probability density functions, so each must �ve unit integral. In the example shown in Fig. 4.8, we have chosen P1 = 0.35, P2 = 0.65, and L* = 0.7L, 80 the condition L* < P2L is satisfied. We have then used equations (4.73) and (4.74) to obtain y(l,s) and Y(l) as shown in Fig. 4.8 (a) and (b) respectively. If we chose different forms ofJ(P) and bi(l), the two-dimensional sizeibirth-size distribution y(l,s) as in Fig. 4.8(a) would still consist of a set of "inverse square" lines along the characteristics up to 1 = L *, then a more rapid decrease down to y(L, s) = O. The dashed line in Fig. 4.8(b) is the line Y = A1L/P, which closely matches the solution Y(L) on a region larger thanjust f.R2L,L* ). 4.2.3 Cell Division Size Precisely Determined by Birth-size The Simplifled Kernel In the previous section, we considered the case in which the hazard rate is completely independent of initial cell size. Here we consider the opposite extreme, where cell division size is precisely determined by initial size, i.e. where all cells born at one size will divide at the same larger size. This means that the minimum and maximium sizes for cell division, L*(s) and L(s) (Fig. 4.4), are made to coincide, 80 we set 8,(/, s) = 1 -H(l - L(s)), 83 (4.78) where H is the unit step function, and b/(/,s) = 6(1 - L(s» , (4.79) where 6 is the Dirac delta function. Thus, cells born at size 3 grow without division until they reach 1 = L(s) and then divide. Under these conditions, we can substitute (4.79) into the kernel of the Fredholm integral (4.37) to obtain so the equation for ;(/) is 3 1 K(/,3) = (L(3» 1I(L(S» ' (4.80) (4.81 ) Once again, ; must satisfy the integral condition (4.33), and we note that/(ljL(s» = 0 for L(s) < I. Despite its apparent simplicity (4.81 ) is not easy to deal with in general. In the follOwing, we choose a realistic form for L(s) together with the simplest possible form of/(P) in order to obtain a series solution. A Special Case: A Power Law L(s), and a Uniform Distribution for f(P) We consider here the special case L(s) = Ll-TST, 0 < r < 1 , (4.82) where L is now a parameter with units oflengtb. We make 0 < r < 1 so that L(s) satisfies the requirements set out in section 4.1 .6, that L(s) is monotonic non-decreasing on (O,L) and L(s)js is monotonic decreasing. The kernel, (4.80), of the Fredholm integral equation then becomes ;-2T I K(l,s) = L2(1-r1(o-rsT ). For /(P), we choose the uniform distribution on (0, 2a-1 ), f(P) = ;H(2a-1 -p), 84 (4.83) (4.84) so that the mean of the distribution is a-I as required by (4.7). This choice of / is unrealistic, but it allows us to illustrate the use of analytical solution techniques. Note that for this form of/(P) it is necessary that a � 2, to prevent formation of daughter cells larger than parent cells. Size- I Figure 4.9: The feasible region for cells when L(s) = Ll -rsr, and/(p) = H(2a-1 - p)a/2. For a = 2, Smu = L, but in this figure a > 2, so from (4.85) Smu < L. For this form of/(P), PI = 0, so from (4.45) we have Smin = O. Substituting (4.82) and P2 = 2a-1 into (4.46), we obtain (4.85) which gives the feasible region for cells shaded in Fig. 4.9. Note that if a = 2, we would have Smu = L 80 the shading would continue up to the point where the line I = S meets the line I = Ll-rsr. In the shaded region, substitution of (4.83) into (4.36) gives a2 [- 2 I sl-2r ;(1) = 2" 0 H(;; - Ll-rsr )L2(I-r) ;(s)ds which can be written as 2 1'-;(l) = 2(I-r) ( 1 )1 ;- 2r ;(s)ds. Smu iiiiU " - Differentiating both sides of (4.87) we obtain ;' (I) =-(_1 ) � 2.f2mu; ( _I ) t smu) . Smu r{3 Smu (4.86) (4.87) (4.88) '1b turn this into a simpler dimensionless form, we make the substitution I = smue-"', with M(m) = ;(1), so M'(m) = ;'(I).!!!.. = -smue-"';'(I) din 85 (4.89) and (4.90) Substitution into (4.87) then gives after some simplification the functional differential equation (4.91) Earlier success applying Laplace transform techniques to functional differential equations involving terms like M(r-1m) (see Chapter 2) encouraged us to try this approach here. Taking Laplace transforms of both sides of(4.91), we obtain pM(P) = 2M«(P - 2)r + 2), which we write in the iterative form Iterating N + 1 times, we get - 2 -M(P) = pM«(P - 2)r + 2). (4.92) (4.93) M(P) = �[(� - 1 )r + 1][(� - 1 )�2 + 1] ... [(i - 1 ),-N + I]M«(P - 2)?+1 + 2), (4.94) which taking limits N - 00 gives (4.95) This limit exists, because r < 1 . Expressing this product as a sum of partial fractions, where N 1 N QII(2(1 - r-II» Ill] [(i - 1 )ya + 1 ] = II� [(i - 1 )ya + 1 ) ' (4.96) 1 QII(P) = � [(� - 1 )r + 1 )] ... [(� _ 1 )ya-1 + 1][(� - 1 )ya+l + 1 ] .. . [(j - 1 ),-N + 1] ' (4.97) so that QII(2(1 - r-II» = (1 _ r-II) ••• (l _ r-1�(1 _ r) . . . (l _ rN-II) · (4.98) Hence in the limit as N - 00, we have _It (-1 )" 1 QII(2(1 - r » = (r-II - 1 ) . .. (r-1 - 1 ) K(r) , 86 (4.99) where 00 K(y) = IT (1 - r"). (4.100) 11=1 The constant K referred to in Chapter 2 is in fact K(a), where the function K is given by (4.100), 80 for a discussion of this function, including alternative expressions and series, the reader is referred to Section (2.4.1). We now write (4.95) as M - M(2) � ( -1)" [ 1 ] (P) K(y) to (y-l -1 ) . . . (y-II - 1 ) (i - 1 )r" + 1 (4.101) so - 2M(2) 00 (-1 )"y-II r _ 1 ] M(P) = K(y) II� (y-l - 1 ) . . . (y-II - 1) LP + 2(y-" - 1 ) , (4.102) where we are using the convention that the empty product in the denominator is given the value 1 when n = O. Taking inverse Laplace transforms term by term we obtain 00 ( 1 )" -11 M(m) - A � - y e-2m(r---1 ) - to (y-l - 1 ) . .. (y-II - 1 ) (4.103) where A is a constant to be evaluated from (4.33). Reverting to the original coordinates, we substitute e-Ift = I/smax and M(m) = ;(1), giving 00 (_l )"y-" ( I )2(r---1) ;(1) = A II� (y-l - 1 ) .. . (y-II - 1 ) Smax . (4.104) 1b find the unknown constant A, we apply the integral condition, (4.33), using Smu as the upper limit of integration. This gives -1 00 (-l )"r-1I 1 [ ( ) ] -1 A = (1 - a )smax 1 (y-l - 1 ) . .. (r-II - 1) 2y-1I - 1 (4.105) The sizelbirth-size distribution in the region shaded in Fig. 4.9 is found by substituting of (4.104) into (4.24), setting BI(I, s) = 1 so s 00 (_l )"y-1 ( S )2(r---1 ) y(l,s) = A P. II� (y-l - 1 ) . . . (y-II _ 1 ) Smu . 87 (4.106) (a) (b) Figure 4.10: (a) The steady sizelbirth-size pdfy(/,s) and (b) the steady size pdfY(/) when f(P) is the uniform distribution on (0, 1 ) and L(s) = Ll-rsr with r = 0.25. Fig. 4.10(a) shows the steady sizelbirth-size distribution ](/,s) obtained with a = 2 and the arbitrary r = 0.25. For this value of r, the series in (4.10ijconverges very rapidly, even for s = Smax. The singu1aJ;ity at (0, 0) appears as a consequence ofchoosingf(P) as a uniform distribution extending top = O. As I, s -- 0 along the line I = s, y(/,s) -- 00, but as I, s -- 0 along the line I = Ll-rsr, y(/,s) -- 0 . 1b find the steady size distribution Y(/), we substitute (4.104) into (4.25) and integrate term by term between a lower limit of s = P./rLl-l/r and an upper limit of I (for 1 < smax) or Smax (for I � Smax). This gives where Y(/) _ As2max � (-1 r F (I) - 2P ,..;0 (r-1 - 1 ) . . . (r-" - 1 ) ,. (4.107) (4.108) Fig. 4.1O(b) shows the form of Y(/) in the case where a = 2 (so Smax = L) and r = 0 .25. This one dimensional size distribution has a finite non-zero value at I = 0, and the series in (4.107) converges rapidly for all values of I. The general form of the corresponding size/age distribution can be found by substituting (4.104) into (4.28). In the particular case where the growth rate is constant, with r(t) = g, 88 we obtain 00 (-I )"r-1I ( I ) 2(r---1 ) y(l, a) = Ae-2Ia L ( -1 - 1 ) ( -11 _ 1 ) -e-IG 'j ' 11=0 r . . . r Smu: 4.2.4 Numerical Solution , (4.109) In earlier sections, we demonstrated methods for finding solutions to the homogeneous Fredholm integral equation (4.36) with a kernel given by (4.37) or (4.38), in some special cases, and with the solution, ;(1), normalised so that it satisfied the integral condition (4.33). In the general case, the eigenfunction ;(1) corresponding to the principal eigenvalue a-1 must be found numerically. Here we scale all cell sizes so that Smu � 1 , then use the straightforward iteration ,(0) (1) = 1 , O < I � I , ,(11+1) (1) = a 101 K(I,s),(II)(s)tU, (4.110) (4.111 ) to estimate ,(1) = limll_OO ,(11)(1). Once such an eigenfunction is found, we can normalise using (4.33) to obtain ;(1) and, as before, apply (4.24), (4.25), and (4.29) to find y(l,s), Y(I), and Yil(l, a). Figs. 4.11 and 4.12 illustrate two examples of solutions using this technique. In each case we obtained convergence to within four decimal places in ,(1) within eight integral iterations. Varying ,(0)(1) to other forms such as triangular or normal distributions did not affect the final solution ,(1), and actually improved the rate of convergence. Figs. 4.11(a) and (b) show the steady size/birtb-size and steady size distributions in the case where cell division size is precisely determined by cell birth-size, as in Section 4.2.3, with the cell division size given by (4.82) and the kernel given by (4.83), but withf(P) as a normal distribution, (ra)1 1 - f(P) = � e ." 21r� (4.112) with mean � = 0.5 (so that from (4.2), a = 2) and standard deviation � = 0.1 . This is more realistic than the uniform distribution for f used in Section 4.2.3. We choose r = 0.25, so the feasible region for cells is the same as it was for the solution given Fig. 4.10(a). AB there is almost no overlap between the size distribution of dividing cells and the size distribution ofnewbom cells, th� dashed line Y = Lfo s;(s)ds]/p. in Fig. 4.11(b) touches the line Y = Y(I). 89 (a) (b) Si2Ie, l Figure 4.1 1 : (a) The steady size/birth-size pdf, y(l, s), and (b) the steady size pdf, Y(I), where/(p) is given by (4.112), and all cells divide on the line I = L(s) with L(s) given by (4.82). The dashed line in (b) is the line Y = ut s;(s)ds]/(l. (a) (b) Size, 1 Figure 4.1 2: (a) The steady sizelbirth-size pdfy(l, s) and (b), the steady size pdf Y(I) where /(P) is given by (4.112), and the cell division size is variable, with b/(I, s) given by (4.11 3). The dashed line in (b) is once again the line Y = ucf s;(s)ds]/(l. 90 As a second example, we again choose f(P) to be the normal distribution (4.11 2) with J.I{ = 0.5 and � = 0.1 , but rather than force a cell of initial size s to divide exactly at size L(s) given by (4.82) we allow some variation about this. We do this by setting b,(/, s), for any s, to be a normal distribution with mean L(s) and standard deviation dependent on s, 80 that 1 (I-y.� b,(l, s) = ..f2i e -I(�.} ) 2#0),(s) We choose a,,(s), the standard deviation of this distribution, to be ( ) � (L(s) - s)(L - L(s» 0), s = Jb (L - s) (4.113) (4.114) where/b is a constant. Fig. 4.1 2(a) and (b) show the steady sizelbirth-size and steady size distributions respectively when L(s) is given by (4.82) with r = 0 .25 as before, and with fb = 0.3 in (4.114). 4.2.5 Discussion In all the models described in this chapter, we have assumed the following: • (Is), - The cell population is spatially homogeneous, so we may ignore the location of cells in space. • (I,) -Cell growth rate is proportional to current cell size, but may vary with time, 80 g(l,s, t) = r(t)/. • (18) - The probability that a cell initially of size s will survive without division to size 1 is independent of time, i.e. B, = B,(l, s). • (//) - The pdf describing the proportions into which a parent cell may divide is independent of both time and current or initial size of the parent cell, i.e. / = /(p). If we follow a cohort of cells as in Fig. 4.3 as it moves along a leaf or root, it passes first through the meristem where both cell division and cell growth are occuring, then through a growth-only region, then finally passes into mature tissue. If we choose a cohort which covers just a short length of root or leaf, it is possible to satisfy assumption (Is) fairly closely, 80 as discussed in Section 4.1 .3, symplastic growth assures us that assumption (I,) 91 also holds. The time-invariance implied by assumptions (IB) and (If) will hold only in parts of the meristem in which B/ and! do not change significantly with spatial position, but provided that the rate of change of B/ and! with position is slow compared with the mean distance covered by a cell during its lifetime, steady size and sizelbirth-size distributions calculated using the models described here should be close to the size and sizelbirth-size distributions observed. Assumption (//) has been made for the sake of simplicity, but there is little data available to support or contradict this. Of course, the models here may be applied to cell growth in situations other than one­ dimensional plant tissues, but we have not addressed this question specifically. We have taken cell size to be synonymous with cell length, but it could equally well be any other property which is conserved when cells divide, such as cell mass, volume, or DNA content. Some types of bacterial cells have been observed to grow approximately exponentially in size (for example, see Schaechter et al. [48] or Bell and Anderson [6]), so assumption (I,) may be made in these situations. Plant tissue growth is symplastic whether the cells are growing in one, two, or three dimensions, so in all these cases, provided that we consider a region of tissue in which (Is) holds, we could expect (I,) to hold too. Our specific choices of forms for B/ and ! have been guided in earlier sections by a desire for simplification, but it is interesting to note that the birth-size independence postulated in Section 4.2.2 has been assumed in a number of papers on bacterial cell size distributions (Koch and Schaechter [30], Collins and Richmond [7], Tyson and DiekmaDD [61]). This assumption is made because it is very difficult to measure cell birth-size (or alternatively cell age), so simple one-dimensional models for Y(l) are the only ones which can be compared easily with existing data. However, a comparison of Figs. 4.8(b), 4.11 (b), and 4.12(b) shows that very similar shapes for Y(l) may be obtained from what are fundamentally quite different models, so one should be wary asserting that a particular model is correct on the grounds of the predicted one-dimensional size distributions. The two-dimensional sizelbirth-size distributions shown in Figs. 4.8(a), 4.11(a), and 4.12(a) show greater variation, so these (or size/age distributions) might be a better basis for comparison of models if measurement difficulties could be overcome. It is unclear whether the most general type of model with distributions like those in Fig. 4.12 provides a significant improvement in accuracy over more deterministic models with distributions like those in Fig. 4.11. The "power law" used for L(s) in Sections 4.2.3 and 4.2.4 is the simplest form which satisfies the requirements for "realism" outlined in Section 4.1.6, but of course other forms 92 are also possible. If we look at Figs. 4.11 or 4.12, we see that only the central section of the line I = L(s) actually affects the solution, so the precise form of L(s) over all s is unimportant. A normal distribution for f(P), which we used in the examples shown in Figs. 4.6, 4.11 , and 4.12, seems a reasonable assumption, but its standard deviation must be chosen to be sufficiently small so that fJ f(P)dp = 1 is a reasonable approximation. While there is little experimental evidence for any particular choice off, the assumption of a normal distribution is certainly more realistic than the assumption that all cells divide exactly in half, as in the models of Hall and Wake [19] and Bell and Anderson [6]. Tyson and Diekmann [61] point out that in the case of exponential growth of cell size and division of cells exactly in half, the steady size distribution is not asymptotically stable, whereas Heijmans [23] shows that if cells are allowed to divide unevenly, the steady size distribution in the birth-size-independent case is approached asymptotically. 93 Chapter 5 Existence" Uniqueness, and Stability of SSDs The objective of this chapter is to examjne the existence, uniqueness, and stability of the steady sizelbirth-size distributions discussed in Chapter 4, where the growth rate of cells is given by g,(/) = r(t)/. 5.1 Existence and Uniqueness Consider the Fredholm integral equation (4.36) for the birth-size distribution ;(/) in Chapter 4, with the kernel given by (4.37) or (4.38). In Section 4.1 .5 we noted that there is just one non-negative solution ;(/) which satisfies both the integral equation (4.36) and the normalising condition (4.33), provided that the kernel is both square-integrable and non-factorable. Let the support off be -r(f) = (pt ,P2]. where 0 � 1'1 < P2 � 1 , as in Section (4.1.6). We show here that the kernel will be square-integrable and non-factorable provided that: • (Ai) There exist continuous functions L*(s) � s and L(s) � L*(s) such that b,(/, s) > O forL*(s) < 1 < L(s) and b,(/,s) = Oelsewhere, andofcourse.rt. 0 such that (5.1 ) • (..4".) There exists a mjnjmum cell birth-size Smin, and a maximum cell birth-size Smax, such that PtL * (Smin) = Smin, (5.2) and (5.3) with 0 $ Smin < Smax. Further, we require that for all S E (Smin,smax), we must have p?,L(s) $ Smax and PtL*(s) � Smin. This enables us to restrict the domain of the integral operator to functions ;(s) such that the support of; is contained in [Smin,Smax], that is (5.4) • (.A.r) For S E (Smin,smax), we require that PtL*(s) < S and p2L(s) > s. In fact, in order to prove that the kernel is non-factorable, we need to strengthen these inequalities, and require that for all S E (Smin,smax), "lE > 0, 36 > 0 : Vs � Smin + E, PtL*(s) $ (1 - 6)s, (5.5) and "lE > 0, 36 > 0 : Vs $ Smax - E, P2L(S) � (1 + 6)s. (5.6) In the particular case where Smin = 0, which from (5.2) can occur either because Pt = 0 or L*(s) = s, one extra condition is required in order to ensure that the kernel is square­ integrable: • (Ao) If Smin = 0 then there must exist constants Mb > 0 and u < 1 such that Vs E (0, smax], (5.7) One consequence of (5.7) is that the expected value of the inverse square of cell division size, J,L(I) �dr, must grow more slowly than Mb/s2 as S .... o. 95 5.LI Square-integrability of the Kernel Consider the kernel (4.37), which as/(P) = 0 for p > 1 may be written in the form K(I,s) = s l°O /(� ) b,(r, s) fir. (5.8) mu(l .. ) r r2 From (5.1 ) we have so as K(I,s) � sM, foo b,(r, s) dr, lmu(l .. ) r2 1 1 _2 < ( )'" r � max(l, s), ,- - max I, s • and b,(r,s) is a probability density function in r for all s, we have sM, K(I, s) � max(l,s)2 . Hence provided that Smtn > 0 we have so sm.uM, K(I,s) � �. ' JDlD 1"1IIU [- (sm.u - s_,_)2s! MJ IK(I,s)12dldf � IIIUI mu f , "miD "miD .r'min and the kernel is square-integrable as required. If Smin = 0, then using (5.7) and (5.9), we find that the kernel must satisfy (5.9) (5.10) (5.11) (5.12) (5.13) K(I s) < sM,Mb ( ) , - scr [max(l, s)]C7 ' 5.14 so 1o"mu fo"lIIU 1K(I, s) 12dldf � M!Mf 1o"mu J:- [m:;�:)]2cr dlds. We now split the inner integral into two parts, I > s and I < s, so 1o"mu 1o"mu IK(I, s) 12dldf � M!Mf [10'-- J: s2-4C7dlds + J:DIaZ l"mu s2-2crr2crdlds] , then changing the order of integration in the second term we obtain 1o"m- 1o"mu IK(/,s) l:ldldf � M!M: [10"- s2�C7 J: dlds + J:DIaZ r2C7 1o' s2-2C7dfdl] . As eT < 1 we can carry out the integrations to obtain {'O ('mu :I MJ�S!!!;C7) (2 - CT) 10 10 IK(I, s) 1 dlds � 2(1 - CT)(3 - 2eT) , so the kernel (5.8) is square-integrable as required. 96 (5.15) (5.16) (5.17) . (5.18) s.u Non-factorabillty of the Kernel Zabreyko et al [69] define a square-integrable kernel K(/, s) on a bounded domain 0 to be non-factorable if and only if for any square-integrable function �(/) � 0, �(/) � 0, which is not identically zero, and any 1 e O, there exists an iteration K(N)(/, s) of the kernel, K(N)(I,s) = la la · · ·la K(/, sdK(st , S2 ) . . . K(SN_t , s)dsN_l . . . ds2dst , (5.19) such that In K(N) (/,s);(s)ds > O. (5.20) Heijmans [24] uses the term TWn-8upporting instead ofnon-factorable. Given a non-negative ;(s) such that ;(s) > 0 for s e (SI ,S.,) C (Smin , smax), and l e (Smin, smax), then if 1 E (SI,S.,), (5.20) is satisfied for all N � O. Otherwise, we must have either 1 e (Smin,s,] or 1 E [S." sma). Here we consider only the first of these cases, as the proof in the second case is almost identical. Given 1 E (Smin ,SI], choose £ = 1 - Smin and a corresponding 6 from (5.5). Becausef(P) > 0 for P E 0 for L(s) < r < L*(s), the kernel (5.8) must satisfy K(l', s) > 0, Pl L*(s) < l' < P'JL(s). Hence given ;(s) > 0, S, < s < s." we must have 1'-K(l', s);(s)ds > 0, PlL*(sl) < l' < P'JL(s.,), '- so applying (5.5) and (5.6), we have r:: K(l', s);(s)ds > 0, (1 - 6)s1 < r < s.,. Now we choose the smallest integer N such that that is the smallest positive integer such that log(f) N > log(l � 6) . Iterating the argument leading to (5.23) N times then gives 1'- K(N)(l',s)�(s)ds > 0, (1 - 6,/,sl < l' < s." '- 97 (5.21 ) (5.22) (5.23) (5.24) (5.25) (5.26) so from (5.24) we have rma K(N) (I, s);(s)ds > o. l'fIdn (5.27) A very similar argument can be applied if l E [S." ,,_), 80 we conclude that the kernel (5.8) is non-factorable. 98 5.2 The Time-dependent Problem 5.2.1 Formulation and Simplification We consider here the time-dependent problem giving rise to the steady sizelbirth-size distributions discussed in Chapter 4. The time-dependent number density is the solution of (4.3), {} {} {}l(l, s, t) + mg,(I, s, t)n(l, s, t) = -h,(l, s, t)g,(l, s, t)n(l, s, t), subject to an initial condition n(l, s, O ) = noel, s), and the boundary or renewal condition given by (4.5), [00 {' s dl gl(S, S, t)n(s, s, t) = 6 10 ht(l, r, t)gl(l, r, t)a(l, r, t)f(l ' l, r, t)n(l, r, t)drZ . (5.28) (5.29) (5.30) The notation used here is exactly that used in Chapter 4. As in that chapter, we consider the special case arising when and g,(I, s, t) = r(t)l, "'(I, s, t) = ht(l, s) , f(P, l, r, t) = f(P). The probability distributionf(p) must have support -r(f) � [0 , 1] and satisfy 10 1 f(P)dp = 1 , and 101 pf(P)tlp = a-I , and a must be a constant. As in Chapter 4, we define B,(l, s) = e- t. II(r)dr, 99 (5.31 ) (5.32) (5.33) (5.34) (5.35) (5.36) and the probability density function oB, b,(/,s) = - m = h,(/,s) B,(/, s). Making the appropriate substitutions into (5.28) gives the time-dependent PDE o 0 otn(/,s, t) + m[r(t)ln(/,s, t)] = -h,(/,s)r(t)ln(/,s, t). with a boundary condition given from (5.30) by sn(s, s, t) = 100 10' h,(/, r)a!(l)n(/, r, t)drdl, and an initial condition given by (5.29). (5.37) (5.38) (5.39) In order to avoid the explicit dependence of the PDE (5.38) on time t, we replace t by the dimensionless time-like variable t' = R(t), (5.40) where R(t) = lot r(-r)d-r. (5.41 ) Note that if r(t) = 1 then t' = t. Also, for any r(t), t' is the log of the factor by which the total biomass has grown since t = O. We can ensure that t' is positive for t > 0 and that (5.40) is invertible, ifwe require that r(t) � rmin > 0, where rmin is a constant. From (5.40) we have on on ot = r(t) lJt" so substituting (5.40) into (5.38) we have on on Ot + 1 {)l = -(IJa,(/,s) + 1 )n(l, s, t) , (5.42) (5.43) (5.44) where we have replaced t' by t in order to avoid excessive notational compexity in what follows. Note that the boundary condition (5.39) and the initial condition (5.29) remain unchanged despite the change in the nature of the time variable t. Equation (5.40) has transformed (5.38) into (5.44), which represents the case where r(t) = 1, so from (5.31) the growth rate of cells may be regarded as just g,(I, s, t) = I. (5.45) In the following, we consider just this special case, on the understanding that we can transform to the more general case using the inverse of (5.40). 100 5.2.2 Solution using Characteristics 1b integrate along characteristics, we choose as new variables in (5.44) 11 = Ie-t and -r = I, so I = 11e�, to obtain lJn lJ-r = -( 11e�h,( 11e�,s) + 1 )n. (5.46) This first order PDE has the solution (5.47) where \If is an "arbitrary" function to be determined by either the boundary condition (5.39) or the initial condition (5.29), whichever is appropriate. Reverting to the original variables, we obtain (5.48) We now define a "birth rate" function cz,(I,S) such that cz,(I, s)dstlt gives the number of cells born in a short time tit in the size range s to s + tb. Then cz,(I,S) = = = = dl n(s,s, t) dtl,=.r n(s,s, l)g,(S, s, I) sn(s,s, l) se-t ",(se-t ,s), because from (5.45) we have g,(S,S,I) = .1, and from (5.36), B,(s,s) = 1 . Hence so (5.48) may be written as (I ) _ �(- IOg(f), s) '" , s - I B,(l s) I n(l,s, l) = z' cz,(t - log(s)'S). 5.2.3 The Region In1luenced Directly by the Initial Condition For the region 101 (5.49) (5.50) (5.51) (5.52) we substitute (5.51 ) into the initial condition (5.29) which gives 80 and I I �(- log(;), s) = Bl(I,S)no(I,S) se-I �(t, S) = B ( 1 )no(se -', s), I se- ,s I le-I �(t - log(;),s) = Bl(le-l, s)no(le-' , S). Substituting this into (5.51 ) then gives (I ) BI(I, s) -I (le-' ) n ,s, t = Bl(le-' , S) e no ,s . (5.53) (5.54) (5.55) (5.56) Note that the fraction at the start of this expression cannot exceed 1 , 80 provided that no is bounded, n(l,s, t) must decay exponentially in the supremum norm in this region. 5.2.4 The Region Influenced by the Boundary Condition For the region t > log (i) we substitute (5.51 ) into the boundary condition (5.39) which leads to tx) (' s I dl �(t,s) = 10 10 at(l)bl(l, r)�(t - log (,J, r)drT' (5.57) (5.58) Changing the order of integration, then making the change of variable a = log( �) in the inner integral, gives 100 100 S �(t, s) = at( -e-)b,(rtf, r)tl>(t - a, r)dadr. o 0 r (5.59) We need to split the inner integral in (5.59) into two parts, with separate integrals for the regions a > t and a < t. In the region a > t, �(t - a, r) is given by (5.54) as rtf-I �(t - a, r) = B ( I )no(rtf -' , r) ,a > t. I re"- , r (5.60) 102 Hence (5.59) becomes where 1000 10 ' S tl»(t,s) = af( -e-)b,(rt!' , r)4l(t - a, r)dadr + �O(t,S), o 0 r 1000 100 S rtfl-I �o(t, s) = a/(-e-lJ)b,(rt!', r)B ( -I )no(rell-', r)dadr. o 1 r I rt!' , r The form of�(t,s) 8S t - 00 is discussed in the following section. 103 (5.61 ) (5.62) 5.3 Stability of the SizelBirth-size Distributions We are interested in the properties ofn(l,s, t), and therefore ofcz,(t,s), as t --+ 00. In order to show that that the steady sizelbirth-size distributions calculated in Chapter 4 are stable, we need to show that for large t, we may write for some ;(s), cz,(t,s) = eI( ;(s) + 0(1 » . (5.63) Because we .have transformed to a dimensionless time variable, we anticipate that the coefficient of t in the exponent is simply one. If we can show that cz,(t, s) takes the form (5.63), then as in Chapter 4 we may interpret ;(s) as the steady size distribution ofnewbom cells. In order to establish (5.63), we follow the approach of Heijmans [24], and look at the properties of a family of operators on the Laplace transform of cz, rather than study the properties of the integral operator in (5.61 ) directly. All of the conditions (Ai), (At), (.A...), (.A.r), and ifnecessary (..40), specified in Section (5.1 ), are assumed to hold throughout this section. We also assume that the initial cell siZ&'i>irth-size distribution, no(l, s), is zero outside the region Smin < s < SJJW[, s < 1 < L(s). 5.3.1 The Laplace Transform of the Operator We now proceed to take Laplace transforms of both sides of equation (5.61 ), denoting the Laplace transform of cz, by � and the Laplace transform of � by cito. Then we have (5.64) The Laplace transform ofcz,o(t, s), �o(A., s), exists and is analytic for A. > 0 provided that no is bounded and has as its support the "feasible region" shaded in Fig. 4.4. Now we make the change of variable a = t - a, giving �(A., s) = rx> (00 (t af(�e-(t-'»)bl(re(t-/) , r)cz,(a', r)e-Alda'dtdr + �o(A., s), (5.65) Jo Jo Jo r then swap the order of the two inner integrals and change variables again using t' = t -a, so that (5.66) 104 Now replace the dummy variables a' and t' by t and a respectively, then change the order of the inner integrals again to obtain cP(A., s) = 1000 [1000 a!(;e-tl)bl(rtl' , r)e-Aada] cP(A., r)dr + cPo(A., s). (5.67) This is in the usual form of an inhomogenous Fredholm integral equation, cP(A., I) = 1000 k,,(I, s)cP(A.,s)ds + cPo(A.,l), (5.68) where (5.69) Alternatively, we can substitute r = sell to obtain (5.70) More formally, using operator terminology we may write (5.71 ) where k(A.) is the operator defined by (5.72) Inverting (5.71 ), we have cP(A.) = (/ - k(A.))-lci»o(.t), (5.73) where I is the identity operator. We define an operator R(.t) = (I -k(.t))-l . (5.74) Now in order to evaluate the inverse Laplace transform of cP( A.) precisely, we would want to be able to find all the singularities of R(.t) and their orders. However, as we are interested only in that part of �(t,s) which grows most rapidly at large time, it is sufficient to find just the .t with the largest real part for which R(.t) is singular. It is clear that R(.t) can be singular only if k(A.) has an eigenvalue of one, 80 we now proceed to find the maximum value of A. for which this is possible. 105 5.8.2 Eigenvalues of the Operator K(A.) .. . " . . � ", I t' • d • .. . " "'. The kernel of the operator k given by (5.70) and the kernel K given by (5.8) both have the same support, so the argument used in Section (5.1 .2) to show that the kernel K(/,s) is non-factorable is equally applicable to the kernel t..t. Hence the kernel tA. is non-factorable for alLt > O. The kernel kA. is clearly square-integrable for ,\ � 1 , as in this case we have kA.(/,s) $ aK(I, s), Vs, 1 E [Smin,smaxJ, and K(/, s) has been shown to be square-integrable in section (5.1 .1 ). For '\ < 1 , the situation is somewhat more complex. Provided that Smin > 0, we can use arguments as in Section 5.1 .1 to obtain [- 1"- " 2 a2(Smu - smin)2s2max.utJ IkA.(/,s) 1 dlds $ �+2X ' IJBiD laiD min (5.75) so kA.(/, s) is square-integrable (et: 5.13). If,\ < 1 and Smin = 0, then we need to assume (5.7), and consider only ,\ > 20' - 1 . Then from (5.70), noting the restrictions tof(P) and b,(/, s) applied in Section 5.1 , (5.76) so as L(s) is finite we may write lL(l) L(sf-A. kA.(/,s ) $ a�MI b,(r, s) r2 tIr. max(',r) (5.77) Then applying (5.7) and an argument identical to that used in deriving (5.18), we obtain -- Ima A 2 a 2MJ��-A.)s�..t-2cr)(l + ,\ - 0') [.m L.m IkA.(/, 8) 1 dlds $ (1 + ,\ _ 0') (1 + 2,\ _ 20') ' (5.78) where Lsup = sup(L(s», O $ 8 $ Smax. Hence the kernel (5.70) is square-integrable if Smin = 0 provided '\ > 20' - 1 , which is not too restrictive because 20' - 1 < 0' < 1. Hence iff(P) and b,(l,s) satisfy the conditions outlined in Section 5.1 , there exists a constant Ao < 1 such that for all '\ > Ao, the kernel kA.(/,s) given by (5.70) is both square-integrable and non-factorable. From the above, we may choose Ao = 0 if Smin > 0, or Ao = 20' - 1 if Smin = O. As pointed out in Section (4.1 .5), this means that the operator K('\) must have a largest positive eigenvalue which is simple, and all other eigenvalues of K('\) have real part strictly less than this. 106 Let hax be the largest positive eigenvalue corresponding to the kernel for some A. > Ao, and let �(l) be a multiple of the corresponding eigenfunction chosen such that ; is non-negative, and not identically zero. Then 7maz;(I) = 1000 k.t(l, s);(s)ds, so substituting from (5.70), multiplying by f- and integrating gives 1iuax 1000 r;(l)dl = 1000 r (1000 k.t(l, s);(s)ds) dl (5.79) _ loo af {OO � lOO/(! )b,(r,s) 11+.tdr;(s)dsdl. (5.80) 10 10 • r r Changing the order of integration, we have (5.81) because/(p) = 0 for p > 1 . Hence as J. b,(r,s)dr = 1 , Vs, and ;(I) is non-negative and not identically zero we have (5.82) • If A. = 1 , then 1S:nax = 1 is an eigenvalue of the operator K{A.), and all othereigenvalues of K(A.) have real part strictly less than one. • If A. > 1 , then as � is monotonic decreasing with A. for p E [0,11, we must have 7max < 1 . This means that 'Y = 1 cannot be an eigenvalue of K{A.). • If A. < 1 , then the monotonicity of p1 implies 1mu > 1 , and all other eigenvalues of K(A.) must have real part strictly less than 1mu. There must be some � say such that if A.t < A. < 1 , then all eigenvalues apart from 1iuax (which is greater than one) must have real part less than one. Hence for � < .t < 1 , it is impossible K(A.) to have an eigenvalue of one. If we now let % = max{Ao,�) < 1 , we can say that the only value of .t with real part greater than A.l for which K(A.) has an eigenvalue of one is A. = 1 . 107 5.3.3 The Inverse Laplace Transform It follows immediately from the previous section that the operator R(A) given by (5.74) has a singularity at .1. = 1, and no other singu1arities in the region .1. > At . Now the operator K(A) is analytic in a neighbourhood of .1. = 1 , so in this neighbourhood we may use the Taylor series expansion 00 K(A) = L (.1. - 1 )"K", (5.83) 11=0 where the XII operators are not functions of .1.. If the order of the pole of R( .1.) at .1. = 1 is p � 1 , then we may write a Laurent series expansion for R( .1.), 00 R(A) = L (.1. - l)"R,,, (5.84) II=-p where R_p :/; o. Now from the definition ofR(A), (5.74), R(A)(I - K(A)) = (1 - K(A»)R(A) = I, (5.85) so substituting the series (5.83) and (5.84), then multiplying through by (.1. - 1)P gives 00 00 ( L (.1. - 1 )"R,._p)(1 - L (A. - 1 fK,.) 11=0 11=0 Equating coefficients of (A. - 1 ) 0 then gives 00 00 = (/ - L (.1. - 1 )"K")( L (.1. - 1 )"R,._p) 11::0 11::0 = (A - 1rI. R_p(1 -Ko) = (I -Ko)R_p = 0, and equating coefficients of (.1. - 1 )1 we have - R-,xl + Rl_p(I -Ko) = ' { 0 p > 1 I, p = 1 . (5.86) (5.87) (5.88) Let us assume that p > 1 and look for a contradiction. Then from (5.87) and (5.88) it follows that (5.89) Now in (5.87), we have Ko = K(l) and R_p :/; 0, 80 the right hand equality can hold only if for all yt, R_pyt is an eigenfunction corresponding to the eigenvalue r = 1 of K(l ). Also Kl = ixK(A.)I.1.::b which is a strictly negative operator, 80 R_pX1R_p cannot be the zero operator. Hencep = 1 and R(A.) has a simple pole at .1. = 1. 108 We now proceed to find the inverse Laplace transform of �(,t), which is given by 1 lc+iOO 4 Cl>(t) = 21fi c-ioo �Cl>(,t)d,t 1 lc+ioo = 21fi c-ioo �R(,t)ci>o(A.)d,t (5.90) where c > 1 as shown in Fig. 5.1 . Consider the integral around the closed path shown in Fig. 5.1 , in the limit as e -+ 00, with � < d < 1. It follows from the previous section that this path encloses just the one simple pole at ,t = 1 , 80 the total integral around the path is equal to 21fi times the residue of the pole at ,t = 1. That is, (5.91) As both � and ci-o(,t) are analytic in a neighbourhood of,t = 1, and R(,t) has a simple pole at A. = 1 , we have (5.92) Now as noted above, R-l is the operator which maps any non-zero function directly to an eigenfunction of the operator K(I ). Provided that Cl>o(t, s) is not identically zero, �o(l ) is not identically zero either. Hence we have (5.93) where ;(s) is a normalised eigenfunction of K(I), and C is a non-zero constant dependent on the ci-o(l ) and therefore on 4»0. We note that K(I ) is the integral operator with a kernel given by (5.8) or (4.37), so ;(s) is a multiple of the steady birth-size distribution as discussed earlier. Ifwe choose to normalise ;(s) using (4.33), then the function ; referred to in (5.93) is exactly the same function ; as in Chapter 4. We note that in the limit as e -+ 00, the integrals along the upper and lower edges of the path in Fig. 5.1 approach zero, 80 from (5.90) we have 1 [ A t+iOO 4 1 Cl>(t, s) = -2 . (21fi) (residue(�R(,t)4to(,t,s)) I.l=l ) + . �R(,t)Cl>o(,t,s)d,t . 1f1 a-lOO Now the residue at ,t = 1 can be found from (5.92) and (5.93) to be residue (e.1lR(,t�o(,t)) l.l=l = tc;(s), and 1 2 1 . r�iOO �R(A.�o(,t)dA.1 5: M(s)eD, 1fl 1d-100 109 (5.94) (5.95) (5.96) I I d�e�--------�--------1C� I C R(A) d • 1 d-ie�--------------""",c-ie r Figure 5.1 : The integration path in the complex plane, r, used to estimate «I»(t,s). We choose c > 1 , A.I < d < 1 , and let e -+ 00. 110 where M(s) = 2 1 11 00 R(d - il1)d>O(d - il1,s)d771. K -00 Substituting (5.95) and (5.96) into (5.94), we get fZ>(t,s) - C;(s)t � M(s�, VAt < d < 1 , so for large t it is possible to write fZ>(t,s) = t(C;(s) + 0(1 ) ) (5.97) (5.98) (5.99) as required, because i(l -d)t must decay with t for d < 1 . Hence, given the conditions outlined in Section 5.1 , the steady birth-size distribution can be regarded as stable, 80 the steady size and steady sizelbirth-size distributions will also be stable. We note that (5.99) is expressed in terms of dimensionless time. Ifwe transform back to real time using (5.40), and let the birth rate function in terms of real time be f%>'(t,s) = cz,(t',s), then we have simply f%>'(t,s) = e'l(I)(C;(S) + 0(1 )), (5.100) so the stability in terms of dimensionless time implies stability in real time. 111 5.4 Examples and Counter-examples Section 5.3 shows that the sufficient conditions for unique steady size distributions (SSDs) and steady size/birth-size distributions to exist, discussed in Section 5.1 , are also sufficient to ensure that they will be stable. By this we mean that given an initial cell size distribution no(l, s) such that no is not identically zero, no is bounded, and no = 0 outside the feasible region shaded in Fig. 4.4, then after a sufficiently large time we will have n(l, s, t) = I'(I)(Cy(I,s) + 0(1 » , (5.101) where y(l,s) is the steady sizelbirth-size distribution given by (4.24). 5.4.1 Some Stable Cases We consider here the cases solved analytically in Chapter 4, assnming that in each case the initial size distribution no(/,s) is zero outside the region defined by Smin < S < Smu, S < 1 < L(s). The simple case discussed in Section 4.2.1 , where all cells divide at the same size, with bl(l, s) given by (4.47), is stable provided that!(p) is bounded. As L*(s) = L(s) = L, we have Smin = fJ1L from (5.2), and Smu = fJ2L from (5.3), which satisfies (AI) and (...4...). 'Ib check that (A.,) is satisfied, let us consider first (5.5). We have for s � Smin + e, e PlL*(s) = fJ1L = Smin :$ s(1 - s), (5.102) so (5.5) is satisfied with 6 = e/smax. We can show that (5.6) holds in a similar manner, 80 (A.,) is satisfied. In the particular case where SmiD = 0, we have for s :$ I < L, 1L "'(r, s) dr - .!. I r2 - L2 ' (5.103) so we can satisfy (Ao) by choosing a = 0 and Mb = 1 /L2. If we also choose a bounded!, 80 that (Af) is satisfied, then it follows that the case where all cells divide at the same size is stable. In the case where the hazard rate is independent of initial cell size, discussed in Section 4.2.2, we have L* and L as constants with L* < L, and also Smin = PlL* and Smax = p,p: (Fig. 4.7). Conditions (At) and (Am) are clearly satisfied, and (Ar) follows in much the same 112 way as in the case where all cells divide at the same size. AB Pl > 0 so Smin > 0, we do not need to check (..40), so we conclude that the steady sizelbirth-size distribution will be stable provided that/(p) is bounded (so (.AI) is satisfied). In Section 4.2.3, we discussed the case where b,(l,s) = 6(1 - L(s», with L*(s) = L(s) = Lt-rsT, and/(p) = H(l - p) SO Pl = 0 and p2 = 1 . Conditions (A,) and (At) are clearly satisfied, and from Section 4.2.3 we have "miD = 0 and Smu given by (4.85), so (A".) holds. 'lb show that (�) is satisfied, we note first that PlL*(s) = 0 so (5.5) is satisfied, then we consider (5.6). Given E > 0, we have for s 5 Smu - E, (5.104) As in Section 4.2.3, we assume that a � 2, so from (4.85) we have Smax ::; L. Hence S < L - E and the factor (�)I -r is strictly greater than one. Thus there exists a 6 such that p2L(s) � s(l + 6), so (A,) is satisfied. Finally, as b,(l, s) = 6(l - L(s» , with L*(s) = L(s) = Ll -rST, we have for all s ::; I ::; L(s), lL(.f) b,(r,s) 1 1 I r2 dr = L(s)2 = L(s)O'L(s)2-O' · (5.105) Now if we choose a = r'11 < 1 , and note that L(s)O' � zcr and L(s)2-O' = (Lt-rST)2-O', we obtain after some simplification lL(.f) b,(r,s) 1 I r2 dr ::; L2(I-O')(ls)O" (5.106) so (..40) is satisfied with Mb = L -2(1 -0'). Hence the case dealt with in Section 4.2.3 is stable. 5.4.2 Some Unstable Cases We consider here some examples in which one or more of the conditions specified in Section 5.1 are violated in some way, and show using elementary arguments in each case that a unique stable sizeA>irth-size distribution will not in general be approached at large time. Multiple solutions for Smin and Smax Let us choose L*(s), L(s), Smin, and Saura: such that conditions (A,) and (A".) are satisfied, and choose/(P) to be bounded (so (At) holds) with Pl > 0 (so (.Ao) is unnecessary). However, 113 we violate (�) by allowing there to be points Smin,; E (Smin, smax) such that PlL*(Smin,;) = Smin,;, and points Smax,; E (Smin, smax) such that s 8min,2 8 . 1 -. 8mos,2 (5.107) (5.108) 1 Figure 5.2: A case where the steady size distribution is not unique because (�) fails. We have smin/L*(smin) = Smin� /L*(smin; d = smin:J,/L*(smin:J, = Pl and smax/L(smax) = smax� /L(smax� ) = smax:J,/L(smax:J, = P2. All cells divide between the lines 1 = L*(s) and 1 = L(s). The example we consider is shown in Fig. 5.2, where there are two points, Smin� , and Smin,2, satisfying (5.107), and two points, smax� and smax:J" satisfying (5.108), with Smax� < Smax,2 < Smin� < Smin,2 as shown. If we restricted the domain of possible cell birth-sizes to (Smin, smax� ), by choosing no(/,s) to be zero outside the lower shaded region in Fig. 5.2, then a steady birth-size distribution would develop with support (Smin , smax,l ), say � (s). Similarly, If we restricted the domain of possible cell birth-sizes to (Smin,2 ,Smax), by choosing no(/,s) to be zero outside the upper shaded region in Fig. 5.2, then a steady birth-size distribution would develop With support (Smin,2,Smax), say �(s). Hence in this case a steady cell birth-size distribution may develop given any no(/,s), but on (Smin, smax) this distribution will not be unique, as changing no{l,s) will change the linear combination of � (s) and �(s) which will make up the final steady birth-size distribution. 114 Mean cell size always decreasina We first consider the way in which the mean cell size changes with the generation-by­ generation iteration ft+1 (/) = a looo K(/, s)�(s)ds. The mean cell size after the iteration is - foOO 1;1,:+ldl a foOO 1 1000 K(/, s);';(s)dsdl �+1 = foOO ;';+ldl = foOO fO K(/, s);I,:(s)dsdl ' (5.109) (5.110) which, after manipulations simi1ar to those used in deriving (4.43) from (4.40), gives - fOO s;';(s) log(3), Jl(a,x) = O. (5.120) (5.121) (5.122) Here a stands for cell age and x for cell size, and b(a,x)dt gives the probability that a cell of age a and size x divides in a short time interval dt (in the notation of Section 1.8, b(a, x) = h/(x, a)g(x, a» . Then the minjmum age at which a cell can divide is ao = log(3), and all ofHeijmans' assumptions (A,), (Ab), (Ap), and (Ad) are satisfied, as well his Assumption 6.4. However, it is a simple matter to show that a cell which is born at size x will reach size 3x + 2 at age ao, so its division size must be greater than size 3x + 2, and both daughter cells will have size greater than 1 .5x + 1 , which is strictly greater than x. Hence all daughter cells will be born at a strictly larger size than the parent cells, 80 no stable size distribution is possible. f(P) a delta function Here we consider the the stability of the case dealt with in Chapter 3, where cells divide into a daughters all of the same size, so f(P) = 6(p - a-I ) , (5.123) where 6 is the Dirac delta function. A number of authors have pointed that a population of cells growing exponentially in size and dividing exactly in half will not, in general, approach a stable size distribution (e.g. Trucco [59], Hannsgen et al [21], Tyson and Diekmann [61 ], and Heijmans [24]). However, as we found in Chapter 3, a steady but unstable size distribution can exist, and we note the justification given a footnote in 1Yson and Diekmann [61] for studying the steady size distribution in this unstable case: 117 "Our rationale is that any small aberration from either true exponential growth or exactly symmetric division yields stable distributions which are close to the ones we calculate." Clearly, iff(P) is given by (5.123) then assumption (AI) cannot be satisfied, 80 we would not necessarily expect this case to be stable even if all the other assumptions were satisfied. A simple "mind experiment" shows that a stable size distribution cannot be approached because there is no "spreading" of cells into wider size-classes. Consider a group of cells whose sizes at time t = 0 all lie between 11 and 12 (Fig. 5.4(a», but which may have any distribution of ages or birth-sizes, and for simplicity choose r(t) = 1 80 the growth rate is simply g,(l) = I. Any offspring produced by these cells in the next instant of time will lie between sizes a-lit and ·a-112, as shown. A short time 'r later, we have the situation shown in Fig. 5.4(b), where the original cells now lie between sizes 11 e� and 12e� and their daughter cells have all grown to sizes between a-lit e� and a-lI2e�. But any cells newborn at t = 'r will fall exactly into the latter size range, i.e. (a-lit e� , a-1 12e�), and all cells born between time t = 0 and t = 'r will also lie in this size range at time t = 'r. Fig. 5.4(c) shows the situation after a large time, with cells from a number of generations present. The cells from each generation are confined to a small range of sizes, with the size limits increasing exponentially with time until all cells of that generation have divided. Hence no stable size distribution is approached. 118 (a) t = o size ,I (b) t = T size ,I (c) t -+ 00 Figure 5.4: Time series of size distributions where cells grow exponentially in size and divide into exactly a equal-sized daughters. (a) At t = 0, an initial size distribution has support (h , 12), with daughter cells being born in the size interval (a-Ih , a-112). (b) A short time -r later, all cells born between t = 0 and t = -r lie in the size range (a-I ll e', a-l/2e'). (c) At large time, as t - 00, each generation of cells still present lies within a range of sizes (lkl , Ik2), where lk2/1ld = 12/h . JI 9 Appendix A MODlents of a LiDliting Distribution In Chapter 2, we showed that in the limit as a -+ 1 +, the solution of the functional differential equation y'(x) = -ay(x) + aay(ax), x > 0, subject to the normalising condition fooo y(x)dx = 1 , tends to a normal distribution. H this is the case, it follows that all the moments of the y(x) given by (2.42) should tend to the corresponding moments of the normal distribution as a -+ 1 + . A.t Formulation as a Combinatoric Limit Rewriting the expression given by (2.24) for the nth moment about the origin in terms of q = a-I , we obtain It! E{x"] = a"(1 - q)(1 - if) . . . (1 - if) ' (At) for all positive integers It, where the notation EO indicates the expected value. We write the mean and variance (from (2.25) and (2.26» as 1 Jl = a(1 _ q) (A2) and 1 a" = 02(1 _ q2) (A3) 120 respectively. Then, using the binomial expansion of (x - JJ)", we obtain so that 11 n I ll_A: k!( -1 ) . ( ) " -le E[(x - JJ)"] = k k ( a(l - q» al(l - q)(l - q2 ) . . . (l - qi) ' ri-le X - JJ _ , _ .!l i � (-1 ) E[( U )"] - n.(l q ) � (n _ k)!(l _ q)"-A:(l _ q)(l _ q2) . . . (1 _ qi) ' Borrowing some notation from combinatorics, we may write this as x - JJ 11 (1 + q) I 11 .n-k n! E[(--q) ] = 1 - q �(-1) (n - k)!k!q ' where we define klq to be A: 1 - tI 1-1 . k!q = TI -1- = TI (l + q + " + . . . + q') , i=1 - q i=O (AA) (AS) (AS) (A.7) as in Goulden and Jackson [14]. Note that in the limit as q - 1 - , klq approaches the usual factorial, so k!q is called a q-analog factorial. lim k!q = kt, q ..... l- (AB) Now if the SSD, y(x), tends to a normal distribution as q - 1 - (i.e. as a - I +), then the distribution of the variable (x - JJ)/u must tend to the unit normal distribution N(O, 1 ). The nth moment of the normal distribution in given by JI:(O�) _ { 0 , n odd - (n - I )! ! , n even ' where (n - 1 )!! is the so-called "double factorial", which for n even is given by n' (n - 1 )!! = 1 .3.5 . . . (n - 1 ) = 21 � " I' (A9) (A10) Comparing (A. b) with (A9), for the limit of the distribution of (x - JJ)/ u as q - 1- (i.e. as a - I +) to be a unit normal distribution, we must have lim [(1 + q) I ± ( 1 t n! 1 { 0 , n odd q ..... l- 1 - q A:=O - (n - k)!klq = (n - I )!! , n even 1 21 (A.ll While there are some results in the literature which are superficially similar to (All) (Goulden and Jackson [14]), this result itself does not appear, nor apparently any other results from which (All) could easily be deduced. The following proof starts from what could be regarded as "first principles" and is therefore rather long! A.2 Formulation Involving Polynomial Coefticients From (A 7), we have for le, n E 1+ and le < n, n!q = k!q(1 + q + . . . + t/) . . . (1 + q + . . . + tt-I ), so the left hand side of (All ) can be written as lim q ( I t n. [ (1 + ) 1 11 ' ] q-+l- 1 - q lo - (n - k)!k!q (AI2) = lim [ (1 + q)l ] Il:o( -118(1 + q + . . . + tf) . . . (1 + q + . . . + (,-1 ) . (A13) 9_1- n!q (1 - q)1 The limit as q - 1 of the fraction in square brackets in this expression is �. Expressing the right-hand fraction in terms of p = 1 - q gives where lim [ (1 + q) t I, (-I )1 n! ] _ 21 lim !Jil q-l- 1 - q 1::0 (n - k)!klq - n! ,-.0+ pi ' A n' 11-1 I f(P) = k�( -1)" (n -·k) !!] �(1 -pi· For (All ) to hold, it follows from (A14) that . [f(P)] { 0 ,l!..r:.. pi = 1I!(�i)!! , n odd , n even. (Al.) (AtS) (A16) Now f(P) is a polynomial in p, so the limit in (A16) can exist only if all coefficients of p' in f(P) are zero for r < i. The coefficient ofp' with r = i (n even only) must be chosen 80 that (All) holds, but coefficients ofP' for r > i are irrelevant. Hence (All) will hold provided that 122 { 0 , r < i [p'1f(P) = 1I!(�"t)!! , r = i, n even ' where the notation flf1f(P) is used to indicate the coefficient ofp' inf(P). A.S An expression for fJfV(P) Now (l -pi = I,(-l Y( i )JI j=O 1 so after a change of order of snmmation we have ±(l -pi = ±(-1YJI± ( i ) . .=0 }=o .=} 1 We can now apply the combinatorial identity ± ( � ) = ( � + 1 ) , i=j J J + 1 which can be proved easily by induction on I, to obtain ±(1 -pi = ±(-1 Y ( 1 + 1 )JI. i=O j=O 1 + 1 Substituting this result into (A.15) we obtain For each k, consider 11 n! 11-1 I ( I + 1 ) i 1(P) = L(-1t (n _ k), fI � i + 1 (-p) . 1=0 • l=k .=o (A.l7) (A.l8) (A.19) (A.20) (A.21 ) (A22) (A23) as a product of (n - k) terms, which we number according to the I value, from k to (n - l ). The constant term in the product is then simply the product of the constant terms in each of the factors, (A24) 123 so that [p0lf(P) = f, ( -1 )t n! n! = n! t ( -ll ( n ) . 1=0 (n - k)! k! t::O k We can then apply the identity so for all n > 0 we have f,(-lt( n ) = { 0 1::0 k 1 , n > 0 , n = 0 (A.25) (A.26) Hence in order to establish (A.17) we may now restrict our attention to 0 < r � j. Note, in particular, that the term in (A.22) with k = n is a constant, so we ignore that term from here on and sum over all k up to (n - 1 ) only. Note also that we may assume n � 2r. (A. 27) Now for any r > 0, we choose some partition ! = (it , i2 , . . . , i",) ofr, such that r = it +i2+ . . . i"" and ij > 0, Vi. Note that therefore the order of the partition, m, must satisfy m � r, (A.28) and equality can occur only in the case where ! = (1 , 1 , . . . , 1 ). For the chosen partition b we choose the constant term from all factors in (A.23) except that we choose the term involving JIt from the it th factor, the term involving Jia from the hth factor, and in general the term involving pit from the if th factor, with it < 12 < . . . < 1",. Note that if m > n - k, this would not be possible, so such a partition could not give rise to terms in p' for the chosen k value. We then calculate the coefficient of ( -p)' corresponding to the selection i = Ut. ,12 , · . . 1",) of factors and the partition b by taking the expression (A.24) for [po] and adjusting it for each of the factors if' so· where 11-1 I ( I + 1 ) . n ' '" [( -P)']ij n L . + 1 (-p)' = k; n Jr(q), - - Id i:::O I • 9=1 124 (A.29) This notation ensures that if there is no term involving p� in the chosen factor, then 1f(q) = O. Since [(-p)'Y(P) = (-1 ),rP'1f(P), summing over all partitions ! (which implicitly sums over all m as well), and summing within each over all possible factor sequences � leads to [JlY(P) = n!(-l)' '':£ (-l )k ( n ) I k=O " ! "f 6(1 ) /I-I-1 ) 6(2)... /If 1f(q)... 1: 6(m - 1) /If 1f(m). (A 30) it =k h=it +1 jf=jf-l +1 j.-l =i_s+1 j.=j. - l +1 If le > n - m, then any partitions of order m should not contribute to the sum. This can be achieved by ensuring at all stages that if the lower limit of a sum is greater than the upper limit, a contribution of zero is made to any resulting total. We now shift the summation with respect to ! to the left, and reverse the order of all the other sums to get A.4 Proof of the Main Result (A31) In order to progress step-wise through the multiple sum in (A 31 ), we introduce the abbreviation Sq to represent the inner sum in this expression back to summation with respect tojq. That is, i.+l -1 i. -1 .is-1 it ( n ) Sq == I 6(q) . . . I 6(2) I 6(1 )n!( -1 'f I ( _1 )k . j.=IfI-q .hi =1 it::O k=O " (A32) Note that if q = m, we replace the upper limit of the outer sum by n - 1 so that S". gives the total coefficient of If for the particular partition !. Given n > 0, it can be proved by induction on m that (A33) 125 80 the sum with respect to it in (A31) is SI = L K(I )n!( -1 y( -1)'1 , il-l . ( n - 1 ) i1=O }1 (AM) We now substitu� Cor K(I ) using (A30), and use the Cact that K(I ) = 0 Cor it < it to raise the lower limit of this sum from it = 0 to it = 11 , A change of summation variable from h to it - il then yields (-I)it .;.-it-l . (n - I )! SI = n!(-I )' C l )' 1:, (-1)'1 ( j , 1)" '1 + . h::/j n - 1 - '1 - , (A35) which can be written as i (n - I )! b-it -1 . ( n - (il + 1 ) ) SI = n!( -1 n -1 ) 1 C 1 )'( (i 1 » ' 1:, (-1)'1 , . , + , n - 1 + ' h::/j 11 (A 36) In the next section we show by induction on q that Sf can be written as i,+1 -6(f) . ( n _ S(q» ) 6(!t-l i,+l-C . ( n - c ) Sq = 'I'(q) L (-I Y' + 2, ac L (-1 Y' , i,:O if t=f i,:O }q (A37) where the Dc can be any constants, and Now s(q) = t(ic: + 1 ), (A 38) c=1 '1'( ) = n'(-1 )'(-1 )�1 � (n - 1 )! 1 (A q , fi1Ll(i.t + l )!](n - s(q» !m.:f(n - s(d» ' 39) rP'"Y(P) = 1:,S", i (A40) so substituting q = m and s(m) = I:'=l (ic + 1 ) = r + m into (A37) we have where the Dc are a set oC constants. 126 In (A41 ), e � r + m - I , 80 e � 2r - 1 from (A28) and, using (A.27), e < n. Hence n - e > 0 so from (A.26) the last sum in (A41 ) must be zero for all e, so all terms involving the constants tic are zero. Also, from (A26), the first sum in (A41) is zero for m < r because n � 2r. However, as pointed out after (A28), the only partition of r for which m I. r is the partition ! = (1 ,1 , ... , 1 ), so this is the only partition we need retain in (A41). Hence, substituting iq = 1 , Vq into (A39) and (A41) we get, for n � 21', , (n - I )! 1 a-2r . ( n - 2r ) fp']f(P) = n!( -1 ) (-1 Y (2')(n _ 2r)! (n _ 2)(n _ 4) ... 4.2 i-'I:o (-1 Y- j", . Since this sum is zero for n > 2r, and one for n = 2r, we get finally [p"lf(P) = { '���1l11 which is equivalent to (A.17), as required. , n > 2r , n = 2r ' A.5 Proof By Induction of the Expression for Sq (A42) (A.43) We know from (A36) that the expression for Sf' (A37), is correct when q = 1 . We now assume (A37) is true for some q and show it must then be true when q is replaced by q + 1 and hence true for all q up to n. Then we have, from the definition (A32), if+I-1 Sq+! = I "(q + 1 ) if+t=q [,¥(q)it+t( f)(_I-;' ( n -.s(q» ) + '(£"1 t1cif+fC(_I '/t ( n � e )] . jt=O }II t=q k=O }q (A«) If we apply (A33) to both the sums with respect to }q, and raise the lower limit on the outer sum in each case to take account of the fact that the inner sum is zero if its lower limit exceeds the upper, we obtain (A45) 127 Now we note that for any c > 0 and i � }, it is possible to choose a set of constants bd such that 1qUq - 1 ) ... Uq - iq + 1 ) = Uq - c)Uq - (c + 1 ) ... Uq - (c + iq) + 1) i,-1 + I, bdUq - c)Uq - (c + 1 )) ... Uq - (c + d) + 1) (A.46) d::O The coefficients, bd, could be evaluated by equating coefficients of powers of}, but we avoid this step as we are not interested in the actual values. Hence, substituting (A.30) into (A45) and rearranging, we can write for some new set of constants bct/, .�) �+1 -1 jf+S-1 + I, 1: bct/ 1: c=q d=O jf+1 =+4 [Uq+! - C)Uq+l - (C.+ 1 )) ... Uq+l - (C + d) + 1 ) (_1 Yt+1 ( n - c - 1 )] . ('q+l + 1 )! iq+l - c (A.47) Note that the lower limits of the outer sums have been raised again because the product is zero in the intermediate cases, and that the unknown constants in the expansion of the first sum have been included in the second, which is why the upper limit on the sum with respect to d has been raised to s( q). We now note that the last sum in (A. 4 7) can be rearranged to give jf!1 [Uq+l - C)Uq+l - (C + 1 )) . . . �q+l - (C + d) + 1 ) (_1 Yt+l ( n - c - 1 )] jf+l=+4 (Iq+l + 1 ). iq+l - C = (n - (c + 1 ))(n - (c + 2)) . .. (n _ (c + d)) jf r 1 (_1 yf+l ( n - (c + d) - 1 ) (A.48) jt+l=+4 1q+l - (c + d) and that a similar rearrangement of the first sum with respect to 1q+l can be made, with c replaced by s(q) and d replaced by iq+l . We now change the summation variable from}q+! to 1q+l - (s(q) + iq+l ) in the first sum, and from 1.+1 to }q+! - (c + d) in the last, and also replace the summation variable c by c + 1 . Then, grouping "'like binomials", we have for a new set of constants llc, [ (I .:(� 1 ) , ( -1 )* )(n - .r(q) + 1 » ... (n - (s(q) + I<+t l) ( -1 )'(f) ( -1 )'] 128 (A49) where the upper limit of the sum with respect to c has been increased from s(q) + 1 to s(q) + i.+l by setting the appropriate constants Ilc to zero. Now s(q) + (if+l + 1 ) = s(q + 1 ) and the first expression in square brackets of (A49) simplifies to 'I'(q + 1 ), so we have finally as required. 'I'(q + 1 ) i f+I�f +l ) [(-1 yf+l ( 11 - s(q + 1 ) ) ] if+l::O if+l + L tic L (_1yf+l . - .r(HI )-1 i.-+.-c [ ( " c )] C=f+1 it+l::O )f+1 129 (A50) Appe ndix B The case b(x) · = b with g(x) = gxk Here we consider again the special case described in Section 1 .9, where the product of the hazard rate, hi, and the time-independent part of the growth rate growth rate, g, is assumed to be constant «1 .79) or (1 .80». As pointed out in Section 1 .9, this means that the probability of cell division in a short time interval ell is independent of cell size I and cell age a, but is proportional to the time-dependent part of the growth rate r(t). We assume also that all cells divide into exactly a daughters all of the same size, but instead of assuming that the time-independent part of the growth rate g is a constant as in Chapter 2, we assume g(/) = gf-k, t > 0, (B.1) where the g on the right hand side is a constant. Substituting (B.1) into (1.82), the functional differential equation for the SSD y(/) becomes Now if we substitute and into (B.2) we obtain d dlgr -ly(/) = -bay(/) + ba2y(al). Z(/) = g(I)y(/) = g/l-ly(l) ba a = -g Z'(l) = r-1 ( -aZ(/) + a�Z(al). 130 (B.2) (B.3) (BA) (B.5) As pointed out in Section 3.1 , this is exactly of the form (3.20), the solution of which is given by (3.25), _ 00 (-1 )"ab' -tJ.;t Z(l) - C! (al - 1) ... (abt _ l )e (B.6) 1b find the constant C in this case, we use the normalisation condition obtained by substituting (B.3) into (1 .43), which is [00 \JIt-l 10 Z(l� dl = g. Substituting the solution (B.6) into (B.7), we obtain 00 00 (-1 )"Ir" -1 -tJ.;t [ _� ] -1 C = g 10 ! (ak _ 1) . . . (abt _ 1 )r e dl (B.7) (B.B) As k > 0, we may now simply interchange the order of summation and integration, and perform the integration term by term to obtain _ 1: (-1 )"ab' [00 r-1 e-� dl ,,::0 (a! - l ) . . . (aAII - 1 ) 10 = I, (-Ira: 1 ,,::0 (ak - l) . . . (a - l ) aaU - !. I, a! (-1 rAIl ' (B.9) a ,,::0 ( - l ) . . . (a - 1 ) Now if we use the definition of K( a,/J) given by (3.38), the sum on the right hand side is simply K( at, 1 ), so from (B.B), C is given by ga C = K(aA,l ) , Hence the solution of(B.5) subject to (B.7) may be written as ga � (-1)"ab' _� Z(l) = K(ak,l ) ,.� (al - 1 ) . . . (abt _ l )e , so the SSD y(x) is given from (B.3) as _ alk-1 00 (-1)"ab' -.;t y(l) - K(ak,l )! (al - 1 ) . .. (abt - 1/ . 131 (B.IO) (B.ll ) (B.12) Appendix C The General Equation Y' (x) = -p(x)y(x) + q(x)y( ax) Equation (1 .55) gives the general form of the functional differential equation arising in the case of a population structured on size only, where all cells divide into exactly a pieces all of the same size. Here we develop a general series solution for equations of this form, assuming that a value for Q is given. Consider the functional differential equation of the general form y(x) = -p(x)y(x) + q(x)y(ax), (C.1 ) where a > 1 and we are interested only in the region x > 0. We wish to find an expression for the solution to this equation, assuming for the moment that it exists. Firstly, let us assume that there is a maximum value of x, x* say, such that y(x) == O,x > x*. (C.2) We then write y(x) as { 0, x � x* y(x) = YA:(x), x* a-A: $ x < x* a-(A:-l ) . (C.3) We will refer to the range of x values x*a-l $ x < x*a-(l-l) as region k, so the regions are labelled in the same way as in Fig. 3.3. The solution in region k we will call Yl(X). 132 In region 0, we have simply yo(x) = o. In region 1 , from (CA) we have y(ax) = 0, 80 (C.1 ) gives >1 (x) = -p(x)Yl (x), so ifwe let P(x) be an anti-derivative ofp(x) we have Yl (x) = Ce-P(z), where C is an arbitrary constant. In region k + 1 , k > 0, (C.I ) gives Yt+l (X) = -p(X)yl+l (X) + q(X)Y1(ax), which if we let for all k simplifies to W1+1 (x) = -af(ax)Wk(ax), where (CA) (C.5) (C.6) (C.7) (C.B) (C.9) (C.IO) We now assume that y(x), and therefore W(x), is continuous at x = x*a-1 so (C.9) can be written in integral form as �a1� W1+l (X) = W1(x*a-1) + lax /(S)W1(s)ds. (C.lI) From (C.6) we know Wl(X) = C, so it is possible to apply (C.lI) repeatedly. By induction, we can then show so the solution to (C.I) can be written as y(x) = Ce-P(z) 1 + L 1 /(s,,) 1 /(s,,-t } . . . I /(st}dst . . . ds,.-lds,. , [ 1-1 �a;l-. x-al-· � 1 11=1 ax a:.r. J cut (C.12) x*a-1 � x < x*a-(k-l) (C.13) 133 in the case where (C.2) holds. Note that for all x > 0, only a finite number of integrals are required, but if x· is particularly large or x is particularly small, the number of integrals required could be very large. In applications where the x is cell size, there will normally be a minimum size X() say below which it is known that y(x) = 0, so for all "feasible" cell sizes x in the range X() $ x $ x*, the number of terms in the sum in (C.l3) will not exceed log(x* /Xo)/ log(a). Now if we choose some x and let x* -+ 00, then it is apparent that k -+ 00 also. Hence we can find the solution to (C.l) without the restriction (C.2) by taking limits of (C.l3) as x· -+ 00 and k -+ 00. The result, independent of the order in which these limits are taken and the value of x, is -P(:t) [ 00 100 100 100 1 y(x) = Ce 1 + l: I(s,,) I(s,,-I ) . . . l(st)t/.rJ. . . . ds"_lds,, , 11=1 ca cu. cut (C.14) provided of course that the integrals exist and the sum converges. A necessary condition for all the integrals to exist is clearly that lim:t_oo!(x) = 0, but it appears that the most straightforward way to check whether the series solution (C.14) exists is to calculate the multiple integrals and substitute to check for convergence of the series. All solutions given in this thesis for the case where cells divide into exactly a pieces of equal size can be regarded as special cases of either the infinite series (C.l4) or the finite series (C.l3). For example, in the case considered in Section 3.3, in (3.20) we have p(x) = axk-1 (C.l5) so P(x) = ixl, (C.l6) and q(x) = �xl-l . (C.l7) From (C.lO) we then have (C.l8) so (C.l9) 134 Once again using induction, we can show that 100 100 100 (-1 yaaMe-(czM-l)fzi z I(s,,) I . I(s,,-d · · · 12 I(st }dst . . . ds,._lds. = (at _ 1 )( a21 _ 1 ) . . . (aM _ 1 r (C.20) The solution then becomes _ -Ix' [ 00 (-1 )"aMe-(aM-l)fzi 1 y(x) - Ce 1 + .� (at - 1 )(a21 - 1 ) . . . (aM - 1 ) ' This is just a slight rearrangement of the solution for Z(x) given in (3.25). 135 (C.21) Appe ndix D Generalised Solutions of the Functional Differential Equation We consider here the possible generalised solutions of (3.8), which can arise at size x = 0 because the growth rate is zero at size zero. Firstly, we show that the Dirac delta function 6(x) is a solution of(3.8) for any given b(x). Intuitively, this corresponds to the idea that if the growth rate is zero at size zero (as it is for g(x) = gx), then given'an initial population all of size 0, the whole population will remain at size 0 for all time! The following elementary results hold for any function b(x): and we note also that 1000 b(x)6(x)dx = b(O), b(x)6(x) = b(0)6(x), b(ax)6(ax) = b(O) 6�) , xcS'(x) = -6(x). Substitution ofy(x) = 6(x) into (3.8) then leads to the trivial identity -g6(x) = -[b(O) + (a - 1 )b(O) + g]6(x) + a2b(O) 6(x) , a 136 (D.l) (D.2) (D.3) (D.4) (D.S) so y(x) = 6(x) is always a solution of (3.8). Note that 6(x) also satisfies the unit integral condition (3.4) and the non-negativity condition (3.5), so y(x) = 6(x) satisfies all the conditions required to be a SSD. In general, 6(x) is NOT the only generalised solution with point support the origin. 1b see whether other generalised solutions are possible, we substitute the most general possible solution with point support at the origin, 00 y(x) = L C" cS<")(X), (D.6) 11=0 into (3.8), then equate coefficients of each derivative of the delta function. Using the identities above, we obtain a set of equations of the form (gn - S)CII = ( -1 )"(1 - aI-") � ( j ) bU-II) (O)(-I '/Cj J=ll n where Substituting n = 0 leads to 00 S = (a - 1) L (-1 '/bV1(0) j=O -SCo = -S so either S = 0 or Co = 1 , then substituting n = 1 gives so either S = g or Cl = O. (D.7) (D.8) (D.9) (D.10) For n > 1 the equation contains as many Cj terms as there are non-zero derivatives of b(x) at x = 0, so the situation becomes very complex. However, it is clear that if b{IR) (0) <> 0 but bV1 (O) = 0 for all j > m, then it is possible to choose C2 up to CIR+l arbitrarily, then all C/s for j > m + 1 can be expressed in terms of these. The term ( � ) on the right hand side seems to ensure that the coefficients will get smaller as n increases. In other words, it is possible to in general choose solutions of the form 00 y(x) = 6(x) + L CII6(1I)(X), (D.ll ) 11=2 where the C.'s are non-zero. 137 For example, if b'(O) <> 0 but bV) = 0 for 1 > 1 , then S = (a - 1 )b(0) and we get as a solution 00 y(x) = 6(x) + A [6(2) (x) + L All 6(11) (x)] 11=3 where A is arbitrary and A _ 2 «a - a-I )b(0) - 2g) . . . « a - a2-II)b(O) - (11 - l )g) 11 - [(1 _ a-I ) . . . (1 _ a2-)]II!(b'(O»"-2 (D.l2) (D.tS) Note that because all derivatives of the delta function have zero integral, solutions of the form (D.ll) satisfy the unit integral condition (S.4). We now define non-negativity of a generalised function t(x) in the natural way, such that • A generalised function t(x) is non-negative iffV; E S (where S is the set of "smooth­ test functions) such that ;(x) � 0, "Ix, J�oo t(x);(x)dx � O. Then we have (D.l4) so from (D.l2), i: y(x);(x)dx = ;(0) + ! CII( -1 )";(")(0). (D.t5) Now for any 11 > 1 , it is possible to choose a ; such that -CII( -1 )";(")(0) > ;(0) for any finite CII' and ;v) = 0, "11 :; 11, j � 2. Therefore if CII '1 0 for any 11 � 2, y(x) cannot be said to be non-negative. Hence y(x) = 6(x) is the only non-negative generalised solution of(S.8). 1 38 Bibliography [1 ] [2] [3] [4] [5] [6] [7] [8] [9] [1 0] Abramowitz,M. and Stegun,I.A, CJIandbook of mathematical functions", Dover, N.Y. (1965) And.rews,G.E., "Number theory", W.B.Saunders, Phil. (1971 ) Anderson,E.C., Bell,G.I., PetersonJ).F., and Thbey,R.A, "Cell growth and division IV. Determination of volume growth rate and division probability", Biophys.J. 9, 246-263. (1969) Anderson,E.C., and Peterson,D.F., "Cell growth and division 11. Exper­ imental studies of cell volume distributions in mammalian suspension cultures", Biophys.J. 7, 353-364. (1967) Bell,G.l, "Cell Growth and Division ill. Conditions for balanced exponen­ tial growth in a mathematical model", Biophys.J. 8, 431-444. (1968) Bell,G.I., Anderson,E.C., "Cell growth and division. I. A mathematical model with applications to cell volume distributions in mammalian SUB­ pension cultures", Biophys. J. 7,329-351. (1967) Collins,J.F. and Richmond,M.H., -&ate of growth of Bacillus cereus be­ tween divisions", J. gen. MicrobioL 28, 1 5-33. (1962) Cox,D.R., and Miller,H.D., "The theory ofparticulate processes", Methuen, wndon. (1965) Erickson,R.O., and Sax,K.B., "Rates of cell division and cell elongation in the growth of the primary root of Zea mays", Proc. Am. Phil. Soc. 100, 499-514. (1956) Errington,F.P., Powell,E.O., and Thompson,N., "Growth characteristics of some gram-negative bacteria" J. gen Microbio1 89, 109-123. (1965) 139 [11 ] [12] [13] Feller,W., "On the integral equation of renewal theorY', Ann.Math.Stat. 12, 243-267. (1941) Gompertz,B.J., -on the nature of the function expressive of the law of human mortality and on a new mode of determining the value of life contingencies", Phil. 7hmB. Roy. Soc. 115, 513-585. (1825) Grover,N.B., Woldringh,C.L., and Koppes,L.J., "Elongation and surface extension of individual cells of Escherichia coli B/r: Comparison of the­ oretical and experimental size distributions", J.Theor.Biol. 129, 337-348. (1987) [14] Goulden,I.P., Jackson,D.M., "Combinatorial enumeration", Wiley, New York. (1983) [15] Gurtin,M.E., and MacCamy,R.C., "Nonlinear age-dependent population dynamics", Arch. Rat. Mech. Anal. 5, 281-300. (1974) [16] Gyllenberg,M., "The size and scar distributions of the yeast Saccharomyces cerevisae", J.Math.Biol. 24, 81-101 . (1986) [17] Gyllenberg,M., and Webb,G.F., "Age-size structure in populations with quiescence", Math.Biosci. 86, 67-95. (1987) [18] Hall,AJ. and Wake,G.C., "A functional differential equation arising in the modelling of cell growth", J. Austral. Math. Soc. Ser. B. 30, 424-435. (1989) [19] Hall,A.J. and Wake,G.C., "Functional differential equations determining steady size distributions for populations of cells growing exponentially", J. Austral. Math. Soc. Ser. B. 31, 434-453. (1990) [20] Hall,AJ., Wake,G.C., and Gandar,P.W., ·Steady size distributions for cells in one-dimensional plant tissues", J. Math. Bio. (in press). (1991) [21] Hannsgen,K.B., Tyson,J.J., and Watson,L. T., "Steady-state size distribu­ tions in probabilistic models of the cell division cycle", SIAM J. Appl. Math. 45, 523-540. [22] Harvey,R.J., Marr,A.G., and Painter,P.R., "Kinetics ofgrowtb of individual cells of Escheria coli and Azobacter agilis", J.Bacteriol. 93(2), 605-617. (1967) [23] Heijmans,H.J.A.M., "On the stable size distribution of populations repro­ ducing by fission into two unequal parts", Math. Biosc. 72, 19-50. (1984) 140 [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] . [36] [37] Heijmans,H.J.AM., "The dynamical behaviour of the age-size-distribution ofa ce1l population", in: Metz,J.A..l., and Diekmann,O. (ed), "The dynamics of physiologically structured populations", pp. 185-202. Springer-Verlag, Berlin. (1986) Hochstadt,H., "Integral equations·, John Wuey and Sons (1973) Jolley,L.B.W., "Summation of series", Dover, New York. (1961 ) Kato,T. and McLeod,J.B. "The functional differential equation y'(x) = ay(Ax) + by(x)", BuU. Amer. Math. Soc. 77, 891-937. (1971) Keyfitz,N., "Applied mathematical demography", Wuey, New York. (1977) Kirkpatrick,M. , "Demographic models based on size, not age, for organ­ isms with indeterminate growth·, Ecology 65, 1874-1884. (1984) Koch,A.L. and Schaechter,M. "A model for statistics of the cell division process", J. gen. Microbiol. 29, 435-454. (1962) Koppes,L.J., Woldringh,C.L., and Grover,N.B., "Predicted steady-state cell size distributions for various growth models", J.theor.Biol. 129, 325-335. (1987) Laird,A.K, "Dynamics of tumour growth. Comparison of growth rate and extrapolation of growth curve to one cell", Brit. J. Cancer 29, 278-291. (1965) La80ta,A., and Mackey.M.C., "Globally asymptotic properties ofproliferat­ ing cell populations", J.Math.BioL 19, 43-62. (1984) wtka,AJ., -&lation between birth rates and death rates", Science N.S. 26, 21-22. (1907) Malthus,T.R., "An essay on the principle of population", St. Paul's, London (1798); reprinted in: Malthus,T.R., "An essay on the principle of population and A summary view of the principle of population", Harmondsworth, England Penguin. (1970) M'Kendrick,A.G., 1I Applications olmathematics to medical problems", hoc. Edinburgh Math. Soc .u, 98-130. (1926) Metz,J.AJ., and Diekmann,O., -rile dynamics of physiologically struc­ tured populations", Springer-VerJag, Berlin. (1986) 141 [38] [39] [40] [41 ] [42] [43] [44] [45] [46] [47] [48] [49] [50] Oster,G., "Lectures in population dynamics", in: DiPrima,R.C. (ed), "Mod­ em modelling of continuum phenomena", pp. 149-190. Am. Math. Soc., Providence, Rhode Is. (1 977) Oster,G., -rile dynamics of nonJinear models with age structure", in: Levin,S.A (ed), "Studies in mathematical biology', Vo1 16, Mathematical Association of America, pp. 411-438. (1978) Oster,G., and Takabasbi,Y., "Models for age-specific interactions in a periodic environment", EcoI.Monogr. 44, 483-501 . (1974) Painter,P.R., and Marr,AG., "Mathematics of microbial populatioU', Ann.R.Microbiol. 22, 519-548. (1968) E.O.Powell, CIA note on Koch & Schaechter's hypothesis about growth and fission of bacteria", J. gen. MicrobioL 37, 231-249. (1964) Randolph,AD., and I.arsen,M.A, -rheory of particulate processes", Aca­ demic Press, New York. (1971 ) Ricciardi,L.M., "Diffusion processes and related topics in biology', Lecture notes in biomathematics 14, Springer-Verlag, Berlin. (1 977) Rotenberg,M. , "rransport theory for growing cell populations", J. theor. Bioi. 103, 181-199. (1983) Rubinow,S.I., "Age-structured equations in the theory ofcell populatioU', in: Levin,s.A (ed), "Studies in mathematical biology', Vo1 16, Mathemati­ cal Association of America, pp. 389-410. (1978) Rubinow,S.I., "A maturity-time representation for cell populations", B;' phYB. J. 8, 1055-1061 . (1968) Schaechter, M., WiUiam80n,J.P., Hood,J.R.Jr., and Koch,AL., "Growth, cell and nuclear divisions in some bacteria", J. gen Microbiol 29, 421-4M. (1962) Scherbaum,O., and Rasch,G., "Ce1l si.ze distribution and single cell growth in Tetrahymenapyriformitl GL",ActaPathol. Microbiol. Scand. 41, 161-182. (1957) Sharpe,F.R., and Lotka,A.J., CIA problem in age distributions", PhiL Mag. 21, 435-438. (1911) 142 [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] Simpson-Herren,L., and Lloyd,H.H., "Kinetic parameters and growth curves for experimental tumour systems", Cancer Chemother. Rep 54, 143-174. (1970) Sinko,J.W., and Streifer,W., "A new model for age-size structure of a population", Ecology 48, 910-918. (1967) Sinko,J.W. and Streifer,W., "A model for populations reproducing by fis­ sion", Ecology 52, 330-335. (1971) Smith,D., and Keyfitz,N., "Mathematical demography (selected papers)", Biomathematics Vol. 6, Springer, Berlin. (1977) Streifer,W., "Realistic models in population ecology", Adv.Ecol.Res. 8, 199- 266. (1974) Sundareshan,M.K., and Fundakowski,R., "On the equivalence of math­ ematical models for cell proliferation kinetics", Cell 7issue Kinet. 17, 609-618. (1984) Trucco,E., "Mathematical equations for cellular systems. The von Foerster equation. Part 1", Bull.Math.Biophys. 27, 285-303. (1965) Trucco,E., "Mathematical equations for cellular systems. The von Foerster equation. Part 11", Bull.Math.Biophys. 27, 449-471 . (1965) Trucco,E., "On the average cellular volume in sychronized cell popula­ tions", Bull. Math. Biophy. S2, 459-473. (1970) Tyson,J.J. , "Size control of cell division", J.theor.Biol. 128, 381-391 . (1987) Tyson,J.J. and Die1nn8Dn,O., "Sloppy size control of the cell division cycle", J. theor. Biol. 118, 405-426. (1986) Tyson,J.J., and Hannsgen,K.B., "The distributions of cell size and genera­ tion time in a model of the cell cycle incorporating size control and random transitions", J.theor.biol. 113, 29-62. (1985) Tyson,J.J., and Hannsgen,K.B., "Cell growth and division: a determinis­ ticlprobabilistic model of the cell cycle", J.Math.Biol. 23, 231-246. (1986) VanSickle,J., "Analysis ofa distributed-parameter population model based on physiological age", J.theor.Biol. 64, 571-586. (1977) 143 [65] [66] [67] [68] [69] Verhulst,P.F., "Notice sur la loi que la population suit dans son accroise­ ment", Correspondence mathematique et physique publiee par A. Quetelet 10, 112-121 . (1838) von Foerster,H., "Some remarks on changing populations", in: Stohlman,­ F.Jr. (ed), "The kinetics of cellular proliferation", pp. 382-407. Grune & Stratton, New York. (1959) Webb,G.F., "Random transitions, size control, and inheritance in cell population dynamics", Math.Biosci. 85, 71-91 . (1987) Weiss,G.H., "Equations for the age structure of growing populations", Bull.Math.Biophys. 30, 427-435. (1968) Zabreyko,P.P., Koshelev,A.I., Krasnosel'skii,M.A, Mikhlin,8.G., Rakovsh­ chik,L.S., and Stet'senko,v.Ya., "Integral equations - a reference text", Noordhoff, Leyden. (1975) 144