Stochastic branching at the edge: individual-based modeling of tumor cell proliferation

An individual-based model of stochastic branching is proposed and studied, in which point particles drift in R¯+:=[0,+∞)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{\mathbb {R}}_{+}:=[0,+\infty )$$\end{document} toward the origin (edge) with unit speed, where each of them splits into two particles that instantly appear in R¯+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{\mathbb {R}}_{+}$$\end{document} at random positions. During their drift, the particles are subject to a random disappearance (death). The model is intended to capture the main features of the proliferation of tumor cells, in which trait x∈R¯+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x\in \bar{\mathbb {R}}_{+}$$\end{document} of a given cell is time to its division and the death is caused by therapeutic factors. The main result of the paper is proving the existence of an honest evolution of this kind and finding a condition that involves the death rate and cell cycle distribution parameters, under which the mean size of the population remains bounded in time.


Introduction
One of the most natural applications of semigroups of bounded positive operators in L 1 -like spaces [1,2,16] is the description of the evolution of probability densities, widely used in population biology, genetics, medical sciences, etc. Among the processes which have multiple applications, one might distinguish branching [6]. In this paper, we propose and study an individual-based model that can describe the proliferation of tumor cells, cf. [15,Sect. 2.2]. Here, 'individual-based' means that each single member of the population is taken into account in an explicit way. This is in contrast to macroscopic models yielding a 'coarse-grained' picture, where populations are described in terms of aggregate parameters such as density, cf., e.g., [13,18]. In the proposed model, each member of a finite population of particles (assuming cells) is characterized by a random trait x ∈ [0, +∞)-time to its branching (fission or division). Then, one of the basic acts of the dynamics is aging-diminishing of the traits with unit speed. At point x = 0, the particle divides into two progenies with randomly distributed traits-the second basic act of its dynamics. Finally, during the whole lifetime each particle is subject to a random death, which in the case of tumor cells is caused by therapeutic factors. The main questions concerning this model which we address here are: (a) can one expect (and under which conditions) that the population dynamics is honest; (b) what might be a condition for the boundedness in time of the population mean size. The mentioned honesty of the dynamics means that the population remains in time almost surely finite (no explosion). In the language of semigroups of positive operators, cf. [1,2], this means that the semigroup of operators mapping the initial probability distribution of the population traits on those corresponding to t > 0 is stochastic. The mentioned boundedness condition ought to involve the death rate and cell cycle distribution parameters. Its practical meaning might be estimating at which level of the therapeutic pressure the tumor cell population stops growing ad infinitum. This aspect of the theory is indeed practical since, for various kinds of tumors, the cell cycle distributions and their parameters are known, see [4,15] and also [5,14,17,19].
The rest of the paper is organized as follows. In Section 2, we recall necessary mathematics and then introduce two models, of which the second one is the principal model mentioned above. The introduction of this model is preceded by a careful investigation of its 'mild' version, in which the particles instead of fission just disappear at the edge. This turns useful in the subsequent study of the principal model. In Section 2, we also formulate the main result-Theorem 1. Its proof is performed in Section 3, whereas concluding remarks are placed in Section 4.

The model and the result
We begin by providing necessary notions and facts. Then, we introduce an auxiliary model, the advantage of which is that it is soluble. This allows us to clarify a number of properties of the principal model introduced afterward. Next, we formulate the result as Theorem 1.

Preliminaries
In this work, we use the following standard notations: R + = (0, +∞),R + = [0, +∞), N 0 = N ∪ {0}, N stands for the set of positive integers. For a Banach space, (E, · E ), with a cone of positive elements, E + , a C 0 -semigroup S = {S(t)} t≥0 of bounded linear operators S(t) : E → E is called sub-stochastic (resp. stochastic) if, By Γ , we denote the set of all finite subset ofR + . Its elements are finite configurations. This set is equipped with the weak topology, see [3], which is metrizable in such a way that the corresponding metric space is separable and complete. Namely, a sequence, {γ n } n∈N ⊂ Γ , is convergent in this topology to some γ ∈ Γ if is measurable if and only if there exists a family of symmetric Borel  functions Note that f (n) are defined up to their values at points of coincidence x i = x j . However, this makes no problem as we are going to deal with L 1 -like spaces, elements of which are defined up to sets of (Lebesgue)-measure zero.
To simplify our notations, in expressions like γ ∪ x, x ∈R + we consider x as a single-element configuration {x}. The Lebesgue-Poisson measure λ on (Γ, B(Γ )) is defined by the integrals holding for suitable measurable f : Γ → R. Such integrals have the following property which we will use throughout the whole paper.
Here, we just renamed the variables and used the assumed symmetry of the function.
Let h : Γ → R + be separated away from zero and such that holding for all f ∈ X := L 1 (Γ, dλ). Then, we set X h = L 1 (Γ, hdλ) and equip it with the norm defined in (2.4). By X + , X + h we will denote the cones of positive elements of X and X h , respectively. Note that where → denotes continuous embedding. Clearly, X h and X + h are dense in X and X + , respectively.
For a given n ∈ N, by W 1,1 (R n + ) we denote the standard Sobolev space [9], whereas W 1,1 s (R n + ) will stand for its subset consisting of all symmetric u, i.e., such that u(x 1 , . . . , x n ) = u(x σ (1) , . . . x σ (n) ) holding for all permutations σ ∈ Σ n . Remark 1. It is known, cf. [9, Theorem 1, page 4], that each element of W 1,1 s (R n + )-as an equivalence class-contains a unique (symmetric) u :R n + → R such that (a) for Lebesgue-almost all (x 1 , . . . , x n−1 ), the mapR + y → u(y, x 1 . . . , x n−1 ) is continuous and its restriction to R + is absolutely continuous; (b) the following holds In the sequel, we will mean this function u when speaking of a given element of Then, k u ∈ W 1,1 (R + ). As mentioned above, see (2.1), each measurable f : Γ → R defines symmetric Borel functions f (n) :R n + → R, n ∈ N 0 . Let f ∈ X be such that each f (n) belongs to the corresponding W 1,1 s (R n + ). Set (2.8) Definition 1. By W we denote the subset of X consisting of all those f for which f (n) ∈ W 1,1 s (R n ) and the following holds Note that the key point here is the convergence of the series.
For γ ∈ Γ , set γ t = {x + t : x ∈ γ }, i.e., γ t is a shift of γ . Obviously, γ t ∈ Γ for t > 0. In the sequel, we will use such shifts also with negative t in the situations where x + t ≥ 0. By (2.8) and (2.9), for f ∈ W, we have For f ∈ W, we then set Proof. The continuity of the embedding W → X follows by (2.11). Take f ∈ X and ε > 0. Then, pick n ε ∈ N 0 such that, see (2.4), Let f (n) be its limit, which exists as W 1,1 (R n ) is complete. Clearly, f (n) is symmetric, i.e., f (n) ∈ W 1,1 s (R n ). Since { f m } m∈N is a Cauchy sequence in X , it converges there to some f ∈ X such that its f (n) are the limits of the sequences { f (n) m } m∈N as just discussed. By Remark 1, these f (n) satisfy (2.6); hence, for all N ∈ N, the following holds Let us then show that the sequence { f N } N ∈N is bounded, and thus f lies in W. Set C = sup m D f m . Then, for each m ∈ N, by the triangle inequality we have that The second inequality holds for a fixed N and m > m N for an appropriate m N . This yields the proof of the first part of the statement. Then, the closedness follows by the fact that (2.11) is exactly the graph norm of (D, W).
Note that the action of the semigroup S 0 can be visualized in the form

A soluble model
As mentioned above, our principal model is a modification of another model, the main advantage of which is that it is in a sense soluble. We will use this fact in studying the principal model below. This soluble model describes the following process. A finite cloud of point particles is distributed overR + . Each particle in the cloud moves toward the origin with unit speed, and disappears at x = 0. The states of the cloud are probability measures on X , which we assume to be absolutely continuous with respect to λ defined in (2.2).
Let f t be the density (Radon-Nikodym derivative) of the state at time t, and f be the density of the initial state. As described above, the cloud undergoes the evolution according to the following formula Here, and in the sequel, by 1 we denote the corresponding indicator. Let us show that { f t } t≥0 has the flow property To this end, we write express f s in the right-hand side of (2.14) by (2.13), and then obtain, see (2.3), Let |γ | denote the number of points in γ ∈ Γ . For f as in (2.13), we define Note that N 1 is just the expected number of points in the cloud as time t = 0. Note also that N 0 = 1 in view of the normalization of f . Let us prove that where the last equality follows by Then, the sum in the last line of (2.16) has only one nonzero term corresponding to ξ = γ t 1 . That is, that yields (2.15). The latter yields also N 0 (t) = 1, i.e., the map f → f t defined in (2.13) preserves the norm. Notably, (γ t 2 ) −t is the part of the initial cloud that remains after time t, shifted toward the origin. Its expected cardinality thus cannot be bigger than that of the initial cloud, that is reflected in the latter estimate. Proof. Clearly, it is enough to prove the statement for positive f only. By (2.13), we have (2.17) Here, χ(ξ ) = 0 whenever ξ = ∅, and χ(ξ ) = 1 otherwise. By Corollary 1, we have that I 1 (t) → 0 as t → 0 + . To estimate I 2 (t), we proceed similarly as in deriving (2.16). That is, To proceed further, we write .
As in the proof of Proposition 2, we may take f m such that f (n) m = 0 for n > n m with sufficiently large n m . Fix some m > m ε and then write which by (2.3) yields in (2.19) the following . As all such k u are continuous, one may pick t ε so that I 3 (t) < ε/2 for t ≤ t ε which by (2.19) yields I 2 (t) < ε for such t. This completes the proof.

Corollary 2.
There exists a unique stochastic semigroup, S 0 = {S 0 (t)} t≥0 , on X such that, for each f ∈ X , f t = S 0 (t) f , where f t and f are the same as in (2.13). The semigroup S 0 leaves invariant each X l , l ∈ N.
Proof. The semigroup property of S 0 follows by (2.14). Its strong continuity follows by Proposition 3, whereas the property S 0 (t) : X l → X l follows by (2.15).
For each f ∈ W and λ-almost all γ ∈ Γ , we know that the map x → f (γ ∪ x) is continuous onR + and absolutely continuous on R + , see Remark 1. Set (2.20) Note that, for each x ∈ R + and f ∈ V, with an appropriate C f > 0, independent of x. Indeed, by (2.9) we have Then, the proof of (2.21) follows by the definition of V and the triangle inequality. For f ∈ V, let f (n) be as in (2.1) and k f (n) , see Remark 1, be as in (2.7) for this f (n) . In view of (2.21), one can define is the expected number of points with traits in the interval [a, b] in the corresponding state. A priory k f need not be integrable on the wholeR + . Define Clearly, L 0 : V → X , in view of which we introduce the following norm In the statement below, we will use the set V consisting of all those f ∈ V which have the following two properties: (a) for each Let us prove that where the closure is taken in · V . For f ∈ V and m ∈ N, let f m be such that f All the three terms of the right-hand side are the remainders of convergent series, that eventually yields (2.26). Proof. First we prove that, for all f ∈ V, it follows that Clearly, it is enough to show this for f ∈ V + := V ∩ X + only. Similarly as in (2.17) and then (2.18), for such f we obtain for some τ ∈ [0, t), see (2.10). Then, to prove (2.27) it suffices to show that holding for positive f ∈ V , see (2.26). By (2.28) we then have (2.30) To estimate J (t), we proceed as follows, cf. (2.16), that can be obtained for a positive f ∈ V similarly as (2.21). Since γ → f (γ ∪ 0) is in W, the first term in the first line of (2.30) also disappears in the limit t → 0 + , which finally yields (2.29).
Let us prove now that V = W. By Corollary 2 (semigroup property and strong continuity) and by (2.27), it follows that f t is differentiable in t at all t ≥ 0 whenever f ∈ V. Then, f t ∈ W, see (2.10). Let us prove that also f t ∈ V in this case. Indeed, by (2.13) and (2.21) for f ∈ V + , we have from which we then obtain By Corollary 2, we know that f t = f , which by (2.31) implies and hence, see (2.24),  ). Then, the graph of (L 0 , W) is closed in the graph norm, which yields the closedness and hence the whole proof.

The model
Our principal model is a modification of the soluble model just described. Its main new aspect is that each particle by reaching the edge divides into two progenies with randomly distributed traits x, y ∈R + . In addition, we assume here that the particles can disappear (die) at random also outside of the origin. The Fokker-Planck-Kolmogorov equation corresponding to this our model is defined by the Kolmogorov operator L that has the following form, cf. (2.24), Here, m ≥ 0 is the mortality rate and b is a symmetric probability density which hereby has the property Our assumption concerning the cell cycle probability density is that holding with some σ ≥ 3 and b * > 0. Then, (2.34) with L given in (2.35) describes a drift of the particles toward the origin (with unit speed) subject to a random death that occurs at x ∈R + with constant rate m. At the origin, the particle produces two progenies whose initial traits (times to their division) are random. According to (2.35), the dynamics of the considered model is characterized by the following competing processes: (a) disappearance of the existing particles at the edge x = 0 and due to the mentioned random death; (b) appearance of new particles in the course of division. It is quite clear that, for m = 0, the branching is supercritical and thus the population will grow ad infinitum. Among our aims in this work is to find a trade-off condition for these two processes that secures the boundedness in time of the population mean size.
As mentioned above, our model is intended to capture the basic aspects of the dynamics of a population of tumor cells consisting in the following: (a) malfunctioning of regulatory mechanisms and hence uncontrolled proliferation with random cycle length; (b) increased mortality caused by therapeutics; (c) death occurring at random with no inter-cell dependence. In the model, aspect (a) corresponds to the independent division with random cycle length, for a given particle equal to its trait x at the moment of its appearance. Aspects (b) and (c) are taken into account in the second and third terms of L, see (2.35). The choice of the model parameters is based on the following reasons: (a) we believe that the therapeutic effect on a cell is nearly independent of its age (phase of mitosis); (b) b(x, y) is often modeled as the product of two Γ -densities x k e −αx , cf. [5,19], which clearly satisfies (2.38). See also [4] for more on cell cycle modeling and Section 4 for further comments.
Let us now define L as an operator in X . Set  3),  For positive ς and α, we set Next, assuming (2.38) holding with σ ≥ 3, we introduce The proof of this theorem will be performed in Sect. 3. Here, we make some comments to its results. The condition m > m 1 -under which the described evolution is honest and hence the system of cells remains almost surely finite-is only sufficient. At this moment, it is unclear how far it may be from the necessary one. The last part of Theorem 1 yields a balance condition between the disappearance of the particles and the appearance of their progenies. Indeed, the expected number of particles at time t is N 1 (t), see (2.15). By (2.45) and Theorem 1, we then have (2.47) and thus N 1 (t) remains bounded if the mortality rate m is bigger than a certain quantity, explicitly computable in terms of the cell cycle parameters, see (3.6). Another conclusion of this sort is that the evolution described by Theorem 1 is honest, cf. [1], since the norm of f t is preserved, i.e., f t = 1. This, in particular, means that the system of particles remains almost surely finite for all t > 0. Indeed, since f t is the Radon-Nikodym derivative of the state at time t, the fact that f t < 1 would mean that the population is finite with probability strictly less than one, and hence the estimate in (2.47) holds provided the system is finite. Since (perhaps) Reuter's seminal paper [12], the evolution of this kind is called dishonest. More on the honesty theory can be found in [1]. The backward Kolmogorov equation is dual to (2.34) in the sense that Here, F t : Γ → R is an observable and that additionally clarifies the nature of the dynamics described by L and its dual L * .

The proof
The proof will be divided into two parts. First we construct a C 0 semigroup S = {S(t)} t≥0 such that the solution in question is obtained in the form f t = S(t) f 0 , for all t ≥ 0 and initial f 0 belonging to the domain of the generator of S. A special attention here will be paid to proving that S is stochastic. In the second part, we prove the stated boundedness that implies (2.47).

The stochastic semigroup
The construction of the mentioned semigroup S is based on a perturbation technique, developed in [16], and some aspects of the honesty theory [1,10]. Its adaptation to the present context is given in the following three statements. Therein, we deal with a Banach space E equipped with a cone of positive elements, E + , that have the following property. There exists a positive linear functional, ϕ E , such that u E = ϕ E (u) whenever u ∈ E + . Thereby, the norm · E is additive on E + .  In this case, T 1 (t)u E < u E , that is, the evolution is dishonest. In order to establish the honesty of T 1 , one has to get additional information on its properties. The first statement in this direction is a simple consequence of Theorem 3.5 and Corollary 3.6 of [1].

Proposition 6. The semigroup T 1 mentioned in Proposition 5 is honest if and only if its generator is the closure of (A + B, D A ).
A more specific fact-applicable in L 1 spaces-is provided by the following statement. Proof. The operator (A, D(A)) generates a substochastic semigroup, S 0 , with see ( see (2.24) and (2.36). Since B is positive, this yields that, for each ε ∈ (0, 1) and f ∈ D + (A), the following holds ϕ(L ε f ) ≤ 0. Then, (L ε , D(A)), see (3.1), generates S ε as stated, and the semigroup S is obtained in accordance with Proposition 5.

Lemma 2.
Let m 1 be as in (2.46). Then, for m > m 1 , the semigroup S constructed in Lemma 1 is generated by the closure of (A + B, D(A)) and hence is honest therefore.
Proof. Here, we employ Proposition 7. To this end, we introduce v ∈ D(A) by the following expression where φ σ is as in (2.37 that yields v ∈ D(A). Thus, to apply Proposition 7 we have to show that For γ = ∅, both terms on the left-hand side of (3.2) vanish. For γ = {x}, Now by (2.38), we obtain Then, If b(x, y) is such that b * ≤ 2σ , we take m 1 = 0 and obtain (3.2). For b * ≤ 2σ , we use the fact that φ σ +1 (x) ≤ φ σ (x), x ≥ 0, and then get from (3.3) the following where the latter inequality holds in view of the assumed m ≥ m 1 , see (2.46).

The boundedness
To prove the boundedness which yields (2.47), we are going to employ another statement of [16]. Thus, in the context of Proposition 5 we further impose the following.

Assumption 1.
There exists a linear subspace, E ⊂ E, which has the following properties: (ii) There exists a norm, · E , on E that makes it a Banach space and the embedding E into E is continuous. (iii) E + := E ∩ E + is a generating cone in E. The norm · E is additive on E + and hence there exists a linear functional, ϕ E , such that u E = ϕ E (u) whenever u ∈ E. (iv) The cone E + is dense in E + . To this end, we employ Proposition 8, where as E we take X h ς,α = X h m . Note that the latter means that they are equal as sets, and that m > m 1 is positive even if m 1 = 0, see (2.46). Clearly, X h ς,α has all the properties as in Assumption 1, cf. (2.5). Moreover, in this case By direct inspection, one checks that both (a) and (b) assumed in Proposition 8 are satisfied for this choice of E. Now we prove that (c) also holds for m ≥ m 2 where the latter has to be found. For f ∈ D A ∩ X + h ς,α , we write, cf. (2.48), (2.49), where ϒ(ς, α) = ς − 1 +β(α), Since β is integrable, by the Riemann-Lebesgue lemma it follows thatβ(α) → 0 as α → +∞. Then, one finds α > 0 such thatβ(α) < 1. For this α, ς = 1 −β(α) is positive and thus can be used in (3.4), where it yields ϒ(ς, α) ≤ 0. Now, for this α and m 1 as in (2.46), we set m 2 = min{m 1 ; α}, (3.6) that, for m ≥ m 2 , yields LHS(3.4) ≤ 0, which then by Proposition 8 yields in turn This completes the proof of Theorem 1 with m 2 defined in (3.6).

Summary and concluding remarks
We begin by making a brief summary of the aspects of the theory presented here, understandable also for non-mathematicians, see also [8] for an extended explanations and analysis. Then, we discuss some aspects of this work, as well as outline its possible continuation.

Summary
In cancer biology, it is well established that cancer cells proliferate wildly by repeated, uncontrolled mitosis. 'Unlike normal cells, cancer cells ignore the usual density-dependent inhibition of growth...piling up until all nutrients are exhausted 1 ' Therefore, to model populations of cancer cells one can use 'particles' that undergo independent branching into two new 'particles' after some random time. Being unharmed populations of such 'particles' grow ad infinitum since the branching number is two. Therapeutic pressure causes disappearance of some of them from the population before branching, the result of which may be restricting the population growth. The effect of the treatment is proportional to its intensity and to the mean length of the inter-mitosis period, during which it acts. Then, the paramount problem of modeling of such populations is to find qualitative relations between the treatment intensity and probabilistic parameters of the cell cycle processes in a given population. In this work, we find this relation in the form m ≥ m 2 with m 2 > 0 defined in (3.6) and (2.46).

The model and its study
The proposed model seems to be the simplest individual-based model that takes into account the basic aspects of the phenomenon which we intended to describe: (a) essential mortality caused by external factors and independent of the interactions inside the population; (b) randomly distributed lifetimes of the population members, at the end of which each of them branches into two progenies; (c) branching independent of the interactions inside the population. The main difficulty of its mathematical study stemmed from the presence of the gradient in the Kolmogorov operator L in (2.35), which is typical for transport problems [10]. A more general version of the proposed model instead of the last summand in (2.35) could contain that corresponds to branching into a 'cloud' ξ with possibly random number of progenies. In fact, this might be done in the present context at the expense of a modification of the bound in (2.38). Our choice was motivated by the reasons of simplicity and practical applications-mitosis with two progenies. Noteworthy, in our model the lifetimes of siblings are in general dependent as random variables. The independent case would correspond to the choice b(x, y) = β(x)β(y) with β as in (3.5). Note, however, that the definition of m 2 in (3.6) remains the same in this case. In order to take into account also dependence like 'parent-progeny,' cf. [4], one would make the trait more complex by including the corresponding parameter. For example, instead of R + one may takeR 2 + consisting of pairsx = (x, y) in which x is still time to fission whereas y is responsible for the mentioned dependence. This additional trait can be used to model, e.g., further mutations of the tumor cells.

The practical meaning
As mentioned above, we believe that the proposed theory can have direct practical applications for the following reasons. There exists a rich literature on modelingparameter fitting including-of various types of cancer, see, e.g., [4,5,17,19] and the sources quoted in these publications. This means that, in a concrete situation, one can calculate m 2 by means of (3.6) and (2.46), which then can be used to estimate the corresponding therapeutic dose.

Further development
Along with the modifications of the model already mentioned above in this section, we plan to consider also its version describing infinite populations. Here, we plan to employ methods developed in [7], of which studying finite populations is a part. We also plan to develop a mesoscopic theory of this model by means of scaling techniques and Poisson approximations, see also [7]. This would allow for connecting the microscopic theory developed in this way to a description based on aggregate parameters, similar to that is [13,18].