Cancer, tumor growth and invasion

Gliomas are very aggressive brain tumors, in which tumor cells gain the ability to invade the surrounding normal tissue. This type of brain tumor comprises an extremely interesting paradigm of invasive tumor, where its dynamics and mechanisms are not yet fully understood. Our work is motivated by the migration/proliferation dichotomy (so-called Go-or-Grow) hypothesis, which might play a central role in the biology of these tumors.

We propose first a “Go-or-Rest” model and describe cell migration as a velocity-jump process including resting phases. By using scaling arguments we derive a continuum (macroscopic) model that provides anomalous diffusion, which is further analyzed. In particular, we show that sub- and super-diffusion regimes can be obtained, and are governed by a parameter describing intrinsic migratory properties of cells. We demonstrate the potential of the model to explain some in vitro data of glioma tumor expansion.

When proliferation is included, the framework previously used is again the base of our extended model. In particular, we developed a lattice-gas cellular automaton, which provides a discrete, stochastic description of the afore-mentioned mesoscopic framework. Our goal is to test hypotheses of cellular mechanisms involved in glioma tumor invasion. To this end, we use again data of in vitro glioma cell cultures that allow for the analysis of the invasive behavior of gliomas. The main observations of these experiments are: (i) different core and invasive radii speeds and (ii) a high radial persistency of cell motion nearby to the core. Our analysis shows that the migration/proliferation dichotomy assumption plays a central role in the expansion of glioma tumor. Morever, we find that a cell-cell repulsion is required to explain the observed radial persistency of glioma cells near to dense areas. The combination of these two mechanisms is a sufficient and necessary condition for the faithful reproduction of the experimentally observed behavior [see corresponding paper] .

Cell migration in the extracellular matrix

Cell migration through tissues is an essential feature of either normal or pathological phenomena in biology, such as embryonic development, wound healing or angio/vasculogenesis. Cells (cancer metastases, fibroblasts, endothelial cells) interact both with other cells and with the surrounding tissue (the ExtraCellular Matrix) which provides them a natural complex scaffold to adhere to during motion.

Recently, much attention has been devoted to the description of the mechanics of motions of cells as a result of their interactions with the ECM. This ECM is a complex set of macromolecules essentially composed by collagen. So it appears as a fibrous environnement whose each fibre owns its proper orientation.

The aim of the model is to account for the  individual interaction mechanisms with the collagen fibres, including both a preferential movement of cells along the collagen fibres (contact guidance) and a randomly reorientation of cell motion due to the interactions among cells. The model also includes the possibility that a chemical substance induces an internal mechanism within the cell, capable to affect its motion. The net effect of chemotaxis is thus modelled as an external force depending on a given chemical profile and on its gradient.

The starting point of the model is a Boltzmann-like transport equation at the mesoscopic scale. We describe the various interactions at the microscopic scale in order to derive a system of moment equations for the macroscopic mean variables (mass and momentum). As a way to close this system, we propose the diffusion limit of the initial transport equation and obtain a single non-linear advection-diffusion equation for the cell density. This equation is capable to account for the motion related to heterogeneity and anisotropy of the ECM.

I developed a numerical solver for the problem. This is, a symmetric splitting operator scheme of order two has been used, for space and operators, respectively. For each part of the equation, I used the finite volume method with a specific scheme related to the type of operator. A high resolution wave-propagation algorithm for spatially varying flux has been implemented for the non-linear hyperbolic part. The non-linear parabolic part is solved using a Crank-Nicholson scheme in which the implicit non-linear term is treated by a Beam and Warming scheme.

In a first step [see corresponding paper] the ECM has been considered as a stationary fixed medium that has the ability to influence the (ameboid) cell motion. In a second step, we derived an equation of evolution for the collagen fibres [see corresponding paper] that includes the proteolitic activity of moving cells to open their own way through the matrix (mesenchymal migration).

Numerical investigations have shown the ability of the model to account for heterogeneity and anisotropy of the ECM and their influcence on a spreading cell population.

We also extended the model to include taxis processes in a more biologically relevant way. There, taxis processes result in a bias of the main cell mouvement. We derived the corresponding advective term appearing in the macroscopic diffusion equation. This term generalizes classical taxis representations by accounting for anisotropy of the ECM [see corresponding book chapter].
We finally performed the modeling of a phenotypic switch bewteen two cell populations with different migratory properties. Starting again from a microscopic description of the interactions between the different populations involved (e.g. resting cells in the core of a tumor, invasive tumor cells, ECM) we proposed a macroscopic PDE system that models phenotype changes triggered by the environmental density [see corresponding paper].

Kinetic theory for multicellular systems

The scientific community is aware that the great revolution of this century is going to be the mathematical formalization of phenomena belonging to the living matter, while the revolution of the past two centuries was the same approach applied to the inert matter. In the ''Mathematical Physics group'' headed by the Professor N. Bellomo, within the Department of Mathematics of the Politecnico of Torino, we develop models for large systems of interacting entities. Mathematical methods of the kinetic theory are particularly well adapted to the description of the collective behaviour starting from a description of microscopic interactions. In the framework of the generalized kinetic theory, the microscopic state of each living entities includes, in addition to mechanical variables (typically position and velocity), also an additional biological variable called activity. This last one is related to a somehow organized and even "intelligent" behaviour, and is able to modify laws of classical mechanics. The system is thus described by a probability distribution over the microscopic state. The evolution equation is obtained, similarly to the case of the Boltzmann equation, as a balance of particles in an elementary state volume. The first difficulty consists of the modelling of the biological interactions. The second one is to describe the biological specificities of each interacting specie, accounting as well for proliferating and/or destructive phenomena.

The first step of my work consisted of the generalization of the existing mathematical framework to a discretization of the biological variable [see corresponding paper].

Then we developed an application related to tumoral expansion, and more specifically to the competition between the tumoral growth and the immune response. This is a very complex process involving tumor cells and the various populations of the immune system. When tumor cells are recognized by immune cells, a competition starts and may end up either with the destruction of tumor cells, or with the inhibition and depression of the immune system.

We reduced the complexity of the system by identifying only three different cell populations:
- the endothelial cells, which are the thin cells in the interior face of blood vessels and can be found in the circulating blood;
- the immune cells whose many species are gathered;
- a specific proteins family, called cytokines, which are a kind of messenger whose function consists of their ability to participate to the activation of the immune cells.

We described the following phenomenology: The endothelial cells have a natural trend to modify their biological state from normal to progressing (abnormal) which is the initial stage of the clonal expansion: abnormal cells may loose their programmed death ability (apoptosis) to start progressing towards metastatic states. Therefore, the biological state of these cells can be defined as normal or abnormal. Immune cells, if active, have the ability to contrast the presence of tumor cells (abnormal endothelial cells) by reducing their progression; if inhibited, they loose this ability. As already mentioned, active cytokines are a kind of messenger able to active immune inhibited cells, but may loose this ability after having interacted with these last ones. The biological state of each element of these interacting populations (cells or cytokines) may thus be described by two discrete values of the biological variable.

In the very early stage of the immune competition, proliferating and destructive phenomena may not play a relevant role. Therefore, at this stage, the mathematical modelling of the immune competition between tumor and immune cells, considering the action of cytokines, simply deals with the variation, due to interactions, of the biological state of each cell population. Taking into account only conservative encounters, the aim of this particular model is the ability to describe the onset of progressing cells as a transition from the normal state of endothelial cells. At a later stage, encounters generating proliferation or destruction processes occur. This perturbation of the cellular dynamics presented above is then modelled by non-conservative encounters able to describe this new phenomenology.

The theoretical aspect and the simplicity of the model lead to the greatest prudence with regards to the conclusions to be made. Nevertheless, despite the above simplicity, the model has shown [see corresponding paper] to be able to describe various aspects of interest in  the action of the immune system against tumor growth.

Fluid Dynamics

Mathematical and numerical studies of exact solutions of the Navier-Stokes equations

  • Final year internship for my Degree in Mechanics

Several mathematical frameworks developed for the study of atmospheric tornadoes are based on axisymmetric exact solutions of the stationary incompressible Navier-Stokes equations. The aim of these phenomenological models is essentially to describe the movement and evolution of flows within developed tornadoes.

The approach we adopted finds its roots in these analytical works. Our model describes a flow of viscous fluid in a half-space (Z positive), with the plan Z=0 corresponding to the solid ground. Moreover, in order to provide a realistic model for these vortices, the prescribed boundary conditions were endowed with a typical adhesion constraint at the ground, and we assumed a behaviour of a free vortex near the vertical axis of the modelled tornado. We extended the model first developed by J. Serrin [The swirling vortex (1972), Phil. Trans. Roy. Soc. London A 271], so as to include the possible appearance of a line of sink or source of mass at the central axis. The mathematical analysis of this particular problem enabled us to model five different families of flows [see corresponding paper], therefore providing us with a wider range of possible behaviours, whose characteristics are in accordance with meteorological observations.

  •  Master thesis

The first part of my project was dedicated to the extension of the results obtained during my final year degree internship. A mathematical analysis of the solutions, arising from problems which involve a source/sink line, enabled us to identify in details the local behaviour of velocity fields near the axis and the horizontal boundary plan [see corresponding paper]. These results were thereafter shown to correlate with the numerical analysis of that same mathematical problem. In fact, this work gave rise to the development of numerical simulations for modelling various phenomenological situations involving swirling flows, as described by, among others, J. Serrin and M. Goldshtik [Collapse in conical viscous flows (1990), J. Fluid Mech. 218]. Moreover, this provided an adequate tool for visualizing such swirling flows.

In the second part of my Master thesis, we focussed our attention on the identification of the various factors that could be at the origin of atmospheric tornadoes. Our investigation enables us to conclude that one of the mechanisms by which tornadoes are created relates to the thermomechanical coupling between the ground and the atmosphere, generated by a high temperature gradient. So as to describe this type of phenomenon, our model provided exact solutions of the Boussinesq equations. Regarding the thermal boundary conditions, various choices for the prescribed heat and temperature fluxes were studied. The resulting solutions consist of a family of velocity and temperature fields involving five parameters that we were able to test numerically.

Contributions to the modelling of a swirling half-line vortex: Application to atmospheric tornado

Following the works I performed during my Master project, which was concerned with the mathematical analysis of stationary and axisymmetric exact solutions of the Boussinesq equations, the mathematical framework that we retained was the one based on the stream-function/vorticity method. The model we adopted enabled us to describe the particular situation in which a swirling vortex of cold air descends along the axis, while some air heated by the ground rolls up around the cold cone [see corresponding paper and proceeding (in French)]. This case provides a qualitative description, of phenomenological nature, of an atmospheric tornado. Moreover, the presence of a source/sink line, by accentuating the singularity of the solution near the axis, remains an appropriate manner of modelling energy supplies, such as air condensation.

It is worth emphasizing that other problems related to the creation and overall stability of atmospheric tornadoes have given rise, over the years, to additional investigations, such as ones dealing with non-stationary solutions of the Navier-Stokes equations or non-axisymmetric solutions of these same equations.

It is worth emphasizing that other problems related to the creation and overall stability of atmospheric tornadoes have given rise, over the years, to additional investigations, such as ones dealing with non-stationary solutions of the Navier-Stokes equations or non-axisymmetric solutions of these same equations.

  • Non-stationary models

This is a generalization of the stationary models previously studied. We enlightened the necessity of introducing additional terms in the Navier-Stokes equations for modelling the non-stationarity. In order for the equations to be homogeneous, the terms were embodied in a framework in which the evolution equations in space and in time were decoupled. We subsequently tried to implement coupling mechanisms that could lead to the formation of atmospheric tornadoes, such as the thermomechanical coupling introduced in the stationary model. This mathematical framework generated a model in which the main constraint is linked to the no-slip boundary condition at the ground, and it did not enable to predict the time evolution of the swirling flow.

  • Non-axisymmetric models

Various observations on axisymmetric breaks of atmospheric tornadoesLoss of symmetry near the ground yielded the subsequent development of  non-axisymmetric model. Due to the complexity of the model obtained, we try and develop a simplified formulation of the problem based on an extension of the model of Aristov [Three-dimensional conical viscous incompressible fluid flows (1998), Fluid Dynamics 33]. That is, we modelled a vortex whose axis was not orthogonal to the horizontal boundary plan. In the specific case in which the symmetry crack is small, the matched asymptotic expansions method predicted analytical solutions for high values of Reynolds number [see corresponding paper]. These solutions were limited to flows directed downwards along the axis. Nevertheless, the model we presented will allow for the study of a wider range of flows, such as swirling vortices flowing upwards along a leaning axis.

  • High Reynolds number vortices

3D trajectoriesSo as to reach realistic orders of values for the upwards and swirl velocities within atmospheric tornadoes, it is essential to consider high values of Reynolds number. However, it should be noted that most models developed by specialists from this particular field of application of fluid mechanics do not enable for such upwards flow with high velocity to verify the no-slip boundary condition at the ground. The principal reason for this lies in the fact that these models cannot account for a discontinuity within the flow. In order to overcome this difficulty, we developed an axisymmetric model which describes a swirling flow descending along the axis and an ascending flow rolling up around it. Importantly, these two swirling flows were to satisfy, at the interface, continuity conditions for the velocity and pressure only. With this approach, we were able to obtain new families of solutions characterizing a wider range of behaviours [see corresponding paper and proceeding (in French)]. Furthermore, the model adopted enabled us to obtain, close to the axis, local solutions involving very high values of Reynolds number, thus reflecting what has been previously observed in meteorology.

  • Potential expansions 

An essential idea arising from this work is the importance of providing an adequate model for the decline of a swirling vortex. That is, when the intensity of a tornado decreases, the azimuthal symmetry disappears, in a similar way as when a drop spindle loses velocity. The approach I propose consists of modelling the weakening of a tornado through three distinct phases, namely:

- A mature phase, involving an axisymmetric geometry and a very high swirl velocity;
- An intermediate weakening phase, still axisymmetric, but involving a considerably lowered swirl velocity;
- A final phase, in which the symmetry disappears.

Each of these states may be described by the models presented in my PhD thesis. Starting from a nonlinear perturbation of the initial velocity field, the mathematical modelling of each phase would generate a velocity field that may adequately describe the final state. This would enable to obtain a model that suitably describes the loss of intensity of the rotation. The variation in the geometry of the vortex, which results from the destabilization of its axis, would ultimately leads to the death of the atmospheric tornado.

  • Azimuthal symmetry break in swirling flows

In the framework of the expansion of my PhD work, concerning a possible destabilization of the vortex axis, I developed a general model of swirling flow whose rotation axis is not orthogonal to the boundary plan representing the ground. This mean there is no restriction related to the values of Reynolds number. It led to a very complex mathematical structure resulting from the particular shape of the exact solutions I used. So I perturbed the axisymmetric model by introducing a small angular deviation from the initial position of the vortex axis. The use of the strained coordinates method enabled to obtain an analytical solution for non-swirling flows either ascending or descending along the axis. This may technically be generalized to swirling flows as well. I also showed how a loss of symmetry of the vertical shear near the ground could trigger off the rotation of an initially non-swirling flow.

Two-dimensional numerical study of two co-rotating merging vortices

The phenomenon of vortex merging is encountered in two- and three-dimensional turbulence. In the simplest (2D) case, such interaction is known to be the principal ingredient of the inverse cascade. This process is also important for the dispersion of passive scalars in meteorology or geophysics. The purpose of the study is to perform a quantitative analysis of the simplest case: a pair of two-dimensional co-rotating vortices. More specifically, the idea was to determine the dependence of the process with respect to the Reynolds number, with the advantage that numerical experiments let all fields to be available compared to experimental approaches.

Before merging, when vortices are far apart, the vortex dynamics may be described by a rotation of the two vortices around their center, similarly to point vortices, at a distance which remains constant. The vorticity contours are elliptical deformed by the local strain induced by each vortex on its companion. In addition to these two essentially non-viscous effects, an increasing radius of each vortex as well as a decreasing distance between them is observed, which leads finally to the merging. Then a central zone is created where major part of the initial vorticity is concentrated while the other part is ejected in an outer region in the form of two filamentary spiral structures. At the end of merging, a unique and bigger vortex is recovered. Direct numerical simulations of the Navier-Stokes equations, at high Reynolds number, have been performed using a pseudo-spectral algorithm based on the two-dimensional Fast Fourier Transform. This enabled us to validate the method.

Simulations for Re=4000 at t=0, t=3.05, t=7.7, t=10.54, t=72.7

Potential expansions (works performed by C. Josserand): The precise study of the loss of vorticity through the filamentary structures requires a numerical parallel approach which will allow to solve the equations on a very small grid. Once developed, the computation time will be short enough to initiate a parametric study with respect to the Reynolds number, and check out the influence of this last one on every stages of the merging which have been described above.