Evolution of a structured cell population endowed with plasticity of traits under constraints on and between the traits

Confronted with the biological problem of managing plasticity in cell populations, which is in particular responsible for transient and reversible drug resistance in cancer, we propose a rationale consisting of an integro-differential and a reaction-advection-diffusion equation, the properties of which are studied theoretically and numerically. By using a constructive finite volume method, we show the existence and uniqueness of a weak solution and illustrate by numerical approximations and their simulations the capacity of the model to exhibit divergence of traits. This feature may be theoretically interpreted as describing a physiological step towards multicellularity in animal evolution and, closer to present-day clinical challenges in oncology, as a possible representation of bet hedging in cancer cell populations.


Introduction: biological background
One of the main theories explaining the origins of cancer states that the hallmark capabilities of cancer are based on latent functions already existing in the genome of normal human cells, and that cancer represents a reversion to a less differentiated and less cooperative cellular behavior. This theory is usually called the atavistic model of cancer (Davies and Lineweaver 2011). It is opposed to the more commonly admitted somatic mutation theory (Lineweaver and Davies 2020), that states that cancer originated in a single cell, following a catastrophic sequence of stochastic tumorigenic mutations, and from which Darwinian selection produced by divisions from this unicellular basis a more or less organised society of cheating cells, more and more escaping controls by the host organism, i.e., a cancer. In the atavistic theory, accompanied or induced by blockade of differentiation or reverse differentiation of normally maturing cells, societies of cells in a multicellular organism (cancer is always a disease of multicellular organisms) somehow, in some location of the organism, escape the fine control under which they are normally placed and revert to a previous, coarse and disorganised state of multicellularity (Davies and Lineweaver 2011). This may be understood as a process of "deDarwinisation", through which cancer cells gain a state of plasticity (Clairambault 2020;Shen and Clairambault 2020;Hanahan 2022;Yuan et al. 2019) representative of a former state in the evolution of multicellularity.
The atavistic theory has thus deep connections with the theory of evolution to multicellularity. The passage from unicellular organisms to multicellular ones led to the regulation of capabilities, resulting in controlled proliferation and differentiation of cells leading to specialisation and cooperation between specialised cells. The role of environment-driven cellular stress in this process of specialisation has recently been stressed by various authors (Nedelcu and Michod 2020;Wagner et al. 2019). The new genes responsible for these regulations became tumour suppressors. The atavistic theory states that if these new suppressors become damaged for some reason, then latent genes, associated with functions from unicellular organisms, will reappear and dominate the scenery, thus resulting in the unconstrained proliferation and the lack of cooperation with the other cells of the host organism, as actually found in tumours.
Understanding how evolution led to the emergence of multicellularity then becomes a problem closely related to that of understanding cancer. Firstly because unravelling in detail what is lost (cohesion) and what is gained (plasticity) in this reverse evolutionary process may help understand the nature of cancer from a functional point of view. Secondly because it is reasonable to assume that Darwinian selection in tumours starts from a primitive state of multicellularity in which cells are very plastic with respect to their phenotypes, which sends us back to states in evolution close to the emergence of multicellularity. The division of work through cell differentiation achieved by coherent multicellularity (i.e., designing a stable, cohesive multicellular organism) (Müller 2001) was of vital importance in order to evolve into more complex and more functionally efficient organisms. In a first stage, this differentiation was very likely reversible, due to the high plasticity that cells were endowed with in a primitive state of multicellularity. Under these conditions, one can reasonably assume that these primitive organisms adopted bet hedging strategies, i.e., common risk-diversifying strategies in unpredictably changing and often aggressive environments, in order to maximise their phenotypic fitness (Philippi and Seger 1989;Seger 1987).
In more detail, bet hedging in cell populations may be defined as a diversification of phenotypes in a cohesive (or at least bound by a common fate community, usually of genetic nature) cell population to optimise its fitness, in particular a minima to ensure its survival in a life-threatening environment, in other words to design a fail-safe strategy for the preservation of a propagating element; one cell may be enough to achieve this goal. Enhancing the ability to (quickly) diversify phenotypes by changing differentiation paths (by reversal, i.e., vertical de-differentiation, or by horizontal trandsdifferentiation in a differentiation tree) may in particular be achieved at the chromatin level by means of epigenetic enzymes, molecular instances of cell plasticity at work (Shen and Clairambault 2020).
Among such commonly described strategies of living organisms (unicellular or multicellular) meant to ensure survival in changing environments have been classically described fright, fight and flight. Fright (or freeze) is not likely to induce phenotype evolution. Fight (establishing barriers, secreting poisons, gathering in colonies) and flight (motility to escape unbeatable predators) can. Differentiation between somatic and germinal cells is also a major step in evolution. Bet hedging strategies were not only present at elementary stages of evolution. They are a common adaptive tool that can still be found in nature at different levels of complexity, from prokaryotic organisms to vertebrate ones. In between, tumour cells, thanks to their high plasticity, in the presence of an aggressive environment provided by immune response of the host body or of any anti-cancer treatment, may adopt bet hedging as a strategy to guarantee a prolonged survival of their colony. The wide presence of bet hedging in nature as an evolutionary mechanism, and its many links to the development of cancer is what motivates us in the present attempt towards a mathematical model representing some of the factors that influence this phenomenon (natural selection, epimutations and environmental stress).
In our modelling choices, we have privileged, for the sake of biological interpretation, as in Chisholm et al. (2015), two phenotypes that are often identified as such in theoretical ecology: viability (potential to resist a deadly insult: the elephant strategy) and fecundity (potential to proliferate, or in a way, surviving by numbers, even under hard environmental conditions: the rat strategy). They may influence cell behaviour in opposed directions, as two different choices, incompatible for the same cell, in the same way as fecundity and motility are notably incompatible (cells that divide do not move, and vice versa). To take into account the faculty of rapidly changing phenotypes in case of a life-threatening insult, a capacity reported about many cancer cells (e.g., epithelial to mesenchymal transition or drug-induced drug persistence) (Shen and Clairambault 2020), here we have added plasticity as a complementary structuring phenotype (or trait) of cells. Among questions at stake we will in particular deal with is the optimisation of fitness strategy by concentration of phenotypes, either in two or more ranges surrounding fixed points (multimodality of traits, i.e., bet hedging) or around a central optimal point (unimodality of trait, i.e. no bet hedging when it is not favourable).

The model
We consider a population (not necessarily of tumour cells) in which each individual has three defining traits: viability associated with the variable x ∈ [0, 1] which reflects the potential to resist deadly insults, fecundity associated with the variable y ∈ [0, 1] representing the potential to proliferate and plasticity associated with the variable θ ∈ [0, 1] which represents the potential to continue to differentiate (or de-differentiate, or transdifferentiate, as long as it has not been fixed at its lower bound θ =0) within a differentiation tree. We assume furthermore that for a certain regular function C : R 2 → R and a positive constant K , (x, y) ∈ := {C(x, y) ≤ K }, so z = (x, y, θ) ranges over the set D := { × [0, 1]}. We then consider the evolution problem (1), (2), (3) on the density of population n = n(t, z) ≥ 0.
In the above equation, chosen for the sake of simplicity as diagonal, the matrix gives the speed at which non-genetic epimutations occur, otherwise said it is a minimally simple representation of how the internal plasticity trait θ affects the non-genetic instability of traits x and y, by tuning the diffusion term; the function represents the sensitivity, meant here as the force of external evolutionary pressure, of the population to abrupt changes in the environment; stands for the total amount of individuals in the population at time t. From a biological point of view, the matrix A(θ ) represents random epimutations leading to non genetic changes, increasing with plasticity θ , in the traits viability x and fecundity y, while V stands for an external factor (abrupt biophysical changes in the ecosystem, exposure to life-threatening drugs in the case of cancer cell populations) inducing cellular stress in the population. Throughout the paper, we refer to ∇ · (V n) as advection term, or drift term, indistinctly.
We chose the term d(z)ρ(t) here as the simplest way to represent a death term in the evolution of the population. This multiplicative representation is classic, but most of the results we show in what follows can be extended to more general reaction terms. Throughout our work we assume (H1) For some p > 2, the initial population density n 0 (z) belongs to L p (D). (H2) The intrinsic growth rate in absence of competition r (z) and the death rate d (z) due to competition for individuals with trait z are positive bounded functions that satisfy 0 < r − ≤ r (z) ≤ r + and 0 < d − ≤ d(z) ≤ d + . (H3) The diffusion parameters a 11 (θ ), a 22 (θ ) and a 33 are strictly positive, with a 11 (θ ) and a 22 (θ ) being non decreasing with respect to θ . Hence, the matrix A(θ ) is elliptic for all values of θ . (H4) The function V (t, z) is continuously differentiable for all values of t > 0 and z ∈ D.
Under these hypotheses we are able to prove the existence and uniqueness of a solution for the problem (1-3) using the finite volume method in order to obtain a convergent semi-discrete scheme.
The problem (1-3) underlying hypotheses (H1) to (H4) sets a structured population model of evolution that, as mentioned before, takes into account some of the factors that might lead to the occurrence of "bet hedging". Amongst the first works in this topic, we can find (Cohen 1966), where is studied the fraction of seeds that germinate and the fraction that remains dormant, in order to maximise the long term expectation of growth. In the same work, the similarities of this phenomenon with economic decision making under risk (so-called "fail-safe strategies") are noted. Other early works on the subject are Philippi and Seger (1989) and Seger (1987). The reaction term in (1) is a simple way of modelling the selection principle. This term and more general ones are used in Perthame (2006) in order to study some basic properties of structured populations undergoing this type of behaviour. The same reaction term is also used in Pouchol and Trélat (2018), where are analysed the global asymptotic stability properties for integro-differential systems of N species structured by different sets of traits. A similar competition term is used in Desvillettes et al. (2008) to provide results about the long time behaviour of such reaction models. The diffusion term here models non-genetic instabilities (also known as epimutations), which constitute the drift of phenotype without alteration of the genotype. In Perthame (2006) an integral operator in order to model mutations arising during reproduction is used, and something similar could be done for the epimutations. The second order operator used in (1) can be obtained then after re-scaling the time variable. The effect of this phenomenon in cancer development is discussed in Huang (2013) from a biological point of view. Epimutations can also occur because of external stress, and this is represented in (1) by means of the advection terms. A biological example for a population changing phenotype due to external stress can be found in Matuła et al. (2019), where the effect of physical stress on the shape and the cell wall thickness of E.coli bacterias is discussed. Two different models are used in Chisholm et al. (2015) to conclude that the three mechanisms described above might reversibly push an actively-proliferating, and drug-sensitive, cell population to transition into a weakly-proliferative and drug-tolerant state, which will eventually facilitate the emergence of more potent, proliferating and drug-tolerant cells. One of such models is an integro-differential model very similar to (1-3), but without including the effect of plasticity on the evolution of the population and without assuming the existence of a constraint between the traits x and y.
The main results of this paper involve the variational formulation of (1-3), which we now introduce. Denote H = L 2 (D), with (·, ·) H the usual scalar product in that space, and V = H 1 (D) with ·, · = ·, · V ×V being the duality product in V. For any given n 0 ∈ H , T > 0, we say that is a variational solution of the problem (1-3) if it is a solution in the following weak sense where for any ϕ ∈ X T . We say that n is a global solution if it is a solution on [0, T ] for any T > 0.
We focus on giving a proof for this theorem using a discretised version of problem (1-3) after applying the Finite Volume Method to it. For this purpose, we define a set D h ⊃ D, that can be covered by the union of N disjoint cubic cells, denoted as D j , of side length h. After integrating the Eq. (1) over each of the cells D j we derive the system of first-order differential equations where ν j is an approximation of the average of the solution n(t, z) over D j , N j is the set of indexes corresponding to the neighbours of D j and the coefficients M j and B jl are functions of V (t, z), A(θ ), r (z) and d(z). A full detailed derivation of the scheme is given in Sect. 3. We can then introduce the following result involving the solution for this system: Theorem 2 For all non-negative n 0 ∈ L p (D), p > 2, there exists a unique nonnegative solution for problem (5)(6). Furthermore, the functionñ h (t, z) defined bỹ converges in L 2 (D T ) to the unique non-negative weak solution of (1-3) as h goes to zero.
The existence and non-negativity of the solution for (5-6) results from the Cauchy-Lipschitz theorem, while the convergence ofñ h (t, z) is the consequence of the compactness of the sequence. The existence result in Theorem 1 will be treated in Sect. 3, but the uniqueness can be directly obtained from the variational formulation with the help of some a posteriori estimates. Let us proceed then to prove the uniqueness before addressing the existence. Assume that there exists a non-negative weak solution n for (1-3). Assume as well that there exist 0 ≤ t 0 < t 1 such that ρ(t) ≤ max{ρ (0) ))1 t∈[t 0 ,t 0 +ε] + 1 t>t 0 +ε . This leads to the equality Taking the limit when ε goes to 0, we obtain that, for all t > t 0 where ρ(t) = D n(t, z)dz is a continuous function due to the fact that n(t, z) ∈ C([0, T ], H ). We can write the previous relation as which, thanks to the hypothesis on ρ(t) over (t 0 , t 1 ), leads to for all t ∈ (t 0 , t 1 ). This represents a contradiction, and implies that for all values of t.
Taking now ϕ = n on the variational formulation, using standard arguments to bound the linear terms from Q[n] and the estimate (7) for the non-linear part together with the Gronwall Lemma, we obtain the relation for some real positive number a. Using Gronwall's lemma, this relation implies that for all values of t. Finally assume the existence of two non-negative weak solutions n 1 and n 2 for the same initial data n 0 . Taking the difference between their respective variational formulations, choosing ϕ = n 1 − n 2 , we get the equality Once again, the linear part of Q[n 1 ] − Q[n 2 ] can be easily bounded using standard methods, leading to

H ds
Consequently, thanks to Gronwall's lemma, n 1 − n 2 2 H = 0 for all t, therefore, the solution is unique. As stated before, the proof of existence of a weak solution will be carried on in Sect. 3. Following the ideas from Carrillo et al. (2020), the semi-discrete numerical scheme (5-6) is developed and the convergence of its solution to the solution of (1-3) is demonstrated. A fully discrete numerical scheme is obtained starting from the semidiscrete one and its convergence is also proved. Section 4 is devoted to the numerical simulations, starting with some numerical computations of the approximation error by comparing the results with an exact solution. Finally, the solutions of some examples with biological meaning is presented.

Existence of a weak solution and numerical approximation
We aim to use the Finite Volume method in order to find a sequence of problems whose solutions converge to the solution of (4).

Preliminaries on the finite volume method
Consider h = 1/M where M is a natural number and define the mesh

Now introduce the sets
For simplicity, we define N := |M| as the amount of elements in M, and denote each of its elements as D j , for j = 1, . . . , N . For each D j , we denote its centre of mass as z j := (x j , y j , θ j ). For each j, define N j as the set of indexes l such that D l and D j share a common boundary. Denote such common boundary as jl , its centre of mass as z jl and n jl the outer normal vector of D j , in the direction of D l . We remark that the distance between the centres of two neighbouring cells D j and D l will be equal to |z l − z j |. Having cubes as the mesh cells guarantees that the condition is fulfilled. It is important to remark that if D has a regular enough boundary (for example: smooth or polygonal), then the area of D h \ D, which we denote as |D h \ D|, will converge to zero as h vanishes. The approximated problem (1-3) is then given by We propose a classical finite volume method based on local averages of the unknown density over cell grids defined by Assume that the coefficients and the solution from Eq. (9) are smooth. Then, integrating it over a cell D j yields the equality After integrating by parts and using the boundary conditions (10), we get For a real function f (t), define the positive and negative part of f as respectively. Using an upwind approximation technique for the advection term, we conclude where u jl (t) = V (t, z jl ) · n jl . On the other hand, we have As A(θ ) is a diagonal matrix and n jl is either one of the vectors from the euclidean canonical base, or one of their opposites, we have the relation where A jl = (A(θ jl )n jl · n jl ). Notice that, thanks to hypothesis (H3), there exists α > 0 such that A jl ≥ α, for all j and l. So that, together with the expression of the normal vectors, (8) this implies Due to the approximation of the gradient Taking into account that the reaction terms can be easily approximated where we have adopted the notation r j := r (z j ) and d j := d(z j ). Finally, using again (12), we get Consequently, collecting (13), (14) and (15), and getting rid of the approximation orders we obtain the semi-discrete scheme where This system of equations can be complemented with the set of initial data where n 0 (z) is the extension by 0 of n 0 (z) to all of R 3 .

Global existence, uniqueness, positivity and boundedness of the solution for the semi-discrete scheme
We prove the local existence and uniqueness of the solution for the problem (16-17), by using the Cauchy-Lipschitz theorem. Then, such solution is proved to be nonnegative and, as a consequence, bounded independently of t. Finally the boundedness property is used to prove the global existence of the solution.

Proposition 3 (Local existence of solution) For all sets of initial data {ν
Proof The RHS term in (16) is Lipschitz continuous for all values of t and ν j , therefore the existence and uniqueness of solution over a certain interval [0, T * ) is a direct consequence of the Cauchy-Lipschitz theorem.
On the other hand, consider a strictly positive set of initial values ν 0 j and define the continuous function f (t) = min j ν j (t). If f (t) ≥ 0 for all t < T * , then the solution remains positive at all times. If f (t) < 0 for some t ∈ (0, T * ), then there exists t 0 > 0 such that f (t 0 ) = 0 and f (t) ≥ 0 for t < t 0 . This implies the existence of j 0 such that ν j 0 (t 0 ) = 0 with ν j 0 (t 0 ) ≤ 0. If ν l (t 0 ) > 0 for some l ∈ N j 0 , then, thanks to (16) we have which is a contradiction with the previously established fact that ν j 0 (t 0 ) ≤ 0. Consequently ν l (t 0 ) = 0 for all l ∈ N j 0 . Furthermore, from the definitions of f (t) and t 0 , we also have that ν l (t 0 ) ≤ 0 for all l ∈ N j 0 . We can iterate the previous argument in order to obtain that ν j (t 0 ) = 0 for all j and consequently, thanks to the uniqueness of the solution, ν j (t) ≡ 0 for all j and all t. This is a contradiction with the assumption that f (t) is negative for some value of t and consequently f (t) ≥ 0 for all t.
Finally, for non-negative initial values ν 0 j and ε small enough, we define ν ε j as Thanks to the previous step, the solution of (16) associated to ν ε j remains non-negative for all t and all ε. Hence, thanks to the continuous dependence of the solution of a system with respect to its initial data, we conclude that ν 0 where r − , r + , d − and d + are the bounds given in (H1) for r (z) and d(z) respectively.
Proof Multiplying (16) by h 3 for each j, adding up all the equations, recalling that The non-negativity of the solution implies then These differential inequalities directly imply the bounds overρ h (t).
Notice that the L 1 bounds are independent of t. Additionally, the upper bound also implies that  (18) which, thanks to the Lipschitz continuity of the right hand side of (16) allows to extend to solution to a certain interval [T h , T * h ), contradicting this way the maximality of T h .

Discrete gradient, L 2 norm estimate and compactness result
In this section we introduce some piecewise constant functions depending on the solution of (16-17) together with some estimates related to such functions. Then, some compactness properties will be proved in order to ensure that such functions converge to some function that will be proved to be a weak solution of (1-2), and their derivatives, respectively. We first introduce n h (t, z) defined as Notice that n h L 1 =ρ h (t). Now, for each l ∈ N j , define the polygonal subsets of D j , denoted D jl , having jl as the common side and z j as a vertex. The subsets D jl are pyramids of area s jl = This function can be regarded as a discrete gradient for n h (t).
Proposition 6 (L 2 bound) For each value of h, define the space H h := L 2 (D h ). Then, there exists positive constants a and b, independent of h, such that the functions n h (t, z) and v h (t, z) satisfy the following estimate Proof Multiplying Eq. (16) by h 3 ν j (t) for each j, and adding them up, we obtain the relation where Let us proceed to estimate each of these terms. In order to simplify the notation, we define and write Consequently, knowing that −u + jl = u − l j ≤ 0, we have where we used that |z l − z j | = 2|z jl − z j |. From the definition of u jl , we conclude that Then, for all ε > 0, Young's inequality implies On the other hand The ellipticity and boundedness of the matrix A(θ ) imply that there exist positive constants α, α such that α ≤ A jl ≤ α for all j and l. From this and the value of s jl we deduce that and consequently notice that the last identity is true since v h is a piecewise constant function. Using the bounds over r (z) and the positiveness of d(z) andρ h (t) we see that Using (21), (22) and (23) in (20), with the relation which, after taking ε = α 3V and using Gronwall's lemma, leads to the estimate (19).
The result of Proposition 6 can be easily generalised if instead of multiplying each Eq. (16) by h 3 ν j (t), we multiply by h 3 ν satisfy the following estimate Before stating the compactness result that will allow us to extract a convergent subsequence from n h we give two important results that are a consequence from estimate (19).
Proof The proof of this proposition follows the same steps as Lemma 18.3 and Remark 18.8 in Eymard et al. (2000). In particular, thanks to the discrete trace inequality given in the Lemma 10.5 within the same reference, we can prove the existence of a constant C, independent of h and η, such that and we get the result from Proposition 8 thanks to (19).
From now on we introduce the notation D T := (0, T ) × D.
Proposition 9 There exists a constant C 1 , independent of h, such that for all ψ ∈ D(D T ), Proof Using the scheme (16), we have that, for all ψ ∈ D(D) ψdz is the mean value of ψ over D j . Reordering the terms in this equality we have We can bound each term on this equality as follows: From the definition of u jl (t) we get that |u ± jl | ≤ V L ∞ (D T ) and consequently Using the previously established relation A jl ≤ α we get Finally, using the boundedness of r (z), d(z) andρ h we see that Putting everything together, we obtain the existence of a constant C independent of h, such that Finally, using the inclusions H 3 (D) ⊂ C 1, 1 2 (D) ⊂ C 1 (D) that hold true in any open subset of R 3 with smooth enough boundary thanks to the Sobolev inequalities, integrating over (0, T ), using the Cauchy-Schwartz inequality and estimate (19), we get to the result stated in the proposition.
Proposition 10 If n 0 ∈ L p (D) for some p > 2, then there exists a function n ∈ L 2 (0, T ; V) such that, up to the extraction of a sub-sequence, the sequence of functions n h strongly converges to n in L 2 (D T ) and v h weakly converges to ∇n in L 2 (D T ).
Proof Propositions 8 and 9 give all the necessary tools in order to ensure the relative compactness of the set of functions {n h } h in L 2 (D T ). In order to prove that such set is indeed compact, we will use a variant adapted to our purposes of the proof of Moussa (2016, Theorem 3) to show that {n h } h is a Cauchy sequence in L 2 (D T ). We consider a sequence of mollifiers ε (z) = ε −3 (ε −1 z) for a positive, symmetric function ∈ D(B(0, 1)) satisfying R 3 (z)dz = 1.
Step 1: We claim that lim ε→0 ε * n h − n h L 2 ((0,T )×R 3 ) = 0, uniformly on h. We have , thanks to the Cauchy-Schwarz inequality and the value of the integral of ε (y). Consequently and due to Proposition 8, we obtain the strong convergence of ε * n h to n h in L 2 ((0, T ) × R 3 ), uniformly in h.
Step 2: We prove that for every fixed ε and any compact ω ⊂ D, the sequence ε * n h is uniformly bounded in H 1 ((0, T ) × ω). Thanks to Young's inequality, for a fixed ε, Furthermore, for any compact ω ⊂ D, ψ ∈ D((0, T ) × ω) and ε small enough, we have ε * ψ ∈ D(D T ) and consequently, using Proposition 9 and Young's inequality for the convolution product, we get Then, by duality, for a fixed ε, d dt ε * n h is uniformly bounded in L ∞ ((0, T ) × ω) and we conclude that for each fixed ε, ε * n h is uniformly bounded in H 1 ((0, T ) × ω).
Step 4: The sequence n h is a Cauchy sequence in L 2 (D T ).
Taking q > 1 and using Hölder's inequality we have Hence, for q = p/2, thanks to (25) goes to 0 when m goes to infinity. Consequently, for all η > 0, we can fix m big enough so that the second term in (26) is smaller than η/2 and then, thanks to Step 3, we can choose h 1 and h 2 in such a way that the first term is also smaller than η/2, which proves that {n h } h is a Cauchy sequence in L 2 (D T ).
Step 5: Let n be the limit of a suitable subsequence of n h and v := (v 1 , v 2 , v 3 ) the weak limit of v h (such limit exists due to the fact that v h is bounded in (L 2 (D T )) 3 ). We claim that n ∈ L 2 (0, T ; V) and v = ∇n. For all φ(t, z) ∈ C([0, T ]; D(D)), the convergence of n h to n implies that where Therefore, we obtain D n h (t, z)∇φ (t, z) Recalling that z jl is the center of mass of jl , we may estimate the integral over jl as as done in (12) with the integral over D j . Using the fact that | jl | = h 2 and the amount of cells on the mesh is Reordering the terms of this last sum, we have where This function strongly converges in L 2 (D T ) to φ. Substituting (29) in (28) and then in (27), and using the weak convergence of v h to v, we obtain that, for all φ ∈

The separability of D(D) finally implies that
As v(t, z) ∈ L 2 (D) for almost all t ∈ [0, T ], n(t, z) belongs to V for almost all t ∈ [0, T ], which proves the statements of the Proposition.
An immediate consequence of Proposition 10, together with estimate (19) is that Noticing that the definition of D h and Proposition 10 implies that the sequence of functions ρ h (t) strongly converges to ρ(t) := D n(t, z)dz in L 2 ((0, T )).

Existence of weak solution
This section is devoted to prove that the function n is a weak solution of problem (1-2).

Proposition 11
The function n satisfies , and for each j, multiply Eq. (16) by h 3 ϕ j (t) := h 3 ϕ(t, z j ), and add them up for all j, obtaining the relation Reordering the terms from A ϕ (t), we get where Thanks to the boundedness of u jl and the regularity of ϕ, this term satisfies Recalling the definition of u jl and the property (8), we get from (32) Defining V j (t) := V (t, z j ) allows us to write Thanks to the regularity of ϕ and V , we have Finally, we write Again, the regularity of ϕ and the boundedness of V imply And this together with (33), (34) and estimate (19) ensures that there exists a constant On the other hand, we can rewrite the remaining term as In Eymard et al. (2006), it was proven that Reordering as well the terms from D ϕ , we get Thanks to the boundedness of the coefficients A jl and the regularity of ϕ we get Recalling the definition of A jl , the fact that A(θ ) is a diagonal matrix and the normal vectors to the boundary of D j are elements of the canonical euclidean basis or their opposites, we have that, for all j and l A jl ∇ϕ(z jl ) · (z jl − z j ) = A(θ jl )∇ϕ(z jl ) · (z jl − z j ).
Consequently, we infer that . The regularity of ϕ, the boundedness The sequence of functions belongs to L ∞ (D T ) and strongly converges in Finally, we conclude that Putting this result together with (35), (37) and (38), we get that, for all ϕ ∈ As The estimate (30) implies that Q[n] ∈ L 2 ((0, T ); V ), consequently, n belongs to H 1 ((0, T ); V ).

Proposition 14 The function n is a weak solution of problem (1-2).
Proof We define ϕ ε (t) = ϕ * t ε for a mollifier ε with compact support included in (0, ∞) so that ϕ ε ∈ C 1 c ([0, T ); V) for any ε > 0 small enough and Writing the Eq. (31) for ϕ ε and passing to the limit ε → 0 we get that (31) also holds true for ϕ. 1, 0)), and such that the integral of χ is −1. For example, fix δ < 1/2 and define Passing to the limit when ε goes to 0 we obtain that n is a solution for the variational formulation.

A discrete implicit scheme
Once we have established the proof of existence of a solution for problem (16-17), together with the obtention of a semi-discrete scheme in order to approximate such solution, we proceed to derive an implicit discrete scheme starting from problem (16-17), and to prove its convergence. Consider a natural number K and define t = T K and t k = k t, ν k j := ν j (t k ), k = 1, . . . , K . Using a forward difference approximation for the time derivative in (16) we get the implicit scheme where M k+1 Theorem 15 Let ν 0 j be non-negative initial data with mass ρ 0 = N j=1 h 3 ν 0 j and assume that t < 1 then there exists a unique non-negative solution ν k j , k = 1, . . . , N to scheme (40). Furthermore, for each h, the sequence of piecewise constant functions strongly converges to the solution of (16-17) in (L 2 ((0, T ))) N when t goes to 0.
Proof For all t satisfying (41), there exists λ < 1 such that where ρ λ = ρ 1−λ . Consider the set and assume ν k ∈ R M to be the solution of (40) for a previous iteration, having all non-negative components and satisfying j h 3 ν n j ≤ ρ. Define the operator ν = F(η) : X → R M as the solution of the linear system where the components of matrix M(η) are defined as From the definition of X and the choice of t we have that P jl (η) is positive for all η if j = l, and non-positive if j = l. Furthermore P(η) is a column-dominant matrix, so that we may conclude that P(η) is an M-matrix, which implies that its inverse exists and has only non-negative entries. As ν k has all non-negative components, then ν = F(η) also has non-negative components for all η ∈ X . Multiplying system (43) by h 3 and adding up all equations, thanks to (42) we obtain This implies that N j=1 h 3 ν j < ρ λ and consequently, F(η) is a continuous application going from X to itself, and thanks to Brouwer's fixed point theorem, F(η) has at least a fixed point on X . Furthermore, a fixed point of F(η) will satisfy h 3 r j − d j ν 1 ν j which in turn implies ν 1 ≤ ρ. As a consequence, the implicit Euler scheme satisfies For the uniqueness, assume there exists two different solutions ν and μ to scheme (40). Let us denote the sign of β ∈ R as sign(β). Taking the difference for each equation, multiplying by h 3 sign(ν j − μ j ), adding up all the equations and recalling that strongly converges to M j (t, ρ h (t))ν * j + l∈N j B jl (t)ν * l , so, taking the limit in (45), we which implies that ν * j is in C 1 and is a solution for (16). Furthermore, as ν * j is the point-wise limit of μ j t , it also satisfies the initial conditions (17).

Simulations
The first part of this section is devoted to the numerical analysis of the approximation error. For certain values of the coefficients of problem (1-3), it is possible to obtain an analytical solution, which we will use in order to compare with our numerical approximation.
Assume that r (x, y, θ) and d (x, y, θ) are constants, that V (t, x, y, θ) is independent of t and that exists W (x, y, θ) such that Assume as well that Then, the solution for problem (1-3) is The existence of an analytic solution allows us to compare its values with those obtained from solving (16) for different values of h, and this way, numerically establish the error order of the method. Choosing For a grid of points {(t k , z j )} with t k = k t, k = 1, . . . , K , t > 0 and j = 1, . . . , N , we define the discrete L 2 (D T ) error for the semi-discrete scheme (16) as where ν(t) is the solution of the scheme for the functions introduced above. We set t = 0.01 and in Fig. 1 we show the dependence in log-log scale of E h = E(0.01, h) with respect to the inverse of the cell size M = 1/h.
In the same way, we define the discrete L 2 (D T ) error for the discrete scheme (40) as where ν k j is the solution of (40). We set h = 1/50 and plot the dependence of E 2 t := E 2 ( t, 0.02) with respect to t = 1/M 1 , again in log-log scale (Fig. 2).

Phenotypic dimorphism
We present now several examples illustrating the effect of the environment on a population and how the plasticity trait plays a role in surviving effects.

Monomorphic population
On the first place, we show the evolution of a population which is under the effects of natural selection and non-genetic epimuations, but without considering plasticity as a trait nor accounting for the environmental pressure. Specifically, we solve the problem (1-2) over the domain = {(x, y) ∈ [0, 1] 2 : (x − 1) 2 + (y − 1) 2 > 1}. This choice of domain represents the existence of a trade-off between traits and it is evidenced here by noticing that the individuals of the population which are close to the maximal value of one of the traits (x = 1 or y = 1), must forcibly be close as well to the minimal value of the other trait (y = 0 or x = 0). We take the growth rate as r (x, y) = e −(x−0.1) 2 −(y−0.1) 2 , the death rate as d(x, y) = 0.5, the diffusion parameters a 11 (θ ) = a 22 (θ ) = 10 −6 , and the drift terms V (t, x, y) = (0, 0). We take an initial condition (3) given by the expression . We choose the value of a in such a way that ρ 0 = n 0 (x, y)dxdy = 1. With this choice of n 0 we intend to represent a "strongly" monomorphic population. This is, a population where most of the individuals are concentrated around a single set of traits.
We observe in Fig. 3 that the dominant phenotype slowly converges to the point which maximises the fitness, in this case, the point (0.1, 0.1), which maximises the growth rate r (x, y).

Dimorphism due to the effect of the environment
For a second example we keep the same parameters, but add a drift term accounting for the effect of the environment (biologically, a "cellular stress"). Specifically we choose Notice that we still are not including plasticity as a trait in our analysis. Also notice that the function V (t, x, y) is not continuous, while the results presented in previous sections needed V (t, z) to be smooth in order to ensure the existence and uniqueness of solution for the problem, as well as the convergence of the finite volume method. This is not a big issue, because the conditions of our problem allow us to use a density argument in order to extend our results to any V (t, z) ∈ C([0, T ), L 2 (D)), and in any case all numerical approximations are smoothing approximations of the drift V . The choice of V can be seen (and experimentally replicated) as a certain type of "training": all the individuals of the population are "pushed" in the direction where they show their largest potential. We can appreciate in Fig. 4 that adding a drift term resulted in the appearance of dimorphism: the final population is concentrated around two different trait configura-  tions. This evolution into dimorphism due the action of environment, contrasts with the results obtained in Lorenzi and Pouchol (2020), where it was proved that dimorphism can occur in the absence of a drift term, given that the growth rate r (x, y) has several maximum point satisfying certain conditions.

Plasticity, environmental effect and dimorphism
We will now consider plasticity as a trait and modify the parameters from the previous examples accordingly. We first consider the growth rate as This choice of V is similar to the one shown before, only that now differentiation imposes a cost on adaptability: the more specialised you are, the harder it gets to adapt to new situations. Notice that a higher plasticity increases the effect of non-genetic epimutations (given by the diffusion term) and stress induced mutations (given by the drift term). For the initial data (3) we take a function of the form  We observe in Fig. 5 that for this set of parameters, the sub-population with the lowest plasticity quickly gets extinct (notice the scale of the density values), while the emergence of dimorphism can be appreciated for the one with higher plasticity.
This kind of behaviour persists up to the final stages of the evolution, when we can observe in Fig. 6 that the low plasticity sub-population is completely extinct while the high plasticity one has completed the differentiation process.

Another example of emergence of dimorphism
We present now a different example of emergence of dimorphism, but this time, not as a result of a response to the effect of the environment, but as a consequence of the existence of two maximum points for the growth rate. We will observe how a population initially concentrated around a single phenotype configuration, will evolve with time into a dimorphic population, in which each upcoming sub-population is more specialised and less plastic that in the initial configuration. For this purpose, over the domain D = × [0, 1] we consider an initial density given by the expression Fig. 6 Evolution of two sub-populations with different levels of plasticity: final stages with f (z) = z−z 0 2 (0.025) 2 , where z 0 = (0.25, 0.25, 0.5) and · is the euclidean norm. We choose the value of a in such a way that ρ 0 = D n 0 (z) = 1. We set the growth rate and the death rate as We choose the diffusion matrix and finally the drift term This time, the "push" towards specialisation imposed by V is inversely proportional to the current set of traits (individuals with traits (x, y) are specialising with a rate proportional to (−y, −x). The growth rate was chosen in such a way that it satisfies the sufficient conditions given in Lorenzi and Pouchol (2020) in order to guarantee the appearance of phenotypic polymorphism.
We show in Fig. 7 that initially the population is concentrated around the phenotype z 0 = (0.25, 0.25, 0.5), and gradually starts differentiating while loosing plasticity.
As the two new sub-populations become more and more differentiated, the loss in plasticity becomes more evident, and we see in Fig. 8 that most of the mass is migrating towards less plastic states, while the differentiation process continues.
Finally we observe in Fig. 9 that once the sub-populations are fully specialised, the concentration process continues and at the final stage t = 1000 we have a dimorphic population which is more specialised but less plastic that the initial one.
Overall, in the absence of a drift term, the effects of the growth and death rate, with small (or zero) diffusion has been widely studied, for example, see (Perthame 2006;Pouchol and Trélat 2018;Lorenzi and Pouchol 2020). In general, the points where the population concentrates are entirely determined by the reaction part, while the diffusion coefficients determine how concentrated the population is. The shape of the domain directly affects where the fitness function (which depends on the birth and death rates) attains its maxima, therefore, affecting as well where the concentration phenomena occurs.
The choice of initial state does not appear to affect the final configuration, but rather the dynamics of the population during the evolution in time. We have shown that if we start with an already concentrated population, we witness the continuous "migration" of the the population towards the fittest state.
Finally, if we consider the presence of the advection term, we observe that it may affect not only the amount of selected sets of traits, but their positions and the dynamics of convergence towards them as well.

Concluding remarks
The validity of the model we constructed is strengthened by the different evolutionary mechanisms described in Wagner et al. (2019), where the authors focus their attention on Stress-Induced Evolutionary Innovation, and compare it to plasticity-based models, in particular the Plasticity First Hypothesis. Quoting the authors: "SIEI and PFH are not competing models but explain different kinds of evolutionary processes that are some- Fig. 9 Final stages of the population density for different values of θ : Around t = 900 (bottom left) the differentiation process is over and most of the population has reached the plasticity level θ = 0.25. At t = 1000 (bottom right) we observe that the population concentrated around any other level of plasticity is almost extinct, and only the one around θ = 0.25 survives times distinct and sometimes combined over evolutionary time". Similar mechanisms were taken into consideration in the construction of our model, environmental stress (aka environmental pressure, biologically, at the single-cell level, "cellular stress") in the form of an advection term and and mutations thanks to plasticity in the form of a diffusion term, both accompanied by natural selection in the form of reaction term.
It is important to highlight the novelty that represents the inclusion of plasticity as a trait, which has not been considered before (with the exception of Bouin et al. 2012;Bouin and Calvez 2014 and the subsequent works on the cane toad spreading rate, where a similar structural variable is used, but to denote spatial diffusivity, and not adaptability potential, as in our case). Another novelty is the inclusion of constraints between traits modelled through a certain relation between the structural variables. Most of the previous work, and in particular (Chisholm et al. 2015), have considered the variable space as the entire unit square, or all of R d , d ∈ N, disregarding the existence of constraints between different traits, and the possible effects on the dynamics of the population this might have.
The finite volume method offers a powerful tool in order to numerically approximate the solution of integro-differential or reaction-diffusion equations, such as the one treated in the present paper. The preservation of the structure of the original problem at a semi-discrete level and the excellent approximation for the non-local terms are just two of the reasons why we chose this method. This way, we were able to obtain two numerical schemes in order to approximate the solution for an evolution problem modelling bet hedging strategies.
We proved the existence and uniqueness of solutions for such schemes, and constructed sequences of functions converging to the solution of the original problem. We approximated the convergence error by establishing a comparison with an exact solution. It is worth mentioning that the constructive character of the proofs may provide new and interesting tools in order to obtain further theoretical results.
After simulating various situations, we observed different ways in which a population can respond to external stress, depending on the plasticity levels of its individuals. A highly plastic sub-population can quickly adapt to its surrounding environment, guaranteeing this way its survival, while a less plastic sub-population might go extinct under the same external factor. Another strategy consists in "trading" some of the plasticity by a higher differentiation level.
Furthermore, the emergence of dimorphism as a consequence of external stress, not only is an interesting alternative to the previously established results from Lorenzi and Pouchol (2020), but also shows that bet-hedging strategies are a suitable response to (abrupt) external changes in the environment, and, at the same time, a possible way to survive them. It is fair mentioning that, throughout all the simulations, the symmetry hypothesis required in the reference (Lorenzi and Pouchol 2020) in order to observe dimorphism were respected. It remains to establish what are the essential conditions that will lead to the appearance of dimorphism when an advection term is present.
We thus provided here a rigorous model for the study of the emergence of dimorphism, an event that is likely to have been at the evolutionary origin of multicellularity by divergence of phenotypes and may thus provide a rationale for a renewed conception of animal evolution towards multicellular organisms, and, more pragmatically and consistently with the atavistic theory of cancer, for a possible origin of phenotype bet hedging in cancer cell populations.
Bet hedging in cancer cell populations is indeed a strategy susceptible to yield maximal probabilities of survival to a plastic cell population exposed to life-threatening insults such as by drugs or other deadly therapies. The modelling setting presented here may thus help in the future to test and optimise combined anticancer therapies involving chemotherapies, targeted therapies, and -what is likely still ahead of us for the present time -possible control of cell plasticity by epigenetic drugs.