## Section: New Results

### Probabilistic numerical methods, stochastic modelling and applications

Participants : Mireille Bossy, Nicolas Champagnat, Madalina Deaconu, Angela Ganz, Inès Henchiri, Samuel Herrmann, Pierre-Emmanuel Jabin, Antoine Lejay, Pierre-Louis Lions, Sylvain Maire, Denis Talay, Etienne Tanré, Samih Zein.

#### Monte Carlo methods for diffusion processes and PDE related problems

Keywords : Monte Carlo methods, divergence form operator, discontinuous media, skew Brownian motion, random walk on spheres/rectangles, backward stochastic differential equations.

S. Zein started his post-doc in September 2007 on the application of optimization algorithm for the method of random walk on rectangles using importance sampling techniques [Oops!] . The idea is then to look for optimal parameters through an automatic procedure in order to reduce the variance of the Monte Carlo estimator, or to simulate rare events.

A. Lejay and S. Maire have worked on Monte Carlo computations of the first eigenvalue of linear operators like the Laplace operator or the neutron transport operator in a bounded domain [Oops!] . This method relies on a branching mechanism and improves our previous algorithm [60] , [58] : The method is more robust, more accurate and provides the first eigenelement of the adjoint operator.

A. Lejay and S. Maire are developing a Monte Carlo method to simulate a diffusion process in discontinuous media. This method relies on a kinetic approximation of the diffusion operator.

All these methods are motivated by problems in physics and finance: estimation of the electrical conductivity in the brain (MEG problem), simulation of diffusions in heterogeneous media, computation of effective coefficients for PDEs at small scales, computation of option prices, ...

#### Statistical estimators of diffusion processes

In numerous applications, one chooses to model a complex dynamical phenomenon by stochastic differential equations or, more generally, by semimartingales, either because random forces excite a mechanical system, or because time–dependent uncertainties disturb a deterministic trend, or because one aims to reduce the dimension of a large scale system by considering that some components contribute stochastically to the evolution of the system. Examples of applications respectively concern mechanical oscillators submitted to random loading, prices of financial assets, molecular dynamics.

Of course the calibration of the model is a crucial issue. A huge literature deals with the statistics of stochastic processes, particularly of diffusion processes: However, somewhat astonishingly, it seems to us that most of the papers consider that the dimension of the noise is known
by the observer. This hypothesis is often questionable: there is no reason to
*a priori* fix this dimension when one observes a basket of assets, or a complex mechanical structure in a random environment. Actually we were motivated to study this question by modelling and simulation issues related to the pricing of contracts based on baskets of energy prices
within our past collaboration with Gaz de France.

Jointly with J. Jacod (Université Paris 6), A. Lejay and D. Talay have tackled the question of estimating the Brownian dimension of an Itô process from the observation of one trajectory during a finite time interval. More precisely, we aim to build estimators which provide an
“explicative Brownian dimension”
r_{B} : a model driven by a
r_{B} Brownian motion satisfyingly fits the information conveyed by the observed path, whereas increasing the Brownian dimension does not allow to fit the data any better. Stated this way, the problem is obviously ill posed, hence our first step consisted in defining a reasonable framework
to develop our study. We developed several estimators and obtained partial theoretical results on their convergence. We also studied these estimators numerically owing to simulated observations of various models.

#### New simulations schemes for the approximation of Feynman-Kac representations

The Feynman-Kac formula is a well-known tool to achieve stochastic representations of the pointwise solution of numerous partial differential equations like diffusion or transport equations. In [55] , the authors introduced sequential Monte Carlo algorithms to compute the solution of the Poisson equation with a great accuracy using spectral methods. In [Oops!] , S. Maire and E. Tanré have proposed new variants of the Euler scheme and of the walk on spheres method for the Monte Carlo computation of Feynman-Kac representations. We optimize these variants using quantization for both source and boundary terms. Numerical tests are given on basic examples and on Monte Carlo versions of spectral methods for the Poisson equation.

#### Exact simulation of SDE

Beskos et al. [44] , [43] proposed an exact simulation algorithm for one dimensional SDE with constant diffusion coefficient and restrictive assumptions on the drift coefficient. E. Tanré and V. Reutenauer (Calyon) have adapted the results. We extend this algorithm to more general SDEs and we apply it to SDEs with non Lipschitz coefficients (for instance, the CIR model for interest rates). During her internship supervised by E. Tanré, I. Henchiri has performed numerical tests on this subject.

#### Coagulation-fragmentation equations

Keywords : Chemical kinetics, Smoluchowski equation.

The phenomenon of coagulation was first described by the Polish physicist Marian von Smoluchowski in 1916, in a paper on the precipitation in colloidal suspensions. It applies and is studied in many domains as: the formation of stars and planets, the behaviour of fuels in combustion engines, the polymerisation phenomenon, etc. Our expertise in this domain is now developed in a new direction.

More precisely, M. Deaconu and E. Tanré are trying to adapt their methodology for a problem arising in the Chilean copper industry. Copper ore is ground in crushers by using steel balls, the aim being to fragment the ore by using the least amount of energy. A stochastic algorithm that evaluates the time spent in the crusher has been already constructed.

A new study is underway, aiming to construct a model which effectively illustrates the phenomenon of fragmentation and accounts for new parameters, such as the position of particles in the crusher, the crusher geometry, the shape, size and number of steel balls used as projectiles and the damage factor. This factor takes into account the effect of low-impact collisions, which may occur before particle fracture. A balance must be found between the model's complexity and its ability to approximate the copper industry crushers. A successful outcome will allow us to optimize crushers efficiency at a minimal cost.

E. Tanré participated to the supervision of the Master internship of Tomás Neumann (Pontificia Universidad Catòlica de Chile) on this subject, defended in June 2007.

#### On the approximation of equilibrium measure of McKean-Vlasov equations

Keywords : Equilibrium measure, McKean Vlasov equation, p-Wasserstein metrics, particle systems, Euler scheme.

In her Ph.D. thesis supervised by M. Bossy and D. Talay, A. Ganz studies stochastic particle numerical methods to approximate the limit, when t goes to infinity, of the McKean Vlasov partial differential equation in dimension 1.

We suppose that b satisfies a monotonicity condition, and that is a non-linear interaction kernel and a small perturbation of b .

Previously, we have used the method introduced in
[48] and written
V(
x,
t) as a cumulative distribution function of a nonlinear stochastic process
X_{t} . We have also established the existence and uniqueness of the equilibrium measure of
X_{t} , and obtained an exponential convergence rate to the stationary measure in the 2-Wasserstein distance. We have estimated the convergence rate of
to the distribution function of the equilibrium measure, where
is the discretization by Euler scheme for the particle system associated to
X_{t} . We also obtained a bound for the
norm of the approximation error.

In 2007 we completed the study of the dimension 1: we established the convergence to the equilibrium measure in the norm and obtained a bound for the approximation error in the norm. In and norms, our approximation converges with a rate proportional to the discretization time step, inversely proportional to , and with an exponential decay in T .

#### Evolutionary branching in a jump process of evolution

In [51] , under the assumptions of rare mutations and large population, N. Champagnat obtained a microscopic interpretation of a jump process of evolution where the evolution proceeds by a sequence of mutant invasions in the population, in the specific case where only one type can survive at any time. This process, called “trait substitution sequence” (TSS) of adaptive dynamics, is the basis of the theory of adaptive dynamics [61] . This theory provides in particular an interpretation of the phenomenon of evolutionary branching, where the distribution of the phenotypic traits under consideration in the population, initially concentrated around a single value, is driven by selection to divide into two (or more) clusters corresponding to sub-populations with different phenotypes, that stably coexist and are still in competition.

In collaboration with M. Benaïm (Université de Neuchatel, Switzerland) and S. Méléard (École Polytechnique, Palaiseau), we continued the study of a similar scaling limit of an individual-based evolutionary system, in the case where coexistence is possible. We were able to justify the convergence to a generalized TSS under the assumption of symmetric competition, thanks to results from dynamical systems. From this generalized TSS, we were also able to mathematically justify the evolutionary branching criterion proposed by biologists in the limit of small mutation, using results on the asymptotic dynamics of 3-dimensional competitive Lotka-Volterra systems [65] . We are preparing a publication.

#### Time scales in adaptive dynamics

In collaboration with A. Bovier (WIAS, Berlin), N. Champagnat is studying the different time scales of evolution as a function of the parameters of a similar model as the one of the previous paragraph. Three parameters are involved: the population size, the mutation probability and the size of mutations. When the first parameter goes to infinity and the two others go to zero, three distinct time scales appear. The first one corresponds to directional evolution, the second one to fast evolutionary branching and the last one to slow evolutionary branching. In the first phase, the number of coexisting types is (essentially) constant and in particular cannot increase. The types evolve by maximizing their relative “fitness” according to an ODE known as the “canonical equation of adaptive dynamics” [54] . Once they reach an equilibrium, evolutionary branching occurs on the second time scale similarly as in the previous paragraph. In particular we are able to compute the asymptotic branching time. Finally, on the last time scale, non-local evolutionary branching occurs. Here again, we are able to compute the asymptotic behaviour of the branching using large deviations arguments. All these time scales have been mathematically justified when the initial population is monotype in a simplified model. Some technical difficulties, linked with stability properties of dynamical systems and with convergence of stochastic processes on long time scales, still remain for the generalization of our results to more realistic models and to multitype populations.

#### The canonical diffusion of adaptive dynamics

N. Champagnat and Amaury Lambert (Laboratoire d'écologie, ENS Paris) have recently obtained [Oops!] the counterpart of the TSS described above when the population size is finite. A diffusion process has also been obtained when the size of the jumps in the TSS vanishes, together with analytical expressions of its coefficients as series. From a biological viewpoint, this diffusion process quantifies the interplay between genetic drift and directional selection, and allows us to give an interpretation of the phenomenon of punctualism, where the evolution of the population is an alternance of long phases of stasis (constant state) and short phases of quick evolution.

We are now working on the numerical computation of the coefficients of the diffusion. They involve the first-order derivatives of the probability of invasion of a mutant with respect to the mutant birth, death and competition rates. We were able to decompose these derivatives into five fundamental components that allow one to determine which mutant type is more likely to invade the population (a problem related to notion of robustness of the population). This numerical study was initiated in [Oops!] and allowed to obtain the convergence of the canonical diffusion to the canonical equation of adaptive dynamics mentioned above when the population size goes to infinity. We also plan to apply our numerical results to a bistable ecological model with two habitats.

#### Probabilistic interpretation of PDEs of evolution and evolutionary branching

In collaboration with P. Del Moral (INRIA Bordeaux) and L. Miclo (CNRS, Marseille), N. Champagnat and P.-E. Jabin are working on nonlinear Feynman-Kac formulæ giving a probabilistic interpretation of PDEs of the form

The Laplacian corresponds to mutations in the phenotype space, and the nonlinear term corresponds to birth and death with interaction. The parameter scales the mutation size. Barles and Perthame [39] showed that, in very specific cases, when 0 , converges to a single Dirac mass which support is characterized by a Hamilton-Jacobi PDE with constraint on its maximum. Using large deviations results, our probabilistic interpretation gives a variational characterization of the solution of this Hamilton-Jacobi equation, that could allow one to extend the results of Barles and Perthame to cases where several phenotypes can coexist in the population, i.e. where evolutionary branching can occur.

#### Two-types birth and death processes conditioned on non-extinction

In collaboration with P. Diaconis (Stanford University and Université de Nice – Sophia Antipolis) and L. Miclo (CNRS, Marseille), N. Champagnat is studying two-types birth and death processes conditioned on non-extinction. The motivation is to study the competition between two populations (for example a resident one and a mutant one) and to study which population survives when the whole population does not go extinct for a long time. We obtained criterions on the eigenvalues of sub-matrices of the transition matrix that characterizes the cases where one population goes extinct or when both populations survive. In the symmetric cases, where both populations are identical, we obtained a full description of the eigenvalues of the transition matrix in terms of the monotype transition matrix. We next aim at studying small deviations from the symmetric case to identify the largest eigenvalue and the corresponding asymptotic behaviour in the neighbourhood of symmetry.

#### Large population scalings of individual-based models with spatial and phenotypic interactions

In [Oops!] , N. Champagnat and S. Méléard (École Polytechnique, Palaiseau) have studied some large population scalings of individual-based ecological birth-death-mutation-competition models where each particle moves in the (physical) space according to a diffusion reflected at the boundary of a fixed domain. A particular care is put on the spatial interaction range, in order to recover integro-PDE models with local spatial interactions proposed by Desvillettes et al. [53] . Several examples are given, illustrating spatial and phenotypic clustering and the importance of such models, that combine space and evolution, for the study of invasion phenomena.

#### Systems of differential equations with singular force fields

P.-E. Jabin is mainly working on systems of differential equations with singular force fields. He was recently able to obtain quantitative estimates of compactness for solutions to ODE's with a BV field and a sort of compressibility condition (weaker than bounded divergence). More precisely it is possible to show that solutions to

_{t}X( t, x) = b( X( t, x)), X(0, x) = x,

satisfy a uniform bound on a quantity like

provided that the field is BV and the Jacobian of the flow is bounded. The above quantity was already introduced in a paper by Crippa and DeLellis
[52] , but with a
W^{1,
p} condition on the field. This bound also implies well posedness for the system and it improves the well-known result of Ambrosio
[36] which required a bounded divergence.

Moreover, in collaboration with Julien Barré (Université de Nice – Sophia Antipolis) and Maxime Hauray (Université Paris 6), we used this idea of quantitative estimates for the study of interacting particles. Considering a number of particles interacting with each other through a potential (electrostatic for example), the question is to obtain the limit when the number of particles goes to infinity. This limit is classically a kinetic equation in the phase space but a rigorous proof was not known whenever the potential was singular. We prove that in a suitable sense, for almost all initial conditions, the trajectories of almost all particles are close to the trajectories predicted by the mean field limit.

#### Poisson-Boltzmann equation in molecular dynamics

M. Bossy, N. Champagnat, P.-E. Jabin, S. Maire and D. Talay are working in collaboration with T. Malliavin (Institut Pasteur, Paris) on a probabilistic interpretation of the Poisson-Boltzmann equation of molecular dynamics that allows one to compute the electrostatic potential around a bio-molecular assembly (for example a protein). We also aim at developing specific Monte Carlo methods linked to the singularities of the model (operator with a divergence form with discontinuous coefficients, singular source terms, and discontinuous non-linear terms).

#### Lower bounds for densities of hypo-elliptic diffusions

Motivated by applications to risk analysis in Finance and reliability problems in Random Mechanics, P.-L. Lions and D. Talay are getting precise estimates on lower bounds of certain hypo-elliptic diffusion processes. Lower bounds of Gaussian type are well known for strongly elliptic diffusions. The techniques are very different in the hypo-elliptic case.