Relaxation of the Cahn-Hilliard equation for Biology

Get Complete Project Material File(s) Now! »

Mathematical representation of living tissues

First steps of mathematical modelling

We present the formulation of a prototypal model of N N∗ components from which we can derive models that include precise effects. We define the domain Ω(t) Rd as the tissue or a part of it, and where d = 1, 2, 3 is the dimension and t is time. The boundary of Ω(t) is denoted by ∂Ω(t) and is assumed to be sufficiently regular. In the following, we adopt the Eulerian point of view: for a fixed point in space, we observe what happens without moving x. The second point of view (that we will not consider in this manuscript) is the Lagrangian or material point of view, where the observer follows the flow i.e. the material point X Ω(t) of the domain moves with respect to the flow given by a certain velocity. To describe how the different forces act on the tissue, we take an arbitrary volume V (t) Ω(t). We define Br(x) as the ball of center x and of radius r 0. We define the mass density of the i-th component at (t, x) R+ Ω in Eulerian where Ji(t, x) = Ji(ρ, ρ, v) is the net flow of the components across the boundary of the arbitrary volume V (t). This flux is defined as a function of the density of the i-th component, of its gradient and of a velocity field vi(t, x) which needs to be defined using a constitutive relation or given by an equation. However, to be able to represent the Physics of the tissue, one needs to use laws of mechanics and thermodynamics to derive a set of physically relevant equations and explain the different constitutive relations. In (1.1), the term Si(t, x) represents the growth of the tissue but also its degradation. This proliferation effect depends on many aspects occurring both at the scale of the cells and the tissue, such as the stress in the tissue, the amount of available nutrients, or the available space in the neighborhood of dividing cells.
The integral form of the continuity equation is useful to give a description of the different forces acting on the arbitrary volume V (t). However, since the equation (1.1) is satisfied for any volume V (t), we use the local form of this equation ∂ρi(t, x) + ∇ · J (t, x) = S (t, x), (1.2)
In the following, we present the derivation of two essential models of living tissues—the Cahn- Hilliard model for Biology and the Keller-Segel model. The first one is used to represent tissues as a multiphase fluid and describes the attractive and repulsive interactions of cells. As a result, the model can reproduce the formation of patterns and is used to represent tumors as fluids. The second is modeling chemotaxis: a type of movement that cells exhibit in nature. Indeed, cells have the capacity to sense their micro-environment, and their migration is driven by signals. The Keller-Segel model represents cells’ movement toward zones of a large concentration of certain chemicals called chemoattractants and due to random Brownian motion.

The Cahn-Hilliard model for modelling of tissues and tumours

Derivation. A detailed description of the derivation of the Cahn-Hilliard equation (CH in short) from the different mechanics and thermodynamics laws requires lengthy calculations. For the sake of simplicity, we give here a simple description of the Cahn-Hilliard equation for the mixture of two incompressible fluids. We also assume that the tissue under investigation is not growing nor degrading. This model is relatively simple and focuses only on the organization of two types of cells. However, as we will see in the following of this manuscript, it already gives a good description of tissues (especially tumors) and is at the center of many research pieces.
To satisfy the previous assumptions, we set N = 2 and Si(t, x) = 0 (for i = 1, 2). For the representation of tissues, we can assume that cells inside the phase i = 1 are cells constituting the tissue of interest (or the tumor), and the phase i = 2 is used to represent the rest of the cells of the micro-environment. From the incompressibility assumption, the continuity equation (1.2) can be simplified and written to follow the evolution of the mass fraction ci or the volume fraction ni. We present the case where we focus on volume fractions, and since these quantities satisfy n1 + n2 = 1, where we have used the fact that J in (1.2) was equal to J = ρihi along with the fact that the fluids are incompressible, i.e. div (ρi) = 0 for all i = 1, 2. Before giving the constitutive relation for the flux h, we first describe the energy associated to this model and use some of basic thermodynamics quantities.
Indeed, the formulation of the free energy associated to the system is often considered the starting point to give a simple derivation of the Cahn-Hilliard model. As in its original descrip- tion, the Ginzburg-Landau free energy is given by E[n](t) := ∫ γ |∇n|2 + ψ(n) dx.
The free energy density is the sum of two important terms. The surface tension γ n 2 is a force occurring at the interface between the two phases. This terms has the effect to penalize large gradient of the order parameter and tends to make the interface between the two phases smooth. Hence, the length of this diffuse interface is given by √γ. The second term ψ(n) is related to the mechanical interactions between the cells, and is called the homogeneous free energy. Attractive and repulsive forces for the two cell types are represented by this term. Therefore, attraction occurs when ψ′′(n) < 0 and repulsion when ψ′′(n) > 0. For volume fractions such that ψ′′(n) = 0, we say that the mixture is at equilibrium, i.e. attractive and repulsive forces balance out.
We define µ = µ(n, n) the chemical potential as the variational derivative of the free energy with respect to the order parameter µ = δE
Going back to the definition of the net flux h, and using a generalized Fick’s law, we obtain h = −b(n)∇µ, where b(n) is a mobility coefficient to represent the active movement of cells. Altogether, and assuming zero-flux boundary conditions, the original Cahn-Hilliard equation for diphasic fluids reads
∂n = div (b(n)∇ (−γ∆n + ψ′(n))) , t > 0, x ∈ Ω,
∂µ = ∂n = 0, t > 0, x ∈ ∂Ω,
where ν is the outward normal vector to the boundary ∂Ω.
The Cahn-Hilliard for Biology. This equation found its original application in the context of material sciences. Cahn and Hilliard [48, 47] first proposed the equation to represent the separation of phases occurring in binary alloys during a sudden cooling, assuming isotropy and constant temperature. Indeed, the equation can model the different stages of phase separation: from the formation of microscopic structures (i.e. spinodal decomposition) to the coarsening of them to form large arrangements. Later, the Cahn-Hilliard equation has been used to represent many phenomena in Physics, such as dealloying occurring in corrosion [81], thin films [191], image processing [53], or even about the formation of the rings of Saturn [196]. The previous references do not cover the extensive literature that exists for each of these subjects. However, the interested reader can find more references in the review book of Miranville [147].
In this manuscript, we are interested in its application to Biology, especially to represent tissues and tumors. Due to qualitative similarities between the formation of patterns in nature and the processes of phase separation occurring in materials, researchers started to use the Cahn- Hilliard equation as a phenomenological model for biological applications, such as population dynamics [62], wound healing [126], tumor growth [203, 90] or even the organization of mussel banks [138]. In the previous references, a source term g(n) is often considered inside the equation to represent growth of the tissue or its degradation, leading to the generalized Cahn-Hilliard equation [59]
∂n = div (b(n)∇ (−γ∆n + ψ′(n))) + g(n), t > 0, x ∈ Ω. (1.5)
Both [62] and [126] used the Cahn-Hilliard to produce results close to observed phenomenon with- out focusing on the consistency of the model with their system thermodynamics. The first work that proposed a mechanically and thermodynamically consistent derivation of the CH equation for Biology is [203]. The authors derived a system of CH-type equations using thermodynamics laws to derive constraints for the constitutive relations of their model. Their model is capable of representing multi-species systems, and they considered as a test case the representation of a system with four species: viable tumor cells, dead tumor cells (necrotic core), healthy cells, and the rest of the micro-environment. In the previous references that propose a CH model for living tissues, the mobility coefficient (i.e. b(n) in (1.5)) is taken to be a constant and the potential ψ(n) is double-welled and logarithmic, i.e. it describes the segregation of the two different types of cells. For this application, a thermodynamically consistent potential used in many research pieces is ψ(n) = 1 n ln n + (1 − n) ln(1 − n) − (n − 1 )2, but is often approximated by a polynomial function for simplicity.
Later, other forms for the mobility and potential have been proposed to give a better rep- resentation of the mechanics of cells inside tissues. Considering a constant mobility in the CH model does not seem biologically relevant. Indeed, to represent the clustering of cells, a degener- ate mobility is more relevant as pointed in [8, 6], and consists in taking the mobility to be zero in pure phase i.e. in zones where cells are too overcrowded (n = 1 and n = 0). A possible example is b(n) = n(1 − n)2.
Featuring this kind of mobility the CH is often referred to as the degenerate Cahn-Hilliard model (DCH in short). Another modification that renders the model biologically relevant for certain applications, concerns the potential ψ(n). As said previously, this potential is often taken to be a double-well logarithmic function with singularity points at the pure phases. However, as proposed in [46], a single-well logarithmic potential seems to be more relevant for the modeling of cancer. Indeed, tumor cells have the capacity to spontaneously form aggregates, and are often the only active cells in the experiments. For example, modelling the formation of in-vitro spheroids of cancerous cells, one can consider the two phases of the fluid to be the tumor cells and the other phase is the inactive gel that serves as the culture medium. For this kind of study, a phenomenological potential is ψ(n) = −(1 − n ) ln(1 − n) − 3 − (1 − n ) 2 − (1 − n )n, in which 0 < n٨ 1 is a parameter used to represent the maximal density at which aggregates are stable, i.e. attraction and repulsion forces balance out.
The DCH with a single-well logarithmic homogeneous free energy has been used for the modelling of skin cancer [55, 56, 61], and of glioblastoma [9]. In both of these studies, the results of numerical simulations of the model were in good agreement with biological observations and experimental results. Furthermore, in [10], the authors used their mathematical and algorithm framework, which involves the DCH model in a concrete study with data coming from a patient suffering from glioblastoma. The numerical simulations were able to represent the evolution of the tumor of the patient even under treatment. The mathematical model was able to match the volume and boundaries of the tumor observed in images realized by Magnetic Resonance Imaging (MRI).
Although the previous Cahn-Hilliard framework provides a good representation of a tumor in a healthy tissue or in vitro, as a diphasic fluid, it is necessary to consider a larger number of cell types and include particular mechanical effects in some applications. To this end, in [98] the authors derived a system of equations that comprises a CH-type equation to represent healthy and tumor cells, coupled to an equation for the diffusion of the nutrients in the tissue. This model is derived based on simple thermodynamics principles and represents nutrient diffusion, chemotaxis, active transport, adhesion between cells, apoptosis (i.e., death of the cells), and proliferation. Therefore, in the equation for the chemical potential, the effect of nutrients is taken into account to represent the chemotactic movement. A velocity field is also taken into account to represent non-active movement (i.e., advection). It is given by Darcy’s law which states that the movement is towards zones of lower pressure. The extension of the model for multispecies systems has been proposed in [97]. Numerical simulations of this model give a good qualitative representation of biological phenomena, such as the emergence of a necrotic core inside the tumor. In the two previous models, the transport due to the fluid movement was given by Darcy’s law. However, to account for the fluid viscosity, a Cahn-Hilliard-Brinkman model has been proposed for tumor growth [72].

READ  Microstructural foundations of volatility properties

The Keller-Segel model and the volume-filling approach

Derivation. As for the original proposition of the Cahn-Hilliard equation, the Keller-Segel model (KS in short) has first been introduced as a phenomenological model to describe cells’ movement toward a specific signal called the chemoattractant. This type of motion for cells and bacteria is called chemotaxis. Starting from the prototypal model (1.2), we define the flux for each population by Ji = −Di(ρi)∇ρi + χ(ρi)vi, (1.6)
where Di(ρi) is a density-dependent diffusion coefficient, χ(ρi) a function which gives the strength of the chemotactic movement. Therefore, chemotaxis is described mathematically by a specific velocity field vi. Hence, vi is a function of one or many other component j = i. For a single cell population, if we consider a flux of the form (1.6) and the chemoattractant is given by a known function F (x), the model (1.2) reduces to the Fokker-Planck equation. Generally, the concentration of chemoattratant is given by a diffusion equation. To simplify the derivation, we consider in the following the biological situation of a unique cell population i = 1, and a single chemoattractant i = 2. The movement of the cells is a combination of random Brownian motion and chemotaxis. We designate by u = ρ1 the density of cells and by c = ρ2 the concentration of the chemoattractant that diffuses in the domain. The two fluxes (1.6) are given by
J1 = −D1(u, c)∇u + χ(u, c)∇c,
J2 = −D2(u, c)∇c.
Therefore, the general Keller-Segel model reads
∂u = div (D1(u, c)∇u) − div (χ(u, c)∇c) + S1(u, c),
∂c = div (D2(u, c)∇c) + S2(u, c),
and is often supplemented by zero-flux boundary conditions
∂u ∂c = ∂ν ∂ν = 0. (1.8)
Function S1(u, c) describes the growth of the proliferating cell population that depends on the local cell density (since the proliferation of existing cells gives growth) and, possibly, on the chemoattractant concentration. The second source term S2(u, c) represents the production and decay of the chemoattractant. Both functions D1(u, c) and D2(u, c) are diffusion coefficients of respectively the cells and the chemoattractant. These two latter functions are often considered to be constant coefficients. Particular attention will be given in the following to the function χ(u, c) referred to as the chemosensitivity and describes the strength of chemotaxis. This function can depend both on the density of cells (as for the mobility in the CH model case that we have seen before) and on the concentration of chemoattractant. In the form (1.7), the KS model is a system of two parabolic equations and is often referred to as the parabolic-parabolic Keller-Segel model. If the diffusion of the chemoattractant is very fast compared to the Brownian motion of the cells i.e. D1(u,c) ≡ 0, then the parabolic-parabolic KS model is approximated by ∂u = div (D1(u, c)∇u) + div (χ(u, c)∇c) + S1(u, c), (1.9) 0 = div (D2(u, c)∇c) + S2(u, c), and is called the parabolic-elliptic Keller-Segel model (since the second equation is now of elliptic type). Further approximations of this system are possible, but are not covered in this manuscript.
The Keller-Segel model. In the article [164], Patlak studied the effect of an external bias on the movement of particles. From a model that describes the cells individually, Patlak obtained a modified Fokker-Planck equation. Independently, Keller and Segel, in the 70’s [123, 122, 124] worked on the modeling of the aggregation of cellular slime mold called acrasiales. Indeed, biolo- gists observed a long-range effect in the morphogenesis of these aggregates depending on chemical cues present in the environment. Therefore, incorporating the different reactions for their appli- cation, they wrote a four equations system that they reduced to the general system (1.7). This was the first time this model was proposed in this form. Since then, mathematicians have started to analyze its analytical properties, and find new ways to simulate it efficiently.
One of the simplest version of the general model (1.7), is the minimal Keller-Segel model, and is obtained by taking D1(u, c) = D1 0, D2(u, c) = D2 0, and χ(u, c) = u. Altogether, the model reads
∂u = D1∆u + div (u∇c) ,
∂t = D2∆c + u − c,
(1.10)
where τ = D1 describes how fast the chemoattractant diffuses compared to the Brownian motion of the cells. This relatively simple model is capable to represent the aggregation of a constant mass of cells due to a chemical signal c produced by the cells themselves. The properties of the solutions of this model have been studied by many authors. One of the most important results concerns the potential blow-up of the solution u in finite time. Indeed, from [158] we know that if d = 1, the model (1.10) has a global weak solution, and u remains bounded. However, for dimension d ≥ 2, it exists a critical mass M such that a global weak solution exists if ∫ u(0, x) dx ≤ M.
The precise value of this critical mass has been found for d = 2 [158], and is M = 4π. For d ≥ 3 [51], if u0 L 2 (Ω) , ∇c0 Ld (Ω) ≤ ε where ε > 0, it exists a global and bounded weak solution u, c . From numerical simulations, we observe that the system forms patterns that appear in Turing-type instabilities. Mainly, from a uniform initial cell density, spikes of large cell density will form, and if the simulation is run long enough, all the cells aggregates in a single sharp peak. However, it is easy to understand that this behavior is not biologically relevant. Indeed, the formation of overcrowded zones is a scenario that cells tend to avoid. To solve this issue and represent chemotaxis in living organisms, variants of the KS model (1.7) have been proposed.
To avoid the blow-up of the solution, particular forms of the functions D1(u, c), D2(u, c), and χ(u, c) can be chosen. Indeed, in [162], Painter and Hillen proposed modifying the functions for chemosensitivity to take into account the effects of « volume-filling » and « quorum-sensing ». The first effect is motivated by the fact that cells have a finite size. Hence, they cannot aggregate indefinitely at a certain point, i.e. the space available for new cells to move at a specific location decreases as the density increases. Therefore, taking into account this assumption, cells tend to form aggregates that have a finite saturation value. The second effect, which is « quorum-sensing », captures how cells behave to achieve homeostasis. Indeed, biological tissues organize themselves in such a way as to avoid excessive cell densities that may result in a depletion of important nutrients and, hence, necrosis. To model this effect, a supplementary chemical w is introduced in the system to allow the cells to sense if the zones are overcrowded and change, consequently the chemosensitivity. In [162], the authors started from a model with continuous-time and discrete space. The domain is divided into discrete locations, and cells have a probability to jump to a neighboring location depending on the chemoattractant. The authors added a set of rules to take into account the effect of volume-filling and quorum-sensing. Then, they formally derived a PDE model and retrieved a Keller-Segel system with a particular form for the chemosensitivity.
Indeed, a simple modification of the original system (1.10) can prevent the blow of the so- lution. Taking χ(u) = χcu(1 − u/u) (where χc is a positive constant and u is the saturation

Table of contents :

1 Introduction 
1.1 Motivations
1.2 Mathematical representation of living tissues
1.2.1 First steps of mathematical modelling
1.2.2 The Cahn-Hilliard model for modelling of tissues and tumours
1.2.3 The Keller-Segel model and the volume-filling approach
1.3 General assumptions and preliminaries
1.4 Numerical simulation: foundations
1.5 Summary of the thesis
1.5.1 Towards an efficient numerical scheme for the degenerate Cahn-Hilliard model for Biology
1.5.2 Modelling of specific scenarios in Biology
1.5.3 Structure-preserving numerical method for nonlinear models
1.6 Discussion and perspectives
1.6.1 Simulation of the relaxed-degenerate Cahn-Hilliard model and effect of the relaxation
1.6.2 Support a deeper understanding of key mechanisms in tumor progression
I The Cahn-Hilliard equation for Biology 
2 Relaxation of the Cahn-Hilliard equation for Biology
2.1 Introduction
2.2 The regularized problem
2.2.1 Regularization procedure
2.2.2 Existence for the regularized problem
2.2.3 Energy, entropy and a priori estimates
2.2.4 Inequalities
2.3 Existence: convergence as ! 0
2.4 Convergence as ! 0
2.5 Long-time behavior
2.6 Conclusion
3 Structure-preserving numerical method for the relaxed-degenerate Cahn- Hilliard model 
3.1 Introduction
3.2 Notations
3.3 Definition of the regularized problem
3.4 Nonlinear semi-implicit scheme
3.4.1 Description of the nonlinear numerical scheme
3.4.2 Well-posedness of the regularized problem and stability bounds
3.4.3 Well-posedness of the non regularized problem and stability
3.4.4 Convergence analysis
3.5 Non-linear semi-implicit multi-dimensional upwind numerical scheme
3.6 Linearized semi-implicit numerical scheme
3.7 Numerical simulations
3.7.1 Numerical results: test cases
3.7.2 Effect of the relaxation parameter
3.8 Conclusion
3.A Proof of M-matrix properties in the 1D and 2D cases
II Modification of existing nonlinear PDE models, numerical simulation, and application in Biology.
4 Treatment-induced shrinking of tumour aggregates: A nonlinear volumefilling chemotactic approach 
4.1 Introduction
4.2 Description of the experiments
4.3 Mathematical model
4.3.1 Volume-filling approach for chemotaxis: first part P1
4.3.2 PDE system including the treatment: Part P2
4.4 Linear stability analysis and pattern formation
4.4.1 Dimensionless model
4.4.2 First part: Formation of the aggregates
4.4.3 Second part: Treatment
4.5 Numerical simulations
4.5.1 Biological relevance of the model parameters
4.5.2 Numerical results for a one dimensional case
4.5.3 Numerical results for a two dimensional case
4.6 Discussion of results and perspectives
4.A Derivation of the general model
4.B Stability analysis
4.C Description of the numerics
4.D One dimensional numerical results
5 Compressible Navier-Stokes-Cahn-Hilliard model for the modelling of tumor invasion in healthy tissue
5.1 Introduction
5.2 Derivation of the model
5.2.1 Notation and definitions
5.2.2 Mass balance equations
5.2.3 Balance of linear momentum
5.2.4 Energy balance
5.2.5 Entropy balance and Clausius-Duhem inequality
5.2.6 Constitutive assumptions and model equations
5.2.7 Summary of the model equations
5.3 General assumptions and biologically relevant choice of the model functions
5.3.1 General forms and assumptions
5.3.2 Biologically consistent choice of functions
5.3.3 Non-dimensionalized model
5.4 Large friction hypothesis
5.5 Finite volume numerical scheme
III Structure-preserving numerical method for nonlinear PDEs 
6 The Scalar Auxiliary Variable method for the volume-filling Keller-Segel model. 
6.1 Introduction
6.2 Numerical scheme
6.2.1 Finite element framework
6.2.2 Fully discrete scheme
6.2.3 Matrix formulation
6.2.4 Upwind stabilization
6.2.5 Solving Algorithm
6.3 Existence of a non-negative solution and stability bound
6.3.1 Existence of a discrete non-negative solution
6.3.2 Discrete energy a priori estimate
6.4 Numerical results
6.4.1 1D numerical results
6.5 Conclusion
7 Conservation properties and long time behavior of the Scalar Auxiliary Variable method for nonlinear dispersive equations. 
7.1 Introduction
7.2 Numerical scheme
7.2.1 Time and space discretisation of the SAV model
7.2.2 The fully discrete SAV scheme
7.3 Conservation properties and inequalities
7.4 Convergence analysis
7.4.1 Notations
7.4.2 Convergence theorem
7.5 Error analysis
7.6 Numerical experiments
7.6.1 First test case: cubic nonlinearity
7.6.2 Second test case: cubic nonlinearity with non-smooth initial condition
7.6.3 Third test case: non-integer exponent
7.6.4 Computing ground states
7.A Gradient flow with discrete normalization for computing ground state
Bibliography

GET THE COMPLETE PROJECT

Related Posts