Mean field limit with proliferation

An interacting particle system with long range interaction is considered. 
Particles, in addition to the interaction, proliferate with a rate depending 
on the empirical measure. We prove convergence of the empirical measure to the 
solution of a parabolic equation with non-local nonlinear transport term and 
proliferation term of logistic type.

1. Introduction. We are concerned with a macroscopic limit result for a system of particles which interact by a mean field potential and may proliferate; we want to prove the convergence of the empirical measure to a parabolic equation with nonlinear non-local transport term, like in classical mean field models, plus a growth term corresponding to the proliferation. Our starting motivation for this research came from Mathematical Oncology, where cells interact by random dynamics and proliferate, and one would like to discover appropriate macroscopic limits (PDE), for instance because a cell-level simulation of a tumor is too demanding, being involved a number of particles of the order of 10 9 .
1.1. The microscopic model. The microscopic system, defined on a filtered probability space (Ω, F, F t , P ), is composed of particles at positions X a,N t ∈ R d . We label particles by a multi-index a = (k, i 1 , ..., i n ) with i 1 , ..., i n ∈ {1, 2} and k = 1, ..., N , where N ∈ N is the number of particles at time t = 0. The particles already alive at time t = 0 are those with label a = (k), k = 1, ..., N . Their descendants require the additional labeling; setting (a, −) := (k, i 1 , ..., i n−1 ) if a = (k, i 1 , ..., i n ), we may say that a is a descendant of (a, −). Each particle "lives" only during a random time interval: particle with label a lives on In the random time interval I a,N , particle with label a interacts with all other living particles through a potential V ∈ C 2 c R d ; the dynamics of X a,N t is described by the gradient system where B a t are independent Brownian motions in R d and σ > 0. Each particle proliferates at a random rate λ a,N t which, in most applications, depends on the density of particles it has in the neighborhood. We prescribe a general structure of the form where the properties of the measurable functionals F N : will be specified below (M + R d is the set of Borel finite positive measures on R d ).
The precise meaning of the previous sentences is that we have standard Poisson processes N 0,a t , independent between themselves and with respect to the Brownian motions and initial conditions X R d the subset of the Sobolev space W 1,2 R d made of all non negative functions and by C β b R d , β ∈ (0, 1), the space of β-Hölder continuous functions; the β-Hölder seminorm of f will be denoted by [f ] β := sup x =y |f (x) − f (y)| / |x − y| β . We say that a map F : W 1,2 + R d → C β b R d satisfies the mild Lipschitz conditions if it is Lipschitz in the C b R d -norm and has linear growth in the C β b R d -norm, namely for every u, v ∈ W 1,2 The macroscopic limit result below requires a number of natural assumptions that we list now, plus the more critical assumption (5) that we emphasize in the statement of the theorem. Let {θ N } be a classical family of compact support smooth mollifiers of the form θ N (x) = −d N θ −1 N x ; let us introduce the mollified empirical measure (the theoretical analog of the numerical method of kernel smoothing) For technical reasons we assume sup N ∈N −d−2 N N < ∞; moreover, at least in the case of example (8)-(9) below, the role of θ N is auxiliary, it does not appear in the model, hence the restriction does not have biological relevance. Concerning the interaction potential V , assume that V ∈ C 2 c R d . Concerning the initial conditions, assume that u 0 ∈ L 2 R d ∩ L 1 R d and S N 0 , φ be convergent in probability to u 0 , φ , as N → ∞, for every φ ∈ C ∞ c R d . Moreover, assume the following technical condition: lim A simple and relevant sufficient condition for this property is given by Proposition 1 below. Finally, the definition of weak solution of the PDE (7) is given below in Section 3.1.
Theorem 1.1. Assume that, for some β ∈ (0, 1), there is a map F : which satisfies the mild Lipschitz condition (2)-(3) above and there is a sequence of positive real numbers {α N } converging to zero such that Assume moreover that there exists a constant C F > 0 such that Then the process h N t (x) converges in probability, as N → ∞, to the unique weak solution of the PDE The topologies of convergence of Having in mind Fisher-Kolmogorov-Petroskii-Piskunov equation and in general the concept of logistic growth, the most natural example of functional F (u) is F (u) (x) = (1 − u (x)), which means that proliferation rate decreases when we approach the threshold u = 1 (the value 1 is conventional). More precisely, due to the term div ((∇V * u t ) u t ), the solution u t may overcome any threshold, so one has to correct the classical logistic term and consider F (u) (x) = (1 − u (x)) + (proliferation is completely inhibited when the threshold is crossed). Our theory covers this case only in dimension d = 1, where W 1,2 1). In this case we take and assumption (5) is obviously satisfied; the others are elementary. In dimension d > 1 this example is not covered but we can treat the case where W is a Lipschitz continuous compact support probability density. In this case we simply take The validity of assumptions (2)-(3)-(5)- (6) in this case is shown in Lemma 3.10 below.
A variety of macroscopic limit results has been proved in the literature; let us quote, related to the present one, [11], [12], [14], [13], [8], [9], [10] and references therein, from which we have taken several elements of inspiration. However, a result of mean field type with proliferation, in the sense described here, is not treated in the previous references. Let us also mention a different class of macroscopic limit results, for instance [18], [17], which require very different techniques. In the more applied literature, two examples of works related to our problem are [1], [3].
A feature of our approach, shared by some of the previously quoted references, is that we do not use only results of tightness of measure-valued processes but of processes. A more distinguished feature, probably shared only by [9] (which however is very different), is that we use typical tools of the theory of stochastic Partial Differential Equations, for instance the tightness criterion used for stochastic Navier-Stokes equations in [6], [2], see also [5], [4].

1.3.
Motivations from mathematical oncology. Although the aim of this paper is mostly theoretical, we have been inspired by lectures about the emerging field of Mathematical Oncology in the choice of the problem and of some details. In this area, roughly speaking, models are classified as macroscopic, when described by partial differential equations, or microscopic, when described by stochastic ordinary differential equations or, even more often, cellular automata and other discrete stochastic models -in addition there are multiscale models with mixture of the previous two cases. Macroscopic models look at the tumor at the tissue level, microscopic ones at the cellular level. The link between the two descriptions is of interest for various reasons, in particular because a precise justification of the macroscopic models is difficult from general arguments based on fluid dynamics or mechanical, a biological tissue being something different. At the cellular level it is easier to be more realistic and thus results on the macroscopic limit of cellular systems is a way to justify or improve the macroscopic models; and also to have different interpretations of the constants appearing in the models -characterizing these constants is a major problem for the real applications. After these general comments, which clarify why we investigate the macroscopic limit, let us say that most models used nowadays in Mathematical Oncology are more complex than our one, since they involve different cell types -e.g. normoxic and hypoxic tumor cells, or cells with different genetic mutations, or cells of the extracellular matrix -molecular fields like oxygen and growth factors, and possibly objects related to the angiogenic cascade. See for instance [16] as an example of complex model. Our model, chosen as a starting point, captures only a few features of such complexity: i) the interaction between cells which may incorporate for instance a certain degree of repulsion resulting from the fact that cells cannot press each other too much -see the term called "crowding effect" in [16], different from our one but corresponding to a similar mechanism; ii) cancer cell proliferation, always present in each model.
Concerning proliferation in detail, most often in the literature of Mathematical Oncology it is taken of logistic form u (1 − u) (which corresponds to F (u) (x) = (1 − u (x)) + above). Such choice of F simply corresponds to the fact that cells proliferate better when the neighbor is not so crowded. However, other forms of F may be interesting as well. A phenomenon observed in vitro is that certain isolated cancer cells, even if embedded in a liquid that is rich of nutrients, do not proliferate: they need to adhere to other cells to proliferate. A functional F which charge -in the sense that decreases proliferation rate -not only the excessive presence of cells in the neighbor but also the opposite case, an excessive isolation, may be more realistic: cells which separate from the main cloud by random motion will continue their travel to meet blood vessels and lead to metastasis, but along the trip they do not proliferate so often as the cells close to the main tumor body. This, although vague, could be a motivation for investigating a general proliferation mechanism of the form λ a,N t = F N S N t , X a,N t above; more modeling work is necessary and is part of our research program.

Preparation. Starting from this section
to simplify notations. Let δ denote a point outside R d , the so called grave state, where we assume the processes X a t live when t / ∈ I a . Hence, whenever a particle proliferates and therefore dies, it stays forever in the grave state δ. In the sequel, the test functions φ are assumed to be defined over R d × {δ} and be such that φ (δ) = 0. Using Itô formula over random time intervals, one can With a few computations, one can see that the empirical measure S N t satisfies for every φ ∈ C 2 b R d and where The total relative mass FRANCO FLANDOLI AND MATTI LEIMBACH plays a central role. Since, in our model, the number of particles may only increase, we have the inequality that we shall use very often. The quantity S N T is, moreover, exponentially integrable, uniformly in N , see Lemma 3.11 below. Sometimes, however, having S N T in the inequalities could spoil properties associated to adaptedness, so we also introduce, for every R > 0, the stopping time (τ N R = +∞ if the set is empty). We also repeatedly use the identity which follows from Fubini Theorem. Other simple rules of calculus we often use are for every bounded measurable f : R d → R. Finally, we often use the inequalities which holds with a suitable constant C > 0. The bound on the first term comes from the assumption that θ has compact support and the assumption sup N −d−2 N /N < ∞; the bound on the second term is similar.
3. Tightness. In this section, about F N , we only use (6).
the left limit of f at t, Further, let (X t ) t≥0 be a stochastic process. We define its quadratic variation, when the limit exists and is independent of the partitions, where the maximal distance of two consequent sites in the partition {0 = t k 0 < t k 1 < · · · < t k n = T } converges to 0 as k → ∞. We denote the continuous part of the quadratic variation by [X] c , i.e. [X] In the following we need a lemma, also known as generalized Itô formula, see for example [15, page 245].
We use the generalized Itô formula to derive the following energy estimate.
Proof. For every x ∈ R d , we apply the generalized Itô formula to ((h(x)) 2 ) t≥0 and integrate over x ∈ R d : in ordinary Lebesgue integrals). It remains to understand the last line of the previous formula. At every jump time r, we have In the next lemmata, we generically denote by C > 0 any constant depending only on σ, T , Proof. By Hölder inequality we have and then we handle the second term by means of the bound (using also (11)) (14). The left-hand-side of the second inequality of the lemma is bounded above by (14) and (6).
Proof. Since we have (using the change of variable x → x − X a s in the Lebesgue integral over R d ) We conclude the first estimate of the lemma using Lemma 3.11 and (15). Similarly, since we deduce the second estimate of the lemma.
Proof. Obviously the expected value on the left-hand-side above is bounded by

FRANCO FLANDOLI AND MATTI LEIMBACH
Concerning the M 2,N s -term, we have is a martingale with respect to the filtration G t := F Λt , hence the last expression is bounded by Since the jumps of N a s and N b s , for a = b, never occur at the same time, we have Hence the last expression is equal to It is known that As above for g N s (y), using (15)  The result of the lemma follows by estimating S N s by R, being the integral in s only up to T ∧ τ N R . Proof.
Step 1. From lemmata 3.3 and 3.4 we have Hence, writing χ N R = 1 h N 0 2 L 2 ≤R , from Lemmata 3.5 and 3.6 we have

ds.
By Gronwall's lemma, we get, with C (R) := R 2 + 2C exp 2C R 2 + R + 1 , Moreover, for the same reasons, where C 1 (R) = R 2 /2 + C + C R 2 + R + 1 T C (R) and we have used the previous bound in the last term.
Step 2. For every R 1 > 0, the probability P sup t∈[0,T ] h N t 2 L 2 ≥ R 1 is bounded above by We may now find functions R 1 → R (R 1 ) such that lim R1→∞ R (R 1 ) = +∞, lim R1→∞ C (R (R 1 )) /R 1 = 0, where the function C(R) has been defined in step 1. We deduce, from assumption (4), that This proves one half of the claim of the corollary.
Step 3. For every R 1 > 0, by similar arguments P As above, we conclude that the second half of the claim of the corollary holds true.
In order to show tightness of the family of the functions {h N } N , in addition to the previous bound which shows a regularity in space, we also need a regularity in time. See the compactness criteria below. Proof.
Step 1. We need to estimate h N t − h N s 2 W −1,2 in such a way that it cancels with the singularity in the denominator at t = s. First, we have and thus by the Hölder inequality Notice that L 2 ⊂ W −1,2 with continuous embedding, namely there exists a constant C > 0 such that f W −1,2 ≤ C f L 2 for all f ∈ L 2 . Moreover, the linear operator div is bounded from L 2 to W −1,2 , namely div f W −1,2 ≤ C f L 2 and the operator ∆ is bounded from W 1,2 to W −1,2 , namely ∆f W −1,2 ≤ C div f W 1,2 . Therefore (we denote by C > 0 any constant independent of N , h N . , t, s) and now using (14), assumption (6) and t s ≤ T 0 in all terms,

FRANCO FLANDOLI AND MATTI LEIMBACH
Accordingly, we split the estimate of P dsdt > R in four more elementary estimates, that now we handle separately; the final result will be a consequence of them. The number C α = T 0 T 0 1 |t−s| 2α dsdt is finite, hence the first term is bounded by (renaming the constant C) and both these terms are, uniformly in N , small for large R, due to Lemma 3.11 and the estimates proved in the tightness part. The second addend, the one with Step 2. Concerning the martingale terms, we now prove that and the proof will be complete. For notational convenience, we abbreviate, for i = 1, 2, Note, that for every x ∈ R d the processes M 1,a (x) and M 2,a (x) are martingales. It follows, with computations similar to those of Lemma 3.6, for t ≥ s where in the last inequality we have used (15) and Lemma 3.11. Similarly, for the second martingale, A version of Aubin-Lions lemma, see [7], [6], [2], states that when E 0 ⊂ E ⊂ E 1 are three Banach spaces with continuous dense embedding, E 0 , E 1 reflexive, with E 0 compactly embedded into E, given p, q ∈ (1, ∞) and α ∈ (0, 1), the space L q (0, T ; E 0 ) ∩ W α,p (0, T ; E 1 ) is compactly embedded into L q (0, T ; E). We use this lemma with E = L 2 (D), E 0 = W 1,2 (D) and E 1 = W −1,2 R d , where D is a regular bounded domain, and with p = q = 2, α ∈ (0, 1/2). The lemma states that L 2 0, T ; W 1,2 (D) ∩ W α,2 0, T ; W −1,2 R d is compactly embedded into L 2 0, T ; L 2 (D) .
Notice that for αp > 1, the space W α,p (0, T ; E 1 ) is embedded into C ([0, T ] ; E 1 ), so it is not suitable for our purposes since we have to deal with discontinuous processes. However, for αp < 1 the space W α,p (0, T ; E 1 ) includes piecewise constant functions, as one can easily check. Therefore it is a good space for càdlàg processes. Now, consider the space Using the Fréchet topology on L 2 one concludes that L 2 0, T ; W 1,2 R d ∩ W α,2 0, T ; W −1,2 R d is compactly embedded into L 2 loc [0, T ] × R d (the proof is elementary, using that if a set is compact in L 2 0, T ; L 2 (B (0, n)) for every n then it is compact in L 2 loc [0, T ] × R d with this topology; see a similar result in [5]). Denoting by L ∞ w * 0, T ; L 2 R d and L 2 w 0, T ; W 1,2 R d the spaces L ∞ 0, T ; L 2 R d and L 2 0, T ; W 1,2 R d endowed respectively with the weak star and weak topology, we have that Y 0 is compactly embedded into Denote by Q N N ∈N the laws of h N N ∈N on Y 0 . From the "boundedness in probability" of the family Q N N ∈N , in Y 0 , stated by Corollary 1 and Lemma 3.7, it follows that the family Q N N ∈N is tight in Y . Hence, by Prohorov theorem, from every subsequence of Q N N ∈N it is possible to extract a further subsequence which converges to a probability measure Q on Y . We shall prove that every such limit measure Q is a delta Dirac Q = δ u concentrated to the same element u ∈ Y , hence the whole sequence Q N N ∈N converges to δ u ; and also the processes h N N ∈N converge in probability to u.
3.1. Auxiliary results. Let V , F satisfy the assumptions of Theorem 1.1 above.
Definition 3.8. Given u 0 ∈ L 2 R d ∩L 1 R d , u 0 ≥ 0, by weak solution of equation (7) we mean a function u ≥ 0 of class

FRANCO FLANDOLI AND MATTI LEIMBACH
such that for almost every t ∈ [0, T ] and for all φ ∈ W 1,2 R d .
Notice that, due to u ∈ L ∞ 0, T ; L 2 R d , and the assumption that ∇V is bounded and compact supported, we have (∇V * u r ) bounded, hence the first integral in the weak equation is well defined. Moreover, since F is uniformly bounded (see (6)), the last term is also well defined.

Remark 1.
If u is a solution in the sense of the definition, then there exists a (unique) re-presentative of u of class C w [0, T ] ; L 2 R d . Indeed, given φ ∈ W 1,2 R d , from identity (18) we deduce that the a.e. defined function t → u t , φ is a.s. equal to a continuous function g φ . Let u : From this it follows that t → u (t) , φ is uniformly continuous on Υ hence uniquely extendible to a continuous function t → L t (φ) on [0, T ]. From this it is easy to extract a re-definition of u (t) for t / ∈ Υ so that t → u (t) , φ is continuous on [0, T ]. Finally, it is not difficult to show that identity (18) holds for all t ∈ [0, T ] for the representative of class C w [0, T ] ; L 2 R d .
s v s , ∇v s ds the terms on the second line can be bounded by The terms on the last line, using assumptions (2)-(6), can be bounded by It is then sufficient to apply Gronwall lemma to deduce v = 0.
Step 2. To complete the proof, we check v ∈ W 1,2 0, T ; W −1,2 R d , which is needed to apply the chain rule. To this end, we check that any weak solution u of equation (7) given by the definition has this regularity property. The term u 0 is in W 1,2 0, T ; W −1,2 R d . The function r → (∇V * u r ) u r is of class L 2 0, T ; L 2 R d , since ∇V * u r is bounded, as remarked after the definition; hence the function r → div ((∇V * u r ) u r ) is of class L 2 0, T ; W −1,2 R d and thus the function t → t 0 div ((∇V * u r ) u r ) dr is of class W 1,2 0, T ; W −1,2 R d . The function r → ∆u r is of class L 2 0, T ; W −1,2 R d because u ∈ L 2 0, T ; W 1,2 R d and thus the function t → t 0 ∆u r dr is of class W 1,2 0, T ; W −1,2 R d . Finally the function r → F (u r ) u r is of class L 2 0, T ; L 2 R d , since F (u r ) is bounded, as remarked after the definition; hence the function t → t 0 div ((∇V * u r ) u r ) dr is of class W 1,2 0, T ; W −1,2 R d . It follows that the function But it easily coincides with the function t → u (t), see also Remark 1. Proof. First, it is obvious that F maps W 1,2 + R d (in fact even L 2 + R d ) into C β b R d for any β ∈ (0, 1) and (6) (2)-(3) are true for F by composition with the Lipschitz bounded function r → (1 − r) + (we also use the fact that W * u ≥ 0). Finally, which implies (5).
Lemma 3.11. There exists a γ > 0 such that sup N E e γ[S N T ] < ∞.
Proof. On the same probability space (Ω, F, P ) one can construct two processes, X a t ; a ∈ A N , which is equal in law to our particle system, and Y a t ; a ∈ A N , such that S (X) N t , 1 ≤ S (Y ) N t , 1 a.s. (we distinguish the objects of the different particle systems by adding "(X)" or "(Y )"), where the branching rate of Y a t ; a ∈ A N is constant and dominates the branching rate of the "X" system, The proof is an easy exercise, similar to a classical computation about the weak law of large numbers. 4. Passage to the limit. Given χ : [0, T ] → R of class C 1 , with χ T = 0 and given ψ ∈ W 1,2 R d , defined φ t := χ t ψ, we have Denote by Q N the law of h N t and assume a subsequence Q N k weakly converges, in the topology of the space Y defined by (17), to a probability measure Q. Given φ t := χ t ψ as above, the functional We shall prove that this lim inf is zero. It will follow Q (u : |Ψ φ (u)| > ) = 0. Since this holds for every > 0, we will deduce Q (u : Ψ φ (u) = 0) = 1. From this fact we shall prove the following result.
Lemma 4.1. Q is supported on the set of weak solutions of equation (7).
Proof. It remains to prove that lim inf k→∞ P Ψ φ h N k · > = 0, that the assertion of the lemma follows from the fact that Q (u ∈ Y : Ψ φ (u) = 0) = 1 for every φ of the form φ t = χ t ψ as above, and that Q is concentrated on non negative functions. We prove these claims in three different steps.
Step 1. We show that lim inf k→∞ P Ψ φ h N k · > = 0. Let us write N in place of N k for simplicity of notation. We have Using the identity satisfied by h N t against the test function φ t , we have where the constant C depends only on the quantities described below, before Lemma 4.2. In order to prove lim N →∞ P Ψ φ h N · > ε = 0, it is sufficient to prove the same result for each one of the previous terms.
We have lim N →∞ P C ( N + α N ) S N T 2 > ε = 0 from Chebyshev inequality and Lemma 3.11. The same applies to the terms C S N T β N and θ − N * φ 0 − φ 0 S N 0 . The term u 0 − S N 0 , φ 0 is obvious by the assumption of convergence in probability on S N 0 . The two martingale terms can be treated again by Chebyshev inequality and Lemma 4.4 below. Finally The first term goes to zero by Chebyshev inequality and Lemma 3.11. The second one goes to zero, as N → ∞, by Corollary 1, by a simple argument based on the definition of limit.
Step 2. We prove that the assertion of the lemma follows from the fact that Q (u ∈ Y : Ψ φ (u) = 0) = 1 for every φ of the form φ t = χ t ψ as above. Let {χ n } be a sequence of functions χ n : [0, T ] → R of class C 1 , with χ n T = 0, which is dense in L 2 (0, T ). Let {ψ m } be a dense sequence in W 1,2 R d . Set φ n,m t := χ n t ψ m ; we have Q (u ∈ Y : Ψ φ n,m (u) = 0, ∀n, m ∈ N) = 1. Let us prove that the set A := {u ∈ Y : Ψ φ n,m (u) = 0, ∀n, m ∈ N} is contained in the set of weak solutions.
If u ∈ A, then Ψ χtψ m (u) = 0 for every m ∈ N and every Lipschitz continuous χ : [0, T ] → R with χ T = 0. The proof of this claim by approximation with the sequence χ n t is not difficult and we omit the details. Given u ∈ Y , there exists a set Υ ⊂ [0, T ) of full measure such that for every t 0 ∈ Υ and every m ∈ N. Given t 0 ∈ Υ, take the new sequence χ n t defined (at least for n large enough) as χ n t = 0 for t > t 0 + 1 n , χ n t = −1 for t < t 0 , χ n t = −1 + n (t − t 0 ) for t ∈ t 0 , t 0 + 1 n . We have We deduce, from Ψ χ n ψ m (u) = 0, ∀n, m ∈ N, that identity (18) holds at time t 0 , for all φ = ψ m , m ∈ N. Therefore it is true, for each φ = ψ m , a.s. in time. By density of {ψ m } and the regularity properties of u it is easy to deduce that (18) holds for every φ ∈ W 1,2 R d .
Step 3. By the Portmanteau Theorem, the weak convergence of Q N k to Q implies that Q (A) ≥ lim sup Thus, Q is concentrated on non-negative functions, which completes the proof.
In the next three lemmata we denote by C any constant depending only on σ 2 , T , C F , D 2 V ∞ , θ