Coagulation, fragmentation and growth processes in a size structured population. Discrete and Continuous Dynamical Systems

An equation describing the dynamical behaviour of phytoplankton cells is considered in which the eﬀects of cell division and aggregration are incorporated by coupling the coagulation-fragmentation equation with the McKendrick-von Foerster renewal model of an age-structured population. Under appropriate conditions on the model parameters, the associated initial-boundary value problem is shown to be well posed in a physically relevant Banach space.

Lamb, W. and Banasiak, Jacek (2009) Coagulation, fragmentation and growth processes in a size structured population. Discrete and Continuous Dynamical Systems, B (11). pp. 1-23. ISSN 1078-0947 http://strathprints.strath.ac.uk/13318/ This is an author produced version of a paper published in Discrete and Continuous Dynamical Systems, B (11). pp. 1-23. ISSN 1078-0947. This version has been peer-reviewed but does not include the final publisher proof corrections, published layout or pagination.
Strathprints is designed to allow users to access the research output of the University of Strathclyde. Copyright © and Moral Rights for the papers on this site are retained by the individual authors and/or other copyright owners. You may not engage in further distribution of the material for any profitmaking activities or any commercial gain. You may freely distribute both the url (http://strathprints.strath.ac.uk) and the content of this paper for research or study, educational, or not-for-profit purposes without prior permission or charge. You may freely distribute the url (http://strathprints.strath.ac.uk) of the Strathprints website.

Introduction
Phytoplankton is a vital ingredient in the majority of freshwater and oceanic food chains as it is the only food available for many species of fish in their larval stage. An important observation is that phytoplankton cells tend to form aggregates; that is, groups of cells living together. Since larvae do not move on their own, they survive only if they are in the vicinity of the aggregates. 1 Hence, being able to determine the spatial distribution of the aggregates is of importance in the study of fish recruitment.
The distribution of aggregates can be studied at different levels. Individual-based models, which can be thought of as providing 'microscopic' models, track the random motion and division of individual particles; see [7,27]. A 'macroscopic' description is provided by advection-diffusionreaction equations describing the concentrations (densities) of individual particles; see [21]. The model which we study in this paper can be considered as lying somewhere in between, on a 'mesoscopic' scale, in that it recognizes the role played by the phytoplankton aggregates which are now individual building blocks labelled by their size. Thus, we describe the phytoplankton using the aggregate density function u(x, t) which gives the number density of aggregates of size x; that is, 'consisting' of 'x' individual building blocks (note that here x is a continuous variable) at time t. Such a model, fitting into the broad coagulation-fragmentation theory, can be obtained as a limiting case of individual-based models [27] or derived from first principles, as in [18,8,1].
Our starting point is the model considered by A.S. Ackleh and B.G. Fitzpatrick in [1]. To introduce this model we begin by noting that, though effects of aggregation on the dynamics of algal communities had been studied prior to [1], incorporating cell division in the aggregates had presented many difficulties. Thus modelling efforts tended to focus on some special mechanisms for cell division; for example, in [18] the author assumes that either the cells in the aggregate are dead (and thus do not divide), or all daughter cells remain in the aggregate after division of the mother cell. In [1], drawing on some earlier work, the authors assume that all daughter cells fall off the aggregates and, joining the single cell population, enter the system as new aggregates, leaving the size of the original aggregate unchanged by cell division. The resulting problem, analyzed in papers [1,4], can be regarded as a combination of the classical coagulation equation with the McKendrick-von Foerster renewal model of an age-structured population : Here 0 ≤ x 0 < x 1 ≤ ∞ are the minimum and, respectively, maximum size of particles. We note that though only the case 0 < x 0 < x 1 < ∞ seems to make a biological sense, it is difficult to put precise bounds for the size of the aggregates and hence it is customary to consider systems 2 with aggregates of arbitrary size. Also, phase transitions observed in fragmentation-coagulation systems, such as shattering and gelation, require the setting with x 0 = 0 and x 1 = ∞. We also note that when x 1 < +∞, then most of Section 2.1 becomes superfluous. The growth part, however, is more interesting in this case, see the discussion of the boundary condition at x 1 .
We note, however, that here we have a size-structured population and the interpretation of the coefficients r and β differs from that in the age-structured models. In (1.1) the function β represents the rate at which daughter cells enter the single cell population and the function r represents the rate of increase of a size x aggregate due to daughter cells remaining in the aggregate. In other words, the model describes divisions of cells which happen simultaneously with the splitting of an aggregate into an aggregate of the same size as the original one, and a single cell, which then enters the system as a new aggregate. The coefficient µ represents the removal term which can include various mechanisms such as death of cells or sinking of aggregates.
Next, the integral terms describe the standard coagulation process with the 'stickiness function' k(x, y) representing the rate at which an aggregate of size x sticks to an aggregate of size y. The precise form of k can be found in [1]; we note in particular that it is a bounded function and this is the only assumption required in our analysis. We note that if 0 < x 0 < x 1 < ∞, the coalescence of aggregrates can occur only if x 1 > 2x 0 . In this case, in (1.1) (and in (1.3) below), the two coagulation terms should strictly be interpreted as where χ U and χ V are the characteristic functions of the intervals U and V given by The main role in the process of coagulation of phytoplankton is played by TEP (Transport Exopolymer Particles). TEP are a by-product of the growth of phytoplankton and their stickiness make the cells remain together upon contact. On the other hand, low levels of concentration of TEP make the aggregates susceptible to being split by external forces, such as currents or turbulence, [8,14,18,23]. Fragmentation terms also appear in a natural way in individual based models of phytoplankton, [27,28]. Thus, here we consider the full fragmentation-coagulation equation 3 coupled with the McKendrick-von Foerster renewal model: with the initial and boundary conditions of (1.1). As in classical fragmentation models, the coefficient a describes the rate of fragmentation; that is, the number of splitting events per unit time.
The function b(x|y) gives the distribution of x-sized daughter aggregates originating from a parent of size y. We assume that the number of daughter particles is bounded 4) and also that  [2], however the author focuses on bounded coefficients on a finite interval and binary fragmentation; the existence is proved locally in time.
The results of [1] were extended in [4] to cover time-dependent coefficients r and µ with x 0 = 0 and x 1 = ∞. However, the global-in-time existence results depend upon the existence of appropriate upper and lower solutions which are constructed so that a large class of admissible initial conditions is covered. We note a number of other generalizations of the model to include nonlinear effects such as light shading in the column of water, [3], and interactions between species, [5,6], which at this moment are beyond the range of applicability of the theory presented here.
The papers [8,27] are concerned with a problem of the form (1.3) (even with diffusion in space) but with the growth coefficient r not requiring any boundary condition at x = x 0 = 0; also the coagulation term has a different form which makes the the global existence result less involved.
The authors work in the physically relevant L 1 space (where the norm gives the total number of 4 aggregates in the ensemble) but, as their main interest is the limit passage from the individual-based model to the mesoscopic one, the mathematical analysis of the latter lacks some rigour.
The work presented in this paper is carried out in the space where the weight is given by the first eigenfunction of the adjoint problem, possibly up to a suitable normalization, see [24,26] and [25,Chapter 3]. This space is very well suited to studying such problems yielding easy existence and long time behaviour results, either through entropy estimates or the theory of partially integral Markov semigroups. The study of (1.3) in this context is a subject of current research. However, since most applied papers on fragmentation and coagulation models are concerned with an analysis of them in either the norm is a measurable characteristic of the system, our results also are presented in the more traditional setting.
This choice, moreover, has an additional advantage in that the fragmentation operator is known to behave well in the space L 1 ([x 0 , x 1 ), xdx) whereas the transport term has nice properties in provides a convenient framework for a comprehensive analysis of the linear part of the problem (1.3) by means of the theory of substochastic semigroups, [9], yielding satisfactory results even for unbounded coefficients.
The bilinear structure of the coagulation terms makes it straightforward to establish Lipschitz continuity and Fréchet differentiability properties if we work again in the space L 1 ([x 0 , x 1 ), (1 + x)dx).
Thus, when (1.3) is expressed as an abstract Cauchy problem, standard results on semigroup theory lead immediately to the existence and uniqueness of a strongly differentiable, local (in time) solution The advantage of the solution being strongly differentiable is that no further justification has to be provided for pushing a time-derivative through the L 1 norm; see the proof of Theorem 1.1. Moreover, the semigroup approach that we adopt enables us to borrow a simple but elegant trick from kinetic theory (see [13,Chapter 8]) to prove the positivity, and, subsequently, global existence, of solutions, extending earlier results of [10,22]. As far as we are aware this technique has not previously been exploited in studies of coagulation-fragmentation processes. It should be noted that other strategies, based on weak compactness arguments, can be used to establish the existence of solutions to coagulation-fragmentation equations. In the case of continuous models, the first results obtained in this manner were presented by Stewart in [29], and later developments can be found, for example, in [20]. Although these compactness-based approaches can deal with certain unbounded coagulation kernels, they do not lead to strongly differentiable solutions. Uniqueness results also have to be established separately; see [30].
General notation. We introduce the following notation for the spaces: for k = 0, 1 Accordingly, we denote by · k and · the natural norms in, respectively, X k and X. Further, we introduce the dual X ∞ to X consisting of measurable functions f for which With this identification, the duality pairing is the integral Conversely, if f ∈ X ∞ , then clearly it is bounded (a.e.) by an affine function.
The main result of the paper is the existence of global strict solutions to (1.3).
Theorem 1.1 Let the coefficients a, µ, b, β, r and k satisfy the following constraints: (i) a is measurable, supp a ⊂ [x 0 , x 1 ) and there exist non-negative constants P a and Q a such that 0 ≤ a(x) < P a x + Q a a.e.; (ii) µ ∈ L ∞ ([x 0 , x 1 )) and µ(x) + a(x) ≥ 0 a.e.; (iii) b satisfies (1.4) and (1.5); (iv) β is non-negative and β( (for suitably small δ) and r( is non-negative and symmetric. Then the abstract Cauchy problem (1.3) has a unique, global, non-negative strict solution u for each where K is the generator of the linear fragmentation semigroup.
The proof of this theorem, given in Section 3, is based on a detailed analysis of the linear part of the problem and, in particular, a precise characterization of the generator K, which is the subject of Section 2. We emphasize that this characterization is crucial in obtaining a priori estimates in

The Linear Part
In this section we provide a comprehensive theory of the linear part of the problem including The first step is to deal with the transport part and we shall find it convenient to consider: with the same boundary and initial conditions, where q(x) = µ(x) + a(x) ≥ 0. We note that though the theory of first-order equations has been well developed in e.g. [11,12,17], the assumptions of op. cit. are not particularly well-suited for coefficients having singularities, which may appear in our model. Since we require a detailed control of the estimates, it proves more convenient to perform directly on (2.1) rather than try to adapt the general theory to cover our particular case.
As indicated in Theorem 1.1, we assume that µ and a are measurable functions with µ essentially bounded and a satisfying We further assume that supp a ⊂ [x 0 , x 1 ) so that no splitting of particles larger than x 1 (should they exist) can influence the dynamics on (x 0 , x 1 ).
It is well-known, [15,Lemma VI.4.2], that the solution of this problem can be obtained via the solutions of the corresponding problem with a homogeneous boundary condition. We shall recover this result in the process as we work in a slightly different setting.
Before we proceed to examine the well-posedness of the initial/boundary value problem (2.2), we first discuss in detail the assumptions and the meaning of the boundary condition so that the proper definition of the domain can be stated. The positive growth rate r is assumed to be absolutely continuous on (x 0 , x 1 ), denoted by r ∈ AC((x 0 , x 1 )); (absolute continuity on an open set means that r is absolutely continuous in the standard sense on each compact subinterval). To ensure global existence of characteristics if x 1 = ∞, we also assume that r ∈ X ∞ ((x 0 , x 1 )). In addition, by imposing the restriction that (where henceforth the symbols α + and β − will indicate integrals in some right (respectively left) neighbourhood of α (respectively β)) we can define and it follows automatically that We point out that assumption (2.4) allows for r(x) becoming zero at x 0 ; for example, we may have Let us consider the boundary condition at x = x 0 . The simplest clearly is the homogeneous condition which, due to the form of the growth term, should be written as Clearly, if r is continuous at x 0 with non-zero limit, this is the same as saying that u(x 0 , t) = 0, but (2.6) can also allow r to be infinite at x 0 . However, (2.6) has to be modified for a process that includes fragmentation events in which large particles split creating an array of smaller particles since some of these smaller particles may have size x 0 . Such x 0 -sized particles are created by fragmentation at the rate x 0 a(y)b(x 0 |y)u(y, t)dy and enter the population x 0 > 0 at the rate In general, we consider the following boundary condition which covers all these cases 8 where 0 ≤ β ∈ X ∞ . One other feature of the problem that has to be considered is the behaviour of r at x 1 if the latter is finite. There are two ways of looking at the system. One is that we are looking only at the 'window' (x 0 , x 1 ) and ignore what happens beyond the maximal size x 1 . In this case particles can grow bigger than x 1 but will then move outside the scope of our model.
This does not require any assumption on r(x 1 ) but we have to impose a cut-off condition on the fragmentation rate for sizes bigger than x 1 so that they do not influence particles in the interval (x 0 , x 1 ). Consequently the system of particles of size x ∈ (x 0 , x 1 ) will not be closed but the evolution of the system will be governed only by what happens in (x 0 , x 1 ).
The second way of modelling the system is to prevent the particles from growing beyond x 1 . If we formally integrate the equation (2.2) with respect to the measure (1 + x)dx we obtain d dt and it is clear that to prevent any escape through x 1 it is enough to assume that r(x 1 )u(x 1 ) = 0.
We note that it is possible to introduce other mechanisms to prevent aggregates growing beyond a certain limit, such as breakage, see [16,31].
Let us denote by T the formal differential expression and by T max the realization of this expression on the maximal domain D(T max ) := {u ∈ X : qu ∈ X, ru ∈ AC((x 0 , x 1 )), (ru) x ∈ X}.
Our main interest is T max restricted to a domain ensuring that appropriate boundary conditions are satisfied. Thus, we define T β as T max restricted to which takes care of both the homogeneous (β = 0) and the 'renewal' boundary condition for either the unbounded domain (x 1 = ∞) or with x 1 < ∞ but without imposing any restriction on the behaviour of the system at x 1 . We note that accordingly T 0 will denote the realization of T max with zero boundary condition (at x = x 0 ).
If we want to prevent growth beyond x 1 then we further restrict T β to T β,1 defined on When it does not lead to any misunderstanding, to simplify notation, we shall denote by (T, D(T )) the operator T max restricted to any domain with the boundary conditions discussed above.
The general 'formal' solution of the resolvent equation is given by where Q(x) = x x 0 q(s)/r(s) ds is well defined due to (2.4).
There are some identities and estimates which appear throughout the paper. We collect them in the following lemma.
Proof. Consider first the case x 1 < +∞. Then, on integrating by parts, we obtain which in turn gives yielding part (a). For (b), since the right-hand side of (2.15) does not depend on x ′′ , and x 1 = ∞, Since the first inequality in (2.13) is obvious, part (b) is proved.
To prove (c), first we note that Hence, integrating by parts with respect to (1 + x)dx, we obtain The solution to the problem λu−T 0 u = f , that is, to (2.10) with homogeneous boundary conditions, is given by (2.11) with C = 0: x 0 e λR(y)+Q(y) f (y)dy .

(2.18)
In fact, standard estimates give the following lemma.

Lemma 2.2
Under the adopted assumptions, if λ > r ∞ , then R(λ, T 0 )f = u λ defines the resolvent of (T 0 , D(T 0 )) and satisfies the estimate The following lemma establishes that the resolvent of T 0,1 is also defined by (2.18). Then u λ ∈ D(T 0,1 ).
Proof. The only thing to check is that x 0 e λR(y)+Q(y) f (y)dy.
Due to the assumption, lim x→x − 1 R(x) = +∞ and, since Q is positive, the first factor tends to 0 whereas the second may diverge to ∞. Let f ∈ C ∞ 0 ((x 0 , x 1 )). Then lim x→x − Consider next an arbitrary sequence (a n ) n∈N ⊂ (x 0 , x 1 ) converging to x 1 and define a sequence of linear functionals on X by x 1 x 0 e λR(y)+Q(y) f (y)dy = 0, and, in general, there will be an out-flux through x 1 creating particles of sizes larger than x 1 .
Next we turn our attention to the problem with β = 0.
The statement for T β,1 follows since, due to (2.20), we have where u is defined by (2.22). 2 is bounded on X. Indeed,
Hence, the remaining part of this section is non-trivial only for x 1 = ∞ with unbounded a (satisfying (2.3)). However, for the sake of completeness, we formulate all results for the general case.
Proof. First we note that if u ∈ D(T ) ⊂ D(A), where D(A) = {u ∈ X : au ∈ X}, then Bu ∈ X.

Indeed, as in (2.26),
Bu ≤ x 1 x 0 |u(y)|a(y)(y + n(y))(1 + y) −1 (1 + y)dy ≤ M au . (2.28) Thus, we can separate the terms T u and Bu in (2.27). Also, since D(T ) ⊂ D(A) and µ is bounded, we can concentrate on the term (ru) x . Thus, taking x 0 < x ′ < x ′′ < x 1 and using the fact that ru ∈ AC((x 0 , x 1 )) we obtain Now, by the definition of the domain, the left hand side converges to shows that ru is integrable on (x 0 , x 1 ) and so the integral on the right hand side converges to x 1 x 0 r(x)u(x)dx. It follows that r(x ′′ )u(x ′′ )(1 + x ′′ ) converges to a limit ℓ as x ′′ → x 1 . Now, if x 1 < ∞ and u ∈ D(T β ) then clearly we obtain the existence of the limit r(x)u(x)| x=x 1 . If x 1 < +∞ and u ∈ D(T β,1 ), then r(x)u(x)| x=x 1 = 0. Finally, if x 1 = ∞ and ℓ = 0, then r(x)u(x) ≥ a ′ (1 + x) −1 for some a ′ > 0 for x large enough, contradicting the integrability of Combining all estimates we obtain (2.27). 2 Now, we have Defining Using these two concepts we can define an operator that can be thought of as the maximal extension of T + B in X: For the extension of the resolvent of T β we recall that Typically, the extension of R(λ, T β ) is obtained first for positive functions f ∈ E for which the expression on the right hand side of (2.32) defines a function in E f ; that is, finite a.e. Now, we note that the right hand side of (2.32) is the sum of two non-negative quantities so that each of them must be finite a.e. We start by defining (recall that an integral of a positive measurable function is always defined). The finiteness of the second term in (2.32) requires βL 0,λ f ∈ X 0 = L 1 ((x 0 , ∞), dx) so that we define We note that the second term in the definition of L β,λ is the same as in R(λ, T β ) since the extended term L 0,λ f enters the expression in a constant multiplier.
Denoting for a moment byT β andL β,λ the analogous extensions of the operatorT and its resolvent R(λ,T β ), respectively, we see thatL β,λ−b = L β,λ , λ >b, and hence we can repeat the proof of [9, Lemma 9.9] to show thatK ⊂T +B. However, as bothK andT differ from K and T , respectively, by the operator of multiplication byb, we obtain immediately that Further, we can establish the following results.
Lemma 2.6 (a) If f ∈ D(L β,λ ) then f ∈ L 1 ((x 0 , N )) for any N < ∞, L β,λ f ∈ C((x 0 , ∞)) and Proof. From the definition of D(L β,λ ) it is enough to consider f ≥ 0. Since L 0,λ f is finite almost everywhere and e λ /r is finite and non-zero for any finite x > x 0 , we see that, for any N ∈ (x 0 , ∞), x ′ x 0 e λR(y)+Q(y) f (y)dy < +∞ for some x ′ > N so that x ′ x 0 e λR(y)+Q(y) f (y)dy < +∞ . Since e λR(y)+Q(y) → 1 as y → x 0 , we arrive at the first assertion made in (a). The other two follow immediately.
To prove (b), we have, by (2.14), Proof. The proof hinges upon two results on the characterization of the generator (K, D(K)).
Next, consider the decomposition By Lemma 2.6(b), BL ′ β g = BL ′ β g ∈ X and thus L 0,λ g ∈ L 1 ((x 0 , N ), (1 + x)dx). Hence where each of the first three terms on the right-hand side is in L 1 ((x 0 , N ), (1 + x)dx) and the last two are both in X. Since Ku ∈ X, we can write where the limit on the right-hand side exists. Since the integral over (x 0 , N ) of each term within this limit exists, we can evaluate Using (2.14), we evaluate Hence N x 0 Let u 0 := [L 0,λ g]. Since u 0 = u − const e λ /r, we see that u 0 ∈ X. Thus µu 0 ∈ X and both ru 0 and (n − 1)au are in L 1 ((x 0 , ∞), dx). Furthermore, there exists a sequence N k → ∞ as k → ∞ for which r(N k )u 0 (N k )(1 + N k ) → 0. Indeed, otherwise r(x)u 0 (x)(1 + x) ≥ ǫ > 0 for some ǫ and all sufficiently large x. But then r(x)u 0 (x) ≥ ǫ(1 + x) −1 which would contradict the integrability of To deal with the last two terms in (2.37), we note that g enters the expression through a constant scalar multiplier and hence first we evaluate e −λR(y)−Q(y) dy, on using (2.14) with (2.29). Hence we obtain µ(y)L ′ β g(y)(y + 1))dy

20
Combining the above results, we see that there is a sequence (N k ) k∈N for such that which proves the thesis.

Global Solutions of the Growth C-F Equation
We now consider the combined coagulation and mass-growth fragmentation equation (1.3). Following the definition (1.2), we define the associated coagulation operator N on X by where, for u, v ∈ X, We recall that the coagulation kernel is a non-negative, bounded and symmetric function. Hence, for all u, v ∈ X, Similarly, N 2 [u, v] ≤ k L∞ u v . It follows that N (X) ⊂ X with N u ≤ 2 k L∞ u 2 for all u ∈ X . Moreover, for a given u 0 ∈ X, if u, v ∈ B(u 0 , ρ) := {w ∈ X : w − u 0 ≤ ρ} then 1) and so N is locally Lipschitz on X. The bilinearity of N also leads to from which we deduce that N is Fréchet differentiable at each u ∈ X, with Fréchet derivative Consequently, N u v ≤ q ρ,u 0 v , ∀v ∈ X, u ∈ B(u 0 , ρ) . Also, for u 1 , u 2 , g ∈ X, Hence, the Fréchet derivative is continuous with respect to u.
These results establish the existence of positive constants ρ 0 , t 0 and a strongly differentiable func- where X + = {u ∈ X : u ≥ 0 a.e. on (x 0 , x 1 )} and K = T + B is the infinitesimal generator of the growth-fragmentation positive semigroup {G K (t)} t≥0 ; see [13,Lemma 14].
To show that this local (in time ) solution is in X + for all t ∈ [0, t 0 ), we adopt the argument used in [13,Chapter 8]. First we note that the solution u of (3.2) is also the unique strongly differentiable for any α ∈ R. Hence u is the unique solution of the integral equation 4) where N α := N + αI .
Proof. This can be established by using the argument presented in the last part of the proof of [19,Theorem 4.3]. 2 To prove the global (in time) existence of a strict non-negative solution to (3.2) we shall apply Gronwall's inequality to establish that the local solution cannot blow up in finite time. The following preliminary lemmas will be required.
where C 1 is a positive constant.