A cell-cell repulsion model on a hyperbolic Keller-Segel equation
Xiaoming Fu, Quentin Griette, Pierre Magal
TL;DR
This paper introduces a hyperbolic Keller-Segel model for two-cell populations, proving solution existence and uniqueness, and explores how initial distributions influence short-term and long-term cell population dynamics through numerical simulations.
Contribution
The work develops a novel hyperbolic Keller-Segel model for cell populations, proving mathematical properties and analyzing the effects of initial conditions on population dynamics.
Findings
Model exhibits competitive exclusion differing from ODE models.
Initial total cell number influences population ratio.
Fast dispersion provides short-term advantage.
Abstract
In this work, we discuss a cell-cell repulsion population dynamic model based on a hyperbolic Keller-Segel equation with two populations. This model can well describe the cell growth and dispersion in the cell co-culture experiment in the work of Pasquier et al. \cite{Pasquier2011}. With the notion of solutions integrated along the characteristics, we prove the existence and uniqueness of the solution and the segregation property of the two species. From a numerical perspective, we can also observe that the model admits a competitive exclusion (the results are different from the corresponding ODE model). More importantly, our model shows the complexity of the short term (6 days) co-cultured cell distribution depending on the initial distribution of each species. Through numerical simulations, the impact of the initial distribution on the population ratio lies in the initial total cell…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 16
Figure 17
Figure 18
Figure 19
Figure 20
Figure 21
Figure 22
Figure 23
Figure 24
Figure 25
Figure 26
Figure 27
Figure 28
Figure 29
Figure 30
Figure 31
Figure 32
Figure 33
Figure 34
Figure 35
Figure 36
Figure 37
Figure 38
Figure 39
Figure 40| Symbol | Interpretation | Value | Unit | Method | Dimensionless value |
|---|---|---|---|---|---|
| time | 1 | - | 1 | ||
| inner radius of the dish | 2.62 | [31] | 1 | ||
| cell total number at | [31] | 0.01 | |||
| growth rate of cell | 0.6420 | fitted | 0.6420 | ||
| growth rate of cell | 0.6359 | fitted | 0.6359 | ||
| intraspecific competition of | fitted | ||||
| intraspecific competition of | fitted | ||||
| dispersion coefficient of | 13.73 | fitted | 2 | ||
| dispersion coefficient of | 13.73 | fitted | 2 | ||
| sensing coefficient | fitted |
| Parameters | Relations | ||||
|---|---|---|---|---|---|
| Case (ii) | 0.4 | 0 | 1 | 1 | |
| Case (iii) | 0.4 | 0 | 0.2 | 5 | |
| Case (iv) | 0.4 | 0 | 1 | 5 |
| Case (i) | Case (ii) | Case (iii) | Case (iv) | |
|---|---|---|---|---|
| Global attractor in ODE | Coexistence | Region dependent | ||
| Stable steady state in PDE |
| Parameters | ||||
|---|---|---|---|---|
| scenario 1: | 0 | 0 | ||
| scenario 2: | 0 | 0 |
| Parameters | ||||
|---|---|---|---|---|
| scenario 3: | 0.1 |
| Drug concentration () | 0 | 0.1 | 0.3 | 1 | 3 | 10 |
|---|---|---|---|---|---|---|
| Extra mortality (day-1) | 0 | 0.6619 | 0.8109 | 1.0118 | 1.5585 | 1.9545 |
| Extra mortality (day-1) | 0 | 0 | 0 | 0.0246 | 0.0569 | 0.2192 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsMathematical Biology Tumor Growth · Mathematical and Theoretical Epidemiology and Ecology Models · Gene Regulatory Network Analysis
A cell-cell repulsion model on a hyperbolic Keller-Segel equation
Xiaoming Fu , Quentin Griette and Pierre Magal
Univ. Bordeaux, IMB, UMR 5251, F-33400 Talence, France
*CNRS, IMB, UMR 5251, F-33400 Talence, France.
The research of this author is supported by China Scholarship Council.
Abstract
In this work, we discuss a cell-cell repulsion population dynamic model based on a hyperbolic Keller-Segel equation with two populations. This model can well describe the cell growth and dispersion in the cell co-culture experiment in the work of Pasquier et al. [31]. With the notion of solutions integrated along the characteristics, we prove the existence and uniqueness of the solution and the segregation property of the two species. From a numerical perspective, we can also observe that the model admits a competitive exclusion (the results are different from the corresponding ODE model). More importantly, our model shows the complexity of the short term (6 days) co-cultured cell distribution depending on the initial distribution of each species. Through numerical simulations, the impact of the initial distribution on the population ratio lies in the initial total cell number and our study shows that the population ratio is not impacted by the law of initial distribution. We also find that a fast dispersion rate gives a short-term advantage while the vital dynamic contributes to a long-term population advantage.
1 Introduction
In many recent biological experiments, the co-culture of multiple types of cells has been used for a better understanding of cell-cell interactions. This is a typical case in the context of studying cancer cells where the interaction between cancer cells and normal cells plays a crucial role in tumor development as well as in the resistance of cells to chemotherapeutic drugs. The goal of this work is to introduce a mathematical model taking care of the cell growth together with the spatial segregation property between two types of cells. Such a phenomenon was observed by Pasquier et al. [32]. They studied the protein transfer between two types of human breast cancer cell. Over a 7-day cell co-culture, the spatial competition was observed between these two types of cells and a clear boundary was formed between them on day 7 (see Figure 1). Segregation property in cell co-culture was also studied recently by Taylor et al. [37]. They compared the experimental results with an individual-based model. They found the heterotypic repulsion and homotypic cohesion can account for the cell segregation and the border formation. A similar segregation property is also found in the mosaic pattern between nections and cadherins in the experiments of Katsunuma et al. [21].
The early attempts to explain the segregation property by continuum equations date back to 1970s. Shigesada, Kawasaki and Teramoto [34] studied segregation with a nonlinear diffusion model and they found the spatial segregation acts to stabilize the coexistence of two similar species, relaxing the competition among different species. Lou and Ni [25] generalized the model and studied the steady state problem for the self/cross-diffusion model. For the nonlinear diffusion model, Bertsch et al. [5] in their work proved the existence of segregated solutions when the reaction term is of Lotka-Volterra type.
Here instead of using nonlinear diffusion models, we will focus on a (hyperbolic) Keller-Segel model. Such models have been used to describe the attraction and the repulsion of cell populations known as chemotaxis models. Theoretical and mathematical modeling of chemotaxis can date to the pioneering works of Patlak [30] in the 1950s and Keller and Segel [22] in the 1970s. It has become an important model in the description of tumor growth or embryonic development. We refer to the review papers of Horstmann [20] and Hillen and Painter [18] for a detailed introduction about the Keller–Segel model. To our best knowledge a model taking care of segregation property and cell-cell repulsion of Keller-Segel type has not been studied.
As we will explain in the paper, our model can also be regarded as a nonlocal advection model. Recently, implementing nonlocal advection models for the study of cell-cell adhesion and repulsion has attracted a lot of attention. As pointed out by many biologists, cell-cell interactions do not only exist in a local scope, but a long-range interaction should be taken into account to guide the mathematical modeling. Armstrong, Painter and Sherratt [2] in their early work purposed a model (APS model) under the principle of the local diffusion plus the nonlocal attraction driven by the adhesion forces to describe the phenomenon of cell mixing, full/partial engulfment and complete sorting in the cell sorting problem. Based on the APS model, Murakawa and Togashi [28] thought that the population pressure should come from the cell volume size instead of the linear diffusion. Therefore, they changed the linear diffusion term into a nonlinear diffusion in order to capture the sharp fronts and the segregation in cell co-culture. Carrillo et al. [8] recently purposed a new assumption on the adhesion velocity field and their model showed a good agreement in the experiments in the work of Katsunuma et al. [21]. The idea of the long-range attraction and short-range repulsion can also be seen in the work of Leverentz, Topaz and Bernoff [24]. They considered a nonlocal advection model to study the asymptotic behavior of the solution. By choosing a Morse-type kernel which follows the attractive-repulsive interactions, they found the solution can asymptotically spread, contract (blow-up), or reach a steady-state. Burger, Fetecau and Huang [6] considered a similar nonlocal adhesion model with nonlinear diffusion, they studied the well-posedness of the model and proved the existence of a compact supported, non-constant steady state. Dyson et al. [12] established the local existence of a classical solution for a nonlocal cell-cell adhesion model in spaces of uniformly continuous functions. For Turing and Turing-Hopf bifurcation due to the nonlocal effect, we refer to Ducrot et al. [10] and Song et al. [35]. We also refer the readers to Mogliner et al. [26], Eftimie et al. [13], Ducrot and Magal [11], Fu and Magal [16] for more topics about nonlocal advection equations. For the derivation of such models, readers can refer to the work of Bellomo et al. [4] and Morale, Capasso and Oelschläger [27].
In this article, we consider a two-dimensional bounded domain (a flat circular petri dish). We use the notion of solution integrated along the characteristics. Thanks to the appropriate boundary condition of the pressure equation (see Equation (2.2)), we deduce that the characteristics stay in the domain for any positive time. The positivity of solutions, the segregation property and a conservation law follow from the notion of solution as well. The main goal in this article is to investigate the complexity of the short-term (6 days) co-cultured cell distribution depending on the initial distribution of each species. Through the numerical simulations, we investigate the impact of the initial population number (as well as the law of initial distributions) on the population ratio. In the above mentioned literature, the numerical simulation are restricted to a rectangular domain with periodic boundary conditions. It is worth mentioning that here the domain is circular with no flux boundary condition for the pressure which requires a finite volume method (see Appendix 5.4).
The plan of the paper is the following. In Section 2, we present the model for the single-species case and we prove the local existence and uniqueness of solutions by considering the solution integrated along the characteristics and we prove a conservation law. Section 3 is devoted to the numerical analysis of the model. In Section 3.1, we consider the model homogeneous in space which corresponds to an ODE model that has been previously studied by Zeeman [39]. In Section 3.2, we investigated the competitive exclusion principle and the impact of the initial distribution on population ratio. The spatial competition due to the dispersion coefficients and cell kinetics is considered in Section 3.3. Section 4 is devoted to discussion and conclusion and all the mathematical details are presented in the Appendix.
2 Mathematical modeling
2.1 Single species model
Let us consider the following one-species model
[TABLE]
where satisfies the following elliptic equation
[TABLE]
We denote to be the unit open disk centered at with radius , i.e., . Here is the outward normal vector, is the dispersion coefficient, is the sensing coefficient. The divergence, gradient and Laplacian are taken with respect to . System (2.1)-(2.2) can be regarded as a hyperbolic Keller-Segel equation (with chemotactic repulsion) on a bounded domain.
Remark 2.1**.**
Equation (2.2) can be derived from the following parabolic equation (which is the classical case in the Keller-Segel equation [20]) as goes to [math]:
[TABLE]
The process of letting corresponds to the assumption that the dynamics of the chemorepellent is fast compared to the evolution of the cell density. In the case of chemoattractant a variant of such a model was considered by Perthame and Dalibard [33], Calvez and Dolak-Struß[7].
Remark 2.2**.**
As we mentioned in the introduction, Equation (2.2) can be regarded as a non-local integral equation by using Green’s representation
[TABLE]
2.1.1 The invariance of domain and the well-posedness of the model
We remark that in System (2.1)-(2.2) we do not impose any boundary condition directly on . Instead, the boundary condition here is induced by . If we consider the associated characteristics flow of (2.1)-(2.2)
[TABLE]
We can prove (see Appendix 5.1) the characteristics can not leave the domain (see Figure 2 for an illustration). In fact, we can prove for any , the mapping is a bijection from to itself (see Lemma 2.10). We consider the solution along the characteristics
[TABLE]
Taking any , there exists such that , and since
[TABLE]
we can reconstruct the solution from and on .
Assumption 2.3**.**
Assume the vector field is continuous in and lipschitzian with respect to in .
Remark 2.4**.**
Assumption 2.3 is a sufficient condition for the existence and uniqueness of the characteristic flow in (2.4).
Definition 2.5**.**
For any bounded open domain . If is bounded and continuous, we write
[TABLE]
For any , the –Hölder norm of is
[TABLE]
where
[TABLE]
The Hölder space consists of all functions for which the norm
[TABLE]
is finite.
Lemma 2.6**.**
[17, Theorem 6.30 and 6.31]** Let to be a unit open disk. Consider the following elliptic equation
[TABLE]
where is the outward unit normal vector on . Then for all , the elliptic problem (2.5) has a unique solution . Moreover,
[TABLE]
where .
The following theorem tells us if we choose our initial value to be smooth enough, then Assumption 2.3 can be satisfied and the existence and the uniqueness of solutions follow.
Theorem 2.7** (Existence and uniqueness of the solution along the characteristic).**
Let . Then for some there exists a unique non-negative solution to (2.1)-(2.2) which satisfies . Moreover for any , we have and .
The proof the above theorem will be detailed in Appendix 5.2.
Remark 2.8**.**
Since for any and for any , we have , we deduce from Lemma 2.6 that . Therefore, is continuous (since ) and lipchitzian with respect to which implies Assumption 2.3.
2.1.2 Conservation law on a volume
If the reaction term in System (2.1)-(2.2), the boundary condition implies the conservation law for . This can be seen through the solution along the characteristics. In fact, we have the following conservation law.
Theorem 2.9**.**
For each volume and each we have
[TABLE]
In particular, if we have , then for any
[TABLE]
This means the total number of cell in the volume is constant along the volumes .
Before proving Theorem 2.9, we need the following lemma.
Lemma 2.10**.**
Let and to be the characteristic flow generated by (2.4). Then the map is continuously differentiable and one has the determinant of Jacobi matrix:
[TABLE]
where is the Jacobian matrix of with respect to at point .
Proof.
From Theorem 2.7 and Remark 2.8, the mapping is and for any if . This ensures the characteristics is continuously differentiable. Taking the partial derivative of Equation (2.4) with respect to yields
[TABLE]
where is the Jacobian matrix of with respect to at point . For any matrix-valued function , the Jacobian formula reads as follows
[TABLE]
Hence, we obtain
[TABLE]
and since . Therefore, we have
[TABLE]
Therefore the result follows. ∎
Proof of Theorem 2.9.
Let to be the characteristic flow generated by (2.4). Given any measurable set and any , we integrate over the volume with respect to
[TABLE]
where we changed the variable to on the right-hand-side.
For the right-hand-side, we will prove in Appendix 5.2 that
[TABLE]
Combining with (2.6) we obtain that
[TABLE]
Substituting the above equation into (2.7) gives us
[TABLE]
which is equivalent to
[TABLE]
The result follows. ∎
2.2 Multi-species model
2.2.1 Multi-species ODE model
Let us consider the corresponding two species model without the spatial variable that is for
[TABLE]
We adopt the Lotka-Volterra model by letting
[TABLE]
where are the growth rate, represent the mutual competition between the species, is the competition among the same species and is the additional mortality rate caused by drug treatment. In Section 2.2.1 we will always assume for (when , one can regard as a whole). If we consider (2.8) in the absence of the other species, we can rewrite (2.9) as
[TABLE]
We always assume that for each , meaning that each species alone exhibits logistic growth. This model has been considered by many authors (for example, see [29, 39]). Here we give a short summary of some qualitative properties of the solution to (2.8) in order to compare with the PDE model.
Equilibrium and stability for (2.8)-(2.9)
One can easily compute the system has the following equilibrium
[TABLE]
where
[TABLE]
The solution is only of relevance when and is strictly positive, which is equivalent to say
[TABLE]
A scheme of the qualitative behavior of the phase trajectory is given in Figure 3 by numerical simulations.
We adapt the main stability results from Zeeman [39] where the author considered a general –species extinction case, Murray [29, Chapter 3.5] and Hirsch [19, Chapter 11] to system (2.8)-(2.9) for the following fours cases (i)-(iv) and discuss their biological implications.
Proposition 2.11**.**
For system (2.8)-(2.9), suppose for each , and for any . Let be the equilibrium for each species alone and assume the initial value lies strictly in the first quadrant that is and . Then for the following four cases we have
- (i).
* This case corresponds to Figure 3 (a). The system (2.8) has four positive equilibrium, namely and . In such case, only is global globally asymptotic stable in the region .*
- (ii).
* This case corresponds to Figure 3 (b). The system (2.8) has three positive equilibrium, namely and . Only is globally stable in the positive quadrant excepted for the axis .*
- (iii).
* This case corresponds to Figure 3 (c). The analysis of the stability is similar to the case (ii). Only is globally stable in the positive quadrant excepted for the axis .*
- (iv).
* This case corresponds to Figure 3 (d). In this case, system (2.8) has four equilibrium, where and are stable while is a saddle point. The steady states and have two non-overlapping domains of attraction, separated by the stable manifold S of equilibria .*
Remark 2.12**.**
Although among the four cases, (ii) and (iii) always lead to the principle of exclusion and so do (iv) due to the natural perturbation in population levels, we still have the case (i) where the two species can coexist in the long term. As we further develop our PDE model for (2.8), we can show numerically that the principle of exclusion dominates even when case (i) is satisfied and this is a evident difference compared to ODE model (2.8).
2.2.2 Multi-species PDE model
We study a two species population dynamic model on a unit open disk given as follows
[TABLE]
where is the outward normal vector, is the dispersion coefficient, is the sensing coefficient. Recall the function is of form
[TABLE]
System (2.11) is supplemented with initial distribution
[TABLE]
2.2.3 Segregation property
From the mono-layer cell populations co-culture experiments, we can see that once the two cell populations confront each other, they will stop growing, thus, forming the separated islets. We can prove that our model (2.11) preserves such segregation property.
Theorem 2.13**.**
Suppose is the solution of (2.11)-(2.12) and assume in (2.11). Then for any initial distribution with for all , we have for any and .
Proof.
We argue by contradiction, assume there exist such that
[TABLE]
Suppose the characteristic flow satisfies the following equation
[TABLE]
Since is invertible from to itself, there exists some such that . Then for any we have
[TABLE]
which implies
[TABLE]
This is a contradiction. ∎
For the one dimensional case , suppose are solutions to (2.11)-(2.12), we give an illustration (see Figure 4) of the segregation for the solutions integrated along the characteristics for . In fact, if there exists for some such that for Then from Equation (2.13) we obtain
[TABLE]
Therefore, the characteristics forms a segregation barrier for the two cell populations.
2.2.4 Conservation law on a volume
If we assume that in system (2.11), we have the following similar conservation law for two species case. Suppose volume and each :
[TABLE]
Therefore, if we have for any
[TABLE]
This means the total cell number of the species is constant along the characteristics starting from the volume .
3 Numerical simulations
3.1 Impact of the segregation on the competitive exclusion
We set to be the total number at time
[TABLE]
We give the parameter values used in the simulations and their interpretations in Table 1. The parameter fitting for the growth rate and the intraspecific coefficients are detailed in Appendix 5.3.
The goal of our simulations is to compare the various cases discussed in Proposition 2.11 (ODE case) with our PDE model with segregation. As we will see in the numerical simulations, the model with spatial structure can present totally different results compared to the previous ODE model. To that aim, we firstly consider the case where the drug (doxorubicine) concentration is low in the cell co-culture for MCF-7 and MCF-7/Doxo. The drug treatment causes an additional mortality to the sensitive population MCF-7 represented by while no extra mortality to the resistant population MCF-7/Doxo represented by (MCF-7/Doxo is resistant to a small quantity of drug treatment see Table 6 in Appendix 5.3).
Now since we consider the presence of the drug, the equilibrium (2.10) in the ODE case should be rewritten as
[TABLE]
Moreover we assume the drug concentration is low such that and , therefore we have
[TABLE]
The case when is similar and will be discussed in the end of this section.
Case (i): By using (3.2), the condition in Case (i) can be interpreted by
[TABLE]
Since we have and , if the coefficients and are small, then Case (i) holds. We give a possible set of parameters satisfying Case (i) :
[TABLE]
We assume that for each species , the initial distribution follows the uniform distribution on a disk with initial cell clusters (represented by the red/green dots in Figure 5 (a)). The initial total cell number is in (3.1) for each species and we assume each cluster contains the same quantity of cells. We present its numerical simulation in Figure 5 from day [math] to day . We also plot the relative cell numbers in Figure 5 (f) where we define the relative cell number for species i as
[TABLE]
In Case (i) of the ODE system (2.8), Proposition 2.11 shows that the two species can coexist with the equilibrium
[TABLE]
However, as shown in Figure 5, we can see the population density tends to 0 and tends to 1. Next, we consider the Cases (ii)-(iv) in Proposition 2.11 by choosing the parameters in each case as follows.
With the simulations in Figure 5, Figure 6 and the results in Proposition 2.11, by setting
[TABLE]
we can compare the stability between ODE model (2.8) and PDE model (2.11) under four different cases.
The numerical simulation strongly indicates that the stable steady states only depend on the relation between and . If (resp. ), the population (resp. ) will dominate and the other species will die out. We also did the four cases when , the results showed that is the only stable steady state, which verifies our conjecture. Since the results are similar we omit the numerical simulations.
One can notice that unlike the ODE system (2.8), the segregation property for the PDE model implies that it is impossible for the two species to coexist at a same position . Moreover, through the numerical simulations we observed that the PDE model (2.11) always undergoes a competitive exclusion principle, unless the equilibrium in (3.2).
3.2 Impact of the initial distribution on the population ratio
In the previous section, we considered the competitive exclusion principle for the two species. By studying the relative proportions of and , we presented the relation of the interspecific competition in our numerical simulation. Moreover, we can discover in Figure 5 (f) and in Figure 6 (a)-(c) that the increase of the proportion of the dominant population (red curve) is varying with time. It is evident to see from day 0 to day 2 the increase of the dominant population is faster than the increase from day 4 to day 6. If we further study the spatial-temporal evolution of the cell co-culture presented in Figure 5 (a)-(e), we can observe that from day 0 to day 2 the competition between the two groups is mainly expressed in the competition for space resources. However, from day 4 to day 6, when the surface of the dish is almost fully occupied by cells, the reaction term in the equation begins to play a major role influencing the change in the number of cells. In order to explore the major factors in cell competition, we consider the impact on the initial distribution. We will mainly focus on two factors, namely the initial cell total number and the law of initial distribution, which might influence the proportions for and . To that aim, we set the following parameters
[TABLE]
and the other parameters are given in Table 1.
3.2.1 Dependency on the initial total cell number
In cell culture, the initial number of cell cluster is an important factor. Bailey et al. [3] study the sphere-forming efficiency of MCF-7 human breast cancer cell by comparing the cell culture with different initial numbers of the cell cluster. Here we consider the impact of the initial cluster number on the final proportion of species. To that aim, we assume the initial distribution follows the uniform distribution on a disk. We consider two sets of initial condition, that is
[TABLE]
where and are defined in (3.1) and (respectively ) is the initial number of cell clusters of species (respectively species ).
The above initial conditions correspond to different types of seeding in the experiment, namely cells are sparsely seeded or densely seeded. We assume the total cell number is proportional to the initial number of cell cluster, meaning the dilution procedure adopted in the experiment is the same, thus the number of cells in each cell cluster is a constant.
In Figure 7, we first give a numerical simulation for the cell growth with sparse seeding.
In Figure 8, we present the simulation under the dense seeding condition, tracking from day 0 to day 6.
In Figure 9 we plot the evolution of the total number and its proportion for species and over 6 days of the simulation.
From Figure 9 (a)-(b), since is resistant to the drug, the number of population is much greater than . However, we can also observe a difference in their cell growth curves. In Figure (a) we can see that both cells are in the period of exponential growth from day 0 to day 6 (a base-10 log scale is used for the y-axis). Conversely, in Figure (b) the growth curves for both cells are converging to a constant from day 4 to day 6, implying that the cell co-culture is reaching a saturation stage. More importantly, in Figure (c), we observe a significant difference in the development of population ratios. In fact, since the spatial competition is still the dominant factor in the first two days, we can hardly see any difference between the dashed lines and solid lines. The proportion of the dominant population grows almost linearly. However, the proportion of the densely seeded group changed much slower after day 4, while the sparsely seeded group still grows linearly. This shows that although the growth rate of is at a competitive disadvantage, due to the sufficient number of cluster in the initial stage, and due to the segregation principle, does not die out in a short time in the competition. Although the competitive exclusion applies in this case, the time for the extinction of will be very long.
3.2.2 Dependency on the law of the initial distribution
In the experiment, the size of the cell dish can be a factor to determine the law of the initial distribution for the cell. In general, under the same total cell number, a small size cell dish will lead to a biased initial distribution and cells are more likely to aggregate at the border. While a big size cell dish will make the cell distribution more homogeneous, thus the initial distribution follows a uniform distribution. Therefore, in this section, we study whether the population ratio can be affected by the law of initial distribution. We will choose the beta distribution for the choice of the radius and the angle follows the uniform distribution on , that is
[TABLE]
The coordinate transformation of the initial distribution to a unit disk is as follows
[TABLE]
In Figure 10, we plot the density function of the beta distribution for different
[TABLE]
where is a normalization constant to ensure that the total integral is .
Our simulation will mainly compare the following two cases
[TABLE]
We plot the initial distributions of the two different cases in Figure 11 where we choose cell clusters (i.e., and in (3.6)) for species and species .
Suppose the initial cell clusters and cell total number , which is equally distributed in each cell cluster. Typical numerical solutions are shown in Figure 12 when and in Figure 13 when .
Now we plot the evolution of the total number for species and over 6 days.
From Figure 14 we can see that the law of initial distribution has almost no influence on the final proportion of species. We also tried different scenarios when the cell total numbers are 20, 50 and 100 or with different extra mortality rate and , the results are similar. Thus we can deduce that the final relative proportion is stable under the variation of the law of the initial distribution.
Combining the above numerical experiments in Section 3.2.1 and Section 3.2.2, we can see that under the competitive principle, the difference in the spatial resources can change the competition induced by the cell dynamics. To be more precise, under the case of sufficient spatial resources, the competitive mechanism can be more sufficiently expressed than in the case of less spatial resources. In Section 3.2.2, although we changed the law of the initial distribution of the cell seeding, as for the overall spatial resources, it is the same for both species. Therefore, the result of the competitive principle is almost the same in terms of the total number and the population ratio of the two populations.
3.3 Impact of the dispersion coefficient on the population ratio
In Section 3.2, when the parameters of the model are the same, the competition induced by the cell dynamics can be reflected by the difference in the spatial resource. Now we assume the spatial resource is the same and we investigate the role of the dispersion coefficient in the evolution of the species.
To that aim, we let the initial distribution of the two species follow the same uniform distribution and they are sparsely seeded on the dish. Furthermore, we let the cell dynamics for the two population to be almost the same, the only variable we control here is the dispersion coefficient for the population. We take the same uniform initial distribution at day 0, with the same number of initial cluster and the same amount of cell total number, i.e.,
[TABLE]
We compare the following two scenarios in Table 4 where the only difference is the dispersion parameters.
In scenario 1, the dispersion coefficients of the two species are the same, while in scenario 2 we suppose the species has an advantage in the spatial competition over its competitor .
Now we plot the evolution of the total number and the population ratios for species and over 6 days.
The main result from Figure 16 is that the dispersion coefficient can have a great impact on the population ratio after 6 days.
Next, we consider the following scenario where has the advantage in dispersion coefficient but is at a disadvantage induced by drug treatment. Therefore
By including now a drug treatment, we can see from Figure 17 and Figure 18 that between day [math] and day , the population dominates over thanks to a larger dispersion rate. After day , since the drug is killing the cell from species while the drug has no effect on the species , the species finally takes over the species . It leads to a gradual increase in its proportion of the population ratio.
In the numerical simulations for the scenarios and in Table 4, we let the cell dynamics of the two species be almost equal. Thus the competition due to the cell dynamics is almost negligible. We have shown the dispersion coefficient of populations can have a great impact on the population ratio after 6 days.
In the simulation for scenario 3 in Table 5, we can observe that despite the competitive exclusion, a larger dispersion coefficient can lead to a short-term advantage in the population. In the long term, the competitive exclusion principle still dominates.
4 Conclusion and Discussion
From the experimental data in the work of Pasquier et al. [31], we modeled the mono-layer cell co-culture by a hyperbolic Keller-Segel equation (2.11). We proved the local existence and uniqueness of solutions by using the notion of the solution integrated along the characteristics in Theorem 2.7 and proved the conservation law in Theorem 2.9. For the asymptotic behavior, we analyzed the problem numerically in Section 3.
In Section 3.1 we discussed the competitive exclusion principle, indicating that the asymptotic behavior of the population depends only on the relationship between the steady states and (see (3.2) for definition) which is different from the ODE case. We found that except for the case , the model with spatial segregation always exhibits an exclusion principle.
Even though the long term dynamics of cell density is decided by the relative values of the equilibrium, the short term behavior need a more delicate description. We studied two factors which may influence the population ratios. The first factor is the initial cell distribution, as measured by the initial total cell number and the law of initial distribution. We found that the impact of the initial distribution on the population ratio lies in the initial total cell number but not in the law of initial distribution.
The second factor influencing the population ratio is the cell movement in space, as measured by the dispersion coefficient . We can conclude that in the case of sparsely seeded initial distribution by the dispersion rate and cell dynamics. In the transient stage (i.e. before the dish is saturated), the dispersion rate are the dominant factor. Once the surface of the dish is saturated by cells, cell dynamics becomes the key factor. Note that the coefficients do not play any role in the competition because of the segregation principle.
We can briefly summarize the following main factors that can influence the population ratio in cell culture for model (2.11):
- (a).
The difference of cell dynamics in the two species (internal factor): if the equilibrium (see (3.2) for definition), then will dominate, will die out (and vice-versa when ) (see Figures 5-6 and Table 3); 2. (b).
If cells are densely seeded at the beginning, despite of competitive advantage, the dominant species can not take out its competitor in a short time (see Figures 7-8). We also concluded that the law of initial distribution has almost no influence on the population ratio (see Figures 12-13); 3. (c).
If cells are sparsely seeded at the beginning, we need to distinguish the period of time needed for the cell to occupy the surface of the dish and the time needed for each species to reach a saturation stage (see Figure 16 and Figure 18).
5 Appendix
5.1 Invariance of domain
In this section, we prove the invariance of domain for the characteristic equation.
Assumption 5.1**.**
Let be an open bounded subset with of class .
Since is a bounded domain of class , there exists a neighborhood of the boundary such that the distance function restricted to has the regularity (see Foote [15, Theorem 1]). Furthermore, by Foote [15, Theorem 1] and Ambrosio [1, Theorem 1 p.11], we have the following properties for .
Lemma 5.2**.**
Let Assumption 5.1 be satisfied. Then
- (i).
There exists a small neighborhood of with such that, for every there is a unique projection satisfying . 2. (ii).
The distance function is on . 3. (iii).
For any , where is the outward normal vector.
We consider the following non-autonomous differential equation on
[TABLE]
Assumption 5.3**.**
The vector field is continuous and satisfies
[TABLE]
Moreover, for any , there exists a constant such that vector field satisfies
[TABLE]
By (5.3), we have the existence and uniqueness of the solutions of (5.1) and the solutions may eventually reach the boundary in finite time. We will prove that (5.2) implies that the solutions of (5.1) actually stay in and can not attain boundary in finite time under Assumption 5.1.
Theorem 5.4**.**
Let Assumption 5.1 and 5.3 be satisfied. For any , let be the solution of (5.1) on . Then for any .
Proof.
We prove this theorem by contradiction. Let be the first time when reaches boundary , i.e.,
[TABLE]
We can find a such that, for any . Since is , the mapping is on . By Lemma 5.2 (iii), we have
[TABLE]
where is the outward normal vector and is the unique projection of onto . By assumption (5.2), we have
[TABLE]
Hence (5.4) becomes
[TABLE]
which yields
[TABLE]
and implies which contradicts our assumption . ∎
5.2 Proof of Theorem 2.7
Solution integrated along the characteristics. Let us temporarily suppose , we can rewrite the first equation in (2.1) as
[TABLE]
Moreover, if we differentiate the solution along the characteristic with respect to then
[TABLE]
The solution along the characteristics can be written as
[TABLE]
For the simplicity of notation, we let in our following discussion and define . We construct the following Banach fixed point problem for the pair . For each , we let
[TABLE]
where we set for any and we define
[TABLE]
where is the resolvent of the Laplacian operator with Neumann boundary condition.
We define
[TABLE]
where are two constants to be fixed later. We also set
[TABLE]
Notice is a complete metric space for the distance induced by the norm . For simplicity, we denote and for .
Theorem 5.5** (Existence and uniqueness of solutions).**
For any initial value and , for any , large enough in (5.8), there exists such that the mapping has a unique fixed point in .
Proof.
For any positive initial value and , we fix to be a constant such that and is a constant defined in (5.19) later in the proof.
We also denote
[TABLE]
and let be the closed ball centered at with radius in with usual product norm
[TABLE]
and we set
[TABLE]
Suppose , we need to prove that there exits a small enough such that the following properties hold
- (a).
For any , in (5.6) and (5.7) belong to and their norms satisfy
[TABLE] 2. (b).
Moreover, we have
[TABLE]
Moreover, we plan to show that the mapping is a contraction: there exists a such that for any \begin{pmatrix}\tilde{w}\\ \tilde{P}\end{pmatrix},\begin{pmatrix}w\\ P\end{pmatrix}\in\overline{B_{\tilde{Z}^{\tau}}}\bigg{(}\begin{pmatrix}w^{0}\\ P^{0}\end{pmatrix},r\bigg{)} we have
[TABLE]
Step 1. We show that there exists a small enough such that for any then
[TABLE]
where is defined in (5.6).
Indeed, since is Lipschitz continuous, then is also Lipschitz continuous. Since maps into , we have
[TABLE]
For any , we let . By the definition of in (5.6), we have
[TABLE]
Next we estimate . We have for any
[TABLE]
Since is the unit open disk, . We can obtain the following estimate
[TABLE]
Thanks to Grönwall’s inequality, we have
[TABLE]
Substituting the (5.15) into (5.14) yields
[TABLE]
Since , we can choose and we obtain
[TABLE]
Thus, Equation (5.9) holds.
Let us now check that satisfies (5.11). Let , we remark that for all . We have
[TABLE]
where . From (5.17) we have
[TABLE]
Since , it suffice to take small enough to ensure (5.11).
Step 2. Next we verify (5.10) and (5.12) for where is defined as the second component of (5.7). We show that there exists small enough such that for any
[TABLE]
Thanks to the Schauder estimate [17, Theorem 6.30], there exists a constant depending only on such that
[TABLE]
Recalling as a consequence of (5.15), we have
[TABLE]
We can now define
[TABLE]
which only depends on and . Finally, we let and we have
[TABLE]
In particular, we have shown (5.10).
Next we prove (5.12). Since is a two-dimensional unit disk, using Morrey’s inequality [14, Chapter 5. Theorem 6], we have
[TABLE]
where is a constant depending only on . For the sake of simplicity, we use the same notation for a universal constant depending only on in the following estimates. Moreover, by the classical elliptic estimates we have
[TABLE]
This implies that
[TABLE]
where we have used (5.18) for the last inequality . We can conclude
[TABLE]
Thus, it suffice to take small enough to ensure the neighborhood condition (5.12).
Step 3. Contraction mapping In order to verify (5.13), we let . We observe that
[TABLE]
Due to the classical inequality which holds for any , we deduce
[TABLE]
To estimate in (5.20), we claim that
[TABLE]
Indeed, we can obtain that
[TABLE]
This leads to
[TABLE]
Again due to Grönwall’s inequality, we conclude that (5.21) holds.
Inserting (5.21) into (5.20) we have
[TABLE]
where we set
[TABLE]
and as .
Next we prove the contraction property for . As before, applying the same argument of Morrey’s inequality and the classical elliptic estimates, we can deduce
[TABLE]
where we used (5.21) in the last inequality and is a constant depending only on . Defining and together with (5.22) we obtain
[TABLE]
Combing with (5.22) and (5.23) we deduce
[TABLE]
where as . If is small enough, this implies (5.13) for some . Since is complete metric space for the distance induced by the norm in , the result follows by the classical Banach fixed point theorem. ∎
5.3 Parameter fitting
From the work in [31], MCF-7 and MCF-7/Doxo cells are cultured at initial cell number separately in mm cell dish with or without doxorubicine. We use the cell proliferation data followed every 12 hours during six days to fit the parameters of the following ordinary differential equation
[TABLE]
Here we use to represent the MCF-7 (sensitive to drug) and to represent the MCF-7/Doxo (resistant to drug) and is the growth rate is the extra mortality rate caused by drug (doxorubicine) treatment and is a coefficient which controls the number of saturation.
In the work [36] cell proliferation kinetics for MCF-7 is studied over 11 days in 150 flask. Following an inoculation of cells at day [math], a maximum cell density of to cells/flask was reached at day . Therefore, we assume the saturation number for each species in mm (surface of ) dish satisfies
[TABLE]
By fixing the saturation number, we first estimate the growth rate of each species under zero drug concentration, namely . We divide the cell number by (the initial cell number) and rescale the parameters as follows
[TABLE]
As seen in Figure 19, without treatment, MCF-7 and MCF-7/Doxo displayed very similar growth rates, 0.6420 and 0.6359 per day, respectively.
By fixing the parameters
[TABLE]
we consider different scenarios with the drug concentration varies from to (see Figure 20) and we estimate the extra mortality rate for each population due to doxorubicine (see Table 6).
5.4 Numerical Scheme
For simplicity, we give the numerical scheme for the following one species and one dimensional model
[TABLE]
The numerical method used is based on finite volume method. We refer to [23, 38] for more results about this subject. Our numerical scheme reads as follows
[TABLE]
with the flux defined as
[TABLE]
and
[TABLE]
where we define
[TABLE]
where is a constant and is the usual linear diffusion matrix with Neumann boundary condition. Therefore, since the Neumann boundary condition corresponds to a no flux boundary condition, we impose
[TABLE]
which corresponds to and .
The numerical scheme at the boundary becomes
[TABLE]
By this boundary condition, we have the conservation of mass for Equation (5.28) when the reaction term .
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] L. Ambrosio, Geometric evolution problems, distance function and viscosity solutions, Calc. Var. Partial Differential Equations Springer, Berlin, Heidelberg, (2000), 5-93.
- 2[2] N. J. Armstrong, K. J., Painter and J. A. Sherratt, A continuum approach to modelling cell-cell adhesion, J. Theoret. Biol., 243 (2006), 98-113.
- 3[3] P. C. Bailey, R. M. Lee, M. I. Vitolo, S. J. Pratt, E. Ory, K. Chakrabarti, C. J. Lee, K. N. Thompson and S. S. Martin, Single-Cell Tracking of Breast Cancer Cells Enables Prediction of Sphere Formation from Early Cell Divisions, i Science , 8 , (2018) 29-39.
- 4[4] N. Bellomo, A. Bellouquid, J. Nieto and J. Soler, On the asymptotic theory from microscopic to macroscopic growing tissue models: An overview with perspectives, Math. Models Methods Appl. Sci. , 22(01) (2012), 1130001.
- 5[5] M. Bertsch, D. Hilhorst, H. Izuhara and M. Mimura, A nonlinear parabolichyperbolic system for contact inhibition of cell-growth, Diff. Equ. Appl. 4 (2012), 137-157.
- 6[6] M. Burger, R. Fetecau and Y. Huang, Stationary states and asymptotic behavior of aggregation models with nonlinear local repulsion, SIAM J. Appl. Dyn. Syst. , 13(1) (2014), 397-424.
- 7[7] V. Calvez and Y. Dolak-Struß, Asymptotic behavior of a two-dimensional Keller-Segel model with and without density control, In Mathematical Modeling of Biological Systems, Volume II (2008) (pp. 323-337). Birkhäuser Boston.
- 8[8] J. A. Carrillo, H. Murakawa, M. Sato, H. Togashi and O. Trush, A population dynamics model of cell-cell adhesion incorporating population pressure and density saturation, J. Theor. Biol. ， 474 (2019) 14-24.