A singular limit for an age structured mutation problem

The spread of a particular trait in a cell population often is modelled by an appropriate system of ordinary differential equations describing how the sizes of subpopulations of the cells with the same genome change in time. On the other hand, it is recognized that cells have their own vital dynamics and mutations, leading to changes in their genome, mostly occurring during the cell division at the end of its life cycle. In this context, the process is described by a system of McKendrick type equations which resembles a network transport problem. In this paper we show that, under an appropriate scaling of the latter, these two descriptions are asymptotically equivalent.


1.
Introduction. An important problem related to mutations is to understand how a particular trait spreads in a population. One of the simplest ways to model this is to describe the change in the sizes of subpopulations having this trait. It can be done by the standard balancing argument: the rate of change of the number of, say, cells with a particular genome γ is equal to the rate of recruitment of cells with other genomes that change, due to mutations, to γ, minus the rate at which the mutations cause the cells with genome γ to move to subpopulations with another genome. This balance equation typically is supplemented by terms describing the death and proliferation of the cells. Models of this type were extensively studied, see for instance [5,11,15,16,23], either in the context of the development of drug resistance in cancer cells, or in describing the micro-satellite repeats. In the former, a population of cells was divided into smaller subpopulations according to the number of the drug resistant gene, while in the latter the feature of interest was the number of the micro-satellite repeats. In both cases the authors only allowed for changes between neighbouring populations, which resulted in the birth-and-death system with proliferation, u ′ 0 = a 0 u 0 + d 1 u 1 , u ′ 1 = a 1 u 1 + d 2 u 2 , u ′ n = a n u n + b n−1 u n−1 + d n+1 u n+1 , n ≥ 2, where u n is the number of cells belonging to the subpopulation n, n = 0, 1, 2 . . . (e.g. having n drug-resistant genes), d n+1 , b n−1 are the rates of recruitment from the populations n + 1 and n − 1 into the population n and a n is the net growth rate of the population n which incorporates birth, death and loss to other populations of cells of type n. We note that the particular form of (1), in which the first equation is decoupled from the rest, is due to the assumption adopted in op.cit. that any object in the state 0 can generate only objects in the same state. This assumption, however, is of no significance in our considerations. It is clear that, in principle, jumps between arbitrary populations can occur and thus there is no need to restrict our attention to tridiagonal matrices. We can consider a general model where u = (u i ) i∈N and L = (l ij ) i,j∈N , N ⊆ N may be infinite and L is a positive off-diagonal matrix, where the off-diagonal term l ij is the rate at which the cells are recruited into the population i, characterized by a particular genotype, from the population j. The diagonal entries represent the loss rates from corresponding classes. At the same time it is recognized that the cells have their own vital dynamics that should be taken into account if a more detailed model of the evolution of the whole population is to be built. Also, the mutations can be divided into various groups. Here, we distinguish two types of mutations: those that are due to the replication errors and occur when the cell divides, and others, due to external factors (mutagenes), that may happen at any moment of the cell's life cycle. This results in a model of the form where we consider a population of cells described by their density u j (x) in each class j ∈ N, and where x is a parameter describing the maturity of the cell, typically its age or size. In general, cells in each class can divide at a different age (or size), say τ j , so that we normalize the age of division to 1 by introducing the maturation Cells are born at x = 0 and divide at x = 1 producing, due to mutations during mitosis, daughter cells of an arbitrary class, the distribution of which is governed by a nonnegative matrix K = (k ij ) i,j∈N . Though typical mitosis produces two daughter particles, this is not always the case. For instance, cancer cells can produce up to five daughter cells, [24], so that we shall not place any further restrictions on K. Further, M = diag(µ j ) j∈N gives the death rates due to external causes and R = (r ij ) i,j∈N describes redistribution of cells caused by mutations due to external factors (mutagenes). Finally, the nonnegative vectorů describes the initial population. A classical example of this type is the discrete Lebowitz-Rubinow-Rotenberg model in which the cells' populations are distinguished precisely by the maturation velocities, [21].
The main objective of this paper is to determine under what conditions can solutions to (3) be approximated by the solutions of (2). There could be several ways to approach this problem and the answer may be not unique. Our approach is to assume that the maturation velocities are very large or, in other words, the cells divide many times in the reference unit of time, while the deaths and mutations due to external causes remain at fixed, independent of the maturation velocity, levels.
To balance the fact that there is a large number of cell divisions in the unit time, we assume that the daughter cells have a tendency to be of the same genotype as the mother (see e.g. [18, p. 19]), which is represented by splitting the boundary operator as K = I + ǫB, ǫ ≪ 1. Thus, we will consider the singularly perturbed problem where I is the identity matrix. By l 1 N we denote R N equipped with the l 1 norm if N is finite or l 1 (the space of absolutely summable sequences) if N is infinite (thus in the latter case we identify N with N). If N = N, we assume that all matrices in (4) represent bounded operators from l 1 N to l 1 N . Alongside (4), we consider the simplified problem We will be working in X = L 1 ([0, 1], l 1 N ). We define A ǫ as the realization of the differential expression ; u(0) = (I + ǫB)u(1)}, see [12]. In line with the interpretation of the model, we assume that I + ǫB ≥ 0. Further, A 0,ǫ is the realization of A 0,ǫ = diag{−ǫ −1 v j ∂ x } j∈N on the same domain. We assume that if N is infinite, then there are numbers v min and v max such that If N is finite, the fact that for each ǫ > 0 the operator (A 0,ǫ , D(A ǫ )) generates a semigroup, denoted by (e tA0,ǫ ) t≥0 , follows from [7, Theorem 3.1]. However, the argument used in op.cit is purely norm dependent and can be repeated for countable N as long as B is bounded and (6) is satisfied, see also [12]. Then the boundedness of M and R allows for the application of the Bounded Perturbation Theorem, [14, Theorem III 1.3], which gives the generation of (e tAǫ ) t≥0 by (A ǫ , D(A ǫ )). The paper is organized as follows. In Section 2 we show the convergence of the resolvents of (4) as ǫ → 0 + and thus, by the Trotter-Kato theorem, the convergence of semigroups (e tAǫ ) t≥0 to the age independent semigroup generated by VB − M + R, but only for initial values that are constant for each j ∈ N. Such a convergence is often referred to as a regular convergence of semigroups, [9]. This result is not fully satisfactory as often singularly perturbed semigroups, possibly corrected by initial or boundary layers, converge to the limit semigroup for all initial conditions, see [4,8]. In Section 3 we show that in this case, in general it is impossible to achieve such a result. However, if we restrict our attention to the macroscopic characteristic of the model; that is, the observable total size of each subpopulation j, obtained from the solution to (5), then the situation changes. In Section 4 we prove that if the maturation times in each subpopulation are natural multiples of some fixed reference time, then, for arbitrary initial valuesů, the population sizes obtained from the solution to (5) by integration over (0, 1), converge to the solution to (2) that emanates from the initial condition 1 0ů (x)dx. This extends a similar result from [8] obtained for velocities independent of j ∈ N.
In the general case, we rewrite the main formulae from the proof of [7, Theorem 3.1], specified for (5). First, we solve where c ǫ = (c ǫ,j ) j∈N is an arbitrary vector and E ǫλ (s) = diag e If ǫ is fixed, E ǫλ (1) can be made arbitrarily small by choosing large λ. Then c ǫ is uniquely defined by the Neumann series and hence the resolvent of A 0,ǫ exists. Since I + ǫB ≥ 0, we can work with f ≥ 0. Adding together the rows in (9) we obtain which is equivalent to the standard norm by (6). Then, by integrating (8), we obtain see [7, Theorem 3.1] for details. The case when b j ≤ 0 for all j ∈ N leads to a contractive semigroup for each ǫ > 0 and the equiboundedness is obvious. Otherwise, as in op. cit., we can assume that each b j is nonnegative. Then the generation is achieved by the application of the Arendt-Batty-Robinson theorem [2,3]. However, to prove equiboundedness of (e tA0,ǫ ) t≥0 we need to understand how the Hille-Yosida estimates are obtained in this theorem. Using, for instance, the proof given in [3,Theorem 3.39], a densely defined resolvent positive operator T on a Banach lattice X generates a positive semigroup if there are λ 0 > s(T ) (the spectral bound of T ) and c > 0 such that for any nonnegative x ∈ X we have R(λ 0 , T )x ≥ c x . In the proof one defines S = T − ωI, where s(T ) < ω ≤ λ 0 and then the Hille-Yosida estimate for S is obtained as Hence, we have to show that s := sup ǫ>0 {s(A 0,ǫ )} < ∞ (and thus, by (11), we can take c = λ −1 for any fixed λ > s) and that for some ω > s the family { R(ω, A 0,ǫ ) } ǫ>0 is equibounded. For this, let us return to (9) and estimate the Neumann series for provided λ > B v max . Next, using l'Hôspital's rule, we find and hence ǫ (I − (I + ǫB)E ǫλ (1)) −1 ≤ L 1 for some constant L 1 (depending on λ), yielding for some constant L 2 . Let us fix ω > v max B . Using (8) and (14), we see that there is a constant L, independent of ǫ such that Then using (11), the equivalence of the norms · v and · and R(λ, A 0,ǫ − ωI) = R(λ + ω, A 0,ǫ ), we write (12) as Theorem 2.2. If λ > ω + Q Lω −1 , where ω and L are defined in (15) and (16), in the uniform operator topology.
Proof. By the previous proof, the resolvent of A 0,ǫ is given by where We observe that for any α ∈ [0, 1] the matrix E ǫλ (α) has the following Taylor expansion where, using the integral form of the reminders, we find First we see that and hence the last term in (19) converges to zero as ǫ → 0 + . Next, using the second equality in (21) with α = 1, we have for λ > v max B Next, using the first equation in (21) we see that

A SINGULAR LIMIT FOR A MUTATION PROBLEM 7
Finally, treating c → E ǫλ (x)c as the operator from l 1 N to X, we estimate This actually shows that the operators converge in the uniform operator norm. Combining all estimates and using the projection operator P we find that From (16) and (15) we have, in particular, for some fixed ω > v max B and arbitrary λ > ω. Since Q = −M + R is bounded, and the series converges uniformly in ǫ for λ > ω + Q Lω −1 . Since the operators B, V, Q are independent of x, they commute with P and, by P 2 = P, we have (QR(λ, VB)P) n = R(λ, VB + Q)P.
Proof. According to the version of the Trotter-Kato approximation theorem given in [10,Theorem 8.4.3], if the resolvents of the generators of an equibounded family of semigroups (strongly) converge to an operator R λ , then R λ is the resolvent of the generator of a semigroup on the closure of its range. Here, the limit operator is the resolvent of the bounded operator VB + Q composed with the projection P : X → l 1 N . Thus, restricted to l 1 N , the range of R(λ, VB + Q)P equals the range of R(λ, VB + Q) which is l 1 N . Hence, (27) follows from the Trotter-Kato theorem.

A counterexample.
This result is not very satisfactory. In asymptotic theory, [4,9,8], the convergence obtained in Corollary 1 is referred to as regular convergence. However, typically it is possible to extend the convergence to initial data from the whole space, albeit at the cost of losing the convergence at t = 0, or adding necessary initial or boundary layers. The following example shows that this is impossible to achieve in our context. Consider the scalar problem where b ∈ R. It is easy to see that if u ǫ (x, t) = [e tA0,ǫů ](x), then for t satisfying nǫ ≤ t < (n + 1)ǫ, so that n = ⌊t/ǫ⌋, we have Then where we used The above is obvious for t = 0 and for t > 0 it follows from n n + 1 ≤ t ǫ ǫ t ≤ 1 and n → ∞ with ǫ → 0. At the same time, let us consider t = 1 and ǫ = 1/k. Then we obtain

A version of irregular convergence.
As demonstrated in Section 3, we should not expect the convergence of (e tA0,ǫ ) t≥0 for initial conditions which are not constant for each population. However, under certain additional assumptions we can derive a stronger version of convergence than that of Corollary 1, which is of relevance to the problem at hand. In fact, we are interested in approximating the solutions to (4) (or (5)) by solutions to the system of ODEs (2). The solution of (2) gives the total size of each subpopulation, while the solution to (4) (or (5)) gives the density in each subpopulation. The total number of individuals in each subpopulation can be calculated by integrating the density. Hence, we can ask how well the total population in each subpopulation, calculated from the solution of (4) ((or (5)) can be approximated by the solution to (2). Using the notation of Section 2, the problem can be expressed as follows: does hold for some some matrix H?
In this section we shall focus on problem (5) and adopt the assumption from [17] that the speeds v j , j ∈ N, are linearly dependent over the field of rational numbers Q or, in other words, In our interpretation of the model, this corresponds to the situation that the maturation times τ j = 1/v j in each subpopulation are natural multiples of some fixed reference maturation time τ = 1/v. We observe that (6) implies that the set of different velocities is finite. Condition (31) allows the problem to be transformed into an analogous problem with unit velocities. Such a transformation appeared in [17] (and in a more detailed version in [20]) in the context of transport on networks and thus, even though (5) is not necessarily related to the network transport, see [6], its interpretation as a network problem allows for a better description of the construction.
Using the graph theoretical terminology, we identify the jth subpopulation with an oriented edge e j parametrized by x ∈ (0, 1), with ageing corresponding to moving from the tail of e j , associated with 0, to the head associated with 1, and denote G = {e j } j∈N . Let Q be the graph with the set of vertices, V (Q), equal to G and the adjacency matrix I + ǫB. We note that, in general, Q is not the line graph as I + ǫB allows for the construction of the vertices connecting the edges e j so that G becomes a graph only in very special cases, [6].
To proceed with the construction, first we re-scale time as τ = vt and, for each j ∈ N, we re-parameterize each edge e j by y = l j x. This converts (5) into a problem with the unit velocity for each j ∈ N, but with the equations defined on (0, l j ). However, since each l j is a natural number, we subdivide each interval (0, l j ) into l j intervals of unit length. The kth subinterval in (0, l j ) is identified with the edge e j,k . This creates from G the set G 1 = {{e j,k } 1≤k≤lj } j∈N consisting of, say, M edges. Then, as in the previous paragraph, we define Q 1 to be the graph with V (Q 1 ) = G 1 . Its adjacency matrix must take into account the connections between the old edges e j , determined by I + ǫB, and the connections across the points subdividing the old edges. More precisely, consider a function f on G that is continuous on each e j . This function is transformed into a function ϕ on G 1 , the components of which are required to have the same values at the head of e j,k and the tail of e j,k+1 , k = 1, . . . , l j − 1, j ∈ N. Then it is easy to see that the adjacency matrix of Q 1 satisfying these conditions is given by Here, T j is a l j × l j dimensional matrix which is either a scalar 1 if l j = 1, or a cyclic matrix, and C ij is an l i × l j matrix, given, respectively, by Summarizing, we converted (5) into Since (32) has the same structure as (5), there is a semigroup (e tA0,ǫ ) t≥0 solving it.
To be more precise, the above construction defines an operator f → ϕ = Sf , where ϕ = (φ j,s ) j∈N,1≤s≤lj , f = (f j ) j∈N and, for s ∈ {1, . . . , l j }, It is easy to see that the inverse f = S −1 ϕ is defined by By direct calculation, see also [20], S : X → X := L 1 ([0, 1], l 1 M ) is an isomorphism such that we have the similarity relation The motivation behind (35), see [12] and [20, Proposition 4.5.1], is the fact that for any ϕ ∈ X with (e tvA0,ǫ ϕ)(x) = ϕ(x − vt/ǫ) for vt/ǫ ≤ x < 1. Let l = lcm{l j } j∈N . We observe that T l = I and For any K ⊂ N, let P K be the projection from L 1 ([0, 1], l 1 K ) to l 1 K , given by (17). The next result is not strictly necessary but it relates operators on G to those on G 1 and introduces in a natural way the projection Π that plays an essential role in the main theorem. where Proof. The estimate (13) carries over to this case by (35), which also gives and, by Theorem 2.2 and the continuity of S, To find the limit resolvent in terms ofC, we modify the calculations from Theorem 2.2. The only different part is related to the calculation of where here E ǫλ (s) = e −ǫλs I. We have, for λ > v max B , Hence, as in (22) with V = l −1 I, Since T is block diagonal with a finite number of finite dimensional blocks T j (different from 1), we can use the finite dimensional theory to claim, by [19, p. 633 by periodicity. This, combined with (40), ends the proof. Using Proof. Let us write f ∈ X as where P N f ∈ l 1 N and P N w = 0. By Corollary 1 and the continuity of P N , lim ǫ→0 + P N [e tA0,ǫ P N f ] = P N e tVB P N f = e tVB P N f . Hence, by linearity, it suffices to show that lim ǫ→0 + P N [e tA0,ǫ w] = 0,