RESEARCH

Axis 1 – Multiscale study of interacting particle systems

Graph Limit for Interacting Particle Systems on Weighted Random Graphs

In 28, we studied the large-population limit of interacting particle systems posed on weighted random graphs. In that aim, we introduced a general framework for the construction of weighted random graphs, generalizing the concept of graphons. We proved quantitative convergence results in probability, as the number of particles tends to infinity, of the finite-dimensional system towards the solution of a deterministic graph-limit equation. In this limit equation, the graphon prescribing the interaction is given by the first moment of the weighted random graph law. We also studied interacting particle systems posed on switching weighted random graphs, which are obtained by resetting the weighted random graph at regular time intervals. We revealed the interplay between the large-population limit and the switching time. In particular, we show that for a fixed switching time, these systems converge as the number of particles tends to infinity to the same graph-limit equation, in which the interaction is prescribed by the constant-in-time first moment of the weighted random graph law.

Localization limit for multispecies cross-diffusions models

In 34, we studied the localization limit of models describing the long-range interaction between individuals, of particular interest for living systems. The starting point models are quadratic, written under the form of transport equations with a nonlocal self-generated drift. We established the localisation limit, that is the convergence of nonlocal to local systems, when the range of interaction tends to 0. These theoretical results are sustained by numerical simulations. The major new feature in our analysis is that we do not need diffusion to gain compactness, at odd with the existing literature. The central compactness result is provided by a full rank assumption on the interaction kernels. In turn, we prove existence of weak solutions for the resulting system, a cross-diffusion system of quadratic type.

Asymptotic preserving schemes for nonlinear kinetic equations

In 40, we studied the diffusion limit of nonlinear kinetic equations modeling the pattern formation in bacterial populations moving by a chemotaxis process. We first proved, by formal arguments, that the diffusion limit of these equations, when both the transport term and the turning operator are density-dependent, lead to volume-exclusion chemotactic equations. We generalise an asymptotic preserving scheme for such nonlinear kinetic equations based on a micromacro decomposition. By properly discretizing the nonlinear term implicitly-explicitly in an upwind manner, the scheme produces accurate approximations also in the case of strong chemosensitivity. We showed, via detailed calculations, that the scheme presents the following properties: asymptotic preserving, positivity preserving and energy dissipation, which are essential for practical applications. We extend this scheme to two dimensional kinetic models and we validate its efficiency by means of 1D and 2D numerical experiments of pattern formation in biological systems.

Derivation of the degenerate Cahn-Hilliard equation

The degenerate Cahn-Hilliard equation can be derived from interacting particle systems with short and long range interactions. These scales lead to consider non-local versions of the equation as well as the local limit, 16. The degenerate Cahn-Hilliard equation describes surface tension between two types of cells and this produces a pressure jump at the interface in the compressible limit which is computed in 15.

Axis 2 – Stochastic models for biological and chemical systems

Hawkes Processes

In 19 a general study of Hawkes processes on $ℝ$ has been achieved. A Hawkes process is a point process whose intensity function at time $t$ is a functional of its past activity before time $t$ . It is defined by its activation function $\Phi$ and its memory function $h$ . In this work, the Hawkes property is expressed as an operator on the sub-space of non-negative sequences associated to distances between its points. By using the classical correspondence between a stationary point process and its Palm measure, we establish a characterization of the corresponding Palm measure as an invariant distribution of a Markovian kernel. We prove that if $\Phi$ is continuous and its growth rate is at most linear with a rate below some constant, then there exists a stationary Hawkes point process. The classical Lipschitz condition of the literature for an unbounded function $\Phi$ is relaxed. Our proofs rely on a combination of coupling methods, monotonicity properties of linear Hawkes processes and classical results on Palm distributions. An investigation of the Hawkes process starting from the null measure, the empty state, on ${ℝ}_{–}$ plays also an important role. The linear case of Hawkes and Oakes is revisited at this occasion.

If the memory function $h$ is an exponential function, under a weak condition it is shown that there exists a unique stationary Hawkes point process. In this case, its Palm measure is expressed in terms of the invariant distribution of a one-dimensional Harris ergodic Markov chain. When the activation function is a polynomial $\Phi$ with degree $>1$ , there does not exist a stationary Hawkes process and if the Hawkes process starts from the empty state, a scaling result for the accumulation of its points is obtained.

Stochastic Chemical Networks

The general goal of this PhD work, started in September 2020, is of developing a scaling approach to analyze stochastic models of chemical networks.

Stability of Stochastic Chemical Reaction Networks

In 47, we have investigated the asymptotic properties of Markov processes associated to stochastic chemical reaction networks (CRNs) driven by the kinetics of the law of mass action. Their transition rates exhibit a polynomial dependence on the state variable, with possible discontinuities of the dynamics along the boundary of the state space. As a natural choice to study stability properties of CRNs, the scaling parameter considered in this paper is the norm of the initial state. Compared to existing scalings of the literature, this scaling does not change neither the topology of a CRN, nor its reactions constants. Functional limit theorems with this scaling parameter can be used to prove positive recurrence of the Markov process. This scaling approach also gives interesting insights on the transient behavior of these networks, to describe how multiple time scales drive the time evolution of their sample paths for example. General stability criteria are presented as well as a possible framework for scaling analyses. Several simple examples of CRNs are investigated with this approach. A detailed stability and scaling analyses of a CRN with slow and fast timescales is worked out.

Allocation of Resources in Prokaryotic Cells

The information flow from DNA genes to proteins is a fundamental process that is common to all living organisms. We analyze it in the context of a bacterial cell. The production of proteins is the most important process that takes place inside a bacterium consuming around 90% of its resources.

Bacterium adapts very quickly to changes in the environment. It succeeds in making a balance between supply allocating the resources available in the environment for its growth. All of this is achieved via regulatory mechanisms in the bacterium.

A Model of Regulation of Transcription

In 17, we studied an important global regulation mechanism of transcription of biological cells using specific macro-molecules, 6SRNAs. The functional property of 6SRNA is of blocking the transcription of RNAs when the environment of the cell is not favorable. We investigate the efficiency of this mechanism with a scaling analysis of a stochastic model. The evolution equations of our model are driven by the law of mass action and the total number of polymerases is used as a scaling parameter. Two regimes are analyzed: exponential phase when the environment of the cell is favorable to its growth, and the stationary phase when resources are scarce. In both regimes, by defining properly occupation measures of the model, we prove an averaging principle for the associated multi-dimensional Markov process on a convenient timescale, as well as convergence results for “fast” variables of the system. An analytical expression of the asymptotic fraction of sequestrated polymerases in stationary phase is in particular obtained. The consequences of these results are discussed.

A Stochastic Analysis of Particle Systems with Pairing

In 41, motivated by a general principle governing regulation mechanisms in biological cells, we investigate a general interaction scheme between different populations of particles and specific particles, referred to as agents. Assuming that each particle follows a random path in the medium, when a particle and an agent meet, they may bind and form a pair which has some specific functional properties. Such a pair is also subject to random events and it splits after some random amount of time. In a stochastic context, using a Markovian model for the vector of the number of paired particles, and by taking the total number of particles as a scaling parameter, we study the asymptotic behavior of the time evolution of the number of paired particles. Two scenarios are investigated: one with a large but fixed number of agents, and the other one, the dynamic case, when agents are created at a bounded rate and may die after some time when they are not paired. A first order limit theorem is established for the time evolution of the system in both cases. The proof of an averaging principle of the dynamic case is one of the main contributions of the paper. Limit theorems for fluctuations are obtained in the case of a fixed number agents. The impact of dynamical arrivals of agents on the level of pairing of the system is discussed.

Axis 3 – Theoretical analysis of nonlinear partial differential equations (PDE) modelling various structured population dynamics

Applications to cancer-immune interactions and to evolutionary biology in multicellularity and in cancer.

Mentioned in the Section ‘Adaptive phenotype-structured cell population dynamics’ of the research program, the PhD theses of Zineb Kaid, on interactions of a tumour cell population with populations of immune cells, NK- and T-lymphocytes, and of Frank Ernesto Alvarez Borges, on phenotype divergence in early development of multicellular animals and in tumours, have been defended, respectively at Tlemcen University (Algeria) in October, and at Paris-Dauphine University, in December, both under the mentoring of Jean Clairambault. It is expected that the collaborative research work with the latter will be pursued, developing more in-depth the idea, proposed in 27, that specialisation and cooperation between cell lines, and later between tissues, is the result of an optimisation process that leads to division of work in the development of organised multicellularity in animals, a process which is partially, nevertheless in an instable and reversible way, conserved in tumours, in which cooperation has indeed been evidenced between subpopulations of cells.

Long-time behaviour of an advection-selection equation

In 43 we studied the long-time behaviour of the advection-selection equation

with $\rho \left(t\right)={\int }_{{ℝ}^d}n(t,x)dx$ and with an initial condition $n(0,·)=n^0$ . In the field of adaptive dynamics, this equation typically describes the evolution of a phenotype-structured population over time. In this case, $x\mapsto n(t,x)$ represents the density of the population characterised by a phenotypic trait $x$ , the advection term ` $\nabla ·\left(f\left(x\right)n(t,x)\right)$ ‘ a cell differentiation phenomenon driving the individuals toward specific regions, and the selection term ` $\left(r\left(x\right)–\rho \left(t\right)\right)n(t,x)$ ‘ the growth of the population, which is of logistic type through the total population size $\rho \left(t\right)={\int }_{{ℝ}^d}n(t,x)dx$ .

In the one-dimensional case $x\in ℝ$ , we prove that the solution to this equation can either converge to a weighted Dirac mass or to a function in $L^1$ . Depending on the parameters $n^0$ , $f$ and $r$ , we determine which of these two regimes of convergence occurs, and we specify the weight and the point where the Dirac mass is supported, or the expression of the $L^1$ -function which is reached. This work was part of the PhD thesis of Jules Guilberteau, co-directed with C. Pouchol.

Infinite times renewal equation

In neuroscience, the time elapsed since the last discharge has been used to predict the probability of the next discharge. Such predictions can be improved taking into account the last discharge times. Such multi-time processes arise in many other areas and there is no universal limitation on the number of times to be used. This observation leads us to study the infinite-times renewal equation as a simple model to understand the meaning and properties of such partial differential equations depending on an infinite number of variables. We define notions of solutions, prove existence and uniqueness of solutions and prove the long time convergence, with exponential rate, to the steady state in different, strong or weak, topologies depending on the coefficients. See 32.

Work with Xuan Dou, Chenjiayue Qi, Delphine Salort and Zhennan Zhou.

Wasserstein contraction for the stochastic Morris-Lecar neuron model

Neuron models have attracted a lot of attention recently, both in mathematics and neuroscience. In 45 we study long-time and large-population emerging properties in a simplified toy model. From a mathematical perspective, this amounts to study the long-time behaviour of a degenerate reflected diffusion process. Using coupling arguments, the flow is proven to be a contraction of the Wasserstein distance for long times, which implies the exponential relaxation toward a (non-explicit) unique globally attractive equilibrium distribution. The result is extended to a McKean-Vlasov type non-linear variation of the model, when the mean-field interaction is sufficiently small. The ergodicity of the process results from a combination of deterministic contraction properties and local diffusion, the noise being sufficient to drive the system away from non-contractive domains.

Work with Maxime Herda and Pierre Monmarché.

Axis 4 – Mathematical epidemiology

Biological control of vectors

Sex-structured model for the design of sex-biased release strategies

Laboratory experiments as well as some field essays have revealed that the intracellular bacterium Wolbachia , deliberately introduced in Aedes spp female mosquitoes, drastically reduces their vector competence for dengue virus and other pathogens. However, female mosquitoes infected with Wolbachia still need to ingest human blood while male mosquitoes, either wild or Wolbachia -carrying, do not bite people. Moreover, Wolbachia -carrying females may transmit the virus to people during blood-feeding, even though with a far less probability than the wild ones. Therefore, massive releases of Wolbachia -carrying females may increase both the nuisance and the epidemiological risk among human residents. In the paper 11, we proposed a sex-structured model of Wolbachia invasion that brings forward the possibility of developing male-biased release strategies of Wolbachia -carriers leading to Wolbachia invasion. Thanks to this model, we studied the minimal amount of mosquitoes necessary to complete this task, according to the relative sex-ratio of the released mosquitoes and to the release schedule. We also paid attention to the estimate of the time needed to achieve the ultimate population replacement.

Estimating disturbing effect of mosquitoes migration in Wolbachia infection strategies

In 26, we investigated an initial-boundary-value problem of a reaction-diffusion equation in a bounded domain with a Robin boundary condition and introduce some particular parameters to consider the non-zero flux on the boundary. This problem arises in the study of mosquito populations under the intervention of the population replacement method, where the boundary condition takes into account the inflow and outflow of individuals through the boundary. Using phase-plane analysis, we studied the existence and properties of non-constant steady-state solutions depending on several parameters. We then used the principle of linearized stability to prove some sufficient conditions for their stability. We showed that the long-time efficiency of this control method depends strongly on the size of the treated zone and the migration rate. To illustrate these theoretical results, some numerical simulations were provided in the framework of mosquito population control.

Control strategy for Sterile Insect Techniques using exponentially decreasing releases to avoid the hair-trigger effect

In 48, we introduced a control strategy for applying the Sterile Insect Technique (SIT) to eliminate the population of Aedes mosquitoes which are vectors of various deadly diseases like dengue, zika, chikungunya… in a wide area. We used a system of reaction-diffusion equations to model the mosquito population and study the effect of releasing sterile males. Without any human intervention, and due to the so-called hair-trigger effect, the introduction of only a few individuals (eggs or fertilized females) can lead to the invasion of mosquitoes in the whole region after some time. To avoid this phenomenon, our strategy has been to keep releasing a small number of sterile males in the treated zone and move this release forward with a negative forcing speed $c$ to push back the invasive front of wild mosquitoes. By using traveling wave analysis, we showed that the strategy succeeds in repulsing the population while consuming a finite amount of mosquitoes in any finite time interval even though we treat a moving half-space $\{x>ct\}$ . Moreover, we succeeded in constructing a ‘forced’ traveling wave for our system moving at the same speed as the releases. We also provided some numerical illustrations for our results.

Optimal outbreak control via instant releases

Optimal release strategies to control vector-borne diseases, such as dengue, Zika, chikungunya and malaria have been studied in 25. Two techniques were considered: the sterile insect one (SIT), which consists in releasing sterilized males among wild vectors in order to perturb their reproduction, and the Wolbachia one (presently used mainly for mosquitoes), which consists in releasing vectors, that are infected with a bacterium limiting their vector capacity, in order to replace the wild population by one with reduced vector capacity. In each case, the time dynamics of the vector population is modeled by a system of ordinary differential equations in which the releases are represented by linear combinations of Dirac measures with positive coefficients determining their intensity. We introduced optimal control problems that we solve numerically using ad-hoc algorithms, based on writing first-order optimality conditions characterizing the best combination of Dirac measures. We then discussed the results obtained, focusing in particular on the complexity and efficiency of optimal controls and comparing the strategies obtained.

Control of infectious diseases

Behavioural epidemiology

Tiered social distancing policies have been adopted by many governments to mitigate the harmful consequences of COVID-19. Such policies have a number of well-established features i.e., they are short-term, adaptive (to the changing epidemiological conditions), and based on a multiplicity of indicators of the prevailing epidemic activity. In 10, we used ideas from Behavioural Epidemiology to represent tiered policies in an SEIRS model by using a composite information index including multiple indicators of current and past epidemic activity mimicking those used by governments during the COVID-19 pandemic, such as transmission intensity, infection incidence and hospitals’ occupancy. In its turn, the dynamics of the information index is assumed to endogenously inform the governmental social distancing interventions. The resulting model is described by a hereditary system showing a noteworthy property i.e., a dependency of the endemic levels of epidemiological variables from initial conditions. This is a consequence of the need to normalise the different indicators to pool them into a single index. Simulations suggested a rich spectrum of possible results. These include policy suggestions and identify pitfalls and undesired outcomes, such as a worsening of epidemic control, that can arise following such types of approaches to epidemic responses.

Computation and properties of the epidemic final size

The final infection size is defined as the total number of individuals that become infected throughout an epidemic. Despite its importance for predicting the fraction of the population that will end infected, it does not capture which part of the infected population will present symptoms. Knowing this information is relevant because it is related to the severity of the epidemics. The objective of the work 9 has been to give a formula for the total number of symptomatic cases throughout an epidemic. Specifically, we focused on different types of structured SIR epidemic models (in which infected individuals can possibly become symptomatic before recovering), and we computed the accumulated number of symptomatic cases when time goes to infinity using a probabilistic approach. The methodology behind the strategy followed is relatively independent of the details of the model.

Optimal control of social distancing

We revisited in 21 the problem of minimizing the epidemic final size in the SIR model through social distancing. In the existing literature, this problem has been considered by imposing a priori interval structure on the time period during which interventions are enforced. We showed that when considering the more general class of controls with an L1 constraint on the confinement effort to reduce the infection rate, the support of the optimal control is still a single time interval. In other words, for the problem of minimizing the epidemic final size in the SIR model through social distancing, there is no benefit in splitting interventions on several disjoint time periods. The techniques deployed are different from the ones used in the literature, and could be potentially applied to other problems.

This paper was one of the six papers of an Invited session dedicated to Control theory perspectives on mathematical epidemiology that we organized with Alain Rapaport (Inrae) at the 2023 IFAC World Congress. See 20.

Observability and identifiability issues for models with reinfections

Observation and identification are important issues for the practical use of compartmental models of epidemic dynamics. Usually, the state and parameters of the epidemic model are evaluated based on the number of infected individuals (the prevalence) or the newly infected cases (the incidence). Other estimation techniques, for example, based on the exploitation of the proportion of primo-infected individuals (easily retrievable data), are rarely considered. We have been thus interested in a general question: may the measure of the number of primo-infected individuals and the prevalence improve simultaneous state and parameter estimation? In the paper 23, we designed a nonlinear adaptive observer for a simple infection model with waning immunity and consequent reinfections to answer this question. The practical asymptotic stability of the estimation errors has been proved using the Lyapunov function method, and the convergence of the observer has been illustrated by simulations.

Axis 5 – Development and analysis of mathematical models for biological tissues confronted to experimental data

An integrative phenotype-structured partial differential equation model for the population dynamics of epithelial-mesenchymal transition

Phenotypic heterogeneity along the epithelial-mesenchymal (E-M) axis contributes to cancer metastasis and drug resistance. Recent experimental efforts have collated detailed time-course data on the emergence and dynamics of E-M heterogeneity in a cell population. However, it remains unclear how different possible processes interplay in shaping the dynamics of E-M heterogeneity: a) intracellular regulatory interaction among biomolecules, b) cell division and death, and c) stochastic cell-state transition (biochemical reaction noise and asymmetric cell division). In 42, we propose a Cell Population Balance (Partial Differential Equation (PDE)) based model that captures the dynamics of cell population density along the E-M phenotypic axis due to abovementioned multi-scale cellular processes. We demonstrate how population distribution resulting from intracellular regulatory networks driving cell-state transition gets impacted by stochastic fluctuations in E-M regulatory biomolecules, differences in growth rates among cell subpopulations, and initial population distribution. Further, we reveal that a linear dependence of the cell growth rate on the population heterogeneity is sufficient to recapitulate the faster in vivo growth of orthotopic injected heterogeneous E-M subclones reported before experimentally. Overall, our model contributes to the combined understanding of intracellular and cell-population levels dynamics in the emergence of E-M heterogeneity in a cell population. This work was part of the PhD thesis of Jules Guilberteau, and was conducted in collaboration with Mohit Kumar Jolly and Paras Jain (Indian Institute of Science, Bangalore) and Camille Pouchol (Université Paris Cité).

Mathematical modeling of axolotl spinal cord regeneration

Axolotls are uniquely able to completely regenerate the spinal cord after amputation. The underlying governing mechanisms of this regenerative response have not yet been fully elucidated. We previously found that spinal cord regeneration is mainly driven by cell cycle acceleration of ependymal cells, recruited by a hypothetical signal propagating from the injury. However, the nature of the signal and its propagation remain unknown. In 30, we took a step forward and proposed a theoretical study in which we investigated whether the regeneration-inducing signal could follow a reaction-diffusion process. We developed a computational model, validated it with experimental data and showed that the signal dynamics can be understood in terms of reaction-diffusion mechanism. By developing a theory of the regenerating outgrowth in the limit of fast reaction-diffusion, we demonstrated that control of regenerative response solely relies on cell-to-signal sensitivity and the signal reaction-diffusion characteristic length. This study lays foundations for further identification of the signal controlling regeneration of the spinal cord.

Development and analaysis of 3D dynamical network models

The Extra-Cellular-Matrix (ECM) is a complex interconnected 3D network that provides structural support for the cells and tissues and defines organ architecture key for their healthy functioning. However, the intimate mechanisms by which ECM acquire their 3D architecture are still largely unknown. In 12, we addressed this question by means of a 3D individual based model of interacting fibers able to spontaneously crosslink or unlink to each other and align at the crosslinks. We showed that such systems are able to spontaneously generate different types of architectures, and performed a thorough analysis of the emerging structures by an exhaustive parametric analysis and the use of appropriate visualization tools and quantifiers in 3D. The most striking result is that the emergence of ordered structures can be fully explained by a single emerging variable : the proportion of crosslinks in the network. This simple variable becomes an important putative target to control and predict the structuring of biological tissues, to suggest possible new therapeutic strategies to restore tissue functions after disruption, and to help in the development of collagen-based scaffolds for tissue engineering. Moreover, the model revealed that the emergence of architecture is a spatially homogeneous process following a unique evolutionary path, and highlights the essential role of dynamical crosslinking in tissue structuring.