François Fages
email: francois.fages@inria.fr

Sylvain Soliman
email: sylvain.soliman@inria.fr

On Robustness Computation and Optimization in BIOCHAM-4

BIOCHAM-4 is a tool for modeling, analyzing and synthesizing biochemical reaction networks with respect to some formal, yet possibly imprecise, specification of their behavior. We focus here on one new capability of this tool to optimize the robustness of a parametric model with respect to a specification of its dynamics in quantitative temporal logic. More precisely, we present two complementary notions of robustness: the statistical notion of model robustness to parameter perturbations, defined as its mean functionality, and a metric notion of formula satisfaction robustness, defined as the penetration depth in the validity domain of the temporal logic constraints. We show how the formula robustness can be used in BIOCHAM-4 with no extra cost as an objective function in the parameter optimization procedure, to actually improve the model robustness. We illustrate these unique features with a classical example of the hybrid systems community and provide some performance figures on a model of MAPK signalling with 37 parameters.

Introduction

Computational systems biology aims at gaining a system level understanding of high-level biological processes from their biochemical realm. Formal methods from Computer Science have been soon introduced at the heart of this effort to go beyond mathematical modeling and master the complexity of cell processes. In particular, model-checking techniques have been used to analyze Boolean gene regulatory networks and signaling or control protein reaction networks. They have also been generalized, in particular in the hybrid systems community, to quantitative temporal logics such as MTL, MITL, STL, or in our case FO-LTL(R lin ) [START_REF] Rizk | Continuous valuations of temporal logic specifications with applications to parameter optimization and robustness measures[END_REF] for dealing with numerical constraints, fitting quantitative models to the time series data observed with increasing accuracy in biological experiments [16,[START_REF] Heitzler | Competing G protein-coupled receptor kinases balance G protein and β-arrestin signaling[END_REF], controling cell processes in real-time, or designing circuits in synthetic biology [START_REF] Courbet | Computer-aided biochemical programming of synthetic microreactors as diagnostic devices[END_REF][START_REF] Rizk | A general computational method for robustness analysis with applications to synthetic gene networks[END_REF].

One striking feature of natural biological processes is their robustness with respect to both external perturbations and their intrinsic stochasticity. Measuring and optimizing robustness is thus an important topic in biological modeling and a useful yet rare feature of modeling tools in this domain. In [START_REF] Kitano | Towards a theory of biological robustness[END_REF], Kitano gives a definition of robustness of a biological system with respect to a dynamical property and a parameter perturbation law, as the mean functionality of the system. Such a statistical definition of robustness can be evaluated in a computational model by sampling the parameters' space according to their distribution laws, and checking the property on simulation traces, i.e. by runtime verification or monitoring. Since the capability of generating simulation traces is the only requirement, this can be done in a very general setting of non-linear hybrid systems.

In [START_REF] Fainekos | Robustness of temporal logic specifications[END_REF][START_REF] Rizk | Continuous valuations of temporal logic specifications with applications to parameter optimization and robustness measures[END_REF], different notions of robustness are proposed with the definition of a continuous degree of satisfaction of a temporal logic property on a trace. Such a satisfaction degree gives quantitative information on the satisfaction of the formula, as the distance to, and penetration depth within, the validity domain of the formula. In the simplest example of a threshold formula x < c, that satisfaction degree can be defined for instance as the value c -x, with negative values representing distance to satisfaction, and positive values the margins achieved for robust satisfaction.

In this paper, we present some unique features of BIOCHAM-41 for defining hybrid systems by reaction networks with rates and events, and focus on the use of quantitative temporal logic language to specify the dynamical properties of the system, observed or wished, verify them in the model, compute parameter sensitivity indices, measure robustness, synthesize model parameters and optimize robustness. We deal with the two complementary notions of robustness mentioned above, i.e. the statistical notion of model robustness to parameter change, and the formula satisfaction robustness. We show how the formula robustness can be used in BIOCHAM-4 with no extra cost as an objective function in the parameter optimization procedure to actually improve the robustness of the model with respect to parameter perturbations. We illustrate these features with a classical example of non-linear hybrid system, the bouncing ball example, and provide some performance figures on the Huang and Ferrell's model of MAPK signaling [START_REF] Huang | Ultrasensitivity in the mitogen-activated protein kinase cascade[END_REF] with 37 parameters.

BIOCHAM Models

A BIOCHAM model is composed of a (multi)set of reactions with rate functions, and/or influences with forces, plus possibly events. Such models can be interpreted in a hierarchy of differential, stochastic, Petri net and Boolean semantics [START_REF] Fages | Influence networks compared with reaction networks: Semantics, expressivity and attractors[END_REF]. We will focus here on the differential semantics. They can be imported from model repositories such as BioModels using the BIOCHAM interface to SBML (SBML-qual for influence models), or from ODE models [START_REF] Fages | Inferring reaction systems from ordinary differential equations[END_REF] through an interface to the XPP format.

This article comes with a set of examples (MAPK signaling, cell cycle, bouncing ball) available online under the form of BIOCHAM-4 Jupyter notebooks2 .

3 Behavior Specifications in FO-LTL(R lin )

Validity Domains of FO-LTL(R lin ) Constraints

In BIOCHAM, quantitative temporal properties of the behavior of a system can be formally specified in a first-order version of Linear Time Logic with free variables, linear constraints over R, and quantifiers, named FO-LTL(R lin ). The grammar is:

φ ::= c | ¬φ | φ ⇒ ψ | φ ∧ φ | φ ∨ φ | ∃x φ | ∀x φ | X φ | F φ | G φ | φ U φ | φ W φ
where c denotes linear constraints over state variables (including Time) and free variables. The validity domain D π,φ ∈ R k of an FO-LTL(R lin ) formula φ containing k ≥ 1 variables on a finite trace π = (s 0 , . . . , s n ) can be defined by recursively bottom-up from the validity domain of the basic constraints to the complex FO-LTL(R lin ) formulae by intersection, union and complementation [START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF]. An FO-LTL(R lin ) formula is false if one validity domain is empty, valid if the validity domains of all variables are R, and satisfiable otherwise.

For instance, the formula F(A ≥ 0.2) where A is a state variable expresses that the concentration of molecule A gets greater than 0.2 at some time point in the trace (F). If needed, the precise time values where the concentration of A gets greater than the threshold value can be expressed by introducing a free variable t with an equality constraint to the real time variable, F(A ≥ 0.2 ∧ t = Time). Constraints between time variables can then relate the time of different events. The maximum value of a state variable A can be specified and set in a variable v by the formula

G(A ≤ v) ∧ F(A ≥ v), a local maximum (or plateau) by F(A ≤ v ∧ X(A = v ∧ X A ≤ v))
. FO-LTL(R lin ) formulae are very expressive and can be used to define complex oscillation properties, with pseudo-period constraints defined by delay constraints between the local maxima, and phase constraints defined by delays between the peaks of different state variables [START_REF] Fages | Temporal logic modeling of dynamical behaviors: Firstorder patterns and solvers[END_REF].

Implementation

The recursive definition of the validity domain of an FO-LTL(R lin ) formula on a finite trace is implemented in BIOCHAM-4 by generating a C program that implements the necessary loops, starting from the last time point to the first, and considering the subformulae in the bottom-up order, i.e. first from the linear constraints at the leaves, to the root of the syntactic tree. The call to the C compiler is responsible for slower response time than the previous interpreted implementation on small examples, but faster on large examples. The Parma Polyhedra Library (PPL)3 is used to solve the linear constraints and simplify the representation of validity domains by finite lists of polyhedra.

Bound constraints, i.e. constraints of the form x ≤ c or x ≥ c where x is a variable and c a constant, define boxes as a particular kind of polyhedra. In that case, the validity domains are finite union domains of boxes, since they are obtained by intersection, union, complementation and projection of boxes. However, it is worth noticing that even in the case of bound constraints, the validity domain of a temporal formula can contain an exponential number of polyhedra in the number of free variables in the FO-LTL(R lin ) formula [START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF]. The issue of trace simplification, such as keeping in the trace only the time points that are local extrema for one state variable, is also important to reduce computation time and justified in a number of practical cases [START_REF] Traynard | Trace simplifications preserving temporal logic formulae with case study in a coupled model of the cell cycle and the circadian clock[END_REF].

Formula Robustness and Satisfaction Degree

In [START_REF] Rizk | A general computational method for robustness analysis with applications to synthetic gene networks[END_REF], the continuous satisfaction degree in the interval [0, 1], of an FO-LTL(R lin ) formula φ on a trace π, was defined as the distance between the validity domain of the free variables x 1 , . . . , x k in φ and some objective values in R k given as a valuation σ of the variables x 1 , . . . , x k , i.e. a vector v σ = (v 1 , . . . , v k ). This definition enforces that the violation degree is in [0, +∞) (0 when the formula σ(φ) is satisfied) and the satisfaction degree in (0, 1] (1 when the formula is satisfied). Those notions can be generalized in order to take into account how robustly a formula is satisfied, by taking into account the penetration depth of the objective values in the validity domain, similarly to the space-time robustness defined in [START_REF] Donzé | Robust satisfaction of temporal logic over real-valued signals[END_REF] for STL: Definition 1. The violation degree vd(π, φ, σ) ∈ (-1, +∞) of an FO-LTL(R lin ) formula φ in a numerical trace π with respect to an objective valuation σ is

min v∈D π,φ d(v, v σ ) if v σ ∈ D π,φ , 1 1 + min v ∈D π,φ d(v, v σ ) -1 if v σ ∈ D π,φ .
The satisfaction degree sd(π, φ, σ) ∈ (0, +∞) of φ in π w.r.t. σ is

1 1 + vd(π, φ, σ)
The first case is the same as [START_REF] Rizk | Continuous valuations of temporal logic specifications with applications to parameter optimization and robustness measures[END_REF], i.e. the distance in [0, +∞) between v σ and the domain D π,φ . In the second case, we get a negative number related to the distance between v σ and the outside of the domain D π,φ . The satisfaction degree is defined as in [START_REF] Rizk | Continuous valuations of temporal logic specifications with applications to parameter optimization and robustness measures[END_REF] but now provides a notion of formula robustness when its value is greater than 1. Indeed sd(π, φ, σ) is bounded by the radius of D π,φ and describes in the space of the variables x 1 , . . . , x k by how much one can change the objective v σ while keeping σ(φ) satisfied on π.

In our implementation in BIOCHAM-4, the distance to, and penetration depth within, validity domains are not computed exactly, but approximated using the notions of generators and constraints in PPL. This is indeed sufficient to guide the search of parameter values and optimize formula robustness, using the Covariance Matrix Adaptive Evolutionary Strategy CMA-ES [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] as black-box continuous optimization tool with satisfaction degree as objective function.

Model Robustness and Parameter Sensitivity

The more classical notion of model robustness defined in [START_REF] Kitano | Towards a theory of biological robustness[END_REF] as the mean functionality with respect to a set of parameter perturbations P can be instanciated in our setting for a FO-LTL(R lin ) formula φ and some objective σ as

R P,φ,σ = p∈P min(1, sd(π(p), φ, σ))prob(p)dp
In BIOCHAM this integral is evaluated by sampling (log-)normally-distributed parameter values given by the user. This sampling stops when a user-given number of samples is reached, or when the relative sample-standard-deviation becomes smaller than a given threshold. However, the computation may be computationally expensive with tens of numerical simulations needed and is generally not usable inside a search for model parameters.

BIOCHAM commands

The bouncing ball is a classical example of the hybrid systems community with which we can illustrate our approach. One single simulation is already quite time consuming (about 2s) because of the checking of the bouncing event condition. computes a value (0.884 in 160s) that quantifies how robustly is satisfied the formula with an objective for h of 3.5 when parameters x 0 , K, D are perturbed (with default Gaussian distribution). On the other hand, the formula robustness necessitates a single simulation to get the value sd(π, φ, σ) = 1.4664. This low computational cost makes it possible to use formula robustness as an optimization criterion within the parameter search procedure. One can ask BIOCHAM-4 to search for parameter values such that the violation degree vd is below -0.5 (i.e., sd ≥ 2): It is worth noting however that, in theory, increasing the formula robustness does not necessarily improve the model robustness. Indeed, the formula robustness may increase with a parameter set that approaches a frontier where the formula is no longer satisfied, in which case the model robustness with respect to parameter perturbation may decrease. In many practical cases however, though the systems we tackle are highly non-linear, improving the formula robustness also did improve the model robustness.

Evaluation on MAPK Signaling Network

In [START_REF] Huang | Ultrasensitivity in the mitogen-activated protein kinase cascade[END_REF], the authors present the dynamics of a signaling network, namely the Mitogen Activated Protein Kinase (MAPK) cascade, with a model encompassing 22 species and fully relying on Mass-Action kinetics with 37 parameters. This model remains today a reference for the MAPK cascade as it properly describes the 3-level global structure and the 2-step nature of each level. Though it was supposed that explicit feedback reactions were necessary for obtaining oscillations in this model, it was demonstrated by Qiao et al. in 2007 [12] that this is actually not necessary. The authors explored blindly the space of 36 parameters and then established a bifurcation diagram in the 37th parameter in order to find oscillations.

The FO-LTL(R lin ) formula used in [START_REF] Rizk | Continuous valuations of temporal logic specifications with applications to parameter optimization and robustness measures[END_REF] to measure model robustness is

F(PPK ≥ up ∧ F(PPK ≤ up -amplitude))
, where PPK represents the output of the signaling cascade, i.e. the concentration of the doubly-phosphorylated kinase at the third level. By using that formula with an objective of for instance 0.5 for amplitude, one can actually search parameter values in a directed way, guided by continuous satisfaction degree. In order to ignore simple overshoots, one can consider the refined formula

F(PPK ≥ up ∧ F(PPK ≤ up -a 1 ∧ F(PPK ≥ up -a 1 + 0.5 * a 2 )))
, with objectives of 0.5 and 1 for respectively a 1 and a 2 . BIOCHAM-4 is able to find a solution for the 37 parameters in a few minutes (averaged to 5 minutes over 100 runs on a 3.6GHz Intel Core i7 machine) and to optimize formula robustness with no extra cost. This is also shown on a cell cycle model in the companion notebook mentionned above.

Conclusion

Building models that are robust to parameter variations is a key issue in systems biology and synthetic biology, because of both their fluctuating environment and the stochastic nature of biochemical reactions. BIOCHAM-4 is one of the very few tools to implement a metric notion of robustness for dynamical properties and optimize it by the parameter search procedure with no extra cost. Although not true in all generality for non-linear hybrid systems, optimizing formula robustness does improve in practice the statistical notion of model robustness with respect to parameter variations. The integration of formula robustness in the notion of satisfaction degree thus provides a simple and effective approach to the design of robust parameterization of reaction networks.

Model robustness to parameter perturbations is estimated by sampling the parameter space. This is time consuming since it involves making one simulation for each parameter set.

The FO-LTL formula used in this example defines the amplitude range of x in the domain of the free variable h The amplitude 3.96 exceeds here the objective of 3.5, the satisfaction degree is thus greater than one showing some formula robustness.

However the model robustness w.r.t. (default) parameter perturbations is below one (0.88) showing that some parameter perturbations destroy the amplitude objective.

biocham: validity_domain(F(x < h /\ F(x > h))). h<3.96646/\h> -1.0e-6 biocham: satisfaction_degree(F(x < h /\ F(x > h)), [h -> 3.5]). Let us search for parameter values exhibiting oscillations, using the same range of values as Qiao et al. [START_REF] Qiao | Bistability and oscillations in the huang-ferrell model of mapk signaling[END_REF] as given in their Table S1. We also add a range for E1 tot which was not searched by Qiao

  The FO-LTL(R lin ) formula F(x<h /\ F(x>h)) gives the height of the first bounce in the free variable h (∼ 4 in the companion notebook). Computing the validity domain with the command validity_domain(F(x<h / F(x>h))) and satisfaction degree satisfaction_degree(F(x<h / F(x>h)), [h -> 3.5]) on such a simulation trace is slightly faster (about 0.1s). Now, the command robustness(F(x<h / F(x>h)), [x0, K, D], [h->3.5]).

  search_parameters(F(x<h /\ F(x>h)), [5<=x0<=10, 0.5<=K<=1.0, -0.1<=D<=0.0], [h->3.5], cmaes_stop_fitness: -0.5). BIOCHAM-4 finds in 2.5s a solution with x 0 = 9.555, K = 0.801, D = -0.0585 such that vd = -0.572 and sd = 2.339. Computing again the model robustness as before we now get (in 101s) an improved model robustness of 0.930.

  (F(x < h /\ F(x > h)), [x0, K, D], [h -> 3.5]).Time: 158.027 s Robustness degree: 0.88414 With a satisfaction degree greater than 1 the formula is already satisfied with some margin. However, formula robustness can be further optimized by parameter optimization with formula robustness as objective function with no extra cost and much faster computation time than for estimating model robustness.It gives an amplitude of 4.84 which obviously improves the satisfaction degree of the formula. This is shown to also improve the model robustness (0.93) to parameter perturbations as expected in many examples. biocham: search_parameters( F(x < h /\ F(x > h)), [5 <= x0 <= 10, 0.5 <= K <= 1.0 , -0.1 <= D <= 0.0], [h -> 3.5], cmaes_stop_fitness: -0.5). Time: 27.491 s Stopping reason: Fitness: function value -5.72e-01 <= stopFitness Best satisfaction degree: 2.339136 [0] parameter(x0=9.555198250463423) [1] parameter(K=0.8014006370116868) [2] parameter(D= -0.058468364144163905) biocham: numerical_simulation. plot. biocham: validity_domain(F(x < h /\ F(x > h))). h<4.83916/\h> -1.0e-6 biocham: satisfaction_degree(F(x < h /\ F(x > h)), [h -> 3.5]). 2.339120 biocham: robustness(F(x < h /\ F(x > h)), [x0, K, D], [h -> 3.5]). Time: 100.736 s Robustness degree: 0.930547Appendix B Cell Cycle ExampleWe consider here the simple model of cell division cycle of Tyson[17] and the search of parameter values optimizing the robustness of some amplitude property of the oscillations. biocham: load(library:examples/cell_cycle/Tyson_1991.bc). biocham: option(time:150). biocham: parameter(k8=100). biocham: numerical_simulation. plot. biocham: satisfaction_degree(exists(max, exists(min, F(Cdc2-Cyclin~{p1} > max /\ F(Cdc2-Cyclin~{p1} < min /\ F(Cdc2-Cyclin~{p1} > max /\ F(Cdc2-Cyclin~{p1} < min)))) /\ max -min > amp)), [amp -> 0.19]). 0.974450 biocham: seed(0). biocham: robustness(exists(max, exists(min, F(Cdc2-Cyclin~{p1} > max /\ F(Cdc2-Cyclin~{p1} < min /\ F(Cdc2-Cyclin~{p1} > max /\ F(Cdc2-Cyclin~{p1} < min)))) /\ max -min > amp)), [k4, k6], [amp -> 0.19]). Time: 7.494 s

  et al. but used instead as parameter in a bifurcation diagram. biocham: timeout 600 Timeout set to 600 biocham: seed(0). biocham: search_parameters( exists(up, F(PP_K >= up /\ F(PP_K <= up -amplitude /\ F(PP_K >= up)))), [ 6.0e-4 <= KKK_tot <= 1.5e-2, 0.24 <= KK_tot <= 6, 0.24 <= K_tot <= 6, 6.0e-5 <= E2_tot <= 1.5e-3, 6.0e-5 <= KKPase_tot <= 1.5e-3, 2.4e-2 <= KPase_tot <= 0.6, biocham: seed(0). biocham: robustness( exists(up, F(PP_K >= up /\ F(PP_K <= up -amplitude /\ F(PP_K >= up)))), [ a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10 ], [amplitude -> 0.5]). Time: 10.38 s Robustness degree: 0.837798 biocham: seed(0). biocham: search_parameters( exists(up, F(PP_K >= up /\ F(PP_K <= up -amplitude /\ F(PP_K >= up)))), [ 6.0e-4 <= KKK_tot <= 1.5e-2,

http://lifeware.inria.fr/biocham4

https://lifeware.inria.fr/wiki/Main/Software#CMSB18

http://www.bugseng.com/ppl

<= a1 <= 5000, 200 <= a2 <= 5000, 200 <= a3 <= 5000, 200 <= a4 <= 5000, 200 <= a5 <= 5000, 200 <= a6 <= 5000, 200 <= a7 <= 5000, 200 <= a8 <= 5000, 200 <= a9 <= 5000, 200 <= a10 <= 5000, 30 <= d1 <= 750, 30 <= d2 <= 750, 30 <= d3 <= 750,

Acknowledgment: This work benefited from partial support from the ANR project HYCLOCK contract DS0401.

Appendix A Bouncing Ball Example

The bouncing ball with air friction is a classical example of the hybrid system community. It can be modeled in BIOCHAM by one first reaction with mass action (MA) law kinetics where the velocity catalyzes the position, a second reaction where the gravity force catalyzes the velocity, a third reaction for air friction where the velocity autocalyzes its degradation, and one event for elastic bouncing by reversing the sign of one parameter:

Appendix C MAPK Signaling Network Example

Here we consider the model of Huang and Ferrell of the Mitogen Activated Protein Kinase (MAPK) signaling network [START_REF] Huang | Ultrasensitivity in the mitogen-activated protein kinase cascade[END_REF]. We illustrate the search of parameter values exhibiting oscillations with some apmlitude properties and the optimization of the robustness of those properties. biocham: load(huang_ferrell.bc).


Jean-Charles Croix

Nicolas Durrande

Mauricio A Álvarez

Bayesian inversion of a diffusion evolution equation with application to Biology

A common task in experimental sciences is to fit mathematical models to real-world measurements to improve understanding of natural phenomenon (reverse-engineering or inverse modeling). When complex dynamical systems are considered, such as partial differential equations, this task may become challenging and ill-posed. In this work, a linear parabolic equation is considered where the objective is to estimate both the differential operator coefficients and the source term at once. The Bayesian methodology for inverse problems provides a form of regularization while quantifying uncertainty as the solution is a probability measure taking in account data. This posterior distribution, which is non-Gaussian and infinite dimensional, is then summarized through a mode and sampled using a state-of-the-art Markov-Chain Monte-Carlo algorithm based on a geometric approach. After a rigorous analysis, this methodology is applied on a dataset of the post-transcriptional regulation of Kni gap gene in the early development of Drosophila Melanogaster where mRNA concentration and both diffusion and depletion rates are inferred from noisy measurement of the protein concentration.

Introduction

The problem of diffusion in a porous media, which is ubiquitous in physics, engineering and biology, is usually represented by the following partial differential 1 equation: ∂y ∂t (x, t) + λ(x, t)y(x, t) -D(x, t)∆y(x, t) = f (x, t), ∀(x, t) ∈ Ω×]0, T ],

y(x, t) = 0, ∀(x, t) ∈ Ω × {t = 0}, y(x, t) = 0, ∀(x, t) ∈ ∂Ω×]0, T ].

(

) 1
where the spatial domain is Ω ⊂ R n (n ≤ 3) and the final time is T ∈ R + (other initial and boundary conditions are possible). In real world applications, the quantity of interest y (hereafter called the solution of equation 1) is typically the concentration of some chemical and evolves from a null initial state under three distinct mechanisms: a) direct variation in concentration, given by the source f , b) diffusion at rate D, c) production or depletion at a rate λ. Different hypotheses on the parameters lead to a well-posedness of this solution, which will be detailed later on. Besides the traditional computation of the solution from the parameters, one can use this model for the determination of an optimal control (e.g. source leading to the minimization of a particular cost functional) or the identification of parameters from partial knowledge of the solution in an inverse setting. This is the problem that will be of interest in this paper. The motivation comes from a a challenging identification problem in Biology where the objective is to infer jointly the decay rate λ, the diffusion rate D and the source f given a limited number of noisy observations of the solution y. Note that given their physical interpretation, the two rates and the solution must be positive so this constrain has to be taken into account in our inference scheme. This problem has already been solved by different approaches, under distinct sets of hypotheses. In [START_REF] Becker | Reverse-Engineering Post-Transcriptional Regulation of Gap Genes in Drosophilia melanogaster[END_REF], the authors use a system of ordinary differential equations instead of equation 1, and minimizes a discrete version of a leastsquare type functional, while confidence intervals on parameters are given by bootstrapping. Alternative methods to solve this problem, in a Bayesian setting, is to use Latent Force Models [START_REF] Alvarez | Latent Forces Models using Gaussian Processes[END_REF][START_REF] Särkkä | Gaussian Process Latent Force Models for Learning and Stochastic Control of Physical Systems[END_REF]. These approaches assume that unknown physical quantities can be modelled with Gaussian Processes [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]. In particular, if f is Gaussian and if the decay and the diffusion are constant, y is Gaussian as well (linear transformation of a Gaussian process). The two constant parameters λ and D can then be estimated through likelihood maximization. The main drawback of this approach is the difficulty to ensure positiveness of the source function f for instance.

In this work, we apply a more general methodology, based on the recent advent of Bayesian Inverse Problems [START_REF] Stuart | Inverse problems: A Bayesian perspective[END_REF] for infinite dimensional spaces. In a sense, it has the advantage of dealing with the ill-posedness while fully integrating the quantification of uncertainties. Moreover, the possibility of taking naturally physical constraints (such as positivity) will be particularly useful. This paper is organized as follows: section 2 presents all the mathematical analysis of the forward model (mapping equation parameters to its solution) and the Bayesian methodology applied in our particular setting. Section 3 focuses on a particular Markov chain Monte-Carlo algorithm that is adapted to the problem and robust to discretization. Finally, section 4 contains all implemen-tation details and the numerical results obtained on a dataset associated to the developmental biology of the Drosophila Melanogaster.

Bayesian inversion

As previously announced, our goal in this work is to infer a source term f (mRNA concentration) jointly with rates of diffusion D and decay λ (i.e. the parameter u = (λ, D, f )) from noisy and partial measurements of the solution y (gap protein concentration). This problem is ill-posed for multiple reasons: a) the parameter u is infinite dimensional and only finite data are available, b) the solution map is not injective and c) observations are noisy. The typical approach to alleviate this issue is to regularize the problem, usually adding constraints with Tikhonov-Philips functionals, to ensure uniqueness and continuity w.r.t. the observations [START_REF] Isakov | Inverse Problems for Partial Differential Equations[END_REF][START_REF] Schuster | Regularization Methods in Banach Spaces[END_REF]. Doing so, the regularized problem's solution will be compatible with the dataset, but possibly very different from the real parameter (if there is such thing). Additionally, a particularly valuable information is a representation of all parameters u that would lead to similar data, giving precise statement on how the dataset is informative [START_REF] Ghanem | Handbook of Uncertainty Quantification[END_REF][START_REF] Biegler | Large-Scale Inverse Problems and Quantification of Uncertainty[END_REF][START_REF] Sullivan | Introduction to uncertainty quantification[END_REF]. One approach consists in treating these 2 objectives sequentially, first regularizing then quantifying uncertainty. However, the Bayesian methodology for inverse problems ( [START_REF] Stuart | Inverse problems: A Bayesian perspective[END_REF] and more recently [START_REF] Dashti | The Bayesian Approach To Inverse Problems[END_REF]) is precisely tailored to complete both tasks at once in an elegant manner. One particularity of these recent contributions is to tackle inverse problems directly in function spaces, postponing discretization at the very end for implementation purposes, leading to algorithms robust to the discretization dimension. Indeed, finite approximations of probability measures may be absolutely continuous while their infinite counterparts are mutually singular. This become particularly problematic in MCMC sampling for instance [START_REF] Cotter | MCMC Methods for Functions: Modifying Old Algorithms to Make Them Faster[END_REF].

In essence, instead of searching for one particular parameter value that would solve the regularized problem, this approach considers the probability distribution of the parameter given the data. Namely, given a prior distribution (see subsection 2.2) and few technical conditions on the forward problem (parameter to data map, see subsection 2.1), Bayes theorem applies and exhibits a unique posterior distribution (see subsection 2.3), which is continuous in the data (w.r.t. Hellinger metric). Finally, one may summarize information from this posterior distribution such as expected value or modes (subsection 2.4). We detail in the next sections how this methodology can be applied to the problem at hand.

Forward model analysis

The first step is to detail precisely the required regularity of the solution map from equation 1. Using common variational techniques from PDE theory (see [START_REF] Evans | Partial differential equations[END_REF] or [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]), one can show that this equation has a unique weak solution (see theorem 2.1 which proof is in the appendix) given u = (λ, D, f ) in a domain U that will be precised later on. Moreover, this solution evolves smoothly when the parameter varies. Without loss of generality and keeping in mind the Biological application, the underlying physical domain will be Ω =]0, L[ with L ∈ R + . Theorem 2.1. Let P = R×]0, +∞[×L 2 ([0, T ], H -1 (Ω)), then for all u ∈ P, equation 1 has a unique weak solution, defining a map:

y : u ∈ P → y(u) ∈ W (0, T, L 2 , H 1 0
). Moreover, this map has the following properties:

1. it satisfies the following estimate ∀u ∈ P,

y(u) W (0,T,L 2 ,H 1 0 ) ≤ C √ D f L 2 ([0,T ],H -1 (Ω))
with C > 0 constant independent of u, 2. it is locally Lipschitz, ∀u ∈ P, ∀r > 0 such that B(u, r) ⊂ P, ∃L(r) > 0,

∀(u 1 , u 2 ) ∈ B(u, r), y(u 1 ) -y(u 2 ) W (0,T,L 2 ,H 1 0 ) ≤ L(r) u 1 -u 2 P 3. it is twice Fréchet differentiable on P.
The proof is given in the appendix. Let us now justify why the properties given in theorem 2.1 are important for the Bayesian inversion (most of them used in section 2.3):

1. The energy estimate will be critical to establish the continuity of the posterior w.r. In the rest of this work, we will restrict ourselves to the subset

U := R + ×]0, +∞[×C([0, T ] × [0, L], R),
which is implicitly equipped with the norm

u U = |λ| + |D| + f ∞ .
Since U ⊂ P, the solution map is well defined on this subset and keeps all its smoothness properties. Moreover, one can show that a weak solution of equation 1 for u ∈ U is also a strong solution [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF][START_REF] Evans | Partial differential equations[END_REF], but we don't need these regularity results here.

Choice of prior distribution

The second step is to choose a prior probability distribution on U, encoding all knowledge on the physics at hand, while being simple enough to keep analysis tractable. Here are the constraints given by the biological application:

• (λ, D) must be positive (which imposes decay, not production),

• f must be positive and continuous at all time and position (it is a concentration).

Starting with the depletion and diffusion parameters (λ, D), we choose respectively Borel prior distributions µ λ 0 and µ D 0 on R + with densities w.r.t. Lebesgue's measure. Now, since f must be positive, we re-parametrize the problem with

f * = exp(f ), (2)
where

f ∈ C([0, T ] × [0, L], R)
. By selecting a Borel probability measure µ f 0 on the Banach space C([0, T ] × [0, L], R), both continuity and positivity of the new source f * will be ensured almost-surely. In this work, we choose µ f 0 as a Gaussian distribution with covariance operator C with continuous realizations (see [START_REF] Igorevich | [END_REF] for a presentation of infinite dimensional Gaussian measures). Finally, we assume independence between the three components, leading to the following prior distribution:

µ 0 (du) := µ λ 0 (dλ) ⊗ µ D 0 (dD) ⊗ µ f 0 (df ). (3)
These choices clearly ensure positivity of u (in the previous sense) and µ 0 (U) = 1. The exponential map in equation 2 can be replaced with any sufficiently differentiable function from R to R + (to keep the second order Fréchet differentiability of the solution map). Alternative distributions are possible for f (Besov measure from [START_REF] Dashti | MAP estimators and their consistency in Bayesian nonparametric inverse problems[END_REF] or more general convex measures from [START_REF] Hosseini | Well-Posed Bayesian Inverse Problems: Priors with Exponential Tails[END_REF]) for the regularization. In practice however, our choice is also motivated by the fact that one can find a Gaussian reference measure µ ref in the form

µ ref = µ λ ref ⊗µ D ref ⊗µ f ref = N (λ ref , σ 2 λ )⊗N (D ref , σ 2 D )⊗N (0, C) = N (u ref , C ref )
where

u ref = (λ ref , D ref , 0) and
C ref : u ∈ R 2 × L 2 ([0, L]) → (σ 2 λ λ, σ 2 D D, Cf ) ∈ R 2 × L 2 ([0, L]), such that µ 0 << µ ref . Indeed, choose (λ ref , D ref ) ∈ R 2 and σ 2 λ , σ 2 D > 0 then µ 0 << µ ref with dµ 0 dµ ref (u) = dµ λ 0 dµ λ ref (λ) dµ D 0 dµ D ref (D).
This reference will be critical for modes analysis (section 2.4) and MCMC sampling (section 3.2).

Posterior distribution

The third and last step in the theoretical analysis of Bayesian inversion is to show that this particular setting (forward model and prior distribution) leads to a well defined posterior measure using Bayes theorem. This is the purpose of theorem 2.2 which is a direct application of the theory initially developed in [START_REF] Stuart | Inverse problems: A Bayesian perspective[END_REF]. Consider a dataset z = (z i ) i∈ [1,n] which corresponds to observations at different times and locations (t i , x i ) i∈ [1,n] and assume they are produced from the following model (in vector notations):

z = G(u) + η, (4)
where η ∼ N (0, σ 2 η I n ) (I n being the identity matrix of dimension n) and G : U → R n is the observation operator, mapping directly the PDE parameter u to the value of the associated solution y at measurement locations (composition of solution map y with Diracs). The following theorem, which is proved in the appendix, establishes the existence of a posterior probability measure µ z (the solution of our inverse problem), expressing how observations z changed prior beliefs µ 0 on the parameter u.

Theorem 2.2. Let G be the observation operator defined in equation 4 and µ 0 the probability measure defined in equation 3 (satisfying µ 0 (U) = 1), then there exists a unique posterior measure µ z for u|z. It is characterized by the following Radon-Nikodym density w.r.t. µ 0 :

dµ z dµ 0 (u) = 1 Z(z) exp (-φ(u; z)) , with φ(u; z) = 1 2σ 2 η z -G(u) 2 ,
and

Z(z) = U exp(-φ(u; z)µ 0 (dz).
Furthermore, the two following integrability conditions:

• E µ D 0 D -1 2 < +∞, • ∃κ > 0, E µ f 0 [exp(κ f ∞ )] < +∞,
imply the continuity of µ z in the data w.r.t. Hellinger distance.

Theorem 2.2 gives two distinct results: a) the existence and uniqueness of a posterior (as long as µ 0 is Radon and µ 0 (U) = 1), b) well-posedness of the Bayesian inverse problem under additional integrability conditions of certain functions. In particular, the need for an exponential moment under µ f 0 comes from the re-parametrization in equation 2. If one chooses a different map between f * and f , this condition may be considerably relaxed (using a polynomial map for instance).

Maximum a posteriori

In the previous section, we proved the well-posedness of the Bayesian inverse problem under specific integrability conditions. However, the posterior distribution is only known up to a multiplicative constant, through its density w.r.t. µ 0 . In our application, we will need to summarize µ z , which is usually done by the selection of elements of interest in U, such as the conditional mean (CM) or its modes. We saw in theorem 2.2 that the posterior expected value is continuous in the data (consequence of Hellinger continuity, see [START_REF] Dashti | The Bayesian Approach To Inverse Problems[END_REF]) but its optimality properties (under the frequentist paradigm) are not yet well-understood in infinite dimension to the best of our knowledge. This is why posterior modes (or Maximum a Posteriori) are more and more considered instead. Indeed, they provide a clear link with Tikhonov-Philips regularization (see [START_REF] Dashti | MAP estimators and their consistency in Bayesian nonparametric inverse problems[END_REF][START_REF] Helin | Maximum a posteriori probability estimates in infinite-dimensional Bayesian inverse problems[END_REF][START_REF] Agapiou | Sparsity-promoting and edge-preserving maximum a posteriori estimators in non-parametric Bayesian inverse problems[END_REF]) and a practical optimization problem (in case of Gaussian or Besov priors) which can be solved numerically, see theorem 2.3 (which proof is given in the appendix).

Theorem 2.3. Let µ 0 be the prior probability measure defined in equation 3 and µ ref the Gaussian reference measure from equation 2.2. Suppose additionally that

u ∈ U → ln dµ 0 dµ ref (u) ∈ R
is locally Lipschitz, then the modes of µ z are exactly the minimizers of the following (generalized) Onsager-Mashlup functional:

I(u) := Φ(u; z) + 1 2 u -u ref 2 µ ref -ln dµ 0 dµ ref (u) ,
where . µ ref and u ref are respectively the norm of the Cameron-Martin space and the mean of µ. A minimizer will be noted u MAP = (λ MAP , D MAP , f MAP ).

The precise application of this theorem to our biological setting is done in section 4.1.

Metropolis-Hastings algorithm

As it was previously announced, our motivation for the Bayesian methodology is the quantification of uncertainty, which will be done by simulation. Among the vast catalogue of methods for probability distributions sampling (Sequential Monte-Carlo, Approximate Bayesian Computations, Transport Maps, etc...), Markov chain Monte-Carlo is very popular (MCMC, see [START_REF] Brooks | Handbook of Markov Chain Monte Carlo[END_REF]) and well defined on function spaces [START_REF] Tierney | A note on Metropolis-Hastings kernels for general state spaces[END_REF] even though ergodicity analysis of such algorithm is still in its infancy [START_REF] Hairer | Analysis of SPDEs arising in path sampling part I: The Gaussian case[END_REF][START_REF] Hairer | Analysis of SPDEs arising in path sampling part II: The nonlinear case[END_REF][START_REF] Hairer | Spectral gaps for a Metropolis-Hastings algorithm in infinite dimensions[END_REF][START_REF] Rudolf | On a generalization of the preconditioned Crank-Nicolson Metropolis algorithm[END_REF]. After a short presentation of the functional Metropolis-Hastings algorithm (section 3.1), we will focus on a state-of-the-art Markov kernel designed to sample from Gaussian measures (section 3.2) and adapt it when the prior is not Gaussian, but absolutely continuous with respect to a Gaussian reference.

Metropolis-Hastings on function spaces

The Metropolis-Hastings algorithm (MH) is a very general [START_REF] Tierney | A note on Metropolis-Hastings kernels for general state spaces[END_REF] method to design Markov chains to sample from a given probability measure. It is based on a twostep process on each iteration: This algorithm provides a sample distributed under a predefined probability measure µ, if one selects α and Q in a specific way (see [START_REF] Tierney | A note on Metropolis-Hastings kernels for general state spaces[END_REF] for a discussion in general state spaces). For instance, let ν(du, dv) = µ(du)Q(u, dv) and ν T (du, dv) = µ(dv)Q(v, du), the Metropolis-Hastings algorithm typically considers the following acceptance probability:

α MH (u, v) = min 1, dν T dν (u, v) , (5)
which, in particular, requires the absolute continuity of ν T w.r.t. ν (detailed balance of the Markov chain). Contrary to finite dimensional situations, this condition may be difficult to satisfy and a common way to overcome this situation in Bayesian Inverse problems (see [START_REF] Dashti | The Bayesian Approach To Inverse Problems[END_REF][START_REF] Mark | Riemann manifold Langevin and Hamiltonian Monte Carlo methods[END_REF][START_REF] Beskos | Geometric MCMC for infinite-dimensional inverse problems[END_REF][START_REF] Cotter | MCMC Methods for Functions: Modifying Old Algorithms to Make Them Faster[END_REF][START_REF] Hairer | Spectral gaps for a Metropolis-Hastings algorithm in infinite dimensions[END_REF]) is to select Q revertible w.r.t. µ 0 . Indeed, in this case (with ν 0 (du, dv) = µ 0 (du)Q(u, dv)):

dν T dν (u, v) = dν T dν T 0 (u, v) dν dν0 (u, v) = dµz dµ0 (v) dµz dµ0 (u) = exp (φ(u; z) -φ(v; z)) . (6)
In theory, the MH algorithm may be implemented with a large family of proposal kernels Q. In practice however, they need to be as efficient as possible and thus adapted to the problem at hand. Two common desirable properties for Q are:

• Adjust the proposal to locally mimic the target distribution µ z ,

• Include a step size to tune acceptance probability.

These two properties may be used to trade-off self-correlation, acceptance rates and convergence speed to high interest areas of the parameter space. The next section presents an algorithm with both properties adapted to Gaussian priors which will then be applied to the problem at hand.

Geometric MCMC under Gaussian reference

du dt = - 1 2 K(u) C -1 ref (u -u ref ) + ∇ u φ(u; z) + K(u) dW dt , (7)
where K(u) is a (possibly position-dependent) preconditioner, W a cylindrical Brownian motion and ∇ u φ(u; y) the gradient in u of the negative log-likelihood.

According to [START_REF] Beskos | Geometric MCMC for infinite-dimensional inverse problems[END_REF], a semi-implicit discretization of equation 7 leads to a Markov chain with the following kernel:

Q(u, dv) = N ρ(u -u ref ) + u ref + 1 -ρ 2 √ h 2 g(u), K(u) , (8)
where h > 0 is a step-size parameter, ρ = 1-h 1+h and:

g(u) = -K(u) (C -1 ref -K(u) -1 )(u -u ref ) + ∇ u φ(u; z) .
This dynamic explores the parameter space with a balance between Newtontype descent to zones of high density and Gaussian exploration. The philosophy behind this kernel is to use alternative Gaussian reference measures locally adapted to the posterior distribution, since it has been recently showed that higher efficiency is obtained from operator weighted proposals ( [START_REF] Kody | Proposals which speed up function-space MCMC[END_REF] and later generalized in [START_REF] Beskos | Geometric MCMC for infinite-dimensional inverse problems[END_REF] and [START_REF] Cui | Dimensionindependent likelihood-informed MCMC[END_REF]). Indeed, highly informative datasets may result in a posterior measure significantly different from the prior in likelihood-informed directions and non-geometric kernels (such as Independent sampler or preconditioned Crank-Nicholson) are ineffective in this case. However, the infinite dimensional manifold Modified Adjusted Langevin Algorithm (∞-mMALA) considers a specific preconditioner:

K(u) = C -1 ref + H Φ (u) -1
, where H Φ (u) is the Gauss-Newton Hessian matrix of φ, which locally adapts to the posterior. This kernel does not preserve the distribution µ z but is shown to be absolutely continuous w.r.t. the reference measure µ ref , almost-surely in u (under technical assumptions regarding K(u)) and the Radon-Nikodym density is:

dQ(u, dv) dµ ref (v) = dN √ h 2 g(u), K(u) dN (0, C) v -ρ(u -u ref ) -u ref 1 -ρ 2 ,
and noting w =

v-ρ(u-u ref )-u ref √ 1-ρ 2
as it is done in [START_REF] Beskos | Geometric MCMC for infinite-dimensional inverse problems[END_REF], it finally comes:

dQ(u, dv) dµ ref (u, v) = exp - h 8 |K(u) -1 2 g(u)| 2 + √ h 2 K(u) -1 2 g(u), K(u) -1 2 w - 1 2 w, H Φ (u)w K 1 2 K(u) -1 2 .
Finally, the acceptance probability associated to the Markov kernel from equation 8 is

α(u, v) = min 1, dQ dµ ref (v, u) dµz dµ0 (v) dQ dµ ref (u, v) dµz dµ0 (u)
.

This algorithm is well-defined on function spaces (reversibility is ensured w.r.t. µ 0 ), thus it is robust to discretization as required. The ∞-mMALA proposal may be computationally expensive, as it requires to compute both gradient ∇Φ, Gauss-Newton Hessian H Φ and the Cholesky decomposition of K(u) -1 at each step. However, different dimension reduction techniques can be used (split in [START_REF] Beskos | Geometric MCMC for infinite-dimensional inverse problems[END_REF] or likelihood-informed in [START_REF] Cui | Dimensionindependent likelihood-informed MCMC[END_REF]) to reduce the computational burden. A second alternative is to choose a constant preconditioner, located at a posterior mode for instance (similar to HMALA in [START_REF] Cui | Dimensionindependent likelihood-informed MCMC[END_REF] and gpCN in [START_REF] Rudolf | On a generalization of the preconditioned Crank-Nicolson Metropolis algorithm[END_REF]).

Numerical application

We now turn to the practical implementation of the previous methodology on the problem of reverse-engineering for post-transcriptional gap-gene in Drosophila Melanogaster. First, we precise our choice of distributions for the parameters compatible with previous assumptions, and give a random series representation for f and precise the generalized Onsager-Mashlup functional. Then, we provide quantitative results on the dataset taken from [START_REF] Becker | Reverse-Engineering Post-Transcriptional Regulation of Gap Genes in Drosophilia melanogaster[END_REF], consisting in protein concentration measurements irregularly spread in space and time.

Choice of measures

We will now specify our choice of prior measure µ 0 with justifications:

• µ λ 0 , µ D 0 . Concerning the decay parameter, the only constraint given in the problem so far is positivity and Lebesgue density. However, the diffusion must also satisfy an integrability condition from theorem 2.2, which is clearly the case for uniform distributions. Finally, we choose

µ λ 0 = U([0, λ m ]) and µ D 0 = U([0, D m ]
), maximum parameters λ m and D m being tuned to 0.5 to be sufficiently large w.r.t. previous estimations from [START_REF] Becker | Reverse-Engineering Post-Transcriptional Regulation of Gap Genes in Drosophilia melanogaster[END_REF].

• µ f 0 . The prior measure µ f 0 will be chosen as a centred Gaussian measure on

L 2 ([0, T ] × [0, L]), with covariance operator on L 2 ([0, T ] × [0, L]): C = (-∆) -α ,
where α ∈ R + is a smoothness parameter tuned to ensure almost-surely continuity of the samples. The precise eigen-decomposition is given as follows i 1 , i 2 ≥ 1:

ϕ i1,i2 (x, t) = 1 √ LT sin i 1 π L x sin i 2 π T t , λ i1,i2 = i 1 πL 2 + i 2 πT 2 -α .
• With these choices done, the Radon-Nikodym density of the prior distribution w.r.t. the reference measure is

dµ 0 dµ ref (u) = 2πσ λ σ D λ m D m exp (λ -λ ref ) 2 2σ 2 λ + (λ -D ref ) 2 2σ 2 D ½ [0,λm] (λ)½ [0,Dm] (D),
which is locally Lipschitz (mean-value theorem). We now tune the parameters of µ ref , by simply choosing

λ ref = λ 2 , σ 2 λ = λ 2 m 12 , D ref = Dm 2 and σ 2 D = D 2 m 12
(minimizing Kullback-Leibler divergence). Finally, we are capable of specifying the exact form of the generalized Onsager-Mashlup functional using theorem 2.3:

I(u) = 1 2σ 2 η z -G(u) 2 + 1 2 f 2 µ f 0 .
In particular, parameters λ and D are only influenced by their range and contribution to the likelihood (the uniform prior is non-informative), contrary to the source f .

Prior and solution map discretization

The analysis conducted in all previous sections happens to be valid for infinite dimensional quantities. In practice however, one needs to discretize for numerical experiments. In this work, the solution map is approximated using finite elements in space (FEniCS library in Python, see [START_REF] Alnaes | The FEniCS Project Version 1.5[END_REF] and [START_REF] Petter | Solving PDEs in Python -The FEniCS Tutorial Volume I[END_REF]) and finite differences in time. We use 100 basis functions and 30 time steps on a desktop computer (Intel i7-3770 with 8Gb of RAM memory) 1 . We set L = 100 and T = 100 is the final time. Concerning the prior measure, we use a truncated Karhunen-Loeve basis of f to simulate from it

f = 1≤i1,i2≤N λ i1,i2 ξ i ϕ i1,i2 , (9)
where (ξ i1,i2 ) 1≤i1,i2≤N are i.i.d. N (0, 1) random variables. We consider N = 10 (100 basis functions) thus ũ = (λ, D, f ) is of dimension 102. All quantities related to negative log-likelihood derivatives (Gradient and Gauss-Newton Hessian matrix) are numerically computed using discrete adjoint methods (see [START_REF] Hinze | Optimization with PDE Constraints[END_REF] or [START_REF] Heinkenschloss | Numerical solution to implicitely constrained optimization problems[END_REF]) to keep scalability in N . The initial point in the chain is chosen at the MAP location, obtained by minimization of the functional in equation 4.1 (Prior based initialization results in long burnin phase). Practical optimization is done using L-BFGS-B algorithm from the Scipy library [START_REF] Byrd | A Limited Memory Algorithm for Bound Constrained Optimization[END_REF].

Results

We now turn to our main objective, the inversion and uncertainty quantification of gap-gene protein concentration from [START_REF] Becker | A Quantitative Study of Translational Regulation in Drosophilia Segment Determination[END_REF]. The dataset consists of 508 different measures which are non-uniformly spread in time and space (precise repartition can be seen in figure 4). The noise variance parameter is estimated using the following routine (10 iterations, 3 multi-start each):

1. Find u MAP minimizing I using the current noise level σ 2 η .

2. Update the current variance estimation

σ 2 η = 2 n-1 z -G(u MAP ) 2 . 3. Go back to 1.
The resulting estimated noise level is σ 2 η = 11, 77. With this estimated value, we compute our initial MAP estimate (numerical minimization of the Onsager-Machlup functional) and use it as initial point in the MCMC sampling. The Markov chain is ran for 11000 total iterations and the resulting traceplot is given in figure 1 for negative log-likelihood, decay, diffusion and first three components of f . The first thousand iterations are used as burnin and according to the hundred as posterior sample (thinning). From this, we compute both posterior mean and MAP estimates, the precise values of decay, diffusion, negative loglikelihood and Onsager-Machlup functional being given in table to the estimated values, one can also look at the marginal distribution on figure 3. Concerning the MAP estimator (figure 4), we recover both events described in [START_REF] Becker | A Quantitative Study of Translational Regulation in Drosophilia Segment Determination[END_REF] and [START_REF] Becker | Reverse-Engineering Post-Transcriptional Regulation of Gap Genes in Drosophilia melanogaster[END_REF], that is 2 pikes of protein concentrations. The first happens on the anterior part of the embryo in the early experiment (x = 35, t = 35). The second is much more intense and happens in the posterior part during the second half of the experiment. The estimated source explains these with an intense and localized increase in concentration. Finally, the uncertainty on both source and solution around data seems to be really low, which provides a good confidence on the level of mRNA at this time and part of the embryo (see figure 5). However, the point-wise variance on the solution y important before the first observations.

Conclusion

In this work, we applied the Bayesian inverse problem methodology from [START_REF] Stuart | Inverse problems: A Bayesian perspective[END_REF] to a practical Biological dynamical system. Doing so, we provide a rigorous and detailed analysis of the forward model, existence and continuity of the posterior measure, characterization of maximum a posteriori (MAP) estimates and state-of-the-art MCMC methodology. Because the forward MAP is non-linear, the unicity of posterior modes is unclear and it appears that local maximas are present. Nevertheless, the Bayesian methodology provides both a regularized solution to the problem, while giving a quantification of uncertainty. However, the estimation of prior hyper-parameters is still out of reach, giving poor confidence in the estimated variance. This direction requires further research, that we will try to address in a future work.

Aknowledgments

This work has been supported by Colciencias and ECOS-Nord under the grant C15M04. The authors would like to thank Dr. Xavier Bay for his precious advice and proof reading.

A Proofs

Proof of theorem 2.1. Using standard notations in PDE theory, let

Ω =]0, L[, H = L 2 (Ω), V = H 1 0 (Ω), V * = H -1 (Ω), P = R×]0, ∞[×L 2 ([0, T ], V * ) and W ([0, T ], H, V ) = y ∈ L 2 ([0, T ], H), y ′ ∈ L 2 ([0, T ], V * ) , equipped with the norm y W ([0,T ],H,V ) = y 2 L 2 ([0,T ],V ) + y ′ 2 L 2 ([0,T ],V * ) 1 2 ,
then the weak form associated to equation 1 has the following reduced form:

Find y ∈ W ([0, T ], H, V ), F (y, u), v = 0, ∀v ∈ L 2 ([0, T ], V ), with F (y, u), v = T 0 y t (t) -f (t), v(t) V * ,V + λ y(t), v(t) H + D y x (t), v x (t) H dt, = T 0 y t (t) -f (t), v(t) V * ,V + B(y(t), v(t), λ, D)dt.
• Existence of a solution map. The PDE operator is uniformly parabolic whenever D > 0, thus there exists a unique weak solution of equation 1 in W (0, T, H, V ) for every u ∈ P (see chapter 7 in [START_REF] Evans | Partial differential equations[END_REF] for instance). Moreover, we have the following energy estimate:

y(u) W (0,T,H,V ) ≤ C √ D f L 2 ([0,T ],V * )
with C a constant independent of u.

• Second order Fréchet differentiability of the solution map. Let u, h u ∈ P such that u + h u ∈ P, and h y ∈ W (0, T, H, V ) then ∀v ∈ L 2 ([0, T ], V ):

F (y + h y , u + h u ) -F (y, u), v = DF (y, u)[h y , h u ], v + c(h u , h y ) with |c(h u , h y )| = h λ T 0 h y , v H dt + h D T 0 h y x , v x H dt ≤ C v L 2 ([0,T ],V ) h 2 U ,
with C an other constant independent from u and:

DF (y, u)[h y , h u ], v = T 0 h y t , v V * ,V dt + T 0 h λ y, v H dt + T 0 h D y x , v x H dt + T 0 λ h y , v H dt + D T 0 h y x , v x V dt - T 0 f, v V * ,V dt.
Moreover, we have:

| DF (y, u)[h y , h u ], v | ≤ (h u , h y ) v
thus DF is bounded, which shows that F is Fréchet-differentiable on P.

Consider F y the partial derivative of F w.r.t. its first variable:

F y (y, u)h y , v = T 0 h y t , v + λ h y , v + D h y x , v x dt.
which defines a unique solution h ∈ W (0, T, H, V ) whenever D > 0 (using same arguments than previously). Here, F -1 y is clearly bounded. Because F is differentiable and F -1 y exists and is bounded, the implicit function theorem applies and y is differentiable on P. The second order differentiability uses the same arguments.

• Local Lipschitz continuity of the solution map. Let u ∈ P, r > 0 such that B(u, r) ⊂ P and (u 1 , u 2 ) ∈ B(u, r). There exists two unique solution with respect to (u 1 , u 2 ) satisfying ∀i ∈ {1, 2}, ∀v ∈ L 2 ([0, T ], V ):

y ′ i , v V * ,V + B(y i , v, λ i , D i ) = f i , v V * ,V
for almost every t ∈ [0, T ], which leads to ∀v ∈ L 2 ([0, T ], V ) and almost-every t ∈ [0, T ]:

y ′ 1 -y ′ 2 , v V * ,V + B(y 1 , v, λ 1 , D 1 ) -B(y 2 , v, λ 2 , D 2 ) = f 1 -f 2 , v V * ,V
. and equivalently:

y ′ 1 -y ′ 2 , v V * ,V + B(y 1 -y 2 , v, λ 1 , D 1 ) = f 1 -f 2 , v V * ,V + B(y 2 , v, λ 1 , D 1 ) -B(y 2 , v, λ 2 , D 2 ).
Now, letting v = y 1y 2 and using the coercivity of B we get:

y ′ 1 -y ′ 2 , y 1 -y 2 V * ,V + D 1 y 1 -y 2 2 V ≤ f 1 -f 2 , y 1 -y 2 V * ,V + B(y 2 , y 1 -y 2 , λ 1 , D 1 ) -B(y 2 , y 1 -y 2 , λ 2 , D 2 ). (10)
Now, we use the identity

d dt 1 2 y 1 -y 2 2 H = y ′ 1 -y ′ 2 , y 1 -y 2 V * ,V
and drop the positive term in the left hand side to get

d dt 1 2 y 1 -y 2 2 H ≤ f 1 -f 2 , y 1 -y 2 V * ,V + B(y 2 , y 1 -y 2 , λ 1 , D 1 ) -B(y 2 , y 1 -y 2 , λ 2 , D 2 ).
We integrate between 0 and t, use Cauchy-Schwartz and Poincarré's inequalities and y 1 (0) = y 2 (0) = 0 to obtain the following estimation

1 2 y 1 -y 2 2 H ≤ C y 1 -y 2 L 2 ([0,T ],V ) u 1 -u 2 , t -a.e.
The same reasoning applies now to equation 10 to get

y 1 -y 2 L 2 ([0,T ],V ) ≤ C u 1 -u 2 .
We start back from equation A to get:

y ′ 1 -y ′ 2 , v V * ,V = -B(y 1 -y 2 , v, λ 1 , D 1 ) + f 1 -f 2 , v V * ,V + B(y 2 , v, λ 1 , D 1 ) -B(y 2 , v, λ 2 , D 2 ),
and taking the supremum on the unit ball of V before integrating:

y ′ 1 -y ′ 2 L 2 ([0,T ],V * ) ≤ C 2 u 1 -u 2
which completes the proof.

Proof of theorem 2.2. Following the program in [START_REF] Dashti | The Bayesian Approach To Inverse Problems[END_REF], let us first show that µ z exists and is unique and then the continuity of the posterior with respect to the Hellinger distance.

• (Existence and unicity of µ z ). Let µ 0 be a probability measure defined as in equation 3 and consider the following Gaussian negative log-likelihood:

Φ(u; z) = 1 2σ 2 η z -G(u) 2 .
Since G is differentiable (theorem 2.1), Φ is measurable w.r.t. µ 0 for all z ∈ R n . Consider now the following integral Z(z) = U exp(-Φ(u; z))µ 0 (du).

The negative log-likelihood being positive, the integral is finite. Furthermore, we have ∀(u, z) ∈ U × R n , Φ(u; z) < +∞ thus Z(z) > 0 for every z ∈ R n . In other words, the following function exp(-Φ(u; z)) Z(z) defines a probability density function w.r.t. µ 0 , and the associated measure is µ z , the (unique) posterior distribution of u|z.

• (Continuity in z). It remains to show that the posterior distribution is continuous in z w.r.t. the Hellinger distance. Following the lines of [START_REF] Stuart | Inverse problems: A Bayesian perspective[END_REF][START_REF] Dashti | The Bayesian Approach To Inverse Problems[END_REF] we proove the following sufficient condition |Φ(u; z 1 ) -Φ(u; z 2 )| ≤ g(u) z 1z 2 , with g ∈ L 1 (µ 0 ).

|Φ(u; z 1 ) -Φ(u; z 2 )| = 1 2σ 2 ǫ z 1 2 -z 2 2 + 2 z 2 -z 1 , G(u) , = 1 2σ 2 ǫ |( z 1 -z 2 )( z 1 + z 2 ) + 2 z 2 -z 1 , G(u) | , ≤ 1 σ 2 η (r + G(u) ) z 1 -z 2 .
Now, using continuous injection between spaces, we have

G(u) ≤ n y(u) ∞ ≤ nC 1 y(u) W ≤ nC 2 f L 2 D ,
the last inequality is given by the energy estimate of theorem 2.1. It remains to see that

g(u) := 1 σ 2 η r + f D is in L 1 (µ 0 ), that is E µ0 [g(u)] = E µ D 0 1 D E µ f 0 [ f L 2 ].
The first term in the right hand side product is finite by assumption on µ D 0 . The second term is finite as well since

exp(f * ) L 2 ≤ exp(f * ) ∞ ≤ exp( f * ∞ )
which is µ f 0 integrable by the assumption on µ f 0 . Finally, Fubini's theorem gives the result on the integrability of g. Now, it remains to follow the proof from [START_REF] Stuart | Inverse problems: A Bayesian perspective[END_REF] as is.

Proof of theorem 2.3. The posterior measure µ z is absolutely continuous w.r.t. µ with the following Radon-Nikodym density: and since G(u 1 ) -G(u 2 ) ≤ C y(u 1 )y(u 2 ) , we conclude that Φ and thus Φ are locally Lipschitz. Now, we follow the lines of [START_REF] Dashti | MAP estimators and their consistency in Bayesian nonparametric inverse problems[END_REF] it comes: ≤ exp (-I(u 1 ) + I(u 2 )) .

J δ (u 1 ) J δ (u 2 ) = B δ (
A similar argument leads to lim inf δ→0 J δ (u 1 ) J δ (u 2 ) ≥ exp (-I(u 1 ) + I(u 2 )) .

We conclude that lim δ→0 J δ (u1)

J δ (u2) = exp (-I(u 1 ) + I(u 2 )). For a fixed value u 2 , this quantity is maximized when u 1 is a minimizer of I, the proof is then complete.

1 . 2 .

 12 Given a current state u ∈ U, propose a new candidate v according to a proposal Markov kernel Q(u, dv) (it is a probability distribution on U for almost any u ∈ U), Accept the new state v with probability α(u, v) or remain at u.

  We are now going to detail a specific Markov proposal kernel Q, tailored to sample distributions having a density w.r.t. a Gaussian measure µ ref . Most of recent work on infinite dimensional MCMC methods are based on the following Langevin stochastic differential equation:

2 Figure 1 :Figure 2 :

 212 Figure 1: Trace plots of Φ(u; z), λ, D, ξ 0 , ξ 1 and ξ 2 (the first 1000 iterations are burned).

Figure 3 :

 3 Figure 3: Marginal posterior distributions for λ|z and D|z.

Figure 4 :

 4 Figure 4: MAP estimator from MCMC sampling. Left: Estimated source term f MAP , right: estimated solution y(u MAP ) with absolute error at data locations. Grey bar represents the error level between data and y(u MAP ).

Figure 5 :

 5 Figure 5: Posterior point-wise variance of source (left figure) and solution (right figure).

2 ηz 2 ) + 2 z 2σ 2 ηG(u 1 )

 2221 exp(-Φ(u; z)).where Φ(u; z) = Φ(u; z)ln dµ0 dµ ref (u) . Let us now show that Φ is locally Lispchitz in its first argument:|Φ(u 1 ; z) -Φ(u 2 ; z)| = 1 2σ -G(u 1 ) 2z -G(u 2 ) 2 , 1 ) 2 -G(u 2 ) 2 + 2 z, G(u 2 ) -G(u 1 ) , ≤ G(u 1 ) + G(u -G(u 2 ) ,

  u1) exp(-Φ(u; z))µ ref (du)B δ (u2) exp(-Φ(v; z))µ ref (dv) , = B δ (u1) exp(-Φ(u; z) + Φ(u 1 ; z)) exp(-Φ(u 1 ; z))µ ref (du) B δ (u2) exp(-Φ(v; z) + Φ(u 2 ; z)) exp(-Φ(u 2 ; z))µ ref (dv). Now,J δ (u 1 ) J δ (u 2 ) ≤ exp δC -Φ(u 1 , z) + Φ(u 2 ; z) B δ (u1) µ ref (du)B δ (u2) µ ref (dv) and finally lim sup δ→0 J δ (u 1 ) J δ (u 2 )

Table 1 :

 1 10 -1 3.23 * 10 -2 246.72 460.73 Conditional mean (MCMC) 3.97 * 10 -1 3.30 * 10 -2 243.54 422.41 Values of decay, diffusion, negative log-likelihood and Onsager-Machlup functional for different estimators.

	1. Additionally

All codes are available online at https://github.com/JeanCharlesCroix/2018_Bayesian_estimation .